FV3DYCORE  Version1.0.0
fv_tracer2d.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 
24 
25 ! Modules Included:
26 ! <table>
27 ! <tr>
28 ! <th>Module Name</th>
29 ! <th>Functions Included</th>
30 ! </tr>
31 ! <table>
32 ! <tr>
33 ! <td>boundary_mod</td>
34 ! <td>nested_grid_BC_apply_intT</td>
35 ! </tr>
36 ! <tr>
37 ! <td>fv_arrays_mod</td>
38 ! <td>fv_grid_type, fv_nest_type, fv_atmos_type, fv_grid_bounds_type</td>
39 ! </tr>
40 ! <tr>
41 ! <tr>
42 ! <td>fv_mp_mod</td>
43 ! <td>mp_reduce_max, ng, mp_gather, is_master, group_halo_update_type,
44 ! start_group_halo_update, complete_group_halo_update</td>
45 ! </tr>
46 ! <tr>
47 ! <td>fv_timing_mod</td>
48 ! <td>timing_on, timing_off</td>
49 ! </tr>
50 ! <tr>
51 ! <td>mpp_mod</td>
52 ! <td>mpp_error, FATAL, mpp_broadcast, mpp_send, mpp_recv, mpp_sum, mpp_max</td>
53 ! </tr>
54 ! <tr>
55 ! <td>mpp_domains_mod</td>
56 ! <td>mpp_update_domains, CGRID_NE, domain2d</td>
57 ! </tr>
58 ! <tr>
59 ! <td>tp_core_mod</td>
60 ! <td>fv_tp_2d, copy_corners</td>
61 ! </tr>
62 ! </table>
63 
64 module fv_tracer2d_mod
65  use tp_core_mod, only: fv_tp_2d, copy_corners
66  use fv_mp_mod, only: mp_reduce_max
67  use fv_mp_mod, only: ng, mp_gather, is_master
68  use fv_mp_mod, only: group_halo_update_type
69  use fv_mp_mod, only: start_group_halo_update, complete_group_halo_update
70  use mpp_domains_mod, only: mpp_update_domains, cgrid_ne, domain2d
71  use fv_timing_mod, only: timing_on, timing_off
72  use boundary_mod, only: nested_grid_bc_apply_intt
73  use fv_regional_mod, only: regional_boundary_update
74  use fv_regional_mod, only: current_time_in_seconds
75  use fv_arrays_mod, only: fv_grid_type, fv_nest_type, fv_atmos_type, fv_grid_bounds_type
76  use mpp_mod, only: mpp_error, fatal, mpp_broadcast, mpp_send, mpp_recv, mpp_sum, mpp_max
77 
78 implicit none
79 private
80 
81 public :: tracer_2d, tracer_2d_nested, tracer_2d_1l
82 
83 real, allocatable, dimension(:,:,:) :: nest_fx_west_accum, nest_fx_east_accum, nest_fx_south_accum, nest_fx_north_accum
84 
85 contains
86 
92 subroutine tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
93  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, regional)
94 
95  type(fv_grid_bounds_type), intent(IN) :: bd
96  integer, intent(IN) :: npx
97  integer, intent(IN) :: npy
98  integer, intent(IN) :: npz
99  integer, intent(IN) :: nq
100  integer, intent(IN) :: hord, nord_tr
101  integer, intent(IN) :: q_split
102  integer, intent(IN) :: id_divg
103  real , intent(IN) :: dt, trdm
104  real , intent(IN) :: lim_fac
105  logical, intent(IN) :: regional
106  type(group_halo_update_type), intent(inout) :: q_pack
107  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq)
108  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
109  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz)
110  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz)
111  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
112  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz)
113  type(fv_grid_type), intent(IN), target :: gridstruct
114  type(domain2d), intent(INOUT) :: domain
115 
116 ! Local Arrays
117  real :: qn2(bd%isd:bd%ied,bd%jsd:bd%jed,nq)
118  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
119  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
120  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
121  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
122  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
123  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
124  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
125  real :: cmax(npz)
126  real :: frac
127  integer :: nsplt
128  integer :: i,j,k,it,iq
129 
130  real, pointer, dimension(:,:) :: area, rarea
131  real, pointer, dimension(:,:,:) :: sin_sg
132  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
133 
134  integer :: is, ie, js, je
135  integer :: isd, ied, jsd, jed
136 
137  is = bd%is
138  ie = bd%ie
139  js = bd%js
140  je = bd%je
141  isd = bd%isd
142  ied = bd%ied
143  jsd = bd%jsd
144  jed = bd%jed
145 
146  area => gridstruct%area
147  rarea => gridstruct%rarea
148 
149  sin_sg => gridstruct%sin_sg
150  dxa => gridstruct%dxa
151  dya => gridstruct%dya
152  dx => gridstruct%dx
153  dy => gridstruct%dy
154 
155 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
156 !$OMP sin_sg,cy,yfx,dya,dx,cmax)
157  do k=1,npz
158  do j=jsd,jed
159  do i=is,ie+1
160  if (cx(i,j,k) > 0.) then
161  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
162  else
163  xfx(i,j,k) = cx(i,j,k)*dxa(i, j)*dy(i,j)*sin_sg(i, j,1)
164  endif
165  enddo
166  enddo
167  do j=js,je+1
168  do i=isd,ied
169  if (cy(i,j,k) > 0.) then
170  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
171  else
172  yfx(i,j,k) = cy(i,j,k)*dya(i,j )*dx(i,j)*sin_sg(i,j, 2)
173  endif
174  enddo
175  enddo
176 
177  cmax(k) = 0.
178  if ( k < npz/6 ) then
179  do j=js,je
180  do i=is,ie
181  cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) )
182  enddo
183  enddo
184  else
185  do j=js,je
186  do i=is,ie
187  cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) )
188  enddo
189  enddo
190  endif
191  enddo ! k-loop
192 
193  call mp_reduce_max(cmax,npz)
194 
195 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx, &
196 !$OMP cy,yfx,mfx,mfy,cmax) &
197 !$OMP private(nsplt, frac)
198  do k=1,npz
199 
200  nsplt = int(1. + cmax(k))
201  if ( nsplt > 1 ) then
202  frac = 1. / real(nsplt)
203  do j=jsd,jed
204  do i=is,ie+1
205  cx(i,j,k) = cx(i,j,k) * frac
206  xfx(i,j,k) = xfx(i,j,k) * frac
207  enddo
208  enddo
209  do j=js,je
210  do i=is,ie+1
211  mfx(i,j,k) = mfx(i,j,k) * frac
212  enddo
213  enddo
214  do j=js,je+1
215  do i=isd,ied
216  cy(i,j,k) = cy(i,j,k) * frac
217  yfx(i,j,k) = yfx(i,j,k) * frac
218  enddo
219  enddo
220  do j=js,je+1
221  do i=is,ie
222  mfy(i,j,k) = mfy(i,j,k) * frac
223  enddo
224  enddo
225  endif
226 
227  enddo
228  call timing_on('COMM_TOTAL')
229  call timing_on('COMM_TRACER')
230  call complete_group_halo_update(q_pack, domain)
231  call timing_off('COMM_TRACER')
232  call timing_off('COMM_TOTAL')
233 
234 ! Begin k-independent tracer transport; can not be OpenMPed because the mpp_update call.
235  do k=1,npz
236 
237 !$OMP parallel do default(none) shared(k,is,ie,js,je,isd,ied,jsd,jed,xfx,area,yfx,ra_x,ra_y)
238  do j=jsd,jed
239  do i=is,ie
240  ra_x(i,j) = area(i,j) + xfx(i,j,k) - xfx(i+1,j,k)
241  enddo
242  if ( j>=js .and. j<=je ) then
243  do i=isd,ied
244  ra_y(i,j) = area(i,j) + yfx(i,j,k) - yfx(i,j+1,k)
245  enddo
246  endif
247  enddo
248 
249  nsplt = int(1. + cmax(k))
250  do it=1,nsplt
251 
252 !$OMP parallel do default(none) shared(k,is,ie,js,je,rarea,mfx,mfy,dp1,dp2)
253  do j=js,je
254  do i=is,ie
255  dp2(i,j) = dp1(i,j,k) + (mfx(i,j,k)-mfx(i+1,j,k)+mfy(i,j,k)-mfy(i,j+1,k))*rarea(i,j)
256  enddo
257  enddo
258 
259 !$OMP parallel do default(none) shared(k,nsplt,it,is,ie,js,je,isd,ied,jsd,jed,npx,npy,cx,xfx,hord,trdm, &
260 !$OMP nord_tr,nq,gridstruct,bd,cy,yfx,mfx,mfy,qn2,q,ra_x,ra_y,dp1,dp2,rarea,lim_fac,regional) &
261 !$OMP private(fx,fy)
262  do iq=1,nq
263  if ( nsplt /= 1 ) then
264  if ( it==1 ) then
265  do j=jsd,jed
266  do i=isd,ied
267  qn2(i,j,iq) = q(i,j,k,iq)
268  enddo
269  enddo
270  endif
271  call fv_tp_2d(qn2(isd,jsd,iq), cx(is,jsd,k), cy(isd,js,k), &
272  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
273  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
274  if ( it < nsplt ) then ! not last call
275  do j=js,je
276  do i=is,ie
277  qn2(i,j,iq) = (qn2(i,j,iq)*dp1(i,j,k)+(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j))/dp2(i,j)
278  enddo
279  enddo
280  else
281  do j=js,je
282  do i=is,ie
283  q(i,j,k,iq) = (qn2(i,j,iq)*dp1(i,j,k)+(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j))/dp2(i,j)
284  enddo
285  enddo
286  endif
287  else
288  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
289  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
290  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
291  do j=js,je
292  do i=is,ie
293  q(i,j,k,iq) = (q(i,j,k,iq)*dp1(i,j,k)+(fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j))/dp2(i,j)
294  enddo
295  enddo
296  endif
297  enddo ! tracer-loop
298 
299  if ( it < nsplt ) then ! not last call
300  do j=js,je
301  do i=is,ie
302  dp1(i,j,k) = dp2(i,j)
303  enddo
304  enddo
305  call timing_on('COMM_TOTAL')
306  call timing_on('COMM_TRACER')
307  call mpp_update_domains(qn2, domain)
308  call timing_off('COMM_TRACER')
309  call timing_off('COMM_TOTAL')
310  endif
311  enddo ! time-split loop
312  enddo ! k-loop
313 
314 end subroutine tracer_2d_1l
315 
317 subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
318  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, regional)
319 
320  type(fv_grid_bounds_type), intent(IN) :: bd
321  integer, intent(IN) :: npx
322  integer, intent(IN) :: npy
323  integer, intent(IN) :: npz
324  integer, intent(IN) :: nq
325  integer, intent(IN) :: hord, nord_tr
326  integer, intent(IN) :: q_split
327  integer, intent(IN) :: id_divg
328  real , intent(IN) :: dt, trdm
329  real , intent(IN) :: lim_fac
330  logical, intent(IN) :: regional
331  type(group_halo_update_type), intent(inout) :: q_pack
332  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq)
333  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
334  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz)
335  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz)
336  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
337  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz)
338  type(fv_grid_type), intent(IN), target :: gridstruct
339  type(domain2d), intent(INOUT) :: domain
340 
341 ! Local Arrays
342  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
343  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
344  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
345  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
346  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
347  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
348  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
349  real :: cmax(npz)
350  real :: c_global
351  real :: frac, rdt
352  integer :: ksplt(npz)
353  integer :: nsplt
354  integer :: i,j,k,it,iq
355 
356  real, pointer, dimension(:,:) :: area, rarea
357  real, pointer, dimension(:,:,:) :: sin_sg
358  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
359 
360  integer :: is, ie, js, je
361  integer :: isd, ied, jsd, jed
362 
363  is = bd%is
364  ie = bd%ie
365  js = bd%js
366  je = bd%je
367  isd = bd%isd
368  ied = bd%ied
369  jsd = bd%jsd
370  jed = bd%jed
371 
372  area => gridstruct%area
373  rarea => gridstruct%rarea
374 
375  sin_sg => gridstruct%sin_sg
376  dxa => gridstruct%dxa
377  dya => gridstruct%dya
378  dx => gridstruct%dx
379  dy => gridstruct%dy
380 
381 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
382 !$OMP sin_sg,cy,yfx,dya,dx,cmax,q_split,ksplt)
383  do k=1,npz
384  do j=jsd,jed
385  do i=is,ie+1
386  if (cx(i,j,k) > 0.) then
387  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
388  else
389  xfx(i,j,k) = cx(i,j,k)*dxa(i,j)*dy(i,j)*sin_sg(i,j,1)
390  endif
391  enddo
392  enddo
393  do j=js,je+1
394  do i=isd,ied
395  if (cy(i,j,k) > 0.) then
396  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
397  else
398  yfx(i,j,k) = cy(i,j,k)*dya(i,j)*dx(i,j)*sin_sg(i,j,2)
399  endif
400  enddo
401  enddo
402 
403  if ( q_split == 0 ) then
404  cmax(k) = 0.
405  if ( k < npz/6 ) then
406  do j=js,je
407  do i=is,ie
408  cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) )
409  enddo
410  enddo
411  else
412  do j=js,je
413  do i=is,ie
414  cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) )
415  enddo
416  enddo
417  endif
418  endif
419  ksplt(k) = 1
420 
421  enddo
422 
423 !--------------------------------------------------------------------------------
424 
425 ! Determine global nsplt:
426  if ( q_split == 0 ) then
427  call mp_reduce_max(cmax,npz)
428 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
429  c_global = cmax(1)
430  if ( npz /= 1 ) then ! if NOT shallow water test case
431  do k=2,npz
432  c_global = max(cmax(k), c_global)
433  enddo
434  endif
435  nsplt = int(1. + c_global)
436  if ( is_master() .and. nsplt > 4 ) write(*,*) 'Tracer_2d_split=', nsplt, c_global
437  else
438  nsplt = q_split
439  endif
440 
441 !--------------------------------------------------------------------------------
442 
443  if( nsplt /= 1 ) then
444 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,mfx,cy,yfx,mfy,cmax,nsplt,ksplt) &
445 !$OMP private( frac )
446  do k=1,npz
447 
448 #ifdef GLOBAL_CFL
449  ksplt(k) = nsplt
450 #else
451  ksplt(k) = int(1. + cmax(k))
452 #endif
453  frac = 1. / real(ksplt(k))
454 
455  do j=jsd,jed
456  do i=is,ie+1
457  cx(i,j,k) = cx(i,j,k) * frac
458  xfx(i,j,k) = xfx(i,j,k) * frac
459  enddo
460  enddo
461  do j=js,je
462  do i=is,ie+1
463  mfx(i,j,k) = mfx(i,j,k) * frac
464  enddo
465  enddo
466 
467  do j=js,je+1
468  do i=isd,ied
469  cy(i,j,k) = cy(i,j,k) * frac
470  yfx(i,j,k) = yfx(i,j,k) * frac
471  enddo
472  enddo
473  do j=js,je+1
474  do i=is,ie
475  mfy(i,j,k) = mfy(i,j,k) * frac
476  enddo
477  enddo
478 
479  enddo
480  endif
481 
482  do it=1,nsplt
483  call timing_on('COMM_TOTAL')
484  call timing_on('COMM_TRACER')
485  call complete_group_halo_update(q_pack, domain)
486  call timing_off('COMM_TRACER')
487  call timing_off('COMM_TOTAL')
488 
489 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq,ksplt,&
490 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm,lim_fac,regional) &
491 !$OMP private(dp2, ra_x, ra_y, fx, fy)
492  do k=1,npz
493 
494  if ( it .le. ksplt(k) ) then
495 
496  do j=js,je
497  do i=is,ie
498  dp2(i,j) = dp1(i,j,k) + (mfx(i,j,k)-mfx(i+1,j,k)+mfy(i,j,k)-mfy(i,j+1,k))*rarea(i,j)
499  enddo
500  enddo
501 
502  do j=jsd,jed
503  do i=is,ie
504  ra_x(i,j) = area(i,j) + xfx(i,j,k) - xfx(i+1,j,k)
505  enddo
506  enddo
507  do j=js,je
508  do i=isd,ied
509  ra_y(i,j) = area(i,j) + yfx(i,j,k) - yfx(i,j+1,k)
510  enddo
511  enddo
512 
513  do iq=1,nq
514  if ( it==1 .and. trdm>1.e-4 ) then
515  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
516  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
517  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k), &
518  mass=dp1(isd,jsd,k), nord=nord_tr, damp_c=trdm)
519  else
520  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
521  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
522  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
523  endif
524  do j=js,je
525  do i=is,ie
526  q(i,j,k,iq) = ( q(i,j,k,iq)*dp1(i,j,k) + &
527  (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j) )/dp2(i,j)
528  enddo
529  enddo
530  enddo
531 
532  if ( it /= nsplt ) then
533  do j=js,je
534  do i=is,ie
535  dp1(i,j,k) = dp2(i,j)
536  enddo
537  enddo
538  endif
539 
540  endif ! ksplt
541 
542  enddo ! npz
543 
544  if ( it /= nsplt ) then
545  call timing_on('COMM_TOTAL')
546  call timing_on('COMM_TRACER')
547  call start_group_halo_update(q_pack, q, domain)
548  call timing_off('COMM_TRACER')
549  call timing_off('COMM_TOTAL')
550  endif
551 
552  enddo ! nsplt
553 
554 
555 end subroutine tracer_2d
556 
557 
558 subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
559  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, &
560  k_split, neststruct, parent_grid, lim_fac, regional)
561 
562  type(fv_grid_bounds_type), intent(IN) :: bd
563  integer, intent(IN) :: npx
564  integer, intent(IN) :: npy
565  integer, intent(IN) :: npz
566  integer, intent(IN) :: nq
567  integer, intent(IN) :: hord, nord_tr
568  integer, intent(IN) :: q_split, k_split
569  integer, intent(IN) :: id_divg
570  real , intent(IN) :: dt, trdm
571  real , intent(IN) :: lim_fac
572  logical, intent(IN) :: regional
573  type(group_halo_update_type), intent(inout) :: q_pack
574  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq)
575  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
576  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz)
577  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz)
578  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
579  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz)
580  type(fv_grid_type), intent(IN), target :: gridstruct
581  type(fv_nest_type), intent(INOUT) :: neststruct
582  type(fv_atmos_type), pointer, intent(IN) :: parent_grid
583  type(domain2d), intent(INOUT) :: domain
584 
585 ! Local Arrays
586  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
587  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
588  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
589  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
590  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
591  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
592  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
593  real :: cmax(npz)
594  real :: cmax_t
595  real :: c_global
596  real :: frac, rdt
597  real :: recip_nsplt,reg_bc_update_time
598  integer :: nsplt, nsplt_parent, msg_split_steps = 1
599  integer :: i,j,k,it,iq
600 
601  real, pointer, dimension(:,:) :: area, rarea
602  real, pointer, dimension(:,:,:) :: sin_sg
603  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
604 
605  integer :: is, ie, js, je
606  integer :: isd, ied, jsd, jed
607 
608  is = bd%is
609  ie = bd%ie
610  js = bd%js
611  je = bd%je
612  isd = bd%isd
613  ied = bd%ied
614  jsd = bd%jsd
615  jed = bd%jed
616 
617  area => gridstruct%area
618  rarea => gridstruct%rarea
619 
620  sin_sg => gridstruct%sin_sg
621  dxa => gridstruct%dxa
622  dya => gridstruct%dya
623  dx => gridstruct%dx
624  dy => gridstruct%dy
625 
626 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
627 !$OMP sin_sg,cy,yfx,dya,dx)
628  do k=1,npz
629  do j=jsd,jed
630  do i=is,ie+1
631  if (cx(i,j,k) > 0.) then
632  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
633  else
634  xfx(i,j,k) = cx(i,j,k)*dxa(i,j)*dy(i,j)*sin_sg(i,j,1)
635  endif
636  enddo
637  enddo
638  do j=js,je+1
639  do i=isd,ied
640  if (cy(i,j,k) > 0.) then
641  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
642  else
643  yfx(i,j,k) = cy(i,j,k)*dya(i,j)*dx(i,j)*sin_sg(i,j,2)
644  endif
645  enddo
646  enddo
647  enddo
648 
649 !--------------------------------------------------------------------------------
650  if ( q_split == 0 ) then
651 ! Determine nsplt
652 
653 !$OMP parallel do default(none) shared(is,ie,js,je,npz,cmax,cx,cy,sin_sg) &
654 !$OMP private(cmax_t )
655  do k=1,npz
656  cmax(k) = 0.
657  if ( k < 4 ) then
658 ! Top layers: C < max( abs(c_x), abs(c_y) )
659  do j=js,je
660  do i=is,ie
661  cmax_t = max( abs(cx(i,j,k)), abs(cy(i,j,k)) )
662  cmax(k) = max( cmax_t, cmax(k) )
663  enddo
664  enddo
665  else
666  do j=js,je
667  do i=is,ie
668  cmax_t = max(abs(cx(i,j,k)), abs(cy(i,j,k))) + 1.-sin_sg(i,j,5)
669  cmax(k) = max( cmax_t, cmax(k) )
670  enddo
671  enddo
672  endif
673  enddo
674  call mp_reduce_max(cmax,npz)
675 
676 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
677  c_global = cmax(1)
678  if ( npz /= 1 ) then ! if NOT shallow water test case
679  do k=2,npz
680  c_global = max(cmax(k), c_global)
681  enddo
682  endif
683  nsplt = int(1. + c_global)
684  if ( is_master() .and. nsplt > 3 ) write(*,*) 'Tracer_2d_split=', nsplt, c_global
685  else
686  nsplt = q_split
687  if (gridstruct%nested .and. neststruct%nestbctype > 1) msg_split_steps = max(q_split/parent_grid%flagstruct%q_split,1)
688  endif
689 
690 !--------------------------------------------------------------------------------
691 
692  frac = 1. / real(nsplt)
693  recip_nsplt = 1. / real(nsplt)
694 
695  if( nsplt /= 1 ) then
696 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,frac,xfx,mfx,cy,yfx,mfy)
697  do k=1,npz
698  do j=jsd,jed
699  do i=is,ie+1
700  cx(i,j,k) = cx(i,j,k) * frac
701  xfx(i,j,k) = xfx(i,j,k) * frac
702  enddo
703  enddo
704  do j=js,je
705  do i=is,ie+1
706  mfx(i,j,k) = mfx(i,j,k) * frac
707  enddo
708  enddo
709 
710  do j=js,je+1
711  do i=isd,ied
712  cy(i,j,k) = cy(i,j,k) * frac
713  yfx(i,j,k) = yfx(i,j,k) * frac
714  enddo
715  enddo
716 
717  do j=js,je+1
718  do i=is,ie
719  mfy(i,j,k) = mfy(i,j,k) * frac
720  enddo
721  enddo
722  enddo
723  endif
724 
725 
726  do it=1,nsplt
727  if ( gridstruct%nested ) then
728  neststruct%tracer_nest_timestep = neststruct%tracer_nest_timestep + 1
729  end if
730  call timing_on('COMM_TOTAL')
731  call timing_on('COMM_TRACER')
732  call complete_group_halo_update(q_pack, domain)
733  call timing_off('COMM_TRACER')
734  call timing_off('COMM_TOTAL')
735 
736  if (gridstruct%nested) then
737  do iq=1,nq
738  call nested_grid_bc_apply_intt(q(isd:ied,jsd:jed,:,iq), &
739  0, 0, npx, npy, npz, bd, &
740  real(neststruct%tracer_nest_timestep)+real(nsplt*k_split), real(nsplt*k_split), &
741  neststruct%q_bc(iq), bctype=neststruct%nestbctype )
742  enddo
743  endif
744 
745  if (regional) then
746  reg_bc_update_time=current_time_in_seconds+(it-1)*recip_nsplt*dt
747  do iq=1,nq
748  call regional_boundary_update(q(:,:,:,iq), 'q', &
749  isd, ied, jsd, jed, npz, &
750  is, ie, js, je, &
751  isd, ied, jsd, jed, &
752  reg_bc_update_time, &
753  iq )
754  enddo
755  endif
756 
757 
758 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq, &
759 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm,lim_fac,regional) &
760 !$OMP private(dp2, ra_x, ra_y, fx, fy)
761  do k=1,npz
762 
763  do j=js,je
764  do i=is,ie
765  dp2(i,j) = dp1(i,j,k) + (mfx(i,j,k)-mfx(i+1,j,k)+mfy(i,j,k)-mfy(i,j+1,k))*rarea(i,j)
766  enddo
767  enddo
768 
769  do j=jsd,jed
770  do i=is,ie
771  ra_x(i,j) = area(i,j) + xfx(i,j,k) - xfx(i+1,j,k)
772  enddo
773  enddo
774  do j=js,je
775  do i=isd,ied
776  ra_y(i,j) = area(i,j) + yfx(i,j,k) - yfx(i,j+1,k)
777  enddo
778  enddo
779 
780  do iq=1,nq
781  if ( it==1 .and. trdm>1.e-4 ) then
782  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
783  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
784  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k), &
785  mass=dp1(isd,jsd,k), nord=nord_tr, damp_c=trdm)
786  else
787  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
788  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
789  gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
790  endif
791  do j=js,je
792  do i=is,ie
793  q(i,j,k,iq) = ( q(i,j,k,iq)*dp1(i,j,k) + &
794  (fx(i,j)-fx(i+1,j)+fy(i,j)-fy(i,j+1))*rarea(i,j) )/dp2(i,j)
795  enddo
796  enddo
797  enddo
798  enddo ! npz
799 
800  if ( it /= nsplt ) then
801  call timing_on('COMM_TOTAL')
802  call timing_on('COMM_TRACER')
803  call start_group_halo_update(q_pack, q, domain)
804  call timing_off('COMM_TRACER')
805  call timing_off('COMM_TOTAL')
806  endif
807  !Apply nested-grid BCs
808  if ( gridstruct%nested ) then
809  do iq=1,nq
810 
811 
812  call nested_grid_bc_apply_intt(q(isd:ied,jsd:jed,:,iq), &
813  0, 0, npx, npy, npz, bd, &
814  real(neststruct%tracer_nest_timestep), real(nsplt*k_split), &
815  neststruct%q_bc(iq), bctype=neststruct%nestbctype )
816 
817  end do
818  end if
819 
820 ! BCs for q at the current time were applied above for the regional mode.
821 
822  enddo ! nsplt
823 
824  if ( id_divg > 0 ) then
825  rdt = 1./(frac*dt)
826 
827 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dp1,xfx,yfx,rarea,rdt)
828  do k=1,npz
829  do j=js,je
830  do i=is,ie
831  dp1(i,j,k) = (xfx(i+1,j,k)-xfx(i,j,k) + yfx(i,j+1,k)-yfx(i,j,k))*rarea(i,j)*rdt
832  enddo
833  enddo
834  enddo
835  endif
836 
837  end subroutine tracer_2d_nested
838 
839 end module fv_tracer2d_mod
fv_regional_mod::regional_boundary_update
subroutine, public regional_boundary_update(array, bc_vbl_name, lbnd_x, ubnd_x, lbnd_y, ubnd_y, ubnd_z, is, ie, js, je, isd, ied, jsd, jed, fcst_time, index4)
Definition: fv_regional_bc.F90:3523
fv_regional_mod::current_time_in_seconds
real, public current_time_in_seconds
Definition: fv_regional_bc.F90:108
fv_mp_mod
The module &#39;fv_mp_mod&#39; is a single program multiple data (SPMD) parallel decompostion/communication m...
Definition: fv_mp_mod.F90:24
fv_timing_mod::timing_off
subroutine timing_off(blk_name)
The subroutine &#39;timing_off&#39; stops a timer.
Definition: fv_timing.F90:180
fv_arrays_mod::fv_grid_type
The type &#39;fv_grid_type&#39; is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
Definition: fv_arrays.F90:120
fv_arrays_mod::fv_grid_bounds_type
Definition: fv_arrays.F90:1137
fv_tracer2d_mod::nest_fx_west_accum
real, dimension(:,:,:), allocatable nest_fx_west_accum
Definition: fv_tracer2d.F90:83
fv_tracer2d_mod::nest_fx_south_accum
real, dimension(:,:,:), allocatable nest_fx_south_accum
Definition: fv_tracer2d.F90:83
fv_timing_mod
The module &#39;fv_timing&#39; contains FV3 timers.
Definition: fv_timing.F90:24
boundary_mod
The module &#39;boundary&#39; contains utility routines for grid nesting and boundary conditions.
Definition: boundary.F90:25
tp_core_mod
The module &#39;tp_core&#39; is a collection of routines to support FV transport.
Definition: tp_core.F90:24
fv_tracer2d_mod
The module &#39;fv_tracer2d.F90&#39; performs sub-cycled tracer advection.
Definition: fv_tracer2d.F90:64
tp_core_mod::copy_corners
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
Definition: tp_core.F90:248
fv_arrays_mod
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
fv_tracer2d_mod::tracer_2d
subroutine, public tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, regional)
The subroutine &#39;tracer_2d&#39; is the standard routine for sub-cycled tracer advection.
Definition: fv_tracer2d.F90:319
fv_regional_mod
Definition: fv_regional_bc.F90:22
fv_timing_mod::timing_on
subroutine timing_on(blk_name)
The subroutine &#39;timing_on&#39; starts a timer.
Definition: fv_timing.F90:116
fv_mp_mod::ng
integer, parameter, public ng
Definition: fv_mp_mod.F90:2716
fv_tracer2d_mod::nest_fx_east_accum
real, dimension(:,:,:), allocatable nest_fx_east_accum
Definition: fv_tracer2d.F90:83
tp_core_mod::fv_tp_2d
subroutine, public fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, lim_fac, regional, mfx, mfy, mass, nord, damp_c)
The subroutine &#39;fv_tp_2d&#39; contains the FV advection scheme .
Definition: tp_core.F90:110
fv_arrays_mod::fv_atmos_type
Definition: fv_arrays.F90:1166
fv_tracer2d_mod::tracer_2d_1l
subroutine, public tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, lim_fac, regional)
The subroutine &#39;tracer_2d_1L&#39; performs 2-D horizontal-to-lagrangian transport.
Definition: fv_tracer2d.F90:94
boundary_mod::nested_grid_bc_apply_intt
subroutine, public nested_grid_bc_apply_intt(var_nest, istag, jstag, npx, npy, npz, bd, step, split, BC, bctype)
The subroutine &#39;nested_grid_BC_apply_intT&#39; performs linear interpolation or extrapolation in time for...
Definition: boundary.F90:1684
fv_tracer2d_mod::nest_fx_north_accum
real, dimension(:,:,:), allocatable nest_fx_north_accum
Definition: fv_tracer2d.F90:83
fv_arrays_mod::fv_nest_type
Definition: fv_arrays.F90:1047
fv_tracer2d_mod::tracer_2d_nested
subroutine, public tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, k_split, neststruct, parent_grid, lim_fac, regional)
Definition: fv_tracer2d.F90:561