FV3 Bundle
fv_tracer2d_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
22  use fv_mp_nlm_mod, only: mp_reduce_max
23  use fv_mp_nlm_mod, only: ng, mp_gather, is_master
24  use fv_mp_nlm_mod, only: group_halo_update_type
25  use fv_mp_nlm_mod, only: start_group_halo_update, complete_group_halo_update
31 
32 implicit none
33 private
34 
36 
38 
39 contains
40 
41 !-----------------------------------------------------------------------
42 ! !ROUTINE: Perform 2D horizontal-to-lagrangian transport
43 !-----------------------------------------------------------------------
44 
45 
46 
47 subroutine tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
48  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, dpA)
49 
50  type(fv_grid_bounds_type), intent(IN) :: bd
51  integer, intent(IN) :: npx
52  integer, intent(IN) :: npy
53  integer, intent(IN) :: npz
54  integer, intent(IN) :: nq ! number of tracers to be advected
55  integer, intent(IN) :: hord, nord_tr
56  integer, intent(IN) :: q_split
57  integer, intent(IN) :: id_divg
58  real , intent(IN) :: dt, trdm
59  type(group_halo_update_type), intent(inout) :: q_pack
60  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) ! Tracers
61  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) ! DELP before dyn_core
62  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) ! Mass Flux X-Dir
63  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) ! Mass Flux Y-Dir
64  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir
65  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) ! Courant Number Y-Dir
66  real , optional, intent(OUT) :: dpa(bd%is:bd%ie,bd%js:bd%je) ! DELP after advection
67  type(fv_grid_type), intent(IN), target :: gridstruct
68  type(domain2d), intent(INOUT) :: domain
69 
70 ! Local Arrays
71  real :: qn2(bd%isd:bd%ied,bd%jsd:bd%jed,nq) ! 3D tracers
72  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
73  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
74  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
75  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
76  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
77  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
78  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
79  real :: cmax(npz)
80  real :: frac
81  integer :: nsplt
82  integer :: i,j,k,it,iq
83 
84  real, pointer, dimension(:,:) :: area, rarea
85  real, pointer, dimension(:,:,:) :: sin_sg
86  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
87 
88  integer :: is, ie, js, je
89  integer :: isd, ied, jsd, jed
90 
91  is = bd%is
92  ie = bd%ie
93  js = bd%js
94  je = bd%je
95  isd = bd%isd
96  ied = bd%ied
97  jsd = bd%jsd
98  jed = bd%jed
99 
100  area => gridstruct%area
101  rarea => gridstruct%rarea
102 
103  sin_sg => gridstruct%sin_sg
104  dxa => gridstruct%dxa
105  dya => gridstruct%dya
106  dx => gridstruct%dx
107  dy => gridstruct%dy
108 
109 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
110 !$OMP sin_sg,cy,yfx,dya,dx,cmax)
111  do k=1,npz
112  do j=jsd,jed
113  do i=is,ie+1
114  if (cx(i,j,k) > 0.) then
115  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
116  else
117  xfx(i,j,k) = cx(i,j,k)*dxa(i, j)*dy(i,j)*sin_sg(i, j,1)
118  endif
119  enddo
120  enddo
121  do j=js,je+1
122  do i=isd,ied
123  if (cy(i,j,k) > 0.) then
124  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
125  else
126  yfx(i,j,k) = cy(i,j,k)*dya(i,j )*dx(i,j)*sin_sg(i,j, 2)
127  endif
128  enddo
129  enddo
130 
131  cmax(k) = 0.
132  if ( k < npz/6 ) then
133  do j=js,je
134  do i=is,ie
135  cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) )
136  enddo
137  enddo
138  else
139  do j=js,je
140  do i=is,ie
141  cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) )
142  enddo
143  enddo
144  endif
145  enddo ! k-loop
146 
147  call mp_reduce_max(cmax,npz)
148 
149 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx, &
150 !$OMP cy,yfx,mfx,mfy,cmax) &
151 !$OMP private(nsplt, frac)
152  do k=1,npz
153 
154  nsplt = int(1. + cmax(k))
155  if ( nsplt > 1 ) then
156  frac = 1. / real(nsplt)
157  do j=jsd,jed
158  do i=is,ie+1
159  cx(i,j,k) = cx(i,j,k) * frac
160  xfx(i,j,k) = xfx(i,j,k) * frac
161  enddo
162  enddo
163  do j=js,je
164  do i=is,ie+1
165  mfx(i,j,k) = mfx(i,j,k) * frac
166  enddo
167  enddo
168  do j=js,je+1
169  do i=isd,ied
170  cy(i,j,k) = cy(i,j,k) * frac
171  yfx(i,j,k) = yfx(i,j,k) * frac
172  enddo
173  enddo
174  do j=js,je+1
175  do i=is,ie
176  mfy(i,j,k) = mfy(i,j,k) * frac
177  enddo
178  enddo
179  endif
180 
181  enddo
182  call timing_on('COMM_TOTAL')
183  call timing_on('COMM_TRACER')
184  call complete_group_halo_update(q_pack, domain)
185  call timing_off('COMM_TRACER')
186  call timing_off('COMM_TOTAL')
187 
188 ! Begin k-independent tracer transport; can not be OpenMPed because the mpp_update call.
189  do k=1,npz
190 
191 !$OMP parallel do default(none) shared(k,is,ie,js,je,isd,ied,jsd,jed,xfx,area,yfx,ra_x,ra_y)
192  do j=jsd,jed
193  do i=is,ie
194  ra_x(i,j) = area(i,j) + (xfx(i,j,k) - xfx(i+1,j,k))
195  enddo
196  if ( j>=js .and. j<=je ) then
197  do i=isd,ied
198  ra_y(i,j) = area(i,j) + (yfx(i,j,k) - yfx(i,j+1,k))
199  enddo
200  endif
201  enddo
202 
203  nsplt = int(1. + cmax(k))
204  do it=1,nsplt
205 
206 !$OMP parallel do default(none) shared(k,is,ie,js,je,rarea,mfx,mfy,dp1,dp2)
207  do j=js,je
208  do i=is,ie
209  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)
210  enddo
211  enddo
212 
213 !$OMP parallel do default(none) shared(k,nsplt,it,is,ie,js,je,isd,ied,jsd,jed,npx,npy,cx,xfx,hord,trdm, &
214 !$OMP nord_tr,nq,gridstruct,bd,cy,yfx,mfx,mfy,qn2,q,ra_x,ra_y,dp1,dp2,rarea) &
215 !$OMP private(fx,fy)
216  do iq=1,nq
217  if ( nsplt /= 1 ) then
218  if ( it==1 ) then
219  do j=jsd,jed
220  do i=isd,ied
221  qn2(i,j,iq) = q(i,j,k,iq)
222  enddo
223  enddo
224  endif
225  call fv_tp_2d(qn2(isd,jsd,iq), cx(is,jsd,k), cy(isd,js,k), &
226  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
227  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
228  if ( it < nsplt ) then ! not last call
229  do j=js,je
230  do i=is,ie
231  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)
232  enddo
233  enddo
234  else
235  do j=js,je
236  do i=is,ie
237  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)
238  enddo
239  enddo
240  endif
241  else
242  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
243  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
244  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
245  do j=js,je
246  do i=is,ie
247  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)
248  enddo
249  enddo
250  endif
251  enddo ! tracer-loop
252 
253  if ( it < nsplt ) then ! not last call
254  do j=js,je
255  do i=is,ie
256  dp1(i,j,k) = dp2(i,j)
257  enddo
258  enddo
259  call timing_on('COMM_TOTAL')
260  call timing_on('COMM_TRACER')
261  call mpp_update_domains(qn2, domain)
262  call timing_off('COMM_TRACER')
263  call timing_off('COMM_TOTAL')
264  endif
265  enddo ! time-split loop
266  enddo ! k-loop
267 
268  if (present(dpa)) then
269  dpa=dp2
270  endif
271 
272 end subroutine tracer_2d_1l
273 
274 
275 subroutine tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
276  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, dpA)
278  type(fv_grid_bounds_type), intent(IN) :: bd
279  integer, intent(IN) :: npx
280  integer, intent(IN) :: npy
281  integer, intent(IN) :: npz
282  integer, intent(IN) :: nq ! number of tracers to be advected
283  integer, intent(IN) :: hord, nord_tr
284  integer, intent(IN) :: q_split
285  integer, intent(IN) :: id_divg
286  real , intent(IN) :: dt, trdm
287  type(group_halo_update_type), intent(inout) :: q_pack
288  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) ! Tracers
289  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) ! DELP before dyn_core
290  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) ! Mass Flux X-Dir
291  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) ! Mass Flux Y-Dir
292  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir
293  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) ! Courant Number Y-Dir
294  real , optional, intent(OUT) :: dpa(bd%is:bd%ie,bd%js:bd%je,npz)! DELP after advection
295  type(fv_grid_type), intent(IN), target :: gridstruct
296  type(domain2d), intent(INOUT) :: domain
297 
298 ! Local Arrays
299  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
300  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
301  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
302  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
303  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
304  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
305  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
306  real :: cmax(npz)
307  real :: c_global
308  real :: frac, rdt
309  integer :: ksplt(npz)
310  integer :: nsplt
311  integer :: i,j,k,it,iq
312 
313  real, pointer, dimension(:,:) :: area, rarea
314  real, pointer, dimension(:,:,:) :: sin_sg
315  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
316 
317  integer :: is, ie, js, je
318  integer :: isd, ied, jsd, jed
319 
320  is = bd%is
321  ie = bd%ie
322  js = bd%js
323  je = bd%je
324  isd = bd%isd
325  ied = bd%ied
326  jsd = bd%jsd
327  jed = bd%jed
328 
329  area => gridstruct%area
330  rarea => gridstruct%rarea
331 
332  sin_sg => gridstruct%sin_sg
333  dxa => gridstruct%dxa
334  dya => gridstruct%dya
335  dx => gridstruct%dx
336  dy => gridstruct%dy
337 
338 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
339 !$OMP sin_sg,cy,yfx,dya,dx,cmax,q_split,ksplt)
340  do k=1,npz
341  do j=jsd,jed
342  do i=is,ie+1
343  if (cx(i,j,k) > 0.) then
344  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
345  else
346  xfx(i,j,k) = cx(i,j,k)*dxa(i,j)*dy(i,j)*sin_sg(i,j,1)
347  endif
348  enddo
349  enddo
350  do j=js,je+1
351  do i=isd,ied
352  if (cy(i,j,k) > 0.) then
353  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
354  else
355  yfx(i,j,k) = cy(i,j,k)*dya(i,j)*dx(i,j)*sin_sg(i,j,2)
356  endif
357  enddo
358  enddo
359 
360  if ( q_split == 0 ) then
361  cmax(k) = 0.
362  if ( k < npz/6 ) then
363  do j=js,je
364  do i=is,ie
365  cmax(k) = max( cmax(k), abs(cx(i,j,k)), abs(cy(i,j,k)) )
366  enddo
367  enddo
368  else
369  do j=js,je
370  do i=is,ie
371  cmax(k) = max( cmax(k), max(abs(cx(i,j,k)),abs(cy(i,j,k)))+1.-sin_sg(i,j,5) )
372  enddo
373  enddo
374  endif
375  endif
376  ksplt(k) = 1
377 
378  enddo
379 
380 !--------------------------------------------------------------------------------
381 
382 ! Determine global nsplt:
383  if ( q_split == 0 ) then
384  call mp_reduce_max(cmax,npz)
385 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
386  c_global = cmax(1)
387  if ( npz /= 1 ) then ! if NOT shallow water test case
388  do k=2,npz
389  c_global = max(cmax(k), c_global)
390  enddo
391  endif
392  nsplt = int(1. + c_global)
393  if ( is_master() .and. nsplt > 4 ) write(*,*) 'Tracer_2d_split=', nsplt, c_global
394  else
395  nsplt = q_split
396  endif
397 
398 !--------------------------------------------------------------------------------
399 
400  if( nsplt /= 1 ) then
401 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,mfx,cy,yfx,mfy,cmax,nsplt,ksplt) &
402 !$OMP private( frac )
403  do k=1,npz
404 
405 #ifdef GLOBAL_CFL
406  ksplt(k) = nsplt
407 #else
408  ksplt(k) = int(1. + cmax(k))
409 #endif
410  frac = 1. / real(ksplt(k))
411 
412  do j=jsd,jed
413  do i=is,ie+1
414  cx(i,j,k) = cx(i,j,k) * frac
415  xfx(i,j,k) = xfx(i,j,k) * frac
416  enddo
417  enddo
418  do j=js,je
419  do i=is,ie+1
420  mfx(i,j,k) = mfx(i,j,k) * frac
421  enddo
422  enddo
423 
424  do j=js,je+1
425  do i=isd,ied
426  cy(i,j,k) = cy(i,j,k) * frac
427  yfx(i,j,k) = yfx(i,j,k) * frac
428  enddo
429  enddo
430  do j=js,je+1
431  do i=is,ie
432  mfy(i,j,k) = mfy(i,j,k) * frac
433  enddo
434  enddo
435 
436  enddo
437  endif
438 
439  do it=1,nsplt
440  call timing_on('COMM_TOTAL')
441  call timing_on('COMM_TRACER')
442  call complete_group_halo_update(q_pack, domain)
443  call timing_off('COMM_TRACER')
444  call timing_off('COMM_TOTAL')
445 
446 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq,ksplt,&
447 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
448 !$OMP private(dp2, ra_x, ra_y, fx, fy)
449  do k=1,npz
450 
451  if ( it .le. ksplt(k) ) then
452 
453  do j=js,je
454  do i=is,ie
455  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)
456  enddo
457  enddo
458 
459  do j=jsd,jed
460  do i=is,ie
461  ra_x(i,j) = area(i,j) + (xfx(i,j,k) - xfx(i+1,j,k))
462  enddo
463  enddo
464  do j=js,je
465  do i=isd,ied
466  ra_y(i,j) = area(i,j) + (yfx(i,j,k) - yfx(i,j+1,k))
467  enddo
468  enddo
469 
470  do iq=1,nq
471  if ( it==1 .and. trdm>1.e-4 ) then
472  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
473  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
474  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k), &
475  mass=dp1(isd,jsd,k), nord=nord_tr, damp_c=trdm)
476  else
477  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
478  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
479  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
480  endif
481  do j=js,je
482  do i=is,ie
483  q(i,j,k,iq) = ( q(i,j,k,iq)*dp1(i,j,k) + &
484  ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*rarea(i,j) )/dp2(i,j)
485  enddo
486  enddo
487  enddo
488 
489  if ( it /= nsplt ) then
490  do j=js,je
491  do i=is,ie
492  dp1(i,j,k) = dp2(i,j)
493  enddo
494  enddo
495  endif
496 
497  endif ! ksplt
498 
499  enddo ! npz
500 
501  if ( it /= nsplt ) then
502  call timing_on('COMM_TOTAL')
503  call timing_on('COMM_TRACER')
504  call start_group_halo_update(q_pack, q, domain)
505  call timing_off('COMM_TRACER')
506  call timing_off('COMM_TOTAL')
507  endif
508 
509  enddo ! nsplt
510 
511 
512  if (present(dpa)) then
513  dpa=dp1(bd%is:bd%ie,bd%js:bd%je,1:npz)
514  endif
515 
516 end subroutine tracer_2d
517 
518 subroutine tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain, npx, npy, npz, &
519  nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, &
520  k_split, neststruct, parent_grid)
522  type(fv_grid_bounds_type), intent(IN) :: bd
523  integer, intent(IN) :: npx
524  integer, intent(IN) :: npy
525  integer, intent(IN) :: npz
526  integer, intent(IN) :: nq ! number of tracers to be advected
527  integer, intent(IN) :: hord, nord_tr
528  integer, intent(IN) :: q_split, k_split
529  integer, intent(IN) :: id_divg
530  real , intent(IN) :: dt, trdm
531  type(group_halo_update_type), intent(inout) :: q_pack
532  real , intent(INOUT) :: q(bd%isd:bd%ied,bd%jsd:bd%jed,npz,nq) ! Tracers
533  real , intent(INOUT) :: dp1(bd%isd:bd%ied,bd%jsd:bd%jed,npz) ! DELP before dyn_core
534  real , intent(INOUT) :: mfx(bd%is:bd%ie+1,bd%js:bd%je, npz) ! Mass Flux X-Dir
535  real , intent(INOUT) :: mfy(bd%is:bd%ie ,bd%js:bd%je+1,npz) ! Mass Flux Y-Dir
536  real , intent(INOUT) :: cx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir
537  real , intent(INOUT) :: cy(bd%isd:bd%ied,bd%js :bd%je +1,npz) ! Courant Number Y-Dir
538  type(fv_grid_type), intent(IN), target :: gridstruct
539  type(fv_nest_type), intent(INOUT) :: neststruct
540  type(fv_atmos_type), intent(INOUT) :: parent_grid
541  type(domain2d), intent(INOUT) :: domain
542 
543 ! Local Arrays
544  real :: dp2(bd%is:bd%ie,bd%js:bd%je)
545  real :: fx(bd%is:bd%ie+1,bd%js:bd%je )
546  real :: fy(bd%is:bd%ie , bd%js:bd%je+1)
547  real :: ra_x(bd%is:bd%ie,bd%jsd:bd%jed)
548  real :: ra_y(bd%isd:bd%ied,bd%js:bd%je)
549  real :: xfx(bd%is:bd%ie+1,bd%jsd:bd%jed ,npz)
550  real :: yfx(bd%isd:bd%ied,bd%js: bd%je+1, npz)
551  real :: cmax(npz)
552  real :: cmax_t
553  real :: c_global
554  real :: frac, rdt
555  integer :: nsplt, nsplt_parent, msg_split_steps = 1
556  integer :: i,j,k,it,iq
557 
558  real, pointer, dimension(:,:) :: area, rarea
559  real, pointer, dimension(:,:,:) :: sin_sg
560  real, pointer, dimension(:,:) :: dxa, dya, dx, dy
561 
562  integer :: is, ie, js, je
563  integer :: isd, ied, jsd, jed
564 
565  is = bd%is
566  ie = bd%ie
567  js = bd%js
568  je = bd%je
569  isd = bd%isd
570  ied = bd%ied
571  jsd = bd%jsd
572  jed = bd%jed
573 
574  area => gridstruct%area
575  rarea => gridstruct%rarea
576 
577  sin_sg => gridstruct%sin_sg
578  dxa => gridstruct%dxa
579  dya => gridstruct%dya
580  dx => gridstruct%dx
581  dy => gridstruct%dy
582 
583 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
584 !$OMP sin_sg,cy,yfx,dya,dx)
585  do k=1,npz
586  do j=jsd,jed
587  do i=is,ie+1
588  if (cx(i,j,k) > 0.) then
589  xfx(i,j,k) = cx(i,j,k)*dxa(i-1,j)*dy(i,j)*sin_sg(i-1,j,3)
590  else
591  xfx(i,j,k) = cx(i,j,k)*dxa(i,j)*dy(i,j)*sin_sg(i,j,1)
592  endif
593  enddo
594  enddo
595  do j=js,je+1
596  do i=isd,ied
597  if (cy(i,j,k) > 0.) then
598  yfx(i,j,k) = cy(i,j,k)*dya(i,j-1)*dx(i,j)*sin_sg(i,j-1,4)
599  else
600  yfx(i,j,k) = cy(i,j,k)*dya(i,j)*dx(i,j)*sin_sg(i,j,2)
601  endif
602  enddo
603  enddo
604  enddo
605 
606 !--------------------------------------------------------------------------------
607  if ( q_split == 0 ) then
608 ! Determine nsplt
609 
610 !$OMP parallel do default(none) shared(is,ie,js,je,npz,cmax,cx,cy,sin_sg) &
611 !$OMP private(cmax_t )
612  do k=1,npz
613  cmax(k) = 0.
614  if ( k < 4 ) then
615 ! Top layers: C < max( abs(c_x), abs(c_y) )
616  do j=js,je
617  do i=is,ie
618  cmax_t = max( abs(cx(i,j,k)), abs(cy(i,j,k)) )
619  cmax(k) = max( cmax_t, cmax(k) )
620  enddo
621  enddo
622  else
623  do j=js,je
624  do i=is,ie
625  cmax_t = max(abs(cx(i,j,k)), abs(cy(i,j,k))) + 1.-sin_sg(i,j,5)
626  cmax(k) = max( cmax_t, cmax(k) )
627  enddo
628  enddo
629  endif
630  enddo
631  call mp_reduce_max(cmax,npz)
632 
633 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
634  c_global = cmax(1)
635  if ( npz /= 1 ) then ! if NOT shallow water test case
636  do k=2,npz
637  c_global = max(cmax(k), c_global)
638  enddo
639  endif
640  nsplt = int(1. + c_global)
641  if ( is_master() .and. nsplt > 3 ) write(*,*) 'Tracer_2d_split=', nsplt, c_global
642  else
643  nsplt = q_split
644  if (gridstruct%nested .and. neststruct%nestbctype > 1) msg_split_steps = max(q_split/parent_grid%flagstruct%q_split,1)
645  endif
646 
647 !--------------------------------------------------------------------------------
648 
649  frac = 1. / real(nsplt)
650 
651  if( nsplt /= 1 ) then
652 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,frac,xfx,mfx,cy,yfx,mfy)
653  do k=1,npz
654  do j=jsd,jed
655  do i=is,ie+1
656  cx(i,j,k) = cx(i,j,k) * frac
657  xfx(i,j,k) = xfx(i,j,k) * frac
658  enddo
659  enddo
660  do j=js,je
661  do i=is,ie+1
662  mfx(i,j,k) = mfx(i,j,k) * frac
663  enddo
664  enddo
665 
666  do j=js,je+1
667  do i=isd,ied
668  cy(i,j,k) = cy(i,j,k) * frac
669  yfx(i,j,k) = yfx(i,j,k) * frac
670  enddo
671  enddo
672 
673  do j=js,je+1
674  do i=is,ie
675  mfy(i,j,k) = mfy(i,j,k) * frac
676  enddo
677  enddo
678  enddo
679  endif
680 
681 
682  do it=1,nsplt
683  if ( gridstruct%nested ) then
684  neststruct%tracer_nest_timestep = neststruct%tracer_nest_timestep + 1
685  end if
686  call timing_on('COMM_TOTAL')
687  call timing_on('COMM_TRACER')
688  call complete_group_halo_update(q_pack, domain)
689  call timing_off('COMM_TRACER')
690  call timing_off('COMM_TOTAL')
691 
692  if (gridstruct%nested) then
693  do iq=1,nq
694  call nested_grid_bc_apply_intt(q(isd:ied,jsd:jed,:,iq), &
695  0, 0, npx, npy, npz, bd, &
696  real(neststruct%tracer_nest_timestep)+real(nsplt*k_split), real(nsplt*k_split), &
697  neststruct%q_bc(iq), bctype=neststruct%nestbctype )
698  enddo
699  endif
700 
701 
702 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq, &
703 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
704 !$OMP private(dp2, ra_x, ra_y, fx, fy)
705  do k=1,npz
706 
707  do j=js,je
708  do i=is,ie
709  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)
710  enddo
711  enddo
712 
713  do j=jsd,jed
714  do i=is,ie
715  ra_x(i,j) = area(i,j) + (xfx(i,j,k) - xfx(i+1,j,k))
716  enddo
717  enddo
718  do j=js,je
719  do i=isd,ied
720  ra_y(i,j) = area(i,j) + (yfx(i,j,k) - yfx(i,j+1,k))
721  enddo
722  enddo
723 
724  do iq=1,nq
725  if ( it==1 .and. trdm>1.e-4 ) then
726  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
727  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
728  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k), &
729  mass=dp1(isd,jsd,k), nord=nord_tr, damp_c=trdm)
730  else
731  call fv_tp_2d(q(isd,jsd,k,iq), cx(is,jsd,k), cy(isd,js,k), &
732  npx, npy, hord, fx, fy, xfx(is,jsd,k), yfx(isd,js,k), &
733  gridstruct, bd, ra_x, ra_y, mfx=mfx(is,js,k), mfy=mfy(is,js,k))
734  endif
735  do j=js,je
736  do i=is,ie
737  q(i,j,k,iq) = ( q(i,j,k,iq)*dp1(i,j,k) + &
738  ((fx(i,j)-fx(i+1,j))+(fy(i,j)-fy(i,j+1)))*rarea(i,j) )/dp2(i,j)
739  enddo
740  enddo
741  enddo
742  enddo ! npz
743 
744  if ( it /= nsplt ) then
745  call timing_on('COMM_TOTAL')
746  call timing_on('COMM_TRACER')
747  call start_group_halo_update(q_pack, q, domain)
748  call timing_off('COMM_TRACER')
749  call timing_off('COMM_TOTAL')
750  endif
751  !Apply nested-grid BCs
752  if ( gridstruct%nested ) then
753  do iq=1,nq
754 
755 
756  call nested_grid_bc_apply_intt(q(isd:ied,jsd:jed,:,iq), &
757  0, 0, npx, npy, npz, bd, &
758  real(neststruct%tracer_nest_timestep), real(nsplt*k_split), &
759  neststruct%q_bc(iq), bctype=neststruct%nestbctype )
760 
761  end do
762  end if
763 
764 
765  enddo ! nsplt
766 
767  if ( id_divg > 0 ) then
768  rdt = 1./(frac*dt)
769 
770 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dp1,xfx,yfx,rarea,rdt)
771  do k=1,npz
772  do j=js,je
773  do i=is,ie
774  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
775  enddo
776  enddo
777  enddo
778  endif
779 
780  end subroutine tracer_2d_nested
781 
782 subroutine offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, &
783  gridstruct, flagstruct, bd, domain, &
784  ak, bk, ptop, npx, npy, npz, &
785  nq, hord, kord, q_split, k_split, dt, z_tracer, fill)
787  use fv_mapz_nlm_mod, only: map1_q2
788  use fv_fill_nlm_mod, only: fillz
789 
790  integer, intent(IN) :: npx
791  integer, intent(IN) :: npy
792  integer, intent(IN) :: npz
793  integer, intent(IN) :: nq ! number of tracers to be advected
794  integer, intent(IN) :: hord
795  integer, intent(IN) :: kord
796  integer, intent(IN) :: q_split
797  integer, intent(IN) :: k_split
798  logical, intent(IN) :: z_tracer
799  logical, intent(IN) :: fill
800  type(fv_grid_bounds_type), intent(IN ) :: bd
801  type(fv_flags_type), intent(INOUT) :: flagstruct
802  type(fv_grid_type), intent(IN), target :: gridstruct
803  type(domain2d), intent(INOUT) :: domain
804 
805  real, intent(IN ) :: dt
806  real, intent(IN ) ::ple0(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP before dyn_core
807  real, intent(INOUT) ::ple1(bd%is:bd%ie,bd%js:bd%je,npz+1) ! DELP after dyn_core
808  real, intent(IN ) :: cx(bd%is:bd%ie,bd%js:bd%je,npz) ! Courant Number X-Dir
809  real, intent(IN ) :: cy(bd%is:bd%ie,bd%js:bd%je,npz) ! Courant Number Y-Dir
810  real, intent(IN ) :: mfx(bd%is:bd%ie,bd%js:bd%je,npz) ! Mass Flux X-Dir
811  real, intent(IN ) :: mfy(bd%is:bd%ie,bd%js:bd%je,npz) ! Mass Flux Y-Dir
812  real, intent(INOUT) :: q(bd%is:bd%ie,bd%js:bd%je,npz,nq) ! Tracers
813  real, intent(IN ) :: ak(npz+1) ! AK for remapping
814  real, intent(IN ) :: bk(npz+1) ! BK for remapping
815  real, intent(IN ) :: ptop
816 ! Local Arrays
817  real :: xl(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) ! X-Dir for MPP Updates
818  real :: yl(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) ! Y-Dir for MPP Updates
819  real :: cxl(bd%is :bd%ie +1,bd%jsd:bd%jed ,npz) ! Courant Number X-Dir
820  real :: cyl(bd%isd:bd%ied ,bd%js :bd%je +1,npz) ! Courant Number Y-Dir
821  real :: mfxl(bd%is :bd%ie +1,bd%js :bd%je ,npz) ! Mass Flux X-Dir
822  real :: mfyl(bd%is :bd%ie ,bd%js :bd%je +1,npz) ! Mass Flux Y-Dir
823  real :: dpl(bd%is :bd%ie ,bd%js :bd%je ,npz) ! Pressure Thickness
824  real :: dpa(bd%is :bd%ie ,bd%js :bd%je ,npz) ! Pressure Thickness
825 ! Local Tracer Arrays
826  real :: q1(bd%is:bd%ie ,bd%js:bd%je, npz )! 2D Tracers
827  real :: q2(bd%isd:bd%ied ,bd%jsd:bd%jed ,nq)! 2D Tracers
828  real :: q3(bd%isd:bd%ied ,bd%jsd:bd%jed, npz,nq)! 3D Tracers
829 ! Local Buffer Arrarys
830  real :: wbuffer(bd%js:bd%je,npz)
831  real :: sbuffer(bd%is:bd%ie,npz)
832  real :: ebuffer(bd%js:bd%je,npz)
833  real :: nbuffer(bd%is:bd%ie,npz)
834 ! Local Remap Arrays
835  real pe1(bd%is:bd%ie,npz+1)
836  real pe2(bd%is:bd%ie,npz+1)
837  real dp2(bd%is:bd%ie,bd%js:bd%je,npz)
838 
839 ! Local indices
840  integer :: i,j,k,n,iq
841 
842  real :: scalingfactor
843 
844  type(group_halo_update_type), save :: i_pack
845 
846  integer :: is, ie, js, je
847  integer :: isd, ied, jsd, jed
848 
849  is = bd%is
850  ie = bd%ie
851  js = bd%js
852  je = bd%je
853  isd = bd%isd
854  ied = bd%ied
855  jsd = bd%jsd
856  jed = bd%jed
857 
858 ! Time-step
859 ! Fill CX/CY C-Grid boundaries and update ghost regions
860  xl(is:ie,js:je,:) = cx(:,:,:)
861  yl(is:ie,js:je,:) = cy(:,:,:)
862  call mpp_get_boundary(xl, yl, domain, &
863  wbufferx=wbuffer, ebufferx=ebuffer, &
864  sbuffery=sbuffer, nbuffery=nbuffer, &
865  gridtype=cgrid_ne )
866  xl(ie+1,js:je,:) = ebuffer
867  yl(is:ie,je+1,:) = nbuffer
868  call mpp_update_domains( xl, yl, domain, gridtype=cgrid_ne, complete=.true.)
869  cxl(is:ie+1,jsd:jed,:) = xl(is:ie+1,jsd:jed,:)
870  cyl(isd:ied,js:je+1,:) = yl(isd:ied,js:je+1,:)
871 
872 ! Fill MFX/MFY C-Grid boundaries
873  xl(is:ie,js:je,:) = mfx(:,:,:)
874  yl(is:ie,js:je,:) = mfy(:,:,:)
875  call mpp_get_boundary(xl, yl, domain, &
876  wbufferx=wbuffer, ebufferx=ebuffer, &
877  sbuffery=sbuffer, nbuffery=nbuffer, &
878  gridtype=cgrid_ne )
879  xl(ie+1,js:je,:) = ebuffer
880  yl(is:ie,je+1,:) = nbuffer
881  mfxl(is:ie+1,js:je,:) = xl(is:ie+1,js:je,:)
882  mfyl(is:ie,js:je+1,:) = yl(is:ie,js:je+1,:)
883 
884 ! Fill local tracers and pressure thickness
885  dpl(:,:,:) = ple0(:,:,2:npz+1) - ple0(:,:,1:npz)
886  q3(is:ie,js:je,:,:) = q(is:ie,js:je,:,:)
887 
888  if ( z_tracer ) then
889 !$omp parallel do default(shared) private(q2)
890  do k=1,npz
891  do iq=1,nq
892  do j=js,je
893  do i=is,ie ! To_do list:
894  q2(i,j,iq) = q3(i,j,k,iq) ! The data copying can be avoided if q is
895  ! re-dimensioned as q(i,j,nq,k)
896  enddo
897  enddo
898  enddo
899  call start_group_halo_update(i_pack, q2, domain)
900  call tracer_2d_1l(q2, dpl(is,js,k), mfxl(is,js,k), mfyl(is,js,k), cxl(is,js,k), cyl(is,js,k), &
901  gridstruct, bd, domain, npx, npy, npz, nq, &
902  flagstruct%hord_tr, q_split, dt, 0, i_pack, &
903  flagstruct%nord_tr, flagstruct%trdm2, dpa=dpa)
904  do iq=1,nq
905  do j=js,je
906  do i=is,ie ! To_do list:
907  q3(i,j,k,iq) = q2(i,j,iq) ! The data copying can be avoided if q is
908  ! re-dimensioned as q(i,j,nq,k)
909  enddo
910  enddo
911  enddo
912  enddo
913  else
914  call start_group_halo_update(i_pack, q3, domain)
915  call tracer_2d(q3, dpl, mfxl, mfyl, cxl, cyl, gridstruct, bd, domain, npx, npy, npz, nq, &
916  flagstruct%hord_tr, q_split, dt, 0, i_pack, &
917  flagstruct%nord_tr, flagstruct%trdm2, dpa=dpa)
918  endif
919 
920 !------------------------------------------------------------------
921 ! Re-Map constituents
922 ! Do remapping one tracer at a time; seems to be faster
923 ! It requires less memory than mapn_ppm
924 !------------------------------------------------------------------
925 
926  do iq=1,nq
927  do j=js,je
928  ! pressures mapping from (dpA is new delp after tracer_2d)
929  pe1(:,1) = ptop
930  do k=2,npz+1
931  pe1(:,k) = pe1(:,k-1) + dpa(:,j,k-1)
932  enddo
933  ! pressures mapping to
934  pe2(:,1) = ptop
935  pe2(:,npz+1) = pe1(:,npz+1)
936  do k=2,npz
937  pe2(: ,k) = ak(k) + bk(k)*pe1(:,npz+1)
938  enddo
939  do k=1,npz
940  dp2(:,j,k) = pe2(:,k+1) - pe2(:,k)
941  enddo
942  call map1_q2(npz, pe1, q3(isd,jsd,1,iq), &
943  npz, pe2, q1(:,j,:), dp2(:,j,:), &
944  is, ie, 0, kord, j, &
945  isd, ied, jsd, jed, 0.)
946  if (fill) call fillz(ie-is+1, npz, 1, q1(:,j,:), dp2(:,j,:))
947  enddo
948  ! Rescale tracers based on ple1 at destination timestep
949  !------------------------------------------------------
950 
951  scalingfactor = calcscalingfactor(q1, dp2, ple1, npx, npy, npz, gridstruct, bd)
952  !scalingFactors = computeScalingFactors(q1, dp2, ple1, npx, npy, npz)
953 
954  ! Return tracers
955  !---------------
956  q(is:ie,js:je,1:npz,iq) = q1(is:ie,js:je,1:npz) * scalingfactor
957  !do k =1,npz
958  !do j = js,je
959  !do i = is,ie
960  !q(i,j,k,iq) = q1(i,j,k)
961  !enddo
962  !enddo
963  !enddo
964 
965  enddo
966 
967 end subroutine offline_tracer_advection
968 
969 !------------------------------------------------------------------------------------
970 
971  function calcscalingfactor(q1, dp2, ple1, npx, npy, npz, gridstruct, bd) result(scaling)
972  use mpp_mod, only: mpp_sum
973  integer, intent(in) :: npx
974  integer, intent(in) :: npy
975  integer, intent(in) :: npz
976  real, intent(in) :: q1(:,:,:)
977  real, intent(in) :: dp2(:,:,:)
978  real, intent(in) :: ple1(:,:,:)
979  type(fv_grid_type), intent(IN ) :: gridstruct
980  type(fv_grid_bounds_type), intent(IN ) :: bd
981  real :: scaling
982 
983  integer :: k
984  real :: partialsums(2,npz), globalsums(2)
985  real, parameter :: tiny_denominator = tiny(1.0)
986 
987  !-------
988  ! Compute partial sum on local array first to minimize communication.
989  ! This algorithm will not be strongly repdroducible under changes do domain
990  ! decomposition, but uses far less communication bandwidth (and memory BW)
991  ! then the preceding implementation.
992  !-------
993  do k = 1, npz
994  ! numerator
995  partialsums(1,k) = sum(q1(:,:,k)*dp2(:,:,k)*gridstruct%area(bd%is:bd%ie,bd%js:bd%je))
996  ! denominator
997  partialsums(2,k) = sum(q1(:,:,k)*(ple1(:,:,k+1)-ple1(:,:,k))*gridstruct%area(bd%is:bd%ie,bd%js:bd%je))
998  end do
999 
1000  globalsums(1) = sum(partialsums(1,:))
1001  globalsums(2) = sum(partialsums(2,:))
1002 
1003  call mpp_sum(globalsums, 2)
1004 
1005  if (globalsums(2) > tiny_denominator) then
1006  scaling = globalsums(1) / globalsums(2)
1007  !#################################################################
1008  ! This line was added to ensure strong reproducibility of the code
1009  !#################################################################
1010  scaling = REAL(scaling, kind=kind(1.00))
1011  else
1012  scaling = 1.d0
1013  end if
1014 
1015  end function calcscalingfactor
1016 
1017 end module fv_tracer2d_nlm_mod
real, dimension(:,:,:), allocatable nest_fx_east_accum
subroutine, public nested_grid_bc_apply_intt(var_nest, istag, jstag, npx, npy, npz, bd, step, split, BC, bctype)
real, dimension(:,:,:), allocatable nest_fx_west_accum
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, dpA)
real function calcscalingfactor(q1, dp2, ple1, npx, npy, npz, gridstruct, bd)
subroutine, public fillz(im, km, nq, q, dp)
Definition: fv_fill_nlm.F90:33
Definition: mpp.F90:39
real, dimension(:,:,:), allocatable nest_fx_north_accum
real, dimension(:,:,:), allocatable nest_fx_south_accum
subroutine, public map1_q2(km, pe1, q1, kn, pe2, q2, dp2, i1, i2, iv, kord, j, ibeg, iend, jbeg, jend, q_min)
integer, parameter, public ng
subroutine timing_on(blk_name)
#define max(a, b)
Definition: mosaic_util.h:33
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
subroutine, public offline_tracer_advection(q, ple0, ple1, mfx, mfy, cx, cy, gridstruct, flagstruct, bd, domain, ak, bk, ptop, npx, npy, npz, nq, hord, kord, q_split, k_split, dt, z_tracer, fill)
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, dpA)
subroutine, public fv_tp_2d(q, crx, cry, npx, npy, hord, fx, fy, xfx, yfx, gridstruct, bd, ra_x, ra_y, mfx, mfy, mass, nord, damp_c)
Definition: tp_core_nlm.F90:80
subroutine timing_off(blk_name)
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)