Line data Source code
1 :
2 : !**************************************************************************************
3 : !
4 : ! fv_mapz contains vertical remapping algorithms that come from the FV3 dycore.
5 : ! They have been minimally modified for use in CAM.
6 : !
7 : ! The following license statement is from the original code.
8 : !
9 : !**************************************************************************************
10 :
11 : !***********************************************************************
12 : !* GNU Lesser General Public License
13 : !*
14 : !* This file is part of the FV3 dynamical core.
15 : !*
16 : !* The FV3 dynamical core is free software: you can redistribute it
17 : !* and/or modify it under the terms of the
18 : !* GNU Lesser General Public License as published by the
19 : !* Free Software Foundation, either version 3 of the License, or
20 : !* (at your option) any later version.
21 : !*
22 : !* The FV3 dynamical core is distributed in the hope that it will be
23 : !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
24 : !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
25 : !* See the GNU General Public License for more details.
26 : !*
27 : !* You should have received a copy of the GNU Lesser General Public
28 : !* License along with the FV3 dynamical core.
29 : !* If not, see <http://www.gnu.org/licenses/>.
30 : !***********************************************************************
31 : module fv_mapz
32 :
33 : use shr_kind_mod, only: r8=>shr_kind_r8
34 : use cam_abortutils, only: endrun
35 :
36 : implicit none
37 :
38 : public map_scalar, map1_ppm, mapn_tracer
39 :
40 : real(kind=r8), parameter:: r3 = 1._r8/3._r8, r23 = 2._r8/3._r8, r12 = 1._r8/12._r8
41 : contains
42 :
43 0 : subroutine map_scalar( km, pe1, q1, qs, &
44 0 : kn, pe2, q2, i1, i2, &
45 : j, ibeg, iend, jbeg, jend, iv, kord, q_min)
46 : ! iv=1
47 : integer, intent(in) :: i1 !< Starting longitude
48 : integer, intent(in) :: i2 !< Finishing longitude
49 : integer, intent(in) :: iv !< Mode: 0 == constituents 1 == temp 2 == remap temp with cs scheme
50 : integer, intent(in) :: kord !< Method order
51 : integer, intent(in) :: j !< Current latitude
52 : integer, intent(in) :: ibeg, iend, jbeg, jend
53 : integer, intent(in) :: km !< Original vertical dimension
54 : integer, intent(in) :: kn !< Target vertical dimension
55 : real(kind=r8), intent(in) :: qs(i1:i2) !< bottom BC
56 : real(kind=r8), intent(in) :: pe1(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the original vertical coordinate
57 : real(kind=r8), intent(in) :: pe2(i1:i2,kn+1) !< pressure at layer edges from model top to bottom surface in the new vertical coordinate
58 : real(kind=r8), intent(in) :: q1(ibeg:iend,jbeg:jend,km) !< Field input
59 : ! INPUT/OUTPUT PARAMETERS:
60 : real(kind=r8), intent(inout):: q2(ibeg:iend,jbeg:jend,kn) !< Field output
61 : real(kind=r8), intent(in):: q_min
62 :
63 : ! DESCRIPTION:
64 : ! IV = 0: constituents
65 : ! pe1: pressure at layer edges (from model top to bottom surface)
66 : ! in the original vertical coordinate
67 : ! pe2: pressure at layer edges (from model top to bottom surface)
68 : ! in the new vertical coordinate
69 : ! LOCAL VARIABLES:
70 0 : real(kind=r8) dp1(i1:i2,km)
71 0 : real(kind=r8) q4(4,i1:i2,km)
72 : real(kind=r8) pl, pr, qsum, dp, esl
73 : integer i, k, l, m, k0
74 :
75 0 : do k=1,km
76 0 : do i=i1,i2
77 0 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
78 0 : q4(1,i,k) = q1(i,j,k)
79 : enddo
80 : enddo
81 :
82 : ! Compute vertical subgrid distribution
83 0 : if ( kord >7 ) then
84 0 : call scalar_profile( qs, q4, dp1, km, i1, i2, iv, kord, q_min )
85 : else
86 0 : call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
87 : endif
88 :
89 0 : do i=i1,i2
90 : k0 = 1
91 0 : do 555 k=1,kn
92 0 : do l=k0,km
93 : ! locate the top edge: pe2(i,k)
94 0 : if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
95 0 : pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
96 0 : if( pe2(i,k+1) <= pe1(i,l+1) ) then
97 : ! entire new grid is within the original grid
98 0 : pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
99 0 : q2(i,j,k) = q4(2,i,l) + 0.5_r8*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
100 0 : *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
101 0 : k0 = l
102 0 : goto 555
103 : else
104 : ! Fractional area...
105 : qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5_r8*(q4(4,i,l)+ &
106 : q4(3,i,l)-q4(2,i,l))*(1._r8+pl)-q4(4,i,l)* &
107 0 : (r3*(1._r8+pl*(1._r8+pl))))
108 0 : do m=l+1,km
109 : ! locate the bottom edge: pe2(i,k+1)
110 0 : if( pe2(i,k+1) > pe1(i,m+1) ) then
111 : ! Whole layer
112 0 : qsum = qsum + dp1(i,m)*q4(1,i,m)
113 : else
114 0 : dp = pe2(i,k+1)-pe1(i,m)
115 0 : esl = dp / dp1(i,m)
116 : qsum = qsum + dp*(q4(2,i,m)+0.5_r8*esl* &
117 0 : (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1._r8-r23*esl)))
118 0 : k0 = m
119 0 : goto 123
120 : endif
121 : enddo
122 : goto 123
123 : endif
124 : endif
125 : enddo
126 0 : 123 q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
127 0 : 555 continue
128 : enddo
129 0 : end subroutine map_scalar
130 :
131 :
132 0 : subroutine mapn_tracer(nq, km, pe1, pe2, q1, dp2, kord, j, &
133 : i1, i2, isd, ied, jsd, jed, q_min, fill)
134 : ! INPUT PARAMETERS:
135 : integer, intent(in):: km !< vertical dimension
136 : integer, intent(in):: j, nq, i1, i2
137 : integer, intent(in):: isd, ied, jsd, jed
138 : integer, intent(in):: kord(nq)
139 : real(kind=r8), intent(in):: pe1(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the original vertical coordinate
140 : real(kind=r8), intent(in):: pe2(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the new vertical coordinate
141 : real(kind=r8), intent(in):: dp2(i1:i2,km)
142 : real(kind=r8), intent(in):: q_min
143 : logical, intent(in):: fill
144 : real(kind=r8), intent(inout):: q1(isd:ied,jsd:jed,km,nq) ! Field input
145 : ! LOCAL VARIABLES:
146 0 : real(kind=r8):: q4(4,i1:i2,km,nq)
147 0 : real(kind=r8):: q2(i1:i2,km,nq) !< Field output
148 0 : real(kind=r8):: qsum(nq)
149 0 : real(kind=r8):: dp1(i1:i2,km)
150 0 : real(kind=r8):: qs(i1:i2)
151 : real(kind=r8):: pl, pr, dp, esl, fac1, fac2
152 : integer:: i, k, l, m, k0, iq
153 :
154 0 : do k=1,km
155 0 : do i=i1,i2
156 0 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
157 : enddo
158 : enddo
159 :
160 0 : do iq=1,nq
161 0 : do k=1,km
162 0 : do i=i1,i2
163 0 : q4(1,i,k,iq) = q1(i,j,k,iq)
164 : enddo
165 : enddo
166 0 : call scalar_profile( qs, q4(1,i1,1,iq), dp1, km, i1, i2, 0, kord(iq), q_min )
167 : enddo
168 : ! Mapping
169 0 : do 1000 i=i1,i2
170 : k0 = 1
171 0 : do 555 k=1,km
172 0 : do 100 l=k0,km
173 : ! locate the top edge: pe2(i,k)
174 0 : if(pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1)) then
175 0 : pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
176 0 : if(pe2(i,k+1) <= pe1(i,l+1)) then
177 : ! entire new grid is within the original grid
178 0 : pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
179 0 : fac1 = pr + pl
180 0 : fac2 = r3*(pr*fac1 + pl*pl)
181 0 : fac1 = 0.5_r8*fac1
182 0 : do iq=1,nq
183 0 : q2(i,k,iq) = q4(2,i,l,iq) + (q4(4,i,l,iq)+q4(3,i,l,iq)-q4(2,i,l,iq))*fac1 &
184 0 : - q4(4,i,l,iq)*fac2
185 : enddo
186 : k0 = l
187 : goto 555
188 : else
189 : ! Fractional area...
190 0 : dp = pe1(i,l+1) - pe2(i,k)
191 0 : fac1 = 1.0_r8 + pl
192 0 : fac2 = r3*(1.0_r8+pl*fac1)
193 0 : fac1 = 0.5_r8*fac1
194 0 : do iq=1,nq
195 0 : qsum(iq) = dp*(q4(2,i,l,iq) + (q4(4,i,l,iq)+ &
196 0 : q4(3,i,l,iq) - q4(2,i,l,iq))*fac1 - q4(4,i,l,iq)*fac2)
197 : enddo
198 0 : do m=l+1,km
199 : ! locate the bottom edge: pe2(i,k+1)
200 0 : if(pe2(i,k+1) > pe1(i,m+1) ) then
201 : ! Whole layer..
202 0 : do iq=1,nq
203 0 : qsum(iq) = qsum(iq) + dp1(i,m)*q4(1,i,m,iq)
204 : enddo
205 : else
206 0 : dp = pe2(i,k+1)-pe1(i,m)
207 0 : esl = dp / dp1(i,m)
208 0 : fac1 = 0.5_r8*esl
209 0 : fac2 = 1.0_r8-r23*esl
210 0 : do iq=1,nq
211 0 : qsum(iq) = qsum(iq) + dp*( q4(2,i,m,iq) + fac1*( &
212 0 : q4(3,i,m,iq)-q4(2,i,m,iq)+q4(4,i,m,iq)*fac2 ) )
213 : enddo
214 : k0 = m
215 : goto 123
216 : endif
217 : enddo
218 : goto 123
219 : endif
220 : endif
221 0 : 100 continue
222 : 123 continue
223 0 : do iq=1,nq
224 0 : q2(i,k,iq) = qsum(iq) / dp2(i,k)
225 : enddo
226 0 : 555 continue
227 0 : 1000 continue
228 :
229 0 : if (fill) call fillz(i2-i1+1, km, nq, q2, dp2)
230 :
231 0 : do iq=1,nq
232 : ! if (fill) call fillz(i2-i1+1, km, 1, q2(i1,1,iq), dp2)
233 0 : do k=1,km
234 0 : do i=i1,i2
235 0 : q1(i,j,k,iq) = q2(i,k,iq)
236 : enddo
237 : enddo
238 : enddo
239 :
240 0 : end subroutine mapn_tracer
241 :
242 :
243 171428400 : subroutine map1_ppm( km, pe1, q1, qs, &
244 171428400 : kn, pe2, q2, i1, i2, &
245 : j, ibeg, iend, jbeg, jend, iv, kord)
246 : integer, intent(in) :: i1 !< Starting longitude
247 : integer, intent(in) :: i2 !< Finishing longitude
248 : integer, intent(in) :: iv !< Mode: 0 == constituents 1 == ??? 2 == remap temp with cs scheme
249 : integer, intent(in) :: kord !< Method order
250 : integer, intent(in) :: j !< Current latitude
251 : integer, intent(in) :: ibeg, iend, jbeg, jend
252 : integer, intent(in) :: km !< Original vertical dimension
253 : integer, intent(in) :: kn !< Target vertical dimension
254 : real(kind=r8), intent(in) :: qs(i1:i2) !< bottom BC
255 : real(kind=r8), intent(in) :: pe1(i1:i2,km+1) !< pressure at layer edges from model top to bottom surface in the original vertical coordinate
256 : real(kind=r8), intent(in) :: pe2(i1:i2,kn+1) !< pressure at layer edges from model top to bottom surface in the new vertical coordinate
257 : real(kind=r8), intent(in) :: q1(ibeg:iend,jbeg:jend,km) !< Field input
258 : ! INPUT/OUTPUT PARAMETERS:
259 : real(kind=r8), intent(inout):: q2(ibeg:iend,jbeg:jend,kn) !< Field output
260 :
261 : ! DESCRIPTION:
262 : ! IV = 0: constituents
263 : ! pe1: pressure at layer edges (from model top to bottom surface)
264 : ! in the original vertical coordinate
265 : ! pe2: pressure at layer edges (from model top to bottom surface)
266 : ! in the new vertical coordinate
267 :
268 : ! LOCAL VARIABLES:
269 342856800 : real(kind=r8) dp1(i1:i2,km)
270 342856800 : real(kind=r8) q4(4,i1:i2,km)
271 : real(kind=r8) pl, pr, qsum, dp, esl
272 : integer i, k, l, m, k0
273 :
274 4628566800 : do k=1,km
275 21241537200 : do i=i1,i2
276 16612970400 : dp1(i,k) = pe1(i,k+1) - pe1(i,k)
277 21070108800 : q4(1,i,k) = q1(i,j,k)
278 : enddo
279 : enddo
280 :
281 : ! Compute vertical subgrid distribution
282 171428400 : if ( kord >7 ) then
283 171428400 : call cs_profile( qs, q4, dp1, km, i1, i2, iv, kord )
284 : else
285 0 : call ppm_profile( q4, dp1, km, i1, i2, iv, kord )
286 : endif
287 :
288 810388800 : do i=i1,i2
289 : k0 = 1
290 17423359200 : do 555 k=1,kn
291 16612970400 : do l=k0,km
292 : ! locate the top edge: pe2(i,k)
293 16612970400 : if( pe2(i,k) >= pe1(i,l) .and. pe2(i,k) <= pe1(i,l+1) ) then
294 16612970400 : pl = (pe2(i,k)-pe1(i,l)) / dp1(i,l)
295 16612970400 : if( pe2(i,k+1) <= pe1(i,l+1) ) then
296 : ! entire new grid is within the original grid
297 2383214901 : pr = (pe2(i,k+1)-pe1(i,l)) / dp1(i,l)
298 4766429802 : q2(i,j,k) = q4(2,i,l) + 0.5_r8*(q4(4,i,l)+q4(3,i,l)-q4(2,i,l)) &
299 2383214901 : *(pr+pl)-q4(4,i,l)*r3*(pr*(pr+pl)+pl**2)
300 2383214901 : k0 = l
301 2383214901 : goto 555
302 : else
303 : ! Fractional area...
304 : qsum = (pe1(i,l+1)-pe2(i,k))*(q4(2,i,l)+0.5_r8*(q4(4,i,l)+ &
305 : q4(3,i,l)-q4(2,i,l))*(1.0_r8+pl)-q4(4,i,l)* &
306 14229755499 : (r3*(1.0_r8+pl*(1.0_r8+pl))))
307 15974010000 : do m=l+1,km
308 : ! locate the bottom edge: pe2(i,k+1)
309 15974010000 : if( pe2(i,k+1) > pe1(i,m+1) ) then
310 : ! Whole layer
311 1744254501 : qsum = qsum + dp1(i,m)*q4(1,i,m)
312 : else
313 14229755499 : dp = pe2(i,k+1)-pe1(i,m)
314 14229755499 : esl = dp / dp1(i,m)
315 : qsum = qsum + dp*(q4(2,i,m)+0.5_r8*esl* &
316 14229755499 : (q4(3,i,m)-q4(2,i,m)+q4(4,i,m)*(1.0_r8-r23*esl)))
317 14229755499 : k0 = m
318 14229755499 : goto 123
319 : endif
320 : enddo
321 : goto 123
322 : endif
323 : endif
324 : enddo
325 14229755499 : 123 q2(i,j,k) = qsum / ( pe2(i,k+1) - pe2(i,k) )
326 638960400 : 555 continue
327 : enddo
328 :
329 171428400 : end subroutine map1_ppm
330 :
331 0 : subroutine ppm_profile(a4, delp, km, i1, i2, iv, kord)
332 :
333 : ! INPUT PARAMETERS:
334 : integer, intent(in):: iv !< iv =-1: winds iv = 0: positive definite scalars iv = 1: others iv = 2: temp (if remap_t) and w (iv=-2)
335 : integer, intent(in):: i1 !< Starting longitude
336 : integer, intent(in):: i2 !< Finishing longitude
337 : integer, intent(in):: km !< Vertical dimension
338 : integer, intent(in):: kord !< Order (or more accurately method no.):
339 : !
340 : real(kind=r8) , intent(in):: delp(i1:i2,km) !< Layer pressure thickness
341 :
342 : ! !INPUT/OUTPUT PARAMETERS:
343 : real(kind=r8) , intent(inout):: a4(4,i1:i2,km) !< Interpolated values
344 :
345 : ! DESCRIPTION:
346 : !
347 : ! Perform the piecewise parabolic reconstruction
348 : !
349 : ! !REVISION HISTORY:
350 : ! S.-J. Lin revised at GFDL 2007
351 : !-----------------------------------------------------------------------
352 : ! local arrays:
353 0 : real(kind=r8) dc(i1:i2,km)
354 0 : real(kind=r8) h2(i1:i2,km)
355 0 : real(kind=r8) delq(i1:i2,km)
356 0 : real(kind=r8) df2(i1:i2,km)
357 0 : real(kind=r8) d4(i1:i2,km)
358 :
359 : ! local scalars:
360 : integer i, k, km1, lmt, it
361 : real(kind=r8) fac
362 : real(kind=r8) a1, a2, c1, c2, c3, d1, d2
363 : real(kind=r8) qm, dq, lac, qmp, pmp
364 :
365 0 : km1 = km - 1
366 0 : it = i2 - i1 + 1
367 :
368 0 : do k=2,km
369 0 : do i=i1,i2
370 0 : delq(i,k-1) = a4(1,i,k) - a4(1,i,k-1)
371 0 : d4(i,k ) = delp(i,k-1) + delp(i,k)
372 : enddo
373 : enddo
374 :
375 0 : do k=2,km1
376 0 : do i=i1,i2
377 0 : c1 = (delp(i,k-1)+0.5_r8*delp(i,k))/d4(i,k+1)
378 0 : c2 = (delp(i,k+1)+0.5_r8*delp(i,k))/d4(i,k)
379 : df2(i,k) = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / &
380 0 : (d4(i,k)+delp(i,k+1))
381 : dc(i,k) = sign( min(abs(df2(i,k)), &
382 : max(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))-a4(1,i,k), &
383 0 : a4(1,i,k)-min(a4(1,i,k-1),a4(1,i,k),a4(1,i,k+1))), df2(i,k) )
384 : enddo
385 : enddo
386 :
387 : !-----------------------------------------------------------
388 : ! 4th order interpolation of the provisional cell edge value
389 : !-----------------------------------------------------------
390 :
391 0 : do k=3,km1
392 0 : do i=i1,i2
393 0 : c1 = delq(i,k-1)*delp(i,k-1) / d4(i,k)
394 0 : a1 = d4(i,k-1) / (d4(i,k) + delp(i,k-1))
395 0 : a2 = d4(i,k+1) / (d4(i,k) + delp(i,k))
396 : a4(2,i,k) = a4(1,i,k-1) + c1 + 2.0_r8/(d4(i,k-1)+d4(i,k+1)) * &
397 : ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - &
398 0 : delp(i,k-1)*a1*dc(i,k ) )
399 : enddo
400 : enddo
401 :
402 : ! if(km>8 .and. kord>4) call steepz(i1, i2, km, a4, df2, dc, delq, delp, d4)
403 :
404 : ! Area preserving cubic with 2nd deriv. = 0 at the boundaries
405 : ! Top
406 0 : do i=i1,i2
407 0 : d1 = delp(i,1)
408 0 : d2 = delp(i,2)
409 0 : qm = (d2*a4(1,i,1)+d1*a4(1,i,2)) / (d1+d2)
410 0 : dq = 2.0_r8*(a4(1,i,2)-a4(1,i,1)) / (d1+d2)
411 0 : c1 = 4.0_r8*(a4(2,i,3)-qm-d2*dq) / ( d2*(2.0_r8*d2*d2+d1*(d2+3.0_r8*d1)) )
412 0 : c3 = dq - 0.5_r8*c1*(d2*(5.0_r8*d1+d2)-3.0_r8*d1*d1)
413 0 : a4(2,i,2) = qm - 0.25_r8*c1*d1*d2*(d2+3.0_r8*d1)
414 : ! Top edge:
415 : !-------------------------------------------------------
416 0 : a4(2,i,1) = d1*(2.0_r8*c1*d1**2-c3) + a4(2,i,2)
417 : !-------------------------------------------------------
418 : ! a4(2,i,1) = (12./7.)*a4(1,i,1)-(13./14.)*a4(1,i,2)+(3./14.)*a4(1,i,3)
419 : !-------------------------------------------------------
420 : ! No over- and undershoot condition
421 0 : a4(2,i,2) = max( a4(2,i,2), min(a4(1,i,1), a4(1,i,2)) )
422 0 : a4(2,i,2) = min( a4(2,i,2), max(a4(1,i,1), a4(1,i,2)) )
423 0 : dc(i,1) = 0.5_r8*(a4(2,i,2) - a4(1,i,1))
424 : enddo
425 :
426 : ! Enforce monotonicity within the top layer
427 :
428 0 : if( iv==0 ) then
429 0 : do i=i1,i2
430 0 : a4(2,i,1) = max(0.0_r8, a4(2,i,1))
431 0 : a4(2,i,2) = max(0.0_r8, a4(2,i,2))
432 : enddo
433 0 : elseif( iv==-1 ) then
434 0 : do i=i1,i2
435 0 : if ( a4(2,i,1)*a4(1,i,1) <= 0.0_r8 ) a4(2,i,1) = 0.0_r8
436 : enddo
437 0 : elseif( abs(iv)==2 ) then
438 0 : do i=i1,i2
439 0 : a4(2,i,1) = a4(1,i,1)
440 0 : a4(3,i,1) = a4(1,i,1)
441 : enddo
442 : endif
443 :
444 : ! Bottom
445 : ! Area preserving cubic with 2nd deriv. = 0 at the surface
446 0 : do i=i1,i2
447 0 : d1 = delp(i,km)
448 0 : d2 = delp(i,km1)
449 0 : qm = (d2*a4(1,i,km)+d1*a4(1,i,km1)) / (d1+d2)
450 0 : dq = 2.0_r8*(a4(1,i,km1)-a4(1,i,km)) / (d1+d2)
451 0 : c1 = (a4(2,i,km1)-qm-d2*dq) / (d2*(2.0_r8*d2*d2+d1*(d2+3.0_r8*d1)))
452 0 : c3 = dq - 2.0_r8*c1*(d2*(5._r8*d1+d2)-3.0_r8*d1*d1)
453 0 : a4(2,i,km) = qm - c1*d1*d2*(d2+3.0_r8*d1)
454 : ! Bottom edge:
455 : !-----------------------------------------------------
456 0 : a4(3,i,km) = d1*(8.0_r8*c1*d1**2-c3) + a4(2,i,km)
457 : ! dc(i,km) = 0.5*(a4(3,i,km) - a4(1,i,km))
458 : !-----------------------------------------------------
459 : ! a4(3,i,km) = (12./7.)*a4(1,i,km)-(13./14.)*a4(1,i,km-1)+(3./14.)*a4(1,i,km-2)
460 : ! No over- and under-shoot condition
461 0 : a4(2,i,km) = max( a4(2,i,km), min(a4(1,i,km), a4(1,i,km1)) )
462 0 : a4(2,i,km) = min( a4(2,i,km), max(a4(1,i,km), a4(1,i,km1)) )
463 0 : dc(i,km) = 0.5_r8*(a4(1,i,km) - a4(2,i,km))
464 : enddo
465 :
466 :
467 : ! Enforce constraint on the "slope" at the surface
468 :
469 : #ifdef BOT_MONO
470 : do i=i1,i2
471 : a4(4,i,km) = 0
472 : if( a4(3,i,km) * a4(1,i,km) <= 0.0_r8 ) a4(3,i,km) = 0.0_r8
473 : d1 = a4(1,i,km) - a4(2,i,km)
474 : d2 = a4(3,i,km) - a4(1,i,km)
475 : if ( d1*d2 < 0.0_r8 ) then
476 : a4(2,i,km) = a4(1,i,km)
477 : a4(3,i,km) = a4(1,i,km)
478 : else
479 : dq = sign(min(abs(d1),abs(d2),0.5_r8*abs(delq(i,km-1))), d1)
480 : a4(2,i,km) = a4(1,i,km) - dq
481 : a4(3,i,km) = a4(1,i,km) + dq
482 : endif
483 : enddo
484 : #else
485 0 : if( iv==0 ) then
486 0 : do i=i1,i2
487 0 : a4(2,i,km) = max(0.0_r8,a4(2,i,km))
488 0 : a4(3,i,km) = max(0.0_r8,a4(3,i,km))
489 : enddo
490 0 : elseif( iv<0 ) then
491 0 : do i=i1,i2
492 0 : if( a4(1,i,km)*a4(3,i,km) <= 0.0_r8 ) a4(3,i,km) = 0.0_r8
493 : enddo
494 : endif
495 : #endif
496 :
497 0 : do k=1,km1
498 0 : do i=i1,i2
499 0 : a4(3,i,k) = a4(2,i,k+1)
500 : enddo
501 : enddo
502 :
503 : !-----------------------------------------------------------
504 : ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
505 : !-----------------------------------------------------------
506 : ! Top 2 and bottom 2 layers always use monotonic mapping
507 0 : do k=1,2
508 0 : do i=i1,i2
509 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
510 : enddo
511 0 : call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
512 : enddo
513 :
514 0 : if(kord >= 7) then
515 : !-----------------------
516 : ! Huynh's 2nd constraint
517 : !-----------------------
518 0 : do k=2,km1
519 0 : do i=i1,i2
520 : ! Method#1
521 : ! h2(i,k) = delq(i,k) - delq(i,k-1)
522 : ! Method#2 - better
523 0 : h2(i,k) = 2.0_r8*(dc(i,k+1)/delp(i,k+1) - dc(i,k-1)/delp(i,k-1)) &
524 : / ( delp(i,k)+0.5_r8*(delp(i,k-1)+delp(i,k+1)) ) &
525 0 : * delp(i,k)**2
526 : ! Method#3
527 : !!! h2(i,k) = dc(i,k+1) - dc(i,k-1)
528 : enddo
529 : enddo
530 :
531 0 : fac = 1.5_r8 ! original quasi-monotone
532 :
533 0 : do k=3,km-2
534 0 : do i=i1,i2
535 : ! Right edges
536 : ! qmp = a4(1,i,k) + 2.0*delq(i,k-1)
537 : ! lac = a4(1,i,k) + fac*h2(i,k-1) + 0.5*delq(i,k-1)
538 : !
539 0 : pmp = 2.0_r8*dc(i,k)
540 0 : qmp = a4(1,i,k) + pmp
541 0 : lac = a4(1,i,k) + fac*h2(i,k-1) + dc(i,k)
542 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), qmp, lac)), &
543 0 : max(a4(1,i,k), qmp, lac) )
544 : ! Left edges
545 : ! qmp = a4(1,i,k) - 2.0*delq(i,k)
546 : ! lac = a4(1,i,k) + fac*h2(i,k+1) - 0.5*delq(i,k)
547 : !
548 0 : qmp = a4(1,i,k) - pmp
549 0 : lac = a4(1,i,k) + fac*h2(i,k+1) - dc(i,k)
550 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), qmp, lac)), &
551 0 : max(a4(1,i,k), qmp, lac))
552 : !-------------
553 : ! Recompute A6
554 : !-------------
555 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
556 : enddo
557 : ! Additional constraint to ensure positivity when kord=7
558 0 : if (iv == 0 .and. kord >= 6 ) &
559 0 : call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 2)
560 : enddo
561 :
562 : else
563 :
564 0 : lmt = kord - 3
565 0 : lmt = max(0, lmt)
566 0 : if (iv == 0) lmt = min(2, lmt)
567 :
568 0 : do k=3,km-2
569 0 : if( kord /= 4) then
570 0 : do i=i1,i2
571 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
572 : enddo
573 : endif
574 0 : if(kord/=6) call ppm_limiters(dc(i1,k), a4(1,i1,k), it, lmt)
575 : enddo
576 : endif
577 :
578 0 : do k=km1,km
579 0 : do i=i1,i2
580 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
581 : enddo
582 0 : call ppm_limiters(dc(i1,k), a4(1,i1,k), it, 0)
583 : enddo
584 :
585 0 : end subroutine ppm_profile
586 :
587 0 : subroutine ppm_limiters(dm, a4, itot, lmt)
588 :
589 : ! INPUT PARAMETERS:
590 : real(kind=r8) , intent(in):: dm(*) !< Linear slope
591 : integer, intent(in) :: itot !< Total Longitudes
592 : integer, intent(in) :: lmt !< 0: Standard PPM constraint 1: Improved full monotonicity constraint
593 : !< (Lin) 2: Positive definite constraint
594 : !< 3: do nothing (return immediately)
595 : ! INPUT/OUTPUT PARAMETERS:
596 : real(kind=r8) , intent(inout) :: a4(4,*) !< PPM array AA <-- a4(1,i) AL <-- a4(2,i) AR <-- a4(3,i) A6 <-- a4(4,i)
597 : ! LOCAL VARIABLES:
598 : real(kind=r8) qmp
599 : real(kind=r8) da1, da2, a6da
600 : real(kind=r8) fmin
601 : integer i
602 :
603 : ! Developer: S.-J. Lin
604 :
605 0 : if ( lmt == 3 ) return
606 :
607 0 : if(lmt == 0) then
608 : ! Standard PPM constraint
609 0 : do i=1,itot
610 0 : if(dm(i) == 0.0_r8) then
611 0 : a4(2,i) = a4(1,i)
612 0 : a4(3,i) = a4(1,i)
613 0 : a4(4,i) = 0.0_r8
614 : else
615 0 : da1 = a4(3,i) - a4(2,i)
616 0 : da2 = da1**2
617 0 : a6da = a4(4,i)*da1
618 0 : if(a6da < -da2) then
619 0 : a4(4,i) = 3.0_r8*(a4(2,i)-a4(1,i))
620 0 : a4(3,i) = a4(2,i) - a4(4,i)
621 0 : elseif(a6da > da2) then
622 0 : a4(4,i) = 3.0_r8*(a4(3,i)-a4(1,i))
623 0 : a4(2,i) = a4(3,i) - a4(4,i)
624 : endif
625 : endif
626 : enddo
627 :
628 0 : elseif (lmt == 1) then
629 :
630 : ! Improved full monotonicity constraint (Lin 2004)
631 : ! Note: no need to provide first guess of A6 <-- a4(4,i)
632 0 : do i=1, itot
633 0 : qmp = 2.0_r8*dm(i)
634 0 : a4(2,i) = a4(1,i)-sign(min(abs(qmp),abs(a4(2,i)-a4(1,i))), qmp)
635 0 : a4(3,i) = a4(1,i)+sign(min(abs(qmp),abs(a4(3,i)-a4(1,i))), qmp)
636 0 : a4(4,i) = 3.0_r8*( 2.0_r8*a4(1,i) - (a4(2,i)+a4(3,i)) )
637 : enddo
638 :
639 0 : elseif (lmt == 2) then
640 :
641 : ! Positive definite constraint
642 0 : do i=1,itot
643 0 : if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
644 0 : fmin = a4(1,i)+0.25_r8*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12
645 0 : if( fmin < 0.0_r8 ) then
646 0 : if(a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i)) then
647 0 : a4(3,i) = a4(1,i)
648 0 : a4(2,i) = a4(1,i)
649 0 : a4(4,i) = 0.0_r8
650 0 : elseif(a4(3,i) > a4(2,i)) then
651 0 : a4(4,i) = 3.0_r8*(a4(2,i)-a4(1,i))
652 0 : a4(3,i) = a4(2,i) - a4(4,i)
653 : else
654 0 : a4(4,i) = 3.0_r8*(a4(3,i)-a4(1,i))
655 0 : a4(2,i) = a4(3,i) - a4(4,i)
656 : endif
657 : endif
658 : endif
659 : enddo
660 :
661 : endif
662 :
663 : end subroutine ppm_limiters
664 :
665 :
666 0 : subroutine scalar_profile(qs, a4, delp, km, i1, i2, iv, kord, qmin)
667 : ! Optimized vertical profile reconstruction:
668 : ! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL
669 : integer, intent(in):: i1, i2
670 : integer, intent(in):: km !< vertical dimension
671 : integer, intent(in):: iv !< iv =-1: winds iv = 0: positive definite scalars iv = 1: others
672 : integer, intent(in):: kord
673 : real(kind=r8), intent(in) :: qs(i1:i2)
674 : real(kind=r8), intent(in) :: delp(i1:i2,km) !< Layer pressure thickness
675 : real(kind=r8), intent(inout):: a4(4,i1:i2,km) !< Interpolated values
676 : real(kind=r8), intent(in):: qmin
677 : !-----------------------------------------------------------------------
678 0 : logical, dimension(i1:i2,km):: extm, ext5, ext6
679 0 : real(kind=r8) gam(i1:i2,km)
680 0 : real(kind=r8) q(i1:i2,km+1)
681 0 : real(kind=r8) d4(i1:i2)
682 : real(kind=r8) bet, a_bot, grat
683 : real(kind=r8) pmp_1, lac_1, pmp_2, lac_2, x0, x1
684 : integer i, k, im
685 :
686 0 : if ( iv .eq. -2 ) then
687 0 : do i=i1,i2
688 0 : gam(i,2) = 0.5_r8
689 0 : q(i,1) = 1.5_r8*a4(1,i,1)
690 : enddo
691 0 : do k=2,km-1
692 0 : do i=i1, i2
693 0 : grat = delp(i,k-1) / delp(i,k)
694 0 : bet = 2.0_r8 + grat + grat - gam(i,k)
695 0 : q(i,k) = (3.0_r8*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
696 0 : gam(i,k+1) = grat / bet
697 : enddo
698 : enddo
699 0 : do i=i1,i2
700 0 : grat = delp(i,km-1) / delp(i,km)
701 0 : q(i,km) = (3.0_r8*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / &
702 0 : (2.0_r8 + grat + grat - gam(i,km))
703 0 : q(i,km+1) = qs(i)
704 : enddo
705 0 : do k=km-1,1,-1
706 0 : do i=i1,i2
707 0 : q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
708 : enddo
709 : enddo
710 : else
711 0 : do i=i1,i2
712 0 : grat = delp(i,2) / delp(i,1) ! grid ratio
713 0 : bet = grat*(grat+0.5_r8)
714 0 : q(i,1) = ( (grat+grat)*(grat+1.0_r8)*a4(1,i,1) + a4(1,i,2) ) / bet
715 0 : gam(i,1) = ( 1.0_r8 + grat*(grat+1.5_r8) ) / bet
716 : enddo
717 :
718 0 : do k=2,km
719 0 : do i=i1,i2
720 0 : d4(i) = delp(i,k-1) / delp(i,k)
721 0 : bet = 2.0_r8 + d4(i) + d4(i) - gam(i,k-1)
722 0 : q(i,k) = ( 3.0_r8*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
723 0 : gam(i,k) = d4(i) / bet
724 : enddo
725 : enddo
726 :
727 0 : do i=i1,i2
728 0 : a_bot = 1.0_r8 + d4(i)*(d4(i)+1.5_r8)
729 0 : q(i,km+1) = (2.0_r8*d4(i)*(d4(i)+1.0_r8)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
730 0 : / ( d4(i)*(d4(i)+0.5_r8) - a_bot*gam(i,km) )
731 : enddo
732 :
733 0 : do k=km,1,-1
734 0 : do i=i1,i2
735 0 : q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
736 : enddo
737 : enddo
738 : endif
739 :
740 : !----- Perfectly linear scheme --------------------------------
741 0 : if ( abs(kord) > 16 ) then
742 0 : do k=1,km
743 0 : do i=i1,i2
744 0 : a4(2,i,k) = q(i,k )
745 0 : a4(3,i,k) = q(i,k+1)
746 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
747 : enddo
748 : enddo
749 : return
750 : endif
751 : !----- Perfectly linear scheme --------------------------------
752 : !------------------
753 : ! Apply constraints
754 : !------------------
755 0 : im = i2 - i1 + 1
756 :
757 : ! Apply *large-scale* constraints
758 0 : do i=i1,i2
759 0 : q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) )
760 0 : q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) )
761 : enddo
762 :
763 0 : do k=2,km
764 0 : do i=i1,i2
765 0 : gam(i,k) = a4(1,i,k) - a4(1,i,k-1)
766 : enddo
767 : enddo
768 :
769 : ! Interior:
770 0 : do k=3,km-1
771 0 : do i=i1,i2
772 0 : if ( gam(i,k-1)*gam(i,k+1)>0.0_r8 ) then
773 : ! Apply large-scale constraint to ALL fields if not local max/min
774 0 : q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) )
775 0 : q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) )
776 : else
777 0 : if ( gam(i,k-1) > 0.0_r8 ) then
778 : ! There exists a local max
779 0 : q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k)))
780 : else
781 : ! There exists a local min
782 0 : q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k)))
783 0 : if ( iv==0 ) q(i,k) = max(0.0_r8, q(i,k))
784 : endif
785 : endif
786 : enddo
787 : enddo
788 :
789 : ! Bottom:
790 0 : do i=i1,i2
791 0 : q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) )
792 0 : q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) )
793 : enddo
794 :
795 0 : do k=1,km
796 0 : do i=i1,i2
797 0 : a4(2,i,k) = q(i,k )
798 0 : a4(3,i,k) = q(i,k+1)
799 : enddo
800 : enddo
801 :
802 0 : do k=1,km
803 0 : if ( k==1 .or. k==km ) then
804 0 : do i=i1,i2
805 0 : extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0.0_r8
806 : enddo
807 : else
808 0 : do i=i1,i2
809 0 : extm(i,k) = gam(i,k)*gam(i,k+1) < 0.0_r8
810 : enddo
811 : endif
812 0 : if ( abs(kord) > 9 ) then
813 0 : do i=i1,i2
814 0 : x0 = 2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))
815 0 : x1 = abs(a4(2,i,k)-a4(3,i,k))
816 0 : a4(4,i,k) = 3.0_r8*x0
817 0 : ext5(i,k) = abs(x0) > x1
818 0 : ext6(i,k) = abs(a4(4,i,k)) > x1
819 : enddo
820 : endif
821 : enddo
822 :
823 : !---------------------------
824 : ! Apply subgrid constraints:
825 : !---------------------------
826 : ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
827 : ! Top 2 and bottom 2 layers always use monotonic mapping
828 :
829 0 : if ( iv==0 ) then
830 0 : do i=i1,i2
831 0 : a4(2,i,1) = max(0.0_r8, a4(2,i,1))
832 : enddo
833 0 : elseif ( iv==-1 ) then
834 0 : do i=i1,i2
835 0 : if ( a4(2,i,1)*a4(1,i,1) <= 0.0_r8 ) a4(2,i,1) = 0.0_r8
836 : enddo
837 0 : elseif ( iv==2 ) then
838 0 : do i=i1,i2
839 0 : a4(2,i,1) = a4(1,i,1)
840 0 : a4(3,i,1) = a4(1,i,1)
841 0 : a4(4,i,1) = 0.0_r8
842 : enddo
843 : endif
844 :
845 0 : if ( iv/=2 ) then
846 0 : do i=i1,i2
847 0 : a4(4,i,1) = 3.0_r8*(2.0_r8*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1)))
848 : enddo
849 0 : call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1)
850 : endif
851 :
852 : ! k=2
853 0 : do i=i1,i2
854 0 : a4(4,i,2) = 3.0_r8*(2.0_r8*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2)))
855 : enddo
856 0 : call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2)
857 :
858 : !-------------------------------------
859 : ! Huynh's 2nd constraint for interior:
860 : !-------------------------------------
861 0 : do k=3,km-2
862 0 : if ( abs(kord)<9 ) then
863 0 : do i=i1,i2
864 : ! Left edges
865 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
866 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
867 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
868 0 : max(a4(1,i,k), pmp_1, lac_1) )
869 : ! Right edges
870 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
871 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
872 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
873 0 : max(a4(1,i,k), pmp_2, lac_2) )
874 :
875 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
876 : enddo
877 :
878 : elseif ( abs(kord)==9 ) then
879 0 : do i=i1,i2
880 0 : if ( extm(i,k) .and. extm(i,k-1) ) then
881 : ! grid-scale 2-delta-z wave detected
882 0 : a4(2,i,k) = a4(1,i,k)
883 0 : a4(3,i,k) = a4(1,i,k)
884 0 : a4(4,i,k) = 0.0_r8
885 0 : else if ( extm(i,k) .and. extm(i,k+1) ) then
886 : ! grid-scale 2-delta-z wave detected
887 0 : a4(2,i,k) = a4(1,i,k)
888 0 : a4(3,i,k) = a4(1,i,k)
889 0 : a4(4,i,k) = 0.0_r8
890 0 : else if ( extm(i,k) .and. a4(1,i,k)<qmin ) then
891 : ! grid-scale 2-delta-z wave detected
892 0 : a4(2,i,k) = a4(1,i,k)
893 0 : a4(3,i,k) = a4(1,i,k)
894 0 : a4(4,i,k) = 0.0_r8
895 : else
896 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
897 : ! Check within the smooth region if subgrid profile is non-monotonic
898 0 : if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
899 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
900 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
901 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
902 0 : max(a4(1,i,k), pmp_1, lac_1) )
903 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
904 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
905 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
906 0 : max(a4(1,i,k), pmp_2, lac_2) )
907 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
908 : endif
909 : endif
910 : enddo
911 : elseif ( abs(kord)==10 ) then
912 0 : do i=i1,i2
913 0 : if( ext5(i,k) ) then
914 0 : if( ext5(i,k-1) .or. ext5(i,k+1) ) then
915 0 : a4(2,i,k) = a4(1,i,k)
916 0 : a4(3,i,k) = a4(1,i,k)
917 0 : elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
918 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
919 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
920 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
921 0 : max(a4(1,i,k), pmp_1, lac_1) )
922 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
923 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
924 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
925 0 : max(a4(1,i,k), pmp_2, lac_2) )
926 : endif
927 0 : elseif( ext6(i,k) ) then
928 0 : if( ext5(i,k-1) .or. ext5(i,k+1) ) then
929 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
930 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
931 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
932 0 : max(a4(1,i,k), pmp_1, lac_1) )
933 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
934 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
935 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
936 0 : max(a4(1,i,k), pmp_2, lac_2) )
937 : endif
938 : endif
939 : enddo
940 0 : do i=i1,i2
941 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
942 : enddo
943 : elseif ( abs(kord)==12 ) then
944 0 : do i=i1,i2
945 0 : if( extm(i,k) ) then
946 0 : a4(2,i,k) = a4(1,i,k)
947 0 : a4(3,i,k) = a4(1,i,k)
948 0 : a4(4,i,k) = 0.0_r8
949 : else ! not a local extremum
950 0 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
951 : ! Check within the smooth region if subgrid profile is non-monotonic
952 0 : if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
953 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
954 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
955 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
956 0 : max(a4(1,i,k), pmp_1, lac_1) )
957 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
958 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
959 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
960 0 : max(a4(1,i,k), pmp_2, lac_2) )
961 0 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
962 : endif
963 : endif
964 : enddo
965 : elseif ( abs(kord)==13 ) then
966 0 : do i=i1,i2
967 0 : if( ext6(i,k) ) then
968 0 : if ( ext6(i,k-1) .and. ext6(i,k+1) ) then
969 : ! grid-scale 2-delta-z wave detected
970 0 : a4(2,i,k) = a4(1,i,k)
971 0 : a4(3,i,k) = a4(1,i,k)
972 : endif
973 : endif
974 : enddo
975 0 : do i=i1,i2
976 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
977 : enddo
978 : elseif ( abs(kord)==14 ) then
979 :
980 0 : do i=i1,i2
981 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
982 : enddo
983 :
984 : elseif ( abs(kord)==15 ) then ! Revised abs(kord)=9 scheme
985 0 : do i=i1,i2
986 0 : if ( ext5(i,k) .and. ext5(i,k-1) ) then
987 0 : a4(2,i,k) = a4(1,i,k)
988 0 : a4(3,i,k) = a4(1,i,k)
989 0 : else if ( ext5(i,k) .and. ext5(i,k+1) ) then
990 0 : a4(2,i,k) = a4(1,i,k)
991 0 : a4(3,i,k) = a4(1,i,k)
992 0 : else if ( ext5(i,k) .and. a4(1,i,k)<qmin ) then
993 0 : a4(2,i,k) = a4(1,i,k)
994 0 : a4(3,i,k) = a4(1,i,k)
995 0 : elseif( ext6(i,k) ) then
996 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
997 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
998 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
999 0 : max(a4(1,i,k), pmp_1, lac_1) )
1000 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1001 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1002 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1003 0 : max(a4(1,i,k), pmp_2, lac_2) )
1004 : endif
1005 : enddo
1006 0 : do i=i1,i2
1007 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1008 : enddo
1009 : elseif ( abs(kord)==16 ) then
1010 0 : do i=i1,i2
1011 0 : if( ext5(i,k) ) then
1012 0 : if ( ext5(i,k-1) .or. ext5(i,k+1) ) then
1013 0 : a4(2,i,k) = a4(1,i,k)
1014 0 : a4(3,i,k) = a4(1,i,k)
1015 0 : elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
1016 : ! Left edges
1017 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1018 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1019 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1020 0 : max(a4(1,i,k), pmp_1, lac_1) )
1021 : ! Right edges
1022 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1023 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1024 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1025 0 : max(a4(1,i,k), pmp_2, lac_2) )
1026 : endif
1027 : endif
1028 : enddo
1029 0 : do i=i1,i2
1030 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1031 : enddo
1032 : else ! kord = 11, 13
1033 0 : do i=i1,i2
1034 0 : if ( ext5(i,k) .and. (ext5(i,k-1).or.ext5(i,k+1).or.a4(1,i,k)<qmin) ) then
1035 : ! Noisy region:
1036 0 : a4(2,i,k) = a4(1,i,k)
1037 0 : a4(3,i,k) = a4(1,i,k)
1038 0 : a4(4,i,k) = 0.0_r8
1039 : else
1040 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1041 : endif
1042 : enddo
1043 : endif
1044 :
1045 : ! Additional constraint to ensure positivity
1046 0 : if ( iv==0 ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 0)
1047 :
1048 : enddo ! k-loop
1049 :
1050 : !----------------------------------
1051 : ! Bottom layer subgrid constraints:
1052 : !----------------------------------
1053 0 : if ( iv==0 ) then
1054 0 : do i=i1,i2
1055 0 : a4(3,i,km) = max(0.0_r8, a4(3,i,km))
1056 : enddo
1057 0 : elseif ( iv .eq. -1 ) then
1058 0 : do i=i1,i2
1059 0 : if ( a4(3,i,km)*a4(1,i,km) <= 0.0_r8 ) a4(3,i,km) = 0.0_r8
1060 : enddo
1061 : endif
1062 :
1063 0 : do k=km-1,km
1064 0 : do i=i1,i2
1065 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1066 : enddo
1067 0 : if(k==(km-1)) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 2)
1068 0 : if(k== km ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1)
1069 : enddo
1070 :
1071 : end subroutine scalar_profile
1072 :
1073 :
1074 171428400 : subroutine cs_profile(qs, a4, delp, km, i1, i2, iv, kord)
1075 : ! Optimized vertical profile reconstruction:
1076 : ! Latest: Apr 2008 S.-J. Lin, NOAA/GFDL
1077 : integer, intent(in):: i1, i2
1078 : integer, intent(in):: km !< vertical dimension
1079 : integer, intent(in):: iv !< iv =-1: winds
1080 : !< iv = 0: positive definite scalars
1081 : !< iv = 1: others
1082 : integer, intent(in):: kord
1083 : real(kind=r8), intent(in) :: qs(i1:i2)
1084 : real(kind=r8), intent(in) :: delp(i1:i2,km) !< layer pressure thickness
1085 : real(kind=r8), intent(inout):: a4(4,i1:i2,km) !< Interpolated values
1086 : !-----------------------------------------------------------------------
1087 342856800 : logical, dimension(i1:i2,km):: extm, ext5, ext6
1088 342856800 : real(kind=r8) gam(i1:i2,km)
1089 342856800 : real(kind=r8) q(i1:i2,km+1)
1090 171428400 : real(kind=r8) d4(i1:i2)
1091 : real(kind=r8) bet, a_bot, grat
1092 : real(kind=r8) pmp_1, lac_1, pmp_2, lac_2, x0, x1
1093 : integer i, k, im
1094 :
1095 171428400 : if ( iv .eq. -2 ) then
1096 0 : do i=i1,i2
1097 0 : gam(i,2) = 0.5_r8
1098 0 : q(i,1) = 1.5_r8*a4(1,i,1)
1099 : enddo
1100 0 : do k=2,km-1
1101 0 : do i=i1, i2
1102 0 : grat = delp(i,k-1) / delp(i,k)
1103 0 : bet = 2.0_r8 + grat + grat - gam(i,k)
1104 0 : q(i,k) = (3.0_r8*(a4(1,i,k-1)+a4(1,i,k)) - q(i,k-1))/bet
1105 0 : gam(i,k+1) = grat / bet
1106 : enddo
1107 : enddo
1108 0 : do i=i1,i2
1109 0 : grat = delp(i,km-1) / delp(i,km)
1110 0 : q(i,km) = (3.0_r8*(a4(1,i,km-1)+a4(1,i,km)) - grat*qs(i) - q(i,km-1)) / &
1111 0 : (2.0_r8 + grat + grat - gam(i,km))
1112 0 : q(i,km+1) = qs(i)
1113 : enddo
1114 0 : do k=km-1,1,-1
1115 0 : do i=i1,i2
1116 0 : q(i,k) = q(i,k) - gam(i,k+1)*q(i,k+1)
1117 : enddo
1118 : enddo
1119 : else
1120 810388800 : do i=i1,i2
1121 638960400 : grat = delp(i,2) / delp(i,1) ! grid ratio
1122 638960400 : bet = grat*(grat+0.5_r8)
1123 638960400 : q(i,1) = ( (grat+grat)*(grat+1.0_r8)*a4(1,i,1) + a4(1,i,2) ) / bet
1124 810388800 : gam(i,1) = ( 1.0_r8 + grat*(grat+1.5_r8) ) / bet
1125 : enddo
1126 :
1127 4457138400 : do k=2,km
1128 20431148400 : do i=i1,i2
1129 15974010000 : d4(i) = delp(i,k-1) / delp(i,k)
1130 15974010000 : bet = 2.0_r8 + d4(i) + d4(i) - gam(i,k-1)
1131 15974010000 : q(i,k) = ( 3.0_r8*(a4(1,i,k-1)+d4(i)*a4(1,i,k)) - q(i,k-1) )/bet
1132 20259720000 : gam(i,k) = d4(i) / bet
1133 : enddo
1134 : enddo
1135 :
1136 810388800 : do i=i1,i2
1137 638960400 : a_bot = 1.0_r8 + d4(i)*(d4(i)+1.5_r8)
1138 1277920800 : q(i,km+1) = (2.0_r8*d4(i)*(d4(i)+1.0_r8)*a4(1,i,km)+a4(1,i,km-1)-a_bot*q(i,km)) &
1139 1449349200 : / ( d4(i)*(d4(i)+0.5_r8) - a_bot*gam(i,km) )
1140 : enddo
1141 :
1142 4628566800 : do k=km,1,-1
1143 21241537200 : do i=i1,i2
1144 21070108800 : q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
1145 : enddo
1146 : enddo
1147 : endif
1148 : !----- Perfectly linear scheme --------------------------------
1149 171428400 : if ( abs(kord) > 16 ) then
1150 0 : do k=1,km
1151 0 : do i=i1,i2
1152 0 : a4(2,i,k) = q(i,k )
1153 0 : a4(3,i,k) = q(i,k+1)
1154 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1155 : enddo
1156 : enddo
1157 : return
1158 : endif
1159 : !----- Perfectly linear scheme --------------------------------
1160 :
1161 : !------------------
1162 : ! Apply constraints
1163 : !------------------
1164 171428400 : im = i2 - i1 + 1
1165 :
1166 : ! Apply *large-scale* constraints
1167 810388800 : do i=i1,i2
1168 638960400 : q(i,2) = min( q(i,2), max(a4(1,i,1), a4(1,i,2)) )
1169 810388800 : q(i,2) = max( q(i,2), min(a4(1,i,1), a4(1,i,2)) )
1170 : enddo
1171 :
1172 4457138400 : do k=2,km
1173 20431148400 : do i=i1,i2
1174 20259720000 : gam(i,k) = a4(1,i,k) - a4(1,i,k-1)
1175 : enddo
1176 : enddo
1177 :
1178 : ! Interior:
1179 4114281600 : do k=3,km-1
1180 18810370800 : do i=i1,i2
1181 18638942400 : if ( gam(i,k-1)*gam(i,k+1)>0.0_r8 ) then
1182 : ! Apply large-scale constraint to ALL fields if not local max/min
1183 10799865931 : q(i,k) = min( q(i,k), max(a4(1,i,k-1),a4(1,i,k)) )
1184 10799865931 : q(i,k) = max( q(i,k), min(a4(1,i,k-1),a4(1,i,k)) )
1185 : else
1186 3896223269 : if ( gam(i,k-1) > 0.0_r8 ) then
1187 : ! There exists a local max
1188 1905969455 : q(i,k) = max(q(i,k), min(a4(1,i,k-1),a4(1,i,k)))
1189 : else
1190 : ! There exists a local min
1191 1990253814 : q(i,k) = min(q(i,k), max(a4(1,i,k-1),a4(1,i,k)))
1192 1990253814 : if ( iv==0 ) q(i,k) = max(0.0_r8, q(i,k))
1193 : endif
1194 : endif
1195 : enddo
1196 : enddo
1197 :
1198 : ! Bottom:
1199 810388800 : do i=i1,i2
1200 638960400 : q(i,km) = min( q(i,km), max(a4(1,i,km-1), a4(1,i,km)) )
1201 810388800 : q(i,km) = max( q(i,km), min(a4(1,i,km-1), a4(1,i,km)) )
1202 : enddo
1203 :
1204 4628566800 : do k=1,km
1205 21241537200 : do i=i1,i2
1206 16612970400 : a4(2,i,k) = q(i,k )
1207 21070108800 : a4(3,i,k) = q(i,k+1)
1208 : enddo
1209 : enddo
1210 :
1211 4628566800 : do k=1,km
1212 4457138400 : if ( k==1 .or. k==km ) then
1213 1620777600 : do i=i1,i2
1214 1620777600 : extm(i,k) = (a4(2,i,k)-a4(1,i,k)) * (a4(3,i,k)-a4(1,i,k)) > 0.0_r8
1215 : enddo
1216 : else
1217 19449331200 : do i=i1,i2
1218 19449331200 : extm(i,k) = gam(i,k)*gam(i,k+1) < 0.0_r8
1219 : enddo
1220 : endif
1221 4628566800 : if ( abs(kord) > 9 ) then
1222 0 : do i=i1,i2
1223 0 : x0 = 2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k))
1224 0 : x1 = abs(a4(2,i,k)-a4(3,i,k))
1225 0 : a4(4,i,k) = 3.0_r8*x0
1226 0 : ext5(i,k) = abs(x0) > x1
1227 0 : ext6(i,k) = abs(a4(4,i,k)) > x1
1228 : enddo
1229 : endif
1230 : enddo
1231 :
1232 : !---------------------------
1233 : ! Apply subgrid constraints:
1234 : !---------------------------
1235 : ! f(s) = AL + s*[(AR-AL) + A6*(1-s)] ( 0 <= s <= 1 )
1236 : ! Top 2 and bottom 2 layers always use monotonic mapping
1237 :
1238 171428400 : if ( iv==0 ) then
1239 498700800 : do i=i1,i2
1240 498700800 : a4(2,i,1) = max(0.0_r8, a4(2,i,1))
1241 : enddo
1242 62337600 : elseif ( iv==-1 ) then
1243 207792000 : do i=i1,i2
1244 207792000 : if ( a4(2,i,1)*a4(1,i,1) <= 0.0_r8 ) a4(2,i,1) = 0.0_r8
1245 : enddo
1246 20779200 : elseif ( iv==2 ) then
1247 0 : do i=i1,i2
1248 0 : a4(2,i,1) = a4(1,i,1)
1249 0 : a4(3,i,1) = a4(1,i,1)
1250 0 : a4(4,i,1) = 0.0_r8
1251 : enddo
1252 : endif
1253 :
1254 171428400 : if ( iv/=2 ) then
1255 810388800 : do i=i1,i2
1256 810388800 : a4(4,i,1) = 3.0_r8*(2.0_r8*a4(1,i,1) - (a4(2,i,1)+a4(3,i,1)))
1257 : enddo
1258 171428400 : call cs_limiters(im, extm(i1,1), a4(1,i1,1), 1)
1259 : endif
1260 :
1261 810388800 : do i=i1,i2
1262 810388800 : a4(4,i,2) = 3.0_r8*(2.0_r8*a4(1,i,2) - (a4(2,i,2)+a4(3,i,2)))
1263 : enddo
1264 171428400 : call cs_limiters(im, extm(i1,2), a4(1,i1,2), 2)
1265 :
1266 : !-------------------------------------
1267 : ! Huynh's 2nd constraint for interior:
1268 : !-------------------------------------
1269 3942853200 : do k=3,km-2
1270 3771424800 : if ( abs(kord)<9 ) then
1271 0 : do i=i1,i2
1272 : ! Left edges
1273 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1274 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1275 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1276 0 : max(a4(1,i,k), pmp_1, lac_1) )
1277 : ! Right edges
1278 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1279 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1280 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1281 0 : max(a4(1,i,k), pmp_2, lac_2) )
1282 :
1283 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1284 : enddo
1285 :
1286 : elseif ( abs(kord)==9 ) then
1287 17828553600 : do i=i1,i2
1288 17828553600 : if ( extm(i,k) .and. extm(i,k-1) ) then ! c90_mp122
1289 : ! grid-scale 2-delta-z wave detected
1290 842272345 : a4(2,i,k) = a4(1,i,k)
1291 842272345 : a4(3,i,k) = a4(1,i,k)
1292 842272345 : a4(4,i,k) = 0.0_r8
1293 13214856455 : else if ( extm(i,k) .and. extm(i,k+1) ) then ! c90_mp122
1294 : ! grid-scale 2-delta-z wave detected
1295 605102218 : a4(2,i,k) = a4(1,i,k)
1296 605102218 : a4(3,i,k) = a4(1,i,k)
1297 605102218 : a4(4,i,k) = 0.0_r8
1298 : else
1299 12609754237 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
1300 : ! Check within the smooth region if subgrid profile is non-monotonic
1301 12609754237 : if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
1302 5119002489 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1303 5119002489 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1304 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1305 5119002489 : max(a4(1,i,k), pmp_1, lac_1) )
1306 5119002489 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1307 5119002489 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1308 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1309 5119002489 : max(a4(1,i,k), pmp_2, lac_2) )
1310 5119002489 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
1311 : endif
1312 : endif
1313 : enddo
1314 : elseif ( abs(kord)==10 ) then
1315 0 : do i=i1,i2
1316 0 : if( ext5(i,k) ) then
1317 0 : if( ext5(i,k-1) .or. ext5(i,k+1) ) then
1318 0 : a4(2,i,k) = a4(1,i,k)
1319 0 : a4(3,i,k) = a4(1,i,k)
1320 0 : elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
1321 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1322 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1323 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1324 0 : max(a4(1,i,k), pmp_1, lac_1) )
1325 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1326 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1327 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1328 0 : max(a4(1,i,k), pmp_2, lac_2) )
1329 : endif
1330 0 : elseif( ext6(i,k) ) then
1331 0 : if( ext5(i,k-1) .or. ext5(i,k+1) ) then
1332 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1333 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1334 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1335 0 : max(a4(1,i,k), pmp_1, lac_1) )
1336 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1337 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1338 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1339 0 : max(a4(1,i,k), pmp_2, lac_2) )
1340 : endif
1341 : endif
1342 : enddo
1343 0 : do i=i1,i2
1344 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1345 : enddo
1346 : elseif ( abs(kord)==12 ) then
1347 0 : do i=i1,i2
1348 0 : if( extm(i,k) ) then
1349 : ! grid-scale 2-delta-z wave detected
1350 0 : a4(2,i,k) = a4(1,i,k)
1351 0 : a4(3,i,k) = a4(1,i,k)
1352 0 : a4(4,i,k) = 0.0_r8
1353 : else ! not a local extremum
1354 0 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
1355 : ! Check within the smooth region if subgrid profile is non-monotonic
1356 0 : if( abs(a4(4,i,k)) > abs(a4(2,i,k)-a4(3,i,k)) ) then
1357 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1358 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1359 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1360 0 : max(a4(1,i,k), pmp_1, lac_1) )
1361 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1362 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1363 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1364 0 : max(a4(1,i,k), pmp_2, lac_2) )
1365 0 : a4(4,i,k) = 6.0_r8*a4(1,i,k) - 3.0_r8*(a4(2,i,k)+a4(3,i,k))
1366 : endif
1367 : endif
1368 : enddo
1369 : elseif ( abs(kord)==13 ) then
1370 0 : do i=i1,i2
1371 0 : if( ext6(i,k) ) then
1372 0 : if ( ext6(i,k-1) .and. ext6(i,k+1) ) then
1373 : ! grid-scale 2-delta-z wave detected
1374 0 : a4(2,i,k) = a4(1,i,k)
1375 0 : a4(3,i,k) = a4(1,i,k)
1376 : endif
1377 : endif
1378 : enddo
1379 0 : do i=i1,i2
1380 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1381 : enddo
1382 : elseif ( abs(kord)==14 ) then
1383 :
1384 0 : do i=i1,i2
1385 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1386 : enddo
1387 :
1388 : elseif ( abs(kord)==15 ) then ! revised kord=9 scehem
1389 0 : do i=i1,i2
1390 0 : if ( ext5(i,k) ) then ! c90_mp122
1391 0 : if ( ext5(i,k-1) .or. ext5(i,k+1) ) then ! c90_mp122
1392 : ! grid-scale 2-delta-z wave detected
1393 0 : a4(2,i,k) = a4(1,i,k)
1394 0 : a4(3,i,k) = a4(1,i,k)
1395 : endif
1396 0 : elseif( ext6(i,k) ) then
1397 : ! Check within the smooth region if subgrid profile is non-monotonic
1398 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1399 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1400 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1401 0 : max(a4(1,i,k), pmp_1, lac_1) )
1402 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1403 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1404 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1405 0 : max(a4(1,i,k), pmp_2, lac_2) )
1406 : endif
1407 : enddo
1408 0 : do i=i1,i2
1409 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1410 : enddo
1411 : elseif ( abs(kord)==16 ) then
1412 0 : do i=i1,i2
1413 0 : if( ext5(i,k) ) then
1414 0 : if ( ext5(i,k-1) .or. ext5(i,k+1) ) then
1415 0 : a4(2,i,k) = a4(1,i,k)
1416 0 : a4(3,i,k) = a4(1,i,k)
1417 0 : elseif ( ext6(i,k-1) .or. ext6(i,k+1) ) then
1418 : ! Left edges
1419 0 : pmp_1 = a4(1,i,k) - 2.0_r8*gam(i,k+1)
1420 0 : lac_1 = pmp_1 + 1.5_r8*gam(i,k+2)
1421 : a4(2,i,k) = min(max(a4(2,i,k), min(a4(1,i,k), pmp_1, lac_1)), &
1422 0 : max(a4(1,i,k), pmp_1, lac_1) )
1423 : ! Right edges
1424 0 : pmp_2 = a4(1,i,k) + 2.0_r8*gam(i,k)
1425 0 : lac_2 = pmp_2 - 1.5_r8*gam(i,k-1)
1426 : a4(3,i,k) = min(max(a4(3,i,k), min(a4(1,i,k), pmp_2, lac_2)), &
1427 0 : max(a4(1,i,k), pmp_2, lac_2) )
1428 : endif
1429 : endif
1430 : enddo
1431 0 : do i=i1,i2
1432 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1433 : enddo
1434 : else ! kord = 11
1435 0 : do i=i1,i2
1436 0 : if ( ext5(i,k) .and. (ext5(i,k-1) .or. ext5(i,k+1)) ) then
1437 : ! Noisy region:
1438 0 : a4(2,i,k) = a4(1,i,k)
1439 0 : a4(3,i,k) = a4(1,i,k)
1440 0 : a4(4,i,k) = 0.0_r8
1441 : else
1442 0 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1443 : endif
1444 : enddo
1445 : endif
1446 :
1447 : ! Additional constraint to ensure positivity
1448 3942853200 : if ( iv==0 ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 0)
1449 :
1450 : enddo ! k-loop
1451 :
1452 : !----------------------------------
1453 : ! Bottom layer subgrid constraints:
1454 : !----------------------------------
1455 171428400 : if ( iv==0 ) then
1456 498700800 : do i=i1,i2
1457 498700800 : a4(3,i,km) = max(0.0_r8, a4(3,i,km))
1458 : enddo
1459 62337600 : elseif ( iv .eq. -1 ) then
1460 207792000 : do i=i1,i2
1461 207792000 : if ( a4(3,i,km)*a4(1,i,km) <= 0.0_r8 ) a4(3,i,km) = 0.0_r8
1462 : enddo
1463 : endif
1464 :
1465 514285200 : do k=km-1,km
1466 1620777600 : do i=i1,i2
1467 1620777600 : a4(4,i,k) = 3.0_r8*(2.0_r8*a4(1,i,k) - (a4(2,i,k)+a4(3,i,k)))
1468 : enddo
1469 342856800 : if(k==(km-1)) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 2)
1470 514285200 : if(k== km ) call cs_limiters(im, extm(i1,k), a4(1,i1,k), 1)
1471 : enddo
1472 :
1473 : end subroutine cs_profile
1474 :
1475 3085711200 : subroutine cs_limiters(im, extm, a4, iv)
1476 : integer, intent(in) :: im
1477 : integer, intent(in) :: iv
1478 : logical, intent(in) :: extm(im)
1479 : real(kind=r8) , intent(inout) :: a4(4,im) !< PPM array
1480 : ! LOCAL VARIABLES:
1481 : real(kind=r8) da1, da2, a6da
1482 : integer i
1483 :
1484 3085711200 : if ( iv==0 ) then
1485 : ! Positive definite constraint
1486 10971417600 : do i=1,im
1487 10971417600 : if( a4(1,i)<=0.0_r8) then
1488 1567798095 : a4(2,i) = a4(1,i)
1489 1567798095 : a4(3,i) = a4(1,i)
1490 1567798095 : a4(4,i) = 0.0_r8
1491 : else
1492 7003621905 : if( abs(a4(3,i)-a4(2,i)) < -a4(4,i) ) then
1493 1258062329 : if( (a4(1,i)+0.25_r8*(a4(3,i)-a4(2,i))**2/a4(4,i)+a4(4,i)*r12) < 0.0_r8 ) then
1494 : ! local minimum is negative
1495 275120694 : if( a4(1,i)<a4(3,i) .and. a4(1,i)<a4(2,i) ) then
1496 891701 : a4(3,i) = a4(1,i)
1497 891701 : a4(2,i) = a4(1,i)
1498 891701 : a4(4,i) = 0.0_r8
1499 274228993 : elseif( a4(3,i) > a4(2,i) ) then
1500 152306596 : a4(4,i) = 3.0_r8*(a4(2,i)-a4(1,i))
1501 152306596 : a4(3,i) = a4(2,i) - a4(4,i)
1502 : else
1503 121922397 : a4(4,i) = 3.0_r8*(a4(3,i)-a4(1,i))
1504 121922397 : a4(2,i) = a4(3,i) - a4(4,i)
1505 : endif
1506 : endif
1507 : endif
1508 : endif
1509 : enddo
1510 685713600 : elseif ( iv==1 ) then
1511 1620777600 : do i=1,im
1512 1620777600 : if( (a4(1,i)-a4(2,i))*(a4(1,i)-a4(3,i))>=0.0_r8 ) then
1513 462717819 : a4(2,i) = a4(1,i)
1514 462717819 : a4(3,i) = a4(1,i)
1515 462717819 : a4(4,i) = 0.0_r8
1516 : else
1517 815202981 : da1 = a4(3,i) - a4(2,i)
1518 815202981 : da2 = da1**2
1519 815202981 : a6da = a4(4,i)*da1
1520 815202981 : if(a6da < -da2) then
1521 174978013 : a4(4,i) = 3.0_r8*(a4(2,i)-a4(1,i))
1522 174978013 : a4(3,i) = a4(2,i) - a4(4,i)
1523 640224968 : elseif(a6da > da2) then
1524 65366476 : a4(4,i) = 3.0_r8*(a4(3,i)-a4(1,i))
1525 65366476 : a4(2,i) = a4(3,i) - a4(4,i)
1526 : endif
1527 : endif
1528 : enddo
1529 : else
1530 : ! Standard PPM constraint
1531 1620777600 : do i=1,im
1532 1620777600 : if( extm(i) ) then
1533 175280777 : a4(2,i) = a4(1,i)
1534 175280777 : a4(3,i) = a4(1,i)
1535 175280777 : a4(4,i) = 0.0_r8
1536 : else
1537 1102640023 : da1 = a4(3,i) - a4(2,i)
1538 1102640023 : da2 = da1**2
1539 1102640023 : a6da = a4(4,i)*da1
1540 1102640023 : if(a6da < -da2) then
1541 205743392 : a4(4,i) = 3.0_r8*(a4(2,i)-a4(1,i))
1542 205743392 : a4(3,i) = a4(2,i) - a4(4,i)
1543 896896631 : elseif(a6da > da2) then
1544 191993177 : a4(4,i) = 3.0_r8*(a4(3,i)-a4(1,i))
1545 191993177 : a4(2,i) = a4(3,i) - a4(4,i)
1546 : endif
1547 : endif
1548 : enddo
1549 : endif
1550 3085711200 : end subroutine cs_limiters
1551 :
1552 :
1553 0 : subroutine fillz(im, km, nq, q, dp)
1554 : integer, intent(in):: im !< No. of longitudes
1555 : integer, intent(in):: km !< No. of levels
1556 : integer, intent(in):: nq !< Total number of tracers
1557 : real(kind=r8) , intent(in):: dp(im,km) !< pressure thickness
1558 : real(kind=r8) , intent(inout) :: q(im,km,nq) !< tracer mixing ratio
1559 : ! LOCAL VARIABLES:
1560 0 : logical:: zfix(im)
1561 0 : real(kind=r8) :: dm(km)
1562 : integer i, k, ic, k1
1563 : real(kind=r8) qup, qly, dup, dq, sum0, sum1, fac
1564 :
1565 0 : do ic=1,nq
1566 : #ifdef DEV_GFS_PHYS
1567 : ! Bottom up:
1568 : do k=km,2,-1
1569 : k1 = k-1
1570 : do i=1,im
1571 : if( q(i,k,ic) < 0.0_r8 ) then
1572 : q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
1573 : q(i,k ,ic) = 0.0_r8
1574 : endif
1575 : enddo
1576 : enddo
1577 : ! Top down:
1578 : do k=1,km-1
1579 : k1 = k+1
1580 : do i=1,im
1581 : if( q(i,k,ic) < 0.0_r8 ) then
1582 : q(i,k1,ic) = q(i,k1,ic) + q(i,k,ic)*dp(i,k)/dp(i,k1)
1583 : q(i,k ,ic) = 0.0_r8
1584 : endif
1585 : enddo
1586 : enddo
1587 : #else
1588 : ! Top layer
1589 0 : do i=1,im
1590 0 : if( q(i,1,ic) < 0.0_r8 ) then
1591 0 : q(i,2,ic) = q(i,2,ic) + q(i,1,ic)*dp(i,1)/dp(i,2)
1592 0 : q(i,1,ic) = 0.0_r8
1593 : endif
1594 : enddo
1595 :
1596 : ! Interior
1597 0 : zfix(:) = .false.
1598 0 : do k=2,km-1
1599 0 : do i=1,im
1600 0 : if( q(i,k,ic) < 0.0_r8 ) then
1601 0 : zfix(i) = .true.
1602 0 : if ( q(i,k-1,ic) > 0.0_r8 ) then
1603 : ! Borrow from above
1604 0 : dq = min ( q(i,k-1,ic)*dp(i,k-1), -q(i,k,ic)*dp(i,k) )
1605 0 : q(i,k-1,ic) = q(i,k-1,ic) - dq/dp(i,k-1)
1606 0 : q(i,k ,ic) = q(i,k ,ic) + dq/dp(i,k )
1607 : endif
1608 0 : if ( q(i,k,ic)<0.0_r8 .and. q(i,k+1,ic)>0.0_r8 ) then
1609 : ! Borrow from below:
1610 0 : dq = min ( q(i,k+1,ic)*dp(i,k+1), -q(i,k,ic)*dp(i,k) )
1611 0 : q(i,k+1,ic) = q(i,k+1,ic) - dq/dp(i,k+1)
1612 0 : q(i,k ,ic) = q(i,k ,ic) + dq/dp(i,k )
1613 : endif
1614 : endif
1615 : enddo
1616 : enddo
1617 :
1618 : ! Bottom layer
1619 0 : k = km
1620 0 : do i=1,im
1621 0 : if( q(i,k,ic)<0.0_r8 .and. q(i,k-1,ic)>0.0_r8) then
1622 0 : zfix(i) = .true.
1623 : ! Borrow from above
1624 0 : qup = q(i,k-1,ic)*dp(i,k-1)
1625 0 : qly = -q(i,k ,ic)*dp(i,k )
1626 0 : dup = min(qly, qup)
1627 0 : q(i,k-1,ic) = q(i,k-1,ic) - dup/dp(i,k-1)
1628 0 : q(i,k, ic) = q(i,k, ic) + dup/dp(i,k )
1629 : endif
1630 : enddo
1631 :
1632 : ! Perform final check and non-local fix if needed
1633 0 : do i=1,im
1634 0 : if ( zfix(i) ) then
1635 : sum0 = 0.0_r8
1636 0 : do k=2,km
1637 0 : dm(k) = q(i,k,ic)*dp(i,k)
1638 0 : sum0 = sum0 + dm(k)
1639 : enddo
1640 :
1641 0 : if ( sum0 > 0.0_r8 ) then
1642 : sum1 = 0.0_r8
1643 0 : do k=2,km
1644 0 : sum1 = sum1 + max(0.0_r8, dm(k))
1645 : enddo
1646 0 : fac = sum0 / sum1
1647 0 : do k=2,km
1648 0 : q(i,k,ic) = max(0.0_r8, fac*dm(k)/dp(i,k))
1649 : enddo
1650 : endif
1651 :
1652 : endif
1653 : enddo
1654 : #endif
1655 :
1656 : enddo
1657 0 : end subroutine fillz
1658 : end module fv_mapz
|
