Line data Source code
1 : #if defined( UNICOSMP ) || defined ( NEC_SX )
2 : #define VECTORIZE
3 : #endif
4 : module tp_core
5 : !BOP
6 : !
7 : ! !MODULE: tp_core --- Utilities for the transport core
8 : !
9 : ! !USES:
10 : use shr_kind_mod, only : r8 => shr_kind_r8
11 :
12 : !
13 : ! !PUBLIC MEMBER FUNCTIONS:
14 : public tp2c, tp2d, xtp, xtpv, fxppm, xmist, steepx, lmppm
15 : public huynh, ytp, ymist, fyppm, tpcc, ycc
16 : !
17 : ! !DESCRIPTION:
18 : !
19 : ! This module provides
20 : !
21 : ! \begin{tabular}{|l|l|} \hline \hline
22 : ! tp2c & \\ \hline
23 : ! tp2d & \\ \hline
24 : ! xtp & \\ \hline
25 : ! fxppm & \\ \hline
26 : ! xmist & \\ \hline
27 : ! steepx & \\ \hline
28 : ! lmppm & \\ \hline
29 : ! huynh & \\ \hline
30 : ! ytp & \\ \hline
31 : ! ymist & \\ \hline
32 : ! fyppm & \\ \hline
33 : ! tpcc & \\ \hline
34 : ! ycc & \\ \hline
35 : ! \hline
36 : ! \end{tabular}
37 : !
38 : ! !REVISION HISTORY:
39 : ! 01.01.15 Lin Routines coalesced into this module
40 : ! 01.03.26 Sawyer Additional ProTeX documentation
41 : ! 03.11.19 Sawyer Merged in CAM changes by Mirin
42 : ! 04.10.07 Sawyer ompinner now from dynamics_vars
43 : ! 05.03.25 Todling shr_kind_r8 can only be referenced once (MIPSpro-7.4.2)
44 : ! 05.05.25 Sawyer Merged CAM and GEOS5 versions (mostly CAM)
45 : ! 06.09.06 Sawyer Turned "magic numbers" into F90 parameters
46 : !
47 : !EOP
48 : !-----------------------------------------------------------------------
49 :
50 : ! Magic numbers used in this module
51 :
52 : private
53 : real(r8), parameter :: D0_0 = 0.0_r8
54 : real(r8), parameter :: D0_05 = 0.05_r8
55 : real(r8), parameter :: D0_25 = 0.25_r8
56 : real(r8), parameter :: D0_5 = 0.5_r8
57 : real(r8), parameter :: D1_0 = 1.0_r8
58 : real(r8), parameter :: D2_0 = 2.0_r8
59 : real(r8), parameter :: D3_0 = 3.0_r8
60 : real(r8), parameter :: D4_0 = 4.0_r8
61 : real(r8), parameter :: D8_0 = 8.0_r8
62 : real(r8), parameter :: D12_0 = 12.0_r8
63 : real(r8), parameter :: D24_0 = 24.0_r8
64 :
65 : CONTAINS
66 :
67 : !-----------------------------------------------------------------------
68 : !BOP
69 : ! !IROUTINE: tp2c --- Perform transport on a C grid
70 : !
71 : ! !INTERFACE:
72 42663936 : subroutine tp2c(dh, va, h, crx, cry, im, jm, &
73 42663936 : iord, jord, ng, fx, fy, ffsl, &
74 42663936 : rcap, acosp, xfx, yfx, cosp, id, jfirst, jlast)
75 : !-----------------------------------------------------------------------
76 :
77 : implicit none
78 :
79 : ! !INPUT PARAMETERS:
80 : integer im, jm ! Dimensions
81 : integer jfirst, jlast ! Latitude strip
82 : integer iord, jord ! Interpolation order in x,y
83 : integer ng ! Max. NS dependencies
84 : integer id ! density (0) (mfx = C)
85 : real (r8) rcap ! Ask S.-J. (polar constant?)
86 : real (r8) acosp(jm) ! Ask S.-J. (difference to cosp??)
87 : logical ffsl(jm) ! Use flux-form semi-Lagrangian trans.?
88 : ! (N*NG S*NG)
89 : real (r8) cosp(jm) ! Critical angle
90 : real (r8) va(im,jfirst:jlast) ! Courant (unghosted)
91 : real (r8) h(im,jfirst-ng:jlast+ng) ! Pressure ( N*NG S*NG )
92 : real (r8) crx(im,jfirst-ng:jlast+ng) ! Ask S.-J. ( N*NG S*NG )
93 : real (r8) cry(im,jfirst:jlast+1) ! Ask S.-J. ( N like FY )
94 : real (r8) xfx(im,jfirst:jlast) ! Ask S.-J. ( unghosted like FX )
95 : real (r8) yfx(im,jfirst:jlast+1) ! Ask S.-J. ( N like FY )
96 :
97 : ! !OUTPUT PARAMETERS:
98 : real (r8) dh(im,jfirst:jlast) ! Ask S.-J. ( unghosted )
99 : real (r8) fx(im,jfirst:jlast) ! Flux in x ( unghosted )
100 : real (r8) fy(im,jfirst:jlast+1) ! Flux in y ( N, see tp2c )
101 :
102 : ! !DESCRIPTION:
103 : ! Perform transport on a C grid. The number of ghost
104 : ! latitudes (NG) depends on what method (JORD) will be used
105 : ! subsequentally. NG is equal to MIN(ABS(JORD),3).
106 : ! Ask S.-J. how exactly this differs from TP2C.
107 : !
108 : ! !REVISION HISTORY:
109 : !
110 : !EOP
111 : !-----------------------------------------------------------------------
112 : !BOC
113 : integer i, j, js2g0, jn2g0
114 : real (r8) sum1
115 :
116 42663936 : js2g0 = max(2,jfirst) ! No ghosting
117 42663936 : jn2g0 = min(jm-1,jlast) ! No ghosting
118 :
119 : call tp2d(va, h, crx, cry, im, jm, iord, jord, ng,fx, fy, ffsl, &
120 42663936 : xfx, yfx, cosp, id, jfirst, jlast)
121 :
122 169322496 : do j=js2g0,jn2g0
123 36477665280 : do i=1,im-1
124 36477665280 : dh(i,j) = fx(i,j) - fx(i+1,j) + (fy(i,j)-fy(i,j+1))*acosp(j)
125 : enddo
126 169322496 : dh(im,j) = fx(im,j) - fx(1,j) + (fy(im,j)-fy(im,j+1))*acosp(j)
127 : enddo
128 :
129 : ! Poles
130 42663936 : if ( jfirst == 1 ) then
131 : ! sum1 = - SUM( fy(1:im, 2) ) * rcap
132 666624 : sum1 = D0_0
133 192654336 : do i=1,im
134 192654336 : sum1 = sum1 + fy(i,2)
135 : enddo
136 666624 : sum1 = -sum1*rcap
137 192654336 : do i=1,im
138 192654336 : dh(i, 1) = sum1
139 : enddo
140 : endif
141 :
142 42663936 : if ( jlast == jm ) then
143 : ! sum1 = SUM( fy(1:im,jm) ) * rcap
144 666624 : sum1 = D0_0
145 192654336 : do i=1,im
146 192654336 : sum1 = sum1 + fy(i,jm)
147 : enddo
148 666624 : sum1 = sum1*rcap
149 192654336 : do i=1,im
150 192654336 : dh(i,jm) = sum1
151 : enddo
152 : endif
153 42663936 : return
154 : !EOC
155 : end subroutine tp2c
156 : !-----------------------------------------------------------------------
157 :
158 : !-----------------------------------------------------------------------
159 : !BOP
160 : ! !IROUTINE: tp2d --- Perform transport on a D grid
161 : !
162 : ! !INTERFACE:
163 43008000 : subroutine tp2d(va, q, crx, cry, im, jm, iord, jord, ng, fx, fy, &
164 43008000 : ffsl, xfx, yfx, cosp, id, jfirst, jlast)
165 : !-----------------------------------------------------------------------
166 : ! !USES:
167 :
168 : implicit none
169 :
170 : ! !INPUT PARAMETERS:
171 : integer im, jm ! Dimensions
172 : integer jfirst, jlast ! Latitude strip
173 : integer iord, jord ! Interpolation order in x,y
174 : integer ng ! Max. NS dependencies
175 : integer id ! density (0) (mfx = C)
176 : ! mixing ratio (1) (mfx = mass flux)
177 : logical ffsl(jm) ! Use flux-form semi-Lagrangian trans.?
178 : ! ghosted N*ng S*ng
179 : real (r8) cosp(jm) ! Critical angle
180 : real (r8) va(im,jfirst:jlast) ! Courant (unghosted)
181 : real (r8) q(im,jfirst-ng:jlast+ng) ! transported scalar ( N*NG S*NG )
182 : real (r8) crx(im,jfirst-ng:jlast+ng) ! Ask S.-J. ( N*NG S*NG )
183 : real (r8) cry(im,jfirst:jlast+1) ! Ask S.-J. ( N like FY )
184 : real (r8) xfx(im,jfirst:jlast) ! Ask S.-J. ( unghosted like FX )
185 : real (r8) yfx(im,jfirst:jlast+1) ! Ask S.-J. ( N like FY )
186 :
187 : ! !OUTPUT PARAMETERS:
188 : real (r8) fx(im,jfirst:jlast) ! Flux in x ( unghosted )
189 : real (r8) fy(im,jfirst:jlast+1) ! Flux in y ( N, see tp2c )
190 :
191 : ! !DESCRIPTION:
192 : ! Perform transport on a D grid. The number of ghost
193 : ! latitudes (NG) depends on what method (JORD) will be used
194 : ! subsequentally. NG is equal to MIN(ABS(JORD),3).
195 : !
196 : !
197 : ! !REVISION HISTORY:
198 : ! WS 99.04.13: Added jfirst:jlast concept
199 : ! 99.04.21: Removed j1 and j2 (j1=2, j2=jm-1 consistently)
200 : ! 99.04.27: Removed dc, wk2 as arguments (local to YTP)
201 : ! 99.04.27: Removed adx as arguments (local here)
202 : ! SJL 99.07.26: ffsl flag added
203 : ! WS 99.09.07: Restructuring, cleaning, documentation
204 : ! WS 99.10.22: NG now argument; arrays pruned
205 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
206 : !
207 : !EOP
208 : !---------------------------------------------------------------------
209 : !BOC
210 :
211 : ! Local:
212 : integer i, j, iad, jp, js2g0, js2gng, jn2g0, jn2gng
213 86016000 : real (r8) adx(im,jfirst-ng:jlast+ng)
214 86016000 : real (r8) wk1v(im,jfirst-ng:jlast+ng)
215 86016000 : real (r8) dm(-im/3:im+im/3)
216 86016000 : real (r8) qtmpv(-im/3:im+im/3,jfirst-ng:jlast+ng)
217 86016000 : real (r8) al(-im/3:im+im/3)
218 86016000 : real (r8) ar(-im/3:im+im/3)
219 86016000 : real (r8) a6(-im/3:im+im/3)
220 :
221 : ! Number of ghost latitudes
222 43008000 : js2g0 = max(2,jfirst) ! No ghosting
223 43008000 : js2gng = max(2,jfirst-ng) ! Number needed on S
224 43008000 : jn2g0 = min(jm-1,jlast) ! No ghosting
225 43008000 : jn2gng = min(jm-1,jlast+ng) ! Number needed on N
226 43008000 : iad = 1
227 :
228 : call xtpv(im, ffsl, wk1v, q, crx, iad, crx, &
229 : cosp, 0, dm, qtmpv, al, ar, a6, &
230 : jfirst, jlast, js2gng, jn2gng, jm, &
231 : 1, jm, jfirst-ng, jlast+ng, &
232 : jfirst-ng, jlast+ng, jfirst-ng, jlast+ng, &
233 43008000 : jfirst-ng, jlast+ng, jfirst-ng, jlast+ng)
234 :
235 422026752 : do j=js2gng,jn2gng ! adx needed on N*ng S*ng
236 :
237 >10915*10^7 : do i=1,im-1
238 >21755*10^7 : adx(i,j) = q(i,j) + D0_5 * &
239 >32671*10^7 : (wk1v(i,j)-wk1v(i+1,j) + q(i,j)*(crx(i+1,j)-crx(i,j)))
240 : enddo
241 758037504 : adx(im,j) = q(im,j) + D0_5 * &
242 1180064256 : (wk1v(im,j)-wk1v(1,j) + q(im,j)*(crx(1,j)-crx(im,j)))
243 : enddo
244 :
245 : ! WS 99.09.07 : Split up north and south pole
246 :
247 43008000 : if ( jfirst-ng <= 1 ) then
248 385308672 : do i=1,im
249 385308672 : adx(i, 1) = q(i,1)
250 : enddo
251 : endif
252 43008000 : if ( jlast+ng >= jm ) then
253 385308672 : do i=1,im
254 385308672 : adx(i,jm) = q(i,jm)
255 : enddo
256 : endif
257 :
258 43008000 : call ytp(im,jm,fy, adx,cry,yfx,ng,jord,0,jfirst,jlast)
259 :
260 170688000 : do j=js2g0,jn2g0
261 36942528000 : do i=1,im
262 36771840000 : jp = j-va(i,j)
263 36899520000 : wk1v(i,j) = q(i,j) +D0_5*va(i,j)*(q(i,jp)-q(i,jp+1))
264 : enddo
265 : enddo
266 :
267 : call xtpv(im, ffsl, fx, wk1v, crx, iord, xfx, &
268 : cosp, id, dm, qtmpv, al, ar, a6, &
269 : jfirst, jlast, js2g0, jn2g0, jm, &
270 : 1, jm, jfirst, jlast, &
271 : jfirst-ng, jlast+ng, jfirst-ng, jlast+ng, &
272 43008000 : jfirst, jlast, jfirst-ng, jlast+ng)
273 :
274 43008000 : return
275 : !EOC
276 : end subroutine tp2d
277 : !-----------------------------------------------------------------------
278 :
279 : #ifndef VECTORIZE
280 : !-----------------------------------------------------------------------
281 : !BOP
282 : ! !IROUTINE: xtpv
283 : !
284 : ! !INTERFACE:
285 87392256 : subroutine xtpv(im, ffslv, fxv, qv, cv, iord, mfxv, &
286 87392256 : cosav, id, dmw, qtmpv, alw, arw, a6w, &
287 : jfirst, jlast, jlow, jhigh, jm, &
288 : jl2, jh2, jl3, jh3, &
289 : jl4, jh4, jl5, jh5, &
290 : jl7, jh7, jl11, jh11)
291 : !-----------------------------------------------------------------------
292 :
293 : implicit none
294 :
295 : ! !INPUT PARAMETERS:
296 : integer id ! ID = 0: density (mfx = C)
297 : ! ID = 1: mixing ratio (mfx is mass flux)
298 :
299 : integer im ! Total longitudes
300 : integer iord
301 : integer jfirst, jlast, jlow, jhigh, jm
302 : integer jl2, jh2, jl3, jh3, jl4, jh4, jl5, jh5
303 : integer jl7, jh7, jl11, jh11
304 : real (r8) cv(im,jl5:jh5) ! Courant numbers
305 : real (r8) qv(im,jl4:jh4)
306 : real (r8) mfxv(im,jl7:jh7)
307 : logical ffslv(jl2:jh2)
308 : real (r8) cosav(jm)
309 :
310 : ! !INPUT/OUTPUT PARAMETERS:
311 : real (r8) qtmpv(-im/3:im+im/3,jl11:jh11) ! Input work arrays:
312 : real (r8) dmw(-im/3:im+im/3)
313 : real (r8) alw(-im/3:im+im/3)
314 : real (r8) arw(-im/3:im+im/3)
315 : real (r8) a6w(-im/3:im+im/3)
316 :
317 : ! !OUTPUT PARAMETERS:
318 : real (r8) fxv(im,jl3:jh3)
319 :
320 : ! !DESCRIPTION:
321 : !
322 : !
323 : ! !REVISION HISTORY:
324 : ! 99.01.01 Lin Creation
325 : ! 01.03.27 Sawyer Additional ProTeX documentation
326 : !
327 : !EOP
328 : !-----------------------------------------------------------------------
329 : !BOC
330 :
331 : ! Local:
332 : real (r8) cos_upw !critical cosine for upwind
333 : real (r8) cos_van !critical cosine for van Leer
334 : real (r8) cos_ppm !critical cosine for ppm
335 :
336 : parameter (cos_upw = D0_05) !roughly at 87 deg.
337 : parameter (cos_van = D0_25) !roughly at 75 deg.
338 : parameter (cos_ppm = D0_25)
339 :
340 : integer i, imp, j
341 : real (r8) qmax, qmin
342 : real (r8) rut, tmp
343 : integer iu, itmp, ist
344 174784512 : integer isave(im)
345 : integer iuw, iue
346 174784512 : real (r8) dm(-im/3:im+im/3)
347 174784512 : real (r8) al(-im/3:im+im/3)
348 174784512 : real (r8) ar(-im/3:im+im/3)
349 174784512 : real (r8) a6(-im/3:im+im/3)
350 :
351 87392256 : imp = im + 1
352 :
353 599801664 : do j = jlow, jhigh
354 :
355 >14808*10^7 : do i=1,im
356 >14808*10^7 : qtmpv(i,j) = qv(i,j)
357 : enddo
358 :
359 599801664 : if( ffslv(j) ) then
360 : ! Flux-Form Semi-Lagrangian transport
361 :
362 : ! Figure out ghost zone for the western edge:
363 15378664 : iuw = -cv(1,j)
364 15378664 : iuw = min(0, iuw)
365 :
366 30757328 : do i=iuw, 0
367 30757328 : qtmpv(i,j) = qv(im+i,j)
368 : enddo
369 :
370 : ! Figure out ghost zone for the eastern edge:
371 15378664 : iue = im - cv(im,j)
372 15378664 : iue = max(imp, iue)
373 :
374 33246870 : do i=imp, iue
375 33246870 : qtmpv(i,j) = qv(i-im,j)
376 : enddo
377 :
378 15378664 : if( iord == 1 .or. cosav(j) < cos_upw) then
379 3926388391 : do i=1,im
380 3912802272 : iu = cv(i,j)
381 3912802272 : if(cv(i,j) .le. D0_0) then
382 1558489343 : itmp = i - iu
383 1558489343 : isave(i) = itmp - 1
384 : else
385 2354312929 : itmp = i - iu - 1
386 2354312929 : isave(i) = itmp + 1
387 : endif
388 3926388391 : fxv(i,j) = (cv(i,j)-iu) * qtmpv(itmp,j)
389 : enddo
390 : else
391 :
392 518045505 : do i=1,im
393 : ! 2nd order slope
394 516252960 : tmp = D0_25*(qtmpv(i+1,j) - qtmpv(i-1,j))
395 516252960 : qmax = max(qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j)) - qtmpv(i,j)
396 516252960 : qmin = qtmpv(i,j) - min(qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j))
397 518045505 : dm(i) = sign(min(abs(tmp),qmax,qmin), tmp)
398 : enddo
399 :
400 :
401 3585090 : do i=iuw, 0
402 3585090 : dm(i) = dm(im+i)
403 : enddo
404 :
405 3585090 : do i=imp, iue
406 3585090 : dm(i) = dm(i-im)
407 : enddo
408 :
409 1792545 : if(iord .ge. 3 .and. cosav(j) .gt. cos_ppm) then
410 : call fxppm(im, cv(:,j), mfxv(:,j), qtmpv(:,j), dm, fxv(:,j), iord, al, ar, a6, &
411 262560 : iuw, iue, ffslv(j), isave)
412 : else
413 442165665 : do i=1,im
414 440635680 : iu = cv(i,j)
415 440635680 : rut = cv(i,j) - iu
416 442165665 : if(cv(i,j) .le. D0_0) then
417 124549908 : itmp = i - iu
418 124549908 : isave(i) = itmp - 1
419 124549908 : fxv(i,j) = rut*(qtmpv(itmp,j)-dm(itmp)*(D1_0+rut))
420 : else
421 316085772 : itmp = i - iu - 1
422 316085772 : isave(i) = itmp + 1
423 316085772 : fxv(i,j) = rut*(qtmpv(itmp,j)+dm(itmp)*(D1_0-rut))
424 : endif
425 : enddo
426 : endif
427 :
428 : endif
429 :
430 4444433896 : do i=1,im
431 4444433896 : if(cv(i,j) .ge. D1_0) then
432 2863285989 : do ist = isave(i),i-1
433 2863285989 : fxv(i,j) = fxv(i,j) + qtmpv(ist,j)
434 : enddo
435 3352868707 : elseif(cv(i,j) .le. -D1_0) then
436 1632720550 : do ist = i,isave(i)
437 1632720550 : fxv(i,j) = fxv(i,j) - qtmpv(ist,j)
438 : enddo
439 : endif
440 : enddo
441 :
442 15378664 : if(id .ne. 0) then
443 1241335920 : do i=1,im
444 1241335920 : fxv(i,j) = fxv(i,j)*mfxv(i,j)
445 : enddo
446 : endif
447 :
448 : else
449 : ! Regular PPM (Eulerian without FFSL extension)
450 :
451 497030744 : qtmpv(imp,j) = qv(1,j)
452 497030744 : qtmpv( 0,j) = qv(im,j)
453 :
454 497030744 : if(iord == 1 .or. cosav(j) < cos_upw) then
455 >10756*10^7 : do i=1,im
456 >10719*10^7 : iu = real(i,r8) - cv(i,j)
457 >10756*10^7 : fxv(i,j) = mfxv(i,j)*qtmpv(iu,j)
458 : enddo
459 : else
460 :
461 124819311 : qtmpv(-1,j) = qv(im-1,j)
462 124819311 : qtmpv(imp+1,j) = qv(2,j)
463 :
464 124819311 : if(iord > 0 .or. cosav(j) < cos_van) then
465 121554735 : call xmist(im, qtmpv(:,j), dm, 2)
466 : else
467 3264576 : call xmist(im, qtmpv(:,j), dm, iord)
468 : endif
469 :
470 124819311 : dm(0) = dm(im)
471 :
472 124819311 : if( abs(iord).eq.2 .or. cosav(j) .lt. cos_van ) then
473 5374849455 : do i=1,im
474 5356251360 : iu = real(i,r8) - cv(i,j)
475 5374849455 : fxv(i,j) = mfxv(i,j)*(qtmpv(iu,j)+dm(iu)*(sign(D1_0,cv(i,j))-cv(i,j)))
476 :
477 : ! if(cv(i,j) .le. 0.) then
478 : ! fxv(i,j) = qtmpv(i,j) - dm(i)*(1.+cv(i,j))
479 : ! else
480 : ! fxv(i,j) = qtmpv(i-1,j) + dm(i-1)*(1.-cv(i,j))
481 : ! endif
482 : ! fxv(i,j) = fxv(i,j)*mfxv(i,j)
483 :
484 : enddo
485 : else
486 : call fxppm(im, cv(:,j), mfxv(:,j), qtmpv(:,j), dm, fxv(:,j), iord, al, ar, a6, &
487 106221216 : iuw, iue, ffslv(j), isave)
488 : endif
489 : endif
490 :
491 : endif
492 :
493 : enddo
494 :
495 87392256 : return
496 : !EOC
497 : end subroutine xtpv
498 : !-----------------------------------------------------------------------
499 :
500 : !-----------------------------------------------------------------------
501 : !BOP
502 : ! !IROUTINE: xmist
503 : !
504 : ! !INTERFACE:
505 124819311 : subroutine xmist(im, q, dm, id)
506 : !-----------------------------------------------------------------------
507 :
508 : implicit none
509 :
510 : ! !INPUT PARAMETERS:
511 : integer im ! Total number of longitudes
512 : integer id ! ID = 0: density (mfx = C)
513 : ! ID = 1: mixing ratio (mfx is mass flux)
514 : real(r8) q(-im/3:im+im/3) ! Input latitude
515 :
516 : ! !OUTPUT PARAMETERS:
517 : real(r8) dm(-im/3:im+im/3) !
518 :
519 : ! !DESCRIPTION:
520 : !
521 : !
522 : ! !REVISION HISTORY:
523 : ! 99.01.01 Lin Creation
524 : ! 01.03.27 Sawyer Additional ProTeX documentation
525 : !
526 : !EOP
527 : !-----------------------------------------------------------------------
528 : !BOC
529 :
530 : real(r8) r24
531 : parameter( r24 = D1_0/D24_0)
532 :
533 : integer i
534 : real(r8) qmin, qmax
535 :
536 124819311 : if(id .le. 2) then
537 36072780879 : do i=1,im
538 36072780879 : dm(i) = r24*(D8_0*(q(i+1) - q(i-1)) + q(i-2) - q(i+2))
539 : enddo
540 : else
541 0 : do i=1,im
542 0 : dm(i) = D0_25*(q(i+1) - q(i-1))
543 : enddo
544 : endif
545 :
546 124819311 : if( id < 0 ) return
547 :
548 : ! Apply monotonicity constraint (Lin et al. 1994, MWR)
549 35129318415 : do i=1,im
550 35007763680 : qmax = max( q(i-1), q(i), q(i+1) ) - q(i)
551 35007763680 : qmin = q(i) - min( q(i-1), q(i), q(i+1) )
552 35129318415 : dm(i) = sign( min(abs(dm(i)), qmax, qmin), dm(i) )
553 : enddo
554 : return
555 : !EOC
556 : end subroutine xmist
557 : !-----------------------------------------------------------------------
558 :
559 : !-----------------------------------------------------------------------
560 : !BOP
561 : ! !IROUTINE: fxppm
562 : !
563 : ! !INTERFACE:
564 106483776 : subroutine fxppm(im, c, mfx, p, dm, fx, iord, al, ar, a6, &
565 106483776 : iuw, iue, ffsl, isave)
566 : !-----------------------------------------------------------------------
567 : !
568 : ! !USES:
569 : implicit none
570 :
571 : ! !INPUT PARAMETERS:
572 : integer im, iord
573 : real (r8) c(im)
574 : real (r8) p(-im/3:im+im/3)
575 : real (r8) dm(-im/3:im+im/3)
576 : real (r8) mfx(im)
577 : integer iuw, iue
578 : logical ffsl
579 : integer isave(im)
580 :
581 : ! !INPUT/OUTPUT PARAMETERS:
582 : real (r8) al(-im/3:im+im/3)
583 : real (r8) ar(-im/3:im+im/3)
584 : real (r8) a6(-im/3:im+im/3)
585 :
586 : ! !OUTPUT PARAMETERS:
587 : real (r8) fx(im)
588 :
589 : ! !DESCRIPTION:
590 : !
591 : !
592 : ! !REVISION HISTORY:
593 : ! 99.01.01 Lin Creation
594 : ! 01.03.27 Sawyer Additional ProTeX documentation
595 : !
596 : !EOP
597 : !-----------------------------------------------------------------------
598 : !BOC
599 : !
600 : ! !LOCAL VARIABLES:
601 : real (r8) r3, r23
602 : parameter ( r3 = D1_0/D3_0, r23 = D2_0/D3_0 )
603 :
604 : integer i, lmt
605 : integer iu, itmp
606 : real (r8) ru
607 : logical steep
608 :
609 106483776 : if( iord == 6 ) then
610 : steep = .true.
611 : else
612 106483776 : steep = .false.
613 : endif
614 :
615 30773811264 : do i=1,im
616 30773811264 : al(i) = D0_5*(p(i-1)+p(i)) + (dm(i-1) - dm(i))*r3
617 : enddo
618 :
619 106483776 : if( steep ) call steepx( im, p, al(1), dm )
620 :
621 30667327488 : do i=1,im-1
622 30667327488 : ar(i) = al(i+1)
623 : enddo
624 106483776 : ar(im) = al(1)
625 :
626 106483776 : if(iord == 7) then
627 : call huynh(im, ar(1), al(1), p(1), a6(1), dm(1))
628 : else
629 106483776 : if(iord .eq. 3 .or. iord .eq. 5) then
630 0 : do i=1,im
631 0 : a6(i) = D3_0*(p(i)+p(i) - (al(i)+ar(i)))
632 : enddo
633 : endif
634 106483776 : lmt = iord - 3
635 : call lmppm( dm(1), a6(1), ar(1), al(1), p(1), im, lmt )
636 : endif
637 :
638 106483776 : if( ffsl ) then
639 :
640 525120 : do i=iuw, 0
641 262560 : al(i) = al(im+i)
642 262560 : ar(i) = ar(im+i)
643 525120 : a6(i) = a6(im+i)
644 : enddo
645 :
646 525120 : do i=im+1, iue
647 262560 : al(i) = al(i-im)
648 262560 : ar(i) = ar(i-im)
649 525120 : a6(i) = a6(i-im)
650 : enddo
651 :
652 75879840 : do i=1,im
653 75617280 : iu = c(i)
654 75617280 : ru = c(i) - iu
655 75879840 : if(c(i) .gt. D0_0) then
656 75500160 : itmp = i - iu - 1
657 75500160 : isave(i) = itmp + 1
658 75500160 : fx(i) = ru*(ar(itmp)+D0_5*ru*(al(itmp)-ar(itmp) + &
659 151000320 : a6(itmp)*(D1_0-r23*ru)) )
660 : else
661 117120 : itmp = i - iu
662 117120 : isave(i) = itmp - 1
663 117120 : fx(i) = ru*(al(itmp)-D0_5*ru*(ar(itmp)-al(itmp) + &
664 234240 : a6(itmp)*(D1_0+r23*ru)) )
665 : endif
666 : enddo
667 :
668 : else
669 106221216 : al(0) = al(im)
670 106221216 : ar(0) = ar(im)
671 106221216 : a6(0) = a6(im)
672 30697931424 : do i=1,im
673 30591710208 : if(c(i) .gt. D0_0) then
674 18555299150 : fx(i) = ar(i-1) + D0_5*c(i)*(al(i-1) - ar(i-1) + &
675 18555299150 : a6(i-1)*(D1_0-r23*c(i)) )
676 : else
677 12036411058 : fx(i) = al(i) - D0_5*c(i)*(ar(i) - al(i) + &
678 12036411058 : a6(i)*(D1_0+r23*c(i)))
679 : endif
680 30697931424 : fx(i) = mfx(i) * fx(i)
681 : enddo
682 : endif
683 106483776 : return
684 : !EOC
685 : end subroutine fxppm
686 : !-----------------------------------------------------------------------
687 :
688 : !-----------------------------------------------------------------------
689 : !BOP
690 : ! !IROUTINE: steepx
691 : !
692 : ! !INTERFACE:
693 0 : subroutine steepx(im, p, al, dm)
694 : !-----------------------------------------------------------------------
695 :
696 : ! !USES:
697 : implicit none
698 :
699 : ! !INPUT PARAMETERS:
700 : integer im
701 : real (r8) p(-im/3:im+im/3)
702 : real (r8) dm(-im/3:im+im/3)
703 :
704 : ! !INPUT/OUTPUT PARAMETERS:
705 : real (r8) al(im)
706 :
707 : ! !DESCRIPTION:
708 : !
709 : !
710 : ! !REVISION HISTORY:
711 : ! 99.01.01 Lin Creation
712 : ! 01.03.27 Sawyer Additional ProTeX documentation
713 : !
714 : !EOP
715 : !-----------------------------------------------------------------------
716 : !BOC
717 : !
718 : ! !LOCAL VARIABLES:
719 : integer i
720 : real (r8) r3
721 : parameter ( r3 = D1_0/D3_0 )
722 :
723 0 : real (r8) dh(0:im)
724 0 : real (r8) d2(0:im+1)
725 0 : real (r8) eta(0:im)
726 : real (r8) xxx, bbb, ccc
727 :
728 0 : do i=0,im
729 0 : dh(i) = p(i+1) - p(i)
730 : enddo
731 :
732 : ! Needs dh(0:im)
733 0 : do i=1,im
734 0 : d2(i) = dh(i) - dh(i-1)
735 : enddo
736 0 : d2(0) = d2(im)
737 0 : d2(im+1) = d2(1)
738 :
739 : ! needs p(-1:im+2), d2(0:im+1)
740 0 : do i=1,im
741 0 : if( d2(i+1)*d2(i-1).lt.D0_0 .and. p(i+1).ne.p(i-1) ) then
742 0 : xxx = D1_0 - D0_5 * ( p(i+2) - p(i-2) ) / ( p(i+1) - p(i-1) )
743 0 : eta(i) = max(D0_0, min(xxx, D0_5) )
744 : else
745 0 : eta(i) = D0_0
746 : endif
747 : enddo
748 :
749 0 : eta(0) = eta(im)
750 :
751 : ! needs eta(0:im), dh(0:im-1), dm(0:im)
752 0 : do i=1,im
753 0 : bbb = ( D2_0*eta(i ) - eta(i-1) ) * dm(i-1)
754 0 : ccc = ( D2_0*eta(i-1) - eta(i ) ) * dm(i )
755 0 : al(i) = al(i) + D0_5*( eta(i-1) - eta(i)) * dh(i-1) + (bbb - ccc) * r3
756 : enddo
757 0 : return
758 : !EOC
759 : end subroutine steepx
760 : !-----------------------------------------------------------------------
761 :
762 : !-----------------------------------------------------------------------
763 : !BOP
764 : ! !IROUTINE: lmppm
765 : !
766 : ! !INTERFACE:
767 148975680 : subroutine lmppm(dm, a6, ar, al, p, im, lmt)
768 : !-----------------------------------------------------------------------
769 :
770 : implicit none
771 :
772 : ! !INPUT PARAMETERS:
773 : integer im ! Total longitudes
774 : integer lmt ! LMT = 0: full monotonicity
775 : ! LMT = 1: Improved and simplified full monotonic constraint
776 : ! LMT = 2: positive-definite constraint
777 : ! LMT = 3: Quasi-monotone constraint
778 : real(r8) p(im)
779 : real(r8) dm(im)
780 :
781 : ! !OUTPUT PARAMETERS:
782 : real(r8) a6(im)
783 : real(r8) ar(im)
784 : real(r8) al(im)
785 :
786 : ! !DESCRIPTION:
787 : !
788 : !
789 : ! !REVISION HISTORY:
790 : ! 99.01.01 Lin Creation
791 : ! 01.03.27 Sawyer Additional ProTeX documentation
792 : !
793 : !EOP
794 : !-----------------------------------------------------------------------
795 : !BOC
796 : !
797 : ! !LOCAL VARIABLES:
798 : real (r8) r12
799 : parameter ( r12 = D1_0/D12_0 )
800 :
801 : real (r8) da1, da2, fmin, a6da
802 : real (r8) dr, dl
803 :
804 : integer i
805 :
806 : ! LMT = 0: full monotonicity
807 : ! LMT = 1: Improved and simplified full monotonic constraint
808 : ! LMT = 2: positive-definite constraint
809 : ! LMT = 3: Quasi-monotone constraint
810 :
811 148975680 : if( lmt == 0 ) then
812 :
813 : ! Full constraint
814 0 : do i=1,im
815 0 : if(dm(i) .eq. D0_0) then
816 0 : ar(i) = p(i)
817 0 : al(i) = p(i)
818 0 : a6(i) = D0_0
819 : else
820 0 : da1 = ar(i) - al(i)
821 0 : da2 = da1**2
822 0 : a6da = a6(i)*da1
823 0 : if(a6da .lt. -da2) then
824 0 : a6(i) = D3_0*(al(i)-p(i))
825 0 : ar(i) = al(i) - a6(i)
826 0 : elseif(a6da .gt. da2) then
827 0 : a6(i) = D3_0*(ar(i)-p(i))
828 0 : al(i) = ar(i) - a6(i)
829 : endif
830 : endif
831 : enddo
832 :
833 148975680 : elseif( lmt == 1 ) then
834 :
835 : ! Improved (Lin 2001?) full constraint
836 91622217792 : do i=1,im
837 91473242112 : da1 = dm(i) + dm(i)
838 91473242112 : dl = sign(min(abs(da1),abs(al(i)-p(i))), da1)
839 91473242112 : dr = sign(min(abs(da1),abs(ar(i)-p(i))), da1)
840 91473242112 : ar(i) = p(i) + dr
841 91473242112 : al(i) = p(i) - dl
842 91622217792 : a6(i) = D3_0*(dl-dr)
843 : enddo
844 :
845 0 : elseif( lmt == 2 ) then
846 : ! Positive definite constraint
847 0 : do 250 i=1,im
848 0 : if(abs(ar(i)-al(i)) .ge. -a6(i)) go to 250
849 0 : fmin = p(i) + D0_25*(ar(i)-al(i))**2/a6(i) + a6(i)*r12
850 0 : if(fmin.ge.D0_0) go to 250
851 0 : if(p(i).lt.ar(i) .and. p(i).lt.al(i)) then
852 0 : ar(i) = p(i)
853 0 : al(i) = p(i)
854 0 : a6(i) = D0_0
855 0 : elseif(ar(i) .gt. al(i)) then
856 0 : a6(i) = D3_0*(al(i)-p(i))
857 0 : ar(i) = al(i) - a6(i)
858 : else
859 0 : a6(i) = D3_0*(ar(i)-p(i))
860 0 : al(i) = ar(i) - a6(i)
861 : endif
862 0 : 250 continue
863 :
864 0 : elseif(lmt .eq. 3) then
865 : ! Quasi-monotone constraint
866 0 : do i=1,im
867 0 : da1 = D4_0*dm(i)
868 0 : dl = sign(min(abs(da1),abs(al(i)-p(i))), da1)
869 0 : dr = sign(min(abs(da1),abs(ar(i)-p(i))), da1)
870 0 : ar(i) = p(i) + dr
871 0 : al(i) = p(i) - dl
872 0 : a6(i) = D3_0*(dl-dr)
873 : enddo
874 : endif
875 148975680 : return
876 : !EOC
877 : end subroutine lmppm
878 : !-----------------------------------------------------------------------
879 :
880 : !-----------------------------------------------------------------------
881 : !BOP
882 : ! !IROUTINE: huynh --- Enforce Huynh's 2nd constraint in 1D periodic domain
883 : !
884 : ! !INTERFACE:
885 0 : subroutine huynh(im, ar, al, p, d2, d1)
886 : !-----------------------------------------------------------------------
887 :
888 : ! !USES:
889 :
890 : implicit none
891 :
892 : ! !INPUT PARAMETERS:
893 : integer im
894 : real(r8) p(im)
895 :
896 : ! !OUTPUT PARAMETERS:
897 : real(r8) ar(im)
898 : real(r8) al(im)
899 : real(r8) d2(im)
900 : real(r8) d1(im)
901 :
902 : ! !DESCRIPTION:
903 : !
904 : ! Enforce Huynh's 2nd constraint in 1D periodic domain
905 : !
906 : ! !REVISION HISTORY:
907 : ! 99.01.01 Lin Creation
908 : ! 01.03.27 Sawyer Additional ProTeX documentation
909 : !
910 : !EOP
911 : !-----------------------------------------------------------------------
912 : !BOC
913 : !
914 : ! !LOCAL VARIABLES:
915 : integer i
916 : real(r8) pmp
917 : real(r8) lac
918 : real(r8) pmin
919 : real(r8) pmax
920 :
921 : ! Compute d1 and d2
922 0 : d1(1) = p(1) - p(im)
923 0 : do i=2,im
924 0 : d1(i) = p(i) - p(i-1)
925 : enddo
926 :
927 0 : do i=1,im-1
928 0 : d2(i) = d1(i+1) - d1(i)
929 : enddo
930 0 : d2(im) = d1(1) - d1(im)
931 :
932 : ! Constraint for AR
933 : ! i = 1
934 0 : pmp = p(1) + D2_0 * d1(1)
935 0 : lac = p(1) + D0_5 * (d1(1)+d2(im)) + d2(im)
936 0 : pmin = min(p(1), pmp, lac)
937 0 : pmax = max(p(1), pmp, lac)
938 0 : ar(1) = min(pmax, max(ar(1), pmin))
939 :
940 0 : do i=2, im
941 0 : pmp = p(i) + D2_0*d1(i)
942 0 : lac = p(i) + D0_5*(d1(i)+d2(i-1)) + d2(i-1)
943 0 : pmin = min(p(i), pmp, lac)
944 0 : pmax = max(p(i), pmp, lac)
945 0 : ar(i) = min(pmax, max(ar(i), pmin))
946 : enddo
947 :
948 : ! Constraint for AL
949 0 : do i=1, im-1
950 0 : pmp = p(i) - D2_0*d1(i+1)
951 0 : lac = p(i) + D0_5*(d2(i+1)-d1(i+1)) + d2(i+1)
952 0 : pmin = min(p(i), pmp, lac)
953 0 : pmax = max(p(i), pmp, lac)
954 0 : al(i) = min(pmax, max(al(i), pmin))
955 : enddo
956 :
957 : ! i=im
958 0 : i = im
959 0 : pmp = p(im) - D2_0*d1(1)
960 0 : lac = p(im) + D0_5*(d2(1)-d1(1)) + d2(1)
961 0 : pmin = min(p(im), pmp, lac)
962 0 : pmax = max(p(im), pmp, lac)
963 0 : al(im) = min(pmax, max(al(im), pmin))
964 :
965 : ! compute A6 (d2)
966 0 : do i=1, im
967 0 : d2(i) = D3_0*(p(i)+p(i) - (al(i)+ar(i)))
968 : enddo
969 0 : return
970 : !EOC
971 : end subroutine huynh
972 : !-----------------------------------------------------------------------
973 : #endif
974 :
975 : !-----------------------------------------------------------------------
976 : !BOP
977 : ! !IROUTINE: ytp
978 : !
979 : ! !INTERFACE:
980 43352064 : subroutine ytp(im, jm, fy, q, c, yfx, ng, jord, iv, jfirst, jlast)
981 : !-----------------------------------------------------------------------
982 :
983 : ! !USES:
984 : implicit none
985 :
986 : ! !INPUT PARAMETERS:
987 : integer im, jm ! Dimensions
988 : integer jfirst, jlast ! Latitude strip
989 : integer ng ! Max. NS dependencies
990 : integer jord ! order of subgrid dist
991 : integer iv ! Scalar=0, Vector=1
992 : real (r8) q(im,jfirst-ng:jlast+ng) ! advected scalar N*jord S*jord
993 : real (r8) c(im,jfirst:jlast+1) ! Courant N (like FY)
994 : real (r8) yfx(im,jfirst:jlast+1) ! Backgrond mass flux
995 :
996 : ! !OUTPUT PARAMETERS:
997 : real (r8) fy(im,jfirst:jlast+1) ! Flux N (see tp2c)
998 :
999 : ! !DESCRIPTION:
1000 : ! This routine calculates the flux FX. The method chosen
1001 : ! depends on the order of the calculation JORD (currently
1002 : ! 1, 2 or 3).
1003 : !
1004 : ! !CALLED FROM:
1005 : ! cd_core
1006 : ! tp2d
1007 : !
1008 : ! !REVISION HISTORY:
1009 : !
1010 : ! SJL 99.04.13: Delivery
1011 : ! WS 99.04.13: Added jfirst:jlast concept
1012 : ! WS 99.04.21: Removed j1 and j2 (j1=2, j2=jm-1 consistently)
1013 : ! removed a6,ar,al from argument list
1014 : ! WS 99.04.27: DM made local to this routine
1015 : ! WS 99.09.09: Documentation; indentation; cleaning
1016 : ! WS 99.10.22: Added NG as argument; pruned arrays
1017 : ! SJL 99.12.24: Revised documentation; optimized for better cache usage
1018 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
1019 : !
1020 : !EOP
1021 : !---------------------------------------------------------------------
1022 : !BOC
1023 : !
1024 : ! !LOCAL VARIABLES:
1025 : integer i, j, jt
1026 : integer js2g0, jn1g1
1027 :
1028 : ! work arrays (should pass in eventually for performance enhancement):
1029 86704128 : real (r8) dm(im,jfirst-ng:jlast+ng)
1030 :
1031 : ! real (r8) ar(im,jfirst-1:jlast+1) ! AR needs to be ghosted on NS
1032 : ! real (r8) al(im,jfirst-1:jlast+2) ! AL needs to be ghosted on N2S
1033 : ! real (r8) a6(im,jfirst-1:jlast+1) ! A6 needs to be ghosted on NS
1034 :
1035 43352064 : js2g0 = max(2,jfirst) ! No ghosting
1036 43352064 : jn1g1 = min(jm,jlast+1) ! Ghost N*1
1037 :
1038 43352064 : if(jord == 1) then
1039 641088 : do j=js2g0,jn1g1
1040 148115520 : do i=1,im
1041 147474432 : jt = real(j,r8) - c(i,j)
1042 147986496 : fy(i,j) = q(i,jt)
1043 : enddo
1044 : enddo
1045 : else
1046 :
1047 : !
1048 : ! YMIST requires q on NS; Only call to YMIST here
1049 : !
1050 43223040 : call ymist(im, jm, q, dm, ng, jord, iv, jfirst, jlast)
1051 :
1052 43223040 : if( abs(jord) .ge. 3 ) then
1053 :
1054 42491904 : call fyppm(c,q,dm,fy,im,jm,ng,jord,iv,jfirst,jlast)
1055 :
1056 : else
1057 : !
1058 : ! JORD can either have the value 2 or -2 at this point
1059 : !
1060 3632832 : do j=js2g0,jn1g1
1061 839321280 : do i=1,im
1062 835688448 : jt = real(j,r8) - c(i,j)
1063 838590144 : fy(i,j) = q(i,jt) + (sign(D1_0,c(i,j))-c(i,j))*dm(i,jt)
1064 : enddo
1065 : enddo
1066 : endif
1067 : endif
1068 :
1069 215405568 : do j=js2g0,jn1g1
1070 49766814720 : do i=1,im
1071 49723462656 : fy(i,j) = fy(i,j)*yfx(i,j)
1072 : enddo
1073 : enddo
1074 43352064 : return
1075 : !EOC
1076 : end subroutine ytp
1077 : !-----------------------------------------------------------------------
1078 :
1079 : !-----------------------------------------------------------------------
1080 : !BOP
1081 : ! !IROUTINE: ymist
1082 : !
1083 : ! !INTERFACE:
1084 43223040 : subroutine ymist(im, jm, q, dm, ng, jord, iv, jfirst, jlast)
1085 : !-----------------------------------------------------------------------
1086 :
1087 : ! !USES:
1088 : implicit none
1089 :
1090 : ! !INPUT PARAMETERS:
1091 : integer im, jm ! Dimensions
1092 : integer jfirst, jlast ! Latitude strip
1093 : integer ng ! NS dependencies
1094 : integer jord ! order of subgrid distribution
1095 : integer iv ! Scalar (==0) Vector (==1)
1096 : real (r8) q(im,jfirst-ng:jlast+ng) ! transported scalar N*ng S*ng
1097 :
1098 : ! !OUTPUT PARAMETERS:
1099 : real (r8) dm(im,jfirst-ng:jlast+ng) ! Slope only N*(ng-1) S*(ng-1) used
1100 :
1101 : ! !DESCRIPTION:
1102 : ! Calculate the slope of the pressure. The number of ghost
1103 : ! latitudes (NG) depends on what method (JORD) will be used
1104 : ! subsequentally. NG is equal to MIN(ABS(JORD),3).
1105 : !
1106 : ! !CALLED FROM:
1107 : ! ytp
1108 : !
1109 : ! !REVISION HISTORY:
1110 : !
1111 : ! SJL 99.04.13: Delivery
1112 : ! WS 99.04.13: Added jfirst:jlast concept
1113 : ! WS 99.09.09: Documentation; indentation; cleaning
1114 : ! SJL 00.01.06: Documentation
1115 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
1116 : !
1117 : !EOP
1118 : !---------------------------------------------------------------------
1119 : !BOC
1120 :
1121 : ! Local variables
1122 :
1123 : integer i, j, jm1, im2, js2gng1, jn2gng1
1124 : real (r8) qmax, qmin, tmp
1125 :
1126 43223040 : js2gng1 = max(2, jfirst-ng+1) ! Number needed on S
1127 43223040 : jn2gng1 = min(jm-1,jlast+ng-1) ! Number needed on N
1128 :
1129 43223040 : jm1 = jm - 1
1130 43223040 : im2 = im / 2
1131 :
1132 340546752 : do j=js2gng1,jn2gng1
1133 85969775808 : do i=1,im
1134 85926552768 : dm(i,j) = D0_25*(q(i,j+1) - q(i,j-1))
1135 : enddo
1136 : enddo
1137 :
1138 43223040 : if( iv == 0 ) then
1139 :
1140 42889728 : if ( jfirst-ng <= 1 ) then
1141 : ! S pole
1142 192979920 : do i=1,im2
1143 191649024 : tmp = D0_25*(q(i,2)-q(i+im2,2))
1144 191649024 : qmax = max(q(i,2),q(i,1), q(i+im2,2)) - q(i,1)
1145 191649024 : qmin = q(i,1) - min(q(i,2),q(i,1), q(i+im2,2))
1146 192979920 : dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp)
1147 : enddo
1148 :
1149 192979920 : do i=im2+1,im
1150 192979920 : dm(i, 1) = - dm(i-im2, 1)
1151 : enddo
1152 : endif
1153 :
1154 42889728 : if ( jlast+ng >= jm ) then
1155 : ! N pole
1156 192979920 : do i=1,im2
1157 191649024 : tmp = D0_25*(q(i+im2,jm1)-q(i,jm1))
1158 191649024 : qmax = max(q(i+im2,jm1),q(i,jm), q(i,jm1)) - q(i,jm)
1159 191649024 : qmin = q(i,jm) - min(q(i+im2,jm1),q(i,jm), q(i,jm1))
1160 192979920 : dm(i,jm) = sign(min(abs(tmp),qmax,qmin),tmp)
1161 : enddo
1162 :
1163 192979920 : do i=im2+1,im
1164 192979920 : dm(i,jm) = - dm(i-im2,jm)
1165 : enddo
1166 : endif
1167 :
1168 : else
1169 :
1170 333312 : if ( jfirst-ng <= 1 ) then
1171 : ! South
1172 1510320 : do i=1,im2
1173 1499904 : tmp = D0_25*(q(i,2)+q(i+im2,2))
1174 1499904 : qmax = max(q(i,2),q(i,1), -q(i+im2,2)) - q(i,1)
1175 1499904 : qmin = q(i,1) - min(q(i,2),q(i,1),-q(i+im2,2))
1176 1510320 : dm(i,1) = sign(min(abs(tmp),qmax,qmin),tmp)
1177 : enddo
1178 :
1179 1510320 : do i=im2+1,im
1180 1510320 : dm(i, 1) = dm(i-im2, 1)
1181 : enddo
1182 : endif
1183 :
1184 333312 : if ( jlast+ng >= jm ) then
1185 : ! North
1186 1510320 : do i=1,im2
1187 1499904 : tmp = -D0_25*(q(i+im2,jm1)+q(i,jm1))
1188 1499904 : qmax = max(-q(i+im2,jm1),q(i,jm), q(i,jm1)) - q(i,jm)
1189 1499904 : qmin = q(i,jm) - min(-q(i+im2,jm1),q(i,jm), q(i,jm1))
1190 1510320 : dm(i,jm) = sign(min(abs(tmp),qmax,qmin),tmp)
1191 : enddo
1192 :
1193 1510320 : do i=im2+1,im
1194 1510320 : dm(i,jm) = dm(i-im2,jm)
1195 : enddo
1196 : endif
1197 :
1198 : endif
1199 :
1200 43223040 : if( jord > 0 ) then
1201 : !
1202 : ! Applies monotonic slope constraint (off if jord less than zero)
1203 : !
1204 336971712 : do j=js2gng1,jn2gng1
1205 85109997504 : do i=1,im
1206 84773025792 : qmax = max(q(i,j-1),q(i,j),q(i,j+1)) - q(i,j)
1207 84773025792 : qmin = q(i,j) - min(q(i,j-1),q(i,j),q(i,j+1))
1208 85067376576 : dm(i,j) = sign(min(abs(dm(i,j)),qmin,qmax),dm(i,j))
1209 : enddo
1210 : enddo
1211 : endif
1212 43223040 : return
1213 : !EOC
1214 : end subroutine ymist
1215 : !-----------------------------------------------------------------------
1216 :
1217 : !-----------------------------------------------------------------------
1218 : !BOP
1219 : ! !IROUTINE: fyppm
1220 : !
1221 : ! !INTERFACE:
1222 42491904 : subroutine fyppm(c, q, dm, flux, im, jm, ng, jord, iv, jfirst, jlast)
1223 : !-----------------------------------------------------------------------
1224 :
1225 : ! !USES:
1226 : implicit none
1227 :
1228 : ! !INPUT PARAMETERS:
1229 : integer im, jm ! Dimensions
1230 : integer jfirst, jlast ! Latitude strip
1231 : integer ng ! Max. NS dependencies
1232 : integer jord ! Approximation order
1233 : integer iv ! Scalar=0, Vector=1
1234 : real (r8) q(im,jfirst-ng:jlast+ng) ! mean value needed only N*2 S*2
1235 : real (r8) dm(im,jfirst-ng:jlast+ng) ! Slope needed only N*2 S*2
1236 : real (r8) c(im,jfirst:jlast+1) ! Courant N (like FLUX)
1237 :
1238 : ! !INPUT/OUTPUT PARAMETERS:
1239 84983808 : real (r8) ar(im,jfirst-1:jlast+1) ! AR needs to be ghosted on NS
1240 84983808 : real (r8) al(im,jfirst-1:jlast+2) ! AL needs to be ghosted on N2S
1241 84983808 : real (r8) a6(im,jfirst-1:jlast+1) ! A6 needs to be ghosted on NS
1242 :
1243 : ! !OUTPUT PARAMETERS:
1244 : real (r8) flux(im,jfirst:jlast+1) ! Flux N (see tp2c)
1245 :
1246 : ! !DESCRIPTION:
1247 : !
1248 : ! NG is passed from YTP for convenience -- it is actually 1 more in NS
1249 : ! than the actual number of latitudes needed here. But in the shared-memory
1250 : ! case it becomes 0, which is much cleaner.
1251 : !
1252 : ! !CALLED FROM:
1253 : ! ytp
1254 : !
1255 : ! !REVISION HISTORY:
1256 : !
1257 : ! SJL 99.04.13: Delivery
1258 : ! WS 99.04.19: Added jfirst:jlast concept; FYPPM only called from YTP
1259 : ! WS 99.04.21: Removed j1, j2 (j1=2, j2=jm-1 consistently)
1260 : ! removed a6,ar,al from argument list
1261 : ! WS 99.09.09: Documentation; indentation; cleaning
1262 : ! WS 99.10.22: Added ng as argument; Pruned arrays
1263 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
1264 : !
1265 : !EOP
1266 : !---------------------------------------------------------------------
1267 : !BOC
1268 : real (r8) r3, r23
1269 : parameter ( r3 = D1_0/D3_0, r23 = D2_0/D3_0 )
1270 : integer i, j, imh, jm1, lmt
1271 : integer js1g1, js2g0, js2g1, jn1g2, jn1g1, jn2g1
1272 : integer jan, jlow, jhigh, ilow, ihigh
1273 : integer ja(jlast-jfirst+3)
1274 : ! logical steep
1275 :
1276 : ! if(jord .eq. 6) then
1277 : ! steep = .true.
1278 : ! else
1279 : ! steep = .false.
1280 : ! endif
1281 :
1282 42491904 : imh = im / 2
1283 42491904 : jm1 = jm - 1
1284 :
1285 42491904 : js1g1 = max(1,jfirst-1) ! Ghost S*1
1286 42491904 : js2g0 = max(2,jfirst) ! No ghosting
1287 42491904 : js2g1 = max(2,jfirst-1) ! Ghost S*1
1288 42491904 : jn1g1 = min(jm,jlast+1) ! Ghost N*1
1289 42491904 : jn1g2 = min(jm,jlast+2) ! Ghost N*2
1290 42491904 : jn2g1 = min(jm-1,jlast+1) ! Ghost N*1
1291 :
1292 294787584 : do j=js2g1,jn1g2 ! AL needed N2S
1293 72955943424 : do i=1,im ! P, dm ghosted N2S2 (at least)
1294 72913451520 : al(i,j) = D0_5*(q(i,j-1)+q(i,j)) + r3*(dm(i,j-1) - dm(i,j))
1295 : enddo
1296 : enddo
1297 :
1298 : ! Yeh's steepening procedure; to be implemented
1299 : ! if(steep) call steepy(im, jm, jfirst, jlast, &
1300 : ! ng, q, al, dm )
1301 :
1302 252959616 : do j=js1g1,jn2g1 ! AR needed NS
1303 60867660672 : do i=1,im
1304 60825168768 : ar(i,j) = al(i,j+1) ! AL ghosted N2S
1305 : enddo
1306 : enddo
1307 :
1308 : ! WS 990726 : Added condition to decide if poles are on this processor
1309 :
1310 : ! Poles:
1311 :
1312 42491904 : if( iv == 0 ) then
1313 :
1314 42190848 : if ( jfirst == 1 ) then
1315 95588640 : do i=1,imh
1316 94929408 : al(i, 1) = al(i+imh,2)
1317 95588640 : al(i+imh,1) = al(i, 2)
1318 : enddo
1319 : endif
1320 :
1321 42190848 : if ( jlast == jm ) then
1322 95588640 : do i=1,imh
1323 94929408 : ar(i, jm) = ar(i+imh,jm1)
1324 95588640 : ar(i+imh,jm) = ar(i, jm1)
1325 : enddo
1326 : endif
1327 :
1328 : else
1329 :
1330 301056 : if ( jfirst == 1 ) then
1331 682080 : do i=1,imh
1332 677376 : al(i, 1) = -al(i+imh,2)
1333 682080 : al(i+imh,1) = -al(i, 2)
1334 : enddo
1335 : endif
1336 :
1337 301056 : if ( jlast == jm ) then
1338 682080 : do i=1,imh
1339 677376 : ar(i, jm) = -ar(i+imh,jm1)
1340 682080 : ar(i+imh,jm) = -ar(i, jm1)
1341 : enddo
1342 : endif
1343 :
1344 : endif
1345 :
1346 42491904 : if( jord == 3 .or. jord == 5 ) then
1347 0 : do j=js1g1,jn1g1 ! A6 needed NS
1348 0 : do i=1,im
1349 0 : a6(i,j) = D3_0*(q(i,j)+q(i,j) - (al(i,j)+ar(i,j)))
1350 : enddo
1351 : enddo
1352 : endif
1353 :
1354 42491904 : lmt = jord - 3
1355 :
1356 : ! do j=js1g1,jn1g1 ! A6, AR, AL needed NS
1357 : ! call lmppm(dm(1,j),a6(1,j),ar(1,j),al(1,j),q(1,j),im,lmt)
1358 : ! enddo
1359 :
1360 : #ifdef VECTORIZE
1361 : jan = 1
1362 : ja(1) = 1
1363 : ilow = 1
1364 : ihigh = im*(jn1g1-js1g1+1)
1365 : jlow = 1
1366 : jhigh = 1
1367 : call lmppmv(dm(1,js1g1), a6(1,js1g1), ar(1,js1g1), &
1368 : al(1,js1g1), q(1,js1g1), im*(jn1g1-js1g1+1), lmt, &
1369 : jan, ja, ilow, ihigh, jlow, jhigh, jlow, jhigh)
1370 : #else
1371 42491904 : call lmppm(dm(1,js1g1), a6(1,js1g1), ar(1,js1g1), &
1372 84983808 : al(1,js1g1), q(1,js1g1), im*(jn1g1-js1g1+1), lmt)
1373 : #endif
1374 :
1375 211131648 : do j=js2g0,jn1g1 ! flux needed N
1376 48779377920 : do i=1,im
1377 48736886016 : if(c(i,j).gt.D0_0) then
1378 23504592994 : flux(i,j) = ar(i,j-1) + D0_5*c(i,j)*(al(i,j-1) - ar(i,j-1) + &
1379 47009185988 : a6(i,j-1)*(D1_0-r23*c(i,j)) )
1380 : else
1381 25063653278 : flux(i,j) = al(i,j) - D0_5*c(i,j)*(ar(i,j) - al(i,j) + &
1382 25063653278 : a6(i,j)*(D1_0+r23*c(i,j)))
1383 : endif
1384 : enddo
1385 : enddo
1386 42491904 : return
1387 : !EOC
1388 : end subroutine fyppm
1389 : !-----------------------------------------------------------------------
1390 :
1391 : !-----------------------------------------------------------------------
1392 : !BOP
1393 : ! !IROUTINE: tpcc
1394 : !
1395 : ! !INTERFACE:
1396 344064 : subroutine tpcc(va, ymass, q, crx, cry, im, jm, ng_c, ng_d, &
1397 344064 : iord, jord, fx, fy, ffsl, cose, jfirst, jlast, &
1398 344064 : dm, qtmp, al, ar, a6 )
1399 : !-----------------------------------------------------------------------
1400 :
1401 : ! !USES:
1402 : implicit none
1403 :
1404 : ! !INPUT PARAMETERS:
1405 : integer im, jm ! Dimensions
1406 : integer ng_c !
1407 : integer ng_d !
1408 : integer jfirst, jlast ! Latitude strip
1409 : integer iord, jord ! Interpolation order in x,y
1410 : logical ffsl(jm) ! Flux-form semi-Lagrangian transport?
1411 : real (r8) cose(jm) ! Critical cosine (replicated)
1412 : real (r8) va(im,jfirst:jlast) ! Courant (unghosted like FX)
1413 : real (r8) q(im,jfirst-ng_d:jlast+ng_d) !
1414 : real (r8) crx(im,jfirst-ng_c:jlast+ng_c)
1415 : real (r8) cry(im,jfirst:jlast) ! Courant # (ghosted like FY)
1416 : real (r8) ymass(im,jfirst:jlast) ! Background y-mass-flux (ghosted like FY)
1417 :
1418 : ! Input 1D work arrays:
1419 : real (r8) dm(-im/3:im+im/3)
1420 : real (r8) qtmp(-im/3:im+im/3)
1421 : real (r8) al(-im/3:im+im/3)
1422 : real (r8) ar(-im/3:im+im/3)
1423 : real (r8) a6(-im/3:im+im/3)
1424 :
1425 : ! !OUTPUT PARAMETERS:
1426 : real (r8) fx(im,jfirst:jlast) ! Flux in x (unghosted)
1427 : real (r8) fy(im,jfirst:jlast) ! Flux in y (unghosted since iv==0)
1428 :
1429 : ! !DESCRIPTION:
1430 : ! In this routine the number
1431 : ! of north ghosted latitude min(abs(jord),2), and south ghosted
1432 : ! latitudes is XXXX
1433 : !
1434 : ! !CALLED FROM:
1435 : ! cd_core
1436 : !
1437 : ! !REVISION HISTORY:
1438 : ! SJL 99.04.13: Delivery
1439 : ! WS 99.04.13: Added jfirst:jlast concept
1440 : ! WS 99.05.10: Replaced JNP with JM, JMR with JM-1, IMR with IM
1441 : ! WS 99.05.10: Removed fvcore.h and JNP, IMH, IML definitions
1442 : ! WS 99.10.20: Pruned arrays
1443 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
1444 : !
1445 : !EOP
1446 : !-----------------------------------------------------------------------
1447 : !BOC
1448 : !
1449 : ! !LOCAL VARIABLES:
1450 :
1451 688128 : real (r8) adx(im,jfirst-1:jlast+2)
1452 : integer north, south
1453 : integer i, j, jp, im2, js2g0, js2gs, jn2g0, jn1g0, jn1gn
1454 688128 : real (r8) wk1v(im,jfirst-1:jlast+2)
1455 688128 : real (r8) fx1(im)
1456 688128 : real (r8) qtmpv(-im/3:im+im/3,jfirst-1:jlast+2)
1457 :
1458 344064 : im2 = im/2
1459 344064 : north = min(2,abs(jord)) ! north == 1 or 2
1460 344064 : south = north-1 ! south == 0 or 1
1461 344064 : js2g0 = max(2,jfirst)
1462 344064 : js2gs = max(2,jfirst-south)
1463 344064 : jn2g0 = min(jm-1,jlast)
1464 344064 : jn1gn = min(jm,jlast+north)
1465 344064 : jn1g0 = min(jm,jlast)
1466 :
1467 : ! This loop must be ghosted N*NG, S*NG
1468 :
1469 : call xtpv( im, ffsl, wk1v, q, crx, 1, crx, &
1470 : cose, 0, dm, qtmpv, al, ar, a6, &
1471 : jfirst, jlast, js2gs, jn1gn, jm, &
1472 : 1, jm, jfirst-1, jlast+2, &
1473 : jfirst-ng_d, jlast+ng_d, jfirst-ng_c, jlast+ng_c, &
1474 344064 : jfirst-ng_c, jlast+ng_c, jfirst-1, jlast+2)
1475 :
1476 2302272 : do j=js2gs,jn1gn
1477 :
1478 563963904 : do i=1,im-1
1479 1124011392 : adx(i,j) = q(i,j) + D0_5 * &
1480 1687975296 : (wk1v(i,j)-wk1v(i+1,j) + q(i,j)*(crx(i+1,j)-crx(i,j)))
1481 : enddo
1482 :
1483 3916416 : adx(im,j) = q(im,j) + D0_5 * &
1484 6218688 : (wk1v(im,j)-wk1v(1,j) + q(im,j)*(crx(1,j)-crx(im,j)))
1485 : enddo
1486 :
1487 344064 : call ycc(im, jm, fy, adx, cry, ymass, jord, 0,jfirst,jlast)
1488 :
1489 : ! For Scalar only!!!
1490 344064 : if ( jfirst == 1 ) then ! ( jfirst -ng_d <= 1 ) fails when
1491 : ! ng_d=3, ng_c=2, jlast-jfirst+1 = 3
1492 779520 : do i=1,im2
1493 779520 : q(i,1) = q(i+im2, 2)
1494 : enddo
1495 779520 : do i=im2+1,im
1496 779520 : q(i,1) = q(i-im2, 2)
1497 : enddo
1498 : endif
1499 :
1500 344064 : if ( jlast == jm ) then
1501 779520 : do i=1,im2
1502 779520 : fx1(i) = q(i+im2,jm)
1503 : enddo
1504 779520 : do i=im2+1,im
1505 779520 : fx1(i) = q(i-im2,jm)
1506 : enddo
1507 :
1508 1553664 : do i=1,im
1509 1553664 : if(va(i,jm) .gt. D0_0) then
1510 770789 : adx(i,jm) = q(i,jm) + D0_5*va(i,jm)*(q(i,jm-1)-q(i,jm))
1511 : else
1512 777499 : adx(i,jm) = q(i,jm) + D0_5*va(i,jm)*(q(i,jm)-fx1(i))
1513 : endif
1514 : enddo
1515 : endif
1516 :
1517 1365504 : do j=js2g0,jn2g0
1518 295540224 : do i=1,im
1519 294174720 : jp = j-va(i,j)
1520 : ! jp = j if va < 0
1521 : ! jp = j -1 if va < 0
1522 : ! q needed max(1, jfirst-1)
1523 295196160 : adx(i,j) = q(i,j) + D0_5*va(i,j)*(q(i,jp)-q(i,jp+1))
1524 : enddo
1525 : enddo
1526 :
1527 : call xtpv( im, ffsl, fx, adx, crx, iord, crx, &
1528 : cose, 0, dm, qtmpv, al, ar, a6, &
1529 : jfirst, jlast, js2g0, jn1g0, jm, &
1530 : 1, jm, jfirst, jlast, &
1531 : jfirst-1, jlast+2,jfirst-ng_c, jlast+ng_c, &
1532 344064 : jfirst-ng_c, jlast+ng_c, jfirst-1, jlast+2)
1533 :
1534 344064 : return
1535 : !EOC
1536 : end subroutine tpcc
1537 : !-----------------------------------------------------------------------
1538 :
1539 : !-----------------------------------------------------------------------
1540 : !BOP
1541 : ! !IROUTINE: ycc
1542 : !
1543 : ! !INTERFACE:
1544 688128 : subroutine ycc(im, jm, fy, q, vc, ymass, jord, iv, jfirst, jlast)
1545 : !-----------------------------------------------------------------------
1546 :
1547 : ! !USES:
1548 : implicit none
1549 :
1550 : ! !INPUT PARAMETERS:
1551 : integer im, jm ! Dimensions
1552 : integer jfirst, jlast ! Latitude strip
1553 : integer jord ! Approximation order
1554 : integer iv ! Scalar=0, Vector=1
1555 : real (r8) q(im,jfirst-1-iv:jlast+2) ! Field (N*2 S*(iv+1))
1556 : real (r8) vc(im,jfirst-iv:jlast) ! Courant (like FY)
1557 : real (r8) ymass(im,jfirst-iv:jlast) ! background mass flux
1558 :
1559 : ! !OUTPUT PARAMETERS:
1560 : real (r8) fy(im,jfirst-iv:jlast) ! Flux (S if iv=1)
1561 :
1562 : ! !DESCRIPTION:
1563 : ! Will Sawyer's note: In this routine the number
1564 : ! of ghosted latitudes NG is min(abs(jord),2). The scalar/vector
1565 : ! flag determines whether the flux FY needs to be ghosted on the
1566 : ! south. If called from CD\_CORE (iv==1) then it does, if called
1567 : ! from TPCC (iv==0) it does not.
1568 : !
1569 : ! !CALLED FROM:
1570 : ! cd_core
1571 : ! tpcc
1572 : !
1573 : ! !REVISION HISTORY:
1574 : !
1575 : ! SJL 99.04.13: Delivery
1576 : ! WS 99.04.19: Added jfirst:jlast concept
1577 : ! WS 99.04.27: DC removed as argument (local to this routine); DC on N
1578 : ! WS 99.05.10: Replaced JNP with JM, JMR with JM-1, IMR with IM
1579 : ! WS 99.05.10: Removed fvcore.h
1580 : ! WS 99.07.27: Built in tests for SP or NP
1581 : ! WS 99.09.09: Documentation; indentation; cleaning; pole treatment
1582 : ! WS 99.09.14: Loop limits
1583 : ! WS 00.05.14: Renamed ghost indices as per Kevin's definitions
1584 : !
1585 : !EOP
1586 : !---------------------------------------------------------------------
1587 : !BOC
1588 :
1589 : ! !LOCAL VARIABLES:
1590 1376256 : real (r8) dc(im,jfirst-iv:jlast+1)
1591 : real (r8) qmax, qmin
1592 : integer i, j, jt, im2, js2giv, js3giv, jn2g1, jn2g0
1593 :
1594 :
1595 688128 : im2 = im/2
1596 :
1597 688128 : js2giv = max(2,jfirst-iv)
1598 688128 : js3giv = max(3,jfirst-iv)
1599 688128 : jn2g1 = min(jm-1,jlast+1)
1600 688128 : jn2g0 = min(jm-1,jlast)
1601 :
1602 688128 : if(jord == 1) then
1603 383712 : do j=js2giv,jn2g0 ! FY needed on S*iv
1604 86120160 : do i=1,im
1605 : ! jt=j if vc > 0; jt=j+1 if vc <=0
1606 85736448 : jt = real(j+1,r8) - vc(i,j) ! VC ghosted like fy
1607 86034144 : fy(i,j) = q(i,jt)*ymass(i,j) ! ymass ghosted like fy
1608 : enddo ! q ghosted N*1, S*iv
1609 : enddo
1610 :
1611 : else
1612 :
1613 3269280 : do j=js3giv,jn2g1 ! dc needed N*1, S*iv
1614 771413664 : do i=1,im
1615 770811552 : dc(i,j) = D0_25*(q(i,j+1)-q(i,j-1)) ! q ghosted N*2, S*(iv+1)
1616 : enddo
1617 : enddo
1618 :
1619 602112 : if(iv.eq.0) then
1620 : ! Scalar.
1621 :
1622 : ! WS 99.07.27 : Split loops in SP and NP regions, added SP/NP tests
1623 :
1624 301056 : if ( jfirst-iv <= 2 ) then
1625 682080 : do i=1,im2
1626 682080 : dc(i, 2) = D0_25 * ( q(i,3) - q(i+im2,2) )
1627 : enddo
1628 :
1629 682080 : do i=im2+1,im
1630 682080 : dc(i, 2) = D0_25 * ( q(i,3) - q(i-im2,2) )
1631 : enddo
1632 : endif
1633 :
1634 301056 : if ( jlast == jm ) then
1635 682080 : do i=1,im2
1636 682080 : dc(i,jm) = D0_25 * ( q(i+im2,jm) - q(i,jm-1) )
1637 : enddo
1638 :
1639 682080 : do i=im2+1,im
1640 682080 : dc(i,jm) = D0_25 * ( q(i-im2,jm) - q(i,jm-1) )
1641 : enddo
1642 : endif
1643 :
1644 : else
1645 : ! Vector winds
1646 :
1647 : ! WS 99.07.27 : Split loops in SP and NP regions, added SP/NP tests
1648 :
1649 301056 : if ( jfirst-iv <= 2 ) then
1650 682080 : do i=1,im2
1651 682080 : dc(i, 2) = D0_25 * ( q(i,3) + q(i+im2,2) )
1652 : enddo
1653 :
1654 682080 : do i=im2+1,im
1655 682080 : dc(i, 2) = D0_25 * ( q(i,3) + q(i-im2,2) )
1656 : enddo
1657 : endif
1658 :
1659 301056 : if ( jlast == jm ) then
1660 682080 : do i=1,im2
1661 682080 : dc(i,jm) = -D0_25 * ( q(i,jm-1) + q(i+im2,jm) )
1662 : enddo
1663 :
1664 682080 : do i=im2+1,im
1665 682080 : dc(i,jm) = -D0_25 * ( q(i,jm-1) + q(i-im2,jm) )
1666 : enddo
1667 : endif
1668 :
1669 : endif
1670 :
1671 602112 : if( jord > 0 ) then
1672 : ! Monotonic constraint
1673 0 : do j=js3giv,jn2g1 ! DC needed N*1, S*iv
1674 0 : do i=1,im ! P ghosted N*2, S*(iv+1)
1675 0 : qmax = max(q(i,j-1),q(i,j),q(i,j+1)) - q(i,j)
1676 0 : qmin = q(i,j) - min(q(i,j-1),q(i,j),q(i,j+1))
1677 0 : dc(i,j) = sign(min(abs(dc(i,j)),qmin,qmax),dc(i,j))
1678 : enddo
1679 : enddo
1680 : !
1681 : ! WS 99.08.03 : Following loop split into SP and NP part
1682 : !
1683 0 : if ( jfirst-iv <= 2 ) then
1684 0 : do i=1,im
1685 0 : dc(i, 2) = D0_0
1686 : enddo
1687 : endif
1688 0 : if ( jlast == jm ) then
1689 0 : do i=1,im
1690 0 : dc(i,jm) = D0_0
1691 : enddo
1692 : endif
1693 : endif
1694 :
1695 2685984 : do j=js2giv,jn2g0 ! fy needed S*iv
1696 602841120 : do i=1,im
1697 600155136 : jt = real(j+1,r8) - vc(i,j) ! vc, ymass ghosted like fy
1698 602239008 : fy(i,j) = (q(i,jt)+(sign(D1_0,vc(i,j))-vc(i,j))*dc(i,jt))*ymass(i,j)
1699 : enddo
1700 : enddo
1701 : endif
1702 688128 : return
1703 : !EOC
1704 : end subroutine ycc
1705 : !-----------------------------------------------------------------------
1706 :
1707 : #ifdef VECTORIZE
1708 : !-----------------------------------------------------------------------
1709 : !BOP
1710 : ! !IROUTINE: xtpv
1711 : !
1712 : ! !INTERFACE:
1713 : subroutine xtpv(im, ffslv, fxv, qv, cv, iord, mfxv, &
1714 : cosav, id, dm, qtmpv, al, ar, a6, &
1715 : jfirst, jlast, jlow, jhigh, jm, &
1716 : jl2, jh2, jl3, jh3, &
1717 : jl4, jh4, jl5, jh5, &
1718 : jl7, jh7, jl11, jh11)
1719 : !-----------------------------------------------------------------------
1720 :
1721 : implicit none
1722 :
1723 : ! !INPUT PARAMETERS:
1724 : integer id ! ID = 0: density (mfx = C)
1725 : ! ID = 1: mixing ratio (mfx is mass flux)
1726 :
1727 : integer im ! Total longitudes
1728 : real (r8) cv(im,jl5:jh5) ! Courant numbers
1729 : real (r8) qv(im,jl4:jh4)
1730 : real (r8) mfxv(im,jl7:jh7)
1731 : logical ffslv(jl2:jh2)
1732 : integer iord
1733 : integer jfirst, jlast, jlow, jhigh, jm
1734 : integer jl2, jh2, jl3, jh3, jl4, jh4, jl5, jh5
1735 : integer jl7, jh7, jl11, jh11
1736 : real (r8) cosav(jm)
1737 :
1738 : ! !INPUT/OUTPUT PARAMETERS:
1739 : real (r8) qtmpv(-im/3:im+im/3,jl11:jh11) ! Input work arrays:
1740 : real (r8) dm(-im/3:im+im/3)
1741 : real (r8) al(-im/3:im+im/3)
1742 : real (r8) ar(-im/3:im+im/3)
1743 : real (r8) a6(-im/3:im+im/3)
1744 :
1745 : ! !OUTPUT PARAMETERS:
1746 : real (r8) fxv(im,jl3:jh3)
1747 :
1748 : ! !DESCRIPTION:
1749 : !
1750 : !
1751 : ! !REVISION HISTORY:
1752 : ! 99.01.01 Lin Creation
1753 : ! 01.03.27 Sawyer Additional ProTeX documentation
1754 : !
1755 : !EOP
1756 : !-----------------------------------------------------------------------
1757 : !BOC
1758 :
1759 : ! Local:
1760 : real (r8) cos_upw !critical cosine for upwind
1761 : real (r8) cos_van !critical cosine for van Leer
1762 : real (r8) cos_ppm !critical cosine for ppm
1763 :
1764 : parameter (cos_upw = D0_05) !roughly at 87 deg.
1765 : parameter (cos_van = D0_25) !roughly at 75 deg.
1766 : parameter (cos_ppm = D0_25)
1767 :
1768 : real (r8) r24
1769 : parameter (r24 = D1_0/D24_0)
1770 :
1771 : integer i, imp, j
1772 : real (r8) qmax, qmin
1773 : real (r8) rut, tmp
1774 : real (r8) dmv(-im/3:im+im/3,jlow:jhigh)
1775 : integer iu, itmp, ist
1776 : integer isave(im,jlow:jhigh)
1777 : integer iuwv(jlow:jhigh), iuev(jlow:jhigh)
1778 :
1779 : integer jatn, jafn, ja
1780 : integer jat(jhigh-jlow+1), jaf(jhigh-jlow+1)
1781 : integer jattn, jatfn, jaftn, jaffn
1782 : integer jatt(jhigh-jlow+1), jatf(jhigh-jlow+1)
1783 : integer jaft(jhigh-jlow+1), jaff(jhigh-jlow+1)
1784 : integer jatftn, jatffn
1785 : integer jatft(jhigh-jlow+1), jatff(jhigh-jlow+1)
1786 : integer jafftn1, jafffn1
1787 : integer jafft1(jhigh-jlow+1), jafff1(jhigh-jlow+1)
1788 : integer jafftn2, jafffn2
1789 : integer jafft2(jhigh-jlow+1), jafff2(jhigh-jlow+1)
1790 : real (r8) qsum((-im/3)-1:im+im/3,jlow:jhigh) ! work arrays
1791 :
1792 :
1793 : jatn = 0
1794 : jafn = 0
1795 : jattn = 0
1796 : jatfn = 0
1797 : jaftn = 0
1798 : jaffn = 0
1799 : jatftn = 0
1800 : jatffn = 0
1801 : jafftn1 = 0
1802 : jafffn1 = 0
1803 : jafftn2 = 0
1804 : jafffn2 = 0
1805 : !call ftrace_region_begin("xtpv_1")
1806 : do j = jlow, jhigh
1807 : if (ffslv(j)) then
1808 : jatn = jatn + 1
1809 : jat(jatn) = j
1810 : if( iord == 1 .or. cosav(j) < cos_upw) then
1811 : jattn = jattn + 1
1812 : jatt(jattn) = j
1813 : else
1814 : jatfn = jatfn + 1
1815 : jatf(jatfn) = j
1816 : if(iord .ge. 3 .and. cosav(j) .gt. cos_ppm) then
1817 : jatftn = jatftn + 1
1818 : jatft(jatftn) = j
1819 : else
1820 : jatffn = jatffn + 1
1821 : jatff(jatffn) = j
1822 : endif
1823 : endif
1824 : else
1825 : jafn = jafn + 1
1826 : jaf(jafn) = j
1827 : if( iord == 1 .or. cosav(j) < cos_upw) then
1828 : jaftn = jaftn + 1
1829 : jaft(jaftn) = j
1830 : else
1831 : jaffn = jaffn + 1
1832 : jaff(jaffn) = j
1833 : if(iord > 0 .or. cosav(j) < cos_van) then
1834 : jafftn1 = jafftn1 + 1
1835 : jafft1(jafftn1) = j
1836 : else
1837 : jafffn1 = jafffn1 + 1
1838 : jafff1(jafffn1) = j
1839 : endif
1840 : if( abs(iord).eq.2 .or. cosav(j) .lt. cos_van ) then
1841 : jafftn2 = jafftn2 + 1
1842 : jafft2(jafftn2) = j
1843 : else
1844 : jafffn2 = jafffn2 + 1
1845 : jafff2(jafffn2) = j
1846 : endif
1847 : endif
1848 : endif
1849 : enddo
1850 : !call ftrace_region_end("xtpv_1")
1851 :
1852 : imp = im + 1
1853 :
1854 : do j = jlow, jhigh
1855 : do i=1,im
1856 : qtmpv(i,j) = qv(i,j)
1857 : enddo
1858 : enddo
1859 :
1860 : ! Flux-Form Semi-Lagrangian transport
1861 :
1862 : !call ftrace_region_begin("xtpv_2")
1863 : do ja = 1, jatn
1864 : j = jat(ja)
1865 :
1866 : ! Figure out ghost zone for the western edge:
1867 : iuwv(j) = -cv(1,j)
1868 : iuwv(j) = min(0, iuwv(j))
1869 :
1870 : do i=iuwv(j), 0
1871 : qtmpv(i,j) = qv(im+i,j)
1872 : enddo
1873 :
1874 : ! Figure out ghost zone for the eastern edge:
1875 : iuev(j) = im - cv(im,j)
1876 : iuev(j) = max(imp, iuev(j))
1877 :
1878 : do i=imp, iuev(j)
1879 : qtmpv(i,j) = qv(i-im,j)
1880 : enddo
1881 :
1882 : enddo
1883 : !call ftrace_region_end("xtpv_2")
1884 :
1885 : !call ftrace_region_begin("xtpv_3")
1886 : do ja = 1, jattn
1887 : j = jatt(ja)
1888 :
1889 : do i=1,im
1890 : iu = cv(i,j)
1891 : if(cv(i,j) .le. D0_0) then
1892 : itmp = i - iu
1893 : isave(i,j) = itmp - 1
1894 : else
1895 : itmp = i - iu - 1
1896 : isave(i,j) = itmp + 1
1897 : endif
1898 : fxv(i,j) = (cv(i,j)-iu) * qtmpv(itmp,j)
1899 : enddo
1900 :
1901 : enddo
1902 : !call ftrace_region_end("xtpv_3")
1903 :
1904 : !call ftrace_region_begin("xtpv_4")
1905 : do ja = 1, jatfn
1906 : j = jatf(ja)
1907 :
1908 : do i=1,im
1909 : ! 2nd order slope
1910 : tmp = D0_25*(qtmpv(i+1,j) - qtmpv(i-1,j))
1911 : qmax = max(qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j)) - qtmpv(i,j)
1912 : qmin = qtmpv(i,j) - min(qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j))
1913 : dmv(i,j) = sign(min(abs(tmp),qmax,qmin), tmp)
1914 : enddo
1915 :
1916 : do i=iuwv(j), 0
1917 : dmv(i,j) = dmv(im+i,j)
1918 : enddo
1919 :
1920 : do i=imp, iuev(j)
1921 : dmv(i,j) = dmv(i-im,j)
1922 : enddo
1923 :
1924 : enddo
1925 : !call ftrace_region_end("xtpv_4")
1926 :
1927 : call fxppmv(im, cv, mfxv, qtmpv, dmv, fxv, iord, &
1928 : iuwv, iuev, ffslv, isave, jatftn, jatft, jlow, jhigh, &
1929 : jl2, jh2, jl3, jh3, jl5, jh5, jl7, jh7, jl11, jh11)
1930 :
1931 : !call ftrace_region_begin("xtpv_5")
1932 : do ja = 1, jatffn
1933 : j = jatff(ja)
1934 :
1935 : do i=1,im
1936 : iu = cv(i,j)
1937 : rut = cv(i,j) - iu
1938 : if(cv(i,j) .le. D0_0) then
1939 : itmp = i - iu
1940 : isave(i,j) = itmp - 1
1941 : fxv(i,j) = rut*(qtmpv(itmp,j)-dmv(itmp,j)*(D1_0+rut))
1942 : else
1943 : itmp = i - iu - 1
1944 : isave(i,j) = itmp + 1
1945 : fxv(i,j) = rut*(qtmpv(itmp,j)+dmv(itmp,j)*(D1_0-rut))
1946 : endif
1947 : enddo
1948 :
1949 : enddo
1950 : !call ftrace_region_end("xtpv_5")
1951 :
1952 : !call ftrace_region_begin("xtpv_6")
1953 : do ja = 1, jatn
1954 : j = jat(ja)
1955 : qsum(iuwv(j)-1,j) = D0_0
1956 : do i = iuwv(j), iuev(j)
1957 : qsum(i,j) = qsum(i-1,j) + qtmpv(i,j)
1958 : end do
1959 :
1960 : !
1961 : ! The boolean terms:
1962 : ! a) .and. (isave(i,j) < i)
1963 : ! b) .and. (i <= isave(i,j))
1964 : ! are needed in the IF statements below because I cannot prove to myself
1965 : ! that the relationship between i and isave are such to guarantee that
1966 : ! there is always at least one term from qsum (qtmpv,j) contributed to fxv.
1967 : !
1968 :
1969 : do i=1,im
1970 : if(cv(i,j) >= D1_0 .and. (isave(i,j) < i) ) then
1971 : fxv(i,j) = fxv(i,j) + (qsum(i-1,j) - qsum(isave(i,j) - 1,j))
1972 : else if (cv(i,j) <= -D1_0 .and. (i <= isave(i,j)) ) then
1973 : fxv(i,j) = fxv(i,j) - (qsum(isave(i,j),j) - qsum(i-1,j))
1974 : end if
1975 : end do
1976 :
1977 : if(id .ne. 0) then
1978 : do i=1,im
1979 : fxv(i,j) = fxv(i,j)*mfxv(i,j)
1980 : enddo
1981 : endif
1982 :
1983 : enddo
1984 : !call ftrace_region_end("xtpv_6")
1985 :
1986 : ! Regular PPM (Eulerian without FFSL extension)
1987 :
1988 : !call ftrace_region_begin("xtpv_7")
1989 : do ja = 1, jafn
1990 : j = jaf(ja)
1991 :
1992 : qtmpv(imp,j) = qv(1,j)
1993 : qtmpv( 0,j) = qv(im,j)
1994 : enddo
1995 :
1996 : do ja = 1, jaftn
1997 : j = jaft(ja)
1998 :
1999 : do i=1,im
2000 : iu = real(i,r8) - cv(i,j)
2001 : fxv(i,j) = mfxv(i,j)*qtmpv(iu,j)
2002 : enddo
2003 : enddo
2004 :
2005 : do ja = 1, jaffn
2006 : j = jaff(ja)
2007 :
2008 : qtmpv(-1,j) = qv(im-1,j)
2009 : qtmpv(imp+1,j) = qv(2,j)
2010 :
2011 : enddo
2012 : !call ftrace_region_end("xtpv_7")
2013 :
2014 : !call ftrace_region_begin("xtpv_8")
2015 : do ja = 1, jafftn1
2016 : j = jafft1(ja)
2017 :
2018 : ! In-line xmist
2019 :
2020 : do i=1,im
2021 : dmv(i,j) = r24*(D8_0*(qtmpv(i+1,j) - qtmpv(i-1,j)) + qtmpv(i-2,j) - qtmpv(i+2,j))
2022 : enddo
2023 :
2024 : ! Apply monotonicity constraint (Lin et al. 1994, MWR)
2025 : do i=1,im
2026 : qmax = max( qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j) ) - qtmpv(i,j)
2027 : qmin = qtmpv(i,j) - min( qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j) )
2028 : dmv(i,j) = sign( min(abs(dmv(i,j)), qmax, qmin), dmv(i,j) )
2029 : enddo
2030 :
2031 : enddo
2032 : !call ftrace_region_end("xtpv_8")
2033 :
2034 : !call ftrace_region_begin("xtpv_9")
2035 : do ja = 1, jafffn1
2036 : j = jafff1(ja)
2037 :
2038 : ! In-line xmist
2039 :
2040 : if(iord .le. 2) then
2041 : do i=1,im
2042 : dmv(i,j) = r24*(D8_0*(qtmpv(i+1,j) - qtmpv(i-1,j)) + qtmpv(i-2,j) - qtmpv(i+2,j))
2043 : enddo
2044 : else
2045 : do i=1,im
2046 : dmv(i,j) = D0_25*(qtmpv(i+1,j) - qtmpv(i-1,j))
2047 : enddo
2048 : endif
2049 :
2050 : if( iord >= 0 ) then
2051 :
2052 : ! Apply monotonicity constraint (Lin et al. 1994, MWR)
2053 : do i=1,im
2054 : qmax = max( qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j) ) - qtmpv(i,j)
2055 : qmin = qtmpv(i,j) - min( qtmpv(i-1,j), qtmpv(i,j), qtmpv(i+1,j) )
2056 : dmv(i,j) = sign( min(abs(dmv(i,j)), qmax, qmin), dmv(i,j) )
2057 : enddo
2058 : endif
2059 :
2060 : enddo
2061 : !call ftrace_region_end("xtpv_9")
2062 :
2063 : !call ftrace_region_begin("xtpv_10")
2064 : do ja = 1, jaffn
2065 : j = jaff(ja)
2066 :
2067 : dmv(0,j) = dmv(im,j)
2068 :
2069 : enddo
2070 : !call ftrace_region_end("xtpv_10")
2071 :
2072 : !call ftrace_region_begin("xtpv_11")
2073 : do ja = 1, jafftn2
2074 : j = jafft2(ja)
2075 :
2076 : do i=1,im
2077 : iu = real(i,r8) - cv(i,j)
2078 : fxv(i,j) = mfxv(i,j)*(qtmpv(iu,j)+dmv(iu,j)*(sign(D1_0,cv(i,j))-cv(i,j)))
2079 : enddo
2080 :
2081 : enddo
2082 : !call ftrace_region_end("xtpv_11")
2083 :
2084 : call fxppmv(im, cv, mfxv, qtmpv, dmv, fxv, iord, &
2085 : iuwv, iuev, ffslv, isave, jafffn2, jafff2, jlow, jhigh, &
2086 : jl2, jh2, jl3, jh3, jl5, jh5, jl7, jh7, jl11, jh11)
2087 :
2088 : return
2089 : !EOC
2090 : end subroutine xtpv
2091 : !-----------------------------------------------------------------------
2092 :
2093 : !-----------------------------------------------------------------------
2094 : !BOP
2095 : ! !IROUTINE: fxppmv
2096 : !
2097 : ! !INTERFACE:
2098 :
2099 : subroutine fxppmv(im, c, mfx, p, dm, fx, iord, &
2100 : iuw, iue, ffsl, isave, jan, ja, jlow, jhigh, &
2101 : jl2, jh2, jl3, jh3, jl5, jh5, jl7, jh7, jl11, jh11)
2102 : !-----------------------------------------------------------------------
2103 : !
2104 : ! !USES:
2105 : implicit none
2106 :
2107 : ! !INPUT PARAMETERS:
2108 : integer jan, ja(jan), jlow, jhigh, jj, j
2109 : integer jl2, jh2, jl3, jh3, jl5, jh5, jl7, jh7, jl11, jh11
2110 : integer im, iord
2111 : real (r8) c(im,jl5:jh5)
2112 : real (r8) p(-im/3:im+im/3,jl11:jh11)
2113 : real (r8) dm(-im/3:im+im/3,jlow:jhigh)
2114 : real (r8) mfx(im,jl7:jh7)
2115 : integer iuw(jlow:jhigh), iue(jlow:jhigh)
2116 : logical ffsl(jl2:jh2)
2117 : integer isave(im,jlow:jhigh)
2118 :
2119 : ! !OUTPUT PARAMETERS:
2120 : real (r8) fx(im,jl3:jh3)
2121 :
2122 : ! !DESCRIPTION:
2123 : !
2124 : !
2125 : ! !REVISION HISTORY:
2126 : ! 99.01.01 Lin Creation
2127 : ! 01.03.27 Sawyer Additional ProTeX documentation
2128 : !
2129 : !EOP
2130 : !-----------------------------------------------------------------------
2131 : !BOC
2132 : !
2133 : ! !LOCAL VARIABLES:
2134 : real (r8) r3, r23
2135 : parameter ( r3 = D1_0/D3_0, r23 = D2_0/D3_0 )
2136 :
2137 : integer i, lmt
2138 : integer iu, itmp
2139 : real (r8) ru
2140 : logical steep
2141 : real (r8) al(-im/3:im+im/3,jlow:jhigh)
2142 : real (r8) ar(-im/3:im+im/3,jlow:jhigh)
2143 : real (r8) a6(-im/3:im+im/3,jlow:jhigh)
2144 :
2145 : integer jbtn, jbfn
2146 : integer jbt(jan), jbf(jan)
2147 : integer ilow, ihigh
2148 :
2149 : ilow = -im/3
2150 : ihigh = im + im/3
2151 :
2152 : if( iord == 6 ) then
2153 : steep = .true.
2154 : else
2155 : steep = .false.
2156 : endif
2157 :
2158 : do jj = 1, jan
2159 : j = ja(jj)
2160 :
2161 : do i=1,im
2162 : al(i,j) = D0_5*(p(i-1,j)+p(i,j)) + (dm(i-1,j) - dm(i,j))*r3
2163 : enddo
2164 :
2165 : enddo
2166 :
2167 : if (steep) then
2168 :
2169 : call steepxv( im, p, al, dm, jan, ja, jlow, jhigh, jl11, jh11 )
2170 :
2171 : endif
2172 :
2173 : do jj = 1, jan
2174 : j = ja(jj)
2175 :
2176 : do i=1,im-1
2177 : ar(i,j) = al(i+1,j)
2178 : enddo
2179 : ar(im,j) = al(1,j)
2180 :
2181 : enddo
2182 :
2183 : if(iord == 7) then
2184 :
2185 : call huynhv(im, ar, al, p, a6, dm, jan, ja, jlow, jhigh, jl11, jh11 )
2186 :
2187 : else
2188 :
2189 : if(iord .eq. 3 .or. iord .eq. 5) then
2190 :
2191 : do jj = 1, jan
2192 : j = ja(jj)
2193 :
2194 : do i=1,im
2195 : a6(i,j) = D3_0*(p(i,j)+p(i,j) - (al(i,j)+ar(i,j)))
2196 : enddo
2197 :
2198 : enddo
2199 : endif
2200 :
2201 : lmt = iord - 3
2202 :
2203 : call lmppmv( dm, a6, ar, al, p, im, lmt, jan, ja, ilow, ihigh, &
2204 : jlow, jhigh, jl11, jh11 )
2205 :
2206 : endif
2207 :
2208 : jbtn = 0
2209 : jbfn = 0
2210 : do jj = 1, jan
2211 : j = ja(jj)
2212 : if( ffsl(j) ) then
2213 : jbtn = jbtn + 1
2214 : jbt(jbtn) = j
2215 : else
2216 : jbfn = jbfn + 1
2217 : jbf(jbfn) = j
2218 : endif
2219 : enddo
2220 :
2221 : do jj = 1, jbtn
2222 : j = jbt(jj)
2223 :
2224 : do i=iuw(j), 0
2225 : al(i,j) = al(im+i,j)
2226 : ar(i,j) = ar(im+i,j)
2227 : a6(i,j) = a6(im+i,j)
2228 : enddo
2229 :
2230 : do i=im+1, iue(j)
2231 : al(i,j) = al(i-im,j)
2232 : ar(i,j) = ar(i-im,j)
2233 : a6(i,j) = a6(i-im,j)
2234 : enddo
2235 :
2236 : do i=1,im
2237 : iu = c(i,j)
2238 : ru = c(i,j) - iu
2239 : if(c(i,j) .gt. D0_0) then
2240 : itmp = i - iu - 1
2241 : isave(i,j) = itmp + 1
2242 : fx(i,j) = ru*(ar(itmp,j)+D0_5*ru*(al(itmp,j)-ar(itmp,j) + &
2243 : a6(itmp,j)*(D1_0-r23*ru)) )
2244 : else
2245 : itmp = i - iu
2246 : isave(i,j) = itmp - 1
2247 : fx(i,j) = ru*(al(itmp,j)-D0_5*ru*(ar(itmp,j)-al(itmp,j) + &
2248 : a6(itmp,j)*(D1_0+r23*ru)) )
2249 : endif
2250 : enddo
2251 :
2252 : enddo
2253 :
2254 : do jj = 1, jbfn
2255 : j = jbf(jj)
2256 :
2257 : al(0,j) = al(im,j)
2258 : ar(0,j) = ar(im,j)
2259 : a6(0,j) = a6(im,j)
2260 : do i=1,im
2261 : if(c(i,j) .gt. D0_0) then
2262 : fx(i,j) = ar(i-1,j) + D0_5*c(i,j)*(al(i-1,j) - ar(i-1,j) + &
2263 : a6(i-1,j)*(D1_0-r23*c(i,j)) )
2264 : else
2265 : fx(i,j) = al(i,j) - D0_5*c(i,j)*(ar(i,j) - al(i,j) + &
2266 : a6(i,j)*(D1_0+r23*c(i,j)))
2267 : endif
2268 : fx(i,j) = mfx(i,j) * fx(i,j)
2269 : enddo
2270 :
2271 : enddo
2272 :
2273 : return
2274 : !EOC
2275 : end subroutine fxppmv
2276 : !-----------------------------------------------------------------------
2277 :
2278 : !-----------------------------------------------------------------------
2279 : !BOP
2280 : ! !IROUTINE: steepxv
2281 : !
2282 : ! !INTERFACE:
2283 : subroutine steepxv(im, p, al, dm, jan, ja, jlow, jhigh, jl11, jh11 )
2284 : !-----------------------------------------------------------------------
2285 :
2286 : ! !USES:
2287 : implicit none
2288 :
2289 : ! !INPUT PARAMETERS:
2290 : integer im
2291 : integer jan, ja(jan), jlow, jhigh, jl11, jh11
2292 : real (r8) p(-im/3:im+im/3,jl11:jh11)
2293 : real (r8) dm(-im/3:im+im/3,jlow:jhigh)
2294 :
2295 : ! !INPUT/OUTPUT PARAMETERS:
2296 : real (r8) al(im,jlow:jhigh)
2297 :
2298 : ! !DESCRIPTION:
2299 : !
2300 : !
2301 : ! !REVISION HISTORY:
2302 : ! 99.01.01 Lin Creation
2303 : ! 01.03.27 Sawyer Additional ProTeX documentation
2304 : !
2305 : !EOP
2306 : !-----------------------------------------------------------------------
2307 : !BOC
2308 : !
2309 : ! !LOCAL VARIABLES:
2310 : integer i, jj, j
2311 : real (r8) r3
2312 : parameter ( r3 = D1_0/D3_0 )
2313 :
2314 : real (r8) dh(0:im,jlow:jhigh)
2315 : real (r8) d2(0:im+1,jlow:jhigh)
2316 : real (r8) eta(0:im,jlow:jhigh)
2317 : real (r8) xxx, bbb, ccc
2318 :
2319 : do jj = 1, jan
2320 : j = ja(jj)
2321 :
2322 : do i=0,im
2323 : dh(i,j) = p(i+1,j) - p(i,j)
2324 : enddo
2325 :
2326 : ! Needs dh(0:im,j)
2327 : do i=1,im
2328 : d2(i,j) = dh(i,j) - dh(i-1,j)
2329 : enddo
2330 : d2(0,j) = d2(im,j)
2331 : d2(im+1,j) = d2(1,j)
2332 :
2333 : ! needs p(-1:im+2,j), d2(0:im+1,j)
2334 : do i=1,im
2335 : if( d2(i+1,j)*d2(i-1,j).lt.D0_0 .and. p(i+1,j).ne.p(i-1,j) ) then
2336 : xxx = D1_0 - D0_5 * ( p(i+2,j) - p(i-2,j) ) / ( p(i+1,j) - p(i-1,j) )
2337 : eta(i,j) = max(D0_0, min(xxx, D0_5) )
2338 : else
2339 : eta(i,j) = D0_0
2340 : endif
2341 : enddo
2342 :
2343 : eta(0,j) = eta(im,j)
2344 :
2345 : ! needs eta(0:im,j), dh(0:im-1,j), dm(0:im,j)
2346 : do i=1,im
2347 : bbb = ( D2_0*eta(i,j ) - eta(i-1,j) ) * dm(i-1,j)
2348 : ccc = ( D2_0*eta(i-1,j) - eta(i,j ) ) * dm(i,j )
2349 : al(i,j) = al(i,j) + D0_5*( eta(i-1,j) - eta(i,j)) * dh(i-1,j) + (bbb - ccc) * r3
2350 : enddo
2351 :
2352 : enddo
2353 :
2354 : return
2355 : !EOC
2356 : end subroutine steepxv
2357 : !-----------------------------------------------------------------------
2358 :
2359 : !-----------------------------------------------------------------------
2360 : !BOP
2361 : ! !IROUTINE: huynhv --- Enforce Huynh's 2nd constraint in 1D periodic domain
2362 : !
2363 : ! !INTERFACE:
2364 : subroutine huynhv(im, ar, al, p, d2, d1, jan, ja, jlow, jhigh, jl11, jh11)
2365 : !-----------------------------------------------------------------------
2366 :
2367 : ! !USES:
2368 :
2369 : implicit none
2370 :
2371 : ! !INPUT PARAMETERS:
2372 : integer im
2373 : integer jan, ja(jan), jlow, jhigh, jl11, jh11
2374 : real(r8) p(im,jl11:jh11)
2375 :
2376 : ! !OUTPUT PARAMETERS:
2377 : real(r8) ar(im,jlow:jhigh)
2378 : real(r8) al(im,jlow:jhigh)
2379 : real(r8) d2(im,jlow:jhigh)
2380 : real(r8) d1(im,jlow:jhigh)
2381 :
2382 : ! !DESCRIPTION:
2383 : !
2384 : ! Enforce Huynh's 2nd constraint in 1D periodic domain
2385 : !
2386 : ! !REVISION HISTORY:
2387 : ! 99.01.01 Lin Creation
2388 : ! 01.03.27 Sawyer Additional ProTeX documentation
2389 : !
2390 : !EOP
2391 : !-----------------------------------------------------------------------
2392 : !BOC
2393 : !
2394 : ! !LOCAL VARIABLES:
2395 : integer i, jj, j
2396 : real(r8) pmp
2397 : real(r8) lac
2398 : real(r8) pmin
2399 : real(r8) pmax
2400 :
2401 : do jj = 1, jan
2402 : j = ja(jj)
2403 :
2404 : ! Compute d1 and d2
2405 : d1(1,j) = p(1,j) - p(im,j)
2406 : do i=2,im
2407 : d1(i,j) = p(i,j) - p(i-1,j)
2408 : enddo
2409 :
2410 : do i=1,im-1
2411 : d2(i,j) = d1(i+1,j) - d1(i,j)
2412 : enddo
2413 : d2(im,j) = d1(1,j) - d1(im,j)
2414 :
2415 : ! Constraint for AR
2416 : ! i = 1
2417 : pmp = p(1,j) + D2_0 * d1(1,j)
2418 : lac = p(1,j) + D0_5 * (d1(1,j)+d2(im,j)) + d2(im,j)
2419 : pmin = min(p(1,j), pmp, lac)
2420 : pmax = max(p(1,j), pmp, lac)
2421 : ar(1,j) = min(pmax, max(ar(1,j), pmin))
2422 :
2423 : do i=2, im
2424 : pmp = p(i,j) + D2_0*d1(i,j)
2425 : lac = p(i,j) + D0_5*(d1(i,j)+d2(i-1,j)) + d2(i-1,j)
2426 : pmin = min(p(i,j), pmp, lac)
2427 : pmax = max(p(i,j), pmp, lac)
2428 : ar(i,j) = min(pmax, max(ar(i,j), pmin))
2429 : enddo
2430 :
2431 : ! Constraint for AL
2432 : do i=1, im-1
2433 : pmp = p(i,j) - D2_0*d1(i+1,j)
2434 : lac = p(i,j) + D0_5*(d2(i+1,j)-d1(i+1,j)) + d2(i+1,j)
2435 : pmin = min(p(i,j), pmp, lac)
2436 : pmax = max(p(i,j), pmp, lac)
2437 : al(i,j) = min(pmax, max(al(i,j), pmin))
2438 : enddo
2439 :
2440 : ! i=im
2441 : i = im
2442 : pmp = p(im,j) - D2_0*d1(1,j)
2443 : lac = p(im,j) + D0_5*(d2(1,j)-d1(1,j)) + d2(1,j)
2444 : pmin = min(p(im,j), pmp, lac)
2445 : pmax = max(p(im,j), pmp, lac)
2446 : al(im,j) = min(pmax, max(al(im,j), pmin))
2447 :
2448 : ! compute A6 (d2)
2449 : do i=1, im
2450 : d2(i,j) = D3_0*(p(i,j)+p(i,j) - (al(i,j)+ar(i,j)))
2451 : enddo
2452 :
2453 : enddo
2454 :
2455 : return
2456 : !EOC
2457 : end subroutine huynhv
2458 : !-----------------------------------------------------------------------
2459 :
2460 : !-----------------------------------------------------------------------
2461 : !BOP
2462 : ! !IROUTINE: lmppmv
2463 : !
2464 : ! !INTERFACE:
2465 : subroutine lmppmv(dm, a6, ar, al, p, im, lmt, jan, ja, &
2466 : ilow, ihigh, jlow, jhigh, jl11, jh11)
2467 : !-----------------------------------------------------------------------
2468 :
2469 : implicit none
2470 :
2471 : ! !INPUT PARAMETERS:
2472 : integer im ! Total longitudes
2473 : integer jan, ja(jan), ilow, ihigh, jlow, jhigh, jl11, jh11
2474 : integer lmt ! LMT = 0: full monotonicity
2475 : ! LMT = 1: Improved and simplified full monotonic constraint
2476 : ! LMT = 2: positive-definite constraint
2477 : ! LMT = 3: Quasi-monotone constraint
2478 : real(r8) p(ilow:ihigh,jl11:jh11)
2479 : real(r8) dm(ilow:ihigh,jlow:jhigh)
2480 :
2481 : ! !OUTPUT PARAMETERS:
2482 : real(r8) a6(ilow:ihigh,jlow:jhigh)
2483 : real(r8) ar(ilow:ihigh,jlow:jhigh)
2484 : real(r8) al(ilow:ihigh,jlow:jhigh)
2485 :
2486 : ! !DESCRIPTION:
2487 : !
2488 : !
2489 : ! !REVISION HISTORY:
2490 : ! 99.01.01 Lin Creation
2491 : ! 01.03.27 Sawyer Additional ProTeX documentation
2492 : !
2493 : !EOP
2494 : !-----------------------------------------------------------------------
2495 : !BOC
2496 : !
2497 : ! !LOCAL VARIABLES:
2498 : real (r8) r12
2499 : parameter ( r12 = D1_0/D12_0 )
2500 :
2501 : real (r8) da1, da2, fmin, a6da
2502 : real (r8) dr, dl
2503 :
2504 : integer i, jj, j
2505 :
2506 : ! LMT = 0: full monotonicity
2507 : ! LMT = 1: Improved and simplified full monotonic constraint
2508 : ! LMT = 2: positive-definite constraint
2509 : ! LMT = 3: Quasi-monotone constraint
2510 :
2511 : if( lmt == 0 ) then
2512 :
2513 : ! Full constraint
2514 :
2515 : do jj = 1, jan
2516 : j = ja(jj)
2517 :
2518 : do i=1,im
2519 : if(dm(i,j) .eq. D0_0) then
2520 : ar(i,j) = p(i,j)
2521 : al(i,j) = p(i,j)
2522 : a6(i,j) = D0_0
2523 : else
2524 : da1 = ar(i,j) - al(i,j)
2525 : da2 = da1**2
2526 : a6da = a6(i,j)*da1
2527 : if(a6da .lt. -da2) then
2528 : a6(i,j) = D3_0*(al(i,j)-p(i,j))
2529 : ar(i,j) = al(i,j) - a6(i,j)
2530 : elseif(a6da .gt. da2) then
2531 : a6(i,j) = D3_0*(ar(i,j)-p(i,j))
2532 : al(i,j) = ar(i,j) - a6(i,j)
2533 : endif
2534 : endif
2535 : enddo
2536 :
2537 : enddo
2538 :
2539 : elseif( lmt == 1 ) then
2540 :
2541 : ! Improved (Lin 2001?) full constraint
2542 :
2543 : do jj = 1, jan
2544 : j = ja(jj)
2545 :
2546 : do i=1,im
2547 : da1 = dm(i,j) + dm(i,j)
2548 : dl = sign(min(abs(da1),abs(al(i,j)-p(i,j))), da1)
2549 : dr = sign(min(abs(da1),abs(ar(i,j)-p(i,j))), da1)
2550 : ar(i,j) = p(i,j) + dr
2551 : al(i,j) = p(i,j) - dl
2552 : a6(i,j) = D3_0*(dl-dr)
2553 : enddo
2554 :
2555 : enddo
2556 :
2557 : elseif( lmt == 2 ) then
2558 :
2559 : ! Positive definite constraint
2560 :
2561 : do jj = 1, jan
2562 : j = ja(jj)
2563 :
2564 : do i=1,im
2565 : if(abs(ar(i,j)-al(i,j)) .lt. -a6(i,j)) then
2566 : fmin = p(i,j) + D0_25*(ar(i,j)-al(i,j))**2/a6(i,j) + a6(i,j)*r12
2567 : if(fmin.lt.D0_0) then
2568 : if(p(i,j).lt.ar(i,j) .and. p(i,j).lt.al(i,j)) then
2569 : ar(i,j) = p(i,j)
2570 : al(i,j) = p(i,j)
2571 : a6(i,j) = D0_0
2572 : elseif(ar(i,j) .gt. al(i,j)) then
2573 : a6(i,j) = D3_0*(al(i,j)-p(i,j))
2574 : ar(i,j) = al(i,j) - a6(i,j)
2575 : else
2576 : a6(i,j) = D3_0*(ar(i,j)-p(i,j))
2577 : al(i,j) = ar(i,j) - a6(i,j)
2578 : endif
2579 : endif
2580 : endif
2581 : enddo
2582 :
2583 : enddo
2584 :
2585 : elseif(lmt .eq. 3) then
2586 :
2587 : ! Quasi-monotone constraint
2588 :
2589 : do jj = 1, jan
2590 : j = ja(jj)
2591 :
2592 : do i=1,im
2593 : da1 = D4_0*dm(i,j)
2594 : dl = sign(min(abs(da1),abs(al(i,j)-p(i,j))), da1)
2595 : dr = sign(min(abs(da1),abs(ar(i,j)-p(i,j))), da1)
2596 : ar(i,j) = p(i,j) + dr
2597 : al(i,j) = p(i,j) - dl
2598 : a6(i,j) = D3_0*(dl-dr)
2599 : enddo
2600 :
2601 : enddo
2602 :
2603 : endif
2604 : return
2605 : !EOC
2606 : end subroutine lmppmv
2607 : !-----------------------------------------------------------------------
2608 : #endif
2609 :
2610 : end module tp_core
|
