FV3 Bundle
external_ic_nlm.F90
Go to the documentation of this file.
1 module external_ic_nlm_mod
2 
3 #ifdef MAPL_MODE
4 #define DEALLOCGLOB_(A) if(associated(A)) then;A=0;call MAPL_DeAllocNodeArray(A,rc=STATUS);if(STATUS==MAPL_NoShm) deallocate(A,stat=STATUS);NULLIFY(A);endif
5 #endif
6 
7 #ifndef DYCORE_SOLO
8  use amip_interp_mod, only: i_sst, j_sst, sst_ncep
9 #endif
10  use fv_arrays_nlm_mod, only: real4, real8, fvprc
11  use fms_mod, only: file_exist, read_data, field_exist
12  use fms_io_mod, only: get_tile_string, field_size
13  use mpp_mod, only: mpp_error, fatal, note, mpp_broadcast,mpp_npes
14  use mpp_parameter_mod, only: agrid_param=>agrid
15  use mpp_domains_mod, only: mpp_get_tile_id, domain2d, mpp_update_domains, mpp_get_boundary, dgrid_ne
16  use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index
17  use field_manager_mod, only: model_atmos
18 
19 #ifndef MAPL_MODE
20  use constants_mod, only: pi=>pi_8, omega, grav, kappa, rdgas, rvgas, cp_air
21 #else
22  use mapl_constantsmod, only: mapl_pi_r8, mapl_omega, mapl_grav, &
23  mapl_kappa, mapl_rgas, mapl_rvap, &
24  mapl_cp
25 ! use MAPL_IOMod
26 ! use MAPL_ShmemMod
27 #endif
28  use, intrinsic :: iso_fortran_env, only: real64, real32
29 
30  use fv_arrays_nlm_mod, only: fv_atmos_type, fv_grid_type, fv_grid_bounds_type, fvprc, real4, real8
31  use fv_diagnostics_nlm_mod,only: prt_maxmin
32  use fv_mp_nlm_mod, only: is_master, ng, mp_barrier, mp_gather, mp_bcst, &
33  is,js,ie,je, isd,jsd,ied,jed, fill_corners, ydir
34  use fv_grid_utils_nlm_mod, only: ptop_min
35  use fv_grid_utils_nlm_mod, only: normalize_vect
36  use fv_mapz_nlm_mod, only: mappm
37  use fv_surf_map_nlm_mod, only: surfdrv
38  use fv_timing_nlm_mod, only: timing_on, timing_off
39  use init_hydro_nlm_mod, only: p_var
40 #ifndef MAPL_MODE
41  use fv_io_nlm_mod, only: fv_io_read_tracers
42  use sim_nc_nlm_mod, only: open_ncfile, close_ncfile, get_ncdim1, get_var1_double, get_var2_real, &
43  get_var3_r4
44 #endif
45  use fv_timing_nlm_mod, only: timing_on, timing_off
46  use mpp_mod, only: mpp_pe
47  use memutils_mod, only: print_memuse_stats
48  use fv_nwp_nudge_nlm_mod, only: t_is_tv
49 
50  implicit none
51 
52 #include "mpif.h"
53 
54  private
55 
56 #ifdef MAPL_MODE
57  real(FVPRC), parameter :: PI = mapl_pi_r8
58  real(FVPRC), parameter :: OMEGA = mapl_omega
59  real(FVPRC), parameter :: GRAV = mapl_grav
60  real(FVPRC), parameter :: KAPPA = mapl_kappa
61  real(FVPRC), parameter :: RDGAS = mapl_rgas
62  real(FVPRC), parameter :: RVGAS = mapl_rvap
63  real(FVPRC), parameter :: CP_AIR = mapl_cp
64 #endif
65 
66  real(FVPRC), parameter:: zvir = rvgas/rdgas - 1.
67  real(FVPRC) :: deg2rad
68 
69  public get_external_ic,get_cubed_sphere_terrain
70 
71  integer, SAVE :: tile, npes_x, npes_y
72 
73  integer :: status
74  integer :: iunit=15
75  integer :: ounit=16
76 
77 contains
78 
79  subroutine get_external_ic( Atm, fv_domain, use_geos_latlon_restart, use_geos_cubed_restart, ntracers )
80 
81  type(fv_atmos_type), intent(inout) :: atm(:)
82  type(domain2d), intent(inout) :: fv_domain
83  logical, optional, intent(in) :: use_geos_latlon_restart
84  logical, optional, intent(in) :: use_geos_cubed_restart
85  integer, optional, intent(in) :: ntracers(3)
86  real(FVPRC):: alpha = 0.
87  integer i,j,k,nq
88  logical :: cubed_sphere,fv_diag_ic
89 
90 ! * Initialize coriolis param:
91 
92  do j=jsd,jed+1
93  do i=isd,ied+1
94  atm(1)%gridstruct%fc(i,j) = 2.*omega*( -1.*cos(atm(1)%gridstruct%grid(i,j,1))*cos(atm(1)%gridstruct%grid(i,j,2))*sin(alpha) + &
95  sin(atm(1)%gridstruct%grid(i,j,2))*cos(alpha) )
96  enddo
97  enddo
98 
99  do j=jsd,jed
100  do i=isd,ied
101  atm(1)%gridstruct%f0(i,j) = 2.*omega*( -1.*cos(atm(1)%gridstruct%agrid(i,j,1))*cos(atm(1)%gridstruct%agrid(i,j,2))*sin(alpha) + &
102  sin(atm(1)%gridstruct%agrid(i,j,2))*cos(alpha) )
103  enddo
104  enddo
105 
106  cubed_sphere=.true.
107  fv_diag_ic=.false.
108  call mpp_update_domains( atm(1)%gridstruct%f0, atm(1)%domain )
109  if ( cubed_sphere ) call fill_corners(atm(1)%gridstruct%f0, atm(1)%npx, atm(1)%npy, ydir)
110 
111 ! Read in cubed_sphere terrain
112  if ( atm(1)%flagstruct%mountain ) then
113  call get_cubed_sphere_terrain(atm, fv_domain)
114  else
115  atm(1)%phis = 0.
116  endif
117 
118 ! Read in the specified external dataset and do all the needed transformation
119  if ( present(use_geos_latlon_restart)) then
120  if (allocated(atm(1)%q)) deallocate( atm(1)%q )
121  allocate ( atm(1)%q(isd:ied,jsd:jed,atm(1)%npz,atm(1)%ncnst) )
122  !call get_geos_latlon_ic( Atm, fv_domain, Atm(1)%ncnst, ntracers )
123  elseif ( present(use_geos_cubed_restart)) then
124  if (allocated(atm(1)%q)) deallocate( atm(1)%q )
125  allocate ( atm(1)%q(isd:ied,jsd:jed,atm(1)%npz,atm(1)%ncnst) )
126  !call get_geos_cubed_ic( Atm, fv_domain, Atm(1)%ncnst, ntracers )
127  else
128  if ( atm(1)%flagstruct%ncep_ic ) then
129  nq = 1
130  call timing_on('NCEP_IC')
131  call get_ncep_ic( atm, fv_domain, nq )
132  call timing_off('NCEP_IC')
133 #ifndef MAPL_MODE
134 #ifndef NO_FV_TRACERS
135  call fv_io_read_tracers( fv_domain, atm )
136  if(is_master()) write(6,*) 'All tracers except sphum replaced by FV IC'
137 #endif
138 #endif
139  elseif ( fv_diag_ic ) then
140 ! Interpolate/remap diagnostic output from a FV model diagnostic output file on uniform lat-lon A grid:
141  nq = size(atm(1)%q,4)
142 ! Needed variables: lon, lat, pfull(dim), zsurf, ps, ucomp, vcomp, w, temp, and all q
143 ! delz not implemnetd yet; set make_nh = .true.
144  call get_diag_ic( atm, fv_domain, nq )
145  else
146 ! The following is to read in lagacy lat-lon FV core restart file
147 ! is Atm%q defined in all cases?
148 
149  nq = size(atm(1)%q,4)
150  call get_fv_ic( atm, fv_domain, nq )
151  endif
152  endif
153 
154  call prt_maxmin('T', atm(1)%pt, is, ie, js, je, ng, atm(1)%npz, 1.0_fvprc)
155 
156  call p_var(atm(1)%npz, is, ie, js, je, atm(1)%ak(1), ptop_min, &
157  atm(1)%delp, atm(1)%delz, atm(1)%pt, atm(1)%ps, &
158  atm(1)%pe, atm(1)%peln, atm(1)%pk, atm(1)%pkz, &
159  kappa, atm(1)%q, ng, atm(1)%ncnst, dble(atm(1)%gridstruct%area),atm(1)%flagstruct%dry_mass, &
160  atm(1)%flagstruct%adjust_dry_mass, atm(1)%flagstruct%mountain, atm(1)%flagstruct%moist_phys, &
161  atm(1)%flagstruct%hydrostatic, atm(1)%flagstruct%nwat, atm(1)%domain,atm(1)%flagstruct%make_nh)
162 
163  end subroutine get_external_ic
164 
165 
166 
167  subroutine get_cubed_sphere_terrain( Atm, fv_domain )
168  type(fv_atmos_type), intent(inout) :: atm(:)
169  type(domain2d), intent(inout) :: fv_domain
170  integer :: ntileme
171  integer, allocatable :: tile_id(:)
172  character(len=64) :: fname
173  integer :: n
174  integer :: jbeg, jend
175  real(FVPRC) ftop
176 
177 #ifndef MAPL_MODE
178 
179 ! ntileMe = size(Atm(:)) ! This will have to be modified for mult tiles per PE
180 !! always one at this point
181 !
182 ! allocate( tile_id(ntileMe) )
183 ! tile_id = mpp_get_tile_id( fv_domain )
184 !
185 ! do n=1,ntileMe
186 !
187 ! call get_tile_string(fname, 'INPUT/fv_core.res.tile', tile_id(n), '.nc' )
188 !
189 ! if( file_exist(fname) ) then
190 ! call read_data(fname, 'phis', Atm(n)%phis(is:ie,js:je), &
191 ! domain=fv_domain, tile_count=n)
192 ! else
193 ! call surfdrv( Atm(n)%npx, Atm(n)%npy, grid, agrid, &
194 ! area, dx, dy, dxc, dyc, Atm(n)%phis, is_master())
195 ! call mpp_error(NOTE,'terrain datasets generated using USGS data')
196 ! endif
197 !
198 ! end do
199 !
200 ! call mpp_update_domains( Atm(1)%phis, domain )
201 ! ftop = g_sum(Atm(1)%phis(is:ie,js:je), is, ie, js, je, ng, area, 1)
202 !
203 ! call prt_maxmin('ZS', Atm(1)%phis, is, ie, js, je, ng, 1, 1./grav)
204 ! if(is_master()) write(6,*) 'mean terrain height (m)=', ftop/grav
205 !
206  deallocate( tile_id )
207 
208 #endif
209 
210  end subroutine get_cubed_sphere_terrain
211 
212 
213  subroutine get_diag_ic( Atm, fv_domain, nq )
214  type(fv_atmos_type), intent(inout) :: Atm(:)
215  type(domain2d), intent(inout) :: fv_domain
216  integer, intent(in):: nq
217 ! local:
218  character(len=128) :: fname, tracer_name
219  real(kind=4), allocatable:: wk1(:), wk2(:,:), wk3(:,:,:)
220  real(FVPRC), allocatable:: tp(:,:,:), qp(:,:,:,:)
221  real(FVPRC), allocatable:: ua(:,:,:), va(:,:,:), wa(:,:,:)
222  real(FVPRC), allocatable:: lat(:), lon(:), ak0(:), bk0(:)
223  real(FVPRC):: s2c(is:ie,js:je,4)
224  integer, dimension(is:ie,js:je):: id1, id2, jdc
225  real(FVPRC) psc(is:ie,js:je)
226  real(FVPRC) gzc(is:ie,js:je)
227  integer:: i, j, k, im, jm, km, npz, npt
228  integer:: i1, i2, j1, ncid
229  integer:: jbeg, jend
230  integer tsize(4), tr_ind
231  logical:: found
232 
233  integer sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
234 #ifndef MAPL_MODE
235 
236 ! deg2rad = pi/180.
237 !
238 ! npz = Atm(1)%npz
239 !
240 !! Zero out all initial tracer fields:
241 ! Atm(1)%q = 0.
242 !
243 ! fname = Atm(1)%res_latlon_dynamics
244 !
245 ! if( file_exist(fname) ) then
246 ! call open_ncfile( fname, ncid ) ! open the file
247 ! call get_ncdim1( ncid, 'lon', tsize(1) )
248 ! call get_ncdim1( ncid, 'lat', tsize(2) )
249 ! call get_ncdim1( ncid, 'pfull', tsize(3) )
250 !
251 ! im = tsize(1); jm = tsize(2); km = tsize(3)
252 !
253 ! if(is_master()) write(*,*) fname, ' FV_diag IC dimensions:', tsize
254 !
255 ! allocate ( lon(im) )
256 ! allocate ( lat(jm) )
257 !
258 ! call get_var1_double (ncid, 'lon', im, lon )
259 ! call get_var1_double (ncid, 'lat', jm, lat )
260 !
261 !! Convert to radian
262 ! do i=1,im
263 ! lon(i) = lon(i) * deg2rad ! lon(1) = 0.
264 ! enddo
265 ! do j=1,jm
266 ! lat(j) = lat(j) * deg2rad
267 ! enddo
268 !
269 ! allocate ( ak0(km+1) )
270 ! allocate ( bk0(km+1) )
271 !! npz
272 ! if ( npz /= km ) then
273 ! call mpp_error(FATAL,'==>Error in get_diag_ic: vertical dim must be the same')
274 ! else
275 ! ak0(:) = Atm(1)%ak(:)
276 ! bk0(:) = Atm(1)%bk(:)
277 ! endif
278 ! else
279 ! call mpp_error(FATAL,'==> Error in get_diag_ic: Expected file '//trim(fname)//' does not exist')
280 ! endif
281 !
282 !! Initialize lat-lon to Cubed bi-linear interpolation coeff:
283 ! call remap_coef( im, jm, lon, lat, id1, id2, jdc, s2c, Atm(1)%gridstruct%agrid, Atm(1)%bd )
284 !
285 !! Find bounding latitudes:
286 ! jbeg = jm-1; jend = 2
287 ! do j=js,je
288 ! do i=is,ie
289 ! j1 = jdc(i,j)
290 ! jbeg = min(jbeg, j1)
291 ! jend = max(jend, j1+1)
292 ! enddo
293 ! enddo
294 !
295 !! remap surface pressure and height:
296 ! allocate ( wk2(im,jbeg:jend) )
297 ! call get_var3_r4( ncid, 'ps', 1,im, jbeg,jend, 1,1, wk2 )
298 ! do j=js,je
299 ! do i=is,ie
300 ! i1 = id1(i,j)
301 ! i2 = id2(i,j)
302 ! j1 = jdc(i,j)
303 ! psc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
304 ! s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
305 ! enddo
306 ! enddo
307 !
308 ! call get_var3_r4( ncid, 'zsurf', 1,im, jbeg,jend, 1,1, wk2 )
309 ! do j=js,je
310 ! do i=is,ie
311 ! i1 = id1(i,j)
312 ! i2 = id2(i,j)
313 ! j1 = jdc(i,j)
314 ! gzc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
315 ! s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
316 ! enddo
317 ! enddo
318 ! deallocate ( wk2 )
319 !
320 !! Read in temperature:
321 ! allocate ( wk3(1:im,jbeg:jend, 1:km) )
322 ! call get_var3_r4( ncid, 'temp', 1,im, jbeg,jend, 1,km, wk3 )
323 ! allocate ( tp(is:ie,js:je,km) )
324 ! do k=1,km
325 ! do j=js,je
326 ! do i=is,ie
327 ! i1 = id1(i,j)
328 ! i2 = id2(i,j)
329 ! j1 = jdc(i,j)
330 ! tp(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
331 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
332 ! enddo
333 ! enddo
334 ! enddo
335 !
336 !! Read in all tracers:
337 ! allocate ( qp(is:ie,js:je,km,nq) )
338 ! qp = 0.
339 ! do tr_ind=1, nq
340 ! call get_tracer_names(MODEL_ATMOS, tr_ind, tracer_name)
341 ! if (field_exist(fname,tracer_name)) then
342 ! call get_var3_r4( ncid, tracer_name, 1,im, jbeg,jend, 1,km, wk3 )
343 ! do k=1,km
344 ! do j=js,je
345 ! do i=is,ie
346 ! i1 = id1(i,j)
347 ! i2 = id2(i,j)
348 ! j1 = jdc(i,j)
349 ! qp(i,j,k,tr_ind) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
350 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
351 ! enddo
352 ! enddo
353 ! enddo
354 ! endif
355 ! enddo
356 ! call remap_scalar(im, jm, km, npz, nq, nq, ak0, bk0, psc, gzc, tp, qp, Atm(1))
357 ! deallocate ( tp )
358 ! deallocate ( qp )
359 !
360 !! Winds:
361 ! call get_var3_r4( ncid, 'ucomp', 1,im, jbeg,jend, 1,km, wk3 )
362 ! allocate ( ua(is:ie,js:je,km) )
363 ! do k=1,km
364 ! do j=js,je
365 ! do i=is,ie
366 ! i1 = id1(i,j)
367 ! i2 = id2(i,j)
368 ! j1 = jdc(i,j)
369 ! ua(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
370 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
371 ! enddo
372 ! enddo
373 ! enddo
374 !
375 ! call get_var3_r4( ncid, 'vcomp', 1,im, jbeg,jend, 1,km, wk3 )
376 ! allocate ( va(is:ie,js:je,km) )
377 ! do k=1,km
378 ! do j=js,je
379 ! do i=is,ie
380 ! i1 = id1(i,j)
381 ! i2 = id2(i,j)
382 ! j1 = jdc(i,j)
383 ! va(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
384 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
385 ! enddo
386 ! enddo
387 ! enddo
388 ! call remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
389 !
390 ! deallocate ( ua )
391 ! deallocate ( va )
392 !
393 ! if ( .not. Atm(1)%hydrostatic ) then
394 ! if (field_exist(fname,'w')) then
395 ! allocate ( wa(is:ie,js:je,km) )
396 ! call get_var3_r4( ncid, 'w', 1,im, jbeg,jend, 1,km, wk3 )
397 ! do k=1,km
398 ! do j=js,je
399 ! do i=is,ie
400 ! i1 = id1(i,j)
401 ! i2 = id2(i,j)
402 ! j1 = jdc(i,j)
403 ! wa(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
404 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
405 ! enddo
406 ! enddo
407 ! enddo
408 ! call remap_wz(im, jm, km, npz, ng, ak0, bk0, psc, wa, Atm(1)%w, Atm)
409 ! deallocate ( wa )
410 ! else ! use "w = - fac * omega" ?
411 ! Atm(1)%w(:,:,:) = 0.
412 ! endif
413 !! delz:
414 ! if (field_exist(fname,'delz')) then
415 ! allocate ( wa(is:ie,js:je,km) )
416 ! call get_var3_r4( ncid, 'delz', 1,im, jbeg,jend, 1,km, wk3 )
417 ! do k=1,km
418 ! do j=js,je
419 ! do i=is,ie
420 ! i1 = id1(i,j)
421 ! i2 = id2(i,j)
422 ! j1 = jdc(i,j)
423 ! wa(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
424 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
425 ! enddo
426 ! enddo
427 ! enddo
428 ! call remap_wz(im, jm, km, npz, 0, ak0, bk0, psc, wa, Atm(1)%delz, Atm)
429 ! deallocate ( wa )
430 ! else ! Force make = T
431 ! Atm(1)%make_nh = .true.
432 ! endif
433 !
434 ! endif ! hydrostatic test
435 !
436 ! call close_ncfile ( ncid )
437 ! deallocate ( wk3 )
438 ! deallocate ( ak0 )
439 ! deallocate ( bk0 )
440 ! deallocate ( lat )
441 ! deallocate ( lon )
442 #endif
443 
444  end subroutine get_diag_ic
445 
446 ! subroutine get_geos_cubed_ic( Atm, fv_domain, nq, ntracers )
447 ! use GHOST_CUBSPH_mod, only : A_grid, ghost_cubsph_update
448 ! use CUB2CUB_mod, only: get_c2c_weight, &
449 ! interpolate_c2c
450 ! type(fv_atmos_type), intent(inout) :: Atm(:)
451 ! type(domain2d), intent(inout) :: fv_domain
452 ! integer, intent(in):: nq, ntracers(3)
453 !
454 ! character(len=128) :: fname, fname1
455 ! real(FVPRC), allocatable:: pkz0(:,:)
456 ! real(FVPRC), allocatable:: ps0(:,:), gz0(:,:), t0(:,:,:), q0(:,:,:)
457 ! real(FVPRC), allocatable:: u0(:,:,:), v0(:,:,:)
458 ! real(FVPRC), allocatable:: ak0(:), bk0(:)
459 ! integer :: i, j, k, l, iq, j1, j2, im, jm, km, npts, npx, npy, npz
460 ! integer :: iq_moist0 , iq_moist1
461 ! integer :: iq_gocart0, iq_gocart1
462 ! integer :: iq_pchem0 , iq_pchem1
463 ! integer :: ntiles=6
464 ! logical found
465 ! integer :: header(6)
466 ! character (len=8) :: imc, jmc
467 !
468 ! real(FVPRC) :: qtmp
469 !
470 ! integer:: i1, i2
471 ! integer:: ic, jc, jc0
472 ! real(FVPRC) psc(is:ie,js:je)
473 ! real(FVPRC) gzc(is:ie,js:je)
474 ! real(FVPRC), allocatable:: tp(:,:,:), qp(:,:,:,:)
475 ! real(FVPRC), allocatable:: ua(:,:,:), va(:,:,:)
476 !
477 ! real(REAL64), allocatable :: akbk_r8(:)
478 !
479 ! real(REAL64), dimension(:,:,:,:), allocatable :: corner_in, corner_out, weight_c2c
480 ! integer, dimension(:,:,:,:), allocatable :: index_c2c
481 !
482 ! integer :: is_i,ie_i, js_i,je_i
483 ! integer :: isd_i,ied_i, jsd_i,jed_i
484 ! integer :: ng_i = 5
485 ! type(domain2d) :: domain_i
486 !
487 ! real(FVPRC), allocatable :: ebuffer(:,:)
488 ! real(FVPRC), allocatable :: nbuffer(:,:)
489 ! real(FVPRC), allocatable :: wbuffer(:,:)
490 ! real(FVPRC), allocatable :: sbuffer(:,:)
491 !
492 ! integer (kind=MPI_OFFSET_KIND) :: slice_2d
493 ! integer (kind=MPI_OFFSET_KIND) :: offset
494 ! character(len=64) :: strTxt
495 ! integer :: nmoist, ngocart, npchem, nqmap
496 !
497 ! integer :: filetype
498 ! logical :: isNC4
499 ! type(MAPL_NCIO) :: ncio
500 ! integer :: nDims, nVars, ivar, dimSizes(3)
501 ! character(len=128) :: vname
502 ! real(FVPRC), allocatable :: gslice_r4(:,:)
503 ! real*8, allocatable :: gslice_r8(:,:)
504 ! integer :: tileoff,lvar_cnt
505 !
506 !!bma added
507 ! character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/)
508 !
509 ! iq_moist0 = 1
510 ! iq_moist1 = ntracers(1)
511 ! iq_gocart0=iq_moist1 +1
512 ! iq_gocart1=iq_moist1 +ntracers(2)
513 ! iq_pchem0 =iq_gocart1+1
514 ! iq_pchem1 =iq_gocart1+ntracers(3)
515 ! nmoist = ntracers(1)
516 ! ngocart = ntracers(2)
517 ! npchem = ntracers(3)
518 !
519 ! npx = Atm(1)%npx
520 ! npy = Atm(1)%npy
521 ! npz = Atm(1)%npz
522 !
523 !! Zero out all initial tracer fields:
524 ! Atm(1)%q = 0.
525 !
526 !! Read input FV core restart file
527 ! fname = "fvcore_internal_restart_in"
528 !
529 ! if( file_exist(fname) ) then
530 !
531 ! call MAPL_NCIOGetFileType(fname,filetype)
532 ! if (filetype >=0 ) then
533 ! isNC4 = .true.
534 ! else
535 ! isNC4 = .false.
536 ! end if
537 !
538 ! if (isNC4) then
539 !
540 ! NCIO = MAPL_NCIOOpen(fname)
541 ! call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km)
542 ! allocate(gslice_r8(im,jm))
543 !
544 ! else
545 !
546 ! open(IUNIT,file=fname ,access='sequential',form='unformatted',status='old')
547 ! read (IUNIT, IOSTAT=status) header
548 ! if (is_master()) print*, header
549 ! read (IUNIT, IOSTAT=status) header(1:5)
550 ! if (is_master()) print*, header(1:5)
551 !
552 ! im=header(1)
553 ! jm=header(2)
554 ! km=header(3)
555 !
556 ! end if
557 !
558 ! if(is_master()) write(*,*) 'Using GEOS restart:', fname
559 !
560 ! if ( file_exist(fname) ) then
561 ! if(is_master()) write(*,*) 'External IC dimensions:', im , jm , km
562 ! if(is_master()) write(*,*) 'Interpolating to :', npx-1, (npy-1)*6, npz
563 ! else
564 ! call mpp_error(FATAL,'==> Error from get_external_ic: &
565 ! & field not found')
566 ! endif
567 !
568 !!--------------------------------------------------------------------!
569 !! setup input cubed-sphere domain !
570 !!--------------------------------------------------------------------!
571 ! call mpp_domain_decomp(domain_i,im+1,(jm/ntiles)+1,ntiles,ng_i,ntiles, &
572 ! is_i,ie_i,js_i,je_i,isd_i,ied_i,jsd_i,jed_i,tile)
573 !!--------------------------------------------------------------------!
574 !! initialize cubed sphere grid: in !
575 !!--------------------------------------------------------------------!
576 ! allocate(corner_in(2,is_i:ie_i+1,js_i:je_i+1,tile:tile))
577 ! call init_cubsph_grid(im+1, is_i,ie_i, js_i,je_i, ntiles, corner_in)
578 ! call print_memuse_stats('get_geos_cubed_ic: init corner_in')
579 !!--------------------------------------------------------------------!
580 !! initialize cubed sphere grid: out !
581 !!--------------------------------------------------------------------!
582 ! allocate(corner_out(2,is:ie+1,js:je+1,tile:tile))
583 ! do l=tile,tile
584 ! do j=js,je+1
585 ! do i=is,ie+1
586 ! corner_out(1,i,j,l) = Atm(1)%gridstruct%grid(i,j,1)
587 ! corner_out(2,i,j,l) = Atm(1)%gridstruct%grid(i,j,2)
588 ! enddo
589 ! enddo
590 ! enddo
591 ! call print_memuse_stats('get_geos_cubed_ic: init corner_out')
592 !!--------------------------------------------------------------------!
593 !! calculate weights and indices from bilinear interpolation !
594 !! from grid_in to grid_out !
595 !!--------------------------------------------------------------------!
596 ! allocate(index_c2c (3, is:ie, js:je, tile:tile))
597 ! allocate(weight_c2c(4, is:ie, js:je, tile:tile))
598 ! call get_c2c_weight(ntiles, im+1, (jm/6)+1, &
599 ! is_i, ie_i, js_i, je_i, isd_i, ied_i, jsd_i, jed_i, &
600 ! corner_in(:,is_i:ie_i+1,js_i:je_i+1,tile), &
601 ! npx, npy, is,ie, js,je, tile,tile, &
602 ! corner_out(:,is:ie+1,js:je+1,tile:tile), &
603 ! index_c2c, weight_c2c, domain_i)
604 ! npts = im
605 ! call print_memuse_stats('get_geos_cubed_ic: get_c2c_weight')
606 !
607 ! allocate ( ak0(km+1) )
608 ! allocate ( bk0(km+1) )
609 ! allocate ( akbk_r8(km+1) )
610 ! if (isNC4) then
611 ! call MAPL_VarRead(NCIO,"AK",akbk_r8)
612 ! else
613 ! read (IUNIT, IOSTAT=status) akbk_r8
614 ! end if
615 ! ak0 = akbk_r8
616 ! if (isNC4) then
617 ! call MAPL_VarRead(NCIO,"BK",akbk_r8)
618 ! else
619 ! read (IUNIT, IOSTAT=status) akbk_r8
620 ! end if
621 ! bk0 = akbk_r8
622 ! deallocate ( akbk_r8 )
623 ! call print_memuse_stats('get_geos_cubed_ic: read ak/bk')
624 ! close (IUNIT)
625 !
626 !! Read U
627 ! allocate ( u0(isd_i:ied_i,jsd_i:jed_i+1,km) )
628 ! u0(:,:,:) = 0.0
629 ! if (isNC4) then
630 ! tileoff = (tile-1)*(jm/ntiles)
631 ! do k=1,km
632 ! call MAPL_VarRead(NCIO,"U",gslice_r8,lev=k)
633 ! u0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i)
634 ! enddo
635 ! else
636 !!offset = sequential access: 4 + INT(6) + 8 + INT(5) + 8 + DBL(NPZ+1) + 8 + DBL(NPZ+1) + 8
637 ! offset = 4 + 24 + 8 + 20 + 8 + (km+1)*8 + 8 + (km+1)*8 + 8
638 ! if (is_master()) print*, offset
639 ! call parallel_read_file_r8(fname, npts, is_i,ie_i, js_i,je_i, km, offset, u0(is_i:ie_i,js_i:je_i,:))
640 ! end if
641 ! call print_memuse_stats('get_geos_cubed_ic: read U')
642 !! Read V
643 ! allocate ( v0(isd_i:ied_i+1,jsd_i:jed_i,km) )
644 ! v0(:,:,:) = 0.0
645 ! if (isNC4) then
646 ! tileoff = (tile-1)*(jm/ntiles)
647 ! do k=1,km
648 ! call MAPL_VarRead(NCIO,"V",gslice_r8,lev=k)
649 ! v0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i)
650 ! enddo
651 ! else
652 ! if (is_master()) print*, offset
653 ! call parallel_read_file_r8(fname, npts, is_i,ie_i, js_i,je_i, km, offset, v0(is_i:ie_i,js_i:je_i,:))
654 ! end if
655 ! call print_memuse_stats('get_geos_cubed_ic: read V')
656 ! allocate ( sbuffer(is_i:ie_i,km) )
657 ! allocate ( wbuffer(js_i:je_i,km) )
658 ! allocate ( nbuffer(is_i:ie_i,km) )
659 ! allocate ( ebuffer(js_i:je_i,km) )
660 ! call mpp_get_boundary(u0, v0, domain_i, &
661 ! wbuffery=wbuffer, ebuffery=ebuffer, &
662 ! sbufferx=sbuffer, nbufferx=nbuffer, &
663 ! gridtype=DGRID_NE )
664 ! do k=1,km
665 ! do i=is_i,ie_i
666 ! u0(i,je_i+1,k) = nbuffer(i,k)
667 ! enddo
668 ! do j=js_i,je_i
669 ! v0(ie_i+1,j,k) = ebuffer(j,k)
670 ! enddo
671 ! enddo
672 ! deallocate ( sbuffer )
673 ! deallocate ( wbuffer )
674 ! deallocate ( nbuffer )
675 ! deallocate ( ebuffer )
676 ! call mpp_update_domains( u0, v0, domain_i, gridtype=DGRID_NE, complete=.true. )
677 ! call prt_maxmin(' U_geos', u0, is_i, ie_i, js_i, je_i, ng_i, km, 1.0_FVPRC)
678 ! call prt_maxmin(' V_geos', v0, is_i, ie_i, js_i, je_i, ng_i, km, 1.0_FVPRC)
679 ! allocate ( ua(is:ie,js:je,km) )
680 ! allocate ( va(is:ie,js:je,km) )
681 ! call interp_c2c_vect(npts, npts, npx-1, npy-1, km, ntiles, domain_i, &
682 ! is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i, &
683 ! u0, v0, ua, va, index_c2c, weight_c2c, corner_in(:,is_i:ie_i+1,js_i:je_i+1,tile), corner_out, Atm(1)%gridstruct)
684 ! deallocate ( v0 )
685 ! deallocate ( u0 )
686 ! deallocate ( corner_in )
687 ! deallocate ( corner_out )
688 !! Read T
689 ! allocate ( t0(isd_i:ied_i,jsd_i:jed_i,km) )
690 ! t0(:,:,:) = 0.0
691 ! if (isNC4) then
692 ! tileoff = (tile-1)*(jm/ntiles)
693 ! do k=1,km
694 ! call MAPL_VarRead(NCIO,"PT",gslice_r8,lev=k)
695 ! t0(is_i:ie_i,js_i:je_i,k) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i)
696 ! enddo
697 ! else
698 ! if (is_master()) print*, offset
699 ! call parallel_read_file_r8(fname, npts, is_i,ie_i, js_i,je_i, km, offset, t0(is_i:ie_i,js_i:je_i,:))
700 ! end if
701 ! call print_memuse_stats('get_geos_cubed_ic: read T')
702 !! Read PE at Surface only
703 ! allocate ( ps0(isd_i:ied_i,jsd_i:jed_i) )
704 ! ps0(:,:) = 0.0
705 ! if (isNC4) then
706 ! tileoff = (tile-1)*(jm/ntiles)
707 ! call MAPL_VarRead(NCIO,"PE",gslice_r8,lev=km+1)
708 ! ps0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i)
709 ! else
710 ! slice_2d = npts*npts*ntiles
711 ! do k=1,km
712 ! offset = offset + slice_2d*8 + 8 ! skip first KM levels of Edge Pressure to find surface pressure
713 ! enddo
714 ! if (is_master()) print*, offset, slice_2d
715 ! call parallel_read_file_r8(fname, npts, is_i,ie_i, js_i,je_i, 1, offset, ps0(is_i:ie_i,js_i:je_i))
716 ! end if
717 ! call mpp_update_domains(ps0, domain_i)
718 !! Read PKZ
719 ! allocate ( pkz0(isd_i:ied_i,jsd_i:jed_i) )
720 ! pkz0(:,:) = 0.0
721 ! if (isNC4) then
722 ! tileoff = (tile-1)*(jm/ntiles)
723 ! do k=1,km
724 ! call MAPL_VarRead(NCIO,"PKZ",gslice_r8,lev=k)
725 ! pkz0(is_i:ie_i,js_i:je_i) = gslice_r8(is_i:ie+i,tileoff+js_i:tileoff+je_i)
726 ! t0(is_i:ie_i,js_i:je_i,k) = t0(is_i:ie_i,js_i:je_i,k)*pkz0(is_i:ie_i,js_i:je_i)
727 ! enddo
728 ! else
729 ! if (is_master()) print*, offset
730 ! do k=1,km
731 ! call parallel_read_file_r8(fname, npts, is_i,ie_i, js_i,je_i, 1, offset, pkz0(is_i:ie_i,js_i:je_i))
732 !! t0 needs to be just temperature with no virtual effect
733 ! t0(is_i:ie_i,js_i:je_i,k) = t0(is_i:ie_i,js_i:je_i,k)*pkz0(is_i:ie_i,js_i:je_i)
734 ! enddo
735 ! end if
736 ! call print_memuse_stats('get_geos_cubed_ic: converted T')
737 ! deallocate ( pkz0 )
738 !
739 ! if (isNC4) then
740 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
741 ! deallocate(gslice_r8)
742 ! end if
743 !
744 ! allocate ( gz0(isd_i:ied_i,jsd_i:jed_i) )
745 ! gz0(:,:) = 0.0
746 !
747 ! write(imc, "(i8)") im
748 ! write(jmc, "(i8)") jm
749 ! imc = adjustl(imc)
750 ! jmc = adjustl(jmc)
751 !
752 ! write(fname1, "('topo_DYN_ave_',a,'x',a,'.data')") trim(imc), trim(jmc)
753 ! if (.not. file_exist(fname1)) then
754 ! call mpp_error(FATAL,'get_geos_cubed_ic: cannot find topo_DYN_ave file')
755 ! endif
756 ! call print_memuse_stats('get_geos_cubed_ic: '//TRIM(fname1)//' being read')
757 ! offset = 4
758 ! call parallel_read_file_r4(fname1, npts, is_i,ie_i, js_i,je_i, 1, offset, gz0(is_i:ie_i,js_i:je_i))
759 ! call mpp_update_domains(gz0, domain_i)
760 ! gz0 = gz0*grav
761 !
762 !! Read cubed-sphere phis from file since IMPORT is not ready yet
763 ! offset = 4
764 ! call parallel_read_file_r4('topo_dynave.data', Atm(1)%npx-1, is,ie, js,je, 1, offset, Atm(1)%phis(is:ie,js:je))
765 ! call mpp_update_domains(Atm(1)%phis, Atm(1)%domain)
766 ! Atm(1)%phis = Atm(1)%phis*grav
767 ! call print_memuse_stats('get_geos_cubed_ic: phis')
768 !
769 !! Horiz Interp for surface pressure
770 ! call prt_maxmin('PS_geos', ps0, is_i, ie_i, js_i, je_i, ng_i, 1, 1.0_FVPRC)
771 ! do j=js,je
772 ! do i=is,ie
773 ! ic=index_c2c(1,i,j,tile)
774 ! jc=index_c2c(2,i,j,tile)
775 ! psc(i,j)=weight_c2c(1,i,j,tile)*ps0(ic ,jc ) &
776 ! +weight_c2c(2,i,j,tile)*ps0(ic ,jc+1) &
777 ! +weight_c2c(3,i,j,tile)*ps0(ic+1,jc+1) &
778 ! +weight_c2c(4,i,j,tile)*ps0(ic+1,jc )
779 ! enddo
780 ! enddo
781 ! deallocate ( ps0 )
782 !! Horiz Interp for surface height
783 ! call prt_maxmin('GZ_geos', gz0, is_i, ie_i, js_i, je_i, ng_i, 1, 1.0/grav)
784 ! do j=js,je
785 ! do i=is,ie
786 ! ic=index_c2c(1,i,j,tile)
787 ! jc=index_c2c(2,i,j,tile)
788 ! gzc(i,j)=weight_c2c(1,i,j,tile)*gz0(ic ,jc ) &
789 ! +weight_c2c(2,i,j,tile)*gz0(ic ,jc+1) &
790 ! +weight_c2c(3,i,j,tile)*gz0(ic+1,jc+1) &
791 ! +weight_c2c(4,i,j,tile)*gz0(ic+1,jc )
792 ! enddo
793 ! enddo
794 ! deallocate ( gz0 )
795 !
796 !! Horiz Interp for Q
797 ! allocate ( q0(isd_i:ied_i,jsd_i:jed_i,km) )
798 ! allocate ( qp(is:ie,js:je,km,nq) )
799 ! q0(:,:,:) = 0.0
800 ! qp(:,:,:,:) = 0.0
801 !
802 !! Horiz Interp for moist tracers
803 !! is there a moist restart file to interpolate?
804 !! Read in tracers: only sphum at this point
805 ! if( file_exist("moist_internal_restart_in") .and. ntracers(1) > 0 ) then
806 ! if (is_master()) print*, 'Trying to interpolate moist_internal_restart_in'
807 !
808 ! call MAPL_NCIOGetFileType("moist_internal_restart_in",filetype)
809 !
810 ! if (filetype /= 0) then
811 ! offset=4
812 ! else
813 ! lvar_cnt = 0
814 ! allocate(gslice_r4(im,jm))
815 ! NCIO = MAPL_NCIOOpen("moist_internal_restart_in")
816 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
817 ! if (nVars /= iq_moist1-iq_moist0+1) call mpp_error(FATAL,'Wrong number of variables in moist file')
818 ! tileoff = (tile-1)*(jm/ntiles)
819 ! end if
820 !
821 ! do ivar=iq_moist0,iq_moist1
822 ! if (filetype ==0) lvar_cnt=lvar_cnt+1
823 ! do k=1,km
824 ! if (filetype /= 0) then
825 ! call parallel_read_file_r4('moist_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k))
826 ! call mpp_update_domains(q0(:,:,k), domain_i)
827 ! else
828 ! vname = trim(moist_order(lvar_cnt))
829 ! call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k)
830 ! q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i)
831 ! end if
832 ! call mpp_update_domains(q0(:,:,k), domain_i)
833 ! do j=js,je
834 ! do i=is,ie
835 ! ic=index_c2c(1,i,j,tile)
836 ! jc=index_c2c(2,i,j,tile)
837 ! qp(i,j,k,ivar)=weight_c2c(1,i,j,tile)*q0(ic ,jc ,k) &
838 ! +weight_c2c(2,i,j,tile)*q0(ic ,jc+1,k) &
839 ! +weight_c2c(3,i,j,tile)*q0(ic+1,jc+1,k) &
840 ! +weight_c2c(4,i,j,tile)*q0(ic+1,jc ,k)
841 ! enddo
842 ! enddo
843 !
844 ! enddo
845 ! if (ivar == 1) t0(is_i:ie_i,js_i:je_i,:) = (t0(is_i:ie_i,js_i:je_i,:)/(1.0 + zvir*q0(is_i:ie_i,js_i:je_i,:)))
846 ! call prt_maxmin( 'Q_geos_moist', q0, is_i, ie_i, js_i, je_i, ng_i, km, 1._FVPRC)
847 ! enddo
848 !
849 ! if (filetype == 0) then
850 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
851 ! deallocate(gslice_r4)
852 ! end if
853 !
854 ! end if
855 !
856 !! Horiz Interp for GOCART tracers
857 !! is there a gocart restart file to interpolate?
858 !! Read in tracers: only sphum at this point
859 ! if( file_exist("gocart_internal_restart_in") .and. ntracers(2) > 0 ) then
860 ! if (is_master()) print*, 'Trying to interpolate gocart_internal_restart_in'
861 !
862 ! call MAPL_NCIOGetFileType("gocart_internal_restart_in",filetype)
863 !
864 ! if (filetype /= 0) then
865 ! offset=4
866 ! else
867 ! lvar_cnt = 0
868 ! allocate(gslice_r4(im,jm))
869 ! NCIO = MAPL_NCIOOpen("gocart_internal_restart_in")
870 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
871 ! if (nVars /= iq_gocart1-iq_gocart0+1) call mpp_error(FATAL,'Wrong number of variables in gocart file')
872 ! tileoff = (tile-1)*(jm/ntiles)
873 ! end if
874 !
875 ! do ivar=iq_gocart0,iq_gocart1
876 ! if (filetype ==0) lvar_cnt=lvar_cnt+1
877 ! do k=1,km
878 ! if (filetype /= 0) then
879 ! call parallel_read_file_r4('gocart_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k))
880 ! call mpp_update_domains(q0(:,:,k), domain_i)
881 ! else
882 ! call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname)
883 ! call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k)
884 ! q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i)
885 ! end if
886 ! call mpp_update_domains(q0(:,:,k), domain_i)
887 ! do j=js,je
888 ! do i=is,ie
889 ! ic=index_c2c(1,i,j,tile)
890 ! jc=index_c2c(2,i,j,tile)
891 ! qp(i,j,k,ivar)=weight_c2c(1,i,j,tile)*q0(ic ,jc ,k) &
892 ! +weight_c2c(2,i,j,tile)*q0(ic ,jc+1,k) &
893 ! +weight_c2c(3,i,j,tile)*q0(ic+1,jc+1,k) &
894 ! +weight_c2c(4,i,j,tile)*q0(ic+1,jc ,k)
895 ! enddo
896 ! enddo
897 !
898 ! enddo
899 ! call prt_maxmin( 'Q_geos_gocart', q0, is_i, ie_i, js_i, je_i, ng_i, km, 1._FVPRC)
900 ! enddo
901 !
902 ! if (filetype == 0) then
903 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
904 ! deallocate(gslice_r4)
905 ! end if
906 !
907 ! end if
908 !! Horiz Interp for pchem tracers
909 !! is there a gocart restart file to interpolate?
910 !! Read in tracers: only sphum at this point
911 ! if( file_exist("pchem_internal_restart_in") .and. ntracers(3) > 0 ) then
912 ! if (is_master()) print*, 'Trying to interpolate pchem_internal_restart_in'
913 !
914 ! call MAPL_NCIOGetFileType("pchem_internal_restart_in",filetype)
915 !
916 ! if (filetype /= 0) then
917 ! offset=4
918 ! else
919 ! lvar_cnt = 0
920 ! allocate(gslice_r4(im,jm))
921 ! NCIO = MAPL_NCIOOpen("pchem_internal_restart_in")
922 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
923 ! if (nVars /= iq_pchem1-iq_pchem0+1) call mpp_error(FATAL,'Wrong number of variables in pchem file')
924 ! tileoff = (tile-1)*(jm/ntiles)
925 ! end if
926 !
927 ! do ivar=iq_pchem0,iq_pchem1
928 ! if (filetype == 0) lvar_cnt=lvar_cnt+1
929 ! do k=1,km
930 ! if (filetype /= 0) then
931 ! call parallel_read_file_r4('pchem_internal_restart_in', npts, is_i,ie_i, js_i,je_i, 1, offset, q0(is_i:ie_i,js_i:je_i,k))
932 ! call mpp_update_domains(q0(:,:,k), domain_i)
933 ! else
934 ! call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname)
935 ! call MAPL_VarRead(NCIO,vname,gslice_r4,lev=k)
936 ! q0(is_i:ie_i,js_i:je_i,k)=gslice_r4(is_i:ie+i,tileoff+js_i:tileoff+je_i)
937 ! end if
938 ! call mpp_update_domains(q0(:,:,k), domain_i)
939 ! do j=js,je
940 ! do i=is,ie
941 ! ic=index_c2c(1,i,j,tile)
942 ! jc=index_c2c(2,i,j,tile)
943 ! qp(i,j,k,ivar)=weight_c2c(1,i,j,tile)*q0(ic ,jc ,k) &
944 ! +weight_c2c(2,i,j,tile)*q0(ic ,jc+1,k) &
945 ! +weight_c2c(3,i,j,tile)*q0(ic+1,jc+1,k) &
946 ! +weight_c2c(4,i,j,tile)*q0(ic+1,jc ,k)
947 ! enddo
948 ! enddo
949 !
950 ! enddo
951 ! call prt_maxmin( 'Q_geos_pchem', q0, is_i, ie_i, js_i, je_i, ng_i, km, 1._FVPRC)
952 ! enddo
953 !
954 ! if (filetype == 0) then
955 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
956 ! deallocate(gslice_r4)
957 ! end if
958 !
959 ! end if
960 !
961 !! Horiz Interp for T
962 ! deallocate ( q0 )
963 ! call mpp_update_domains(t0, domain_i)
964 ! call prt_maxmin( 'T_geos', t0, is_i, ie_i, js_i, je_i, ng_i, km, 1.0_FVPRC)
965 ! allocate ( tp(is:ie,js:je,km) )
966 ! do k=1,km
967 ! do j=js,je
968 ! do i=is,ie
969 ! ic=index_c2c(1,i,j,tile)
970 ! jc=index_c2c(2,i,j,tile)
971 ! tp(i,j,k)=weight_c2c(1,i,j,tile)*t0(ic ,jc ,k) &
972 ! +weight_c2c(2,i,j,tile)*t0(ic ,jc+1,k) &
973 ! +weight_c2c(3,i,j,tile)*t0(ic+1,jc+1,k) &
974 ! +weight_c2c(4,i,j,tile)*t0(ic+1,jc ,k)
975 ! enddo
976 ! enddo
977 ! enddo
978 ! deallocate ( t0 )
979 ! deallocate( index_c2c )
980 ! deallocate( weight_c2c )
981 !
982 !! Horz/Vert remap for scalars
983 ! nqmap = nmoist + ngocart + npchem
984 !
985 ! call remap_scalar(im, jm, km, npz, nqmap, nqmap, ak0, bk0, psc, gzc, tp, qp, Atm(1))
986 !
987 ! deallocate ( tp )
988 ! deallocate ( qp )
989 ! call print_memuse_stats('get_geos_cubed_ic: remap_scalar')
990 !! Horz/Vert remap for U/V
991 ! call remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
992 ! deallocate ( ua )
993 ! deallocate ( va )
994 ! call print_memuse_stats('get_geos_cubed_ic: remap_winds')
995 !
996 ! else
997 ! call mpp_error(FATAL,'==> Error from get_external_ic: &
998 ! & Expected file '//trim(fname)//' does not exist')
999 ! endif
1000 !
1001 ! if (allocated(bk0)) deallocate ( bk0 )
1002 ! if (allocated(ak0)) deallocate ( ak0 )
1003 !
1004 !! Finished, let's check the results !
1005 !
1006 ! call prt_maxmin('GZ_model', Atm(1)%phis, is, ie, js, je, ng, 1, 1.0/grav)
1007 ! call prt_maxmin('PS_model', Atm(1)%ps, is, ie, js, je, ng, 1, 0.01_FVPRC)
1008 ! call prt_maxmin('DP_model', Atm(1)%delp, is, ie, js, je, ng, npz, 1.0_FVPRC)
1009 ! call prt_maxmin(' U_model', Atm(1)%u, is, ie, js, je, ng, npz, 1.0_FVPRC)
1010 ! call prt_maxmin(' V_model', Atm(1)%v, is, ie, js, je, ng, npz, 1.0_FVPRC)
1011 ! call prt_maxmin('PT_model', Atm(1)%pt, is, ie, js, je, ng, npz, 1.0_FVPRC)
1012 !! Range check the MOIST tracers
1013 ! do iq=iq_moist0,iq_moist1
1014 ! do k=1,npz
1015 ! do j=js,je
1016 ! do i=is,ie
1017 ! Atm(1)%q(i,j,k,iq) = MIN(Atm(1)%q(i,j,k,iq),1.d0)
1018 ! Atm(1)%q(i,j,k,iq) = MAX(Atm(1)%q(i,j,k,iq),0.d0)
1019 ! enddo
1020 ! enddo
1021 ! enddo
1022 ! enddo
1023 ! do iq=1,nmoist+npchem+ngocart
1024 ! call prt_maxmin('QP_model', Atm(1)%q(is:ie,js:je,1:npz,iq), is, ie, js, je, 0, npz, 1._FVPRC)
1025 ! enddo
1026 !
1027 ! end subroutine get_geos_cubed_ic
1028 !
1029 ! subroutine get_geos_latlon_ic( Atm, fv_domain, nq, ntracers)
1030 ! type(fv_atmos_type), intent(inout) :: Atm(:)
1031 ! type(domain2d), intent(inout) :: fv_domain
1032 ! integer, intent(in):: nq,ntracers(3)
1033 !
1034 ! character(len=128) :: fname, fname1
1035 ! real(FVPRC), allocatable:: pkz0(:,:)
1036 ! real(FVPRC), allocatable:: ps0(:,:), gz0(:,:), t0(:,:,:), q0(:,:,:)
1037 ! real(FVPRC), allocatable:: u0(:,:,:), v0(:,:,:), ua0(:,:), va0(:,:)
1038 ! real(FVPRC), allocatable:: lat(:), lon(:), ak0(:), bk0(:)
1039 ! integer :: i, j, k, l, iq, j1, j2, im, jm, km, npz
1040 ! logical found
1041 ! real(FVPRC) dak, dbk
1042 ! integer :: header(6)
1043 ! character (len=8) :: imc, jmc
1044 !
1045 ! integer:: i1, i2, nmoist, ngocart, npchem, nqmap, offset
1046 ! real(FVPRC):: s2c(is:ie,js:je,4)
1047 ! integer, dimension(is:ie,js:je):: id1, id2, jdc
1048 ! real(FVPRC) psc(is:ie,js:je)
1049 ! real(FVPRC) gzc(is:ie,js:je)
1050 ! real(FVPRC), allocatable:: tp(:,:,:), qp(:,:,:,:)
1051 ! real(FVPRC), allocatable:: ua(:,:,:), va(:,:,:)
1052 !
1053 ! real(REAL4), allocatable :: phis_r4(:,:)
1054 ! real(REAL64), allocatable :: r8latlon(:,:)
1055 ! real(REAL4), allocatable :: r4latlon(:,:)
1056 ! real(REAL64), allocatable :: akbk_r8(:)
1057 !
1058 ! character(len=128) :: tag
1059 ! character(len=3) :: cid
1060 !
1061 ! integer :: filetype
1062 ! logical :: isNC4
1063 ! type(MAPL_NCIO) :: ncio
1064 ! integer :: nDims, nVars, ivar, dimSizes(3)
1065 ! character(len=128) :: vname
1066 ! integer :: iq_moist0 , iq_moist1
1067 ! integer :: iq_gocart0, iq_gocart1
1068 ! integer :: iq_pchem0 , iq_pchem1
1069 ! integer :: lvar_cnt
1070 !!bma added
1071 ! character(len=128) :: moist_order(9) = (/"Q ","QLLS","QLCN","CLLS","CLCN","QILS","QICN","NCPL","NCPI"/)
1072 !
1073 ! iq_moist0 = 1
1074 ! iq_moist1 = ntracers(1)
1075 ! iq_gocart0=iq_moist1 +1
1076 ! iq_gocart1=iq_moist1 +ntracers(2)
1077 ! iq_pchem0 =iq_gocart1+1
1078 ! iq_pchem1 =iq_gocart1+ntracers(3)
1079 ! nmoist = ntracers(1)
1080 ! ngocart = ntracers(2)
1081 ! npchem = ntracers(3)
1082 !
1083 ! npz = Atm(1)%npz
1084 !
1085 !! Zero out all initial tracer fields:
1086 ! Atm(1)%q = 0.
1087 !
1088 !! Read in lat-lon FV core restart file
1089 ! fname = "fvcore_internal_restart_in"
1090 !
1091 ! if( file_exist(fname) ) then
1092 !
1093 !
1094 ! call MAPL_NCIOGetFileType(fname,filetype)
1095 ! if (filetype >=0 ) then
1096 ! isNC4 = .true.
1097 ! else
1098 ! isNC4 = .false.
1099 ! end if
1100 !
1101 ! if (isNC4) then
1102 !
1103 ! NCIO = MAPL_NCIOOpen(fname)
1104 ! call MAPL_NCIOGetDimSizes(NCIO,lon=im,lat=jm,lev=km)
1105 !
1106 ! else
1107 !
1108 ! open(IUNIT,file=fname ,access='sequential',form='unformatted',status='old')
1109 ! read (IUNIT, IOSTAT=status) header
1110 ! if (is_master()) print*, header
1111 ! read (IUNIT, IOSTAT=status) header(1:5)
1112 ! if (is_master()) print*, header(1:5)
1113 !
1114 ! im=header(1)
1115 ! jm=header(2)
1116 ! km=header(3)
1117 !
1118 ! end if
1119 !
1120 ! if(is_master()) write(*,*) 'Using GEOS restart:', fname
1121 ! if(is_master()) write(*,*) 'External IC dimensions:', im, jm, km
1122 !
1123 ! allocate ( lon(im) )
1124 ! do i=1,im
1125 ! lon(i) = (0.5 + real(i-1)) * 2.*pi/real(im)
1126 ! enddo
1127 ! allocate ( lat(jm) )
1128 ! do j=1,jm
1129 ! lat(j) = -0.5*pi + real(j-1)*pi/real(jm-1) ! SP to NP
1130 ! enddo
1131 !
1132 ! call remap_coef( im, jm, lon, lat, id1, id2, jdc, s2c , Atm(1)%gridstruct%agrid, Atm(1)%bd)
1133 !
1134 ! allocate ( ak0(km+1) )
1135 ! allocate ( bk0(km+1) )
1136 ! allocate ( akbk_r8(km+1) )
1137 ! if (isNC4) then
1138 ! call MAPL_VarRead(NCIO,"AK",akbk_r8)
1139 ! else
1140 ! read (IUNIT, IOSTAT=status) akbk_r8
1141 ! end if
1142 ! ak0 = akbk_r8
1143 ! if (isNC4) then
1144 ! call MAPL_VarRead(NCIO,"BK",akbk_r8)
1145 ! else
1146 ! read (IUNIT, IOSTAT=status) akbk_r8
1147 ! end if
1148 ! bk0 = akbk_r8
1149 ! deallocate ( akbk_r8 )
1150 !
1151 ! call print_memuse_stats('get_geos_latlon_ic: read ak/bk')
1152 !
1153 ! allocate ( r8latlon(im,jm) )
1154 !! Read U
1155 ! allocate ( u0(im,jm,km) )
1156 ! do k=1,km
1157 ! if (isNC4) then
1158 ! call MAPL_VarRead(NCIO,"U",r8latlon,lev=k)
1159 ! else
1160 ! read (IUNIT, IOSTAT=status) r8latlon
1161 ! end if
1162 !! Regrid from -180:180 to 0:360
1163 ! u0(1 :im/2,:,k) = r8latlon(im/2 + 1 :im , :)
1164 ! u0(im/2 + 1:im ,:,k) = r8latlon(1 :im/2, :)
1165 ! enddo
1166 ! call print_memuse_stats('get_geos_latlon_ic: read u')
1167 !! Read V
1168 ! allocate ( v0(im,jm,km) )
1169 ! do k=1,km
1170 ! if (isNC4) then
1171 ! call MAPL_VarRead(NCIO,"V",r8latlon,lev=k)
1172 ! else
1173 ! read (IUNIT, IOSTAT=status) r8latlon
1174 ! end if
1175 !! Regrid from -180:180 to 0:360
1176 ! v0(1 :im/2,:,k) = r8latlon(im/2 + 1 :im , :)
1177 ! v0(im/2 + 1:im ,:,k) = r8latlon(1 :im/2, :)
1178 ! enddo
1179 ! call print_memuse_stats('get_geos_latlon_ic: read v')
1180 ! if(is_master()) call pmaxmin( 'U_geos', u0(:,2:jm,:), im*(jm-1), km, 1.0_FVPRC)
1181 ! if(is_master()) call pmaxmin( 'V_geos', v0, im*jm , km, 1.0_FVPRC)
1182 ! allocate ( ua(is:ie,js:je,km) )
1183 ! allocate ( va(is:ie,js:je,km) )
1184 ! allocate ( ua0(im,jm) )
1185 ! allocate ( va0(im,jm) )
1186 ! do k=1,km
1187 !! Move latlon D winds to cell centers (A-grid)
1188 ! call d2a3d(u0(:,:,k), v0(:,:,k), ua0(:,:), va0(:,:), im, jm, 1, lon)
1189 !! Horiz Interp for U
1190 ! do j=js,je
1191 ! do i=is,ie
1192 ! i1 = id1(i,j)
1193 ! i2 = id2(i,j)
1194 ! j1 = jdc(i,j)
1195 ! ua(i,j,k) = s2c(i,j,1)*ua0(i1,j1 ) + s2c(i,j,2)*ua0(i2,j1 ) + &
1196 ! s2c(i,j,3)*ua0(i2,j1+1) + s2c(i,j,4)*ua0(i1,j1+1)
1197 ! enddo
1198 ! enddo
1199 !! Horiz Interp for V
1200 ! do j=js,je
1201 ! do i=is,ie
1202 ! i1 = id1(i,j)
1203 ! i2 = id2(i,j)
1204 ! j1 = jdc(i,j)
1205 ! va(i,j,k) = s2c(i,j,1)*va0(i1,j1 ) + s2c(i,j,2)*va0(i2,j1 ) + &
1206 ! s2c(i,j,3)*va0(i2,j1+1) + s2c(i,j,4)*va0(i1,j1+1)
1207 ! enddo
1208 ! enddo
1209 ! enddo
1210 ! call print_memuse_stats('get_geos_latlon_ic: d2a3d')
1211 ! deallocate ( v0 )
1212 ! deallocate ( u0 )
1213 ! deallocate ( ua0 )
1214 ! deallocate ( va0 )
1215 !! Read T
1216 ! allocate ( t0(im,jm,km) )
1217 ! do k=1,km
1218 ! if (isNC4) then
1219 ! call MAPL_VarRead(NCIO,"PT",r8latlon,lev=k)
1220 ! else
1221 ! read (IUNIT, IOSTAT=status) r8latlon
1222 ! end if
1223 !! Regrid from -180:180 to 0:360
1224 ! t0(1 :im/2,:,k) = r8latlon(im/2 + 1 :im , :)
1225 ! t0(im/2 + 1:im ,:,k) = r8latlon(1 :im/2, :)
1226 ! enddo
1227 ! call print_memuse_stats('get_geos_latlon_ic: read t')
1228 !! Read PE
1229 ! do k=1,km+1
1230 ! if (isNC4) then
1231 ! call MAPL_VarRead(NCIO,"PE",r8latlon,lev=k)
1232 ! else
1233 ! read (IUNIT, IOSTAT=status) r8latlon
1234 ! end if
1235 ! enddo
1236 !! Regrid from -180:180 to 0:360
1237 ! allocate ( ps0(im,jm) )
1238 ! ps0(1 :im/2,:) = r8latlon(im/2 + 1 :im , :)
1239 ! ps0(im/2 + 1:im ,:) = r8latlon(1 :im/2, :)
1240 ! allocate ( pkz0(im,jm) )
1241 ! do k=1,km
1242 ! if (isNC4) then
1243 ! call MAPL_VarRead(NCIO,"PKZ",r8latlon,lev=k)
1244 ! else
1245 ! read (IUNIT, IOSTAT=status) r8latlon
1246 ! end if
1247 !! Regrid from -180:180 to 0:360
1248 ! pkz0(1 :im/2,:) = r8latlon(im/2 + 1 :im , :)
1249 ! pkz0(im/2 + 1:im ,:) = r8latlon(1 :im/2, :)
1250 !! t0 needs to be just temperature with no virtual effect
1251 ! t0(:,:,k) = t0(:,:,k)*pkz0
1252 ! enddo
1253 ! deallocate ( r8latlon )
1254 ! call print_memuse_stats('get_geos_latlon_ic: converted T')
1255 ! deallocate ( pkz0 )
1256 ! if (isNC4) then
1257 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
1258 ! else
1259 ! close (IUNIT)
1260 ! end if
1261 !
1262 ! write(imc, "(i8)") im
1263 ! write(jmc, "(i8)") jm
1264 ! imc = adjustl(imc)
1265 ! jmc = adjustl(jmc)
1266 !
1267 ! write(fname1, "('topo_DYN_ave_',a,'x',a,'_DC.data')") trim(imc), trim(jmc)
1268 ! if (.not. file_exist(fname1)) then
1269 ! CALL mpp_error(FATAL,'get_geos_latlon_ic: cannot find topo_DYN_ave file')
1270 ! endif
1271 ! call print_memuse_stats('get_geos_latlon_ic: '//TRIM(fname1)//' being read')
1272 ! allocate ( r4latlon(im,jm) )
1273 ! open(IUNIT,file=fname1,form='unformatted',status='old')
1274 ! read(IUNIT) r4latlon
1275 ! close(IUNIT)
1276 !! Regrid from -180:180 to 0:360
1277 ! allocate ( gz0(im,jm) )
1278 ! gz0(1 :im/2,:) = r4latlon(im/2 + 1 :im , :)
1279 ! gz0(im/2 + 1:im ,:) = r4latlon(1 :im/2, :)
1280 ! gz0 = gz0*grav
1281 ! deallocate ( r4latlon )
1282 !
1283 !! Read cubed-sphere phis from file since IMPORT is not ready yet
1284 ! write(imc, "(i8)") Atm(1)%npx-1
1285 ! write(jmc, "(i8)") 6*(Atm(1)%npy-1)
1286 ! imc = adjustl(imc)
1287 ! jmc = adjustl(jmc)
1288 !
1289 ! write(fname1, "('topo_DYN_ave_',a,'x',a,'.data')") trim(imc), trim(jmc)
1290 ! if (.not. file_exist(fname1)) then
1291 ! call mpp_error(FATAL,'get_geos_latlon_ic: cannot find topo_DYN_ave file')
1292 ! endif
1293 ! allocate( phis_r4(Atm(1)%npx-1,6*(Atm(1)%npy-1)) )
1294 ! open(IUNIT,file=fname1,form='unformatted',status='old')
1295 ! read(IUNIT) phis_r4
1296 ! close(IUNIT)
1297 ! Atm(1)%phis(is:ie,js:je) = phis_r4(is:ie,js+(tile-1)*(Atm(1)%npy-1):je+(tile-1)*(Atm(1)%npy-1))*grav
1298 ! call mpp_update_domains(Atm(1)%phis, Atm(1)%domain)
1299 ! deallocate( phis_r4 )
1300 ! call print_memuse_stats('get_geos_latlon_ic: phis')
1301 !
1302 !! Horiz Interp for surface pressure
1303 ! if(is_master()) call pmaxmin( 'PS_geos', ps0, im, jm, 0.01_FVPRC)
1304 ! do j=js,je
1305 ! do i=is,ie
1306 ! i1 = id1(i,j)
1307 ! i2 = id2(i,j)
1308 ! j1 = jdc(i,j)
1309 ! psc(i,j) = s2c(i,j,1)*ps0(i1,j1 ) + s2c(i,j,2)*ps0(i2,j1 ) + &
1310 ! s2c(i,j,3)*ps0(i2,j1+1) + s2c(i,j,4)*ps0(i1,j1+1)
1311 ! enddo
1312 ! enddo
1313 ! deallocate ( ps0 )
1314 !! Horiz Interp for surface height
1315 ! if(is_master()) call pmaxmin( 'ZS_geos', gz0, im, jm, 1.0/grav)
1316 ! do j=js,je
1317 ! do i=is,ie
1318 ! i1 = id1(i,j)
1319 ! i2 = id2(i,j)
1320 ! j1 = jdc(i,j)
1321 ! gzc(i,j) = s2c(i,j,1)*gz0(i1,j1 ) + s2c(i,j,2)*gz0(i2,j1 ) + &
1322 ! s2c(i,j,3)*gz0(i2,j1+1) + s2c(i,j,4)*gz0(i1,j1+1)
1323 ! enddo
1324 ! enddo
1325 ! deallocate ( gz0 )
1326 !
1327 !! Horiz Interp for MOIST
1328 !! ----------------------
1329 ! allocate ( q0(im,jm,km) )
1330 ! allocate ( qp(is:ie,js:je,km,nq) )
1331 ! qp = 0.0
1332 !
1333 !! Horiz Interp for moist tracers
1334 !! is there a moist restart file to interpolate?
1335 !! Read in tracers: only sphum at this point
1336 ! if( file_exist("moist_internal_restart_in") .and. ntracers(1) > 0 ) then
1337 ! if (is_master()) print*, 'Trying to interpolate moist_internal_restart_in'
1338 ! allocate ( r4latlon(im,jm) )
1339 !
1340 ! call MAPL_NCIOGetFileType("moist_internal_restart_in",filetype)
1341 !
1342 ! if (filetype /= 0) then
1343 ! open(IUNIT,file="moist_internal_restart_in" ,access='sequential',form='unformatted',status='old')
1344 ! else
1345 ! lvar_cnt = 0
1346 ! NCIO = MAPL_NCIOOpen("moist_internal_restart_in")
1347 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
1348 ! if (nVars /= iq_moist1-iq_moist0+1) call mpp_error(FATAL,'Wrong number of variables in moist file')
1349 ! end if
1350 !
1351 ! do ivar=iq_moist0,iq_moist1
1352 ! if (filetype ==0) lvar_cnt=lvar_cnt+1
1353 ! do k=1,km
1354 ! if (filetype /= 0) then
1355 ! read (IUNIT, IOSTAT=status) r4latlon
1356 ! else
1357 ! vname = trim(moist_order(lvar_cnt))
1358 ! call MAPL_VarRead(NCIO,vname,r4latlon,lev=k)
1359 ! end if
1360 ! q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360
1361 ! q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360
1362 ! do j=js,je
1363 ! do i=is,ie
1364 ! i1 = id1(i,j)
1365 ! i2 = id2(i,j)
1366 ! j1 = jdc(i,j)
1367 ! qp(i,j,k,ivar) = s2c(i,j,1)*q0(i1,j1 ,k) + s2c(i,j,2)*q0(i2,j1 ,k) + &
1368 ! s2c(i,j,3)*q0(i2,j1+1,k) + s2c(i,j,4)*q0(i1,j1+1,k)
1369 ! enddo
1370 ! enddo
1371 ! enddo
1372 ! if (ivar == 1) t0 = (t0/(1.0 + zvir*q0(:,:,:)))
1373 ! if (is_master()) call pmaxmin( 'MOIST_Q_', q0(:,:,:), im*jm, km, 1.0_FVPRC)
1374 ! enddo
1375 ! if (filetype == 0) then
1376 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
1377 ! else
1378 ! close(IUNIT)
1379 ! end if
1380 ! deallocate(r4latlon)
1381 !
1382 ! end if
1383 !
1384 !! Horiz Interp for GOCART tracers
1385 !! is there a gocart restart file to interpolate?
1386 !! Read in tracers: only sphum at this point
1387 ! if( file_exist("gocart_internal_restart_in") .and. ntracers(2) > 0 ) then
1388 ! if (is_master()) print*, 'Trying to interpolate gocart_internal_restart_in'
1389 ! allocate ( r4latlon(im,jm) )
1390 ! call MAPL_NCIOGetFileType("gocart_internal_restart_in",filetype)
1391 !
1392 ! if (filetype /= 0) then
1393 ! open(IUNIT,file="gocart_internal_restart_in" ,access='sequential',form='unformatted',status='old')
1394 ! else
1395 ! lvar_cnt = 0
1396 ! NCIO = MAPL_NCIOOpen("gocart_internal_restart_in")
1397 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
1398 ! if (nVars /= iq_gocart1-iq_gocart0+1) call mpp_error(FATAL,'Wrong number of variables in gocart file')
1399 ! end if
1400 !
1401 ! do ivar=iq_gocart0,iq_gocart1
1402 ! if (filetype ==0) lvar_cnt=lvar_cnt+1
1403 ! do k=1,km
1404 ! if (filetype /= 0) then
1405 ! read (IUNIT, IOSTAT=status) r4latlon
1406 ! else
1407 ! call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname)
1408 ! call MAPL_VarRead(NCIO,vname,r4latlon,lev=k)
1409 ! end if
1410 ! q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360
1411 ! q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360
1412 ! do j=js,je
1413 ! do i=is,ie
1414 ! i1 = id1(i,j)
1415 ! i2 = id2(i,j)
1416 ! j1 = jdc(i,j)
1417 ! qp(i,j,k,ivar) = s2c(i,j,1)*q0(i1,j1 ,k) + s2c(i,j,2)*q0(i2,j1 ,k) + &
1418 ! s2c(i,j,3)*q0(i2,j1+1,k) + s2c(i,j,4)*q0(i1,j1+1,k)
1419 ! enddo
1420 ! enddo
1421 !
1422 ! enddo
1423 ! if (is_master()) call pmaxmin( 'GOCART_Q_', q0(:,:,:), im*jm, km, 1.0_FVPRC)
1424 ! enddo
1425 !
1426 ! if (filetype == 0) then
1427 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
1428 ! else
1429 ! close(IUNIT)
1430 ! end if
1431 ! deallocate(r4latlon)
1432 ! end if
1433 !
1434 !! Horiz Interp for pchem tracers
1435 !! is there a pchem restart file to interpolate?
1436 !! Read in tracers: only sphum at this point
1437 ! if( file_exist("pchem_internal_restart_in") .and. ntracers(3) > 0 ) then
1438 ! if (is_master()) print*, 'Trying to interpolate pchem_internal_restart_in'
1439 ! allocate ( r4latlon(im,jm) )
1440 ! call MAPL_NCIOGetFileType("pchem_internal_restart_in",filetype)
1441 !
1442 ! if (filetype /= 0) then
1443 ! open(IUNIT,file="pchem_internal_restart_in" ,access='sequential',form='unformatted',status='old')
1444 ! else
1445 ! lvar_cnt = 0
1446 ! NCIO = MAPL_NCIOOpen("pchem_internal_restart_in")
1447 ! call MAPL_NCIOGetDimSizes(NCIO,nVars=nVars)
1448 ! if (nVars /= iq_pchem1-iq_pchem0+1) call mpp_error(FATAL,'Wrong number of variables in pchem file')
1449 ! end if
1450 !
1451 ! do ivar=iq_pchem0,iq_pchem1
1452 ! if (filetype ==0) lvar_cnt=lvar_cnt+1
1453 ! do k=1,km
1454 ! if (filetype /= 0) then
1455 ! read (IUNIT, IOSTAT=status) r4latlon
1456 ! else
1457 ! call MAPL_NCIOGetVarName(NCIO,lvar_cnt,vname)
1458 ! call MAPL_VarRead(NCIO,vname,r4latlon,lev=k)
1459 ! end if
1460 ! q0(1 :im/2,:,k) = r4latlon(im/2 + 1 :im , :) ! Regrid from -180:180 to 0:360
1461 ! q0(im/2 + 1:im ,:,k) = r4latlon(1 :im/2, :) ! Regrid from -180:180 to 0:360
1462 ! do j=js,je
1463 ! do i=is,ie
1464 ! i1 = id1(i,j)
1465 ! i2 = id2(i,j)
1466 ! j1 = jdc(i,j)
1467 ! qp(i,j,k,ivar) = s2c(i,j,1)*q0(i1,j1 ,k) + s2c(i,j,2)*q0(i2,j1 ,k) + &
1468 ! s2c(i,j,3)*q0(i2,j1+1,k) + s2c(i,j,4)*q0(i1,j1+1,k)
1469 ! enddo
1470 ! enddo
1471 !
1472 ! enddo
1473 ! if (is_master()) call pmaxmin( 'PCHEM_Q_', q0(:,:,:), im*jm, km, 1.0_FVPRC)
1474 ! enddo
1475 !
1476 ! if (filetype == 0) then
1477 ! call MAPL_NCIOClose(NCIO,destroy=.true.)
1478 ! else
1479 ! close(IUNIT)
1480 ! end if
1481 ! deallocate(r4latlon)
1482 ! end if
1483 !
1484 ! call print_memuse_stats('get_geos_latlon_ic: remap_tracers')
1485 ! deallocate ( q0 )
1486 !
1487 !! Horiz Interp for T
1488 ! if(is_master()) call pmaxmin( 'T_geos', t0, im*jm, km, 1.0_FVPRC)
1489 ! allocate ( tp(is:ie,js:je,km) )
1490 ! do k=1,km
1491 ! do j=js,je
1492 ! do i=is,ie
1493 ! i1 = id1(i,j)
1494 ! i2 = id2(i,j)
1495 ! j1 = jdc(i,j)
1496 ! tp(i,j,k) = s2c(i,j,1)*t0(i1,j1 ,k) + s2c(i,j,2)*t0(i2,j1 ,k) + &
1497 ! s2c(i,j,3)*t0(i2,j1+1,k) + s2c(i,j,4)*t0(i1,j1+1,k)
1498 ! enddo
1499 ! enddo
1500 ! enddo
1501 ! deallocate ( t0 )
1502 ! call print_memuse_stats('get_geos_latlon_ic: remap_t')
1503 !
1504 !! Horz/Vert remap for MOIST, GOCART, and PCHEM scalars (Assuming Total Number is divisible by KM)
1505 !! -----------------------------------------------------------------------------------------------
1506 ! nqmap = nmoist + ngocart + npchem
1507 !
1508 ! call remap_scalar(im, jm, km, npz, nqmap, nqmap, ak0, bk0, psc, gzc, tp, qp, Atm(1))
1509 !
1510 ! deallocate ( tp )
1511 ! deallocate ( qp )
1512 ! call print_memuse_stats('get_geos_latlon_ic: remap_scalar')
1513 !
1514 !! Horz/Vert remap for U/V
1515 ! call remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
1516 ! deallocate ( ua )
1517 ! deallocate ( va )
1518 ! call print_memuse_stats('get_geos_latlon_ic: remap_winds')
1519 !
1520 ! else
1521 ! call mpp_error(FATAL,'==> Error from get_external_ic: &
1522 ! & Expected file '//trim(fname)//' does not exist')
1523 ! endif
1524 !
1525 ! if (allocated(bk0)) deallocate ( bk0 )
1526 ! if (allocated(ak0)) deallocate ( ak0 )
1527 ! if (allocated(lat)) deallocate ( lat )
1528 ! if (allocated(lon)) deallocate ( lon )
1529 !
1530 !! Finished, let's check the results !
1531 !
1532 ! call prt_maxmin('GZ_model', Atm(1)%phis, is, ie, js, je, ng, 1, 1.0/grav)
1533 ! call prt_maxmin('PS_model', Atm(1)%ps , is, ie, js, je, ng, 1, 0.01_FVPRC)
1534 ! call prt_maxmin('DP_model', Atm(1)%delp, is, ie, js, je, ng, npz, 1.0_FVPRC)
1535 ! call prt_maxmin(' U_model', Atm(1)%u , is, ie, js, je, ng, npz, 1.0_FVPRC)
1536 ! call prt_maxmin(' V_model', Atm(1)%v , is, ie, js, je, ng, npz, 1.0_FVPRC)
1537 ! call prt_maxmin('PT_model', Atm(1)%pt , is, ie, js, je, ng, npz, 1.0_FVPRC)
1538 !
1539 !! Range check the MOIST tracers
1540 ! do iq=iq_moist0,iq_moist1
1541 ! do k=1,npz
1542 ! do j=js,je
1543 ! do i=is,ie
1544 ! Atm(1)%q(i,j,k,iq) = MIN(Atm(1)%q(i,j,k,iq),1.d0)
1545 ! Atm(1)%q(i,j,k,iq) = MAX(Atm(1)%q(i,j,k,iq),0.d0)
1546 ! enddo
1547 ! enddo
1548 ! enddo
1549 ! enddo
1550 !
1551 ! if (is_master()) print*
1552 ! do iq=iq_moist0,iq_moist1
1553 ! call prt_maxmin('QP_MOIST_Q', Atm(1)%q(is:ie,js:je,1:npz,iq), is, ie, js, je, 0, npz, 1._FVPRC)
1554 ! enddo
1555 ! if (is_master()) print*
1556 ! do iq=iq_gocart0,iq_gocart1
1557 ! call prt_maxmin('QP_GOCART_Q', Atm(1)%q(is:ie,js:je,1:npz,iq), is, ie, js, je, 0, npz, 1._FVPRC)
1558 ! enddo
1559 ! if (is_master()) print*
1560 ! do iq=iq_pchem0,iq_pchem1
1561 ! call prt_maxmin('QP_PCHEM_Q', Atm(1)%q(is:ie,js:je,1:npz,iq), is, ie, js, je, 0, npz, 1._FVPRC)
1562 ! enddo
1563 !
1564 ! end subroutine get_geos_latlon_ic
1565 
1566  subroutine get_ncep_ic( Atm, fv_domain, nq )
1567  type(fv_atmos_type), intent(inout) :: Atm(:)
1568  type(domain2d), intent(inout) :: fv_domain
1569  integer, intent(in):: nq
1570 ! local:
1571  character(len=128) :: fname
1572  real(kind=4), allocatable:: wk1(:), wk2(:,:), wk3(:,:,:)
1573  real(FVPRC), allocatable:: tp(:,:,:), qp(:,:,:)
1574  real(FVPRC), allocatable:: ua(:,:,:), va(:,:,:)
1575  real(FVPRC), allocatable:: lat(:), lon(:), ak0(:), bk0(:)
1576  real(FVPRC):: s2c(is:ie,js:je,4)
1577  integer, dimension(is:ie,js:je):: id1, id2, jdc
1578  real(FVPRC) psc(is:ie,js:je)
1579  real(FVPRC) gzc(is:ie,js:je)
1580  real(FVPRC) tmean
1581  integer:: i, j, k, im, jm, km, npz, npt
1582  integer:: i1, i2, j1, ncid
1583  integer:: jbeg, jend
1584  integer tsize(4)
1585  logical:: read_ts = .true.
1586  logical:: land_ts = .false.
1587  logical:: found
1588 #ifndef MAPL_MODE
1589 
1590 ! deg2rad = pi/180.
1591 !
1592 ! npz = Atm(1)%npz
1593 !
1594 !! Zero out all initial tracer fields:
1595 ! Atm(1)%q = 0.
1596 !
1597 ! fname = Atm(1)%res_latlon_dynamics
1598 !
1599 ! if( file_exist(fname) ) then
1600 ! call open_ncfile( fname, ncid ) ! open the file
1601 ! call get_ncdim1( ncid, 'lon', tsize(1) )
1602 ! call get_ncdim1( ncid, 'lat', tsize(2) )
1603 ! call get_ncdim1( ncid, 'lev', tsize(3) )
1604 !
1605 ! im = tsize(1); jm = tsize(2); km = tsize(3)
1606 !
1607 ! if(is_master()) write(*,*) fname, ' NCEP IC dimensions:', tsize
1608 !
1609 ! allocate ( lon(im) )
1610 ! allocate ( lat(jm) )
1611 !
1612 ! call get_var1_double (ncid, 'lon', im, lon )
1613 ! call get_var1_double (ncid, 'lat', jm, lat )
1614 !
1615 !! Convert to radian
1616 ! do i=1,im
1617 ! lon(i) = lon(i) * deg2rad ! lon(1) = 0.
1618 ! enddo
1619 ! do j=1,jm
1620 ! lat(j) = lat(j) * deg2rad
1621 ! enddo
1622 !
1623 ! allocate ( ak0(km+1) )
1624 ! allocate ( bk0(km+1) )
1625 ! call get_var1_double (ncid, 'hyai', km+1, ak0, found )
1626 ! if ( .not. found ) ak0(:) = 0.
1627 !
1628 ! call get_var1_double (ncid, 'hybi', km+1, bk0 )
1629 !
1630 !! Note: definition of NCEP hybrid is p(k) = a(k)*1.E5 + b(k)*ps
1631 ! ak0(:) = ak0(:) * 1.E5
1632 !
1633 !! Limiter to prevent NAN at top during remapping
1634 ! if ( bk0(1) < 1.E-9 ) ak0(1) = max(1.e-9, ak0(1))
1635 ! else
1636 ! call mpp_error(FATAL,'==> Error in get_external_ic: Expected file '//trim(fname)//' does not exist')
1637 ! endif
1638 !
1639 !! Initialize lat-lon to Cubed bi-linear interpolation coeff:
1640 ! call remap_coef( im, jm, lon, lat, id1, id2, jdc, s2c, Atm(1)%gridstruct%agrid, Atm(1)%bd )
1641 !
1642 !! Find bounding latitudes:
1643 ! jbeg = jm-1; jend = 2
1644 ! do j=js,je
1645 ! do i=is,ie
1646 ! j1 = jdc(i,j)
1647 ! jbeg = min(jbeg, j1)
1648 ! jend = max(jend, j1+1)
1649 ! enddo
1650 ! enddo
1651 !
1652 !! remap surface pressure and height:
1653 !
1654 ! allocate ( wk2(im,jbeg:jend) )
1655 ! call get_var3_r4( ncid, 'PS', 1,im, jbeg,jend, 1,1, wk2 )
1656 !
1657 ! do j=js,je
1658 ! do i=is,ie
1659 ! i1 = id1(i,j)
1660 ! i2 = id2(i,j)
1661 ! j1 = jdc(i,j)
1662 ! psc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1663 ! s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1664 ! enddo
1665 ! enddo
1666 !
1667 ! call get_var3_r4( ncid, 'PHIS', 1,im, jbeg,jend, 1,1, wk2 )
1668 ! do j=js,je
1669 ! do i=is,ie
1670 ! i1 = id1(i,j)
1671 ! i2 = id2(i,j)
1672 ! j1 = jdc(i,j)
1673 ! gzc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1674 ! s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1675 ! enddo
1676 ! enddo
1677 !
1678 ! deallocate ( wk2 )
1679 ! allocate ( wk2(im,jm) )
1680 !
1681 ! if ( read_ts ) then ! read skin temperature; could be used for SST
1682 !
1683 ! call get_var2_real( ncid, 'TS', im, jm, wk2 )
1684 !
1685 ! if ( .not. land_ts ) then
1686 ! allocate ( wk1(im) )
1687 !
1688 ! do j=1,jm
1689 !! Read NCEP ORO (1; land; 0: ocean; 2: sea_ice)
1690 ! call get_var3_r4( ncid, 'ORO', 1,im, j,j, 1,1, wk1 )
1691 ! tmean = 0.
1692 ! npt = 0
1693 ! do i=1,im
1694 ! if( abs(wk1(i)-1.) > 0.99 ) then ! ocean or sea ice
1695 ! tmean = tmean + wk2(i,j)
1696 ! npt = npt + 1
1697 ! endif
1698 ! enddo
1699 !!------------------------------------------------------
1700 !! Replace TS over interior land with zonal mean SST/Ice
1701 !!------------------------------------------------------
1702 ! if ( npt /= 0 ) then
1703 ! tmean= tmean / real(npt)
1704 ! do i=1,im
1705 ! if( abs(wk1(i)-1.) <= 0.99 ) then ! Land points
1706 ! if ( i==1 ) then
1707 ! i1 = im; i2 = 2
1708 ! elseif ( i==im ) then
1709 ! i1 = im-1; i2 = 1
1710 ! else
1711 ! i1 = i-1; i2 = i+1
1712 ! endif
1713 ! if ( abs(wk1(i2)-1.)>0.99 ) then ! east side has priority
1714 ! wk2(i,j) = wk2(i2,j)
1715 ! elseif ( abs(wk1(i1)-1.)>0.99 ) then ! west side
1716 ! wk2(i,j) = wk2(i1,j)
1717 ! else
1718 ! wk2(i,j) = tmean
1719 ! endif
1720 ! endif
1721 ! enddo
1722 ! endif
1723 ! enddo ! j-loop
1724 ! deallocate ( wk1 )
1725 ! endif !(.not.land_ts)
1726 !
1727 ! do j=js,je
1728 ! do i=is,ie
1729 ! i1 = id1(i,j)
1730 ! i2 = id2(i,j)
1731 ! j1 = jdc(i,j)
1732 ! Atm(1)%ts(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1733 ! s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1734 ! enddo
1735 ! enddo
1736 ! call prt_maxmin('SST_model', Atm(1)%ts, is, ie, js, je, 0, 1, 1._FVPRC)
1737 !
1738 !! Perform interp to FMS SST format/grid
1739 !#ifndef DYCORE_SOLO
1740 ! call ncep2fms(im, jm, lon, lat, wk2)
1741 ! if( is_master() ) then
1742 ! write(*,*) 'External_ic_mod: i_sst=', i_sst, ' j_sst=', j_sst
1743 ! call pmaxmin( 'SST_ncep_fms', sst_ncep, i_sst, j_sst, 1._FVPRC)
1744 ! endif
1745 !#endif
1746 ! endif !(read_ts)
1747 !
1748 ! deallocate ( wk2 )
1749 !
1750 !! Read in temperature:
1751 ! allocate ( wk3(1:im,jbeg:jend, 1:km) )
1752 ! call get_var3_r4( ncid, 'T', 1,im, jbeg,jend, 1,km, wk3 )
1753 !
1754 ! allocate ( tp(is:ie,js:je,km) )
1755 ! do k=1,km
1756 ! do j=js,je
1757 ! do i=is,ie
1758 ! i1 = id1(i,j)
1759 ! i2 = id2(i,j)
1760 ! j1 = jdc(i,j)
1761 ! tp(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1762 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1763 ! enddo
1764 ! enddo
1765 ! enddo
1766 !
1767 !! Read in tracers: only sphum at this point
1768 ! call get_var3_r4( ncid, 'Q', 1,im, jbeg,jend, 1,km, wk3 )
1769 !
1770 ! allocate ( qp(is:ie,js:je,km) )
1771 ! do k=1,km
1772 ! do j=js,je
1773 ! do i=is,ie
1774 ! i1 = id1(i,j)
1775 ! i2 = id2(i,j)
1776 ! j1 = jdc(i,j)
1777 ! qp(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1778 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1779 ! enddo
1780 ! enddo
1781 ! enddo
1782 !
1783 ! call remap_scalar(im, jm, km, npz, nq, nq, ak0, bk0, psc, gzc, tp, qp, Atm)
1784 ! deallocate ( tp )
1785 ! deallocate ( qp )
1786 !
1787 !! Winds:
1788 ! call get_var3_r4( ncid, 'U', 1,im, jbeg,jend, 1,km, wk3 )
1789 !
1790 ! allocate ( ua(is:ie,js:je,km) )
1791 ! do k=1,km
1792 ! do j=js,je
1793 ! do i=is,ie
1794 ! i1 = id1(i,j)
1795 ! i2 = id2(i,j)
1796 ! j1 = jdc(i,j)
1797 ! ua(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1798 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1799 ! enddo
1800 ! enddo
1801 ! enddo
1802 !
1803 ! call get_var3_r4( ncid, 'V', 1,im, jbeg,jend, 1,km, wk3 )
1804 ! call close_ncfile ( ncid )
1805 !
1806 ! allocate ( va(is:ie,js:je,km) )
1807 ! do k=1,km
1808 ! do j=js,je
1809 ! do i=is,ie
1810 ! i1 = id1(i,j)
1811 ! i2 = id2(i,j)
1812 ! j1 = jdc(i,j)
1813 ! va(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1814 ! s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1815 ! enddo
1816 ! enddo
1817 ! enddo
1818 ! deallocate ( wk3 )
1819 !
1820 ! call remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
1821 !
1822 ! deallocate ( ua )
1823 ! deallocate ( va )
1824 !
1825 ! deallocate ( ak0 )
1826 ! deallocate ( bk0 )
1827 ! deallocate ( lat )
1828 ! deallocate ( lon )
1829 #endif
1830 
1831  end subroutine get_ncep_ic
1832 
1833 
1834 
1835  subroutine get_fv_ic( Atm, fv_domain, nq )
1836  type(fv_atmos_type), intent(inout) :: Atm(:)
1837  type(domain2d), intent(inout) :: fv_domain
1838  integer, intent(in):: nq
1839 
1840  character(len=128) :: fname, tracer_name
1841  real(FVPRC), allocatable:: ps0(:,:), gz0(:,:), u0(:,:,:), v0(:,:,:), t0(:,:,:), dp0(:,:,:), q0(:,:,:,:)
1842  real(FVPRC), allocatable:: ua(:,:,:), va(:,:,:)
1843  real(FVPRC), allocatable:: lat(:), lon(:), ak0(:), bk0(:)
1844  integer :: i, j, k, im, jm, km, npz, tr_ind
1845  integer tsize(4)
1846 ! integer sphum, liq_wat, ice_wat, cld_amt ! GFDL AM2 physics
1847  logical found
1848 
1849 #ifndef MAPL_MODE
1850 ! npz = Atm(1)%npz
1851 !
1852 !! Zero out all initial tracer fields:
1853 ! Atm(1)%q = 0.
1854 !
1855 !! Read in lat-lon FV core restart file
1856 ! fname = Atm(1)%res_latlon_dynamics
1857 !
1858 ! if( file_exist(fname) ) then
1859 ! call field_size(fname, 'T', tsize, field_found=found)
1860 ! if(is_master()) write(*,*) 'Using lat-lon FV restart:', fname
1861 !
1862 ! if ( found ) then
1863 ! im = tsize(1); jm = tsize(2); km = tsize(3)
1864 ! if(is_master()) write(*,*) 'External IC dimensions:', tsize
1865 ! else
1866 ! call mpp_error(FATAL,'==> Error in get_external_ic: field not found')
1867 ! endif
1868 !
1869 !! Define the lat-lon coordinate:
1870 ! allocate ( lon(im) )
1871 ! allocate ( lat(jm) )
1872 !
1873 ! do i=1,im
1874 ! lon(i) = (0.5 + real(i-1)) * 2.*pi/real(im)
1875 ! enddo
1876 !
1877 ! do j=1,jm
1878 ! lat(j) = -0.5*pi + real(j-1)*pi/real(jm-1) ! SP to NP
1879 ! enddo
1880 !
1881 ! allocate ( ak0(1:km+1) )
1882 ! allocate ( bk0(1:km+1) )
1883 ! allocate ( ps0(1:im,1:jm) )
1884 ! allocate ( gz0(1:im,1:jm) )
1885 ! allocate ( u0(1:im,1:jm,1:km) )
1886 ! allocate ( v0(1:im,1:jm,1:km) )
1887 ! allocate ( t0(1:im,1:jm,1:km) )
1888 ! allocate ( dp0(1:im,1:jm,1:km) )
1889 !
1890 ! call read_data (fname, 'ak', ak0)
1891 ! call read_data (fname, 'bk', bk0)
1892 ! call read_data (fname, 'Surface_geopotential', gz0)
1893 ! call read_data (fname, 'U', u0)
1894 ! call read_data (fname, 'V', v0)
1895 ! call read_data (fname, 'T', t0)
1896 ! call read_data (fname, 'DELP', dp0)
1897 !
1898 !! Share the load
1899 ! if(is_master()) call pmaxmin( 'ZS_data', gz0, im, jm, 1./grav)
1900 ! if(is_master()) call pmaxmin( 'U_data', u0, im*jm, km, 1._FVPRC)
1901 ! if(is_master()) call pmaxmin( 'V_data', v0, im*jm, km, 1._FVPRC)
1902 ! if(is_master()) call pmaxmin( 'T_data', t0, im*jm, km, 1._FVPRC)
1903 ! if(is_master()) call pmaxmin( 'DEL-P', dp0, im*jm, km, 0.01_FVPRC)
1904 !
1905 !
1906 ! else
1907 ! call mpp_error(FATAL,'==> Error in get_external_ic: Expected file '//trim(fname)//' does not exist')
1908 ! endif
1909 !
1910 !! Read in tracers: only AM2 "physics tracers" at this point
1911 ! fname = Atm(1)%res_latlon_tracers
1912 !
1913 ! if( file_exist(fname) ) then
1914 ! if(is_master()) write(*,*) 'Using lat-lon tracer restart:', fname
1915 !
1916 ! allocate ( q0(im,jm,km,Atm(1)%ncnst) )
1917 ! q0 = 0.
1918 !
1919 ! do tr_ind = 1, nq
1920 ! call get_tracer_names(MODEL_ATMOS, tr_ind, tracer_name)
1921 ! if (field_exist(fname,tracer_name)) then
1922 ! call read_data(fname, tracer_name, q0(1:im,1:jm,1:km,tr_ind))
1923 ! call mpp_error(NOTE,'==> Have read tracer '//trim(tracer_name)//' from '//trim(fname))
1924 ! cycle
1925 ! endif
1926 ! enddo
1927 ! else
1928 ! call mpp_error(FATAL,'==> Error in get_external_ic: Expected file '//trim(fname)//' does not exist')
1929 ! endif
1930 !
1931 !! D to A transform on lat-lon grid:
1932 ! allocate ( ua(im,jm,km) )
1933 ! allocate ( va(im,jm,km) )
1934 !
1935 ! call d2a3d(u0, v0, ua, va, im, jm, km, lon)
1936 !
1937 ! deallocate ( u0 )
1938 ! deallocate ( v0 )
1939 !
1940 ! if(is_master()) call pmaxmin( 'UA', ua, im*jm, km, 1._FVPRC)
1941 ! if(is_master()) call pmaxmin( 'VA', va, im*jm, km, 1._FVPRC)
1942 !
1943 ! do j=1,jm
1944 ! do i=1,im
1945 ! ps0(i,j) = ak0(1)
1946 ! enddo
1947 ! enddo
1948 !
1949 ! do k=1,km
1950 ! do j=1,jm
1951 ! do i=1,im
1952 ! ps0(i,j) = ps0(i,j) + dp0(i,j,k)
1953 ! enddo
1954 ! enddo
1955 ! enddo
1956 !
1957 ! if (is_master()) call pmaxmin( 'PS_data (mb)', ps0, im, jm, 0.01_FVPRC)
1958 !
1959 !! Horizontal interpolation to the cubed sphere grid center
1960 !! remap vertically with terrain adjustment
1961 !
1962 ! call remap_xyz( im, 1, jm, jm, km, npz, nq, Atm(1)%ncnst, lon, lat, ak0, bk0, &
1963 ! ps0, gz0, ua, va, t0, q0, Atm )
1964 !
1965 ! deallocate ( ak0 )
1966 ! deallocate ( bk0 )
1967 ! deallocate ( ps0 )
1968 ! deallocate ( gz0 )
1969 ! deallocate ( t0 )
1970 ! deallocate ( q0 )
1971 ! deallocate ( dp0 )
1972 ! deallocate ( ua )
1973 ! deallocate ( va )
1974 ! deallocate ( lat )
1975 ! deallocate ( lon )
1976 #endif
1977 
1978  end subroutine get_fv_ic
1979 
1980 
1981 #ifndef DYCORE_SOLO
1982  subroutine ncep2fms(im, jm, lon, lat, wk)
1983 
1984  integer, intent(in):: im, jm
1985  real(FVPRC), intent(in):: lon(im), lat(jm)
1986  real(kind=4), intent(in):: wk(im,jm)
1987 ! local:
1988  real(FVPRC) :: rdlon(im)
1989  real(FVPRC) :: rdlat(jm)
1990  real(FVPRC):: a1, b1
1991  real(FVPRC):: delx, dely
1992  real(FVPRC):: xc, yc ! "data" location
1993  real(FVPRC):: c1, c2, c3, c4
1994  integer i,j, i1, i2, jc, i0, j0, it, jt
1995 
1996  do i=1,im-1
1997  rdlon(i) = 1. / (lon(i+1) - lon(i))
1998  enddo
1999  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
2000 
2001  do j=1,jm-1
2002  rdlat(j) = 1. / (lat(j+1) - lat(j))
2003  enddo
2004 
2005 ! * Interpolate to "FMS" 1x1 SST data grid
2006 ! lon: 0.5, 1.5, ..., 359.5
2007 ! lat: -89.5, -88.5, ... , 88.5, 89.5
2008 
2009  delx = 360./real(i_sst)
2010  dely = 180./real(j_sst)
2011 
2012  jt = 1
2013  do 5000 j=1,j_sst
2014 
2015  yc = (-90. + dely * (0.5+real(j-1))) * deg2rad
2016  if ( yc<lat(1) ) then
2017  jc = 1
2018  b1 = 0.
2019  elseif ( yc>lat(jm) ) then
2020  jc = jm-1
2021  b1 = 1.
2022  else
2023  do j0=jt,jm-1
2024  if ( yc>=lat(j0) .and. yc<=lat(j0+1) ) then
2025  jc = j0
2026  jt = j0
2027  b1 = (yc-lat(jc)) * rdlat(jc)
2028  go to 222
2029  endif
2030  enddo
2031  endif
2032 222 continue
2033  it = 1
2034 
2035  do i=1,i_sst
2036  xc = delx * (0.5+real(i-1)) * deg2rad
2037  if ( xc>lon(im) ) then
2038  i1 = im; i2 = 1
2039  a1 = (xc-lon(im)) * rdlon(im)
2040  elseif ( xc<lon(1) ) then
2041  i1 = im; i2 = 1
2042  a1 = (xc+2.*pi-lon(im)) * rdlon(im)
2043  else
2044  do i0=it,im-1
2045  if ( xc>=lon(i0) .and. xc<=lon(i0+1) ) then
2046  i1 = i0; i2 = i0+1
2047  it = i0
2048  a1 = (xc-lon(i1)) * rdlon(i0)
2049  go to 111
2050  endif
2051  enddo
2052  endif
2053 111 continue
2054 
2055  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
2056  write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
2057  endif
2058 
2059  c1 = (1.-a1) * (1.-b1)
2060  c2 = a1 * (1.-b1)
2061  c3 = a1 * b1
2062  c4 = (1.-a1) * b1
2063 ! Interpolated surface pressure
2064  sst_ncep(i,j) = c1*wk(i1,jc ) + c2*wk(i2,jc ) + &
2065  c3*wk(i2,jc+1) + c4*wk(i1,jc+1)
2066  enddo !i-loop
2067 5000 continue ! j-loop
2068 
2069  end subroutine ncep2fms
2070 #endif
2071 
2072 
2073  subroutine remap_coef( im, jm, lon, lat, id1, id2, jdc, s2c, agrid, bd )
2074 
2075  type(fv_grid_bounds_type), intent(IN) :: bd
2076  integer, intent(in):: im, jm
2077  real(FVPRC), intent(in):: lon(im), lat(jm)
2078  real(FVPRC), intent(out):: s2c(bd%is:bd%ie,bd%js:bd%je,4)
2079  integer, intent(out), dimension(bd%is:bd%ie,bd%js:bd%je):: id1, id2, jdc
2080  real(FVPRC), intent(in):: agrid(bd%isd:bd%ied,bd%jsd:bd%jed,2)
2081 ! local:
2082  real(FVPRC) :: rdlon(im)
2083  real(FVPRC) :: rdlat(jm)
2084  real(FVPRC):: a1, b1
2085  integer i,j, i1, i2, jc, i0, j0
2086 
2087  integer :: is, ie, js, je
2088  integer :: isd, ied, jsd, jed
2089 
2090  is = bd%is
2091  ie = bd%ie
2092  js = bd%js
2093  je = bd%je
2094  isd = bd%isd
2095  ied = bd%ied
2096  jsd = bd%jsd
2097  jed = bd%jed
2098  do i=1,im-1
2099  rdlon(i) = 1. / (lon(i+1) - lon(i))
2100  enddo
2101  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
2102 
2103  do j=1,jm-1
2104  rdlat(j) = 1. / (lat(j+1) - lat(j))
2105  enddo
2106 
2107 ! * Interpolate to cubed sphere cell center
2108  do 5000 j=js,je
2109 
2110  do i=is,ie
2111 
2112  if ( agrid(i,j,1)>lon(im) ) then
2113  i1 = im; i2 = 1
2114  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
2115  elseif ( agrid(i,j,1)<lon(1) ) then
2116  i1 = im; i2 = 1
2117  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
2118  else
2119  do i0=1,im-1
2120  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
2121  i1 = i0; i2 = i0+1
2122  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
2123  go to 111
2124  endif
2125  enddo
2126  endif
2127 111 continue
2128  if ( agrid(i,j,2)<lat(1) ) then
2129  jc = 1
2130  b1 = 0.
2131  elseif ( agrid(i,j,2)>lat(jm) ) then
2132  jc = jm-1
2133  b1 = 1.
2134  else
2135  do j0=1,jm-1
2136  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
2137  jc = j0
2138  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
2139  go to 222
2140  endif
2141  enddo
2142  endif
2143 222 continue
2144 
2145  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
2146  write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
2147  endif
2148 
2149  s2c(i,j,1) = (1.-a1) * (1.-b1)
2150  s2c(i,j,2) = a1 * (1.-b1)
2151  s2c(i,j,3) = a1 * b1
2152  s2c(i,j,4) = (1.-a1) * b1
2153  id1(i,j) = i1
2154  id2(i,j) = i2
2155  jdc(i,j) = jc
2156  enddo !i-loop
2157 5000 continue ! j-loop
2158 
2159  end subroutine remap_coef
2160 
2161  subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, Atm)
2162  type(fv_atmos_type), intent(inout) :: Atm
2163  integer, intent(in):: im, jm, km, npz, nq, ncnst
2164  real(FVPRC), intent(in):: ak0(km+1), bk0(km+1)
2165  real(FVPRC), intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: psc, gzc
2166  real(FVPRC), intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: ta
2167  real(FVPRC), intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km,ncnst):: qa
2168 ! local:
2169  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,km):: tp
2170  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0, pn0
2171  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1
2172  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1, pn1
2173  real(FVPRC) pt0(km), gz(km+1), pk0(km+1)
2174  real(FVPRC) qp(Atm%bd%is:Atm%bd%ie,km,ncnst)
2175  real(FVPRC) pst, p1, p2, alpha, rdg
2176  integer i,j,k, iq
2177  integer sphum
2178  integer :: is, ie, js, je
2179  integer :: isd, ied, jsd, jed
2180 
2181  is = atm%bd%is
2182  ie = atm%bd%ie
2183  js = atm%bd%js
2184  je = atm%bd%je
2185  isd = atm%bd%isd
2186  ied = atm%bd%ied
2187  jsd = atm%bd%jsd
2188  jed = atm%bd%jed
2189 
2190 
2191 #ifdef MAPL_MODE
2192  sphum = 1
2193 #else
2194  sphum = get_tracer_index(model_atmos, 'sphum')
2195 #endif
2196  if ( sphum/=1 ) then
2197  call mpp_error(fatal,'SPHUM must be 1st tracer')
2198  endif
2199 
2200  do 5000 j=js,je
2201  do i=is,ie
2202 
2203  do iq=1,ncnst
2204  do k=1,km
2205  qp(i,k,iq) = qa(i,j,k,iq)
2206  enddo
2207  enddo
2208 
2209  if ( t_is_tv ) then
2210 ! The "T" field in NCEP analysis is actually virtual temperature (Larry H. post processing)
2211 ! BEFORE 20051201
2212  do k=1,km
2213  tp(i,k) = ta(i,j,k)
2214  enddo
2215  else
2216  do k=1,km
2217  tp(i,k) = ta(i,j,k)*(1.+zvir*qp(i,k,sphum))
2218  enddo
2219  endif
2220 ! Tracers:
2221 
2222  do k=1,km+1
2223  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2224  pn0(i,k) = log(pe0(i,k))
2225  pk0(k) = pe0(i,k)**kappa
2226  enddo
2227 
2228 #ifdef USE_DATA_ZS
2229  atm% ps(i,j) = psc(i,j)
2230  atm%phis(i,j) = gzc(i,j)
2231 #else
2232 
2233 ! * Adjust interpolated ps to model terrain
2234  gz(km+1) = gzc(i,j)
2235  do k=km,1,-1
2236  gz(k) = gz(k+1) + rdgas*tp(i,k)*(pn0(i,k+1)-pn0(i,k))
2237  enddo
2238 ! Only lowest layer potential temp is needed
2239  pt0(km) = tp(i,km)/(pk0(km+1)-pk0(km))*(kappa*(pn0(i,km+1)-pn0(i,km)))
2240  if( atm%phis(i,j)>gzc(i,j) ) then
2241  do k=km,1,-1
2242  if( atm%phis(i,j) < gz(k) .and. &
2243  atm%phis(i,j) >= gz(k+1) ) then
2244  pst = pk0(k) + (pk0(k+1)-pk0(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
2245  go to 123
2246  endif
2247  enddo
2248  else
2249 ! Extrapolation into the ground
2250  pst = pk0(km+1) + (gzc(i,j)-atm%phis(i,j))/(cp_air*pt0(km))
2251  endif
2252 
2253 123 atm%ps(i,j) = pst**(1./kappa)
2254 #endif
2255  enddo !i-loop
2256 
2257 
2258  do i=is,ie
2259  pe1(i,1) = atm%ak(1)
2260  pn1(i,1) = log(pe1(i,1))
2261  enddo
2262  do k=2,npz+1
2263  do i=is,ie
2264  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
2265  pn1(i,k) = log(pe1(i,k))
2266  enddo
2267  enddo
2268 
2269 ! * Compute delp
2270  do k=1,npz
2271  do i=is,ie
2272  atm%delp(i,j,k) = pe1(i,k+1) - pe1(i,k)
2273  enddo
2274  enddo
2275 
2276 !---------------
2277 ! map tracers
2278 !----------------
2279  do iq=1,ncnst
2280  call mappm(km, pe0, qp(is,1,iq), npz, pe1, qn1, is,ie, 0, 11, atm%ptop)
2281  if ( iq==sphum .and. atm%flagstruct%ncep_ic ) then
2282  p1 = 200.e2
2283  p2 = 75.e2
2284 ! Blend model sphum with NCEP data
2285  do k=1,npz
2286  do i=is,ie
2287  pst = 0.5*(pe1(i,k)+pe1(i,k+1))
2288  if ( pst > p1 ) then
2289  atm%q(i,j,k,iq) = qn1(i,k)
2290  elseif( pst > p2 ) then ! p2 < pst < p1
2291  alpha = (pst-p2)/(p1-p2)
2292  atm%q(i,j,k,1) = qn1(i,k)*alpha + atm%q(i,j,k,1)*(1.-alpha)
2293  endif
2294  enddo
2295  enddo
2296  else
2297  do k=1,npz
2298  do i=is,ie
2299  atm%q(i,j,k,iq) = qn1(i,k)
2300  enddo
2301  enddo
2302  endif
2303  enddo
2304 
2305 !-------------------------------------------------------------
2306 ! map virtual temperature using geopotential conserving scheme.
2307 !-------------------------------------------------------------
2308  call mappm(km, pn0, tp, npz, pn1, qn1, is,ie, 1, 9, atm%ptop)
2309  do k=1,npz
2310  do i=is,ie
2311  atm%pt(i,j,k) = qn1(i,k)/(1.+zvir*atm%q(i,j,k,sphum))
2312  enddo
2313  enddo
2314 
2315  if ( .not. atm%flagstruct%hydrostatic .and. atm%flagstruct%ncep_ic ) then
2316 ! Replace delz with NCEP hydrostatic state
2317  rdg = -rdgas / grav
2318  do k=1,npz
2319  do i=is,ie
2320  atm%delz(i,j,k) = rdg*qn1(i,k)*(pn1(i,k+1)-pn1(i,k))
2321  enddo
2322  enddo
2323  endif
2324 
2325 5000 continue
2326 
2327  call prt_maxmin('PS_model', atm%ps, is, ie, js, je, ng, 1, 0.01_fvprc)
2328 
2329  if (is_master()) write(*,*) 'done remap_scalar'
2330 
2331  end subroutine remap_scalar
2332 
2333  subroutine remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
2334  type(fv_atmos_type), intent(inout) :: Atm(:)
2335  integer, intent(in):: im, jm, km, npz
2336  real(FVPRC), intent(in):: ak0(km+1), bk0(km+1)
2337  real(FVPRC), intent(in):: psc(is:ie,js:je)
2338  real(FVPRC), intent(in), dimension(is:ie,js:je,km):: ua, va
2339 ! local:
2340  real(FVPRC), dimension(isd:ied,jsd:jed,npz):: ut, vt ! winds
2341  real(FVPRC), dimension(is:ie, km+1):: pe0
2342  real(FVPRC), dimension(is:ie,npz+1):: pe1
2343  real(FVPRC), dimension(is:ie,npz):: qn1
2344  integer i,j,k
2345 
2346  ut = 0.0
2347  vt = 0.0
2348  do 5000 j=js,je
2349 
2350  do k=1,km+1
2351  do i=is,ie
2352  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2353  enddo
2354  enddo
2355 
2356  do k=1,npz+1
2357  do i=is,ie
2358  pe1(i,k) = atm(1)%ak(k) + atm(1)%bk(k)*atm(1)%ps(i,j)
2359  enddo
2360  enddo
2361 
2362 !------
2363 ! map u
2364 !------
2365  call mappm(km, pe0, ua(is:ie,j,1:km), npz, pe1, qn1, is,ie, -1, 4, atm(1)%ptop)
2366  do k=1,npz
2367  do i=is,ie
2368  ut(i,j,k) = qn1(i,k)
2369  enddo
2370  enddo
2371 !------
2372 ! map v
2373 !------
2374  call mappm(km, pe0, va(is:ie,j,1:km), npz, pe1, qn1, is,ie, -1, 4, atm(1)%ptop)
2375  do k=1,npz
2376  do i=is,ie
2377  vt(i,j,k) = qn1(i,k)
2378  enddo
2379  enddo
2380 
2381 5000 continue
2382 
2383  call prt_maxmin('UT', ut, is, ie, js, je, ng, npz, 1.0_fvprc)
2384  call prt_maxmin('VT', vt, is, ie, js, je, ng, npz, 1.0_fvprc)
2385 
2386 !----------------------------------------------
2387 ! winds: lat-lon ON A to Cubed-D transformation:
2388 !----------------------------------------------
2389  call cubed_a2d(atm(1)%npx, atm(1)%npy, npz, ut, vt, atm(1)%u, atm(1)%v, atm(1)%gridstruct, &
2390  atm(1)%domain, atm(1)%bd )
2391 
2392  if (is_master()) write(*,*) 'done remap_winds'
2393 
2394  end subroutine remap_winds
2395 
2396 
2397  subroutine remap_wz(im, jm, km, npz, mg, ak0, bk0, psc, wa, wz, Atm)
2398  type(fv_atmos_type), intent(inout) :: Atm(:)
2399  integer, intent(in):: im, jm, km, npz
2400  integer, intent(in):: mg ! mg = 0 for delz; mg=3 for w
2401  real(FVPRC), intent(in):: ak0(km+1), bk0(km+1)
2402  real(FVPRC), intent(in):: psc(is:ie,js:je)
2403  real(FVPRC), intent(in), dimension(is:ie,js:je,km):: wa
2404  real(FVPRC), intent(out):: wz(is-mg:ie+mg,js-mg:je+mg,npz)
2405 ! local:
2406  real(FVPRC), dimension(is:ie, km+1):: pe0
2407  real(FVPRC), dimension(is:ie,npz+1):: pe1
2408  real(FVPRC), dimension(is:ie,npz):: qn1
2409  integer i,j,k
2410 
2411  do 5000 j=js,je
2412 
2413  do k=1,km+1
2414  do i=is,ie
2415  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2416  enddo
2417  enddo
2418 
2419  do k=1,npz+1
2420  do i=is,ie
2421  pe1(i,k) = atm(1)%ak(k) + atm(1)%bk(k)*atm(1)%ps(i,j)
2422  enddo
2423  enddo
2424 
2425 !------
2426 ! map w
2427 !------
2428  call mappm(km, pe0, wa(is:ie,j,1:km), npz, pe1, qn1, is,ie, -1, 4, atm(1)%ptop)
2429  do k=1,npz
2430  do i=is,ie
2431  wz(i,j,k) = qn1(i,k)
2432  enddo
2433  enddo
2434 
2435 5000 continue
2436 
2437 ! call prt_maxmin('WZ', wz, is, ie, js, je, mg, npz, 1._FVPRC, is_master())
2438 ! if (is_master()) write(*,*) 'done remap_wz'
2439 
2440  end subroutine remap_wz
2441 
2442  subroutine remap_xyz( im, jbeg, jend, jm, km, npz, nq, ncnst, lon, lat, ak0, bk0, ps0, gz0, &
2443  ua, va, ta, qa, Atm )
2444 
2445  type(fv_atmos_type), intent(inout), target :: Atm
2446  integer, intent(in):: im, jm, km, npz, nq, ncnst
2447  integer, intent(in):: jbeg, jend
2448  real(FVPRC), intent(in):: lon(im), lat(jm), ak0(km+1), bk0(km+1)
2449  real(FVPRC), intent(in):: gz0(im,jbeg:jend), ps0(im,jbeg:jend)
2450  real(FVPRC), intent(in), dimension(im,jbeg:jend,km):: ua, va, ta
2451  real(FVPRC), intent(in), dimension(im,jbeg:jend,km,ncnst):: qa
2452 
2453  real(FVPRC), pointer, dimension(:,:,:) :: agrid
2454 
2455 ! local:
2456  real(FVPRC), dimension(Atm%bd%isd:Atm%bd%ied,Atm%bd%jsd:Atm%bd%jed,npz):: ut, vt ! winds
2457  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,km):: up, vp, tp
2458  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0, pn0
2459  real(FVPRC) pt0(km), gz(km+1), pk0(km+1)
2460  real(FVPRC) qp(Atm%bd%is:Atm%bd%ie,km,ncnst)
2461  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1
2462  real(FVPRC), dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1, pn1
2463  real(FVPRC) :: rdlon(im)
2464  real(FVPRC) :: rdlat(jm)
2465  real(FVPRC):: a1, b1, c1, c2, c3, c4
2466  real(FVPRC):: gzc, psc, pst
2467  integer i,j,k, i1, i2, jc, i0, j0, iq
2468 ! integer sphum, liq_wat, ice_wat, cld_amt
2469  integer sphum
2470  integer :: is, ie, js, je
2471  integer :: isd, ied, jsd, jed
2472 
2473  is = atm%bd%is
2474  ie = atm%bd%ie
2475  js = atm%bd%js
2476  je = atm%bd%je
2477  isd = atm%bd%isd
2478  ied = atm%bd%ied
2479  jsd = atm%bd%jsd
2480  jed = atm%bd%jed
2481 
2482  !!NOTE: Only Atm is used in this routine.
2483  agrid => atm%gridstruct%agrid
2484 
2485  sphum = get_tracer_index(model_atmos, 'sphum')
2486 ! liq_wat = get_tracer_index(MODEL_ATMOS, 'liq_wat')
2487 ! ice_wat = get_tracer_index(MODEL_ATMOS, 'ice_wat')
2488 ! cld_amt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
2489 
2490  if ( sphum/=1 ) then
2491  call mpp_error(fatal,'SPHUM must be 1st tracer')
2492  endif
2493  pk0(1) = ak0(1)**kappa
2494 
2495  do i=1,im-1
2496  rdlon(i) = 1. / (lon(i+1) - lon(i))
2497  enddo
2498  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
2499 
2500  do j=1,jm-1
2501  rdlat(j) = 1. / (lat(j+1) - lat(j))
2502  enddo
2503 
2504 ! * Interpolate to cubed sphere cell center
2505  do 5000 j=js,je
2506 
2507  do i=is,ie
2508  pe0(i,1) = ak0(1)
2509  pn0(i,1) = log(ak0(1))
2510  enddo
2511 
2512  do i=is,ie
2513 
2514  if ( agrid(i,j,1)>lon(im) ) then
2515  i1 = im; i2 = 1
2516  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
2517  elseif ( agrid(i,j,1)<lon(1) ) then
2518  i1 = im; i2 = 1
2519  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
2520  else
2521  do i0=1,im-1
2522  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
2523  i1 = i0; i2 = i0+1
2524  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
2525  go to 111
2526  endif
2527  enddo
2528  endif
2529 
2530 111 continue
2531  if ( agrid(i,j,2)<lat(1) ) then
2532  jc = 1
2533  b1 = 0.
2534  elseif ( agrid(i,j,2)>lat(jm) ) then
2535  jc = jm-1
2536  b1 = 1.
2537  else
2538  do j0=1,jm-1
2539  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
2540  jc = j0
2541  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
2542  go to 222
2543  endif
2544  enddo
2545  endif
2546 222 continue
2547 
2548 #ifndef DEBUG_REMAP
2549  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
2550  write(*,*) i,j,a1, b1
2551  endif
2552 #endif
2553  c1 = (1.-a1) * (1.-b1)
2554  c2 = a1 * (1.-b1)
2555  c3 = a1 * b1
2556  c4 = (1.-a1) * b1
2557 
2558 ! Interpolated surface pressure
2559  psc = c1*ps0(i1,jc ) + c2*ps0(i2,jc ) + &
2560  c3*ps0(i2,jc+1) + c4*ps0(i1,jc+1)
2561 
2562 ! Interpolated surface geopotential
2563  gzc = c1*gz0(i1,jc ) + c2*gz0(i2,jc ) + &
2564  c3*gz0(i2,jc+1) + c4*gz0(i1,jc+1)
2565 
2566 ! 3D fields:
2567  do iq=1,ncnst
2568 ! if ( iq==sphum .or. iq==liq_wat .or. iq==ice_wat .or. iq==cld_amt ) then
2569  do k=1,km
2570  qp(i,k,iq) = c1*qa(i1,jc, k,iq) + c2*qa(i2,jc, k,iq) + &
2571  c3*qa(i2,jc+1,k,iq) + c4*qa(i1,jc+1,k,iq)
2572  enddo
2573 ! endif
2574  enddo
2575  do k=1,km
2576  up(i,k) = c1*ua(i1,jc, k) + c2*ua(i2,jc, k) + &
2577  c3*ua(i2,jc+1,k) + c4*ua(i1,jc+1,k)
2578  vp(i,k) = c1*va(i1,jc, k) + c2*va(i2,jc, k) + &
2579  c3*va(i2,jc+1,k) + c4*va(i1,jc+1,k)
2580  tp(i,k) = c1*ta(i1,jc, k) + c2*ta(i2,jc, k) + &
2581  c3*ta(i2,jc+1,k) + c4*ta(i1,jc+1,k)
2582 ! Virtual effect:
2583  tp(i,k) = tp(i,k)*(1.+zvir*qp(i,k,sphum))
2584  enddo
2585 ! Tracers:
2586 
2587  do k=2,km+1
2588  pe0(i,k) = ak0(k) + bk0(k)*psc
2589  pn0(i,k) = log(pe0(i,k))
2590  pk0(k) = pe0(i,k)**kappa
2591  enddo
2592 
2593 #ifdef USE_DATA_ZS
2594  atm% ps(i,j) = psc
2595  atm%phis(i,j) = gzc
2596 #else
2597 
2598 ! * Adjust interpolated ps to model terrain
2599  gz(km+1) = gzc
2600  do k=km,1,-1
2601  gz(k) = gz(k+1) + rdgas*tp(i,k)*(pn0(i,k+1)-pn0(i,k))
2602  enddo
2603 ! Only lowest layer potential temp is needed
2604  pt0(km) = tp(i,km)/(pk0(km+1)-pk0(km))*(kappa*(pn0(i,km+1)-pn0(i,km)))
2605  if( atm%phis(i,j)>gzc ) then
2606  do k=km,1,-1
2607  if( atm%phis(i,j) < gz(k) .and. &
2608  atm%phis(i,j) >= gz(k+1) ) then
2609  pst = pk0(k) + (pk0(k+1)-pk0(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
2610  go to 123
2611  endif
2612  enddo
2613  else
2614 ! Extrapolation into the ground
2615  pst = pk0(km+1) + (gzc-atm%phis(i,j))/(cp_air*pt0(km))
2616  endif
2617 
2618 123 atm%ps(i,j) = pst**(1./kappa)
2619 #endif
2620  enddo !i-loop
2621 
2622 
2623 ! * Compute delp from ps
2624  do i=is,ie
2625  pe1(i,1) = atm%ak(1)
2626  pn1(i,1) = log(pe1(i,1))
2627  enddo
2628  do k=2,npz+1
2629  do i=is,ie
2630  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
2631  pn1(i,k) = log(pe1(i,k))
2632  enddo
2633  enddo
2634 
2635  do k=1,npz
2636  do i=is,ie
2637  atm%delp(i,j,k) = pe1(i,k+1) - pe1(i,k)
2638  enddo
2639  enddo
2640 
2641 ! Use kord=9 for winds; kord=11 for tracers
2642 !------
2643 ! map u
2644 !------
2645  call mappm(km, pe0, up, npz, pe1, qn1, is,ie, -1, 9, atm%ptop)
2646  do k=1,npz
2647  do i=is,ie
2648  ut(i,j,k) = qn1(i,k)
2649  enddo
2650  enddo
2651 !------
2652 ! map v
2653 !------
2654  call mappm(km, pe0, vp, npz, pe1, qn1, is,ie, -1, 9, atm%ptop)
2655  do k=1,npz
2656  do i=is,ie
2657  vt(i,j,k) = qn1(i,k)
2658  enddo
2659  enddo
2660 
2661 !---------------
2662 ! map tracers
2663 !----------------
2664  do iq=1,ncnst
2665 ! Note: AM2 physics tracers only
2666 ! if ( iq==sphum .or. iq==liq_wat .or. iq==ice_wat .or. iq==cld_amt ) then
2667  call mappm(km, pe0, qp(is,1,iq), npz, pe1, qn1, is,ie, 0, 11, atm%ptop)
2668  do k=1,npz
2669  do i=is,ie
2670  atm%q(i,j,k,iq) = qn1(i,k)
2671  enddo
2672  enddo
2673 ! endif
2674  enddo
2675 
2676 !-------------------------------------------------------------
2677 ! map virtual temperature using geopotential conserving scheme.
2678 !-------------------------------------------------------------
2679  call mappm(km, pn0, tp, npz, pn1, qn1, is,ie, 1, 9, atm%ptop)
2680  do k=1,npz
2681  do i=is,ie
2682  atm%pt(i,j,k) = qn1(i,k)/(1.+zvir*atm%q(i,j,k,sphum))
2683  enddo
2684  enddo
2685 
2686 5000 continue
2687 
2688  call prt_maxmin('PS_model', atm%ps, is, ie, js, je, ng, 1, 0.01_fvprc)
2689  call prt_maxmin('UT', ut, is, ie, js, je, ng, npz, 1._fvprc)
2690  call prt_maxmin('VT', vt, is, ie, js, je, ng, npz, 1._fvprc)
2691 
2692 !----------------------------------------------
2693 ! winds: lat-lon ON A to Cubed-D transformation:
2694 !----------------------------------------------
2695  call cubed_a2d(atm%npx, atm%npy, npz, ut, vt, atm%u, atm%v, atm%gridstruct, atm%domain, atm%bd )
2696 
2697  if (is_master()) write(*,*) 'done remap_xyz'
2698 
2699  end subroutine remap_xyz
2700 
2701 
2702 
2703 ! subroutine init_cubsph_grid(npts, is,ie, js,je, ntiles, sph_corner) !------------------------------------------------------------------! ! read/generate cubed sphere grid ! ! calculate cell center from cell corners !
2704 !! !
2705 !! input: !
2706 !! npts, is,ie, js,je, ntiles number of grid points and tiles !
2707 !! !
2708 !! output: !
2709 !! sph_corner cell corners in spherical coor !
2710 !!------------------------------------------------------------------!
2711 ! use GHOST_CUBSPH_mod, only: B_grid, ghost_cubsph_update
2712 ! use fv_grid_utils_nlm_mod, only : gnomonic_grids
2713 ! use fv_grid_tools_nlm_mod, only : mirror_grid
2714 !
2715 ! integer, intent(in) :: npts, is,ie, js,je, ntiles
2716 ! real*8, dimension(2,is:ie+1,js:je+1), intent(out) :: sph_corner
2717 !!------------------------------------------------------------------!
2718 !! local variables !
2719 !!------------------------------------------------------------------!
2720 ! integer :: i, j, l, n
2721 ! real*8, pointer :: xs(:,:), ys(:,:)
2722 ! real*8, pointer :: grid_in(:,:,:,:)
2723 ! integer :: grid_type = 0
2724 !!------------------------------------------------------------------!
2725 !! create sph_corner !
2726 !!------------------------------------------------------------------!
2727 !#ifdef SMEM_MAPL_MODE
2728 !! allocate global arrays (preferable in shared memory)
2729 ! if(MAPL_ShmInitialized) then
2730 ! if (is_master()) write(*,*) 'Using MAPL_Shmem in external_ic: init_cubsph_grid'
2731 ! call MAPL_AllocNodeArray(grid_in,Shp=(/npts,npts,2,ntiles/),rc=STATUS)
2732 ! else
2733 ! if (is_master()) write(*,*) 'WARNING... in external_ic: Global grid allocate'
2734 ! allocate( grid_in(npts,npts,2,ntiles) )
2735 ! endif
2736 ! if (is_master()) then
2737 ! allocate( xs(npts,npts) )
2738 ! allocate( ys(npts,npts) )
2739 ! call gnomonic_grids(grid_type, npts-1, xs, ys)
2740 ! do j=1,npts
2741 ! do i=1,npts
2742 ! grid_in(i,j,1,1) = xs(i,j)
2743 ! grid_in(i,j,2,1) = ys(i,j)
2744 ! enddo
2745 ! enddo
2746 ! deallocate ( xs )
2747 ! deallocate ( ys )
2748 ! endif
2749 ! if(MAPL_ShmInitialized) then
2750 ! call MAPL_SyncSharedMemory(rc=STATUS)
2751 ! call MAPL_BroadcastToNodes( grid_in, N=size(grid_in), ROOT=masterproc, RC=status)
2752 ! call MAPL_SyncSharedMemory(rc=STATUS)
2753 ! else
2754 ! call mpp_broadcast(grid_in, size(grid_in), masterproc)
2755 ! endif
2756 !#else
2757 ! allocate( xs(npts,npts) )
2758 ! allocate( ys(npts,npts) )
2759 ! allocate( grid_in(npts,npts,2,ntiles) )
2760 ! call gnomonic_grids(grid_type, npts-1, xs, ys)
2761 ! do j=1,npts
2762 ! do i=1,npts
2763 ! grid_in(i,j,1,1) = xs(i,j)
2764 ! grid_in(i,j,2,1) = ys(i,j)
2765 ! enddo
2766 ! enddo
2767 ! deallocate ( xs )
2768 ! deallocate ( ys )
2769 !#endif
2770 !
2771 !! mirror_grid assumes that the tile=1 is centered on equator and greenwich meridian Lon[-pi,pi]
2772 ! call mirror_grid(grid_in, 0, npts, npts, 2, ntiles)
2773 ! sph_corner(1,is:ie+1,js:je+1) = grid_in(is:ie+1,js:je+1,1,tile)
2774 ! sph_corner(2,is:ie+1,js:je+1) = grid_in(is:ie+1,js:je+1,2,tile)
2775 ! do j=js,je+1
2776 ! do i=is,ie+1
2777 !!---------------------------------
2778 !! Shift the corner away from Japan
2779 !!---------------------------------
2780 !! This will result in the corner close to east coast of China
2781 ! sph_corner(1,i,j) = sph_corner(1,i,j) - pi/18.
2782 ! if ( sph_corner(1,i,j) < 0. ) &
2783 ! sph_corner(1,i,j) = sph_corner(1,i,j) + 2.*pi
2784 ! if (ABS(sph_corner(1,i,j)) < 1.e-10) sph_corner(1,i,j) = 0.0
2785 ! if (ABS(sph_corner(2,i,j)) < 1.e-10) sph_corner(2,i,j) = 0.0
2786 ! enddo
2787 ! enddo
2788 !#ifdef SMEM_MAPL_MODE
2789 ! call MAPL_SyncSharedMemory(rc=STATUS)
2790 ! DEALLOCGLOB_(grid_in)
2791 ! call MAPL_SyncSharedMemory(rc=STATUS)
2792 !#else
2793 ! deallocate ( grid_in )
2794 !#endif
2795 !!------------------------------------------------------------------!
2796 !! do halo update !
2797 !!------------------------------------------------------------------!
2798 !
2799 ! end subroutine init_cubsph_grid
2800 
2801 ! subroutine interp_c2c_vect(npx_in, npy_in, npx_out, npy_out, npz, ntiles, domain_i, &
2802 ! is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i, &
2803 ! u_in, v_in, u_out, v_out, index_c2c, weight_c2c, corner_in, corner_out, gridstruct)
2804 ! use GRID_UTILS_mod, only: latlon2xyz
2805 ! use GRID_UTILS_mod, only: get_dx, get_dxa, get_dy, get_dya, &
2806 ! get_center_vect, get_west_vect, &
2807 ! get_south_vect, get_cosa_center
2808 ! use FLOW_PROJ_mod, only: d2a, d2a_vect, a2d_vect
2809 ! use GHOST_CUBSPH_mod, only : A_grid, ghost_cubsph_update
2810 ! use CUB2CUB_mod, only: do_c2c_interpolation
2811 ! integer, intent(IN) :: npx_in, npy_in, npx_out, npy_out, npz, ntiles
2812 ! type(domain2d), intent(INOUT) :: domain_i
2813 ! integer, intent(IN) :: is,ie, js,je, isd_i,ied_i, jsd_i,jed_i, is_i,ie_i, js_i,je_i
2814 ! integer, intent(IN) :: index_c2c(3, is:ie, js:je )
2815 ! real(REAL64), intent(IN) :: weight_c2c(4, is:ie, js:je )
2816 ! real(REAL64), intent(IN) :: corner_in(2,is_i:ie_i+1,js_i:je_i+1)
2817 ! real(REAL64), intent(IN) :: corner_out(2,is:ie+1,js:je+1)
2818 ! real(FVPRC), intent(IN) :: u_in(isd_i:ied_i,jsd_i:jed_i+1,npz)
2819 ! real(FVPRC), intent(IN) :: v_in(isd_i:ied_i+1,jsd_i:jed_i,npz)
2820 ! real(FVPRC), intent(OUT):: u_out(is:ie,js:je,npz)
2821 ! real(FVPRC), intent(OUT):: v_out(is:ie,js:je,npz)
2822 ! type(fv_grid_type), intent(IN), target :: gridstruct
2823 !
2824 ! real(FVPRC) :: tmp(isd_i:ied_i,jsd_i:jed_i)
2825 !
2826 ! integer :: i,j,l,k,n
2827 ! real(REAL64), dimension(:,:,:), allocatable :: vxyz_in, vxyz_out
2828 ! real(REAL64), dimension(:,:,:), allocatable :: ec1, ec2, ew1, ew2, es1, es2
2829 ! real(REAL64), dimension(:,:) , allocatable :: dx, dy, dxa, dya, rdxa, rdya, cosa_s, sina_s
2830 ! real(FVPRC) :: u1, v1, vx, vy, vz
2831 ! real(REAL64) :: pc1(3), pc2(3)
2832 ! integer :: ic, jc, lc
2833 !
2834 ! real(REAL64), dimension(:,:,:), allocatable :: xyz_corner_in, xyz_corner_out
2835 !
2836 ! character(len=64) :: strTxt
2837 !
2838 !!------------------------------------------------------------------!
2839 !! calculate xyz cell corners and cell centers !
2840 !!------------------------------------------------------------------!
2841 ! allocate(xyz_corner_in (3, isd_i:ied_i+1, jsd_i:jed_i+1), &
2842 ! xyz_corner_out(3, is :ie +1, js :je +1))
2843 ! do j=js_i,je_i+1
2844 ! do i=is_i,ie_i+1
2845 ! call latlon2xyz(corner_in(:,i,j), xyz_corner_in(:,i,j))
2846 ! enddo
2847 ! enddo
2848 ! do j=js,je+1
2849 ! do i=is,ie+1
2850 ! call latlon2xyz(corner_out(:,i,j), xyz_corner_out(:,i,j))
2851 ! enddo
2852 ! enddo
2853 ! call print_memuse_stats('interp_c2c_vect: GRIDS')
2854 !
2855 !!----------------------------------------------------------!
2856 !! allocate horizontal flow variables !
2857 !!----------------------------------------------------------!
2858 ! allocate(vxyz_in(3,isd_i:ied_i,jsd_i:jed_i))
2859 ! allocate(dx(isd_i:ied_i,jsd_i:jed_i+1), dxa(isd_i:ied_i,jsd_i:jed_i), rdxa(isd_i:ied_i,jsd_i:jed_i))
2860 ! allocate(dy(isd_i:ied_i+1,jsd_i:jed_i), dya(isd_i:ied_i,jsd_i:jed_i), rdya(isd_i:ied_i,jsd_i:jed_i))
2861 ! allocate(ec1(3,isd_i:ied_i,jsd_i:jed_i), ec2(3,isd_i:ied_i,jsd_i:jed_i))
2862 ! allocate(cosa_s(isd_i:ied_i,jsd_i:jed_i), sina_s(isd_i:ied_i,jsd_i:jed_i))
2863 ! allocate(vxyz_out(3,is:ie,js:je))
2864 !
2865 !!-------------------------------------------------------!
2866 !! geometrical properties of input grid !
2867 !!-------------------------------------------------------!
2868 ! call get_dx (xyz_corner_in(:,isd_i:ied_i+1,jsd_i:jed_i+1), &
2869 ! isd_i,ied_i ,jsd_i,jed_i, &
2870 ! is_i ,ie_i ,js_i ,je_i , dx)
2871 ! call get_dxa(xyz_corner_in(:,isd_i:ied_i+1,jsd_i:jed_i+1), &
2872 ! isd_i,ied_i ,jsd_i,jed_i, &
2873 ! is_i ,ie_i ,js_i ,je_i , dxa, rdxa=rdxa)
2874 ! call get_dy (xyz_corner_in(:,isd_i:ied_i+1,jsd_i:jed_i+1), &
2875 ! isd_i,ied_i ,jsd_i,jed_i, &
2876 ! is_i ,ie_i ,js_i ,je_i , dy)
2877 ! call get_dya(xyz_corner_in(:,isd_i:ied_i+1,jsd_i:jed_i+1), &
2878 ! isd_i,ied_i ,jsd_i,jed_i, &
2879 ! is_i ,ie_i ,js_i ,je_i , dya, rdya=rdya)
2880 ! call get_center_vect(xyz_corner_in(:,isd_i:ied_i+1,jsd_i:jed_i+1), &
2881 ! isd_i,ied_i ,jsd_i,jed_i, &
2882 ! is_i ,ie_i ,js_i ,je_i , ec1, ec2)
2883 ! call get_cosa_center(ec1, ec2, isd_i,ied_i ,jsd_i,jed_i, &
2884 ! is_i ,ie_i ,js_i ,je_i , cosa_s, sina_s)
2885 !
2886 !! Flow interpolation for U and V components
2887 ! do k=1,npz
2888 !!-------------------------------------------------------!
2889 !! calculate flow vector for a-grid !
2890 !!-------------------------------------------------------!
2891 ! call d2a_vect(DBLE(u_in(:,:,k)), DBLE(v_in(:,:,k)), DBLE(dx), DBLE(dy), DBLE(rdxa), DBLE(rdya), DBLE(cosa_s), DBLE(ec1), DBLE(ec2), &
2892 ! isd_i, ied_i, jsd_i, jed_i, 1, 1, &
2893 ! is_i , ie_i , js_i , je_i , 1, 1, &
2894 ! vxyz_in(:,isd_i:ied_i,jsd_i:jed_i))
2895 ! write(strTxt,'(A,i3.3)') 'interp_c2c_vect: INPUT D2A level:', k
2896 ! if (k==npz) call print_memuse_stats(strTxt)
2897 !!----------------------------------------------------------!
2898 !! ghost cell update of vxyz_in !
2899 !!----------------------------------------------------------!
2900 ! do n=1,3
2901 ! tmp(is_i:ie_i,js_i:je_i) = vxyz_in(n,is_i:ie_i,js_i:je_i)
2902 ! call mpp_update_domains(tmp, domain_i)
2903 ! vxyz_in(n,:,:) = tmp
2904 ! enddo
2905 ! do n=1,3
2906 ! do j=js,je
2907 ! do i=is,ie
2908 !!----------------------------------------------------------!
2909 !! do interpolation of flow vector on A-Grid !
2910 !!----------------------------------------------------------!
2911 ! ic=index_c2c(1,i,j)
2912 ! jc=index_c2c(2,i,j)
2913 ! vxyz_out(n,i,j)=weight_c2c(1,i,j)*vxyz_in(n,ic ,jc ) &
2914 ! +weight_c2c(2,i,j)*vxyz_in(n,ic ,jc+1) &
2915 ! +weight_c2c(3,i,j)*vxyz_in(n,ic+1,jc+1) &
2916 ! +weight_c2c(4,i,j)*vxyz_in(n,ic+1,jc )
2917 ! enddo
2918 ! enddo
2919 ! enddo
2920 ! do j=js,je
2921 ! do i=is,ie
2922 ! vx = vxyz_out(1,i,j)
2923 ! vy = vxyz_out(2,i,j)
2924 ! vz = vxyz_out(3,i,j)
2925 !!----------------------------------------------------------!
2926 !! convert flow vector to wind vectors on A-Grid !
2927 !!----------------------------------------------------------!
2928 ! pc1(:)=xyz_corner_out(:,i+1,j)+xyz_corner_out(:,i+1,j+1) &
2929 ! -xyz_corner_out(:,i ,j)-xyz_corner_out(:,i ,j+1)
2930 ! pc2(:)=xyz_corner_out(:,i,j+1)+xyz_corner_out(:,i+1,j+1) &
2931 ! -xyz_corner_out(:,i,j )-xyz_corner_out(:,i+1,j )
2932 ! call normalize_vect(pc1(:))
2933 ! call normalize_vect(pc2(:))
2934 ! u1 = vx*pc1(1) + vy*pc1(2) + vz*pc1(3)
2935 ! v1 = vx*pc2(1) + vy*pc2(2) + vz*pc2(3)
2936 !!----------------------------------------------------------!
2937 !! rotate wind vectors from cubed to latlon orientation !
2938 !!----------------------------------------------------------!
2939 ! u_out(i,j,k) = 2.0*(gridstruct%a11(i,j)*u1 + gridstruct%a12(i,j)*v1)
2940 ! v_out(i,j,k) = 2.0*(gridstruct%a21(i,j)*u1 + gridstruct%a22(i,j)*v1)
2941 ! enddo
2942 ! enddo
2943 ! write(strTxt,'(A,i3.3)') 'interp_c2c_vect: OUTPUT A-grid C2L level:', k
2944 ! if (k==npz) call print_memuse_stats(strTxt)
2945 ! enddo ! npz
2946 !
2947 ! deallocate(dx, dy, dxa, dya, rdxa, rdya, ec1, ec2, cosa_s, sina_s)
2948 ! deallocate(vxyz_in)
2949 ! deallocate(vxyz_out)
2950 ! deallocate ( xyz_corner_in, xyz_corner_out )
2951 !
2952 ! end subroutine interp_c2c_vect
2953 
2954  Function reverse(A) Result(B)
2955  real(FVPRC), Intent(In) :: a(:,:)
2956  real(FVPRC) :: b(size(a,1),size(a,2))
2957 
2958  Integer :: i, n
2959 
2960  n = Size(a, 1)
2961 
2962  Do i = 1, n
2963  b(i,:) = a(1+n-i,:)
2964  End Do
2965 
2966  End Function reverse
2967 
2968  subroutine cubed_a2d( npx, npy, npz, ua, va, u, v, gridstruct, fv_domain, bd )
2969 
2970 ! Purpose; Transform wind on A grid to D grid
2971 
2972  use mpp_domains_mod, only: mpp_update_domains
2973 
2974  type(fv_grid_bounds_type), intent(IN) :: bd
2975  integer, intent(in):: npx, npy, npz
2976  real(FVPRC), intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va
2977  real(FVPRC), intent(out):: u(bd%isd:bd%ied, bd%jsd:bd%jed+1,npz)
2978  real(FVPRC), intent(out):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
2979  type(fv_grid_type), intent(IN), target :: gridstruct
2980  type(domain2d), intent(INOUT) :: fv_domain
2981 ! local:
2982  real(FVPRC) v3(3,bd%is-1:bd%ie+1,bd%js-1:bd%je+1)
2983  real(FVPRC) ue(3,bd%is-1:bd%ie+1,bd%js:bd%je+1) ! 3D winds at edges
2984  real(FVPRC) ve(3,bd%is:bd%ie+1,bd%js-1:bd%je+1) ! 3D winds at edges
2985  real(FVPRC), dimension(bd%is:bd%ie):: ut1, ut2, ut3
2986  real(FVPRC), dimension(bd%js:bd%je):: vt1, vt2, vt3
2987  integer i, j, k, im2, jm2
2988 
2989  real(REAL64), pointer, dimension(:,:,:) :: vlon, vlat
2990  real(REAL64), pointer, dimension(:) :: edge_vect_w, edge_vect_e, edge_vect_s, edge_vect_n
2991  real(REAL64), pointer, dimension(:,:,:,:) :: ew, es
2992 
2993  integer :: is, ie, js, je
2994  integer :: isd, ied, jsd, jed
2995 
2996  is = bd%is
2997  ie = bd%ie
2998  js = bd%js
2999  je = bd%je
3000  isd = bd%isd
3001  ied = bd%ied
3002  jsd = bd%jsd
3003  jed = bd%jed
3004  vlon => gridstruct%vlon
3005  vlat => gridstruct%vlat
3006 
3007  edge_vect_w => gridstruct%edge_vect_w
3008  edge_vect_e => gridstruct%edge_vect_e
3009  edge_vect_s => gridstruct%edge_vect_s
3010  edge_vect_n => gridstruct%edge_vect_n
3011 
3012  ew => gridstruct%ew
3013  es => gridstruct%es
3014 
3015  call mpp_update_domains(ua, fv_domain, complete=.false.)
3016  call mpp_update_domains(va, fv_domain, complete=.true.)
3017 
3018  im2 = (npx-1)/2
3019  jm2 = (npy-1)/2
3020  do k=1, npz
3021 ! Compute 3D wind on A grid
3022  do j=js-1,je+1
3023  do i=is-1,ie+1
3024  !v3(1,i,j) = ua(i,j,k)*vlon(1,i,j) + va(i,j,k)*vlat(1,i,j)
3025  !v3(2,i,j) = ua(i,j,k)*vlon(2,i,j) + va(i,j,k)*vlat(2,i,j)
3026  !v3(3,i,j) = ua(i,j,k)*vlon(3,i,j) + va(i,j,k)*vlat(3,i,j)
3027  v3(1,i,j) = ua(i,j,k)*vlon(i,j,1) + va(i,j,k)*vlat(i,j,1)
3028  v3(2,i,j) = ua(i,j,k)*vlon(i,j,2) + va(i,j,k)*vlat(i,j,2)
3029  v3(3,i,j) = ua(i,j,k)*vlon(i,j,3) + va(i,j,k)*vlat(i,j,3)
3030  enddo
3031  enddo
3032 
3033 ! A --> D
3034 ! Interpolate to cell edges
3035  do j=js,je+1
3036  do i=is-1,ie+1
3037  ue(1,i,j) = 0.5*(v3(1,i,j-1) + v3(1,i,j))
3038  ue(2,i,j) = 0.5*(v3(2,i,j-1) + v3(2,i,j))
3039  ue(3,i,j) = 0.5*(v3(3,i,j-1) + v3(3,i,j))
3040  enddo
3041  enddo
3042 
3043  do j=js-1,je+1
3044  do i=is,ie+1
3045  ve(1,i,j) = 0.5*(v3(1,i-1,j) + v3(1,i,j))
3046  ve(2,i,j) = 0.5*(v3(2,i-1,j) + v3(2,i,j))
3047  ve(3,i,j) = 0.5*(v3(3,i-1,j) + v3(3,i,j))
3048  enddo
3049  enddo
3050 
3051 ! --- E_W edges (for v-wind):
3052  if (.not. gridstruct%nested) then
3053  if ( is==1) then
3054  i = 1
3055  do j=js,je
3056  if ( j>jm2 ) then
3057  vt1(j) = edge_vect_w(j)*ve(1,i,j-1)+(1.-edge_vect_w(j))*ve(1,i,j)
3058  vt2(j) = edge_vect_w(j)*ve(2,i,j-1)+(1.-edge_vect_w(j))*ve(2,i,j)
3059  vt3(j) = edge_vect_w(j)*ve(3,i,j-1)+(1.-edge_vect_w(j))*ve(3,i,j)
3060  else
3061  vt1(j) = edge_vect_w(j)*ve(1,i,j+1)+(1.-edge_vect_w(j))*ve(1,i,j)
3062  vt2(j) = edge_vect_w(j)*ve(2,i,j+1)+(1.-edge_vect_w(j))*ve(2,i,j)
3063  vt3(j) = edge_vect_w(j)*ve(3,i,j+1)+(1.-edge_vect_w(j))*ve(3,i,j)
3064  endif
3065  enddo
3066  do j=js,je
3067  ve(1,i,j) = vt1(j)
3068  ve(2,i,j) = vt2(j)
3069  ve(3,i,j) = vt3(j)
3070  enddo
3071  endif
3072 
3073  if ( (ie+1)==npx ) then
3074  i = npx
3075  do j=js,je
3076  if ( j>jm2 ) then
3077  vt1(j) = edge_vect_e(j)*ve(1,i,j-1)+(1.-edge_vect_e(j))*ve(1,i,j)
3078  vt2(j) = edge_vect_e(j)*ve(2,i,j-1)+(1.-edge_vect_e(j))*ve(2,i,j)
3079  vt3(j) = edge_vect_e(j)*ve(3,i,j-1)+(1.-edge_vect_e(j))*ve(3,i,j)
3080  else
3081  vt1(j) = edge_vect_e(j)*ve(1,i,j+1)+(1.-edge_vect_e(j))*ve(1,i,j)
3082  vt2(j) = edge_vect_e(j)*ve(2,i,j+1)+(1.-edge_vect_e(j))*ve(2,i,j)
3083  vt3(j) = edge_vect_e(j)*ve(3,i,j+1)+(1.-edge_vect_e(j))*ve(3,i,j)
3084  endif
3085  enddo
3086  do j=js,je
3087  ve(1,i,j) = vt1(j)
3088  ve(2,i,j) = vt2(j)
3089  ve(3,i,j) = vt3(j)
3090  enddo
3091  endif
3092 
3093 ! N-S edges (for u-wind):
3094  if ( js==1 ) then
3095  j = 1
3096  do i=is,ie
3097  if ( i>im2 ) then
3098  ut1(i) = edge_vect_s(i)*ue(1,i-1,j)+(1.-edge_vect_s(i))*ue(1,i,j)
3099  ut2(i) = edge_vect_s(i)*ue(2,i-1,j)+(1.-edge_vect_s(i))*ue(2,i,j)
3100  ut3(i) = edge_vect_s(i)*ue(3,i-1,j)+(1.-edge_vect_s(i))*ue(3,i,j)
3101  else
3102  ut1(i) = edge_vect_s(i)*ue(1,i+1,j)+(1.-edge_vect_s(i))*ue(1,i,j)
3103  ut2(i) = edge_vect_s(i)*ue(2,i+1,j)+(1.-edge_vect_s(i))*ue(2,i,j)
3104  ut3(i) = edge_vect_s(i)*ue(3,i+1,j)+(1.-edge_vect_s(i))*ue(3,i,j)
3105  endif
3106  enddo
3107  do i=is,ie
3108  ue(1,i,j) = ut1(i)
3109  ue(2,i,j) = ut2(i)
3110  ue(3,i,j) = ut3(i)
3111  enddo
3112  endif
3113 
3114  if ( (je+1)==npy ) then
3115  j = npy
3116  do i=is,ie
3117  if ( i>im2 ) then
3118  ut1(i) = edge_vect_n(i)*ue(1,i-1,j)+(1.-edge_vect_n(i))*ue(1,i,j)
3119  ut2(i) = edge_vect_n(i)*ue(2,i-1,j)+(1.-edge_vect_n(i))*ue(2,i,j)
3120  ut3(i) = edge_vect_n(i)*ue(3,i-1,j)+(1.-edge_vect_n(i))*ue(3,i,j)
3121  else
3122  ut1(i) = edge_vect_n(i)*ue(1,i+1,j)+(1.-edge_vect_n(i))*ue(1,i,j)
3123  ut2(i) = edge_vect_n(i)*ue(2,i+1,j)+(1.-edge_vect_n(i))*ue(2,i,j)
3124  ut3(i) = edge_vect_n(i)*ue(3,i+1,j)+(1.-edge_vect_n(i))*ue(3,i,j)
3125  endif
3126  enddo
3127  do i=is,ie
3128  ue(1,i,j) = ut1(i)
3129  ue(2,i,j) = ut2(i)
3130  ue(3,i,j) = ut3(i)
3131  enddo
3132  endif
3133 
3134  endif ! .not. nested
3135  do j=js,je+1
3136  do i=is,ie
3137  u(i,j,k) = ue(1,i,j)*es(1,i,j,1) + &
3138  ue(2,i,j)*es(2,i,j,1) + &
3139  ue(3,i,j)*es(3,i,j,1)
3140  enddo
3141  enddo
3142  do j=js,je
3143  do i=is,ie+1
3144  v(i,j,k) = ve(1,i,j)*ew(1,i,j,2) + &
3145  ve(2,i,j)*ew(2,i,j,2) + &
3146  ve(3,i,j)*ew(3,i,j,2)
3147  enddo
3148  enddo
3149 
3150  enddo ! k-loop
3151 
3152  end subroutine cubed_a2d
3153 
3154 
3155  subroutine d2a3d(u, v, ua, va, im, jm, km, lon)
3156  integer, intent(in):: im, jm, km ! Dimensions
3157  real(FVPRC), intent(in ) :: lon(im)
3158  real(FVPRC), intent(in ), dimension(im,jm,km):: u, v
3159  real(FVPRC), intent(out), dimension(im,jm,km):: ua, va
3160 ! local
3161  real(FVPRC) :: coslon(im),sinlon(im) ! Sine and cosine in longitude
3162  integer i, j, k
3163  integer imh
3164  real(FVPRC) un, vn, us, vs
3165 
3166  integer :: ks, ke
3167 
3168  imh = im/2
3169 
3170  do i=1,im
3171  sinlon(i) = sin(lon(i))
3172  coslon(i) = cos(lon(i))
3173  enddo
3174 
3175  do k=1,km
3176  do j=2,jm-1
3177  do i=1,im
3178  ua(i,j,k) = 0.5*(u(i,j,k) + u(i,j+1,k))
3179  enddo
3180  enddo
3181 
3182  do j=2,jm-1
3183  do i=1,im-1
3184  va(i,j,k) = 0.5*(v(i,j,k) + v(i+1,j,k))
3185  enddo
3186  va(im,j,k) = 0.5*(v(im,j,k) + v(1,j,k))
3187  enddo
3188 
3189 ! Projection at SP
3190  us = 0.
3191  vs = 0.
3192  do i=1,imh
3193  us = us + (ua(i+imh,2,k)-ua(i,2,k))*sinlon(i) &
3194  + (va(i,2,k)-va(i+imh,2,k))*coslon(i)
3195  vs = vs + (ua(i+imh,2,k)-ua(i,2,k))*coslon(i) &
3196  + (va(i+imh,2,k)-va(i,2,k))*sinlon(i)
3197  enddo
3198  us = us/im
3199  vs = vs/im
3200  do i=1,imh
3201  ua(i,1,k) = -us*sinlon(i) - vs*coslon(i)
3202  va(i,1,k) = us*coslon(i) - vs*sinlon(i)
3203  ua(i+imh,1,k) = -ua(i,1,k)
3204  va(i+imh,1,k) = -va(i,1,k)
3205  enddo
3206 
3207 ! Projection at NP
3208  un = 0.
3209  vn = 0.
3210  do i=1,imh
3211  un = un + (ua(i+imh,jm-1,k)-ua(i,jm-1,k))*sinlon(i) &
3212  + (va(i+imh,jm-1,k)-va(i,jm-1,k))*coslon(i)
3213  vn = vn + (ua(i,jm-1,k)-ua(i+imh,jm-1,k))*coslon(i) &
3214  + (va(i+imh,jm-1,k)-va(i,jm-1,k))*sinlon(i)
3215  enddo
3216 
3217  un = un/im
3218  vn = vn/im
3219  do i=1,imh
3220  ua(i,jm,k) = -un*sinlon(i) + vn*coslon(i)
3221  va(i,jm,k) = -un*coslon(i) - vn*sinlon(i)
3222  ua(i+imh,jm,k) = -ua(i,jm,k)
3223  va(i+imh,jm,k) = -va(i,jm,k)
3224  enddo
3225  enddo
3226 
3227  end subroutine d2a3d
3228 
3229 
3230 
3231  subroutine pmaxmin( qname, a, im, jm, fac )
3232 
3233  integer, intent(in):: im, jm
3234  character(len=*) :: qname
3235  integer i, j
3236  real(FVPRC) a(im,jm)
3237 
3238  real(FVPRC) qmin(jm), qmax(jm)
3239  real(FVPRC) pmax, pmin
3240  real(FVPRC) fac ! multiplication factor
3241 
3242  do j=1,jm
3243  pmax = a(1,j)
3244  pmin = a(1,j)
3245  do i=2,im
3246  pmax = max(pmax, a(i,j))
3247  pmin = min(pmin, a(i,j))
3248  enddo
3249  qmax(j) = pmax
3250  qmin(j) = pmin
3251  enddo
3252 !
3253 ! Now find max/min of amax/amin
3254 !
3255  pmax = qmax(1)
3256  pmin = qmin(1)
3257  do j=2,jm
3258  pmax = max(pmax, qmax(j))
3259  pmin = min(pmin, qmin(j))
3260  enddo
3261 
3262  write(*,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac
3263 
3264  end subroutine pmaxmin
3265 
3266  subroutine pmaxmin4d( qname, a, im, jm, km, lm, fac )
3267 
3268  character*(*) qname
3269  integer, intent(in):: im, jm, km, lm
3270  integer i, j, k, l
3271  real(FVPRC) a(im,jm,km,lm)
3272 
3273  real(FVPRC) qmin(jm), qmax(jm)
3274  real(FVPRC) pmax, pmin
3275  real(FVPRC) fac ! multiplication factor
3276 
3277  pmax = a(1,1,1,1)
3278  pmin = a(1,1,1,1)
3279  do l=1,lm
3280  do k=1,km
3281  do j=1,jm
3282  do i=1,im
3283  pmax = max(pmax, a(i,j,k,l))
3284  pmin = min(pmin, a(i,j,k,l))
3285  enddo
3286  enddo
3287  enddo
3288  enddo
3289 
3290  write(*,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac
3291 
3292  end subroutine pmaxmin4d
3293 
3294 
3295  subroutine mpp_domain_decomp(domain,npx,npy,nregions,ng,grid_type, &
3296  is,ie,js,je,isd,ied,jsd,jed,tile)
3297  use mpp_domains_mod, only: mpp_domains_init, mpp_domain_time, &
3298  mpp_define_mosaic, mpp_get_compute_domain, mpp_get_data_domain, &
3299  mpp_domains_set_stack_size, mpp_define_layout
3300  use mpp_mod, only : mpp_pe
3301  type(domain2D), intent(OUT) :: domain
3302  integer, intent(IN) :: npx,npy,nregions,ng,grid_type
3303  integer, intent(OUT) :: is,ie,js,je,isd,ied,jsd,jed,tile
3304 
3305  integer :: layout(2)
3306  integer, allocatable :: pe_start(:), pe_end(:)
3307 
3308  integer :: num_contact, ntiles, npes_per_tile
3309  integer, allocatable, dimension(:) :: npes_tile, tile1, tile2
3310  integer, allocatable, dimension(:) :: istart1, iend1, jstart1, jend1
3311  integer, allocatable, dimension(:) :: istart2, iend2, jstart2, jend2
3312  integer, allocatable, dimension(:,:) :: layout2D, global_indices
3313 
3314  character*80 :: evalue
3315  integer :: ios,nx,ny,n,num_alloc
3316  character(len=32) :: type = "unknown"
3317  logical :: is_symmetry
3318  logical :: debug=.false.
3319  integer :: npes
3320 
3321  npes = mpp_npes()
3322 
3323  nx = npx-1
3324  ny = npy-1
3325 
3326  call print_memuse_stats('external_ic:mpp_domain_decomp: top')
3327 
3328  call mpp_domains_init(mpp_domain_time)
3329 
3330  call print_memuse_stats('external_ic:mpp_domain_decomp: mpp_domains_init')
3331 
3332 ! call mpp_domains_set_stack_size(10000)
3333 ! call mpp_domains_set_stack_size(900000)
3334 ! call mpp_domains_set_stack_size(1500000)
3335  call mpp_domains_set_stack_size(3000000)
3336 
3337  select case(nregions)
3338  case ( 1 ) ! Lat-Lon "cyclic"
3339 
3340  select case (grid_type)
3341  case (3) ! Lat-Lon "cyclic"
3342  type="Lat-Lon: cyclic"
3343  ntiles = 4
3344  num_contact = 8
3345  if( mod(npes,ntiles) .NE. 0 ) then
3346  call mpp_error(note,'TEST_MPP_DOMAINS: for Cyclic mosaic, npes should be multiple of ntiles. ' // &
3347  'No test is done for Cyclic mosaic. ' )
3348  return
3349  end if
3350  npes_per_tile = npes/ntiles
3351  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3352  layout = (/1,npes_per_tile/) ! force decomp only in Lat-Direction
3353  case (4) ! Cartesian, double periodic
3354  type="Cartesian: double periodic"
3355  ntiles = 1
3356  num_contact = 2
3357  npes_per_tile = npes/ntiles
3358  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3359  case (5) ! latlon patch
3360  type="Lat-Lon: patch"
3361  ntiles = 1
3362  num_contact = 0
3363  npes_per_tile = npes/ntiles
3364  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3365  case (6) ! latlon strip
3366  type="Lat-Lon: strip"
3367  ntiles = 1
3368  num_contact = 1
3369  npes_per_tile = npes/ntiles
3370  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3371  case (7) ! Cartesian, channel
3372  type="Cartesian: channel"
3373  ntiles = 1
3374  num_contact = 1
3375  npes_per_tile = npes/ntiles
3376  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3377  end select
3378 
3379  case ( 6 ) ! Cubed-Sphere
3380  type="Cubic: cubed-sphere"
3381  ntiles = 6
3382  num_contact = 12
3383 !--- cubic grid always have six tiles, so npes should be multiple of 6
3384  if( mod(npes,ntiles) .NE. 0 .OR. npx-1 .NE. npy-1) then
3385  call mpp_error(note,'mpp_domain_decomp: for Cubic_grid mosaic, npes should be multiple of ntiles(6) ' // &
3386  'and npx-1 should equal npy-1, mpp_domain_decomp is NOT done for Cubic-grid mosaic. ' )
3387  return
3388  end if
3389  npes_per_tile = npes/ntiles
3390  call mpp_define_layout( (/1,npx-1,1,npy-1/), npes_per_tile, layout )
3391  npes_x = layout(1)
3392  npes_y = layout(2)
3393  if ( (npx/npes_x < ng) .or. (npy/npes_y < ng) ) then
3394  write(*,310) npes_x, npes_y, npx/npes_x, npy/npes_y
3395 310 format('Invalid layout, NPES_X:',i4.4,'NPES_Y:',i4.4,'ncells_X:',i4.4,'ncells_Y:',i4.4)
3396  call mpp_error(fatal, 'mpp_domain_decomp: bad decomp')
3397  endif
3398  case default
3399  call mpp_error(fatal, 'mpp_domain_decomp: no such test: '//type)
3400  end select
3401 
3402  call print_memuse_stats('external_ic:mpp_domain_decomp: mpp_define_layout')
3403 
3404  allocate(layout2d(2,ntiles), global_indices(4,ntiles), npes_tile(ntiles) )
3405  allocate(pe_start(ntiles),pe_end(ntiles))
3406  npes_tile = npes_per_tile
3407  do n = 1, ntiles
3408  global_indices(:,n) = (/1,npx-1,1,npy-1/)
3409  layout2d(:,n) = layout
3410  pe_start(n) = (n-1)*layout(1)*layout(2)
3411  pe_end(n) = pe_start(n) + layout(1)*layout(2) -1
3412  end do
3413  num_alloc=max(1,num_contact)
3414  allocate(tile1(num_alloc), tile2(num_alloc) )
3415  allocate(istart1(num_alloc), iend1(num_alloc), jstart1(num_alloc), jend1(num_alloc) )
3416  allocate(istart2(num_alloc), iend2(num_alloc), jstart2(num_alloc), jend2(num_alloc) )
3417 
3418  is_symmetry = .true.
3419 
3420  call print_memuse_stats('external_ic:mpp_domain_decomp: allocates 1')
3421 
3422  select case(nregions)
3423  case ( 1 )
3424 
3425  select case (grid_type)
3426  case (3) ! Lat-Lon "cyclic"
3427 !--- Contact line 1, between tile 1 (EAST) and tile 2 (WEST)
3428  tile1(1) = 1; tile2(1) = 2
3429  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
3430  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
3431 !--- Contact line 2, between tile 1 (SOUTH) and tile 3 (NORTH) --- cyclic
3432  tile1(2) = 1; tile2(2) = 3
3433  istart1(2) = 1; iend1(2) = nx; jstart1(2) = 1; jend1(2) = 1
3434  istart2(2) = 1; iend2(2) = nx; jstart2(2) = ny; jend2(2) = ny
3435 !--- Contact line 3, between tile 1 (WEST) and tile 2 (EAST) --- cyclic
3436  tile1(3) = 1; tile2(3) = 2
3437  istart1(3) = 1; iend1(3) = 1; jstart1(3) = 1; jend1(3) = ny
3438  istart2(3) = nx; iend2(3) = nx; jstart2(3) = 1; jend2(3) = ny
3439 !--- Contact line 4, between tile 1 (NORTH) and tile 3 (SOUTH)
3440  tile1(4) = 1; tile2(4) = 3
3441  istart1(4) = 1; iend1(4) = nx; jstart1(4) = ny; jend1(4) = ny
3442  istart2(4) = 1; iend2(4) = nx; jstart2(4) = 1; jend2(4) = 1
3443 !--- Contact line 5, between tile 2 (SOUTH) and tile 4 (NORTH) --- cyclic
3444  tile1(5) = 2; tile2(5) = 4
3445  istart1(5) = 1; iend1(5) = nx; jstart1(5) = 1; jend1(5) = 1
3446  istart2(5) = 1; iend2(5) = nx; jstart2(5) = ny; jend2(5) = ny
3447 !--- Contact line 6, between tile 2 (NORTH) and tile 4 (SOUTH)
3448  tile1(6) = 2; tile2(6) = 4
3449  istart1(6) = 1; iend1(6) = nx; jstart1(6) = ny; jend1(6) = ny
3450  istart2(6) = 1; iend2(6) = nx; jstart2(6) = 1; jend2(6) = 1
3451 !--- Contact line 7, between tile 3 (EAST) and tile 4 (WEST)
3452  tile1(7) = 3; tile2(7) = 4
3453  istart1(7) = nx; iend1(7) = nx; jstart1(7) = 1; jend1(7) = ny
3454  istart2(7) = 1; iend2(7) = 1; jstart2(7) = 1; jend2(7) = ny
3455 !--- Contact line 8, between tile 3 (WEST) and tile 4 (EAST) --- cyclic
3456  tile1(8) = 3; tile2(8) = 4
3457  istart1(8) = 1; iend1(8) = 1; jstart1(8) = 1; jend1(8) = ny
3458  istart2(8) = nx; iend2(8) = nx; jstart2(8) = 1; jend2(8) = ny
3459  is_symmetry = .false.
3460  case (4) ! Cartesian, double periodic
3461 !--- Contact line 1, between tile 1 (EAST) and tile 1 (WEST)
3462  tile1(1) = 1; tile2(1) = 1
3463  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
3464  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
3465 !--- Contact line 2, between tile 1 (SOUTH) and tile 1 (NORTH) --- cyclic
3466  tile1(2) = 1; tile2(2) = 1
3467  istart1(2) = 1; iend1(2) = nx; jstart1(2) = 1; jend1(2) = 1
3468  istart2(2) = 1; iend2(2) = nx; jstart2(2) = ny; jend2(2) = ny
3469  case (5) ! latlon patch
3470 
3471  case (6) !latlon strip
3472 !--- Contact line 1, between tile 1 (EAST) and tile 1 (WEST)
3473  tile1(1) = 1; tile2(1) = 1
3474  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
3475  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
3476  case (7) ! Cartesian, channel
3477 !--- Contact line 1, between tile 1 (EAST) and tile 1 (WEST)
3478  tile1(1) = 1; tile2(1) = 1
3479  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
3480  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
3481  end select
3482 
3483  case ( 6 ) ! Cubed-Sphere
3484 !--- Contact line 1, between tile 1 (EAST) and tile 2 (WEST)
3485  tile1(1) = 1; tile2(1) = 2
3486  istart1(1) = nx; iend1(1) = nx; jstart1(1) = 1; jend1(1) = ny
3487  istart2(1) = 1; iend2(1) = 1; jstart2(1) = 1; jend2(1) = ny
3488 !--- Contact line 2, between tile 1 (NORTH) and tile 3 (WEST)
3489  tile1(2) = 1; tile2(2) = 3
3490  istart1(2) = 1; iend1(2) = nx; jstart1(2) = ny; jend1(2) = ny
3491  istart2(2) = 1; iend2(2) = 1; jstart2(2) = ny; jend2(2) = 1
3492 !--- Contact line 3, between tile 1 (WEST) and tile 5 (NORTH)
3493  tile1(3) = 1; tile2(3) = 5
3494  istart1(3) = 1; iend1(3) = 1; jstart1(3) = 1; jend1(3) = ny
3495  istart2(3) = nx; iend2(3) = 1; jstart2(3) = ny; jend2(3) = ny
3496 !--- Contact line 4, between tile 1 (SOUTH) and tile 6 (NORTH)
3497  tile1(4) = 1; tile2(4) = 6
3498  istart1(4) = 1; iend1(4) = nx; jstart1(4) = 1; jend1(4) = 1
3499  istart2(4) = 1; iend2(4) = nx; jstart2(4) = ny; jend2(4) = ny
3500 !--- Contact line 5, between tile 2 (NORTH) and tile 3 (SOUTH)
3501  tile1(5) = 2; tile2(5) = 3
3502  istart1(5) = 1; iend1(5) = nx; jstart1(5) = ny; jend1(5) = ny
3503  istart2(5) = 1; iend2(5) = nx; jstart2(5) = 1; jend2(5) = 1
3504 !--- Contact line 6, between tile 2 (EAST) and tile 4 (SOUTH)
3505  tile1(6) = 2; tile2(6) = 4
3506  istart1(6) = nx; iend1(6) = nx; jstart1(6) = 1; jend1(6) = ny
3507  istart2(6) = nx; iend2(6) = 1; jstart2(6) = 1; jend2(6) = 1
3508 !--- Contact line 7, between tile 2 (SOUTH) and tile 6 (EAST)
3509  tile1(7) = 2; tile2(7) = 6
3510  istart1(7) = 1; iend1(7) = nx; jstart1(7) = 1; jend1(7) = 1
3511  istart2(7) = nx; iend2(7) = nx; jstart2(7) = ny; jend2(7) = 1
3512 !--- Contact line 8, between tile 3 (EAST) and tile 4 (WEST)
3513  tile1(8) = 3; tile2(8) = 4
3514  istart1(8) = nx; iend1(8) = nx; jstart1(8) = 1; jend1(8) = ny
3515  istart2(8) = 1; iend2(8) = 1; jstart2(8) = 1; jend2(8) = ny
3516 !--- Contact line 9, between tile 3 (NORTH) and tile 5 (WEST)
3517  tile1(9) = 3; tile2(9) = 5
3518  istart1(9) = 1; iend1(9) = nx; jstart1(9) = ny; jend1(9) = ny
3519  istart2(9) = 1; iend2(9) = 1; jstart2(9) = ny; jend2(9) = 1
3520 !--- Contact line 10, between tile 4 (NORTH) and tile 5 (SOUTH)
3521  tile1(10) = 4; tile2(10) = 5
3522  istart1(10) = 1; iend1(10) = nx; jstart1(10) = ny; jend1(10) = ny
3523  istart2(10) = 1; iend2(10) = nx; jstart2(10) = 1; jend2(10) = 1
3524 !--- Contact line 11, between tile 4 (EAST) and tile 6 (SOUTH)
3525  tile1(11) = 4; tile2(11) = 6
3526  istart1(11) = nx; iend1(11) = nx; jstart1(11) = 1; jend1(11) = ny
3527  istart2(11) = nx; iend2(11) = 1; jstart2(11) = 1; jend2(11) = 1
3528 !--- Contact line 12, between tile 5 (EAST) and tile 6 (WEST)
3529  tile1(12) = 5; tile2(12) = 6
3530  istart1(12) = nx; iend1(12) = nx; jstart1(12) = 1; jend1(12) = ny
3531  istart2(12) = 1; iend2(12) = 1; jstart2(12) = 1; jend2(12) = ny
3532  end select
3533 
3534  call mpp_define_mosaic(global_indices, layout2d, domain, ntiles, num_contact, tile1, tile2, &
3535  istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2, &
3536  pe_start=pe_start, pe_end=pe_end, symmetry=is_symmetry, &
3537  shalo = ng, nhalo = ng, whalo = ng, ehalo = ng, name = type)
3538  call print_memuse_stats('external_ic:mpp_domain_decomp: mpp_define_mosaic')
3539 
3540  deallocate(pe_start,pe_end)
3541 
3542 !--- find the tile number
3543  tile = mpp_pe()/npes_per_tile+1
3544  call mpp_get_compute_domain( domain, is, ie, js, je )
3545  call mpp_get_data_domain ( domain, isd, ied, jsd, jed )
3546 
3547  call print_memuse_stats('external_ic:mpp_domain_decomp: mpp_get domains')
3548 
3549  end subroutine mpp_domain_decomp
3550 !
3551 ! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ !
3552 !-------------------------------------------------------------------------------
3553 
3554  subroutine parallel_read_file_r8(fname, npts, is,ie, js,je, km, offset, var)
3555  character(len=*), intent(IN) :: fname
3556  integer, intent(IN) :: npts, is,ie, js,je, km
3557  integer (kind=MPI_OFFSET_KIND), intent(INOUT) :: offset
3558  real(FVPRC), intent(INOUT) :: var(is:ie, js:je, km)
3559 
3560  integer :: ntiles=6
3561  real(REAL64) :: var_r8(is:ie, js:je)
3562  integer :: k
3563 
3564  integer :: MUNIT=17
3565  integer :: lsize, gsizes(2), distribs(2), dargs(2), psizes(2)
3566  integer :: filetype
3567  integer :: mcol, mrow, irow, jcol, mpiio_rank
3568  integer :: rank, total_pes
3569  integer :: mpistatus(MPI_STATUS_SIZE)
3570  integer (kind=MPI_OFFSET_KIND) :: slice_2d
3571 
3572  real(FVPRC) :: xmod, ymod
3573  character(128) :: strErr
3574 
3575  xmod = mod(npts,npes_x)
3576  write(strerr, "(i4.4,' not evenly divisible by ',i4.4)") npts, npes_x
3577  if (xmod /= 0) call mpp_error(fatal, strerr)
3578  ymod = mod(npts*6,npes_y)
3579  write(strerr, "(i4.4,' not evenly divisible by ',i4.4)") npts*6, npes_y
3580  if (ymod /= 0) call mpp_error(fatal, strerr)
3581 
3582  call mpi_file_open(mpi_comm_world, fname, mpi_mode_rdonly, mpi_info_null, munit, status)
3583  gsizes(1) = npts
3584  gsizes(2) = npts * 6
3585  distribs(1) = mpi_distribute_block
3586  distribs(2) = mpi_distribute_block
3587  dargs(1) = mpi_distribute_dflt_darg
3588  dargs(2) = mpi_distribute_dflt_darg
3589  psizes(1) = npes_x
3590  psizes(2) = npes_y * 6
3591  call mpi_comm_size(mpi_comm_world, total_pes, status)
3592  call mpi_comm_rank(mpi_comm_world, rank, status)
3593  mcol = npes_x
3594  mrow = npes_y*ntiles
3595  irow = rank/mcol !! logical row number
3596  jcol = mod(rank, mcol) !! logical column number
3597  mpiio_rank = jcol*mrow + irow
3598  call mpi_type_create_darray(total_pes, mpiio_rank, 2, gsizes, distribs, dargs, psizes, mpi_order_fortran, mpi_double_precision, filetype, status)
3599  call mpi_type_commit(filetype, status)
3600  lsize = (ie-is+1)*(je-js+1)
3601  slice_2d = npts*npts*ntiles
3602  do k=1,km
3603  call mpi_file_set_view(munit, offset, mpi_double_precision, filetype, "native", mpi_info_null, status)
3604  call mpi_file_read_all(munit, var_r8, lsize, mpi_double_precision, mpistatus, status)
3605  var(:,:,k) = var_r8
3606  offset = offset + slice_2d*8 + 8
3607  enddo
3608  call mpi_file_close(munit, status)
3609 
3610  end subroutine parallel_read_file_r8
3611 
3612  subroutine parallel_read_file_r4(fname, npts, is,ie, js,je, km, offset, var)
3613  character(len=*), intent(IN) :: fname
3614  integer, intent(IN) :: npts, is,ie, js,je, km
3615  integer (kind=MPI_OFFSET_KIND), intent(INOUT) :: offset
3616  real(FVPRC), intent(INOUT) :: var(is:ie, js:je, km)
3617 
3618  integer :: ntiles=6
3619  real(REAL4) :: var_r4(is:ie, js:je)
3620  integer :: k
3621 
3622  integer :: MUNIT=17
3623  integer :: lsize, gsizes(2), distribs(2), dargs(2), psizes(2)
3624  integer :: filetype
3625  integer :: mcol, mrow, irow, jcol, mpiio_rank
3626  integer :: rank, total_pes
3627  integer :: mpistatus(MPI_STATUS_SIZE)
3628  integer (kind=MPI_OFFSET_KIND) :: slice_2d
3629 
3630  real(FVPRC) :: xmod, ymod
3631  character(128) :: strErr
3632 
3633  xmod = mod(npts,npes_x)
3634  write(strerr, "(i4.4,' not evenly divisible by ',i4.4)") npts, npes_x
3635  if (xmod /= 0) call mpp_error(fatal, strerr)
3636  ymod = mod(npts*6,npes_y)
3637  write(strerr, "(i4.4,' not evenly divisible by ',i4.4)") npts*6, npes_y
3638  if (ymod /= 0) call mpp_error(fatal, strerr)
3639 
3640  call mpi_file_open(mpi_comm_world, fname, mpi_mode_rdonly, mpi_info_null, munit, status)
3641  gsizes(1) = npts
3642  gsizes(2) = npts * 6
3643  distribs(1) = mpi_distribute_block
3644  distribs(2) = mpi_distribute_block
3645  dargs(1) = mpi_distribute_dflt_darg
3646  dargs(2) = mpi_distribute_dflt_darg
3647  psizes(1) = npes_x
3648  psizes(2) = npes_y * 6
3649  call mpi_comm_size(mpi_comm_world, total_pes, status)
3650  call mpi_comm_rank(mpi_comm_world, rank, status)
3651  mcol = npes_x
3652  mrow = npes_y*ntiles
3653  irow = rank/mcol !! logical row number
3654  jcol = mod(rank, mcol) !! logical column number
3655  mpiio_rank = jcol*mrow + irow
3656  call mpi_type_create_darray(total_pes, mpiio_rank, 2, gsizes, distribs, dargs, psizes, mpi_order_fortran, mpi_real, filetype, status)
3657  call mpi_type_commit(filetype, status)
3658  lsize = (ie-is+1)*(je-js+1)
3659  slice_2d = npts*npts*ntiles
3660  do k=1,km
3661  call mpi_file_set_view(munit, offset, mpi_real, filetype, "native", mpi_info_null, status)
3662  call mpi_file_read_all(munit, var_r4, lsize, mpi_real, mpistatus, status)
3663  var(:,:,k) = var_r4
3664  offset = offset + slice_2d*4 + 8
3665  enddo
3666  call mpi_file_close(munit, status)
3667 
3668  end subroutine parallel_read_file_r4
3669 
3670  end module external_ic_nlm_mod
3671 
memutils_mod::print_memuse_stats
subroutine, public print_memuse_stats(text, unit, always)
Definition: memutils.F90:282
fms_mod
Definition: fms.F90:20
external_ic_nlm_mod::iunit
integer iunit
Definition: external_ic_nlm.F90:74
tracer_manager_mod::get_tracer_index
Definition: tracer_manager.F90:124
fv_mp_nlm_mod
Definition: fv_mp_nlm_mod.F90:24
mapl_constantsmod::mapl_omega
real, parameter, public mapl_omega
Definition: MAPL_Constants.F90:20
constants_mod
Definition: constants.F90:24
fv_arrays_nlm_mod::real8
integer, parameter real8
Definition: fv_arrays_nlm.F90:34
mpp_parameter_mod
Definition: mpp_parameter.F90:19
field_manager_mod::model_atmos
integer, parameter, public model_atmos
Definition: field_manager.F90:299
init_hydro_nlm_mod::p_var
subroutine, public p_var(km, ifirst, ilast, jfirst, jlast, ptop, ptop_min, delp, delz, pt, ps, pe, peln, pk, pkz, cappa, q, ng, nq, area, dry_mass, adjust_dry_mass, mountain, moist_phys, hydrostatic, nwat, domain, make_nh)
Definition: init_hydro_nlm.F90:47
constants_mod::omega
real, parameter, public omega
Rotation rate of the Earth [1/s].
Definition: constants.F90:75
sim_nc_nlm_mod::open_ncfile
subroutine, public open_ncfile(iflnm, ncid)
Definition: sim_nc_nlm_mod.F90:46
external_ic_nlm_mod::mpp_domain_decomp
subroutine mpp_domain_decomp(domain, npx, npy, nregions, ng, grid_type, is, ie, js, je, isd, ied, jsd, jed, tile)
Definition: external_ic_nlm.F90:3297
field_manager_mod
Definition: field_manager.F90:20
fv_grid_utils_nlm_mod::ptop_min
real, parameter, public ptop_min
Definition: fv_grid_utils_nlm.F90:54
sim_nc_nlm_mod::get_var2_real
subroutine, public get_var2_real(ncid, var_name, im, jm, var2)
Definition: sim_nc_nlm_mod.F90:130
constants_mod::pi_8
real(kind=8), parameter, public pi_8
Ratio of circle circumference to diameter [N/A].
Definition: constants.F90:73
fv_arrays_nlm_mod::fv_grid_type
Definition: fv_arrays_nlm.F90:115
external_ic_nlm_mod::remap_wz
subroutine remap_wz(im, jm, km, npz, mg, ak0, bk0, psc, wa, wz, Atm)
Definition: external_ic_nlm.F90:2398
external_ic_nlm_mod::parallel_read_file_r8
subroutine parallel_read_file_r8(fname, npts, is, ie, js, je, km, offset, var)
Definition: external_ic_nlm.F90:3555
mpp_domains_mod::mpp_update_domains
Definition: mpp_domains.F90:1068
external_ic_nlm_mod::npes_y
integer, save npes_y
Definition: external_ic_nlm.F90:71
mapl_constantsmod::mapl_rvap
real, parameter, public mapl_rvap
Definition: MAPL_Constants.F90:34
fv_nwp_nudge_nlm_mod
Definition: fv_nudge_nlm.F90:26
external_ic_nlm_mod::ncep2fms
subroutine ncep2fms(im, jm, lon, lat, wk)
Definition: external_ic_nlm.F90:1983
sim_nc_nlm_mod
Definition: sim_nc_nlm_mod.F90:20
mapl_constantsmod::mapl_pi_r8
real(kind=mapl_r8), parameter, public mapl_pi_r8
Definition: MAPL_Constants.F90:13
external_ic_nlm_mod::d2a3d
subroutine d2a3d(u, v, ua, va, im, jm, km, lon)
Definition: external_ic_nlm.F90:3156
external_ic_nlm_mod::status
integer status
Definition: external_ic_nlm.F90:73
sim_nc_nlm_mod::get_var3_r4
subroutine, public get_var3_r4(ncid, var3_name, is, ie, js, je, ks, ke, var3, time_slice)
Definition: sim_nc_nlm_mod.F90:227
constants_mod::rdgas
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
Definition: constants.F90:77
external_ic_nlm_mod::remap_scalar
subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, Atm)
Definition: external_ic_nlm.F90:2162
tracer_manager_mod::get_number_tracers
subroutine, public get_number_tracers(model, num_tracers, num_prog, num_diag, num_family)
Definition: tracer_manager.F90:641
fv_arrays_nlm_mod::fv_atmos_type
Definition: fv_arrays_nlm.F90:611
fms_io_mod
Definition: fms_io.F90:20
mapl_constantsmod::mapl_cp
real, parameter, public mapl_cp
Definition: MAPL_Constants.F90:46
fv_arrays_nlm_mod::real4
integer, parameter real4
Definition: fv_arrays_nlm.F90:33
external_ic_nlm_mod::cubed_a2d
subroutine cubed_a2d(npx, npy, npz, ua, va, u, v, gridstruct, fv_domain, bd)
Definition: external_ic_nlm.F90:2969
mpp_mod::mpp_broadcast
Definition: mpp.F90:1124
external_ic_nlm_mod::get_diag_ic
subroutine get_diag_ic(Atm, fv_domain, nq)
Definition: external_ic_nlm.F90:214
external_ic_nlm_mod::deg2rad
real(fvprc) deg2rad
Definition: external_ic_nlm.F90:67
mpp_mod
Definition: mpp.F90:39
external_ic_nlm_mod::get_external_ic
subroutine, public get_external_ic(Atm, fv_domain, use_geos_latlon_restart, use_geos_cubed_restart, ntracers)
Definition: external_ic_nlm.F90:80
mpp_domains_mod::mpp_define_layout
Definition: mpp_domains.F90:742
fv_grid_utils_nlm_mod
Definition: fv_grid_utils_nlm.F90:20
fv_timing_nlm_mod
Definition: fv_timing_nlm.F90:20
external_ic_nlm_mod::get_fv_ic
subroutine get_fv_ic(Atm, fv_domain, nq)
Definition: external_ic_nlm.F90:1836
external_ic_nlm_mod
Definition: external_ic_nlm.F90:1
external_ic_nlm_mod::parallel_read_file_r4
subroutine parallel_read_file_r4(fname, npts, is, ie, js, je, km, offset, var)
Definition: external_ic_nlm.F90:3613
mpp_domains_mod
Definition: mpp_domains.F90:128
external_ic_nlm_mod::get_cubed_sphere_terrain
subroutine, public get_cubed_sphere_terrain(Atm, fv_domain)
Definition: external_ic_nlm.F90:168
mpp_parameter_mod::agrid
integer, parameter, public agrid
Definition: mpp_parameter.F90:101
external_ic_nlm_mod::ounit
integer ounit
Definition: external_ic_nlm.F90:75
fv_nwp_nudge_nlm_mod::t_is_tv
logical, public t_is_tv
Definition: fv_nudge_nlm.F90:114
mpp_domains_mod::mpp_get_data_domain
Definition: mpp_domains.F90:2617
mapl_constantsmod::mapl_kappa
real, parameter, public mapl_kappa
Definition: MAPL_Constants.F90:38
field_exist
int field_exist(const char *file, const char *name)
Definition: read_mosaic.c:69
external_ic_nlm_mod::pmaxmin4d
subroutine pmaxmin4d(qname, a, im, jm, km, lm, fac)
Definition: external_ic_nlm.F90:3267
mpp_domains_mod::domain2d
Definition: mpp_domains.F90:378
constants_mod::rvgas
real, parameter, public rvgas
Gas constant for water vapor [J/kg/deg].
Definition: constants.F90:78
fv_mp_nlm_mod::ng
integer, parameter, public ng
Definition: fv_mp_nlm_mod.F90:2465
fms_io_mod::field_size
subroutine, public field_size(filename, fieldname, siz, field_found, domain, no_domain)
Definition: fms_io.F90:4941
fv_timing_nlm_mod::timing_on
subroutine timing_on(blk_name)
Definition: fv_timing_nlm.F90:101
amip_interp_mod
Definition: amip_interp.F90:21
constants_mod::cp_air
real, parameter, public cp_air
Specific heat capacity of dry air at constant pressure [J/kg/deg].
Definition: constants.F90:83
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
fv_surf_map_nlm_mod
Definition: fv_surf_map_nlm.F90:20
external_ic_nlm_mod::npes_x
integer, save npes_x
Definition: external_ic_nlm.F90:71
fv_surf_map_nlm_mod::surfdrv
subroutine, public surfdrv(npx, npy, grid, agrid, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, phis, stretch_fac, nested, npx_global, domain, grid_number, bd)
Definition: fv_surf_map_nlm.F90:91
external_ic_nlm_mod::reverse
real(fvprc) function, dimension(size(a, 1), size(a, 2)) reverse(A)
Definition: external_ic_nlm.F90:2955
sim_nc_nlm_mod::close_ncfile
subroutine, public close_ncfile(ncid)
Definition: sim_nc_nlm_mod.F90:58
init_hydro_nlm_mod
Definition: init_hydro_nlm.F90:21
mpp_domains_mod::mpp_get_boundary
Definition: mpp_domains.F90:1900
external_ic_nlm_mod::zvir
real(fvprc), parameter zvir
Definition: external_ic_nlm.F90:66
amip_interp_mod::sst_ncep
real, dimension(:,:), allocatable, public sst_ncep
Definition: amip_interp.F90:120
constants_mod::grav
real, parameter, public grav
Acceleration due to gravity [m/s^2].
Definition: constants.F90:76
tracer_manager_mod
Definition: tracer_manager.F90:20
fv_mapz_nlm_mod
Definition: fv_mapz_nlm.F90:23
sim_nc_nlm_mod::get_var1_double
subroutine, public get_var1_double(ncid, var1_name, im, var1, var_exist)
Definition: sim_nc_nlm_mod.F90:86
fv_mapz_nlm_mod::mappm
subroutine, public mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
Definition: fv_mapz_nlm.F90:3225
sim_nc_nlm_mod::get_ncdim1
subroutine, public get_ncdim1(ncid, var1_name, im)
Definition: sim_nc_nlm_mod.F90:69
mpp_mod::mpp_error
Definition: mpp.F90:402
normalize_vect
void normalize_vect(double *e)
Definition: mosaic_util.c:679
size
************************************************************************GNU Lesser General Public License **This file is part of the GFDL Flexible Modeling System(FMS). ! *! *FMS is free software without even the implied warranty of MERCHANTABILITY or *FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License *for more details **You should have received a copy of the GNU Lesser General Public *License along with FMS If see< http:! ***********************************************************************subroutine READ_RECORD_CORE_(unit, field, nwords, data, start, axsiz) integer, intent(in) ::unit type(fieldtype), intent(in) ::field integer, intent(in) ::nwords MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in) ::start(:), axsiz(:) integer(SHORT_KIND) ::i2vals(nwords)!rab used in conjunction with transfer intrinsic to determine size of a variable integer(KIND=1) ::one_byte(8) integer ::word_sz!#ifdef __sgi integer(INT_KIND) ::ivals(nwords) real(FLOAT_KIND) ::rvals(nwords)!#else! integer ::ivals(nwords)! real ::rvals(nwords)!#endif real(DOUBLE_KIND) ::r8vals(nwords) pointer(ptr1, i2vals) pointer(ptr2, ivals) pointer(ptr3, rvals) pointer(ptr4, r8vals) if(mpp_io_stack_size< nwords) call mpp_io_set_stack_size(nwords) call mpp_error(FATAL, 'MPP_READ currently requires use_netCDF option') end subroutine READ_RECORD_CORE_ subroutine READ_RECORD_(unit, field, nwords, data, time_level, domain, position, tile_count, start_in, axsiz_in)!routine that is finally called by all mpp_read routines to perform the read!a non-netCDF record contains:! field ID! a set of 4 coordinates(is:ie, js:je) giving the data subdomain! a timelevel and a timestamp(=NULLTIME if field is static)! 3D real data(stored as 1D)!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above!in a global direct access file, record position on PE is given by %record.!Treatment of timestamp:! We assume that static fields have been passed without a timestamp.! Here that is converted into a timestamp of NULLTIME.! For non-netCDF fields, field is treated no differently, but is written! with a timestamp of NULLTIME. There is no check in the code to prevent! the user from repeatedly writing a static field. integer, intent(in) ::unit, nwords type(fieldtype), intent(in) ::field MPP_TYPE_, intent(inout) ::data(nwords) integer, intent(in), optional ::time_level type(domain2D), intent(in), optional ::domain integer, intent(in), optional ::position, tile_count integer, intent(in), optional ::start_in(:), axsiz_in(:) integer, dimension(size(field%axes(:))) ::start, axsiz integer ::tlevel !, subdomain(4) integer ::i, error, is, ie, js, je, isg, ieg, jsg, jeg type(domain2d), pointer ::io_domain=> tlevel if(PRESENT(start_in) .AND. PRESENT(axsiz_in)) then if(size(start(! the data domain and compute domain must refer to the subdomain being passed ! In this ! since that attempts to gather all data on PE size(field%axes(:)) axsiz(i)
amip_interp_mod::j_sst
integer, public j_sst
Definition: amip_interp.F90:117
max
#define max(a, b)
Definition: mosaic_util.h:33
fms_io_mod::get_tile_string
subroutine, public get_tile_string(str_out, str_in, tile, str2_in)
Definition: fms_io.F90:7735
amip_interp_mod::i_sst
integer, public i_sst
Definition: amip_interp.F90:116
tracer_manager_mod::get_tracer_names
subroutine, public get_tracer_names(model, n, name, longname, units, err_msg)
Definition: tracer_manager.F90:948
fv_io_nlm_mod::fv_io_read_tracers
subroutine, public fv_io_read_tracers(fv_domain, Atm)
Definition: fv_io_nlm.F90:187
mpp_domains_mod::mpp_get_compute_domain
Definition: mpp_domains.F90:2578
external_ic_nlm_mod::tile
integer, save tile
Definition: external_ic_nlm.F90:71
mapl_constantsmod::mapl_rgas
real, parameter, public mapl_rgas
Definition: MAPL_Constants.F90:45
mapl_constantsmod
Definition: MAPL_Constants.F90:1
mapl_constantsmod::mapl_grav
real, parameter, public mapl_grav
Definition: MAPL_Constants.F90:18
external_ic_nlm_mod::remap_coef
subroutine remap_coef(im, jm, lon, lat, id1, id2, jdc, s2c, agrid, bd)
Definition: external_ic_nlm.F90:2074
fv_diagnostics_nlm_mod::prt_maxmin
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
Definition: fv_diagnostics_nlm.F90:2771
external_ic_nlm_mod::remap_xyz
subroutine remap_xyz(im, jbeg, jend, jm, km, npz, nq, ncnst, lon, lat, ak0, bk0, ps0, gz0, ua, va, ta, qa, Atm)
Definition: external_ic_nlm.F90:2444
external_ic_nlm_mod::pmaxmin
subroutine pmaxmin(qname, a, im, jm, fac)
Definition: external_ic_nlm.F90:3232
memutils_mod
Definition: memutils.F90:19
external_ic_nlm_mod::remap_winds
subroutine remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
Definition: external_ic_nlm.F90:2334
min
#define min(a, b)
Definition: mosaic_util.h:32
external_ic_nlm_mod::get_ncep_ic
subroutine get_ncep_ic(Atm, fv_domain, nq)
Definition: external_ic_nlm.F90:1567
fv_arrays_nlm_mod::fvprc
integer, parameter fvprc
Definition: fv_arrays_nlm.F90:38
fv_diagnostics_nlm_mod
Definition: fv_diagnostics_nlm.F90:20
constants_mod::kappa
real, parameter, public kappa
RDGAS / CP_AIR [dimensionless].
Definition: constants.F90:82
fv_arrays_nlm_mod::fv_grid_bounds_type
Definition: fv_arrays_nlm.F90:601
fv_io_nlm_mod
Definition: fv_io_nlm.F90:24
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86
fv_timing_nlm_mod::timing_off
subroutine timing_off(blk_name)
Definition: fv_timing_nlm.F90:175