FV3 Bundle
fv_tracer2d_tlm.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 !***********************************************************************
23  use fv_mp_nlm_mod, only: mp_reduce_max
24  use fv_mp_nlm_mod, only: ng, mp_gather, is_master
25  use fv_mp_nlm_mod, only: group_halo_update_type
35 
36 implicit none
37 private
38 
41 
43 
44 CONTAINS
45 ! Differentiation of tracer_2d_1l in forward (tangent) mode:
46 ! variations of useful results: q dp1
47 ! with respect to varying inputs: q dp1 mfx mfy cx cy
48 !-----------------------------------------------------------------------
49 ! !ROUTINE: Perform 2D horizontal-to-lagrangian transport
50 !-----------------------------------------------------------------------
51  SUBROUTINE tracer_2d_1l_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, &
52 & mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz&
53 & , nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, &
54 & nord_tr_pert, trdm_pert, split_damp_tr, dpa)
55  IMPLICIT NONE
56  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
57  INTEGER, INTENT(IN) :: npx
58  INTEGER, INTENT(IN) :: npy
59  INTEGER, INTENT(IN) :: npz
60 ! number of tracers to be advected
61  INTEGER, INTENT(IN) :: nq
62  INTEGER, INTENT(IN) :: hord, nord_tr
63  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
64  LOGICAL, INTENT(IN) :: split_damp_tr
65  INTEGER, INTENT(IN) :: q_split
66  INTEGER, INTENT(IN) :: id_divg
67  REAL, INTENT(IN) :: dt, trdm, trdm_pert
68  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
69 ! Tracers
70  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
71  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
72 ! DELP before dyn_core
73  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
74  REAL, INTENT(INOUT) :: dp1_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
75 ! Mass Flux X-Dir
76  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
77  REAL, INTENT(INOUT) :: mfx_tl(bd%is:bd%ie+1, bd%js:bd%je, npz)
78 ! Mass Flux Y-Dir
79  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
80  REAL, INTENT(INOUT) :: mfy_tl(bd%is:bd%ie, bd%js:bd%je+1, npz)
81 ! Courant Number X-Dir
82  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
83  REAL, INTENT(INOUT) :: cx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
84 ! Courant Number Y-Dir
85  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
86  REAL, INTENT(INOUT) :: cy_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
87 ! DELP after advection
88  REAL, OPTIONAL, INTENT(OUT) :: dpa(bd%is:bd%ie, bd%js:bd%je)
89  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
90  TYPE(domain2d), INTENT(INOUT) :: domain
91 ! Local Arrays
92 ! 3D tracers
93  REAL :: qn2(bd%isd:bd%ied, bd%jsd:bd%jed, nq)
94  REAL :: qn2_tl(bd%isd:bd%ied, bd%jsd:bd%jed, nq)
95  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
96  REAL :: dp2_tl(bd%is:bd%ie, bd%js:bd%je)
97  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
98  REAL :: fx_tl(bd%is:bd%ie+1, bd%js:bd%je)
99  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
100  REAL :: fy_tl(bd%is:bd%ie, bd%js:bd%je+1)
101  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
102  REAL :: ra_x_tl(bd%is:bd%ie, bd%jsd:bd%jed)
103  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
104  REAL :: ra_y_tl(bd%isd:bd%ied, bd%js:bd%je)
105  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
106  REAL :: xfx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
107  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
108  REAL :: yfx_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
109  REAL :: cmax(npz)
110  REAL :: frac
111  INTEGER :: nsplt
112  INTEGER :: i, j, k, it, iq
113  REAL, DIMENSION(:, :), POINTER :: area, rarea
114  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
115  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
116  INTEGER :: is, ie, js, je
117  INTEGER :: isd, ied, jsd, jed
118  INTRINSIC abs
119  INTRINSIC max
120  INTRINSIC int
121  INTRINSIC real
122  INTRINSIC PRESENT
123  REAL :: max1
124  REAL :: x1
125  REAL :: z1
126  REAL :: y3
127  REAL :: y2
128  REAL :: y1
129  is = bd%is
130  ie = bd%ie
131  js = bd%js
132  je = bd%je
133  isd = bd%isd
134  ied = bd%ied
135  jsd = bd%jsd
136  jed = bd%jed
137  area => gridstruct%area
138  rarea => gridstruct%rarea
139  sin_sg => gridstruct%sin_sg
140  dxa => gridstruct%dxa
141  dya => gridstruct%dya
142  dx => gridstruct%dx
143  dy => gridstruct%dy
144  xfx_tl = 0.0
145  yfx_tl = 0.0
146 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
147 !$OMP sin_sg,cy,yfx,dya,dx,cmax)
148  DO k=1,npz
149  DO j=jsd,jed
150  DO i=is,ie+1
151  IF (cx(i, j, k) .GT. 0.) THEN
152  xfx_tl(i, j, k) = dxa(i-1, j)*dy(i, j)*sin_sg(i-1, j, 3)*&
153 & cx_tl(i, j, k)
154  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
155 & j, 3)
156  ELSE
157  xfx_tl(i, j, k) = dxa(i, j)*dy(i, j)*sin_sg(i, j, 1)*cx_tl(i&
158 & , j, k)
159  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
160 & )
161  END IF
162  END DO
163  END DO
164  DO j=js,je+1
165  DO i=isd,ied
166  IF (cy(i, j, k) .GT. 0.) THEN
167  yfx_tl(i, j, k) = dya(i, j-1)*dx(i, j)*sin_sg(i, j-1, 4)*&
168 & cy_tl(i, j, k)
169  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
170 & 1, 4)
171  ELSE
172  yfx_tl(i, j, k) = dya(i, j)*dx(i, j)*sin_sg(i, j, 2)*cy_tl(i&
173 & , j, k)
174  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
175 & )
176  END IF
177  END DO
178  END DO
179  cmax(k) = 0.
180  IF (k .LT. npz/6) THEN
181  DO j=js,je
182  DO i=is,ie
183  IF (cx(i, j, k) .GE. 0.) THEN
184  y1 = cx(i, j, k)
185  ELSE
186  y1 = -cx(i, j, k)
187  END IF
188  IF (cy(i, j, k) .GE. 0.) THEN
189  z1 = cy(i, j, k)
190  ELSE
191  z1 = -cy(i, j, k)
192  END IF
193  IF (cmax(k) .LT. y1) THEN
194  IF (y1 .LT. z1) THEN
195  cmax(k) = z1
196  ELSE
197  cmax(k) = y1
198  END IF
199  ELSE IF (cmax(k) .LT. z1) THEN
200  cmax(k) = z1
201  ELSE
202  cmax(k) = cmax(k)
203  END IF
204  END DO
205  END DO
206  ELSE
207  DO j=js,je
208  DO i=is,ie
209  IF (cx(i, j, k) .GE. 0.) THEN
210  x1 = cx(i, j, k)
211  ELSE
212  x1 = -cx(i, j, k)
213  END IF
214  IF (cy(i, j, k) .GE. 0.) THEN
215  y3 = cy(i, j, k)
216  ELSE
217  y3 = -cy(i, j, k)
218  END IF
219  IF (x1 .LT. y3) THEN
220  max1 = y3
221  ELSE
222  max1 = x1
223  END IF
224  y2 = max1 + 1. - sin_sg(i, j, 5)
225  IF (cmax(k) .LT. y2) THEN
226  cmax(k) = y2
227  ELSE
228  cmax(k) = cmax(k)
229  END IF
230  END DO
231  END DO
232  END IF
233  END DO
234 ! k-loop
235  CALL mp_reduce_max(cmax, npz)
236 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx, &
237 !$OMP cy,yfx,mfx,mfy,cmax) &
238 !$OMP private(nsplt, frac)
239  DO k=1,npz
240  nsplt = int(1. + cmax(k))
241  IF (nsplt .GT. 1) THEN
242  frac = 1./REAL(nsplt)
243  DO j=jsd,jed
244  DO i=is,ie+1
245  cx_tl(i, j, k) = frac*cx_tl(i, j, k)
246  cx(i, j, k) = cx(i, j, k)*frac
247  xfx_tl(i, j, k) = frac*xfx_tl(i, j, k)
248  xfx(i, j, k) = xfx(i, j, k)*frac
249  END DO
250  END DO
251  DO j=js,je
252  DO i=is,ie+1
253  mfx_tl(i, j, k) = frac*mfx_tl(i, j, k)
254  mfx(i, j, k) = mfx(i, j, k)*frac
255  END DO
256  END DO
257  DO j=js,je+1
258  DO i=isd,ied
259  cy_tl(i, j, k) = frac*cy_tl(i, j, k)
260  cy(i, j, k) = cy(i, j, k)*frac
261  yfx_tl(i, j, k) = frac*yfx_tl(i, j, k)
262  yfx(i, j, k) = yfx(i, j, k)*frac
263  END DO
264  END DO
265  DO j=js,je+1
266  DO i=is,ie
267  mfy_tl(i, j, k) = frac*mfy_tl(i, j, k)
268  mfy(i, j, k) = mfy(i, j, k)*frac
269  END DO
270  END DO
271  END IF
272  END DO
273  CALL timing_on('COMM_TOTAL')
274  CALL timing_on('COMM_TRACER')
275  CALL complete_group_halo_update(q_pack, domain)
276  CALL timing_off('COMM_TRACER')
277  CALL timing_off('COMM_TOTAL')
278  dp2_tl = 0.0
279  qn2_tl = 0.0
280  ra_x_tl = 0.0
281  ra_y_tl = 0.0
282  fx_tl = 0.0
283  fy_tl = 0.0
284 ! Begin k-independent tracer transport; can not be OpenMPed because the mpp_update call.
285  DO k=1,npz
286 !$OMP parallel do default(none) shared(k,is,ie,js,je,isd,ied,jsd,jed,xfx,area,yfx,ra_x,ra_y)
287  DO j=jsd,jed
288  DO i=is,ie
289  ra_x_tl(i, j) = xfx_tl(i, j, k) - xfx_tl(i+1, j, k)
290  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
291  END DO
292  IF (j .GE. js .AND. j .LE. je) THEN
293  DO i=isd,ied
294  ra_y_tl(i, j) = yfx_tl(i, j, k) - yfx_tl(i, j+1, k)
295  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
296  END DO
297  END IF
298  END DO
299  nsplt = int(1. + cmax(k))
300  DO it=1,nsplt
301 !$OMP parallel do default(none) shared(k,is,ie,js,je,rarea,mfx,mfy,dp1,dp2)
302  DO j=js,je
303  DO i=is,ie
304  dp2_tl(i, j) = dp1_tl(i, j, k) + rarea(i, j)*(mfx_tl(i, j, k&
305 & )-mfx_tl(i+1, j, k)+mfy_tl(i, j, k)-mfy_tl(i, j+1, k))
306  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(mfy&
307 & (i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
308  END DO
309  END DO
310 !$OMP parallel do default(none) shared(k,nsplt,it,is,ie,js,je,isd,ied,jsd,jed,npx,npy,cx,xfx,hord,trdm, &
311 !$OMP nord_tr,nq,gridstruct,bd,cy,yfx,mfx,mfy,qn2,q,ra_x,ra_y,dp1,dp2,rarea) &
312 !$OMP private(fx,fy)
313  DO iq=1,nq
314  IF (nsplt .NE. 1) THEN
315  IF (it .EQ. 1) THEN
316  DO j=jsd,jed
317  DO i=isd,ied
318  qn2_tl(i, j, iq) = q_tl(i, j, k, iq)
319  qn2(i, j, iq) = q(i, j, k, iq)
320  END DO
321  END DO
322  END IF
323  IF (hord .EQ. hord_pert) THEN
324  CALL fv_tp_2d_tlm(qn2(isd:ied, jsd:jed, iq), qn2_tl(isd&
325 & :ied, jsd:jed, iq), cx(is:ie+1, jsd:jed, k)&
326 & , cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, &
327 & js:je+1, k), cy_tl(isd:ied, js:je+1, k), &
328 & npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx(&
329 & is:ie+1, jsd:jed, k), xfx_tl(is:ie+1, jsd:&
330 & jed, k), yfx(isd:ied, js:je+1, k), yfx_tl(&
331 & isd:ied, js:je+1, k), gridstruct, bd, ra_x&
332 & , ra_x_tl, ra_y, ra_y_tl, mfx=mfx(is:ie+1, &
333 & js:je, k), mfx_tl=mfx_tl(is:ie+1, js:je, k)&
334 & , mfy=mfy(is:ie, js:je+1, k), mfy_tl=mfy_tl&
335 & (is:ie, js:je+1, k))
336  ELSE
337  CALL fv_tp_2d_tlm(qn2(isd:ied, jsd:jed, iq), qn2_tl(isd:&
338 & ied, jsd:jed, iq), cx(is:ie+1, jsd:jed, k), &
339 & cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, js:je+&
340 & 1, k), cy_tl(isd:ied, js:je+1, k), npx, npy, &
341 & hord_pert, fx, fx_tl, fy, fy_tl, xfx(is:ie+1, &
342 & jsd:jed, k), xfx_tl(is:ie+1, jsd:jed, k), yfx(&
343 & isd:ied, js:je+1, k), yfx_tl(isd:ied, js:je+1&
344 & , k), gridstruct, bd, ra_x, ra_x_tl, ra_y, &
345 & ra_y_tl, mfx=mfx(is:ie+1, js:je, k), mfx_tl=&
346 & mfx_tl(is:ie+1, js:je, k), mfy=mfy(is:ie, js:&
347 & je+1, k), mfy_tl=mfy_tl(is:ie, js:je+1, k))
348  call fv_tp_2d(qn2(isd:ied,jsd:jed,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
349  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
350  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
351  END IF
352  IF (it .LT. nsplt) THEN
353 ! not last call
354  DO j=js,je
355  DO i=is,ie
356  qn2_tl(i, j, iq) = ((qn2_tl(i, j, iq)*dp1(i, j, k)+qn2&
357 & (i, j, iq)*dp1_tl(i, j, k)+rarea(i, j)*(fx_tl(i, j)-&
358 & fx_tl(i+1, j)+fy_tl(i, j)-fy_tl(i, j+1)))*dp2(i, j)-&
359 & (qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)-fx(i+1, j)+(fy&
360 & (i, j)-fy(i, j+1)))*rarea(i, j))*dp2_tl(i, j))/dp2(i&
361 & , j)**2
362  qn2(i, j, iq) = (qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)-&
363 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i&
364 & , j)
365  END DO
366  END DO
367  ELSE
368  DO j=js,je
369  DO i=is,ie
370  q_tl(i, j, k, iq) = ((qn2_tl(i, j, iq)*dp1(i, j, k)+&
371 & qn2(i, j, iq)*dp1_tl(i, j, k)+rarea(i, j)*(fx_tl(i, &
372 & j)-fx_tl(i+1, j)+fy_tl(i, j)-fy_tl(i, j+1)))*dp2(i, &
373 & j)-(qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)-fx(i+1, j)+&
374 & (fy(i, j)-fy(i, j+1)))*rarea(i, j))*dp2_tl(i, j))/&
375 & dp2(i, j)**2
376  q(i, j, k, iq) = (qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)&
377 & -fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(&
378 & i, j)
379  END DO
380  END DO
381  END IF
382  ELSE
383  IF (hord .EQ. hord_pert) THEN
384  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:&
385 & ied, jsd:jed, k, iq), cx(is:ie+1, jsd:jed, &
386 & k), cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied&
387 & , js:je+1, k), cy_tl(isd:ied, js:je+1, k), &
388 & npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx(&
389 & is:ie+1, jsd:jed, k), xfx_tl(is:ie+1, jsd:&
390 & jed, k), yfx(isd:ied, js:je+1, k), yfx_tl(&
391 & isd:ied, js:je+1, k), gridstruct, bd, ra_x&
392 & , ra_x_tl, ra_y, ra_y_tl, mfx=mfx(is:ie+1, &
393 & js:je, k), mfx_tl=mfx_tl(is:ie+1, js:je, k)&
394 & , mfy=mfy(is:ie, js:je+1, k), mfy_tl=mfy_tl&
395 & (is:ie, js:je+1, k))
396  ELSE
397  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:ied&
398 & , jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k), &
399 & cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, js:je+&
400 & 1, k), cy_tl(isd:ied, js:je+1, k), npx, npy, &
401 & hord_pert, fx, fx_tl, fy, fy_tl, xfx(is:ie+1, &
402 & jsd:jed, k), xfx_tl(is:ie+1, jsd:jed, k), yfx(&
403 & isd:ied, js:je+1, k), yfx_tl(isd:ied, js:je+1&
404 & , k), gridstruct, bd, ra_x, ra_x_tl, ra_y, &
405 & ra_y_tl, mfx=mfx(is:ie+1, js:je, k), mfx_tl=&
406 & mfx_tl(is:ie+1, js:je, k), mfy=mfy(is:ie, js:&
407 & je+1, k), mfy_tl=mfy_tl(is:ie, js:je+1, k))
408  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
409  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
410  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
411  END IF
412  DO j=js,je
413  DO i=is,ie
414  q_tl(i, j, k, iq) = ((q_tl(i, j, k, iq)*dp1(i, j, k)+q(i&
415 & , j, k, iq)*dp1_tl(i, j, k)+rarea(i, j)*(fx_tl(i, j)-&
416 & fx_tl(i+1, j)+fy_tl(i, j)-fy_tl(i, j+1)))*dp2(i, j)-(q&
417 & (i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-fx(i+1, j)+(fy(i&
418 & , j)-fy(i, j+1)))*rarea(i, j))*dp2_tl(i, j))/dp2(i, j)&
419 & **2
420  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-&
421 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, &
422 & j)
423  END DO
424  END DO
425  END IF
426  END DO
427 ! tracer-loop
428  IF (it .LT. nsplt) THEN
429 ! not last call
430  DO j=js,je
431  DO i=is,ie
432  dp1_tl(i, j, k) = dp2_tl(i, j)
433  dp1(i, j, k) = dp2(i, j)
434  END DO
435  END DO
436  CALL timing_on('COMM_TOTAL')
437  CALL timing_on('COMM_TRACER')
438  CALL mpp_update_domains_tlm(qn2, qn2_tl, domain)
439  CALL timing_off('COMM_TRACER')
440  CALL timing_off('COMM_TOTAL')
441  END IF
442  END DO
443  END DO
444 ! time-split loop
445 ! k-loop
446  IF (PRESENT(dpa)) dpa = dp2
447  END SUBROUTINE tracer_2d_1l_tlm
448 !-----------------------------------------------------------------------
449 ! !ROUTINE: Perform 2D horizontal-to-lagrangian transport
450 !-----------------------------------------------------------------------
451  SUBROUTINE tracer_2d_1l(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, &
452 & domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, &
453 & nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, &
454 & dpa)
455  IMPLICIT NONE
456  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
457  INTEGER, INTENT(IN) :: npx
458  INTEGER, INTENT(IN) :: npy
459  INTEGER, INTENT(IN) :: npz
460 ! number of tracers to be advected
461  INTEGER, INTENT(IN) :: nq
462  INTEGER, INTENT(IN) :: hord, nord_tr
463  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
464  LOGICAL, INTENT(IN) :: split_damp_tr
465  INTEGER, INTENT(IN) :: q_split
466  INTEGER, INTENT(IN) :: id_divg
467  REAL, INTENT(IN) :: dt, trdm, trdm_pert
468  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
469 ! Tracers
470  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
471 ! DELP before dyn_core
472  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
473 ! Mass Flux X-Dir
474  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
475 ! Mass Flux Y-Dir
476  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
477 ! Courant Number X-Dir
478  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
479 ! Courant Number Y-Dir
480  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
481 ! DELP after advection
482  REAL, OPTIONAL, INTENT(OUT) :: dpa(bd%is:bd%ie, bd%js:bd%je)
483  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
484  TYPE(domain2d), INTENT(INOUT) :: domain
485 ! Local Arrays
486 ! 3D tracers
487  REAL :: qn2(bd%isd:bd%ied, bd%jsd:bd%jed, nq)
488  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
489  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
490  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
491  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
492  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
493  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
494  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
495  REAL :: cmax(npz)
496  REAL :: frac
497  INTEGER :: nsplt
498  INTEGER :: i, j, k, it, iq
499  REAL, DIMENSION(:, :), POINTER :: area, rarea
500  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
501  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
502  INTEGER :: is, ie, js, je
503  INTEGER :: isd, ied, jsd, jed
504  INTRINSIC abs
505  INTRINSIC max
506  INTRINSIC int
507  INTRINSIC real
508  INTRINSIC PRESENT
509  REAL :: max1
510  REAL :: x1
511  REAL :: z1
512  REAL :: y3
513  REAL :: y2
514  REAL :: y1
515  is = bd%is
516  ie = bd%ie
517  js = bd%js
518  je = bd%je
519  isd = bd%isd
520  ied = bd%ied
521  jsd = bd%jsd
522  jed = bd%jed
523  area => gridstruct%area
524  rarea => gridstruct%rarea
525  sin_sg => gridstruct%sin_sg
526  dxa => gridstruct%dxa
527  dya => gridstruct%dya
528  dx => gridstruct%dx
529  dy => gridstruct%dy
530 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
531 !$OMP sin_sg,cy,yfx,dya,dx,cmax)
532  DO k=1,npz
533  DO j=jsd,jed
534  DO i=is,ie+1
535  IF (cx(i, j, k) .GT. 0.) THEN
536  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
537 & j, 3)
538  ELSE
539  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
540 & )
541  END IF
542  END DO
543  END DO
544  DO j=js,je+1
545  DO i=isd,ied
546  IF (cy(i, j, k) .GT. 0.) THEN
547  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
548 & 1, 4)
549  ELSE
550  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
551 & )
552  END IF
553  END DO
554  END DO
555  cmax(k) = 0.
556  IF (k .LT. npz/6) THEN
557  DO j=js,je
558  DO i=is,ie
559  IF (cx(i, j, k) .GE. 0.) THEN
560  y1 = cx(i, j, k)
561  ELSE
562  y1 = -cx(i, j, k)
563  END IF
564  IF (cy(i, j, k) .GE. 0.) THEN
565  z1 = cy(i, j, k)
566  ELSE
567  z1 = -cy(i, j, k)
568  END IF
569  IF (cmax(k) .LT. y1) THEN
570  IF (y1 .LT. z1) THEN
571  cmax(k) = z1
572  ELSE
573  cmax(k) = y1
574  END IF
575  ELSE IF (cmax(k) .LT. z1) THEN
576  cmax(k) = z1
577  ELSE
578  cmax(k) = cmax(k)
579  END IF
580  END DO
581  END DO
582  ELSE
583  DO j=js,je
584  DO i=is,ie
585  IF (cx(i, j, k) .GE. 0.) THEN
586  x1 = cx(i, j, k)
587  ELSE
588  x1 = -cx(i, j, k)
589  END IF
590  IF (cy(i, j, k) .GE. 0.) THEN
591  y3 = cy(i, j, k)
592  ELSE
593  y3 = -cy(i, j, k)
594  END IF
595  IF (x1 .LT. y3) THEN
596  max1 = y3
597  ELSE
598  max1 = x1
599  END IF
600  y2 = max1 + 1. - sin_sg(i, j, 5)
601  IF (cmax(k) .LT. y2) THEN
602  cmax(k) = y2
603  ELSE
604  cmax(k) = cmax(k)
605  END IF
606  END DO
607  END DO
608  END IF
609  END DO
610 ! k-loop
611  CALL mp_reduce_max(cmax, npz)
612 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx, &
613 !$OMP cy,yfx,mfx,mfy,cmax) &
614 !$OMP private(nsplt, frac)
615  DO k=1,npz
616  nsplt = int(1. + cmax(k))
617  IF (nsplt .GT. 1) THEN
618  frac = 1./REAL(nsplt)
619  DO j=jsd,jed
620  DO i=is,ie+1
621  cx(i, j, k) = cx(i, j, k)*frac
622  xfx(i, j, k) = xfx(i, j, k)*frac
623  END DO
624  END DO
625  DO j=js,je
626  DO i=is,ie+1
627  mfx(i, j, k) = mfx(i, j, k)*frac
628  END DO
629  END DO
630  DO j=js,je+1
631  DO i=isd,ied
632  cy(i, j, k) = cy(i, j, k)*frac
633  yfx(i, j, k) = yfx(i, j, k)*frac
634  END DO
635  END DO
636  DO j=js,je+1
637  DO i=is,ie
638  mfy(i, j, k) = mfy(i, j, k)*frac
639  END DO
640  END DO
641  END IF
642  END DO
643  CALL timing_on('COMM_TOTAL')
644  CALL timing_on('COMM_TRACER')
645  CALL complete_group_halo_update(q_pack, domain)
646  CALL timing_off('COMM_TRACER')
647  CALL timing_off('COMM_TOTAL')
648 ! Begin k-independent tracer transport; can not be OpenMPed because the mpp_update call.
649  DO k=1,npz
650 !$OMP parallel do default(none) shared(k,is,ie,js,je,isd,ied,jsd,jed,xfx,area,yfx,ra_x,ra_y)
651  DO j=jsd,jed
652  DO i=is,ie
653  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
654  END DO
655  IF (j .GE. js .AND. j .LE. je) THEN
656  DO i=isd,ied
657  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
658  END DO
659  END IF
660  END DO
661  nsplt = int(1. + cmax(k))
662  DO it=1,nsplt
663 !$OMP parallel do default(none) shared(k,is,ie,js,je,rarea,mfx,mfy,dp1,dp2)
664  DO j=js,je
665  DO i=is,ie
666  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(mfy&
667 & (i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
668  END DO
669  END DO
670 !$OMP parallel do default(none) shared(k,nsplt,it,is,ie,js,je,isd,ied,jsd,jed,npx,npy,cx,xfx,hord,trdm, &
671 !$OMP nord_tr,nq,gridstruct,bd,cy,yfx,mfx,mfy,qn2,q,ra_x,ra_y,dp1,dp2,rarea) &
672 !$OMP private(fx,fy)
673  DO iq=1,nq
674  IF (nsplt .NE. 1) THEN
675  IF (it .EQ. 1) THEN
676  DO j=jsd,jed
677  DO i=isd,ied
678  qn2(i, j, iq) = q(i, j, k, iq)
679  END DO
680  END DO
681  END IF
682  IF (hord .EQ. hord_pert) THEN
683  CALL fv_tp_2d(qn2(isd:ied, jsd:jed, iq), cx(is:ie+1, &
684 & jsd:jed, k), cy(isd:ied, js:je+1, k), npx, npy&
685 & , hord, fx, fy, xfx(is:ie+1, jsd:jed, k), yfx(&
686 & isd:ied, js:je+1, k), gridstruct, bd, ra_x, &
687 & ra_y, mfx=mfx(is:ie+1, js:je, k), mfy=mfy(is:ie&
688 & , js:je+1, k))
689  ELSE
690  call fv_tp_2d(qn2(isd:ied,jsd:jed,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
691  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
692  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
693  END IF
694  IF (it .LT. nsplt) THEN
695 ! not last call
696  DO j=js,je
697  DO i=is,ie
698  qn2(i, j, iq) = (qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)-&
699 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i&
700 & , j)
701  END DO
702  END DO
703  ELSE
704  DO j=js,je
705  DO i=is,ie
706  q(i, j, k, iq) = (qn2(i, j, iq)*dp1(i, j, k)+(fx(i, j)&
707 & -fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(&
708 & i, j)
709  END DO
710  END DO
711  END IF
712  ELSE
713  IF (hord .EQ. hord_pert) THEN
714  CALL fv_tp_2d(q(isd:ied, jsd:jed, k, iq), cx(is:ie+1, &
715 & jsd:jed, k), cy(isd:ied, js:je+1, k), npx, npy&
716 & , hord, fx, fy, xfx(is:ie+1, jsd:jed, k), yfx(&
717 & isd:ied, js:je+1, k), gridstruct, bd, ra_x, &
718 & ra_y, mfx=mfx(is:ie+1, js:je, k), mfy=mfy(is:ie&
719 & , js:je+1, k))
720  ELSE
721  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
722  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
723  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
724  END IF
725  DO j=js,je
726  DO i=is,ie
727  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-&
728 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, &
729 & j)
730  END DO
731  END DO
732  END IF
733  END DO
734 ! tracer-loop
735  IF (it .LT. nsplt) THEN
736 ! not last call
737  DO j=js,je
738  DO i=is,ie
739  dp1(i, j, k) = dp2(i, j)
740  END DO
741  END DO
742  CALL timing_on('COMM_TOTAL')
743  CALL timing_on('COMM_TRACER')
744  CALL mpp_update_domains(qn2, domain)
745  CALL timing_off('COMM_TRACER')
746  CALL timing_off('COMM_TOTAL')
747  END IF
748  END DO
749  END DO
750 ! time-split loop
751 ! k-loop
752  IF (PRESENT(dpa)) dpa = dp2
753  END SUBROUTINE tracer_2d_1l
754 ! Differentiation of tracer_2d in forward (tangent) mode:
755 ! variations of useful results: q dp1
756 ! with respect to varying inputs: q dp1 mfx mfy cx cy
757  SUBROUTINE tracer_2d_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, &
758 & mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz&
759 & , nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, &
760 & nord_tr_pert, trdm_pert, split_damp_tr, dpa)
761  IMPLICIT NONE
762  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
763  INTEGER, INTENT(IN) :: npx
764  INTEGER, INTENT(IN) :: npy
765  INTEGER, INTENT(IN) :: npz
766 ! number of tracers to be advected
767  INTEGER, INTENT(IN) :: nq
768  INTEGER, INTENT(IN) :: hord, nord_tr
769  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
770  LOGICAL, INTENT(IN) :: split_damp_tr
771  INTEGER, INTENT(IN) :: q_split
772  INTEGER, INTENT(IN) :: id_divg
773  REAL, INTENT(IN) :: dt, trdm, trdm_pert
774  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
775 ! Tracers
776  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
777  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
778 ! DELP before dyn_core
779  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
780  REAL, INTENT(INOUT) :: dp1_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
781 ! Mass Flux X-Dir
782  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
783  REAL, INTENT(INOUT) :: mfx_tl(bd%is:bd%ie+1, bd%js:bd%je, npz)
784 ! Mass Flux Y-Dir
785  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
786  REAL, INTENT(INOUT) :: mfy_tl(bd%is:bd%ie, bd%js:bd%je+1, npz)
787 ! Courant Number X-Dir
788  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
789  REAL, INTENT(INOUT) :: cx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
790 ! Courant Number Y-Dir
791  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
792  REAL, INTENT(INOUT) :: cy_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
793 ! DELP after advection
794  REAL, OPTIONAL, INTENT(OUT) :: dpa(bd%is:bd%ie, bd%js:bd%je, npz)
795  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
796  TYPE(domain2d), INTENT(INOUT) :: domain
797 ! Local Arrays
798  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
799  REAL :: dp2_tl(bd%is:bd%ie, bd%js:bd%je)
800  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
801  REAL :: fx_tl(bd%is:bd%ie+1, bd%js:bd%je)
802  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
803  REAL :: fy_tl(bd%is:bd%ie, bd%js:bd%je+1)
804  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
805  REAL :: ra_x_tl(bd%is:bd%ie, bd%jsd:bd%jed)
806  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
807  REAL :: ra_y_tl(bd%isd:bd%ied, bd%js:bd%je)
808  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
809  REAL :: xfx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
810  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
811  REAL :: yfx_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
812  REAL :: cmax(npz)
813  REAL :: c_global
814  REAL :: frac, rdt
815  INTEGER :: ksplt(npz)
816  INTEGER :: nsplt
817  INTEGER :: i, j, k, it, iq
818  REAL, DIMENSION(:, :), POINTER :: area, rarea
819  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
820  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
821  INTEGER :: is, ie, js, je
822  INTEGER :: isd, ied, jsd, jed
823  INTRINSIC abs
824  INTRINSIC max
825  INTRINSIC int
826  INTRINSIC real
827  INTRINSIC PRESENT
828  REAL :: max1
829  REAL :: x1
830  REAL :: z1
831  REAL :: y3
832  REAL :: y2
833  REAL :: y1
834  is = bd%is
835  ie = bd%ie
836  js = bd%js
837  je = bd%je
838  isd = bd%isd
839  ied = bd%ied
840  jsd = bd%jsd
841  jed = bd%jed
842  area => gridstruct%area
843  rarea => gridstruct%rarea
844  sin_sg => gridstruct%sin_sg
845  dxa => gridstruct%dxa
846  dya => gridstruct%dya
847  dx => gridstruct%dx
848  dy => gridstruct%dy
849  xfx_tl = 0.0
850  yfx_tl = 0.0
851 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
852 !$OMP sin_sg,cy,yfx,dya,dx,cmax,q_split,ksplt)
853  DO k=1,npz
854  DO j=jsd,jed
855  DO i=is,ie+1
856  IF (cx(i, j, k) .GT. 0.) THEN
857  xfx_tl(i, j, k) = dxa(i-1, j)*dy(i, j)*sin_sg(i-1, j, 3)*&
858 & cx_tl(i, j, k)
859  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
860 & j, 3)
861  ELSE
862  xfx_tl(i, j, k) = dxa(i, j)*dy(i, j)*sin_sg(i, j, 1)*cx_tl(i&
863 & , j, k)
864  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
865 & )
866  END IF
867  END DO
868  END DO
869  DO j=js,je+1
870  DO i=isd,ied
871  IF (cy(i, j, k) .GT. 0.) THEN
872  yfx_tl(i, j, k) = dya(i, j-1)*dx(i, j)*sin_sg(i, j-1, 4)*&
873 & cy_tl(i, j, k)
874  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
875 & 1, 4)
876  ELSE
877  yfx_tl(i, j, k) = dya(i, j)*dx(i, j)*sin_sg(i, j, 2)*cy_tl(i&
878 & , j, k)
879  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
880 & )
881  END IF
882  END DO
883  END DO
884  IF (q_split .EQ. 0) THEN
885  cmax(k) = 0.
886  IF (k .LT. npz/6) THEN
887  DO j=js,je
888  DO i=is,ie
889  IF (cx(i, j, k) .GE. 0.) THEN
890  y1 = cx(i, j, k)
891  ELSE
892  y1 = -cx(i, j, k)
893  END IF
894  IF (cy(i, j, k) .GE. 0.) THEN
895  z1 = cy(i, j, k)
896  ELSE
897  z1 = -cy(i, j, k)
898  END IF
899  IF (cmax(k) .LT. y1) THEN
900  IF (y1 .LT. z1) THEN
901  cmax(k) = z1
902  ELSE
903  cmax(k) = y1
904  END IF
905  ELSE IF (cmax(k) .LT. z1) THEN
906  cmax(k) = z1
907  ELSE
908  cmax(k) = cmax(k)
909  END IF
910  END DO
911  END DO
912  ELSE
913  DO j=js,je
914  DO i=is,ie
915  IF (cx(i, j, k) .GE. 0.) THEN
916  x1 = cx(i, j, k)
917  ELSE
918  x1 = -cx(i, j, k)
919  END IF
920  IF (cy(i, j, k) .GE. 0.) THEN
921  y3 = cy(i, j, k)
922  ELSE
923  y3 = -cy(i, j, k)
924  END IF
925  IF (x1 .LT. y3) THEN
926  max1 = y3
927  ELSE
928  max1 = x1
929  END IF
930  y2 = max1 + 1. - sin_sg(i, j, 5)
931  IF (cmax(k) .LT. y2) THEN
932  cmax(k) = y2
933  ELSE
934  cmax(k) = cmax(k)
935  END IF
936  END DO
937  END DO
938  END IF
939  END IF
940  ksplt(k) = 1
941  END DO
942 !--------------------------------------------------------------------------------
943 ! Determine global nsplt:
944  IF (q_split .EQ. 0) THEN
945  CALL mp_reduce_max(cmax, npz)
946 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
947  c_global = cmax(1)
948  IF (npz .NE. 1) THEN
949 ! if NOT shallow water test case
950  DO k=2,npz
951  IF (cmax(k) .LT. c_global) THEN
952  c_global = c_global
953  ELSE
954  c_global = cmax(k)
955  END IF
956  END DO
957  END IF
958  nsplt = int(1. + c_global)
959  IF (is_master() .AND. nsplt .GT. 4) WRITE(*, *) 'Tracer_2d_split='&
960 & , nsplt, c_global
961  ELSE
962  nsplt = q_split
963  END IF
964 !--------------------------------------------------------------------------------
965  IF (nsplt .NE. 1) THEN
966 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,mfx,cy,yfx,mfy,cmax,nsplt,ksplt) &
967 !$OMP private( frac )
968  DO k=1,npz
969  ksplt(k) = int(1. + cmax(k))
970  frac = 1./REAL(ksplt(k))
971  DO j=jsd,jed
972  DO i=is,ie+1
973  cx_tl(i, j, k) = frac*cx_tl(i, j, k)
974  cx(i, j, k) = cx(i, j, k)*frac
975  xfx_tl(i, j, k) = frac*xfx_tl(i, j, k)
976  xfx(i, j, k) = xfx(i, j, k)*frac
977  END DO
978  END DO
979  DO j=js,je
980  DO i=is,ie+1
981  mfx_tl(i, j, k) = frac*mfx_tl(i, j, k)
982  mfx(i, j, k) = mfx(i, j, k)*frac
983  END DO
984  END DO
985  DO j=js,je+1
986  DO i=isd,ied
987  cy_tl(i, j, k) = frac*cy_tl(i, j, k)
988  cy(i, j, k) = cy(i, j, k)*frac
989  yfx_tl(i, j, k) = frac*yfx_tl(i, j, k)
990  yfx(i, j, k) = yfx(i, j, k)*frac
991  END DO
992  END DO
993  DO j=js,je+1
994  DO i=is,ie
995  mfy_tl(i, j, k) = frac*mfy_tl(i, j, k)
996  mfy(i, j, k) = mfy(i, j, k)*frac
997  END DO
998  END DO
999  END DO
1000  dp2_tl = 0.0
1001  ra_x_tl = 0.0
1002  ra_y_tl = 0.0
1003  fx_tl = 0.0
1004  fy_tl = 0.0
1005  ELSE
1006  dp2_tl = 0.0
1007  ra_x_tl = 0.0
1008  ra_y_tl = 0.0
1009  fx_tl = 0.0
1010  fy_tl = 0.0
1011  END IF
1012  DO it=1,nsplt
1013  CALL timing_on('COMM_TOTAL')
1014  CALL timing_on('COMM_TRACER')
1015  CALL complete_group_halo_update(q_pack, domain)
1016  CALL timing_off('COMM_TRACER')
1017  CALL timing_off('COMM_TOTAL')
1018 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq,ksplt,&
1019 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
1020 !$OMP private(dp2, ra_x, ra_y, fx, fy)
1021  DO k=1,npz
1022 ! ksplt
1023  IF (it .LE. ksplt(k)) THEN
1024  DO j=js,je
1025  DO i=is,ie
1026  dp2_tl(i, j) = dp1_tl(i, j, k) + rarea(i, j)*(mfx_tl(i, j&
1027 & , k)-mfx_tl(i+1, j, k)+mfy_tl(i, j, k)-mfy_tl(i, j+1, k)&
1028 & )
1029  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(&
1030 & mfy(i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
1031  END DO
1032  END DO
1033  DO j=jsd,jed
1034  DO i=is,ie
1035  ra_x_tl(i, j) = xfx_tl(i, j, k) - xfx_tl(i+1, j, k)
1036  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
1037  END DO
1038  END DO
1039  DO j=js,je
1040  DO i=isd,ied
1041  ra_y_tl(i, j) = yfx_tl(i, j, k) - yfx_tl(i, j+1, k)
1042  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
1043  END DO
1044  END DO
1045  DO iq=1,nq
1046  IF (it .EQ. 1 .AND. trdm .GT. 1.e-4) THEN
1047  IF (hord .EQ. hord_pert) THEN
1048  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(&
1049 & isd:ied, jsd:jed, k, iq), cx(is:ie+1, jsd&
1050 & :jed, k), cx_tl(is:ie+1, jsd:jed, k), cy(&
1051 & isd:ied, js:je+1, k), cy_tl(isd:ied, js:&
1052 & je+1, k), npx, npy, hord, fx, fx_tl, fy, &
1053 & fy_tl, xfx(is:ie+1, jsd:jed, k), xfx_tl(&
1054 & is:ie+1, jsd:jed, k), yfx(isd:ied, js:je+&
1055 & 1, k), yfx_tl(isd:ied, js:je+1, k), &
1056 & gridstruct, bd, ra_x, ra_x_tl, ra_y, &
1057 & ra_y_tl, mfx=mfx(is:ie+1, js:je, k), &
1058 & mfx_tl=mfx_tl(is:ie+1, js:je, k), mfy=mfy&
1059 & (is:ie, js:je+1, k), mfy_tl=mfy_tl(is:ie&
1060 & , js:je+1, k), mass=dp1(isd:ied, jsd:jed&
1061 & , k), mass_tl=dp1_tl(isd:ied, jsd:jed, k)&
1062 & , nord=nord_tr, damp_c=trdm)
1063  ELSE
1064  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:&
1065 & ied, jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k&
1066 & ), cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, &
1067 & js:je+1, k), cy_tl(isd:ied, js:je+1, k), npx&
1068 & , npy, hord_pert, fx, fx_tl, fy, fy_tl, xfx(&
1069 & is:ie+1, jsd:jed, k), xfx_tl(is:ie+1, jsd:&
1070 & jed, k), yfx(isd:ied, js:je+1, k), yfx_tl(&
1071 & isd:ied, js:je+1, k), gridstruct, bd, ra_x, &
1072 & ra_x_tl, ra_y, ra_y_tl, mfx=mfx(is:ie+1, js:&
1073 & je, k), mfx_tl=mfx_tl(is:ie+1, js:je, k), &
1074 & mfy=mfy(is:ie, js:je+1, k), mfy_tl=mfy_tl(is&
1075 & :ie, js:je+1, k), mass=dp1(isd:ied, jsd:jed&
1076 & , k), mass_tl=dp1_tl(isd:ied, jsd:jed, k), &
1077 & nord=nord_tr_pert, damp_c=trdm_pert)
1078  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1079  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1080  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k), &
1081  mass=dp1(isd:ied,jsd:jed,k), nord=nord_tr, damp_c=trdm)
1082  END IF
1083  ELSE IF (hord .EQ. hord_pert) THEN
1084  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:&
1085 & ied, jsd:jed, k, iq), cx(is:ie+1, jsd:jed, &
1086 & k), cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied&
1087 & , js:je+1, k), cy_tl(isd:ied, js:je+1, k), &
1088 & npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx(&
1089 & is:ie+1, jsd:jed, k), xfx_tl(is:ie+1, jsd:&
1090 & jed, k), yfx(isd:ied, js:je+1, k), yfx_tl(&
1091 & isd:ied, js:je+1, k), gridstruct, bd, ra_x&
1092 & , ra_x_tl, ra_y, ra_y_tl, mfx=mfx(is:ie+1, &
1093 & js:je, k), mfx_tl=mfx_tl(is:ie+1, js:je, k)&
1094 & , mfy=mfy(is:ie, js:je+1, k), mfy_tl=mfy_tl&
1095 & (is:ie, js:je+1, k))
1096  ELSE
1097  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:ied&
1098 & , jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k), &
1099 & cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, js:je+&
1100 & 1, k), cy_tl(isd:ied, js:je+1, k), npx, npy, &
1101 & hord_pert, fx, fx_tl, fy, fy_tl, xfx(is:ie+1, &
1102 & jsd:jed, k), xfx_tl(is:ie+1, jsd:jed, k), yfx(&
1103 & isd:ied, js:je+1, k), yfx_tl(isd:ied, js:je+1&
1104 & , k), gridstruct, bd, ra_x, ra_x_tl, ra_y, &
1105 & ra_y_tl, mfx=mfx(is:ie+1, js:je, k), mfx_tl=&
1106 & mfx_tl(is:ie+1, js:je, k), mfy=mfy(is:ie, js:&
1107 & je+1, k), mfy_tl=mfy_tl(is:ie, js:je+1, k))
1108  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1109  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1110  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
1111  END IF
1112  DO j=js,je
1113  DO i=is,ie
1114  q_tl(i, j, k, iq) = ((q_tl(i, j, k, iq)*dp1(i, j, k)+q(i&
1115 & , j, k, iq)*dp1_tl(i, j, k)+rarea(i, j)*(fx_tl(i, j)-&
1116 & fx_tl(i+1, j)+fy_tl(i, j)-fy_tl(i, j+1)))*dp2(i, j)-(q&
1117 & (i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-fx(i+1, j)+(fy(i&
1118 & , j)-fy(i, j+1)))*rarea(i, j))*dp2_tl(i, j))/dp2(i, j)&
1119 & **2
1120  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-&
1121 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, &
1122 & j)
1123  END DO
1124  END DO
1125  END DO
1126  IF (it .NE. nsplt) THEN
1127  DO j=js,je
1128  DO i=is,ie
1129  dp1_tl(i, j, k) = dp2_tl(i, j)
1130  dp1(i, j, k) = dp2(i, j)
1131  END DO
1132  END DO
1133  END IF
1134  END IF
1135  END DO
1136 ! npz
1137  IF (it .NE. nsplt) THEN
1138  CALL timing_on('COMM_TOTAL')
1139  CALL timing_on('COMM_TRACER')
1140  CALL start_group_halo_update_tlm(q_pack, q, q_tl, domain)
1141  CALL timing_off('COMM_TRACER')
1142  CALL timing_off('COMM_TOTAL')
1143  END IF
1144  END DO
1145 ! nsplt
1146  IF (PRESENT(dpa)) dpa = dp1(bd%is:bd%ie, bd%js:bd%je, 1:npz)
1147  END SUBROUTINE tracer_2d_tlm
1148  SUBROUTINE tracer_2d(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, domain&
1149 & , npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, &
1150 & trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
1151  IMPLICIT NONE
1152  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
1153  INTEGER, INTENT(IN) :: npx
1154  INTEGER, INTENT(IN) :: npy
1155  INTEGER, INTENT(IN) :: npz
1156 ! number of tracers to be advected
1157  INTEGER, INTENT(IN) :: nq
1158  INTEGER, INTENT(IN) :: hord, nord_tr
1159  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
1160  LOGICAL, INTENT(IN) :: split_damp_tr
1161  INTEGER, INTENT(IN) :: q_split
1162  INTEGER, INTENT(IN) :: id_divg
1163  REAL, INTENT(IN) :: dt, trdm, trdm_pert
1164  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
1165 ! Tracers
1166  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
1167 ! DELP before dyn_core
1168  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1169 ! Mass Flux X-Dir
1170  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
1171 ! Mass Flux Y-Dir
1172  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
1173 ! Courant Number X-Dir
1174  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1175 ! Courant Number Y-Dir
1176  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1177 ! DELP after advection
1178  REAL, OPTIONAL, INTENT(OUT) :: dpa(bd%is:bd%ie, bd%js:bd%je, npz)
1179  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
1180  TYPE(domain2d), INTENT(INOUT) :: domain
1181 ! Local Arrays
1182  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
1183  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
1184  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
1185  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
1186  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
1187  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1188  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1189  REAL :: cmax(npz)
1190  REAL :: c_global
1191  REAL :: frac, rdt
1192  INTEGER :: ksplt(npz)
1193  INTEGER :: nsplt
1194  INTEGER :: i, j, k, it, iq
1195  REAL, DIMENSION(:, :), POINTER :: area, rarea
1196  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
1197  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
1198  INTEGER :: is, ie, js, je
1199  INTEGER :: isd, ied, jsd, jed
1200  INTRINSIC abs
1201  INTRINSIC max
1202  INTRINSIC int
1203  INTRINSIC real
1204  INTRINSIC PRESENT
1205  REAL :: max1
1206  REAL :: x1
1207  REAL :: z1
1208  REAL :: y3
1209  REAL :: y2
1210  REAL :: y1
1211  is = bd%is
1212  ie = bd%ie
1213  js = bd%js
1214  je = bd%je
1215  isd = bd%isd
1216  ied = bd%ied
1217  jsd = bd%jsd
1218  jed = bd%jed
1219  area => gridstruct%area
1220  rarea => gridstruct%rarea
1221  sin_sg => gridstruct%sin_sg
1222  dxa => gridstruct%dxa
1223  dya => gridstruct%dya
1224  dx => gridstruct%dx
1225  dy => gridstruct%dy
1226 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
1227 !$OMP sin_sg,cy,yfx,dya,dx,cmax,q_split,ksplt)
1228  DO k=1,npz
1229  DO j=jsd,jed
1230  DO i=is,ie+1
1231  IF (cx(i, j, k) .GT. 0.) THEN
1232  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
1233 & j, 3)
1234  ELSE
1235  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
1236 & )
1237  END IF
1238  END DO
1239  END DO
1240  DO j=js,je+1
1241  DO i=isd,ied
1242  IF (cy(i, j, k) .GT. 0.) THEN
1243  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
1244 & 1, 4)
1245  ELSE
1246  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
1247 & )
1248  END IF
1249  END DO
1250  END DO
1251  IF (q_split .EQ. 0) THEN
1252  cmax(k) = 0.
1253  IF (k .LT. npz/6) THEN
1254  DO j=js,je
1255  DO i=is,ie
1256  IF (cx(i, j, k) .GE. 0.) THEN
1257  y1 = cx(i, j, k)
1258  ELSE
1259  y1 = -cx(i, j, k)
1260  END IF
1261  IF (cy(i, j, k) .GE. 0.) THEN
1262  z1 = cy(i, j, k)
1263  ELSE
1264  z1 = -cy(i, j, k)
1265  END IF
1266  IF (cmax(k) .LT. y1) THEN
1267  IF (y1 .LT. z1) THEN
1268  cmax(k) = z1
1269  ELSE
1270  cmax(k) = y1
1271  END IF
1272  ELSE IF (cmax(k) .LT. z1) THEN
1273  cmax(k) = z1
1274  ELSE
1275  cmax(k) = cmax(k)
1276  END IF
1277  END DO
1278  END DO
1279  ELSE
1280  DO j=js,je
1281  DO i=is,ie
1282  IF (cx(i, j, k) .GE. 0.) THEN
1283  x1 = cx(i, j, k)
1284  ELSE
1285  x1 = -cx(i, j, k)
1286  END IF
1287  IF (cy(i, j, k) .GE. 0.) THEN
1288  y3 = cy(i, j, k)
1289  ELSE
1290  y3 = -cy(i, j, k)
1291  END IF
1292  IF (x1 .LT. y3) THEN
1293  max1 = y3
1294  ELSE
1295  max1 = x1
1296  END IF
1297  y2 = max1 + 1. - sin_sg(i, j, 5)
1298  IF (cmax(k) .LT. y2) THEN
1299  cmax(k) = y2
1300  ELSE
1301  cmax(k) = cmax(k)
1302  END IF
1303  END DO
1304  END DO
1305  END IF
1306  END IF
1307  ksplt(k) = 1
1308  END DO
1309 !--------------------------------------------------------------------------------
1310 ! Determine global nsplt:
1311  IF (q_split .EQ. 0) THEN
1312  CALL mp_reduce_max(cmax, npz)
1313 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
1314  c_global = cmax(1)
1315  IF (npz .NE. 1) THEN
1316 ! if NOT shallow water test case
1317  DO k=2,npz
1318  IF (cmax(k) .LT. c_global) THEN
1319  c_global = c_global
1320  ELSE
1321  c_global = cmax(k)
1322  END IF
1323  END DO
1324  END IF
1325  nsplt = int(1. + c_global)
1326  IF (is_master() .AND. nsplt .GT. 4) WRITE(*, *) 'Tracer_2d_split='&
1327 & , nsplt, c_global
1328  ELSE
1329  nsplt = q_split
1330  END IF
1331 !--------------------------------------------------------------------------------
1332  IF (nsplt .NE. 1) THEN
1333 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,mfx,cy,yfx,mfy,cmax,nsplt,ksplt) &
1334 !$OMP private( frac )
1335  DO k=1,npz
1336  ksplt(k) = int(1. + cmax(k))
1337  frac = 1./REAL(ksplt(k))
1338  DO j=jsd,jed
1339  DO i=is,ie+1
1340  cx(i, j, k) = cx(i, j, k)*frac
1341  xfx(i, j, k) = xfx(i, j, k)*frac
1342  END DO
1343  END DO
1344  DO j=js,je
1345  DO i=is,ie+1
1346  mfx(i, j, k) = mfx(i, j, k)*frac
1347  END DO
1348  END DO
1349  DO j=js,je+1
1350  DO i=isd,ied
1351  cy(i, j, k) = cy(i, j, k)*frac
1352  yfx(i, j, k) = yfx(i, j, k)*frac
1353  END DO
1354  END DO
1355  DO j=js,je+1
1356  DO i=is,ie
1357  mfy(i, j, k) = mfy(i, j, k)*frac
1358  END DO
1359  END DO
1360  END DO
1361  END IF
1362  DO it=1,nsplt
1363  CALL timing_on('COMM_TOTAL')
1364  CALL timing_on('COMM_TRACER')
1365  CALL complete_group_halo_update(q_pack, domain)
1366  CALL timing_off('COMM_TRACER')
1367  CALL timing_off('COMM_TOTAL')
1368 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq,ksplt,&
1369 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
1370 !$OMP private(dp2, ra_x, ra_y, fx, fy)
1371  DO k=1,npz
1372 ! ksplt
1373  IF (it .LE. ksplt(k)) THEN
1374  DO j=js,je
1375  DO i=is,ie
1376  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(&
1377 & mfy(i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
1378  END DO
1379  END DO
1380  DO j=jsd,jed
1381  DO i=is,ie
1382  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
1383  END DO
1384  END DO
1385  DO j=js,je
1386  DO i=isd,ied
1387  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
1388  END DO
1389  END DO
1390  DO iq=1,nq
1391  IF (it .EQ. 1 .AND. trdm .GT. 1.e-4) THEN
1392  IF (hord .EQ. hord_pert) THEN
1393  CALL fv_tp_2d(q(isd:ied, jsd:jed, k, iq), cx(is:ie+1&
1394 & , jsd:jed, k), cy(isd:ied, js:je+1, k), npx, &
1395 & npy, hord, fx, fy, xfx(is:ie+1, jsd:jed, k), &
1396 & yfx(isd:ied, js:je+1, k), gridstruct, bd, &
1397 & ra_x, ra_y, mfx(is:ie+1, js:je, k), mfy(is:ie&
1398 & , js:je+1, k), dp1(isd:ied, jsd:jed, k), &
1399 & nord_tr, trdm)
1400  ELSE
1401  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1402  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1403  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k), &
1404  mass=dp1(isd:ied,jsd:jed,k), nord=nord_tr, damp_c=trdm)
1405  END IF
1406  ELSE IF (hord .EQ. hord_pert) THEN
1407  CALL fv_tp_2d(q(isd:ied, jsd:jed, k, iq), cx(is:ie+1, &
1408 & jsd:jed, k), cy(isd:ied, js:je+1, k), npx, npy&
1409 & , hord, fx, fy, xfx(is:ie+1, jsd:jed, k), yfx(&
1410 & isd:ied, js:je+1, k), gridstruct, bd, ra_x, &
1411 & ra_y, mfx=mfx(is:ie+1, js:je, k), mfy=mfy(is:ie&
1412 & , js:je+1, k))
1413  ELSE
1414  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1415  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1416  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
1417  END IF
1418  DO j=js,je
1419  DO i=is,ie
1420  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-&
1421 & fx(i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, &
1422 & j)
1423  END DO
1424  END DO
1425  END DO
1426  IF (it .NE. nsplt) THEN
1427  DO j=js,je
1428  DO i=is,ie
1429  dp1(i, j, k) = dp2(i, j)
1430  END DO
1431  END DO
1432  END IF
1433  END IF
1434  END DO
1435 ! npz
1436  IF (it .NE. nsplt) THEN
1437  CALL timing_on('COMM_TOTAL')
1438  CALL timing_on('COMM_TRACER')
1439  CALL start_group_halo_update(q_pack, q, domain)
1440  CALL timing_off('COMM_TRACER')
1441  CALL timing_off('COMM_TOTAL')
1442  END IF
1443  END DO
1444 ! nsplt
1445  IF (PRESENT(dpa)) dpa = dp1(bd%is:bd%ie, bd%js:bd%je, 1:npz)
1446  END SUBROUTINE tracer_2d
1447 ! Differentiation of tracer_2d_nested in forward (tangent) mode:
1448 ! variations of useful results: q dp1
1449 ! with respect to varying inputs: q dp1 mfx mfy cx cy
1450  SUBROUTINE tracer_2d_nested_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy&
1451 & , mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, &
1452 & npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, k_split&
1453 & , neststruct, parent_grid, hord_pert, nord_tr_pert, trdm_pert, &
1454 & split_damp_tr)
1455  IMPLICIT NONE
1456  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
1457  INTEGER, INTENT(IN) :: npx
1458  INTEGER, INTENT(IN) :: npy
1459  INTEGER, INTENT(IN) :: npz
1460 ! number of tracers to be advected
1461  INTEGER, INTENT(IN) :: nq
1462  INTEGER, INTENT(IN) :: hord, nord_tr
1463  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
1464  LOGICAL, INTENT(IN) :: split_damp_tr
1465  INTEGER, INTENT(IN) :: q_split, k_split
1466  INTEGER, INTENT(IN) :: id_divg
1467  REAL, INTENT(IN) :: dt, trdm, trdm_pert
1468  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
1469 ! Tracers
1470  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
1471  REAL, INTENT(INOUT) :: q_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
1472 ! DELP before dyn_core
1473  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1474  REAL, INTENT(INOUT) :: dp1_tl(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1475 ! Mass Flux X-Dir
1476  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
1477  REAL, INTENT(INOUT) :: mfx_tl(bd%is:bd%ie+1, bd%js:bd%je, npz)
1478 ! Mass Flux Y-Dir
1479  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
1480  REAL, INTENT(INOUT) :: mfy_tl(bd%is:bd%ie, bd%js:bd%je+1, npz)
1481 ! Courant Number X-Dir
1482  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1483  REAL, INTENT(INOUT) :: cx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1484 ! Courant Number Y-Dir
1485  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1486  REAL, INTENT(INOUT) :: cy_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1487  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
1488  TYPE(fv_nest_type), INTENT(INOUT) :: neststruct
1489  TYPE(fv_atmos_type), INTENT(INOUT) :: parent_grid
1490  TYPE(domain2d), INTENT(INOUT) :: domain
1491 ! Local Arrays
1492  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
1493  REAL :: dp2_tl(bd%is:bd%ie, bd%js:bd%je)
1494  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
1495  REAL :: fx_tl(bd%is:bd%ie+1, bd%js:bd%je)
1496  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
1497  REAL :: fy_tl(bd%is:bd%ie, bd%js:bd%je+1)
1498  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
1499  REAL :: ra_x_tl(bd%is:bd%ie, bd%jsd:bd%jed)
1500  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
1501  REAL :: ra_y_tl(bd%isd:bd%ied, bd%js:bd%je)
1502  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1503  REAL :: xfx_tl(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1504  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1505  REAL :: yfx_tl(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1506  REAL :: cmax(npz)
1507  REAL :: cmax_t
1508  REAL :: c_global
1509  REAL :: frac, rdt
1510  INTEGER :: nsplt, nsplt_parent
1511  INTEGER, SAVE :: msg_split_steps=1
1512  INTEGER :: i, j, k, it, iq
1513  REAL, DIMENSION(:, :), POINTER :: area, rarea
1514  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
1515  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
1516  INTEGER :: is, ie, js, je
1517  INTEGER :: isd, ied, jsd, jed
1518  INTRINSIC abs
1519  INTRINSIC max
1520  INTRINSIC int
1521  INTRINSIC real
1522  REAL :: max1
1523  REAL :: x2
1524  REAL :: x1
1525  REAL :: y2
1526  REAL :: y1
1527  is = bd%is
1528  ie = bd%ie
1529  js = bd%js
1530  je = bd%je
1531  isd = bd%isd
1532  ied = bd%ied
1533  jsd = bd%jsd
1534  jed = bd%jed
1535  area => gridstruct%area
1536  rarea => gridstruct%rarea
1537  sin_sg => gridstruct%sin_sg
1538  dxa => gridstruct%dxa
1539  dya => gridstruct%dya
1540  dx => gridstruct%dx
1541  dy => gridstruct%dy
1542  xfx_tl = 0.0
1543  yfx_tl = 0.0
1544 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
1545 !$OMP sin_sg,cy,yfx,dya,dx)
1546  DO k=1,npz
1547  DO j=jsd,jed
1548  DO i=is,ie+1
1549  IF (cx(i, j, k) .GT. 0.) THEN
1550  xfx_tl(i, j, k) = dxa(i-1, j)*dy(i, j)*sin_sg(i-1, j, 3)*&
1551 & cx_tl(i, j, k)
1552  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
1553 & j, 3)
1554  ELSE
1555  xfx_tl(i, j, k) = dxa(i, j)*dy(i, j)*sin_sg(i, j, 1)*cx_tl(i&
1556 & , j, k)
1557  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
1558 & )
1559  END IF
1560  END DO
1561  END DO
1562  DO j=js,je+1
1563  DO i=isd,ied
1564  IF (cy(i, j, k) .GT. 0.) THEN
1565  yfx_tl(i, j, k) = dya(i, j-1)*dx(i, j)*sin_sg(i, j-1, 4)*&
1566 & cy_tl(i, j, k)
1567  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
1568 & 1, 4)
1569  ELSE
1570  yfx_tl(i, j, k) = dya(i, j)*dx(i, j)*sin_sg(i, j, 2)*cy_tl(i&
1571 & , j, k)
1572  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
1573 & )
1574  END IF
1575  END DO
1576  END DO
1577  END DO
1578 !--------------------------------------------------------------------------------
1579  IF (q_split .EQ. 0) THEN
1580 ! Determine nsplt
1581 !$OMP parallel do default(none) shared(is,ie,js,je,npz,cmax,cx,cy,sin_sg) &
1582 !$OMP private(cmax_t )
1583  DO k=1,npz
1584  cmax(k) = 0.
1585  IF (k .LT. 4) THEN
1586 ! Top layers: C < max( abs(c_x), abs(c_y) )
1587  DO j=js,je
1588  DO i=is,ie
1589  IF (cx(i, j, k) .GE. 0.) THEN
1590  x1 = cx(i, j, k)
1591  ELSE
1592  x1 = -cx(i, j, k)
1593  END IF
1594  IF (cy(i, j, k) .GE. 0.) THEN
1595  y1 = cy(i, j, k)
1596  ELSE
1597  y1 = -cy(i, j, k)
1598  END IF
1599  IF (x1 .LT. y1) THEN
1600  cmax_t = y1
1601  ELSE
1602  cmax_t = x1
1603  END IF
1604  IF (cmax_t .LT. cmax(k)) THEN
1605  cmax(k) = cmax(k)
1606  ELSE
1607  cmax(k) = cmax_t
1608  END IF
1609  END DO
1610  END DO
1611  ELSE
1612  DO j=js,je
1613  DO i=is,ie
1614  IF (cx(i, j, k) .GE. 0.) THEN
1615  x2 = cx(i, j, k)
1616  ELSE
1617  x2 = -cx(i, j, k)
1618  END IF
1619  IF (cy(i, j, k) .GE. 0.) THEN
1620  y2 = cy(i, j, k)
1621  ELSE
1622  y2 = -cy(i, j, k)
1623  END IF
1624  IF (x2 .LT. y2) THEN
1625  max1 = y2
1626  ELSE
1627  max1 = x2
1628  END IF
1629  cmax_t = max1 + 1. - sin_sg(i, j, 5)
1630  IF (cmax_t .LT. cmax(k)) THEN
1631  cmax(k) = cmax(k)
1632  ELSE
1633  cmax(k) = cmax_t
1634  END IF
1635  END DO
1636  END DO
1637  END IF
1638  END DO
1639  CALL mp_reduce_max(cmax, npz)
1640 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
1641  c_global = cmax(1)
1642  IF (npz .NE. 1) THEN
1643 ! if NOT shallow water test case
1644  DO k=2,npz
1645  IF (cmax(k) .LT. c_global) THEN
1646  c_global = c_global
1647  ELSE
1648  c_global = cmax(k)
1649  END IF
1650  END DO
1651  END IF
1652  nsplt = int(1. + c_global)
1653  IF (is_master() .AND. nsplt .GT. 3) WRITE(*, *) 'Tracer_2d_split='&
1654 & , nsplt, c_global
1655  ELSE
1656  nsplt = q_split
1657  IF (gridstruct%nested .AND. neststruct%nestbctype .GT. 1) THEN
1658  IF (q_split/parent_grid%flagstruct%q_split .LT. 1) THEN
1659  msg_split_steps = 1
1660  ELSE
1661  msg_split_steps = q_split/parent_grid%flagstruct%q_split
1662  END IF
1663  END IF
1664  END IF
1665 !--------------------------------------------------------------------------------
1666  frac = 1./REAL(nsplt)
1667  IF (nsplt .NE. 1) THEN
1668 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,frac,xfx,mfx,cy,yfx,mfy)
1669  DO k=1,npz
1670  DO j=jsd,jed
1671  DO i=is,ie+1
1672  cx_tl(i, j, k) = frac*cx_tl(i, j, k)
1673  cx(i, j, k) = cx(i, j, k)*frac
1674  xfx_tl(i, j, k) = frac*xfx_tl(i, j, k)
1675  xfx(i, j, k) = xfx(i, j, k)*frac
1676  END DO
1677  END DO
1678  DO j=js,je
1679  DO i=is,ie+1
1680  mfx_tl(i, j, k) = frac*mfx_tl(i, j, k)
1681  mfx(i, j, k) = mfx(i, j, k)*frac
1682  END DO
1683  END DO
1684  DO j=js,je+1
1685  DO i=isd,ied
1686  cy_tl(i, j, k) = frac*cy_tl(i, j, k)
1687  cy(i, j, k) = cy(i, j, k)*frac
1688  yfx_tl(i, j, k) = frac*yfx_tl(i, j, k)
1689  yfx(i, j, k) = yfx(i, j, k)*frac
1690  END DO
1691  END DO
1692  DO j=js,je+1
1693  DO i=is,ie
1694  mfy_tl(i, j, k) = frac*mfy_tl(i, j, k)
1695  mfy(i, j, k) = mfy(i, j, k)*frac
1696  END DO
1697  END DO
1698  END DO
1699  dp2_tl = 0.0
1700  ra_x_tl = 0.0
1701  ra_y_tl = 0.0
1702  fx_tl = 0.0
1703  fy_tl = 0.0
1704  ELSE
1705  dp2_tl = 0.0
1706  ra_x_tl = 0.0
1707  ra_y_tl = 0.0
1708  fx_tl = 0.0
1709  fy_tl = 0.0
1710  END IF
1711  DO it=1,nsplt
1712  IF (gridstruct%nested) neststruct%tracer_nest_timestep = &
1713 & neststruct%tracer_nest_timestep + 1
1714  CALL timing_on('COMM_TOTAL')
1715  CALL timing_on('COMM_TRACER')
1716  CALL complete_group_halo_update(q_pack, domain)
1717  CALL timing_off('COMM_TRACER')
1718  CALL timing_off('COMM_TOTAL')
1719  IF (gridstruct%nested) THEN
1720  DO iq=1,nq
1721  CALL nested_grid_bc_apply_intt_tlm(q(isd:ied, jsd:jed, :, iq)&
1722 & , q_tl(isd:ied, jsd:jed, :, iq), &
1723 & 0, 0, npx, npy, npz, bd, REAL(& & neststruct%tracer_nest_timestep) &
1724 & + REAL(nsplt*k_split), REAL(nsplt&
1725 & *k_split), neststruct%q_bc(iq), &
1726 & bctype=neststruct%nestbctype)
1727  end do
1728  END IF
1729 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq, &
1730 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
1731 !$OMP private(dp2, ra_x, ra_y, fx, fy)
1732  DO k=1,npz
1733  DO j=js,je
1734  DO i=is,ie
1735  dp2_tl(i, j) = dp1_tl(i, j, k) + rarea(i, j)*(mfx_tl(i, j, k&
1736 & )-mfx_tl(i+1, j, k)+mfy_tl(i, j, k)-mfy_tl(i, j+1, k))
1737  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(mfy&
1738 & (i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
1739  END DO
1740  END DO
1741  DO j=jsd,jed
1742  DO i=is,ie
1743  ra_x_tl(i, j) = xfx_tl(i, j, k) - xfx_tl(i+1, j, k)
1744  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
1745  END DO
1746  END DO
1747  DO j=js,je
1748  DO i=isd,ied
1749  ra_y_tl(i, j) = yfx_tl(i, j, k) - yfx_tl(i, j+1, k)
1750  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
1751  END DO
1752  END DO
1753  DO iq=1,nq
1754  IF (it .EQ. 1 .AND. trdm .GT. 1.e-4) THEN
1755  IF (hord .EQ. hord_pert) THEN
1756  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:&
1757 & ied, jsd:jed, k, iq), cx(is:ie+1, jsd:jed, &
1758 & k), cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied&
1759 & , js:je+1, k), cy_tl(isd:ied, js:je+1, k), &
1760 & npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx(&
1761 & is:ie+1, jsd:jed, k), xfx_tl(is:ie+1, jsd:&
1762 & jed, k), yfx(isd:ied, js:je+1, k), yfx_tl(&
1763 & isd:ied, js:je+1, k), gridstruct, bd, ra_x&
1764 & , ra_x_tl, ra_y, ra_y_tl, mfx=mfx(is:ie+1, &
1765 & js:je, k), mfx_tl=mfx_tl(is:ie+1, js:je, k)&
1766 & , mfy=mfy(is:ie, js:je+1, k), mfy_tl=mfy_tl&
1767 & (is:ie, js:je+1, k), mass=dp1(isd:ied, jsd:&
1768 & jed, k), mass_tl=dp1_tl(isd:ied, jsd:jed, k&
1769 & ), nord=nord_tr, damp_c=trdm)
1770  ELSE
1771  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:ied&
1772 & , jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k), &
1773 & cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, js:je+&
1774 & 1, k), cy_tl(isd:ied, js:je+1, k), npx, npy, &
1775 & hord_pert, fx, fx_tl, fy, fy_tl, xfx(is:ie+1, &
1776 & jsd:jed, k), xfx_tl(is:ie+1, jsd:jed, k), yfx(&
1777 & isd:ied, js:je+1, k), yfx_tl(isd:ied, js:je+1&
1778 & , k), gridstruct, bd, ra_x, ra_x_tl, ra_y, &
1779 & ra_y_tl, mfx=mfx(is:ie+1, js:je, k), mfx_tl=&
1780 & mfx_tl(is:ie+1, js:je, k), mfy=mfy(is:ie, js:&
1781 & je+1, k), mfy_tl=mfy_tl(is:ie, js:je+1, k), &
1782 & mass=dp1(isd:ied, jsd:jed, k), mass_tl=dp1_tl(&
1783 & isd:ied, jsd:jed, k), nord=nord_tr_pert, &
1784 & damp_c=trdm_pert)
1785  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1786  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1787  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k), &
1788  mass=dp1(isd:ied,jsd:jed,k), nord=nord_tr, damp_c=trdm)
1789  END IF
1790  ELSE IF (hord .EQ. hord_pert) THEN
1791  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:&
1792 & ied, jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k)&
1793 & , cx_tl(is:ie+1, jsd:jed, k), cy(isd:ied, js:&
1794 & je+1, k), cy_tl(isd:ied, js:je+1, k), npx, &
1795 & npy, hord, fx, fx_tl, fy, fy_tl, xfx(is:ie+1&
1796 & , jsd:jed, k), xfx_tl(is:ie+1, jsd:jed, k), &
1797 & yfx(isd:ied, js:je+1, k), yfx_tl(isd:ied, js:&
1798 & je+1, k), gridstruct, bd, ra_x, ra_x_tl, ra_y&
1799 & , ra_y_tl, mfx=mfx(is:ie+1, js:je, k), mfx_tl&
1800 & =mfx_tl(is:ie+1, js:je, k), mfy=mfy(is:ie, js&
1801 & :je+1, k), mfy_tl=mfy_tl(is:ie, js:je+1, k))
1802  ELSE
1803  CALL fv_tp_2d_tlm(q(isd:ied, jsd:jed, k, iq), q_tl(isd:ied, &
1804 & jsd:jed, k, iq), cx(is:ie+1, jsd:jed, k), cx_tl(&
1805 & is:ie+1, jsd:jed, k), cy(isd:ied, js:je+1, k), &
1806 & cy_tl(isd:ied, js:je+1, k), npx, npy, hord_pert&
1807 & , fx, fx_tl, fy, fy_tl, xfx(is:ie+1, jsd:jed, k)&
1808 & , xfx_tl(is:ie+1, jsd:jed, k), yfx(isd:ied, js:&
1809 & je+1, k), yfx_tl(isd:ied, js:je+1, k), &
1810 & gridstruct, bd, ra_x, ra_x_tl, ra_y, ra_y_tl, &
1811 & mfx=mfx(is:ie+1, js:je, k), mfx_tl=mfx_tl(is:ie+&
1812 & 1, js:je, k), mfy=mfy(is:ie, js:je+1, k), mfy_tl&
1813 & =mfy_tl(is:ie, js:je+1, k))
1814  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
1815  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
1816  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
1817  END IF
1818  DO j=js,je
1819  DO i=is,ie
1820  q_tl(i, j, k, iq) = ((q_tl(i, j, k, iq)*dp1(i, j, k)+q(i, &
1821 & j, k, iq)*dp1_tl(i, j, k)+rarea(i, j)*(fx_tl(i, j)-fx_tl&
1822 & (i+1, j)+fy_tl(i, j)-fy_tl(i, j+1)))*dp2(i, j)-(q(i, j, &
1823 & k, iq)*dp1(i, j, k)+(fx(i, j)-fx(i+1, j)+(fy(i, j)-fy(i&
1824 & , j+1)))*rarea(i, j))*dp2_tl(i, j))/dp2(i, j)**2
1825  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-fx&
1826 & (i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, j)
1827  END DO
1828  END DO
1829  END DO
1830  END DO
1831 ! npz
1832  IF (it .NE. nsplt) THEN
1833  CALL timing_on('COMM_TOTAL')
1834  CALL timing_on('COMM_TRACER')
1835  CALL start_group_halo_update_tlm(q_pack, q, q_tl, domain)
1836  CALL timing_off('COMM_TRACER')
1837  CALL timing_off('COMM_TOTAL')
1838  END IF
1839 !Apply nested-grid BCs
1840  IF (gridstruct%nested) THEN
1841  DO iq=1,nq
1842  CALL nested_grid_bc_apply_intt_tlm(q(isd:ied, jsd:jed, :, iq)&
1843 & , q_tl(isd:ied, jsd:jed, :, iq), &
1844 & 0, 0, npx, npy, npz, bd, REAL(&
1845 & neststruct%tracer_nest_timestep)&
1846 & , REAL(nsplt*k_split), neststruct&
1847 & %q_bc(iq), bctype=neststruct%&
1848 & nestbctype)
1849  end do
1850  END IF
1851  END DO
1852 ! nsplt
1853  IF (id_divg .GT. 0) THEN
1854  rdt = 1./(frac*dt)
1855 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dp1,xfx,yfx,rarea,rdt)
1856  DO k=1,npz
1857  DO j=js,je
1858  DO i=is,ie
1859  dp1_tl(i, j, k) = rarea(i, j)*rdt*(xfx_tl(i+1, j, k)-xfx_tl(&
1860 & i, j, k)+yfx_tl(i, j+1, k)-yfx_tl(i, j, k))
1861  dp1(i, j, k) = (xfx(i+1, j, k)-xfx(i, j, k)+(yfx(i, j+1, k)-&
1862 & yfx(i, j, k)))*rarea(i, j)*rdt
1863  END DO
1864  END DO
1865  END DO
1866  END IF
1867  END SUBROUTINE tracer_2d_nested_tlm
1868  SUBROUTINE tracer_2d_nested(q, dp1, mfx, mfy, cx, cy, gridstruct, bd, &
1869 & domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, &
1870 & nord_tr, trdm, k_split, neststruct, parent_grid, hord_pert, &
1871 & nord_tr_pert, trdm_pert, split_damp_tr)
1872  IMPLICIT NONE
1873  TYPE(fv_grid_bounds_type), INTENT(IN) :: bd
1874  INTEGER, INTENT(IN) :: npx
1875  INTEGER, INTENT(IN) :: npy
1876  INTEGER, INTENT(IN) :: npz
1877 ! number of tracers to be advected
1878  INTEGER, INTENT(IN) :: nq
1879  INTEGER, INTENT(IN) :: hord, nord_tr
1880  INTEGER, INTENT(IN) :: hord_pert, nord_tr_pert
1881  LOGICAL, INTENT(IN) :: split_damp_tr
1882  INTEGER, INTENT(IN) :: q_split, k_split
1883  INTEGER, INTENT(IN) :: id_divg
1884  REAL, INTENT(IN) :: dt, trdm, trdm_pert
1885  TYPE(group_halo_update_type), INTENT(INOUT) :: q_pack
1886 ! Tracers
1887  REAL, INTENT(INOUT) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq)
1888 ! DELP before dyn_core
1889  REAL, INTENT(INOUT) :: dp1(bd%isd:bd%ied, bd%jsd:bd%jed, npz)
1890 ! Mass Flux X-Dir
1891  REAL, INTENT(INOUT) :: mfx(bd%is:bd%ie+1, bd%js:bd%je, npz)
1892 ! Mass Flux Y-Dir
1893  REAL, INTENT(INOUT) :: mfy(bd%is:bd%ie, bd%js:bd%je+1, npz)
1894 ! Courant Number X-Dir
1895  REAL, INTENT(INOUT) :: cx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1896 ! Courant Number Y-Dir
1897  REAL, INTENT(INOUT) :: cy(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1898  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
1899  TYPE(fv_nest_type), INTENT(INOUT) :: neststruct
1900  TYPE(fv_atmos_type), INTENT(INOUT) :: parent_grid
1901  TYPE(domain2d), INTENT(INOUT) :: domain
1902 ! Local Arrays
1903  REAL :: dp2(bd%is:bd%ie, bd%js:bd%je)
1904  REAL :: fx(bd%is:bd%ie+1, bd%js:bd%je)
1905  REAL :: fy(bd%is:bd%ie, bd%js:bd%je+1)
1906  REAL :: ra_x(bd%is:bd%ie, bd%jsd:bd%jed)
1907  REAL :: ra_y(bd%isd:bd%ied, bd%js:bd%je)
1908  REAL :: xfx(bd%is:bd%ie+1, bd%jsd:bd%jed, npz)
1909  REAL :: yfx(bd%isd:bd%ied, bd%js:bd%je+1, npz)
1910  REAL :: cmax(npz)
1911  REAL :: cmax_t
1912  REAL :: c_global
1913  REAL :: frac, rdt
1914  INTEGER :: nsplt, nsplt_parent
1915  INTEGER, SAVE :: msg_split_steps=1
1916  INTEGER :: i, j, k, it, iq
1917  REAL, DIMENSION(:, :), POINTER :: area, rarea
1918  REAL, DIMENSION(:, :, :), POINTER :: sin_sg
1919  REAL, DIMENSION(:, :), POINTER :: dxa, dya, dx, dy
1920  INTEGER :: is, ie, js, je
1921  INTEGER :: isd, ied, jsd, jed
1922  INTRINSIC abs
1923  INTRINSIC max
1924  INTRINSIC int
1925  INTRINSIC real
1926  REAL :: max1
1927  REAL :: x2
1928  REAL :: x1
1929  REAL :: y2
1930  REAL :: y1
1931  is = bd%is
1932  ie = bd%ie
1933  js = bd%js
1934  je = bd%je
1935  isd = bd%isd
1936  ied = bd%ied
1937  jsd = bd%jsd
1938  jed = bd%jed
1939  area => gridstruct%area
1940  rarea => gridstruct%rarea
1941  sin_sg => gridstruct%sin_sg
1942  dxa => gridstruct%dxa
1943  dya => gridstruct%dya
1944  dx => gridstruct%dx
1945  dy => gridstruct%dy
1946 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,xfx,dxa,dy, &
1947 !$OMP sin_sg,cy,yfx,dya,dx)
1948  DO k=1,npz
1949  DO j=jsd,jed
1950  DO i=is,ie+1
1951  IF (cx(i, j, k) .GT. 0.) THEN
1952  xfx(i, j, k) = cx(i, j, k)*dxa(i-1, j)*dy(i, j)*sin_sg(i-1, &
1953 & j, 3)
1954  ELSE
1955  xfx(i, j, k) = cx(i, j, k)*dxa(i, j)*dy(i, j)*sin_sg(i, j, 1&
1956 & )
1957  END IF
1958  END DO
1959  END DO
1960  DO j=js,je+1
1961  DO i=isd,ied
1962  IF (cy(i, j, k) .GT. 0.) THEN
1963  yfx(i, j, k) = cy(i, j, k)*dya(i, j-1)*dx(i, j)*sin_sg(i, j-&
1964 & 1, 4)
1965  ELSE
1966  yfx(i, j, k) = cy(i, j, k)*dya(i, j)*dx(i, j)*sin_sg(i, j, 2&
1967 & )
1968  END IF
1969  END DO
1970  END DO
1971  END DO
1972 !--------------------------------------------------------------------------------
1973  IF (q_split .EQ. 0) THEN
1974 ! Determine nsplt
1975 !$OMP parallel do default(none) shared(is,ie,js,je,npz,cmax,cx,cy,sin_sg) &
1976 !$OMP private(cmax_t )
1977  DO k=1,npz
1978  cmax(k) = 0.
1979  IF (k .LT. 4) THEN
1980 ! Top layers: C < max( abs(c_x), abs(c_y) )
1981  DO j=js,je
1982  DO i=is,ie
1983  IF (cx(i, j, k) .GE. 0.) THEN
1984  x1 = cx(i, j, k)
1985  ELSE
1986  x1 = -cx(i, j, k)
1987  END IF
1988  IF (cy(i, j, k) .GE. 0.) THEN
1989  y1 = cy(i, j, k)
1990  ELSE
1991  y1 = -cy(i, j, k)
1992  END IF
1993  IF (x1 .LT. y1) THEN
1994  cmax_t = y1
1995  ELSE
1996  cmax_t = x1
1997  END IF
1998  IF (cmax_t .LT. cmax(k)) THEN
1999  cmax(k) = cmax(k)
2000  ELSE
2001  cmax(k) = cmax_t
2002  END IF
2003  END DO
2004  END DO
2005  ELSE
2006  DO j=js,je
2007  DO i=is,ie
2008  IF (cx(i, j, k) .GE. 0.) THEN
2009  x2 = cx(i, j, k)
2010  ELSE
2011  x2 = -cx(i, j, k)
2012  END IF
2013  IF (cy(i, j, k) .GE. 0.) THEN
2014  y2 = cy(i, j, k)
2015  ELSE
2016  y2 = -cy(i, j, k)
2017  END IF
2018  IF (x2 .LT. y2) THEN
2019  max1 = y2
2020  ELSE
2021  max1 = x2
2022  END IF
2023  cmax_t = max1 + 1. - sin_sg(i, j, 5)
2024  IF (cmax_t .LT. cmax(k)) THEN
2025  cmax(k) = cmax(k)
2026  ELSE
2027  cmax(k) = cmax_t
2028  END IF
2029  END DO
2030  END DO
2031  END IF
2032  END DO
2033  CALL mp_reduce_max(cmax, npz)
2034 ! find global max courant number and define nsplt to scale cx,cy,mfx,mfy
2035  c_global = cmax(1)
2036  IF (npz .NE. 1) THEN
2037 ! if NOT shallow water test case
2038  DO k=2,npz
2039  IF (cmax(k) .LT. c_global) THEN
2040  c_global = c_global
2041  ELSE
2042  c_global = cmax(k)
2043  END IF
2044  END DO
2045  END IF
2046  nsplt = int(1. + c_global)
2047  IF (is_master() .AND. nsplt .GT. 3) WRITE(*, *) 'Tracer_2d_split='&
2048 & , nsplt, c_global
2049  ELSE
2050  nsplt = q_split
2051  IF (gridstruct%nested .AND. neststruct%nestbctype .GT. 1) THEN
2052  IF (q_split/parent_grid%flagstruct%q_split .LT. 1) THEN
2053  msg_split_steps = 1
2054  ELSE
2055  msg_split_steps = q_split/parent_grid%flagstruct%q_split
2056  END IF
2057  END IF
2058  END IF
2059 !--------------------------------------------------------------------------------
2060  frac = 1./REAL(nsplt)
2061  IF (nsplt .NE. 1) THEN
2062 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,cx,frac,xfx,mfx,cy,yfx,mfy)
2063  DO k=1,npz
2064  DO j=jsd,jed
2065  DO i=is,ie+1
2066  cx(i, j, k) = cx(i, j, k)*frac
2067  xfx(i, j, k) = xfx(i, j, k)*frac
2068  END DO
2069  END DO
2070  DO j=js,je
2071  DO i=is,ie+1
2072  mfx(i, j, k) = mfx(i, j, k)*frac
2073  END DO
2074  END DO
2075  DO j=js,je+1
2076  DO i=isd,ied
2077  cy(i, j, k) = cy(i, j, k)*frac
2078  yfx(i, j, k) = yfx(i, j, k)*frac
2079  END DO
2080  END DO
2081  DO j=js,je+1
2082  DO i=is,ie
2083  mfy(i, j, k) = mfy(i, j, k)*frac
2084  END DO
2085  END DO
2086  END DO
2087  END IF
2088  DO it=1,nsplt
2089  IF (gridstruct%nested) neststruct%tracer_nest_timestep = &
2090 & neststruct%tracer_nest_timestep + 1
2091  CALL timing_on('COMM_TOTAL')
2092  CALL timing_on('COMM_TRACER')
2093  CALL complete_group_halo_update(q_pack, domain)
2094  CALL timing_off('COMM_TRACER')
2095  CALL timing_off('COMM_TOTAL')
2096  IF (gridstruct%nested) THEN
2097  DO iq=1,nq
2098  CALL nested_grid_bc_apply_intt(q(isd:ied, jsd:jed, :, iq), 0, &
2099 & 0, npx, npy, npz, bd, REAL(neststruct& & %tracer_nest_timestep) + REAL(nsplt*& & k_split), REAL(nsplt*k_split), &
2100 & neststruct%q_bc(iq), neststruct%&
2101 & nestbctype)
2102  end do
2103  END IF
2104 !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,mfx,mfy,rarea,nq, &
2105 !$OMP area,xfx,yfx,q,cx,cy,npx,npy,hord,gridstruct,bd,it,nsplt,nord_tr,trdm) &
2106 !$OMP private(dp2, ra_x, ra_y, fx, fy)
2107  DO k=1,npz
2108  DO j=js,je
2109  DO i=is,ie
2110  dp2(i, j) = dp1(i, j, k) + (mfx(i, j, k)-mfx(i+1, j, k)+(mfy&
2111 & (i, j, k)-mfy(i, j+1, k)))*rarea(i, j)
2112  END DO
2113  END DO
2114  DO j=jsd,jed
2115  DO i=is,ie
2116  ra_x(i, j) = area(i, j) + (xfx(i, j, k)-xfx(i+1, j, k))
2117  END DO
2118  END DO
2119  DO j=js,je
2120  DO i=isd,ied
2121  ra_y(i, j) = area(i, j) + (yfx(i, j, k)-yfx(i, j+1, k))
2122  END DO
2123  END DO
2124  DO iq=1,nq
2125  IF (it .EQ. 1 .AND. trdm .GT. 1.e-4) THEN
2126  IF (hord .EQ. hord_pert) THEN
2127  CALL fv_tp_2d(q(isd:ied, jsd:jed, k, iq), cx(is:ie+1, &
2128 & jsd:jed, k), cy(isd:ied, js:je+1, k), npx, npy&
2129 & , hord, fx, fy, xfx(is:ie+1, jsd:jed, k), yfx(&
2130 & isd:ied, js:je+1, k), gridstruct, bd, ra_x, &
2131 & ra_y, mfx(is:ie+1, js:je, k), mfy(is:ie, js:je+&
2132 & 1, k), dp1(isd:ied, jsd:jed, k), nord_tr, trdm)
2133  ELSE
2134  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
2135  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
2136  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k), &
2137  mass=dp1(isd:ied,jsd:jed,k), nord=nord_tr, damp_c=trdm)
2138  END IF
2139  ELSE IF (hord .EQ. hord_pert) THEN
2140  CALL fv_tp_2d(q(isd:ied, jsd:jed, k, iq), cx(is:ie+1, jsd&
2141 & :jed, k), cy(isd:ied, js:je+1, k), npx, npy, hord&
2142 & , fx, fy, xfx(is:ie+1, jsd:jed, k), yfx(isd:ied, &
2143 & js:je+1, k), gridstruct, bd, ra_x, ra_y, mfx=mfx(&
2144 & is:ie+1, js:je, k), mfy=mfy(is:ie, js:je+1, k))
2145  ELSE
2146  call fv_tp_2d(q(isd:ied,jsd:jed,k,iq), cx(is:ie+1,jsd:jed,k), cy(isd:ied,js:je+1,k), &
2147  npx, npy, hord, fx, fy, xfx(is:ie+1,jsd:jed,k), yfx(isd:ied,js:je+1,k), &
2148  gridstruct, bd, ra_x, ra_y, mfx=mfx(is:ie+1,js:je,k), mfy=mfy(is:ie,js:je+1,k))
2149  END IF
2150  DO j=js,je
2151  DO i=is,ie
2152  q(i, j, k, iq) = (q(i, j, k, iq)*dp1(i, j, k)+(fx(i, j)-fx&
2153 & (i+1, j)+(fy(i, j)-fy(i, j+1)))*rarea(i, j))/dp2(i, j)
2154  END DO
2155  END DO
2156  END DO
2157  END DO
2158 ! npz
2159  IF (it .NE. nsplt) THEN
2160  CALL timing_on('COMM_TOTAL')
2161  CALL timing_on('COMM_TRACER')
2162  CALL start_group_halo_update(q_pack, q, domain)
2163  CALL timing_off('COMM_TRACER')
2164  CALL timing_off('COMM_TOTAL')
2165  END IF
2166 !Apply nested-grid BCs
2167  IF (gridstruct%nested) THEN
2168  DO iq=1,nq
2169  CALL nested_grid_bc_apply_intt(q(isd:ied, jsd:jed, :, iq), 0, &
2170 & 0, npx, npy, npz, bd, REAL(neststruct& & %tracer_nest_timestep), REAL(nsplt*& & k_split), neststruct%q_bc(iq), &
2171 & neststruct%nestbctype)
2172  end do
2173  END IF
2174  END DO
2175 ! nsplt
2176  IF (id_divg .GT. 0) THEN
2177  rdt = 1./(frac*dt)
2178 !$OMP parallel do default(none) shared(is,ie,js,je,npz,dp1,xfx,yfx,rarea,rdt)
2179  DO k=1,npz
2180  DO j=js,je
2181  DO i=is,ie
2182  dp1(i, j, k) = (xfx(i+1, j, k)-xfx(i, j, k)+(yfx(i, j+1, k)-&
2183 & yfx(i, j, k)))*rarea(i, j)*rdt
2184  END DO
2185  END DO
2186  END DO
2187  END IF
2188  END SUBROUTINE tracer_2d_nested
2189 
2190 end module fv_tracer2d_tlm_mod
2191 
subroutine, public nested_grid_bc_apply_intt(var_nest, istag, jstag, npx, npy, npz, bd, step, split, bc, bctype)
subroutine, public fv_tp_2d_tlm(q, q_tl, crx, crx_tl, cry, cry_tl, npx, npy, hord, fx, fx_tl, fy, fy_tl, xfx, xfx_tl, yfx, yfx_tl, gridstruct, bd, ra_x, ra_x_tl, ra_y, ra_y_tl, mfx, mfx_tl, mfy, mfy_tl, mass, mass_tl, nord, damp_c)
real, dimension(:,:,:), allocatable nest_fx_south_accum
subroutine, public tracer_2d_1l_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
subroutine, public tracer_2d_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
Definition: mpp.F90:39
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, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr)
subroutine, public complete_group_halo_update(group, group_tl, domain)
Definition: fv_mp_tlm.F90:815
real, dimension(:,:,:), allocatable nest_fx_north_accum
real, dimension(:,:,:), allocatable nest_fx_west_accum
integer, parameter, public ng
subroutine timing_on(blk_name)
real, dimension(:,:,:), allocatable nest_fx_east_accum
subroutine, public copy_corners(q, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
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, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, dpa)
#define max(a, b)
Definition: mosaic_util.h:33
subroutine, public nested_grid_bc_apply_intt_tlm(var_nest, var_nest_tl, istag, jstag, npx, npy, npz, bd, step, split, bc, bctype)
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, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr, 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_tlm.F90:85
subroutine, public copy_corners_tlm(q, q_tl, npx, npy, dir, nested, bd, sw_corner, se_corner, nw_corner, ne_corner)
subroutine, public tracer_2d_nested_tlm(q, q_tl, dp1, dp1_tl, mfx, mfx_tl, mfy, mfy_tl, cx, cx_tl, cy, cy_tl, gridstruct, bd, domain, npx, npy, npz, nq, hord, q_split, dt, id_divg, q_pack, nord_tr, trdm, k_split, neststruct, parent_grid, hord_pert, nord_tr_pert, trdm_pert, split_damp_tr)
subroutine timing_off(blk_name)