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