FV3DYCORE  Version 1.1.0
external_ic.F90
Go to the documentation of this file.
1 
2 !***********************************************************************
3 !* GNU Lesser General Public License
4 !*
5 !* This file is part of the FV3 dynamical core.
6 !*
7 !* The FV3 dynamical core is free software: you can redistribute it
8 !* and/or modify it under the terms of the
9 !* GNU Lesser General Public License as published by the
10 !* Free Software Foundation, either version 3 of the License, or
11 !* (at your option) any later version.
12 !*
13 !* The FV3 dynamical core is distributed in the hope that it will be
14 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
15 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16 !* See the GNU General Public License for more details.
17 !*
18 !* You should have received a copy of the GNU Lesser General Public
19 !* License along with the FV3 dynamical core.
20 !* If not, see <http://www.gnu.org/licenses/>.
21 !***********************************************************************
22 
23 #ifdef OVERLOAD_R4
24 #define _GET_VAR1 get_var1_real
25 #else
26 #define _GET_VAR1 get_var1_double
27 #endif
28 
31 
33 
34 ! <table>
35 ! <tr>
36 ! <th>Module Name</th>
37 ! <th>Functions Included</th>
38 ! </tr>
39 ! <tr>
40 ! <td>constants_mod</td>
41 ! <td>pi=>pi_8, omega, grav, kappa, rdgas, rvgas, cp_air</td>
42 ! </tr>
43 ! <tr>
44 ! <td>external_sst_mod</td>
45 ! <td>i_sst, j_sst, sst_ncep</td>
46 ! </tr>
47 ! <tr>
48 ! <td>field_manager_mod</td>
49 ! <td>MODEL_ATMOS</td>
50 ! </tr>
51 ! <tr>
52 ! <td>fms_mod</td>
53 ! <td>file_exist, read_data, field_exist, write_version_number,
54 ! open_namelist_file, check_nml_error, close_file,
55 ! get_mosaic_tile_file, read_data, error_mesg</td>
56 ! </tr>
57 ! <tr>
58 ! <td>fms_io_mod</td>
59 ! <td>get_tile_string, field_size, free_restart_type,
60 ! restart_file_type, register_restart_field,
61 ! save_restart, restore_state</td>
62 ! </tr>
63 ! <tr>
64 ! <td>fv_arrays_mod</td>
65 ! <td>fv_atmos_type, fv_grid_type, fv_grid_bounds_type, R_GRID</td>
66 ! </tr>
67 ! <tr>
68 ! <td>fv_control_mod</td>
69 ! <td>fv_init, fv_end, ngrids</td>
70 ! </tr>
71 ! <tr>
72 ! <td>fv_diagnostics_mod</td>
73 ! <td>prt_maxmin, prt_gb_nh_sh, prt_height</td>
74 ! </tr>
75 ! <tr>
76 ! <td>fv_eta_mod</td>
77 ! <td>set_eta, set_external_eta</td>
78 ! </tr>
79 ! <tr>
80 ! <td>fv_fill_mod</td>
81 ! <td>fillz</td>
82 ! </tr>
83 ! <tr>
84 ! <td>fv_grid_utils_mod</td>
85 ! <td>ptop_min, g_sum,mid_pt_sphere,get_unit_vect2,
86 ! get_latlon_vector,inner_prod</td>
87 ! </tr>
88 ! <tr>
89 ! <td>fv_io_mod</td>
90 ! <td>fv_io_read_tracers</td>
91 ! </tr>
92 ! <tr>
93 ! <td>fv_mp_mod</td>
94 ! <td>ng, is_master, fill_corners, YDir, mp_reduce_min, mp_reduce_max</td>
95 ! </tr>
96 ! <tr>
97 ! <td>fv_mapz_mod</td>
98 ! <td>mappm</td>
99 ! </tr>
100 ! <tr>
101 ! <td>fv_nwp_nudge_mod</td>
102 ! <td>T_is_Tv</td>
103 ! </tr>
104 ! <tr>
105 ! <td>fv_surf_map_mod</td>
106 ! <td>surfdrv, FV3_zs_filter,sgh_g, oro_g,del2_cubed_sphere, del4_cubed_sphere</td>
107 ! </tr>
108 ! <tr>
109 ! <td>fv_timing_mod</td>
110 ! <td>timing_on, timing_off</td>
111 ! </tr>
112 ! <tr>
113 ! <td>fv_update_phys_mod</td>
114 ! <td>fv_update_phys</td>
115 ! </tr>
116 ! <tr>
117 ! <td>init_hydro_mod</td>
118 ! <td>p_var</td>
119 ! </tr>
120 ! <tr>
121 ! <td>mpp_mod</td>
122 ! <td>mpp_error, FATAL, NOTE, mpp_pe, mpp_root_pe,stdlog, input_nml_file</td>
123 ! </tr>
124 ! <tr>
125 ! <td>mpp_domains_mod</td>
126 ! <td>mpp_get_tile_id, domain2d, mpp_update_domains, NORTH, EAST</td>
127 ! </tr>
128 ! <tr>
129 ! <td>mpp_parameter_mod</td>
130 ! <td>AGRID_PARAM=>AGRID</td>
131 ! </tr>
132 ! <tr>
133 ! <td>sim_nc_mod</td>
134 ! <td>open_ncfile, close_ncfile, get_ncdim1, get_var1_double, get_var2_real,
135 ! get_var3_r4, get_var2_r4, get_var1_real, get_var_att_double</td>
136 ! </tr>
137 ! <tr>
138 ! <td>tracer_manager_mod</td>
139 ! <td>get_tracer_names, get_number_tracers, get_tracer_index, set_tracer_profile</td>
140 ! </tr>
141 ! <tr>
142 ! <td>test_cases_mod</td>
143 ! <td>checker_tracers</td>
144 ! </tr>
145 ! </table>
146 
147  use netcdf
148  use external_sst_mod, only: i_sst, j_sst, sst_ncep
149  use fms_mod, only: file_exist, read_data, field_exist, write_version_number
150  use fms_mod, only: open_namelist_file, check_nml_error, close_file
151  use fms_mod, only: get_mosaic_tile_file, read_data, error_mesg
152  use fms_io_mod, only: get_tile_string, field_size, free_restart_type
153  use fms_io_mod, only: restart_file_type, register_restart_field
154  use fms_io_mod, only: save_restart, restore_state
155  use mpp_mod, only: mpp_error, fatal, note, mpp_pe, mpp_root_pe
156  use mpp_mod, only: stdlog, input_nml_file
157  use mpp_parameter_mod, only: agrid_param=>agrid
158  use mpp_domains_mod, only: mpp_get_tile_id, domain2d, mpp_update_domains, north, east
159  use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index
160  use tracer_manager_mod, only: set_tracer_profile
161  use field_manager_mod, only: model_atmos
162 
163  use constants_mod, only: pi=>pi_8, omega, grav, kappa, rdgas, rvgas, cp_air
167  use fv_io_mod, only: fv_io_read_tracers
168  use fv_mapz_mod, only: mappm
169 
171  use fv_mp_mod, only: ng, is_master, fill_corners, ydir, mp_reduce_min, mp_reduce_max
174  use fv_surf_map_mod, only: sgh_g, oro_g
177  use init_hydro_mod, only: p_var
178  use fv_fill_mod, only: fillz
182  use fv_nwp_nudge_mod, only: t_is_tv
184 
185 ! The "T" field in NCEP analysis is actually virtual temperature (Larry H. post processing)
186 ! BEFORE 20051201
187 
189  use mpp_domains_mod, only: mpp_get_data_domain, mpp_get_global_domain, mpp_get_compute_domain
190 
191 #ifdef MULTI_GASES
192  use multi_gases_mod, only: virq, virqd, vicpqd
193 #endif
194 
195  implicit none
196  private
197 
198  real, parameter:: zvir = rvgas/rdgas - 1.
199  real(kind=R_GRID), parameter :: cnst_0p20=0.20d0
200  real :: deg2rad
201  character (len = 80) :: source ! This tells what the input source was for the data
203 
204 ! version number of this module
205 ! Include variable "version" to be written to log file.
206 #include<file_version.h>
207 
208 contains
209 
210  subroutine get_external_ic( Atm, fv_domain, cold_start )
212  type(fv_atmos_type), intent(inout), target :: Atm(:)
213  type(domain2d), intent(inout) :: fv_domain
214  logical, intent(IN) :: cold_start
215  real:: alpha = 0.
216  real rdg
217  integer i,j,k,nq
218 
219  real, pointer, dimension(:,:,:) :: grid, agrid
220  real, pointer, dimension(:,:) :: fC, f0
221 
222  integer :: is, ie, js, je
223  integer :: isd, ied, jsd, jed
224  integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel
225 #ifdef CCPP
226  integer :: liq_aero, ice_aero
227 #endif
228 #ifdef MULTI_GASES
229  integer :: spfo, spfo2, spfo3
230 #else
231  integer :: o3mr
232 #endif
233 
234  is = atm(1)%bd%is
235  ie = atm(1)%bd%ie
236  js = atm(1)%bd%js
237  je = atm(1)%bd%je
238  isd = atm(1)%bd%isd
239  ied = atm(1)%bd%ied
240  jsd = atm(1)%bd%jsd
241  jed = atm(1)%bd%jed
242 
243  grid => atm(1)%gridstruct%grid
244  agrid => atm(1)%gridstruct%agrid
245 
246  fc => atm(1)%gridstruct%fC
247  f0 => atm(1)%gridstruct%f0
248 
249 ! * Initialize coriolis param:
250 
251  do j=jsd,jed+1
252  do i=isd,ied+1
253  fc(i,j) = 2.*omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(alpha) + &
254  sin(grid(i,j,2))*cos(alpha) )
255  enddo
256  enddo
257 
258  do j=jsd,jed
259  do i=isd,ied
260  f0(i,j) = 2.*omega*( -1.*cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(alpha) + &
261  sin(agrid(i,j,2))*cos(alpha) )
262  enddo
263  enddo
264 
265  call mpp_update_domains( f0, fv_domain )
266  if ( atm(1)%gridstruct%cubed_sphere .and. (.not. (atm(1)%neststruct%nested .or. atm(1)%flagstruct%regional)))then
267  call fill_corners(f0, atm(1)%npx, atm(1)%npy, ydir)
268  endif
269 
270 ! Read in cubed_sphere terrain
271  if ( atm(1)%flagstruct%mountain ) then
272  call get_cubed_sphere_terrain(atm, fv_domain)
273  else
274  if (.not. atm(1)%neststruct%nested) atm(1)%phis = 0.
275  endif
276 
277 ! Read in the specified external dataset and do all the needed transformation
278  if ( atm(1)%flagstruct%ncep_ic ) then
279  nq = 1
280  call timing_on('NCEP_IC')
281  call get_ncep_ic( atm, fv_domain, nq )
282  call timing_off('NCEP_IC')
283 #ifdef FV_TRACERS
284  if (.not. cold_start) then
285  call fv_io_read_tracers( fv_domain, atm )
286  if(is_master()) write(*,*) 'All tracers except sphum replaced by FV IC'
287  endif
288 #endif
289  elseif ( atm(1)%flagstruct%nggps_ic ) then
290  call timing_on('NGGPS_IC')
291  call get_nggps_ic( atm, fv_domain )
292  call timing_off('NGGPS_IC')
293  elseif ( atm(1)%flagstruct%ecmwf_ic ) then
294  if( is_master() ) write(*,*) 'Calling get_ecmwf_ic'
295  call timing_on('ECMWF_IC')
296  call get_ecmwf_ic( atm, fv_domain )
297  call timing_off('ECMWF_IC')
298  else
299 ! The following is to read in legacy lat-lon FV core restart file
300 ! is Atm%q defined in all cases?
301  nq = size(atm(1)%q,4)
302  call get_fv_ic( atm, fv_domain, nq )
303  endif
304 
305  call prt_maxmin('PS', atm(1)%ps, is, ie, js, je, ng, 1, 0.01)
306  call prt_maxmin('T', atm(1)%pt, is, ie, js, je, ng, atm(1)%npz, 1.)
307  if (.not.atm(1)%flagstruct%hydrostatic) call prt_maxmin('W', atm(1)%w, is, ie, js, je, ng, atm(1)%npz, 1.)
308  call prt_maxmin('SPHUM', atm(1)%q(:,:,:,1), is, ie, js, je, ng, atm(1)%npz, 1.)
309  if ( atm(1)%flagstruct%nggps_ic ) then
310  call prt_maxmin('TS', atm(1)%ts, is, ie, js, je, 0, 1, 1.)
311  endif
312  if ( atm(1)%flagstruct%nggps_ic .or. atm(1)%flagstruct%ecmwf_ic ) then
313  sphum = get_tracer_index(model_atmos, 'sphum')
314  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
315  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
316  rainwat = get_tracer_index(model_atmos, 'rainwat')
317  snowwat = get_tracer_index(model_atmos, 'snowwat')
318  graupel = get_tracer_index(model_atmos, 'graupel')
319 #ifdef MULTI_GASES
320  spfo = get_tracer_index(model_atmos, 'spfo')
321  spfo2 = get_tracer_index(model_atmos, 'spfo2')
322  spfo3 = get_tracer_index(model_atmos, 'spfo3')
323 #else
324  o3mr = get_tracer_index(model_atmos, 'o3mr')
325 #endif
326 #ifdef CCPP
327  liq_aero = get_tracer_index(model_atmos, 'liq_aero')
328  ice_aero = get_tracer_index(model_atmos, 'ice_aero')
329 #endif
330 
331  if ( liq_wat > 0 ) &
332  call prt_maxmin('liq_wat', atm(1)%q(:,:,:,liq_wat), is, ie, js, je, ng, atm(1)%npz, 1.)
333  if ( ice_wat > 0 ) &
334  call prt_maxmin('ice_wat', atm(1)%q(:,:,:,ice_wat), is, ie, js, je, ng, atm(1)%npz, 1.)
335  if ( rainwat > 0 ) &
336  call prt_maxmin('rainwat', atm(1)%q(:,:,:,rainwat), is, ie, js, je, ng, atm(1)%npz, 1.)
337  if ( snowwat > 0 ) &
338  call prt_maxmin('snowwat', atm(1)%q(:,:,:,snowwat), is, ie, js, je, ng, atm(1)%npz, 1.)
339  if ( graupel > 0 ) &
340  call prt_maxmin('graupel', atm(1)%q(:,:,:,graupel), is, ie, js, je, ng, atm(1)%npz, 1.)
341 #ifdef MULTI_GASES
342  if ( spfo > 0 ) &
343  call prt_maxmin('SPFO', atm(1)%q(:,:,:,spfo), is, ie, js, je, ng, atm(1)%npz, 1.)
344  if ( spfo2 > 0 ) &
345  call prt_maxmin('SPFO2', atm(1)%q(:,:,:,spfo2), is, ie, js, je, ng, atm(1)%npz, 1.)
346  if ( spfo3 > 0 ) &
347  call prt_maxmin('SPFO3', atm(1)%q(:,:,:,spfo3), is, ie, js, je, ng, atm(1)%npz, 1.)
348 #else
349  if ( o3mr > 0 ) &
350  call prt_maxmin('O3MR', atm(1)%q(:,:,:,o3mr), is, ie, js, je, ng, atm(1)%npz, 1.)
351 #endif
352 #ifdef CCPP
353  if ( liq_aero > 0) &
354  call prt_maxmin('liq_aero',atm(1)%q(:,:,:,liq_aero),is, ie, js, je, ng, atm(1)%npz, 1.)
355  if ( ice_aero > 0) &
356  call prt_maxmin('ice_aero',atm(1)%q(:,:,:,ice_aero),is, ie, js, je, ng, atm(1)%npz, 1.)
357 #endif
358  endif
359 
360  call p_var(atm(1)%npz, is, ie, js, je, atm(1)%ak(1), ptop_min, &
361  atm(1)%delp, atm(1)%delz, atm(1)%pt, atm(1)%ps, &
362  atm(1)%pe, atm(1)%peln, atm(1)%pk, atm(1)%pkz, &
363  kappa, atm(1)%q, ng, atm(1)%ncnst, atm(1)%gridstruct%area_64, atm(1)%flagstruct%dry_mass, &
364  atm(1)%flagstruct%adjust_dry_mass, atm(1)%flagstruct%mountain, atm(1)%flagstruct%moist_phys, &
365  atm(1)%flagstruct%hydrostatic, atm(1)%flagstruct%nwat, atm(1)%domain, atm(1)%flagstruct%make_nh)
366 
367  end subroutine get_external_ic
368 
369 
370 !------------------------------------------------------------------
371  subroutine get_cubed_sphere_terrain( Atm, fv_domain )
372  type(fv_atmos_type), intent(inout), target :: Atm(:)
373  type(domain2d), intent(inout) :: fv_domain
374  integer :: ntileMe
375  integer, allocatable :: tile_id(:)
376  character(len=64) :: fname
377  character(len=7) :: gn
378  integer :: n
379  integer :: jbeg, jend
380  real ftop
381  real, allocatable :: g_dat2(:,:,:)
382  real, allocatable :: pt_coarse(:,:,:)
383  integer isc_p, iec_p, jsc_p, jec_p, isg, ieg, jsg,jeg
384 
385  integer :: is, ie, js, je
386  integer :: isd, ied, jsd, jed
387 
388  is = atm(1)%bd%is
389  ie = atm(1)%bd%ie
390  js = atm(1)%bd%js
391  je = atm(1)%bd%je
392  isd = atm(1)%bd%isd
393  ied = atm(1)%bd%ied
394  jsd = atm(1)%bd%jsd
395  jed = atm(1)%bd%jed
396 
397  if (atm(1)%grid_number > 1) then
398  !write(gn,'(A2, I1)') ".g", Atm(1)%grid_number
399  write(gn,'(A5, I2.2)') ".nest", atm(1)%grid_number
400  else
401  gn = ''
402  end if
403 
404  ntileme = size(atm(:)) ! This will have to be modified for mult tiles per PE
405  ! ASSUMED always one at this point
406 
407  allocate( tile_id(ntileme) )
408  tile_id = mpp_get_tile_id( fv_domain )
409  do n=1,ntileme
410 
411  call get_tile_string(fname, 'INPUT/fv_core.res'//trim(gn)//'.tile', tile_id(n), '.nc' )
412  if (mpp_pe() == mpp_root_pe()) print*, 'external_ic: looking for ', fname
413 
414 
415  if( file_exist(fname) ) then
416  call read_data(fname, 'phis', atm(n)%phis(is:ie,js:je), &
417  domain=fv_domain, tile_count=n)
418  else
419  call surfdrv( atm(n)%npx, atm(n)%npy, atm(n)%gridstruct%grid_64, atm(n)%gridstruct%agrid_64, &
420  atm(n)%gridstruct%area_64, atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, &
421  atm(n)%gridstruct%dxa, atm(n)%gridstruct%dya, &
422  atm(n)%gridstruct%dxc, atm(n)%gridstruct%dyc, atm(n)%gridstruct%sin_sg, &
423  atm(n)%phis, atm(n)%flagstruct%stretch_fac, &
424  atm(n)%neststruct%nested, atm(n)%neststruct%npx_global, atm(n)%domain, &
425  atm(n)%flagstruct%grid_number, atm(n)%bd, atm(n)%flagstruct%regional )
426  call mpp_error(note,'terrain datasets generated using USGS data')
427  endif
428 
429  end do
430 
431 ! Needed for reproducibility. DON'T REMOVE THIS!!
432  call mpp_update_domains( atm(1)%phis, atm(1)%domain )
433  ftop = g_sum(atm(1)%domain, atm(1)%phis(is:ie,js:je), is, ie, js, je, ng, atm(1)%gridstruct%area_64, 1)
434 
435  call prt_maxmin('ZS', atm(1)%phis, is, ie, js, je, ng, 1, 1./grav)
436  if(is_master()) write(*,*) 'mean terrain height (m)=', ftop/grav
437 
438  deallocate( tile_id )
439 
440  end subroutine get_cubed_sphere_terrain
441 
445  subroutine get_nggps_ic (Atm, fv_domain)
466 #ifdef __PGI
467  use gfs_restart, only : gfs_restart_type
468 
469  implicit none
470 #endif
471 
472  type(fv_atmos_type), intent(inout) :: Atm(:)
473  type(domain2d), intent(inout) :: fv_domain
474 ! local:
475  real, dimension(:), allocatable:: ak, bk
476  real, dimension(:,:), allocatable:: wk2, ps, oro_g
477  real, dimension(:,:,:), allocatable:: ud, vd, u_w, v_w, u_s, v_s, omga, temp
478  real, dimension(:,:,:), allocatable:: zh(:,:,:) ! 3D height at 65 edges
479  real, dimension(:,:,:,:), allocatable:: q
480  real, dimension(:,:), allocatable :: phis_coarse ! lmh
481  real rdg, wt, qt, m_fac
482  integer:: n, npx, npy, npz, itoa, nt, ntprog, ntdiag, ntracers, ntrac, iq
483  integer :: is, ie, js, je
484  integer :: isd, ied, jsd, jed
485  integer :: ios, ierr, unit, id_res
486  type(restart_file_type) :: oro_restart, sfc_restart, gfs_restart
487  character(len=6) :: gn, stile_name
488  character(len=64) :: tracer_name
489  character(len=64) :: fn_gfs_ctl = 'gfs_ctrl.nc'
490  character(len=64) :: fn_gfs_ics = 'gfs_data.nc'
491  character(len=64) :: fn_sfc_ics = 'sfc_data.nc'
492  character(len=64) :: fn_oro_ics = 'oro_data.nc'
493  ! DH* character(len=64) :: fn_aero_ics = 'aero_data.nc' *DH
494  logical :: remap
495  logical :: filtered_terrain = .true.
496  logical :: gfs_dwinds = .true.
497  integer :: levp = 64
498  logical :: checker_tr = .false.
499  integer :: nt_checker = 0
500  real(kind=R_GRID), dimension(2):: p1, p2, p3
501  real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
502  integer:: i,j,k,nts, ks
503  integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, ntclamt
504  namelist /external_ic_nml/ filtered_terrain, levp, gfs_dwinds, &
505  checker_tr, nt_checker
506 #ifdef GFSL64
507  real, dimension(65):: ak_sj, bk_sj
508  data ak_sj/20.00000, 68.00000, 137.79000, &
509  221.95800, 318.26600, 428.43400, &
510  554.42400, 698.45700, 863.05803, &
511  1051.07995, 1265.75194, 1510.71101, &
512  1790.05098, 2108.36604, 2470.78817, &
513  2883.03811, 3351.46002, 3883.05187, &
514  4485.49315, 5167.14603, 5937.04991, &
515  6804.87379, 7780.84698, 8875.64338, &
516  9921.40745, 10760.99844, 11417.88354, &
517  11911.61193, 12258.61668, 12472.89642, &
518  12566.58298, 12550.43517, 12434.26075, &
519  12227.27484, 11938.39468, 11576.46910, &
520  11150.43640, 10669.41063, 10142.69482, &
521  9579.72458, 8989.94947, 8382.67090, &
522  7766.85063, 7150.91171, 6542.55077, &
523  5948.57894, 5374.81094, 4825.99383, &
524  4305.79754, 3816.84622, 3360.78848, &
525  2938.39801, 2549.69756, 2194.08449, &
526  1870.45732, 1577.34218, 1313.00028, &
527  1075.52114, 862.90778, 673.13815, &
528  504.22118, 354.22752, 221.32110, &
529  103.78014, 0./
530  data bk_sj/0.00000, 0.00000, 0.00000, &
531  0.00000, 0.00000, 0.00000, &
532  0.00000, 0.00000, 0.00000, &
533  0.00000, 0.00000, 0.00000, &
534  0.00000, 0.00000, 0.00000, &
535  0.00000, 0.00000, 0.00000, &
536  0.00000, 0.00000, 0.00000, &
537  0.00000, 0.00000, 0.00000, &
538  0.00179, 0.00705, 0.01564, &
539  0.02749, 0.04251, 0.06064, &
540  0.08182, 0.10595, 0.13294, &
541  0.16266, 0.19492, 0.22950, &
542  0.26615, 0.30455, 0.34435, &
543  0.38516, 0.42656, 0.46815, &
544  0.50949, 0.55020, 0.58989, &
545  0.62825, 0.66498, 0.69987, &
546  0.73275, 0.76351, 0.79208, &
547  0.81845, 0.84264, 0.86472, &
548  0.88478, 0.90290, 0.91923, &
549  0.93388, 0.94697, 0.95865, &
550  0.96904, 0.97826, 0.98642, &
551  0.99363, 1./
552 #else
553 ! The following L63 setting is the same as NCEP GFS's L64 except the top layer
554  real, dimension(64):: ak_sj, bk_sj
555  data ak_sj/64.247, 137.790, 221.958, &
556  318.266, 428.434, 554.424, &
557  698.457, 863.05803, 1051.07995, &
558  1265.75194, 1510.71101, 1790.05098, &
559  2108.36604, 2470.78817, 2883.03811, &
560  3351.46002, 3883.05187, 4485.49315, &
561  5167.14603, 5937.04991, 6804.87379, &
562  7780.84698, 8875.64338, 10100.20534, &
563  11264.35673, 12190.64366, 12905.42546, &
564  13430.87867, 13785.88765, 13986.77987, &
565  14047.96335, 13982.46770, 13802.40331, &
566  13519.33841, 13144.59486, 12689.45608, &
567  12165.28766, 11583.57006, 10955.84778, &
568  10293.60402, 9608.08306, 8910.07678, &
569  8209.70131, 7516.18560, 6837.69250, &
570  6181.19473, 5552.39653, 4955.72632, &
571  4394.37629, 3870.38682, 3384.76586, &
572  2937.63489, 2528.37666, 2155.78385, &
573  1818.20722, 1513.68173, 1240.03585, &
574  994.99144, 776.23591, 581.48797, &
575  408.53400, 255.26520, 119.70243, 0. /
576 
577  data bk_sj/0.00000, 0.00000, 0.00000, &
578  0.00000, 0.00000, 0.00000, &
579  0.00000, 0.00000, 0.00000, &
580  0.00000, 0.00000, 0.00000, &
581  0.00000, 0.00000, 0.00000, &
582  0.00000, 0.00000, 0.00000, &
583  0.00000, 0.00000, 0.00000, &
584  0.00000, 0.00000, 0.00000, &
585  0.00201, 0.00792, 0.01755, &
586  0.03079, 0.04751, 0.06761, &
587  0.09097, 0.11746, 0.14690, &
588  0.17911, 0.21382, 0.25076, &
589  0.28960, 0.32994, 0.37140, &
590  0.41353, 0.45589, 0.49806, &
591  0.53961, 0.58015, 0.61935, &
592  0.65692, 0.69261, 0.72625, &
593  0.75773, 0.78698, 0.81398, &
594  0.83876, 0.86138, 0.88192, &
595  0.90050, 0.91722, 0.93223, &
596  0.94565, 0.95762, 0.96827, &
597  0.97771, 0.98608, 0.99347, 1./
598 #endif
599 
600 #ifdef TEMP_GFSPLV
601  real, dimension(64):: ak_sj, bk_sj
602  data ak_sj/64.247, 137.79, 221.958, &
603  318.266, 428.434, 554.424, &
604  698.457, 863.058, 1051.08, &
605  1265.752, 1510.711, 1790.051, &
606  2108.366, 2470.788, 2883.038, &
607  3351.46, 3883.052, 4485.493, &
608  5167.146, 5937.05, 6804.874, &
609  7777.15, 8832.537, 9936.614, &
610  11054.85, 12152.94, 13197.07, &
611  14154.32, 14993.07, 15683.49, &
612  16197.97, 16511.74, 16611.6, &
613  16503.14, 16197.32, 15708.89, &
614  15056.34, 14261.43, 13348.67, &
615  12344.49, 11276.35, 10171.71, &
616  9057.051, 7956.908, 6893.117, &
617  5884.206, 4945.029, 4086.614, &
618  3316.217, 2637.553, 2051.15, &
619  1554.789, 1143.988, 812.489, &
620  552.72, 356.223, 214.015, &
621  116.899, 55.712, 21.516, &
622  5.741, 0.575, 0., 0. /
623 
624  data bk_sj/0.00000, 0.00000, 0.00000, &
625  0.00000, 0.00000, 0.00000, &
626  0.00000, 0.00000, 0.00000, &
627  0.00000, 0.00000, 0.00000, &
628  0.00000, 0.00000, 0.00000, &
629  0.00000, 0.00000, 0.00000, &
630  0.00000, 0.00000, 0.00000, &
631  0.00003697, 0.00043106, 0.00163591, &
632  0.00410671, 0.00829402, 0.01463712, &
633  0.02355588, 0.03544162, 0.05064684, &
634  0.06947458, 0.09216691, 0.1188122, &
635  0.1492688, 0.1832962, 0.2205702, &
636  0.2606854, 0.3031641, 0.3474685, &
637  0.3930182, 0.4392108, 0.4854433, &
638  0.5311348, 0.5757467, 0.6187996, &
639  0.659887, 0.6986829, 0.7349452, &
640  0.7685147, 0.7993097, 0.8273188, &
641  0.8525907, 0.8752236, 0.895355, &
642  0.913151, 0.9287973, 0.9424911, &
643  0.9544341, 0.9648276, 0.9738676, &
644  0.9817423, 0.9886266, 0.9946712, 1./
645 #endif
646 
647  call mpp_error(note,'Using external_IC::get_nggps_ic which is valid only for data which has been &
648  &horizontally interpolated to the current cubed-sphere grid')
649 #ifdef INTERNAL_FILE_NML
650  read (input_nml_file,external_ic_nml,iostat=ios)
651  ierr = check_nml_error(ios,'external_ic_nml')
652 #else
653  unit=open_namelist_file()
654  read (unit,external_ic_nml,iostat=ios)
655  ierr = check_nml_error(ios,'external_ic_nml')
656  call close_file(unit)
657 #endif
658 
659  unit = stdlog()
660  call write_version_number ( 'EXTERNAL_IC_mod::get_nggps_ic', version )
661  write(unit, nml=external_ic_nml)
662 
663  remap = .true.
664  if (atm(1)%flagstruct%external_eta) then
665  if (filtered_terrain) then
666  call mpp_error(note,'External_IC::get_nggps_ic - use externally-generated, filtered terrain &
667  &and NCEP pressure levels (no vertical remapping)')
668  else if (.not. filtered_terrain) then
669  call mpp_error(note,'External_IC::get_nggps_ic - use externally-generated, raw terrain &
670  &and NCEP pressure levels (no vertical remapping)')
671  endif
672  else ! (.not.external_eta)
673  if (filtered_terrain) then
674  call mpp_error(note,'External_IC::get_nggps_ic - use externally-generated, filtered terrain &
675  &and FV3 pressure levels (vertical remapping)')
676  else if (.not. filtered_terrain) then
677  call mpp_error(note,'External_IC::get_nggps_ic - use externally-generated, raw terrain &
678  &and FV3 pressure levels (vertical remapping)')
679  endif
680  endif
681 
682  is = atm(1)%bd%is
683  ie = atm(1)%bd%ie
684  js = atm(1)%bd%js
685  je = atm(1)%bd%je
686  isd = atm(1)%bd%isd
687  ied = atm(1)%bd%ied
688  jsd = atm(1)%bd%jsd
689  jed = atm(1)%bd%jed
690  npz = atm(1)%npz
691  write(*,22001)is,ie,js,je,isd,ied,jsd,jed
692 22001 format(' enter get_nggps_ic is=',i4,' ie=',i4,' js=',i4,' je=',i4,' isd=',i4,' ied=',i4,' jsd=',i4,' jed=',i4)
693  call get_number_tracers(model_atmos, num_tracers=ntracers, num_prog=ntprog)
694  ntdiag = ntracers-ntprog
695 
696 !--- test for existence of the GFS control file
697  if (.not. file_exist('INPUT/'//trim(fn_gfs_ctl), no_domain=.true.)) then
698  call mpp_error(fatal,'==> Error in External_ic::get_nggps_ic: file '//trim(fn_gfs_ctl)//' for NGGPS IC does not exist')
699  endif
700  call mpp_error(note,'==> External_ic::get_nggps_ic: using control file '//trim(fn_gfs_ctl)//' for NGGPS IC')
701 
702 !--- read in the number of tracers in the NCEP NGGPS ICs
703  call read_data ('INPUT/'//trim(fn_gfs_ctl), 'ntrac', ntrac, no_domain=.true.)
704  if (ntrac > ntracers) call mpp_error(fatal,'==> External_ic::get_nggps_ic: more NGGPS tracers &
705  &than defined in field_table '//trim(fn_gfs_ctl)//' for NGGPS IC')
706 
707 !--- read in ak and bk from the gfs control file using fms_io read_data ---
708  allocate (wk2(levp+1,2))
709  allocate (ak(levp+1))
710  allocate (bk(levp+1))
711  call read_data('INPUT/'//trim(fn_gfs_ctl),'vcoord',wk2, no_domain=.true.)
712  ak(1:levp+1) = wk2(1:levp+1,1)
713  bk(1:levp+1) = wk2(1:levp+1,2)
714  deallocate (wk2)
715 
716  if (.not. file_exist('INPUT/'//trim(fn_oro_ics), domain=atm(1)%domain)) then
717  call mpp_error(fatal,'==> Error in External_ic::get_nggps_ic: tiled file '//trim(fn_oro_ics)//' for NGGPS IC does not exist')
718  endif
719  call mpp_error(note,'==> External_ic::get_nggps_ic: using tiled data file '//trim(fn_oro_ics)//' for NGGPS IC')
720 
721  if (.not. file_exist('INPUT/'//trim(fn_sfc_ics), domain=atm(1)%domain)) then
722  call mpp_error(fatal,'==> Error in External_ic::get_nggps_ic: tiled file '//trim(fn_sfc_ics)//' for NGGPS IC does not exist')
723  endif
724  call mpp_error(note,'==> External_ic::get_nggps_ic: using tiled data file '//trim(fn_sfc_ics)//' for NGGPS IC')
725 
726  if (.not. file_exist('INPUT/'//trim(fn_gfs_ics), domain=atm(1)%domain)) then
727  call mpp_error(fatal,'==> Error in External_ic::get_nggps_ic: tiled file '//trim(fn_gfs_ics)//' for NGGPS IC does not exist')
728  endif
729  call mpp_error(note,'==> External_ic::get_nggps_ic: using tiled data file '//trim(fn_gfs_ics)//' for NGGPS IC')
730 !
731  call get_data_source(source,atm(1)%flagstruct%regional)
732 !
733  allocate (zh(is:ie,js:je,levp+1)) ! SJL
734  allocate (ps(is:ie,js:je))
735  allocate (omga(is:ie,js:je,levp))
736  allocate (q(is:ie,js:je,levp,ntracers))
737  allocate ( u_w(is:ie+1, js:je, 1:levp) )
738  allocate ( v_w(is:ie+1, js:je, 1:levp) )
739  allocate ( u_s(is:ie, js:je+1, 1:levp) )
740  allocate ( v_s(is:ie, js:je+1, 1:levp) )
741  allocate (temp(is:ie,js:je,levp))
742 
743  do n = 1,size(atm(:))
744 
745  !!! If a nested grid, save the filled coarse-grid topography for blending
746  if (atm(n)%neststruct%nested) then
747  allocate(phis_coarse(isd:ied,jsd:jed))
748  do j=jsd,jed
749  do i=isd,ied
750  phis_coarse(i,j) = atm(n)%phis(i,j)
751  enddo
752  enddo
753  endif
754 
755 !--- read in surface temperature (k) and land-frac
756  ! surface skin temperature
757  id_res = register_restart_field(sfc_restart, fn_sfc_ics, 'tsea', atm(n)%ts, domain=atm(n)%domain)
758 
759  ! terrain surface height -- (needs to be transformed into phis = zs*grav)
760  if (filtered_terrain) then
761  id_res = register_restart_field(oro_restart, fn_oro_ics, 'orog_filt', atm(n)%phis, domain=atm(n)%domain)
762  elseif (.not. filtered_terrain) then
763  id_res = register_restart_field(oro_restart, fn_oro_ics, 'orog_raw', atm(n)%phis, domain=atm(n)%domain)
764  endif
765 
766  if ( atm(n)%flagstruct%full_zs_filter) then
767  allocate (oro_g(isd:ied,jsd:jed))
768  oro_g = 0.
769  ! land-frac
770  id_res = register_restart_field(oro_restart, fn_oro_ics, 'land_frac', oro_g, domain=atm(n)%domain)
771  call mpp_update_domains(oro_g, atm(n)%domain)
772  if (atm(n)%neststruct%nested) then
773  call extrapolation_bc(oro_g, 0, 0, atm(n)%npx, atm(n)%npy, atm(n)%bd, .true.)
774  endif
775  endif
776 
777  if ( atm(n)%flagstruct%fv_land ) then
778  ! stddev
779  id_res = register_restart_field(oro_restart, fn_oro_ics, 'stddev', atm(n)%sgh, domain=atm(n)%domain)
780  ! land-frac
781  id_res = register_restart_field(oro_restart, fn_oro_ics, 'land_frac', atm(n)%oro, domain=atm(n)%domain)
782  endif
783 
784  ! surface pressure (Pa)
785  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'ps', ps, domain=atm(n)%domain)
786 
787  ! D-grid west face tangential wind component (m/s)
788  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'u_w', u_w, domain=atm(n)%domain,position=east)
789  ! D-grid west face normal wind component (m/s)
790  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'v_w', v_w, domain=atm(n)%domain,position=east)
791  ! D-grid south face tangential wind component (m/s)
792  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'u_s', u_s, domain=atm(n)%domain,position=north)
793  ! D-grid south face normal wind component (m/s)
794  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'v_s', v_s, domain=atm(n)%domain,position=north)
795 
796  ! vertical velocity 'omega' (Pa/s)
797  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'w', omga, domain=atm(n)%domain)
798  ! GFS grid height at edges (including surface height)
799  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'ZH', zh, domain=atm(n)%domain)
800  ! real temperature (K)
801  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 't', temp, mandatory=.false., &
802  domain=atm(n)%domain)
803  ! prognostic tracers
804  do nt = 1, ntracers
805  call get_tracer_names(model_atmos, nt, tracer_name)
806  ! DH* if aerosols are in separate file, need to test for indices liq_aero and ice_aero and change fn_gfs_ics to fn_aero_ics *DH
807  id_res = register_restart_field(gfs_restart, fn_gfs_ics, trim(tracer_name), q(:,:,:,nt), &
808  mandatory=.false.,domain=atm(n)%domain)
809  enddo
810 
811  ! initialize all tracers to default values prior to being input
812  do nt = 1, ntprog
813  call get_tracer_names(model_atmos, nt, tracer_name)
814  ! set all tracers to an initial profile value
815  call set_tracer_profile (model_atmos, nt, atm(n)%q(:,:,:,nt) )
816  enddo
817  do nt = ntprog+1, ntracers
818  call get_tracer_names(model_atmos, nt, tracer_name)
819  ! set all tracers to an initial profile value
820  call set_tracer_profile (model_atmos, nt, atm(n)%qdiag(:,:,:,nt) )
821  enddo
822 
823  ! read in the restart
824  call restore_state (oro_restart)
825  call restore_state (sfc_restart)
826  call restore_state (gfs_restart)
827 
828  ! free the restart type to be re-used by the nest
829  call free_restart_type(oro_restart)
830  call free_restart_type(sfc_restart)
831  call free_restart_type(gfs_restart)
832 
833  ! multiply NCEP ICs terrain 'phis' by gravity to be true geopotential
834  atm(n)%phis = atm(n)%phis*grav
835 
836  ! set the pressure levels and ptop to be used
837  if (atm(1)%flagstruct%external_eta) then
838  itoa = levp - npz + 1
839  atm(n)%ptop = ak(itoa)
840  atm(n)%ak(1:npz+1) = ak(itoa:levp+1)
841  atm(n)%bk(1:npz+1) = bk(itoa:levp+1)
842  call set_external_eta (atm(n)%ak, atm(n)%bk, atm(n)%ptop, atm(n)%ks)
843  else
844  if ( npz <= 64 ) then
845  atm(n)%ak(:) = ak_sj(:)
846  atm(n)%bk(:) = bk_sj(:)
847  atm(n)%ptop = atm(n)%ak(1)
848  else
849  call set_eta(npz, ks, atm(n)%ptop, atm(n)%ak, atm(n)%bk)
850  endif
851  endif
852  ! call vertical remapping algorithms
853  if(is_master()) write(*,*) 'GFS ak =', ak,' FV3 ak=',atm(n)%ak
854  ak(1) = max(1.e-9, ak(1))
855 
856 !*** For regional runs read in each of the BC variables from the NetCDF boundary file
857 !*** and remap in the vertical from the input levels to the model integration levels.
858 !*** Here in the initialization we begn by allocating the regional domain's boundary
859 !*** objects. Then we need to read the first two regional BC files so the integration
860 !*** can begin interpolating between those two times as the forecast proceeds.
861 
862  if (n==1.and.atm(1)%flagstruct%regional) then
863 
864  call start_regional_cold_start(atm(1), ak, bk, levp, &
865  is, ie, js, je, &
866  isd, ied, jsd, jed )
867  endif
868 
869 !
870 !*** Remap the variables in the compute domain.
871 !
872  call remap_scalar_nggps(atm(n), levp, npz, ntracers, ak, bk, ps, temp, q, omga, zh)
873 
874  allocate ( ud(is:ie, js:je+1, 1:levp) )
875  allocate ( vd(is:ie+1,js:je, 1:levp) )
876 
877 !$OMP parallel do default(none) shared(is,ie,js,je,levp,Atm,ud,vd,u_s,v_s,u_w,v_w) &
878 !$OMP private(p1,p2,p3,e1,e2,ex,ey)
879  do k=1,levp
880  do j=js,je+1
881  do i=is,ie
882  p1(:) = atm(1)%gridstruct%grid(i, j,1:2)
883  p2(:) = atm(1)%gridstruct%grid(i+1,j,1:2)
884  call mid_pt_sphere(p1, p2, p3)
885  call get_unit_vect2(p1, p2, e1)
886  call get_latlon_vector(p3, ex, ey)
887  ud(i,j,k) = u_s(i,j,k)*inner_prod(e1,ex) + v_s(i,j,k)*inner_prod(e1,ey)
888  enddo
889  enddo
890  do j=js,je
891  do i=is,ie+1
892  p1(:) = atm(1)%gridstruct%grid(i,j ,1:2)
893  p2(:) = atm(1)%gridstruct%grid(i,j+1,1:2)
894  call mid_pt_sphere(p1, p2, p3)
895  call get_unit_vect2(p1, p2, e2)
896  call get_latlon_vector(p3, ex, ey)
897  vd(i,j,k) = u_w(i,j,k)*inner_prod(e2,ex) + v_w(i,j,k)*inner_prod(e2,ey)
898  enddo
899  enddo
900  enddo
901  deallocate ( u_w )
902  deallocate ( v_w )
903  deallocate ( u_s )
904  deallocate ( v_s )
905 
906  call remap_dwinds(levp, npz, ak, bk, ps, ud, vd, atm(n))
907 
908  deallocate ( ud )
909  deallocate ( vd )
910 
911  if (atm(n)%neststruct%nested) then
912  if (is_master()) write(*,*) 'Blending nested and coarse grid topography'
913  npx = atm(n)%npx
914  npy = atm(n)%npy
915  do j=jsd,jed
916  do i=isd,ied
917  wt = max(0.,min(1.,real(5 - min(i,j,npx-i,npy-j,5))/5. ))
918  atm(n)%phis(i,j) = (1.-wt)*atm(n)%phis(i,j) + wt*phis_coarse(i,j)
919  enddo
920  enddo
921  endif
922 
923 
924  !!! Perform terrain smoothing, if desired
925  if ( atm(n)%flagstruct%full_zs_filter ) then
926 
927  call mpp_update_domains(atm(n)%phis, atm(n)%domain)
928 
929  call fv3_zs_filter( atm(n)%bd, isd, ied, jsd, jed, npx, npy, atm(n)%neststruct%npx_global, &
930  atm(n)%flagstruct%stretch_fac, atm(n)%neststruct%nested, atm(n)%domain, &
931  atm(n)%gridstruct%area_64, atm(n)%gridstruct%dxa, atm(n)%gridstruct%dya, &
932  atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, atm(n)%gridstruct%dxc, &
933  atm(n)%gridstruct%dyc, atm(n)%gridstruct%grid_64, atm(n)%gridstruct%agrid_64, &
934  atm(n)%gridstruct%sin_sg, atm(n)%phis, oro_g, atm(n)%flagstruct%regional)
935  deallocate(oro_g)
936  endif
937 
938 
939  if ( atm(n)%flagstruct%n_zs_filter > 0 ) then
940 
941  if ( atm(n)%flagstruct%nord_zs_filter == 2 ) then
942  call del2_cubed_sphere(atm(n)%npx, atm(n)%npy, atm(n)%phis, &
943  atm(n)%gridstruct%area_64, atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, &
944  atm(n)%gridstruct%dxc, atm(n)%gridstruct%dyc, atm(n)%gridstruct%sin_sg, &
945  atm(n)%flagstruct%n_zs_filter, cnst_0p20*atm(n)%gridstruct%da_min, &
946  .false., oro_g, atm(n)%neststruct%nested, atm(n)%domain, atm(n)%bd, atm(n)%flagstruct%regional)
947  if ( is_master() ) write(*,*) 'Warning !!! del-2 terrain filter has been applied ', &
948  atm(n)%flagstruct%n_zs_filter, ' times'
949  else if( atm(n)%flagstruct%nord_zs_filter == 4 ) then
950  call del4_cubed_sphere(atm(n)%npx, atm(n)%npy, atm(n)%phis, atm(n)%gridstruct%area_64, &
951  atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, &
952  atm(n)%gridstruct%dxc, atm(n)%gridstruct%dyc, atm(n)%gridstruct%sin_sg, &
953  atm(n)%flagstruct%n_zs_filter, .false., oro_g, atm(n)%neststruct%nested, &
954  atm(n)%domain, atm(n)%bd, atm(n)%flagstruct%regional)
955  if ( is_master() ) write(*,*) 'Warning !!! del-4 terrain filter has been applied ', &
956  atm(n)%flagstruct%n_zs_filter, ' times'
957  endif
958 
959  endif
960 
961  if ( atm(n)%neststruct%nested .and. ( atm(n)%flagstruct%n_zs_filter > 0 .or. atm(n)%flagstruct%full_zs_filter ) ) then
962  npx = atm(n)%npx
963  npy = atm(n)%npy
964  do j=jsd,jed
965  do i=isd,ied
966  wt = max(0.,min(1.,real(5 - min(i,j,npx-i,npy-j,5))/5. ))
967  atm(n)%phis(i,j) = (1.-wt)*atm(n)%phis(i,j) + wt*phis_coarse(i,j)
968  enddo
969  enddo
970  deallocate(phis_coarse)
971  endif
972 
973  call mpp_update_domains( atm(n)%phis, atm(n)%domain, complete=.true. )
974  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
975  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
976  rainwat = get_tracer_index(model_atmos, 'rainwat')
977  snowwat = get_tracer_index(model_atmos, 'snowwat')
978  graupel = get_tracer_index(model_atmos, 'graupel')
979  ntclamt = get_tracer_index(model_atmos, 'cld_amt')
980  if (trim(source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
981  do k=1,npz
982  do j=js,je
983  do i=is,ie
984  wt = atm(n)%delp(i,j,k)
985  if ( atm(n)%flagstruct%nwat == 6 ) then
986  qt = wt*(1. + atm(n)%q(i,j,k,liq_wat) + &
987  atm(n)%q(i,j,k,ice_wat) + &
988  atm(n)%q(i,j,k,rainwat) + &
989  atm(n)%q(i,j,k,snowwat) + &
990  atm(n)%q(i,j,k,graupel))
991  else ! all other values of nwat
992  qt = wt*(1. + sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
993  endif
994  atm(n)%delp(i,j,k) = qt
995  if (ntclamt > 0) atm(n)%q(i,j,k,ntclamt) = 0.0 ! Moorthi
996  enddo
997  enddo
998  enddo
999  else
1000 !--- Add cloud condensate from GFS to total MASS
1001 ! 20160928: Adjust the mixing ratios consistently...
1002  do k=1,npz
1003  do j=js,je
1004  do i=is,ie
1005  wt = atm(n)%delp(i,j,k)
1006  if ( atm(n)%flagstruct%nwat == 6 ) then
1007  qt = wt*(1. + atm(n)%q(i,j,k,liq_wat) + &
1008  atm(n)%q(i,j,k,ice_wat) + &
1009  atm(n)%q(i,j,k,rainwat) + &
1010  atm(n)%q(i,j,k,snowwat) + &
1011  atm(n)%q(i,j,k,graupel))
1012  else ! all other values of nwat
1013  qt = wt*(1. + sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
1014  endif
1015  m_fac = wt / qt
1016  do iq=1,ntracers
1017  atm(n)%q(i,j,k,iq) = m_fac * atm(n)%q(i,j,k,iq)
1018  enddo
1019  atm(n)%delp(i,j,k) = qt
1020  if (ntclamt > 0) atm(n)%q(i,j,k,ntclamt) = 0.0 ! Moorthi
1021  enddo
1022  enddo
1023  enddo
1024  endif !end trim(source) test
1025 
1026 !--- reset the tracers beyond condensate to a checkerboard pattern
1027  if (checker_tr) then
1028  nts = ntracers - nt_checker+1
1029  call checker_tracers(is,ie, js,je, isd,ied, jsd,jed, nt_checker, &
1030  npz, atm(n)%q(:,:,:,nts:ntracers), &
1031  atm(n)%gridstruct%agrid_64(is:ie,js:je,1), &
1032  atm(n)%gridstruct%agrid_64(is:ie,js:je,2), 9., 9.)
1033  endif
1034  enddo ! n-loop
1035 
1036  atm(1)%flagstruct%make_nh = .false.
1037 
1038  deallocate (ak)
1039  deallocate (bk)
1040  deallocate (ps)
1041  deallocate (q )
1042  deallocate (temp)
1043  deallocate (omga)
1044 
1045  end subroutine get_nggps_ic
1046 !------------------------------------------------------------------
1047 !------------------------------------------------------------------
1049  subroutine get_ncep_ic( Atm, fv_domain, nq )
1050  type(fv_atmos_type), intent(inout) :: Atm(:)
1051  type(domain2d), intent(inout) :: fv_domain
1052  integer, intent(in):: nq
1053 ! local:
1054 #ifdef HIWPP_ETA
1055  real :: ak_HIWPP(65), bk_HIWPP(65)
1056  data ak_hiwpp/ &
1057  0, 0.00064247, 0.0013779, 0.00221958, 0.00318266, 0.00428434, &
1058  0.00554424, 0.00698457, 0.00863058, 0.0105108, 0.01265752, 0.01510711, &
1059  0.01790051, 0.02108366, 0.02470788, 0.02883038, 0.0335146, 0.03883052, &
1060  0.04485493, 0.05167146, 0.0593705, 0.06804874, 0.0777715, 0.08832537, &
1061  0.09936614, 0.1105485, 0.1215294, 0.1319707, 0.1415432, 0.1499307, &
1062  0.1568349, 0.1619797, 0.1651174, 0.166116, 0.1650314, 0.1619731, &
1063  0.1570889, 0.1505634, 0.1426143, 0.1334867, 0.1234449, 0.1127635, &
1064  0.1017171, 0.09057051, 0.07956908, 0.06893117, 0.05884206, 0.04945029, &
1065  0.04086614, 0.03316217, 0.02637553, 0.0205115, 0.01554789, 0.01143988, &
1066  0.00812489, 0.0055272, 0.00356223, 0.00214015, 0.00116899, 0.00055712, &
1067  0.00021516, 5.741e-05, 5.75e-06, 0, 0 /
1068 
1069  data bk_hiwpp/ &
1070  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1071  3.697e-05, 0.00043106, 0.00163591, 0.00410671, 0.00829402, 0.01463712, &
1072  0.02355588, 0.03544162, 0.05064684, 0.06947458, 0.09216691, 0.1188122, &
1073  0.1492688, 0.1832962, 0.2205702, 0.2606854, 0.3031641, 0.3474685, &
1074  0.3930182, 0.4392108, 0.4854433, 0.5311348, 0.5757467, 0.6187996, &
1075  0.659887, 0.6986829, 0.7349452, 0.7685147, 0.7993097, 0.8273188, &
1076  0.8525907, 0.8752236, 0.895355, 0.913151, 0.9287973, 0.9424911, &
1077  0.9544341, 0.9648276, 0.9738676, 0.9817423, 0.9886266, 0.9946712, 1 /
1078 #endif
1079  character(len=128) :: fname
1080  real(kind=4), allocatable:: wk1(:), wk2(:,:), wk3(:,:,:)
1081  real, allocatable:: tp(:,:,:), qp(:,:,:)
1082  real, allocatable:: ua(:,:,:), va(:,:,:)
1083  real, allocatable:: lat(:), lon(:), ak0(:), bk0(:)
1084  real:: s2c(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je,4)
1085  integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je):: id1, id2, jdc
1086  real psc(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je)
1087  real gzc(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je)
1088  real tmean
1089  integer:: i, j, k, im, jm, km, npz, npt
1090  integer:: i1, i2, j1, ncid
1091  integer:: jbeg, jend
1092  integer tsize(3)
1093  logical:: read_ts = .true.
1094  logical:: land_ts = .false.
1095  logical:: found
1096  integer :: is, ie, js, je
1097  integer :: isd, ied, jsd, jed
1098 
1099  is = atm(1)%bd%is
1100  ie = atm(1)%bd%ie
1101  js = atm(1)%bd%js
1102  je = atm(1)%bd%je
1103  isd = atm(1)%bd%isd
1104  ied = atm(1)%bd%ied
1105  jsd = atm(1)%bd%jsd
1106  jed = atm(1)%bd%jed
1107 
1108  deg2rad = pi/180.
1109 
1110  npz = atm(1)%npz
1111 
1112 ! Zero out all initial tracer fields:
1113 ! SJL: 20110716
1114 ! Atm(1)%q = 0.
1115 
1116  fname = atm(1)%flagstruct%res_latlon_dynamics
1117 
1118  if( file_exist(fname) ) then
1119  call open_ncfile( fname, ncid ) ! open the file
1120  call get_ncdim1( ncid, 'lon', tsize(1) )
1121  call get_ncdim1( ncid, 'lat', tsize(2) )
1122  call get_ncdim1( ncid, 'lev', tsize(3) )
1123 
1124  im = tsize(1); jm = tsize(2); km = tsize(3)
1125 
1126  if(is_master()) write(*,*) fname
1127  if(is_master()) write(*,*) ' NCEP IC dimensions:', tsize
1128 
1129  allocate ( lon(im) )
1130  allocate ( lat(jm) )
1131 
1132  call _get_var1(ncid, 'lon', im, lon )
1133  call _get_var1(ncid, 'lat', jm, lat )
1134 
1135 ! Convert to radian
1136  do i=1,im
1137  lon(i) = lon(i) * deg2rad ! lon(1) = 0.
1138  enddo
1139  do j=1,jm
1140  lat(j) = lat(j) * deg2rad
1141  enddo
1142 
1143  allocate ( ak0(km+1) )
1144  allocate ( bk0(km+1) )
1145 
1146 #ifdef HIWPP_ETA
1147 ! The HIWPP data from Jeff does not contain (ak,bk)
1148  do k=1, km+1
1149  ak0(k) = ak_hiwpp(k)
1150  bk0(k) = bk_hiwpp(k)
1151  enddo
1152 #else
1153  call _get_var1(ncid, 'hyai', km+1, ak0, found )
1154  if ( .not. found ) ak0(:) = 0.
1155 
1156  call _get_var1(ncid, 'hybi', km+1, bk0 )
1157 #endif
1158  if( is_master() ) then
1159  do k=1,km+1
1160  write(*,*) k, ak0(k), bk0(k)
1161  enddo
1162  endif
1163 
1164 ! Note: definition of NCEP hybrid is p(k) = a(k)*1.E5 + b(k)*ps
1165  ak0(:) = ak0(:) * 1.e5
1166 
1167 ! Limiter to prevent NAN at top during remapping
1168  if ( bk0(1) < 1.e-9 ) ak0(1) = max(1.e-9, ak0(1))
1169 
1170  else
1171  call mpp_error(fatal,'==> Error in get_external_ic: Expected file '//trim(fname)//' for NCEP IC does not exist')
1172  endif
1173 
1174 ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
1175  call remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
1176  im, jm, lon, lat, id1, id2, jdc, s2c , atm(1)%gridstruct%agrid)
1177 
1178 ! Find bounding latitudes:
1179  jbeg = jm-1; jend = 2
1180  do j=js,je
1181  do i=is,ie
1182  j1 = jdc(i,j)
1183  jbeg = min(jbeg, j1)
1184  jend = max(jend, j1+1)
1185  enddo
1186  enddo
1187 
1188 ! remap surface pressure and height:
1189 
1190  allocate ( wk2(im,jbeg:jend) )
1191  call get_var3_r4( ncid, 'PS', 1,im, jbeg,jend, 1,1, wk2 )
1192 
1193  do j=js,je
1194  do i=is,ie
1195  i1 = id1(i,j)
1196  i2 = id2(i,j)
1197  j1 = jdc(i,j)
1198  psc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1199  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1200  enddo
1201  enddo
1202 
1203  call get_var3_r4( ncid, 'PHIS', 1,im, jbeg,jend, 1,1, wk2 )
1204  do j=js,je
1205  do i=is,ie
1206  i1 = id1(i,j)
1207  i2 = id2(i,j)
1208  j1 = jdc(i,j)
1209  gzc(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1210  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1211  enddo
1212  enddo
1213 
1214  deallocate ( wk2 )
1215  allocate ( wk2(im,jm) )
1216 
1217  if ( read_ts ) then ! read skin temperature; could be used for SST
1218 
1219  call get_var2_real( ncid, 'TS', im, jm, wk2 )
1220 
1221  if ( .not. land_ts ) then
1222  allocate ( wk1(im) )
1223 
1224  do j=1,jm
1225 ! Read NCEP ORO (1; land; 0: ocean; 2: sea_ice)
1226  call get_var3_r4( ncid, 'ORO', 1,im, j,j, 1,1, wk1 )
1227  tmean = 0.
1228  npt = 0
1229  do i=1,im
1230  if( abs(wk1(i)-1.) > 0.99 ) then ! ocean or sea ice
1231  tmean = tmean + wk2(i,j)
1232  npt = npt + 1
1233  endif
1234  enddo
1235 !------------------------------------------------------
1236 ! Replace TS over interior land with zonal mean SST/Ice
1237 !------------------------------------------------------
1238  if ( npt /= 0 ) then
1239  tmean= tmean / real(npt)
1240  do i=1,im
1241  if( abs(wk1(i)-1.) <= 0.99 ) then ! Land points
1242  if ( i==1 ) then
1243  i1 = im; i2 = 2
1244  elseif ( i==im ) then
1245  i1 = im-1; i2 = 1
1246  else
1247  i1 = i-1; i2 = i+1
1248  endif
1249  if ( abs(wk1(i2)-1.)>0.99 ) then ! east side has priority
1250  wk2(i,j) = wk2(i2,j)
1251  elseif ( abs(wk1(i1)-1.)>0.99 ) then ! west side
1252  wk2(i,j) = wk2(i1,j)
1253  else
1254  wk2(i,j) = tmean
1255  endif
1256  endif
1257  enddo
1258  endif
1259  enddo ! j-loop
1260  deallocate ( wk1 )
1261  endif !(.not.land_ts)
1262 
1263  do j=js,je
1264  do i=is,ie
1265  i1 = id1(i,j)
1266  i2 = id2(i,j)
1267  j1 = jdc(i,j)
1268  atm(1)%ts(i,j) = s2c(i,j,1)*wk2(i1,j1 ) + s2c(i,j,2)*wk2(i2,j1 ) + &
1269  s2c(i,j,3)*wk2(i2,j1+1) + s2c(i,j,4)*wk2(i1,j1+1)
1270  enddo
1271  enddo
1272  call prt_maxmin('SST_model', atm(1)%ts, is, ie, js, je, 0, 1, 1.)
1273 
1274 ! Perform interp to FMS SST format/grid
1275 #ifndef DYCORE_SOLO
1276  call ncep2fms(im, jm, lon, lat, wk2)
1277  if( is_master() ) then
1278  write(*,*) 'External_ic_mod: i_sst=', i_sst, ' j_sst=', j_sst
1279  call pmaxmin( 'SST_ncep_fms', sst_ncep, i_sst, j_sst, 1.)
1280  endif
1281 #endif
1282  endif !(read_ts)
1283 
1284  deallocate ( wk2 )
1285 
1286 ! Read in temperature:
1287  allocate ( wk3(1:im,jbeg:jend, 1:km) )
1288  call get_var3_r4( ncid, 'T', 1,im, jbeg,jend, 1,km, wk3 )
1289 
1290  allocate ( tp(is:ie,js:je,km) )
1291  do k=1,km
1292  do j=js,je
1293  do i=is,ie
1294  i1 = id1(i,j)
1295  i2 = id2(i,j)
1296  j1 = jdc(i,j)
1297  tp(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1298  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1299  enddo
1300  enddo
1301  enddo
1302 
1303 ! Read in tracers: only sphum at this point
1304  call get_var3_r4( ncid, 'Q', 1,im, jbeg,jend, 1,km, wk3 )
1305 
1306  allocate ( qp(is:ie,js:je,km) )
1307  do k=1,km
1308  do j=js,je
1309  do i=is,ie
1310  i1 = id1(i,j)
1311  i2 = id2(i,j)
1312  j1 = jdc(i,j)
1313  qp(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1314  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1315  enddo
1316  enddo
1317  enddo
1318 
1319  call remap_scalar(im, jm, km, npz, nq, nq, ak0, bk0, psc, gzc, tp, qp, atm(1))
1320  deallocate ( tp )
1321  deallocate ( qp )
1322 
1323 ! Winds:
1324  call get_var3_r4( ncid, 'U', 1,im, jbeg,jend, 1,km, wk3 )
1325 
1326  allocate ( ua(is:ie,js:je,km) )
1327  do k=1,km
1328  do j=js,je
1329  do i=is,ie
1330  i1 = id1(i,j)
1331  i2 = id2(i,j)
1332  j1 = jdc(i,j)
1333  ua(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1334  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1335  enddo
1336  enddo
1337  enddo
1338 
1339  call get_var3_r4( ncid, 'V', 1,im, jbeg,jend, 1,km, wk3 )
1340  call close_ncfile ( ncid )
1341 
1342  allocate ( va(is:ie,js:je,km) )
1343  do k=1,km
1344  do j=js,je
1345  do i=is,ie
1346  i1 = id1(i,j)
1347  i2 = id2(i,j)
1348  j1 = jdc(i,j)
1349  va(i,j,k) = s2c(i,j,1)*wk3(i1,j1 ,k) + s2c(i,j,2)*wk3(i2,j1 ,k) + &
1350  s2c(i,j,3)*wk3(i2,j1+1,k) + s2c(i,j,4)*wk3(i1,j1+1,k)
1351  enddo
1352  enddo
1353  enddo
1354  deallocate ( wk3 )
1355  call remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, atm(1))
1356 
1357  deallocate ( ua )
1358  deallocate ( va )
1359 
1360  deallocate ( ak0 )
1361  deallocate ( bk0 )
1362  deallocate ( lat )
1363  deallocate ( lon )
1364 
1365  end subroutine get_ncep_ic
1366 
1370  subroutine get_ecmwf_ic( Atm, fv_domain )
1372 #ifdef __PGI
1373  use gfs_restart, only : gfs_restart_type
1374 
1375  implicit none
1376 #endif
1377 
1378  type(fv_atmos_type), intent(inout) :: Atm(:)
1379  type(domain2d), intent(inout) :: fv_domain
1380 ! local:
1381  real :: ak_ec(138), bk_ec(138)
1382  data ak_ec/ 0.000000, 2.000365, 3.102241, 4.666084, 6.827977, 9.746966, &
1383  13.605424, 18.608931, 24.985718, 32.985710, 42.879242, 54.955463, &
1384  69.520576, 86.895882, 107.415741, 131.425507, 159.279404, 191.338562, &
1385  227.968948, 269.539581, 316.420746, 368.982361, 427.592499, 492.616028, &
1386  564.413452, 643.339905, 729.744141, 823.967834, 926.344910, 1037.201172, &
1387  1156.853638, 1285.610352, 1423.770142, 1571.622925, 1729.448975, 1897.519287, &
1388  2076.095947, 2265.431641, 2465.770508, 2677.348145, 2900.391357, 3135.119385, &
1389  3381.743652, 3640.468262, 3911.490479, 4194.930664, 4490.817383, 4799.149414, &
1390  5119.895020, 5452.990723, 5798.344727, 6156.074219, 6526.946777, 6911.870605, &
1391  7311.869141, 7727.412109, 8159.354004, 8608.525391, 9076.400391, 9562.682617, &
1392  10065.978516, 10584.631836, 11116.662109, 11660.067383, 12211.547852, 12766.873047, &
1393  13324.668945, 13881.331055, 14432.139648, 14975.615234, 15508.256836, 16026.115234, &
1394  16527.322266, 17008.789063, 17467.613281, 17901.621094, 18308.433594, 18685.718750, &
1395  19031.289063, 19343.511719, 19620.042969, 19859.390625, 20059.931641, 20219.664063, &
1396  20337.863281, 20412.308594, 20442.078125, 20425.718750, 20361.816406, 20249.511719, &
1397  20087.085938, 19874.025391, 19608.572266, 19290.226563, 18917.460938, 18489.707031, &
1398  18006.925781, 17471.839844, 16888.687500, 16262.046875, 15596.695313, 14898.453125, &
1399  14173.324219, 13427.769531, 12668.257813, 11901.339844, 11133.304688, 10370.175781, &
1400  9617.515625, 8880.453125, 8163.375000, 7470.343750, 6804.421875, 6168.531250, &
1401  5564.382813, 4993.796875, 4457.375000, 3955.960938, 3489.234375, 3057.265625, &
1402  2659.140625, 2294.242188, 1961.500000, 1659.476563, 1387.546875, 1143.250000, &
1403  926.507813, 734.992188, 568.062500, 424.414063, 302.476563, 202.484375, &
1404  122.101563, 62.781250, 22.835938, 3.757813, 0.000000, 0.000000 /
1405 
1406  data bk_ec/ 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1407  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1408  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1409  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1410  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1411  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1412  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1413  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1414  0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, &
1415  0.000000, 0.000007, 0.000024, 0.000059, 0.000112, 0.000199, &
1416  0.000340, 0.000562, 0.000890, 0.001353, 0.001992, 0.002857, &
1417  0.003971, 0.005378, 0.007133, 0.009261, 0.011806, 0.014816, &
1418  0.018318, 0.022355, 0.026964, 0.032176, 0.038026, 0.044548, &
1419  0.051773, 0.059728, 0.068448, 0.077958, 0.088286, 0.099462, &
1420  0.111505, 0.124448, 0.138313, 0.153125, 0.168910, 0.185689, &
1421  0.203491, 0.222333, 0.242244, 0.263242, 0.285354, 0.308598, &
1422  0.332939, 0.358254, 0.384363, 0.411125, 0.438391, 0.466003, &
1423  0.493800, 0.521619, 0.549301, 0.576692, 0.603648, 0.630036, &
1424  0.655736, 0.680643, 0.704669, 0.727739, 0.749797, 0.770798, &
1425  0.790717, 0.809536, 0.827256, 0.843881, 0.859432, 0.873929, &
1426  0.887408, 0.899900, 0.911448, 0.922096, 0.931881, 0.940860, &
1427  0.949064, 0.956550, 0.963352, 0.969513, 0.975078, 0.980072, &
1428  0.984542, 0.988500, 0.991984, 0.995003, 0.997630, 1.000000 /
1429 
1430 ! The following L63 will be used in the model
1431 ! The setting is the same as NCEP GFS's L64 except the top layer
1432  real, dimension(64):: ak_sj, bk_sj
1433  data ak_sj/64.247, 137.790, 221.958, &
1434  318.266, 428.434, 554.424, &
1435  698.457, 863.05803, 1051.07995, &
1436  1265.75194, 1510.71101, 1790.05098, &
1437  2108.36604, 2470.78817, 2883.03811, &
1438  3351.46002, 3883.05187, 4485.49315, &
1439  5167.14603, 5937.04991, 6804.87379, &
1440  7780.84698, 8875.64338, 10100.20534, &
1441  11264.35673, 12190.64366, 12905.42546, &
1442  13430.87867, 13785.88765, 13986.77987, &
1443  14047.96335, 13982.46770, 13802.40331, &
1444  13519.33841, 13144.59486, 12689.45608, &
1445  12165.28766, 11583.57006, 10955.84778, &
1446  10293.60402, 9608.08306, 8910.07678, &
1447  8209.70131, 7516.18560, 6837.69250, &
1448  6181.19473, 5552.39653, 4955.72632, &
1449  4394.37629, 3870.38682, 3384.76586, &
1450  2937.63489, 2528.37666, 2155.78385, &
1451  1818.20722, 1513.68173, 1240.03585, &
1452  994.99144, 776.23591, 581.48797, &
1453  408.53400, 255.26520, 119.70243, 0. /
1454 
1455  data bk_sj/0.00000, 0.00000, 0.00000, &
1456  0.00000, 0.00000, 0.00000, &
1457  0.00000, 0.00000, 0.00000, &
1458  0.00000, 0.00000, 0.00000, &
1459  0.00000, 0.00000, 0.00000, &
1460  0.00000, 0.00000, 0.00000, &
1461  0.00000, 0.00000, 0.00000, &
1462  0.00000, 0.00000, 0.00000, &
1463  0.00201, 0.00792, 0.01755, &
1464  0.03079, 0.04751, 0.06761, &
1465  0.09097, 0.11746, 0.14690, &
1466  0.17911, 0.21382, 0.25076, &
1467  0.28960, 0.32994, 0.37140, &
1468  0.41353, 0.45589, 0.49806, &
1469  0.53961, 0.58015, 0.61935, &
1470  0.65692, 0.69261, 0.72625, &
1471  0.75773, 0.78698, 0.81398, &
1472  0.83876, 0.86138, 0.88192, &
1473  0.90050, 0.91722, 0.93223, &
1474  0.94565, 0.95762, 0.96827, &
1475  0.97771, 0.98608, 0.99347, 1./
1476 
1477  character(len=128) :: fname
1478  real, allocatable:: wk2(:,:)
1479  real(kind=4), allocatable:: wk2_r4(:,:)
1480  real, dimension(:,:,:), allocatable:: ud, vd
1481  real, allocatable:: wc(:,:,:)
1482  real(kind=4), allocatable:: uec(:,:,:), vec(:,:,:), tec(:,:,:), wec(:,:,:)
1483  real(kind=4), allocatable:: psec(:,:), zsec(:,:), zhec(:,:,:), qec(:,:,:,:)
1484  real(kind=4), allocatable:: psc(:,:)
1485  real(kind=4), allocatable:: sphumec(:,:,:)
1486  real, allocatable:: psc_r8(:,:), zhc(:,:,:), qc(:,:,:,:)
1487  real, allocatable:: lat(:), lon(:), ak0(:), bk0(:)
1488  real, allocatable:: pt_c(:,:,:), pt_d(:,:,:)
1489  real:: s2c(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je,4)
1490  real:: s2c_c(atm(1)%bd%is:atm(1)%bd%ie+1,atm(1)%bd%js:atm(1)%bd%je,4)
1491  real:: s2c_d(atm(1)%bd%is:atm(1)%bd%ie,atm(1)%bd%js:atm(1)%bd%je+1,4)
1492  integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je):: &
1493  id1, id2, jdc
1494  integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie+1,Atm(1)%bd%js:Atm(1)%bd%je):: &
1495  id1_c, id2_c, jdc_c
1496  integer, dimension(Atm(1)%bd%is:Atm(1)%bd%ie,Atm(1)%bd%js:Atm(1)%bd%je+1):: &
1497  id1_d, id2_d, jdc_d
1498  real:: utmp, vtmp
1499  integer:: i, j, k, n, im, jm, km, npz, npt
1500  integer:: i1, i2, j1, ncid
1501  integer:: jbeg, jend, jn
1502  integer tsize(3)
1503  logical:: read_ts = .true.
1504  logical:: land_ts = .false.
1505  logical:: found
1506  integer :: is, ie, js, je
1507  integer :: isd, ied, jsd, jed
1508  integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel
1509 #ifdef MULTI_GASES
1510  integer :: spfo, spfo2, spfo3
1511 #else
1512  integer :: o3mr
1513 #endif
1514  real:: wt, qt, m_fac
1515  real(kind=8) :: scale_value, offset, ptmp
1516  real(kind=R_GRID), dimension(2):: p1, p2, p3
1517  real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
1518  real, allocatable:: ps_gfs(:,:), zh_gfs(:,:,:)
1519 #ifdef MULTI_GASES
1520  real, allocatable:: spfo_gfs(:,:,:), spfo2_gfs(:,:,:), spfo3_gfs(:,:,:)
1521 #else
1522  real, allocatable:: o3mr_gfs(:,:,:)
1523 #endif
1524  real, allocatable:: ak_gfs(:), bk_gfs(:)
1525  integer :: id_res, ntprog, ntracers, ks, iq, nt
1526  character(len=64) :: tracer_name
1527  integer :: levp_gfs = 64
1528  type(restart_file_type) :: oro_restart, gfs_restart
1529  character(len=64) :: fn_oro_ics = 'oro_data.nc'
1530  character(len=64) :: fn_gfs_ics = 'gfs_data.nc'
1531  character(len=64) :: fn_gfs_ctl = 'gfs_ctrl.nc'
1532  logical :: filtered_terrain = .true.
1533  namelist /external_ic_nml/ filtered_terrain
1534 
1535  is = atm(1)%bd%is
1536  ie = atm(1)%bd%ie
1537  js = atm(1)%bd%js
1538  je = atm(1)%bd%je
1539  isd = atm(1)%bd%isd
1540  ied = atm(1)%bd%ied
1541  jsd = atm(1)%bd%jsd
1542  jed = atm(1)%bd%jed
1543 
1544  deg2rad = pi/180.
1545 
1546  npz = atm(1)%npz
1547  call get_number_tracers(model_atmos, num_tracers=ntracers, num_prog=ntprog)
1548  if(is_master()) write(*,*) 'ntracers = ', ntracers, 'ntprog = ',ntprog
1549 
1550  sphum = get_tracer_index(model_atmos, 'sphum')
1551  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
1552  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
1553  rainwat = get_tracer_index(model_atmos, 'rainwat')
1554  snowwat = get_tracer_index(model_atmos, 'snowwat')
1555  graupel = get_tracer_index(model_atmos, 'graupel')
1556 #ifdef MULTI_GASES
1557  spfo = get_tracer_index(model_atmos, 'spfo')
1558  spfo2 = get_tracer_index(model_atmos, 'spfo2')
1559  spfo3 = get_tracer_index(model_atmos, 'spfo3')
1560 #else
1561  o3mr = get_tracer_index(model_atmos, 'o3mr')
1562 #endif
1563 
1564  if (is_master()) then
1565  print *, 'sphum = ', sphum
1566  print *, 'liq_wat = ', liq_wat
1567  if ( atm(1)%flagstruct%nwat .eq. 6 ) then
1568  print *, 'rainwat = ', rainwat
1569  print *, 'iec_wat = ', ice_wat
1570  print *, 'snowwat = ', snowwat
1571  print *, 'graupel = ', graupel
1572  endif
1573 #ifdef MULTI_GASES
1574  print *, ' spfo3 = ', spfo3
1575  print *, ' spfo = ', spfo
1576  print *, ' spfo2 = ', spfo2
1577 #else
1578  print *, ' o3mr = ', o3mr
1579 #endif
1580  endif
1581 
1582 
1583 ! Set up model's ak and bk
1584  if ( npz <= 64 ) then
1585  atm(1)%ak(:) = ak_sj(:)
1586  atm(1)%bk(:) = bk_sj(:)
1587  atm(1)%ptop = atm(1)%ak(1)
1588  else
1589  call set_eta(npz, ks, atm(1)%ptop, atm(1)%ak, atm(1)%bk)
1590  endif
1591 
1592 !! Read in model terrain from oro_data.tile?.nc
1593  if (filtered_terrain) then
1594  id_res = register_restart_field(oro_restart, fn_oro_ics, 'orog_filt', atm(1)%phis, domain=atm(1)%domain)
1595  elseif (.not. filtered_terrain) then
1596  id_res = register_restart_field(oro_restart, fn_oro_ics, 'orog_raw', atm(1)%phis, domain=atm(1)%domain)
1597  endif
1598  call restore_state (oro_restart)
1599  call free_restart_type(oro_restart)
1600  atm(1)%phis = atm(1)%phis*grav
1601  if(is_master()) write(*,*) 'done reading model terrain from oro_data.nc'
1602  call mpp_update_domains( atm(1)%phis, atm(1)%domain )
1603 
1604 !! Read in o3mr, ps and zh from GFS_data.tile?.nc
1605 #ifdef MULTI_GASES
1606  allocate (spfo3_gfs(is:ie,js:je,levp_gfs))
1607  allocate ( spfo_gfs(is:ie,js:je,levp_gfs))
1608  allocate (spfo2_gfs(is:ie,js:je,levp_gfs))
1609 #else
1610  allocate (o3mr_gfs(is:ie,js:je,levp_gfs))
1611 #endif
1612  allocate (ps_gfs(is:ie,js:je))
1613  allocate (zh_gfs(is:ie,js:je,levp_gfs+1))
1614 
1615 #ifdef MULTI_GASES
1616  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'spfo3', spfo3_gfs, &
1617  mandatory=.false.,domain=atm(1)%domain)
1618  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'spfo', spfo_gfs, &
1619  mandatory=.false.,domain=atm(1)%domain)
1620  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'spfo2', spfo2_gfs, &
1621  mandatory=.false.,domain=atm(1)%domain)
1622 #else
1623  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'o3mr', o3mr_gfs, &
1624  mandatory=.false.,domain=atm(1)%domain)
1625 #endif
1626  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'ps', ps_gfs, domain=atm(1)%domain)
1627  id_res = register_restart_field(gfs_restart, fn_gfs_ics, 'ZH', zh_gfs, domain=atm(1)%domain)
1628  call restore_state (gfs_restart)
1629  call free_restart_type(gfs_restart)
1630 
1631 
1632  ! Get GFS ak, bk for o3mr vertical interpolation
1633  allocate (wk2(levp_gfs+1,2))
1634  allocate (ak_gfs(levp_gfs+1))
1635  allocate (bk_gfs(levp_gfs+1))
1636  call read_data('INPUT/'//trim(fn_gfs_ctl),'vcoord',wk2, no_domain=.true.)
1637  ak_gfs(1:levp_gfs+1) = wk2(1:levp_gfs+1,1)
1638  bk_gfs(1:levp_gfs+1) = wk2(1:levp_gfs+1,2)
1639  deallocate (wk2)
1640 
1641  if ( bk_gfs(1) < 1.e-9 ) ak_gfs(1) = max(1.e-9, ak_gfs(1))
1642 
1643 #ifdef MULTI_GASES
1644  iq = spfo3
1645  if(is_master()) write(*,*) 'Reading spfo3 from GFS_data.nc:'
1646  if(is_master()) write(*,*) 'spfo3 =', iq
1647  call remap_scalar_single(atm(1), levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo3_gfs, zh_gfs, iq)
1648  iq = spfo
1649  if(is_master()) write(*,*) 'Reading spfo from GFS_data.nc:'
1650  if(is_master()) write(*,*) 'spfo =', iq
1651  call remap_scalar_single(atm(1), levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo_gfs, zh_gfs, iq)
1652  iq = spfo2
1653  if(is_master()) write(*,*) 'Reading spfo2 from GFS_data.nc:'
1654  if(is_master()) write(*,*) 'spfo2 =', iq
1655  call remap_scalar_single(atm(1), levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, spfo2_gfs, zh_gfs, iq)
1656 #else
1657  iq = o3mr
1658  if(is_master()) write(*,*) 'Reading o3mr from GFS_data.nc:'
1659  if(is_master()) write(*,*) 'o3mr =', iq
1660  call remap_scalar_single(atm(1), levp_gfs, npz, ak_gfs, bk_gfs, ps_gfs, o3mr_gfs, zh_gfs, iq)
1661 #endif
1662 
1663  deallocate (ak_gfs, bk_gfs)
1664  deallocate (ps_gfs, zh_gfs)
1665 #ifdef MULTI_GASES
1666  deallocate (spfo3_gfs)
1667  deallocate ( spfo_gfs)
1668  deallocate (spfo2_gfs)
1669 #else
1670  deallocate (o3mr_gfs)
1671 #endif
1672 
1673 !! Start to read EC data
1674  fname = atm(1)%flagstruct%res_latlon_dynamics
1675 
1676  if( file_exist(fname) ) then
1677  call open_ncfile( fname, ncid ) ! open the file
1678 
1679  call get_ncdim1( ncid, 'longitude', tsize(1) )
1680  call get_ncdim1( ncid, 'latitude', tsize(2) )
1681  call get_ncdim1( ncid, 'level', tsize(3) )
1682 
1683  im = tsize(1); jm = tsize(2); km = tsize(3)
1684 
1685  if(is_master()) write(*,*) fname
1686  if(is_master()) write(*,*) ' ECMWF IC dimensions:', tsize
1687 
1688  allocate ( lon(im) )
1689  allocate ( lat(jm) )
1690 
1691  call _get_var1(ncid, 'longitude', im, lon )
1692  call _get_var1(ncid, 'latitude', jm, lat )
1693 
1694 !! Convert to radian
1695  do i = 1, im
1696  lon(i) = lon(i) * deg2rad ! lon(1) = 0.
1697  enddo
1698  do j = 1, jm
1699  lat(j) = lat(j) * deg2rad
1700  enddo
1701 
1702  allocate ( ak0(km+1) )
1703  allocate ( bk0(km+1) )
1704 
1705 ! The ECMWF data from does not contain (ak,bk)
1706  do k=1, km+1
1707  ak0(k) = ak_ec(k)
1708  bk0(k) = bk_ec(k)
1709  enddo
1710 
1711  if( is_master() ) then
1712  do k=1,km+1
1713  write(*,*) k, ak0(k), bk0(k)
1714  enddo
1715  endif
1716 
1717 ! Limiter to prevent NAN at top during remapping
1718  if ( bk0(1) < 1.e-9 ) ak0(1) = max(1.e-9, ak0(1))
1719 
1720  else
1721  call mpp_error(fatal,'==> Error in get_external_ic: Expected file '//trim(fname)//' for NCEP IC does not exist')
1722  endif
1723 
1724 ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
1725  call remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
1726  im, jm, lon, lat, id1, id2, jdc, s2c , atm(1)%gridstruct%agrid )
1727 
1728 ! Find bounding latitudes:
1729  jbeg = jm-1; jend = 2
1730  do j=js,je
1731  do i=is,ie
1732  j1 = jdc(i,j)
1733  jbeg = min(jbeg, j1)
1734  jend = max(jend, j1+1)
1735  enddo
1736  enddo
1737 
1738  if(is_master()) write(*,*) 'jbeg, jend = ', jbeg, jend
1739 ! read in surface pressure and height:
1740  allocate ( psec(im,jbeg:jend) )
1741  allocate ( zsec(im,jbeg:jend) )
1742  allocate ( wk2_r4(im,jbeg:jend) )
1743 
1744  call get_var2_r4( ncid, 'lnsp', 1,im, jbeg,jend, wk2_r4 )
1745  call get_var_att_double ( ncid, 'lnsp', 'scale_factor', scale_value )
1746  call get_var_att_double ( ncid, 'lnsp', 'add_offset', offset )
1747  psec(:,:) = exp(wk2_r4(:,:)*scale_value + offset)
1748  if(is_master()) write(*,*) 'done reading psec'
1749 
1750  call get_var2_r4( ncid, 'z', 1,im, jbeg,jend, wk2_r4 )
1751  call get_var_att_double ( ncid, 'z', 'scale_factor', scale_value )
1752  call get_var_att_double ( ncid, 'z', 'add_offset', offset )
1753  zsec(:,:) = (wk2_r4(:,:)*scale_value + offset)/grav
1754  if(is_master()) write(*,*) 'done reading zsec'
1755 
1756  deallocate ( wk2_r4 )
1757 
1758 ! Read in temperature:
1759  allocate ( tec(1:im,jbeg:jend, 1:km) )
1760 
1761  call get_var3_r4( ncid, 't', 1,im, jbeg,jend, 1,km, tec )
1762  call get_var_att_double ( ncid, 't', 'scale_factor', scale_value )
1763  call get_var_att_double ( ncid, 't', 'add_offset', offset )
1764  tec(:,:,:) = tec(:,:,:)*scale_value + offset
1765  if(is_master()) write(*,*) 'done reading tec'
1766 
1767 ! read in specific humidity:
1768  allocate ( sphumec(1:im,jbeg:jend, 1:km) )
1769 
1770  call get_var3_r4( ncid, 'q', 1,im, jbeg,jend, 1,km, sphumec(:,:,:) )
1771  call get_var_att_double ( ncid, 'q', 'scale_factor', scale_value )
1772  call get_var_att_double ( ncid, 'q', 'add_offset', offset )
1773  sphumec(:,:,:) = sphumec(:,:,:)*scale_value + offset
1774  if(is_master()) write(*,*) 'done reading sphum ec'
1775 
1776 ! Read in other tracers from EC data and remap them into cubic sphere grid:
1777  allocate ( qec(1:im,jbeg:jend,1:km,5) )
1778 
1779  do n = 1, 5
1780  if (n == sphum) then
1781  qec(:,:,:,sphum) = sphumec(:,:,:)
1782  deallocate ( sphumec )
1783  else if (n == liq_wat) then
1784  call get_var3_r4( ncid, 'clwc', 1,im, jbeg,jend, 1,km, qec(:,:,:,liq_wat) )
1785  call get_var_att_double ( ncid, 'clwc', 'scale_factor', scale_value )
1786  call get_var_att_double ( ncid, 'clwc', 'add_offset', offset )
1787  qec(:,:,:,liq_wat) = qec(:,:,:,liq_wat)*scale_value + offset
1788  if(is_master()) write(*,*) 'done reading clwc ec'
1789  else if (n == rainwat) then
1790  call get_var3_r4( ncid, 'crwc', 1,im, jbeg,jend, 1,km, qec(:,:,:,rainwat) )
1791  call get_var_att_double ( ncid, 'crwc', 'scale_factor', scale_value )
1792  call get_var_att_double ( ncid, 'crwc', 'add_offset', offset )
1793  qec(:,:,:,rainwat) = qec(:,:,:,rainwat)*scale_value + offset
1794  if(is_master()) write(*,*) 'done reading crwc ec'
1795  else if (n == ice_wat) then
1796  call get_var3_r4( ncid, 'ciwc', 1,im, jbeg,jend, 1,km, qec(:,:,:,ice_wat) )
1797  call get_var_att_double ( ncid, 'ciwc', 'scale_factor', scale_value )
1798  call get_var_att_double ( ncid, 'ciwc', 'add_offset', offset )
1799  qec(:,:,:,ice_wat) = qec(:,:,:,ice_wat)*scale_value + offset
1800  if(is_master()) write(*,*) 'done reading ciwc ec'
1801  else if (n == snowwat) then
1802  call get_var3_r4( ncid, 'cswc', 1,im, jbeg,jend, 1,km, qec(:,:,:,snowwat) )
1803  call get_var_att_double ( ncid, 'cswc', 'scale_factor', scale_value )
1804  call get_var_att_double ( ncid, 'cswc', 'add_offset', offset )
1805  qec(:,:,:,snowwat) = qec(:,:,:,snowwat)*scale_value + offset
1806  if(is_master()) write(*,*) 'done reading cswc ec'
1807  else
1808  if(is_master()) write(*,*) 'nq is more then 5!'
1809  endif
1810 
1811  enddo
1812 
1813 
1814 !!!! Compute height on edges, zhec [ use psec, zsec, tec, sphum]
1815  allocate ( zhec(1:im,jbeg:jend, km+1) )
1816  jn = jend - jbeg + 1
1817 
1818  call compute_zh(im, jn, km, ak0, bk0, psec, zsec, tec, qec, 5, zhec )
1819  if(is_master()) write(*,*) 'done compute zhec'
1820 
1821 ! convert zhec, psec, zsec from EC grid to cubic grid
1822  allocate (psc(is:ie,js:je))
1823  allocate (psc_r8(is:ie,js:je))
1824 
1825 #ifdef LOGP_INTP
1826  do j=jbeg,jend
1827  do i=1,im
1828  psec(i,j) = log(psec(i,j))
1829  enddo
1830  enddo
1831 #endif
1832  do j=js,je
1833  do i=is,ie
1834  i1 = id1(i,j)
1835  i2 = id2(i,j)
1836  j1 = jdc(i,j)
1837 #ifdef LOGP_INTP
1838  ptmp = s2c(i,j,1)*psec(i1,j1 ) + s2c(i,j,2)*psec(i2,j1 ) + &
1839  s2c(i,j,3)*psec(i2,j1+1) + s2c(i,j,4)*psec(i1,j1+1)
1840  psc(i,j) = exp(ptmp)
1841 #else
1842  psc(i,j) = s2c(i,j,1)*psec(i1,j1 ) + s2c(i,j,2)*psec(i2,j1 ) + &
1843  s2c(i,j,3)*psec(i2,j1+1) + s2c(i,j,4)*psec(i1,j1+1)
1844 #endif
1845  enddo
1846  enddo
1847  deallocate ( psec )
1848  deallocate ( zsec )
1849 
1850  allocate (zhc(is:ie,js:je,km+1))
1851 !$OMP parallel do default(none) shared(is,ie,js,je,km,s2c,id1,id2,jdc,zhc,zhec) &
1852 !$OMP private(i1,i2,j1)
1853  do k=1,km+1
1854  do j=js,je
1855  do i=is,ie
1856  i1 = id1(i,j)
1857  i2 = id2(i,j)
1858  j1 = jdc(i,j)
1859  zhc(i,j,k) = s2c(i,j,1)*zhec(i1,j1 ,k) + s2c(i,j,2)*zhec(i2,j1 ,k) + &
1860  s2c(i,j,3)*zhec(i2,j1+1,k) + s2c(i,j,4)*zhec(i1,j1+1,k)
1861  enddo
1862  enddo
1863  enddo
1864  deallocate ( zhec )
1865 
1866  if(is_master()) write(*,*) 'done interpolate psec/zsec/zhec into cubic grid psc/zhc!'
1867 
1868 ! Read in other tracers from EC data and remap them into cubic sphere grid:
1869  allocate ( qc(is:ie,js:je,km,6) )
1870 
1871  do n = 1, 5
1872 !$OMP parallel do default(none) shared(n,is,ie,js,je,km,s2c,id1,id2,jdc,qc,qec) &
1873 !$OMP private(i1,i2,j1)
1874  do k=1,km
1875  do j=js,je
1876  do i=is,ie
1877  i1 = id1(i,j)
1878  i2 = id2(i,j)
1879  j1 = jdc(i,j)
1880  qc(i,j,k,n) = s2c(i,j,1)*qec(i1,j1 ,k,n) + s2c(i,j,2)*qec(i2,j1 ,k,n) + &
1881  s2c(i,j,3)*qec(i2,j1+1,k,n) + s2c(i,j,4)*qec(i1,j1+1,k,n)
1882  enddo
1883  enddo
1884  enddo
1885  enddo
1886 
1887  qc(:,:,:,graupel) = 0. ! note Graupel must be tracer #6
1888 
1889  deallocate ( qec )
1890  if(is_master()) write(*,*) 'done interpolate tracers (qec) into cubic (qc)'
1891 
1892 ! Read in vertical wind from EC data and remap them into cubic sphere grid:
1893  allocate ( wec(1:im,jbeg:jend, 1:km) )
1894  allocate ( wc(is:ie,js:je,km))
1895 
1896  call get_var3_r4( ncid, 'w', 1,im, jbeg,jend, 1,km, wec )
1897  call get_var_att_double ( ncid, 'w', 'scale_factor', scale_value )
1898  call get_var_att_double ( ncid, 'w', 'add_offset', offset )
1899  wec(:,:,:) = wec(:,:,:)*scale_value + offset
1900  !call p_maxmin('wec', wec, 1, im, jbeg, jend, km, 1.)
1901 
1902 !$OMP parallel do default(none) shared(is,ie,js,je,km,id1,id2,jdc,s2c,wc,wec) &
1903 !$OMP private(i1,i2,j1)
1904  do k=1,km
1905  do j=js,je
1906  do i=is,ie
1907  i1 = id1(i,j)
1908  i2 = id2(i,j)
1909  j1 = jdc(i,j)
1910  wc(i,j,k) = s2c(i,j,1)*wec(i1,j1 ,k) + s2c(i,j,2)*wec(i2,j1 ,k) + &
1911  s2c(i,j,3)*wec(i2,j1+1,k) + s2c(i,j,4)*wec(i1,j1+1,k)
1912  enddo
1913  enddo
1914  enddo
1915  !call p_maxmin('wc', wc, is, ie, js, je, km, 1.)
1916 
1917  deallocate ( wec )
1918  if(is_master()) write(*,*) 'done reading and interpolate vertical wind (w) into cubic'
1919 
1920 ! remap tracers
1921  psc_r8(:,:) = psc(:,:)
1922  deallocate ( psc )
1923 
1924  call remap_scalar_ec(atm(1), km, npz, 6, ak0, bk0, psc_r8, qc, wc, zhc )
1925  call mpp_update_domains(atm(1)%phis, atm(1)%domain)
1926  if(is_master()) write(*,*) 'done remap_scalar_ec'
1927 
1928  deallocate ( zhc )
1929  deallocate ( wc )
1930  deallocate ( qc )
1931 
1932 !! Winds:
1933  ! get lat/lon values of pt_c and pt_d from grid data (pt_b)
1934  allocate (pt_c(isd:ied+1,jsd:jed ,2))
1935  allocate (pt_d(isd:ied ,jsd:jed+1,2))
1936  allocate (ud(is:ie , js:je+1, km))
1937  allocate (vd(is:ie+1, js:je , km))
1938 
1939  call get_staggered_grid( is, ie, js, je, &
1940  isd, ied, jsd, jed, &
1941  atm(1)%gridstruct%grid, pt_c, pt_d)
1942 
1943  !------ pt_c part ------
1944  ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
1945  call remap_coef( is, ie+1, js, je, isd, ied+1, jsd, jed, &
1946  im, jm, lon, lat, id1_c, id2_c, jdc_c, s2c_c, pt_c)
1947 
1948  ! Find bounding latitudes:
1949  jbeg = jm-1; jend = 2
1950  do j=js,je
1951  do i=is,ie+1
1952  j1 = jdc_c(i,j)
1953  jbeg = min(jbeg, j1)
1954  jend = max(jend, j1+1)
1955  enddo
1956  enddo
1957 
1958  ! read in EC wind data
1959  allocate ( uec(1:im,jbeg:jend, 1:km) )
1960  allocate ( vec(1:im,jbeg:jend, 1:km) )
1961 
1962  call get_var3_r4( ncid, 'u', 1,im, jbeg,jend, 1,km, uec )
1963  call get_var_att_double ( ncid, 'u', 'scale_factor', scale_value )
1964  call get_var_att_double ( ncid, 'u', 'add_offset', offset )
1965  do k=1,km
1966  do j=jbeg, jend
1967  do i=1,im
1968  uec(i,j,k) = uec(i,j,k)*scale_value + offset
1969  enddo
1970  enddo
1971  enddo
1972  if(is_master()) write(*,*) 'first time done reading uec'
1973 
1974  call get_var3_r4( ncid, 'v', 1,im, jbeg,jend, 1,km, vec )
1975  call get_var_att_double ( ncid, 'v', 'scale_factor', scale_value )
1976  call get_var_att_double ( ncid, 'v', 'add_offset', offset )
1977  do k=1,km
1978  do j=jbeg, jend
1979  do i=1,im
1980  vec(i,j,k) = vec(i,j,k)*scale_value + offset
1981  enddo
1982  enddo
1983  enddo
1984 
1985  if(is_master()) write(*,*) 'first time done reading vec'
1986 
1987 !$OMP parallel do default(none) shared(is,ie,js,je,km,s2c_c,id1_c,id2_c,jdc_c,uec,vec,Atm,vd) &
1988 !$OMP private(i1,i2,j1,p1,p2,p3,e2,ex,ey,utmp,vtmp)
1989  do k=1,km
1990  do j=js,je
1991  do i=is,ie+1
1992  i1 = id1_c(i,j)
1993  i2 = id2_c(i,j)
1994  j1 = jdc_c(i,j)
1995  p1(:) = atm(1)%gridstruct%grid(i,j ,1:2)
1996  p2(:) = atm(1)%gridstruct%grid(i,j+1,1:2)
1997  call mid_pt_sphere(p1, p2, p3)
1998  call get_unit_vect2(p1, p2, e2)
1999  call get_latlon_vector(p3, ex, ey)
2000  utmp = s2c_c(i,j,1)*uec(i1,j1 ,k) + &
2001  s2c_c(i,j,2)*uec(i2,j1 ,k) + &
2002  s2c_c(i,j,3)*uec(i2,j1+1,k) + &
2003  s2c_c(i,j,4)*uec(i1,j1+1,k)
2004  vtmp = s2c_c(i,j,1)*vec(i1,j1 ,k) + &
2005  s2c_c(i,j,2)*vec(i2,j1 ,k) + &
2006  s2c_c(i,j,3)*vec(i2,j1+1,k) + &
2007  s2c_c(i,j,4)*vec(i1,j1+1,k)
2008  vd(i,j,k) = utmp*inner_prod(e2,ex) + vtmp*inner_prod(e2,ey)
2009  enddo
2010  enddo
2011  enddo
2012 
2013  deallocate ( uec, vec )
2014 
2015  !------ pt_d part ------
2016  ! Initialize lat-lon to Cubed bi-linear interpolation coeff:
2017  call remap_coef( is, ie, js, je+1, isd, ied, jsd, jed+1, &
2018  im, jm, lon, lat, id1_d, id2_d, jdc_d, s2c_d, pt_d)
2019  deallocate ( pt_c, pt_d )
2020 
2021  ! Find bounding latitudes:
2022  jbeg = jm-1; jend = 2
2023  do j=js,je+1
2024  do i=is,ie
2025  j1 = jdc_d(i,j)
2026  jbeg = min(jbeg, j1)
2027  jend = max(jend, j1+1)
2028  enddo
2029  enddo
2030 
2031  ! read in EC wind data
2032  allocate ( uec(1:im,jbeg:jend, 1:km) )
2033  allocate ( vec(1:im,jbeg:jend, 1:km) )
2034 
2035  call get_var3_r4( ncid, 'u', 1,im, jbeg,jend, 1,km, uec )
2036  call get_var_att_double ( ncid, 'u', 'scale_factor', scale_value )
2037  call get_var_att_double ( ncid, 'u', 'add_offset', offset )
2038  uec(:,:,:) = uec(:,:,:)*scale_value + offset
2039  if(is_master()) write(*,*) 'second time done reading uec'
2040 
2041  call get_var3_r4( ncid, 'v', 1,im, jbeg,jend, 1,km, vec )
2042  call get_var_att_double ( ncid, 'v', 'scale_factor', scale_value )
2043  call get_var_att_double ( ncid, 'v', 'add_offset', offset )
2044  vec(:,:,:) = vec(:,:,:)*scale_value + offset
2045  if(is_master()) write(*,*) 'second time done reading vec'
2046 
2047 !$OMP parallel do default(none) shared(is,ie,js,je,km,id1_d,id2_d,jdc_d,s2c_d,uec,vec,Atm,ud) &
2048 !$OMP private(i1,i2,j1,p1,p2,p3,e1,ex,ey,utmp,vtmp)
2049  do k=1,km
2050  do j=js,je+1
2051  do i=is,ie
2052  i1 = id1_d(i,j)
2053  i2 = id2_d(i,j)
2054  j1 = jdc_d(i,j)
2055  p1(:) = atm(1)%gridstruct%grid(i, j,1:2)
2056  p2(:) = atm(1)%gridstruct%grid(i+1,j,1:2)
2057  call mid_pt_sphere(p1, p2, p3)
2058  call get_unit_vect2(p1, p2, e1)
2059  call get_latlon_vector(p3, ex, ey)
2060  utmp = s2c_d(i,j,1)*uec(i1,j1 ,k) + &
2061  s2c_d(i,j,2)*uec(i2,j1 ,k) + &
2062  s2c_d(i,j,3)*uec(i2,j1+1,k) + &
2063  s2c_d(i,j,4)*uec(i1,j1+1,k)
2064  vtmp = s2c_d(i,j,1)*vec(i1,j1 ,k) + &
2065  s2c_d(i,j,2)*vec(i2,j1 ,k) + &
2066  s2c_d(i,j,3)*vec(i2,j1+1,k) + &
2067  s2c_d(i,j,4)*vec(i1,j1+1,k)
2068  ud(i,j,k) = utmp*inner_prod(e1,ex) + vtmp*inner_prod(e1,ey)
2069  enddo
2070  enddo
2071  enddo
2072  deallocate ( uec, vec )
2073 
2074  call remap_dwinds(km, npz, ak0, bk0, psc_r8, ud, vd, atm(1))
2075  deallocate ( ud, vd )
2076 
2077 #ifndef COND_IFS_IC
2078 ! Add cloud condensate from IFS to total MASS
2079 ! Adjust the mixing ratios consistently...
2080  do k=1,npz
2081  do j=js,je
2082  do i=is,ie
2083  wt = atm(1)%delp(i,j,k)
2084  if ( atm(1)%flagstruct%nwat .eq. 2 ) then
2085  qt = wt*(1.+atm(1)%q(i,j,k,liq_wat))
2086  elseif ( atm(1)%flagstruct%nwat .eq. 6 ) then
2087  qt = wt*(1. + atm(1)%q(i,j,k,liq_wat) + &
2088  atm(1)%q(i,j,k,ice_wat) + &
2089  atm(1)%q(i,j,k,rainwat) + &
2090  atm(1)%q(i,j,k,snowwat) + &
2091  atm(1)%q(i,j,k,graupel))
2092  endif
2093  m_fac = wt / qt
2094  do iq=1,ntracers
2095  atm(1)%q(i,j,k,iq) = m_fac * atm(1)%q(i,j,k,iq)
2096  enddo
2097  atm(1)%delp(i,j,k) = qt
2098  enddo
2099  enddo
2100  enddo
2101 #endif
2102 
2103  deallocate ( ak0, bk0 )
2104 ! deallocate ( psc )
2105  deallocate ( psc_r8 )
2106  deallocate ( lat, lon )
2107 
2108  atm(1)%flagstruct%make_nh = .false.
2109 
2110  end subroutine get_ecmwf_ic
2111 !------------------------------------------------------------------
2112 !------------------------------------------------------------------
2113  subroutine get_fv_ic( Atm, fv_domain, nq )
2114  type(fv_atmos_type), intent(inout) :: Atm(:)
2115  type(domain2d), intent(inout) :: fv_domain
2116  integer, intent(in):: nq
2117 
2118  character(len=128) :: fname, tracer_name
2119  real, allocatable:: ps0(:,:), gz0(:,:), u0(:,:,:), v0(:,:,:), t0(:,:,:), dp0(:,:,:), q0(:,:,:,:)
2120  real, allocatable:: ua(:,:,:), va(:,:,:)
2121  real, allocatable:: lat(:), lon(:), ak0(:), bk0(:)
2122  integer :: i, j, k, im, jm, km, npz, tr_ind
2123  integer tsize(3)
2124 ! integer sphum, liq_wat, ice_wat, cld_amt ! GFDL AM2 physics
2125  logical found
2126 
2127  npz = atm(1)%npz
2128 
2129 ! Zero out all initial tracer fields:
2130  atm(1)%q = 0.
2131 
2132 ! Read in lat-lon FV core restart file
2133  fname = atm(1)%flagstruct%res_latlon_dynamics
2134 
2135  if( file_exist(fname) ) then
2136  call field_size(fname, 'T', tsize, field_found=found)
2137  if(is_master()) write(*,*) 'Using lat-lon FV restart:', fname
2138 
2139  if ( found ) then
2140  im = tsize(1); jm = tsize(2); km = tsize(3)
2141  if(is_master()) write(*,*) 'External IC dimensions:', tsize
2142  else
2143  call mpp_error(fatal,'==> Error in get_external_ic: field not found')
2144  endif
2145 
2146 ! Define the lat-lon coordinate:
2147  allocate ( lon(im) )
2148  allocate ( lat(jm) )
2149 
2150  do i=1,im
2151  lon(i) = (0.5 + real(i-1)) * 2.*pi/real(im)
2152  enddo
2153 
2154  do j=1,jm
2155  lat(j) = -0.5*pi + real(j-1)*pi/real(jm-1) ! SP to NP
2156  enddo
2157 
2158  allocate ( ak0(1:km+1) )
2159  allocate ( bk0(1:km+1) )
2160  allocate ( ps0(1:im,1:jm) )
2161  allocate ( gz0(1:im,1:jm) )
2162  allocate ( u0(1:im,1:jm,1:km) )
2163  allocate ( v0(1:im,1:jm,1:km) )
2164  allocate ( t0(1:im,1:jm,1:km) )
2165  allocate ( dp0(1:im,1:jm,1:km) )
2166 
2167  call read_data (fname, 'ak', ak0)
2168  call read_data (fname, 'bk', bk0)
2169  call read_data (fname, 'Surface_geopotential', gz0)
2170  call read_data (fname, 'U', u0)
2171  call read_data (fname, 'V', v0)
2172  call read_data (fname, 'T', t0)
2173  call read_data (fname, 'DELP', dp0)
2174 
2175 ! Share the load
2176  if(is_master()) call pmaxmin( 'ZS_data', gz0, im, jm, 1./grav)
2177  if(mpp_pe()==1) call pmaxmin( 'U_data', u0, im*jm, km, 1.)
2178  if(mpp_pe()==1) call pmaxmin( 'V_data', v0, im*jm, km, 1.)
2179  if(mpp_pe()==2) call pmaxmin( 'T_data', t0, im*jm, km, 1.)
2180  if(mpp_pe()==3) call pmaxmin( 'DEL-P', dp0, im*jm, km, 0.01)
2181 
2182 
2183  else
2184  call mpp_error(fatal,'==> Error in get_external_ic: Expected file '//trim(fname)//' for dynamics does not exist')
2185  endif
2186 
2187 ! Read in tracers: only AM2 "physics tracers" at this point
2188  fname = atm(1)%flagstruct%res_latlon_tracers
2189 
2190  if( file_exist(fname) ) then
2191  if(is_master()) write(*,*) 'Using lat-lon tracer restart:', fname
2192 
2193  allocate ( q0(im,jm,km,atm(1)%ncnst) )
2194  q0 = 0.
2195 
2196  do tr_ind = 1, nq
2197  call get_tracer_names(model_atmos, tr_ind, tracer_name)
2198  if (field_exist(fname,tracer_name)) then
2199  call read_data(fname, tracer_name, q0(1:im,1:jm,1:km,tr_ind))
2200  call mpp_error(note,'==> Have read tracer '//trim(tracer_name)//' from '//trim(fname))
2201  cycle
2202  endif
2203  enddo
2204  else
2205  call mpp_error(fatal,'==> Error in get_external_ic: Expected file '//trim(fname)//' for tracers does not exist')
2206  endif
2207 
2208 ! D to A transform on lat-lon grid:
2209  allocate ( ua(im,jm,km) )
2210  allocate ( va(im,jm,km) )
2211 
2212  call d2a3d(u0, v0, ua, va, im, jm, km, lon)
2213 
2214  deallocate ( u0 )
2215  deallocate ( v0 )
2216 
2217  if(mpp_pe()==4) call pmaxmin( 'UA', ua, im*jm, km, 1.)
2218  if(mpp_pe()==4) call pmaxmin( 'VA', va, im*jm, km, 1.)
2219 
2220  do j=1,jm
2221  do i=1,im
2222  ps0(i,j) = ak0(1)
2223  enddo
2224  enddo
2225 
2226  do k=1,km
2227  do j=1,jm
2228  do i=1,im
2229  ps0(i,j) = ps0(i,j) + dp0(i,j,k)
2230  enddo
2231  enddo
2232  enddo
2233 
2234  if (is_master()) call pmaxmin( 'PS_data (mb)', ps0, im, jm, 0.01)
2235 
2236 ! Horizontal interpolation to the cubed sphere grid center
2237 ! remap vertically with terrain adjustment
2238 
2239  call remap_xyz( im, 1, jm, jm, km, npz, nq, atm(1)%ncnst, lon, lat, ak0, bk0, &
2240  ps0, gz0, ua, va, t0, q0, atm(1) )
2241 
2242  deallocate ( ak0 )
2243  deallocate ( bk0 )
2244  deallocate ( ps0 )
2245  deallocate ( gz0 )
2246  deallocate ( t0 )
2247  deallocate ( q0 )
2248  deallocate ( dp0 )
2249  deallocate ( ua )
2250  deallocate ( va )
2251  deallocate ( lat )
2252  deallocate ( lon )
2253 
2254  end subroutine get_fv_ic
2255 !------------------------------------------------------------------
2256 !------------------------------------------------------------------
2257 #ifndef DYCORE_SOLO
2258  subroutine ncep2fms(im, jm, lon, lat, wk)
2260  integer, intent(in):: im, jm
2261  real, intent(in):: lon(im), lat(jm)
2262  real(kind=4), intent(in):: wk(im,jm)
2263 ! local:
2264  real :: rdlon(im)
2265  real :: rdlat(jm)
2266  real:: a1, b1
2267  real:: delx, dely
2268  real:: xc, yc ! "data" location
2269  real:: c1, c2, c3, c4
2270  integer i,j, i1, i2, jc, i0, j0, it, jt
2271 
2272  do i=1,im-1
2273  rdlon(i) = 1. / (lon(i+1) - lon(i))
2274  enddo
2275  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
2276 
2277  do j=1,jm-1
2278  rdlat(j) = 1. / (lat(j+1) - lat(j))
2279  enddo
2280 
2281 ! * Interpolate to "FMS" 1x1 SST data grid
2282 ! lon: 0.5, 1.5, ..., 359.5
2283 ! lat: -89.5, -88.5, ... , 88.5, 89.5
2284 
2285  delx = 360./real(i_sst)
2286  dely = 180./real(j_sst)
2287 
2288  jt = 1
2289  do 5000 j=1,j_sst
2290 
2291  yc = (-90. + dely * (0.5+real(j-1))) * deg2rad
2292  if ( yc<lat(1) ) then
2293  jc = 1
2294  b1 = 0.
2295  elseif ( yc>lat(jm) ) then
2296  jc = jm-1
2297  b1 = 1.
2298  else
2299  do j0=jt,jm-1
2300  if ( yc>=lat(j0) .and. yc<=lat(j0+1) ) then
2301  jc = j0
2302  jt = j0
2303  b1 = (yc-lat(jc)) * rdlat(jc)
2304  go to 222
2305  endif
2306  enddo
2307  endif
2308 222 continue
2309  it = 1
2310 
2311  do i=1,i_sst
2312  xc = delx * (0.5+real(i-1)) * deg2rad
2313  if ( xc>lon(im) ) then
2314  i1 = im; i2 = 1
2315  a1 = (xc-lon(im)) * rdlon(im)
2316  elseif ( xc<lon(1) ) then
2317  i1 = im; i2 = 1
2318  a1 = (xc+2.*pi-lon(im)) * rdlon(im)
2319  else
2320  do i0=it,im-1
2321  if ( xc>=lon(i0) .and. xc<=lon(i0+1) ) then
2322  i1 = i0; i2 = i0+1
2323  it = i0
2324  a1 = (xc-lon(i1)) * rdlon(i0)
2325  go to 111
2326  endif
2327  enddo
2328  endif
2329 111 continue
2330 
2331  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
2332  write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
2333  endif
2334 
2335  c1 = (1.-a1) * (1.-b1)
2336  c2 = a1 * (1.-b1)
2337  c3 = a1 * b1
2338  c4 = (1.-a1) * b1
2339 ! Interpolated surface pressure
2340  sst_ncep(i,j) = c1*wk(i1,jc ) + c2*wk(i2,jc ) + &
2341  c3*wk(i2,jc+1) + c4*wk(i1,jc+1)
2342  enddo !i-loop
2343 5000 continue ! j-loop
2344 
2345  end subroutine ncep2fms
2346 #endif
2347 
2348 
2349  subroutine remap_coef( is, ie, js, je, isd, ied, jsd, jed, &
2350  im, jm, lon, lat, id1, id2, jdc, s2c, agrid )
2352  integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
2353  integer, intent(in):: im, jm
2354  real, intent(in):: lon(im), lat(jm)
2355  real, intent(out):: s2c(is:ie,js:je,4)
2356  integer, intent(out), dimension(is:ie,js:je):: id1, id2, jdc
2357  real, intent(in):: agrid(isd:ied,jsd:jed,2)
2358 ! local:
2359  real :: rdlon(im)
2360  real :: rdlat(jm)
2361  real:: a1, b1
2362  integer i,j, i1, i2, jc, i0, j0
2363 
2364  do i=1,im-1
2365  rdlon(i) = 1. / (lon(i+1) - lon(i))
2366  enddo
2367  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
2368 
2369  do j=1,jm-1
2370  rdlat(j) = 1. / (lat(j+1) - lat(j))
2371  enddo
2372 
2373 ! * Interpolate to cubed sphere cell center
2374  do 5000 j=js,je
2375 
2376  do i=is,ie
2377 
2378  if ( agrid(i,j,1)>lon(im) ) then
2379  i1 = im; i2 = 1
2380  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
2381  elseif ( agrid(i,j,1)<lon(1) ) then
2382  i1 = im; i2 = 1
2383  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
2384  else
2385  do i0=1,im-1
2386  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
2387  i1 = i0; i2 = i0+1
2388  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
2389  go to 111
2390  endif
2391  enddo
2392  endif
2393 111 continue
2394 
2395  if ( agrid(i,j,2)<lat(1) ) then
2396  jc = 1
2397  b1 = 0.
2398  elseif ( agrid(i,j,2)>lat(jm) ) then
2399  jc = jm-1
2400  b1 = 1.
2401  else
2402  do j0=1,jm-1
2403  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
2404  jc = j0
2405  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
2406  go to 222
2407  endif
2408  enddo
2409  endif
2410 222 continue
2411 
2412  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
2413  write(*,*) 'gid=', mpp_pe(), i,j,a1, b1
2414  endif
2415 
2416  s2c(i,j,1) = (1.-a1) * (1.-b1)
2417  s2c(i,j,2) = a1 * (1.-b1)
2418  s2c(i,j,3) = a1 * b1
2419  s2c(i,j,4) = (1.-a1) * b1
2420  id1(i,j) = i1
2421  id2(i,j) = i2
2422  jdc(i,j) = jc
2423  enddo !i-loop
2424 5000 continue ! j-loop
2425 
2426  end subroutine remap_coef
2427 
2428 
2429  subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, Atm)
2430  type(fv_atmos_type), intent(inout) :: Atm
2431  integer, intent(in):: im, jm, km, npz, nq, ncnst
2432  real, intent(in):: ak0(km+1), bk0(km+1)
2433  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: psc, gzc
2434  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: ta
2435  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km,ncnst):: qa
2436 ! local:
2437  real, dimension(Atm%bd%is:Atm%bd%ie,km):: tp
2438  real, dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0, pn0
2439  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1
2440  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1, pn1
2441  real(kind=R_GRID), dimension(2*km+1):: gz, pn
2442  real pk0(km+1)
2443  real qp(atm%bd%is:atm%bd%ie,km,ncnst)
2444  real p1, p2, alpha, rdg
2445  real(kind=R_GRID):: pst, pt0
2446 #ifdef MULTI_GASES
2447  integer spfo, spfo2, spfo3
2448 #else
2449  integer o3mr
2450 #endif
2451  integer i,j,k, k2,l, iq
2452  integer sphum, clwmr
2453  integer :: is, ie, js, je
2454  integer :: isd, ied, jsd, jed
2455 
2456  is = atm%bd%is
2457  ie = atm%bd%ie
2458  js = atm%bd%js
2459  je = atm%bd%je
2460  isd = atm%bd%isd
2461  ied = atm%bd%ied
2462  jsd = atm%bd%jsd
2463  jed = atm%bd%jed
2464 
2465  k2 = max(10, km/2)
2466 
2467 ! nq is always 1
2468  sphum = get_tracer_index(model_atmos, 'sphum')
2469 
2470  if (mpp_pe()==1) then
2471  print *, 'sphum = ', sphum, ' ncnst=', ncnst
2472  print *, 'T_is_Tv = ', t_is_tv, ' zvir=', zvir, ' kappa=', kappa
2473  endif
2474 
2475  if ( sphum/=1 ) then
2476  call mpp_error(fatal,'SPHUM must be 1st tracer')
2477  endif
2478 
2479  call prt_maxmin('ZS_FV3', atm%phis, is, ie, js, je, 3, 1, 1./grav)
2480  call prt_maxmin('ZS_GFS', gzc, is, ie, js, je, 0, 1, 1./grav)
2481  call prt_maxmin('PS_Data', psc, is, ie, js, je, 0, 1, 0.01)
2482  call prt_maxmin('T_Data', ta, is, ie, js, je, 0, km, 1.)
2483  call prt_maxmin('q_Data', qa(is:ie,js:je,1:km,1), is, ie, js, je, 0, km, 1.)
2484 
2485  do 5000 j=js,je
2486 
2487  do i=is,ie
2488 
2489  do iq=1,ncnst
2490  do k=1,km
2491  qp(i,k,iq) = qa(i,j,k,iq)
2492  enddo
2493  enddo
2494 
2495  if ( t_is_tv ) then
2496 ! The "T" field in NCEP analysis is actually virtual temperature (Larry H. post processing)
2497 ! BEFORE 20051201
2498  do k=1,km
2499  tp(i,k) = ta(i,j,k)
2500  enddo
2501  else
2502  do k=1,km
2503 #ifdef MULTI_GASES
2504  tp(i,k) = ta(i,j,k)*virq(qp(i,k,:))
2505 #else
2506  tp(i,k) = ta(i,j,k)*(1.+zvir*qp(i,k,sphum))
2507 #endif
2508  enddo
2509  endif
2510 ! Tracers:
2511 
2512  do k=1,km+1
2513  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2514  pn0(i,k) = log(pe0(i,k))
2515  pk0(k) = pe0(i,k)**kappa
2516  enddo
2517 ! gzc is geopotential
2518 
2519 ! Note the following line, gz is actully Z (from Jeff's data).
2520  gz(km+1) = gzc(i,j)
2521  do k=km,1,-1
2522  gz(k) = gz(k+1) + rdgas*tp(i,k)*(pn0(i,k+1)-pn0(i,k))
2523  enddo
2524 
2525  do k=1,km+1
2526  pn(k) = pn0(i,k)
2527  enddo
2528 ! Use log-p for interpolation/extrapolation
2529 ! mirror image method:
2530  do k=km+2, km+k2
2531  l = 2*(km+1) - k
2532  gz(k) = 2.*gz(km+1) - gz(l)
2533  pn(k) = 2.*pn(km+1) - pn(l)
2534  enddo
2535  do k=km+k2-1, 2, -1
2536  if( atm%phis(i,j).le.gz(k) .and. atm%phis(i,j).ge.gz(k+1) ) then
2537  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
2538  go to 123
2539  endif
2540  enddo
2541 123 atm%ps(i,j) = exp(pst)
2542  enddo ! i-loop
2543 
2544  do i=is,ie
2545  pe1(i,1) = atm%ak(1)
2546  pn1(i,1) = log(pe1(i,1))
2547  enddo
2548  do k=2,npz+1
2549  do i=is,ie
2550  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
2551  pn1(i,k) = log(pe1(i,k))
2552  enddo
2553  enddo
2554 
2555 ! * Compute delp
2556  do k=1,npz
2557  do i=is,ie
2558  atm%delp(i,j,k) = pe1(i,k+1) - pe1(i,k)
2559  enddo
2560  enddo
2561 
2562 !---------------
2563 ! map shpum, o3mr, clwmr tracers
2564 !----------------
2565  do iq=1,ncnst
2566  call mappm(km, pe0, qp(is,1,iq), npz, pe1, qn1, is,ie, 0, 11, atm%ptop)
2567  do k=1,npz
2568  do i=is,ie
2569  atm%q(i,j,k,iq) = qn1(i,k)
2570  enddo
2571  enddo
2572  enddo
2573 
2574 !-------------------------------------------------------------
2575 ! map virtual temperature using geopotential conserving scheme.
2576 !-------------------------------------------------------------
2577  call mappm(km, pn0, tp, npz, pn1, qn1, is,ie, 1, 9, atm%ptop)
2578  do k=1,npz
2579  do i=is,ie
2580 #ifdef MULTI_GASES
2581  atm%pt(i,j,k) = qn1(i,k)/virq(atm%q(i,j,k,:))
2582 #else
2583  atm%pt(i,j,k) = qn1(i,k)/(1.+zvir*atm%q(i,j,k,sphum))
2584 #endif
2585  enddo
2586  enddo
2587 
2588  if ( .not. atm%flagstruct%hydrostatic .and. atm%flagstruct%ncep_ic ) then
2589 ! Replace delz with NCEP hydrostatic state
2590  rdg = -rdgas / grav
2591  do k=1,npz
2592  do i=is,ie
2593  atm%delz(i,j,k) = rdg*qn1(i,k)*(pn1(i,k+1)-pn1(i,k))
2594  enddo
2595  enddo
2596  endif
2597 
2598 5000 continue
2599 
2600  call prt_maxmin('PS_model', atm%ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
2601 
2602  if (is_master()) write(*,*) 'done remap_scalar'
2603 
2604  end subroutine remap_scalar
2605 
2606 
2607  subroutine remap_scalar_nggps(Atm, km, npz, ncnst, ak0, bk0, psc, t_in, qa, omga, zh)
2608  type(fv_atmos_type), intent(inout) :: Atm
2609  integer, intent(in):: km, npz, ncnst
2610  real, intent(in):: ak0(km+1), bk0(km+1)
2611  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: psc
2612  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: t_in
2613  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: omga
2614  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km,ncnst):: qa
2615  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km+1):: zh
2616 ! local:
2617  real, dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0
2618  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1, dp2
2619  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1
2620  real qp(atm%bd%is:atm%bd%ie,km)
2621  real wk(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je)
2622  real, dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: z500
2623 !!! High-precision
2624  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pn1
2625  real(kind=R_GRID):: gz_fv(npz+1)
2626  real(kind=R_GRID), dimension(2*km+1):: gz, pn
2627  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,km+1):: pn0
2628  real(kind=R_GRID):: pst
2629 !!! High-precision
2630  integer i,j,k,l,m, k2,iq
2631  integer sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt, liq_aero, ice_aero
2632 #ifdef MULTI_GASES
2633  integer spfo, spfo2, spfo3
2634 #else
2635  integer o3mr
2636 #endif
2637  integer :: is, ie, js, je
2638 
2639  is = atm%bd%is
2640  ie = atm%bd%ie
2641  js = atm%bd%js
2642  je = atm%bd%je
2643 
2644  sphum = get_tracer_index(model_atmos, 'sphum')
2645  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
2646  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
2647  rainwat = get_tracer_index(model_atmos, 'rainwat')
2648  snowwat = get_tracer_index(model_atmos, 'snowwat')
2649  graupel = get_tracer_index(model_atmos, 'graupel')
2650  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
2651 #ifdef MULTI_GASES
2652  spfo = get_tracer_index(model_atmos, 'spfo')
2653  spfo2 = get_tracer_index(model_atmos, 'spfo2')
2654  spfo3 = get_tracer_index(model_atmos, 'spfo3')
2655 #else
2656  o3mr = get_tracer_index(model_atmos, 'o3mr')
2657 #endif
2658  liq_aero = get_tracer_index(model_atmos, 'liq_aero')
2659  ice_aero = get_tracer_index(model_atmos, 'ice_aero')
2660 
2661  k2 = max(10, km/2)
2662 
2663  if (mpp_pe()==1) then
2664  print *, 'sphum = ', sphum
2665  print *, 'clwmr = ', liq_wat
2666 #ifdef MULTI_GASES
2667  print *, 'spfo3 = ', spfo3
2668  print *, ' spfo = ', spfo
2669  print *, 'spfo2 = ', spfo2
2670 #else
2671  print *, ' o3mr = ', o3mr
2672 #endif
2673  print *, 'liq_aero = ', liq_aero
2674  print *, 'ice_aero = ', ice_aero
2675  print *, 'ncnst = ', ncnst
2676  endif
2677 
2678  if ( sphum/=1 ) then
2679  call mpp_error(fatal,'SPHUM must be 1st tracer')
2680  endif
2681 
2682 #ifdef USE_GFS_ZS
2683  atm%phis(is:ie,js:je) = zh(is:ie,js:je,km+1)*grav
2684 #endif
2685 
2686 !$OMP parallel do default(none) &
2687 !$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,liq_aero,ice_aero,source, &
2688 !$OMP cld_amt,ncnst,npz,is,ie,js,je,km,k2,ak0,bk0,psc,t_in,zh,omga,qa,Atm,z500) &
2689 !$OMP private(l,m,pst,pn,gz,pe0,pn0,pe1,pn1,dp2,qp,qn1,gz_fv)
2690  do 5000 j=js,je
2691  do k=1,km+1
2692  do i=is,ie
2693  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2694  pn0(i,k) = log(pe0(i,k))
2695  enddo
2696  enddo
2697 
2698  do i=is,ie
2699  do k=1,km+1
2700  pn(k) = pn0(i,k)
2701  gz(k) = zh(i,j,k)*grav
2702  enddo
2703 ! Use log-p for interpolation/extrapolation
2704 ! mirror image method:
2705  do k=km+2, km+k2
2706  l = 2*(km+1) - k
2707  gz(k) = 2.*gz(km+1) - gz(l)
2708  pn(k) = 2.*pn(km+1) - pn(l)
2709  enddo
2710 
2711  do k=km+k2-1, 2, -1
2712  if( atm%phis(i,j).le.gz(k) .and. atm%phis(i,j).ge.gz(k+1) ) then
2713  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
2714  go to 123
2715  endif
2716  enddo
2717 123 atm%ps(i,j) = exp(pst)
2718 
2719 ! ------------------
2720 ! Find 500-mb height
2721 ! ------------------
2722  pst = log(500.e2)
2723  do k=km+k2-1, 2, -1
2724  if( pst.le.pn(k+1) .and. pst.ge.pn(k) ) then
2725  z500(i,j) = (gz(k+1) + (gz(k)-gz(k+1))*(pn(k+1)-pst)/(pn(k+1)-pn(k)))/grav
2726  go to 124
2727  endif
2728  enddo
2729 124 continue
2730 
2731  enddo ! i-loop
2732 
2733  do i=is,ie
2734  pe1(i,1) = atm%ak(1)
2735  pn1(i,1) = log(pe1(i,1))
2736  enddo
2737  do k=2,npz+1
2738  do i=is,ie
2739  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
2740  pn1(i,k) = log(pe1(i,k))
2741  enddo
2742  enddo
2743 
2744 ! * Compute delp
2745  do k=1,npz
2746  do i=is,ie
2747  dp2(i,k) = pe1(i,k+1) - pe1(i,k)
2748  atm%delp(i,j,k) = dp2(i,k)
2749  enddo
2750  enddo
2751 
2752 ! map tracers
2753  do iq=1,ncnst
2754  do k=1,km
2755  do i=is,ie
2756  qp(i,k) = qa(i,j,k,iq)
2757  enddo
2758  enddo
2759  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 0, 8, atm%ptop)
2760  if ( iq==sphum ) then
2761  call fillq(ie-is+1, npz, 1, qn1, dp2)
2762  else
2763  call fillz(ie-is+1, npz, 1, qn1, dp2)
2764  endif
2765 ! The HiRam step of blending model sphum with NCEP data is obsolete because nggps is always cold starting...
2766  do k=1,npz
2767  do i=is,ie
2768  atm%q(i,j,k,iq) = qn1(i,k)
2769  enddo
2770  enddo
2771  enddo
2772 
2773 !---------------------------------------------------
2774 ! Retrive temperature using GFS geopotential height
2775 !---------------------------------------------------
2776  do i=is,ie
2777 ! Make sure FV3 top is lower than GFS; can not do extrapolation above the top at this point
2778  if ( pn1(i,1) .lt. pn0(i,1) ) then
2779  call mpp_error(fatal,'FV3 top higher than NCEP/GFS')
2780  endif
2781 
2782  do k=1,km+1
2783  pn(k) = pn0(i,k)
2784  gz(k) = zh(i,j,k)*grav
2785  enddo
2786 !-------------------------------------------------
2787  do k=km+2, km+k2
2788  l = 2*(km+1) - k
2789  gz(k) = 2.*gz(km+1) - gz(l)
2790  pn(k) = 2.*pn(km+1) - pn(l)
2791  enddo
2792 !-------------------------------------------------
2793 
2794  gz_fv(npz+1) = atm%phis(i,j)
2795 
2796  m = 1
2797 
2798  do k=1,npz
2799 ! Searching using FV3 log(pe): pn1
2800 #ifdef USE_ISOTHERMO
2801  do l=m,km
2802  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
2803  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
2804  goto 555
2805  elseif ( pn1(i,k) .gt. pn(km+1) ) then
2806 ! Isothermal under ground; linear in log-p extra-polation
2807  gz_fv(k) = gz(km+1) + (gz_fv(npz+1)-gz(km+1))*(pn1(i,k)-pn(km+1))/(pn1(i,npz+1)-pn(km+1))
2808  goto 555
2809  endif
2810  enddo
2811 #else
2812  do l=m,km+k2-1
2813  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
2814  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
2815  goto 555
2816  endif
2817  enddo
2818 #endif
2819 555 m = l
2820  enddo
2821 
2822  do k=1,npz+1
2823  atm%peln(i,k,j) = pn1(i,k)
2824  enddo
2825 
2826 !----------------------------------------------------
2827 ! Compute true temperature using hydrostatic balance
2828 !----------------------------------------------------
2829  if (trim(source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
2830  do k=1,npz
2831 #ifdef MULTI_GASES
2832  atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*virq(atm%q(i,j,k,:)) )
2833 #else
2834  atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*atm%q(i,j,k,sphum)) )
2835 #endif
2836  enddo
2837 !------------------------------
2838 ! Remap input T linearly in p.
2839 !------------------------------
2840  else
2841  do k=1,km
2842  qp(i,k) = t_in(i,j,k)
2843  enddo
2844 
2845  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 2, 4, atm%ptop)
2846 
2847  do k=1,npz
2848  atm%pt(i,j,k) = qn1(i,k)
2849  enddo
2850  endif
2851 
2852  if ( .not. atm%flagstruct%hydrostatic ) then
2853  do k=1,npz
2854  atm%delz(i,j,k) = (gz_fv(k+1) - gz_fv(k)) / grav
2855  enddo
2856  endif
2857 
2858  enddo ! i-loop
2859 
2860 !-----------------------------------------------------------------------
2861 ! seperate cloud water and cloud ice
2862 ! From Jan-Huey Chen's HiRAM code
2863 !-----------------------------------------------------------------------
2864  if (cld_amt .gt. 0) atm%q(:,:,:,cld_amt) = 0.
2865  if (trim(source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
2866  if ( atm%flagstruct%nwat .eq. 6 ) then
2867  do k=1,npz
2868  do i=is,ie
2869  qn1(i,k) = atm%q(i,j,k,liq_wat)
2870  atm%q(i,j,k,rainwat) = 0.
2871  atm%q(i,j,k,snowwat) = 0.
2872  atm%q(i,j,k,graupel) = 0.
2873 ! if (cld_amt .gt. 0) Atm%q(i,j,k,cld_amt) = 0.
2874  if ( atm%pt(i,j,k) > 273.16 ) then ! > 0C all liq_wat
2875  atm%q(i,j,k,liq_wat) = qn1(i,k)
2876  atm%q(i,j,k,ice_wat) = 0.
2877 #ifdef ORIG_CLOUDS_PART
2878  else if ( atm%pt(i,j,k) < 258.16 ) then ! < -15C all ice_wat
2879  atm%q(i,j,k,liq_wat) = 0.
2880  atm%q(i,j,k,ice_wat) = qn1(i,k)
2881  else ! between -15~0C: linear interpolation
2882  atm%q(i,j,k,liq_wat) = qn1(i,k)*((atm%pt(i,j,k)-258.16)/15.)
2883  atm%q(i,j,k,ice_wat) = qn1(i,k) - atm%q(i,j,k,liq_wat)
2884  endif
2885 #else
2886  else if ( atm%pt(i,j,k) < 233.16 ) then ! < -40C all ice_wat
2887  atm%q(i,j,k,liq_wat) = 0.
2888  atm%q(i,j,k,ice_wat) = qn1(i,k)
2889  else
2890  if ( k.eq.1 ) then ! between [-40,0]: linear interpolation
2891  atm%q(i,j,k,liq_wat) = qn1(i,k)*((atm%pt(i,j,k)-233.16)/40.)
2892  atm%q(i,j,k,ice_wat) = qn1(i,k) - atm%q(i,j,k,liq_wat)
2893  else
2894  if (atm%pt(i,j,k)<258.16 .and. atm%q(i,j,k-1,ice_wat)>1.e-5 ) then
2895  atm%q(i,j,k,liq_wat) = 0.
2896  atm%q(i,j,k,ice_wat) = qn1(i,k)
2897  else ! between [-40,0]: linear interpolation
2898  atm%q(i,j,k,liq_wat) = qn1(i,k)*((atm%pt(i,j,k)-233.16)/40.)
2899  atm%q(i,j,k,ice_wat) = qn1(i,k) - atm%q(i,j,k,liq_wat)
2900  endif
2901  endif
2902  endif
2903 #endif
2904  call mp_auto_conversion(atm%q(i,j,k,liq_wat), atm%q(i,j,k,rainwat), &
2905  atm%q(i,j,k,ice_wat), atm%q(i,j,k,snowwat) )
2906  enddo
2907  enddo
2908  endif
2909  endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
2910 
2911 ! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
2912 ! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
2913 ! For FV3GFS NEMSIO input, w is already in m/s (but the code reads in as omga) and just needs to be remapped
2914 !-------------------------------------------------------------
2915 ! map omega
2916 !------- ------------------------------------------------------
2917  if ( .not. atm%flagstruct%hydrostatic ) then
2918  do k=1,km
2919  do i=is,ie
2920  qp(i,k) = omga(i,j,k)
2921  enddo
2922  enddo
2923  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, atm%ptop)
2924  if (trim(source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
2925  do k=1,npz
2926  do i=is,ie
2927  atm%w(i,j,k) = qn1(i,k)
2928  enddo
2929  enddo
2930 
2931  else
2932  do k=1,npz
2933  do i=is,ie
2934  atm%w(i,j,k) = qn1(i,k)/atm%delp(i,j,k)*atm%delz(i,j,k)
2935  enddo
2936  enddo
2937  endif
2938  endif !.not. Atm%flagstruct%hydrostatic
2939 5000 continue
2940 
2941 ! Add some diagnostics:
2942  call p_maxmin('PS_model (mb)', atm%ps(is:ie,js:je), is, ie, js, je, 1, 0.01)
2943  call p_maxmin('PT_model', atm%pt(is:ie,js:je,1:npz), is, ie, js, je, npz, 1.)
2944  do j=js,je
2945  do i=is,ie
2946  wk(i,j) = atm%phis(i,j)/grav - zh(i,j,km+1)
2947  enddo
2948  enddo
2949  call pmaxmn('ZS_diff (m)', wk, is, ie, js, je, 1, 1., atm%gridstruct%area_64, atm%domain)
2950 
2951  if (.not.atm%neststruct%nested) then
2952  call prt_gb_nh_sh('GFS_IC Z500', is,ie, js,je, z500, atm%gridstruct%area_64(is:ie,js:je), atm%gridstruct%agrid_64(is:ie,js:je,2))
2953  if ( .not. atm%flagstruct%hydrostatic ) &
2954  call prt_height('fv3_IC Z500', is,ie, js,je, 3, npz, 500.e2, atm%phis, atm%delz, atm%peln, &
2955  atm%gridstruct%area_64(is:ie,js:je), atm%gridstruct%agrid_64(is:ie,js:je,2))
2956  endif
2957 
2958  do j=js,je
2959  do i=is,ie
2960  wk(i,j) = atm%ps(i,j) - psc(i,j)
2961  enddo
2962  enddo
2963  call pmaxmn('PS_diff (mb)', wk, is, ie, js, je, 1, 0.01, atm%gridstruct%area_64, atm%domain)
2964 
2965  if (is_master()) write(*,*) 'done remap_scalar_nggps'
2966 
2967  end subroutine remap_scalar_nggps
2968 
2969  subroutine remap_scalar_ec(Atm, km, npz, ncnst, ak0, bk0, psc, qa, wc, zh)
2970  type(fv_atmos_type), intent(inout) :: Atm
2971  integer, intent(in):: km, npz, ncnst
2972  real, intent(in):: ak0(km+1), bk0(km+1)
2973  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: psc
2974  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: wc
2975  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km,ncnst):: qa
2976  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km+1):: zh
2977 ! local:
2978  real, dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0
2979  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1, dp2
2980  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1
2981  real qp(atm%bd%is:atm%bd%ie,km)
2982  real wk(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je)
2983 !!! High-precision
2984  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pn1
2985  real(kind=R_GRID):: gz_fv(npz+1)
2986  real(kind=R_GRID), dimension(2*km+1):: gz, pn
2987  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,km+1):: pn0
2988  real(kind=R_GRID):: pst
2989  real, dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: z500
2990 !!! High-precision
2991  integer:: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
2992 #ifdef MULTI_GASES
2993  integer:: spfo, spfo2, spfo3
2994 #else
2995  integer:: o3mr
2996 #endif
2997  integer:: i,j,k,l,m,k2, iq
2998  integer:: is, ie, js, je
2999 
3000  is = atm%bd%is
3001  ie = atm%bd%ie
3002  js = atm%bd%js
3003  je = atm%bd%je
3004 
3005  sphum = get_tracer_index(model_atmos, 'sphum')
3006  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
3007 
3008  if ( atm%flagstruct%nwat .eq. 6 ) then
3009  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
3010  rainwat = get_tracer_index(model_atmos, 'rainwat')
3011  snowwat = get_tracer_index(model_atmos, 'snowwat')
3012  graupel = get_tracer_index(model_atmos, 'graupel')
3013  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
3014  endif
3015  if (cld_amt .gt. 0) atm%q(:,:,:,cld_amt) = 0.
3016 
3017  k2 = max(10, km/2)
3018 
3019  if (mpp_pe()==1) then
3020  print *, 'In remap_scalar_ec:'
3021  print *, 'ncnst = ', ncnst
3022  print *, 'sphum = ', sphum
3023  print *, 'liq_wat = ', liq_wat
3024  if ( atm%flagstruct%nwat .eq. 6 ) then
3025  print *, 'rainwat = ', rainwat
3026  print *, 'ice_wat = ', ice_wat
3027  print *, 'snowwat = ', snowwat
3028  print *, 'graupel = ', graupel
3029  endif
3030  endif
3031 
3032 !$OMP parallel do default(none) shared(sphum,ncnst,npz,is,ie,js,je,km,k2,ak0,bk0,psc,zh,qa,wc,Atm,z500) &
3033 !$OMP private(l,m,pst,pn,gz,pe0,pn0,pe1,pn1,dp2,qp,qn1,gz_fv)
3034  do 5000 j=js,je
3035  do k=1,km+1
3036  do i=is,ie
3037  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
3038  pn0(i,k) = log(pe0(i,k))
3039  enddo
3040  enddo
3041 
3042  do i=is,ie
3043  do k=1,km+1
3044  pn(k) = pn0(i,k)
3045  gz(k) = zh(i,j,k)*grav
3046  enddo
3047 ! Use log-p for interpolation/extrapolation
3048 ! mirror image method:
3049  do k=km+2, km+k2
3050  l = 2*(km+1) - k
3051  gz(k) = 2.*gz(km+1) - gz(l)
3052  pn(k) = 2.*pn(km+1) - pn(l)
3053  enddo
3054 
3055  do k=km+k2-1, 2, -1
3056  if( atm%phis(i,j).le.gz(k) .and. atm%phis(i,j).ge.gz(k+1) ) then
3057  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
3058  go to 123
3059  endif
3060  enddo
3061 123 atm%ps(i,j) = exp(pst)
3062 
3063 ! ------------------
3064 ! Find 500-mb height
3065 ! ------------------
3066  pst = log(500.e2)
3067  do k=km+k2-1, 2, -1
3068  if( pst.le.pn(k+1) .and. pst.ge.pn(k) ) then
3069  z500(i,j) = (gz(k+1) + (gz(k)-gz(k+1))*(pn(k+1)-pst)/(pn(k+1)-pn(k)))/grav
3070  go to 125
3071  endif
3072  enddo
3073 125 continue
3074 
3075  enddo ! i-loop
3076 
3077  do i=is,ie
3078  pe1(i,1) = atm%ak(1)
3079  pn1(i,1) = log(pe1(i,1))
3080  enddo
3081  do k=2,npz+1
3082  do i=is,ie
3083  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
3084  pn1(i,k) = log(pe1(i,k))
3085  enddo
3086  enddo
3087 
3088 ! * Compute delp
3089  do k=1,npz
3090  do i=is,ie
3091  dp2(i,k) = pe1(i,k+1) - pe1(i,k)
3092  atm%delp(i,j,k) = dp2(i,k)
3093  enddo
3094  enddo
3095 
3096 ! map shpum, liq_wat, ice_wat, rainwat, snowwat tracers
3097  do iq=1,ncnst
3098  do k=1,km
3099  do i=is,ie
3100  qp(i,k) = qa(i,j,k,iq)
3101  enddo
3102  enddo
3103  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 0, 8, atm%ptop)
3104  if ( iq==1 ) then
3105  call fillq(ie-is+1, npz, 1, qn1, dp2)
3106  else
3107  call fillz(ie-is+1, npz, 1, qn1, dp2)
3108  endif
3109 ! The HiRam step of blending model sphum with NCEP data is obsolete because nggps is always cold starting...
3110  do k=1,npz
3111  do i=is,ie
3112  atm%q(i,j,k,iq) = qn1(i,k)
3113  enddo
3114  enddo
3115  enddo
3116 !---------------------------------------------------
3117 ! Retrive temperature using EC geopotential height
3118 !---------------------------------------------------
3119  do i=is,ie
3120 ! Make sure FV3 top is lower than GFS; can not do extrapolation above the top at this point
3121  if ( pn1(i,1) .lt. pn0(i,1) ) then
3122  call mpp_error(fatal,'FV3 top higher than ECMWF')
3123  endif
3124 
3125  do k=1,km+1
3126  pn(k) = pn0(i,k)
3127  gz(k) = zh(i,j,k)*grav
3128  enddo
3129 !-------------------------------------------------
3130  do k=km+2, km+k2
3131  l = 2*(km+1) - k
3132  gz(k) = 2.*gz(km+1) - gz(l)
3133  pn(k) = 2.*pn(km+1) - pn(l)
3134  enddo
3135 !-------------------------------------------------
3136  gz_fv(npz+1) = atm%phis(i,j)
3137 
3138  m = 1
3139  do k=1,npz
3140 ! Searching using FV3 log(pe): pn1
3141 #ifdef USE_ISOTHERMO
3142  do l=m,km
3143  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
3144  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
3145  goto 555
3146  elseif ( pn1(i,k) .gt. pn(km+1) ) then
3147 ! Isothermal under ground; linear in log-p extra-polation
3148  gz_fv(k) = gz(km+1) + (gz_fv(npz+1)-gz(km+1))*(pn1(i,k)-pn(km+1))/(pn1(i,npz+1)-pn(km+1))
3149  goto 555
3150  endif
3151  enddo
3152 #else
3153  do l=m,km+k2-1
3154  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
3155  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
3156  goto 555
3157  endif
3158  enddo
3159 #endif
3160 555 m = l
3161  enddo
3162 
3163  do k=1,npz+1
3164  atm%peln(i,k,j) = pn1(i,k)
3165  enddo
3166 
3167 ! Compute true temperature using hydrostatic balance
3168  do k=1,npz
3169 ! qc = 1.-(Atm%q(i,j,k,liq_wat)+Atm%q(i,j,k,rainwat)+Atm%q(i,j,k,ice_wat)+Atm%q(i,j,k,snowwat))
3170 ! Atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))*qc/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*Atm%q(i,j,k,sphum)) )
3171 #ifdef MULTI_GASES
3172  atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*virq(atm%q(i,j,k,:)) )
3173 #else
3174  atm%pt(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*atm%q(i,j,k,sphum)) )
3175 #endif
3176  enddo
3177  if ( .not. atm%flagstruct%hydrostatic ) then
3178  do k=1,npz
3179  atm%delz(i,j,k) = (gz_fv(k+1) - gz_fv(k)) / grav
3180  enddo
3181  endif
3182 
3183  enddo ! i-loop
3184 
3185 !-------------------------------------------------------------
3186 ! map omega
3187 !------- ------------------------------------------------------
3188  if ( .not. atm%flagstruct%hydrostatic ) then
3189  do k=1,km
3190  do i=is,ie
3191  qp(i,k) = wc(i,j,k)
3192  enddo
3193  enddo
3194  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, atm%ptop)
3195  do k=1,npz
3196  do i=is,ie
3197  atm%w(i,j,k) = qn1(i,k)/atm%delp(i,j,k)*atm%delz(i,j,k)
3198  enddo
3199  enddo
3200  endif
3201 
3202 5000 continue
3203 
3204 ! Add some diagnostics:
3205  call p_maxmin('PS_model (mb)', atm%ps(is:ie,js:je), is, ie, js, je, 1, 0.01)
3206  call p_maxmin('PT_model', atm%pt(is:ie,js:je,1:npz), is, ie, js, je, npz, 1.)
3207  call pmaxmn('ZS_model', atm%phis(is:ie,js:je)/grav, is, ie, js, je, 1, 1., atm%gridstruct%area_64, atm%domain)
3208  call pmaxmn('ZS_EC', zh(is:ie,js:je,km+1), is, ie, js, je, 1, 1., atm%gridstruct%area_64, atm%domain)
3209  do j=js,je
3210  do i=is,ie
3211  wk(i,j) = atm%phis(i,j)/grav - zh(i,j,km+1)
3212  ! if ((wk(i,j) > 1800.).or.(wk(i,j)<-1600.)) then
3213  ! print *,' '
3214  ! print *, 'Diff = ', wk(i,j), 'Atm%phis =', Atm%phis(i,j)/grav, 'zh = ', zh(i,j,km+1)
3215  ! print *, 'lat = ', Atm%gridstruct%agrid(i,j,2)/deg2rad, 'lon = ', Atm%gridstruct%agrid(i,j,1)/deg2rad
3216  ! endif
3217  enddo
3218  enddo
3219  call pmaxmn('ZS_diff (m)', wk, is, ie, js, je, 1, 1., atm%gridstruct%area_64, atm%domain)
3220 
3221  if (.not.atm%neststruct%nested) then
3222  call prt_gb_nh_sh('IFS_IC Z500', is,ie, js,je, z500, atm%gridstruct%area_64(is:ie,js:je), atm%gridstruct%agrid_64(is:ie,js:je,2))
3223  if ( .not. atm%flagstruct%hydrostatic ) &
3224  call prt_height('fv3_IC Z500', is,ie, js,je, 3, npz, 500.e2, atm%phis, atm%delz, atm%peln, &
3225  atm%gridstruct%area_64(is:ie,js:je), atm%gridstruct%agrid_64(is:ie,js:je,2))
3226  endif
3227 
3228  do j=js,je
3229  do i=is,ie
3230  wk(i,j) = atm%ps(i,j) - psc(i,j)
3231  enddo
3232  enddo
3233  call pmaxmn('PS_diff (mb)', wk, is, ie, js, je, 1, 0.01, atm%gridstruct%area_64, atm%domain)
3234 
3235  end subroutine remap_scalar_ec
3236 
3237  subroutine remap_scalar_single(Atm, km, npz, ak0, bk0, psc, qa, zh ,iq)
3238  type(fv_atmos_type), intent(inout) :: Atm
3239  integer, intent(in):: km, npz, iq
3240  real, intent(in):: ak0(km+1), bk0(km+1)
3241  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je):: psc
3242  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: qa
3243  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km+1):: zh
3244 ! local:
3245  real, dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0
3246  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1, dp2
3247  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1
3248  real qp(atm%bd%is:atm%bd%ie,km)
3249  real wk(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je)
3250 !!! High-precision
3251  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pn1
3252  real(kind=R_GRID):: gz_fv(npz+1)
3253  real(kind=R_GRID), dimension(2*km+1):: gz, pn
3254  real(kind=R_GRID), dimension(Atm%bd%is:Atm%bd%ie,km+1):: pn0
3255  real(kind=R_GRID):: pst
3256 !!! High-precision
3257  integer i,j,k, k2, l
3258  integer :: is, ie, js, je
3259  real, allocatable:: ps_temp(:,:)
3260 
3261  is = atm%bd%is
3262  ie = atm%bd%ie
3263  js = atm%bd%js
3264  je = atm%bd%je
3265 
3266  k2 = max(10, km/2)
3267 
3268  allocate(ps_temp(is:ie,js:je))
3269 
3270  do 5000 j=js,je
3271  do k=1,km+1
3272  do i=is,ie
3273  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
3274  pn0(i,k) = log(pe0(i,k))
3275  enddo
3276  enddo
3277 
3278  do i=is,ie
3279  do k=1,km+1
3280  pn(k) = pn0(i,k)
3281  gz(k) = zh(i,j,k)*grav
3282  enddo
3283 ! Use log-p for interpolation/extrapolation
3284 ! mirror image method:
3285  do k=km+2, km+k2
3286  l = 2*(km+1) - k
3287  gz(k) = 2.*gz(km+1) - gz(l)
3288  pn(k) = 2.*pn(km+1) - pn(l)
3289  enddo
3290 
3291  do k=km+k2-1, 2, -1
3292  if( atm%phis(i,j).le.gz(k) .and. atm%phis(i,j).ge.gz(k+1) ) then
3293  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
3294  go to 123
3295  endif
3296  enddo
3297 123 ps_temp(i,j) = exp(pst)
3298  enddo ! i-loop
3299 
3300  do i=is,ie
3301  pe1(i,1) = atm%ak(1)
3302  pn1(i,1) = log(pe1(i,1))
3303  enddo
3304  do k=2,npz+1
3305  do i=is,ie
3306  pe1(i,k) = atm%ak(k) + atm%bk(k)*ps_temp(i,j)
3307  pn1(i,k) = log(pe1(i,k))
3308  enddo
3309  enddo
3310 
3311 ! * Compute delp
3312  do k=1,npz
3313  do i=is,ie
3314  dp2(i,k) = pe1(i,k+1) - pe1(i,k)
3315  enddo
3316  enddo
3317 
3318  ! map o3mr
3319  do k=1,km
3320  do i=is,ie
3321  qp(i,k) = qa(i,j,k)
3322  enddo
3323  enddo
3324  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 0, 8, atm%ptop)
3325  if ( iq==1 ) then
3326  call fillq(ie-is+1, npz, 1, qn1, dp2)
3327  else
3328  call fillz(ie-is+1, npz, 1, qn1, dp2)
3329  endif
3330 ! The HiRam step of blending model sphum with NCEP data is obsolete because nggps is always cold starting...
3331  do k=1,npz
3332  do i=is,ie
3333  atm%q(i,j,k,iq) = qn1(i,k)
3334  enddo
3335  enddo
3336 
3337 5000 continue
3338  call p_maxmin('o3mr remap', atm%q(is:ie,js:je,1:npz,iq), is, ie, js, je, npz, 1.)
3339 
3340  deallocate(ps_temp)
3341 
3342  end subroutine remap_scalar_single
3343 
3344 
3345  subroutine mp_auto_conversion(ql, qr, qi, qs)
3346  real, intent(inout):: ql, qr, qi, qs
3347  real, parameter:: qi0_max = 2.0e-3
3348  real, parameter:: ql0_max = 2.5e-3
3349 
3350 ! Convert excess cloud water into rain:
3351  if ( ql > ql0_max ) then
3352  qr = ql - ql0_max
3353  ql = ql0_max
3354  endif
3355 ! Convert excess cloud ice into snow:
3356  if ( qi > qi0_max ) then
3357  qs = qi - qi0_max
3358  qi = qi0_max
3359  endif
3360 
3361  end subroutine mp_auto_conversion
3362 
3363 
3364  subroutine remap_dwinds(km, npz, ak0, bk0, psc, ud, vd, Atm)
3365  type(fv_atmos_type), intent(inout) :: Atm
3366  integer, intent(in):: km, npz
3367  real, intent(in):: ak0(km+1), bk0(km+1)
3368  real, intent(in):: psc(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je)
3369  real, intent(in):: ud(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je+1,km)
3370  real, intent(in):: vd(atm%bd%is:atm%bd%ie+1,atm%bd%js:atm%bd%je,km)
3371 ! local:
3372  real, dimension(Atm%bd%isd:Atm%bd%ied,Atm%bd%jsd:Atm%bd%jed):: psd
3373  real, dimension(Atm%bd%is:Atm%bd%ie+1, km+1):: pe0
3374  real, dimension(Atm%bd%is:Atm%bd%ie+1,npz+1):: pe1
3375  real, dimension(Atm%bd%is:Atm%bd%ie+1,npz):: qn1
3376  integer i,j,k
3377  integer :: is, ie, js, je
3378  integer :: isd, ied, jsd, jed
3379 
3380  is = atm%bd%is
3381  ie = atm%bd%ie
3382  js = atm%bd%js
3383  je = atm%bd%je
3384  isd = atm%bd%isd
3385  ied = atm%bd%ied
3386  jsd = atm%bd%jsd
3387  jed = atm%bd%jed
3388 
3389  if (atm%neststruct%nested .or. atm%flagstruct%regional) then
3390  do j=jsd,jed
3391  do i=isd,ied
3392  psd(i,j) = atm%ps(i,j)
3393  enddo
3394  enddo
3395  else
3396  do j=js,je
3397  do i=is,ie
3398  psd(i,j) = psc(i,j)
3399  enddo
3400  enddo
3401  endif
3402  call mpp_update_domains( psd, atm%domain, complete=.false. )
3403  call mpp_update_domains( atm%ps, atm%domain, complete=.true. )
3404 
3405 !$OMP parallel do default(none) shared(is,ie,js,je,npz,km,ak0,bk0,Atm,psc,psd,ud,vd) &
3406 !$OMP private(pe1,pe0,qn1)
3407  do 5000 j=js,je+1
3408 !------
3409 ! map u
3410 !------
3411  do k=1,km+1
3412  do i=is,ie
3413  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j-1)+psd(i,j))
3414  enddo
3415  enddo
3416  do k=1,npz+1
3417  do i=is,ie
3418  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(atm%ps(i,j-1)+atm%ps(i,j))
3419  enddo
3420  enddo
3421  call mappm(km, pe0(is:ie,1:km+1), ud(is:ie,j,1:km), npz, pe1(is:ie,1:npz+1), &
3422  qn1(is:ie,1:npz), is,ie, -1, 8, atm%ptop)
3423  do k=1,npz
3424  do i=is,ie
3425  atm%u(i,j,k) = qn1(i,k)
3426  enddo
3427  enddo
3428 !------
3429 ! map v
3430 !------
3431  if ( j/=(je+1) ) then
3432 
3433  do k=1,km+1
3434  do i=is,ie+1
3435  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i-1,j)+psd(i,j))
3436  enddo
3437  enddo
3438  do k=1,npz+1
3439  do i=is,ie+1
3440  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(atm%ps(i-1,j)+atm%ps(i,j))
3441  enddo
3442  enddo
3443  call mappm(km, pe0(is:ie+1,1:km+1), vd(is:ie+1,j,1:km), npz, pe1(is:ie+1,1:npz+1), &
3444  qn1(is:ie+1,1:npz), is,ie+1, -1, 8, atm%ptop)
3445  do k=1,npz
3446  do i=is,ie+1
3447  atm%v(i,j,k) = qn1(i,k)
3448  enddo
3449  enddo
3450 
3451  endif
3452 
3453 5000 continue
3454 
3455  if (is_master()) write(*,*) 'done remap_dwinds'
3456 
3457  end subroutine remap_dwinds
3458 
3459 
3460  subroutine remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
3461  type(fv_atmos_type), intent(inout) :: Atm
3462  integer, intent(in):: im, jm, km, npz
3463  real, intent(in):: ak0(km+1), bk0(km+1)
3464  real, intent(in):: psc(atm%bd%is:atm%bd%ie,atm%bd%js:atm%bd%je)
3465  real, intent(in), dimension(Atm%bd%is:Atm%bd%ie,Atm%bd%js:Atm%bd%je,km):: ua, va
3466 ! local:
3467  real, dimension(Atm%bd%isd:Atm%bd%ied,Atm%bd%jsd:Atm%bd%jed,npz):: ut, vt ! winds
3468  real, dimension(Atm%bd%is:Atm%bd%ie, km+1):: pe0
3469  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1
3470  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1
3471  integer i,j,k
3472 
3473  integer :: is, ie, js, je
3474  integer :: isd, ied, jsd, jed
3475 
3476  is = atm%bd%is
3477  ie = atm%bd%ie
3478  js = atm%bd%js
3479  je = atm%bd%je
3480  isd = atm%bd%isd
3481  ied = atm%bd%ied
3482  jsd = atm%bd%jsd
3483  jed = atm%bd%jed
3484 
3485  do 5000 j=js,je
3486 
3487  do k=1,km+1
3488  do i=is,ie
3489  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
3490  enddo
3491  enddo
3492 
3493  do k=1,npz+1
3494  do i=is,ie
3495  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
3496  enddo
3497  enddo
3498 
3499 !------
3500 ! map u
3501 !------
3502  call mappm(km, pe0, ua(is:ie,j,1:km), npz, pe1, qn1, is,ie, -1, 8, atm%ptop)
3503  do k=1,npz
3504  do i=is,ie
3505  ut(i,j,k) = qn1(i,k)
3506  enddo
3507  enddo
3508 !------
3509 ! map v
3510 !------
3511  call mappm(km, pe0, va(is:ie,j,1:km), npz, pe1, qn1, is,ie, -1, 8, atm%ptop)
3512  do k=1,npz
3513  do i=is,ie
3514  vt(i,j,k) = qn1(i,k)
3515  enddo
3516  enddo
3517 
3518 5000 continue
3519 
3520  call prt_maxmin('UT', ut, is, ie, js, je, ng, npz, 1.)
3521  call prt_maxmin('VT', vt, is, ie, js, je, ng, npz, 1.)
3522  call prt_maxmin('UA_top',ut(:,:,1), is, ie, js, je, ng, 1, 1.)
3523 
3524 !----------------------------------------------
3525 ! winds: lat-lon ON A to Cubed-D transformation:
3526 !----------------------------------------------
3527  call cubed_a2d(atm%npx, atm%npy, npz, ut, vt, atm%u, atm%v, atm%gridstruct, atm%domain, atm%bd )
3528 
3529  if (is_master()) write(*,*) 'done remap_winds'
3530 
3531  end subroutine remap_winds
3532 
3533 
3534  subroutine remap_xyz( im, jbeg, jend, jm, km, npz, nq, ncnst, lon, lat, ak0, bk0, ps0, gz0, &
3535  ua, va, ta, qa, Atm )
3537  type(fv_atmos_type), intent(inout), target :: Atm
3538  integer, intent(in):: im, jm, km, npz, nq, ncnst
3539  integer, intent(in):: jbeg, jend
3540  real, intent(in):: lon(im), lat(jm), ak0(km+1), bk0(km+1)
3541  real, intent(in):: gz0(im,jbeg:jend), ps0(im,jbeg:jend)
3542  real, intent(in), dimension(im,jbeg:jend,km):: ua, va, ta
3543  real, intent(in), dimension(im,jbeg:jend,km,ncnst):: qa
3544 
3545  real, pointer, dimension(:,:,:) :: agrid
3546 
3547 ! local:
3548  real, dimension(Atm%bd%isd:Atm%bd%ied,Atm%bd%jsd:Atm%bd%jed,npz):: ut, vt ! winds
3549  real, dimension(Atm%bd%is:Atm%bd%ie,km):: up, vp, tp
3550  real, dimension(Atm%bd%is:Atm%bd%ie,km+1):: pe0, pn0
3551  real pt0(km), gz(km+1), pk0(km+1)
3552  real qp(atm%bd%is:atm%bd%ie,km,ncnst)
3553  real, dimension(Atm%bd%is:Atm%bd%ie,npz):: qn1
3554  real, dimension(Atm%bd%is:Atm%bd%ie,npz+1):: pe1, pn1
3555  real :: rdlon(im)
3556  real :: rdlat(jm)
3557  real:: a1, b1, c1, c2, c3, c4
3558  real:: gzc, psc, pst
3559 #ifdef MULTI_GASES
3560  real:: kappax, pkx
3561 #endif
3562  integer i,j,k, i1, i2, jc, i0, j0, iq
3563 ! integer sphum, liq_wat, ice_wat, cld_amt
3564  integer sphum
3565  integer :: is, ie, js, je
3566  integer :: isd, ied, jsd, jed
3567 
3568  is = atm%bd%is
3569  ie = atm%bd%ie
3570  js = atm%bd%js
3571  je = atm%bd%je
3572  isd = atm%bd%isd
3573  ied = atm%bd%ied
3574  jsd = atm%bd%jsd
3575  jed = atm%bd%jed
3576 
3577  !!NOTE: Only Atm is used in this routine.
3578  agrid => atm%gridstruct%agrid
3579 
3580  sphum = get_tracer_index(model_atmos, 'sphum')
3581 
3582  if ( sphum/=1 ) then
3583  call mpp_error(fatal,'SPHUM must be 1st tracer')
3584  endif
3585 
3586  pk0(1) = ak0(1)**kappa
3587 
3588  do i=1,im-1
3589  rdlon(i) = 1. / (lon(i+1) - lon(i))
3590  enddo
3591  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
3592 
3593  do j=1,jm-1
3594  rdlat(j) = 1. / (lat(j+1) - lat(j))
3595  enddo
3596 
3597 ! * Interpolate to cubed sphere cell center
3598  do 5000 j=js,je
3599 
3600  do i=is,ie
3601  pe0(i,1) = ak0(1)
3602  pn0(i,1) = log(ak0(1))
3603  enddo
3604 
3605 
3606  do i=is,ie
3607 
3608  if ( agrid(i,j,1)>lon(im) ) then
3609  i1 = im; i2 = 1
3610  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
3611  elseif ( agrid(i,j,1)<lon(1) ) then
3612  i1 = im; i2 = 1
3613  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
3614  else
3615  do i0=1,im-1
3616  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
3617  i1 = i0; i2 = i0+1
3618  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
3619  go to 111
3620  endif
3621  enddo
3622  endif
3623 
3624 111 continue
3625 
3626  if ( agrid(i,j,2)<lat(1) ) then
3627  jc = 1
3628  b1 = 0.
3629  elseif ( agrid(i,j,2)>lat(jm) ) then
3630  jc = jm-1
3631  b1 = 1.
3632  else
3633  do j0=1,jm-1
3634  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
3635  jc = j0
3636  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
3637  go to 222
3638  endif
3639  enddo
3640  endif
3641 222 continue
3642 
3643 #ifndef DEBUG_REMAP
3644  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
3645  write(*,*) i,j,a1, b1
3646  endif
3647 #endif
3648  c1 = (1.-a1) * (1.-b1)
3649  c2 = a1 * (1.-b1)
3650  c3 = a1 * b1
3651  c4 = (1.-a1) * b1
3652 
3653 ! Interpolated surface pressure
3654  psc = c1*ps0(i1,jc ) + c2*ps0(i2,jc ) + &
3655  c3*ps0(i2,jc+1) + c4*ps0(i1,jc+1)
3656 
3657 ! Interpolated surface geopotential
3658  gzc = c1*gz0(i1,jc ) + c2*gz0(i2,jc ) + &
3659  c3*gz0(i2,jc+1) + c4*gz0(i1,jc+1)
3660 
3661 ! 3D fields:
3662  do iq=1,ncnst
3663 ! if ( iq==sphum .or. iq==liq_wat .or. iq==ice_wat .or. iq==cld_amt ) then
3664  do k=1,km
3665  qp(i,k,iq) = c1*qa(i1,jc, k,iq) + c2*qa(i2,jc, k,iq) + &
3666  c3*qa(i2,jc+1,k,iq) + c4*qa(i1,jc+1,k,iq)
3667  enddo
3668 ! endif
3669  enddo
3670 
3671  do k=1,km
3672  up(i,k) = c1*ua(i1,jc, k) + c2*ua(i2,jc, k) + &
3673  c3*ua(i2,jc+1,k) + c4*ua(i1,jc+1,k)
3674  vp(i,k) = c1*va(i1,jc, k) + c2*va(i2,jc, k) + &
3675  c3*va(i2,jc+1,k) + c4*va(i1,jc+1,k)
3676  tp(i,k) = c1*ta(i1,jc, k) + c2*ta(i2,jc, k) + &
3677  c3*ta(i2,jc+1,k) + c4*ta(i1,jc+1,k)
3678 ! Virtual effect:
3679 #ifdef MULTI_GASES
3680  tp(i,k) = tp(i,k)*virq(qp(i,k,:))
3681 #else
3682  tp(i,k) = tp(i,k)*(1.+zvir*qp(i,k,sphum))
3683 #endif
3684  enddo
3685 ! Tracers:
3686 
3687  do k=2,km+1
3688  pe0(i,k) = ak0(k) + bk0(k)*psc
3689  pn0(i,k) = log(pe0(i,k))
3690  pk0(k) = pe0(i,k)**kappa
3691  enddo
3692 
3693 #ifdef USE_DATA_ZS
3694  atm% ps(i,j) = psc
3695  atm%phis(i,j) = gzc
3696 #else
3697 
3698 ! * Adjust interpolated ps to model terrain
3699  gz(km+1) = gzc
3700  do k=km,1,-1
3701  gz(k) = gz(k+1) + rdgas*tp(i,k)*(pn0(i,k+1)-pn0(i,k))
3702  enddo
3703 ! Only lowest layer potential temp is needed
3704 #ifdef MULTI_GASES
3705  kappax = virqd(qp(i,km,:))/vicpqd(qp(i,km,:))
3706  pkx = (pk0(km+1)-pk0(km))*(kappa*(pn0(i,km+1)-pn0(i,km)))
3707  pkx = exp( kappax*log(pkx) )
3708  pt0(km) = tp(i,km)/pkx
3709 #else
3710  pt0(km) = tp(i,km)/(pk0(km+1)-pk0(km))*(kappa*(pn0(i,km+1)-pn0(i,km)))
3711 #endif
3712  if( atm%phis(i,j)>gzc ) then
3713  do k=km,1,-1
3714  if( atm%phis(i,j) < gz(k) .and. &
3715  atm%phis(i,j) >= gz(k+1) ) then
3716  pst = pk0(k) + (pk0(k+1)-pk0(k))*(gz(k)-atm%phis(i,j))/(gz(k)-gz(k+1))
3717  go to 123
3718  endif
3719  enddo
3720  else
3721 ! Extrapolation into the ground
3722 #ifdef MULTI_GASES
3723  pst = pk0(km+1) + (gzc-atm%phis(i,j))/(cp_air*pt0(km)*pkx)
3724 #else
3725  pst = pk0(km+1) + (gzc-atm%phis(i,j))/(cp_air*pt0(km))
3726 #endif
3727  endif
3728 
3729 #ifdef MULTI_GASES
3730 123 atm%ps(i,j) = pst**(1./(kappa*kappax))
3731 #else
3732 123 atm%ps(i,j) = pst**(1./kappa)
3733 #endif
3734 #endif
3735  enddo !i-loop
3736 
3737 
3738 ! * Compute delp from ps
3739  do i=is,ie
3740  pe1(i,1) = atm%ak(1)
3741  pn1(i,1) = log(pe1(i,1))
3742  enddo
3743  do k=2,npz+1
3744  do i=is,ie
3745  pe1(i,k) = atm%ak(k) + atm%bk(k)*atm%ps(i,j)
3746  pn1(i,k) = log(pe1(i,k))
3747  enddo
3748  enddo
3749 
3750  do k=1,npz
3751  do i=is,ie
3752  atm%delp(i,j,k) = pe1(i,k+1) - pe1(i,k)
3753  enddo
3754  enddo
3755 
3756 ! Use kord=9 for winds; kord=11 for tracers
3757 !------
3758 ! map u
3759 !------
3760  call mappm(km, pe0, up, npz, pe1, qn1, is,ie, -1, 9, atm%ptop)
3761  do k=1,npz
3762  do i=is,ie
3763  ut(i,j,k) = qn1(i,k)
3764  enddo
3765  enddo
3766 !------
3767 ! map v
3768 !------
3769  call mappm(km, pe0, vp, npz, pe1, qn1, is,ie, -1, 9, atm%ptop)
3770  do k=1,npz
3771  do i=is,ie
3772  vt(i,j,k) = qn1(i,k)
3773  enddo
3774  enddo
3775 
3776 !---------------
3777 ! map tracers
3778 !----------------
3779  do iq=1,ncnst
3780 ! Note: AM2 physics tracers only
3781 ! if ( iq==sphum .or. iq==liq_wat .or. iq==ice_wat .or. iq==cld_amt ) then
3782  call mappm(km, pe0, qp(is,1,iq), npz, pe1, qn1, is,ie, 0, 11, atm%ptop)
3783  do k=1,npz
3784  do i=is,ie
3785  atm%q(i,j,k,iq) = qn1(i,k)
3786  enddo
3787  enddo
3788 ! endif
3789  enddo
3790 
3791 !-------------------------------------------------------------
3792 ! map virtual temperature using geopotential conserving scheme.
3793 !-------------------------------------------------------------
3794  call mappm(km, pn0, tp, npz, pn1, qn1, is,ie, 1, 9, atm%ptop)
3795  do k=1,npz
3796  do i=is,ie
3797 #ifdef MULTI_GASES
3798  atm%pt(i,j,k) = qn1(i,k)/virq(atm%q(i,j,k,:))
3799 #else
3800  atm%pt(i,j,k) = qn1(i,k)/(1.+zvir*atm%q(i,j,k,sphum))
3801 #endif
3802  enddo
3803  enddo
3804 
3805 5000 continue
3806 
3807  call prt_maxmin('PS_model', atm%ps, is, ie, js, je, ng, 1, 0.01)
3808  call prt_maxmin('UT', ut, is, ie, js, je, ng, npz, 1.)
3809  call prt_maxmin('VT', vt, is, ie, js, je, ng, npz, 1.)
3810 
3811 !----------------------------------------------
3812 ! winds: lat-lon ON A to Cubed-D transformation:
3813 !----------------------------------------------
3814  call cubed_a2d(atm%npx, atm%npy, npz, ut, vt, atm%u, atm%v, atm%gridstruct, atm%domain, atm%bd )
3815 
3816  if (is_master()) write(*,*) 'done remap_xyz'
3817 
3818  end subroutine remap_xyz
3819 
3821  subroutine cubed_a2d( npx, npy, npz, ua, va, u, v, gridstruct, fv_domain, bd )
3822  use mpp_domains_mod, only: mpp_update_domains
3823 
3824  type(fv_grid_bounds_type), intent(IN) :: bd
3825  integer, intent(in):: npx, npy, npz
3826  real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va
3827  real, intent(out):: u(bd%isd:bd%ied, bd%jsd:bd%jed+1,npz)
3828  real, intent(out):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
3829  type(fv_grid_type), intent(IN), target :: gridstruct
3830  type(domain2d), intent(INOUT) :: fv_domain
3831 ! local:
3832  real v3(3,bd%is-1:bd%ie+1,bd%js-1:bd%je+1)
3833  real ue(3,bd%is-1:bd%ie+1,bd%js:bd%je+1)
3834  real ve(3,bd%is:bd%ie+1,bd%js-1:bd%je+1)
3835  real, dimension(bd%is:bd%ie):: ut1, ut2, ut3
3836  real, dimension(bd%js:bd%je):: vt1, vt2, vt3
3837  integer i, j, k, im2, jm2
3838 
3839  real(kind=R_GRID), pointer, dimension(:,:,:) :: vlon, vlat
3840  real(kind=R_GRID), pointer, dimension(:) :: edge_vect_w, edge_vect_e, edge_vect_s, edge_vect_n
3841  real(kind=R_GRID), pointer, dimension(:,:,:,:) :: ew, es
3842 
3843  integer :: is, ie, js, je
3844  integer :: isd, ied, jsd, jed
3845 
3846  is = bd%is
3847  ie = bd%ie
3848  js = bd%js
3849  je = bd%je
3850  isd = bd%isd
3851  ied = bd%ied
3852  jsd = bd%jsd
3853  jed = bd%jed
3854 
3855  vlon => gridstruct%vlon
3856  vlat => gridstruct%vlat
3857 
3858  edge_vect_w => gridstruct%edge_vect_w
3859  edge_vect_e => gridstruct%edge_vect_e
3860  edge_vect_s => gridstruct%edge_vect_s
3861  edge_vect_n => gridstruct%edge_vect_n
3862 
3863  ew => gridstruct%ew
3864  es => gridstruct%es
3865 
3866  call mpp_update_domains(ua, fv_domain, complete=.false.)
3867  call mpp_update_domains(va, fv_domain, complete=.true.)
3868 
3869  im2 = (npx-1)/2
3870  jm2 = (npy-1)/2
3871 
3872  do k=1, npz
3873 ! Compute 3D wind on A grid
3874  do j=js-1,je+1
3875  do i=is-1,ie+1
3876  v3(1,i,j) = ua(i,j,k)*vlon(i,j,1) + va(i,j,k)*vlat(i,j,1)
3877  v3(2,i,j) = ua(i,j,k)*vlon(i,j,2) + va(i,j,k)*vlat(i,j,2)
3878  v3(3,i,j) = ua(i,j,k)*vlon(i,j,3) + va(i,j,k)*vlat(i,j,3)
3879  enddo
3880  enddo
3881 
3882 ! A --> D
3883 ! Interpolate to cell edges
3884  do j=js,je+1
3885  do i=is-1,ie+1
3886  ue(1,i,j) = 0.5*(v3(1,i,j-1) + v3(1,i,j))
3887  ue(2,i,j) = 0.5*(v3(2,i,j-1) + v3(2,i,j))
3888  ue(3,i,j) = 0.5*(v3(3,i,j-1) + v3(3,i,j))
3889  enddo
3890  enddo
3891 
3892  do j=js-1,je+1
3893  do i=is,ie+1
3894  ve(1,i,j) = 0.5*(v3(1,i-1,j) + v3(1,i,j))
3895  ve(2,i,j) = 0.5*(v3(2,i-1,j) + v3(2,i,j))
3896  ve(3,i,j) = 0.5*(v3(3,i-1,j) + v3(3,i,j))
3897  enddo
3898  enddo
3899 
3900 ! --- E_W edges (for v-wind):
3901  if (.not. gridstruct%nested) then
3902  if ( is==1) then
3903  i = 1
3904  do j=js,je
3905  if ( j>jm2 ) then
3906  vt1(j) = edge_vect_w(j)*ve(1,i,j-1)+(1.-edge_vect_w(j))*ve(1,i,j)
3907  vt2(j) = edge_vect_w(j)*ve(2,i,j-1)+(1.-edge_vect_w(j))*ve(2,i,j)
3908  vt3(j) = edge_vect_w(j)*ve(3,i,j-1)+(1.-edge_vect_w(j))*ve(3,i,j)
3909  else
3910  vt1(j) = edge_vect_w(j)*ve(1,i,j+1)+(1.-edge_vect_w(j))*ve(1,i,j)
3911  vt2(j) = edge_vect_w(j)*ve(2,i,j+1)+(1.-edge_vect_w(j))*ve(2,i,j)
3912  vt3(j) = edge_vect_w(j)*ve(3,i,j+1)+(1.-edge_vect_w(j))*ve(3,i,j)
3913  endif
3914  enddo
3915  do j=js,je
3916  ve(1,i,j) = vt1(j)
3917  ve(2,i,j) = vt2(j)
3918  ve(3,i,j) = vt3(j)
3919  enddo
3920  endif
3921 
3922  if ( (ie+1)==npx ) then
3923  i = npx
3924  do j=js,je
3925  if ( j>jm2 ) then
3926  vt1(j) = edge_vect_e(j)*ve(1,i,j-1)+(1.-edge_vect_e(j))*ve(1,i,j)
3927  vt2(j) = edge_vect_e(j)*ve(2,i,j-1)+(1.-edge_vect_e(j))*ve(2,i,j)
3928  vt3(j) = edge_vect_e(j)*ve(3,i,j-1)+(1.-edge_vect_e(j))*ve(3,i,j)
3929  else
3930  vt1(j) = edge_vect_e(j)*ve(1,i,j+1)+(1.-edge_vect_e(j))*ve(1,i,j)
3931  vt2(j) = edge_vect_e(j)*ve(2,i,j+1)+(1.-edge_vect_e(j))*ve(2,i,j)
3932  vt3(j) = edge_vect_e(j)*ve(3,i,j+1)+(1.-edge_vect_e(j))*ve(3,i,j)
3933  endif
3934  enddo
3935  do j=js,je
3936  ve(1,i,j) = vt1(j)
3937  ve(2,i,j) = vt2(j)
3938  ve(3,i,j) = vt3(j)
3939  enddo
3940  endif
3941 
3942 ! N-S edges (for u-wind):
3943  if ( js==1 ) then
3944  j = 1
3945  do i=is,ie
3946  if ( i>im2 ) then
3947  ut1(i) = edge_vect_s(i)*ue(1,i-1,j)+(1.-edge_vect_s(i))*ue(1,i,j)
3948  ut2(i) = edge_vect_s(i)*ue(2,i-1,j)+(1.-edge_vect_s(i))*ue(2,i,j)
3949  ut3(i) = edge_vect_s(i)*ue(3,i-1,j)+(1.-edge_vect_s(i))*ue(3,i,j)
3950  else
3951  ut1(i) = edge_vect_s(i)*ue(1,i+1,j)+(1.-edge_vect_s(i))*ue(1,i,j)
3952  ut2(i) = edge_vect_s(i)*ue(2,i+1,j)+(1.-edge_vect_s(i))*ue(2,i,j)
3953  ut3(i) = edge_vect_s(i)*ue(3,i+1,j)+(1.-edge_vect_s(i))*ue(3,i,j)
3954  endif
3955  enddo
3956  do i=is,ie
3957  ue(1,i,j) = ut1(i)
3958  ue(2,i,j) = ut2(i)
3959  ue(3,i,j) = ut3(i)
3960  enddo
3961  endif
3962 
3963  if ( (je+1)==npy ) then
3964  j = npy
3965  do i=is,ie
3966  if ( i>im2 ) then
3967  ut1(i) = edge_vect_n(i)*ue(1,i-1,j)+(1.-edge_vect_n(i))*ue(1,i,j)
3968  ut2(i) = edge_vect_n(i)*ue(2,i-1,j)+(1.-edge_vect_n(i))*ue(2,i,j)
3969  ut3(i) = edge_vect_n(i)*ue(3,i-1,j)+(1.-edge_vect_n(i))*ue(3,i,j)
3970  else
3971  ut1(i) = edge_vect_n(i)*ue(1,i+1,j)+(1.-edge_vect_n(i))*ue(1,i,j)
3972  ut2(i) = edge_vect_n(i)*ue(2,i+1,j)+(1.-edge_vect_n(i))*ue(2,i,j)
3973  ut3(i) = edge_vect_n(i)*ue(3,i+1,j)+(1.-edge_vect_n(i))*ue(3,i,j)
3974  endif
3975  enddo
3976  do i=is,ie
3977  ue(1,i,j) = ut1(i)
3978  ue(2,i,j) = ut2(i)
3979  ue(3,i,j) = ut3(i)
3980  enddo
3981  endif
3982 
3983  endif ! .not. nested
3984 
3985  do j=js,je+1
3986  do i=is,ie
3987  u(i,j,k) = ue(1,i,j)*es(1,i,j,1) + &
3988  ue(2,i,j)*es(2,i,j,1) + &
3989  ue(3,i,j)*es(3,i,j,1)
3990  enddo
3991  enddo
3992  do j=js,je
3993  do i=is,ie+1
3994  v(i,j,k) = ve(1,i,j)*ew(1,i,j,2) + &
3995  ve(2,i,j)*ew(2,i,j,2) + &
3996  ve(3,i,j)*ew(3,i,j,2)
3997  enddo
3998  enddo
3999 
4000  enddo ! k-loop
4001 
4002  end subroutine cubed_a2d
4003 
4004 
4005  subroutine d2a3d(u, v, ua, va, im, jm, km, lon)
4006  integer, intent(in):: im, jm, km ! Dimensions
4007  real, intent(in ) :: lon(im)
4008  real, intent(in ), dimension(im,jm,km):: u, v
4009  real, intent(out), dimension(im,jm,km):: ua, va
4010 ! local
4011  real :: coslon(im),sinlon(im) ! Sine and cosine in longitude
4012  integer i, j, k
4013  integer imh
4014  real un, vn, us, vs
4015 
4016  integer :: ks, ke
4017 
4018  imh = im/2
4019 
4020  do i=1,im
4021  sinlon(i) = sin(lon(i))
4022  coslon(i) = cos(lon(i))
4023  enddo
4024 
4025  do k=1,km
4026  do j=2,jm-1
4027  do i=1,im
4028  ua(i,j,k) = 0.5*(u(i,j,k) + u(i,j+1,k))
4029  enddo
4030  enddo
4031 
4032  do j=2,jm-1
4033  do i=1,im-1
4034  va(i,j,k) = 0.5*(v(i,j,k) + v(i+1,j,k))
4035  enddo
4036  va(im,j,k) = 0.5*(v(im,j,k) + v(1,j,k))
4037  enddo
4038 
4039 ! Projection at SP
4040  us = 0.
4041  vs = 0.
4042  do i=1,imh
4043  us = us + (ua(i+imh,2,k)-ua(i,2,k))*sinlon(i) &
4044  + (va(i,2,k)-va(i+imh,2,k))*coslon(i)
4045  vs = vs + (ua(i+imh,2,k)-ua(i,2,k))*coslon(i) &
4046  + (va(i+imh,2,k)-va(i,2,k))*sinlon(i)
4047  enddo
4048  us = us/im
4049  vs = vs/im
4050  do i=1,imh
4051  ua(i,1,k) = -us*sinlon(i) - vs*coslon(i)
4052  va(i,1,k) = us*coslon(i) - vs*sinlon(i)
4053  ua(i+imh,1,k) = -ua(i,1,k)
4054  va(i+imh,1,k) = -va(i,1,k)
4055  enddo
4056 
4057 ! Projection at NP
4058  un = 0.
4059  vn = 0.
4060  do i=1,imh
4061  un = un + (ua(i+imh,jm-1,k)-ua(i,jm-1,k))*sinlon(i) &
4062  + (va(i+imh,jm-1,k)-va(i,jm-1,k))*coslon(i)
4063  vn = vn + (ua(i,jm-1,k)-ua(i+imh,jm-1,k))*coslon(i) &
4064  + (va(i+imh,jm-1,k)-va(i,jm-1,k))*sinlon(i)
4065  enddo
4066 
4067  un = un/im
4068  vn = vn/im
4069  do i=1,imh
4070  ua(i,jm,k) = -un*sinlon(i) + vn*coslon(i)
4071  va(i,jm,k) = -un*coslon(i) - vn*sinlon(i)
4072  ua(i+imh,jm,k) = -ua(i,jm,k)
4073  va(i+imh,jm,k) = -va(i,jm,k)
4074  enddo
4075  enddo
4076 
4077  end subroutine d2a3d
4078 
4079 
4080  subroutine pmaxmin( qname, a, im, jm, fac )
4082  integer, intent(in):: im, jm
4083  character(len=*) :: qname
4084  integer i, j
4085  real a(im,jm)
4086 
4087  real qmin(jm), qmax(jm)
4088  real pmax, pmin
4089  real fac ! multiplication factor
4090 
4091  do j=1,jm
4092  pmax = a(1,j)
4093  pmin = a(1,j)
4094  do i=2,im
4095  pmax = max(pmax, a(i,j))
4096  pmin = min(pmin, a(i,j))
4097  enddo
4098  qmax(j) = pmax
4099  qmin(j) = pmin
4100  enddo
4101 !
4102 ! Now find max/min of amax/amin
4103 !
4104  pmax = qmax(1)
4105  pmin = qmin(1)
4106  do j=2,jm
4107  pmax = max(pmax, qmax(j))
4108  pmin = min(pmin, qmin(j))
4109  enddo
4110 
4111  write(*,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac
4112 
4113  end subroutine pmaxmin
4114 
4115 subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
4116  character(len=*), intent(in):: qname
4117  integer, intent(in):: is, ie, js, je
4118  integer, intent(in):: km
4119  real, intent(in):: q(is:ie, js:je, km)
4120  real, intent(in):: fac
4121  real(kind=R_GRID), intent(IN):: area(is-3:ie+3, js-3:je+3)
4122  type(domain2d), intent(INOUT) :: domain
4123 !---local variables
4124  real qmin, qmax, gmean
4125  integer i,j,k
4126 
4127  qmin = q(is,js,1)
4128  qmax = qmin
4129  gmean = 0.
4130 
4131  do k=1,km
4132  do j=js,je
4133  do i=is,ie
4134  if( q(i,j,k) < qmin ) then
4135  qmin = q(i,j,k)
4136  elseif( q(i,j,k) > qmax ) then
4137  qmax = q(i,j,k)
4138  endif
4139  enddo
4140  enddo
4141  enddo
4142 
4143  call mp_reduce_min(qmin)
4144  call mp_reduce_max(qmax)
4145 
4146  gmean = g_sum(domain, q(is,js,km), is, ie, js, je, 3, area, 1, reproduce=.true.)
4147  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac, gmean*fac
4148 
4149  end subroutine pmaxmn
4150 
4151  subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
4152  character(len=*), intent(in):: qname
4153  integer, intent(in):: is, ie, js, je, km
4154  real, intent(in):: q(is:ie, js:je, km)
4155  real, intent(in):: fac
4156  real qmin, qmax
4157  integer i,j,k
4158 
4159  qmin = q(is,js,1)
4160  qmax = qmin
4161  do k=1,km
4162  do j=js,je
4163  do i=is,ie
4164  if( q(i,j,k) < qmin ) then
4165  qmin = q(i,j,k)
4166  elseif( q(i,j,k) > qmax ) then
4167  qmax = q(i,j,k)
4168  endif
4169  enddo
4170  enddo
4171  enddo
4172  call mp_reduce_min(qmin)
4173  call mp_reduce_max(qmax)
4174  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac
4175 
4176  end subroutine p_maxmin
4177 
4178  subroutine fillq(im, km, nq, q, dp)
4179  integer, intent(in):: im
4180  integer, intent(in):: km
4181  integer, intent(in):: nq
4182  real , intent(in):: dp(im,km)
4183  real , intent(inout) :: q(im,km,nq)
4184 ! !LOCAL VARIABLES:
4185  integer i, k, ic, k1
4186 
4187  do ic=1,nq
4188 ! Bottom up:
4189  do k=km,2,-1
4190  k1 = k-1
4191  do i=1,im
4192  if( q(i,k,ic) < 0. ) then
4193  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
4194  q(i,k ,ic) = 0.
4195  endif
4196  enddo
4197  enddo
4198 ! Top down:
4199  do k=1,km-1
4200  k1 = k+1
4201  do i=1,im
4202  if( q(i,k,ic) < 0. ) then
4203  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
4204  q(i,k ,ic) = 0.
4205  endif
4206  enddo
4207  enddo
4208 
4209  enddo
4210 
4211  end subroutine fillq
4212 
4213  subroutine compute_zh(im, jm, levp, ak0, bk0, ps, zs, t, q, nq, zh )
4214  implicit none
4215  integer, intent(in):: levp, im,jm, nq
4216  real, intent(in), dimension(levp+1):: ak0, bk0
4217  real(kind=4), intent(in), dimension(im,jm):: ps, zs
4218  real(kind=4), intent(in), dimension(im,jm,levp):: t
4219  real(kind=4), intent(in), dimension(im,jm,levp,nq):: q
4220  real(kind=4), intent(out), dimension(im,jm,levp+1):: zh
4221  ! Local:
4222  real, dimension(im,levp+1):: pe0, pn0
4223 ! real:: qc
4224  integer:: i,j,k
4225 
4226 !$OMP parallel do default(none) shared(im,jm,levp,ak0,bk0,zs,ps,t,q,zh) &
4227 !$OMP private(pe0,pn0)
4228  do j = 1, jm
4229 
4230  do i=1, im
4231  pe0(i,1) = ak0(1)
4232  pn0(i,1) = log(pe0(i,1))
4233  zh(i,j,levp+1) = zs(i,j)
4234  enddo
4235 
4236  do k=2,levp+1
4237  do i=1,im
4238  pe0(i,k) = ak0(k) + bk0(k)*ps(i,j)
4239  pn0(i,k) = log(pe0(i,k))
4240  enddo
4241  enddo
4242 
4243  do k = levp, 1, -1
4244  do i = 1, im
4245 ! qc = 1.-(q(i,j,k,2)+q(i,j,k,3)+q(i,j,k,4)+q(i,j,k,5))
4246  zh(i,j,k) = zh(i,j,k+1)+(t(i,j,k)*(1.+zvir*q(i,j,k,1))*(pn0(i,k+1)-pn0(i,k)))*(rdgas/grav)
4247  enddo
4248  enddo
4249  enddo
4250 
4251  !if(is_master()) call pmaxmin( 'zh levp+1', zh(:,:,levp+1), im, jm, 1.)
4252 
4253  end subroutine compute_zh
4254 
4255  subroutine get_staggered_grid( is, ie, js, je, isd, ied, jsd, jed, pt_b, pt_c, pt_d)
4256  integer, intent(in):: is, ie, js, je, isd, ied, jsd, jed
4257  real, dimension(isd:ied+1,jsd:jed+1,2), intent(in) :: pt_b
4258  real, dimension(isd:ied+1,jsd:jed ,2), intent(out) :: pt_c
4259  real, dimension(isd:ied ,jsd:jed+1,2), intent(out) :: pt_d
4260  ! local
4261  real(kind=R_GRID), dimension(2):: p1, p2, p3
4262  integer :: i, j
4263 
4264  do j=js,je+1
4265  do i=is,ie
4266  p1(:) = pt_b(i, j,1:2)
4267  p2(:) = pt_b(i+1,j,1:2)
4268  call mid_pt_sphere(p1, p2, p3)
4269  pt_d(i,j,1:2) = p3(:)
4270  enddo
4271  enddo
4272 
4273  do j=js,je
4274  do i=is,ie+1
4275  p1(:) = pt_b(i,j ,1:2)
4276  p2(:) = pt_b(i,j+1,1:2)
4277  call mid_pt_sphere(p1, p2, p3)
4278  pt_c(i,j,1:2) = p3(:)
4279  enddo
4280  enddo
4281 
4282  end subroutine get_staggered_grid
4283 
4284  end module external_ic_mod
4285 
subroutine remap_dwinds(km, npz, ak0, bk0, psc, ud, vd, Atm)
subroutine, public mid_pt_sphere(p1, p2, pm)
subroutine, public prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
subroutine, public del2_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, cd, zero_ocean, oro, nested, domain, bd, regional)
pure real function, public vicpqd(q)
subroutine remap_scalar_ec(Atm, km, npz, ncnst, ak0, bk0, psc, qa, wc, zh)
subroutine, public get_var3_r4(ncid, var3_name, is, ie, js, je, ks, ke, var3, time_slice)
Definition: sim_nc_mod.F90:246
real, parameter zvir
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:120
integer, parameter, public h_stagger
subroutine, public mappm(km, pe1, q1, kn, pe2, q2, i1, i2, iv, kord, ptop)
The subroutine &#39;mappm&#39; is a general-purpose routine for remapping one set of vertical levels to anoth...
Definition: fv_mapz.F90:3278
subroutine, public set_eta(km, ks, ptop, ak, bk)
This is the version of set_eta used in fvGFS and AM4.
Definition: fv_eta.F90:420
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
subroutine compute_zh(im, jm, levp, ak0, bk0, ps, zs, t, q, nq, zh)
subroutine, public fv3_zs_filter(bd, isd, ied, jsd, jed, npx, npy, npx_global, stretch_fac, nested, domain, area, dxa, dya, dx, dy, dxc, dyc, grid, agrid, sin_sg, phis, oro, regional)
subroutine, public del4_cubed_sphere(npx, npy, q, area, dx, dy, dxc, dyc, sin_sg, nmax, zero_ocean, oro, nested, domain, bd, regional)
subroutine get_nggps_ic(Atm, fv_domain)
The subroutine &#39;get_nggps_ic&#39; reads in data after it has been preprocessed with NCEP/EMC orography ma...
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function &#39;g_sum&#39; is the fast version of &#39;globalsum&#39;.
real function, public inner_prod(v1, v2)
The module &#39;fv_io&#39; contains restart facilities for FV core.
Definition: fv_io.F90:30
subroutine get_staggered_grid(is, ie, js, je, isd, ied, jsd, jed, pt_b, pt_c, pt_d)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine, public get_unit_vect2(e1, e2, uc)
subroutine remap_xyz(im, jbeg, jend, jm, km, npz, nq, ncnst, lon, lat, ak0, bk0, ps0, gz0, ua, va, ta, qa, Atm)
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
pure real function, public virqd(q)
The module fv_nwp_nudge contains routines for nudging to input analyses. note This module is currentl...
Definition: fv_nudge.F90:32
real(kind=r_grid), parameter cnst_0p20
subroutine, public get_var1_real(ncid, var1_name, im, var1, var_exist)
Definition: sim_nc_mod.F90:116
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
pure real function, public virq(q)
logical, public t_is_tv
Definition: fv_nudge.F90:195
The module &#39;boundary&#39; contains utility routines for grid nesting and boundary conditions.
Definition: boundary.F90:25
subroutine pmaxmin(qname, a, im, jm, fac)
subroutine, public get_var1_double(ncid, var1_name, im, var1, var_exist)
The &#39;get_var&#39; subroutines read in variables from netcdf files.
Definition: sim_nc_mod.F90:93
real, parameter, public ptop_min
subroutine, public fv_io_read_tracers(fv_domain, Atm)
The subroutine &#39;fv_io_read_tracers&#39; reads in only tracers from restart files.
Definition: fv_io.F90:225
The module &#39;fv_mapz&#39; contains the vertical mapping routines .
Definition: fv_mapz.F90:27
subroutine, public fillz(im, km, nq, q, dp)
The subroutine &#39;fillz&#39; is for mass-conservative filling of nonphysical negative values in the tracers...
Definition: fv_fill.F90:52
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
interface &#39;nested_grid_BC&#39; includes subroutines &#39;nested_grid_BC_2d&#39; and &#39;nested_grid_BC_3d&#39; that fetc...
Definition: boundary.F90:88
The module &#39;sim_nc&#39; is a netcdf file reader.
Definition: sim_nc_mod.F90:28
subroutine ncep2fms(im, jm, lon, lat, wk)
subroutine, public set_external_eta(ak, bk, ptop, ks)
The subroutine &#39;set_external_eta&#39; sets &#39;ptop&#39; (model top) and &#39;ks&#39; (first level of pure pressure coor...
Definition: fv_eta.F90:1517
subroutine, public get_var2_real(ncid, var_name, im, jm, var2)
Definition: sim_nc_mod.F90:137
subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
The module &#39;fv_eta&#39; contains routine to set up the reference (Eulerian) pressure coordinate.
Definition: fv_eta.F90:25
subroutine cubed_a2d(npx, npy, npz, ua, va, u, v, gridstruct, fv_domain, bd)
The subroutine &#39;cubed_a2d&#39; transforms the wind from the A Grid to the D Grid.
subroutine, public get_var2_r4(ncid, var2_name, is, ie, js, je, var2, time_slice)
Definition: sim_nc_mod.F90:153
The module &#39;external_ic_mod&#39; contains routines that read in and remap initial conditions.
Definition: external_ic.F90:32
subroutine d2a3d(u, v, ua, va, im, jm, km, lon)
real, dimension(:,:), allocatable, public sgh_g
subroutine remap_scalar_nggps(Atm, km, npz, ncnst, ak0, bk0, psc, t_in, qa, omga, zh)
subroutine, public open_ncfile(iflnm, ncid)
Definition: sim_nc_mod.F90:55
subroutine remap_coef(is, ie, js, je, isd, ied, jsd, jed, im, jm, lon, lat, id1, id2, jdc, s2c, agrid)
subroutine, public get_data_source(source, regional)
subroutine, public get_var_att_double(ncid, var_name, att_name, value)
Definition: sim_nc_mod.F90:399
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
subroutine get_fv_ic(Atm, fv_domain, nq)
subroutine get_ncep_ic(Atm, fv_domain, nq)
The subroutine &#39;get_ncep_ic&#39; reads in the specified NCEP analysis or reanalysis dataset.
integer, parameter, public ng
Definition: fv_mp_mod.F90:2716
subroutine remap_winds(im, jm, km, npz, ak0, bk0, psc, ua, va, Atm)
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, regional)
integer, parameter, public v_stagger
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
subroutine, public prt_gb_nh_sh(qname, is, ie, js, je, a2, area, lat)
real, dimension(:,:), allocatable, public oro_g
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
character(len=80) source
subroutine, public get_cubed_sphere_terrain(Atm, fv_domain)
subroutine, public checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, nq, km, q, lon, lat, nx, ny, rn)
subroutine mp_auto_conversion(ql, qr, qi, qs)
subroutine remap_scalar_single(Atm, km, npz, ak0, bk0, psc, qa, zh, iq)
subroutine, public get_external_ic(Atm, fv_domain, cold_start)
subroutine remap_scalar(im, jm, km, npz, nq, ncnst, ak0, bk0, psc, gzc, ta, qa, Atm)
subroutine get_ecmwf_ic(Atm, fv_domain)
The subroutine &#39;get_ecmwf_ic&#39; reads in initial conditions from ECMWF analyses (EXPERIMENTAL: contact ...
subroutine, public get_ncdim1(ncid, var1_name, im)
Definition: sim_nc_mod.F90:78
subroutine, public start_regional_cold_start(Atm, ak, bk, levp, is, ie, js, je, isd, ied, jsd, jed)
subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
subroutine, public extrapolation_bc(q, istag, jstag, npx, npy, bd, pd_in, debug_in)
The subroutine &#39;extrapolation_BC&#39; performs linear extrapolation into the halo region.
Definition: boundary.F90:118
subroutine fillq(im, km, nq, q, dp)
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)
the subroutine &#39;p_var&#39; computes auxiliary pressure variables for a hydrostatic state.
Definition: init_hydro.F90:87
integer, parameter, public u_stagger
subroutine, public close_ncfile(ncid)
Definition: sim_nc_mod.F90:67