FV3DYCORE  Version 2.0.0
fv_regional_bc.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU Lesser General Public License
3 !*
4 !* This file is part of the FV3 dynamical core.
5 !*
6 !* The FV3 dynamical core is free software: you can redistribute it
7 !* and/or modify it under the terms of the
8 !* GNU Lesser General Public License as published by the
9 !* Free Software Foundation, either version 3 of the License, or
10 !* (at your option) any later version.
11 !*
12 !* The FV3 dynamical core is distributed in the hope that it will be
13 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
14 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
15 !* See the GNU General Public License for more details.
16 !*
17 !* You should have received a copy of the GNU Lesser General Public
18 !* License along with the FV3 dynamical core.
19 !* If not, see <http://www.gnu.org/licenses/>.
20 !***********************************************************************
21 
22 module fv_regional_mod
23 
24  use netcdf
25  use mpp_domains_mod, only: domain2d
26  use mpp_domains_mod, only: domain1d, mpp_get_domain_components, &
27  mpp_get_global_domain, &
28  mpp_get_data_domain, &
29  mpp_get_compute_domain, &
30  north, south, east, west, &
31  center, corner, &
32  mpp_domains_set_stack_size, &
33  mpp_update_domains, mpp_get_neighbor_pe
34  use mpp_mod, only: fatal, input_nml_file, &
35  mpp_error ,mpp_pe, mpp_sync, &
36  mpp_npes, mpp_root_pe, mpp_gather, &
37  mpp_get_current_pelist, note, null_pe
38  use mpp_io_mod
39  use tracer_manager_mod,only: get_tracer_index,get_tracer_names
40  use field_manager_mod, only: model_atmos
41  use time_manager_mod, only: get_time &
42  ,operator(-),operator(/) &
43  ,time_type,time_type_to_real
44  use constants_mod, only: cp_air, cp_vapor, grav, kappa &
45  ,pi=>pi_8,rdgas, rvgas
46  use fv_arrays_mod, only: fv_atmos_type &
47  ,fv_grid_bounds_type &
48  ,fv_regional_bc_bounds_type &
49  ,r_grid &
50  ,fv_nest_bc_type_3d &
51  ,allocate_fv_nest_bc_type
52 
53  use fv_diagnostics_mod,only: prt_gb_nh_sh, prt_height
54  use fv_grid_utils_mod, only: g_sum,mid_pt_sphere,get_unit_vect2 &
55  ,get_latlon_vector,inner_prod &
56  ,cell_center2
57  use fv_mapz_mod, only: mappm, moist_cp, moist_cv
58  use fv_mp_mod, only: is_master, mp_reduce_min, mp_reduce_max
59  use fv_fill_mod, only: fillz
60  use fv_eta_mod, only: get_eta_level
61  use fms_mod, only: check_nml_error,file_exist
62  use fms_io_mod, only: read_data,get_global_att_value
63  use boundary_mod, only: fv_nest_bc_type_3d
64 
65  implicit none
66 
67  private
68 
69  public ak_in, bk_in &
70  ,bc_hour &
71  ,bc_time_interval &
72  ,bc_t0,bc_t1 &
73  ,begin_regional_restart,exch_uv &
74  ,ntimesteps_per_bc_update &
75  ,read_new_bc_data &
76  ,regional_bc_data &
77  ,regional_bc_t1_to_t0 &
78  ,regional_boundary_update &
79  ,next_time_to_read_bcs &
80  ,set_regional_bcs &
81  ,setup_regional_bc &
82  ,start_regional_cold_start &
83  ,start_regional_restart &
84  ,dump_field &
85  ,current_time_in_seconds &
86  ,a_step, p_step, k_step, n_step, get_data_source &
87  ,write_full_fields
88  integer,parameter :: bc_time_interval=3 &
89  ,nhalo_data =4 &
90  ,nhalo_model=3
91 !
92  integer, public, parameter :: h_stagger = 1
93  integer, public, parameter :: u_stagger = 2
94  integer, public, parameter :: v_stagger = 3
95 
96  !These parameters are ONLY used for the dump_field debugging routines
97  real, parameter :: stretch_factor = 1.5
98  real, parameter :: target_lon = -97.5
99  real, parameter :: target_lat = 35.5
100  integer, parameter :: parent_tile = 6
101  integer, parameter :: refine_ratio = 3
102 
103  integer, parameter :: cube_res = 96
104  integer, parameter :: istart_nest = 26
105  integer, parameter :: jstart_nest = 36
106  integer, parameter :: iend_nest = 167
107  integer, parameter :: jend_nest = 165
108 
109 ! integer, parameter :: cube_res = 768
110 ! integer, parameter :: istart_nest = 191
111 ! integer, parameter :: jstart_nest = 327
112 ! integer, parameter :: iend_nest = 1346
113 ! integer, parameter :: jend_nest = 1290
114 
115  integer,parameter :: nvars_core=7 & !<-- # of prognostic variables in core restart file
116  ,ndims_core=6 &
117  ,ndims_tracers=4
118 !
119  real,parameter :: blend_exp1=0.5,blend_exp2=10.
120  ! for prescribed external values in the
121  ! blending rows inside the domain boundary.
122  real :: current_time_in_seconds
123 !
124  integer,save :: isd_mod,ied_mod,jsd_mod,jed_mod
125 !
126  integer,save :: ncid,next_time_to_read_bcs,nfields_tracers &
127  ,npz,ntracers
128 !
129  integer,save :: k_split,n_split
130 !
131  integer,save :: bc_hour, ntimesteps_per_bc_update
132 !
133  integer,save :: cld_amt_index & !<--
134  ,graupel_index & !
135  ,ice_water_index & ! Locations of
136  ,liq_water_index & ! tracer vbls
137  ,o3mr_index & ! in the tracers
138  ,rain_water_index & ! array.
139  ,snow_water_index & !
140  ,sphum_index
141 !
142  integer,save :: lbnd_x_tracers,lbnd_y_tracers & !<-- Local lower bounds of x,y for tracer arrays
143  ,ubnd_x_tracers,ubnd_y_tracers
144 !
145  integer,save :: nrows_blend
146 !
147  real,save :: dt_atmos & !<-- The physics (large) timestep (sec)
148  ,dyn_timestep
149 !
150  real(kind=R_GRID),dimension(:,:,:),allocatable :: agrid_reg & !<-- Lon/lat of cell centers
151  ,grid_reg
152 
153  real,dimension(:,:),allocatable :: phis_reg
154 
155  real,dimension(:),allocatable :: ak_in, bk_in
156 
157  logical,save :: north_bc,south_bc,east_bc,west_bc &
158  ,begin_regional_restart=.true.
159 
160  logical,dimension(:),allocatable,save :: blend_this_tracer
161 
162  character(len=50) :: filename_core='INPUT/fv_core.res.temp.nc'
163  character(len=50) :: filename_core_new='RESTART/fv_core.res.tile1_new.nc'
164  character(len=50) :: filename_tracers='INPUT/fv_tracer.res.temp.nc'
165  character(len=50) :: filename_tracers_new='RESTART/fv_tracer.res.tile1_new.nc'
166 
167  type fv_regional_bc_variables
168  real,dimension(:,:,:),allocatable :: delp_bc, divgd_bc, u_bc, v_bc, uc_bc, vc_bc
169  real,dimension(:,:,:,:),allocatable :: q_bc
170 #ifndef SW_DYNAMICS
171  real,dimension(:,:,:),allocatable :: pt_bc, w_bc, delz_bc
172 #ifdef USE_COND
173  real,dimension(:,:,:),allocatable :: q_con_bc
174 #ifdef MOIST_CAPPA
175  real,dimension(:,:,:),allocatable :: cappa_bc
176 #endif
177 #endif
178 #endif
179  end type fv_regional_bc_variables
180 
181  type fv_domain_sides
182  type(fv_regional_bc_variables) :: north, south, east, west
183  end type fv_domain_sides
184 
185  type single_vbl3d_sides
186  real,dimension(:,:,:),pointer :: north, south, east, west
187  end type single_vbl3d_sides
188 
189  type vars_2d
190  real,dimension(:,:),pointer :: ptr
191  character(len=10) :: name
192  end type vars_2d
193 
194  type vars_3d
195  real,dimension(:,:,:),pointer :: ptr
196  character(len=10) :: name
197  end type vars_3d
198 
199  type(fv_domain_sides),target,save :: bc_t0, bc_t1
200 
201  type(fv_regional_bc_variables),pointer,save :: bc_north_t0 &
202  ,bc_south_t0 &
203  ,bc_west_t0 &
204  ,bc_east_t0 &
205  ,bc_north_t1 &
206  ,bc_south_t1 &
207  ,bc_west_t1 &
208  ,bc_east_t1
209 
210  type(fv_regional_bc_variables),pointer :: bc_side_t0,bc_side_t1
211 
212  type(fv_regional_bc_bounds_type),pointer,save :: regional_bounds
213 
214  type(vars_3d),dimension(:),allocatable :: fields_core &
215  ,fields_tracers
216  type(fv_nest_bc_type_3d), public :: delz_regbc ! lmh
217 
218  type(single_vbl3d_sides) :: delz_auxiliary
219 
220  integer :: ns = 0 ! lmh
221 
222  real,parameter :: tice=273.16 &
223  ,t_i0=15.
224  real, parameter :: c_liq = 4185.5 ! gfdl: heat capacity of liquid at 15 deg c
225  real, parameter :: c_ice = 1972.0 ! gfdl: heat capacity of ice at - 15 deg c
226  real, parameter :: zvir = rvgas/rdgas - 1. &
227  ,cv_air = cp_air - rdgas &
228  ,cv_vap = cp_vapor - rvgas
229 
230  real,dimension(:),allocatable :: dum1d, pref
231 
232  character(len=100) :: grid_data='grid.tile7.halo4.nc' &
233  ,oro_data ='oro_data.tile7.halo4.nc'
234 
235 #ifdef OVERLOAD_R4
236  real, parameter:: real_snan=x'FFBFFFFF'
237 #else
238  real, parameter:: real_snan=x'FFF7FFFFFFFFFFFF'
239 #endif
240  real(kind=R_GRID), parameter:: dbl_snan=x'FFF7FFFFFFFFFFFF'
241 
242  interface dump_field
243  module procedure dump_field_3d
244  module procedure dump_field_2d
245  end interface dump_field
246 
247  integer,save :: bc_update_interval, nrows_blend_user
248 
249  integer :: a_step, p_step, k_step, n_step
250 !
251  character(len=80) :: data_source
252 contains
253 
254 !-----------------------------------------------------------------------
255 !
256  subroutine setup_regional_bc(Atm, dt_atmos &
257  ,isd,ied,jsd,jed &
258  ,npx,npy )
259 !
260 !-----------------------------------------------------------------------
261 !*** Regional boundary data is obtained from the external BC file.
262 !-----------------------------------------------------------------------
263  use netcdf
264 !-----------------------------------------------------------------------
265  implicit none
266 !-----------------------------------------------------------------------
267 !
268 !---------------------
269 !*** Input variables
270 !---------------------
271 !
272  integer,intent(in) :: isd,ied,jsd,jed,npx,npy
273 !
274  real,intent(in) :: dt_atmos
275 !
276  type(fv_atmos_type),target,intent(inout) :: Atm
277 !
278 !--------------------
279 !*** Local variables
280 !--------------------
281 !
282  integer :: dimid,i,i_start,i_end,j,j_start,j_end,klev_out &
283  ,nrows_bc_data,nrows_blend_in_data,sec
284 !
285  real :: ps1
286 !
287  character(len=2) :: char2_1,char2_2
288  character(len=3) :: int_to_char
289  character(len=6) :: fmt='(i3.3)'
290  character(len=50) :: file_name
291 !
292 !-----------------------------------------------------------------------
293 !***********************************************************************
294 !-----------------------------------------------------------------------
295 !
296 !-----------------------------------------------------------------------
297 !*** The boundary data is laid out so that the pieces for the north
298 !*** and south sides span the entire distance from the east side of
299 !*** of the east halo to the west side of the west halo. Therefore
300 !*** there the # of cells in the x direction in the north/south BC
301 !*** data is nx+2*nhalo where nx is the # of cells in the x direction
302 !*** on the compute domain. This means the # of cells spanned in the
303 !*** west/east side BC data is just ny (the # of cells in the y
304 !*** direction on the compute domain) and not ny+2*nhalo since the
305 !*** halo values on the south and north ends of the east/west sides
306 !*** are already part of the BC data on the north/south sides.
307 !-----------------------------------------------------------------------
308 !
309 ! nhalo_model=3
310 !
311 ! |----------- nxp-1 -----------| <-- east/west compute points
312 ! |---------- north BC data ----------|
313 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
314 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
315 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
316 ! --- ooo ---j=1--- ooo --- ---
317 ! | ooo ooo | |
318 ! | ooo |ooo | |
319 ! ooo i=1-->|ooo
320 ! west BC data ooo| |ooo east BC data nyp-1 <-- north/south compute points
321 ! ooo|<--i=isd-nhalo_model ooo
322 ! | ooo| ooo | |
323 ! | ooo ooo | |
324 ! --- ooo ---j=jsd-nhalo_model--- ooo --- ---
325 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
326 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
327 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
328 ! |---------- south BC data ----------|
329 !
330 !-----------------------------------------------------------------------
331 !
332  north_bc=.false.
333  south_bc=.false.
334  east_bc =.false.
335  west_bc =.false.
336 !
337 ! write(0,*)' enter setup_regional_BC isd=',isd,' ied=',ied,' jsd=',jsd,' jed=',jed
338  isd_mod=isd
339  ied_mod=ied
340  jsd_mod=jsd
341  jed_mod=jed
342 !
343 !-----------------------------------------------------------------------
344 !*** Which side(s) of the domain does this task lie on if any?
345 !-----------------------------------------------------------------------
346 !
347  if(jsd<0)then
348  north_bc=.true.
349  endif
350 
351  if(jed>npy-1)then
352  south_bc=.true.
353  endif
354 
355  if(isd<0)then
356  east_bc=.true.
357  endif
358 
359  if(ied>npx-1)then
360  west_bc=.true.
361  endif
362 !
363  bc_update_interval=atm%flagstruct%bc_update_interval
364 !
365  k_split=atm%flagstruct%k_split
366  n_split=atm%flagstruct%n_split
367 !
368  dyn_timestep=dt_atmos/real(k_split*n_split)
369 !
370 !-----------------------------------------------------------------------
371 !*** Is blending row data present in the BC file and if so how many
372 !*** rows of data are there? All blending data that is present will
373 !*** be read even if the user requests fewer rows be applied.
374 !*** Construct the name of the regional BC file to be read.
375 !*** We must know whether this is a standard BC file from chgres
376 !*** or a new BC file generated from DA-updated data from enlarged
377 !*** restart files that include the boundary rows.
378 !-----------------------------------------------------------------------
379 !
380  write(int_to_char,fmt) bc_hour
381  if(.not.atm%flagstruct%regional_bcs_from_gsi)then
382  file_name='INPUT/gfs_bndy.tile7.'//int_to_char//'.nc'
383  else
384  file_name='INPUT/gfs_bndy.tile7.'//int_to_char//'_gsi.nc'
385  endif
386 !
387  if (is_master()) then
388  write(*,20011)trim(file_name)
389 20011 format(' regional_bc_data file_name=',a)
390  endif
391 !-----------------------------------------------------------------------
392 !*** Open the regional BC file.
393 !-----------------------------------------------------------------------
394 !
395  call check(nf90_open(file_name,nf90_nowrite,ncid))
396  if (is_master()) then
397  write(0,*)' opened BC file ',trim(file_name)
398  endif
399 !
400 !-----------------------------------------------------------------------
401 !*** Check if the desired number of blending rows are present in
402 !*** the boundary files.
403 !-----------------------------------------------------------------------
404 !
405  nrows_blend_user=atm%flagstruct%nrows_blend
406 !
407  call check(nf90_inq_dimid(ncid,'halo',dimid))
408  call check(nf90_inquire_dimension(ncid,dimid,len=nrows_bc_data))
409 !
410  nrows_blend_in_data=nrows_bc_data-nhalo_data
411 !
412  if(nrows_blend_user>nrows_blend_in_data)then
413  write(char2_1,'(I2.2)')nrows_blend_user
414  write(char2_2,'(I2.2)')nrows_blend_in_data
415  call mpp_error(fatal,'User wants to use '//char2_1//' blending rows but only '//char2_2//' blending rows are in the BC file!')
416  else
417  nrows_blend=nrows_blend_in_data
418  endif
419 !
420  call check(nf90_close(ncid))
421 !
422 !-----------------------------------------------------------------------
423 !*** Compute the index limits within the boundary region on each
424 !*** side of the domain for both scalars and winds. Since the
425 !*** domain does not move then the computations need to be done
426 !*** only once.
427 !-----------------------------------------------------------------------
428 !
429  call compute_regional_bc_indices(atm%regional_bc_bounds)
430 !
431 !-----------------------------------------------------------------------
432 !
433  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
434  return
435  endif
436 !
437 !-----------------------------------------------------------------------
438 !
439 ! ntracers=Atm%ncnst - Atm%flagstruct%dnats !<-- # of advected tracers
440  ntracers=atm%ncnst
441  npz=atm%npz
442  klev_out=npz
443 !
444  regional_bounds=>atm%regional_bc_bounds
445 !
446 !-----------------------------------------------------------------------
447 !*** Allocate the objects that will hold the boundary variables
448 !*** at the two time levels surrounding each piece of the regional
449 !*** domain's integration. Data is read from the BC files into
450 !*** time level t1 while t0 holds the data from the preceding
451 !*** BC file.
452 !-----------------------------------------------------------------------
453 !*** Point pointers at each side's boundary data for both time levels.
454 !*** Those are needed when the actual update of boundary points is
455 !*** executed.
456 !-----------------------------------------------------------------------
457 !
458  if(north_bc)then
459  call allocate_regional_bc_arrays('north' &
460  ,north_bc,south_bc &
461  ,east_bc,west_bc &
462  ,atm%regional_bc_bounds%is_north &
463  ,atm%regional_bc_bounds%ie_north &
464  ,atm%regional_bc_bounds%js_north &
465  ,atm%regional_bc_bounds%je_north &
466  ,atm%regional_bc_bounds%is_north_uvs &
467  ,atm%regional_bc_bounds%ie_north_uvs &
468  ,atm%regional_bc_bounds%js_north_uvs &
469  ,atm%regional_bc_bounds%je_north_uvs &
470  ,atm%regional_bc_bounds%is_north_uvw &
471  ,atm%regional_bc_bounds%ie_north_uvw &
472  ,atm%regional_bc_bounds%js_north_uvw &
473  ,atm%regional_bc_bounds%je_north_uvw &
474  ,klev_out &
475  ,ntracers &
476  ,bc_t1%north &
477  ,delz_auxiliary%north )
478 !
479  call allocate_regional_bc_arrays('north' &
480  ,north_bc,south_bc &
481  ,east_bc,west_bc &
482  ,atm%regional_bc_bounds%is_north &
483  ,atm%regional_bc_bounds%ie_north &
484  ,atm%regional_bc_bounds%js_north &
485  ,atm%regional_bc_bounds%je_north &
486  ,atm%regional_bc_bounds%is_north_uvs &
487  ,atm%regional_bc_bounds%ie_north_uvs &
488  ,atm%regional_bc_bounds%js_north_uvs &
489  ,atm%regional_bc_bounds%je_north_uvs &
490  ,atm%regional_bc_bounds%is_north_uvw &
491  ,atm%regional_bc_bounds%ie_north_uvw &
492  ,atm%regional_bc_bounds%js_north_uvw &
493  ,atm%regional_bc_bounds%je_north_uvw &
494  ,klev_out &
495  ,ntracers &
496  ,bc_t0%north )
497 !
498  bc_north_t0=>bc_t0%north
499  bc_north_t1=>bc_t1%north
500 !
501  endif
502 
503  if(south_bc)then
504  call allocate_regional_bc_arrays('south' &
505  ,north_bc,south_bc &
506  ,east_bc,west_bc &
507  ,atm%regional_bc_bounds%is_south &
508  ,atm%regional_bc_bounds%ie_south &
509  ,atm%regional_bc_bounds%js_south &
510  ,atm%regional_bc_bounds%je_south &
511  ,atm%regional_bc_bounds%is_south_uvs &
512  ,atm%regional_bc_bounds%ie_south_uvs &
513  ,atm%regional_bc_bounds%js_south_uvs &
514  ,atm%regional_bc_bounds%je_south_uvs &
515  ,atm%regional_bc_bounds%is_south_uvw &
516  ,atm%regional_bc_bounds%ie_south_uvw &
517  ,atm%regional_bc_bounds%js_south_uvw &
518  ,atm%regional_bc_bounds%je_south_uvw &
519  ,klev_out &
520  ,ntracers &
521  ,bc_t1%south &
522  ,delz_auxiliary%south )
523 !
524  call allocate_regional_bc_arrays('south' &
525  ,north_bc,south_bc &
526  ,east_bc,west_bc &
527  ,atm%regional_bc_bounds%is_south &
528  ,atm%regional_bc_bounds%ie_south &
529  ,atm%regional_bc_bounds%js_south &
530  ,atm%regional_bc_bounds%je_south &
531  ,atm%regional_bc_bounds%is_south_uvs &
532  ,atm%regional_bc_bounds%ie_south_uvs &
533  ,atm%regional_bc_bounds%js_south_uvs &
534  ,atm%regional_bc_bounds%je_south_uvs &
535  ,atm%regional_bc_bounds%is_south_uvw &
536  ,atm%regional_bc_bounds%ie_south_uvw &
537  ,atm%regional_bc_bounds%js_south_uvw &
538  ,atm%regional_bc_bounds%je_south_uvw &
539  ,klev_out &
540  ,ntracers &
541  ,bc_t0%south )
542 !
543  bc_south_t0=>bc_t0%south
544  bc_south_t1=>bc_t1%south
545 !
546  endif
547 !
548  if(east_bc)then
549  call allocate_regional_bc_arrays('east ' &
550  ,north_bc,south_bc &
551  ,east_bc,west_bc &
552  ,atm%regional_bc_bounds%is_east &
553  ,atm%regional_bc_bounds%ie_east &
554  ,atm%regional_bc_bounds%js_east &
555  ,atm%regional_bc_bounds%je_east &
556  ,atm%regional_bc_bounds%is_east_uvs &
557  ,atm%regional_bc_bounds%ie_east_uvs &
558  ,atm%regional_bc_bounds%js_east_uvs &
559  ,atm%regional_bc_bounds%je_east_uvs &
560  ,atm%regional_bc_bounds%is_east_uvw &
561  ,atm%regional_bc_bounds%ie_east_uvw &
562  ,atm%regional_bc_bounds%js_east_uvw &
563  ,atm%regional_bc_bounds%je_east_uvw &
564  ,klev_out &
565  ,ntracers &
566  ,bc_t1%east &
567  ,delz_auxiliary%east )
568 !
569  call allocate_regional_bc_arrays('east ' &
570  ,north_bc,south_bc &
571  ,east_bc,west_bc &
572  ,atm%regional_bc_bounds%is_east &
573  ,atm%regional_bc_bounds%ie_east &
574  ,atm%regional_bc_bounds%js_east &
575  ,atm%regional_bc_bounds%je_east &
576  ,atm%regional_bc_bounds%is_east_uvs &
577  ,atm%regional_bc_bounds%ie_east_uvs &
578  ,atm%regional_bc_bounds%js_east_uvs &
579  ,atm%regional_bc_bounds%je_east_uvs &
580  ,atm%regional_bc_bounds%is_east_uvw &
581  ,atm%regional_bc_bounds%ie_east_uvw &
582  ,atm%regional_bc_bounds%js_east_uvw &
583  ,atm%regional_bc_bounds%je_east_uvw &
584  ,klev_out &
585  ,ntracers &
586  ,bc_t0%east )
587 !
588  bc_east_t0=>bc_t0%east
589  bc_east_t1=>bc_t1%east
590 !
591  endif
592 !
593  if(west_bc)then
594  call allocate_regional_bc_arrays('west ' &
595  ,north_bc,south_bc &
596  ,east_bc,west_bc &
597  ,atm%regional_bc_bounds%is_west &
598  ,atm%regional_bc_bounds%ie_west &
599  ,atm%regional_bc_bounds%js_west &
600  ,atm%regional_bc_bounds%je_west &
601  ,atm%regional_bc_bounds%is_west_uvs &
602  ,atm%regional_bc_bounds%ie_west_uvs &
603  ,atm%regional_bc_bounds%js_west_uvs &
604  ,atm%regional_bc_bounds%je_west_uvs &
605  ,atm%regional_bc_bounds%is_west_uvw &
606  ,atm%regional_bc_bounds%ie_west_uvw &
607  ,atm%regional_bc_bounds%js_west_uvw &
608  ,atm%regional_bc_bounds%je_west_uvw &
609  ,klev_out &
610  ,ntracers &
611  ,bc_t1%west &
612  ,delz_auxiliary%west )
613 !
614  call allocate_regional_bc_arrays('west ' &
615  ,north_bc,south_bc &
616  ,east_bc,west_bc &
617  ,atm%regional_bc_bounds%is_west &
618  ,atm%regional_bc_bounds%ie_west &
619  ,atm%regional_bc_bounds%js_west &
620  ,atm%regional_bc_bounds%je_west &
621  ,atm%regional_bc_bounds%is_west_uvs &
622  ,atm%regional_bc_bounds%ie_west_uvs &
623  ,atm%regional_bc_bounds%js_west_uvs &
624  ,atm%regional_bc_bounds%je_west_uvs &
625  ,atm%regional_bc_bounds%is_west_uvw &
626  ,atm%regional_bc_bounds%ie_west_uvw &
627  ,atm%regional_bc_bounds%js_west_uvw &
628  ,atm%regional_bc_bounds%je_west_uvw &
629  ,klev_out &
630  ,ntracers &
631  ,bc_t0%west )
632 !
633  bc_west_t0=>bc_t0%west
634  bc_west_t1=>bc_t1%west
635 !
636  endif
637 
638  call allocate_fv_nest_bc_type(delz_regbc,atm,ns,0,0,.false.)
639 !
640 !-----------------------------------------------------------------------
641 !*** We need regional versions of the arrays for surface elevation,
642 !*** latitude/longitude of grid cell corners, and lat/lon of the
643 !*** cell centers because those variables are needed an extra row
644 !*** beyond FV3's normal bounday region width of nhalo_model rows.
645 !-----------------------------------------------------------------------
646 !
647  allocate(phis_reg(isd-1:ied+1,jsd-1:jed+1)) ; phis_reg=real_snan
648 !
649  allocate(agrid_reg(isd-1:ied+1,jsd-1:jed+1,2)); agrid_reg=dbl_snan
650  allocate(grid_reg(isd-1:ied+2,jsd-1:jed+2,2)) ; grid_reg=dbl_snan
651 !
652 !-----------------------------------------------------------------------
653 !*** From the data holding nhalo_model rows of boundary values
654 !*** read in the lat/lon of the grid cell corners and fill in
655 !*** the values of the grid cell centers. The regional mode needs
656 !*** the extra row of data.
657 !-----------------------------------------------------------------------
658 !
659  call read_regional_lon_lat
660 !
661 !-----------------------------------------------------------------------
662 !*** From the data holding nhalo_model rows of filtered topography
663 !*** read in those values. The regional mode needs the extra row
664 !*** of data.
665 !-----------------------------------------------------------------------
666 !
667  call read_regional_filtered_topo
668 !
669 !-----------------------------------------------------------------------
670 !*** In the init step Atm%phis is given values only in the integration
671 !*** domain but in a regional run values are also needed in the
672 !*** boundary rows. Since the same data is read in the preceding
673 !*** subroutine call as when Atm%phis was first filled, fill its
674 !*** boundary rows now.
675 !-----------------------------------------------------------------------
676 !
677  if(north_bc)then
678  i_start=isd
679  i_end =ied
680  j_start=jsd
681  if(.not.atm%flagstruct%warm_start)then
682  j_end=jsd+nhalo_model-1
683  else
684  j_end=jsd+nhalo_model+1
685  endif
686  do j=j_start,j_end
687  do i=i_start,i_end
688  atm%phis(i,j)=phis_reg(i,j)
689  enddo
690  enddo
691  endif
692 !
693  if(south_bc)then
694  i_start=isd
695  i_end =ied
696  j_end =jed
697  if(.not.atm%flagstruct%warm_start)then
698  j_start=jed-nhalo_model+1
699  else
700  j_start=jed-nhalo_model-1
701  endif
702  do j=j_start,j_end
703  do i=i_start,i_end
704  atm%phis(i,j)=phis_reg(i,j)
705  enddo
706  enddo
707  endif
708  if(east_bc)then
709  i_start=isd
710  j_start=jsd
711  j_end =jed
712  if(.not.atm%flagstruct%warm_start)then
713  i_end=isd+nhalo_model-1
714  else
715  i_end=isd+nhalo_model+1
716  endif
717  do j=j_start,j_end
718  do i=i_start,i_end
719  atm%phis(i,j)=phis_reg(i,j)
720  enddo
721  enddo
722  endif
723  if(west_bc)then
724  i_end =ied
725  j_start=jsd
726  j_end =jed
727  if(.not.atm%flagstruct%warm_start)then
728  i_start=ied-nhalo_model+1
729  else
730  i_start=ied-nhalo_model-1
731  endif
732  do j=j_start,j_end
733  do i=i_start,i_end
734  atm%phis(i,j)=phis_reg(i,j)
735  enddo
736  enddo
737  endif
738 !
739  sphum_index = get_tracer_index(model_atmos, 'sphum')
740  liq_water_index = get_tracer_index(model_atmos, 'liq_wat')
741  ice_water_index = get_tracer_index(model_atmos, 'ice_wat')
742  rain_water_index = get_tracer_index(model_atmos, 'rainwat')
743  snow_water_index = get_tracer_index(model_atmos, 'snowwat')
744  graupel_index = get_tracer_index(model_atmos, 'graupel')
745  cld_amt_index = get_tracer_index(model_atmos, 'cld_amt')
746  o3mr_index = get_tracer_index(model_atmos, 'o3mr')
747 ! write(0,*)' setup_regional_bc'
748 ! write(0,*)' sphum_index=',sphum_index
749 ! write(0,*)' liq_water_index=',liq_water_index
750 ! write(0,*)' ice_water_index=',ice_water_index
751 ! write(0,*)' rain_water_index=',rain_water_index
752 ! write(0,*)' snow_water_index=',snow_water_index
753 ! write(0,*)' graupel_index=',graupel_index
754 ! write(0,*)' cld_amt_index=',cld_amt_index
755 ! write(0,*)' o3mr_index=',o3mr_index
756 !
757 !-----------------------------------------------------------------------
758 !*** When nudging of specific humidity is selected then we need a
759 !*** reference pressure profile. Compute it now.
760 !-----------------------------------------------------------------------
761 !
762  allocate(pref(npz+1))
763  allocate(dum1d(npz+1))
764 !
765  ps1=101325.
766  pref(npz+1)=ps1
767  call get_eta_level(npz,ps1,pref(1),dum1d,atm%ak,atm%bk )
768 !
769 !-----------------------------------------------------------------------
770 
771  contains
772 
773 !-----------------------------------------------------------------------
774 !
775  subroutine compute_regional_bc_indices(regional_bc_bounds)
776 !
777 !-----------------------------------------------------------------------
778 !*** This routine computes the starting and ending indices for
779 !*** working arrays of task subdomains that lie on the edges
780 !*** of the FV3 regional domain. These arrays will hold boundary
781 !*** region values of scalar variables located at the grid cell
782 !*** centers and wind components lying on the east/west sides
783 !*** and north/south sides of each cell. Note that the width
784 !*** of the domain's boundary region (4 rows) is currently
785 !*** greater than the fundamental width of the task subdomains'
786 !*** halo regions (3 rows). The variables isd,ied,jsd,jed are
787 !*** the task subdomain index limits including their halos.
788 !*** The diagram in subroutine regional_bc_data will help to
789 !*** understand these index limits have the values they do.
790 !-----------------------------------------------------------------------
791 !
792 !------------------------
793 !*** Argument variables
794 !------------------------
795 !
796  type(fv_regional_bc_bounds_type),intent(out) :: regional_bc_bounds
797 !
798 !---------------------
799 !*** Local variables
800 !---------------------
801 !
802  integer, parameter :: invalid_index = -99
803  integer :: halo_diff
804 !
805 !-----------------------------------------------------------------------
806 !***********************************************************************
807 !-----------------------------------------------------------------------
808 !
809  regional_bc_bounds%is_north = invalid_index
810  regional_bc_bounds%ie_north = invalid_index
811  regional_bc_bounds%js_north = invalid_index
812  regional_bc_bounds%je_north = invalid_index
813  regional_bc_bounds%is_north_uvs = invalid_index
814  regional_bc_bounds%ie_north_uvs = invalid_index
815  regional_bc_bounds%js_north_uvs = invalid_index
816  regional_bc_bounds%je_north_uvs = invalid_index
817  regional_bc_bounds%is_north_uvw = invalid_index
818  regional_bc_bounds%ie_north_uvw = invalid_index
819  regional_bc_bounds%js_north_uvw = invalid_index
820  regional_bc_bounds%je_north_uvw = invalid_index
821 
822  regional_bc_bounds%is_south = invalid_index
823  regional_bc_bounds%ie_south = invalid_index
824  regional_bc_bounds%js_south = invalid_index
825  regional_bc_bounds%je_south = invalid_index
826  regional_bc_bounds%is_south_uvs = invalid_index
827  regional_bc_bounds%ie_south_uvs = invalid_index
828  regional_bc_bounds%js_south_uvs = invalid_index
829  regional_bc_bounds%je_south_uvs = invalid_index
830  regional_bc_bounds%is_south_uvw = invalid_index
831  regional_bc_bounds%ie_south_uvw = invalid_index
832  regional_bc_bounds%js_south_uvw = invalid_index
833  regional_bc_bounds%je_south_uvw = invalid_index
834 
835  regional_bc_bounds%is_east = invalid_index
836  regional_bc_bounds%ie_east = invalid_index
837  regional_bc_bounds%js_east = invalid_index
838  regional_bc_bounds%je_east = invalid_index
839  regional_bc_bounds%is_east_uvs = invalid_index
840  regional_bc_bounds%ie_east_uvs = invalid_index
841  regional_bc_bounds%js_east_uvs = invalid_index
842  regional_bc_bounds%je_east_uvs = invalid_index
843  regional_bc_bounds%is_east_uvw = invalid_index
844  regional_bc_bounds%ie_east_uvw = invalid_index
845  regional_bc_bounds%js_east_uvw = invalid_index
846  regional_bc_bounds%je_east_uvw = invalid_index
847 
848  regional_bc_bounds%is_west = invalid_index
849  regional_bc_bounds%ie_west = invalid_index
850  regional_bc_bounds%js_west = invalid_index
851  regional_bc_bounds%je_west = invalid_index
852  regional_bc_bounds%is_west_uvs = invalid_index
853  regional_bc_bounds%ie_west_uvs = invalid_index
854  regional_bc_bounds%js_west_uvs = invalid_index
855  regional_bc_bounds%je_west_uvs = invalid_index
856  regional_bc_bounds%is_west_uvw = invalid_index
857  regional_bc_bounds%ie_west_uvw = invalid_index
858  regional_bc_bounds%js_west_uvw = invalid_index
859  regional_bc_bounds%je_west_uvw = invalid_index
860 !
861 !-----------------------------------------------------------------------
862 !*** Scalar BC indices
863 !-----------------------------------------------------------------------
864 !*** These must reach one row beyond nhalo_model since we must
865 !*** surround the wind points on the cell edges with mass points.
866 !
867 !*** NOTE: The value of nrows_blend is the total number of
868 !*** blending rows in the BC files.
869 !-----------------------------------------------------------------------
870 !
871  halo_diff=nhalo_data-nhalo_model
872 !
873 !-----------
874 !*** North
875 !-----------
876 !
877  if (north_bc) then
878  regional_bc_bounds%is_north=isd-1
879  regional_bc_bounds%ie_north=ied+1
880 !
881  regional_bc_bounds%js_north=jsd-1
882  regional_bc_bounds%je_north=nrows_blend
883  endif
884 !
885 !-----------
886 !*** South
887 !-----------
888 !
889  if (south_bc) then
890  regional_bc_bounds%is_south=isd-1
891  regional_bc_bounds%ie_south=ied+1
892 !
893  regional_bc_bounds%js_south=jed-nhalo_model-nrows_blend+1
894  regional_bc_bounds%je_south=jed+1
895  endif
896 !
897 !----------
898 !*** East
899 !----------
900 !
901  if (east_bc) then
902  regional_bc_bounds%is_east=isd-1
903  regional_bc_bounds%ie_east=nrows_blend
904 !
905  regional_bc_bounds%js_east=jsd-1
906  if(north_bc)then
907  regional_bc_bounds%js_east=1
908  endif
909 !
910  regional_bc_bounds%je_east=jed+1
911  if(south_bc)then
912  regional_bc_bounds%je_east=jed-nhalo_model
913  endif
914  endif
915 !
916 !----------
917 !*** West
918 !----------
919 !
920  if (west_bc) then
921  regional_bc_bounds%is_west=ied-nhalo_model-nrows_blend+1
922  regional_bc_bounds%ie_west=ied+1
923 !
924  regional_bc_bounds%js_west=jsd-1
925  if(north_bc)then
926  regional_bc_bounds%js_west=1
927  endif
928 !
929  regional_bc_bounds%je_west=jed+1
930  if(south_bc)then
931  regional_bc_bounds%je_west=jed-nhalo_model
932  endif
933  endif
934 !
935 !-----------------------------------------------------------------------
936 !*** Wind component BC indices
937 !-----------------------------------------------------------------------
938 !
939 !-----------
940 !*** North
941 !-----------
942 !
943  if (north_bc) then
944  regional_bc_bounds%is_north_uvs=isd
945  regional_bc_bounds%ie_north_uvs=ied
946 !
947  regional_bc_bounds%js_north_uvs=jsd
948  regional_bc_bounds%je_north_uvs=nrows_blend+1
949 !
950  regional_bc_bounds%is_north_uvw=isd
951  regional_bc_bounds%ie_north_uvw=ied+1
952 !
953  regional_bc_bounds%js_north_uvw=jsd
954  regional_bc_bounds%je_north_uvw=nrows_blend
955  endif
956 !
957 !-----------
958 !*** South
959 !-----------
960 !
961  if (south_bc) then
962  regional_bc_bounds%is_south_uvs=isd
963  regional_bc_bounds%ie_south_uvs=ied
964 !
965  regional_bc_bounds%js_south_uvs=jed-nhalo_model-nrows_blend+1
966  regional_bc_bounds%je_south_uvs=jed+1
967 !
968  regional_bc_bounds%is_south_uvw=isd
969  regional_bc_bounds%ie_south_uvw=ied+1
970 !
971  regional_bc_bounds%js_south_uvw=jed-nhalo_model-nrows_blend+1
972  regional_bc_bounds%je_south_uvw=jed
973  endif
974 !
975 !----------
976 !*** East
977 !----------
978 !
979  if (east_bc) then
980  regional_bc_bounds%is_east_uvs=isd
981  regional_bc_bounds%ie_east_uvs=nrows_blend
982 !
983  regional_bc_bounds%js_east_uvs=jsd
984  if(north_bc)then
985  regional_bc_bounds%js_east_uvs=1
986  endif
987 !
988  regional_bc_bounds%je_east_uvs=jed+1
989  if(south_bc)then
990  regional_bc_bounds%je_east_uvs=jed-nhalo_model+1
991  endif
992 !
993 ! regional_bc_bounds%is_east_uvw=isd-1
994  regional_bc_bounds%is_east_uvw=isd
995  regional_bc_bounds%ie_east_uvw=nrows_blend
996 !
997  regional_bc_bounds%js_east_uvw=jsd
998  if(north_bc)then
999  regional_bc_bounds%js_east_uvw=1
1000  endif
1001  regional_bc_bounds%je_east_uvw=jed
1002  if(south_bc)then
1003  regional_bc_bounds%je_east_uvw=jed-nhalo_model
1004  endif
1005  endif
1006 !
1007 !----------
1008 !*** West
1009 !----------
1010 !
1011  if (west_bc) then
1012  regional_bc_bounds%is_west_uvs=ied-nhalo_model-nrows_blend+1
1013  regional_bc_bounds%ie_west_uvs=ied
1014 !
1015  regional_bc_bounds%js_west_uvs=jsd
1016  if(north_bc)then
1017  regional_bc_bounds%js_west_uvs=1
1018  endif
1019 !
1020  regional_bc_bounds%je_west_uvs=jed+1
1021  if(south_bc)then
1022  regional_bc_bounds%je_west_uvs=jed-nhalo_model+1
1023  endif
1024 !
1025  regional_bc_bounds%is_west_uvw=ied-nhalo_model-nrows_blend+1
1026  regional_bc_bounds%ie_west_uvw=ied+1
1027 !
1028  regional_bc_bounds%js_west_uvw=jsd
1029  if(north_bc)then
1030  regional_bc_bounds%js_west_uvw=1
1031  endif
1032 !
1033  regional_bc_bounds%je_west_uvw=jed
1034  if(south_bc)then
1035  regional_bc_bounds%je_west_uvw=jed-nhalo_model
1036  endif
1037  endif
1038 !
1039 !-----------------------------------------------------------------------
1040 !
1041  end subroutine compute_regional_bc_indices
1042 !
1043 !-----------------------------------------------------------------------
1044 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1045 !-----------------------------------------------------------------------
1046 !
1047  subroutine read_regional_lon_lat
1048 !
1049 !-----------------------------------------------------------------------
1050 !*** Read the longitude/latitude of the grid cell corners from
1051 !*** the external file holding the additional row of data required
1052 !*** by the regional domain.
1053 !-----------------------------------------------------------------------
1054  use netcdf
1055 !-----------------------------------------------------------------------
1056  implicit none
1057 !-----------------------------------------------------------------------
1058 !
1059 !--------------------
1060 !*** Local variables
1061 !--------------------
1062 !
1063  integer :: i_start_data,istat,j_start_data,n,ncid_grid,var_id
1064 !
1065  character(len=150) :: filename,vname
1066 !
1067 !-----------------------------------------------------------------------
1068 !***********************************************************************
1069 !-----------------------------------------------------------------------
1070 !
1071 !-----------------------------------------------------------------------
1072 !*** Open the grid data file.
1073 !-----------------------------------------------------------------------
1074 !
1075  filename='INPUT/'//trim(grid_data)
1076 !
1077  call check(nf90_open(filename,nf90_nowrite,ncid_grid))
1078 !
1079  call mpp_error(note,' opened grid file '//trim(filename))
1080 !
1081 !-----------------------------------------------------------------------
1082 !*** The longitude and latitude are on the super grid. We need only
1083 !*** the points on each corner of the grid cells which is every other
1084 !*** point on the super grid.
1085 !-----------------------------------------------------------------------
1086 !
1087  i_start_data=2*(isd+nhalo_model)-1
1088  j_start_data=2*(jsd+nhalo_model)-1
1089 !
1090 !---------------
1091 !*** Longitude
1092 !---------------
1093 !
1094  vname='x'
1095  call check(nf90_inq_varid(ncid_grid,vname,var_id))
1096  call check(nf90_get_var(ncid_grid,var_id &
1097  ,grid_reg(isd-1:ied+2,jsd-1:jed+2,1) &
1098  ,start=(/i_start_data,j_start_data/) &
1099  ,stride=(/2,2/) ) )
1100 !
1101 !--------------
1102 !*** Latitude
1103 !--------------
1104 !
1105  vname='y'
1106  call check(nf90_inq_varid(ncid_grid,vname,var_id))
1107  call check(nf90_get_var(ncid_grid,var_id &
1108  ,grid_reg(isd-1:ied+2,jsd-1:jed+2,2) &
1109  ,start=(/i_start_data,j_start_data/) &
1110  ,stride=(/2,2/) ) )
1111 !
1112  call check(nf90_close(ncid_grid))
1113 !
1114 !-----------------------------------------------------------------------
1115 !*** Convert from degrees to radians.
1116 !-----------------------------------------------------------------------
1117 !
1118  do n=1,2
1119  do j=jsd-1,jed+2
1120  do i=isd-1,ied+2
1121  grid_reg(i,j,n)=grid_reg(i,j,n)*pi/180.
1122  enddo
1123  enddo
1124  enddo
1125 !
1126 !-----------------------------------------------------------------------
1127 !*** Compute the longitude/latitude in the cell centers.
1128 !-----------------------------------------------------------------------
1129 !
1130  do j=jsd-1,jed+1
1131  do i=isd-1,ied+1
1132  call cell_center2(grid_reg(i,j, 1:2), grid_reg(i+1,j, 1:2), &
1133  grid_reg(i,j+1,1:2), grid_reg(i+1,j+1,1:2), &
1134  agrid_reg(i,j,1:2) )
1135  enddo
1136  enddo
1137 !
1138 !-----------------------------------------------------------------------
1139 !
1140  end subroutine read_regional_lon_lat
1141 !
1142 !-----------------------------------------------------------------------
1143 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1144 !-----------------------------------------------------------------------
1145 !
1146  subroutine read_regional_filtered_topo
1147 !
1148 !-----------------------------------------------------------------------
1149 !*** Read the filtered topography including the extra outer row.
1150 !-----------------------------------------------------------------------
1151 !-----------------------------------------------------------------------
1152  implicit none
1153 !-----------------------------------------------------------------------
1154 !
1155 !---------------------
1156 !*** Local variables
1157 !---------------------
1158 !
1159  integer :: i,i_start_data,istat,j,j_start_data,ncid_oro,var_id
1160 !
1161  character(len=150) :: filename,vname
1162 !
1163 !-----------------------------------------------------------------------
1164 !***********************************************************************
1165 !-----------------------------------------------------------------------
1166 !
1167 !-----------------------------------------------------------------------
1168 !*** Get the name of the working directory. Open the topography data
1169 !*** file.
1170 !-----------------------------------------------------------------------
1171 !
1172  filename='INPUT/'//trim(oro_data)
1173 
1174  if (is_master()) then
1175  write(*,23421)trim(filename)
1176 23421 format(' topo filename=',a)
1177  endif
1178 !
1179  call check(nf90_open(filename,nf90_nowrite,ncid_oro))
1180 !
1181 !-----------------------------------------------------------------------
1182 !*** Read in the data including the extra outer row.
1183 !-----------------------------------------------------------------------
1184 !
1185  i_start_data=isd+nhalo_model
1186  j_start_data=jsd+nhalo_model
1187 !
1188  vname='orog_filt'
1189  call check(nf90_inq_varid(ncid_oro,vname,var_id))
1190  call check(nf90_get_var(ncid_oro,var_id &
1191  ,phis_reg(isd-1:ied+1,jsd-1:jed+1) &
1192  ,start=(/i_start_data,j_start_data/)))
1193 !
1194  call check(nf90_close(ncid_oro))
1195 !
1196 !-----------------------------------------------------------------------
1197 !*** We want the geopotential.
1198 !-----------------------------------------------------------------------
1199 !
1200  do j=jsd-1,jed+1
1201  do i=isd-1,ied+1
1202  phis_reg(i,j)=phis_reg(i,j)*grav
1203  enddo
1204  enddo
1205 !
1206 !-----------------------------------------------------------------------
1207 !***********************************************************************
1208 !-----------------------------------------------------------------------
1209 !
1210  end subroutine read_regional_filtered_topo
1211 !
1212 !-----------------------------------------------------------------------
1213 !
1214  end subroutine setup_regional_bc
1215 !
1216 !-----------------------------------------------------------------------
1217 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1218 !-----------------------------------------------------------------------
1219 !
1220  subroutine start_regional_cold_start(Atm, dt_atmos, ak, bk, levp &
1221  ,is ,ie ,js ,je &
1222  ,isd,ied,jsd,jed )
1223 !
1224 !-----------------------------------------------------------------------
1225 !*** Prepare the regional run for a cold start.
1226 !-----------------------------------------------------------------------
1227  implicit none
1228 !-----------------------------------------------------------------------
1229 !
1230 !------------------------
1231 !*** Argument variables
1232 !------------------------
1233 !
1234  type(fv_atmos_type),intent(inout) :: Atm
1235 !
1236  integer ,intent(in) :: is ,ie ,js ,je & !<-- Integration limits of task subdomain
1237  ,isd,ied,jsd,jed & !<-- Memory limits of task subdomain
1238  ,levp
1239 !
1240  real,intent(in) :: dt_atmos
1241  real,intent(in) :: ak(1:levp+1), bk(1:levp+1)
1242 !
1243 !---------------------
1244 !*** Local variables
1245 !---------------------
1246 !
1247  integer :: k
1248 !
1249 !-----------------------------------------------------------------------
1250 !***********************************************************************
1251 !-----------------------------------------------------------------------
1252 !
1253 !-----------------------------------------------------------------------
1254 !*** Get the source of the input data
1255 !-----------------------------------------------------------------------
1256 !
1257  call get_data_source(data_source,atm%flagstruct%regional)
1258 !
1259  call setup_regional_bc(atm, dt_atmos &
1260  ,isd, ied, jsd, jed &
1261  ,atm%npx, atm%npy )
1262 !
1263  bc_hour=0
1264  call regional_bc_data(atm, bc_hour &
1265  ,is, ie, js, je &
1266  ,isd, ied, jsd, jed &
1267  ,ak, bk )
1268  call regional_bc_t1_to_t0(bc_t1, bc_t0 & !
1269  ,atm%npz &
1270  ,ntracers &
1271  ,atm%regional_bc_bounds ) !
1272 !
1273  bc_hour=bc_hour+bc_update_interval
1274 !
1275 !-----------------------------------------------------------------------
1276 !*** If this is a DA run and the first BC file was updated by
1277 !*** the GSI then reset the gsi flag so that all subsequent
1278 !*** BC files are read normally.
1279 !-----------------------------------------------------------------------
1280 !
1281  if(atm%flagstruct%regional_bcs_from_gsi)then
1282  atm%flagstruct%regional_bcs_from_gsi=.false.
1283  endif
1284 !
1285  call regional_bc_data(atm, bc_hour &
1286  ,is, ie, js, je & ! from the 2nd time level
1287  ,isd, ied, jsd, jed & ! in the BC file.
1288  ,ak, bk ) !
1289 !
1290  allocate (ak_in(1:levp+1))
1291  allocate (bk_in(1:levp+1)) ! remapping BC updates during the forecast.
1292  do k=1,levp+1
1293  ak_in(k)=ak(k)
1294  bk_in(k)=bk(k)
1295  enddo
1296 !
1297 !-----------------------------------------------------------------------
1298 !*** If the GSI will need restart files that includes the
1299 !*** fields' boundary rows. Those files were already created.
1300 !*** Prepare the objects that hold their variables' names and
1301 !*** values.
1302 !-----------------------------------------------------------------------
1303 !
1304  if(atm%flagstruct%write_restart_with_bcs)then
1305  call prepare_full_fields(atm)
1306  endif
1307 !
1308 !-----------------------------------------------------------------------
1309 !
1310  end subroutine start_regional_cold_start
1311 !
1312 !-----------------------------------------------------------------------
1313 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1314 !-----------------------------------------------------------------------
1315 !
1316  subroutine start_regional_restart(Atm, dt_atmos &
1317  ,isc,iec,jsc,jec &
1318  ,isd,ied,jsd,jed )
1319 !
1320 !-----------------------------------------------------------------------
1321 !*** Prepare the regional forecast for a restart.
1322 !-----------------------------------------------------------------------
1323  implicit none
1324 !-----------------------------------------------------------------------
1325 !
1326 !------------------------
1327 !*** Argument variables
1328 !------------------------
1329 !
1330  type(fv_atmos_type),intent(inout) :: Atm
1331 !
1332  real,intent(in) :: dt_atmos
1333 !
1334  integer ,intent(in) :: isc,iec,jsc,jec & !<-- Integration limits of task subdomain
1335  ,isd,ied,jsd,jed
1336 !
1337 !---------------------
1338 !*** Local variables
1339 !---------------------
1340 !
1341  integer :: ierr, ios
1342  real, allocatable :: wk2(:,:)
1343 !
1344  logical :: filtered_terrain = .true.
1345  logical :: gfs_dwinds = .true.
1346  integer :: levp = 64
1347  logical :: checker_tr = .false.
1348  integer :: nt_checker = 0
1349  namelist /external_ic_nml/ filtered_terrain, levp, gfs_dwinds &
1350  ,checker_tr, nt_checker
1351 !-----------------------------------------------------------------------
1352 !***********************************************************************
1353 !-----------------------------------------------------------------------
1354 !
1355 !-----------------------------------------------------------------------
1356 !*** Read the number of model layers in the external forecast (=levp).
1357 !-----------------------------------------------------------------------
1358 !
1359  read (input_nml_file,external_ic_nml,iostat=ios)
1360  ierr = check_nml_error(ios,'external_ic_nml')
1361  if(ierr/=0)then
1362  write(0,11011)ierr
1363 11011 format(' start_regional_restart failed to read external_ic_nml ierr=',i3)
1364  endif
1365 !
1366 !-----------------------------------------------------------------------
1367 !*** Get the source of the input data.
1368 !-----------------------------------------------------------------------
1369 !
1370  call get_data_source(data_source,atm%flagstruct%regional)
1371 !
1372 !-----------------------------------------------------------------------
1373 !*** Preliminary setup for the forecast.
1374 !-----------------------------------------------------------------------
1375 !
1376  call setup_regional_bc(atm, dt_atmos &
1377  ,isd, ied, jsd, jed &
1378  ,atm%npx, atm%npy )
1379 !
1380  allocate (wk2(levp+1,2))
1381  allocate (ak_in(levp+1))
1382  allocate (bk_in(levp+1)) ! remapping BC updates during the forecast.
1383  call read_data('INPUT/gfs_ctrl.nc','vcoord',wk2, no_domain=.true.)
1384  ak_in(1:levp+1) = wk2(1:levp+1,1)
1385  ak_in(1) = 1.e-9
1386  bk_in(1:levp+1) = wk2(1:levp+1,2)
1387  deallocate(wk2)
1388  bc_hour=nint(current_time_in_seconds/3600.)
1389 !
1390 !-----------------------------------------------------------------------
1391 !*** Fill time level t1 from the BC file at the restart time.
1392 !-----------------------------------------------------------------------
1393 !
1394  call regional_bc_data(atm, bc_hour &
1395  ,isc, iec, jsc, jec &
1396  ,isd, ied, jsd, jed &
1397  ,ak_in, bk_in )
1398 !
1399 !-----------------------------------------------------------------------
1400 !*** If this is a DA run and the first BC file was updated by
1401 !*** the GSI then that file was read differently in the preceding
1402 !*** call to subroutine regional_bc_data. Now reset the gsi
1403 !*** flag so that all subsequent BC files are read normally.
1404 !-----------------------------------------------------------------------
1405 !
1406  if(atm%flagstruct%regional_bcs_from_gsi)then
1407  atm%flagstruct%regional_bcs_from_gsi=.false.
1408  endif
1409 !
1410 !-----------------------------------------------------------------------
1411 !*** If the GSI will need restart files that include the
1412 !*** fields' boundary rows after this forecast or forecast
1413 !*** segment completes then prepare the objects that will
1414 !*** hold their variables' names and values.
1415 !-----------------------------------------------------------------------
1416 !
1417  if(atm%flagstruct%write_restart_with_bcs)then
1418  call prepare_full_fields(atm)
1419  endif
1420 !
1421 !-----------------------------------------------------------------------
1422 !
1423  end subroutine start_regional_restart
1424 !
1425 !-----------------------------------------------------------------------
1426 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1427 !-----------------------------------------------------------------------
1428 !
1429  subroutine read_new_bc_data(Atm, Time, Time_step_atmos, p_split &
1430  ,isd,ied,jsd,jed )
1431 !
1432 !-----------------------------------------------------------------------
1433 !*** When it is time to read new boundary data from the external files
1434 !*** move time level t1 to t0 and then read the data into t1.
1435 !-----------------------------------------------------------------------
1436  implicit none
1437 !-----------------------------------------------------------------------
1438 !
1439 !------------------------
1440 !*** Argument variables
1441 !------------------------
1442 !
1443  type(fv_atmos_type),intent(inout) :: Atm
1444  type(time_type),intent(in) :: Time
1445  type(time_type),intent(in) :: time_step_atmos
1446 
1447  integer,intent(in) :: isd,ied,jsd,jed & !<-- Memory limits of task subdomain
1448  ,p_split
1449 !
1450 !---------------------
1451 !*** Local variables
1452 !---------------------
1453 !
1454  integer :: atmos_time_step, sec
1455  real :: dt_atmos
1456  type(time_type) :: atmos_time
1457 !
1458 !-----------------------------------------------------------------------
1459 !***********************************************************************
1460 !-----------------------------------------------------------------------
1461 !
1462  atmos_time = time - atm%Time_init
1463  atmos_time_step = atmos_time / time_step_atmos
1464  current_time_in_seconds = time_type_to_real( atmos_time )
1465  if (mpp_pe() == 0 .and. atm%flagstruct%fv_debug) write(*,"('current_time_seconds = ',f9.1)")current_time_in_seconds
1466 !
1467  call get_time (time_step_atmos, sec)
1468  dt_atmos = real(sec)
1469 !
1470  if(atmos_time_step==0.or.atm%flagstruct%warm_start)then
1471  ntimesteps_per_bc_update=nint(bc_update_interval*3600./(dt_atmos/real(abs(p_split))))
1472  endif
1473 !
1474  if(atmos_time_step+1>=ntimesteps_per_bc_update.and.mod(atmos_time_step,ntimesteps_per_bc_update)==0 &
1475  .or. &
1476  atm%flagstruct%warm_start.and.begin_regional_restart)then
1477 !
1478  begin_regional_restart=.false.
1479  bc_hour=bc_hour+bc_update_interval
1480 !
1481 !-----------------------------------------------------------------------
1482 !*** Transfer the time level t1 data to t0.
1483 !-----------------------------------------------------------------------
1484 !
1485  call regional_bc_t1_to_t0(bc_t1, bc_t0 &
1486  ,atm%npz &
1487  ,ntracers &
1488  ,atm%regional_bc_bounds )
1489 !
1490 !-----------------------------------------------------------------------
1491 !*** Fill time level t1 from the BC file containing data from
1492 !*** the next time level.
1493 !-----------------------------------------------------------------------
1494 !
1495  call regional_bc_data(atm, bc_hour &
1496  ,atm%bd%is, atm%bd%ie &
1497  ,atm%bd%js, atm%bd%je &
1498  ,isd, ied, jsd, jed &
1499  ,ak_in, bk_in )
1500  endif
1501 !
1502 !-----------------------------------------------------------------------
1503 !
1504  end subroutine read_new_bc_data
1505 !
1506 !-----------------------------------------------------------------------
1507 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1508 !-----------------------------------------------------------------------
1509 !
1510  subroutine regional_bc_data(Atm,bc_hour &
1511  ,is,ie,js,je &
1512  ,isd,ied,jsd,jed &
1513  ,ak,bk )
1514 !
1515 !-----------------------------------------------------------------------
1516 !*** Regional boundary data is obtained from the external BC file.
1517 !-----------------------------------------------------------------------
1518 !-----------------------------------------------------------------------
1519  implicit none
1520 !-----------------------------------------------------------------------
1521 !
1522 !------------------------
1523 !*** Argument variables
1524 !------------------------
1525 !
1526 !-----------
1527 !*** Input
1528 !-----------
1529 !
1530  integer,intent(in) :: bc_hour
1531 !
1532  integer,intent(in) :: is,ie,js,je & !<-- Compute limits of task subdomain
1533  ,isd,ied,jsd,jed
1534 !
1535  real,dimension(:),intent(in) :: ak,bk
1536 !
1537 !-----------------
1538 !*** Input/output
1539 !-----------------
1540 !
1541  type(fv_atmos_type),target,intent(inout) :: Atm
1542 !
1543 !---------------------
1544 !*** Local variables
1545 !---------------------
1546 !
1547  integer :: dimid,i,j,k,klev_in,klev_out,n,nlev
1548 !
1549  integer :: is_north,is_south,is_east,is_west &
1550  ,ie_north,ie_south,ie_east,ie_west &
1551  ,js_north,js_south,js_east,js_west &
1552  ,je_north,je_south,je_east,je_west
1553 !
1554  integer :: is_u,ie_u,js_u,je_u &
1555  ,is_v,ie_v,js_v,je_v
1556 !
1557  integer :: is_input,ie_input,js_input,je_input
1558 !
1559  integer :: i_start,i_end,j_start,j_end
1560 !
1561  integer :: nside,nt,index
1562 !
1563  real,dimension(:,:,:),allocatable :: ud,vd,uc,vc
1564 !
1565  real,dimension(:,:),allocatable :: ps_reg
1566  real,dimension(:,:,:),allocatable :: delp_input,delz_input &
1567  ,ps_input,t_input &
1568  ,w_input,zh_input
1569  real,dimension(:,:,:),allocatable :: u_s_input,v_s_input &
1570  ,u_w_input,v_w_input
1571  real,dimension(:,:,:,:),allocatable :: tracers_input
1572 !
1573  real(kind=R_GRID), dimension(2):: p1, p2, p3, p4
1574  real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
1575 
1576 #undef USE_FMS_READ
1577 #ifdef USE_FMS_READ
1578  integer :: isc2, iec2, jsc2, jec2
1579  real(kind=R_GRID), allocatable, dimension(:,:) :: tmpx, tmpy
1580  integer :: start(4), nread(4)
1581  real(kind=R_GRID), allocatable, dimension(:,:,:) :: reg_grid
1582  real(kind=R_GRID), allocatable, dimension(:,:,:) :: reg_agrid
1583 #endif
1584 !
1585  logical,save :: computed_regional_bc_indices=.false.
1586 !
1587  character(len=3) :: int_to_char
1588  character(len=5) :: side
1589  character(len=6) :: fmt='(i3.3)'
1590 !
1591  character(len=50) :: file_name
1592 !
1593  integer,save :: kount1=0,kount2=0
1594  integer :: istart, iend, jstart, jend
1595  integer :: npx, npy
1596 !
1597  character(len=60) :: var_name_root
1598  logical :: required
1599 !
1600  logical :: call_remap
1601 !
1602 !-----------------------------------------------------------------------
1603 !***********************************************************************
1604 !-----------------------------------------------------------------------
1605 !
1606 !-----------------------------------------------------------------------
1607 !*** Only boundary tasks are needed.
1608 !-----------------------------------------------------------------------
1609 !
1610  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
1611  return
1612  endif
1613 !
1614 !-----------------------------------------------------------------------
1615 !
1616  klev_out=atm%npz
1617 !
1618 !-----------------------------------------------------------------------
1619 !*** Construct the name of the regional BC file to be read.
1620 !*** We must know whether this is a standard BC file from chgres
1621 !*** or a new BC file generated from DA-updated data from enlarged
1622 !*** restart files that include the boundary rows.
1623 !-----------------------------------------------------------------------
1624 !
1625  write(int_to_char,fmt) bc_hour
1626  if(.not.atm%flagstruct%regional_bcs_from_gsi)then
1627  file_name='INPUT/gfs_bndy.tile7.'//int_to_char//'.nc'
1628  else
1629  file_name='INPUT/gfs_bndy.tile7.'//int_to_char//'_gsi.nc'
1630  endif
1631 !
1632  if (is_master()) then
1633  write(*,22211)trim(file_name)
1634 22211 format(' regional_bc_data file_name=',a)
1635  endif
1636 !-----------------------------------------------------------------------
1637 !*** Open the regional BC file.
1638 !*** Find the # of layers (klev_in) in the BC input.
1639 !-----------------------------------------------------------------------
1640 !
1641  call check(nf90_open(file_name,nf90_nowrite,ncid))
1642  if (is_master()) then
1643  write(0,*)' opened BC file ',trim(file_name)
1644  endif
1645 !
1646  call check(nf90_inq_dimid(ncid,'lev',dimid))
1647  call check(nf90_inquire_dimension(ncid,dimid,len=klev_in))
1648 !
1649 !-----------------------------------------------------------------------
1650 !*** Allocate the boundary variables and initialize them to garbage.
1651 !-----------------------------------------------------------------------
1652 !
1653  is_input=is-nhalo_data
1654  ie_input=ie+nhalo_data
1655  js_input=js-nhalo_data
1656  je_input=je+nhalo_data
1657  npx = atm%npx
1658  npy = atm%npy
1659 !
1660  allocate( ps_input(is_input:ie_input,js_input:je_input,1)) ; ps_input=real_snan
1661  allocate( t_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; t_input=real_snan
1662  allocate( w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; w_input=real_snan
1663  allocate(u_s_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; u_s_input=real_snan
1664  allocate(v_s_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; v_s_input=real_snan
1665  allocate(u_w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; u_w_input=real_snan
1666  allocate(v_w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; v_w_input=real_snan
1667 !
1668  if(atm%flagstruct%regional_bcs_from_gsi)then
1669  allocate(delp_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; delp_input=real_snan
1670  allocate(delz_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; delz_input=real_snan
1671  else
1672  allocate( zh_input(is_input:ie_input,js_input:je_input,1:klev_in+1)) ; zh_input=real_snan
1673  endif
1674 !
1675  allocate(tracers_input(is_input:ie_input,js_input:je_input,klev_in,ntracers)) ; tracers_input=real_snan
1676 !
1677 !-----------------------------------------------------------------------
1678 !*** Extract each variable from the regional BC file. The final
1679 !*** argument is the object being filled.
1680 !-----------------------------------------------------------------------
1681 !
1682 !------------------
1683 !*** Sfc pressure
1684 !------------------
1685 !
1686  nlev=1
1687  var_name_root='ps'
1688  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1689  ,nlev &
1690  ,ntracers &
1691 ! ,Atm%regional_bc_bounds &
1692  ,var_name_root &
1693  ,array_3d=ps_input )
1694 !
1695 !-----------------------
1696 !*** Vertical velocity
1697 !-----------------------
1698 !
1699  nlev=klev_in
1700  var_name_root='w'
1701  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1702  ,nlev &
1703  ,ntracers &
1704 ! ,Atm%regional_bc_bounds &
1705  ,var_name_root &
1706  ,array_3d=w_input)
1707 !
1708 !-----------------------
1709 !*** Interface heights
1710 !-----------------------
1711 !
1712  if(.not.atm%flagstruct%regional_bcs_from_gsi)then
1713  nlev=klev_in+1
1714  var_name_root='zh'
1715  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1716  ,nlev &
1717  ,ntracers &
1718  ,var_name_root &
1719  ,array_3d=zh_input)
1720  endif
1721 !
1722 !--------------------------
1723 !*** Sensible temperature
1724 !--------------------------
1725 !
1726  if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
1727  nlev=klev_in
1728  var_name_root='t'
1729  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1730  ,nlev &
1731  ,ntracers &
1732 ! ,Atm%regional_bc_bounds &
1733  ,var_name_root &
1734  ,array_3d=t_input)
1735  endif
1736 !
1737 !-----------------------------
1738 !*** U component south/north
1739 !-----------------------------
1740 !
1741  nlev=klev_in
1742  var_name_root='u_s'
1743  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1744  ,nlev &
1745  ,ntracers &
1746 ! ,Atm%regional_bc_bounds &
1747  ,var_name_root &
1748  ,array_3d=u_s_input)
1749 !
1750 !-----------------------------
1751 !*** V component south/north
1752 !-----------------------------
1753 !
1754  nlev=klev_in
1755  var_name_root='v_s'
1756  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1757  ,nlev &
1758  ,ntracers &
1759 ! ,Atm%regional_bc_bounds &
1760  ,var_name_root &
1761  ,array_3d=v_s_input)
1762 !
1763 !---------------------------
1764 !*** U component east/west
1765 !---------------------------
1766 !
1767  nlev=klev_in
1768  var_name_root='u_w'
1769  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1770  ,nlev &
1771  ,ntracers &
1772 ! ,Atm%regional_bc_bounds &
1773  ,var_name_root &
1774  ,array_3d=u_w_input)
1775 !
1776 !---------------------------
1777 !*** V component east/west
1778 !---------------------------
1779 !
1780  nlev=klev_in
1781  var_name_root='v_w'
1782  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1783  ,nlev &
1784  ,ntracers &
1785 ! ,Atm%regional_bc_bounds &
1786  ,var_name_root &
1787  ,array_3d=v_w_input)
1788 !
1789 !-----------------------------------------------------------------------
1790 !*** If this is a DA-updated BC file then also read in the layer
1791 !*** pressure depths.
1792 !-----------------------------------------------------------------------
1793 !
1794  if(atm%flagstruct%regional_bcs_from_gsi)then
1795  nlev=klev_in
1796  var_name_root='delp'
1797  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1798  ,nlev &
1799  ,ntracers &
1800  ,var_name_root &
1801  ,array_3d=delp_input)
1802  var_name_root='delz'
1803  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1804  ,nlev &
1805  ,ntracers &
1806  ,var_name_root &
1807  ,array_3d=delz_input)
1808  endif
1809 !
1810 !-------------
1811 !*** Tracers
1812 !-------------
1813 
1814  nlev=klev_in
1815 !
1816 !-----------------------------------------------------------------------
1817 !*** Read the tracers specified in the field_table. If they are not
1818 !*** in the input data then print a warning and set them to 0 in the
1819 !*** boundary. Some tracers are mandatory to have, because they are
1820 !*** used later for calculating virtual potential temperature etc.
1821 !-----------------------------------------------------------------------
1822 !
1823  do nt = 1, ntracers
1824  call get_tracer_names(model_atmos, nt, var_name_root)
1825  index= get_tracer_index(model_atmos,trim(var_name_root))
1826  if (index==liq_water_index .or. index==sphum_index) then
1827  required = .true.
1828  else
1829  required = .false.
1830  endif
1831  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1832  ,nlev &
1833  ,ntracers &
1834  ,var_name_root &
1835  ,array_4d=tracers_input &
1836  ,tlev=index &
1837  ,required=required )
1838  enddo
1839 !
1840 !-----------------------------------------------------------------------
1841 !*** For a DA-updated BC file we can simply transfer the data
1842 !*** from the *_input arrays into the model's boundary arrays
1843 !*** since they came out of restart files. Otherwise proceed
1844 !*** with vertical remapping from input layers to model forecast
1845 !*** layers and rotate the winds from geographic lat/lon to the
1846 !*** integration grid.
1847 !-----------------------------------------------------------------------
1848 !
1849  data_to_bc: if(atm%flagstruct%regional_bcs_from_gsi)then
1850 !
1851 !-----------------------------------------------------------------------
1852 !
1853  call fill_bc_for_da
1854 !
1855 !-----------------------------------------------------------------------
1856 !
1857  else
1858 !
1859 !-----------------------------------------------------------------------
1860 !*** One final array needs to be allocated. It is the sfc pressure
1861 !*** in the domain's boundary region that is derived from the input
1862 !*** sfc pressure from the BC files. The derived sfc pressure will
1863 !*** be needed in the vertical remapping of the wind components to
1864 !*** the integration levels.
1865 !-----------------------------------------------------------------------
1866 !
1867  allocate(ps_reg(is_input:ie_input,js_input:je_input)) ; ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed
1868 !
1869 !-----------------------------------------------------------------------
1870 !*** We have the boundary variables from the BC file on the levels
1871 !*** of the input data. Remap the scalars (tracers, vertical
1872 !*** velocity, ozone) to the FV3 domain levels. Scalar remapping
1873 !*** must be done on all four sides before remapping of the winds
1874 !*** since pressures are needed on each side of wind points and so
1875 !*** for a given wind component those pressures could include values
1876 !*** from two different boundary side regions.
1877 !-----------------------------------------------------------------------
1878 !
1879 ! Definitions in this module greatly differ from those in existing nesting
1880 ! code or elsewhere in FMS. North <--> South, East <--> West, and
1881 ! North and South always span [isd-1 , ied+1] while East and West do not
1882 ! go into the outermost corners (so the they span [1, je], always.)
1883 !-----------------------------------------------------------------------
1884  sides_scalars: do nside=1,4
1885 !-----------------------------------------------------------------------
1886 !-----------
1887 !*** North
1888 !-----------
1889 !
1890  call_remap=.false.
1891 !
1892  if(nside==1)then
1893  if(north_bc)then
1894  call_remap=.true.
1895  side='north'
1896  bc_side_t1=>bc_t1%north
1897  bc_side_t0=>bc_t0%north
1898  endif
1899  endif
1900 !
1901  if(nside==2)then
1902  if(south_bc)then
1903  call_remap=.true.
1904  side='south'
1905  bc_side_t1=>bc_t1%south
1906  bc_side_t0=>bc_t0%south
1907  endif
1908  endif
1909 !
1910  if(nside==3)then
1911  if(east_bc)then
1912  call_remap=.true.
1913  side='east '
1914  bc_side_t1=>bc_t1%east
1915  bc_side_t0=>bc_t0%east
1916  endif
1917  endif
1918 !
1919  if(nside==4)then
1920  if(west_bc)then
1921  call_remap=.true.
1922  side='west '
1923  bc_side_t1=>bc_t1%west
1924  bc_side_t0=>bc_t0%west
1925  endif
1926  endif
1927 !
1928  if(call_remap)then
1929  call remap_scalar_nggps_regional_bc(atm &
1930  ,side &
1931 
1932  ,isd,ied,jsd,jed &
1933 
1934  ,is_input &
1935  ,ie_input & ! Input array
1936  ,js_input & ! index limits.
1937  ,je_input &
1938 
1939  ,klev_in, klev_out &
1940  ,ntracers &
1941  ,ak, bk &
1942 
1943  ,ps_input &
1944  ,t_input & ! BC vbls
1945  ,tracers_input & ! on input
1946  ,w_input & ! model levels
1947  ,zh_input &
1948 
1949  ,phis_reg &
1950 
1951  ,ps_reg &
1952 
1953  ,bc_side_t1 )
1954 !
1955  call set_delp_and_tracers(bc_side_t1,atm%npz,atm%flagstruct%nwat)
1956 !
1957  if(nside==1)then
1958  if(north_bc)then
1959  if (is == 1) then
1960  istart = 1
1961  else
1962  istart = isd
1963  endif
1964  if (ie == npx-1) then
1965  iend = npx-1
1966  else
1967  iend = ied
1968  endif
1969 
1970  do k=1,npz
1971  do j=jsd,0
1972  do i=istart,iend
1973  delz_regbc%south_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
1974  delz_regbc%south_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
1975  enddo
1976  enddo
1977  enddo
1978 
1979  ! North, south include all corners
1980  if (is == 1) then
1981  do k=1,npz
1982  do j=jsd,0
1983  do i=isd,0
1984  delz_regbc%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
1985  delz_regbc%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
1986  enddo
1987  enddo
1988  enddo
1989  endif
1990 
1991  if (ie == npx-1) then
1992  do k=1,npz
1993  do j=jsd,0
1994  do i=npx,ied
1995  delz_regbc%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
1996  delz_regbc%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
1997  enddo
1998  enddo
1999  enddo
2000  endif
2001  endif
2002  endif
2003 
2004  if(nside==2)then
2005  if(south_bc)then
2006  if (is == 1) then
2007  istart = 1
2008  else
2009  istart = isd
2010  endif
2011  if (ie == npx-1) then
2012  iend = npx-1
2013  else
2014  iend = ied
2015  endif
2016 
2017  do k=1,npz
2018  do j=npy,jed
2019  do i=istart,iend
2020  delz_regbc%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
2021  delz_regbc%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
2022  enddo
2023  enddo
2024  enddo
2025 
2026  ! North, south include all corners
2027  if (is == 1) then
2028  do k=1,npz
2029  do j=npy,jed
2030  do i=isd,0
2031  delz_regbc%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
2032  delz_regbc%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
2033  enddo
2034  enddo
2035  enddo
2036  endif
2037 
2038 
2039  if (ie == npx-1) then
2040  do k=1,npz
2041  do j=npy,jed
2042  do i=npx,ied
2043  delz_regbc%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
2044  delz_regbc%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
2045  enddo
2046  enddo
2047  enddo
2048  endif
2049  endif
2050  endif
2051 
2052 !
2053 
2054  if(nside==3)then
2055  if(east_bc)then
2056  if (js == 1) then
2057  jstart = 1
2058  else
2059  jstart = jsd
2060  endif
2061  if (je == npy-1) then
2062  jend = je
2063  else
2064  jend = jed
2065  endif
2066 
2067 
2068  do k=1,npz
2069  do j=jstart,jend
2070  do i=isd,0
2071  delz_regbc%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
2072  delz_regbc%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
2073  enddo
2074  enddo
2075  enddo
2076  endif
2077  endif
2078 
2079  if(nside==4)then
2080  if(west_bc)then
2081  if (js == 1) then
2082  jstart = 1
2083  else
2084  jstart = jsd
2085  endif
2086  if (je == npy-1) then
2087  jend = je
2088  else
2089  jend = jed
2090  endif
2091 
2092 
2093  do k=1,npz
2094  do j=jstart,jend
2095  do i=npx,ied
2096  delz_regbc%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k)
2097  delz_regbc%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k)
2098  enddo
2099  enddo
2100  enddo
2101  endif
2102  endif
2103 
2104  endif
2105 !
2106 !-----------------------------------------------------------------------
2107  enddo sides_scalars
2108 !-----------------------------------------------------------------------
2109 !
2110 !-----------------------------------------------------------------------
2111 !*** Now that we have the pressure throughout the boundary region
2112 !*** including a row beyond the boundary winds we are ready to
2113 !*** finalize those winds.
2114 !-----------------------------------------------------------------------
2115 !
2116 !-----------------------------------------------------------------------
2117 !*** Transform the D-grid wind components on the four sides of
2118 !*** the regional domain then remap them from the input levels
2119 !*** to the integration levels.
2120 !-----------------------------------------------------------------------
2121 !
2122 #ifdef USE_FMS_READ
2123  isc2 = 2*(isd-1+nhalo_data)-1
2124  iec2 = 2*(ied+2+nhalo_data)-1
2125  jsc2 = 2*(jsd-1+nhalo_data)-1
2126  jec2 = 2*(jed+2+nhalo_data)-1
2127  allocate(tmpx(isc2:iec2, jsc2:jec2)) ; tmpx=dbl_snan
2128  allocate(tmpy(isc2:iec2, jsc2:jec2)) ; tmpy=dbl_snan
2129  start = 1; nread = 1
2130  start(1) = isc2; nread(1) = iec2 - isc2 + 1
2131  start(2) = jsc2; nread(2) = jec2 - jsc2 + 1
2132  call read_data("INPUT/grid.tile7.halo4.nc", 'x', tmpx, start, nread, no_domain=.true.)
2133  call read_data("INPUT/grid.tile7.halo4.nc", 'y', tmpy, start, nread, no_domain=.true.)
2134 
2135  allocate(reg_grid(isd-1:ied+2,jsd-1:jed+2,1:2)) ; reg_grid=dbl_snan
2136  do j = jsd-1, jed+2
2137  do i = isd-1, ied+2
2138  reg_grid(i,j,1) = tmpx(2*(i+nhalo_data)-1, 2*(j+nhalo_data)-1)*pi/180.
2139  reg_grid(i,j,2) = tmpy(2*(i+nhalo_data)-1, 2*(j+nhalo_data)-1)*pi/180.
2140  if ( reg_grid(i,j,1) /= grid_reg(i,j,1) ) then
2141  write(0,*)' reg_grid(i,j,1) /= grid_reg(i,j,1) ',i,j, reg_grid(i,j,1),grid_reg(i,j,1)
2142  endif
2143  enddo
2144  enddo
2145 
2146  allocate(reg_agrid(isd-1:ied+1,jsd-1:jed+1,1:2)) ; reg_agrid=dbl_snan
2147  do j=jsd-1,jed+1
2148  do i=isd-1,ied+1
2149  call cell_center2(reg_grid(i,j, 1:2), reg_grid(i+1,j, 1:2), &
2150  reg_grid(i,j+1,1:2), reg_grid(i+1,j+1,1:2), &
2151  reg_agrid(i,j,1:2) )
2152  enddo
2153  enddo
2154 #endif
2155 !
2156 !-----------------------------------------------------------------------
2157 !*** Loop through the four sides of the domain.
2158 !-----------------------------------------------------------------------
2159 !
2160 !-----------------------------------------------------------------------
2161  sides_winds: do nside=1,4
2162 !-----------------------------------------------------------------------
2163 !
2164  call_remap=.false.
2165 
2166  if(nside==1)then
2167  if(north_bc)then
2168  call_remap=.true.
2169  bc_side_t1=>bc_t1%north
2170 !
2171  is_u=atm%regional_bc_bounds%is_north_uvs
2172  ie_u=atm%regional_bc_bounds%ie_north_uvs
2173  js_u=atm%regional_bc_bounds%js_north_uvs
2174  je_u=atm%regional_bc_bounds%je_north_uvs
2175 !
2176  is_v=atm%regional_bc_bounds%is_north_uvw
2177  ie_v=atm%regional_bc_bounds%ie_north_uvw
2178  js_v=atm%regional_bc_bounds%js_north_uvw
2179  je_v=atm%regional_bc_bounds%je_north_uvw
2180  endif
2181  endif
2182 !
2183  if(nside==2)then
2184  if(south_bc)then
2185  call_remap=.true.
2186  bc_side_t1=>bc_t1%south
2187 !
2188  is_u=atm%regional_bc_bounds%is_south_uvs
2189  ie_u=atm%regional_bc_bounds%ie_south_uvs
2190  js_u=atm%regional_bc_bounds%js_south_uvs
2191  je_u=atm%regional_bc_bounds%je_south_uvs
2192 !
2193  is_v=atm%regional_bc_bounds%is_south_uvw
2194  ie_v=atm%regional_bc_bounds%ie_south_uvw
2195  js_v=atm%regional_bc_bounds%js_south_uvw
2196  je_v=atm%regional_bc_bounds%je_south_uvw
2197  endif
2198  endif
2199 !
2200  if(nside==3)then
2201  if(east_bc)then
2202  call_remap=.true.
2203  bc_side_t1=>bc_t1%east
2204 !
2205  is_u=atm%regional_bc_bounds%is_east_uvs
2206  ie_u=atm%regional_bc_bounds%ie_east_uvs
2207  js_u=atm%regional_bc_bounds%js_east_uvs
2208  je_u=atm%regional_bc_bounds%je_east_uvs
2209 !
2210  is_v=atm%regional_bc_bounds%is_east_uvw
2211  ie_v=atm%regional_bc_bounds%ie_east_uvw
2212  js_v=atm%regional_bc_bounds%js_east_uvw
2213  je_v=atm%regional_bc_bounds%je_east_uvw
2214  endif
2215  endif
2216 !
2217  if(nside==4)then
2218  if(west_bc)then
2219  call_remap=.true.
2220  bc_side_t1=>bc_t1%west
2221 !
2222  is_u=atm%regional_bc_bounds%is_west_uvs
2223  ie_u=atm%regional_bc_bounds%ie_west_uvs
2224  js_u=atm%regional_bc_bounds%js_west_uvs
2225  je_u=atm%regional_bc_bounds%je_west_uvs
2226 !
2227  is_v=atm%regional_bc_bounds%is_west_uvw
2228  ie_v=atm%regional_bc_bounds%ie_west_uvw
2229  js_v=atm%regional_bc_bounds%js_west_uvw
2230  je_v=atm%regional_bc_bounds%je_west_uvw
2231  endif
2232  endif
2233 !
2234  if(call_remap)then
2235 !
2236  allocate(ud(is_u:ie_u,js_u:je_u,1:nlev)) ; ud=real_snan
2237  allocate(vd(is_v:ie_v,js_v:je_v,1:nlev)) ; vd=real_snan
2238  allocate(vc(is_u:ie_u,js_u:je_u,1:nlev)) ; vc=real_snan
2239  allocate(uc(is_v:ie_v,js_v:je_v,1:nlev)) ; uc=real_snan
2240 !
2241  do k=1,nlev
2242  do j=js_u,je_u
2243  do i=is_u,ie_u
2244  p1(:) = grid_reg(i, j,1:2)
2245  p2(:) = grid_reg(i+1,j,1:2)
2246  call mid_pt_sphere(p1, p2, p3)
2247  call get_unit_vect2(p1, p2, e1)
2248  call get_latlon_vector(p3, ex, ey)
2249  ud(i,j,k) = u_s_input(i,j,k)*inner_prod(e1,ex)+v_s_input(i,j,k)*inner_prod(e1,ey)
2250  p4(:) = agrid_reg(i,j,1:2) ! cell centroid
2251  call get_unit_vect2(p3, p4, e2) !C-grid V-wind unit vector
2252  vc(i,j,k) = u_s_input(i,j,k)*inner_prod(e2,ex)+v_s_input(i,j,k)*inner_prod(e2,ey)
2253  enddo
2254  enddo
2255 !
2256  do j=js_v,je_v
2257  do i=is_v,ie_v
2258  p1(:) = grid_reg(i,j ,1:2)
2259  p2(:) = grid_reg(i,j+1,1:2)
2260  call mid_pt_sphere(p1, p2, p3)
2261  call get_unit_vect2(p1, p2, e2)
2262  call get_latlon_vector(p3, ex, ey)
2263  vd(i,j,k) = u_w_input(i,j,k)*inner_prod(e2,ex)+v_w_input(i,j,k)*inner_prod(e2,ey)
2264  p4(:) = agrid_reg(i,j,1:2) ! cell centroid
2265  call get_unit_vect2(p3, p4, e1) !C-grid U-wind unit vector
2266  uc(i,j,k) = u_w_input(i,j,k)*inner_prod(e1,ex)+v_w_input(i,j,k)*inner_prod(e1,ey)
2267  enddo
2268  enddo
2269  enddo
2270 !
2271  call remap_dwinds_regional_bc(atm &
2272 
2273  ,is_input &
2274  ,ie_input & ! Index limits for scalars
2275  ,js_input & ! at center of north BC region grid cells.
2276  ,je_input &
2277 
2278  ,is_u &
2279  ,ie_u & ! Index limits for u component
2280  ,js_u & ! on north edge of BC region grid cells.
2281  ,je_u &
2282 
2283  ,is_v &
2284  ,ie_v & ! Index limits for v component
2285  ,js_v & ! on north edge of BC region grid cells.
2286  ,je_v &
2287 
2288  ,klev_in, klev_out &
2289  ,ak, bk &
2290 
2291  ,ps_reg &
2292  ,ud ,vd &
2293  ,uc ,vc &
2294  ,bc_side_t1 )
2295 !
2296  deallocate(ud,vd,uc,vc)
2297 !
2298  endif
2299 !
2300 !-----------------------------------------------------------------------
2301  enddo sides_winds
2302 !-----------------------------------------------------------------------
2303 !
2304  endif data_to_bc
2305 !
2306 !-----------------------------------------------------------------------
2307 !*** Close the boundary file.
2308 !-----------------------------------------------------------------------
2309 !
2310  call check(nf90_close(ncid))
2311 !
2312 !-----------------------------------------------------------------------
2313 !*** Deallocate working arrays.
2314 !-----------------------------------------------------------------------
2315 !
2316 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2317  if(allocated(ps_input))then
2318  deallocate(ps_input)
2319  endif
2320  if(allocated(t_input))then
2321  deallocate(t_input)
2322  endif
2323  if(allocated(zh_input))then
2324  deallocate(zh_input)
2325  endif
2326  if(allocated(w_input))then
2327  deallocate(w_input)
2328  endif
2329  if(allocated(tracers_input))then
2330  deallocate(tracers_input)
2331  endif
2332  if(allocated(u_s_input))then
2333  deallocate(u_s_input)
2334  endif
2335  if(allocated(u_w_input))then
2336  deallocate(u_w_input)
2337  endif
2338  if(allocated(v_s_input))then
2339  deallocate(v_s_input)
2340  endif
2341  if(allocated(v_w_input))then
2342  deallocate(v_w_input)
2343  endif
2344  if(allocated(delp_input))then
2345  deallocate(delp_input)
2346  endif
2347  if(allocated(delz_input))then
2348  deallocate(delz_input)
2349  endif
2350 !
2351 !-----------------------------------------------------------------------
2352 !*** Fill the remaining boundary arrays starting with the divergence.
2353 !-----------------------------------------------------------------------
2354 !
2355  call fill_divgd_bc
2356 !
2357 !-----------------------------------------------------------------------
2358 !*** Fill the total condensate in the regional boundary array.
2359 !-----------------------------------------------------------------------
2360 !
2361 #ifdef USE_COND
2362  call fill_q_con_bc
2363 #endif
2364 !
2365 !-----------------------------------------------------------------------
2366 !*** Fill moist kappa in the regional domain boundary array.
2367 !-----------------------------------------------------------------------
2368 !
2369 #ifdef MOIST_CAPPA
2370  call fill_cappa_bc
2371 #endif
2372 !
2373 !-----------------------------------------------------------------------
2374 !*** Convert the boundary region sensible temperature array to
2375 !*** FV3's modified virtual potential temperature.
2376 !-----------------------------------------------------------------------
2377 !
2378  call convert_to_virt_pot_temp(isd,ied,jsd,jed,npz)
2379 !
2380 !-----------------------------------------------------------------------
2381 !*** If nudging of the specific humidity has been selected then
2382 !*** nudge the boundary values in the same way as is done for the
2383 !*** interior.
2384 !-----------------------------------------------------------------------
2385 !
2386  if(atm%flagstruct%nudge_qv)then
2387  call nudge_qv_bc(atm,isd,ied,jsd,jed)
2388  endif
2389 !
2390 !-----------------------------------------------------------------------
2391 
2392  contains
2393 
2394 !-----------------------------------------------------------------------
2395 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2396 !-----------------------------------------------------------------------
2397 !
2398  subroutine fill_bc_for_da
2399 !
2400 !-----------------------------------------------------------------------
2401 !*** Transfer the input boundary data directly into the BC object.
2402 !-----------------------------------------------------------------------
2403 !
2404 !---------------------
2405 !*** Local variables
2406 !---------------------
2407 !
2408  integer :: i,j,k,n
2409 !
2410 !-----------------------------------------------------------------------
2411 !***********************************************************************
2412 !-----------------------------------------------------------------------
2413 !
2414 !-----------------------------------------------------------------------
2415 !*** Since corner tasks are on more than one side we cannot
2416 !*** generalize the transfer of data into a given side's
2417 !*** arrays. Do each side separately.
2418 !
2419 !*** Simply obtain the loop limits from the bounds of one of the
2420 !*** BC arrays to be filled.
2421 !-----------------------------------------------------------------------
2422 !
2423 !-----------
2424 !*** North
2425 !-----------
2426 !
2427  if(north_bc)then
2428 !
2429  is_input=lbound(bc_t1%north%delp_BC,1)
2430  ie_input=ubound(bc_t1%north%delp_BC,1) ! Index limits for
2431  js_input=lbound(bc_t1%north%delp_BC,2) ! mass variables.
2432  je_input=ubound(bc_t1%north%delp_BC,2)
2433 !
2434  do k=1,klev_in
2435  do j=js_input,je_input
2436  do i=is_input,ie_input
2437  bc_t1%north%delp_BC(i,j,k)=delp_input(i,j,k)
2438  bc_t1%north%pt_BC(i,j,k)=t_input(i,j,k)
2439  bc_t1%north%w_BC(i,j,k)=w_input(i,j,k)
2440  bc_t1%north%delz_BC(i,j,k)=delz_input(i,j,k)
2441  enddo
2442  enddo
2443  enddo
2444 !
2445  do n=1,ntracers
2446  do k=1,klev_in
2447  do j=js_input,je_input
2448  do i=is_input,ie_input
2449  bc_t1%north%q_BC(i,j,k,n)=tracers_input(i,j,k,n)
2450  enddo
2451  enddo
2452  enddo
2453  enddo
2454 !
2455  is_input=lbound(bc_t1%north%u_BC,1)
2456  ie_input=ubound(bc_t1%north%u_BC,1) ! Index limits for
2457  js_input=lbound(bc_t1%north%u_BC,2) ! D-grid u and C-grid v.
2458  je_input=ubound(bc_t1%north%u_BC,2)
2459 !
2460  do k=1,klev_in
2461  do j=js_input,je_input
2462  do i=is_input,ie_input
2463  bc_t1%north%u_BC(i,j,k)=u_s_input(i,j,k)
2464  bc_t1%north%vc_BC(i,j,k)=v_s_input(i,j,k)
2465  enddo
2466  enddo
2467  enddo
2468 !
2469  is_input=lbound(bc_t1%north%v_BC,1)
2470  ie_input=ubound(bc_t1%north%v_BC,1) ! Index limits for
2471  js_input=lbound(bc_t1%north%v_BC,2) ! D-grid v and C-grid u.
2472  je_input=ubound(bc_t1%north%v_BC,2)
2473 !
2474  do k=1,klev_in
2475  do j=js_input,je_input
2476  do i=is_input,ie_input
2477  bc_t1%north%v_BC(i,j,k)=v_w_input(i,j,k)
2478  bc_t1%north%uc_BC(i,j,k)=u_w_input(i,j,k)
2479  enddo
2480  enddo
2481  enddo
2482 !
2483  endif
2484 !
2485 !-----------
2486 !*** South
2487 !-----------
2488 !
2489  if(south_bc)then
2490  is_input=lbound(bc_t1%south%delp_BC,1)
2491  ie_input=ubound(bc_t1%south%delp_BC,1) ! Index limits for
2492  js_input=lbound(bc_t1%south%delp_BC,2) ! mass variables.
2493  je_input=ubound(bc_t1%south%delp_BC,2)
2494 !
2495  do k=1,klev_in
2496  do j=js_input,je_input
2497  do i=is_input,ie_input
2498  bc_t1%south%delp_BC(i,j,k)=delp_input(i,j,k)
2499  bc_t1%south%pt_BC(i,j,k)=t_input(i,j,k)
2500  bc_t1%south%w_BC(i,j,k)=w_input(i,j,k)
2501  bc_t1%south%delz_BC(i,j,k)=delz_input(i,j,k)
2502  enddo
2503  enddo
2504  enddo
2505 !
2506  do n=1,ntracers
2507  do k=1,klev_in
2508  do j=js_input,je_input
2509  do i=is_input,ie_input
2510  bc_t1%south%q_BC(i,j,k,n)=tracers_input(i,j,k,n)
2511  enddo
2512  enddo
2513  enddo
2514  enddo
2515 !
2516  is_input=lbound(bc_t1%south%u_BC,1)
2517  ie_input=ubound(bc_t1%south%u_BC,1) ! Index limits for
2518  js_input=lbound(bc_t1%south%u_BC,2) ! D-grid u and C-grid v.
2519  je_input=ubound(bc_t1%south%u_BC,2)
2520 !
2521  do k=1,klev_in
2522  do j=js_input,je_input
2523  do i=is_input,ie_input
2524  bc_t1%south%u_BC(i,j,k)=u_s_input(i,j,k)
2525  bc_t1%south%vc_BC(i,j,k)=v_s_input(i,j,k)
2526  enddo
2527  enddo
2528  enddo
2529 !
2530  is_input=lbound(bc_t1%south%v_BC,1)
2531  ie_input=ubound(bc_t1%south%v_BC,1) ! Index limits for
2532  js_input=lbound(bc_t1%south%v_BC,2) ! D-grid v and C-grid u.
2533  je_input=ubound(bc_t1%south%v_BC,2)
2534 !
2535  do k=1,klev_in
2536  do j=js_input,je_input
2537  do i=is_input,ie_input
2538  bc_t1%south%v_BC(i,j,k)=v_w_input(i,j,k)
2539  bc_t1%south%uc_BC(i,j,k)=u_w_input(i,j,k)
2540  enddo
2541  enddo
2542  enddo
2543 !
2544  endif
2545 !
2546 !----------
2547 !*** East
2548 !----------
2549 !
2550  if(east_bc)then
2551  is_input=lbound(bc_t1%east%delp_BC,1)
2552  ie_input=ubound(bc_t1%east%delp_BC,1) ! Index limits
2553  js_input=lbound(bc_t1%east%delp_BC,2) ! for mass variables.
2554  je_input=ubound(bc_t1%east%delp_BC,2)
2555 !
2556  do k=1,klev_in
2557  do j=js_input,je_input
2558  do i=is_input,ie_input
2559  bc_t1%east%delp_BC(i,j,k)=delp_input(i,j,k)
2560  bc_t1%east%pt_BC(i,j,k)=t_input(i,j,k)
2561  bc_t1%east%w_BC(i,j,k)=w_input(i,j,k)
2562  bc_t1%east%delz_BC(i,j,k)=delz_input(i,j,k)
2563  enddo
2564  enddo
2565  enddo
2566 !
2567  do n=1,ntracers
2568  do k=1,klev_in
2569  do j=js_input,je_input
2570  do i=is_input,ie_input
2571  bc_t1%east%q_BC(i,j,k,n)=tracers_input(i,j,k,n)
2572  enddo
2573  enddo
2574  enddo
2575  enddo
2576 !
2577  is_input=lbound(bc_t1%east%u_BC,1)
2578  ie_input=ubound(bc_t1%east%u_BC,1) ! Index limits for
2579  js_input=lbound(bc_t1%east%u_BC,2) ! D-grid u and C-grid v.
2580  je_input=ubound(bc_t1%east%u_BC,2)
2581 !
2582  do k=1,klev_in
2583  do j=js_input,je_input
2584  do i=is_input,ie_input
2585  bc_t1%east%u_BC(i,j,k)=u_s_input(i,j,k)
2586  bc_t1%east%vc_BC(i,j,k)=v_s_input(i,j,k)
2587  enddo
2588  enddo
2589  enddo
2590 !
2591  is_input=lbound(bc_t1%east%v_BC,1)
2592  ie_input=ubound(bc_t1%east%v_BC,1) ! Index limits for
2593  js_input=lbound(bc_t1%east%v_BC,2) ! D-grid v and C-grid u.
2594  je_input=ubound(bc_t1%east%v_BC,2)
2595 !
2596  do k=1,klev_in
2597  do j=js_input,je_input
2598  do i=is_input,ie_input
2599  bc_t1%east%v_BC(i,j,k)=v_w_input(i,j,k)
2600  bc_t1%east%uc_BC(i,j,k)=u_w_input(i,j,k)
2601  enddo
2602  enddo
2603  enddo
2604 !
2605  endif
2606 !
2607 !----------
2608 !*** West
2609 !----------
2610 !
2611  if(west_bc)then
2612  is_input=lbound(bc_t1%west%delp_BC,1)
2613  ie_input=ubound(bc_t1%west%delp_BC,1) ! Index limits for
2614  js_input=lbound(bc_t1%west%delp_BC,2) ! mass variables.
2615  je_input=ubound(bc_t1%west%delp_BC,2)
2616 !
2617  do k=1,klev_in
2618  do j=js_input,je_input
2619  do i=is_input,ie_input
2620  bc_t1%west%delp_BC(i,j,k)=delp_input(i,j,k)
2621  bc_t1%west%pt_BC(i,j,k)=t_input(i,j,k)
2622  bc_t1%west%w_BC(i,j,k)=w_input(i,j,k)
2623  bc_t1%west%delz_BC(i,j,k)=delz_input(i,j,k)
2624  enddo
2625  enddo
2626  enddo
2627 !
2628  do n=1,ntracers
2629  do k=1,klev_in
2630  do j=js_input,je_input
2631  do i=is_input,ie_input
2632  bc_t1%west%q_BC(i,j,k,n)=tracers_input(i,j,k,n)
2633  enddo
2634  enddo
2635  enddo
2636  enddo
2637 !
2638  is_input=lbound(bc_t1%west%u_BC,1)
2639  ie_input=ubound(bc_t1%west%u_BC,1) ! Index limits for
2640  js_input=lbound(bc_t1%west%u_BC,2) ! D-grid u and C-grid v.
2641  je_input=ubound(bc_t1%west%u_BC,2)
2642 !
2643  do k=1,klev_in
2644  do j=js_input,je_input
2645  do i=is_input,ie_input
2646  bc_t1%west%u_BC(i,j,k)=u_s_input(i,j,k)
2647  bc_t1%west%vc_BC(i,j,k)=v_s_input(i,j,k)
2648  enddo
2649  enddo
2650  enddo
2651 !
2652  is_input=lbound(bc_t1%west%v_BC,1)
2653  ie_input=ubound(bc_t1%west%v_BC,1) ! Index limits for
2654  js_input=lbound(bc_t1%west%v_BC,2) ! D-grid v and C-grid u.
2655  je_input=ubound(bc_t1%west%v_BC,2)
2656 !
2657  do k=1,klev_in
2658  do j=js_input,je_input
2659  do i=is_input,ie_input
2660  bc_t1%west%v_BC(i,j,k)=v_w_input(i,j,k)
2661  bc_t1%west%uc_BC(i,j,k)=u_w_input(i,j,k)
2662  enddo
2663  enddo
2664  enddo
2665 !
2666  endif
2667 !
2668 !-----------------------------------------------------------------------
2669 !
2670  end subroutine fill_bc_for_da
2671 !
2672 !-----------------------------------------------------------------------
2673 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2674 !-----------------------------------------------------------------------
2675 !
2676  subroutine fill_divgd_bc
2677 !
2678 !-----------------------------------------------------------------------
2679 !*** For now fill the boundary divergence with zero.
2680 !-----------------------------------------------------------------------
2681  implicit none
2682 !-----------------------------------------------------------------------
2683 !
2684 !--------------------
2685 !*** Local variables
2686 !--------------------
2687 !
2688  integer :: i,ie0,is0,j,je0,js0,k,nside
2689 !
2690  logical :: call_set
2691 !
2692 !-----------------------------------------------------------------------
2693 !***********************************************************************
2694 !-----------------------------------------------------------------------
2695 !
2696 !-----------------------------------------------------------------------
2697 !*** Loop through the four sides.
2698 !-----------------------------------------------------------------------
2699 !
2700  do nside=1,4
2701 !
2702  call_set=.false.
2703 !
2704  if(nside==1)then
2705  if(north_bc)then
2706  call_set=.true.
2707  bc_side_t1=>bc_t1%north
2708  is0=lbound(bc_t1%north%divgd_BC,1)
2709  ie0=ubound(bc_t1%north%divgd_BC,1)
2710  js0=lbound(bc_t1%north%divgd_BC,2)
2711  je0=ubound(bc_t1%north%divgd_BC,2)
2712  endif
2713  endif
2714 !
2715  if(nside==2)then
2716  if(south_bc)then
2717  call_set=.true.
2718  bc_side_t1=>bc_t1%south
2719  is0=lbound(bc_t1%south%divgd_BC,1)
2720  ie0=ubound(bc_t1%south%divgd_BC,1)
2721  js0=lbound(bc_t1%south%divgd_BC,2)
2722  je0=ubound(bc_t1%south%divgd_BC,2)
2723  endif
2724  endif
2725 !
2726  if(nside==3)then
2727  if(east_bc)then
2728  call_set=.true.
2729  bc_side_t1=>bc_t1%east
2730  is0=lbound(bc_t1%east%divgd_BC,1)
2731  ie0=ubound(bc_t1%east%divgd_BC,1)
2732  js0=lbound(bc_t1%east%divgd_BC,2)
2733  je0=ubound(bc_t1%east%divgd_BC,2)
2734  endif
2735  endif
2736 !
2737  if(nside==4)then
2738  if(west_bc)then
2739  call_set=.true.
2740  bc_side_t1=>bc_t1%west
2741  is0=lbound(bc_t1%west%divgd_BC,1)
2742  ie0=ubound(bc_t1%west%divgd_BC,1)
2743  js0=lbound(bc_t1%west%divgd_BC,2)
2744  je0=ubound(bc_t1%west%divgd_BC,2)
2745  endif
2746  endif
2747 !
2748  if(call_set)then
2749  do k=1,klev_out
2750  do j=js0,je0
2751  do i=is0,ie0
2752  bc_side_t1%divgd_BC(i,j,k)=0.
2753  enddo
2754  enddo
2755  enddo
2756  endif
2757 !
2758  enddo
2759 !
2760 !-----------------------------------------------------------------------
2761 !
2762  end subroutine fill_divgd_bc
2763 !
2764 !-----------------------------------------------------------------------
2765 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2766 !-----------------------------------------------------------------------
2767 !
2768 #ifdef USE_COND
2769  subroutine fill_q_con_bc
2770 !
2771 !-----------------------------------------------------------------------
2772 !*** For now fill the total condensate in the boundary regiona
2773 !*** with only the liquid water content.
2774 !-----------------------------------------------------------------------
2775  implicit none
2776 !-----------------------------------------------------------------------
2777 !
2778 !--------------------
2779 !*** Local variables
2780 !--------------------
2781 !
2782  integer :: i,ie0,is0,j,je0,js0,k,nside
2783 !
2784  logical :: call_set
2785 !
2786 !-----------------------------------------------------------------------
2787 !***********************************************************************
2788 !-----------------------------------------------------------------------
2789 !
2790 !-----------------------------------------------------------------------
2791 !*** Loop through the four sides.
2792 !-----------------------------------------------------------------------
2793 !
2794  do nside=1,4
2795  call_set=.false.
2796 !
2797  if(nside==1)then
2798  if(north_bc)then
2799  call_set=.true.
2800  bc_side_t1=>bc_t1%north
2801  is0=lbound(bc_t1%north%q_con_BC,1)
2802  ie0=ubound(bc_t1%north%q_con_BC,1)
2803  js0=lbound(bc_t1%north%q_con_BC,2)
2804  je0=ubound(bc_t1%north%q_con_BC,2)
2805  endif
2806  endif
2807 !
2808  if(nside==2)then
2809  if(south_bc)then
2810  call_set=.true.
2811  bc_side_t1=>bc_t1%south
2812  is0=lbound(bc_t1%south%q_con_BC,1)
2813  ie0=ubound(bc_t1%south%q_con_BC,1)
2814  js0=lbound(bc_t1%south%q_con_BC,2)
2815  je0=ubound(bc_t1%south%q_con_BC,2)
2816  endif
2817  endif
2818 !
2819  if(nside==3)then
2820  if(east_bc)then
2821  call_set=.true.
2822  bc_side_t1=>bc_t1%east
2823  is0=lbound(bc_t1%east%q_con_BC,1)
2824  ie0=ubound(bc_t1%east%q_con_BC,1)
2825  js0=lbound(bc_t1%east%q_con_BC,2)
2826  je0=ubound(bc_t1%east%q_con_BC,2)
2827  endif
2828  endif
2829 !
2830  if(nside==4)then
2831  if(west_bc)then
2832  call_set=.true.
2833  bc_side_t1=>bc_t1%west
2834  is0=lbound(bc_t1%west%q_con_BC,1)
2835  ie0=ubound(bc_t1%west%q_con_BC,1)
2836  js0=lbound(bc_t1%west%q_con_BC,2)
2837  je0=ubound(bc_t1%west%q_con_BC,2)
2838  endif
2839  endif
2840 !
2841  if(call_set)then
2842  do k=1,klev_out
2843  do j=js0,je0
2844  do i=is0,ie0
2845  bc_side_t1%q_con_BC(i,j,k)=bc_side_t1%q_BC(i,j,k,liq_water_index)
2846  enddo
2847  enddo
2848  enddo
2849  endif
2850 !
2851  enddo
2852 !
2853 !-----------------------------------------------------------------------
2854 !
2855  end subroutine fill_q_con_bc
2856 #endif
2857 !
2858 !-----------------------------------------------------------------------
2859 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2860 !-----------------------------------------------------------------------
2861 !
2862 #ifdef MOIST_CAPPA
2863  subroutine fill_cappa_bc
2864 !
2865 !-----------------------------------------------------------------------
2866 !*** Compute cappa in the regional domain boundary area following
2867 !*** Zhao-Carr microphysics.
2868 !-----------------------------------------------------------------------
2869  implicit none
2870 !-----------------------------------------------------------------------
2871 !
2872 !---------------------
2873 !*** Local variables
2874 !---------------------
2875 !
2876  integer :: i1,i2,j1,j2,nside
2877 !
2878  real,dimension(:,:,:),pointer :: cappa,temp,liq_wat,sphum
2879 !
2880  logical :: call_compute
2881 !
2882 !-----------------------------------------------------------------------
2883 !***********************************************************************
2884 !-----------------------------------------------------------------------
2885 !
2886  do nside=1,4
2887  call_compute=.false.
2888 !
2889  if(nside==1)then
2890  if(north_bc)then
2891  call_compute=.true.
2892  bc_side_t1=>bc_t1%north
2893  endif
2894  endif
2895 !
2896  if(nside==2)then
2897  if(south_bc)then
2898  call_compute=.true.
2899  bc_side_t1=>bc_t1%south
2900  endif
2901  endif
2902 !
2903  if(nside==3)then
2904  if(east_bc)then
2905  call_compute=.true.
2906  bc_side_t1=>bc_t1%east
2907  endif
2908  endif
2909 !
2910  if(nside==4)then
2911  if(west_bc)then
2912  call_compute=.true.
2913  bc_side_t1=>bc_t1%west
2914  endif
2915  endif
2916 !
2917  if(call_compute)then
2918  i1=lbound(bc_side_t1%cappa_BC,1)
2919  i2=ubound(bc_side_t1%cappa_BC,1)
2920  j1=lbound(bc_side_t1%cappa_BC,2)
2921  j2=ubound(bc_side_t1%cappa_BC,2)
2922  cappa =>bc_side_t1%cappa_BC
2923  temp =>bc_side_t1%pt_BC
2924  liq_wat=>bc_side_t1%q_BC(:,:,:,liq_water_index)
2925  sphum =>bc_side_t1%q_BC(:,:,:,sphum_index)
2926  call compute_cappa(i1,i2,j1,j2,cappa,temp,liq_wat,sphum)
2927  endif
2928 !
2929  enddo
2930 !
2931 !-----------------------------------------------------------------------
2932 !
2933  end subroutine fill_cappa_bc
2934 !
2935 !-----------------------------------------------------------------------
2936 !***********************************************************************
2937 !-----------------------------------------------------------------------
2938 !
2939  subroutine compute_cappa(i1,i2,j1,j2,cappa,temp,liq_wat,sphum)
2940 !
2941 !-----------------------------------------------------------------------
2942  implicit none
2943 !-----------------------------------------------------------------------
2944 !
2945 !------------------------
2946 !*** Argument variables
2947 !------------------------
2948 !
2949  integer,intent(in) :: i1,i2,j1,j2
2950 !
2951  real,dimension(i1:i2,j1:j2,1:npz),intent(in) :: temp,liq_wat,sphum
2952  real,dimension(i1:i2,j1:j2,1:npz),intent(inout) :: cappa
2953 !
2954 !---------------------
2955 !*** Local variables
2956 !---------------------
2957 !
2958  integer :: i,ie,is,j,je,js,k
2959 !
2960  real :: cvm,qd,ql,qs,qv
2961 !
2962 !-----------------------------------------------------------------------
2963 !***********************************************************************
2964 !-----------------------------------------------------------------------
2965 !
2966  is=lbound(cappa,1)
2967  ie=ubound(cappa,1)
2968  js=lbound(cappa,2)
2969  je=ubound(cappa,2)
2970 !
2971  do k=1,klev_out
2972  do j=js,je
2973  do i=is,ie
2974  qd=max(0.,liq_wat(i,j,k))
2975  if( temp(i,j,k) > tice )then
2976  qs=0.
2977  elseif( temp(i,j,k) < tice-t_i0 )then
2978  qs=qd
2979  else
2980  qs=qd*(tice-temp(i,j,k))/t_i0
2981  endif
2982  ql=qd-qs
2983  qv=max(0.,sphum(i,j,k))
2984  cvm=(1.-(qv+qd))*cv_air + qv*cv_vap + ql*c_liq + qs*c_ice
2985 !
2986  cappa(i,j,k)=rdgas/(rdgas+cvm/(1.+zvir*sphum(i,j,k)))
2987 !
2988  enddo
2989  enddo
2990  enddo
2991 !
2992 !-----------------------------------------------------------------------
2993 !
2994  end subroutine compute_cappa
2995 #endif
2996 !
2997 !-----------------------------------------------------------------------
2998 !
2999  end subroutine regional_bc_data
3000 
3001 !-----------------------------------------------------------------------
3002 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3003 !-----------------------------------------------------------------------
3004  subroutine read_regional_bc_file(is_input,ie_input &
3005  ,js_input,je_input &
3006  ,nlev &
3007  ,ntracers &
3008  ,var_name_root &
3009  ,array_3d &
3010  ,array_4d &
3011  ,tlev &
3012  ,required )
3013 !-----------------------------------------------------------------------
3014 !*** Read the boundary data from the external file generated by
3015 !*** chgres.
3016 !-----------------------------------------------------------------------
3017 !-----------------------------------------------------------------------
3018  implicit none
3019 !-----------------------------------------------------------------------
3020 !
3021 !------------------------
3022 !*** Argument variables
3023 !------------------------
3024 !
3025 !----------
3026 !*** Input
3027 !----------
3028 !
3029  integer,intent(in) :: is_input,ie_input,js_input,je_input,nlev
3030  integer,intent(in) :: ntracers
3031 !
3032  integer,intent(in),optional :: tlev
3033 !
3034  character(len=*),intent(in) :: var_name_root
3035  logical,intent(in),optional :: required
3036 !
3037 !------------
3038 !*** Output
3039 !------------
3040 !
3041  real,dimension(is_input:ie_input,js_input:je_input,1:nlev),intent(out),optional :: array_3d
3042 !
3043  real,dimension(is_input:ie_input,js_input:je_input,1:nlev,1:ntracers),intent(out),optional :: array_4d
3044 !
3045 !---------------------
3046 !*** Local variables
3047 !---------------------
3048 !
3049  integer :: halo,lat,lev,lon
3050 !
3051  integer :: i_count,i_start_array,i_start_data,i_end_array &
3052  ,j_count,j_start_array,j_start_data,j_end_array
3053 !
3054  integer :: dim_id,nctype,ndims,var_id
3055  integer :: nside,status
3056 !
3057  character(len=5) :: dim_name_x & !<-- Dimension names in
3058  ,dim_name_y ! the BC file
3059 !
3060  character(len=80) :: var_name
3061 !
3062  logical :: call_get_var,is_root_pe
3063  logical :: required_local
3064 !
3065 !-----------------------------------------------------------------------
3066 !***********************************************************************
3067 !-----------------------------------------------------------------------
3068 !
3069 !-----------------------------------------------------------------------
3070 !*** Process optional argument required, default value is .true.
3071 !-----------------------------------------------------------------------
3072 !
3073  if(present(required)) then
3074  required_local=required
3075  else
3076  required_local=.true.
3077  endif
3078 !
3079  is_root_pe=(mpp_pe()==mpp_root_pe())
3080 !
3081 !-----------------------------------------------------------------------
3082 !*** Loop through the four sides of the domain.
3083 !-----------------------------------------------------------------------
3084 !
3085 !-----------------------------------------------------------------------
3086  sides: do nside=1,4
3087 !-----------------------------------------------------------------------
3088 !
3089  call_get_var=.false.
3090 !
3091 !-----------------------------------------------------------------------
3092 !*** Construct the variable's name in the NetCDF file and set
3093 !*** the start locations and point counts for the data file and
3094 !*** for the BC arrays being filled. The input array begins
3095 !*** receiving data at (i_start_array,j_start_array), etc.
3096 !*** The read of the data for the given input array begins at
3097 !*** (i_start_data,j_start_data) and encompasses i_count by
3098 !*** j_count datapoints in each direction.
3099 !-----------------------------------------------------------------------
3100 !
3101 !-----------
3102 !*** North
3103 !-----------
3104 !
3105  if(nside==1)then
3106  if(north_bc)then
3107  call_get_var=.true.
3108 !
3109  var_name=trim(var_name_root)//"_bottom"
3110 !
3111  i_start_array=is_input
3112  i_end_array =ie_input
3113  j_start_array=js_input
3114  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
3115  j_end_array=js_input+nhalo_data+nrows_blend
3116  else
3117  j_end_array=js_input+nhalo_data+nrows_blend-1
3118  endif
3119 !
3120  i_start_data=i_start_array+nhalo_data
3121  i_count=i_end_array-i_start_array+1
3122  j_start_data=1
3123  j_count=j_end_array-j_start_array+1
3124 !
3125  endif
3126  endif
3127 !
3128 !-----------
3129 !*** South
3130 !-----------
3131 !
3132  if(nside==2)then
3133  if(south_bc)then
3134  call_get_var=.true.
3135 !
3136  var_name=trim(var_name_root)//"_top"
3137 !
3138  i_start_array=is_input
3139  i_end_array =ie_input
3140  j_start_array=je_input-nhalo_data-nrows_blend+1
3141  j_end_array =je_input
3142 !
3143  i_start_data=i_start_array+nhalo_data
3144  i_count=i_end_array-i_start_array+1
3145  j_start_data=1
3146  j_count=j_end_array-j_start_array+1
3147 !
3148  endif
3149  endif
3150 !
3151 !----------
3152 !*** East
3153 !----------
3154 !
3155  if(nside==3)then
3156  if(east_bc)then
3157  call_get_var=.true.
3158 !
3159  var_name=trim(var_name_root)//"_left"
3160 !
3161  j_start_array=js_input
3162  j_end_array =je_input
3163 !
3164  i_start_array=is_input
3165 !
3166  if(trim(var_name_root)=='u_w'.or.trim(var_name_root)=='v_w')then
3167  i_end_array=is_input+nhalo_data+nrows_blend
3168  else
3169  i_end_array=is_input+nhalo_data+nrows_blend-1
3170  endif
3171 !
3172  if(north_bc)then
3173  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
3174  j_start_array=js_input+nhalo_data+1
3175  else
3176  j_start_array=js_input+nhalo_data
3177  endif
3178  endif
3179  if(south_bc)then
3180  j_end_array =je_input-nhalo_data
3181  endif
3182 !
3183  i_start_data=1
3184  i_count=i_end_array-i_start_array+1
3185  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
3186  j_start_data=j_start_array-1
3187  else
3188  j_start_data=j_start_array
3189  endif
3190  j_count=j_end_array-j_start_array+1
3191 !
3192  endif
3193  endif
3194 !
3195 !----------
3196 !*** West
3197 !----------
3198 !
3199  if(nside==4)then
3200  if(west_bc)then
3201  call_get_var=.true.
3202 !
3203  var_name=trim(var_name_root)//"_right"
3204 !
3205  j_start_array=js_input
3206  j_end_array =je_input
3207 !
3208  i_start_array=ie_input-nhalo_data-nrows_blend+1
3209  i_end_array=ie_input
3210 !
3211  if(north_bc)then
3212  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
3213  j_start_array=js_input+nhalo_data+1
3214  else
3215  j_start_array=js_input+nhalo_data
3216  endif
3217  endif
3218 !
3219  if(south_bc)then
3220  j_end_array =je_input-nhalo_data
3221  endif
3222 !
3223  i_start_data=1
3224  i_count=i_end_array-i_start_array+1
3225  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
3226  j_start_data=j_start_array-1
3227  else
3228  j_start_data=j_start_array
3229  endif
3230  j_count=j_end_array-j_start_array+1
3231 !
3232  endif
3233  endif
3234 !
3235 !-----------------------------------------------------------------------
3236 !*** Fill this task's subset of boundary data for this 3-D
3237 !*** or 4-D variable. This includes the data in the domain's
3238 !*** halo region as well as the blending region that overlaps
3239 !*** the outer nhalo_blend rows of the integration domain.
3240 !*** If the variable is a tracer then check if it is present
3241 !*** in the input data. If it is not then print a warning
3242 !*** and set it to zero.
3243 !-----------------------------------------------------------------------
3244 !
3245  if(call_get_var)then
3246  if (present(array_4d)) then
3247  status=nf90_inq_varid(ncid,trim(var_name),var_id)
3248  if (required_local) then
3249  call check(status)
3250  endif
3251  if (status /= nf90_noerr) then
3252  if (east_bc.and.is_master()) write(0,*)' WARNING: Tracer ',trim(var_name),' not in input file'
3253  array_4d(:,:,:,tlev)=0.
3254 !
3255  blend_this_tracer(tlev)=.false.
3256 !
3257  else
3258  call check(nf90_get_var(ncid,var_id &
3259  ,array_4d(i_start_array:i_end_array &
3260  ,j_start_array:j_end_array &
3261  ,1:nlev, tlev) &
3262  ,start=(/i_start_data,j_start_data,1,tlev/) &
3263  ,count=(/i_count,j_count,nlev,1/)))
3264 !
3265  endif
3266 !
3267  else
3268  call check(nf90_inq_varid(ncid,trim(var_name),var_id))
3269  call check(nf90_get_var(ncid,var_id &
3270  ,array_3d(i_start_array:i_end_array &
3271  ,j_start_array:j_end_array &
3272  ,1:nlev) &
3273  ,start=(/i_start_data,j_start_data,1/) &
3274  ,count=(/i_count,j_count,nlev/)))
3275 !
3276  endif
3277  endif
3278 !
3279  enddo sides
3280 !
3281 !-----------------------------------------------------------------------
3282 !
3283  end subroutine read_regional_bc_file
3284 !
3285 !-----------------------------------------------------------------------
3286 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3287 !-----------------------------------------------------------------------
3288 !
3289  subroutine check(status)
3290  integer,intent(in) :: status
3291 !
3292  if(status /= nf90_noerr) then
3293  write(0,*)' check netcdf status=',status
3294  call mpp_error(fatal, ' NetCDF error ' // trim(nf90_strerror(status)))
3295  endif
3296 !
3297  end subroutine check
3298 !
3299 !-----------------------------------------------------------------------
3300 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3301 !-----------------------------------------------------------------------
3302 !
3303  subroutine allocate_regional_bc_arrays(side &
3304  ,north_bc,south_bc &
3305  ,east_bc,west_bc &
3306  ,is_0,ie_0,js_0,je_0 &
3307  ,is_sn,ie_sn,js_sn,je_sn &
3308  ,is_we,ie_we,js_we,je_we &
3309  ,klev &
3310  ,ntracers &
3311  ,BC_side &
3312  ,delz_side )
3313 !
3314 !-----------------------------------------------------------------------
3315  implicit none
3316 !-----------------------------------------------------------------------
3317 !
3318 !------------------------
3319 !*** Argument variables
3320 !------------------------
3321 !
3322  integer,intent(in) :: klev,ntracers
3323 !
3324  integer,intent(in) :: is_0,ie_0,js_0,je_0
3325  integer,intent(in) :: is_sn,ie_sn,js_sn,je_sn
3326  integer,intent(in) :: is_we,ie_we,js_we,je_we
3327 !
3328  character(len=5),intent(in) :: side
3329 !
3330  logical,intent(in) :: north_bc,south_bc,east_bc,west_bc
3331 !
3332  type(fv_regional_bc_variables),intent(out) :: BC_side
3333 !
3334  real,dimension(:,:,:),pointer,intent(inout),optional :: delz_side
3335 !
3336 !---------------------------------------------------------------------
3337 !*********************************************************************
3338 !---------------------------------------------------------------------
3339 !
3340  if(allocated(bc_side%delp_BC))then
3341  return
3342  endif
3343 !
3344  allocate(bc_side%delp_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%delp_BC=real_snan
3345  allocate(bc_side%divgd_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%divgd_BC=real_snan
3346 !
3347  allocate(bc_side%q_BC (is_0:ie_0,js_0:je_0,1:klev,1:ntracers)) ; bc_side%q_BC=real_snan
3348 !
3349  if(.not.allocated(blend_this_tracer))then
3350  allocate(blend_this_tracer(1:ntracers))
3351  blend_this_tracer=.true.
3352  endif
3353 !
3354 #ifndef SW_DYNAMICS
3355  allocate(bc_side%pt_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%pt_BC=real_snan
3356  allocate(bc_side%w_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%w_BC=real_snan
3357  allocate(bc_side%delz_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%delz_BC=real_snan
3358  if(present(delz_side))then
3359  if(.not.associated(delz_side))then
3360  allocate(delz_side(is_0:ie_0,js_0:je_0,klev)) ; delz_side=real_snan
3361  endif
3362  endif
3363 #ifdef USE_COND
3364  allocate(bc_side%q_con_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%q_con_BC=real_snan
3365 #ifdef MOIST_CAPPA
3366  allocate(bc_side%cappa_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%cappa_BC=real_snan
3367 #endif
3368 #endif
3369 #endif
3370 !
3371 !--------------------
3372 !*** Wind components
3373 !--------------------
3374 !
3375 !** D-grid u, C-grid v
3376 !
3377  allocate(bc_side%u_BC (is_sn:ie_sn, js_sn:je_sn, klev)) ; bc_side%u_BC=real_snan
3378  allocate(bc_side%vc_BC(is_sn:ie_sn, js_sn:je_sn, klev)) ; bc_side%vc_BC=real_snan
3379 !
3380 !** C-grid u, D-grid v
3381 !
3382  allocate(bc_side%uc_BC(is_we:ie_we, js_we:je_we, klev)) ; bc_side%uc_BC=real_snan
3383  allocate(bc_side%v_BC (is_we:ie_we, js_we:je_we, klev)) ; bc_side%v_BC=real_snan
3384 !
3385 !---------------------------------------------------------------------
3386 !
3387  end subroutine allocate_regional_bc_arrays
3388 !
3389 !---------------------------------------------------------------------
3390 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3391 !---------------------------------------------------------------------
3392 
3393 subroutine remap_scalar_nggps_regional_bc(Atm &
3394  ,side &
3395  ,isd,ied,jsd,jed &
3396  ,is_bc,ie_bc,js_bc,je_bc &
3397  ,km, npz, ncnst, ak0, bk0 &
3398  ,psc, t_in, qa, omga, zh &
3399  ,phis_reg &
3400  ,ps &
3401  ,BC_side )
3402 
3403  type(fv_atmos_type), intent(inout) :: Atm
3404  integer, intent(in):: isd,ied,jsd,jed
3405  integer, intent(in):: is_bc,ie_bc,js_bc,je_bc
3406  integer, intent(in):: km & !<-- # of levels in 3-D input variables
3407  ,npz & !<-- # of levels in final 3-D integration variables
3408  ,ncnst
3409  real, intent(in):: ak0(km+1), bk0(km+1)
3410  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc):: psc
3411  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: t_in
3412  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: omga
3413  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km,ncnst):: qa
3414  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km+1):: zh
3415  real, intent(inout), dimension(isd-1:ied+1,jsd-1:jed+1):: phis_reg
3416  real, intent(out),dimension(is_bc:ie_bc,js_bc:je_bc) :: ps
3417  character(len=5),intent(in) :: side
3418  type(fv_regional_bc_variables),intent(inout) :: BC_side
3419 
3420 ! local:
3421 !
3422  real, dimension(:,:),allocatable :: pe0
3423  real, dimension(:,:),allocatable :: qn1
3424  real, dimension(:,:),allocatable :: dp2
3425  real, dimension(:,:),allocatable :: pe1
3426  real, dimension(:,:),allocatable :: qp
3427 !
3428  real wk(is_bc:ie_bc,js_bc:je_bc)
3429  real, dimension(is_bc:ie_bc,js_bc:je_bc):: phis
3430 
3431 !!! High-precision
3432  real(kind=R_GRID), dimension(is_bc:ie_bc,npz+1):: pn1
3433  real(kind=R_GRID):: gz_fv(npz+1)
3434  real(kind=R_GRID), dimension(2*km+1):: gz, pn
3435  real(kind=R_GRID), dimension(is_bc:ie_bc,km+1):: pn0
3436  real(kind=R_GRID):: pst
3437 !!! High-precision
3438  integer i,ie,is,j,je,js,k,l,m, k2,iq
3439  integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
3440 !
3441 !---------------------------------------------------------------------------------
3442 !
3443  sphum = sphum_index
3444  liq_wat = liq_water_index
3445  ice_wat = ice_water_index
3446  rainwat = rain_water_index
3447  snowwat = snow_water_index
3448  graupel = graupel_index
3449  cld_amt = cld_amt_index
3450  o3mr = o3mr_index
3451 
3452  k2 = max(10, km/2)
3453 
3454  if (mpp_pe()==1) then
3455  print *, 'sphum = ', sphum
3456  print *, 'clwmr = ', liq_wat
3457  print *, ' o3mr = ', o3mr
3458  print *, 'ncnst = ', ncnst
3459  print *, 'ntracers = ', ntracers
3460  endif
3461 
3462  if ( sphum/=1 ) then
3463  call mpp_error(fatal,'SPHUM must be 1st tracer')
3464  endif
3465 !
3466 !---------------------------------------------------------------------------------
3467 !*** First compute over the extended boundary regions with halo=nhalo_data.
3468 !*** This is needed to obtain pressures that will surround the wind points.
3469 !---------------------------------------------------------------------------------
3470 !
3471  is=is_bc
3472  if(side=='west')then
3473  is=ie_bc-nhalo_data-nrows_blend+1
3474  endif
3475 !
3476  ie=ie_bc
3477  if(side=='east')then
3478  ie=is_bc+nhalo_data+nrows_blend-1
3479  endif
3480 !
3481  js=js_bc
3482  if(side=='south')then
3483  js=je_bc-nhalo_data-nrows_blend+1
3484  endif
3485 !
3486  je=je_bc
3487  if(side=='north')then
3488  je=js_bc+nhalo_data+nrows_blend-1
3489  endif
3490 !
3491  allocate(pe0(is:ie,km+1)) ; pe0=real_snan
3492  allocate(qn1(is:ie,npz)) ; qn1=real_snan
3493  allocate(dp2(is:ie,npz)) ; dp2=real_snan
3494  allocate(pe1(is:ie,npz+1)) ; pe1=real_snan
3495  allocate(qp(is:ie,km)) ; qp=real_snan
3496 !
3497 !---------------------------------------------------------------------------------
3498  jloop1: do j=js,je
3499 !---------------------------------------------------------------------------------
3500 !
3501  do k=1,km+1
3502  do i=is,ie
3503  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
3504  pn0(i,k) = log(pe0(i,k))
3505  enddo
3506  enddo
3507 
3508  do i=is,ie
3509  do k=1,km+1
3510  pn(k) = pn0(i,k)
3511  gz(k) = zh(i,j,k)*grav
3512  enddo
3513 ! Use log-p for interpolation/extrapolation
3514 ! mirror image method:
3515  do k=km+2, km+k2
3516  l = 2*(km+1) - k
3517  gz(k) = 2.*gz(km+1) - gz(l)
3518  pn(k) = 2.*pn(km+1) - pn(l)
3519  enddo
3520 
3521  do k=km+k2-1, 2, -1
3522  if( phis_reg(i,j).le.gz(k) .and. phis_reg(i,j).ge.gz(k+1) ) then
3523  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-phis_reg(i,j))/(gz(k)-gz(k+1))
3524  go to 123
3525  endif
3526  enddo
3527  123 ps(i,j) = exp(pst)
3528 
3529  enddo ! i-loop
3530 
3531 !---------------------------------------------------------------------------------
3532  enddo jloop1
3533 !---------------------------------------------------------------------------------
3534 
3535 !---------------------------------------------------------------------------------
3536 !*** Transfer values from the expanded boundary array for sfc pressure into
3537 !*** the Atm object.
3538 !---------------------------------------------------------------------------------
3539 !
3540  is=lbound(atm%ps,1)
3541  ie=ubound(atm%ps,1)
3542  js=lbound(atm%ps,2)
3543  je=ubound(atm%ps,2)
3544 !
3545  do j=js,je
3546  do i=is,ie
3547  atm%ps(i,j)=ps(i,j)
3548  enddo
3549  enddo
3550 !
3551 !---------------------------------------------------------------------------------
3552 !*** Now compute over the normal boundary regions with halo=nhalo_model
3553 !*** extended through nrows_blend rows into the integration domain.
3554 !*** Use the dimensions of one of the permanent BC variables in Atm
3555 !*** as the loop limits so any side of the domain can be addressed.
3556 !---------------------------------------------------------------------------------
3557 !
3558  is=lbound(bc_side%delp_BC,1)
3559  ie=ubound(bc_side%delp_BC,1)
3560  js=lbound(bc_side%delp_BC,2)
3561  je=ubound(bc_side%delp_BC,2)
3562 !
3563 !---------------------------------------------------------------------------------
3564  jloop2: do j=js,je
3565 !---------------------------------------------------------------------------------
3566  do k=1,km+1
3567  do i=is,ie
3568  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
3569  pn0(i,k) = log(pe0(i,k))
3570  enddo
3571  enddo
3572 !
3573  do i=is,ie
3574  pe1(i,1) = atm%ak(1)
3575  pn1(i,1) = log(pe1(i,1))
3576  enddo
3577  do k=2,npz+1
3578  do i=is,ie
3579  pe1(i,k) = atm%ak(k) + atm%bk(k)*ps(i,j)
3580  pn1(i,k) = log(pe1(i,k))
3581  enddo
3582  enddo
3583 
3584 ! * Compute delp
3585  do k=1,npz
3586  do i=is,ie
3587  dp2(i,k) = pe1(i,k+1) - pe1(i,k)
3588  bc_side%delp_BC(i,j,k) = dp2(i,k)
3589  enddo
3590  enddo
3591 
3592 ! map shpum, o3mr, liq_wat tracers
3593  do iq=1,ncnst
3594 ! if (iq == sphum .or. iq == liq_wat .or. iq == o3mr) then ! only remap if the data is already set
3595  if (iq /= cld_amt) then ! don't remap cld_amt
3596  do k=1,km
3597  do i=is,ie
3598  qp(i,k) = qa(i,j,k,iq)
3599  enddo
3600  enddo
3601 
3602  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 0, 8, atm%ptop)
3603 
3604  if ( iq==sphum ) then
3605  call fillq(ie-is+1, npz, 1, qn1, dp2)
3606  else
3607  call fillz(ie-is+1, npz, 1, qn1, dp2)
3608  endif
3609 ! The HiRam step of blending model sphum with NCEP data is obsolete because nggps is always cold starting...
3610  do k=1,npz
3611  do i=is,ie
3612  bc_side%q_BC(i,j,k,iq) = qn1(i,k)
3613  enddo
3614  enddo
3615  endif
3616  enddo
3617 
3618 !---------------------------------------------------
3619 ! Retrieve temperature using GFS geopotential height
3620 !---------------------------------------------------
3621 !
3622  i_loop: do i=is,ie
3623 !
3624 ! Make sure FV3 top is lower than GFS; can not do extrapolation above the top at this point
3625  if ( pn1(i,1) .lt. pn0(i,1) ) then
3626  call mpp_error(fatal,'FV3 top higher than NCEP/GFS')
3627  endif
3628 
3629  do k=1,km+1
3630  pn(k) = pn0(i,k)
3631  gz(k) = zh(i,j,k)*grav
3632  enddo
3633 !-------------------------------------------------
3634  do k=km+2, km+k2
3635  l = 2*(km+1) - k
3636  gz(k) = 2.*gz(km+1) - gz(l)
3637  pn(k) = 2.*pn(km+1) - pn(l)
3638  enddo
3639 !-------------------------------------------------
3640 
3641  gz_fv(npz+1) = phis_reg(i,j)
3642 
3643  m = 1
3644 
3645  do k=1,npz
3646 ! Searching using FV3 log(pe): pn1
3647 #ifdef USE_ISOTHERMO
3648  do l=m,km
3649  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
3650  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
3651  goto 555
3652  elseif ( pn1(i,k) .gt. pn(km+1) ) then
3653 ! Isothermal under ground; linear in log-p extra-polation
3654  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))
3655  goto 555
3656  endif
3657  enddo
3658 #else
3659  do l=m,km+k2-1
3660  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
3661  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
3662  goto 555
3663  endif
3664  enddo
3665 #endif
3666 555 m = l
3667  enddo
3668 
3669 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
3670 !xxx DO WE NEED Atm%peln to have values in the boundary region?
3671 !xxx FOR NOW COMMENT IT OUT.
3672 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
3673 !xxx do k=1,npz+1
3674 !xxx Atm%peln(i,k,j) = pn1(i,k)
3675 !xxx enddo
3676 
3677 ! Compute true temperature using hydrostatic balance if not read from input.
3678 
3679  if (data_source /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
3680  do k=1,npz
3681  bc_side%pt_BC(i,j,k) = (gz_fv(k)-gz_fv(k+1))/( rdgas*(pn1(i,k+1)-pn1(i,k))*(1.+zvir*bc_side%q_BC(i,j,k,sphum)) )
3682  enddo
3683  endif
3684 
3685 
3686  if ( .not. atm%flagstruct%hydrostatic ) then
3687  do k=1,npz
3688  bc_side%delz_BC(i,j,k) = (gz_fv(k+1) - gz_fv(k)) / grav
3689  enddo
3690  endif
3691 
3692  enddo i_loop
3693 
3694 !-----------------------------------------------------------------------
3695 ! separate cloud water and cloud ice
3696 ! From Jan-Huey Chen's HiRAM code
3697 !-----------------------------------------------------------------------
3698 !
3699 ! If the source is FV3GFS GAUSSIAN NEMSIO FILE then all the tracers are in the boundary files
3700 ! and will be read in.
3701 ! If the source is from old GFS or operational GSM then the tracers will be fixed in the boundaries
3702 ! and may not provide a very good result
3703 !
3704  if (cld_amt .gt. 0) bc_side%q_BC(:,:,:,cld_amt) = 0.
3705  if (trim(data_source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
3706  if ( atm%flagstruct%nwat .eq. 6 ) then
3707  do k=1,npz
3708  do i=is,ie
3709  qn1(i,k) = bc_side%q_BC(i,j,k,liq_wat)
3710  bc_side%q_BC(i,j,k,rainwat) = 0.
3711  bc_side%q_BC(i,j,k,snowwat) = 0.
3712  bc_side%q_BC(i,j,k,graupel) = 0.
3713  if ( bc_side%pt_BC(i,j,k) > 273.16 ) then ! > 0C all liq_wat
3714  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)
3715  bc_side%q_BC(i,j,k,ice_wat) = 0.
3716 #ifdef ORIG_CLOUDS_PART
3717  else if ( bc_side%pt_BC(i,j,k) < 258.16 ) then ! < -15C all ice_wat
3718  bc_side%q_BC(i,j,k,liq_wat) = 0.
3719  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
3720  else ! between -15~0C: linear interpolation
3721  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-258.16)/15.)
3722  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3723  endif
3724 #else
3725  else if ( bc_side%pt_BC(i,j,k) < 233.16 ) then ! < -40C all ice_wat
3726  bc_side%q_BC(i,j,k,liq_wat) = 0.
3727  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
3728  else
3729  if ( k.eq.1 ) then ! between [-40,0]: linear interpolation
3730  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-233.16)/40.)
3731  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3732  else
3733  if (bc_side%pt_BC(i,j,k)<258.16 .and. bc_side%q_BC(i,j,k-1,ice_wat)>1.e-5 ) then
3734  bc_side%q_BC(i,j,k,liq_wat) = 0.
3735  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
3736  else ! between [-40,0]: linear interpolation
3737  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-233.16)/40.)
3738  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3739  endif
3740  endif
3741  endif
3742 #endif
3743  call mp_auto_conversion(bc_side%q_BC(i,j,k,liq_wat), bc_side%q_BC(i,j,k,rainwat), &
3744  bc_side%q_BC(i,j,k,ice_wat), bc_side%q_BC(i,j,k,snowwat) )
3745  enddo
3746  enddo
3747  endif
3748  endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
3749 !
3750 ! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
3751 ! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
3752 ! For FV3GFS NEMSIO input, w is already in m/s (but the code reads in as omga) and just needs to be remapped
3753 !-------------------------------------------------------------
3754 ! map omega
3755 !------- ------------------------------------------------------
3756  if ( .not. atm%flagstruct%hydrostatic ) then
3757  do k=1,km
3758  do i=is,ie
3759  qp(i,k) = omga(i,j,k)
3760  enddo
3761  enddo
3762 
3763  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, atm%ptop)
3764 
3765  if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
3766  do k=1,npz
3767  do i=is,ie
3768  bc_side%w_BC(i,j,k) = qn1(i,k)
3769  enddo
3770  enddo
3771 !------------------------------
3772 ! Remap input T linearly in p.
3773 !------------------------------
3774  do k=1,km
3775  do i=is,ie
3776  qp(i,k) = t_in(i,j,k)
3777  enddo
3778  enddo
3779 
3780  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 2, 4, atm%ptop)
3781 
3782  do k=1,npz
3783  do i=is,ie
3784  bc_side%pt_BC(i,j,k) = qn1(i,k)
3785  enddo
3786  enddo
3787 
3788  else
3789  do k=1,npz
3790  do i=is,ie
3791  bc_side%w_BC(i,j,k) = qn1(i,k)/bc_side%delp_BC(i,j,k)*bc_side%delz_BC(i,j,k)
3792  enddo
3793  enddo
3794  endif
3795 
3796  endif !.not. Atm%flagstruct%hydrostatic
3797 
3798  enddo jloop2
3799 
3800 ! Add some diagnostics:
3801 ! call p_maxmin('PS_model (mb)', Atm%ps(is:ie,js:je), is, ie, js, je, 1, 0.01)
3802 ! call p_maxmin('PT_model', Atm%pt(is:ie,js:je,1:npz), is, ie, js, je, npz, 1.)
3803  do j=js,je
3804  do i=is,ie
3805  wk(i,j) = phis_reg(i,j)/grav - zh(i,j,km+1)
3806  enddo
3807  enddo
3808 ! call pmaxmn('ZS_diff (m)', wk, is, ie, js, je, 1, 1., Atm%gridstruct%area_64, Atm%domain)
3809 
3810  do j=js,je
3811  do i=is,ie
3812  wk(i,j) = ps(i,j) - psc(i,j)
3813  enddo
3814  enddo
3815 ! call pmaxmn('PS_diff (mb)', wk, is, ie, js, je, 1, 0.01, Atm%gridstruct%area_64, Atm%domain)
3816  deallocate (pe0,qn1,dp2,pe1,qp)
3817  if (is_master()) write(*,*) 'done remap_scalar_nggps_regional_bc'
3818 !---------------------------------------------------------------------
3819 
3820  end subroutine remap_scalar_nggps_regional_bc
3821 
3822 !---------------------------------------------------------------------
3823 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3824 !---------------------------------------------------------------------
3825 
3826  subroutine remap_dwinds_regional_bc(Atm &
3827  ,is_input,ie_input &
3828  ,js_input,je_input &
3829  ,is_u,ie_u,js_u,je_u &
3830  ,is_v,ie_v,js_v,je_v &
3831  ,km, npz &
3832  ,ak0, bk0 &
3833  ,psc, ud, vd, uc, vc &
3834  ,BC_side )
3835  type(fv_atmos_type), intent(inout) :: Atm
3836  integer, intent(in):: is_input, ie_input, js_input, je_input
3837  integer, intent(in):: is_u,ie_u,js_u,je_u
3838  integer, intent(in):: is_v,ie_v,js_v,je_v
3839  integer, intent(in):: km & !<-- # of levels in 3-D input variables
3840  ,npz
3841  real, intent(in):: ak0(km+1), bk0(km+1)
3842 
3843  real, intent(in) :: psc(is_input:ie_input,js_input:je_input)
3844 
3845  real, intent(in):: ud(is_u:ie_u,js_u:je_u,km)
3846  real, intent(in):: vc(is_u:ie_u,js_u:je_u,km)
3847  real, intent(in):: vd(is_v:ie_v,js_v:je_v,km)
3848  real, intent(in):: uc(is_v:ie_v,js_v:je_v,km)
3849  type(fv_regional_bc_variables),intent(inout) :: BC_side
3850 ! local:
3851  real, dimension(:,:),allocatable :: pe0
3852  real, dimension(:,:),allocatable :: pe1
3853  real, dimension(:,:),allocatable :: qn1_d,qn1_c
3854  integer i,j,k
3855 
3856  allocate(pe0(is_u:ie_u, km+1)) ; pe0=real_snan
3857  allocate(pe1(is_u:ie_u, npz+1)) ; pe1=real_snan
3858  allocate(qn1_d(is_u:ie_u, npz)) ; qn1_d=real_snan
3859  allocate(qn1_c(is_u:ie_u, npz)) ; qn1_c=real_snan
3860 
3861 !----------------------------------------------------------------------------------------------
3862  j_loopu: do j=js_u,je_u
3863 !----------------------------------------------------------------------------------------------
3864 
3865 !------
3866 ! map u
3867 !------
3868  do k=1,km+1
3869  do i=is_u,ie_u
3870  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j))
3871  enddo
3872  enddo
3873  do k=1,npz+1
3874  do i=is_u,ie_u
3875  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j))
3876  enddo
3877  enddo
3878  call mappm(km, pe0(is_u:ie_u,1:km+1), ud(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), &
3879  qn1_d(is_u:ie_u,1:npz), is_u,ie_u, -1, 8, atm%ptop )
3880  call mappm(km, pe0(is_u:ie_u,1:km+1), vc(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), &
3881  qn1_c(is_u:ie_u,1:npz), is_u,ie_u, -1, 8, atm%ptop )
3882  do k=1,npz
3883  do i=is_u,ie_u
3884  bc_side%u_BC(i,j,k) = qn1_d(i,k)
3885  bc_side%vc_BC(i,j,k) = qn1_c(i,k)
3886  enddo
3887  enddo
3888 
3889  enddo j_loopu
3890 
3891  deallocate(pe0)
3892  deallocate(pe1)
3893  deallocate(qn1_d)
3894  deallocate(qn1_c)
3895 
3896  allocate(pe0(is_v:ie_v, km+1)) ; pe0=real_snan
3897  allocate(pe1(is_v:ie_v, npz+1)) ; pe1=real_snan
3898  allocate(qn1_d(is_v:ie_v, npz)) ; qn1_d=real_snan
3899  allocate(qn1_c(is_v:ie_v, npz)) ; qn1_c=real_snan
3900 
3901 !----------------------------------------------------------------------------------------------
3902  j_loopv: do j=js_v,je_v
3903 !----------------------------------------------------------------------------------------------
3904 !
3905 !------
3906 ! map v
3907 !------
3908 
3909  do k=1,km+1
3910  do i=is_v,ie_v
3911  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j))
3912  enddo
3913  enddo
3914  do k=1,npz+1
3915  do i=is_v,ie_v
3916  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j))
3917  enddo
3918  enddo
3919  call mappm(km, pe0(is_v:ie_v,1:km+1), vd(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), &
3920  qn1_d(is_v:ie_v,1:npz), is_v,ie_v, -1, 8, atm%ptop)
3921  call mappm(km, pe0(is_v:ie_v,1:km+1), uc(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), &
3922  qn1_c(is_v:ie_v,1:npz), is_v,ie_v, -1, 8, atm%ptop)
3923  do k=1,npz
3924  do i=is_v,ie_v
3925  bc_side%v_BC(i,j,k) = qn1_d(i,k)
3926  bc_side%uc_BC(i,j,k) = qn1_c(i,k)
3927  enddo
3928  enddo
3929 
3930  enddo j_loopv
3931 
3932  deallocate(pe0)
3933  deallocate(pe1)
3934  deallocate(qn1_d)
3935  deallocate(qn1_c)
3936 
3937  if (is_master()) write(*,*) 'done remap_dwinds'
3938 
3939  end subroutine remap_dwinds_regional_bc
3940 
3941 !---------------------------------------------------------------------
3942 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3943 !---------------------------------------------------------------------
3944 
3945  subroutine set_regional_bcs(delp,delz,w,pt &
3946 #ifdef USE_COND
3947  ,q_con &
3948 #endif
3949 #ifdef MOIST_CAPPA
3950  ,cappa &
3951 #endif
3952  ,q &
3953  ,u,v,uc,vc &
3954  ,bd, nlayers &
3955  ,fcst_time )
3956 !
3957 !---------------------------------------------------------------------
3958 !*** Select the boundary variables' boundary data at the two
3959 !*** bracketing time levels and apply them to the updating
3960 !*** of the variables' boundary regions at the appropriate
3961 !*** forecast time. This is done at the beginning of every
3962 !*** large timestep in fv_dynamics.
3963 !---------------------------------------------------------------------
3964  implicit none
3965 !---------------------------------------------------------------------
3966 !
3967 !--------------------
3968 !*** Input variables
3969 !--------------------
3970 !
3971  integer,intent(in) :: nlayers
3972 !
3973  real,intent(in) :: fcst_time
3974 !
3975  type(fv_grid_bounds_type),intent(in) :: bd
3976 !
3977 !----------------------
3978 !*** Output variables
3979 !----------------------
3980 !
3981  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),intent(out) :: &
3982  delp &
3983  ,pt
3984 !
3985  real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: w
3986  real,dimension(bd%is:,bd%js:,1:),intent(out) :: delz
3987 #ifdef USE_COND
3988  real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: q_con
3989 #endif
3990 
3991 !
3992  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,ntracers),intent(out) :: q
3993 !
3994 #ifdef MOIST_CAPPA
3995  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),intent(out) :: cappa
3996 !#else
3997 ! real,dimension(isd:isd,jsd:jsd,1),intent(out) :: cappa
3998 #endif
3999 !
4000  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz),intent(out) :: u,vc
4001 !
4002  real,dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz),intent(out) :: uc,v
4003 !
4004 !---------------------
4005 !*** Local variables
4006 !---------------------
4007 !
4008  real :: fraction_interval
4009 !
4010 !---------------------------------------------------------------------
4011 !*********************************************************************
4012 !---------------------------------------------------------------------
4013 !
4014 !---------------------------------------------------------------------
4015 !*** The current forecast time is this fraction of the way from
4016 !*** time level 0 to time level 1.
4017 !---------------------------------------------------------------------
4018 !
4019  fraction_interval=mod(fcst_time,(bc_update_interval*3600.)) &
4020  /(bc_update_interval*3600.)
4021 !
4022 !---------------------------------------------------------------------
4023 !
4024  if(north_bc)then
4025  call bc_values_into_arrays(bc_t0%north,bc_t1%north &
4026  ,'north' &
4027  ,bd%isd &
4028  ,bd%ied &
4029  ,bd%jsd &
4030  ,bd%js-1 &
4031  ,bd%isd &
4032  ,bd%ied &
4033  ,bd%jsd &
4034  ,bd%js-1 &
4035  ,bd%isd &
4036  ,bd%ied+1 &
4037  ,bd%jsd &
4038  ,bd%js-1)
4039  endif
4040 !
4041  if(south_bc)then
4042  call bc_values_into_arrays(bc_t0%south,bc_t1%south &
4043  ,'south' &
4044  ,bd%isd &
4045  ,bd%ied &
4046  ,bd%je+1 &
4047  ,bd%jed &
4048  ,bd%isd &
4049  ,bd%ied &
4050  ,bd%je+2 &
4051  ,bd%jed+1 &
4052  ,bd%isd &
4053  ,bd%ied+1 &
4054  ,bd%je+1 &
4055  ,bd%jed )
4056  endif
4057 !
4058  if(east_bc)then
4059  call bc_values_into_arrays(bc_t0%east,bc_t1%east &
4060  ,'east ' &
4061  ,bd%isd &
4062  ,bd%is-1 &
4063  ,bd%js &
4064  ,bd%je &
4065  ,bd%isd &
4066  ,bd%is-1 &
4067  ,bd%js &
4068  ,bd%je+1 &
4069  ,bd%isd &
4070  ,bd%is-1 &
4071  ,bd%js &
4072  ,bd%je )
4073  endif
4074 !
4075  if(west_bc)then
4076  call bc_values_into_arrays(bc_t0%west,bc_t1%west &
4077  ,'west ' &
4078  ,bd%ie+1 &
4079  ,bd%ied &
4080  ,bd%js &
4081  ,bd%je &
4082  ,bd%ie+1 &
4083  ,bd%ied &
4084  ,bd%js &
4085  ,bd%je+1 &
4086  ,bd%ie+2 &
4087  ,bd%ied+1 &
4088  ,bd%js &
4089  ,bd%je )
4090  endif
4091 !
4092 !---------------------------------------------------------------------
4093 
4094  contains
4095 
4096 !---------------------------------------------------------------------
4097 !
4098  subroutine bc_values_into_arrays(side_t0,side_t1 &
4099  ,side &
4100  ,i1,i2,j1,j2 &
4101  ,i1_uvs,i2_uvs,j1_uvs,j2_uvs &
4102  ,i1_uvw,i2_uvw,j1_uvw,j2_uvw )
4103 !
4104 !---------------------------------------------------------------------
4105 !*** Apply boundary values to the prognostic arrays at the
4106 !*** desired time.
4107 !---------------------------------------------------------------------
4108  implicit none
4109 !---------------------------------------------------------------------
4110 !
4111 !---------------------
4112 !*** Input arguments
4113 !---------------------
4114 !
4115  type(fv_regional_bc_variables),intent(in) :: side_t0 &
4116  ,side_t1
4117 !
4118  character(len=*),intent(in) :: side
4119 !
4120  integer,intent(in) :: i1,i2,j1,j2 &
4121  ,i1_uvs,i2_uvs,j1_uvs,j2_uvs &
4122  ,i1_uvw,i2_uvw,j1_uvw,j2_uvw
4123 !
4124 !---------------------
4125 !*** Local arguments
4126 !---------------------
4127 !
4128  integer :: i,ie,j,je,jend,jend_uvs,jend_uvw &
4129  ,jstart,jstart_uvs,jstart_uvw,k,nt,nz
4130 !
4131  real,dimension(:,:,:),pointer :: delz_ptr
4132 !
4133 !---------------------------------------------------------------------
4134 !*********************************************************************
4135 !---------------------------------------------------------------------
4136 !
4137  jstart=j1
4138  jend =j2
4139  jstart_uvs=j1_uvs
4140  jend_uvs =j2_uvs
4141  jstart_uvw=j1_uvw
4142  jend_uvw =j2_uvw
4143  if((trim(side)=='east'.or.trim(side)=='west').and..not.north_bc)then
4144  jstart=j1-nhalo_model
4145  jstart_uvs=j1_uvs-nhalo_model
4146  jstart_uvw=j1_uvw-nhalo_model
4147  endif
4148  if((trim(side)=='east'.or.trim(side)=='west').and..not.south_bc)then
4149  jend=j2+nhalo_model
4150  jend_uvs=j2_uvs+nhalo_model
4151  jend_uvw=j2_uvw+nhalo_model
4152  endif
4153 !
4154  select case (trim(side))
4155  case ('north')
4156  delz_ptr=>delz_auxiliary%north
4157  case ('south')
4158  delz_ptr=>delz_auxiliary%south
4159  case ('east')
4160  delz_ptr=>delz_auxiliary%east
4161  case ('west')
4162  delz_ptr=>delz_auxiliary%west
4163  end select
4164 !
4165  do k=1,nlayers
4166  do j=jstart,jend
4167  do i=i1,i2
4168  delp(i,j,k)=side_t0%delp_BC(i,j,k) &
4169  +(side_t1%delp_BC(i,j,k)-side_t0%delp_BC(i,j,k)) &
4170  *fraction_interval
4171  pt(i,j,k)=side_t0%pt_BC(i,j,k) &
4172  +(side_t1%pt_BC(i,j,k)-side_t0%pt_BC(i,j,k)) &
4173  *fraction_interval
4174 ! delz(i,j,k)=side_t0%delz_BC(i,j,k) &
4175 ! +(side_t1%delz_BC(i,j,k)-side_t0%delz_BC(i,j,k)) &
4176 ! *fraction_interval
4177  delz_ptr(i,j,k)=side_t0%delz_BC(i,j,k) &
4178  +(side_t1%delz_BC(i,j,k)-side_t0%delz_BC(i,j,k)) &
4179  *fraction_interval
4180 #ifdef MOIST_CAPPA
4181  cappa(i,j,k)=side_t0%cappa_BC(i,j,k) &
4182  +(side_t1%cappa_BC(i,j,k)-side_t0%cappa_BC(i,j,k)) &
4183  *fraction_interval
4184 #endif
4185 #ifdef USE_COND
4186  q_con(i,j,k)=side_t0%q_con_BC(i,j,k) &
4187  +(side_t1%q_con_BC(i,j,k)-side_t0%q_con_BC(i,j,k)) &
4188  *fraction_interval
4189 #endif
4190  w(i,j,k)=side_t0%w_BC(i,j,k) &
4191  +(side_t1%w_BC(i,j,k)-side_t0%w_BC(i,j,k)) &
4192  *fraction_interval
4193  enddo
4194  enddo
4195 !
4196  do j=jstart_uvs,jend_uvs
4197  do i=i1_uvs,i2_uvs
4198  u(i,j,k)=side_t0%u_BC(i,j,k) &
4199  +(side_t1%u_BC(i,j,k)-side_t0%u_BC(i,j,k)) &
4200  *fraction_interval
4201  vc(i,j,k)=side_t0%vc_BC(i,j,k) &
4202  +(side_t1%vc_BC(i,j,k)-side_t0%vc_BC(i,j,k)) &
4203  *fraction_interval
4204  enddo
4205  enddo
4206 !
4207  do j=jstart_uvw,jend_uvw
4208  do i=i1_uvw,i2_uvw
4209  v(i,j,k)=side_t0%v_BC(i,j,k) &
4210  +(side_t1%v_BC(i,j,k)-side_t0%v_BC(i,j,k)) &
4211  *fraction_interval
4212  uc(i,j,k)=side_t0%uc_BC(i,j,k) &
4213  +(side_t1%uc_BC(i,j,k)-side_t0%uc_BC(i,j,k)) &
4214  *fraction_interval
4215  enddo
4216  enddo
4217  enddo
4218 !
4219  ie=min(ubound(side_t0%w_BC,1),ubound(w,1))
4220  je=min(ubound(side_t0%w_BC,2),ubound(w,2))
4221  nz=ubound(w,3)
4222 !
4223  do nt=1,ntracers
4224  do k=1,nz
4225  do j=jstart,jend
4226  do i=i1,i2
4227  q(i,j,k,nt)=side_t0%q_BC(i,j,k,nt) &
4228  +(side_t1%q_BC(i,j,k,nt)-side_t0%q_BC(i,j,k,nt)) &
4229  *fraction_interval
4230  enddo
4231  enddo
4232  enddo
4233  enddo
4234 !
4235 !---------------------------------------------------------------------
4236 !
4237  end subroutine bc_values_into_arrays
4238 !
4239 !---------------------------------------------------------------------
4240 !
4241  end subroutine set_regional_bcs
4242 !
4243 !---------------------------------------------------------------------
4244 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4245 !---------------------------------------------------------------------
4246  subroutine regional_boundary_update(array &
4247  ,bc_vbl_name &
4248  ,lbnd_x,ubnd_x &
4249  ,lbnd_y,ubnd_y &
4250  ,ubnd_z &
4251  ,is,ie,js,je &
4252  ,isd,ied,jsd,jed &
4253  ,fcst_time &
4254  ,index4 )
4255 !
4256 !---------------------------------------------------------------------
4257 !*** Select the given variable's boundary data at the two
4258 !*** bracketing time levels and apply them to the updating
4259 !*** of the variable's boundary region at the appropriate
4260 !*** forecast time.
4261 !---------------------------------------------------------------------
4262  implicit none
4263 !---------------------------------------------------------------------
4264 !
4265 !--------------------
4266 !*** Input variables
4267 !--------------------
4268 !
4269  integer,intent(in) :: lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z
4270 !
4271  integer,intent(in) :: is,ie,js,je & !<-- Compute limits
4272  ,isd,ied,jsd,jed
4273 !
4274  integer,intent(in),optional :: index4
4275 !
4276  real,intent(in) :: fcst_time
4277 !
4278  character(len=*),intent(in) :: bc_vbl_name
4279 !
4280 !----------------------
4281 !*** Output variables
4282 !----------------------
4283 !
4284  real,dimension(lbnd_x:ubnd_x,lbnd_y:ubnd_y,1:ubnd_z) &
4285  ,intent(out) :: array
4286 !
4287 !---------------------
4288 !*** Local variables
4289 !---------------------
4290 !
4291  integer :: i1,i2,j1,j2
4292  integer :: i_bc,j_bc
4293  integer :: i1_blend,i2_blend,j1_blend,j2_blend
4294  integer :: lbnd1,ubnd1,lbnd2,ubnd2
4295  integer :: iq
4296  integer :: nside
4297 !
4298  real,dimension(:,:,:),pointer :: bc_t0,bc_t1
4299 !
4300  logical :: blend,call_interp
4301 !
4302 !---------------------------------------------------------------------
4303 !*********************************************************************
4304 !---------------------------------------------------------------------
4305 !
4306  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
4307  return
4308  endif
4309 !
4310  blend=.true.
4311  iq=0
4312  if(present(index4))then
4313  iq=index4
4314  blend=blend_this_tracer(iq)
4315  endif
4316 !
4317 !---------------------------------------------------------------------
4318 !*** Loop through the sides of the domain and find the limits
4319 !*** of the region to update in the boundary.
4320 !---------------------------------------------------------------------
4321 !
4322 !---------------------------------------------------------------------
4323  sides: do nside=1,4
4324 !---------------------------------------------------------------------
4325 !
4326  call_interp=.false.
4327 !
4328 !-----------
4329 !*** North
4330 !-----------
4331 !
4332  if(nside==1)then
4333  if(north_bc)then
4334  call_interp=.true.
4335  bc_side_t0=>bc_north_t0
4336  bc_side_t1=>bc_north_t1
4337 !
4338  i1=isd
4339  i2=ied
4340  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
4341  i2=ied+1
4342  endif
4343 !
4344  j1=jsd
4345  j2=js-1
4346 !
4347  i1_blend=is
4348  i2_blend=ie
4349  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
4350  i2_blend=ie+1
4351  endif
4352  j1_blend=js
4353  j2_blend=js+nrows_blend_user-1
4354  i_bc=-9e9
4355  j_bc=j2
4356 !
4357  endif
4358  endif
4359 !
4360 !-----------
4361 !*** South
4362 !-----------
4363 !
4364  if(nside==2)then
4365  if(south_bc)then
4366  call_interp=.true.
4367  bc_side_t0=>bc_south_t0
4368  bc_side_t1=>bc_south_t1
4369 !
4370  i1=isd
4371  i2=ied
4372  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
4373  i2=ied+1
4374  endif
4375 !
4376  j1=je+1
4377  j2=jed
4378  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4379  j1=je+2
4380  j2=jed+1
4381  endif
4382 !
4383  i1_blend=is
4384  i2_blend=ie
4385  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
4386  i2_blend=ie+1
4387  endif
4388  j2_blend=je
4389  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4390  j2_blend=je+1
4391  endif
4392  j1_blend=j2_blend-nrows_blend_user+1
4393  i_bc=-9e9
4394  j_bc=j1
4395 !
4396  endif
4397  endif
4398 !
4399 !----------
4400 !*** East
4401 !----------
4402 !
4403  if(nside==3)then
4404  if(east_bc)then
4405  call_interp=.true.
4406  bc_side_t0=>bc_east_t0
4407  bc_side_t1=>bc_east_t1
4408 !
4409  j1=jsd
4410  j2=jed
4411 !
4412 ! CHJ --- s ---
4413  if(trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='u')then
4414  j2=jed+1
4415  endif
4416 ! CHJ --- e ---
4417 
4418  i1=isd
4419  i2=is-1
4420 !
4421  if(north_bc)then
4422  j1=js
4423  endif
4424  if(south_bc)then
4425  j2=je
4426  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4427  j2=je+1
4428  endif
4429  endif
4430 !
4431  i1_blend=is
4432  i2_blend=is+nrows_blend_user-1
4433  j1_blend=js
4434  j2_blend=je
4435  if(north_bc)then
4436  j1_blend=js+nrows_blend_user
4437  endif
4438  if(south_bc)then
4439  j2_blend=je-nrows_blend_user
4440  endif
4441  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4442  j2_blend=j2_blend+1
4443  endif
4444  i_bc=i2
4445  j_bc=-9e9
4446 !
4447  endif
4448  endif
4449 !
4450 !----------
4451 !*** West
4452 !----------
4453 !
4454  if(nside==4)then
4455  if(west_bc)then
4456  call_interp=.true.
4457  bc_side_t0=>bc_west_t0
4458  bc_side_t1=>bc_west_t1
4459 !
4460  j1=jsd
4461  j2=jed
4462 
4463 ! CHJ --- s ---
4464  if(trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='u')then
4465  j2=jed+1
4466  endif
4467 ! CHJ --- e ---
4468 
4469  i1=ie+1
4470  i2=ied
4471  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
4472  i1=ie+2
4473  i2=ied+1
4474  endif
4475 !
4476  if(north_bc)then
4477  j1=js
4478  endif
4479  if(south_bc)then
4480  j2=je
4481  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4482  j2=je+1
4483  endif
4484  endif
4485 !
4486  i1_blend=i1-nrows_blend_user
4487  i2_blend=i1-1
4488  j1_blend=js
4489  j2_blend=je
4490  if(north_bc)then
4491  j1_blend=js+nrows_blend_user
4492  endif
4493  if(south_bc)then
4494  j2_blend=je-nrows_blend_user
4495  endif
4496  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
4497  j2_blend=j2_blend+1
4498  endif
4499  i_bc=i1
4500  j_bc=-9e9
4501 !
4502  endif
4503  endif
4504 !
4505 !---------------------------------------------------------------------
4506 !*** Get the pointers pointing at the boundary arrays holding the
4507 !*** two time levels of the given prognostic array's boundary region
4508 !*** then update the boundary points.
4509 !---------------------------------------------------------------------
4510 !
4511  if(call_interp)then
4512 !
4513  call retrieve_bc_variable_data(bc_vbl_name &
4514  ,bc_side_t0,bc_side_t1 &
4515  ,bc_t0,bc_t1 &
4516  ,lbnd1,ubnd1,lbnd2,ubnd2 &
4517  ,iq )
4518 !
4519  call bc_time_interpolation(array &
4520  ,lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z &
4521  ,bc_t0,bc_t1 &
4522  ,lbnd1,ubnd1,lbnd2,ubnd2 &
4523  ,i1,i2,j1,j2 &
4524  ,is,ie,js,je &
4525  ,fcst_time &
4526  ,bc_update_interval &
4527  ,i1_blend,i2_blend,j1_blend,j2_blend &
4528  ,i_bc,j_bc,nside,bc_vbl_name,blend )
4529  endif
4530 !
4531 !---------------------------------------------------------------------
4532 !
4533  enddo sides
4534 !
4535 !---------------------------------------------------------------------
4536 !
4537  end subroutine regional_boundary_update
4538 
4539 !---------------------------------------------------------------------
4540 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4541 !---------------------------------------------------------------------
4542 
4543  subroutine retrieve_bc_variable_data(bc_vbl_name &
4544  ,bc_side_t0,bc_side_t1 &
4545  ,bc_t0,bc_t1 &
4546  ,lbnd1,ubnd1,lbnd2,ubnd2 &
4547  ,iq )
4548 
4549 !---------------------------------------------------------------------
4550 !*** Select the boundary variable associated with the prognostic
4551 !*** array that needs its boundary region to be updated.
4552 !---------------------------------------------------------------------
4553  implicit none
4554 !---------------------------------------------------------------------
4555 !
4556 !---------------------
4557 !*** Input variables
4558 !---------------------
4559 !
4560  integer,intent(in) :: iq
4561 !
4562  character(len=*),intent(in) :: bc_vbl_name
4563 !
4564  type(fv_regional_bc_variables),pointer :: bc_side_t0,bc_side_t1
4565 !
4566 !
4567 !----------------------
4568 !*** Output variables
4569 !----------------------
4570 !
4571  integer,intent(out) :: lbnd1,ubnd1,lbnd2,ubnd2
4572 !
4573  real,dimension(:,:,:),pointer :: bc_t0,bc_t1
4574 !
4575 !---------------------------------------------------------------------
4576 !*********************************************************************
4577 !---------------------------------------------------------------------
4578 !
4579  select case (bc_vbl_name)
4580 !
4581  case ('delp')
4582  bc_t0=>bc_side_t0%delp_BC
4583  bc_t1=>bc_side_t1%delp_BC
4584  case ('delz')
4585  bc_t0=>bc_side_t0%delz_BC
4586  bc_t1=>bc_side_t1%delz_BC
4587  case ('pt')
4588  bc_t0=>bc_side_t0%pt_BC
4589  bc_t1=>bc_side_t1%pt_BC
4590  case ('w')
4591  bc_t0=>bc_side_t0%w_BC
4592  bc_t1=>bc_side_t1%w_BC
4593  case ('divgd')
4594  bc_t0=>bc_side_t0%divgd_BC
4595  bc_t1=>bc_side_t1%divgd_BC
4596 #ifdef MOIST_CAPPA
4597  case ('cappa')
4598  bc_t0=>bc_side_t0%cappa_BC
4599  bc_t1=>bc_side_t1%cappa_BC
4600 #endif
4601 #ifdef USE_COND
4602  case ('q_con')
4603  bc_t0=>bc_side_t0%q_con_BC
4604  bc_t1=>bc_side_t1%q_con_BC
4605 #endif
4606  case ('q')
4607  if(iq<1)then
4608  call mpp_error(fatal,' iq<1 is not a valid index for q_BC array in retrieve_bc_variable_data')
4609  endif
4610  lbnd1=lbound(bc_side_t0%q_BC,1)
4611  lbnd2=lbound(bc_side_t0%q_BC,2)
4612  ubnd1=ubound(bc_side_t0%q_BC,1)
4613  ubnd2=ubound(bc_side_t0%q_BC,2)
4614  bc_t0=>bc_side_t0%q_BC(:,:,:,iq)
4615  bc_t1=>bc_side_t1%q_BC(:,:,:,iq)
4616  case ('u')
4617  bc_t0=>bc_side_t0%u_BC
4618  bc_t1=>bc_side_t1%u_BC
4619  case ('v')
4620  bc_t0=>bc_side_t0%v_BC
4621  bc_t1=>bc_side_t1%v_BC
4622  case ('uc')
4623  bc_t0=>bc_side_t0%uc_BC
4624  bc_t1=>bc_side_t1%uc_BC
4625  case ('vc')
4626  bc_t0=>bc_side_t0%vc_BC
4627  bc_t1=>bc_side_t1%vc_BC
4628 !
4629  end select
4630 !
4631  if(trim(bc_vbl_name)/='q')then
4632  lbnd1=lbound(bc_t0,1)
4633  lbnd2=lbound(bc_t0,2)
4634  ubnd1=ubound(bc_t0,1)
4635  ubnd2=ubound(bc_t0,2)
4636  endif
4637 !
4638 !---------------------------------------------------------------------
4639 !
4640  end subroutine retrieve_bc_variable_data
4641 !
4642 !---------------------------------------------------------------------
4643 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4644 !---------------------------------------------------------------------
4645 !
4646  subroutine bc_time_interpolation(array &
4647  ,lbnd_x, ubnd_x &
4648  ,lbnd_y, ubnd_y &
4649  ,ubnd_z &
4650  ,bc_t0, bc_t1 &
4651  ,lbnd1, ubnd1 &
4652  ,lbnd2, ubnd2 &
4653  ,i1,i2,j1,j2 &
4654  ,is,ie,js,je &
4655  ,fcst_time &
4656  ,bc_update_interval &
4657  ,i1_blend,i2_blend,j1_blend,j2_blend &
4658  ,i_bc,j_bc,nside,bc_vbl_name,blend )
4659 
4660 !---------------------------------------------------------------------
4661 !*** Update the boundary region of the input array at the given
4662 !*** forecast time that is within the interval bracketed by the
4663 !*** two current boundary region states.
4664 !---------------------------------------------------------------------
4665  implicit none
4666 !---------------------------------------------------------------------
4667 !
4668 !---------------------
4669 !*** Input variables
4670 !---------------------
4671 !
4672  integer,intent(in) :: lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z
4673 !
4674  integer,intent(in) :: lbnd1,ubnd1,lbnd2,ubnd2
4675 !
4676  integer,intent(in) :: i1,i2,j1,j2 & !<-- Index limits of the updated boundary region.
4677  ,i_bc,j_bc & !<-- Innermost bndry indices (anchor pts for blending)
4678  ,i1_blend,i2_blend,j1_blend,j2_blend & !<-- Index limits of the updated blending region.
4679  ,nside
4680 !
4681  integer,intent(in) :: is,ie,js,je
4682 !
4683  integer,intent(in) :: bc_update_interval
4684 !
4685  real,intent(in) :: fcst_time
4686 !
4687  real,dimension(lbnd1:ubnd1,lbnd2:ubnd2,1:ubnd_z),intent(in) :: bc_t0 & !<-- Interpolate between these
4688  ,bc_t1 ! two boundary region states.
4689 !
4690  character(len=*),intent(in) :: bc_vbl_name
4691 !
4692  logical,intent(in) :: blend
4693 !
4694 !---------------------
4695 !*** Output variables
4696 !---------------------
4697 !
4698  real,dimension(lbnd_x:ubnd_x,lbnd_y:ubnd_y,1:ubnd_z) &
4699  ,intent(out) :: array
4700 !
4701 !---------------------
4702 !*** Local variables
4703 !---------------------
4704 !
4705  integer :: i,j,k
4706 !
4707  real :: blend_value,factor_dist,fraction_interval,rdenom
4708 !
4709 !---------------------------------------------------------------------
4710 !*********************************************************************
4711 !---------------------------------------------------------------------
4712 !
4713 !---------------------------------------------------------------------
4714 !*** The current forecast time is this fraction of the way from
4715 !*** time level 0 to time level 1.
4716 !---------------------------------------------------------------------
4717 !
4718  fraction_interval=mod(fcst_time,(bc_update_interval*3600.)) &
4719  /(bc_update_interval*3600.)
4720 !
4721 !---------------------------------------------------------------------
4722 !
4723  do k=1,ubnd_z
4724  do j=j1,j2
4725  do i=i1,i2
4726  array(i,j,k)=bc_t0(i,j,k) &
4727  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval
4728  enddo
4729  enddo
4730  enddo
4731 !
4732 !---------------------------------------------------------------------
4733 !*** If this tracer is not in the external BC file then it will not
4734 !*** be blended.
4735 !---------------------------------------------------------------------
4736 !
4737  if(.not.blend)then
4738  return
4739  endif
4740 !
4741 !---------------------------------------------------------------------
4742 !*** Use specified external data to blend with integration values
4743 !*** across nrows_blend rows immediately within the domain's
4744 !*** boundary rows. The weighting of the external data drops
4745 !*** off exponentially.
4746 !---------------------------------------------------------------------
4747 !
4748 !-----------
4749 !*** North
4750 !-----------
4751 !
4752  if(nside==1.and.north_bc)then
4753  rdenom=1./real(j2_blend-j_bc-1)
4754  do k=1,ubnd_z
4755  do j=j1_blend,j2_blend
4756  factor_dist=exp(-(blend_exp1+blend_exp2*(j-j_bc-1)*rdenom))
4757  do i=i1_blend,i2_blend
4758  blend_value=bc_t0(i,j,k) &
4759  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
4760 !
4761  array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
4762  enddo
4763  enddo
4764  enddo
4765  endif
4766 !
4767 !-----------
4768 !*** South
4769 !-----------
4770 !
4771  if(nside==2.and.south_bc)then
4772  rdenom=1./real(j_bc-j1_blend-1)
4773  do k=1,ubnd_z
4774  do j=j1_blend,j2_blend
4775  factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom))
4776  do i=i1_blend,i2_blend
4777  blend_value=bc_t0(i,j,k) &
4778  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
4779  array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
4780  enddo
4781  enddo
4782  enddo
4783  endif
4784 !
4785 !----------
4786 !*** East
4787 !----------
4788 !
4789  if(nside==3.and.east_bc)then
4790  rdenom=1./real(i2_blend-i_bc-1)
4791  do k=1,ubnd_z
4792  do j=j1_blend,j2_blend
4793  do i=i1_blend,i2_blend
4794 !
4795  blend_value=bc_t0(i,j,k) &
4796  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
4797 !
4798  factor_dist=exp(-(blend_exp1+blend_exp2*(i-i_bc-1)*rdenom))
4799 !
4800  array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
4801  enddo
4802  enddo
4803  enddo
4804  endif
4805 !
4806 !----------
4807 !*** West
4808 !----------
4809 !
4810  if(nside==4.and.west_bc)then
4811  rdenom=1./real(i_bc-i1_blend-1)
4812  do k=1,ubnd_z
4813  do j=j1_blend,j2_blend
4814  do i=i1_blend,i2_blend
4815 !
4816  blend_value=bc_t0(i,j,k) &
4817  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval ! between t0 and t1.
4818 !
4819  factor_dist=exp(-(blend_exp1+blend_exp2*(i_bc-i-1)*rdenom))
4820 !
4821  array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
4822  enddo
4823  enddo
4824  enddo
4825  endif
4826 !
4827 !---------------------------------------------------------------------
4828 !
4829  end subroutine bc_time_interpolation
4830 !
4831 !---------------------------------------------------------------------
4832 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4833 !---------------------------------------------------------------------
4834 !
4835  subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 &
4836  ,nlev,ntracers,bnds )
4837 !
4838 !---------------------------------------------------------------------
4839 !*** BC data has been read into the time level t1 object. Now
4840 !*** move the t1 data into the t1 object before refilling t1
4841 !*** with the next data from the BC file.
4842 !---------------------------------------------------------------------
4843  implicit none
4844 !---------------------------------------------------------------------
4845 !
4846 !------------------------
4847 !*** Argument variables
4848 !------------------------
4849 !
4850  integer,intent(in) :: nlev & !<-- # of model layers.
4851  ,ntracers
4852 !
4853  type(fv_regional_bc_bounds_type),intent(in) :: bnds
4854 !
4855  type(fv_domain_sides),target,intent(in) :: BC_t1
4856 !
4857  type(fv_domain_sides),target,intent(inout) :: BC_t0
4858 !
4859 !---------------------
4860 !*** Local variables
4861 !---------------------
4862 !
4863  integer :: i,ie_c,ie_s,ie_w,is_c,is_s,is_w &
4864  ,j,je_c,je_s,je_w,js_c,js_s,js_w &
4865  ,k,n,nside
4866 !
4867  logical :: move
4868 !
4869 !---------------------------------------------------------------------
4870 !*********************************************************************
4871 !---------------------------------------------------------------------
4872 !
4873 !---------------------------------------------------------------------
4874 !*** Loop through the four sides of the domain.
4875 !---------------------------------------------------------------------
4876 !
4877  sides: do nside=1,4
4878 !
4879  move=.false.
4880 !
4881  if(nside==1)then
4882  if(north_bc)then
4883  move=.true.
4884  bc_side_t0=>bc_t0%north
4885  bc_side_t1=>bc_t1%north
4886 !
4887  is_c=bnds%is_north
4888  ie_c=bnds%ie_north ! North BC index limits
4889  js_c=bnds%js_north ! for centers of grid cells.
4890  je_c=bnds%je_north
4891 !
4892  is_s=bnds%is_north_uvs
4893  ie_s=bnds%ie_north_uvs ! North BC index limits
4894  js_s=bnds%js_north_uvs ! for winds on N/S sides of grid cells.
4895  je_s=bnds%je_north_uvs
4896 !
4897  is_w=bnds%is_north_uvw
4898  ie_w=bnds%ie_north_uvw ! North BC index limits
4899  js_w=bnds%js_north_uvw ! for winds on E/W sides of grid cells.
4900  je_w=bnds%je_north_uvw
4901  endif
4902  endif
4903 !
4904  if(nside==2)then
4905  if(south_bc)then
4906  move=.true.
4907  bc_side_t0=>bc_t0%south
4908  bc_side_t1=>bc_t1%south
4909 !
4910  is_c=bnds%is_south
4911  ie_c=bnds%ie_south ! South BC index limits
4912  js_c=bnds%js_south ! for centers of grid cells.
4913  je_c=bnds%je_south
4914 !
4915  is_s=bnds%is_south_uvs
4916  ie_s=bnds%ie_south_uvs ! South BC index limits
4917  js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells.
4918  je_s=bnds%je_south_uvs
4919 !
4920  is_w=bnds%is_south_uvw
4921  ie_w=bnds%ie_south_uvw ! South BC index limits
4922  js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells.
4923  je_w=bnds%je_south_uvw
4924  endif
4925  endif
4926 !
4927  if(nside==3)then
4928  if(east_bc)then
4929  move=.true.
4930  bc_side_t0=>bc_t0%east
4931  bc_side_t1=>bc_t1%east
4932 !
4933  is_c=bnds%is_east
4934  ie_c=bnds%ie_east ! East BC index limits
4935  js_c=bnds%js_east ! for centers of grid cells.
4936  je_c=bnds%je_east
4937 !
4938  is_s=bnds%is_east_uvs
4939  ie_s=bnds%ie_east_uvs ! East BC index limits
4940  js_s=bnds%js_east_uvs ! for winds on N/S sides of grid cells.
4941  je_s=bnds%je_east_uvs
4942 !
4943  is_w=bnds%is_east_uvw
4944  ie_w=bnds%ie_east_uvw ! East BC index limits
4945  js_w=bnds%js_east_uvw ! for winds on E/W sides of grid cells.
4946  je_w=bnds%je_east_uvw
4947  endif
4948  endif
4949 !
4950  if(nside==4)then
4951  if(west_bc)then
4952  move=.true.
4953  bc_side_t0=>bc_t0%west
4954  bc_side_t1=>bc_t1%west
4955 !
4956  is_c=bnds%is_west
4957  ie_c=bnds%ie_west ! West BC index limits
4958  js_c=bnds%js_west ! for centers of grid cells.
4959  je_c=bnds%je_west
4960 !
4961  is_s=bnds%is_west_uvs
4962  ie_s=bnds%ie_west_uvs ! West BC index limits
4963  js_s=bnds%js_west_uvs ! for winds on N/S sides of grid cells.
4964  je_s=bnds%je_west_uvs
4965 !
4966  is_w=bnds%is_west_uvw
4967  ie_w=bnds%ie_west_uvw ! West BC index limits
4968  js_w=bnds%js_west_uvw ! for winds on E/W sides of grid cells.
4969  je_w=bnds%je_west_uvw
4970  endif
4971  endif
4972 !
4973  if(move)then
4974  do k=1,nlev
4975  do j=js_c,je_c
4976  do i=is_c,ie_c
4977  bc_side_t0%delp_BC(i,j,k) =bc_side_t1%delp_BC(i,j,k)
4978  bc_side_t0%divgd_BC(i,j,k)=bc_side_t1%divgd_BC(i,j,k)
4979  enddo
4980  enddo
4981  enddo
4982 !
4983  do n=1,ntracers
4984  do k=1,nlev
4985  do j=js_c,je_c
4986  do i=is_c,ie_c
4987  bc_side_t0%q_BC(i,j,k,n)=bc_side_t1%q_BC(i,j,k,n)
4988  enddo
4989  enddo
4990  enddo
4991  enddo
4992 !
4993  do k=1,nlev
4994  do j=js_c,je_c
4995  do i=is_c,ie_c
4996 #ifndef SW_DYNAMICS
4997  bc_side_t0%w_BC(i,j,k) =bc_side_t1%w_BC(i,j,k)
4998  bc_side_t0%pt_BC(i,j,k) =bc_side_t1%pt_BC(i,j,k)
4999  bc_side_t0%delz_BC(i,j,k) =bc_side_t1%delz_BC(i,j,k)
5000 #ifdef USE_COND
5001  bc_side_t0%q_con_BC(i,j,k)=bc_side_t1%q_con_BC(i,j,k)
5002 #ifdef MOIST_CAPPA
5003  bc_side_t0%cappa_BC(i,j,k)=bc_side_t1%cappa_BC(i,j,k)
5004 #endif
5005 #endif
5006 #endif
5007  enddo
5008  enddo
5009  enddo
5010 !
5011  do k=1,nlev
5012  do j=js_s,je_s
5013  do i=is_s,ie_s
5014  bc_side_t0%u_BC(i,j,k) =bc_side_t1%u_BC(i,j,k)
5015  bc_side_t0%vc_BC(i,j,k)=bc_side_t1%vc_BC(i,j,k)
5016  enddo
5017  enddo
5018  enddo
5019 !
5020  do k=1,nlev
5021  do j=js_w,je_w
5022  do i=is_w,ie_w
5023  bc_side_t0%v_BC(i,j,k) =bc_side_t1%v_BC(i,j,k)
5024  bc_side_t0%uc_BC(i,j,k)=bc_side_t1%uc_BC(i,j,k)
5025  enddo
5026  enddo
5027  enddo
5028 !
5029  endif
5030 !
5031  enddo sides
5032 !
5033 !---------------------------------------------------------------------
5034 !
5035  end subroutine regional_bc_t1_to_t0
5036 !
5037 !---------------------------------------------------------------------
5038 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5039 !---------------------------------------------------------------------
5040 !
5041  subroutine convert_to_virt_pot_temp(isd,ied,jsd,jed,npz)
5042 !
5043 !-----------------------------------------------------------------------
5044 !*** Convert the incoming sensible temperature to virtual potential
5045 !*** temperature.
5046 !-----------------------------------------------------------------------
5047  implicit none
5048 !-----------------------------------------------------------------------
5049 !
5050 !------------------------
5051 !*** Argument variables
5052 !------------------------
5053 !
5054  integer,intent(in) :: isd,ied,jsd,jed,npz
5055 !
5056 !---------------------
5057 !*** Local variables
5058 !---------------------
5059 !
5060  integer :: i1,i2,j1,j2
5061 !
5062  real :: rdg
5063 !
5064  real,dimension(:,:,:),pointer :: delp,delz,pt
5065 #ifdef USE_COND
5066  real,dimension(:,:,:),pointer :: q_con
5067 #endif
5068 #ifdef MOIST_CAPPA
5069  real,dimension(:,:,:),pointer ::cappa
5070 #endif
5071 !
5072  real,dimension(:,:,:,:),pointer :: q
5073 !
5074 !-----------------------------------------------------------------------
5075 !***********************************************************************
5076 !-----------------------------------------------------------------------
5077 !
5078  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
5079  return
5080  endif
5081 !
5082  rdg=-rdgas/grav
5083 !
5084  if(north_bc)then
5085  i1=regional_bounds%is_north
5086  i2=regional_bounds%ie_north
5087  j1=regional_bounds%js_north
5088  j2=regional_bounds%je_north
5089  q =>bc_t1%north%q_BC
5090 #ifdef USE_COND
5091  q_con=>bc_t1%north%q_con_BC
5092 #endif
5093  delp =>bc_t1%north%delp_BC
5094  delz =>bc_t1%north%delz_BC
5095 #ifdef MOIST_CAPPA
5096  cappa=>bc_t1%north%cappa_BC
5097 #endif
5098  pt =>bc_t1%north%pt_BC
5099  call compute_vpt
5100  endif
5101 !
5102  if(south_bc)then
5103  i1=regional_bounds%is_south
5104  i2=regional_bounds%ie_south
5105  j1=regional_bounds%js_south
5106  j2=regional_bounds%je_south
5107  q =>bc_t1%south%q_BC
5108 #ifdef USE_COND
5109  q_con=>bc_t1%south%q_con_BC
5110 #endif
5111  delp =>bc_t1%south%delp_BC
5112  delz =>bc_t1%south%delz_BC
5113 #ifdef MOIST_CAPPA
5114  cappa=>bc_t1%south%cappa_BC
5115 #endif
5116  pt =>bc_t1%south%pt_BC
5117  call compute_vpt
5118  endif
5119 !
5120  if(east_bc)then
5121  i1=regional_bounds%is_east
5122  i2=regional_bounds%ie_east
5123  j1=regional_bounds%js_east
5124  j2=regional_bounds%je_east
5125  q =>bc_t1%east%q_BC
5126 #ifdef USE_COND
5127  q_con=>bc_t1%east%q_con_BC
5128 #endif
5129  delp =>bc_t1%east%delp_BC
5130  delz =>bc_t1%east%delz_BC
5131 #ifdef MOIST_CAPPA
5132  cappa=>bc_t1%east%cappa_BC
5133 #endif
5134  pt =>bc_t1%east%pt_BC
5135  call compute_vpt
5136  endif
5137 !
5138  if(west_bc)then
5139  i1=regional_bounds%is_west
5140  i2=regional_bounds%ie_west
5141  j1=regional_bounds%js_west
5142  j2=regional_bounds%je_west
5143  q =>bc_t1%west%q_BC
5144 #ifdef USE_COND
5145  q_con=>bc_t1%west%q_con_BC
5146 #endif
5147  delp =>bc_t1%west%delp_BC
5148  delz =>bc_t1%west%delz_BC
5149 #ifdef MOIST_CAPPA
5150  cappa=>bc_t1%west%cappa_BC
5151 #endif
5152  pt =>bc_t1%west%pt_BC
5153  call compute_vpt
5154  endif
5155 !
5156 !-----------------------------------------------------------------------
5157 
5158  contains
5159 
5160 !-----------------------------------------------------------------------
5161 !
5162  subroutine compute_vpt
5163 !
5164 !-----------------------------------------------------------------------
5165 !*** Compute the virtual potential temperature as done in fv_dynamics.
5166 !-----------------------------------------------------------------------
5167 !
5168 !---------------------
5169 !*** Local variables
5170 !---------------------
5171 !
5172  integer :: i,j,k
5173 !
5174  real :: cvm,dp1,pkz
5175 !
5176 !-----------------------------------------------------------------------
5177 !***********************************************************************
5178 !-----------------------------------------------------------------------
5179 !
5180  do k=1,npz
5181 !
5182  do j=j1,j2
5183  do i=i1,i2
5184  dp1 = zvir*q(i,j,k,sphum_index)
5185 #ifdef USE_COND
5186 #ifdef MOIST_CAPPA
5187  cvm=(1.-q(i,j,k,sphum_index)+q_con(i,j,k))*cv_air &
5188  +q(i,j,k,sphum_index)*cv_vap+q(i,j,k,liq_water_index)*c_liq
5189  pkz=exp(cappa(i,j,k)*log(rdg*delp(i,j,k)*pt(i,j,k) &
5190  *(1.+dp1)*(1.-q_con(i,j,k))/delz(i,j,k)))
5191 #else
5192  pkz=exp(kappa*log(rdg*delp(i,j,k)*pt(i,j,k) &
5193  *(1.+dp1)*(1.-q_con(i,j,k))/delz(i,j,k)))
5194 #endif
5195  pt(i,j,k)=pt(i,j,k)*(1.+dp1)*(1.-q_con(i,j,k))/pkz
5196 #else
5197  pkz=exp(kappa*log(rdg*delp(i,j,k)*pt(i,j,k) &
5198  *(1.+dp1)/delz(i,j,k)))
5199  pt(i,j,k)=pt(i,j,k)*(1.+dp1)/pkz
5200 #endif
5201  enddo
5202  enddo
5203 !
5204  enddo
5205 !
5206 !-----------------------------------------------------------------------
5207 !
5208  end subroutine compute_vpt
5209 !
5210 !-----------------------------------------------------------------------
5211 !
5212  end subroutine convert_to_virt_pot_temp
5213 !
5214 !-----------------------------------------------------------------------
5215 !
5216 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5217 !---------------------------------------------------------------------
5218 !*** The following four subroutines are exact copies from
5219 !*** external_ic_mod. That module must USE this module therefore
5220 !*** this module cannout USE external_IC_mod to get at those
5221 !*** subroutines. The routines may be moved to their own module.
5222 !---------------------------------------------------------------------
5223 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5224 !---------------------------------------------------------------------
5225  subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
5226  character(len=*), intent(in):: qname
5227  integer, intent(in):: is, ie, js, je, km
5228  real, intent(in):: q(is:ie, js:je, km)
5229  real, intent(in):: fac
5230  real qmin, qmax
5231  integer i,j,k
5232 
5233  qmin = q(is,js,1)
5234  qmax = qmin
5235  do k=1,km
5236  do j=js,je
5237  do i=is,ie
5238  if( q(i,j,k) < qmin ) then
5239  qmin = q(i,j,k)
5240  elseif( q(i,j,k) > qmax ) then
5241  qmax = q(i,j,k)
5242  endif
5243  enddo
5244  enddo
5245  enddo
5246  call mp_reduce_min(qmin)
5247  call mp_reduce_max(qmax)
5248  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac
5249 
5250  end subroutine p_maxmin
5251 
5252 
5253  subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
5254  character(len=*), intent(in):: qname
5255  integer, intent(in):: is, ie, js, je
5256  integer, intent(in):: km
5257  real, intent(in):: q(is:ie, js:je, km)
5258  real, intent(in):: fac
5259  real(kind=R_GRID), intent(IN):: area(is-3:ie+3, js-3:je+3)
5260  type(domain2d), intent(INOUT) :: domain
5261 !---local variables
5262  real qmin, qmax, gmean
5263  integer i,j,k
5264 
5265  qmin = q(is,js,1)
5266  qmax = qmin
5267  gmean = 0.
5268 
5269  do k=1,km
5270  do j=js,je
5271  do i=is,ie
5272  if( q(i,j,k) < qmin ) then
5273  qmin = q(i,j,k)
5274  elseif( q(i,j,k) > qmax ) then
5275  qmax = q(i,j,k)
5276  endif
5277  enddo
5278  enddo
5279  enddo
5280 
5281  call mp_reduce_min(qmin)
5282  call mp_reduce_max(qmax)
5283 
5284  gmean = g_sum(domain, q(is,js,km), is, ie, js, je, 3, area, 1, reproduce=.true.)
5285  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac, gmean*fac
5286 
5287  end subroutine pmaxmn
5288 
5289 
5290  subroutine fillq(im, km, nq, q, dp)
5291  integer, intent(in):: im ! No. of longitudes
5292  integer, intent(in):: km ! No. of levels
5293  integer, intent(in):: nq ! Total number of tracers
5294  real , intent(in):: dp(im,km) ! pressure thickness
5295  real , intent(inout) :: q(im,km,nq) ! tracer mixing ratio
5296 ! !LOCAL VARIABLES:
5297  integer i, k, ic, k1
5298 
5299  do ic=1,nq
5300 ! Bottom up:
5301  do k=km,2,-1
5302  k1 = k-1
5303  do i=1,im
5304  if( q(i,k,ic) < 0. ) then
5305  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
5306  q(i,k ,ic) = 0.
5307  endif
5308  enddo
5309  enddo
5310 ! Top down:
5311  do k=1,km-1
5312  k1 = k+1
5313  do i=1,im
5314  if( q(i,k,ic) < 0. ) then
5315  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
5316  q(i,k ,ic) = 0.
5317  endif
5318  enddo
5319  enddo
5320 
5321  enddo
5322 
5323  end subroutine fillq
5324 
5325  subroutine mp_auto_conversion(ql, qr, qi, qs)
5326  real, intent(inout):: ql, qr, qi, qs
5327  real, parameter:: qi0_max = 2.0e-3
5328  real, parameter:: ql0_max = 2.5e-3
5329 
5330 ! Convert excess cloud water into rain:
5331  if ( ql > ql0_max ) then
5332  qr = ql - ql0_max
5333  ql = ql0_max
5334  endif
5335 ! Convert excess cloud ice into snow:
5336  if ( qi > qi0_max ) then
5337  qs = qi - qi0_max
5338  qi = qi0_max
5339  endif
5340 
5341  end subroutine mp_auto_conversion
5342 
5343 !-----------------------------------------------------------------------
5344 !
5345  subroutine nudge_qv_bc(Atm,isd,ied,jsd,jed)
5346 !
5347 !-----------------------------------------------------------------------
5348 !*** When nudging of specific humidity is selected then we must also
5349 !*** nudge the values in the regional boundary.
5350 !-----------------------------------------------------------------------
5351  implicit none
5352 !-----------------------------------------------------------------------
5353 !
5354 !------------------------
5355 !*** Argument variables
5356 !------------------------
5357 !
5358  integer,intent(in) :: isd,ied,jsd,jed
5359 !
5360  type(fv_atmos_type),target,intent(inout) :: Atm
5361 !
5362 !---------------------
5363 !*** Local variables
5364 !---------------------
5365 !
5366  integer :: i,i_x,ie,is,j,j_x,je,js,k
5367 !
5368  real, parameter:: q1_h2o = 2.2e-6
5369  real, parameter:: q7_h2o = 3.8e-6
5370  real, parameter:: q100_h2o = 3.8e-6
5371  real, parameter:: q1000_h2o = 3.1e-6
5372  real, parameter:: q2000_h2o = 2.8e-6
5373  real, parameter:: q3000_h2o = 3.0e-6
5374  real, parameter:: wt=2., xt=1./(1.+wt)
5375 !
5376  real :: p00,q00
5377 !
5378  type(fv_regional_bc_bounds_type),pointer :: bnds
5379 !
5380 !-----------------------------------------------------------------------
5381 !***********************************************************************
5382 !-----------------------------------------------------------------------
5383 !
5384  bnds=>atm%regional_bc_bounds
5385 !
5386 !-----------
5387 !*** North
5388 !-----------
5389 !
5390  if(north_bc)then
5391  is=lbound(bc_t1%north%q_BC,1)
5392  ie=ubound(bc_t1%north%q_BC,1)
5393  js=lbound(bc_t1%north%q_BC,2)
5394  je=ubound(bc_t1%north%q_BC,2)
5395 !
5396  i_x=isd
5397  j_x=jsd ! this location.
5398 !
5399  p00=atm%ptop
5400 !
5401  n_loopk: do k=1,npz
5402  if(p00<3000.)then
5403  call get_q00
5404  do j=js,je
5405  do i=is,ie
5406  bc_t1%north%q_BC(i,j,k,sphum_index)= &
5407  xt*(bc_t1%north%q_BC(i,j,k,sphum_index)+wt*q00)
5408  enddo
5409  enddo
5410  else
5411  exit n_loopk
5412  endif
5413  p00=p00+bc_t1%north%delp_BC(i_x,j_x,k)
5414  enddo n_loopk
5415  endif
5416 !
5417 !-----------
5418 !*** South
5419 !-----------
5420 !
5421  if(south_bc)then
5422  is=lbound(bc_t1%south%q_BC,1)
5423  ie=ubound(bc_t1%south%q_BC,1)
5424  js=lbound(bc_t1%south%q_BC,2)
5425  je=ubound(bc_t1%south%q_BC,2)
5426 !
5427  i_x=isd
5428  j_x=jed ! this location.
5429 !
5430  p00=atm%ptop
5431 !
5432  s_loopk: do k=1,npz
5433  if(p00<3000.)then
5434  call get_q00
5435  do j=js,je
5436  do i=is,ie
5437  bc_t1%south%q_BC(i,j,k,sphum_index)= &
5438  xt*(bc_t1%south%q_BC(i,j,k,sphum_index)+wt*q00)
5439  enddo
5440  enddo
5441  else
5442  exit s_loopk
5443  endif
5444  p00=p00+bc_t1%south%delp_BC(i_x,j_x,k)
5445  enddo s_loopk
5446  endif
5447 !
5448 !----------
5449 !*** East
5450 !----------
5451 !
5452  if(east_bc)then
5453  is=lbound(bc_t1%east%q_BC,1)
5454  ie=ubound(bc_t1%east%q_BC,1)
5455  js=lbound(bc_t1%east%q_BC,2)
5456  je=ubound(bc_t1%east%q_BC,2)
5457 !
5458  i_x=isd
5459  j_x=jsd+nhalo_model ! this location.
5460 !
5461  p00=atm%ptop
5462 !
5463  e_loopk: do k=1,npz
5464  if(p00<3000.)then
5465  call get_q00
5466  do j=js,je
5467  do i=is,ie
5468  bc_t1%east%q_BC(i,j,k,sphum_index)= &
5469  xt*(bc_t1%east%q_BC(i,j,k,sphum_index)+wt*q00)
5470  enddo
5471  enddo
5472  else
5473  exit e_loopk
5474  endif
5475  p00=p00+bc_t1%east%delp_BC(i_x,j_x,k)
5476  enddo e_loopk
5477  endif
5478 !
5479 !----------
5480 !*** West
5481 !----------
5482 !
5483  if(west_bc)then
5484  is=lbound(bc_t1%west%q_BC,1)
5485  ie=ubound(bc_t1%west%q_BC,1)
5486  js=lbound(bc_t1%west%q_BC,2)
5487  je=ubound(bc_t1%west%q_BC,2)
5488 !
5489  i_x=ied
5490  j_x=jsd+nhalo_model ! this location.
5491 !
5492  p00=atm%ptop
5493 !
5494  w_loopk: do k=1,npz
5495  if(p00<3000.)then
5496  call get_q00
5497  do j=js,je
5498  do i=is,ie
5499  bc_t1%west%q_BC(i,j,k,sphum_index)= &
5500  xt*(bc_t1%west%q_BC(i,j,k,sphum_index)+wt*q00)
5501  enddo
5502  enddo
5503  else
5504  exit w_loopk
5505  endif
5506  p00=p00+bc_t1%west%delp_BC(i_x,j_x,k)
5507  enddo w_loopk
5508  endif
5509 !
5510 !-----------------------------------------------------------------------
5511 !
5512  contains
5513 !
5514 !-----------------------------------------------------------------------
5515 !
5516  subroutine get_q00
5517 !
5518 !-----------------------------------------------------------------------
5519 !*** This is an internal subroutine to subroutine nudge_qv_bc that
5520 !*** computes the climatological contribution to the nudging ot the
5521 !*** input specific humidity.
5522 !-----------------------------------------------------------------------
5523 !***********************************************************************
5524 !-----------------------------------------------------------------------
5525 !
5526  if ( p00 < 30.e2 ) then
5527  if ( p00 < 1. ) then
5528  q00 = q1_h2o
5529  elseif ( p00 <= 7. .and. p00 >= 1. ) then
5530  q00 = q1_h2o + (q7_h2o-q1_h2o)*log(pref(k)/1.)/log(7.)
5531  elseif ( p00 < 100. .and. p00 >= 7. ) then
5532  q00 = q7_h2o + (q100_h2o-q7_h2o)*log(pref(k)/7.)/log(100./7.)
5533  elseif ( p00 < 1000. .and. p00 >= 100. ) then
5534  q00 = q100_h2o + (q1000_h2o-q100_h2o)*log(pref(k)/1.e2)/log(10.)
5535  elseif ( p00 < 2000. .and. p00 >= 1000. ) then
5536  q00 = q1000_h2o + (q2000_h2o-q1000_h2o)*log(pref(k)/1.e3)/log(2.)
5537  else
5538  q00 = q2000_h2o + (q3000_h2o-q2000_h2o)*log(pref(k)/2.e3)/log(1.5)
5539  endif
5540  endif
5541 !
5542 !-----------------------------------------------------------------------
5543 !
5544  end subroutine get_q00
5545 !
5546 !-----------------------------------------------------------------------
5547 !
5548  end subroutine nudge_qv_bc
5549 !
5550 !-----------------------------------------------------------------------
5551 !-----------------------------------------------------------------------
5552 
5553  subroutine dump_field_3d (domain, name, field, isd, ied, jsd, jed, nlev, stag)
5554 
5555 !-----------------------------------------------------------------------
5556 !*** Subroutines dump_field_2d and dump_field_3d are module
5557 !*** procedures with the generic interface 'dump_field'.
5558 !*** Use these routines to write out NetCDF files containing
5559 !*** FULL fields that include the variables' boundary region.
5560 !*** See the following four examples for guidance on how to
5561 !*** call the routines.
5562 !-----------------------------------------------------------------------
5563 ! call dump_field(Atm(1)%domain,"atm_pt", Atm(1)%pt, isd, ied, jsd, jed, Atm(1)%npz, stag=H_STAGGER)
5564 ! call dump_field(Atm(1)%domain,"atm_u", Atm(1)%u, isd, ied, jsd, jed+1, Atm(1)%npz, stag=U_STAGGER)
5565 ! call dump_field(Atm(1)%domain,"atm_v", Atm(1)%v, isd, ied+1, jsd, jed, Atm(1)%npz, stag=V_STAGGER)
5566 ! call dump_field(Atm(1)%domain,"atm_phis", Atm(1)%phis, isd, ied, jsd, jed, stag=H_STAGGER)
5567 
5568  type(domain2d), intent(INOUT) :: domain
5569  character(len=*), intent(IN) :: name
5570  real, dimension(isd:ied,jsd:jed,1:nlev), intent(INOUT) :: field
5571  integer, intent(IN) :: isd, ied, jsd, jed, nlev
5572  integer, intent(IN) :: stag
5573 
5574  integer :: unit
5575  character(len=128) :: fname
5576  type(axistype) :: x, y, z
5577  type(fieldtype) :: f
5578  type(domain1d) :: xdom, ydom
5579  integer :: nz
5580  integer :: is, ie, js, je
5581  integer :: isg, ieg, jsg, jeg, nxg, nyg, npx, npy
5582  integer :: i, j, halo, iext, jext
5583  logical :: is_root_pe
5584  real, allocatable, dimension(:,:,:) :: glob_field
5585  integer, allocatable, dimension(:) :: pelist
5586  character(len=1) :: stagname
5587  integer :: isection_s, isection_e, jsection_s, jsection_e
5588 
5589  write(fname,"(A,A,A,I1.1,A)") "regional_",name,".tile", 7 , ".nc"
5590  write(0,*)'dump_field_3d: file name = |', trim(fname) , '|'
5591 
5592  call mpp_get_domain_components( domain, xdom, ydom )
5593  call mpp_get_compute_domain( domain, is, ie, js, je )
5594  call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=npx, ysize=npy, position=center )
5595 
5596  halo = is - isd
5597  if ( halo /= 3 ) then
5598  write(0,*) 'dusan- halo should be 3 ', halo
5599  endif
5600 
5601  iext = 0
5602  jext = 0
5603  stagname = "h";
5604  if (stag == u_stagger) then
5605  jext = 1
5606  stagname = "u";
5607  endif
5608  if (stag == v_stagger) then
5609  iext = 1
5610  stagname = "v";
5611  endif
5612 
5613  nxg = npx + 2*halo + iext
5614  nyg = npy + 2*halo + jext
5615  nz = size(field,dim=3)
5616 
5617  allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) )
5618 
5619  isection_s = is
5620  isection_e = ie
5621  jsection_s = js
5622  jsection_e = je
5623 
5624  if ( isd < 0 ) isection_s = isd
5625  if ( ied > npx-1 ) isection_e = ied
5626  if ( jsd < 0 ) jsection_s = jsd
5627  if ( jed > npy-1 ) jsection_e = jed
5628 
5629  allocate( pelist(mpp_npes()) )
5630  call mpp_get_current_pelist(pelist)
5631 
5632  is_root_pe = (mpp_pe()==mpp_root_pe())
5633 
5634  call mpp_gather(isection_s,isection_e,jsection_s,jsection_e, nz, &
5635  pelist, field(isection_s:isection_e,jsection_s:jsection_e,:), glob_field, is_root_pe, halo, halo)
5636 
5637  call mpp_open( unit, trim(fname), action=mpp_overwr, form=mpp_netcdf, threading=mpp_single)
5638 
5639  call mpp_write_meta( unit, x, 'grid_xt', 'km', 'X distance', 'X', domain=xdom, data=(/(i*1.0,i=1,nxg)/) )
5640  call mpp_write_meta( unit, y, 'grid_yt', 'km', 'Y distance', 'Y', domain=ydom, data=(/(j*1.0,j=1,nyg)/) )
5641  call mpp_write_meta( unit, z, 'lev', 'km', 'Z distance', data=(/(i*1.0,i=1,nz)/) )
5642 
5643  call mpp_write_meta( unit, f, (/x,y,z/), name, 'unit', name)
5644  call mpp_write_meta( unit, "stretch_factor", rval=stretch_factor )
5645  call mpp_write_meta( unit, "target_lon", rval=target_lon )
5646  call mpp_write_meta( unit, "target_lat", rval=target_lat )
5647  call mpp_write_meta( unit, "cube_res", ival= cube_res)
5648  call mpp_write_meta( unit, "parent_tile", ival=parent_tile )
5649  call mpp_write_meta( unit, "refine_ratio", ival=refine_ratio )
5650  call mpp_write_meta( unit, "istart_nest", ival=istart_nest )
5651  call mpp_write_meta( unit, "jstart_nest", ival=jstart_nest )
5652  call mpp_write_meta( unit, "iend_nest", ival=iend_nest )
5653  call mpp_write_meta( unit, "jend_nest", ival=jend_nest )
5654  call mpp_write_meta( unit, "ihalo_shift", ival=halo )
5655  call mpp_write_meta( unit, "jhalo_shift", ival=halo )
5656  call mpp_write_meta( unit, mpp_get_id(f), "hstagger", cval=stagname )
5657  call mpp_write( unit, x )
5658  call mpp_write( unit, y )
5659  call mpp_write( unit, z )
5660  call mpp_write( unit, f, glob_field )
5661 
5662  call mpp_close( unit )
5663 
5664  end subroutine dump_field_3d
5665 
5666  subroutine dump_field_2d (domain, name, field, isd, ied, jsd, jed, stag)
5667 
5668  type(domain2d), intent(INOUT) :: domain
5669  character(len=*), intent(IN) :: name
5670  real, dimension(isd:ied,jsd:jed), intent(INOUT) :: field
5671  integer, intent(IN) :: isd, ied, jsd, jed
5672  integer, intent(IN) :: stag
5673 
5674  integer :: unit
5675  character(len=128) :: fname
5676  type(axistype) :: x, y
5677  type(fieldtype) :: f
5678  type(domain1d) :: xdom, ydom
5679  integer :: is, ie, js, je
5680  integer :: isg, ieg, jsg, jeg, nxg, nyg, npx, npy
5681  integer :: i, j, halo, iext, jext
5682  logical :: is_root_pe
5683  real, allocatable, dimension(:,:) :: glob_field
5684  integer, allocatable, dimension(:) :: pelist
5685  character(len=1) :: stagname
5686  integer :: isection_s, isection_e, jsection_s, jsection_e
5687 
5688  write(fname,"(A,A,A,I1.1,A)") "regional_",name,".tile", 7 , ".nc"
5689  write(0,*)'dump_field_3d: file name = |', trim(fname) , '|'
5690 
5691  call mpp_get_domain_components( domain, xdom, ydom )
5692  call mpp_get_compute_domain( domain, is, ie, js, je )
5693  call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=npx, ysize=npy, position=center )
5694 
5695  halo = is - isd
5696  if ( halo /= 3 ) then
5697  write(0,*) 'dusan- halo should be 3 ', halo
5698  endif
5699 
5700  iext = 0
5701  jext = 0
5702  stagname = "h";
5703  if (stag == u_stagger) then
5704  jext = 1
5705  stagname = "u";
5706  endif
5707  if (stag == v_stagger) then
5708  iext = 1
5709  stagname = "v";
5710  endif
5711 
5712  nxg = npx + 2*halo + iext
5713  nyg = npy + 2*halo + jext
5714 
5715  allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) )
5716 
5717  isection_s = is
5718  isection_e = ie
5719  jsection_s = js
5720  jsection_e = je
5721 
5722  if ( isd < 0 ) isection_s = isd
5723  if ( ied > npx-1 ) isection_e = ied
5724  if ( jsd < 0 ) jsection_s = jsd
5725  if ( jed > npy-1 ) jsection_e = jed
5726 
5727  allocate( pelist(mpp_npes()) )
5728  call mpp_get_current_pelist(pelist)
5729 
5730  is_root_pe = (mpp_pe()==mpp_root_pe())
5731 
5732  call mpp_gather(isection_s,isection_e,jsection_s,jsection_e, &
5733  pelist, field(isection_s:isection_e,jsection_s:jsection_e), glob_field, is_root_pe, halo, halo)
5734 
5735  call mpp_open( unit, trim(fname), action=mpp_overwr, form=mpp_netcdf, threading=mpp_single)
5736 
5737  call mpp_write_meta( unit, x, 'grid_xt', 'km', 'X distance', 'X', domain=xdom, data=(/(i*1.0,i=1,nxg)/) )
5738  call mpp_write_meta( unit, y, 'grid_yt', 'km', 'Y distance', 'Y', domain=ydom, data=(/(j*1.0,j=1,nyg)/) )
5739 
5740  call mpp_write_meta( unit, f, (/x,y/), name, 'unit', name)
5741  call mpp_write_meta( unit, "stretch_factor", rval=stretch_factor )
5742  call mpp_write_meta( unit, "target_lon", rval=target_lon )
5743  call mpp_write_meta( unit, "target_lat", rval=target_lat )
5744  call mpp_write_meta( unit, "cube_res", ival= cube_res)
5745  call mpp_write_meta( unit, "parent_tile", ival=parent_tile )
5746  call mpp_write_meta( unit, "refine_ratio", ival=refine_ratio )
5747  call mpp_write_meta( unit, "istart_nest", ival=istart_nest )
5748  call mpp_write_meta( unit, "jstart_nest", ival=jstart_nest )
5749  call mpp_write_meta( unit, "iend_nest", ival=iend_nest )
5750  call mpp_write_meta( unit, "jend_nest", ival=jend_nest )
5751  call mpp_write_meta( unit, "ihalo_shift", ival=halo )
5752  call mpp_write_meta( unit, "jhalo_shift", ival=halo )
5753  call mpp_write_meta( unit, mpp_get_id(f), "hstagger", cval=stagname )
5754  call mpp_write( unit, x )
5755  call mpp_write( unit, y )
5756  call mpp_write( unit, f, glob_field )
5757 
5758  call mpp_close( unit )
5759 
5760  end subroutine dump_field_2d
5761 
5762 !-----------------------------------------------------------------------
5763 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5764 !-----------------------------------------------------------------------
5765 !
5766  subroutine prepare_full_fields(Atm)
5767 !
5768 !-----------------------------------------------------------------------
5769 !*** Prepare the objects that will hold the names and values of
5770 !*** the core and tracer fields to be written into the expanded
5771 !*** restart files that include the boundary rows so the GSI
5772 !*** can update both the interior and BCs.
5773 !-----------------------------------------------------------------------
5774 !
5775 !------------------------
5776 !*** Argument variables
5777 !------------------------
5778 !
5779  type(fv_atmos_type),target,intent(inout) :: Atm
5780 !
5781 !---------------------
5782 !*** Local variables
5783 !---------------------
5784 !
5785  integer :: index,istat,n &
5786  ,ncid_core_new &
5787  ,ncid_tracers_new &
5788  ,ndims,nkount,nv_core,nv_tracers &
5789  ,var_id
5790 !
5791  integer :: lbnd1,lbnd2,lbnd3,ubnd1,ubnd2,ubnd3
5792 !
5793  integer,dimension(ndims_core) :: dim_lengths_core
5794 !
5795  integer,dimension(ndims_tracers) :: dim_lengths_tracers
5796 !
5797  integer,dimension(1:4) :: dimids=(/0,0,0,0/)
5798 !
5799  real,dimension(:),allocatable :: dim_values
5800 !
5801  character(len=50) :: att_name,var_name
5802 !
5803  character(len=9),dimension(ndims_core) :: dim_names_core=(/ &
5804  'xaxis_1' &
5805  ,'xaxis_2' &
5806  ,'yaxis_1' &
5807  ,'yaxis_2' &
5808  ,'zaxis_1' &
5809  ,'Time ' &
5810  /)
5811 !
5812  character(len=9),dimension(ndims_tracers) :: dim_names_tracers=(/ &
5813  'xaxis_1' &
5814  ,'yaxis_1' &
5815  ,'zaxis_1' &
5816  ,'Time ' &
5817  /)
5818 !
5819 !-----------------------------------------------------------------------
5820 !***********************************************************************
5821 !-----------------------------------------------------------------------
5822 !
5823 !-----------------------------------------------------------------------
5824 !*** The first file to be handled is the core restart file.
5825 !-----------------------------------------------------------------------
5826 !
5827 !-----------------------------------------------------------------------
5828 !*** All tasks are given pointers into the model data that will
5829 !*** be written to the new restart file. The following are the
5830 !*** prognostic variables in the core restart file. Note that
5831 !*** we must add the halo region back into DZ since we need the
5832 !*** domain boundary points for all the fields.
5833 !-----------------------------------------------------------------------
5834 !
5835  allocate(fields_core(1:nvars_core))
5836 !
5837  fields_core(1)%ptr=>atm%u
5838  fields_core(1)%name='u'
5839 !
5840  fields_core(2)%ptr=>atm%v
5841  fields_core(2)%name='v'
5842 !
5843  fields_core(3)%ptr=>atm%w
5844  fields_core(3)%name='W'
5845 !
5846  lbnd1=lbound(atm%delz,1)
5847  ubnd1=ubound(atm%delz,1)
5848  lbnd2=lbound(atm%delz,2)
5849  ubnd2=ubound(atm%delz,2)
5850  lbnd3=1
5851  ubnd3=ubound(atm%delz,3)
5852  allocate(fields_core(4)%ptr(lbnd1-nhalo_model:ubnd1+nhalo_model &
5853  ,lbnd2-nhalo_model:ubnd2+nhalo_model &
5854  ,lbnd3:ubnd3))
5855  fields_core(4)%name='DZ'
5856 !
5857  fields_core(5)%ptr=>atm%pt
5858  fields_core(5)%name='T'
5859 !
5860  fields_core(6)%ptr=>atm%delp
5861  fields_core(6)%name='delp'
5862 !
5863  allocate(fields_core(7)%ptr(lbound(atm%phis,1):ubound(atm%phis,1) &
5864  ,lbound(atm%phis,2):ubound(atm%phis,2) &
5865  ,1:1))
5866  fields_core(7)%ptr(:,:,1)=atm%phis(:,:)
5867  fields_core(7)%name='phis'
5868 !
5869 !-----------------------------------------------------------------------
5870 !*** We need to point at the tracers in the model's tracer array.
5871 !*** Those tracers depend on the physics that was selected so they
5872 !*** cannot be pre-specified like the variables in the core restart
5873 !*** file were. Read them from the expanded tracer restart file
5874 !*** that was created prior to the start for the forecast.
5875 !-----------------------------------------------------------------------
5876 !
5877  call check(nf90_open(path=filename_tracers_new &
5878  ,mode=nf90_nowrite &
5879  ,ncid=ncid_tracers_new ))
5880 !
5881  call check(nf90_inquire(ncid =ncid_tracers_new &
5882  ,nvariables=nv_tracers ))
5883 !
5884  nfields_tracers=nv_tracers-ndims_tracers
5885  allocate(fields_tracers(1:nfields_tracers),stat=istat)
5886  if(istat/=0)then
5887  call mpp_error(fatal,' Failed to allocate fields_tracers.')
5888  else
5889  if(is_master())then
5890  write(0,33012)nfields_tracers
5891 33012 format(' Allocated fields_tracers(1:',i3,')')
5892  endif
5893  endif
5894  nkount=0
5895 !
5896  do n=1,nv_tracers
5897  var_id=n
5898  call check(nf90_inquire_variable(ncid =ncid_tracers_new &
5899  ,varid=var_id &
5900  ,name =var_name ))
5901 !
5902  if(n>ndims_tracers)then
5903  nkount=nkount+1
5904  fields_tracers(nkount)%name=trim(var_name)
5905  index=get_tracer_index(model_atmos, trim(var_name))
5906  fields_tracers(nkount)%ptr=>atm%q(:,:,:, index)
5907  endif
5908 !
5909  enddo
5910 !
5911 !-----------------------------------------------------------------------
5912 !
5913  call check(nf90_close(ncid_tracers_new))
5914 !
5915 !-----------------------------------------------------------------------
5916 !
5917  end subroutine prepare_full_fields
5918 !
5919 !-----------------------------------------------------------------------
5920 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5921 !-----------------------------------------------------------------------
5922 !
5923  subroutine write_full_fields(Atm)
5924 !
5925 !-----------------------------------------------------------------------
5926 !*** Write out full fields of the primary restart variables
5927 !*** INCLUDING BOUNDARY ROWS so the GSI can include BCs in its
5928 !*** update. This is done in a restart look-alike file.
5929 !-----------------------------------------------------------------------
5930 !
5931  type(fv_atmos_type), intent(inout), target :: Atm
5932 !
5933  integer :: count_i,count_j
5934  integer :: iend,istart,jend,jstart,kend,kstart,nz
5935  integer :: iend_ptr,istart_ptr,jend_ptr,jstart_ptr
5936  integer :: iend_g,istart_g,jend_g,jstart_g
5937  integer :: ieg,iext,isg,jeg,jext,jsg,k
5938  integer :: n,ncid_core_new,ncid_tracers_new,nv,var_id
5939  integer :: halo
5940 !
5941  integer,dimension(:),allocatable :: pelist
5942 !
5943  real,dimension(:,:,:),allocatable :: global_field
5944  real,dimension(:,:,:),pointer :: field_3d
5945 !
5946  character(len=10) :: var_name
5947 !
5948  logical :: is_root_pe
5949 !
5950 !-----------------------------------------------------------------------
5951 !***********************************************************************
5952 !-----------------------------------------------------------------------
5953 !
5954  allocate( pelist(mpp_npes()) )
5955  call mpp_get_current_pelist(pelist)
5956 ! write(0,*)' pelist=',pelist
5957 !
5958  halo=nhalo_model
5959 !
5960  is_root_pe = (mpp_pe()==mpp_root_pe())
5961  if(is_root_pe)then
5962  call check(nf90_open(filename_core_new,nf90_write,ncid_core_new))
5963  write(0,*)' Opened core restart with BCs: ',trim(filename_core_new)
5964  endif
5965 !
5966 !-----------------------------------------------------------------------
5967 !*** Save the global limits of the domain and its vertical extent.
5968 !-----------------------------------------------------------------------
5969 !
5970  call mpp_get_global_domain (atm%domain, isg, ieg, jsg, jeg, position=center )
5971 !
5972 !-----------------------------------------------------------------------
5973 !*** Begin with the core restart file.
5974 !*** Loop through that file's prognostic variables.
5975 !-----------------------------------------------------------------------
5976 !
5977  vbls_core: do nv=1,nvars_core
5978 !
5979  var_name=fields_core(nv)%name
5980  if(is_root_pe)then
5981  call check(nf90_inq_varid(ncid_core_new,var_name,var_id))
5982  endif
5983 !
5984 !-----------------------------------------------------------------------
5985 !*** What is the full domain extent of this variable including
5986 !*** boundary rows?
5987 !-----------------------------------------------------------------------
5988 !
5989  iext=0
5990  jext=0
5991  if(var_name=='u'.or.var_name=='vc')then
5992  jext=1
5993  endif
5994  if(var_name=='v'.or.var_name=='uc')then
5995  iext=1
5996  endif
5997 !
5998  call mpp_get_global_domain (atm%domain, isg, ieg, jsg, jeg, position=center )
5999  istart_g=isg-halo
6000  iend_g =ieg+halo+iext
6001  jstart_g=jsg-halo
6002  jend_g =jeg+halo+jext
6003 !
6004  count_i=iend_g-istart_g+1
6005  count_j=jend_g-jstart_g+1
6006 !
6007  nz=size(fields_core(nv)%ptr,3)
6008 !
6009  allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
6010 !
6011 !-----------------------------------------------------------------------
6012 !*** What is the local extent of the variable on the task subdomain?
6013 !*** We must exclude inner halo data since the data is not updated
6014 !*** there in some of the variables. Of course the outer halo data
6015 !*** around the domain boundary is included.
6016 !-----------------------------------------------------------------------
6017 !
6018  istart=lbound(fields_core(nv)%ptr,1)
6019  if(istart>1)then
6020  istart=istart+halo
6021  endif
6022 !
6023  iend =ubound(fields_core(nv)%ptr,1)
6024  if(iend<ieg-halo)then
6025  iend=iend-halo
6026  endif
6027 !
6028  jstart=lbound(fields_core(nv)%ptr,2)
6029  if(jstart>1)then
6030  jstart=jstart+halo
6031  endif
6032 !
6033  jend =ubound(fields_core(nv)%ptr,2)
6034  if(jend<jeg-halo)then
6035  jend=jend-halo
6036  endif
6037 !
6038 !-----------------------------------------------------------------------
6039 !*** The interior values of the pt array are the sensible
6040 !*** temperature. The halo points though remain as the
6041 !*** special potential temperature used inside the dynamics
6042 !*** since those halo values never needed to be converted
6043 !*** back to sensible. We are now writing out the full
6044 !*** field including boundary rows for the GSI so the domain
6045 !*** halo points need to be converted to sensible temperature.
6046 !*** Each boundary task will now do the conversion before the
6047 !*** values are gathered onto the root task for writing out.
6048 !*** Also since the DZ array no longer has halo points we must
6049 !*** insert values back into the domain boundary rows of the
6050 !*** object holding DZ.
6051 !-----------------------------------------------------------------------
6052 !
6053  if(trim(fields_core(nv)%name)=='T' &
6054  .or. &
6055  trim(fields_core(nv)%name)=='DZ') then
6056 !
6057  call apply_delz_boundary(istart,iend,jstart,jend,nz &
6058  ,atm &
6059  ,fields_core(nv)%name &
6060  ,fields_core(nv)%ptr(istart:iend,jstart:jend,:))
6061  endif
6062 !
6063 !-----------------------------------------------------------------------
6064 !*** Gather onto a single task one layer at a time. That task
6065 !*** writes the full data to the new larger restart file.
6066 !-----------------------------------------------------------------------
6067 !
6068  do k=1,nz
6069  call mpp_gather(istart,iend,jstart,jend &
6070  ,pelist, fields_core(nv)%ptr(istart:iend,jstart:jend,k) &
6071  ,global_field(:,:,k), is_root_pe, halo, halo)
6072 !
6073  if(is_root_pe)then
6074  call check(nf90_put_var(ncid_core_new,var_id &
6075  ,global_field(:,:,k) &
6076  ,start=(/1,1,k/) &
6077  ,count=(/count_i,count_j,1/)))
6078  endif
6079  enddo
6080 !
6081  deallocate(global_field)
6082 !
6083  enddo vbls_core
6084 !
6085  if(is_root_pe)then
6086  call check(nf90_close(ncid_core_new))
6087  endif
6088 !
6089 !-----------------------------------------------------------------------
6090 !*** Now open the new tracer restart file.
6091 !-----------------------------------------------------------------------
6092 !
6093  if(is_root_pe)then
6094  call check(nf90_open(filename_tracers_new,nf90_write,ncid_tracers_new))
6095  write(0,*)' Opened tracer restart with BCs: ',trim(filename_tracers_new)
6096  endif
6097 !
6098 !-----------------------------------------------------------------------
6099 !*** What is the full domain extent of this variable including
6100 !*** boundary rows?
6101 !-----------------------------------------------------------------------
6102 !
6103  call mpp_get_global_domain (atm%domain, isg, ieg, jsg, jeg, position=center )
6104  istart_g=isg-halo
6105  iend_g =ieg+halo
6106  jstart_g=jsg-halo
6107  jend_g =jeg+halo
6108 !
6109  count_i=iend_g-istart_g+1
6110  count_j=jend_g-jstart_g+1
6111  nz=size(fields_tracers(1)%ptr,3)
6112 !
6113  allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
6114 !
6115 !-----------------------------------------------------------------------
6116 !*** What is the local extent of the variable on the task subdomain?
6117 !*** We must exclude inner halo data since the data is not updated
6118 !*** there in some of the variables. Of course the outer halo data
6119 !*** around the domain boundary is included. These values are the
6120 !*** same for all the tracers.
6121 !-----------------------------------------------------------------------
6122 !
6123 !-----------------------------------------------------------------------
6124 !*** The following are local bounds based on the global indexing.
6125 !-----------------------------------------------------------------------
6126 !
6127  lbnd_x_tracers=lbound(atm%q,1)
6128  ubnd_x_tracers=ubound(atm%q,1)
6129  lbnd_y_tracers=lbound(atm%q,2)
6130  ubnd_y_tracers=ubound(atm%q,2)
6131 !
6132  istart=lbnd_x_tracers
6133  if(istart>1)then
6134  istart=istart+halo
6135  endif
6136 !
6137  iend =ubnd_x_tracers
6138  if(iend<ieg-halo)then
6139  iend=iend-halo
6140  endif
6141 !
6142  jstart=lbnd_y_tracers
6143  if(jstart>1)then
6144  jstart=jstart+halo
6145  endif
6146 !
6147  jend =ubnd_y_tracers
6148  if(jend<jeg-halo)then
6149  jend=jend-halo
6150  endif
6151 !
6152 !-----------------------------------------------------------------------
6153 !*** The following are local bounds based on the indexing of
6154 !*** the pointers to the tracer arrays in memory. They are all
6155 !*** relative to 1 since each task pointed at the arrays as
6156 !*** ptr => Atm%q(:,:,:,sphum_index) rather than as was done
6157 !*** for the core arrays which was ptr => Atm%u .
6158 !-----------------------------------------------------------------------
6159 !
6160  istart_ptr=halo+1
6161  iend_ptr =ubnd_x_tracers-lbnd_x_tracers+1-halo
6162  jstart_ptr=halo+1
6163  jend_ptr =ubnd_y_tracers-lbnd_y_tracers+1-halo
6164 !
6165  if(north_bc)then
6166  jstart_ptr=1
6167  endif
6168  if(south_bc)then
6169  jend_ptr=ubnd_y_tracers-lbnd_y_tracers+1
6170  endif
6171  if(east_bc)then
6172  istart_ptr=1
6173  endif
6174  if(west_bc)then
6175  iend_ptr=ubnd_x_tracers-lbnd_x_tracers+1
6176  endif
6177 !
6178 !-----------------------------------------------------------------------
6179 !*** Loop through that file's prognostic tracers.
6180 !-----------------------------------------------------------------------
6181 !
6182  vbls_tracers: do nv=1,nfields_tracers
6183 !
6184  var_name=fields_tracers(nv)%name
6185  if(is_root_pe)then
6186  call check(nf90_inq_varid(ncid_tracers_new,var_name,var_id))
6187  endif
6188 !
6189 !-----------------------------------------------------------------------
6190 !*** Gather onto a single task one layer at a time. That task
6191 !*** writes the full data to the new larger restart file.
6192 !-----------------------------------------------------------------------
6193 !
6194  do k=1,nz
6195  call mpp_gather(istart,iend,jstart,jend &
6196  ,pelist, fields_tracers(nv)%ptr(istart_ptr:iend_ptr,jstart_ptr:jend_ptr,k) &
6197  ,global_field(:,:,k), is_root_pe, halo, halo)
6198 !
6199  if(is_root_pe)then
6200  call check(nf90_put_var(ncid_tracers_new,var_id &
6201  ,global_field(:,:,k) &
6202  ,start=(/1,1,k/) &
6203  ,count=(/count_i,count_j,1/)))
6204  endif
6205  enddo
6206 !
6207  enddo vbls_tracers
6208 !
6209  deallocate(global_field)
6210 !
6211  if(is_root_pe)then
6212  call check(nf90_close(ncid_tracers_new))
6213  endif
6214 !
6215 !---------------------------------------------------------------------
6216  end subroutine write_full_fields
6217 !---------------------------------------------------------------------
6218 !---------------------------------------------------------------------
6219 !
6220  subroutine apply_delz_boundary(istart,iend,jstart,jend,nz &
6221  ,Atm &
6222  ,name &
6223  ,field)
6224 !
6225 !---------------------------------------------------------------------
6226 !*** Use the current boundary values of delz to convert the
6227 !*** boundary potential temperature to sensible temperature
6228 !*** and to fill in the boundary rows of the 3D delz array.
6229 !---------------------------------------------------------------------
6230 !
6231 !------------------------
6232 !*** Argument variables
6233 !------------------------
6234 !
6235  integer,intent(in) :: istart,iend,jstart,jend
6236 !
6237  character(len=*),intent(in) :: name
6238 !
6239  type(fv_atmos_type),intent(inout) :: Atm
6240 !
6241  real,dimension(istart:iend,jstart:jend,1:nz),intent(inout) :: field
6242 !
6243 !---------------------
6244 !*** Local variables
6245 !---------------------
6246 !
6247  integer :: i1,i2,j1,j2,nz
6248  integer :: lbnd1,lbnd2,ubnd1,ubnd2,i,j,k
6249 !
6250  real :: rdg
6251 !
6252  real,dimension(:,:,:),pointer :: delz_ptr
6253 !
6254 !---------------------------------------------------------------------
6255 !*********************************************************************
6256 !---------------------------------------------------------------------
6257 !
6258 !fill interior delz points before dealing with boundaries
6259 !
6260 !---------------------------------------------------------------------
6261 !*** Fill the interior of the full delz array using Atm%delz
6262 !*** which does not have a boundary.
6263 !---------------------------------------------------------------------
6264 !
6265  if (trim(name)=='DZ') then
6266  lbnd1=lbound(atm%delz,1)
6267  ubnd1=ubound(atm%delz,1)
6268  lbnd2=lbound(atm%delz,2)
6269  ubnd2=ubound(atm%delz,2)
6270 !
6271  do k=1,nz
6272  do j=lbnd2,ubnd2
6273  do i=lbnd1,ubnd1
6274  field(i,j,k)=atm%delz(i,j,k)
6275  enddo
6276  enddo
6277  enddo
6278  endif
6279 !
6280  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
6281  return
6282  endif
6283 !
6284  rdg=-rdgas/grav
6285 !
6286 !---------------------------------------------------------------------
6287 !
6288  if(north_bc)then
6289  i1=istart
6290  i2=iend
6291  j1=jstart
6292  j2=jstart+nhalo_model-1
6293  delz_ptr=>delz_auxiliary%north
6294 !
6295  if(trim(name)=='T')then
6296  call compute_halo_t
6297  elseif(trim(name)=='DZ')then
6298  call fill_delz
6299  endif
6300 !
6301  endif
6302 !
6303  if(south_bc)then
6304  i1=istart
6305  i2=iend
6306  j1=jend-nhalo_model+1
6307  j2=jend
6308  delz_ptr=>delz_auxiliary%south
6309 !
6310  if(trim(name)=='T')then
6311  call compute_halo_t
6312  elseif(trim(name)=='DZ')then
6313  call fill_delz
6314  endif
6315 !
6316  endif
6317 !
6318  if(east_bc)then
6319  i1=istart
6320  i2=istart+nhalo_model-1
6321  j1=jstart
6322  j2=jend
6323  if(north_bc)then
6324  j1=jstart+nhalo_model
6325  elseif(south_bc)then
6326  j2=jend-nhalo_model
6327  endif
6328  delz_ptr=>delz_auxiliary%east
6329 !
6330  if(trim(name)=='T')then
6331  call compute_halo_t
6332  elseif(trim(name)=='DZ')then
6333  call fill_delz
6334  endif
6335 !
6336  endif
6337 !
6338  if(west_bc)then
6339  i1=iend-nhalo_model+1
6340  i2=iend
6341  j1=jstart
6342  j2=jend
6343  if(north_bc)then
6344  j1=jstart+nhalo_model
6345  elseif(south_bc)then
6346  j2=jend-nhalo_model
6347  endif
6348  delz_ptr=>delz_auxiliary%west
6349 !
6350  if(trim(name)=='T')then
6351  call compute_halo_t
6352  elseif(trim(name)=='DZ')then
6353  call fill_delz
6354  endif
6355 !
6356  endif
6357 !
6358 !---------------------------------------------------------------------
6359  contains
6360 !---------------------------------------------------------------------
6361 !
6362  subroutine compute_halo_t
6363 !
6364 !---------------------------------------------------------------------
6365 !
6366  integer :: i,j,k
6367 !
6368  real :: cappa,cvm,dp1,part1,part2
6369 !
6370 !---------------------------------------------------------------------
6371 !*********************************************************************
6372 !---------------------------------------------------------------------
6373 !
6374  do k=1,nz
6375  do j=j1,j2
6376  do i=i1,i2
6377  dp1 = zvir*atm%q(i,j,k,sphum_index)
6378  cvm=(1.-atm%q(i,j,k,sphum_index)+atm%q_con(i,j,k))*cv_air &
6379  +atm%q(i,j,k,sphum_index)*cv_vap &
6380  +atm%q(i,j,k,liq_water_index)*c_liq
6381  cappa=rdgas/(rdgas+cvm/(1.+dp1))
6382 !
6383  part1=(1.+dp1)*(1.-atm%q_con(i,j,k))
6384  part2=rdg*atm%delp(i,j,k)*(1.+dp1)*(1.-atm%q_con(i,j,k)) &
6385  /delz_ptr(i,j,k)
6386  field(i,j,k)=exp((log(field(i,j,k))-log(part1)+cappa*log(part2)) &
6387  /(1.-cappa))
6388  enddo
6389  enddo
6390  enddo
6391 !
6392 !---------------------------------------------------------------------
6393  end subroutine compute_halo_t
6394 !---------------------------------------------------------------------
6395 !
6396  subroutine fill_delz
6397 !
6398 !---------------------------------------------------------------------
6399 !
6400  integer :: i,j,k
6401  integer :: lbnd1,lbnd2,ubnd1,ubnd2
6402 !
6403 !---------------------------------------------------------------------
6404 !*********************************************************************
6405 !---------------------------------------------------------------------
6406 !
6407 !---------------------------------------------------------------------
6408 !*** Now fill the boundary rows using data from the BC files.
6409 !---------------------------------------------------------------------
6410 !
6411  do k=1,nz
6412  do j=j1,j2
6413  do i=i1,i2
6414  field(i,j,k)=delz_ptr(i,j,k)
6415  enddo
6416  enddo
6417  enddo
6418 !
6419 !---------------------------------------------------------------------
6420  end subroutine fill_delz
6421 !---------------------------------------------------------------------
6422 !
6423  end subroutine apply_delz_boundary
6424 !
6425 !---------------------------------------------------------------------
6426 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6427 !---------------------------------------------------------------------
6428 
6429  subroutine exch_uv(domain, bd, npz, u, v)
6430  use mpi
6431 
6432  implicit none
6433 
6434  type(domain2d), intent(inout) :: domain
6435  type(fv_grid_bounds_type), intent(in) :: bd
6436  integer, intent(in) :: npz
6437  real, intent(inout) :: u (bd%isd:bd%ied ,bd%jsd:bd%jed+1,1:npz)
6438  real, intent(inout) :: v (bd%isd:bd%ied+1,bd%jsd:bd%jed ,1:npz)
6439 
6440  real, dimension(:), allocatable :: buf1,buf2,buf3,buf4
6441  integer :: ihandle1,ihandle2,ihandle3,ihandle4
6442  integer,dimension(MPI_STATUS_SIZE) :: istat
6443  integer :: ic, i, j, k, is, ie, js, je
6444  integer :: irecv, isend, ierr
6445 
6446  integer :: mype
6447  integer :: north_pe, south_pe, east_pe, west_pe
6448 
6449  mype = mpp_pe()
6450  call mpp_get_neighbor_pe( domain, north, north_pe)
6451  call mpp_get_neighbor_pe( domain, south, south_pe)
6452  call mpp_get_neighbor_pe( domain, west, west_pe)
6453  call mpp_get_neighbor_pe( domain, east, east_pe)
6454 
6455  ! write(0,*) ' north_pe = ', north_pe
6456  ! write(0,*) ' south_pe = ', south_pe
6457  ! write(0,*) ' west_pe = ', west_pe
6458  ! write(0,*) ' east_pe = ', east_pe
6459 
6460  is=bd%is
6461  ie=bd%ie
6462  js=bd%js
6463  je=bd%je
6464 
6465  ! The size of these buffers must match the number of indices
6466  ! required below to send/receive the data. In particular,
6467  ! buf1 and buf4 must be of the same size (sim. for buf2 and buf3).
6468  ! Changes to the code below should be tested with debug flags
6469  ! enabled (out-of-bounds reads/writes).
6470  allocate(buf1(1:24*npz))
6471  allocate(buf2(1:36*npz))
6472  allocate(buf3(1:36*npz))
6473  allocate(buf4(1:24*npz))
6474 
6475 ! FIXME: MPI_COMM_WORLD
6476 
6477 #ifdef OVERLOAD_R4
6478 #define _DYN_MPI_REAL MPI_REAL
6479 #else
6480 #define _DYN_MPI_REAL MPI_DOUBLE_PRECISION
6481 #endif
6482 
6483 ! Receive from north
6484  if( north_pe /= null_pe )then
6485  call mpi_irecv(buf1,size(buf1),_dyn_mpi_real,north_pe,north_pe &
6486  ,mpi_comm_world,ihandle1,irecv)
6487  endif
6488 
6489 ! Receive from south
6490  if( south_pe /= null_pe )then
6491  call mpi_irecv(buf2,size(buf2),_dyn_mpi_real,south_pe,south_pe &
6492  ,mpi_comm_world,ihandle2,irecv)
6493  endif
6494 
6495 ! Send to north
6496  if( north_pe /= null_pe )then
6497  ic=0
6498  do k=1,npz
6499 
6500  do j=je-3+1,je-1+1
6501  do i=is-3,is-1
6502  ic=ic+1
6503  buf3(ic)=u(i,j,k)
6504  enddo
6505  do i=ie+1,ie+3
6506  ic=ic+1
6507  buf3(ic)=u(i,j,k)
6508  enddo
6509  enddo
6510 
6511  do j=je-2,je
6512  do i=is-3,is-1
6513  ic=ic+1
6514  buf3(ic)=v(i,j,k)
6515  enddo
6516  do i=ie+1,ie+3
6517  ic=ic+1
6518  buf3(ic)=v(i,j,k)
6519  enddo
6520  enddo
6521  enddo
6522  if (ic/=size(buf2).or.ic/=size(buf3)) &
6523  call mpp_error(fatal,'Buffer sizes buf2 and buf3 in routine exch_uv do not match actual message size')
6524  call mpi_issend(buf3,size(buf3),_dyn_mpi_real,north_pe,mype &
6525  ,mpi_comm_world,ihandle3,isend)
6526  endif
6527 
6528 ! Send to south
6529  if( south_pe /= null_pe )then
6530  ic=0
6531  do k=1,npz
6532 
6533  do j=js+2,js+3
6534  do i=is-3,is-1
6535  ic=ic+1
6536  buf4(ic)=u(i,j,k)
6537  enddo
6538  do i=ie+1,ie+3
6539  ic=ic+1
6540  buf4(ic)=u(i,j,k)
6541  enddo
6542  enddo
6543 
6544  do j=js+1,js+2
6545  do i=is-3,is-1
6546  ic=ic+1
6547  buf4(ic)=v(i,j,k)
6548  enddo
6549  do i=ie+1,ie+3
6550  ic=ic+1
6551  buf4(ic)=v(i,j,k)
6552  enddo
6553  enddo
6554 
6555  enddo
6556  if (ic/=size(buf1).or.ic/=size(buf4)) &
6557  call mpp_error(fatal,'Buffer sizes buf1 and buf4 in routine exch_uv do not match actual message size')
6558  call mpi_issend(buf4,size(buf4),_dyn_mpi_real,south_pe,mype &
6559  ,mpi_comm_world,ihandle4,isend)
6560  endif
6561 
6562 ! Store from south
6563  if( south_pe /= null_pe )then
6564  ic=0
6565  call mpi_wait(ihandle2,istat,ierr)
6566  do k=1,npz
6567 
6568  do j=js-3,js-1
6569  do i=is-3,is-1
6570  ic=ic+1
6571  u(i,j,k)=buf2(ic)
6572  enddo
6573  do i=ie+1,ie+3
6574  ic=ic+1
6575  u(i,j,k)=buf2(ic)
6576  enddo
6577  enddo
6578 
6579  do j=js-3,js-1
6580  do i=is-3,is-1
6581  ic=ic+1
6582  v(i,j,k)=buf2(ic)
6583  enddo
6584  do i=ie+1,ie+3
6585  ic=ic+1
6586  v(i,j,k)=buf2(ic)
6587  enddo
6588  enddo
6589 
6590  enddo
6591  endif
6592 
6593 ! Store from north
6594  if( north_pe /= null_pe )then
6595  ic=0
6596  call mpi_wait(ihandle1,istat,ierr)
6597  do k=1,npz
6598 
6599  do j=je+2+1,je+3+1
6600  do i=is-3,is-1
6601  ic=ic+1
6602  u(i,j,k)=buf1(ic)
6603  enddo
6604  do i=ie+1,ie+3
6605  ic=ic+1
6606  u(i,j,k)=buf1(ic)
6607  enddo
6608  enddo
6609 
6610  do j=je+2,je+3
6611  do i=is-3,is-1
6612  ic=ic+1
6613  v(i,j,k)=buf1(ic)
6614  enddo
6615  do i=ie+1,ie+3
6616  ic=ic+1
6617  v(i,j,k)=buf1(ic)
6618  enddo
6619  enddo
6620 
6621  enddo
6622  endif
6623 
6624  deallocate(buf1)
6625  deallocate(buf2)
6626  deallocate(buf3)
6627  deallocate(buf4)
6628 
6629  end subroutine exch_uv
6630 
6631 !---------------------------------------------------------------------
6632 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6633 !---------------------------------------------------------------------
6634 
6635  subroutine get_data_source(source,regional)
6636 !
6637 ! This routine extracts the data source information if it is present in the datafile.
6638 !
6639  character (len = 80) :: source
6640  integer :: ncids,sourceLength
6641  logical :: lstatus,regional
6642 !
6643 ! Use the fms call here so we can actually get the return code value.
6644 !
6645  if (regional) then
6646  lstatus = get_global_att_value('INPUT/gfs_data.nc',"source", source)
6647  else
6648  lstatus = get_global_att_value('INPUT/gfs_data.tile1.nc',"source", source)
6649  endif
6650  if (.not. lstatus) then
6651  if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute'
6652  source='No Source Attribute'
6653  endif
6654  end subroutine get_data_source
6655 
6656 !---------------------------------------------------------------------
6657 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6658 !---------------------------------------------------------------------
6659 
6660  subroutine set_delp_and_tracers(BC_side,npz,nwat)
6661 !
6662 ! This routine mimics what is done in external_ic to add mass back to delp
6663 ! and remove it from the tracers
6664 !
6665  integer :: npz,nwat
6666  type(fv_regional_bc_variables),intent(inout) :: BC_side
6667 !
6668 ! local variables
6669 !
6670  integer :: k, j, i, iq, is, ie, js, je
6671  integer :: liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
6672  real :: qt, wt, m_fac
6673 
6674  is=lbound(bc_side%delp_BC,1)
6675  ie=ubound(bc_side%delp_BC,1)
6676  js=lbound(bc_side%delp_BC,2)
6677  je=ubound(bc_side%delp_BC,2)
6678 !
6679  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
6680  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
6681  rainwat = get_tracer_index(model_atmos, 'rainwat')
6682  snowwat = get_tracer_index(model_atmos, 'snowwat')
6683  graupel = get_tracer_index(model_atmos, 'graupel')
6684  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
6685 !
6686  source: if (trim(data_source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
6687 !
6688 ! if (cld_amt > 0) BC_side%q_BC(:,:,:,cld_amt) = 0.0 ! Moorthi
6689  do k=1,npz
6690  do j=js,je
6691  do i=is,ie
6692  wt = bc_side%delp_BC(i,j,k)
6693  if ( nwat == 6 ) then
6694  qt = wt*(1. + bc_side%q_BC(i,j,k,liq_wat) + &
6695  bc_side%q_BC(i,j,k,ice_wat) + &
6696  bc_side%q_BC(i,j,k,rainwat) + &
6697  bc_side%q_BC(i,j,k,snowwat) + &
6698  bc_side%q_BC(i,j,k,graupel))
6699  else ! all other values of nwat
6700  qt = wt*(1. + sum(bc_side%q_BC(i,j,k,2:nwat)))
6701  endif
6702 !--- Adjust delp with tracer mass.
6703  bc_side%delp_BC(i,j,k) = qt
6704  enddo
6705  enddo
6706  enddo
6707 !
6708  else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO FILE
6709 !
6710 ! 20160928: Adjust the mixing ratios consistently...
6711  do k=1,npz
6712  do j=js,je
6713  do i=is,ie
6714  wt = bc_side%delp_BC(i,j,k)
6715  if ( nwat == 6 ) then
6716  qt = wt*(1. + bc_side%q_BC(i,j,k,liq_wat) + &
6717  bc_side%q_BC(i,j,k,ice_wat) + &
6718  bc_side%q_BC(i,j,k,rainwat) + &
6719  bc_side%q_BC(i,j,k,snowwat) + &
6720  bc_side%q_BC(i,j,k,graupel))
6721  else ! all other values of nwat
6722  qt = wt*(1. + sum(bc_side%q_BC(i,j,k,2:nwat)))
6723  endif
6724  m_fac = wt / qt
6725  do iq=1,ntracers
6726  bc_side%q_BC(i,j,k,iq) = m_fac * bc_side%q_BC(i,j,k,iq)
6727  enddo
6728  bc_side%delp_BC(i,j,k) = qt
6729  enddo
6730  enddo
6731  enddo
6732 !
6733  endif source
6734 !
6735  end subroutine set_delp_and_tracers
6736 
6737 !---------------------------------------------------------------------
6738 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6739 !---------------------------------------------------------------------
6740 
6741 end module fv_regional_mod
6742 
6743 !---------------------------------------------------------------------
fv_regional_mod::dump_field_2d
subroutine dump_field_2d(domain, name, field, isd, ied, jsd, jed, stag)
Definition: fv_regional_bc.F90:5667
fv_regional_mod::npz
integer, save npz
Definition: fv_regional_bc.F90:126
fv_regional_mod::filename_tracers_new
character(len=50) filename_tracers_new
Definition: fv_regional_bc.F90:165
fv_regional_mod::jend_nest
integer, parameter jend_nest
Definition: fv_regional_bc.F90:107
fv_regional_mod::iend_nest
integer, parameter iend_nest
Definition: fv_regional_bc.F90:106
fv_grid_utils_mod::mid_pt_sphere
subroutine, public mid_pt_sphere(p1, p2, pm)
Definition: fv_grid_utils.F90:2020
fv_diagnostics_mod::prt_height
subroutine, public prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
Definition: fv_diagnostics.F90:4119
fv_regional_mod::lbnd_x_tracers
integer, save lbnd_x_tracers
Definition: fv_regional_bc.F90:142
fv_regional_mod::k_step
integer, public k_step
Definition: fv_regional_bc.F90:249
fv_regional_mod::c_ice
real, parameter c_ice
Definition: fv_regional_bc.F90:225
fv_regional_mod::dump_field_3d
subroutine dump_field_3d(domain, name, field, isd, ied, jsd, jed, nlev, stag)
Definition: fv_regional_bc.F90:5554
fv_arrays_mod::fv_regional_bc_bounds_type
Definition: fv_arrays.F90:1197
fv_arrays_mod::fv_nest_bc_type_3d
Definition: fv_arrays.F90:1043
fv_regional_mod::begin_regional_restart
logical, save, public begin_regional_restart
Definition: fv_regional_bc.F90:157
fv_regional_mod::ns
integer ns
Definition: fv_regional_bc.F90:220
fv_regional_mod::regional_boundary_update
subroutine, public regional_boundary_update(array, bc_vbl_name, lbnd_x, ubnd_x, lbnd_y, ubnd_y, ubnd_z, is, ie, js, je, isd, ied, jsd, jed, fcst_time, index4)
Definition: fv_regional_bc.F90:4255
fv_regional_mod::current_time_in_seconds
real, public current_time_in_seconds
Definition: fv_regional_bc.F90:122
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
fv_regional_mod::fields_core
type(vars_3d), dimension(:), allocatable fields_core
Definition: fv_regional_bc.F90:214
fv_regional_mod::refine_ratio
integer, parameter refine_ratio
Definition: fv_regional_bc.F90:101
fv_regional_mod::regional_bc_t1_to_t0
subroutine, public regional_bc_t1_to_t0(BC_t1, BC_t0, nlev, ntracers, bnds)
Definition: fv_regional_bc.F90:4837
fv_regional_mod::ied_mod
integer, save ied_mod
Definition: fv_regional_bc.F90:124
fv_regional_mod::blend_this_tracer
logical, dimension(:), allocatable, save blend_this_tracer
Definition: fv_regional_bc.F90:160
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1187
fv_regional_mod::h_stagger
integer, parameter, public h_stagger
Definition: fv_regional_bc.F90:92
fv_regional_mod::cld_amt_index
integer, save cld_amt_index
--
Definition: fv_regional_bc.F90:133
fv_regional_mod::p_step
integer, public p_step
Definition: fv_regional_bc.F90:249
fv_regional_mod::bc_update_interval
integer, save bc_update_interval
Definition: fv_regional_bc.F90:247
fv_regional_mod::setup_regional_bc
subroutine, public setup_regional_bc(Atm, dt_atmos, isd, ied, jsd, jed, npx, npy)
Definition: fv_regional_bc.F90:259
fv_regional_mod::grid_data
character(len=100) grid_data
Definition: fv_regional_bc.F90:232
fv_regional_mod::stretch_factor
real, parameter stretch_factor
Definition: fv_regional_bc.F90:97
fv_mapz_mod::moist_cv
subroutine, public moist_cv(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cvm, t1)
The subroutine &#39;moist_cv&#39; computes the FV3-consistent moist heat capacity under constant volume...
Definition: fv_mapz.F90:3418
fv_mapz_mod::mappm
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
fv_regional_mod::bc_east_t1
type(fv_regional_bc_variables), pointer, save bc_east_t1
Definition: fv_regional_bc.F90:201
fv_regional_mod::regional_bc_data
subroutine, public regional_bc_data(Atm, bc_hour, is, ie, js, je, isd, ied, jsd, jed, ak, bk)
Definition: fv_regional_bc.F90:1514
fv_regional_mod::start_regional_cold_start
subroutine, public start_regional_cold_start(Atm, dt_atmos, ak, bk, levp, is, ie, js, je, isd, ied, jsd, jed)
Definition: fv_regional_bc.F90:1223
fv_regional_mod::nfields_tracers
integer, save nfields_tracers
Definition: fv_regional_bc.F90:126
fv_regional_mod::mp_auto_conversion
subroutine mp_auto_conversion(ql, qr, qi, qs)
Definition: fv_regional_bc.F90:5326
fv_regional_mod::remap_scalar_nggps_regional_bc
subroutine remap_scalar_nggps_regional_bc(Atm, side, isd, ied, jsd, jed, is_bc, ie_bc, js_bc, je_bc, km, npz, ncnst, ak0, bk0, psc, t_in, qa, omga, zh, phis_reg, ps, BC_side)
Definition: fv_regional_bc.F90:3402
fv_regional_mod::cube_res
integer, parameter cube_res
Definition: fv_regional_bc.F90:103
fv_regional_mod::vars_3d
Definition: fv_regional_bc.F90:194
fv_regional_mod::o3mr_index
integer, save o3mr_index
Definition: fv_regional_bc.F90:133
fv_regional_mod::filename_core
character(len=50) filename_core
Definition: fv_regional_bc.F90:162
fv_eta_mod::get_eta_level
subroutine, public get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
The subroutine &#39;get_eta_level&#39; returns the interface and layer-mean pressures for reference...
Definition: fv_eta.F90:1833
fv_regional_mod::tice
real, parameter tice
Definition: fv_regional_bc.F90:222
fv_regional_mod::dump_field
Definition: fv_regional_bc.F90:242
fv_regional_mod::bc_t1
type(fv_domain_sides), target, save, public bc_t1
– Boundary values for all BC variables at successive times from the regional BC file ...
Definition: fv_regional_bc.F90:199
fv_grid_utils_mod::g_sum
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;.
Definition: fv_grid_utils.F90:2968
fv_regional_mod::liq_water_index
integer, save liq_water_index
Definition: fv_regional_bc.F90:133
fv_regional_mod::n_step
integer, public n_step
Definition: fv_regional_bc.F90:249
fv_regional_mod::north_bc
logical, save north_bc
Definition: fv_regional_bc.F90:157
fv_regional_mod::ndims_core
integer, parameter ndims_core
-- # of core restart dimensions
Definition: fv_regional_bc.F90:115
fv_regional_mod::start_regional_restart
subroutine, public start_regional_restart(Atm, dt_atmos, isc, iec, jsc, jec, isd, ied, jsd, jed)
Definition: fv_regional_bc.F90:1319
fv_grid_utils_mod::inner_prod
real function, public inner_prod(v1, v2)
Definition: fv_grid_utils.F90:1039
fv_regional_mod::check
subroutine check(status)
Definition: fv_regional_bc.F90:3290
fv_regional_mod::convert_to_virt_pot_temp
subroutine convert_to_virt_pot_temp(isd, ied, jsd, jed, npz)
Definition: fv_regional_bc.F90:5042
fv_grid_utils_mod::get_latlon_vector
subroutine, public get_latlon_vector(pp, elon, elat)
Definition: fv_grid_utils.F90:3341
fv_grid_utils_mod::get_unit_vect2
subroutine, public get_unit_vect2(e1, e2, uc)
Definition: fv_grid_utils.F90:1887
fv_regional_mod::nudge_qv_bc
subroutine nudge_qv_bc(Atm, isd, ied, jsd, jed)
Definition: fv_regional_bc.F90:5346
fv_regional_mod::jsd_mod
integer, save jsd_mod
Definition: fv_regional_bc.F90:124
fv_arrays_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays.F90:34
fv_regional_mod::ndims_tracers
integer, parameter ndims_tracers
– # of tracer restart dimensions
Definition: fv_regional_bc.F90:115
fv_regional_mod::next_time_to_read_bcs
integer, save, public next_time_to_read_bcs
Definition: fv_regional_bc.F90:126
fv_arrays_mod::allocate_fv_nest_bc_type
&#39;allocate_fv_nest_BC_type&#39; is an interface to subroutines that allocate the &#39;fv_nest_BC_type&#39; structu...
Definition: fv_arrays.F90:1176
fv_regional_mod::parent_tile
integer, parameter parent_tile
Definition: fv_regional_bc.F90:100
read_regional_lon_lat
subroutine read_regional_lon_lat
Definition: fv_regional_bc.F90:1048
compute_halo_t
subroutine compute_halo_t
Definition: fv_regional_bc.F90:6363
fv_regional_mod::bc_south_t1
type(fv_regional_bc_variables), pointer, save bc_south_t1
Definition: fv_regional_bc.F90:201
fv_regional_mod::fv_regional_bc_variables
Definition: fv_regional_bc.F90:167
fv_regional_mod::grid_reg
real(kind=r_grid), dimension(:,:,:), allocatable grid_reg
– Lon/lat of cell corners
Definition: fv_regional_bc.F90:150
fv_regional_mod::filename_core_new
character(len=50) filename_core_new
Definition: fv_regional_bc.F90:163
fv_regional_mod::read_regional_bc_file
subroutine read_regional_bc_file(is_input, ie_input, js_input, je_input, nlev, ntracers, var_name_root, array_3d, array_4d, tlev, required)
Definition: fv_regional_bc.F90:3013
fv_regional_mod::bc_time_interval
integer, parameter, public bc_time_interval
Definition: fv_regional_bc.F90:88
fv_regional_mod::bc_south_t0
type(fv_regional_bc_variables), pointer, save bc_south_t0
Definition: fv_regional_bc.F90:201
fv_regional_mod::bc_t0
type(fv_domain_sides), target, save, public bc_t0
Definition: fv_regional_bc.F90:199
fv_regional_mod::dum1d
real, dimension(:), allocatable dum1d
Definition: fv_regional_bc.F90:230
boundary_mod
The module &#39;boundary&#39; contains utility routines for grid nesting and boundary conditions.
Definition: boundary.F90:25
fv_regional_mod::bc_side_t0
type(fv_regional_bc_variables), pointer bc_side_t0
Definition: fv_regional_bc.F90:210
fv_regional_mod::nhalo_data
integer, parameter nhalo_data
Definition: fv_regional_bc.F90:88
fv_regional_mod::data_source
character(len=80) data_source
Definition: fv_regional_bc.F90:251
fv_regional_mod::k_split
integer, save k_split
Definition: fv_regional_bc.F90:129
fv_regional_mod::real_snan
real, parameter real_snan
Definition: fv_regional_bc.F90:238
fv_regional_mod::prepare_full_fields
subroutine prepare_full_fields(Atm)
Definition: fv_regional_bc.F90:5767
fv_regional_mod::jed_mod
integer, save jed_mod
Definition: fv_regional_bc.F90:124
fv_regional_mod::set_regional_bcs
subroutine, public set_regional_bcs(delp, delz, w, pt ifdef USE_COND
Definition: fv_regional_bc.F90:3947
fv_regional_mod::bc_north_t1
type(fv_regional_bc_variables), pointer, save bc_north_t1
Definition: fv_regional_bc.F90:201
fv_regional_mod::west_bc
logical, save west_bc
Definition: fv_regional_bc.F90:157
fv_regional_mod::bk_in
real, dimension(:), allocatable, public bk_in
Definition: fv_regional_bc.F90:155
fv_mapz_mod
The module &#39;fv_mapz&#39; contains the vertical mapping routines .
Definition: fv_mapz.F90:27
fv_regional_mod::zvir
real, parameter zvir
Definition: fv_regional_bc.F90:226
fv_regional_mod::pref
real, dimension(:), allocatable pref
Definition: fv_regional_bc.F90:230
fv_regional_mod::target_lon
real, parameter target_lon
Definition: fv_regional_bc.F90:98
fv_regional_mod::ncid
integer, save ncid
Definition: fv_regional_bc.F90:126
fv_regional_mod::ice_water_index
integer, save ice_water_index
Definition: fv_regional_bc.F90:133
fv_regional_mod::bc_time_interpolation
subroutine bc_time_interpolation(array, lbnd_x, ubnd_x, lbnd_y, ubnd_y, ubnd_z, bc_t0, bc_t1, lbnd1, ubnd1, lbnd2, ubnd2, i1, i2, j1, j2, is, ie, js, je, fcst_time, bc_update_interval, i1_blend, i2_blend, j1_blend, j2_blend, i_bc, j_bc, nside, bc_vbl_name, blend)
Definition: fv_regional_bc.F90:4659
fv_fill_mod::fillz
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
fv_regional_mod::dt_atmos
real, save dt_atmos
-- The physics (large) timestep (sec)
Definition: fv_regional_bc.F90:147
fv_regional_mod::apply_delz_boundary
subroutine apply_delz_boundary(istart, iend, jstart, jend, nz, Atm, name, field)
Definition: fv_regional_bc.F90:6224
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
get_q00
subroutine get_q00
Definition: fv_regional_bc.F90:5517
fv_regional_mod::allocate_regional_bc_arrays
subroutine allocate_regional_bc_arrays(side, north_bc, south_bc, east_bc, west_bc, is_0, ie_0, js_0, je_0, is_sn, ie_sn, js_sn, je_sn, is_we, ie_we, js_we, je_we, klev, ntracers, BC_side, delz_side)
Definition: fv_regional_bc.F90:3313
fv_regional_mod::p_maxmin
subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
Definition: fv_regional_bc.F90:5226
fv_regional_mod::delz_auxiliary
type(single_vbl3d_sides) delz_auxiliary
– Boundary delz that follows integration through forecast time.
Definition: fv_regional_bc.F90:218
fv_regional_mod::t_i0
real, parameter t_i0
Definition: fv_regional_bc.F90:222
fv_regional_mod::graupel_index
integer, save graupel_index
Definition: fv_regional_bc.F90:133
fv_regional_mod::ubnd_y_tracers
integer, save ubnd_y_tracers
– Local upper bounds of x,y for tracer arrays
Definition: fv_regional_bc.F90:142
fv_grid_utils_mod::cell_center2
subroutine, public cell_center2(q1, q2, q3, q4, e2)
Definition: fv_grid_utils.F90:2726
bc_values_into_arrays
subroutine bc_values_into_arrays(side_t0, side_t1, side, i1, i2, j1, j2, i1_uvs, i2_uvs, j1_uvs, j2_uvs, i1_uvw, i2_uvw, j1_uvw, j2_uvw)
Definition: fv_regional_bc.F90:4103
fv_eta_mod
The module &#39;fv_eta&#39; contains routine to set up the reference (Eulerian) pressure coordinate.
Definition: fv_eta.F90:25
fv_regional_mod::single_vbl3d_sides
Definition: fv_regional_bc.F90:185
fv_regional_mod::oro_data
character(len=100) oro_data
Definition: fv_regional_bc.F90:232
fv_regional_mod::isd_mod
integer, save isd_mod
Definition: fv_regional_bc.F90:124
fv_regional_mod::nvars_core
integer, parameter nvars_core
-- # of prognostic variables in core restart file
Definition: fv_regional_bc.F90:115
fv_regional_mod::fields_tracers
type(vars_3d), dimension(:), allocatable fields_tracers
Definition: fv_regional_bc.F90:214
fv_regional_mod::fv_domain_sides
Definition: fv_regional_bc.F90:181
fv_regional_mod::bc_east_t0
type(fv_regional_bc_variables), pointer, save bc_east_t0
Definition: fv_regional_bc.F90:201
fv_regional_mod
Definition: fv_regional_bc.F90:22
fv_regional_mod::istart_nest
integer, parameter istart_nest
Definition: fv_regional_bc.F90:104
fv_mapz_mod::moist_cp
subroutine, public moist_cp(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cpm, t1)
The subroutine &#39;moist_cp&#39; computes the FV3-consistent moist heat capacity under constant pressure...
Definition: fv_mapz.F90:3539
fv_regional_mod::south_bc
logical, save south_bc
Definition: fv_regional_bc.F90:157
fv_regional_mod::ubnd_x_tracers
integer, save ubnd_x_tracers
Definition: fv_regional_bc.F90:142
fv_regional_mod::delz_regbc
type(fv_nest_bc_type_3d), public delz_regbc
Definition: fv_regional_bc.F90:216
fv_regional_mod::get_data_source
subroutine, public get_data_source(source, regional)
Definition: fv_regional_bc.F90:6636
compute_regional_bc_indices
subroutine compute_regional_bc_indices(regional_bc_bounds)
Definition: fv_regional_bc.F90:776
fv_regional_mod::bc_side_t1
type(fv_regional_bc_variables), pointer bc_side_t1
Definition: fv_regional_bc.F90:210
fv_regional_mod::sphum_index
integer, save sphum_index
–
Definition: fv_regional_bc.F90:133
fv_regional_mod::regional_bounds
type(fv_regional_bc_bounds_type), pointer, save regional_bounds
Definition: fv_regional_bc.F90:212
fill_divgd_bc
subroutine fill_divgd_bc
Definition: fv_regional_bc.F90:2677
fv_regional_mod::phis_reg
real, dimension(:,:), allocatable phis_reg
– Filtered sfc geopotential
Definition: fv_regional_bc.F90:153
fv_regional_mod::set_delp_and_tracers
subroutine set_delp_and_tracers(BC_side, npz, nwat)
Definition: fv_regional_bc.F90:6661
fv_regional_mod::blend_exp1
real, parameter blend_exp1
Definition: fv_regional_bc.F90:119
fv_regional_mod::v_stagger
integer, parameter, public v_stagger
Definition: fv_regional_bc.F90:94
fv_regional_mod::lbnd_y_tracers
integer, save lbnd_y_tracers
-- Local lower bounds of x,y for tracer arrays
Definition: fv_regional_bc.F90:142
fv_regional_mod::c_liq
real, parameter c_liq
Definition: fv_regional_bc.F90:224
fv_diagnostics_mod
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
Definition: fv_diagnostics.F90:24
fv_regional_mod::cv_vap
real, parameter cv_vap
Definition: fv_regional_bc.F90:226
fv_regional_mod::snow_water_index
integer, save snow_water_index
Definition: fv_regional_bc.F90:133
fv_grid_utils_mod
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
Definition: fv_grid_utils.F90:26
fv_regional_mod::dbl_snan
real(kind=r_grid), parameter dbl_snan
Definition: fv_regional_bc.F90:240
fv_regional_mod::bc_hour
integer, save, public bc_hour
Definition: fv_regional_bc.F90:131
fv_regional_mod::dyn_timestep
real, save dyn_timestep
– The dynamics timestep (sec)
Definition: fv_regional_bc.F90:147
fv_regional_mod::bc_west_t0
type(fv_regional_bc_variables), pointer, save bc_west_t0
Definition: fv_regional_bc.F90:201
fv_regional_mod::bc_north_t0
type(fv_regional_bc_variables), pointer, save bc_north_t0
Definition: fv_regional_bc.F90:201
fv_regional_mod::write_full_fields
subroutine, public write_full_fields(Atm)
Definition: fv_regional_bc.F90:5924
fv_diagnostics_mod::prt_gb_nh_sh
subroutine, public prt_gb_nh_sh(qname, is, ie, js, je, a2, area, lat)
Definition: fv_diagnostics.F90:4168
fv_regional_mod::target_lat
real, parameter target_lat
Definition: fv_regional_bc.F90:99
fv_regional_mod::nrows_blend_user
integer, save nrows_blend_user
Definition: fv_regional_bc.F90:247
fv_regional_mod::agrid_reg
real(kind=r_grid), dimension(:,:,:), allocatable agrid_reg
-- Lon/lat of cell centers
Definition: fv_regional_bc.F90:150
fv_regional_mod::remap_dwinds_regional_bc
subroutine remap_dwinds_regional_bc(Atm, is_input, ie_input, js_input, je_input, is_u, ie_u, js_u, je_u, is_v, ie_v, js_v, je_v, km, npz, ak0, bk0, psc, ud, vd, uc, vc, BC_side)
Definition: fv_regional_bc.F90:3835
compute_vpt
subroutine compute_vpt
Definition: fv_regional_bc.F90:5163
fv_arrays_mod::fv_atmos_type
Definition: fv_arrays.F90:1216
fv_regional_mod::vars_2d
Definition: fv_regional_bc.F90:189
fv_regional_mod::jstart_nest
integer, parameter jstart_nest
Definition: fv_regional_bc.F90:105
fv_regional_mod::pmaxmn
subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
Definition: fv_regional_bc.F90:5254
fv_regional_mod::retrieve_bc_variable_data
subroutine retrieve_bc_variable_data(bc_vbl_name, bc_side_t0, bc_side_t1, bc_t0, bc_t1, lbnd1, ubnd1, lbnd2, ubnd2, iq)
Definition: fv_regional_bc.F90:4548
fv_regional_mod::a_step
integer, public a_step
Definition: fv_regional_bc.F90:249
fv_regional_mod::ak_in
real, dimension(:), allocatable, public ak_in
Definition: fv_regional_bc.F90:155
fv_regional_mod::ntimesteps_per_bc_update
integer, save, public ntimesteps_per_bc_update
Definition: fv_regional_bc.F90:131
fv_regional_mod::east_bc
logical, save east_bc
Definition: fv_regional_bc.F90:157
fv_regional_mod::nhalo_model
integer, parameter nhalo_model
Definition: fv_regional_bc.F90:88
fv_regional_mod::exch_uv
subroutine, public exch_uv(domain, bd, npz, u, v)
Definition: fv_regional_bc.F90:6430
fv_regional_mod::nrows_blend
integer, save nrows_blend
– # of blending rows in the BC data files.
Definition: fv_regional_bc.F90:145
read_regional_filtered_topo
subroutine read_regional_filtered_topo
Definition: fv_regional_bc.F90:1147
fv_regional_mod::rain_water_index
integer, save rain_water_index
Definition: fv_regional_bc.F90:133
fill_bc_for_da
subroutine fill_bc_for_da
Definition: fv_regional_bc.F90:2399
fv_regional_mod::ntracers
integer, save ntracers
Definition: fv_regional_bc.F90:126
fv_regional_mod::n_split
integer, save n_split
Definition: fv_regional_bc.F90:129
fv_regional_mod::filename_tracers
character(len=50) filename_tracers
Definition: fv_regional_bc.F90:164
fv_regional_mod::blend_exp2
real, parameter blend_exp2
– Define the exponential dropoff of weights
Definition: fv_regional_bc.F90:119
fv_regional_mod::fillq
subroutine fillq(im, km, nq, q, dp)
Definition: fv_regional_bc.F90:5291
fv_regional_mod::bc_west_t1
type(fv_regional_bc_variables), pointer, save bc_west_t1
Definition: fv_regional_bc.F90:201
fv_regional_mod::cv_air
real, parameter cv_air
Definition: fv_regional_bc.F90:226
fv_fill_mod
Definition: fv_fill.F90:21
fv_regional_mod::read_new_bc_data
subroutine, public read_new_bc_data(Atm, Time, Time_step_atmos, p_split, isd, ied, jsd, jed)
Definition: fv_regional_bc.F90:1431
fill_delz
subroutine fill_delz
Definition: fv_regional_bc.F90:6397
fv_regional_mod::u_stagger
integer, parameter, public u_stagger
Definition: fv_regional_bc.F90:93