Line data Source code
1 : module mapz_module
2 :
3 : use shr_kind_mod, only : r8 => shr_kind_r8
4 : use FVperf_module, only : FVstartclock, FVstopclock
5 : use cam_abortutils, only : endrun
6 : use cam_logfile, only : iulog
7 :
8 : public map1_cubic_te, map1_ppm, mapn_ppm, mapn_ppm_tracer, ppme
9 :
10 : private
11 :
12 : real(r8), parameter :: D0_0 = 0.0_r8
13 : real(r8), parameter :: D1EM14 = 1.0e-14_r8
14 : real(r8), parameter :: D0_125 = 0.125_r8
15 : real(r8), parameter :: D0_1875 = 0.1875_r8
16 : real(r8), parameter :: D0_25 = 0.25_r8
17 : real(r8), parameter :: D0_5 = 0.5_r8
18 : real(r8), parameter :: D1_0 = 1.0_r8
19 : real(r8), parameter :: D1_5 = 1.5_r8
20 : real(r8), parameter :: D2_0 = 2.0_r8
21 : real(r8), parameter :: D3_0 = 3.0_r8
22 : real(r8), parameter :: D4_0 = 4.0_r8
23 : real(r8), parameter :: D5_0 = 5.0_r8
24 : real(r8), parameter :: D8_0 = 8.0_r8
25 : real(r8), parameter :: D12_0 = 12.0_r8
26 :
27 : contains
28 :
29 : !-----------------------------------------------------------------------
30 : !BOP
31 : ! !ROUTINE: map1_cubic_te --- Cubic Interpolation for TE mapping
32 : !
33 : ! !INTERFACE:
34 0 : subroutine map1_cubic_te ( km, pe1, q1, kn, pe2, q2, &
35 : ng_s, ng_n, itot, i1, i2, &
36 : j, jfirst, jlast, iv, kord)
37 : implicit none
38 :
39 : ! !INPUT PARAMETERS:
40 : integer, intent(in) :: i1 ! Starting longitude
41 : integer, intent(in) :: i2 ! Finishing longitude
42 : integer, intent(in) :: itot ! Total latitudes
43 : integer, intent(in) :: iv ! Mode: 0 == constituents 1 == ???
44 : integer, intent(in) :: kord ! Method order
45 : integer, intent(in) :: j ! Current latitude
46 : integer, intent(in) :: jfirst ! Starting latitude
47 : integer, intent(in) :: jlast ! Finishing latitude
48 : integer, intent(in) :: ng_s ! Ghosted latitudes south
49 : integer, intent(in) :: ng_n ! Ghosted latitudes north
50 : integer, intent(in) :: km ! Original vertical dimension
51 : integer, intent(in) :: kn ! Target vertical dimension
52 :
53 : real(r8), intent(in) :: pe1(itot,km+1) ! pressure at layer edges
54 : ! (from model top to bottom surface)
55 : ! in the original vertical coordinate
56 : real(r8), intent(in) :: pe2(itot,kn+1) ! pressure at layer edges
57 : ! (from model top to bottom surface)
58 : ! in the new vertical coordinate
59 :
60 : real(r8), intent(in) :: q1(itot,jfirst-ng_s:jlast+ng_n,km) ! Field input
61 :
62 : ! !INPUT/OUTPUT PARAMETERS:
63 : real(r8), intent(inout):: q2(itot,jfirst-ng_s:jlast+ng_n,kn) ! Field output
64 :
65 : ! !DESCRIPTION:
66 : !
67 : ! Perform Cubic Interpolation a given latitude
68 : ! pe1: pressure at layer edges (from model top to bottom surface)
69 : ! in the original vertical coordinate
70 : ! pe2: pressure at layer edges (from model top to bottom surface)
71 : ! in the new vertical coordinate
72 : !
73 : ! !REVISION HISTORY:
74 : ! 05.11.14 Takacs Initial Code
75 : !
76 : !EOP
77 : !-----------------------------------------------------------------------
78 : !BOC
79 : !
80 : ! !LOCAL VARIABLES:
81 0 : real(r8) qx(i1:i2,km)
82 0 : real(r8) logpl1(i1:i2,km)
83 0 : real(r8) logpl2(i1:i2,kn)
84 0 : real(r8) dlogp1(i1:i2,km)
85 0 : real(r8) vsum1(i1:i2)
86 0 : real(r8) vsum2(i1:i2)
87 : real(r8) am2,am1,ap0,ap1,P,PLP1,PLP0,PLM1,PLM2,DLP0,DLM1,DLM2
88 :
89 : integer i, k, LM2,LM1,LP0,LP1
90 :
91 : ! Initialization
92 : ! --------------
93 0 : do k=1,km
94 0 : qx(:,k) = q1(:,j,k)
95 0 : logpl1(:,k) = log( D0_5*(pe1(:,k)+pe1(:,k+1)) )
96 : enddo
97 0 : do k=1,kn
98 0 : logpl2(:,k) = log( D0_5*(pe2(:,k)+pe2(:,k+1)) )
99 : enddo
100 :
101 0 : do k=1,km-1
102 0 : dlogp1(:,k) = logpl1(:,k+1)-logpl1(:,k)
103 : enddo
104 :
105 : ! Compute vertical integral of Input TE
106 : ! -------------------------------------
107 0 : vsum1(:) = D0_0
108 0 : do i=i1,i2
109 0 : do k=1,km
110 0 : vsum1(i) = vsum1(i) + qx(i,k)*( pe1(i,k+1)-pe1(i,k) )
111 : enddo
112 0 : vsum1(i) = vsum1(i) / ( pe1(i,km+1)-pe1(i,1) )
113 : enddo
114 :
115 : ! Interpolate TE onto target Pressures
116 : ! ------------------------------------
117 0 : do i=i1,i2
118 0 : do k=1,kn
119 0 : LM1 = 1
120 : LP0 = 1
121 0 : do while( logpl1(i,LP0).lt.logpl2(i,k) .and. LP0.le.km )
122 0 : LP0 = LP0+1
123 : enddo
124 0 : LM1 = max(LP0-1,1)
125 0 : LP0 = min(LP0, km)
126 :
127 : ! Extrapolate Linearly in LogP above first model level
128 : ! ----------------------------------------------------
129 0 : if( LM1.eq.1 .and. LP0.eq.1 ) then
130 0 : q2(i,j,k) = qx(i,1) + ( qx(i,2)-qx(i,1) )*( logpl2(i,k)-logpl1(i,1) ) &
131 0 : /( logpl1(i,2)-logpl1(i,1) )
132 :
133 : ! Extrapolate Linearly in LogP below last model level
134 : ! ---------------------------------------------------
135 0 : else if( LM1.eq.km .and. LP0.eq.km ) then
136 0 : q2(i,j,k) = qx(i,km) + ( qx(i,km)-qx(i,km-1) )*( logpl2(i,k )-logpl1(i,km ) ) &
137 0 : /( logpl1(i,km)-logpl1(i,km-1) )
138 :
139 : ! Interpolate Linearly in LogP between levels 1 => 2 and km-1 => km
140 : ! -----------------------------------------------------------------
141 0 : else if( LM1.eq.1 .or. LP0.eq.km ) then
142 0 : q2(i,j,k) = qx(i,LP0) + ( qx(i,LM1)-qx(i,LP0) )*( logpl2(i,k )-logpl1(i,LP0) ) &
143 0 : /( logpl1(i,LM1)-logpl1(i,LP0) )
144 :
145 : ! Interpolate Cubicly in LogP between other model levels
146 : ! ------------------------------------------------------
147 : else
148 0 : LP1 = LP0+1
149 0 : LM2 = LM1-1
150 0 : P = logpl2(i,k)
151 0 : PLP1 = logpl1(i,LP1)
152 0 : PLP0 = logpl1(i,LP0)
153 0 : PLM1 = logpl1(i,LM1)
154 0 : PLM2 = logpl1(i,LM2)
155 0 : DLP0 = dlogp1(i,LP0)
156 0 : DLM1 = dlogp1(i,LM1)
157 0 : DLM2 = dlogp1(i,LM2)
158 :
159 0 : ap1 = (P-PLP0)*(P-PLM1)*(P-PLM2)/( DLP0*(DLP0+DLM1)*(DLP0+DLM1+DLM2) )
160 0 : ap0 = (PLP1-P)*(P-PLM1)*(P-PLM2)/( DLP0* DLM1 *( DLM1+DLM2) )
161 0 : am1 = (PLP1-P)*(PLP0-P)*(P-PLM2)/( DLM1* DLM2 *(DLP0+DLM1 ) )
162 0 : am2 = (PLP1-P)*(PLP0-P)*(PLM1-P)/( DLM2*(DLM1+DLM2)*(DLP0+DLM1+DLM2) )
163 :
164 0 : q2(i,j,k) = ap1*qx(i,LP1) + ap0*qx(i,LP0) + am1*qx(i,LM1) + am2*qx(i,LM2)
165 :
166 : endif
167 :
168 : enddo
169 : enddo
170 :
171 : ! Compute vertical integral of Output TE
172 : ! --------------------------------------
173 0 : vsum2(:) = D0_0
174 0 : do i=i1,i2
175 0 : do k=1,kn
176 0 : vsum2(i) = vsum2(i) + q2(i,j,k)*( pe2(i,k+1)-pe2(i,k) )
177 : enddo
178 0 : vsum2(i) = vsum2(i) / ( pe2(i,kn+1)-pe2(i,1) )
179 : enddo
180 :
181 : ! Adjust Final TE to conserve
182 : ! ---------------------------
183 0 : do i=i1,i2
184 0 : do k=1,kn
185 0 : q2(i,j,k) = q2(i,j,k) + vsum1(i)-vsum2(i)
186 : ! q2(i,j,k) = q2(i,j,k) * vsum1(i)/vsum2(i)
187 : enddo
188 : enddo
189 :
190 0 : return
191 : !EOC
192 : end subroutine map1_cubic_te
193 : !-----------------------------------------------------------------------
194 :
195 : !-----------------------------------------------------------------------
196 : !BOP
197 : ! !ROUTINE: map1_ppm --- Piecewise parabolic mapping, variant 1
198 : !
199 : ! !INTERFACE:
200 577584 : subroutine map1_ppm( km, pe1, q1, kn, pe2, q2, &
201 : ng_s, ng_n, itot, i1, i2, &
202 : j, jfirst, jlast, iv, kord)
203 :
204 : implicit none
205 :
206 : ! !INPUT PARAMETERS:
207 : integer, intent(in) :: i1 ! Starting longitude
208 : integer, intent(in) :: i2 ! Finishing longitude
209 : integer, intent(in) :: itot ! Total latitudes
210 : integer, intent(in) :: iv ! Mode: 0 == constituents 1 == ???
211 : integer, intent(in) :: kord ! Method order
212 : integer, intent(in) :: j ! Current latitude
213 : integer, intent(in) :: jfirst ! Starting latitude
214 : integer, intent(in) :: jlast ! Finishing latitude
215 : integer, intent(in) :: ng_s ! Ghosted latitudes south
216 : integer, intent(in) :: ng_n ! Ghosted latitudes north
217 : integer, intent(in) :: km ! Original vertical dimension
218 : integer, intent(in) :: kn ! Target vertical dimension
219 :
220 : real(r8), intent(in) :: pe1(itot,km+1) ! pressure at layer edges
221 : ! (from model top to bottom surface)
222 : ! in the original vertical coordinate
223 : real(r8), intent(in) :: pe2(itot,kn+1) ! pressure at layer edges
224 : ! (from model top to bottom surface)
225 : ! in the new vertical coordinate
226 : real(r8), intent(in) :: q1(itot,jfirst-ng_s:jlast+ng_n,km) ! Field input
227 :
228 : ! !INPUT/OUTPUT PARAMETERS:
229 : real(r8), intent(inout):: q2(itot,jfirst-ng_s:jlast+ng_n,kn) ! Field output
230 :
231 : ! !DESCRIPTION:
232 : !
233 : ! Perform piecewise parabolic method on a given latitude
234 : ! IV = 0: constituents
235 : ! pe1: pressure at layer edges (from model top to bottom surface)
236 : ! in the original vertical coordinate
237 : ! pe2: pressure at layer edges (from model top to bottom surface)
238 : ! in the new vertical coordinate
239 : !
240 : ! !REVISION HISTORY:
241 : ! 00.04.24 Lin Last modification
242 : ! 01.03.26 Sawyer Added ProTeX documentation
243 : ! 02.04.04 Sawyer incorporated latest FVGCM version
244 : ! 02.06.20 Sawyer made Q2 inout since the args for Q1/Q2 same
245 : ! 03.07.22 Parks Cleaned main loop, removed gotos
246 : ! 05.05.25 Sawyer Merged CAM and GEOS5 versions
247 : !
248 : !EOP
249 : !-----------------------------------------------------------------------
250 : !BOC
251 : !
252 : ! !LOCAL VARIABLES:
253 : real(r8) r3, r23
254 : parameter (r3 = D1_0/D3_0, r23 = D2_0/D3_0)
255 1155168 : real(r8) dp1(i1:i2,km)
256 1155168 : real(r8) q4(4,i1:i2,km)
257 :
258 1155168 : integer i, k, kk, kl, k0(i1:i2,0:kn+1), k0found
259 1155168 : real(r8) pl, pr, qsum, qsumk(i1:i2,kn), delp, esl
260 :
261 19060272 : do k=1,km
262 462644784 : do i=i1,i2
263 443584512 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
264 462067200 : q4(1,i,k) = q1(i,j,k)
265 : enddo
266 : enddo
267 :
268 : ! Mapping
269 :
270 : ! Compute vertical subgrid distribution
271 577584 : call ppm2m( q4, dp1, km, i1, i2, iv, kord )
272 :
273 : ! For each pe2(i,k), determine lowest pe1 interval = smallest k0 (= k0(i,k))
274 : ! such that pe1(i,k0) <= pe2(i,k) <= pe1(i,k0+1)
275 : ! Note that pe2(i,1)==pe1(i,1) and pe2(i,kn+1)==pe1(i,kn+1)
276 : ! Note also that pe1, pe2 are assumed to be monotonically increasing
277 : #if defined( UNICOSMP ) || defined ( NEC_SX )
278 : do kk = km, 1, -1
279 : do k = 1, kn+1
280 : do i = i1, i2
281 : if (pe2(i,k) <= pe1(i,kk+1)) then
282 : k0(i,k) = kk
283 : endif
284 : enddo
285 : enddo
286 : enddo
287 : #else
288 14439600 : do i = i1, i2
289 13862016 : k0(i,0) = 1
290 471886128 : do k = 1, kn+1
291 457446528 : k0found = -1
292 887169024 : do kk = k0(i,k-1), km
293 887169024 : if (pe2(i,k) <= pe1(i,kk+1)) then
294 457446528 : k0(i,k) = kk
295 457446528 : k0found = kk
296 457446528 : exit
297 : endif
298 : enddo
299 471308544 : if (k0found .lt. 0) then
300 0 : write(iulog,*) 'mapz error - k0found i j k (kk,pe1,pe2) = ', &
301 0 : k0found, i, j, k, (kk,pe1(i,kk),pe2(i,kk),kk=1,km+1)
302 0 : call endrun('MAPZ_MODULE')
303 0 : return
304 : endif
305 : enddo
306 : enddo
307 : #endif
308 :
309 : ! Interpolate
310 19060272 : do k = 1, kn
311 :
312 : ! Prepare contribution between pe1(i,ko(i,k)+1) and pe1(i,k0(i,k+1))
313 462067200 : qsumk(:,k) = D0_0
314 462067200 : do i = i1, i2
315 496422010 : do kl = k0(i,k)+1, k0(i,k+1)-1
316 477939322 : qsumk(i,k) = qsumk(i,k) + dp1(i,kl)*q4(1,i,kl)
317 : enddo
318 : enddo
319 :
320 462644784 : do i = i1, i2
321 443584512 : kk = k0(i,k)
322 : ! Consider contribution between pe1(i,kk) and pe2(i,k)
323 443584512 : pl = (pe2(i,k)-pe1(i,kk)) / dp1(i,kk)
324 : ! Check to see if pe2(i,k+1) and pe2(i,k) are in same pe1 interval
325 462067200 : if (k0(i,k+1) == k0(i,k)) then
326 48216826 : pr = (pe2(i,k+1)-pe1(i,kk)) / dp1(i,kk)
327 48216826 : q2(i,j,k) = q4(2,i,kk) + D0_5*(q4(4,i,kk)+q4(3,i,kk)-q4(2,i,kk)) &
328 96433652 : *(pr+pl) - q4(4,i,kk)*r3*(pr*(pr+pl)+pl**2)
329 : else
330 : ! Consider contribution between pe2(i,k) and pe1(i,kk+1)
331 395367686 : qsum = (pe1(i,kk+1)-pe2(i,k))*(q4(2,i,kk)+D0_5*(q4(4,i,kk)+ &
332 : q4(3,i,kk)-q4(2,i,kk))*(D1_0+pl)-q4(4,i,kk)* &
333 395367686 : (r3*(D1_0+pl*(D1_0+pl))))
334 : ! Next consider contribution between pe1(i,kk+1) and pe1(i,k0(i,k+1))
335 395367686 : qsum = qsum + qsumk(i,k)
336 : ! Now consider contribution between pe1(i,k0(i,k+1)) and pe2(i,k+1)
337 395367686 : kl = k0(i,k+1)
338 395367686 : delp = pe2(i,k+1)-pe1(i,kl)
339 395367686 : esl = delp / dp1(i,kl)
340 : qsum = qsum + delp*(q4(2,i,kl)+D0_5*esl* &
341 395367686 : (q4(3,i,kl)-q4(2,i,kl)+q4(4,i,kl)*(D1_0-r23*esl)))
342 395367686 : q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
343 : endif
344 : enddo
345 : enddo
346 :
347 : return
348 : !EOC
349 : end subroutine map1_ppm
350 : !-----------------------------------------------------------------------
351 :
352 : !-----------------------------------------------------------------------
353 : !BOP
354 : ! !ROUTINE: mapn_ppm --- Piecewise parabolic mapping, variant 1
355 : !
356 : ! !INTERFACE:
357 0 : subroutine mapn_ppm(km, pe1, q1, nq, &
358 0 : kn, pe2, q2, ng_s, ng_n, &
359 : itot, i1, i2, j, &
360 : jfirst, jlast, iv, kord)
361 :
362 : ! !USES:
363 : implicit none
364 :
365 : ! !INPUT PARAMETERS:
366 : integer, intent(in) :: i1 ! Starting longitude
367 : integer, intent(in) :: i2 ! Finishing longitude
368 : integer, intent(in) :: itot ! Total latitudes
369 : integer, intent(in) :: iv ! Mode: 0 == constituents 1 == ???
370 : integer, intent(in) :: kord ! Method order
371 : integer, intent(in) :: j ! Current latitude
372 : integer, intent(in) :: jfirst ! Starting latitude
373 : integer, intent(in) :: jlast ! Finishing latitude
374 : integer, intent(in) :: ng_s ! Ghosted latitudes south
375 : integer, intent(in) :: ng_n ! Ghosted latitudes north
376 : integer, intent(in) :: km ! Original vertical dimension
377 : integer, intent(in) :: kn ! Target vertical dimension
378 : integer, intent(in) :: nq ! Number of tracers
379 :
380 : real(r8), intent(in) :: pe1(itot,km+1) ! pressure at layer edges
381 : ! (from model top to bottom surface)
382 : ! in the original vertical coordinate
383 : real(r8), intent(in) :: pe2(itot,kn+1) ! pressure at layer edges
384 : ! (from model top to bottom surface)
385 : ! in the new vertical coordinate
386 : real(r8), intent(in) :: q1(itot,jfirst-ng_s:jlast+ng_n,km,nq) ! Field input
387 : ! !INPUT/OUTPUT PARAMETERS:
388 : real(r8), intent(inout):: q2(itot,jfirst-ng_s:jlast+ng_n,kn,nq) ! Field output
389 :
390 : ! !DESCRIPTION:
391 : !
392 : ! Perform piecewise parabolic method on a given latitude
393 : ! IV = 0: constituents
394 : ! pe1: pressure at layer edges (from model top to bottom surface)
395 : ! in the original vertical coordinate
396 : ! pe2: pressure at layer edges (from model top to bottom surface)
397 : ! in the new vertical coordinate
398 : !
399 : ! !REVISION HISTORY:
400 : ! 02.04.04 Sawyer incorporated latest FVGCM version, ProTeX
401 : ! 02.06.20 Sawyer made Q2 inout since the args for Q1/Q2 same
402 : ! 03.07.22 Parks Cleaned main loop, removed gotos
403 : !
404 : !EOP
405 : !-----------------------------------------------------------------------
406 : !BOC
407 : !
408 : ! !LOCAL VARIABLES:
409 : real(r8) r3, r23
410 : parameter (r3 = D1_0/D3_0, r23 = D2_0/D3_0)
411 0 : real(r8) dp1(i1:i2,km)
412 0 : real(r8) q4(4,i1:i2,km)
413 :
414 0 : integer i, k, kk, kl, k0(i1:i2,0:kn+1), iq
415 0 : real(r8) pl, pr, qsum, qsumk(i1:i2,kn), delp, esl
416 :
417 0 : do k=1,km
418 0 : do i=i1,i2
419 0 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
420 : enddo
421 : enddo
422 :
423 : ! Mapping
424 :
425 : ! For each pe2(i,k), determine lowest pe1 interval = smallest k0 (= k0(i,k))
426 : ! such that pe1(i,k0) <= pe2(i,k) <= pe1(i,k0+1)
427 : ! Note that pe2(i,1)==pe1(i,1) and pe2(i,kn+1)==pe1(i,kn+1)
428 : ! Note also that pe1, pe2 are assumed to be monotonically increasing
429 : #if defined( UNICOSMP ) || defined ( NEC_SX )
430 : do kk = km, 1, -1
431 : do k = 1, kn+1
432 : do i = i1, i2
433 : if (pe2(i,k) <= pe1(i,kk+1)) then
434 : k0(i,k) = kk
435 : endif
436 : enddo
437 : enddo
438 : enddo
439 : #else
440 0 : do i = i1, i2
441 0 : k0(i,0) = 1
442 0 : do k = 1, kn+1
443 0 : do kk = k0(i,k-1), km
444 0 : if (pe2(i,k) <= pe1(i,kk+1)) then
445 0 : k0(i,k) = kk
446 0 : exit
447 : endif
448 : enddo
449 : enddo
450 : enddo
451 : #endif
452 :
453 0 : do iq=1,nq
454 :
455 0 : do k=1,km
456 0 : do i=i1,i2
457 0 : q4(1,i,k) = q1(i,j,k,iq)
458 : enddo
459 : enddo
460 :
461 : ! Compute vertical subgrid distribution
462 0 : call ppm2m( q4, dp1, km, i1, i2, iv, kord )
463 : ! Interpolate
464 0 : do k = 1, kn
465 :
466 : ! Prepare contribution between pe1(i,ko(i,k)+1) and pe1(i,k0(i,k+1))
467 0 : qsumk(:,k) = D0_0
468 0 : do i = i1, i2
469 0 : do kl = k0(i,k)+1, k0(i,k+1)-1
470 0 : qsumk(i,k) = qsumk(i,k) + dp1(i,kl)*q4(1,i,kl)
471 : enddo
472 : enddo
473 :
474 0 : do i = i1, i2
475 0 : kk = k0(i,k)
476 : ! Consider contribution between pe1(i,kk) and pe2(i,k)
477 0 : pl = (pe2(i,k)-pe1(i,kk)) / dp1(i,kk)
478 : ! Check to see if pe2(i,k+1) and pe2(i,k) are in same pe1 interval
479 0 : if (k0(i,k+1) == k0(i,k)) then
480 0 : pr = (pe2(i,k+1)-pe1(i,kk)) / dp1(i,kk)
481 0 : q2(i,j,k,iq) = q4(2,i,kk) + D0_5*(q4(4,i,kk)+q4(3,i,kk)-q4(2,i,kk)) &
482 0 : *(pr+pl) - q4(4,i,kk)*r3*(pr*(pr+pl)+pl**2)
483 : else
484 : ! Consider contribution between pe2(i,k) and pe1(i,kk+1)
485 0 : qsum = (pe1(i,kk+1)-pe2(i,k))*(q4(2,i,kk)+D0_5*(q4(4,i,kk)+ &
486 : q4(3,i,kk)-q4(2,i,kk))*(D1_0+pl)-q4(4,i,kk)* &
487 0 : (r3*(D1_0+pl*(D1_0+pl))))
488 : ! Next consider contribution between pe1(i,kk+1) and pe1(i,k0(i,k+1))
489 0 : qsum = qsum + qsumk(i,k)
490 : ! Now consider contribution between pe1(i,k0(i,k+1)) and pe2(i,k+1)
491 0 : kl = k0(i,k+1)
492 0 : delp = pe2(i,k+1)-pe1(i,kl)
493 0 : esl = delp / dp1(i,kl)
494 : qsum = qsum + delp*(q4(2,i,kl)+D0_5*esl* &
495 0 : (q4(3,i,kl)-q4(2,i,kl)+q4(4,i,kl)*(D1_0-r23*esl)))
496 0 : q2(i,j,k,iq) = qsum / ( pe2(i,k+1) - pe2(i,k) )
497 : endif
498 : enddo
499 : enddo
500 :
501 : enddo
502 :
503 0 : return
504 : !EOC
505 : end subroutine mapn_ppm
506 : !-----------------------------------------------------------------------
507 :
508 :
509 : !-----------------------------------------------------------------------
510 : !BOP
511 : ! !ROUTINE: mapn_ppm_tracer --- Piecewise parabolic mapping, multiple tracers
512 : !
513 : ! !INTERFACE:
514 96768 : subroutine mapn_ppm_tracer(km, pe1, tracer, nq, &
515 96768 : kn, pe2, i1, i2, j, &
516 : ifirst, ilast, jfirst, jlast, iv, kord)
517 :
518 : ! !USES:
519 : implicit none
520 :
521 : ! !INPUT PARAMETERS:
522 : integer, intent(in) :: i1 ! Starting longitude
523 : integer, intent(in) :: i2 ! Finishing longitude
524 : integer, intent(in) :: iv ! Mode: 0 == constituents 1 == ???
525 : integer, intent(in) :: kord ! Method order
526 : integer, intent(in) :: j ! Current latitude
527 : integer, intent(in) :: ifirst ! Starting segment
528 : integer, intent(in) :: ilast ! Finishing segment
529 : integer, intent(in) :: jfirst ! Starting latitude
530 : integer, intent(in) :: jlast ! Finishing latitude
531 : integer, intent(in) :: km ! Original vertical dimension
532 : integer, intent(in) :: kn ! Target vertical dimension
533 : integer, intent(in) :: nq ! Number of tracers
534 :
535 : real(r8), intent(in) :: pe1(ifirst:ilast,km+1) ! pressure at layer edges
536 : ! (from model top to bottom surface)
537 : ! in the original vertical coordinate
538 : real(r8), intent(in) :: pe2(ifirst:ilast,kn+1) ! pressure at layer edges
539 : ! (from model top to bottom surface)
540 : ! in the new vertical coordinate
541 : ! !INPUT/OUTPUT PARAMETERS:
542 : real (r8), intent(inout):: tracer(ifirst:ilast,jfirst:jlast,km,nq) ! Field output
543 :
544 : ! !DESCRIPTION:
545 : !
546 : ! Perform piecewise parabolic method on a given latitude
547 : ! IV = 0: constituents
548 : ! pe1: pressure at layer edges (from model top to bottom surface)
549 : ! in the original vertical coordinate
550 : ! pe2: pressure at layer edges (from model top to bottom surface)
551 : ! in the new vertical coordinate
552 : !
553 : ! !REVISION HISTORY:
554 : ! 05.03.20 Sawyer Created from mapn_ppm
555 : ! 05.04.04 Sawyer Simplified indexing, removed ifirst
556 : ! 05.04.12 Sawyer Added r4/r8 distinction
557 : ! 05.10.12 Worley Made mapn_ppm_tracer vector-friendly
558 : !
559 : !EOP
560 : !-----------------------------------------------------------------------
561 : !BOC
562 : !
563 : ! !LOCAL VARIABLES:
564 :
565 : real(r8) r3, r23
566 : parameter (r3 = D1_0/D3_0, r23 = D2_0/D3_0)
567 :
568 193536 : real(r8) dp1(i1:i2,km)
569 193536 : real(r8) q4(4,i1:i2,km)
570 :
571 193536 : integer i, k, kk, kl, k0(i1:i2,0:kn+1), iq
572 193536 : real(r8) pl, pr, qsum, qsumk(i1:i2,kn), delp, esl
573 :
574 3193344 : do k=1,km
575 77511168 : do i=i1,i2
576 77414400 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
577 : enddo
578 : enddo
579 :
580 : ! Mapping
581 :
582 : ! For each pe2(i,k), determine lowest pe1 interval = smallest k0 (= k0(i,k))
583 : ! such that pe1(i,k0) <= pe2(i,k) <= pe1(i,k0+1)
584 : ! Note that pe2(i,1)==pe1(i,1) and pe2(i,kn+1)==pe1(i,kn+1)
585 : ! Note also that pe1, pe2 are assumed to be monotonically increasing
586 : #if defined( UNICOSMP ) || defined ( NEC_SX )
587 : do kk = km, 1, -1
588 : do k = 1, kn+1
589 : do i = i1, i2
590 : if (pe2(i,k) <= pe1(i,kk+1)) then
591 : k0(i,k) = kk
592 : endif
593 : enddo
594 : enddo
595 : enddo
596 : #else
597 2419200 : do i = i1, i2
598 2322432 : k0(i,0) = 1
599 79059456 : do k = 1, kn+1
600 150958080 : do kk = k0(i,k-1), km
601 148635648 : if (pe2(i,k) <= pe1(i,kk+1)) then
602 76640256 : k0(i,k) = kk
603 76640256 : exit
604 : endif
605 : enddo
606 : enddo
607 : enddo
608 : #endif
609 :
610 23321088 : do iq=1,nq
611 766402560 : do k=1,km
612 18602680320 : do i=i1,i2
613 18579456000 : q4(1,i,k) = tracer(i,j,k,iq)
614 : enddo
615 : enddo
616 :
617 : ! Compute vertical subgrid distribution
618 23224320 : call ppm2m( q4, dp1, km, i1, i2, iv, kord )
619 :
620 : ! Interpolate
621 766499328 : do k = 1, kn
622 :
623 : ! Prepare contribution between pe1(i,ko(i,k)+1) and pe1(i,k0(i,k+1))
624 18579456000 : qsumk(:,k) = D0_0
625 18579456000 : do i = i1, i2
626 19976881440 : do kl = k0(i,k)+1, k0(i,k+1)-1
627 19233703200 : qsumk(i,k) = qsumk(i,k) + dp1(i,kl)*q4(1,i,kl)
628 : enddo
629 : enddo
630 :
631 18602680320 : do i = i1, i2
632 17836277760 : kk = k0(i,k)
633 : ! Consider contribution between pe1(i,kk) and pe2(i,k)
634 17836277760 : pl = (pe2(i,k)-pe1(i,kk)) / dp1(i,kk)
635 : ! Check to see if pe2(i,k+1) and pe2(i,k) are in same pe1 interval
636 18579456000 : if (k0(i,k+1) == k0(i,k)) then
637 1954809120 : pr = (pe2(i,k+1)-pe1(i,kk)) / dp1(i,kk)
638 3909618240 : tracer(i,j,k,iq) = q4(2,i,kk) &
639 : + D0_5*(q4(4,i,kk)+q4(3,i,kk)-q4(2,i,kk)) &
640 5864427360 : *(pr+pl)-q4(4,i,kk)*r3*(pr*(pr+pl)+pl**2)
641 : else
642 : ! Consider contribution between pe2(i,k) and pe1(i,kk+1)
643 15881468640 : qsum = (pe1(i,kk+1)-pe2(i,k))*(q4(2,i,kk)+D0_5*(q4(4,i,kk)+ &
644 : q4(3,i,kk)-q4(2,i,kk))*(D1_0+pl)-q4(4,i,kk)* &
645 15881468640 : (r3*(D1_0+pl*(D1_0+pl))))
646 : ! Next consider contribution between pe1(i,kk+1) and pe1(i,k0(i,k+1))
647 15881468640 : qsum = qsum + qsumk(i,k)
648 : ! Now consider contribution between pe1(i,k0(i,k+1)) and pe2(i,k+1)
649 15881468640 : kl = k0(i,k+1)
650 15881468640 : delp = pe2(i,k+1)-pe1(i,kl)
651 15881468640 : esl = delp / dp1(i,kl)
652 : qsum = qsum + delp*(q4(2,i,kl)+D0_5*esl* &
653 15881468640 : (q4(3,i,kl)-q4(2,i,kl)+q4(4,i,kl)*(D1_0-r23*esl)))
654 15881468640 : tracer(i,j,k,iq) = qsum / ( pe2(i,k+1) - pe2(i,k) )
655 : endif
656 : enddo
657 : enddo
658 :
659 : enddo ! do iq=1,nq
660 :
661 96768 : return
662 : !EOC
663 : end subroutine mapn_ppm_tracer
664 : !-----------------------------------------------------------------------
665 :
666 :
667 : !-----------------------------------------------------------------------
668 : !BOP
669 : ! !ROUTINE: ppm2m --- Piecewise parabolic method for fields
670 : !
671 : ! !INTERFACE:
672 23801904 : subroutine ppm2m(a4, delp, km, i1, i2, iv, kord)
673 :
674 : ! !USES:
675 : implicit none
676 :
677 : ! !INPUT PARAMETERS:
678 : integer, intent(in):: iv ! iv =-1: winds
679 : ! iv = 0: positive definite scalars
680 : ! iv = 1: others
681 : integer, intent(in):: i1 ! Starting longitude
682 : integer, intent(in):: i2 ! Finishing longitude
683 : integer, intent(in):: km ! vertical dimension
684 : integer, intent(in):: kord ! Order (or more accurately method no.):
685 : !
686 : real (r8), intent(in):: delp(i1:i2,km) ! layer pressure thickness
687 :
688 : ! !INPUT/OUTPUT PARAMETERS:
689 : real (r8), intent(inout):: a4(4,i1:i2,km) ! Interpolated values
690 :
691 : ! !DESCRIPTION:
692 : !
693 : ! Perform the piecewise parabolic method
694 : !
695 : ! !REVISION HISTORY:
696 : ! ??.??.?? Lin Creation
697 : ! 02.04.04 Sawyer Newest release from FVGCM
698 : ! 02.04.23 Sawyer Incorporated minor algorithmic change to
699 : ! maintain CAM zero diffs (see comments inline)
700 : !
701 : !EOP
702 : !-----------------------------------------------------------------------
703 : !BOC
704 : !
705 : ! !LOCAL VARIABLES:
706 : ! local arrays:
707 47603808 : real(r8) dc(i1:i2,km)
708 47603808 : real(r8) h2(i1:i2,km)
709 47603808 : real(r8) delq(i1:i2,km)
710 47603808 : real(r8) df2(i1:i2,km)
711 47603808 : real(r8) d4(i1:i2,km)
712 :
713 : ! local scalars:
714 : real(r8) fac
715 : real(r8) a1, a2, c1, c2, c3, d1, d2
716 : real(r8) qmax, qmin, cmax, cmin
717 : real(r8) qm, dq, tmp
718 :
719 : integer i, k, km1, lmt
720 : real(r8) qmp, pmp
721 : real(r8) lac
722 : integer it
723 :
724 23801904 : km1 = km - 1
725 23801904 : it = i2 - i1 + 1
726 :
727 761660928 : do k=2,km
728 18470277504 : do i=i1,i2
729 17708616576 : delq(i,k-1) = a4(1,i,k) - a4(1,i,k-1)
730 18446475600 : d4(i,k ) = delp(i,k-1) + delp(i,k)
731 : enddo
732 : enddo
733 :
734 737859024 : do k=2,km1
735 17875229904 : do i=i1,i2
736 17137370880 : c1 = (delp(i,k-1)+D0_5*delp(i,k))/d4(i,k+1)
737 17137370880 : c2 = (delp(i,k+1)+D0_5*delp(i,k))/d4(i,k)
738 : tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / &
739 17137370880 : (d4(i,k)+delp(i,k+1))
740 17137370880 : qmax = max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1)) - a4(1,i,k)
741 17137370880 : qmin = a4(1,i,k) - min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))
742 17137370880 : dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp)
743 17851428000 : df2(i,k) = tmp
744 : enddo
745 : enddo
746 :
747 : !****6***0*********0*********0*********0*********0*********0**********72
748 : ! 4th order interpolation of the provisional cell edge value
749 : !****6***0*********0*********0*********0*********0*********0**********72
750 :
751 714057120 : do k=3,km1
752 17280182304 : do i=i1,i2
753 16566125184 : c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k)
754 16566125184 : a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1))
755 16566125184 : a2 = d4(i,k+1) / (d4(i,k) + delp(i,k))
756 : a4(2,i,k) = a4(1,i,k-1) + c1 + D2_0/(d4(i,k-1)+d4(i,k+1)) * &
757 : ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - &
758 17256380400 : delp(i,k-1)*a1*dc(i,k ) )
759 : enddo
760 : enddo
761 :
762 23801904 : call steepz(i1, i2, km, a4, df2, dc, delq, delp, d4)
763 :
764 : ! Area preserving cubic with 2nd deriv. = 0 at the boundaries
765 : ! Top
766 595047600 : do i=i1,i2
767 571245696 : d1 = delp(i,1)
768 571245696 : d2 = delp(i,2)
769 571245696 : qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2)
770 571245696 : dq = D2_0*(a4(1,i,2)-a4(1,i,1)) / (d1+d2)
771 571245696 : c1 = D4_0*(a4(2,i,3)-qm-d2*dq) / ( d2*(D2_0*d2*d2+d1*(d2+D3_0*d1)) )
772 571245696 : c3 = dq - D0_5*c1*(d2*(D5_0*d1+d2)-D3_0*d1**2)
773 571245696 : a4(2,i,2) = qm - D0_25*c1*d1*d2*(d2+D3_0*d1)
774 571245696 : a4(2,i,1) = d1*(D2_0*c1*d1**2-c3) + a4(2,i,2)
775 571245696 : dc(i,1) = 0.5_r8*(a4(1,i,1) - a4(2,i,1))
776 : ! No over- and undershoot condition
777 571245696 : cmax = max(a4(1,i,1), a4(1,i,2))
778 571245696 : cmin = min(a4(1,i,1), a4(1,i,2))
779 571245696 : a4(2,i,2) = max(cmin,a4(2,i,2))
780 595047600 : a4(2,i,2) = min(cmax,a4(2,i,2))
781 : enddo
782 :
783 23801904 : if( iv == 0 ) then
784 580608000 : do i=i1,i2
785 : !
786 : ! WS: 02.04.23 Algorithmic difference with FVGCM. FVGCM does this:
787 : !
788 : !!! a4(2,i,1) = a4(1,i,1)
789 : !!! a4(3,i,1) = a4(1,i,1)
790 : !
791 : ! CAM does this:
792 : !
793 557383680 : a4(2,i,1) = max(D0_0,a4(2,i,1))
794 580608000 : a4(2,i,2) = max(D0_0,a4(2,i,2))
795 : enddo
796 577584 : elseif ( iv == -1 ) then
797 : ! Winds:
798 384048 : if( km > 32 ) then
799 0 : do i=i1,i2
800 : ! More dampping: top layer as the sponge
801 0 : a4(2,i,1) = a4(1,i,1)
802 0 : a4(3,i,1) = a4(1,i,1)
803 : enddo
804 : else
805 9601200 : do i=i1,i2
806 9601200 : if( a4(1,i,1)*a4(2,i,1) <= D0_0 ) then
807 420062 : a4(2,i,1) = D0_0
808 : else
809 : a4(2,i,1) = sign(min(abs(a4(1,i,1)), &
810 : abs(a4(2,i,1))), &
811 8797090 : a4(1,i,1) )
812 : endif
813 : enddo
814 : endif
815 : endif
816 :
817 : ! Bottom
818 : ! Area preserving cubic with 2nd deriv. = 0 at the surface
819 595047600 : do i=i1,i2
820 571245696 : d1 = delp(i,km)
821 571245696 : d2 = delp(i,km1)
822 571245696 : qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2)
823 571245696 : dq = D2_0*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2)
824 571245696 : c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(D2_0*d2*d2+d1*(d2+D3_0*d1)))
825 571245696 : c3 = dq - D2_0*c1*(d2*(D5_0*d1+d2)-D3_0*d1**2)
826 571245696 : a4(2,i,km) = qm - c1*d1*d2*(d2+D3_0*d1)
827 571245696 : a4(3,i,km) = d1*(D8_0*c1*d1**2-c3) + a4(2,i,km)
828 571245696 : dc(i,km) = 0.5_r8*(a4(3,i,km) - a4(1,i,km))
829 : ! No over- and under-shoot condition
830 571245696 : cmax = max(a4(1,i,km), a4(1,i,km1))
831 571245696 : cmin = min(a4(1,i,km), a4(1,i,km1))
832 571245696 : a4(2,i,km) = max(cmin,a4(2,i,km))
833 595047600 : a4(2,i,km) = min(cmax,a4(2,i,km))
834 : enddo
835 :
836 : ! Enforce constraint at the surface
837 :
838 23801904 : if ( iv == 0 ) then
839 : ! Positive definite scalars:
840 580608000 : do i=i1,i2
841 580608000 : a4(3,i,km) = max(D0_0, a4(3,i,km))
842 : enddo
843 577584 : elseif ( iv == -1 ) then
844 : ! Winds:
845 9601200 : do i=i1,i2
846 9601200 : if( a4(1,i,km)*a4(3,i,km) <= D0_0 ) then
847 998754 : a4(3,i,km) = D0_0
848 : else
849 : a4(3,i,km) = sign( min(abs(a4(1,i,km)), &
850 : abs(a4(3,i,km))), &
851 8218398 : a4(1,i,km) )
852 : endif
853 : enddo
854 : endif
855 :
856 761660928 : do k=1,km1
857 18470277504 : do i=i1,i2
858 18446475600 : a4(3,i,k) = a4(2,i,k+1)
859 : enddo
860 : enddo
861 :
862 : ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
863 :
864 : ! Top 2 and bottom 2 layers always use monotonic mapping
865 71405712 : do k=1,2
866 1190095200 : do i=i1,i2
867 1190095200 : a4(4,i,k) = D3_0*(D2_0*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
868 : enddo
869 71405712 : call kmppm(dc(i1,k), a4(1,i1,k), it, 0)
870 : enddo
871 :
872 23801904 : if(kord >= 7) then
873 : !****6***0*********0*********0*********0*********0*********0**********72
874 : ! Huynh's 2nd constraint
875 : !****6***0*********0*********0*********0*********0*********0**********72
876 719953920 : do k=2, km1
877 17441464320 : do i=i1,i2
878 : ! Method#1
879 : ! h2(i,k) = delq(i,k) - delq(i,k-1)
880 : ! Method#2
881 : ! h2(i,k) = D2_0*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1))
882 : ! & / ( delp(i,k)+D0_5*(delp(i,k-1)+delp(i,k+1)) )
883 : ! & * delp(i,k)**2
884 : ! Method#3
885 17418240000 : h2(i,k) = dc(i,k+1) - dc(i,k-1)
886 : enddo
887 : enddo
888 :
889 23224320 : if( kord == 7 ) then
890 : fac = D1_5 ! original quasi-monotone
891 : else
892 0 : fac = D0_125 ! full monotone
893 : endif
894 :
895 673505280 : do k=3, km-2
896 16257024000 : do i=i1,i2
897 : ! Right edges
898 : ! qmp = a4(1,i,k) + D2_0*delq(i,k-1)
899 : ! lac = a4(1,i,k) + fac*h2(i,k-1) + D0_5*delq(i,k-1)
900 : !
901 15606743040 : pmp = D2_0*dc(i,k)
902 15606743040 : qmp = a4(1,i,k) + pmp
903 15606743040 : lac = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k)
904 15606743040 : qmin = min(a4(1,i,k), qmp, lac)
905 15606743040 : qmax = max(a4(1,i,k), qmp, lac)
906 15606743040 : a4(3,i,k) = min(max(a4(3,i,k), qmin), qmax)
907 : ! Left edges
908 : ! qmp = a4(1,i,k) - D2_0*delq(i,k)
909 : ! lac = a4(1,i,k) + fac*h2(i,k+1) - D0_5*delq(i,k)
910 : !
911 15606743040 : qmp = a4(1,i,k) - pmp
912 15606743040 : lac = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k)
913 15606743040 : qmin = min(a4(1,i,k), qmp, lac)
914 15606743040 : qmax = max(a4(1,i,k), qmp, lac)
915 15606743040 : a4(2,i,k) = min(max(a4(2,i,k), qmin), qmax)
916 : ! Recompute A6
917 16257024000 : a4(4,i,k) = D3_0*(D2_0*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
918 : enddo
919 : ! Additional constraint to prevent negatives when kord=7
920 673505280 : if (iv == 0 .and. kord == 7) then
921 650280960 : call kmppm(dc(i1,k), a4(1,i1,k), it, 2)
922 : endif
923 : enddo
924 :
925 : else
926 :
927 577584 : lmt = kord - 3
928 577584 : lmt = max(0, lmt)
929 577584 : if (iv == 0) lmt = min(2, lmt)
930 :
931 16749936 : do k=3, km-2
932 16172352 : if( kord /= 4) then
933 0 : do i=i1,i2
934 0 : a4(4,i,k) = D3_0*(D2_0*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
935 : enddo
936 : endif
937 16749936 : call kmppm(dc(i1,k), a4(1,i1,k), it, lmt)
938 : enddo
939 : endif
940 :
941 71405712 : do k=km1,km
942 1190095200 : do i=i1,i2
943 1190095200 : a4(4,i,k) = D3_0*(D2_0*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
944 : enddo
945 71405712 : call kmppm(dc(i1,k), a4(1,i1,k), it, 0)
946 : enddo
947 :
948 23801904 : return
949 : !EOC
950 : end subroutine ppm2m
951 : !-----------------------------------------------------------------------
952 :
953 :
954 : !-----------------------------------------------------------------------
955 : !BOP
956 : ! !ROUTINE: ppme --- PPM scheme at vertical edges
957 : !
958 : ! !INTERFACE:
959 0 : subroutine ppme(p,qe,delp,im,km)
960 : ! !USES:
961 : use shr_kind_mod, only : r8 => shr_kind_r8, i4 => shr_kind_i4
962 : implicit none
963 :
964 : ! !INPUT PARAMETERS:
965 : integer, intent(in) :: im, km
966 : real(r8), intent(in) :: p(im,km), delp(im,km)
967 :
968 : ! !INPUT/OUTPUT PARAMETERS:
969 : real(r8), intent(out) :: qe(im,km+1)
970 :
971 : ! !DESCRIPTION:
972 : !
973 : !
974 : ! !REVISION HISTORY:
975 : ! 05.06.13 Sawyer Inserted file ppme.F90 here, added ProTeX
976 : !
977 : !EOP
978 : !-----------------------------------------------------------------------
979 : !BOC
980 :
981 : integer(i4) km1
982 : integer(i4) i, k
983 : ! local arrays.
984 0 : real(r8) dc(im,km),delq(im,km), a6(im,km)
985 : real(r8) c1, c2, c3, tmp, qmax, qmin
986 : real(r8) a1, a2, s1, s2, s3, s4, ss3, s32, s34, s42
987 : real(r8) a3, b2, sc, dm, d1, d2, f1, f2, f3, f4
988 : real(r8) qm, dq
989 :
990 0 : km1 = km - 1
991 :
992 0 : do 500 k=2,km
993 0 : do 500 i=1,im
994 0 : 500 a6(i,k) = delp(i,k-1) + delp(i,k)
995 :
996 0 : do 1000 k=1,km1
997 0 : do 1000 i=1,im
998 0 : delq(i,k) = p(i,k+1) - p(i,k)
999 0 : 1000 continue
1000 :
1001 0 : do 1220 k=2,km1
1002 0 : do 1220 i=1,im
1003 0 : c1 = (delp(i,k-1)+D0_5*delp(i,k))/a6(i,k+1)
1004 0 : c2 = (delp(i,k+1)+D0_5*delp(i,k))/a6(i,k)
1005 : tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / &
1006 0 : (a6(i,k)+delp(i,k+1))
1007 0 : qmax = max(p(i,k-1),p(i,k),p(i,k+1)) - p(i,k)
1008 0 : qmin = p(i,k) - min(p(i,k-1),p(i,k),p(i,k+1))
1009 0 : dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp)
1010 0 : 1220 continue
1011 :
1012 : !****6***0*********0*********0*********0*********0*********0**********72
1013 : ! 4th order interpolation of the provisional cell edge value
1014 : !****6***0*********0*********0*********0*********0*********0**********72
1015 :
1016 0 : do 12 k=3,km1
1017 0 : do 12 i=1,im
1018 0 : c1 = delq(i,k-1)*delp(i,k-1) / a6(i,k)
1019 0 : a1 = a6(i,k-1) / (a6(i,k) + delp(i,k-1))
1020 0 : a2 = a6(i,k+1) / (a6(i,k) + delp(i,k))
1021 0 : qe(i,k) = p(i,k-1) + c1 + D2_0/(a6(i,k-1)+a6(i,k+1)) * &
1022 : ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - &
1023 0 : delp(i,k-1)*a1*dc(i,k ) )
1024 0 : 12 continue
1025 :
1026 : ! three-cell parabolic subgrid distribution at model top
1027 :
1028 0 : do 10 i=1,im
1029 : ! three-cell PP-distribution
1030 : ! Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
1031 : ! a3 = a / 3
1032 : ! b2 = b / 2
1033 0 : s1 = delp(i,1)
1034 0 : s2 = delp(i,2) + s1
1035 : !
1036 0 : s3 = delp(i,2) + delp(i,3)
1037 0 : s4 = s3 + delp(i,4)
1038 0 : ss3 = s3 + s1
1039 0 : s32 = s3*s3
1040 0 : s42 = s4*s4
1041 0 : s34 = s3*s4
1042 : ! model top
1043 0 : a3 = (delq(i,2) - delq(i,1)*s3/s2) / (s3*ss3)
1044 : !
1045 0 : if(abs(a3) .gt. D1EM14) then
1046 0 : b2 = delq(i,1)/s2 - a3*(s1+s2)
1047 0 : sc = -b2/(D3_0*a3)
1048 0 : if(sc .lt. D0_0 .or. sc .gt. s1) then
1049 0 : qe(i,1) = p(i,1) - s1*(a3*s1 + b2)
1050 : else
1051 0 : qe(i,1) = p(i,1) - delq(i,1)*s1/s2
1052 : endif
1053 : else
1054 : ! Linear
1055 0 : qe(i,1) = p(i,1) - delq(i,1)*s1/s2
1056 : endif
1057 0 : dc(i,1) = p(i,1) - qe(i,1)
1058 : ! compute coef. for the off-centered area preserving cubic poly.
1059 0 : dm = delp(i,1) / (s34*ss3*(delp(i,2)+s3)*(s4+delp(i,1)))
1060 0 : f1 = delp(i,2)*s34 / ( s2*ss3*(s4+delp(i,1)) )
1061 : f2 = (delp(i,2)+s3) * (ss3*(delp(i,2)*s3+s34+delp(i,2)*s4) &
1062 0 : + s42*(delp(i,2)+s3+s32/s2))
1063 : f3 = -delp(i,2)*( ss3*(s32*(s3+s4)/(s4-delp(i,2)) &
1064 : + (delp(i,2)*s3+s34+delp(i,2)*s4)) &
1065 0 : + s42*(delp(i,2)+s3) )
1066 0 : f4 = ss3*delp(i,2)*s32*(delp(i,2)+s3) / (s4-delp(i,2))
1067 0 : qe(i,2) = f1*p(i,1)+(f2*p(i,2)+f3*p(i,3)+f4*p(i,4))*dm
1068 0 : 10 continue
1069 :
1070 : ! Bottom
1071 : ! Area preserving cubic with 2nd deriv. = 0 at the surface
1072 0 : do 15 i=1,im
1073 0 : d1 = delp(i,km)
1074 0 : d2 = delp(i,km1)
1075 0 : qm = (d2*p(i,km)+d1*p(i,km1)) / (d1+d2)
1076 0 : dq = D2_0*(p(i,km1)-p(i,km)) / (d1+d2)
1077 0 : c1 = (qe(i,km1)-qm-d2*dq) / (d2*(D2_0*d2*d2+d1*(d2+D3_0*d1)))
1078 0 : c3 = dq - D2_0*c1*(d2*(D5_0*d1+d2)-D3_0*d1**2)
1079 0 : qe(i,km ) = qm - c1*d1*d2*(d2+D3_0*d1)
1080 0 : qe(i,km+1) = d1*(D8_0*c1*d1**2-c3) + qe(i,km)
1081 0 : 15 continue
1082 0 : return
1083 : !EOC
1084 : end subroutine ppme
1085 : !-----------------------------------------------------------------------
1086 :
1087 : !-----------------------------------------------------------------------
1088 : !BOP
1089 : ! !ROUTINE: kmppm --- Perform piecewise parabolic method in vertical
1090 : !
1091 : ! !INTERFACE:
1092 761660928 : subroutine kmppm(dm, a4, itot, lmt)
1093 :
1094 : ! !USES:
1095 : implicit none
1096 :
1097 : ! !INPUT PARAMETERS:
1098 : real(r8), intent(in):: dm(*) ! ??????
1099 : integer, intent(in) :: itot ! Total Longitudes
1100 : integer, intent(in) :: lmt ! 0: Standard PPM constraint
1101 : ! 1: Improved full monotonicity constraint (Lin)
1102 : ! 2: Positive definite constraint
1103 : ! 3: do nothing (return immediately)
1104 :
1105 : ! !INPUT/OUTPUT PARAMETERS:
1106 : real(r8), intent(inout) :: a4(4,*) ! ???????
1107 : ! AA <-- a4(1,i)
1108 : ! AL <-- a4(2,i)
1109 : ! AR <-- a4(3,i)
1110 : ! A6 <-- a4(4,i)
1111 :
1112 : ! !DESCRIPTION:
1113 : !
1114 : ! Writes a standard set of data to the history buffer.
1115 : !
1116 : ! !REVISION HISTORY:
1117 : ! 00.04.24 Lin Last modification
1118 : ! 01.03.26 Sawyer Added ProTeX documentation
1119 : ! 02.04.04 Sawyer Incorporated newest FVGCM version
1120 : !
1121 : !EOP
1122 : !-----------------------------------------------------------------------
1123 : !BOC
1124 : !
1125 : ! !LOCAL VARIABLES:
1126 :
1127 : real(r8) r12
1128 : parameter (r12 = D1_0/D12_0)
1129 :
1130 : real(r8) qmp
1131 : integer i
1132 : real(r8) da1, da2, a6da
1133 : real(r8) fmin
1134 :
1135 : ! Developer: S.-J. Lin, NASA-GSFC
1136 : ! Last modified: Apr 24, 2000
1137 :
1138 761660928 : if ( lmt == 3 ) return
1139 :
1140 761660928 : if(lmt == 0) then
1141 : ! Standard PPM constraint
1142 2380190400 : do i=1,itot
1143 2380190400 : if(dm(i) == D0_0) then
1144 218917873 : a4(2,i) = a4(1,i)
1145 218917873 : a4(3,i) = a4(1,i)
1146 218917873 : a4(4,i) = D0_0
1147 : else
1148 2066064911 : da1 = a4(3,i) - a4(2,i)
1149 2066064911 : da2 = da1**2
1150 2066064911 : a6da = a4(4,i)*da1
1151 2066064911 : if(a6da < -da2) then
1152 451201981 : a4(4,i) = D3_0*(a4(2,i)-a4(1,i))
1153 451201981 : a4(3,i) = a4(2,i) - a4(4,i)
1154 1614862930 : elseif(a6da > da2) then
1155 419957240 : a4(4,i) = D3_0*(a4(3,i)-a4(1,i))
1156 419957240 : a4(2,i) = a4(3,i) - a4(4,i)
1157 : endif
1158 : endif
1159 : enddo
1160 :
1161 666453312 : elseif (lmt == 1) then
1162 :
1163 : ! Improved full monotonicity constraint (Lin)
1164 : ! Note: no need to provide first guess of A6 <-- a4(4,i)
1165 404308800 : do i=1, itot
1166 388136448 : qmp = D2_0*dm(i)
1167 388136448 : a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp)
1168 388136448 : a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp)
1169 404308800 : a4(4,i) = D3_0*( D2_0*a4(1,i) - (a4(2,i)+a4(3,i)) )
1170 : enddo
1171 :
1172 650280960 : elseif (lmt == 2) then
1173 :
1174 : ! Positive definite constraint
1175 16257024000 : do i=1,itot
1176 16257024000 : if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
1177 5491117131 : fmin = a4(1,i)+D0_25*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12
1178 5491117131 : if( fmin < D0_0 ) then
1179 247317870 : if(a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i)) then
1180 23233987 : a4(3,i) = a4(1,i)
1181 23233987 : a4(2,i) = a4(1,i)
1182 23233987 : a4(4,i) = D0_0
1183 224083883 : elseif(a4(3,i) > a4(2,i)) then
1184 155544247 : a4(4,i) = D3_0*(a4(2,i)-a4(1,i))
1185 155544247 : a4(3,i) = a4(2,i) - a4(4,i)
1186 : else
1187 68539636 : a4(4,i) = D3_0*(a4(3,i)-a4(1,i))
1188 68539636 : a4(2,i) = a4(3,i) - a4(4,i)
1189 : endif
1190 : endif
1191 : endif
1192 : enddo
1193 :
1194 : endif
1195 :
1196 : return
1197 : !EOC
1198 : end subroutine kmppm
1199 : !-----------------------------------------------------------------------
1200 :
1201 : !-----------------------------------------------------------------------
1202 : !BOP
1203 : ! !ROUTINE: steepz --- Calculate attributes for PPM
1204 : !
1205 : ! !INTERFACE:
1206 23801904 : subroutine steepz(i1, i2, km, a4, df2, dm, dq, dp, d4)
1207 :
1208 : ! !USES:
1209 : implicit none
1210 :
1211 : ! !INPUT PARAMETERS:
1212 : integer, intent(in) :: km ! Total levels
1213 : integer, intent(in) :: i1 ! Starting longitude
1214 : integer, intent(in) :: i2 ! Finishing longitude
1215 : real(r8), intent(in) :: dp(i1:i2,km) ! grid size
1216 : real(r8), intent(in) :: dq(i1:i2,km) ! backward diff of q
1217 : real(r8), intent(in) :: d4(i1:i2,km) ! backward sum: dp(k)+ dp(k-1)
1218 : real(r8), intent(in) :: df2(i1:i2,km) ! first guess mismatch
1219 : real(r8), intent(in) :: dm(i1:i2,km) ! monotonic mismatch
1220 :
1221 : ! !INPUT/OUTPUT PARAMETERS:
1222 : real(r8), intent(inout) :: a4(4,i1:i2,km) ! first guess/steepened
1223 :
1224 : !
1225 : ! !DESCRIPTION:
1226 : ! This is complicated stuff related to the Piecewise Parabolic Method
1227 : ! and I need to read the Collela/Woodward paper before documenting
1228 : ! thoroughly.
1229 : !
1230 : ! !REVISION HISTORY:
1231 : ! ??.??.?? Lin? Creation
1232 : ! 01.03.26 Sawyer Added ProTeX documentation
1233 : !
1234 : !EOP
1235 : !-----------------------------------------------------------------------
1236 : !BOC
1237 : !
1238 : ! !LOCAL VARIABLES:
1239 : integer i, k
1240 47603808 : real(r8) alfa(i1:i2,km)
1241 47603808 : real(r8) f(i1:i2,km)
1242 47603808 : real(r8) rat(i1:i2,km)
1243 : real(r8) dg2
1244 :
1245 : ! Compute ratio of dq/dp
1246 761660928 : do k=2,km
1247 18470277504 : do i=i1,i2
1248 18446475600 : rat(i,k) = dq(i,k-1) / d4(i,k)
1249 : enddo
1250 : enddo
1251 :
1252 : ! Compute F
1253 737859024 : do k=2,km-1
1254 17875229904 : do i=i1,i2
1255 51412112640 : f(i,k) = (rat(i,k+1) - rat(i,k)) &
1256 69263540640 : / ( dp(i,k-1)+dp(i,k)+dp(i,k+1) )
1257 : enddo
1258 : enddo
1259 :
1260 690255216 : do k=3,km-2
1261 16685134704 : do i=i1,i2
1262 16661332800 : if(f(i,k+1)*f(i,k-1)<D0_0 .and. df2(i,k)/=D0_0) then
1263 : dg2 = (f(i,k+1)-f(i,k-1))*((dp(i,k+1)-dp(i,k-1))**2 &
1264 6924169699 : + d4(i,k)*d4(i,k+1) )
1265 6924169699 : alfa(i,k) = max(D0_0, min(D0_5, -D0_1875*dg2/df2(i,k)))
1266 : else
1267 9070709789 : alfa(i,k) = D0_0
1268 : endif
1269 : enddo
1270 : enddo
1271 :
1272 666453312 : do k=4,km-2
1273 16090087104 : do i=i1,i2
1274 46270901376 : a4(2,i,k) = (D1_0-alfa(i,k-1)-alfa(i,k)) * a4(2,i,k) + &
1275 : alfa(i,k-1)*(a4(1,i,k)-dm(i,k)) + &
1276 62337186576 : alfa(i,k)*(a4(1,i,k-1)+dm(i,k-1))
1277 : enddo
1278 : enddo
1279 :
1280 23801904 : return
1281 : !EOC
1282 : end subroutine steepz
1283 : !-----------------------------------------------------------------------
1284 :
1285 : end module mapz_module
|
