FV3DYCORE  Version1.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 
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, 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 &
49  ,r_grid
54  use fv_mapz_mod, only: mappm, moist_cp, moist_cv
55  use fv_mp_mod, only: is_master, mp_reduce_min, mp_reduce_max
56  use fv_fill_mod, only: fillz
57  use fv_eta_mod, only: get_eta_level
58  use fms_mod, only: check_nml_error,file_exist
59  use fms_io_mod, only: read_data,get_global_att_value
60 
61  implicit none
62 
63  private
64 
65  public ak_in, bk_in &
66  ,bc_hour &
67  ,bc_t0,bc_t1 &
79  ,dump_field &
82 
83  integer,parameter :: nhalo_data =4 &
84  ,nhalo_model=3
85 
86  integer, public, parameter :: h_stagger = 1
87  integer, public, parameter :: u_stagger = 2
88  integer, public, parameter :: v_stagger = 3
89 
90  real, parameter :: stretch_factor = 1.50
91  real, parameter :: target_lon = -97.5
92  real, parameter :: target_lat = 35.5
93  integer, parameter :: parent_tile = 6
94  integer, parameter :: refine_ratio = 3
95 
96  integer, parameter :: cube_res = 96
97  integer, parameter :: istart_nest = 26
98  integer, parameter :: jstart_nest = 36
99  integer, parameter :: iend_nest = 167
100  integer, parameter :: jend_nest = 165
101 
102 ! integer, parameter :: cube_res = 768
103 ! integer, parameter :: istart_nest = 191
104 ! integer, parameter :: jstart_nest = 327
105 ! integer, parameter :: iend_nest = 1346
106 ! integer, parameter :: jend_nest = 1290
107 
110  integer,save :: liq_water_index,sphum_index
112 
113  real(kind=R_GRID),dimension(:,:,:),allocatable :: agrid_reg & !<-- Lon/lat of cell centers
114  ,grid_reg
115 
116  real,dimension(:,:),allocatable :: phis_reg
117 
118  real,dimension(:),allocatable :: ak_in, bk_in
119 
120  logical,save :: north_bc,south_bc,east_bc,west_bc &
121  ,begin_regional_restart=.true.
122 
124  real,dimension(:,:,:),allocatable :: delp_bc, divgd_bc, u_bc, v_bc, uc_bc, vc_bc
125  real,dimension(:,:,:,:),allocatable :: q_bc
126 #ifndef SW_DYNAMICS
127  real,dimension(:,:,:),allocatable :: pt_bc, w_bc, delz_bc
128 #ifdef USE_COND
129  real,dimension(:,:,:),allocatable :: q_con_bc
130 #ifdef MOIST_CAPPA
131  real,dimension(:,:,:),allocatable :: cappa_bc
132 #endif
133 #endif
134 #endif
135  end type fv_regional_bc_variables
136 
138  type(fv_regional_bc_variables) :: north, south, east, west
139  end type fv_domain_sides
140 
141  type(fv_domain_sides),target,save :: bc_t0, bc_t1
142 
143  type(fv_regional_bc_variables),pointer,save :: bc_north_t0 &
144  ,bc_south_t0 &
145  ,bc_west_t0 &
146  ,bc_east_t0 &
147  ,bc_north_t1 &
148  ,bc_south_t1 &
149  ,bc_west_t1 &
150  ,bc_east_t1
151 
153 
155 
156  real,parameter :: tice=273.16 &
157  ,t_i0=15.
158  real, parameter :: c_liq = 4185.5 ! gfdl: heat capacity of liquid at 15 deg c
159  real, parameter :: c_ice = 1972.0 ! gfdl: heat capacity of ice at - 15 deg c
160  real, parameter :: zvir = rvgas/rdgas - 1. &
161  ,cv_air = cp_air - rdgas &
162  ,cv_vap = cp_vapor - rvgas
163 
164  real,dimension(:),allocatable :: dum1d, pref
165 
166  character(len=100) :: grid_data='grid.tile7.halo4.nc' &
167  ,oro_data ='oro_data.tile7.halo4.nc'
168 
169 #ifdef OVERLOAD_R4
170  real, parameter:: real_snan=x'FFBFFFFF'
171 #else
172  real, parameter:: real_snan=x'FFF7FFFFFFFFFFFF'
173 #endif
174  real(kind=R_GRID), parameter:: dbl_snan=x'FFF7FFFFFFFFFFFF'
175 
176  interface dump_field
177  module procedure dump_field_3d
178  module procedure dump_field_2d
179  end interface dump_field
180 
181  integer,save :: bc_update_interval
182 
183  integer :: a_step, p_step, k_step, n_step
184 !
185  character(len=80) :: data_source
186 contains
187 
188 !-----------------------------------------------------------------------
189 !
190  subroutine setup_regional_bc(Atm &
191  ,isd,ied,jsd,jed &
192  ,npx,npy )
193 !
194 !-----------------------------------------------------------------------
195 !*** Regional boundary data is obtained from the external BC file.
196 !-----------------------------------------------------------------------
197  use netcdf
198 !-----------------------------------------------------------------------
199  implicit none
200 !-----------------------------------------------------------------------
201 !
202 !---------------------
203 !*** Input variables
204 !---------------------
205 !
206  integer,intent(in) :: isd,ied,jsd,jed,npx,npy
207 !
208  type(fv_atmos_type),target,intent(inout) :: atm
209 !
210 !--------------------
211 !*** Local variables
212 !--------------------
213 !
214  integer :: i,i_start,i_end,j,j_start,j_end,klev_out
215 !
216  real :: ps1
217 !
218 !-----------------------------------------------------------------------
219 !***********************************************************************
220 !-----------------------------------------------------------------------
221 !
222 !-----------------------------------------------------------------------
223 !*** The boundary data is laid out so that the pieces for the north
224 !*** and south sides span the entire distance from the east side of
225 !*** of the east halo to the west side of the west halo. Therefore
226 !*** there the # of cells in the x direction in the north/south BC
227 !*** data is nx+2*nhalo where nx is the # of cells in the x direction
228 !*** on the compute domain. This means the # of cells spanned in the
229 !*** west/east side BC data is just ny (the # of cells in the y
230 !*** direction on the compute domain) and not ny+2*nhalo since the
231 !*** halo values on the south and north ends of the east/west sides
232 !*** are already part of the BC data on the north/south sides.
233 !-----------------------------------------------------------------------
234 !
235 ! nhalo_model=3
236 !
237 ! |----------- nxp-1 -----------| <-- east/west compute points
238 ! |---------- north BC data ----------|
239 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
240 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
241 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
242 ! --- ooo ---j=1--- ooo --- ---
243 ! | ooo ooo | |
244 ! | ooo |ooo | |
245 ! ooo i=1-->|ooo
246 ! west BC data ooo| |ooo east BC data nyp-1 <-- north/south compute points
247 ! ooo|<--i=isd-nhalo_model ooo
248 ! | ooo| ooo | |
249 ! | ooo ooo | |
250 ! --- ooo ---j=jsd-nhalo_model--- ooo --- ---
251 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
252 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
253 ! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
254 ! |---------- south BC data ----------|
255 !
256 !-----------------------------------------------------------------------
257 !
258  north_bc=.false.
259  south_bc=.false.
260  east_bc =.false.
261  west_bc =.false.
262 !
263 !-----------------------------------------------------------------------
264 !*** Which side(s) of the domain does this task lie on if any?
265 !-----------------------------------------------------------------------
266 !
267  if(jsd<0)then
268  north_bc=.true.
269  endif
270 
271  if(jed>npy-1)then
272  south_bc=.true.
273  endif
274 
275  if(isd<0)then
276  east_bc=.true.
277  endif
278 
279  if(ied>npx-1)then
280  west_bc=.true.
281  endif
282 !
283  bc_update_interval=atm%flagstruct%bc_update_interval
284 !
285  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
286  return
287  endif
288 !
289 !
290 !-----------------------------------------------------------------------
291 !
292  ntracers=atm%ncnst - atm%flagstruct%dnats
293  npz=atm%npz
294  klev_out=npz
295 !
296  regional_bounds=>atm%regional_bc_bounds
297 !
298 !-----------------------------------------------------------------------
299 !*** Compute the index limits within the boundary region on each
300 !*** side of the domain for both scalars and winds. Since the
301 !*** domain does not move then the computations need to be done
302 !*** only once. Likewise find and save the locations of the
303 !*** available tracers in the tracers array.
304 !-----------------------------------------------------------------------
305 !
306  call compute_regional_bc_indices(atm%regional_bc_bounds)
307 !
308  liq_water_index = get_tracer_index(model_atmos, 'liq_wat')
309  sphum_index = get_tracer_index(model_atmos, 'sphum')
310 !
311 !-----------------------------------------------------------------------
312 !*** Allocate the objects that will hold the boundary variables
313 !*** at the two time levels surrounding each piece of the regional
314 !*** domain's integration. Data is read from the BC files into
315 !*** time level t1 while t0 holds the data from the preceding
316 !*** BC file.
317 !-----------------------------------------------------------------------
318 !*** Point pointers at each side's boundary data for both time levels.
319 !*** Those are needed when the actual update of boundary points is
320 !*** executed.
321 !-----------------------------------------------------------------------
322 !
323  if(north_bc)then
324  call allocate_regional_bc_arrays('north' &
325  ,north_bc,south_bc &
326  ,east_bc,west_bc &
327  ,atm%regional_bc_bounds%is_north &
328  ,atm%regional_bc_bounds%ie_north &
329  ,atm%regional_bc_bounds%js_north &
330  ,atm%regional_bc_bounds%je_north &
331  ,atm%regional_bc_bounds%is_north_uvs &
332  ,atm%regional_bc_bounds%ie_north_uvs &
333  ,atm%regional_bc_bounds%js_north_uvs &
334  ,atm%regional_bc_bounds%je_north_uvs &
335  ,atm%regional_bc_bounds%is_north_uvw &
336  ,atm%regional_bc_bounds%ie_north_uvw &
337  ,atm%regional_bc_bounds%js_north_uvw &
338  ,atm%regional_bc_bounds%je_north_uvw &
339  ,klev_out &
340  ,ntracers &
341  ,bc_t1%north )
342 !
343  call allocate_regional_bc_arrays('north' &
344  ,north_bc,south_bc &
345  ,east_bc,west_bc &
346  ,atm%regional_bc_bounds%is_north &
347  ,atm%regional_bc_bounds%ie_north &
348  ,atm%regional_bc_bounds%js_north &
349  ,atm%regional_bc_bounds%je_north &
350  ,atm%regional_bc_bounds%is_north_uvs &
351  ,atm%regional_bc_bounds%ie_north_uvs &
352  ,atm%regional_bc_bounds%js_north_uvs &
353  ,atm%regional_bc_bounds%je_north_uvs &
354  ,atm%regional_bc_bounds%is_north_uvw &
355  ,atm%regional_bc_bounds%ie_north_uvw &
356  ,atm%regional_bc_bounds%js_north_uvw &
357  ,atm%regional_bc_bounds%je_north_uvw &
358  ,klev_out &
359  ,ntracers &
360  ,bc_t0%north )
361 !
362  bc_north_t0=>bc_t0%north
363  bc_north_t1=>bc_t1%north
364 !
365  endif
366 
367  if(south_bc)then
368  call allocate_regional_bc_arrays('south' &
369  ,north_bc,south_bc &
370  ,east_bc,west_bc &
371  ,atm%regional_bc_bounds%is_south &
372  ,atm%regional_bc_bounds%ie_south &
373  ,atm%regional_bc_bounds%js_south &
374  ,atm%regional_bc_bounds%je_south &
375  ,atm%regional_bc_bounds%is_south_uvs &
376  ,atm%regional_bc_bounds%ie_south_uvs &
377  ,atm%regional_bc_bounds%js_south_uvs &
378  ,atm%regional_bc_bounds%je_south_uvs &
379  ,atm%regional_bc_bounds%is_south_uvw &
380  ,atm%regional_bc_bounds%ie_south_uvw &
381  ,atm%regional_bc_bounds%js_south_uvw &
382  ,atm%regional_bc_bounds%je_south_uvw &
383  ,klev_out &
384  ,ntracers &
385  ,bc_t1%south )
386 !
387  call allocate_regional_bc_arrays('south' &
388  ,north_bc,south_bc &
389  ,east_bc,west_bc &
390  ,atm%regional_bc_bounds%is_south &
391  ,atm%regional_bc_bounds%ie_south &
392  ,atm%regional_bc_bounds%js_south &
393  ,atm%regional_bc_bounds%je_south &
394  ,atm%regional_bc_bounds%is_south_uvs &
395  ,atm%regional_bc_bounds%ie_south_uvs &
396  ,atm%regional_bc_bounds%js_south_uvs &
397  ,atm%regional_bc_bounds%je_south_uvs &
398  ,atm%regional_bc_bounds%is_south_uvw &
399  ,atm%regional_bc_bounds%ie_south_uvw &
400  ,atm%regional_bc_bounds%js_south_uvw &
401  ,atm%regional_bc_bounds%je_south_uvw &
402  ,klev_out &
403  ,ntracers &
404  ,bc_t0%south )
405 !
406  bc_south_t0=>bc_t0%south
407  bc_south_t1=>bc_t1%south
408 !
409  endif
410 !
411  if(east_bc)then
412  call allocate_regional_bc_arrays('east ' &
413  ,north_bc,south_bc &
414  ,east_bc,west_bc &
415  ,atm%regional_bc_bounds%is_east &
416  ,atm%regional_bc_bounds%ie_east &
417  ,atm%regional_bc_bounds%js_east &
418  ,atm%regional_bc_bounds%je_east &
419  ,atm%regional_bc_bounds%is_east_uvs &
420  ,atm%regional_bc_bounds%ie_east_uvs &
421  ,atm%regional_bc_bounds%js_east_uvs &
422  ,atm%regional_bc_bounds%je_east_uvs &
423  ,atm%regional_bc_bounds%is_east_uvw &
424  ,atm%regional_bc_bounds%ie_east_uvw &
425  ,atm%regional_bc_bounds%js_east_uvw &
426  ,atm%regional_bc_bounds%je_east_uvw &
427  ,klev_out &
428  ,ntracers &
429  ,bc_t1%east )
430 !
431  call allocate_regional_bc_arrays('east ' &
432  ,north_bc,south_bc &
433  ,east_bc,west_bc &
434  ,atm%regional_bc_bounds%is_east &
435  ,atm%regional_bc_bounds%ie_east &
436  ,atm%regional_bc_bounds%js_east &
437  ,atm%regional_bc_bounds%je_east &
438  ,atm%regional_bc_bounds%is_east_uvs &
439  ,atm%regional_bc_bounds%ie_east_uvs &
440  ,atm%regional_bc_bounds%js_east_uvs &
441  ,atm%regional_bc_bounds%je_east_uvs &
442  ,atm%regional_bc_bounds%is_east_uvw &
443  ,atm%regional_bc_bounds%ie_east_uvw &
444  ,atm%regional_bc_bounds%js_east_uvw &
445  ,atm%regional_bc_bounds%je_east_uvw &
446  ,klev_out &
447  ,ntracers &
448  ,bc_t0%east )
449 !
450  bc_east_t0=>bc_t0%east
451  bc_east_t1=>bc_t1%east
452 !
453  endif
454 !
455  if(west_bc)then
456  call allocate_regional_bc_arrays('west ' &
457  ,north_bc,south_bc &
458  ,east_bc,west_bc &
459  ,atm%regional_bc_bounds%is_west &
460  ,atm%regional_bc_bounds%ie_west &
461  ,atm%regional_bc_bounds%js_west &
462  ,atm%regional_bc_bounds%je_west &
463  ,atm%regional_bc_bounds%is_west_uvs &
464  ,atm%regional_bc_bounds%ie_west_uvs &
465  ,atm%regional_bc_bounds%js_west_uvs &
466  ,atm%regional_bc_bounds%je_west_uvs &
467  ,atm%regional_bc_bounds%is_west_uvw &
468  ,atm%regional_bc_bounds%ie_west_uvw &
469  ,atm%regional_bc_bounds%js_west_uvw &
470  ,atm%regional_bc_bounds%je_west_uvw &
471  ,klev_out &
472  ,ntracers &
473  ,bc_t1%west )
474 !
475  call allocate_regional_bc_arrays('west ' &
476  ,north_bc,south_bc &
477  ,east_bc,west_bc &
478  ,atm%regional_bc_bounds%is_west &
479  ,atm%regional_bc_bounds%ie_west &
480  ,atm%regional_bc_bounds%js_west &
481  ,atm%regional_bc_bounds%je_west &
482  ,atm%regional_bc_bounds%is_west_uvs &
483  ,atm%regional_bc_bounds%ie_west_uvs &
484  ,atm%regional_bc_bounds%js_west_uvs &
485  ,atm%regional_bc_bounds%je_west_uvs &
486  ,atm%regional_bc_bounds%is_west_uvw &
487  ,atm%regional_bc_bounds%ie_west_uvw &
488  ,atm%regional_bc_bounds%js_west_uvw &
489  ,atm%regional_bc_bounds%je_west_uvw &
490  ,klev_out &
491  ,ntracers &
492  ,bc_t0%west )
493 !
494  bc_west_t0=>bc_t0%west
495  bc_west_t1=>bc_t1%west
496 !
497  endif
498 !
499 !-----------------------------------------------------------------------
500 !*** We need regional versions of the arrays for surface elevation,
501 !*** latitude/longitude of grid cell corners, and lat/lon of the
502 !*** cell centers because those variables are needed an extra row
503 !*** beyond FV3's normal bounday region width of nhalo_model rows.
504 !-----------------------------------------------------------------------
505 !
506  allocate(phis_reg(isd-1:ied+1,jsd-1:jed+1)) ; phis_reg=real_snan
507 !
508  allocate(agrid_reg(isd-1:ied+1,jsd-1:jed+1,2)); agrid_reg=dbl_snan
509  allocate(grid_reg(isd-1:ied+2,jsd-1:jed+2,2)) ; grid_reg=dbl_snan
510 !
511 !-----------------------------------------------------------------------
512 !*** From the data holding nhalo_model rows of boundary values
513 !*** read in the lat/lon of the grid cell corners and fill in
514 !*** the values of the grid cell centers. The regional mode needs
515 !*** the extra row of data.
516 !-----------------------------------------------------------------------
517 !
519 !
520 !-----------------------------------------------------------------------
521 !*** From the data holding nhalo_model rows of filtered topography
522 !*** read in those values. The regional mode needs the extra row
523 !*** of data.
524 !-----------------------------------------------------------------------
525 !
527 !
528 !-----------------------------------------------------------------------
529 !*** In the init step Atm%phis is given values only in the integration
530 !*** domain but in a regional run values are also needed in the
531 !*** boundary rows. Since the same data is read in the preceding
532 !*** subroutine call as when Atm%phis was first filled, fill its
533 !*** boundary rows now.
534 !-----------------------------------------------------------------------
535 !
536  if(north_bc)then
537  i_start=isd
538  i_end =ied
539  j_start=jsd
540  if(.not.atm%flagstruct%warm_start)then
541  j_end=jsd+nhalo_model-1
542  else
543  j_end=jsd+nhalo_model+1
544  endif
545  do j=j_start,j_end
546  do i=i_start,i_end
547  atm%phis(i,j)=phis_reg(i,j)
548  enddo
549  enddo
550  endif
551 !
552  if(south_bc)then
553  i_start=isd
554  i_end =ied
555  j_end =jed
556  if(.not.atm%flagstruct%warm_start)then
557  j_start=jed-nhalo_model+1
558  else
559  j_start=jed-nhalo_model-1
560  endif
561  do j=j_start,j_end
562  do i=i_start,i_end
563  atm%phis(i,j)=phis_reg(i,j)
564  enddo
565  enddo
566  endif
567  if(east_bc)then
568  i_start=isd
569  j_start=jsd
570  j_end =jed
571  if(.not.atm%flagstruct%warm_start)then
572  i_end=isd+nhalo_model-1
573  else
574  i_end=isd+nhalo_model+1
575  endif
576  do j=j_start,j_end
577  do i=i_start,i_end
578  atm%phis(i,j)=phis_reg(i,j)
579  enddo
580  enddo
581  endif
582  if(west_bc)then
583  i_end =ied
584  j_start=jsd
585  j_end =jed
586  if(.not.atm%flagstruct%warm_start)then
587  i_start=ied-nhalo_model+1
588  else
589  i_start=ied-nhalo_model-1
590  endif
591  do j=j_start,j_end
592  do i=i_start,i_end
593  atm%phis(i,j)=phis_reg(i,j)
594  enddo
595  enddo
596  endif
597 !
598 !-----------------------------------------------------------------------
599 !*** When nudging of specific humidity is selected then we need a
600 !*** reference pressure profile. Compute it now.
601 !-----------------------------------------------------------------------
602 !
603  allocate(pref(npz+1))
604  allocate(dum1d(npz+1))
605 !
606  ps1=101325.
607  pref(npz+1)=ps1
608  call get_eta_level(npz,ps1,pref(1),dum1d,atm%ak,atm%bk )
609 !
610 !-----------------------------------------------------------------------
611 
612  contains
613 
614 !-----------------------------------------------------------------------
615 !
616  subroutine compute_regional_bc_indices(regional_bc_bounds)
617 !
618 !-----------------------------------------------------------------------
619 !*** This routine computes the starting and ending indices for
620 !*** working arrays of task subdomains that lie on the edges
621 !*** of the FV3 regional domain. These arrays will hold boundary
622 !*** region values of scalar variables located at the grid cell
623 !*** centers and wind components lying on the east/west sides
624 !*** and north/south sides of each cell. Note that the width
625 !*** of the domain's boundary region (4 rows) is currently
626 !*** greater than the fundamental width of the task subdomains'
627 !*** halo regions (3 rows). The variables isd,ied,jsd,jed are
628 !*** the task subdomain index limits including their halos.
629 !*** The diagram in subroutine regional_bc_data will help to
630 !*** understand these index limits have the values they do.
631 !-----------------------------------------------------------------------
632 !
633 !------------------------
634 !*** Argument variables
635 !------------------------
636 !
637  type(fv_regional_bc_bounds_type),intent(out) :: regional_bc_bounds
638 !
639 !---------------------
640 !*** Local variables
641 !---------------------
642 !
643  integer, parameter :: invalid_index = -99
644  integer :: halo_diff
645 !
646 !-----------------------------------------------------------------------
647 !***********************************************************************
648 !-----------------------------------------------------------------------
649 !
650  regional_bc_bounds%is_north = invalid_index
651  regional_bc_bounds%ie_north = invalid_index
652  regional_bc_bounds%js_north = invalid_index
653  regional_bc_bounds%je_north = invalid_index
654  regional_bc_bounds%is_north_uvs = invalid_index
655  regional_bc_bounds%ie_north_uvs = invalid_index
656  regional_bc_bounds%js_north_uvs = invalid_index
657  regional_bc_bounds%je_north_uvs = invalid_index
658  regional_bc_bounds%is_north_uvw = invalid_index
659  regional_bc_bounds%ie_north_uvw = invalid_index
660  regional_bc_bounds%js_north_uvw = invalid_index
661  regional_bc_bounds%je_north_uvw = invalid_index
662 
663  regional_bc_bounds%is_south = invalid_index
664  regional_bc_bounds%ie_south = invalid_index
665  regional_bc_bounds%js_south = invalid_index
666  regional_bc_bounds%je_south = invalid_index
667  regional_bc_bounds%is_south_uvs = invalid_index
668  regional_bc_bounds%ie_south_uvs = invalid_index
669  regional_bc_bounds%js_south_uvs = invalid_index
670  regional_bc_bounds%je_south_uvs = invalid_index
671  regional_bc_bounds%is_south_uvw = invalid_index
672  regional_bc_bounds%ie_south_uvw = invalid_index
673  regional_bc_bounds%js_south_uvw = invalid_index
674  regional_bc_bounds%je_south_uvw = invalid_index
675 
676  regional_bc_bounds%is_east = invalid_index
677  regional_bc_bounds%ie_east = invalid_index
678  regional_bc_bounds%js_east = invalid_index
679  regional_bc_bounds%je_east = invalid_index
680  regional_bc_bounds%is_east_uvs = invalid_index
681  regional_bc_bounds%ie_east_uvs = invalid_index
682  regional_bc_bounds%js_east_uvs = invalid_index
683  regional_bc_bounds%je_east_uvs = invalid_index
684  regional_bc_bounds%is_east_uvw = invalid_index
685  regional_bc_bounds%ie_east_uvw = invalid_index
686  regional_bc_bounds%js_east_uvw = invalid_index
687  regional_bc_bounds%je_east_uvw = invalid_index
688 
689  regional_bc_bounds%is_west = invalid_index
690  regional_bc_bounds%ie_west = invalid_index
691  regional_bc_bounds%js_west = invalid_index
692  regional_bc_bounds%je_west = invalid_index
693  regional_bc_bounds%is_west_uvs = invalid_index
694  regional_bc_bounds%ie_west_uvs = invalid_index
695  regional_bc_bounds%js_west_uvs = invalid_index
696  regional_bc_bounds%je_west_uvs = invalid_index
697  regional_bc_bounds%is_west_uvw = invalid_index
698  regional_bc_bounds%ie_west_uvw = invalid_index
699  regional_bc_bounds%js_west_uvw = invalid_index
700  regional_bc_bounds%je_west_uvw = invalid_index
701 !
702 !-----------------------------------------------------------------------
703 !*** Scalar BC indices
704 !-----------------------------------------------------------------------
705 !*** These must reach one row beyond nhalo_model since we must
706 !*** surround the wind points on the cell edges with mass points.
707 !-----------------------------------------------------------------------
708 !
709  halo_diff=nhalo_data-nhalo_model
710 !
711 !-----------
712 !*** North
713 !-----------
714 !
715  if (north_bc) then
716  regional_bc_bounds%is_north=isd-1
717  regional_bc_bounds%ie_north=ied+1
718 !
719  regional_bc_bounds%js_north=jsd-1
720  regional_bc_bounds%je_north=0
721  endif
722 !
723 !-----------
724 !*** South
725 !-----------
726 !
727  if (south_bc) then
728  regional_bc_bounds%is_south=isd-1
729  regional_bc_bounds%ie_south=ied+1
730 !
731  regional_bc_bounds%js_south=jed-nhalo_model+1
732  regional_bc_bounds%je_south=jed+1
733  endif
734 !
735 !----------
736 !*** East
737 !----------
738 !
739  if (east_bc) then
740  regional_bc_bounds%is_east=isd-1
741  regional_bc_bounds%ie_east=0
742 !
743  regional_bc_bounds%js_east=jsd-1
744  if(north_bc)then
745  regional_bc_bounds%js_east=1
746  endif
747 !
748  regional_bc_bounds%je_east=jed+1
749  if(south_bc)then
750  regional_bc_bounds%je_east=jed-nhalo_model
751  endif
752  endif
753 !
754 !----------
755 !*** West
756 !----------
757 !
758  if (west_bc) then
759  regional_bc_bounds%is_west=ied-nhalo_model+1
760  regional_bc_bounds%ie_west=ied+1
761 !
762  regional_bc_bounds%js_west=jsd-1
763  if(north_bc)then
764  regional_bc_bounds%js_west=1
765  endif
766 !
767  regional_bc_bounds%je_west=jed+1
768  if(south_bc)then
769  regional_bc_bounds%je_west=jed-nhalo_model
770  endif
771  endif
772 !
773 !-----------------------------------------------------------------------
774 !*** Wind component BC indices
775 !-----------------------------------------------------------------------
776 !
777 !-----------
778 !*** North
779 !-----------
780 !
781  if (north_bc) then
782  regional_bc_bounds%is_north_uvs=isd
783  regional_bc_bounds%ie_north_uvs=ied
784 !
785  regional_bc_bounds%js_north_uvs=jsd
786 !xxxxxx regional_bc_bounds%je_north_uvs=0
787 !xxxxxx regional_bc_bounds%je_north_uvs=1
788  regional_bc_bounds%je_north_uvs=1
789 !
790  regional_bc_bounds%is_north_uvw=isd
791  regional_bc_bounds%ie_north_uvw=ied+1
792 !
793  regional_bc_bounds%js_north_uvw=jsd
794  regional_bc_bounds%je_north_uvw=0
795  endif
796 !
797 !-----------
798 !*** South
799 !-----------
800 !
801  if (south_bc) then
802  regional_bc_bounds%is_south_uvs=isd
803  regional_bc_bounds%ie_south_uvs=ied
804 !
805 !xxxxxregional_bc_bounds%js_south_uvs=jed-nhalo_model+2
806  regional_bc_bounds%js_south_uvs=jed-nhalo_model+1
807  regional_bc_bounds%je_south_uvs=jed+1
808 !
809  regional_bc_bounds%is_south_uvw=isd
810  regional_bc_bounds%ie_south_uvw=ied+1
811 !
812  regional_bc_bounds%js_south_uvw=jed-nhalo_model+1
813  regional_bc_bounds%je_south_uvw=jed
814  endif
815 !
816 !----------
817 !*** East
818 !----------
819 !
820  if (east_bc) then
821  regional_bc_bounds%is_east_uvs=isd
822  regional_bc_bounds%ie_east_uvs=0
823 !
824  regional_bc_bounds%js_east_uvs=jsd
825  if(north_bc)then
826 !xxxx regional_bc_bounds%js_east_uvs=2 !<-- north side of cell at j=2 (north bdry contains north side of j=1)
827  regional_bc_bounds%js_east_uvs=1
828  endif
829 !
830  regional_bc_bounds%je_east_uvs=jed+1
831  if(south_bc)then
832 !xxxx regional_bc_bounds%je_east_uvs=jed-nhalo_model
833  regional_bc_bounds%je_east_uvs=jed-nhalo_model+1
834  endif
835 !
836 ! regional_bc_bounds%is_east_uvw=isd-1
837  regional_bc_bounds%is_east_uvw=isd
838  regional_bc_bounds%ie_east_uvw=0
839 !
840  regional_bc_bounds%js_east_uvw=jsd
841  if(north_bc)then
842  regional_bc_bounds%js_east_uvw=1
843  endif
844  regional_bc_bounds%je_east_uvw=jed
845  if(south_bc)then
846  regional_bc_bounds%je_east_uvw=jed-nhalo_model
847  endif
848  endif
849 !
850 !----------
851 !*** West
852 !----------
853 !
854  if (west_bc) then
855  regional_bc_bounds%is_west_uvs=ied-nhalo_model+1
856  regional_bc_bounds%ie_west_uvs=ied
857 !
858  regional_bc_bounds%js_west_uvs=jsd
859  if(north_bc)then
860 !xxxx regional_bc_bounds%js_west_uvs=2
861  regional_bc_bounds%js_west_uvs=1
862  endif
863 !
864  regional_bc_bounds%je_west_uvs=jed+1
865  if(south_bc)then
866 !xxxx regional_bc_bounds%je_west_uvs=jed-nhalo_model
867  regional_bc_bounds%je_west_uvs=jed-nhalo_model+1
868  endif
869 !
870  regional_bc_bounds%is_west_uvw=ied-nhalo_model+2
871  regional_bc_bounds%ie_west_uvw=ied+1
872 !
873  regional_bc_bounds%js_west_uvw=jsd
874  if(north_bc)then
875  regional_bc_bounds%js_west_uvw=1
876  endif
877 !
878  regional_bc_bounds%je_west_uvw=jed
879  if(south_bc)then
880  regional_bc_bounds%je_west_uvw=jed-nhalo_model
881  endif
882  endif
883 !
884 !-----------------------------------------------------------------------
885 !
886  end subroutine compute_regional_bc_indices
887 !
888 !-----------------------------------------------------------------------
889 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
890 !-----------------------------------------------------------------------
891 !
892  subroutine read_regional_lon_lat
893 !
894 !-----------------------------------------------------------------------
895 !*** Read the longitude/latitude of the grid cell corners from
896 !*** the external file holding the additional row of data required
897 !*** by the regional domain.
898 !-----------------------------------------------------------------------
899  use netcdf
900 !-----------------------------------------------------------------------
901  implicit none
902 !-----------------------------------------------------------------------
903 !
904 !--------------------
905 !*** Local variables
906 !--------------------
907 !
908  integer :: i_start_data,istat,j_start_data,n,ncid_grid,var_id
909 !
910  character(len=150) :: filename,vname
911 !
912 !-----------------------------------------------------------------------
913 !***********************************************************************
914 !-----------------------------------------------------------------------
915 !
916 !-----------------------------------------------------------------------
917 !*** Open the data file.
918 !-----------------------------------------------------------------------
919 !
920  filename='INPUT/'//trim(grid_data)
921 !
922  call check(nf90_open(filename,nf90_nowrite,ncid_grid))
923 !
924 ! write(0,*)' opened grid file',trim(filename)
925 !-----------------------------------------------------------------------
926 !*** The longitude and latitude are on the super grid. We need only
927 !*** the points on each corner of the grid cells which is every other
928 !*** point on the super grid.
929 !-----------------------------------------------------------------------
930 !
931  i_start_data=2*(isd+nhalo_model)-1
932  j_start_data=2*(jsd+nhalo_model)-1
933 !
934 ! write(0,11110)i_start_data,j_start_data
935 11110 format(' i_start_data=',i5,' j_start_data=',i5)
936 !---------------
937 !*** Longitude
938 !---------------
939 !
940  vname='x'
941  call check(nf90_inq_varid(ncid_grid,vname,var_id))
942  call check(nf90_get_var(ncid_grid,var_id &
943  ,grid_reg(isd-1:ied+2,jsd-1:jed+2,1) &
944  ,start=(/i_start_data,j_start_data/) &
945  ,stride=(/2,2/) ) )
946 !
947 !--------------
948 !*** Latitude
949 !--------------
950 !
951  vname='y'
952  call check(nf90_inq_varid(ncid_grid,vname,var_id))
953  call check(nf90_get_var(ncid_grid,var_id &
954  ,grid_reg(isd-1:ied+2,jsd-1:jed+2,2) &
955  ,start=(/i_start_data,j_start_data/) &
956  ,stride=(/2,2/) ) )
957 !
958  call check(nf90_close(ncid_grid))
959 !
960 !-----------------------------------------------------------------------
961 !*** Convert from degrees to radians.
962 !-----------------------------------------------------------------------
963 !
964  do n=1,2
965  do j=jsd-1,jed+2
966  do i=isd-1,ied+2
967  grid_reg(i,j,n)=grid_reg(i,j,n)*pi/180.
968  enddo
969  enddo
970  enddo
971 !
972 !-----------------------------------------------------------------------
973 !*** Compute the longitude/latitude in the cell centers.
974 !-----------------------------------------------------------------------
975 !
976  do j=jsd-1,jed+1
977  do i=isd-1,ied+1
978  call cell_center2(grid_reg(i,j, 1:2), grid_reg(i+1,j, 1:2), &
979  grid_reg(i,j+1,1:2), grid_reg(i+1,j+1,1:2), &
980  agrid_reg(i,j,1:2) )
981  enddo
982  enddo
983 !
984 !-----------------------------------------------------------------------
985 !
986  end subroutine read_regional_lon_lat
987 !
988 !-----------------------------------------------------------------------
989 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
990 !-----------------------------------------------------------------------
991 !
992  subroutine read_regional_filtered_topo
993 !
994 !-----------------------------------------------------------------------
995 !*** Read the filtered topography including the extra outer row.
996 !-----------------------------------------------------------------------
997 !-----------------------------------------------------------------------
998  implicit none
999 !-----------------------------------------------------------------------
1000 !
1001 !---------------------
1002 !*** Local variables
1003 !---------------------
1004 !
1005  integer :: i,i_start_data,istat,j,j_start_data,ncid_oro,var_id
1006 !
1007  character(len=150) :: filename,vname
1008 !
1009 !-----------------------------------------------------------------------
1010 !***********************************************************************
1011 !-----------------------------------------------------------------------
1012 !
1013 !-----------------------------------------------------------------------
1014 !*** Get the name of the working directory. Open the data file.
1015 !-----------------------------------------------------------------------
1016 !
1017  filename='INPUT/'//trim(oro_data)
1018 
1019  if (is_master()) then
1020  write(*,23421)trim(filename)
1021 23421 format(' topo filename=',a)
1022  endif
1023 !
1024  call check(nf90_open(filename,nf90_nowrite,ncid_oro))
1025 !
1026 !-----------------------------------------------------------------------
1027 !*** Read in the data including the extra outer row.
1028 !-----------------------------------------------------------------------
1029 !
1030  i_start_data=isd+nhalo_model
1031  j_start_data=jsd+nhalo_model
1032 !
1033  vname='orog_filt'
1034  call check(nf90_inq_varid(ncid_oro,vname,var_id))
1035  call check(nf90_get_var(ncid_oro,var_id &
1036  ,phis_reg(isd-1:ied+1,jsd-1:jed+1) &
1037  ,start=(/i_start_data,j_start_data/)))
1038 !
1039  call check(nf90_close(ncid_oro))
1040 !
1041 !-----------------------------------------------------------------------
1042 !*** We want the geopotential.
1043 !-----------------------------------------------------------------------
1044 !
1045  do j=jsd-1,jed+1
1046  do i=isd-1,ied+1
1047  phis_reg(i,j)=phis_reg(i,j)*grav
1048  enddo
1049  enddo
1050 !
1051 !-----------------------------------------------------------------------
1052 !***********************************************************************
1053 !-----------------------------------------------------------------------
1054 !
1055  end subroutine read_regional_filtered_topo
1056 !
1057 !-----------------------------------------------------------------------
1058 !
1059  end subroutine setup_regional_bc
1060 !
1061 !-----------------------------------------------------------------------
1062 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1063 !-----------------------------------------------------------------------
1064 !
1065  subroutine start_regional_cold_start(Atm, ak, bk, levp &
1066  ,is ,ie ,js ,je &
1067  ,isd,ied,jsd,jed )
1069 !-----------------------------------------------------------------------
1070 !*** Prepare the regional run for a cold start.
1071 !-----------------------------------------------------------------------
1072  implicit none
1073 !-----------------------------------------------------------------------
1074 !
1075 !------------------------
1076 !*** Argument variables
1077 !------------------------
1078 !
1079  type(fv_atmos_type),intent(inout) :: atm
1080 !
1081  integer ,intent(in) :: is ,ie ,js ,je & !<-- Integration limits of task subdomain
1082  ,isd,ied,jsd,jed & !<-- Memory limits of task subdomain
1083  ,levp
1084 !
1085  real,intent(in) :: ak(1:levp+1), bk(1:levp+1)
1086 !
1087 !---------------------
1088 !*** Local variables
1089 !---------------------
1090 !
1091  integer :: k
1092 !
1093 !-----------------------------------------------------------------------
1094 !***********************************************************************
1095 !-----------------------------------------------------------------------
1096 !
1097 ! get the source of the input data
1098 !
1099  call get_data_source(data_source,atm%flagstruct%regional)
1100 !
1101  call setup_regional_bc(atm &
1102  ,isd, ied, jsd, jed &
1103  ,atm%npx, atm%npy )
1104 !
1105  bc_hour=0
1106  call regional_bc_data(atm, bc_hour &
1107  ,is, ie, js, je &
1108  ,isd, ied, jsd, jed &
1109  ,ak, bk )
1110  call regional_bc_t1_to_t0(bc_t1, bc_t0 & !
1111  ,atm%npz &
1112  ,ntracers &
1113  ,atm%regional_bc_bounds ) !
1114 !
1116 !
1117  call regional_bc_data(atm, bc_hour &
1118  ,is, ie, js, je & ! from the 2nd time level
1119  ,isd, ied, jsd, jed & ! in the BC file.
1120  ,ak, bk ) !
1121 !
1122  allocate (ak_in(1:levp+1))
1123  allocate (bk_in(1:levp+1)) ! remapping BC updates during the forecast.
1124  do k=1,levp+1
1125  ak_in(k)=ak(k)
1126  bk_in(k)=bk(k)
1127  enddo
1128 !
1129 !-----------------------------------------------------------------------
1130 !
1131  end subroutine start_regional_cold_start
1132 !
1133 !-----------------------------------------------------------------------
1134 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1135 !-----------------------------------------------------------------------
1136 !
1137  subroutine start_regional_restart(Atm &
1138  ,isc,iec,jsc,jec &
1139  ,isd,ied,jsd,jed )
1141 !-----------------------------------------------------------------------
1142 !*** Prepare the regional forecast for a restart.
1143 !-----------------------------------------------------------------------
1144  implicit none
1145 !-----------------------------------------------------------------------
1146 !
1147 !------------------------
1148 !*** Argument variables
1149 !------------------------
1150 !
1151  type(fv_atmos_type),intent(inout) :: atm
1152 !
1153  integer ,intent(in) :: isc,iec,jsc,jec & !<-- Integration limits of task subdomain
1154  ,isd,ied,jsd,jed
1155 !
1156 !---------------------
1157 !*** Local variables
1158 !---------------------
1159 !
1160  integer :: ierr, ios
1161  real, allocatable :: wk2(:,:)
1162 !
1163  logical :: filtered_terrain
1164  logical :: gfs_dwinds
1165  integer :: levp
1166  logical :: checker_tr
1167  integer :: nt_checker
1168  namelist /external_ic_nml/ filtered_terrain, levp, gfs_dwinds &
1169  ,checker_tr, nt_checker
1170 !-----------------------------------------------------------------------
1171 !***********************************************************************
1172 !-----------------------------------------------------------------------
1173 !
1174 !-----------------------------------------------------------------------
1175 !*** Read the number of model layers in the external forecast (=levp).
1176 !-----------------------------------------------------------------------
1177 !
1178  read (input_nml_file,external_ic_nml,iostat=ios)
1179  ierr = check_nml_error(ios,'external_ic_nml')
1180  if(ierr/=0)then
1181  write(0,11011)ierr
1182 11011 format(' start_regional_restart failed to read external_ic_nml ierr=',i3)
1183  endif
1184 !
1185 !-----------------------------------------------------------------------
1186 !*** Preliminary setup for the forecast.
1187 !-----------------------------------------------------------------------
1188 !
1189  call setup_regional_bc(atm &
1190  ,isd, ied, jsd, jed &
1191  ,atm%npx, atm%npy )
1192 !
1193  allocate (wk2(levp+1,2))
1194  allocate (ak_in(levp+1))
1195  allocate (bk_in(levp+1)) ! remapping BC updates during the forecast.
1196  call read_data('INPUT/gfs_ctrl.nc','vcoord',wk2, no_domain=.true.)
1197  ak_in(1:levp+1) = wk2(1:levp+1,1)
1198  ak_in(1) = 1.e-9
1199  bk_in(1:levp+1) = wk2(1:levp+1,2)
1200  deallocate(wk2)
1201  bc_hour=nint(current_time_in_seconds/3600.)
1202 !
1203 !-----------------------------------------------------------------------
1204 !*** Fill time level t1 from the BC file at the restart time.
1205 !-----------------------------------------------------------------------
1206 !
1207  call regional_bc_data(atm, bc_hour &
1208  ,isc, iec, jsc, jec &
1209  ,isd, ied, jsd, jed &
1210  ,ak_in, bk_in )
1211 !
1212 !-----------------------------------------------------------------------
1213 !
1214  end subroutine start_regional_restart
1215 !
1216 !-----------------------------------------------------------------------
1217 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1218 !-----------------------------------------------------------------------
1219 !
1220  subroutine read_new_bc_data(Atm, Time, Time_step_atmos, p_split &
1221  ,isd,ied,jsd,jed )
1223 !-----------------------------------------------------------------------
1224 !*** When it is time to read new boundary data from the external files
1225 !*** move time level t1 to t0 and then read the data into t1.
1226 !-----------------------------------------------------------------------
1227  implicit none
1228 !-----------------------------------------------------------------------
1229 !
1230 !------------------------
1231 !*** Argument variables
1232 !------------------------
1233 !
1234  type(fv_atmos_type),intent(inout) :: atm
1235  type(time_type),intent(in) :: time
1236  type(time_type),intent(in) :: time_step_atmos
1237 
1238  integer,intent(in) :: isd,ied,jsd,jed & !<-- Memory limits of task subdomain
1239  ,p_split
1240 !
1241 !---------------------
1242 !*** Local variables
1243 !---------------------
1244 !
1245  integer :: atmos_time_step, sec
1246  real :: dt_atmos
1247  type(time_type) :: atmos_time
1248 !
1249 !-----------------------------------------------------------------------
1250 !***********************************************************************
1251 !-----------------------------------------------------------------------
1252 !
1253  atmos_time = time - atm%Time_init
1254  atmos_time_step = atmos_time / time_step_atmos
1255  current_time_in_seconds = time_type_to_real( atmos_time )
1256  if (mpp_pe() == 0) write(*,"('current_time_seconds = ',f9.1)")current_time_in_seconds
1257 !
1258  call get_time (time_step_atmos, sec)
1259  dt_atmos = real(sec)
1260 !
1261  if(atmos_time_step==0.or.atm%flagstruct%warm_start)then
1262  ntimesteps_per_bc_update=nint(bc_update_interval*3600./(dt_atmos/real(abs(p_split))))
1263  endif
1264 !
1265  if(atmos_time_step+1>=ntimesteps_per_bc_update.and.mod(atmos_time_step,ntimesteps_per_bc_update)==0 &
1266  .or. &
1267  atm%flagstruct%warm_start.and.begin_regional_restart)then
1268 !
1269  begin_regional_restart=.false.
1271 !
1272 !-----------------------------------------------------------------------
1273 !*** Transfer the time level t1 data to t0.
1274 !-----------------------------------------------------------------------
1275 !
1277  ,atm%npz &
1278  ,ntracers &
1279  ,atm%regional_bc_bounds )
1280 !
1281 !-----------------------------------------------------------------------
1282 !*** Fill time level t1 from the BC file containing data from
1283 !*** the next time level.
1284 !-----------------------------------------------------------------------
1285 !
1286  call regional_bc_data(atm, bc_hour &
1287  ,atm%bd%is, atm%bd%ie &
1288  ,atm%bd%js, atm%bd%je &
1289  ,isd, ied, jsd, jed &
1290  ,ak_in, bk_in )
1291  endif
1292 !
1293 !-----------------------------------------------------------------------
1294 !
1295  end subroutine read_new_bc_data
1296 !
1297 !-----------------------------------------------------------------------
1298 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1299 !-----------------------------------------------------------------------
1300 !
1301  subroutine regional_bc_data(Atm,bc_hour &
1302  ,is,ie,js,je &
1303  ,isd,ied,jsd,jed &
1304  ,ak,bk )
1306 !-----------------------------------------------------------------------
1307 !*** Regional boundary data is obtained from the external BC file.
1308 !-----------------------------------------------------------------------
1309 !-----------------------------------------------------------------------
1310  implicit none
1311 !-----------------------------------------------------------------------
1312 !
1313 !------------------------
1314 !*** Argument variables
1315 !------------------------
1316 !
1317 !-----------
1318 !*** Input
1319 !-----------
1320 !
1321  integer,intent(in) :: bc_hour
1322 !
1323  integer,intent(in) :: is,ie,js,je & !<-- Compute limits of task subdomain
1324  ,isd,ied,jsd,jed
1325 !
1326  real,dimension(:),intent(in) :: ak,bk
1327 !
1328 !-----------------
1329 !*** Input/output
1330 !-----------------
1331 !
1332  type(fv_atmos_type),target,intent(inout) :: atm
1333 !
1334 !---------------------
1335 !*** Local variables
1336 !---------------------
1337 !
1338  integer :: dimid,i,j,k,klev_in,klev_out,n,nlev
1339 !
1340  integer :: is_north,is_south,is_east,is_west &
1341  ,ie_north,ie_south,ie_east,ie_west &
1342  ,js_north,js_south,js_east,js_west &
1343  ,je_north,je_south,je_east,je_west
1344 !
1345  integer :: is_u,ie_u,js_u,je_u &
1346  ,is_v,ie_v,js_v,je_v
1347 !
1348  integer :: is_input,ie_input,js_input,je_input
1349 !
1350  integer :: i_start,i_end,j_start,j_end
1351 !
1352  real,dimension(:,:,:),allocatable :: ud,vd,uc,vc
1353 !
1354  real,dimension(:,:),allocatable :: ps_reg
1355  real,dimension(:,:,:),allocatable :: ps_input,t_input &
1356  ,w_input,zh_input
1357  real,dimension(:,:,:),allocatable :: u_s_input,v_s_input &
1358  ,u_w_input,v_w_input
1359  real,dimension(:,:,:,:),allocatable :: tracers_input
1360 !
1361  real(kind=R_GRID), dimension(2):: p1, p2, p3, p4
1362  real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
1363 
1364 #undef USE_FMS_READ
1365 #ifdef USE_FMS_READ
1366  integer :: isc2, iec2, jsc2, jec2
1367  real(kind=R_GRID), allocatable, dimension(:,:) :: tmpx, tmpy
1368  integer :: start(4), nread(4)
1369  real(kind=R_GRID), allocatable, dimension(:,:,:) :: reg_grid
1370  real(kind=R_GRID), allocatable, dimension(:,:,:) :: reg_agrid
1371 #endif
1372 !
1373  logical,save :: computed_regional_bc_indices=.false.
1374 !
1375  character(len=3) :: int_to_char
1376  character(len=5) :: side
1377  character(len=6) :: fmt='(i3.3)'
1378 !
1379  character(len=50) :: file_name
1380 !
1381  integer,save :: kount1=0,kount2=0
1382 !
1383  character(len=60) :: var_name_root
1384  integer :: nside,nt,index
1385  logical :: required
1386 !
1387  logical :: call_remap
1388 !
1389 !-----------------------------------------------------------------------
1390 !***********************************************************************
1391 !-----------------------------------------------------------------------
1392 !
1393 !-----------------------------------------------------------------------
1394 !*** Only boundary tasks are needed.
1395 !-----------------------------------------------------------------------
1396 !
1397  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
1398  return
1399  endif
1400 !if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE')
1401 !
1402 !-----------------------------------------------------------------------
1403 !
1404  klev_out=atm%npz
1405 !
1406 !-----------------------------------------------------------------------
1407 !*** Construct the name of the regional BC file to be read.
1408 !-----------------------------------------------------------------------
1409 !
1410  write(int_to_char,fmt) bc_hour
1411  file_name='INPUT/gfs_bndy.tile7.'//int_to_char//'.nc'
1412 !
1413  if (is_master()) then
1414  write(*,22211)trim(file_name)
1415 22211 format(' regional_bc_data file_name=',a)
1416  endif
1417 !-----------------------------------------------------------------------
1418 !*** Open the regional BC file.
1419 !*** Find the # of layers (klev_in) in the BC input.
1420 !-----------------------------------------------------------------------
1421 !
1422  call check(nf90_open(file_name,nf90_nowrite,ncid))
1423 !
1424  call check(nf90_inq_dimid(ncid,'lev',dimid))
1425  call check(nf90_inquire_dimension(ncid,dimid,len=klev_in))
1426 !
1427 !-----------------------------------------------------------------------
1428 !*** Allocate the boundary variables and initialize them to garbage.
1429 !-----------------------------------------------------------------------
1430 !
1431  is_input=is-nhalo_data
1432  ie_input=ie+nhalo_data
1433  js_input=js-nhalo_data
1434  je_input=je+nhalo_data
1435 !
1436  allocate( ps_input(is_input:ie_input,js_input:je_input,1)) ; ps_input=real_snan
1437  allocate( t_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; t_input=real_snan
1438  allocate( w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; w_input=real_snan
1439  allocate( zh_input(is_input:ie_input,js_input:je_input,1:klev_in+1)) ; zh_input=real_snan
1440  allocate(u_s_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; u_s_input=real_snan
1441  allocate(v_s_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; v_s_input=real_snan
1442  allocate(u_w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; u_w_input=real_snan
1443  allocate(v_w_input(is_input:ie_input,js_input:je_input,1:klev_in)) ; v_w_input=real_snan
1444 !
1445  allocate(tracers_input(is_input:ie_input,js_input:je_input,klev_in,ntracers)) ; tracers_input=real_snan
1446 !
1447 !-----------------------------------------------------------------------
1448 !*** Extract each variable from the regional BC file. The final
1449 !*** argument is the object being filled.
1450 !-----------------------------------------------------------------------
1451 !
1452 !------------------
1453 !*** Sfc pressure
1454 !------------------
1455 !
1456  nlev=1
1457  var_name_root='ps'
1458  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1459  ,nlev &
1460  ,ntracers &
1461 ! ,Atm%regional_bc_bounds &
1462  ,var_name_root &
1463  ,array_3d=ps_input )
1464 !
1465 !-----------------------
1466 !*** Vertical velocity
1467 !-----------------------
1468 !
1469  nlev=klev_in
1470  var_name_root='w'
1471  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1472  ,nlev &
1473  ,ntracers &
1474 ! ,Atm%regional_bc_bounds &
1475  ,var_name_root &
1476  ,array_3d=w_input)
1477 !
1478 !-----------------------
1479 !*** Interface heights
1480 !-----------------------
1481 !
1482  nlev=klev_in+1
1483  var_name_root='zh'
1484  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1485  ,nlev &
1486  ,ntracers &
1487 ! ,Atm%regional_bc_bounds &
1488  ,var_name_root &
1489  ,array_3d=zh_input)
1490 !
1491 !--------------------------
1492 !*** Sensible temperature
1493 !--------------------------
1494 !
1495  if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
1496  nlev=klev_in
1497  var_name_root='t'
1498  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1499  ,nlev &
1500  ,ntracers &
1501 ! ,Atm%regional_bc_bounds &
1502  ,var_name_root &
1503  ,array_3d=t_input)
1504  endif
1505 !
1506 !-----------------------------
1507 !*** U component south/north
1508 !-----------------------------
1509 !
1510  nlev=klev_in
1511  var_name_root='u_s'
1512  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1513  ,nlev &
1514  ,ntracers &
1515 ! ,Atm%regional_bc_bounds &
1516  ,var_name_root &
1517  ,array_3d=u_s_input)
1518 !
1519 !-----------------------------
1520 !*** V component south/north
1521 !-----------------------------
1522 !
1523  nlev=klev_in
1524  var_name_root='v_s'
1525  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1526  ,nlev &
1527  ,ntracers &
1528 ! ,Atm%regional_bc_bounds &
1529  ,var_name_root &
1530  ,array_3d=v_s_input)
1531 !
1532 !---------------------------
1533 !*** U component east/west
1534 !---------------------------
1535 !
1536  nlev=klev_in
1537  var_name_root='u_w'
1538  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1539  ,nlev &
1540  ,ntracers &
1541 ! ,Atm%regional_bc_bounds &
1542  ,var_name_root &
1543  ,array_3d=u_w_input)
1544 !
1545 !---------------------------
1546 !*** V component east/west
1547 !---------------------------
1548 !
1549  nlev=klev_in
1550  var_name_root='v_w'
1551  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1552  ,nlev &
1553  ,ntracers &
1554 ! ,Atm%regional_bc_bounds &
1555  ,var_name_root &
1556  ,array_3d=v_w_input)
1557 !----------------------
1558 !*** tracers
1559 !-----------------------
1560 
1561  nlev=klev_in
1562 !
1563 !-----------------------------------------------------------------------
1564 !*** Read the tracers specified in the field_table. If they are not
1565 !*** in the input data then print a warning and set them to 0 in the
1566 !*** boundary. Some tracers are mandatory to have, because they are
1567 !*** used later for calculating virtual potential temperature etc.
1568 !-----------------------------------------------------------------------
1569 !
1570  do nt = 1, ntracers
1571  call get_tracer_names(model_atmos, nt, var_name_root)
1572  index= get_tracer_index(model_atmos,trim(var_name_root))
1573  if (index==liq_water_index .or. index==sphum_index) then
1574  required = .true.
1575  else
1576  required = .false.
1577  endif
1578  call read_regional_bc_file(is_input,ie_input,js_input,je_input &
1579  ,nlev &
1580  ,ntracers &
1581  ,var_name_root &
1582  ,array_4d=tracers_input &
1583  ,tlev=index &
1584  ,required=required )
1585  enddo
1586 !
1587 !-----------------------------------------------------------------------
1588 !*** We now have the boundary variables from the BC file on the
1589 !*** levels of the input data. Before remapping the 3-D variables
1590 !*** from the input levels to the model integration levels we will
1591 !*** simply copy the 2-D sfc pressure (ps) into the model array.
1592 !-----------------------------------------------------------------------
1593 !
1594 ! do j=jsd,jed
1595 ! do i=isd,ied
1596 ! Atm%ps(i,j)=ps(i,j)
1597 ! enddo
1598 ! enddo
1599 !
1600 ! deallocate(ps%north,ps%south,ps%east,ps%west)
1601 !
1602 !-----------------------------------------------------------------------
1603 !*** One final array needs to be allocated. It is the sfc pressure
1604 !*** in the domain's boundary region that is derived from the input
1605 !*** sfc pressure from the BC files. The derived sfc pressure will
1606 !*** be needed in the vertical remapping of the wind components to
1607 !*** the integration levels.
1608 !-----------------------------------------------------------------------
1609 !
1610  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
1611 !
1612 !-----------------------------------------------------------------------
1613 !*** We have the boundary variables from the BC file on the levels
1614 !*** of the input data. Remap the scalars (tracers, vertical
1615 !*** velocity, ozone) to the FV3 domain levels. Scalar remapping
1616 !*** must be done on all four sides before remapping of the winds
1617 !*** since pressures are needed on each side of wind points and so
1618 !*** for a given wind component those pressures could include values
1619 !*** from two different boundary side regions.
1620 !-----------------------------------------------------------------------
1621 !
1622 !-----------------------------------------------------------------------
1623  sides_scalars: do nside=1,4
1624 !-----------------------------------------------------------------------
1625 !
1626  call_remap=.false.
1627 !
1628  if(nside==1)then
1629  if(north_bc)then
1630  call_remap=.true.
1631  side='north'
1632  bc_side_t1=>bc_t1%north
1633  endif
1634  endif
1635 !
1636  if(nside==2)then
1637  if(south_bc)then
1638  call_remap=.true.
1639  side='south'
1640  bc_side_t1=>bc_t1%south
1641  endif
1642  endif
1643 !
1644  if(nside==3)then
1645  if(east_bc)then
1646  call_remap=.true.
1647  side='east '
1648  bc_side_t1=>bc_t1%east
1649  endif
1650  endif
1651 !
1652  if(nside==4)then
1653  if(west_bc)then
1654  call_remap=.true.
1655  side='west '
1656  bc_side_t1=>bc_t1%west
1657  endif
1658  endif
1659 !
1660  if(call_remap)then
1662  ,side &
1663 
1664  ,isd,ied,jsd,jed &
1665 
1666  ,is_input &
1667  ,ie_input & ! Input array
1668  ,js_input & ! index limits.
1669  ,je_input &
1670 
1671  ,klev_in, klev_out &
1672  ,ntracers &
1673  ,ak, bk &
1674 
1675  ,ps_input &
1676  ,t_input & ! BC vbls
1677  ,tracers_input & ! on input
1678  ,w_input & ! model levels
1679  ,zh_input &
1680 
1681  ,phis_reg &
1682 
1683  ,ps_reg &
1684 
1685  ,bc_side_t1 )
1686 !
1687  call set_delp_and_tracers(bc_side_t1,atm%npz,atm%flagstruct%nwat)
1688 !
1689  endif
1690 !
1691 !-----------------------------------------------------------------------
1692  enddo sides_scalars
1693 !-----------------------------------------------------------------------
1694 !
1695 !-----------------------------------------------------------------------
1696 !*** Now that we have the pressure throughout the boundary region
1697 !*** including a row beyond the boundary winds we are ready to
1698 !*** finalize those winds.
1699 !-----------------------------------------------------------------------
1700 !
1701 !-----------------------------------------------------------------------
1702 !*** Transform the D-grid wind components on the north side of
1703 !*** the regional domain then remap them from the input levels
1704 !*** to the integration levels.
1705 !-----------------------------------------------------------------------
1706 !
1707 #ifdef USE_FMS_READ
1708  isc2 = 2*(isd-1+nhalo_data)-1
1709  iec2 = 2*(ied+2+nhalo_data)-1
1710  jsc2 = 2*(jsd-1+nhalo_data)-1
1711  jec2 = 2*(jed+2+nhalo_data)-1
1712  allocate(tmpx(isc2:iec2, jsc2:jec2)) ; tmpx=dbl_snan
1713  allocate(tmpy(isc2:iec2, jsc2:jec2)) ; tmpy=dbl_snan
1714  start = 1; nread = 1
1715  start(1) = isc2; nread(1) = iec2 - isc2 + 1
1716  start(2) = jsc2; nread(2) = jec2 - jsc2 + 1
1717  call read_data("INPUT/grid.tile7.halo4.nc", 'x', tmpx, start, nread, no_domain=.true.)
1718  call read_data("INPUT/grid.tile7.halo4.nc", 'y', tmpy, start, nread, no_domain=.true.)
1719 
1720  allocate(reg_grid(isd-1:ied+2,jsd-1:jed+2,1:2)) ; reg_grid=dbl_snan
1721  do j = jsd-1, jed+2
1722  do i = isd-1, ied+2
1723  reg_grid(i,j,1) = tmpx(2*(i+nhalo_data)-1, 2*(j+nhalo_data)-1)*pi/180.
1724  reg_grid(i,j,2) = tmpy(2*(i+nhalo_data)-1, 2*(j+nhalo_data)-1)*pi/180.
1725  if ( reg_grid(i,j,1) /= grid_reg(i,j,1) ) then
1726  write(0,*)' reg_grid(i,j,1) /= grid_reg(i,j,1) ',i,j, reg_grid(i,j,1),grid_reg(i,j,1)
1727  endif
1728  enddo
1729  enddo
1730 
1731  allocate(reg_agrid(isd-1:ied+1,jsd-1:jed+1,1:2)) ; reg_agrid=dbl_snan
1732  do j=jsd-1,jed+1
1733  do i=isd-1,ied+1
1734  call cell_center2(reg_grid(i,j, 1:2), reg_grid(i+1,j, 1:2), &
1735  reg_grid(i,j+1,1:2), reg_grid(i+1,j+1,1:2), &
1736  reg_agrid(i,j,1:2) )
1737  enddo
1738  enddo
1739 #endif
1740 !
1741 !-----------------------------------------------------------------------
1742 !*** Loop through the four sides of the domain.
1743 !-----------------------------------------------------------------------
1744 !
1745 !-----------------------------------------------------------------------
1746  sides_winds: do nside=1,4
1747 !-----------------------------------------------------------------------
1748 !
1749  call_remap=.false.
1750 
1751  if(nside==1)then
1752  if(north_bc)then
1753  call_remap=.true.
1754  bc_side_t1=>bc_t1%north
1755 !
1756  is_u=atm%regional_bc_bounds%is_north_uvs
1757  ie_u=atm%regional_bc_bounds%ie_north_uvs
1758  js_u=atm%regional_bc_bounds%js_north_uvs
1759  je_u=atm%regional_bc_bounds%je_north_uvs
1760 !
1761  is_v=atm%regional_bc_bounds%is_north_uvw
1762  ie_v=atm%regional_bc_bounds%ie_north_uvw
1763  js_v=atm%regional_bc_bounds%js_north_uvw
1764  je_v=atm%regional_bc_bounds%je_north_uvw
1765  endif
1766  endif
1767 !
1768  if(nside==2)then
1769  if(south_bc)then
1770  call_remap=.true.
1771  bc_side_t1=>bc_t1%south
1772 !
1773  is_u=atm%regional_bc_bounds%is_south_uvs
1774  ie_u=atm%regional_bc_bounds%ie_south_uvs
1775  js_u=atm%regional_bc_bounds%js_south_uvs
1776  je_u=atm%regional_bc_bounds%je_south_uvs
1777 !
1778  is_v=atm%regional_bc_bounds%is_south_uvw
1779  ie_v=atm%regional_bc_bounds%ie_south_uvw
1780  js_v=atm%regional_bc_bounds%js_south_uvw
1781  je_v=atm%regional_bc_bounds%je_south_uvw
1782  endif
1783  endif
1784 !
1785  if(nside==3)then
1786  if(east_bc)then
1787  call_remap=.true.
1788  bc_side_t1=>bc_t1%east
1789 !
1790  is_u=atm%regional_bc_bounds%is_east_uvs
1791  ie_u=atm%regional_bc_bounds%ie_east_uvs
1792  js_u=atm%regional_bc_bounds%js_east_uvs
1793  je_u=atm%regional_bc_bounds%je_east_uvs
1794 !
1795  is_v=atm%regional_bc_bounds%is_east_uvw
1796  ie_v=atm%regional_bc_bounds%ie_east_uvw
1797  js_v=atm%regional_bc_bounds%js_east_uvw
1798  je_v=atm%regional_bc_bounds%je_east_uvw
1799  endif
1800  endif
1801 !
1802  if(nside==4)then
1803  if(west_bc)then
1804  call_remap=.true.
1805  bc_side_t1=>bc_t1%west
1806 !
1807  is_u=atm%regional_bc_bounds%is_west_uvs
1808  ie_u=atm%regional_bc_bounds%ie_west_uvs
1809  js_u=atm%regional_bc_bounds%js_west_uvs
1810  je_u=atm%regional_bc_bounds%je_west_uvs
1811 !
1812  is_v=atm%regional_bc_bounds%is_west_uvw
1813  ie_v=atm%regional_bc_bounds%ie_west_uvw
1814  js_v=atm%regional_bc_bounds%js_west_uvw
1815  je_v=atm%regional_bc_bounds%je_west_uvw
1816  endif
1817  endif
1818 !
1819  if(call_remap)then
1820 !
1821  allocate(ud(is_u:ie_u,js_u:je_u,1:nlev)) ; ud=real_snan
1822  allocate(vd(is_v:ie_v,js_v:je_v,1:nlev)) ; vd=real_snan
1823  allocate(vc(is_u:ie_u,js_u:je_u,1:nlev)) ; vc=real_snan
1824  allocate(uc(is_v:ie_v,js_v:je_v,1:nlev)) ; uc=real_snan
1825 !
1826  do k=1,nlev
1827  do j=js_u,je_u
1828  do i=is_u,ie_u
1829  p1(:) = grid_reg(i, j,1:2)
1830  p2(:) = grid_reg(i+1,j,1:2)
1831  call mid_pt_sphere(p1, p2, p3)
1832  call get_unit_vect2(p1, p2, e1)
1833  call get_latlon_vector(p3, ex, ey)
1834  ud(i,j,k) = u_s_input(i,j,k)*inner_prod(e1,ex)+v_s_input(i,j,k)*inner_prod(e1,ey)
1835  p4(:) = agrid_reg(i,j,1:2) ! cell centroid
1836  call get_unit_vect2(p3, p4, e2) !C-grid V-wind unit vector
1837  vc(i,j,k) = u_s_input(i,j,k)*inner_prod(e2,ex)+v_s_input(i,j,k)*inner_prod(e2,ey)
1838  enddo
1839  enddo
1840 !
1841  do j=js_v,je_v
1842  do i=is_v,ie_v
1843  p1(:) = grid_reg(i,j ,1:2)
1844  p2(:) = grid_reg(i,j+1,1:2)
1845  call mid_pt_sphere(p1, p2, p3)
1846  call get_unit_vect2(p1, p2, e2)
1847  call get_latlon_vector(p3, ex, ey)
1848  vd(i,j,k) = u_w_input(i,j,k)*inner_prod(e2,ex)+v_w_input(i,j,k)*inner_prod(e2,ey)
1849  p4(:) = agrid_reg(i,j,1:2) ! cell centroid
1850  call get_unit_vect2(p3, p4, e1) !C-grid U-wind unit vector
1851  uc(i,j,k) = u_w_input(i,j,k)*inner_prod(e1,ex)+v_w_input(i,j,k)*inner_prod(e1,ey)
1852  enddo
1853  enddo
1854  enddo
1855 !
1856  call remap_dwinds_regional_bc(atm &
1857 
1858  ,is_input &
1859  ,ie_input & ! Index limits for scalars
1860  ,js_input & ! at center of north BC region grid cells.
1861  ,je_input &
1862 
1863  ,is_u &
1864  ,ie_u & ! Index limits for u component
1865  ,js_u & ! on north edge of BC region grid cells.
1866  ,je_u &
1867 
1868  ,is_v &
1869  ,ie_v & ! Index limits for v component
1870  ,js_v & ! on north edge of BC region grid cells.
1871  ,je_v &
1872 
1873  ,klev_in, klev_out &
1874  ,ak, bk &
1875 
1876  ,ps_reg &
1877  ,ud ,vd &
1878  ,uc ,vc &
1879  ,bc_side_t1 )
1880 
1881 !
1882  deallocate(ud,vd,uc,vc)
1883 !
1884  endif
1885 !
1886 !-----------------------------------------------------------------------
1887  enddo sides_winds
1888 !-----------------------------------------------------------------------
1889 !
1890 !-----------------------------------------------------------------------
1891 !*** Close the boundary file.
1892 !-----------------------------------------------------------------------
1893 !
1894  call check(nf90_close(ncid))
1895 ! write(0,*)' closed BC netcdf file'
1896 !
1897 !-----------------------------------------------------------------------
1898 !*** Deallocate working arrays.
1899 !-----------------------------------------------------------------------
1900 !
1901 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1902  if(allocated(ps_input))then
1903  deallocate(ps_input)
1904  endif
1905  if(allocated(t_input))then
1906  deallocate(t_input)
1907  endif
1908  if(allocated(zh_input))then
1909  deallocate(zh_input)
1910  endif
1911  if(allocated(w_input))then
1912  deallocate(w_input)
1913  endif
1914  if(allocated(tracers_input))then
1915  deallocate(tracers_input)
1916  endif
1917  if(allocated(u_s_input))then
1918  deallocate(u_s_input)
1919  endif
1920  if(allocated(u_w_input))then
1921  deallocate(u_w_input)
1922  endif
1923  if(allocated(v_s_input))then
1924  deallocate(v_s_input)
1925  endif
1926  if(allocated(v_w_input))then
1927  deallocate(v_w_input)
1928  endif
1929 !
1930 !-----------------------------------------------------------------------
1931 !*** Fill the remaining boundary arrays starting with the divergence.
1932 !-----------------------------------------------------------------------
1933 !
1934  call fill_divgd_bc
1935 !
1936 !-----------------------------------------------------------------------
1937 !*** Fill the total condensate in the regional boundary array.
1938 !-----------------------------------------------------------------------
1939 !
1940 #ifdef USE_COND
1941  call fill_q_con_bc
1942 #endif
1943 !
1944 !-----------------------------------------------------------------------
1945 !*** Fill moist kappa in the regional domain boundary array.
1946 !-----------------------------------------------------------------------
1947 !
1948 #ifdef MOIST_CAPPA
1949  call fill_cappa_bc
1950 #endif
1951 !
1952 !-----------------------------------------------------------------------
1953 !*** Convert the boundary region sensible temperature array to
1954 !*** FV3's modified virtual potential temperature.
1955 !-----------------------------------------------------------------------
1956 !
1957  call convert_to_virt_pot_temp(isd,ied,jsd,jed,npz &
1959 !
1960 !-----------------------------------------------------------------------
1961 !*** If nudging of the specific humidity has been selected then
1962 !*** nudge the boundary values in the same way as is done for the
1963 !*** interior.
1964 !-----------------------------------------------------------------------
1965 !
1966  if(atm%flagstruct%nudge_qv)then
1967  call nudge_qv_bc(atm,isd,ied,jsd,jed)
1968  endif
1969 !
1970 !-----------------------------------------------------------------------
1971 
1972  contains
1973 
1974 !-----------------------------------------------------------------------
1975 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1976 !-----------------------------------------------------------------------
1977 !
1978  subroutine fill_divgd_bc
1980 !-----------------------------------------------------------------------
1981 !*** For now fill the boundary divergence with zero.
1982 !-----------------------------------------------------------------------
1983  implicit none
1984 !-----------------------------------------------------------------------
1985 !
1986 !--------------------
1987 !*** Local variables
1988 !--------------------
1989 !
1990  integer :: i,ie0,is0,j,je0,js0,k,nside
1991 !
1992  logical :: call_set
1993 !
1994 !-----------------------------------------------------------------------
1995 !***********************************************************************
1996 !-----------------------------------------------------------------------
1997 !
1998 !-----------------------------------------------------------------------
1999 !*** Loop through the four sides.
2000 !-----------------------------------------------------------------------
2001 !
2002  do nside=1,4
2003 !
2004  call_set=.false.
2005 !
2006  if(nside==1)then
2007  if(north_bc)then
2008  call_set=.true.
2009  bc_side_t1=>bc_t1%north
2010  is0=lbound(bc_t1%north%divgd_BC,1)
2011  ie0=ubound(bc_t1%north%divgd_BC,1)
2012  js0=lbound(bc_t1%north%divgd_BC,2)
2013  je0=ubound(bc_t1%north%divgd_BC,2)
2014  endif
2015  endif
2016 !
2017  if(nside==2)then
2018  if(south_bc)then
2019  call_set=.true.
2020  bc_side_t1=>bc_t1%south
2021  is0=lbound(bc_t1%south%divgd_BC,1)
2022  ie0=ubound(bc_t1%south%divgd_BC,1)
2023  js0=lbound(bc_t1%south%divgd_BC,2)
2024  je0=ubound(bc_t1%south%divgd_BC,2)
2025  endif
2026  endif
2027 !
2028  if(nside==3)then
2029  if(east_bc)then
2030  call_set=.true.
2031  bc_side_t1=>bc_t1%east
2032  is0=lbound(bc_t1%east%divgd_BC,1)
2033  ie0=ubound(bc_t1%east%divgd_BC,1)
2034  js0=lbound(bc_t1%east%divgd_BC,2)
2035  je0=ubound(bc_t1%east%divgd_BC,2)
2036  endif
2037  endif
2038 !
2039  if(nside==4)then
2040  if(west_bc)then
2041  call_set=.true.
2042  bc_side_t1=>bc_t1%west
2043  is0=lbound(bc_t1%west%divgd_BC,1)
2044  ie0=ubound(bc_t1%west%divgd_BC,1)
2045  js0=lbound(bc_t1%west%divgd_BC,2)
2046  je0=ubound(bc_t1%west%divgd_BC,2)
2047  endif
2048  endif
2049 !
2050  if(call_set)then
2051  do k=1,klev_out
2052  do j=js0,je0
2053  do i=is0,ie0
2054  bc_side_t1%divgd_BC(i,j,k)=0.
2055  enddo
2056  enddo
2057  enddo
2058  endif
2059 !
2060  enddo
2061 !
2062 !-----------------------------------------------------------------------
2063 !
2064  end subroutine fill_divgd_bc
2065 !
2066 !-----------------------------------------------------------------------
2067 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2068 !-----------------------------------------------------------------------
2069 !
2070 #ifdef USE_COND
2071  subroutine fill_q_con_bc
2072 !
2073 !-----------------------------------------------------------------------
2074 !*** For now fill the total condensate in the boundary regiona
2075 !*** with only the liquid water content.
2076 !-----------------------------------------------------------------------
2077  implicit none
2078 !-----------------------------------------------------------------------
2079 !
2080 !--------------------
2081 !*** Local variables
2082 !--------------------
2083 !
2084  integer :: i,ie0,is0,j,je0,js0,k,nside
2085 !
2086  logical :: call_set
2087 !
2088 !-----------------------------------------------------------------------
2089 !***********************************************************************
2090 !-----------------------------------------------------------------------
2091 !
2092 !-----------------------------------------------------------------------
2093 !*** Loop through the four sides.
2094 !-----------------------------------------------------------------------
2095 !
2096  do nside=1,4
2097  call_set=.false.
2098 !
2099  if(nside==1)then
2100  if(north_bc)then
2101  call_set=.true.
2102  bc_side_t1=>bc_t1%north
2103  is0=lbound(bc_t1%north%q_con_BC,1)
2104  ie0=ubound(bc_t1%north%q_con_BC,1)
2105  js0=lbound(bc_t1%north%q_con_BC,2)
2106  je0=ubound(bc_t1%north%q_con_BC,2)
2107  endif
2108  endif
2109 !
2110  if(nside==2)then
2111  if(south_bc)then
2112  call_set=.true.
2113  bc_side_t1=>bc_t1%south
2114  is0=lbound(bc_t1%south%q_con_BC,1)
2115  ie0=ubound(bc_t1%south%q_con_BC,1)
2116  js0=lbound(bc_t1%south%q_con_BC,2)
2117  je0=ubound(bc_t1%south%q_con_BC,2)
2118  endif
2119  endif
2120 !
2121  if(nside==3)then
2122  if(east_bc)then
2123  call_set=.true.
2124  bc_side_t1=>bc_t1%east
2125  is0=lbound(bc_t1%east%q_con_BC,1)
2126  ie0=ubound(bc_t1%east%q_con_BC,1)
2127  js0=lbound(bc_t1%east%q_con_BC,2)
2128  je0=ubound(bc_t1%east%q_con_BC,2)
2129  endif
2130  endif
2131 !
2132  if(nside==4)then
2133  if(west_bc)then
2134  call_set=.true.
2135  bc_side_t1=>bc_t1%west
2136  is0=lbound(bc_t1%west%q_con_BC,1)
2137  ie0=ubound(bc_t1%west%q_con_BC,1)
2138  js0=lbound(bc_t1%west%q_con_BC,2)
2139  je0=ubound(bc_t1%west%q_con_BC,2)
2140  endif
2141  endif
2142 !
2143  if(call_set)then
2144  do k=1,klev_out
2145  do j=js0,je0
2146  do i=is0,ie0
2147  bc_side_t1%q_con_BC(i,j,k)=bc_side_t1%q_BC(i,j,k,liq_water_index)
2148  enddo
2149  enddo
2150  enddo
2151  endif
2152 !
2153  enddo
2154 !
2155 !-----------------------------------------------------------------------
2156 !
2157  end subroutine fill_q_con_bc
2158 #endif
2159 !
2160 !-----------------------------------------------------------------------
2161 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2162 !-----------------------------------------------------------------------
2163 !
2164 #ifdef MOIST_CAPPA
2165  subroutine fill_cappa_bc
2166 !
2167 !-----------------------------------------------------------------------
2168 !*** Compute cappa in the regional domain boundary area following
2169 !*** Zhao-Carr microphysics.
2170 !-----------------------------------------------------------------------
2171  implicit none
2172 !-----------------------------------------------------------------------
2173 !
2174 !---------------------
2175 !*** Local variables
2176 !---------------------
2177 !
2178  integer :: i1,i2,j1,j2,nside
2179 !
2180  real,dimension(:,:,:),pointer :: cappa,temp,liq_wat,sphum
2181 !
2182  logical :: call_compute
2183 !
2184 !-----------------------------------------------------------------------
2185 !***********************************************************************
2186 !-----------------------------------------------------------------------
2187 !
2188  do nside=1,4
2189  call_compute=.false.
2190 !
2191  if(nside==1)then
2192  if(north_bc)then
2193  call_compute=.true.
2194  bc_side_t1=>bc_t1%north
2195  endif
2196  endif
2197 !
2198  if(nside==2)then
2199  if(south_bc)then
2200  call_compute=.true.
2201  bc_side_t1=>bc_t1%south
2202  endif
2203  endif
2204 !
2205  if(nside==3)then
2206  if(east_bc)then
2207  call_compute=.true.
2208  bc_side_t1=>bc_t1%east
2209  endif
2210  endif
2211 !
2212  if(nside==4)then
2213  if(west_bc)then
2214  call_compute=.true.
2215  bc_side_t1=>bc_t1%west
2216  endif
2217  endif
2218 !
2219  if(call_compute)then
2220  i1=lbound(bc_side_t1%cappa_BC,1)
2221  i2=ubound(bc_side_t1%cappa_BC,1)
2222  j1=lbound(bc_side_t1%cappa_BC,2)
2223  j2=ubound(bc_side_t1%cappa_BC,2)
2224  cappa =>bc_side_t1%cappa_BC
2225  temp =>bc_side_t1%pt_BC
2226  liq_wat=>bc_side_t1%q_BC(:,:,:,liq_water_index)
2227  sphum =>bc_side_t1%q_BC(:,:,:,sphum_index)
2228  call compute_cappa(i1,i2,j1,j2,cappa,temp,liq_wat,sphum)
2229  endif
2230 !
2231  enddo
2232 !
2233 !-----------------------------------------------------------------------
2234 !
2235  end subroutine fill_cappa_bc
2236 !
2237 !-----------------------------------------------------------------------
2238 !***********************************************************************
2239 !-----------------------------------------------------------------------
2240 !
2241  subroutine compute_cappa(i1,i2,j1,j2,cappa,temp,liq_wat,sphum)
2242 !
2243 !-----------------------------------------------------------------------
2244  implicit none
2245 !-----------------------------------------------------------------------
2246 !
2247 !---------------------
2248 !*** Input variables
2249 !---------------------
2250 !
2251  integer,intent(in) :: i1,i2,j1,j2
2252 !
2253  real,dimension(i1:i2,j1:j2,1:npz) :: cappa,temp,liq_wat,sphum
2254 !
2255 !----------------------
2256 !*** Output variables
2257 !----------------------
2258 !
2259 !---------------------
2260 !*** Local variables
2261 !---------------------
2262 !
2263  integer :: i,ie,is,j,je,js,k
2264 !
2265  real :: cvm,qd,ql,qs,qv
2266 !
2267 !-----------------------------------------------------------------------
2268 !***********************************************************************
2269 !-----------------------------------------------------------------------
2270 !
2271  is=lbound(cappa,1)
2272  ie=ubound(cappa,1)
2273  js=lbound(cappa,2)
2274  je=ubound(cappa,2)
2275 !
2276  do k=1,klev_out
2277  do j=js,je
2278  do i=is,ie
2279  qd=max(0.,liq_wat(i,j,k))
2280  if( temp(i,j,k) > tice )then
2281  qs=0.
2282  elseif( temp(i,j,k) < tice-t_i0 )then
2283  qs=qd
2284  else
2285  qs=qd*(tice-temp(i,j,k))/t_i0
2286  endif
2287  ql=qd-qs
2288  qv=max(0.,sphum(i,j,k))
2289  cvm=(1.-(qv+qd))*cv_air + qv*cv_vap + ql*c_liq + qs*c_ice
2290  !
2291  cappa(i,j,k)=rdgas/(rdgas+cvm/(1.+zvir*sphum(i,j,k)))
2292 !
2293  enddo
2294  enddo
2295  enddo
2296 !
2297 !-----------------------------------------------------------------------
2298 !
2299  end subroutine compute_cappa
2300 #endif
2301 !
2302 !-----------------------------------------------------------------------
2303 !
2304  end subroutine regional_bc_data
2305 
2306 !-----------------------------------------------------------------------
2307 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2308 !-----------------------------------------------------------------------
2309  subroutine read_regional_bc_file(is_input,ie_input &
2310  ,js_input,je_input &
2311  ,nlev &
2312  ,ntracers &
2313  ,var_name_root &
2314  ,array_3d &
2315  ,array_4d &
2316  ,tlev &
2317  ,required )
2318 !-----------------------------------------------------------------------
2319 !*** Read the boundary data from the external file generated by
2320 !*** chgres.
2321 !-----------------------------------------------------------------------
2322 !-----------------------------------------------------------------------
2323  implicit none
2324 !-----------------------------------------------------------------------
2325 !
2326 !------------------------
2327 !*** Argument variables
2328 !------------------------
2329 !
2330 !----------
2331 !*** Input
2332 !----------
2333 !
2334  integer,intent(in) :: is_input,ie_input,js_input,je_input,nlev
2335  integer,intent(in) :: ntracers
2336 !
2337  integer,intent(in),optional :: tlev
2338 !
2339  character(len=*),intent(in) :: var_name_root
2340  logical,intent(in),optional :: required
2341 !
2342 !------------
2343 !*** Output
2344 !------------
2345 !
2346  real,dimension(is_input:ie_input,js_input:je_input,1:nlev),intent(out),optional :: array_3d
2347 !
2348  real,dimension(is_input:ie_input,js_input:je_input,1:nlev,1:ntracers),intent(out),optional :: array_4d
2349 !
2350 !---------------------
2351 !*** Local variables
2352 !---------------------
2353 !
2354  integer :: halo,lat,lev,lon
2355 !
2356  integer :: i_count,i_start_array,i_start_data,i_end_array &
2357  ,j_count,j_start_array,j_start_data,j_end_array
2358 !
2359  integer :: dim_id,nctype,ndims,var_id
2360  integer :: nside,status
2361 !
2362  character(len=5) :: dim_name_x & !<-- Dimension names in
2363  ,dim_name_y ! the BC file
2364 !
2365  character(len=80) :: var_name
2366 !
2367  logical :: call_get_var
2368  logical :: required_local
2369 !
2370 !-----------------------------------------------------------------------
2371 !***********************************************************************
2372 !-----------------------------------------------------------------------
2373 !
2374 !-----------------------------------------------------------------------
2375 !*** Process optional argument required, default value is .true.
2376 !-----------------------------------------------------------------------
2377 !
2378  if(present(required)) then
2379  required_local=required
2380  else
2381  required_local=.true.
2382  endif
2383 !
2384 !-----------------------------------------------------------------------
2385 !*** Loop through the four sides of the domain.
2386 !-----------------------------------------------------------------------
2387 !
2388 !-----------------------------------------------------------------------
2389  sides: do nside=1,4
2390 !-----------------------------------------------------------------------
2391 !
2392  call_get_var=.false.
2393 !
2394 !-----------------------------------------------------------------------
2395 !*** Construct the variable's name in the NetCDF file and set
2396 !*** the start locations and point counts for the data file and
2397 !*** for the BC arrays being filled. The input array begins
2398 !*** receiving data at (i_start_array,j_start_array), etc.
2399 !*** The read of the data for the given input array begins at
2400 !*** (i_start_data,j_start_data) and encompasses i_count by
2401 !*** j_count datapoints in each direction.
2402 !-----------------------------------------------------------------------
2403 !
2404 !-----------
2405 !*** North
2406 !-----------
2407 !
2408  if(nside==1)then
2409  if(north_bc)then
2410  call_get_var=.true.
2411 !
2412  var_name=trim(var_name_root)//"_bottom"
2413 !
2414  i_start_array=is_input
2415  i_end_array =ie_input
2416  j_start_array=js_input
2417  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
2418  j_end_array=js_input+nhalo_data
2419  else
2420  j_end_array=js_input+nhalo_data-1
2421  endif
2422 !
2423  i_start_data=i_start_array+nhalo_data
2424  i_count=i_end_array-i_start_array+1
2425  j_start_data=1
2426  j_count=j_end_array-j_start_array+1
2427  endif
2428  endif
2429 !
2430 !-----------
2431 !*** South
2432 !-----------
2433 !
2434  if(nside==2)then
2435  if(south_bc)then
2436  call_get_var=.true.
2437 !
2438  var_name=trim(var_name_root)//"_top"
2439 !
2440  i_start_array=is_input
2441  i_end_array =ie_input
2442  j_start_array=je_input-nhalo_data+1
2443  j_end_array =je_input
2444 !
2445  i_start_data=i_start_array+nhalo_data
2446  i_count=i_end_array-i_start_array+1
2447  j_start_data=1
2448  j_count=j_end_array-j_start_array+1
2449  endif
2450  endif
2451 !
2452 !----------
2453 !*** East
2454 !----------
2455 !
2456  if(nside==3)then
2457  if(east_bc)then
2458  call_get_var=.true.
2459 !
2460  var_name=trim(var_name_root)//"_left"
2461 !
2462  j_start_array=js_input
2463  j_end_array =je_input
2464 !
2465  i_start_array=is_input
2466 !
2467  if(trim(var_name_root)=='u_w'.or.trim(var_name_root)=='v_w')then
2468  i_end_array=is_input+nhalo_data
2469  else
2470  i_end_array=is_input+nhalo_data-1
2471  endif
2472 !
2473  if(north_bc)then
2474  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
2475  j_start_array=js_input+nhalo_data+1
2476  else
2477  j_start_array=js_input+nhalo_data
2478  endif
2479  endif
2480  if(south_bc)then
2481  j_end_array =je_input-nhalo_data
2482  endif
2483 !
2484  i_start_data=1
2485  i_count=i_end_array-i_start_array+1
2486  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
2487  j_start_data=j_start_array-1
2488  else
2489  j_start_data=j_start_array
2490  endif
2491  j_count=j_end_array-j_start_array+1
2492  endif
2493  endif
2494 !
2495 !----------
2496 !*** West
2497 !----------
2498 !
2499  if(nside==4)then
2500  if(west_bc)then
2501  call_get_var=.true.
2502 !
2503  var_name=trim(var_name_root)//"_right"
2504 !
2505  j_start_array=js_input
2506  j_end_array =je_input
2507 !
2508  i_start_array=ie_input-nhalo_data+1
2509  i_end_array=ie_input
2510 !
2511  if(north_bc)then
2512  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
2513  j_start_array=js_input+nhalo_data+1
2514  else
2515  j_start_array=js_input+nhalo_data
2516  endif
2517  endif
2518 !
2519  if(south_bc)then
2520  j_end_array =je_input-nhalo_data
2521  endif
2522 !
2523  i_start_data=1
2524  i_count=i_end_array-i_start_array+1
2525  if(trim(var_name_root)=='u_s'.or.trim(var_name_root)=='v_s')then
2526  j_start_data=j_start_array-1
2527  else
2528  j_start_data=j_start_array
2529  endif
2530  j_count=j_end_array-j_start_array+1
2531  endif
2532  endif
2533 !
2534 !-----------------------------------------------------------------------
2535 !*** Fill this task's subset of boundary data for this 3-D
2536 !*** or 4-D variable. If the variable is a tracer then
2537 !*** check if it is present in the input data. If it is
2538 !*** not then print a warning and set it to zero.
2539 !-----------------------------------------------------------------------
2540 !
2541  if(call_get_var)then
2542  if (present(array_4d)) then
2543  status=nf90_inq_varid(ncid,trim(var_name),var_id)
2544  if (required_local) then
2545  call check(status)
2546  endif
2547  if (status /= nf90_noerr) then
2548  if (east_bc) write(0,*)' WARNING: Tracer ',trim(var_name),' not in input file'
2549  array_4d(:,:,:,tlev)=0.
2550  else
2551  call check(nf90_get_var(ncid,var_id &
2552  ,array_4d(i_start_array:i_end_array &
2553  ,j_start_array:j_end_array &
2554  ,1:nlev, tlev) &
2555  ,start=(/i_start_data,j_start_data,1,tlev/) &
2556  ,count=(/i_count,j_count,nlev,1/)))
2557  endif
2558 !
2559  else
2560  call check(nf90_inq_varid(ncid,trim(var_name),var_id))
2561  call check(nf90_get_var(ncid,var_id &
2562  ,array_3d(i_start_array:i_end_array &
2563  ,j_start_array:j_end_array &
2564  ,1:nlev) &
2565  ,start=(/i_start_data,j_start_data,1/) &
2566  ,count=(/i_count,j_count,nlev/)))
2567  endif
2568  endif
2569 !
2570  enddo sides
2571 !
2572 !-----------------------------------------------------------------------
2573 !
2574  end subroutine read_regional_bc_file
2575 !
2576 !-----------------------------------------------------------------------
2577 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2578 !-----------------------------------------------------------------------
2579 !
2580  subroutine check(status)
2581  integer,intent(in) :: status
2582 !
2583  if(status /= nf90_noerr) then
2584  write(0,*)' check netcdf status=',status
2585  write(0,10001)trim(nf90_strerror(status))
2586 10001 format(' NetCDF error ',a)
2587  stop "Stopped"
2588  endif
2589  end subroutine check
2590 !
2591 !-----------------------------------------------------------------------
2592 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2593 !-----------------------------------------------------------------------
2594 !
2595  subroutine allocate_regional_bc_arrays(side &
2596  ,north_bc,south_bc &
2597  ,east_bc,west_bc &
2598  ,is_0,ie_0,js_0,je_0 &
2599  ,is_sn,ie_sn,js_sn,je_sn &
2600  ,is_we,ie_we,js_we,je_we &
2601  ,klev &
2602  ,ntracers &
2603  ,BC_side )
2605 !-----------------------------------------------------------------------
2606  implicit none
2607 !-----------------------------------------------------------------------
2608 !
2609 !------------------------
2610 !*** Argument variables
2611 !------------------------
2612 !
2613  integer,intent(in) :: klev,ntracers
2614 !
2615  integer,intent(in) :: is_0,ie_0,js_0,je_0
2616  integer,intent(in) :: is_sn,ie_sn,js_sn,je_sn
2617  integer,intent(in) :: is_we,ie_we,js_we,je_we
2618 !
2619  character(len=5),intent(in) :: side
2620 !
2621  logical,intent(in) :: north_bc,south_bc,east_bc,west_bc
2622 !
2623  type(fv_regional_BC_variables),intent(out) :: BC_side
2624 !
2625 !---------------------------------------------------------------------
2626 !*********************************************************************
2627 !---------------------------------------------------------------------
2628 !
2629  if(allocated(bc_side%delp_BC))then
2630  return
2631  endif
2632 !
2633  allocate(bc_side%delp_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%delp_BC=real_snan
2634  allocate(bc_side%divgd_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%divgd_BC=real_snan
2635 !
2636  allocate(bc_side%q_BC (is_0:ie_0,js_0:je_0,1:klev,1:ntracers)) ; bc_side%q_BC=real_snan
2637 !
2638 #ifndef SW_DYNAMICS
2639  allocate(bc_side%pt_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%pt_BC=real_snan
2640  allocate(bc_side%w_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%w_BC=real_snan
2641  allocate(bc_side%delz_BC (is_0:ie_0,js_0:je_0,klev)) ; bc_side%delz_BC=real_snan
2642 #ifdef USE_COND
2643  allocate(bc_side%q_con_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%q_con_BC=real_snan
2644 #ifdef MOIST_CAPPA
2645  allocate(bc_side%cappa_BC(is_0:ie_0,js_0:je_0,klev)) ; bc_side%cappa_BC=real_snan
2646 #endif
2647 #endif
2648 #endif
2649 !
2650 !--------------------
2651 !*** Wind components
2652 !--------------------
2653 !
2654 !** D-grid u, C-grid v
2655 !
2656  allocate(bc_side%u_BC (is_sn:ie_sn, js_sn:je_sn, klev)) ; bc_side%u_BC=real_snan
2657  allocate(bc_side%vc_BC(is_sn:ie_sn, js_sn:je_sn, klev)) ; bc_side%vc_BC=real_snan
2658 !
2659 !** C-grid u, D-grid v
2660 !
2661  allocate(bc_side%uc_BC(is_we:ie_we, js_we:je_we, klev)) ; bc_side%uc_BC=real_snan
2662  allocate(bc_side%v_BC (is_we:ie_we, js_we:je_we, klev)) ; bc_side%v_BC=real_snan
2663 !
2664 !---------------------------------------------------------------------
2665 !
2666  end subroutine allocate_regional_bc_arrays
2667 !
2668 !---------------------------------------------------------------------
2669 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2670 !---------------------------------------------------------------------
2671 
2672 subroutine remap_scalar_nggps_regional_bc(Atm &
2673  ,side &
2674  ,isd,ied,jsd,jed &
2675  ,is_bc,ie_bc,js_bc,je_bc &
2676  ,km, npz, ncnst, ak0, bk0 &
2677  ,psc, t_in, qa, omga, zh &
2678  ,phis_reg &
2679  ,ps &
2680  ,BC_side )
2682  type(fv_atmos_type), intent(inout) :: Atm
2683  integer, intent(in):: isd,ied,jsd,jed
2684  integer, intent(in):: is_bc,ie_bc,js_bc,je_bc
2685  integer, intent(in):: km & !<-- # of levels in 3-D input variables
2686  ,npz & !<-- # of levels in final 3-D integration variables
2687  ,ncnst
2688  real, intent(in):: ak0(km+1), bk0(km+1)
2689  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc):: psc
2690  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: t_in
2691  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: omga
2692  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km,ncnst):: qa
2693  real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km+1):: zh
2694 !xreal, intent(in), dimension(isd-1:ied+1,jsd-1:jed+1):: phis_reg !<-- Filtered sfc geopotential from preprocessing.
2695  real, intent(inout), dimension(isd-1:ied+1,jsd-1:jed+1):: phis_reg
2696  real, intent(out),dimension(is_bc:ie_bc,js_bc:je_bc) :: ps
2697  character(len=5),intent(in) :: side
2698  type(fv_regional_BC_variables),intent(inout) :: BC_side
2699 
2700 ! local:
2701 !
2702  real, dimension(:,:),allocatable :: pe0
2703  real, dimension(:,:),allocatable :: qn1
2704  real, dimension(:,:),allocatable :: dp2
2705  real, dimension(:,:),allocatable :: pe1
2706  real, dimension(:,:),allocatable :: qp
2707 !
2708  real wk(is_bc:ie_bc,js_bc:je_bc)
2709  real, dimension(is_bc:ie_bc,js_bc:je_bc):: phis
2710 
2711 !!! High-precision
2712  real(kind=R_GRID), dimension(is_bc:ie_bc,npz+1):: pn1
2713  real(kind=R_GRID):: gz_fv(npz+1)
2714  real(kind=R_GRID), dimension(2*km+1):: gz, pn
2715  real(kind=R_GRID), dimension(is_bc:ie_bc,km+1):: pn0
2716  real(kind=R_GRID):: pst
2717 !!! High-precision
2718  integer i,ie,is,j,je,js,k,l,m, k2,iq
2719  integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
2720 !
2721 !---------------------------------------------------------------------------------
2722 !
2723  sphum = get_tracer_index(model_atmos, 'sphum')
2724  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
2725  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
2726  rainwat = get_tracer_index(model_atmos, 'rainwat')
2727  snowwat = get_tracer_index(model_atmos, 'snowwat')
2728  graupel = get_tracer_index(model_atmos, 'graupel')
2729  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
2730  o3mr = get_tracer_index(model_atmos, 'o3mr')
2731 
2732  k2 = max(10, km/2)
2733 
2734  if (mpp_pe()==1) then
2735  print *, 'sphum = ', sphum
2736  print *, 'clwmr = ', liq_wat
2737  print *, ' o3mr = ', o3mr
2738  print *, 'ncnst = ', ncnst
2739  endif
2740 
2741  if ( sphum/=1 ) then
2742  call mpp_error(fatal,'SPHUM must be 1st tracer')
2743  endif
2744 !
2745 !---------------------------------------------------------------------------------
2746 !*** First compute over the extended boundary regions with halo=nhalo_data.
2747 !*** This is needed to obtain pressures that will surround the wind points.
2748 !---------------------------------------------------------------------------------
2749 !
2750  is=is_bc
2751  if(side=='west')then
2752  is=ie_bc-nhalo_data+1
2753  endif
2754 !
2755  ie=ie_bc
2756  if(side=='east')then
2757  ie=is_bc+nhalo_data-1
2758  endif
2759 !
2760  js=js_bc
2761  if(side=='south')then
2762  js=je_bc-nhalo_data+1
2763  endif
2764 !
2765  je=je_bc
2766  if(side=='north')then
2767  je=js_bc+nhalo_data-1
2768  endif
2769 !
2770  allocate(pe0(is:ie,km+1)) ; pe0=real_snan
2771  allocate(qn1(is:ie,npz)) ; qn1=real_snan
2772  allocate(dp2(is:ie,npz)) ; dp2=real_snan
2773  allocate(pe1(is:ie,npz+1)) ; pe1=real_snan
2774  allocate(qp(is:ie,km)) ; qp=real_snan
2775 !
2776 !---------------------------------------------------------------------------------
2777  jloop1: do j=js,je
2778 !---------------------------------------------------------------------------------
2779 !
2780  do k=1,km+1
2781  do i=is,ie
2782  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2783  pn0(i,k) = log(pe0(i,k))
2784  enddo
2785  enddo
2786 
2787  do i=is,ie
2788  do k=1,km+1
2789  pn(k) = pn0(i,k)
2790  gz(k) = zh(i,j,k)*grav
2791  enddo
2792 ! Use log-p for interpolation/extrapolation
2793 ! mirror image method:
2794  do k=km+2, km+k2
2795  l = 2*(km+1) - k
2796  gz(k) = 2.*gz(km+1) - gz(l)
2797  pn(k) = 2.*pn(km+1) - pn(l)
2798  enddo
2799 
2800  do k=km+k2-1, 2, -1
2801  if( phis_reg(i,j).le.gz(k) .and. phis_reg(i,j).ge.gz(k+1) ) then
2802  pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-phis_reg(i,j))/(gz(k)-gz(k+1))
2803  go to 123
2804  endif
2805  enddo
2806  123 ps(i,j) = exp(pst)
2807 
2808  enddo ! i-loop
2809 
2810 !---------------------------------------------------------------------------------
2811  enddo jloop1
2812 !---------------------------------------------------------------------------------
2813 
2814 !---------------------------------------------------------------------------------
2815 !*** Transfer values from the expanded boundary array for sfc pressure into
2816 !*** the Atm object.
2817 !---------------------------------------------------------------------------------
2818 !
2819  is=lbound(atm%ps,1)
2820  ie=ubound(atm%ps,1)
2821  js=lbound(atm%ps,2)
2822  je=ubound(atm%ps,2)
2823 !
2824  do j=js,je
2825  do i=is,ie
2826  atm%ps(i,j)=ps(i,j)
2827  enddo
2828  enddo
2829 !
2830 !---------------------------------------------------------------------------------
2831 !*** Now compute over the normal boundary regions with halo=nhalo_model.
2832 !*** Use the dimensions of one of the permanent BC variables in Atm
2833 !*** as the loop limits so any side of the domain can be addressed.
2834 !---------------------------------------------------------------------------------
2835 !
2836  is=lbound(bc_side%delp_BC,1)
2837  ie=ubound(bc_side%delp_BC,1)
2838  js=lbound(bc_side%delp_BC,2)
2839  je=ubound(bc_side%delp_BC,2)
2840 !
2841 !---------------------------------------------------------------------------------
2842  jloop2: do j=js,je
2843 !---------------------------------------------------------------------------------
2844  do k=1,km+1
2845  do i=is,ie
2846  pe0(i,k) = ak0(k) + bk0(k)*psc(i,j)
2847  pn0(i,k) = log(pe0(i,k))
2848  enddo
2849  enddo
2850 !
2851  do i=is,ie
2852  pe1(i,1) = atm%ak(1)
2853  pn1(i,1) = log(pe1(i,1))
2854  enddo
2855  do k=2,npz+1
2856  do i=is,ie
2857  pe1(i,k) = atm%ak(k) + atm%bk(k)*ps(i,j)
2858  pn1(i,k) = log(pe1(i,k))
2859  enddo
2860  enddo
2861 
2862 ! * Compute delp
2863  do k=1,npz
2864  do i=is,ie
2865  dp2(i,k) = pe1(i,k+1) - pe1(i,k)
2866  bc_side%delp_BC(i,j,k) = dp2(i,k)
2867  enddo
2868  enddo
2869 
2870 ! map shpum, o3mr, liq_wat tracers
2871  do iq=1,ncnst
2872 ! if (iq == sphum .or. iq == liq_wat .or. iq == o3mr) then ! only remap if the data is already set
2873  if (iq /= cld_amt) then ! don't remap cld_amt
2874  do k=1,km
2875  do i=is,ie
2876  qp(i,k) = qa(i,j,k,iq)
2877  enddo
2878  enddo
2879 
2880  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 0, 8, atm%ptop)
2881 
2882  if ( iq==sphum ) then
2883  call fillq(ie-is+1, npz, 1, qn1, dp2)
2884  else
2885  call fillz(ie-is+1, npz, 1, qn1, dp2)
2886  endif
2887 ! The HiRam step of blending model sphum with NCEP data is obsolete because nggps is always cold starting...
2888  do k=1,npz
2889  do i=is,ie
2890  bc_side%q_BC(i,j,k,iq) = qn1(i,k)
2891  enddo
2892  enddo
2893 !jaa endif ! for remapping only the tracers included in the data that is read in
2894  endif ! skip cld_amt in the remap since it is not included in the input
2895  enddo
2896 
2897 !---------------------------------------------------
2898 ! Retrieve temperature using GFS geopotential height
2899 !---------------------------------------------------
2900 !
2901  i_loop: do i=is,ie
2902 !
2903 ! Make sure FV3 top is lower than GFS; can not do extrapolation above the top at this point
2904  if ( pn1(i,1) .lt. pn0(i,1) ) then
2905  call mpp_error(fatal,'FV3 top higher than NCEP/GFS')
2906  endif
2907 
2908  do k=1,km+1
2909  pn(k) = pn0(i,k)
2910  gz(k) = zh(i,j,k)*grav
2911  enddo
2912 !-------------------------------------------------
2913  do k=km+2, km+k2
2914  l = 2*(km+1) - k
2915  gz(k) = 2.*gz(km+1) - gz(l)
2916  pn(k) = 2.*pn(km+1) - pn(l)
2917  enddo
2918 !-------------------------------------------------
2919 
2920  gz_fv(npz+1) = phis_reg(i,j)
2921 
2922  m = 1
2923 
2924  do k=1,npz
2925 ! Searching using FV3 log(pe): pn1
2926 #ifdef USE_ISOTHERMO
2927  do l=m,km
2928  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
2929  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
2930  goto 555
2931  elseif ( pn1(i,k) .gt. pn(km+1) ) then
2932 ! Isothermal under ground; linear in log-p extra-polation
2933  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))
2934  goto 555
2935  endif
2936  enddo
2937 #else
2938  do l=m,km+k2-1
2939  if ( (pn1(i,k).le.pn(l+1)) .and. (pn1(i,k).ge.pn(l)) ) then
2940  gz_fv(k) = gz(l) + (gz(l+1)-gz(l))*(pn1(i,k)-pn(l))/(pn(l+1)-pn(l))
2941  goto 555
2942  endif
2943  enddo
2944 #endif
2945 555 m = l
2946  enddo
2947 
2948 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2949 !xxx DO WE NEED Atm%peln to have values in the boundary region?
2950 !xxx FOR NOW COMMENT IT OUT.
2951 !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2952 !xxx do k=1,npz+1
2953 !xxx Atm%peln(i,k,j) = pn1(i,k)
2954 !xxx enddo
2955 
2956 ! Compute true temperature using hydrostatic balance if not read from input.
2957 
2958  if (data_source /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
2959  do k=1,npz
2960  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)) )
2961  enddo
2962  endif
2963 
2964 
2965  if ( .not. atm%flagstruct%hydrostatic ) then
2966  do k=1,npz
2967  bc_side%delz_BC(i,j,k) = (gz_fv(k+1) - gz_fv(k)) / grav
2968  enddo
2969  endif
2970 
2971  enddo i_loop
2972 
2973 !-----------------------------------------------------------------------
2974 ! seperate cloud water and cloud ice
2975 ! From Jan-Huey Chen's HiRAM code
2976 !-----------------------------------------------------------------------
2977 !
2978 ! If the source is FV3GFS GAUSSIAN NEMSIO FILE then all the tracers are in the boundary files
2979 ! and will be read in.
2980 ! If the source is from old GFS or operational GSM then the tracers will be fixed in the boundaries
2981 ! and may not provide a very good result
2982 !
2983 ! if (cld_amt .gt. 0) BC_side%q_BC(:,:,:,cld_amt) = 0.
2984  if (trim(data_source) /= 'FV3GFS GAUSSIAN NEMSIO FILE') then
2985  if ( atm%flagstruct%nwat .eq. 6 ) then
2986  do k=1,npz
2987  do i=is,ie
2988  qn1(i,k) = bc_side%q_BC(i,j,k,liq_wat)
2989  bc_side%q_BC(i,j,k,rainwat) = 0.
2990  bc_side%q_BC(i,j,k,snowwat) = 0.
2991  bc_side%q_BC(i,j,k,graupel) = 0.
2992  if ( bc_side%pt_BC(i,j,k) > 273.16 ) then ! > 0C all liq_wat
2993  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)
2994  bc_side%q_BC(i,j,k,ice_wat) = 0.
2995 #ifdef ORIG_CLOUDS_PART
2996  else if ( bc_side%pt_BC(i,j,k) < 258.16 ) then ! < -15C all ice_wat
2997  bc_side%q_BC(i,j,k,liq_wat) = 0.
2998  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
2999  else ! between -15~0C: linear interpolation
3000  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-258.16)/15.)
3001  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3002  endif
3003 #else
3004  else if ( bc_side%pt_BC(i,j,k) < 233.16 ) then ! < -40C all ice_wat
3005  bc_side%q_BC(i,j,k,liq_wat) = 0.
3006  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
3007  else
3008  if ( k.eq.1 ) then ! between [-40,0]: linear interpolation
3009  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-233.16)/40.)
3010  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3011  else
3012  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
3013  bc_side%q_BC(i,j,k,liq_wat) = 0.
3014  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k)
3015  else ! between [-40,0]: linear interpolation
3016  bc_side%q_BC(i,j,k,liq_wat) = qn1(i,k)*((bc_side%pt_BC(i,j,k)-233.16)/40.)
3017  bc_side%q_BC(i,j,k,ice_wat) = qn1(i,k) - bc_side%q_BC(i,j,k,liq_wat)
3018  endif
3019  endif
3020  endif
3021 #endif
3022  call mp_auto_conversion(bc_side%q_BC(i,j,k,liq_wat), bc_side%q_BC(i,j,k,rainwat), &
3023  bc_side%q_BC(i,j,k,ice_wat), bc_side%q_BC(i,j,k,snowwat) )
3024  enddo
3025  enddo
3026  endif
3027  endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
3028 !
3029 ! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
3030 ! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
3031 ! For FV3GFS NEMSIO input, w is already in m/s (but the code reads in as omga) and just needs to be remapped
3032 !-------------------------------------------------------------
3033 ! map omega
3034 !------- ------------------------------------------------------
3035  if ( .not. atm%flagstruct%hydrostatic ) then
3036  do k=1,km
3037  do i=is,ie
3038  qp(i,k) = omga(i,j,k)
3039  enddo
3040  enddo
3041 
3042  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, -1, 4, atm%ptop)
3043 
3044  if (data_source == 'FV3GFS GAUSSIAN NEMSIO FILE') then
3045  do k=1,npz
3046  do i=is,ie
3047  bc_side%w_BC(i,j,k) = qn1(i,k)
3048  enddo
3049  enddo
3050 !------------------------------
3051 ! Remap input T linearly in p.
3052 !------------------------------
3053  do k=1,km
3054  do i=is,ie
3055  qp(i,k) = t_in(i,j,k)
3056  enddo
3057  enddo
3058 
3059  call mappm(km, pe0, qp, npz, pe1, qn1, is,ie, 2, 4, atm%ptop)
3060 
3061  do k=1,npz
3062  do i=is,ie
3063  bc_side%pt_BC(i,j,k) = qn1(i,k)
3064  enddo
3065  enddo
3066 
3067  else
3068  do k=1,npz
3069  do i=is,ie
3070  bc_side%w_BC(i,j,k) = qn1(i,k)/bc_side%delp_BC(i,j,k)*bc_side%delz_BC(i,j,k)
3071  enddo
3072  enddo
3073  endif
3074 
3075  endif !.not. Atm%flagstruct%hydrostatic
3076 
3077  enddo jloop2
3078 
3079 ! Add some diagnostics:
3080 !xxxcall p_maxmin('PS_model (mb)', Atm%ps(is:ie,js:je), is, ie, js, je, 1, 0.01)
3081 !xxxcall p_maxmin('PT_model', Atm%pt(is:ie,js:je,1:npz), is, ie, js, je, npz, 1.)
3082  do j=js,je
3083  do i=is,ie
3084  wk(i,j) = phis_reg(i,j)/grav - zh(i,j,km+1)
3085  enddo
3086  enddo
3087 !xxxcall pmaxmn('ZS_diff (m)', wk, is, ie, js, je, 1, 1., Atm%gridstruct%area_64, Atm%domain)
3088 
3089  do j=js,je
3090  do i=is,ie
3091  wk(i,j) = ps(i,j) - psc(i,j)
3092  enddo
3093  enddo
3094 !xxxcall pmaxmn('PS_diff (mb)', wk, is, ie, js, je, 1, 0.01, Atm%gridstruct%area_64, Atm%domain)
3095  deallocate (pe0,qn1,dp2,pe1,qp)
3096  if (is_master()) write(*,*) 'done remap_scalar_nggps_regional_bc'
3097 !---------------------------------------------------------------------
3098 
3099  end subroutine remap_scalar_nggps_regional_bc
3100 
3101 !---------------------------------------------------------------------
3102 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3103 !---------------------------------------------------------------------
3104 
3105  subroutine remap_dwinds_regional_bc(Atm &
3106  ,is_input,ie_input &
3107  ,js_input,je_input &
3108  ,is_u,ie_u,js_u,je_u &
3109  ,is_v,ie_v,js_v,je_v &
3110  ,km, npz &
3111  ,ak0, bk0 &
3112  ,psc, ud, vd, uc, vc &
3113  ,BC_side )
3114  type(fv_atmos_type), intent(inout) :: Atm
3115  integer, intent(in):: is_input, ie_input, js_input, je_input
3116  integer, intent(in):: is_u,ie_u,js_u,je_u
3117  integer, intent(in):: is_v,ie_v,js_v,je_v
3118  integer, intent(in):: km & !<-- # of levels in 3-D input variables
3119  ,npz
3120  real, intent(in):: ak0(km+1), bk0(km+1)
3121 
3122  real, intent(in) :: psc(is_input:ie_input,js_input:je_input)
3123 
3124  real, intent(in):: ud(is_u:ie_u,js_u:je_u,km)
3125  real, intent(in):: vc(is_u:ie_u,js_u:je_u,km)
3126  real, intent(in):: vd(is_v:ie_v,js_v:je_v,km)
3127  real, intent(in):: uc(is_v:ie_v,js_v:je_v,km)
3128  type(fv_regional_BC_variables),intent(inout) :: BC_side
3129 ! local:
3130  real, dimension(:,:),allocatable :: pe0
3131  real, dimension(:,:),allocatable :: pe1
3132  real, dimension(:,:),allocatable :: qn1_d,qn1_c
3133  integer i,j,k
3134 
3135  allocate(pe0(is_u:ie_u, km+1)) ; pe0=real_snan
3136  allocate(pe1(is_u:ie_u, npz+1)) ; pe1=real_snan
3137  allocate(qn1_d(is_u:ie_u, npz)) ; qn1_d=real_snan
3138  allocate(qn1_c(is_u:ie_u, npz)) ; qn1_c=real_snan
3139 
3140 !----------------------------------------------------------------------------------------------
3141  j_loopu: do j=js_u,je_u
3142 !----------------------------------------------------------------------------------------------
3143 
3144 !------
3145 ! map u
3146 !------
3147  do k=1,km+1
3148  do i=is_u,ie_u
3149  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j))
3150  enddo
3151  enddo
3152  do k=1,npz+1
3153  do i=is_u,ie_u
3154  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j))
3155  enddo
3156  enddo
3157  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), &
3158  qn1_d(is_u:ie_u,1:npz), is_u,ie_u, -1, 8, atm%ptop )
3159  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), &
3160  qn1_c(is_u:ie_u,1:npz), is_u,ie_u, -1, 8, atm%ptop )
3161  do k=1,npz
3162  do i=is_u,ie_u
3163  bc_side%u_BC(i,j,k) = qn1_d(i,k)
3164  bc_side%vc_BC(i,j,k) = qn1_c(i,k)
3165  enddo
3166  enddo
3167 
3168  enddo j_loopu
3169 
3170  deallocate(pe0)
3171  deallocate(pe1)
3172  deallocate(qn1_d)
3173  deallocate(qn1_c)
3174 
3175  allocate(pe0(is_v:ie_v, km+1)) ; pe0=real_snan
3176  allocate(pe1(is_v:ie_v, npz+1)) ; pe1=real_snan
3177  allocate(qn1_d(is_v:ie_v, npz)) ; qn1_d=real_snan
3178  allocate(qn1_c(is_v:ie_v, npz)) ; qn1_c=real_snan
3179 
3180 !----------------------------------------------------------------------------------------------
3181  j_loopv: do j=js_v,je_v
3182 !----------------------------------------------------------------------------------------------
3183 !
3184 !------
3185 ! map v
3186 !------
3187 
3188  do k=1,km+1
3189  do i=is_v,ie_v
3190  pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j))
3191  enddo
3192  enddo
3193  do k=1,npz+1
3194  do i=is_v,ie_v
3195  pe1(i,k) = atm%ak(k) + atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j))
3196  enddo
3197  enddo
3198  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), &
3199  qn1_d(is_v:ie_v,1:npz), is_v,ie_v, -1, 8, atm%ptop)
3200  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), &
3201  qn1_c(is_v:ie_v,1:npz), is_v,ie_v, -1, 8, atm%ptop)
3202  do k=1,npz
3203  do i=is_v,ie_v
3204  bc_side%v_BC(i,j,k) = qn1_d(i,k)
3205  bc_side%uc_BC(i,j,k) = qn1_c(i,k)
3206  enddo
3207  enddo
3208 
3209  enddo j_loopv
3210 
3211  deallocate(pe0)
3212  deallocate(pe1)
3213  deallocate(qn1_d)
3214  deallocate(qn1_c)
3215 
3216  if (is_master()) write(*,*) 'done remap_dwinds'
3217 
3218  end subroutine remap_dwinds_regional_bc
3219 
3220 !---------------------------------------------------------------------
3221 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3222 !---------------------------------------------------------------------
3223 
3224  subroutine set_regional_bcs(delp,delz,w,pt &
3225 #ifdef USE_COND
3226  ,q_con &
3227 #endif
3228 #ifdef MOIST_CAPPA
3229  ,cappa &
3230 #endif
3231  ,q &
3232  ,u,v,uc,vc &
3233  ,bd, nlayers &
3234  ,fcst_time )
3235 !
3236 !---------------------------------------------------------------------
3237 !*** Select the given variable's boundary data at the two
3238 !*** bracketing time levels and apply them to the updating
3239 !*** of the variable's boundary region at the appropriate
3240 !*** forecast time.
3241 !---------------------------------------------------------------------
3242  implicit none
3243 !---------------------------------------------------------------------
3244 !
3245 !--------------------
3246 !*** Input variables
3247 !--------------------
3248 !
3249  integer,intent(in) :: nlayers
3250 !
3251  real,intent(in) :: fcst_time
3252 !
3253  type(fv_grid_bounds_type),intent(in) :: bd
3254 !
3255 !----------------------
3256 !*** Output variables
3257 !----------------------
3258 !
3259  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),intent(out) :: &
3260  delp &
3261  ,pt
3262 !
3263  real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: delz,w
3264 #ifdef USE_COND
3265  real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: q_con
3266 #endif
3267 
3268 !
3269  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz,ntracers),intent(out) :: q
3270 !
3271 #ifdef MOIST_CAPPA
3272  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),intent(out) :: cappa
3273 !#else
3274 ! real,dimension(isd:isd,jsd:jsd,1),intent(out) :: cappa
3275 #endif
3276 !
3277  real,dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz),intent(out) :: u,vc
3278 !
3279  real,dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz),intent(out) :: uc,v
3280 !
3281 !---------------------
3282 !*** Local variables
3283 !---------------------
3284 !
3285  real :: fraction_interval
3286 !
3287 !---------------------------------------------------------------------
3288 !*********************************************************************
3289 !---------------------------------------------------------------------
3290 !
3291 !---------------------------------------------------------------------
3292 !*** The current forecast time is this fraction of the way from
3293 !*** time level 0 to time level 1.
3294 !---------------------------------------------------------------------
3295 !
3296  fraction_interval=mod(fcst_time,(bc_update_interval*3600.)) &
3297  /(bc_update_interval*3600.)
3298 !
3299 !---------------------------------------------------------------------
3300 !
3301  if(north_bc)then
3302  call bc_values_into_arrays(bc_t0%north,bc_t1%north &
3303  ,'north' &
3304  ,bd%isd &
3305  ,bd%ied &
3306  ,bd%jsd &
3307  ,bd%js-1 &
3308  ,bd%isd &
3309  ,bd%ied &
3310  ,bd%jsd &
3311  ,bd%js-1 &
3312  ,bd%isd &
3313  ,bd%ied+1 &
3314  ,bd%jsd &
3315  ,bd%js-1)
3316  endif
3317 !
3318  if(south_bc)then
3319  call bc_values_into_arrays(bc_t0%south,bc_t1%south &
3320  ,'south' &
3321  ,bd%isd &
3322  ,bd%ied &
3323  ,bd%je+1 &
3324  ,bd%jed &
3325  ,bd%isd &
3326  ,bd%ied &
3327  ,bd%je+2 &
3328  ,bd%jed+1 &
3329  ,bd%isd &
3330  ,bd%ied+1 &
3331  ,bd%je+1 &
3332  ,bd%jed )
3333  endif
3334 !
3335  if(east_bc)then
3336  call bc_values_into_arrays(bc_t0%east,bc_t1%east &
3337  ,'east ' &
3338  ,bd%isd &
3339  ,bd%is-1 &
3340  ,bd%js &
3341  ,bd%je &
3342  ,bd%isd &
3343  ,bd%is-1 &
3344  ,bd%js &
3345  ,bd%je+1 &
3346  ,bd%isd &
3347  ,bd%is-1 &
3348  ,bd%js &
3349  ,bd%je )
3350  endif
3351 !
3352  if(west_bc)then
3353  call bc_values_into_arrays(bc_t0%west,bc_t1%west &
3354  ,'west ' &
3355  ,bd%ie+1 &
3356  ,bd%ied &
3357  ,bd%js &
3358  ,bd%je &
3359  ,bd%ie+1 &
3360  ,bd%ied &
3361  ,bd%js &
3362  ,bd%je+1 &
3363  ,bd%ie+2 &
3364  ,bd%ied+1 &
3365  ,bd%js &
3366  ,bd%je )
3367  endif
3368 !
3369 !---------------------------------------------------------------------
3370 
3371  contains
3372 
3373 !---------------------------------------------------------------------
3374 !
3375  subroutine bc_values_into_arrays(side_t0,side_t1 &
3376  ,side &
3377  ,i1,i2,j1,j2 &
3378  ,i1_uvs,i2_uvs,j1_uvs,j2_uvs &
3379  ,i1_uvw,i2_uvw,j1_uvw,j2_uvw )
3381 !---------------------------------------------------------------------
3382 !*** Apply boundary values to the prognostic arrays at the
3383 !*** desired time.
3384 !---------------------------------------------------------------------
3385  implicit none
3386 !---------------------------------------------------------------------
3387 !
3388 !---------------------
3389 !*** Input arguments
3390 !---------------------
3391 !
3392  type(fv_regional_BC_variables),intent(in) :: side_t0 &
3393  ,side_t1
3394 !
3395  character(len=*),intent(in) :: side
3396 !
3397  integer,intent(in) :: i1,i2,j1,j2 &
3398  ,i1_uvs,i2_uvs,j1_uvs,j2_uvs &
3399  ,i1_uvw,i2_uvw,j1_uvw,j2_uvw
3400 !
3401 !---------------------
3402 !*** Local arguments
3403 !---------------------
3404 !
3405  integer :: i,ie,j,je,jend,jend_uvs,jend_uvw &
3406  ,jstart,jstart_uvs,jstart_uvw,k,nt,nz
3407 !
3408 !---------------------------------------------------------------------
3409 !*********************************************************************
3410 !---------------------------------------------------------------------
3411 !
3412  jstart=j1
3413  jend =j2
3414  jstart_uvs=j1_uvs
3415  jend_uvs =j2_uvs
3416  jstart_uvw=j1_uvw
3417  jend_uvw =j2_uvw
3418  if((trim(side)=='east'.or.trim(side)=='west').and..not.north_bc)then
3419  jstart=j1-nhalo_model
3420  jstart_uvs=j1_uvs-nhalo_model
3421  jstart_uvw=j1_uvw-nhalo_model
3422  endif
3423  if((trim(side)=='east'.or.trim(side)=='west').and..not.south_bc)then
3424  jend=j2+nhalo_model
3425  jend_uvs=j2_uvs+nhalo_model
3426  jend_uvw=j2_uvw+nhalo_model
3427  endif
3428 !
3429  do k=1,nlayers
3430  do j=jstart,jend
3431  do i=i1,i2
3432  delp(i,j,k)=side_t0%delp_BC(i,j,k) &
3433  +(side_t1%delp_BC(i,j,k)-side_t0%delp_BC(i,j,k)) &
3434  *fraction_interval
3435  pt(i,j,k)=side_t0%pt_BC(i,j,k) &
3436  +(side_t1%pt_BC(i,j,k)-side_t0%pt_BC(i,j,k)) &
3437  *fraction_interval
3438 #ifdef MOIST_CAPPA
3439  cappa(i,j,k)=side_t0%cappa_BC(i,j,k) &
3440  +(side_t1%cappa_BC(i,j,k)-side_t0%cappa_BC(i,j,k)) &
3441  *fraction_interval
3442 #endif
3443  enddo
3444  enddo
3445 !
3446  do j=jstart_uvs,jend_uvs
3447  do i=i1_uvs,i2_uvs
3448  u(i,j,k)=side_t0%u_BC(i,j,k) &
3449  +(side_t1%u_BC(i,j,k)-side_t0%u_BC(i,j,k)) &
3450  *fraction_interval
3451  vc(i,j,k)=side_t0%vc_BC(i,j,k) &
3452  +(side_t1%vc_BC(i,j,k)-side_t0%vc_BC(i,j,k)) &
3453  *fraction_interval
3454  enddo
3455  enddo
3456 !
3457  do j=jstart_uvw,jend_uvw
3458  do i=i1_uvw,i2_uvw
3459  v(i,j,k)=side_t0%v_BC(i,j,k) &
3460  +(side_t1%v_BC(i,j,k)-side_t0%v_BC(i,j,k)) &
3461  *fraction_interval
3462  uc(i,j,k)=side_t0%uc_BC(i,j,k) &
3463  +(side_t1%uc_BC(i,j,k)-side_t0%uc_BC(i,j,k)) &
3464  *fraction_interval
3465  enddo
3466  enddo
3467  enddo
3468 !
3469  ie=min(ubound(side_t0%delz_BC,1),ubound(delz,1))
3470  je=min(ubound(side_t0%delz_BC,2),ubound(delz,2))
3471  nz=ubound(delz,3)
3472 !
3473  do k=1,nz
3474  do j=jstart,jend
3475  do i=i1,ie
3476  delz(i,j,k)=side_t0%delz_BC(i,j,k) &
3477  +(side_t1%delz_BC(i,j,k)-side_t0%delz_BC(i,j,k)) &
3478  *fraction_interval
3479 #ifdef USE_COND
3480  q_con(i,j,k)=side_t0%q_con_BC(i,j,k) &
3481  +(side_t1%q_con_BC(i,j,k)-side_t0%q_con_BC(i,j,k)) &
3482  *fraction_interval
3483 #endif
3484  w(i,j,k)=side_t0%w_BC(i,j,k) &
3485  +(side_t1%w_BC(i,j,k)-side_t0%w_BC(i,j,k)) &
3486  *fraction_interval
3487  enddo
3488  enddo
3489  enddo
3490 !
3491  do nt=1,ntracers
3492  do k=1,nz
3493  do j=jstart,jend
3494  do i=i1,i2
3495  q(i,j,k,nt)=side_t0%q_BC(i,j,k,nt) &
3496  +(side_t1%q_BC(i,j,k,nt)-side_t0%q_BC(i,j,k,nt)) &
3497  *fraction_interval
3498  enddo
3499  enddo
3500  enddo
3501  enddo
3502 !
3503 !---------------------------------------------------------------------
3504 !
3505  end subroutine bc_values_into_arrays
3506 !
3507 !---------------------------------------------------------------------
3508 !
3509  end subroutine set_regional_bcs
3510 !
3511 !---------------------------------------------------------------------
3512 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3513 !---------------------------------------------------------------------
3514  subroutine regional_boundary_update(array &
3515  ,bc_vbl_name &
3516  ,lbnd_x,ubnd_x &
3517  ,lbnd_y,ubnd_y &
3518  ,ubnd_z &
3519  ,is,ie,js,je &
3520  ,isd,ied,jsd,jed &
3521  ,fcst_time &
3522  ,index4 )
3524 !---------------------------------------------------------------------
3525 !*** Select the given variable's boundary data at the two
3526 !*** bracketing time levels and apply them to the updating
3527 !*** of the variable's boundary region at the appropriate
3528 !*** forecast time.
3529 !---------------------------------------------------------------------
3530  implicit none
3531 !---------------------------------------------------------------------
3532 !
3533 !--------------------
3534 !*** Input variables
3535 !--------------------
3536 !
3537  integer,intent(in) :: lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z
3538 !
3539  integer,intent(in) :: is,ie,js,je & !<-- Compute limits
3540  ,isd,ied,jsd,jed
3541 !
3542  integer,intent(in),optional :: index4
3543 !
3544  real,intent(in) :: fcst_time
3545 !
3546  character(len=*),intent(in) :: bc_vbl_name
3547 !
3548 !----------------------
3549 !*** Output variables
3550 !----------------------
3551 !
3552  real,dimension(lbnd_x:ubnd_x,lbnd_y:ubnd_y,1:ubnd_z) &
3553  ,intent(out) :: array
3554 !
3555 !---------------------
3556 !*** Local variables
3557 !---------------------
3558 !
3559  integer :: i1,i2,j1,j2
3560  integer :: lbnd1,ubnd1,lbnd2,ubnd2
3561  integer :: iq
3562  integer :: nside
3563 !
3564  real,dimension(:,:,:),pointer :: bc_t0,bc_t1
3565 !
3566  logical :: call_interp
3567 !
3568 !---------------------------------------------------------------------
3569 !*********************************************************************
3570 !---------------------------------------------------------------------
3571 !
3572  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
3573  return
3574  endif
3575 !
3576  iq=0
3577  if(present(index4))then
3578  iq=index4
3579  endif
3580 !
3581 !---------------------------------------------------------------------
3582 !*** Loop through the sides of the domain and find the limits
3583 !*** of the region to update in the boundary.
3584 !---------------------------------------------------------------------
3585 !
3586 !---------------------------------------------------------------------
3587  sides: do nside=1,4
3588 !---------------------------------------------------------------------
3589 !
3590  call_interp=.false.
3591 !
3592 !-----------
3593 !*** North
3594 !-----------
3595 !
3596  if(nside==1)then
3597  if(north_bc)then
3598  call_interp=.true.
3601 !
3602  i1=isd
3603  i2=ied
3604  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
3605  i2=ied+1
3606  endif
3607 !
3608  j1=jsd
3609  j2=js-1
3610  endif
3611  endif
3612 !
3613 !-----------
3614 !*** South
3615 !-----------
3616 !
3617  if(nside==2)then
3618  if(south_bc)then
3619  call_interp=.true.
3622 !
3623  i1=isd
3624  i2=ied
3625  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
3626  i2=ied+1
3627  endif
3628 !
3629  j1=je+1
3630  j2=jed
3631  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
3632  j1=je+2
3633  j2=jed+1
3634  endif
3635  endif
3636  endif
3637 !
3638 !----------
3639 !*** East
3640 !----------
3641 !
3642  if(nside==3)then
3643  if(east_bc)then
3644  call_interp=.true.
3647 !
3648  j1=jsd
3649  j2=jed
3650 !
3651  i1=isd
3652  i2=is-1
3653 !
3654  if(north_bc)then
3655  j1=js
3656  endif
3657  if(south_bc)then
3658  j2=je
3659  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
3660  j2=je+1
3661  endif
3662  endif
3663  endif
3664  endif
3665 !
3666 !----------
3667 !*** West
3668 !----------
3669 !
3670  if(nside==4)then
3671  if(west_bc)then
3672  call_interp=.true.
3675 !
3676  j1=jsd
3677  j2=jed
3678 !
3679  i1=ie+1
3680  i2=ied
3681  if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v')then
3682  i1=ie+2
3683  i2=ied+1
3684  endif
3685 !
3686  if(north_bc)then
3687  j1=js
3688  endif
3689  if(south_bc)then
3690  j2=je
3691  if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc')then
3692  j2=je+1
3693  endif
3694  endif
3695  endif
3696  endif
3697 !
3698 !---------------------------------------------------------------------
3699 !*** Get the pointers pointing at the boundary arrays holding the
3700 !*** two time levels of the given prognostic array's boundary region
3701 !*** then update the boundary points.
3702 !---------------------------------------------------------------------
3703 !
3704  if(call_interp)then
3705 !
3706  call retrieve_bc_variable_data(bc_vbl_name &
3708  ,bc_t0,bc_t1 &
3709  ,lbnd1,ubnd1,lbnd2,ubnd2 &
3710  ,iq )
3711 !
3712  call bc_time_interpolation(array &
3713  ,lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z &
3714  ,bc_t0,bc_t1 &
3715  ,lbnd1,ubnd1,lbnd2,ubnd2 &
3716  ,i1,i2,j1,j2 &
3717  ,fcst_time &
3719  endif
3720 !
3721 !---------------------------------------------------------------------
3722  enddo sides
3723 !---------------------------------------------------------------------
3724 !
3725 !---------------------------------------------------------------------
3726 
3727  end subroutine regional_boundary_update
3728 
3729 !---------------------------------------------------------------------
3730 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3731 !---------------------------------------------------------------------
3732 
3733  subroutine retrieve_bc_variable_data(bc_vbl_name &
3734  ,bc_side_t0,bc_side_t1 &
3735  ,bc_t0,bc_t1 &
3736  ,lbnd1,ubnd1,lbnd2,ubnd2 &
3737  ,iq )
3739 !---------------------------------------------------------------------
3740 !*** Select the boundary variable associated with the prognostic
3741 !*** array that needs its boundary region to be updated.
3742 !---------------------------------------------------------------------
3743  implicit none
3744 !---------------------------------------------------------------------
3745 !
3746 !---------------------
3747 !*** Input variables
3748 !---------------------
3749 !
3750  integer,intent(in) :: iq
3751 !
3752  character(len=*),intent(in) :: bc_vbl_name
3753 !
3754  type(fv_regional_BC_variables),pointer :: bc_side_t0,bc_side_t1
3755 !
3756 !
3757 !----------------------
3758 !*** Output variables
3759 !----------------------
3760 !
3761  integer,intent(out) :: lbnd1,ubnd1,lbnd2,ubnd2
3762 !
3763  real,dimension(:,:,:),pointer :: bc_t0,bc_t1
3764 !
3765 !---------------------------------------------------------------------
3766 !*********************************************************************
3767 !---------------------------------------------------------------------
3768 !
3769  select case (bc_vbl_name)
3770 !
3771  case ('delp')
3772  bc_t0=>bc_side_t0%delp_BC
3773  bc_t1=>bc_side_t1%delp_BC
3774  case ('delz')
3775  bc_t0=>bc_side_t0%delz_BC
3776  bc_t1=>bc_side_t1%delz_BC
3777  case ('pt')
3778  bc_t0=>bc_side_t0%pt_BC
3779  bc_t1=>bc_side_t1%pt_BC
3780  case ('w')
3781  bc_t0=>bc_side_t0%w_BC
3782  bc_t1=>bc_side_t1%w_BC
3783  case ('divgd')
3784  bc_t0=>bc_side_t0%divgd_BC
3785  bc_t1=>bc_side_t1%divgd_BC
3786 #ifdef MOIST_CAPPA
3787  case ('cappa')
3788  bc_t0=>bc_side_t0%cappa_BC
3789  bc_t1=>bc_side_t1%cappa_BC
3790 #endif
3791 #ifdef USE_COND
3792  case ('q_con')
3793  bc_t0=>bc_side_t0%q_con_BC
3794  bc_t1=>bc_side_t1%q_con_BC
3795 #endif
3796  case ('q')
3797  if(iq<1)then
3798  write(0,101)
3799  101 format(' iq<1 is not a valid index for q_BC array in retrieve_bc_variable_data')
3800  endif
3801  lbnd1=lbound(bc_side_t0%q_BC,1)
3802  lbnd2=lbound(bc_side_t0%q_BC,2)
3803  ubnd1=ubound(bc_side_t0%q_BC,1)
3804  ubnd2=ubound(bc_side_t0%q_BC,2)
3805  bc_t0=>bc_side_t0%q_BC(:,:,:,iq)
3806  bc_t1=>bc_side_t1%q_BC(:,:,:,iq)
3807  case ('u')
3808  bc_t0=>bc_side_t0%u_BC
3809  bc_t1=>bc_side_t1%u_BC
3810  case ('v')
3811  bc_t0=>bc_side_t0%v_BC
3812  bc_t1=>bc_side_t1%v_BC
3813  case ('uc')
3814  bc_t0=>bc_side_t0%uc_BC
3815  bc_t1=>bc_side_t1%uc_BC
3816  case ('vc')
3817  bc_t0=>bc_side_t0%vc_BC
3818  bc_t1=>bc_side_t1%vc_BC
3819 !
3820  end select
3821 !
3822  if(trim(bc_vbl_name)/='q')then
3823  lbnd1=lbound(bc_t0,1)
3824  lbnd2=lbound(bc_t0,2)
3825  ubnd1=ubound(bc_t0,1)
3826  ubnd2=ubound(bc_t0,2)
3827  endif
3828 !
3829 !---------------------------------------------------------------------
3830 !
3831  end subroutine retrieve_bc_variable_data
3832 !
3833 !---------------------------------------------------------------------
3834 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3835 !---------------------------------------------------------------------
3836 !
3837  subroutine bc_time_interpolation(array &
3838  ,lbnd_x, ubnd_x &
3839  ,lbnd_y, ubnd_y &
3840  ,ubnd_z &
3841  ,bc_t0, bc_t1 &
3842  ,lbnd1, ubnd1 &
3843  ,lbnd2, ubnd2 &
3844  ,i1,i2,j1,j2 &
3845  ,fcst_time &
3846  ,bc_update_interval )
3848 !---------------------------------------------------------------------
3849 !*** Update the boundary region of the input array at the given
3850 !*** forecast time that is within the interval bracketed by the
3851 !*** two current boundary region states.
3852 !---------------------------------------------------------------------
3853  implicit none
3854 !---------------------------------------------------------------------
3855 !
3856 !---------------------
3857 !*** Input variables
3858 !---------------------
3859 !
3860  integer,intent(in) :: lbnd_x,ubnd_x,lbnd_y,ubnd_y,ubnd_z
3861 !
3862  integer,intent(in) :: lbnd1,ubnd1,lbnd2,ubnd2
3863 !
3864  integer,intent(in) :: i1,i2,j1,j2
3865 !
3866  integer,intent(in) :: bc_update_interval
3867 !
3868  real,intent(in) :: fcst_time
3869 !
3870  real,dimension(lbnd1:ubnd1,lbnd2:ubnd2,1:ubnd_z) :: bc_t0 & !<-- Interpolate between these
3871  ,bc_t1 ! two boundary region states.
3872 !
3873 !---------------------
3874 !*** Output variables
3875 !---------------------
3876 !
3877  real,dimension(lbnd_x:ubnd_x,lbnd_y:ubnd_y,1:ubnd_z) &
3878  ,intent(out) :: array
3879 !
3880 !---------------------
3881 !*** Local variables
3882 !---------------------
3883 !
3884  integer :: i,j,k
3885 !
3886  real :: fraction_interval
3887 !
3888 !---------------------------------------------------------------------
3889 !*********************************************************************
3890 !---------------------------------------------------------------------
3891 !
3892 !---------------------------------------------------------------------
3893 !*** The current forecast time is this fraction of the way from
3894 !*** time level 0 to time level 1.
3895 !---------------------------------------------------------------------
3896 !
3897  fraction_interval=mod(fcst_time,(bc_update_interval*3600.)) &
3898  /(bc_update_interval*3600.)
3899 !
3900 !---------------------------------------------------------------------
3901 !
3902  do k=1,ubnd_z
3903  do j=j1,j2
3904  do i=i1,i2
3905  array(i,j,k)=bc_t0(i,j,k) &
3906  +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval
3907  enddo
3908  enddo
3909  enddo
3910 !
3911 !---------------------------------------------------------------------
3912 
3913  end subroutine bc_time_interpolation
3914 !---------------------------------------------------------------------
3915 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3916 !---------------------------------------------------------------------
3917 !
3918  subroutine bc_time_interpolation_general(is,ie,js,je &
3919  ,is_s,ie_s,js_s,je_s &
3920  ,is_w,ie_w,js_w,je_w &
3921  ,fraction &
3922  ,t0,t1 &
3923  ,Atm )
3925 !---------------------------------------------------------------------
3926 !*** Execute the linear time interpolation between t0 and t1
3927 !*** generically for any side of the regional domain's boundary
3928 !*** region.
3929 !---------------------------------------------------------------------
3930  implicit none
3931 !---------------------------------------------------------------------
3932 !
3933 !------------------------
3934 !*** Argument variables
3935 !------------------------
3936 !
3937  integer,intent(in) :: is,ie,js,je & !<-- Index limits for centers of grid cells
3938  ,is_s,ie_s,js_s,je_s & !<-- Index limits for south/north edges of grid cells
3939  ,is_w,ie_w,js_w,je_w
3940 !
3941  real,intent(in) :: fraction
3942 !
3943  type(fv_regional_BC_variables),intent(in) :: t0,t1
3944 !
3945  type(fv_atmos_type),intent(inout) :: Atm
3946 !
3947 !---------------------
3948 !*** Local variables
3949 !---------------------
3950 !
3951  integer :: i,j,k,n,nlayers!,ntracers
3952 !
3953 !---------------------------------------------------------------------
3954 !*********************************************************************
3955 !---------------------------------------------------------------------
3956 !
3957  nlayers =atm%npz
3958 ! ntracers=Atm%ncnst !<-- # of advected tracers
3959 !
3960 !---------------------------------------------------------------------
3961 !
3962  k_loop: do k=1,nlayers
3963 !
3964 !---------------------------------------------------------------------
3965 !
3966 !-------------
3967 !*** Scalars
3968 !-------------
3969 !
3970  do j=js,je
3971  do i=is,ie
3972 !
3973  atm%delp(i,j,k)=t0%delp_BC(i,j,k) &
3974  +(t1%delp_BC(i,j,k)-t0%delp_BC(i,j,k)) &
3975  *fraction
3976 !
3977 #ifndef SW_DYNAMICS
3978  atm%delz(i,j,k)=t0%delz_BC(i,j,k) &
3979  +(t1%delz_BC(i,j,k)-t0%delz_BC(i,j,k)) &
3980  *fraction
3981 !
3982  atm%w(i,j,k)=t0%w_BC(i,j,k) &
3983  +(t1%w_BC(i,j,k)-t0%w_BC(i,j,k)) &
3984  *fraction
3985 !
3986  atm%pt(i,j,k)=t0%pt_BC(i,j,k) &
3987  +(t1%pt_BC(i,j,k)-t0%pt_BC(i,j,k)) &
3988  *fraction
3989 #ifdef USE_COND
3990  atm%q_con(i,j,k)=t0%q_con_BC(i,j,k) &
3991  +(t1%q_con_BC(i,j,k)-t0%q_con_BC(i,j,k)) &
3992  *fraction
3993 #ifdef MOIST_CAPPA
3994 ! Atm%cappa(i,j,k)=t0%pt_BC(i,j,k) & !<-- Update cappa.
3995 ! +(t1%cappa_BC(i,j,k)-t0%cappa_BC(i,j,k)) &
3996 ! *fraction
3997 #endif
3998 #endif
3999 #endif
4000 !
4001  enddo
4002  enddo
4003 !
4004  do n=1,ntracers
4005 !
4006  do j=js,je
4007  do i=is,ie
4008  atm%q(i,j,k,n)=t0%q_BC(i,j,k,n) &
4009  +(t1%q_BC(i,j,k,n)-t0%q_BC(i,j,k,n)) &
4010  *fraction
4011  enddo
4012  enddo
4013 !
4014  enddo
4015 !
4016 !-----------
4017 !*** Winds
4018 !-----------
4019 !
4020  do j=js_s,je_s
4021  do i=is_s,ie_s
4022  atm%u(i,j,k)=t0%u_BC(i,j,k) &
4023  +(t1%u_BC(i,j,k)-t0%u_BC(i,j,k)) &
4024  *fraction
4025  atm%vc(i,j,k)=t0%vc_BC(i,j,k) &
4026  +(t1%vc_BC(i,j,k)-t0%vc_BC(i,j,k)) &
4027  *fraction
4028  enddo
4029  enddo
4030 !
4031 !
4032  do j=js_w,je_w
4033  do i=is_w,ie_w
4034  atm%v(i,j,k)=t0%v_BC(i,j,k) &
4035  +(t1%v_BC(i,j,k)-t0%v_BC(i,j,k)) &
4036  *fraction
4037  atm%uc(i,j,k)=t0%uc_BC(i,j,k) &
4038  +(t1%uc_BC(i,j,k)-t0%uc_BC(i,j,k)) &
4039  *fraction
4040  enddo
4041  enddo
4042 !
4043 !---------------------------------------------------------------------
4044 !
4045  enddo k_loop
4046 !
4047 !---------------------------------------------------------------------
4048 !
4049  end subroutine bc_time_interpolation_general
4050 !
4051 !---------------------------------------------------------------------
4052 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4053 !---------------------------------------------------------------------
4054 !
4055  subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 &
4056  ,nlev,ntracers,bnds )
4058 !---------------------------------------------------------------------
4059 !*** BC data has been read into the time level t1 object. Now
4060 !*** move the t1 data into the t1 object before refilling t1
4061 !*** with the next data from the BC file.
4062 !---------------------------------------------------------------------
4063  implicit none
4064 !---------------------------------------------------------------------
4065 !
4066 !------------------------
4067 !*** Argument variables
4068 !------------------------
4069 !
4070  integer,intent(in) :: nlev & !<-- # of model layers.
4071  ,ntracers
4072 !
4073  type(fv_regional_bc_bounds_type),intent(in) :: bnds
4074 !
4075  type(fv_domain_sides),target,intent(in) :: bc_t1
4076 !
4077  type(fv_domain_sides),target,intent(inout) :: bc_t0
4078 !
4079 !---------------------
4080 !*** Local variables
4081 !---------------------
4082 !
4083  integer :: i,ie_c,ie_s,ie_w,is_c,is_s,is_w &
4084  ,j,je_c,je_s,je_w,js_c,js_s,js_w &
4085  ,k,n,nside
4086 !
4087  logical :: move
4088 !
4089 !---------------------------------------------------------------------
4090 !*********************************************************************
4091 !---------------------------------------------------------------------
4092 !
4093 !---------------------------------------------------------------------
4094 !*** Loop through the four sides of the domain.
4095 !---------------------------------------------------------------------
4096 !
4097  sides: do nside=1,4
4098 !
4099  move=.false.
4100 !
4101  if(nside==1)then
4102  if(north_bc)then
4103  move=.true.
4104  bc_side_t0=>bc_t0%north
4105  bc_side_t1=>bc_t1%north
4106 !
4107  is_c=bnds%is_north
4108  ie_c=bnds%ie_north ! North BC index limits
4109  js_c=bnds%js_north ! for centers of grid cells.
4110  je_c=bnds%je_north
4111 !
4112  is_s=bnds%is_north_uvs
4113  ie_s=bnds%ie_north_uvs ! North BC index limits
4114  js_s=bnds%js_north_uvs ! for winds on N/S sides of grid cells.
4115  je_s=bnds%je_north_uvs
4116 !
4117  is_w=bnds%is_north_uvw
4118  ie_w=bnds%ie_north_uvw ! North BC index limits
4119  js_w=bnds%js_north_uvw ! for winds on E/W sides of grid cells.
4120  je_w=bnds%je_north_uvw
4121  endif
4122  endif
4123 !
4124  if(nside==2)then
4125  if(south_bc)then
4126  move=.true.
4127  bc_side_t0=>bc_t0%south
4128  bc_side_t1=>bc_t1%south
4129 !
4130  is_c=bnds%is_south
4131  ie_c=bnds%ie_south ! South BC index limits
4132  js_c=bnds%js_south ! for centers of grid cells.
4133  je_c=bnds%je_south
4134 !
4135  is_s=bnds%is_south_uvs
4136  ie_s=bnds%ie_south_uvs ! South BC index limits
4137  js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells.
4138  je_s=bnds%je_south_uvs
4139 !
4140  is_w=bnds%is_south_uvw
4141  ie_w=bnds%ie_south_uvw ! South BC index limits
4142  js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells.
4143  je_w=bnds%je_south_uvw
4144  endif
4145  endif
4146 !
4147  if(nside==3)then
4148  if(east_bc)then
4149  move=.true.
4150  bc_side_t0=>bc_t0%east
4151  bc_side_t1=>bc_t1%east
4152 !
4153  is_c=bnds%is_east
4154  ie_c=bnds%ie_east ! East BC index limits
4155  js_c=bnds%js_east ! for centers of grid cells.
4156  je_c=bnds%je_east
4157 !
4158  is_s=bnds%is_east_uvs
4159  ie_s=bnds%ie_east_uvs ! East BC index limits
4160  js_s=bnds%js_east_uvs ! for winds on N/S sides of grid cells.
4161  je_s=bnds%je_east_uvs
4162 !
4163  is_w=bnds%is_east_uvw
4164  ie_w=bnds%ie_east_uvw ! East BC index limits
4165  js_w=bnds%js_east_uvw ! for winds on E/W sides of grid cells.
4166  je_w=bnds%je_east_uvw
4167  endif
4168  endif
4169 !
4170  if(nside==4)then
4171  if(west_bc)then
4172  move=.true.
4173  bc_side_t0=>bc_t0%west
4174  bc_side_t1=>bc_t1%west
4175 !
4176  is_c=bnds%is_west
4177  ie_c=bnds%ie_west ! West BC index limits
4178  js_c=bnds%js_west ! for centers of grid cells.
4179  je_c=bnds%je_west
4180 !
4181  is_s=bnds%is_west_uvs
4182  ie_s=bnds%ie_west_uvs ! West BC index limits
4183  js_s=bnds%js_west_uvs ! for winds on N/S sides of grid cells.
4184  je_s=bnds%je_west_uvs
4185 !
4186  is_w=bnds%is_west_uvw
4187  ie_w=bnds%ie_west_uvw ! West BC index limits
4188  js_w=bnds%js_west_uvw ! for winds on E/W sides of grid cells.
4189  je_w=bnds%je_west_uvw
4190  endif
4191  endif
4192 !
4193  if(move)then
4194  do k=1,nlev
4195  do j=js_c,je_c
4196  do i=is_c,ie_c
4197  bc_side_t0%delp_BC(i,j,k) =bc_side_t1%delp_BC(i,j,k)
4198  bc_side_t0%divgd_BC(i,j,k)=bc_side_t1%divgd_BC(i,j,k)
4199  enddo
4200  enddo
4201  enddo
4202 !
4203  do n=1,ntracers
4204  do k=1,nlev
4205  do j=js_c,je_c
4206  do i=is_c,ie_c
4207  bc_side_t0%q_BC(i,j,k,n)=bc_side_t1%q_BC(i,j,k,n)
4208  enddo
4209  enddo
4210  enddo
4211  enddo
4212 !
4213  do k=1,nlev
4214  do j=js_c,je_c
4215  do i=is_c,ie_c
4216 #ifndef SW_DYNAMICS
4217  bc_side_t0%w_BC(i,j,k) =bc_side_t1%w_BC(i,j,k)
4218  bc_side_t0%pt_BC(i,j,k) =bc_side_t1%pt_BC(i,j,k)
4219  bc_side_t0%delz_BC(i,j,k) =bc_side_t1%delz_BC(i,j,k)
4220 #ifdef USE_COND
4221  bc_side_t0%q_con_BC(i,j,k)=bc_side_t1%q_con_BC(i,j,k)
4222 #ifdef MOIST_CAPPA
4223  bc_side_t0%cappa_BC(i,j,k)=bc_side_t1%cappa_BC(i,j,k)
4224 #endif
4225 #endif
4226 #endif
4227  enddo
4228  enddo
4229  enddo
4230 !
4231  do k=1,nlev
4232  do j=js_s,je_s
4233  do i=is_s,ie_s
4234  bc_side_t0%u_BC(i,j,k) =bc_side_t1%u_BC(i,j,k)
4235  bc_side_t0%vc_BC(i,j,k)=bc_side_t1%vc_BC(i,j,k)
4236  enddo
4237  enddo
4238  enddo
4239 !
4240  do k=1,nlev
4241  do j=js_w,je_w
4242  do i=is_w,ie_w
4243  bc_side_t0%v_BC(i,j,k) =bc_side_t1%v_BC(i,j,k)
4244  bc_side_t0%uc_BC(i,j,k)=bc_side_t1%uc_BC(i,j,k)
4245  enddo
4246  enddo
4247  enddo
4248 !
4249  endif
4250 !
4251  enddo sides
4252 !
4253 !---------------------------------------------------------------------
4254 !
4255  end subroutine regional_bc_t1_to_t0
4256 !
4257 !---------------------------------------------------------------------
4258 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4259 !---------------------------------------------------------------------
4260 !
4261  subroutine convert_to_virt_pot_temp(isd,ied,jsd,jed,npz &
4262  ,sphum,liq_wat )
4264 !-----------------------------------------------------------------------
4265 !*** Convert the incoming sensible temperature to virtual potential
4266 !*** temperature.
4267 !-----------------------------------------------------------------------
4268  implicit none
4269 !-----------------------------------------------------------------------
4270 !
4271 !---------------------
4272 !*** Input arguments
4273 !---------------------
4274 !
4275  integer,intent(in) :: isd,ied,jsd,jed,npz
4276 !
4277  integer,intent(in) :: liq_wat,sphum
4278 !
4279 !---------------------
4280 !*** Local variables
4281 !---------------------
4282 !
4283  integer :: i1,i2,j1,j2
4284 !
4285  real :: rdg
4286 !
4287  real,dimension(:,:,:),pointer :: delp,delz,pt
4288 #ifdef USE_COND
4289  real,dimension(:,:,:),pointer :: q_con
4290 #endif
4291 #ifdef MOIST_CAPPA
4292  real,dimension(:,:,:),pointer ::cappa
4293 #endif
4294 !
4295  real,dimension(:,:,:,:),pointer :: q
4296 !
4297 !-----------------------------------------------------------------------
4298 !***********************************************************************
4299 !-----------------------------------------------------------------------
4300 !
4301  if(.not.(north_bc.or.south_bc.or.east_bc.or.west_bc))then
4302  return
4303  endif
4304 !
4305  rdg=-rdgas/grav
4306 !
4307  if(north_bc)then
4308  i1=regional_bounds%is_north
4309  i2=regional_bounds%ie_north
4310  j1=regional_bounds%js_north
4311  j2=regional_bounds%je_north
4312  q =>bc_t1%north%q_BC
4313 #ifdef USE_COND
4314  q_con=>bc_t1%north%q_con_BC
4315 #endif
4316  delp =>bc_t1%north%delp_BC
4317  delz =>bc_t1%north%delz_BC
4318 #ifdef MOIST_CAPPA
4319  cappa=>bc_t1%north%cappa_BC
4320 #endif
4321  pt =>bc_t1%north%pt_BC
4322  call compute_vpt
4323  endif
4324 !
4325  if(south_bc)then
4326  i1=regional_bounds%is_south
4327  i2=regional_bounds%ie_south
4328  j1=regional_bounds%js_south
4329  j2=regional_bounds%je_south
4330  q =>bc_t1%south%q_BC
4331 #ifdef USE_COND
4332  q_con=>bc_t1%south%q_con_BC
4333 #endif
4334  delp =>bc_t1%south%delp_BC
4335  delz =>bc_t1%south%delz_BC
4336 #ifdef MOIST_CAPPA
4337  cappa=>bc_t1%south%cappa_BC
4338 #endif
4339  pt =>bc_t1%south%pt_BC
4340  call compute_vpt
4341  endif
4342 !
4343  if(east_bc)then
4344  i1=regional_bounds%is_east
4345  i2=regional_bounds%ie_east
4346  j1=regional_bounds%js_east
4347  j2=regional_bounds%je_east
4348  q =>bc_t1%east%q_BC
4349 #ifdef USE_COND
4350  q_con=>bc_t1%east%q_con_BC
4351 #endif
4352  delp =>bc_t1%east%delp_BC
4353  delz =>bc_t1%east%delz_BC
4354 #ifdef MOIST_CAPPA
4355  cappa=>bc_t1%east%cappa_BC
4356 #endif
4357  pt =>bc_t1%east%pt_BC
4358  call compute_vpt
4359  endif
4360 !
4361  if(west_bc)then
4362  i1=regional_bounds%is_west
4363  i2=regional_bounds%ie_west
4364  j1=regional_bounds%js_west
4365  j2=regional_bounds%je_west
4366  q =>bc_t1%west%q_BC
4367 #ifdef USE_COND
4368  q_con=>bc_t1%west%q_con_BC
4369 #endif
4370  delp =>bc_t1%west%delp_BC
4371  delz =>bc_t1%west%delz_BC
4372 #ifdef MOIST_CAPPA
4373  cappa=>bc_t1%west%cappa_BC
4374 #endif
4375  pt =>bc_t1%west%pt_BC
4376  call compute_vpt
4377  endif
4378 !
4379 !-----------------------------------------------------------------------
4380 
4381  contains
4382 
4383 !-----------------------------------------------------------------------
4384 !
4385  subroutine compute_vpt
4387 !-----------------------------------------------------------------------
4388 !*** Compute the virtual potential temperature as done in fv_dynamics.
4389 !-----------------------------------------------------------------------
4390 !
4391 !---------------------
4392 !*** Local variables
4393 !---------------------
4394 !
4395  integer :: i,j,k
4396 !
4397  real :: cvm,dp1,pkz
4398 !
4399 !-----------------------------------------------------------------------
4400 !***********************************************************************
4401 !-----------------------------------------------------------------------
4402 !
4403  do k=1,npz
4404 !
4405  do j=j1,j2
4406  do i=i1,i2
4407  dp1 = zvir*q(i,j,k,sphum)
4408 #ifdef USE_COND
4409 #ifdef MOIST_CAPPA
4410  cvm=(1.-q(i,j,k,sphum)+q_con(i,j,k))*cv_air &
4411  +q(i,j,k,sphum)*cv_vap+q(i,j,k,liq_wat)*c_liq
4412  pkz=exp(cappa(i,j,k)*log(rdg*delp(i,j,k)*pt(i,j,k) &
4413  *(1.+dp1)*(1.-q_con(i,j,k))/delz(i,j,k)))
4414 #else
4415  pkz=exp(kappa*log(rdg*delp(i,j,k)*pt(i,j,k) &
4416  *(1.+dp1)*(1.-q_con(i,j,k))/delz(i,j,k)))
4417 #endif
4418  pt(i,j,k)=pt(i,j,k)*(1.+dp1)*(1.-q_con(i,j,k))/pkz
4419 #else
4420  pkz=exp(kappa*log(rdg*delp(i,j,k)*pt(i,j,k) &
4421  *(1.+dp1)/delz(i,j,k)))
4422  pt(i,j,k)=pt(i,j,k)*(1.+dp1)/pkz
4423 #endif
4424  enddo
4425  enddo
4426 !
4427  enddo
4428 !
4429 !-----------------------------------------------------------------------
4430 !
4431  end subroutine compute_vpt
4432 !
4433 !-----------------------------------------------------------------------
4434 !
4435  end subroutine convert_to_virt_pot_temp
4436 !
4437 !-----------------------------------------------------------------------
4438 !
4439 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4440 !---------------------------------------------------------------------
4441 !*** The following four subroutines are exact copies from
4442 !*** external_ic_mod. That module must USE this module therefore
4443 !*** this module cannout USE external_IC_mod to get at those
4444 !*** subroutines. The routines may be moved to their own module.
4445 !---------------------------------------------------------------------
4446 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4447 !---------------------------------------------------------------------
4448  subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
4449  character(len=*), intent(in):: qname
4450  integer, intent(in):: is, ie, js, je, km
4451  real, intent(in):: q(is:ie, js:je, km)
4452  real, intent(in):: fac
4453  real qmin, qmax
4454  integer i,j,k
4455 
4456  qmin = q(is,js,1)
4457  qmax = qmin
4458  do k=1,km
4459  do j=js,je
4460  do i=is,ie
4461  if( q(i,j,k) < qmin ) then
4462  qmin = q(i,j,k)
4463  elseif( q(i,j,k) > qmax ) then
4464  qmax = q(i,j,k)
4465  endif
4466  enddo
4467  enddo
4468  enddo
4469  call mp_reduce_min(qmin)
4470  call mp_reduce_max(qmax)
4471  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac
4472 
4473  end subroutine p_maxmin
4474 
4475 
4476  subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
4477  character(len=*), intent(in):: qname
4478  integer, intent(in):: is, ie, js, je
4479  integer, intent(in):: km
4480  real, intent(in):: q(is:ie, js:je, km)
4481  real, intent(in):: fac
4482  real(kind=R_GRID), intent(IN):: area(is-3:ie+3, js-3:je+3)
4483  type(domain2d), intent(INOUT) :: domain
4484 !---local variables
4485  real qmin, qmax, gmean
4486  integer i,j,k
4487 
4488  qmin = q(is,js,1)
4489  qmax = qmin
4490  gmean = 0.
4491 
4492  do k=1,km
4493  do j=js,je
4494  do i=is,ie
4495  if( q(i,j,k) < qmin ) then
4496  qmin = q(i,j,k)
4497  elseif( q(i,j,k) > qmax ) then
4498  qmax = q(i,j,k)
4499  endif
4500  enddo
4501  enddo
4502  enddo
4503 
4504  call mp_reduce_min(qmin)
4505  call mp_reduce_max(qmax)
4506 
4507  gmean = g_sum(domain, q(is,js,km), is, ie, js, je, 3, area, 1, reproduce=.true.)
4508  if(is_master()) write(6,*) qname, qmax*fac, qmin*fac, gmean*fac
4509 
4510  end subroutine pmaxmn
4511 
4512 
4513  subroutine fillq(im, km, nq, q, dp)
4514  integer, intent(in):: im ! No. of longitudes
4515  integer, intent(in):: km ! No. of levels
4516  integer, intent(in):: nq ! Total number of tracers
4517  real , intent(in):: dp(im,km) ! pressure thickness
4518  real , intent(inout) :: q(im,km,nq) ! tracer mixing ratio
4519 ! !LOCAL VARIABLES:
4520  integer i, k, ic, k1
4521 
4522  do ic=1,nq
4523 ! Bottom up:
4524  do k=km,2,-1
4525  k1 = k-1
4526  do i=1,im
4527  if( q(i,k,ic) < 0. ) then
4528  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
4529  q(i,k ,ic) = 0.
4530  endif
4531  enddo
4532  enddo
4533 ! Top down:
4534  do k=1,km-1
4535  k1 = k+1
4536  do i=1,im
4537  if( q(i,k,ic) < 0. ) then
4538  q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
4539  q(i,k ,ic) = 0.
4540  endif
4541  enddo
4542  enddo
4543 
4544  enddo
4545 
4546  end subroutine fillq
4547 
4548  subroutine mp_auto_conversion(ql, qr, qi, qs)
4549  real, intent(inout):: ql, qr, qi, qs
4550  real, parameter:: qi0_max = 2.0e-3
4551  real, parameter:: ql0_max = 2.5e-3
4552 
4553 ! Convert excess cloud water into rain:
4554  if ( ql > ql0_max ) then
4555  qr = ql - ql0_max
4556  ql = ql0_max
4557  endif
4558 ! Convert excess cloud ice into snow:
4559  if ( qi > qi0_max ) then
4560  qs = qi - qi0_max
4561  qi = qi0_max
4562  endif
4563 
4564  end subroutine mp_auto_conversion
4565 
4566 !-----------------------------------------------------------------------
4567 !
4568  subroutine nudge_qv_bc(Atm,isd,ied,jsd,jed)
4570 !-----------------------------------------------------------------------
4571 !*** When nudging of specific humidity is selected then we must also
4572 !*** nudge the values in the regional boundary.
4573 !-----------------------------------------------------------------------
4574  implicit none
4575 !-----------------------------------------------------------------------
4576 !
4577 !------------------------
4578 !*** Argument variables
4579 !------------------------
4580 !
4581  integer,intent(in) :: isd,ied,jsd,jed
4582 !
4583  type(fv_atmos_type),target,intent(inout) :: Atm
4584 !
4585 !---------------------
4586 !*** Local variables
4587 !---------------------
4588 !
4589  integer :: i,i_x,ie,is,j,j_x,je,js,k
4590 !
4591  real, parameter:: q1_h2o = 2.2e-6
4592  real, parameter:: q7_h2o = 3.8e-6
4593  real, parameter:: q100_h2o = 3.8e-6
4594  real, parameter:: q1000_h2o = 3.1e-6
4595  real, parameter:: q2000_h2o = 2.8e-6
4596  real, parameter:: q3000_h2o = 3.0e-6
4597  real, parameter:: wt=2., xt=1./(1.+wt)
4598 !
4599  real :: p00,q00
4600 !
4601  type(fv_regional_bc_bounds_type),pointer :: bnds
4602 !
4603 !-----------------------------------------------------------------------
4604 !***********************************************************************
4605 !-----------------------------------------------------------------------
4606 !
4607  bnds=>atm%regional_bc_bounds
4608 !
4609 !-----------
4610 !*** North
4611 !-----------
4612 !
4613  if(north_bc)then
4614  is=lbound(bc_t1%north%q_BC,1)
4615  ie=ubound(bc_t1%north%q_BC,1)
4616  js=lbound(bc_t1%north%q_BC,2)
4617  je=ubound(bc_t1%north%q_BC,2)
4618 !
4619  i_x=isd
4620  j_x=jsd ! this location.
4621 !
4622  p00=atm%ptop
4623 !
4624  n_loopk: do k=1,npz
4625  if(p00<3000.)then
4626  call get_q00
4627  do j=js,je
4628  do i=is,ie
4629  bc_t1%north%q_BC(i,j,k,sphum_index)= &
4630  xt*(bc_t1%north%q_BC(i,j,k,sphum_index)+wt*q00)
4631  enddo
4632  enddo
4633  else
4634  exit n_loopk
4635  endif
4636  p00=p00+bc_t1%north%delp_BC(i_x,j_x,k)
4637  enddo n_loopk
4638  endif
4639 !
4640 !-----------
4641 !*** South
4642 !-----------
4643 !
4644  if(south_bc)then
4645  is=lbound(bc_t1%south%q_BC,1)
4646  ie=ubound(bc_t1%south%q_BC,1)
4647  js=lbound(bc_t1%south%q_BC,2)
4648  je=ubound(bc_t1%south%q_BC,2)
4649 !
4650  i_x=isd
4651  j_x=jed ! this location.
4652 !
4653  p00=atm%ptop
4654 !
4655  s_loopk: do k=1,npz
4656  if(p00<3000.)then
4657  call get_q00
4658  do j=js,je
4659  do i=is,ie
4660  bc_t1%south%q_BC(i,j,k,sphum_index)= &
4661  xt*(bc_t1%south%q_BC(i,j,k,sphum_index)+wt*q00)
4662  enddo
4663  enddo
4664  else
4665  exit s_loopk
4666  endif
4667  p00=p00+bc_t1%south%delp_BC(i_x,j_x,k)
4668  enddo s_loopk
4669  endif
4670 !
4671 !----------
4672 !*** East
4673 !----------
4674 !
4675  if(east_bc)then
4676  is=lbound(bc_t1%east%q_BC,1)
4677  ie=ubound(bc_t1%east%q_BC,1)
4678  js=lbound(bc_t1%east%q_BC,2)
4679  je=ubound(bc_t1%east%q_BC,2)
4680 !
4681  i_x=isd
4682  j_x=jsd+nhalo_model ! this location.
4683 !
4684  p00=atm%ptop
4685 !
4686  e_loopk: do k=1,npz
4687  if(p00<3000.)then
4688  call get_q00
4689  do j=js,je
4690  do i=is,ie
4691  bc_t1%east%q_BC(i,j,k,sphum_index)= &
4692  xt*(bc_t1%east%q_BC(i,j,k,sphum_index)+wt*q00)
4693  enddo
4694  enddo
4695  else
4696  exit e_loopk
4697  endif
4698  p00=p00+bc_t1%east%delp_BC(i_x,j_x,k)
4699  enddo e_loopk
4700  endif
4701 !
4702 !----------
4703 !*** West
4704 !----------
4705 !
4706  if(west_bc)then
4707  is=lbound(bc_t1%west%q_BC,1)
4708  ie=ubound(bc_t1%west%q_BC,1)
4709  js=lbound(bc_t1%west%q_BC,2)
4710  je=ubound(bc_t1%west%q_BC,2)
4711 !
4712  i_x=ied
4713  j_x=jsd+nhalo_model ! this location.
4714 !
4715  p00=atm%ptop
4716 !
4717  w_loopk: do k=1,npz
4718  if(p00<3000.)then
4719  call get_q00
4720  do j=js,je
4721  do i=is,ie
4722  bc_t1%west%q_BC(i,j,k,sphum_index)= &
4723  xt*(bc_t1%west%q_BC(i,j,k,sphum_index)+wt*q00)
4724  enddo
4725  enddo
4726  else
4727  exit w_loopk
4728  endif
4729  p00=p00+bc_t1%west%delp_BC(i_x,j_x,k)
4730  enddo w_loopk
4731  endif
4732 !
4733 !-----------------------------------------------------------------------
4734 !
4735  contains
4736 !
4737 !-----------------------------------------------------------------------
4738 !
4739  subroutine get_q00
4741 !-----------------------------------------------------------------------
4742 !*** This is an internal subroutine to subroutine nudge_qv_bc that
4743 !*** computes the climatological contribution to the nudging ot the
4744 !*** input specific humidity.
4745 !-----------------------------------------------------------------------
4746 !***********************************************************************
4747 !-----------------------------------------------------------------------
4748 !
4749  if ( p00 < 30.e2 ) then
4750  if ( p00 < 1. ) then
4751  q00 = q1_h2o
4752  elseif ( p00 <= 7. .and. p00 >= 1. ) then
4753  q00 = q1_h2o + (q7_h2o-q1_h2o)*log(pref(k)/1.)/log(7.)
4754  elseif ( p00 < 100. .and. p00 >= 7. ) then
4755  q00 = q7_h2o + (q100_h2o-q7_h2o)*log(pref(k)/7.)/log(100./7.)
4756  elseif ( p00 < 1000. .and. p00 >= 100. ) then
4757  q00 = q100_h2o + (q1000_h2o-q100_h2o)*log(pref(k)/1.e2)/log(10.)
4758  elseif ( p00 < 2000. .and. p00 >= 1000. ) then
4759  q00 = q1000_h2o + (q2000_h2o-q1000_h2o)*log(pref(k)/1.e3)/log(2.)
4760  else
4761  q00 = q2000_h2o + (q3000_h2o-q2000_h2o)*log(pref(k)/2.e3)/log(1.5)
4762  endif
4763  endif
4764 !
4765 !-----------------------------------------------------------------------
4766 !
4767  end subroutine get_q00
4768 !
4769 !-----------------------------------------------------------------------
4770 !
4771  end subroutine nudge_qv_bc
4772 !
4773 !-----------------------------------------------------------------------
4774 !-----------------------------------------------------------------------
4775 
4776  subroutine dump_field_3d (domain, name, field, isd, ied, jsd, jed, nlev, stag)
4778 ! call dump_field(Atm(1)%domain,"atm_pt", Atm(1)%pt, isd, ied, jsd, jed, Atm(1)%npz, stag=H_STAGGER)
4779 ! call dump_field(Atm(1)%domain,"atm_u", Atm(1)%u, isd, ied, jsd, jed+1, Atm(1)%npz, stag=U_STAGGER)
4780 ! call dump_field(Atm(1)%domain,"atm_v", Atm(1)%v, isd, ied+1, jsd, jed, Atm(1)%npz, stag=V_STAGGER)
4781 ! call dump_field(Atm(1)%domain,"atm_phis", Atm(1)%phis, isd, ied, jsd, jed, stag=H_STAGGER)
4782 
4783  type(domain2d), intent(INOUT) :: domain
4784  character(len=*), intent(IN) :: name
4785  real, dimension(isd:ied,jsd:jed,1:nlev), intent(INOUT) :: field
4786  integer, intent(IN) :: isd, ied, jsd, jed, nlev
4787  integer, intent(IN) :: stag
4788 
4789  integer :: unit
4790  character(len=128) :: fname
4791  type(axistype) :: x, y, z
4792  type(fieldtype) :: f
4793  type(domain1D) :: xdom, ydom
4794  integer :: nz
4795  integer :: is, ie, js, je
4796  integer :: isg, ieg, jsg, jeg, nxg, nyg, npx, npy
4797  integer :: i, j, halo, iext, jext
4798  logical :: is_root_pe
4799  real, allocatable, dimension(:,:,:) :: glob_field
4800  integer, allocatable, dimension(:) :: pelist
4801  character(len=1) :: stagname
4802  integer :: isection_s, isection_e, jsection_s, jsection_e
4803 
4804  write(fname,"(A,A,A,I1.1,A)") "regional_",name,".tile", 7 , ".nc"
4805  write(0,*)'dump_field_3d: file name = |', trim(fname) , '|'
4806 
4807  call mpp_get_domain_components( domain, xdom, ydom )
4808  call mpp_get_compute_domain( domain, is, ie, js, je )
4809  call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=npx, ysize=npy, position=center )
4810 
4811  halo = is - isd
4812  if ( halo /= 3 ) then
4813  write(0,*) 'dusan- halo should be 3 ', halo
4814  endif
4815 
4816  iext = 0
4817  jext = 0
4818  stagname = "h";
4819  if (stag == u_stagger) then
4820  jext = 1
4821  stagname = "u";
4822  endif
4823  if (stag == v_stagger) then
4824  iext = 1
4825  stagname = "v";
4826  endif
4827 
4828  nxg = npx + 2*halo + iext
4829  nyg = npy + 2*halo + jext
4830  nz = size(field,dim=3)
4831 
4832  allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) )
4833 
4834  isection_s = is
4835  isection_e = ie
4836  jsection_s = js
4837  jsection_e = je
4838 
4839  if ( isd < 0 ) isection_s = isd
4840  if ( ied > npx-1 ) isection_e = ied
4841  if ( jsd < 0 ) jsection_s = jsd
4842  if ( jed > npy-1 ) jsection_e = jed
4843 
4844  allocate( pelist(mpp_npes()) )
4845  call mpp_get_current_pelist(pelist)
4846 
4847  is_root_pe = (mpp_pe()==mpp_root_pe())
4848 
4849  call mpp_gather(isection_s,isection_e,jsection_s,jsection_e, nz, &
4850  pelist, field(isection_s:isection_e,jsection_s:jsection_e,:), glob_field, is_root_pe, halo, halo)
4851 
4852  call mpp_open( unit, trim(fname), action=mpp_wronly, form=mpp_netcdf, threading=mpp_single)
4853 
4854  call mpp_write_meta( unit, x, 'grid_xt', 'km', 'X distance', 'X', domain=xdom, data=(/(i*1.0,i=1,nxg)/) )
4855  call mpp_write_meta( unit, y, 'grid_yt', 'km', 'Y distance', 'Y', domain=ydom, data=(/(j*1.0,j=1,nyg)/) )
4856  call mpp_write_meta( unit, z, 'lev', 'km', 'Z distance', data=(/(i*1.0,i=1,nz)/) )
4857 
4858  call mpp_write_meta( unit, f, (/x,y,z/), name, 'unit', name)
4859  call mpp_write_meta( unit, "stretch_factor", rval=stretch_factor )
4860  call mpp_write_meta( unit, "target_lon", rval=target_lon )
4861  call mpp_write_meta( unit, "target_lat", rval=target_lat )
4862  call mpp_write_meta( unit, "cube_res", ival= cube_res)
4863  call mpp_write_meta( unit, "parent_tile", ival=parent_tile )
4864  call mpp_write_meta( unit, "refine_ratio", ival=refine_ratio )
4865  call mpp_write_meta( unit, "istart_nest", ival=istart_nest )
4866  call mpp_write_meta( unit, "jstart_nest", ival=jstart_nest )
4867  call mpp_write_meta( unit, "iend_nest", ival=iend_nest )
4868  call mpp_write_meta( unit, "jend_nest", ival=jend_nest )
4869  call mpp_write_meta( unit, "ihalo_shift", ival=halo )
4870  call mpp_write_meta( unit, "jhalo_shift", ival=halo )
4871  call mpp_write_meta( unit, mpp_get_id(f), "hstagger", cval=stagname )
4872  call mpp_write( unit, x )
4873  call mpp_write( unit, y )
4874  call mpp_write( unit, z )
4875  call mpp_write( unit, f, glob_field )
4876 
4877  call mpp_close( unit )
4878 
4879  end subroutine dump_field_3d
4880 
4881  subroutine dump_field_2d (domain, name, field, isd, ied, jsd, jed, stag)
4883  type(domain2d), intent(INOUT) :: domain
4884  character(len=*), intent(IN) :: name
4885  real, dimension(isd:ied,jsd:jed), intent(INOUT) :: field
4886  integer, intent(IN) :: isd, ied, jsd, jed
4887  integer, intent(IN) :: stag
4888 
4889  integer :: unit
4890  character(len=128) :: fname
4891  type(axistype) :: x, y
4892  type(fieldtype) :: f
4893  type(domain1D) :: xdom, ydom
4894  integer :: is, ie, js, je
4895  integer :: isg, ieg, jsg, jeg, nxg, nyg, npx, npy
4896  integer :: i, j, halo, iext, jext
4897  logical :: is_root_pe
4898  real, allocatable, dimension(:,:) :: glob_field
4899  integer, allocatable, dimension(:) :: pelist
4900  character(len=1) :: stagname
4901  integer :: isection_s, isection_e, jsection_s, jsection_e
4902 
4903  write(fname,"(A,A,A,I1.1,A)") "regional_",name,".tile", 7 , ".nc"
4904 ! write(0,*)'dump_field_3d: file name = |', trim(fname) , '|'
4905 
4906  call mpp_get_domain_components( domain, xdom, ydom )
4907  call mpp_get_compute_domain( domain, is, ie, js, je )
4908  call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=npx, ysize=npy, position=center )
4909 
4910  halo = is - isd
4911  if ( halo /= 3 ) then
4912  write(0,*) 'dusan- halo should be 3 ', halo
4913  endif
4914 
4915  iext = 0
4916  jext = 0
4917  stagname = "h";
4918  if (stag == u_stagger) then
4919  jext = 1
4920  stagname = "u";
4921  endif
4922  if (stag == v_stagger) then
4923  iext = 1
4924  stagname = "v";
4925  endif
4926 
4927  nxg = npx + 2*halo + iext
4928  nyg = npy + 2*halo + jext
4929 
4930  allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) )
4931 
4932  isection_s = is
4933  isection_e = ie
4934  jsection_s = js
4935  jsection_e = je
4936 
4937  if ( isd < 0 ) isection_s = isd
4938  if ( ied > npx-1 ) isection_e = ied
4939  if ( jsd < 0 ) jsection_s = jsd
4940  if ( jed > npy-1 ) jsection_e = jed
4941 
4942  allocate( pelist(mpp_npes()) )
4943  call mpp_get_current_pelist(pelist)
4944 
4945  is_root_pe = (mpp_pe()==mpp_root_pe())
4946 
4947  call mpp_gather(isection_s,isection_e,jsection_s,jsection_e, &
4948  pelist, field(isection_s:isection_e,jsection_s:jsection_e), glob_field, is_root_pe, halo, halo)
4949 
4950  call mpp_open( unit, trim(fname), action=mpp_wronly, form=mpp_netcdf, threading=mpp_single)
4951 
4952  call mpp_write_meta( unit, x, 'grid_xt', 'km', 'X distance', 'X', domain=xdom, data=(/(i*1.0,i=1,nxg)/) )
4953  call mpp_write_meta( unit, y, 'grid_yt', 'km', 'Y distance', 'Y', domain=ydom, data=(/(j*1.0,j=1,nyg)/) )
4954 
4955  call mpp_write_meta( unit, f, (/x,y/), name, 'unit', name)
4956  call mpp_write_meta( unit, "stretch_factor", rval=stretch_factor )
4957  call mpp_write_meta( unit, "target_lon", rval=target_lon )
4958  call mpp_write_meta( unit, "target_lat", rval=target_lat )
4959  call mpp_write_meta( unit, "cube_res", ival= cube_res)
4960  call mpp_write_meta( unit, "parent_tile", ival=parent_tile )
4961  call mpp_write_meta( unit, "refine_ratio", ival=refine_ratio )
4962  call mpp_write_meta( unit, "istart_nest", ival=istart_nest )
4963  call mpp_write_meta( unit, "jstart_nest", ival=jstart_nest )
4964  call mpp_write_meta( unit, "iend_nest", ival=iend_nest )
4965  call mpp_write_meta( unit, "jend_nest", ival=jend_nest )
4966  call mpp_write_meta( unit, "ihalo_shift", ival=halo )
4967  call mpp_write_meta( unit, "jhalo_shift", ival=halo )
4968  call mpp_write_meta( unit, mpp_get_id(f), "hstagger", cval=stagname )
4969  call mpp_write( unit, x )
4970  call mpp_write( unit, y )
4971  call mpp_write( unit, f, glob_field )
4972 
4973  call mpp_close( unit )
4974 
4975  end subroutine dump_field_2d
4976 
4977 !---------------------------------------------------------------------
4978 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
4979 !---------------------------------------------------------------------
4980 
4981  subroutine exch_uv(domain, bd, npz, u, v)
4982  use mpi
4983 
4984  implicit none
4985 
4986  type(domain2d), intent(inout) :: domain
4987  type(fv_grid_bounds_type), intent(in) :: bd
4988  integer, intent(in) :: npz
4989  real, intent(inout) :: u (bd%isd:bd%ied ,bd%jsd:bd%jed+1,1:npz)
4990  real, intent(inout) :: v (bd%isd:bd%ied+1,bd%jsd:bd%jed ,1:npz)
4991 
4992  integer,parameter :: ibufexch=2500000
4993  real,dimension(ibufexch) :: buf1,buf2,buf3,buf4
4994  integer :: ihandle1,ihandle2,ihandle3,ihandle4
4995  integer,dimension(MPI_STATUS_SIZE) :: istat
4996  integer :: ic, i, j, k, is, ie, js, je
4997  integer :: irecv, isend, ierr
4998 
4999  integer :: mype
5000  integer :: north_pe, south_pe, east_pe, west_pe
5001 
5002  mype = mpp_pe()
5003  call mpp_get_neighbor_pe( domain, north, north_pe)
5004  call mpp_get_neighbor_pe( domain, south, south_pe)
5005  call mpp_get_neighbor_pe( domain, west, west_pe)
5006  call mpp_get_neighbor_pe( domain, east, east_pe)
5007 
5008  ! write(0,*) ' north_pe = ', north_pe
5009  ! write(0,*) ' south_pe = ', south_pe
5010  ! write(0,*) ' west_pe = ', west_pe
5011  ! write(0,*) ' east_pe = ', east_pe
5012 
5013  is=bd%is
5014  ie=bd%ie
5015  js=bd%js
5016  je=bd%je
5017 
5018 ! FIXME: MPI_COMM_WORLD
5019 
5020 
5021 ! Receive from north
5022  if( north_pe /= null_pe )then
5023  call mpi_irecv(buf1,ibufexch,mpi_real,north_pe,north_pe &
5024  ,mpi_comm_world,ihandle1,irecv)
5025  endif
5026 
5027 ! Receive from south
5028  if( south_pe /= null_pe )then
5029  call mpi_irecv(buf2,ibufexch,mpi_real,south_pe,south_pe &
5030  ,mpi_comm_world,ihandle2,irecv)
5031  endif
5032 
5033 ! Send to north
5034  if( north_pe /= null_pe )then
5035  ic=0
5036  do k=1,npz
5037 
5038  do j=je-3+1,je-1+1
5039  do i=is-3,is-1
5040  ic=ic+1
5041  buf3(ic)=u(i,j,k)
5042  enddo
5043  do i=ie+1,ie+3
5044  ic=ic+1
5045  buf3(ic)=u(i,j,k)
5046  enddo
5047  enddo
5048 
5049  do j=je-2,je
5050  do i=is-3,is-1
5051  ic=ic+1
5052  buf3(ic)=v(i,j,k)
5053  enddo
5054  do i=ie+1,ie+3
5055  ic=ic+1
5056  buf3(ic)=v(i,j,k)
5057  enddo
5058  enddo
5059 
5060  enddo
5061  call mpi_issend(buf3,ic,mpi_real,north_pe,mype &
5062  ,mpi_comm_world,ihandle3,isend)
5063  endif
5064 
5065 ! Send to south
5066  if( south_pe /= null_pe )then
5067  ic=0
5068  do k=1,npz
5069 
5070  do j=js+2,js+3
5071  do i=is-3,is-1
5072  ic=ic+1
5073  buf4(ic)=u(i,j,k)
5074  enddo
5075  do i=ie+1,ie+3
5076  ic=ic+1
5077  buf4(ic)=u(i,j,k)
5078  enddo
5079  enddo
5080 
5081  do j=js+1,js+2
5082  do i=is-3,is-1
5083  ic=ic+1
5084  buf4(ic)=v(i,j,k)
5085  enddo
5086  do i=ie+1,ie+3
5087  ic=ic+1
5088  buf4(ic)=v(i,j,k)
5089  enddo
5090  enddo
5091 
5092  enddo
5093  call mpi_issend(buf4,ic,mpi_real,south_pe,mype &
5094  ,mpi_comm_world,ihandle4,isend)
5095  endif
5096 
5097 ! Store from south
5098  if( south_pe /= null_pe )then
5099  ic=0
5100  call mpi_wait(ihandle2,istat,ierr)
5101  do k=1,npz
5102 
5103  do j=js-3,js-1
5104  do i=is-3,is-1
5105  ic=ic+1
5106  u(i,j,k)=buf2(ic)
5107  enddo
5108  do i=ie+1,ie+3
5109  ic=ic+1
5110  u(i,j,k)=buf2(ic)
5111  enddo
5112  enddo
5113 
5114  do j=js-3,js-1
5115  do i=is-3,is-1
5116  ic=ic+1
5117  v(i,j,k)=buf2(ic)
5118  enddo
5119  do i=ie+1,ie+3
5120  ic=ic+1
5121  v(i,j,k)=buf2(ic)
5122  enddo
5123  enddo
5124 
5125  enddo
5126  endif
5127 
5128 ! Store from north
5129  if( north_pe /= null_pe )then
5130  ic=0
5131  call mpi_wait(ihandle1,istat,ierr)
5132  do k=1,npz
5133 
5134  do j=je+2+1,je+3+1
5135  do i=is-3,is-1
5136  ic=ic+1
5137  u(i,j,k)=buf1(ic)
5138  enddo
5139  do i=ie+1,ie+3
5140  ic=ic+1
5141  u(i,j,k)=buf1(ic)
5142  enddo
5143  enddo
5144 
5145  do j=je+2,je+3
5146  do i=is-3,is-1
5147  ic=ic+1
5148  v(i,j,k)=buf1(ic)
5149  enddo
5150  do i=ie+1,ie+3
5151  ic=ic+1
5152  v(i,j,k)=buf1(ic)
5153  enddo
5154  enddo
5155 
5156  enddo
5157  endif
5158 
5159  end subroutine exch_uv
5160 
5161  subroutine get_data_source(source,regional)
5163 ! This routine extracts the data source information if it is present in the datafile.
5164 !
5165  character (len = 80) :: source
5166  integer :: ncids,sourcelength
5167  logical :: lstatus,regional
5168 !
5169 ! Use the fms call here so we can actually get the return code value.
5170 !
5171  if (regional) then
5172  lstatus = get_global_att_value('INPUT/gfs_data.nc',"source", source)
5173  else
5174  lstatus = get_global_att_value('INPUT/gfs_data.tile1.nc',"source", source)
5175  endif
5176  if (.not. lstatus) then
5177  if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute'
5178  source='No Source Attribute'
5179  endif
5180  end subroutine get_data_source
5181 
5182  subroutine set_delp_and_tracers(BC_side,npz,nwat)
5184 ! This routine mimics what is done in external_ic to add mass back to delp
5185 ! and remove it from the tracers
5186 !
5187  integer :: npz,nwat
5188  type(fv_regional_BC_variables),intent(inout) :: BC_side
5189 !
5190 ! local variables
5191 !
5192  integer :: k, j, i, iq, is, ie, js, je
5193  integer :: liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
5194  real :: qt, wt, m_fac
5195 
5196  is=lbound(bc_side%delp_BC,1)
5197  ie=ubound(bc_side%delp_BC,1)
5198  js=lbound(bc_side%delp_BC,2)
5199  je=ubound(bc_side%delp_BC,2)
5200 !
5201  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
5202  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
5203  rainwat = get_tracer_index(model_atmos, 'rainwat')
5204  snowwat = get_tracer_index(model_atmos, 'snowwat')
5205  graupel = get_tracer_index(model_atmos, 'graupel')
5206  cld_amt = get_tracer_index(model_atmos, 'cld_amt')
5207 !
5208  source: if (trim(data_source) == 'FV3GFS GAUSSIAN NEMSIO FILE') then
5209 !
5210 ! if (cld_amt > 0) BC_side%q_BC(:,:,:,cld_amt) = 0.0 ! Moorthi
5211  do k=1,npz
5212  do j=js,je
5213  do i=is,ie
5214  wt = bc_side%delp_BC(i,j,k)
5215  if ( nwat == 6 ) then
5216  qt = wt*(1. + bc_side%q_BC(i,j,k,liq_wat) + &
5217  bc_side%q_BC(i,j,k,ice_wat) + &
5218  bc_side%q_BC(i,j,k,rainwat) + &
5219  bc_side%q_BC(i,j,k,snowwat) + &
5220  bc_side%q_BC(i,j,k,graupel))
5221  else ! all other values of nwat
5222  qt = wt*(1. + sum(bc_side%q_BC(i,j,k,2:nwat)))
5223  endif
5224 !--- Adjust delp with tracer mass.
5225  bc_side%delp_BC(i,j,k) = qt
5226  enddo
5227  enddo
5228  enddo
5229 !
5230  else source ! This else block is for all sources other than FV3GFS GAUSSIAN NEMSIO FILE
5231 !
5232 ! 20160928: Adjust the mixing ratios consistently...
5233  do k=1,npz
5234  do j=js,je
5235  do i=is,ie
5236  wt = bc_side%delp_BC(i,j,k)
5237  if ( nwat == 6 ) then
5238  qt = wt*(1. + bc_side%q_BC(i,j,k,liq_wat) + &
5239  bc_side%q_BC(i,j,k,ice_wat) + &
5240  bc_side%q_BC(i,j,k,rainwat) + &
5241  bc_side%q_BC(i,j,k,snowwat) + &
5242  bc_side%q_BC(i,j,k,graupel))
5243  else ! all other values of nwat
5244  qt = wt*(1. + sum(bc_side%q_BC(i,j,k,2:nwat)))
5245  endif
5246  m_fac = wt / qt
5247  do iq=1,ntracers
5248  bc_side%q_BC(i,j,k,iq) = m_fac * bc_side%q_BC(i,j,k,iq)
5249  enddo
5250  bc_side%delp_BC(i,j,k) = qt
5251  enddo
5252  enddo
5253  enddo
5254 !
5255  endif source
5256 !
5257  end subroutine set_delp_and_tracers
5258 
5259 !---------------------------------------------------------------------
5260 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
5261 !---------------------------------------------------------------------
5262 
5263 end module fv_regional_mod
5264 
5265 !---------------------------------------------------------------------
subroutine dump_field_2d(domain, name, field, isd, ied, jsd, jed, stag)
integer, save npz
integer, parameter jend_nest
integer, parameter iend_nest
subroutine, public mid_pt_sphere(p1, p2, pm)
subroutine, public prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
integer, public k_step
real, parameter c_ice
subroutine dump_field_3d(domain, name, field, isd, ied, jsd, jed, nlev, stag)
logical, save, public begin_regional_restart
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)
real, public current_time_in_seconds
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
integer, parameter refine_ratio
subroutine, public regional_bc_t1_to_t0(BC_t1, BC_t0, nlev, ntracers, bnds)
integer, parameter, public h_stagger
integer, public p_step
integer, save bc_update_interval
character(len=100) grid_data
real, parameter stretch_factor
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:3437
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:3300
type(fv_regional_bc_variables), pointer, save bc_east_t1
subroutine, public regional_bc_data(Atm, bc_hour, is, ie, js, je, isd, ied, jsd, jed, ak, bk)
subroutine mp_auto_conversion(ql, qr, qi, qs)
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)
integer, parameter cube_res
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:2588
real, parameter tice
type(fv_domain_sides), target, save, public bc_t1
– Boundary values for all BC variables at successive times from the regional BC file ...
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;.
integer, save liq_water_index
integer, public n_step
logical, save north_bc
real function, public inner_prod(v1, v2)
subroutine check(status)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine, public get_unit_vect2(e1, e2, uc)
subroutine nudge_qv_bc(Atm, isd, ied, jsd, jed)
integer, parameter, public r_grid
Definition: fv_arrays.F90:35
integer, save, public next_time_to_read_bcs
integer, parameter parent_tile
subroutine read_regional_lon_lat
type(fv_regional_bc_variables), pointer, save bc_south_t1
real(kind=r_grid), dimension(:,:,:), allocatable grid_reg
– Lon/lat of cell corners
subroutine read_regional_bc_file(is_input, ie_input, js_input, je_input, nlev, ntracers, var_name_root, array_3d, array_4d, tlev, required)
type(fv_regional_bc_variables), pointer, save bc_south_t0
type(fv_domain_sides), target, save, public bc_t0
real, dimension(:), allocatable dum1d
type(fv_regional_bc_variables), pointer bc_side_t0
integer, parameter nhalo_data
subroutine bc_time_interpolation_general(is, ie, js, je, is_s, ie_s, js_s, je_s, is_w, ie_w, js_w, je_w, fraction, t0, t1, Atm)
character(len=80) data_source
real, parameter real_snan
subroutine, public set_regional_bcs(delp, delz, w, pt ifdef USE_COND
type(fv_regional_bc_variables), pointer, save bc_north_t1
subroutine, public setup_regional_bc(Atm, isd, ied, jsd, jed, npx, npy)
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)
logical, save west_bc
real, dimension(:), allocatable, public bk_in
The module &#39;fv_mapz&#39; contains the vertical mapping routines .
Definition: fv_mapz.F90:27
real, parameter zvir
real, dimension(:), allocatable pref
real, parameter target_lon
integer, save ncid
subroutine, public fillz(im, km, nq, q, dp)
The subroutine &#39;fillz&#39; is for mass-conservative filling of nonphysical negative values in the tracers...
Definition: fv_fill.F90:52
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
subroutine get_q00
subroutine p_maxmin(qname, q, is, ie, js, je, km, fac)
real, parameter t_i0
subroutine, public cell_center2(q1, q2, q3, q4, e2)
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)
The module &#39;fv_eta&#39; contains routine to set up the reference (Eulerian) pressure coordinate.
Definition: fv_eta.F90:25
character(len=100) oro_data
subroutine convert_to_virt_pot_temp(isd, ied, jsd, jed, npz, sphum, liq_wat)
type(fv_regional_bc_variables), pointer, save bc_east_t0
integer, parameter istart_nest
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:3540
logical, save south_bc
subroutine, public get_data_source(source, regional)
subroutine compute_regional_bc_indices(regional_bc_bounds)
type(fv_regional_bc_variables), pointer bc_side_t1
integer, save sphum_index
– Locations of tracer vbls in the tracers array
type(fv_regional_bc_bounds_type), pointer, save regional_bounds
subroutine fill_divgd_bc
real, dimension(:,:), allocatable phis_reg
– Filtered sfc geopotential
subroutine set_delp_and_tracers(BC_side, npz, nwat)
integer, parameter, public v_stagger
real, parameter c_liq
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
real, parameter cv_vap
The module &#39;fv_grid_utils&#39; contains routines for setting up and computing grid-related quantities...
real(kind=r_grid), parameter dbl_snan
integer, save, public bc_hour
type(fv_regional_bc_variables), pointer, save bc_west_t0
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, fcst_time, bc_update_interval)
type(fv_regional_bc_variables), pointer, save bc_north_t0
subroutine, public prt_gb_nh_sh(qname, is, ie, js, je, a2, area, lat)
real, parameter target_lat
real(kind=r_grid), dimension(:,:,:), allocatable agrid_reg
-- Lon/lat of cell centers
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)
subroutine compute_vpt
subroutine, public start_regional_restart(Atm, isc, iec, jsc, jec, isd, ied, jsd, jed)
integer, parameter jstart_nest
subroutine pmaxmn(qname, q, is, ie, js, je, km, fac, area, domain)
subroutine retrieve_bc_variable_data(bc_vbl_name, bc_side_t0, bc_side_t1, bc_t0, bc_t1, lbnd1, ubnd1, lbnd2, ubnd2, iq)
integer, public a_step
real, dimension(:), allocatable, public ak_in
integer, save, public ntimesteps_per_bc_update
logical, save east_bc
integer, parameter nhalo_model
subroutine, public exch_uv(domain, bd, npz, u, v)
subroutine read_regional_filtered_topo
integer, save ntracers
subroutine fillq(im, km, nq, q, dp)
subroutine, public start_regional_cold_start(Atm, ak, bk, levp, is, ie, js, je, isd, ied, jsd, jed)
type(fv_regional_bc_variables), pointer, save bc_west_t1
real, parameter cv_air
subroutine, public read_new_bc_data(Atm, Time, Time_step_atmos, p_split, isd, ied, jsd, jed)
integer, parameter, public u_stagger