UPP (develop)
Loading...
Searching...
No Matches
BNDLYR.f
Go to the documentation of this file.
1
36! ------------------------------------------------------------------------------------------
54
55 SUBROUTINE bndlyr(PBND,TBND,QBND,RHBND,UBND,VBND, &
56 WBND,OMGBND,PWTBND,QCNVBND,LVLBND)
57
58!
59!
60 use vrbls3d, only: pint, q, uh, vh, pmid, t, omga, wh, cwm
61 use masks, only: lmh
62 use params_mod, only: d00, gi, pq0, a2, a3, a4
63 use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, modelname, &
64 jsta_m, jend_m, im, nbnd, spval, ista_2l, iend_2u, ista_m, iend_m, ista, iend
65 use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1
66 use gridspec_mod, only: gridtype
67 use upp_physics, only: fpvsnew
68!
69 implicit none
70!
71! DECLARE VARIABLES.
72!
73 real,PARAMETER :: DPBND=30.e2
74 integer, dimension(ista:iend,jsta:jend,NBND),intent(inout) :: LVLBND
75 real, dimension(ista:iend,jsta:jend,NBND),intent(inout) :: PBND,TBND, &
76 qbnd,rhbnd,ubnd,vbnd,wbnd,omgbnd,pwtbnd,qcnvbnd
77
78 REAL Q1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U),V1D(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U), &
79 u1d(ista_2l:iend_2u,jsta_2l:jend_2u),qcnv1d(ista_2l:iend_2u,jsta_2l:jend_2u)
80!
81 REAL, ALLOCATABLE :: PBINT(:,:,:),QSBND(:,:,:)
82 REAL, ALLOCATABLE :: PSUM(:,:,:), QCNVG(:,:,:)
83 REAL, ALLOCATABLE :: PVSUM(:,:,:),NSUM(:,:,:)
84!
85 integer I,J,L,IE,IW,LL,LV,LBND
86 real DP,QSAT,PV1,PV2,PMV,RPSUM,RPVSUM,PMIN,PM,DELP,PMINV,DELPV
87 real es
88!
89!*****************************************************************************
90! START BNDLYR HERE
91!
92 ALLOCATE (pbint(ista_2l:iend_2u,jsta_2l:jend_2u,nbnd+1))
93 ALLOCATE (qsbnd(ista_2l:iend_2u,jsta_2l:jend_2u,nbnd))
94 ALLOCATE (psum(ista_2l:iend_2u,jsta_2l:jend_2u,nbnd))
95 ALLOCATE (qcnvg(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
96 ALLOCATE (pvsum(ista_2l:iend_2u,jsta_2l:jend_2u,nbnd))
97 ALLOCATE (nsum(ista_2l:iend_2u,jsta_2l:jend_2u,nbnd))
98!
99! LOOP OVER HORIZONTAL GRID. AT EACH MASS POINT COMPUTE
100! PRESSURE AT THE INTERFACE OF EACH BOUNDARY LAYER.
101!
102!$omp parallel do private(i,j)
103 DO j=jsta,jend
104 DO i=ista,iend
105 pbint(i,j,1) = pint(i,j,nint(lmh(i,j))+1)
106 ENDDO
107 ENDDO
108!
109 DO lbnd=2,nbnd+1
110!$omp parallel do private(i,j)
111 DO j=jsta,jend
112 DO i=ista,iend
113 pbint(i,j,lbnd) = pbint(i,j,lbnd-1) - dpbnd
114 ENDDO
115 ENDDO
116 ENDDO
117
118! COMPUTE MOISTURE CONVERGENCE FOR EVERY LEVEL
119 DO l=1,lm
120!$omp parallel do private(i,j)
121 DO j=jsta_2l,jend_2u
122 DO i=ista_2l,iend_2u
123 q1d(i,j) = q(i,j,l)
124 u1d(i,j) = uh(i,j,l)
125 v1d(i,j) = vh(i,j,l)
126 ENDDO
127 ENDDO
128 CALL calmcvg(q1d,u1d,v1d,qcnv1d)
129!$omp parallel do private(i,j)
130 DO j=jsta,jend
131 DO i=ista,iend
132 qcnvg(i,j,l)=qcnv1d(i,j)
133 ENDDO
134 ENDDO
135 ENDDO
136
137!
138! LOOP OVER HORIZONTAL. AT EACH MASS POINT COMPUTE
139! MASS WEIGHTED LAYER MEAN P, T, Q, U, V, OMEGA,
140! WAND PRECIPITABLE WATER IN EACH BOUNDARY LAYER FROM THE SURFACE UP.
141!
142!!$omp+ private(dp,pm,qsat)
143!!$omp parallel do private(i,j,lbnd,l,ie,iw,dp,pm,qsat,pv1,pv2,pmv)
144 DO lbnd=1,nbnd
145!$omp parallel do private(i,j)
146 DO j=jsta,jend
147 DO i=ista,iend
148 pbnd(i,j,lbnd) = d00
149 tbnd(i,j,lbnd) = d00
150 qbnd(i,j,lbnd) = d00
151 qsbnd(i,j,lbnd) = d00
152 rhbnd(i,j,lbnd) = d00
153 ubnd(i,j,lbnd) = d00
154 vbnd(i,j,lbnd) = d00
155 wbnd(i,j,lbnd) = d00
156 omgbnd(i,j,lbnd) = d00
157 lvlbnd(i,j,lbnd) = 0
158 nsum(i,j,lbnd) = 0
159 psum(i,j,lbnd) = d00
160 pvsum(i,j,lbnd) = d00
161 pwtbnd(i,j,lbnd) = d00
162 qcnvbnd(i,j,lbnd)= d00
163 ENDDO
164 ENDDO
165!
166!!$omp parallel do private(i,j,l,dp,pm,es,qsat)
167 DO l=1,lm
168!$omp parallel do private(i,j,dp,pm,es,qsat)
169 DO j=jsta,jend
170 DO i=ista,iend
171!
172 pm = pmid(i,j,l)
173 IF(pm<spval)THEN
174 IF((pbint(i,j,lbnd) >= pm).AND. &
175 (pbint(i,j,lbnd+1) <= pm)) THEN
176 dp = pint(i,j,l+1) - pint(i,j,l)
177 psum(i,j,lbnd) = psum(i,j,lbnd) + dp
178 nsum(i,j,lbnd) = nsum(i,j,lbnd) + 1
179 lvlbnd(i,j,lbnd) = lvlbnd(i,j,lbnd) + l
180 tbnd(i,j,lbnd) = tbnd(i,j,lbnd) + t(i,j,l)*dp
181 qbnd(i,j,lbnd) = qbnd(i,j,lbnd) + q(i,j,l)*dp
182 omgbnd(i,j,lbnd) = omgbnd(i,j,lbnd) + omga(i,j,l)*dp
183 IF(gridtype == 'A')THEN
184 ubnd(i,j,lbnd) = ubnd(i,j,lbnd) + uh(i,j,l)*dp
185 vbnd(i,j,lbnd) = vbnd(i,j,lbnd) + vh(i,j,l)*dp
186 END IF
187 wbnd(i,j,lbnd) = wbnd(i,j,lbnd) + wh(i,j,l)*dp
188 qcnvbnd(i,j,lbnd) = qcnvbnd(i,j,lbnd) + qcnvg(i,j,l)*dp
189 pwtbnd(i,j,lbnd) = pwtbnd(i,j,lbnd) &
190 + ( q(i,j,l)+cwm(i,j,l))*dp*gi
191 IF(modelname == 'GFS')THEN
192 es = min(fpvsnew(t(i,j,l)),pm)
193 qsat = con_eps*es/(pm+con_epsm1*es)
194 ELSE
195 qsat = pq0/pm*exp(a2*(t(i,j,l)-a3)/(t(i,j,l)-a4))
196 END IF
197 qsbnd(i,j,lbnd) = qsbnd(i,j,lbnd) + qsat*dp
198 ENDIF
199 ELSE !undeined grids
200 pbnd(i,j,lbnd)=spval
201 tbnd(i,j,lbnd)=spval
202 ubnd(i,j,lbnd)=spval
203 vbnd(i,j,lbnd)=spval
204 wbnd(i,j,lbnd)=spval
205 omgbnd(i,j,lbnd)=spval
206 qcnvbnd(i,j,lbnd)=spval
207 pwtbnd(i,j,lbnd)=spval
208 qbnd(i,j,lbnd)=spval
209 qsbnd(i,j,lbnd)=spval
210 rhbnd(i,j,lbnd)=spval
211 ENDIF
212 ENDDO
213 ENDDO
214 ENDDO
215
216
217 IF(gridtype=='E')THEN
218 CALL exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,1))
219 DO l=1,lm
220 CALL exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l+1))
221!$omp parallel do private(i,j,ie,iw,dp,pv1,pv2,pmv)
222 DO j=jsta_m,jend_m
223 DO i=ista_m,iend_m
224 ie = i+mod(j,2)
225 iw = i+mod(j,2)-1
226 pv1 = 0.25*(pint(iw,j,l) + pint(ie,j,l) &
227 +pint(i,j+1,l) + pint(i,j-1,l))
228 pv2 = 0.25*(pint(iw,j,l+1) + pint(ie,j,l+1) &
229 +pint(i,j+1,l+1) + pint(i,j-1,l+1))
230 dp = pv2-pv1
231 pmv = 0.5*(pv1+pv2)
232 IF((pbint(iw,j,lbnd)>=pmv).AND. &
233 (pbint(iw,j,lbnd+1)<=pmv)) THEN
234 pvsum(i,j,lbnd) = pvsum(i,j,lbnd) + dp
235 ubnd(i,j,lbnd) = ubnd(i,j,lbnd) + dp* uh(i,j,l)
236 vbnd(i,j,lbnd) = vbnd(i,j,lbnd) + dp*vh(i,j,l)
237 ENDIF
238!
239 ENDDO
240 ENDDO
241 ENDDO
242 ELSE IF (gridtype=='B')THEN
243 CALL exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,1))
244 DO l=1,lm
245 CALL exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l+1))
246!$omp parallel do private(i,j,ie,iw,dp,pv1,pv2,pmv)
247 DO j=jsta_m,jend_m
248 DO i=ista_m,iend_m
249 ie = i+1
250 iw = i
251 pv1 = 0.25*(pint(iw,j,l) + pint(ie,j,l) &
252 +pint(iw,j+1,l) + pint(ie,j+1,l))
253 pv2 = 0.25*(pint(iw,j,l+1) + pint(ie,j,l+1) &
254 +pint(iw,j+1,l+1) + pint(ie,j+1,l+1))
255 dp = pv2-pv1
256 pmv = 0.5*(pv1+pv2)
257 IF((pbint(iw,j,lbnd)>=pmv).AND. &
258 (pbint(iw,j,lbnd+1)<=pmv)) THEN
259 pvsum(i,j,lbnd) = pvsum(i,j,lbnd)+dp
260 ubnd(i,j,lbnd) = ubnd(i,j,lbnd)+uh(i,j,l)*dp
261 vbnd(i,j,lbnd) = vbnd(i,j,lbnd)+vh(i,j,l)*dp
262 ENDIF
263!
264 ENDDO
265 ENDDO
266 ENDDO
267 END IF
268
269 ENDDO ! end of lbnd loop
270!
271!!$omp+ private(rpsum)
272!$omp parallel do private(i,j,lbnd,rpsum,rpvsum)
273 DO lbnd=1,nbnd
274 DO j=jsta,jend
275 DO i=ista,iend
276 IF(psum(i,j,lbnd)/=0..AND.tbnd(i,j,lbnd)<spval)THEN
277 rpsum = 1./psum(i,j,lbnd)
278 lvlbnd(i,j,lbnd)= lvlbnd(i,j,lbnd)/nsum(i,j,lbnd)
279 pbnd(i,j,lbnd) = (pbint(i,j,lbnd)+pbint(i,j,lbnd+1))*0.5
280 tbnd(i,j,lbnd) = tbnd(i,j,lbnd)*rpsum
281 qbnd(i,j,lbnd) = qbnd(i,j,lbnd)*rpsum
282 qsbnd(i,j,lbnd) = qsbnd(i,j,lbnd)*rpsum
283 omgbnd(i,j,lbnd)= omgbnd(i,j,lbnd)*rpsum
284 IF(gridtype=='A')THEN
285 ubnd(i,j,lbnd) = ubnd(i,j,lbnd)*rpsum
286 vbnd(i,j,lbnd) = vbnd(i,j,lbnd)*rpsum
287 END IF
288 wbnd(i,j,lbnd) = wbnd(i,j,lbnd)*rpsum
289 IF(qcnvbnd(i,j,lbnd)<spval) &
290 qcnvbnd(i,j,lbnd) = qcnvbnd(i,j,lbnd)*rpsum
291 ENDIF
292 ENDDO
293 ENDDO
294
295 IF(gridtype=='E' .or. gridtype=='B')THEN
296 DO j=jsta_m,jend_m
297 DO i=ista_m,iend_m
298 IF(pvsum(i,j,lbnd)/=0.)THEN
299 rpvsum = 1./pvsum(i,j,lbnd)
300 ubnd(i,j,lbnd) = ubnd(i,j,lbnd)*rpvsum
301 vbnd(i,j,lbnd) = vbnd(i,j,lbnd)*rpvsum
302 ENDIF
303 ENDDO
304 ENDDO
305 END IF
306 ENDDO
307!
308! IF NO ETA MID LAYER PRESSURES FELL WITHIN A BND LYR,
309! FIND THE CLOSEST LAYER TO THE BND LYR AND ASSIGN THE VALUES THERE
310!
311!!$omp+ private(delp,dp,l,pm,pmin,qsat)
312!$omp parallel do private(i,j,lbnd,l,ll,ie,iw,pminv,delp,dp,pm,pmin,es, &
313!$omp& qsat,pmv,delpv,lv)
314 DO lbnd=1,nbnd
315 DO j=jsta,jend
316 DO i=ista,iend
317 IF(psum(i,j,lbnd)==0..AND.pbnd(i,j,lbnd)<spval)THEN
318 l = lm
319 pmin = 9999999.
320 pbnd(i,j,lbnd) = (pbint(i,j,lbnd)+pbint(i,j,lbnd+1))*0.5
321!
322 DO ll=1,lm
323 pm = pmid(i,j,ll)
324 delp = abs(pm-pbnd(i,j,lbnd))
325 IF(delp<pmin)THEN
326 pmin = delp
327 l = ll
328 ENDIF
329 ENDDO
330!
331 dp = pint(i,j,l+1)-pint(i,j,l)
332 pm = pmid(i,j,l)
333 lvlbnd(i,j,lbnd) = l
334 tbnd(i,j,lbnd) = t(i,j,l)
335 qbnd(i,j,lbnd) = q(i,j,l)
336 IF(gridtype == 'A')THEN
337 ubnd(i,j,lbnd) = uh(i,j,l)
338 vbnd(i,j,lbnd) = vh(i,j,l)
339 END IF
340 wbnd(i,j,lbnd) = wh(i,j,l)
341 qcnvbnd(i,j,lbnd) = qcnvg(i,j,l)
342 IF(modelname == 'GFS' .OR. modelname == 'FV3R')THEN
343 es = fpvsnew(t(i,j,l))
344 es = min(es,pm)
345 qsat = con_eps*es/(pm+con_epsm1*es)
346 ELSE
347 qsat=pq0/pm*exp(a2*(t(i,j,l)-a3)/(t(i,j,l)-a4))
348 END IF
349 qsbnd(i,j,lbnd) = qsat
350 omgbnd(i,j,lbnd) = omga(i,j,l)
351 pwtbnd(i,j,lbnd) = (q(i,j,l)+cwm(i,j,l))*dp*gi
352 ENDIF
353!
354! RH, BOUNDS CHECK
355!
356 IF(qsbnd(i,j,lbnd)/=0..AND.qbnd(i,j,lbnd)<spval)THEN
357 rhbnd(i,j,lbnd) = qbnd(i,j,lbnd)/qsbnd(i,j,lbnd)
358 IF (rhbnd(i,j,lbnd)>1.0) THEN
359 rhbnd(i,j,lbnd) = 1.0
360 qbnd(i,j,lbnd) = rhbnd(i,j,lbnd)*qsbnd(i,j,lbnd)
361 ENDIF
362 IF (rhbnd(i,j,lbnd)<0.01) THEN
363 rhbnd(i,j,lbnd) = 0.01
364 qbnd(i,j,lbnd) = rhbnd(i,j,lbnd)*qsbnd(i,j,lbnd)
365 ENDIF
366 ENDIF
367 ENDDO
368 ENDDO
369!
370 IF(gridtype == 'E')THEN
371 DO j=jsta_m,jend_m
372 DO i=ista_m,iend_m
373 IF(pvsum(i,j,lbnd)==0.)THEN
374 lv = lm
375 pminv = 9999999.
376 ie = i+mod(j,2)
377 iw = i+mod(j,2)-1
378!
379! PINT HALOS UPDATED ALREADY
380!
381 DO ll=1,lm
382 pmv = 0.125*(pint(iw,j,ll) + pint(ie,j,ll) + &
383 pint(i,j+1,ll) + pint(i,j-1,ll) + &
384 pint(iw,j,ll+1) + pint(ie,j,ll+1) + &
385 pint(i,j+1,ll+1) + pint(i,j-1,ll+1))
386 delpv = abs(pmv-pbnd(i,j,lbnd))
387 IF(delpv<pminv)THEN
388 pminv = delpv
389 lv = ll
390 ENDIF
391 ENDDO
392!
393 ubnd(i,j,lbnd) = uh(i,j,lv)
394 vbnd(i,j,lbnd) = vh(i,j,lv)
395 ENDIF
396 ENDDO
397 ENDDO
398! END IF
399
400 ELSE IF(gridtype=='B')THEN
401 DO j=jsta_m,jend_m
402 DO i=ista_m,iend_m
403 IF(pvsum(i,j,lbnd)==0.)THEN
404 lv=lm
405 pminv=9999999.
406 ie=i+1
407 iw=i
408!
409! PINT HALOS UPDATED ALREADY
410!
411 DO ll=1,lm
412 pmv=0.125*(pint(iw,j,ll)+pint(ie,j,ll)+ &
413 pint(iw,j+1,ll)+pint(ie,j+1,ll)+ &
414 pint(iw,j,ll+1)+pint(ie,j,ll+1)+ &
415 pint(iw,j+1,ll+1)+pint(ie,j+1,ll+1))
416 delpv=abs(pmv-pbnd(i,j,lbnd))
417 IF(delpv<pminv)THEN
418 pminv=delpv
419 lv=ll
420 ENDIF
421 ENDDO
422
423 ubnd(i,j,lbnd) = uh(i,j,lv)
424 vbnd(i,j,lbnd) = vh(i,j,lv)
425 ENDIF
426 ENDDO
427 ENDDO
428 END IF
429 ENDDO ! end of lbnd loop
430!
431 DEALLOCATE (pbint, qsbnd, psum, pvsum, qcnvg, nsum)
432!
433! END OF ROUTINE
434!
435 RETURN
436 END
437!
bndlyr
subroutine bndlyr(pbnd, tbnd, qbnd, rhbnd, ubnd, vbnd, wbnd, omgbnd, pwtbnd, qcnvbnd, lvlbnd)
Computes boundary layer fields.
Definition BNDLYR.f:57
calmcvg
subroutine calmcvg(q1d, u1d, v1d, qcnvg)
Subroutine that computes moisture convergence.
Definition CALMCVG.f:44
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27