Line data Source code
1 : !-----------------------------------------------------------------------
2 : !BOP
3 : ! !ROUTINE: trac2d --- Remap Lagrangian to fixed coordinates
4 : !
5 : ! !INTERFACE:
6 32256 : subroutine trac2d( grid, dp1, tracer, cx, cy, &
7 32256 : mfx, mfy, iord, jord, fill, &
8 32256 : nlo, nhi, va, flx )
9 :
10 : ! !USES:
11 :
12 : use shr_kind_mod, only : r8 => shr_kind_r8, r4 => shr_kind_r4
13 : use tp_core, only : tp2c
14 : use fill_module, only : fillxy
15 : use dynamics_vars, only : T_FVDYCORE_GRID
16 : use FVperf_module, only : FVstartclock, FVstopclock, FVbarrierclock
17 :
18 : #if defined( SPMD )
19 : use parutilitiesmodule, only: maxop, parcollective
20 : use mod_comm, only : mp_send4d_ns, mp_recv4d_ns, &
21 : mp_send4d_ns_r4, mp_recv4d_ns_r4, &
22 : mp_send3d_2, mp_recv3d_2
23 : #endif
24 :
25 : implicit none
26 :
27 : ! !INPUT PARAMETERS:
28 :
29 : type (T_FVDYCORE_GRID), intent(inout) :: grid
30 : integer, intent(in):: iord, jord
31 :
32 : logical, intent(in):: fill
33 : integer, intent(in):: nlo, nhi ! Tracer index range
34 :
35 : ! !INPUT/OUTPUT PARAMETERS:
36 : real(r8), intent(inout):: dp1(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
37 : real(r8), intent(inout):: cx(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d,grid%kfirst:grid%klast)
38 : real(r8), intent(inout):: cy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast)
39 : real(r8), intent(inout):: mfx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
40 : real(r8), intent(inout):: mfy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast)
41 : real(r8), intent(inout):: tracer(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast,grid%ntotq)
42 :
43 : ! !OUTPUT PARAMETERS:
44 : real(r8), intent(out):: va(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
45 : real(r8), intent(out):: flx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
46 :
47 : ! !DESCRIPTION:
48 : !
49 : ! Perform large-time-step tracer transport using accumulated Courant
50 : ! numbers (cx, cy) and the mass fluxes (mfx, mfy) within the Lagrangian
51 : ! layers. This routine is 100\% parallel in the vertical direction
52 : ! (with SMP). Merdional Courant number will be further split, if
53 : ! necessary, to ensure stability. Cy <= 1 away from poles; Cy $\le$
54 : ! 1/2 at the latitudes closest to the poles.
55 : !
56 : ! !REVISION HISTORY:
57 : !
58 : ! SJL 99.04.13: Delivery
59 : ! WS 99.05.26: Added jfirst:jlast concept; im, jm, km as parameters
60 : ! replaced IMR, JMR, JNP, NL with IM, JM-1, JM and KM
61 : ! WS 99.09.27: Documentation; indentation; jfirst:jlast
62 : ! WS 99.09.30: Ghosting; loop limits; full parallelization; tested
63 : ! SJL 99.10.15: nsplt migrated to outermost loop to remove bug
64 : ! SJL 99.12.19: Local 2D arrays trimmed!
65 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
66 : ! WS 00.07.13: Changed PILGRIM API
67 : ! AAM 00.08.29: Added kfirst, klast
68 : ! AAM 01.06.27: Added y communicators
69 : ! SJL 30.07.01: MPI optimization/simplification
70 : ! WS 02.04.24: New mod_comm interfaces
71 : ! WS 02.07.04: Fixed 2D decomposition bug dest/src for mp_send3d
72 : ! WS 03.11.19: Merged in CAM changes by Mirin
73 : ! WS 03.12.03: Added GRID as argument, dynamics_vars removed
74 : ! WS 04.08.25: Simplification of interface with GRID
75 : ! WS 04.10.07: Removed dependency on spmd_dyn; info now in GRID
76 : ! WS 05.04.04: Transitioned to type T_TRACERS (supports r4 and r8)
77 : ! WS 05.04.09: Each tracer now ghosted individually (save buffers)
78 : ! WS 05.04.12: Full support for either r4 or r8 tracers
79 : ! WS 05.05.25: Merged CAM and GEOS5, e.g. nsplt(k) opt. in CAM
80 : ! PW 05.10.12: Changes for Cray X1(E), alternative implementation
81 : ! of double buffering logic
82 : ! WS 06.09.08: Magic numbers are now F90 parameters
83 : !
84 : !EOP
85 : !---------------------------------------------------------------------
86 : !BOC
87 :
88 : real(r8), parameter :: D1EM10 = 1.0e-10_r8
89 : real(r8), parameter :: D1_0 = 1.0_r8
90 : real(r8), parameter :: D0_0 = 0.0_r8
91 :
92 : ! Local variables:
93 : ! 2d arrays
94 64512 : real(r8) a2(grid%im,grid%jfirst:grid%jlast)
95 64512 : real(r8) fx(grid%im,grid%jfirst:grid%jlast)
96 64512 : real(r8) fy(grid%im,grid%jfirst:grid%jlast+1)
97 64512 : real(r8) cymax(grid%kfirst:grid%klast)
98 : ! Temporary r8 array for Q
99 : real(r8) :: &
100 64512 : q_r8(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d,grid%kfirst:grid%klast,1:2)
101 :
102 64512 : real(r8) dp2(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
103 64512 : logical ffsl(grid%jm,grid%kfirst:grid%klast)
104 64512 : integer :: nsplt(grid%kfirst:grid%klast)
105 :
106 : integer :: im, jm, km ! Dimensions
107 : integer :: ng ! Max number of ghost latitudes
108 : integer :: jfirst, jlast, kfirst, klast ! YZ decomposition limits
109 : integer :: cur, nxt ! current and next q_r8 buffer indices
110 :
111 : integer i, j, k
112 : integer it, iq, kq, max_nsplt
113 : integer :: k_courant, kend
114 : integer ktot
115 : integer js1gd, js2g0, js2gd, jn2g0,jn2gd,jn1g1,jn1gd
116 : #if defined( SPMD )
117 : integer :: dest, src
118 : #endif
119 :
120 : real(r8) cy_global
121 : real(r8) frac
122 : real(r8) cmax
123 : real(r8) sum1, sum2
124 :
125 32256 : cur = 1
126 32256 : nxt = 2
127 :
128 32256 : im = grid%im
129 32256 : jm = grid%jm
130 32256 : km = grid%km
131 32256 : ng = grid%ng_d
132 :
133 32256 : jfirst = grid%jfirst
134 32256 : jlast = grid%jlast
135 32256 : kfirst = grid%kfirst
136 32256 : klast = grid%klast
137 :
138 32256 : ktot = klast - kfirst + 1
139 32256 : js2g0 = max(2,jfirst)
140 32256 : jn2g0 = min(jm-1,jlast)
141 32256 : jn1g1 = min(jm,jlast+1)
142 32256 : js1gd = max(1,jfirst-ng) ! NG latitudes on S (starting at 1)
143 32256 : js2gd = max(2,jfirst-ng) ! NG latitudes on S (starting at 2)
144 32256 : jn2gd = min(jm-1,jlast+ng) ! NG latitudes on S (ending at jm-1)
145 32256 : jn1gd = min(jm,jlast+ng) ! NG latitudes on N (ending at jm)
146 :
147 : #if defined( SPMD )
148 32256 : call FVstartclock(grid,'---TRAC2D_COMM')
149 : call mp_send4d_ns( grid%commyz, im, jm, km, &
150 : 1, jfirst, jlast, kfirst, klast, &
151 32256 : ng, ng, cx )
152 : ! Send one latitude of both cy and mfy to the south
153 32256 : dest = grid%iam-1
154 32256 : src = grid%iam+1
155 32256 : if ( mod(grid%iam,grid%npr_y) == 0 ) dest = -1
156 32256 : if ( mod(grid%iam+1,grid%npr_y) == 0 ) src = -1
157 : call mp_send3d_2( grid%commyz, dest, src, im, jm, km, &
158 : 1, im, jfirst, jlast+1, kfirst, klast, &
159 32256 : 1, im, jfirst, jfirst, kfirst, klast, cy, mfy)
160 32256 : call FVstopclock(grid,'---TRAC2D_COMM')
161 : #endif
162 :
163 : !$omp parallel do default(shared) private(i,j,k,cmax)
164 118272 : do k=kfirst,klast
165 86016 : cymax(k) = D0_0
166 374976 : do j=js2g0,jlast
167 256704 : cmax = D0_0
168 74187456 : do i=1,im
169 74187456 : cmax = max( abs(cy(i,j,k)), cmax)
170 : enddo
171 342720 : cymax(k) = max(cymax(k), cmax*(D1_0 + grid%sine(j)**16) )
172 : enddo
173 : enddo
174 :
175 : #if defined( SPMD )
176 32256 : call FVstartclock(grid,'---TRAC2D_COMM')
177 : call mp_recv4d_ns( grid%commyz, im, jm, km, &
178 : 1, jfirst, jlast, kfirst, klast, &
179 32256 : ng, ng, cx )
180 : call mp_recv3d_2( grid%commyz, src, im, jm, km, &
181 : 1, im, jfirst, jlast+1, kfirst, klast, &
182 32256 : 1, im, jlast+1, jlast+1, kfirst, klast, cy, mfy)
183 :
184 32256 : call parcollective( grid%comm_y, MAXOP, ktot, cymax )
185 32256 : call FVstopclock(grid,'---TRAC2D_COMM')
186 : #endif
187 :
188 : !---------------------------------------------------------------------
189 : ! Determine the required value of nsplt for each level
190 : !---------------------------------------------------------------------
191 118272 : nsplt(:) = int( D1_0 + cymax(:) )
192 118272 : max_nsplt = maxval( nsplt(:) )
193 : #if defined( SPMD )
194 32256 : call FVstartclock(grid,'---TRAC2D_COMM')
195 32256 : call parcollective( grid%comm_z, MAXOP, max_nsplt ) ! Find global max
196 32256 : call FVstopclock(grid,'---TRAC2D_COMM')
197 : #endif
198 32256 : if (grid%ptop>1._r8) then
199 118272 : nsplt(:) = max_nsplt
200 : endif
201 32256 : do k_courant = klast,kfirst,-1
202 32256 : if( nsplt(k_courant) > 1 ) then
203 : exit
204 : end if
205 : end do
206 32256 : k_courant = max( kfirst,k_courant )
207 : !!! if (max_nsplt /= 1) write(iulog,*) 'trac2d: max_nsplt,k_courant = ', max_nsplt,k_courant
208 : !!! write(iulog,*) "max_nsplt", max_nsplt, "k_cour", k_courant, "nsplt", nsplt(:)
209 :
210 : !$omp parallel do default(shared) private(i,j,k,frac) schedule(dynamic,1)
211 :
212 : #if !defined(USE_OMP)
213 : !CSD$ PARALLEL DO PRIVATE (I, J, K, FRAC)
214 : #endif
215 118272 : do 4000 k=kfirst,klast
216 :
217 86016 : if( nsplt(k) .ne. 1 ) then
218 86016 : frac = D1_0 / nsplt(k)
219 846720 : do j=js2gd,jn2gd
220 219929472 : do i=1,im
221 219843456 : cx(i,j,k) = cx(i,j,k) * frac ! cx ghosted on N*ng S*ng
222 : enddo
223 : enddo
224 :
225 341376 : do j=js2g0,jn2g0
226 73885056 : do i=1,im
227 73799040 : mfx(i,j,k) = mfx(i,j,k) * frac
228 : enddo
229 : enddo
230 :
231 427392 : do j=js2g0,jn1g1
232 98743680 : do i=1,im
233 98316288 : cy(i,j,k) = cy(i,j,k) * frac ! cy ghosted on N
234 98657664 : mfy(i,j,k) = mfy(i,j,k) * frac ! mfy ghosted on N
235 : enddo
236 : enddo
237 : endif
238 :
239 341376 : do j=js2g0,jn2g0
240 73885056 : do i=1,im
241 73799040 : if(cy(i,j,k)*cy(i,j+1,k) > D0_0) then
242 70407564 : if( cy(i,j,k) > D0_0) then
243 34032985 : va(i,j,k) = cy(i,j,k)
244 : else
245 36374579 : va(i,j,k) = cy(i,j+1,k) ! cy ghosted on N
246 : endif
247 : else
248 3136116 : va(i,j,k) = D0_0
249 : endif
250 : enddo
251 : enddo
252 :
253 : ! Check if FFSL extension is needed.
254 :
255 846720 : do j=js2gd,jn2gd ! flux needed on N*ng S*ng
256 760704 : ffsl(j,k) = .false.
257 215595885 : do i=1,im
258 215509869 : if( abs(cx(i,j,k)) > D1_0 ) then ! cx ghosted on N*ng S*ng
259 22729 : ffsl(j,k) = .true.
260 22729 : exit
261 : endif
262 : enddo
263 : enddo
264 :
265 : ! Scale E-W mass fluxes by CX if FFSL
266 341376 : do j=js2g0,jn2g0
267 341376 : if( ffsl(j,k) ) then
268 2576435 : do i=1,im
269 5135040 : flx(i,j,k) = mfx(i,j,k) / sign( max(abs(cx(i,j,k)), D1EM10), &
270 5143955 : cx(i,j,k) )
271 : enddo
272 : else
273 71222605 : do i=1,im
274 71222605 : flx(i,j,k) = mfx(i,j,k)
275 : enddo
276 : endif
277 : enddo
278 32256 : 4000 continue
279 : #if !defined(USE_OMP)
280 : !CSD$ END PARALLEL DO
281 : #endif
282 :
283 32256 : call FVbarrierclock(grid,'sync_trac2d_tracer',grid%commyz)
284 32256 : call FVstartclock(grid,'---TRAC2D_TRACER')
285 :
286 96768 : do 6000 it=1, max_nsplt
287 64512 : if ( it == 1 ) then
288 32256 : kend = klast ! The entire vertical slab needs to be considered
289 : else
290 32256 : kend = k_courant ! Only the subset including courant # > 1 considered
291 : endif
292 : ! WS 05.04.06: send only the first tracer the rest at end of do iq loop
293 : ! NOTE: there is per definition at least one tracer
294 0 : q_r8(1:im,jfirst:jlast,kfirst:kend,1) = &
295 149388288 : tracer(1:im,jfirst:jlast,kfirst:kend,nlo)
296 : #if defined( SPMD )
297 64512 : call FVstartclock(grid,'---TRAC2D_TRACER_COMM')
298 : call mp_send4d_ns( grid%commyz, im, jm, km, &
299 : 1, jfirst, jlast, kfirst, kend, &
300 64512 : ng, ng, q_r8(1,jfirst-ng,kfirst,1) )
301 64512 : call FVstopclock(grid,'---TRAC2D_TRACER_COMM')
302 : #endif
303 :
304 : !$omp parallel do default(shared) private(i,j,k,sum1,sum2)
305 :
306 236544 : do 3000 k=kfirst,kend
307 172032 : if (it <= nsplt(k)) then
308 682752 : do j=js2g0,jn2g0
309 147087360 : do i=1,im-1
310 439729920 : dp2(i,j,k) = dp1(i,j,k) + mfx(i,j,k) - mfx(i+1,j,k) + &
311 586817280 : (mfy(i,j,k) - mfy(i,j+1,k)) * grid%acosp(j)
312 : enddo
313 1021440 : dp2(im,j,k) = dp1(im,j,k) + mfx(im,j,k) - mfx(1,j,k) + &
314 1704192 : (mfy(im,j,k) - mfy(im,j+1,k)) * grid%acosp(j)
315 : enddo
316 :
317 172032 : if ( jfirst == 1 ) then
318 2688 : sum1 = D0_0
319 776832 : do i=1,im
320 776832 : sum1 = sum1 + mfy(i,2,k)
321 : end do
322 :
323 2688 : sum1 = - sum1 * grid%rcap
324 776832 : do i=1,im
325 776832 : dp2(i,1,k) = dp1(i,1,k) + sum1
326 : enddo
327 : endif
328 :
329 172032 : if ( jlast == jm ) then
330 2688 : sum2 = D0_0
331 776832 : do i=1,im
332 776832 : sum2 = sum2 + mfy(i,jm,k)
333 : end do
334 :
335 2688 : sum2 = sum2 * grid%rcap
336 776832 : do i=1,im
337 776832 : dp2(i,jm,k) = dp1(i,jm,k) + sum2
338 : enddo
339 : endif
340 : endif
341 64512 : 3000 continue
342 :
343 15547392 : do iq = nlo, nhi
344 : #if defined( SPMD )
345 15482880 : call FVstartclock(grid,'---TRAC2D_TRACER_COMM')
346 : !
347 : ! The buffer indices are exchanged, so that cur points to the current buffer,
348 : ! while nxt points to the one which will be used next.
349 : !
350 15482880 : if ( mod(iq-nlo+1,2) == 0 ) then
351 : cur = 2
352 : nxt = 1
353 : else
354 7741440 : cur = 1
355 7741440 : nxt = 2
356 : endif
357 : call mp_recv4d_ns( grid%commyz, im, jm, km, &
358 : 1, jfirst, jlast, kfirst, kend, &
359 15482880 : ng, ng, q_r8(1,jfirst-ng,kfirst,cur) )
360 :
361 : !
362 : ! Pre-send the next tracer
363 : !
364 15482880 : if ( iq < nhi ) then
365 15418368 : q_r8(1:im,jfirst:jlast,kfirst:kend,nxt) = &
366 35719219200 : tracer(1:im,jfirst:jlast,kfirst:kend,iq+1)
367 : call mp_send4d_ns( grid%commyz, im, jm, km, &
368 : 1, jfirst, jlast, kfirst, kend, &
369 15418368 : ng, ng, q_r8(1,jfirst-ng,kfirst,nxt) )
370 : endif
371 15482880 : call FVstopclock(grid,'---TRAC2D_TRACER_COMM')
372 : #else
373 : !
374 : ! No message passing -- simply copy the tracer into q_r8
375 : !
376 : q_r8(1:im,jfirst:jlast,kfirst:kend,cur) = &
377 : tracer(1:im,jfirst:jlast,kfirst:kend,iq)
378 : #endif
379 :
380 : #if (!defined USE_OMP)
381 : !CSD$ PARALLEL DO PRIVATE (I, J, K, KQ, FX, FY, A2)
382 : #endif
383 56835072 : do 5000 k=kfirst,kend
384 41287680 : if ( it <= nsplt(k) ) then
385 41287680 : call tp2c(a2, va(1,jfirst,k), q_r8(1:,jfirst-ng:,k,cur), &
386 41287680 : cx(1,jfirst-ng,k) , cy(1,jfirst,k), &
387 : im, jm, iord, jord, ng, &
388 41287680 : fx, fy, ffsl(1,k), grid%rcap, grid%acosp, &
389 : flx(1,jfirst,k), mfy(1,jfirst,k), &
390 123863040 : grid%cosp, 1, jfirst, jlast )
391 :
392 165150720 : do j=jfirst,jlast
393 35837706240 : do i=1,im
394 35796418560 : q_r8(i,j,k,cur) = q_r8(i,j,k,cur)*dp1(i,j,k) + a2(i,j)
395 : enddo
396 : enddo
397 :
398 41287680 : if (fill) call fillxy (q_r8(1:,jfirst:,k,cur), im, jm, jfirst, &
399 41287680 : jlast, grid%acap, grid%cosp, grid%acosp)
400 :
401 165150720 : do j=jfirst,jlast
402 35837706240 : do i=1,im
403 35796418560 : tracer(i,j,k,iq) = q_r8(i,j,k,cur) / dp2(i,j,k)
404 : enddo
405 : enddo
406 : endif
407 15482880 : 5000 continue
408 : #if (!defined USE_OMP)
409 : !CSD$ END PARALLEL DO
410 : #endif
411 :
412 : enddo ! End of do iq=nlo, nhi
413 :
414 : !$omp parallel do private(i, j, k) schedule( dynamic,1 )
415 236544 : do k=kfirst,kend
416 236544 : if ( it < nsplt(k) ) then
417 344064 : do j=jfirst,jlast
418 74661888 : do i=1,im
419 74575872 : dp1(i,j,k) = dp2(i,j,k)
420 : enddo
421 : enddo
422 : endif
423 : enddo
424 :
425 32256 : 6000 continue
426 32256 : call FVstopclock(grid,'---TRAC2D_TRACER')
427 :
428 32256 : return
429 : !EOC
430 32256 : end subroutine trac2d
431 : !-----------------------------------------------------------------------
432 :
433 :
|
