| blob | 8845d5a78b873a452444eb1bde99e598784e8144 |
1 !WRF:model_layer:physics
2 !
3 !
4 !
5 !
6 module module_shcu_grims
7 !
8 integer,parameter :: nxtdp = 5001
9 integer,parameter :: nxthe = 241,nythe = 151
10 integer,parameter :: nxma = 151,nyma = 121
11 integer,parameter :: nxsvp = 7501
12 !
13 real,parameter :: t0c = 2.7315e+2
14 real,parameter :: psat = 6.1078e+2
15 real,parameter :: rd = 2.8705e+2
16 real,parameter :: rv = 4.6150e+2
17 real,parameter :: cp = 1.0046e+3
18 real,parameter :: hvap = 2.5000e+6
19 real,parameter :: cvap = 1.8460e+3
20 real,parameter :: cliq = 4.1855e+3
21 real,parameter :: cice = 2.1060E+3
22 real,parameter :: hsub = 2.8340E+6
23 real,parameter :: terrm = 1.e-4
24 !
25 real,parameter :: rocp=rd/cp
26 real,parameter :: cpor=cp/rd
27 real,parameter :: eps=rd/rv
28 real,parameter :: ttp=t0c+0.01
29 real,parameter :: psatk=psat*1.e-3
30 real,parameter :: psatb=psatk*1.e-2
31 real,parameter :: dldt=cvap-cliq,xa=-dldt/rv,xb=xa+hvap/(rv*ttp)
32 real,parameter :: dldti=cvap-cice,xai=-dldti/rv,xbi=xai+hsub/(rv*ttp)
33 !
34 real,save :: c1xma,c2xma,c1yma,c2yma,c1xpvs,c2xpvs
35 real,save :: c1xtdp,c2xtdp,c1xthe,c2xthe,c1ythe,c2ythe
36 real,save :: tbtdp(nxtdp)
37 real,save :: tbthe(nxthe,nythe)
38 real,save :: tbtma(nxma,nyma), tbqma(nxma,nyma)
39 real,save :: tbpvs(nxsvp)
40 contains
41 !
42 !-------------------------------------------------------------------------------
43 subroutine grims(qv3d,t3d,p3di,p3d,pi3d,z3di, &
44 wstar,hpbl,delta, &
45 rthshten,rqvshten, &
46 dt,g,xlv,rd,rv,rcp,p1000mb, &
47 kpbl2d,znu,raincv, &
48 ids,ide, jds,jde, kds,kde, &
49 ims,ime, jms,jme, kms,kme, &
50 its,ite, jts,jte, kts,kte)
51 !-------------------------------------------------------------------------------
52 implicit none
53 !-------------------------------------------------------------------------------
54 !
55 ! input argument
56 !
57 !-- qv3d 3d specific humidity (kgkg-1)
58 !-- t3d 3d temperature (k)
59 !-- p3di 3d pressure (pa) at interface level
60 !-- p3d 3d pressure (pa)
61 !-- pi3d 3d exner function (dimensionless)
62 !-- z3di 3d z at interface level (m)
63 !-- wstar convective velocity scale (ms-1) from pbl
64 !-- hpbl pbl height (m)
65 !-- delta entrainment layer depth (m)
66 !-- rthshten computed theta tendency due to shallow convection scheme
67 !-- rqvshten computed q_v tendency due to shallow convection scheme
68 !-- dt time step (s)
69 !-- g acceleration due to gravity (m/s^2)
70 !-- xlv latent heat of vaporization (j/kg)
71 !-- rd gas constant for dry air (j/kg/k)
72 !-- rv gas constant for water vapor (j/kg/k)
73 !-- kpbl2d k-index for pbl top
74 !-- raincv time-step precipitation from cumulus convection scheme
75 !-- znu eta values (sigma values)
76 !-- ids start index for i in domain
77 !-- ide end index for i in domain
78 !-- jds start index for j in domain
79 !-- jde end index for j in domain
80 !-- kds start index for k in domain
81 !-- kde end index for k in domain
82 !-- ims start index for i in memory
83 !-- ime end index for i in memory
84 !-- jms start index for j in memory
85 !-- jme end index for j in memory
86 !-- kms start index for k in memory
87 !-- kme end index for k in memory
88 !-- its start index for i in tile
89 !-- ite end index for i in tile
90 !-- jts start index for j in tile
91 !-- jte end index for j in tile
92 !-- kts start index for k in tile
93 !-- kte end index for k in tile
94 !
95 ! output argument
96 !-- rthshten computed theta tendency due to shallow convection scheme
97 !-- rqvshten computed q_v tendency due to shallow convection scheme
98 !-------------------------------------------------------------------------------
99 !
100 ! local
101 !
102 !-- icps cps index, =1 for deep convection
103 !-- pi2di 2d exner function at interface level (dimensionless)
104 !-- delp2di 2d pressuer depth (pa) between interface levels
105 !-- zl 2d z (m)
106 !-- t1 2d temperature (k) will be changed by shallow convection
107 !-- q1 2d specific humidity (kgkg-1) will be changed by shallow convection
108 !-- levshc maximum k-level for shallow convection
109 !
110 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
111 ims,ime, jms,jme, kms,kme, &
112 its,ite, jts,jte, kts,kte
113 !
114 real, intent(in ) :: dt,g,xlv,rd,rv,rcp,p1000mb
115 !
116 integer, dimension( ims:ime, jms:jme ) , &
117 intent(in ) :: kpbl2d
118 !
119 real, dimension( ims:ime, kms:kme, jms:jme ) , &
120 intent(in ) :: qv3d, &
121 t3d, &
122 p3di, &
123 p3d, &
124 pi3d, &
125 z3di
126 !
127 real, dimension( ims:ime, jms:jme ) , &
128 intent(in ) :: hpbl, &
129 wstar, &
130 delta, &
131 raincv
132 !
133 real, dimension( kms:kme ) , &
134 intent(in ) :: znu
135 !
136 real, dimension( ims:ime, kms:kme, jms:jme ) , &
137 optional , &
138 intent(inout) :: rthshten, &
139 rqvshten
140 !
141 ! local variables
142 !
143 integer :: i,j,k,levshc
144 real :: sigshc,rdelt
145 !
146 integer, dimension( its:ite ) :: icps
147 real, dimension( its:ite, kts:kte+1 ) :: pi2di
148 real, dimension( its:ite, kts:kte ) :: delp2di, &
149 zl, &
150 t1, &
151 q1
152 !
153 rdelt = 1.0/dt
154 sigshc = 0.6
155 levshc = 0
156 !
157 do k = kts,kte
158 if(znu(k).gt.sigshc) levshc = k + 1
159 enddo
160 !
161 pi2di(:,:) = 0.0
162 delp2di(:,:) = 0.0
163 zl(:,:) = 0.0
164 t1(:,:) = 0.0
165 q1(:,:) = 0.0
166 !
167 do j = jts,jte ! outmost j loop
168 icps(:) = 0
169 do k = kts,kte
170 do i = its,ite
171 pi2di(i,k) = (p3di(i,k,j)/p1000mb)**rcp
172 enddo
173 enddo
174 do i = its,ite
175 pi2di(i,kte+1) = (p3di(i,kte+1,j)/p1000mb)**rcp
176 enddo
177 !
178 do k = kts,kte
179 do i = its,ite
180 delp2di(i,k) = p3di(i,k,j)-p3di(i,k+1,j)
181 zl(i,k) = 0.5*(z3di(i,k,j)+z3di(i,k+1,j))
182 t1(i,k) = t3d(i,k,j)
183 q1(i,k) = qv3d(i,k,j)
184 enddo
185 enddo
186 !
187 do i = its,ite
188 if(raincv(i,j) .gt. 1.e-30) icps(i)=1
189 enddo
190 !
191 call grims2d(q=q1(its,kts),t=t1(its,kts),prsi=p3di(ims,kms,j), &
192 prsik=pi2di(its,kts),delprsi=delp2di(its,kts), &
193 prsl=p3d(ims,kms,j),prslk=pi3d(ims,kms,j),zl=zl(its,kts), &
194 wstar=wstar(ims,j),hpbl=hpbl(ims,j),delta=delta(ims,j), &
195 dt=dt,cp=cp,g=g,xlv=xlv,rd=rd,rv=rv, &
196 icps=icps(its),kpbl=kpbl2d(ims,j),levshc=levshc, &
197 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
198 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
199 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
200 !
201 if(present(rthshten).and.present(rqvshten)) then
202 do k = kts,kte
203 do i = its,ite
204 rthshten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
205 rqvshten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
206 enddo
207 enddo
208 endif
209 !
210 enddo ! end outmost j loop
211 !
212 end subroutine grims
213 !-------------------------------------------------------------------------------
214 !
215 !-------------------------------------------------------------------------------
216 subroutine grims2d(q,t,prsi,prsik,delprsi,prsl,prslk,zl, &
217 wstar,hpbl,delta, &
218 dt,cp,g,xlv,rd,rv, &
219 icps,kpbl,levshc, &
220 ids,ide, jds,jde, kds,kde, &
221 ims,ime, jms,jme, kms,kme, &
222 its,ite, jts,jte, kts,kte)
223 !-------------------------------------------------------------------------------
224 implicit none
225 !-------------------------------------------------------------------------------
226 !
227 !
228 ! this scheme applies an eddy-diffusion approach within the shallow convective
229 ! layer defined by the moist static energy profile and is coupled to the
230 ! ysu pbl properties. this scheme names after the grims
231 ! shallow convection scheme since it was developed/evaluated in grims.
232 !
233 ! coded by song-you hong (yonsei university; ysu) and
234 ! implemented into wrf by jihyeon jang, hyeyum hailey shin, junhong lee (ysu),
235 ! and wei wang (ncar) winter 2012
236 !
237 ! references:
238 ! hong et al. (2013, manuscript in preparation)
239 ! hong et al. (2013, asia-pacific j. atmos. sci.) the global/regional
240 ! integrated model system (grims)
241 !
242 !-------------------------------------------------------------------------------
243 !
244 integer, intent(in ) :: levshc, &
245 ids,ide, jds,jde, kds,kde, &
246 ims,ime, jms,jme, kms,kme, &
247 its,ite, jts,jte, kts,kte
248 !
249 real, intent(in ) :: dt,cp,g,xlv,rd,rv
250 !
251 integer, dimension( ims:ime ) , &
252 intent(in ) :: kpbl
253 !
254 real, dimension( ims:ime, kms:kme ) , &
255 intent(in ) :: prsi
256 !
257 real, dimension( its:ite, kts:kte+1 ) , &
258 intent(in ) :: prsik
259 !
260 real, dimension( its:ite, kts:kte ) , &
261 intent(in ) :: delprsi, &
262 zl
263 !
264 real, dimension( ims:ime, kms:kme ) , &
265 intent(in ) :: prsl, &
266 prslk
267 !
268 integer, dimension( its:ite ) , &
269 intent(in ) :: icps
270 !
271 real, dimension( its:ite, kts:kte ) , &
272 intent(inout) :: q, &
273 t
274 !
275 real, dimension( ims:ime ) , &
276 intent(in ) :: hpbl, &
277 wstar, &
278 delta
279 !
280 ! profile shape parameter
281 !
282 real,parameter :: pfac = 3.
283 !
284 ! maximum and minimum diffusivity
285 !
286 real,parameter :: xkzmax = 50., xkzmin = 0.001
287 !
288 ! maxium distance of a parcel to lcl (m)
289 !
290 real,parameter :: zdiffcr1 = 1000., zdiffcr2 = 1000.
291 !
292 ! bounds of parcel origin
293 !
294 integer,parameter :: kliftl=2,kliftu=2
295 !
296 ! scale factor for wstar
297 !
298 real,parameter :: wsfac = 1.47
299 !
300 ! local variables and arrays
301 !
302 logical :: lshc(its:ite),flg(ite-its+1)
303 integer :: i,ik,ik1,iku,k,k1,k2,kt,n2
304 integer :: index2(ite-its+1)
305 integer :: klcl(ite-its+1),kbot(ite-its+1)
306 integer :: ktop(ite-its+1)
307 integer :: lmin(ite-its+1)
308 integer :: kb(ite-its+1),kbcon(ite-its+1)
309 real :: eps,epsm1
310 real :: eldq,xkzh,cpdt,rtdls
311 real :: dmse,dtodsu,dtodsl,dsig,dsdz1,dsdz2
312 real :: q2((ite-its+1)*kte)
313 real :: t2((ite-its+1)*kte)
314 real :: al((ite-its+1)*(kte-1))
315 real :: ad((ite-its+1)*kte)
316 real :: au((ite-its+1)*(kte-1))
317 real :: delprsi2((ite-its+1)*kte)
318 real :: prsi2((ite-its+1)*kte),prsik2((ite-its+1)*kte)
319 real :: prsl2((ite-its+1)*kte),prslk2((ite-its+1)*kte)
320 real :: qeso2((ite-its+1)*kte),rh2(ite-its+1)
321 real :: depth(ite-its+1),zdiff1(ite-its+1),zdiff2(ite-its+1)
322 real :: hmin(ite-its+1),hmax(ite-its+1)
323 real :: z(1:(ite-its+1),kts:kte)
324 real :: heo(1:(ite-its+1),kts:kte)
325 real :: heso(1:(ite-its+1),kts:kte)
326 real :: pik,height,xkzfac
327 !-------------------------------------------------------------------------------
328 !
329 eps = rd/rv
330 epsm1 = eps-1
331 !
332 do i = its,ite
333 lshc(i)=.false.
334 enddo
335 !
336 ! check for moist static instability to trigger convection
337 !
338 do k = kts,levshc-1
339 do i = its,ite
340 if(icps(i).eq.0.and.kpbl(i).ge.2) then
341 eldq = xlv*(q(i,k)-q(i,k+1))
342 cpdt = cp*(t(i,k)-t(i,k+1))
343 rtdls = (prsl(i,k)-prsl(i,k+1))/prsi(i,k+1)*rd*0.5*(t(i,k)+t(i,k+1))
344 dmse = eldq+cpdt-rtdls
345 lshc(i) = lshc(i).or.dmse.gt.0.
346 endif
347 enddo
348 enddo
349 do i = its,ite
350 if(wstar(i).lt.0.001) lshc(i)=.false.
351 enddo
352 !
353 ! reset i-dimension for active clouds
354 !
355 n2 = 0
356 do i = its,ite
357 if(lshc(i)) then
358 n2 = n2+1
359 index2(n2) = i
360 endif
361 enddo
362 !
363 if(n2.eq.0) return
364 !
365 ! prepare the variables
366 !
367 do k = kts,levshc
368 do i = 1,n2
369 if(lshc(index2(i))) then
370 ik = (k-1)*n2+i
371 pik = prsl(index2(i),k)
372 q2(ik) = q(index2(i),k)
373 t2(ik) = t(index2(i),k)
374 delprsi2(ik) = delprsi(index2(i),k)
375 prsi2(ik) = prsi(index2(i),k)
376 prsl2(ik) = prsl(index2(i),k)
377 prsik2(ik)= prsik(index2(i),k)
378 prslk2(ik)= prslk(index2(i),k)
379 z(i,k) = zl(index2(i),k)
380 qeso2(ik) = fpvs_pa(t2(ik))
381 qeso2(ik) = eps * qeso2(ik) / (pik + epsm1 * qeso2(ik))
382 qeso2(ik) = max(qeso2(ik),1.E-8)
383 heo(i,k) = g * z(i,k) + cp* t2(ik) + xlv * q2(ik)
384 heso(i,k) = g * z(i,k) + cp* t2(ik) + xlv * qeso2(ik)
385 endif
386 enddo
387 enddo
388 !
389 ! determine largest moist static energy for updraft originating level
390 !
391 do i = 1,n2
392 if(lshc(index2(i))) then
393 hmax(i) = heo(i,1)
394 kb(i) = 1
395 kbcon(i) = levshc
396 flg(i) = .true.
397 zdiff1(i) = -1.0
398 zdiff2(i) = -1.0
399 endif
400 enddo
401 !
402 do k = kts+1,levshc-1
403 do i = 1,n2
404 if(lshc(index2(i))) then
405 if(heo(i,k).gt.hmax(i).and.k.le.kpbl(index2(i))) then
406 kb(i) = k
407 hmax(i) = heo(i,k)
408 endif
409 endif
410 enddo
411 enddo
412 !
413 ! check the buoyancy of a parcel by the distance to lcl and lfc
414 !
415 do k = kts+1,levshc-1
416 do i = 1,n2
417 if(lshc(index2(i)).and.flg(i)) then
418 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
419 flg(i) = .false.
420 kbcon(i) = k
421 endif
422 endif
423 enddo
424 enddo
425 !
426 do i = 1,n2
427 if(lshc(index2(i))) then
428 zdiff1(i) = z(i,kpbl(index2(i)))-z(i,kb(i))
429 zdiff2(i) = z(i,kbcon(i))-z(i,kb(I))
430 if(zdiff1(i).gt.zdiffcr1.or.zdiff1(i).lt.0.) lshc(index2(i)) = .false.
431 if(zdiff2(i).gt.zdiffcr2) lshc(index2(i)) = .false.
432 endif
433 enddo
434 !
435 ! compute moist adiabat and determine cloud top
436 !
437 call phys_moist_adiabat_pa(n2,levshc-1,kliftl,kliftu, &
438 prsl2,prsik2,prslk2,t2,q2, &
439 klcl,kbot,ktop,al,au,rd,rv, &
440 ids,ide, jds,jde, kds,kde, &
441 ims,ime, jms,jme, kms,kme, &
442 its,ite, jts,jte, kts,kte)
443 !
444 do i = 1,n2
445 if(lshc(index2(i))) then
446 kbot(i) = max(kpbl(index2(i)),2)
447 kbot(i) = min(kbot(i),levshc)
448 ktop(i) = ktop(i)+1
449 endif
450 enddo
451 !
452 ! revise the cloud top below minimum moist static energy
453 !
454 do i = 1,n2
455 if(lshc(index2(i))) then
456 hmin(i) = heo(i,kbot(i))
457 lmin(i) = levshc
458 endif
459 enddo
460 !
461 do k = kts,levshc
462 do i = 1,n2
463 if(lshc(index2(i))) then
464 if(heo(i,k).lt.hmin(i).and.k.ge.kbot(i).and.k.le.ktop(i)) then
465 lmin(i) = k + 1
466 hmin(i) = heo(i,k)
467 endif
468 endif
469 enddo
470 enddo
471 !
472 do i = 1,n2
473 if(lshc(index2(i))) then
474 ktop(i) = min(ktop(i),lmin(i))
475 endif
476 enddo
477 !
478 do i = 1,n2
479 if(lshc(index2(i))) then
480 if((ktop(i)-kbot(i)).le.1) then
481 lshc(index2(i)) = .false.
482 endif
483 endif
484 enddo
485 !
486 ! compute diffusion properties
487 !
488 do i = 1,n2
489 if(lshc(index2(i))) then
490 k = kbot(i)-1
491 ik=(k-1)*n2+i
492 rh2(i) = q2(ik)/qeso2(ik)
493 endif
494 enddo
495 !
496 k1 = levshc+1
497 k2 = 0
498 do i = 1,n2
499 if(.not.lshc(index2(i))) then
500 kbot(i) = levshc+1
501 ktop(i) = 0
502 else
503 depth(i) = z(i,ktop(i))-z(i,kbot(i))
504 endif
505 k1 = min(k1,kbot(i))
506 k2 = max(k2,ktop(i))
507 enddo
508 kt = k2-k1+1
509 if(kt.lt.2) return
510 !
511 ! set eddy viscosity coefficient xkzh at sigma interfaces
512 !
513 do i = 1,n2
514 ik = (k1-1)*n2+i
515 ad(ik) = 1.
516 enddo
517 !
518 do k = kts,levshc-1
519 do i = 1,n2
520 if(k.ge.kbot(i).and.k.lt.ktop(i)) then
521 ik = (k-1)*n2+i
522 iku = k*n2+i
523 rtdls = (prsl2(ik)-prsl2(iku))/prsi2(iku)*rd*0.5*(t2(ik)+t2(iku))
524 au(ik) = g/rtdls
525 endif
526 enddo
527 enddo
528 !
529 do k = k1,k2-1
530 do i = 1,n2
531 ik = (k-1)*n2+i
532 iku = k*n2+i
533 dtodsl = 2.*dt/delprsi2(ik)
534 dtodsu = 2.*dt/delprsi2(iku)
535 dsig = prsl2(ik)-prsl2(iku)
536 if(k.ge.kbot(i).and.k.lt.ktop(i)) then
537 height = z(i,k)-z(i,kbot(i))
538 xkzfac = rh2(i)*wsfac*wstar(index2(i))*delta(index2(i))
539 xkzh = min(max(xkzfac*(1.-(height+hpbl(index2(i))) &
540 /(depth(i)+hpbl(index2(i))))**pfac,xkzmin),xkzmax)
541 else
542 xkzh = 0.
543 endif
544 dsdz1 = xkzh*dsig*au(ik)*g/cp
545 dsdz2 = xkzh*dsig*au(ik)*au(ik)
546 au(ik) = -dtodsl*dsdz2
547 al(ik) = -dtodsu*dsdz2
548 ad(ik) = ad(ik)-au(ik)
549 ad(iku) = 1.-al(ik)
550 t2(ik) = t2(ik)+dtodsl*dsdz1
551 t2(iku) = t2(iku)-dtodsu*dsdz1
552 enddo
553 enddo
554 !
555 ! solve tri-diagonal matrix
556 !
557 ik1 = (k1-1)*n2+1
558 call scv_tri_diagonal_grims(n2,n2,kt,al(ik1),ad(ik1),au(ik1), &
559 q2(ik1),t2(ik1),au(ik1),q2(ik1),t2(ik1))
560 !
561 ! feedback to large-scale variables
562 !
563 do k = k1,k2
564 do i = 1,n2
565 if(lshc(index2(i))) then
566 ik = (k-1)*n2+i
567 q(index2(i),k) = q2(ik)
568 t(index2(i),k) = t2(ik)
569 endif
570 enddo
571 enddo
572 !
573 return
574 end subroutine grims2d
575 !-------------------------------------------------------------------------------
576 !
577 !-------------------------------------------------------------------------------
578 subroutine scv_tri_diagonal_grims(lons2,l,n,cl,cm,cu,r1,r2,au,a1,a2)
579 !-------------------------------------------------------------------------------
580 !
581 ! subprogram: scv_tri_diagonal
582 !
583 ! abstract: this routine solves multiple tridiagonal matrix problems
584 ! with 2 right-hand-side and solution vectors for every matrix.
585 ! the solutions are found by eliminating off-diagonal coefficients,
586 ! marching first foreward then backward along the matrix diagonal.
587 ! the computations are vectorized around the number of matrices.
588 ! no checks are made for zeroes on the diagonal or singularity.
589 !
590 ! program history log:
591 ! 1991-05-07 iredell
592 ! 2009-03-01 jung-eun kim fortran 90 and modules
593 !
594 ! usage: call scv_tri_diagonal(l,n,cl,cm,cu,r1,r2,au,a1,a2)
595 !
596 ! input argument list:
597 ! l - integer number of tridiagonal matrices
598 ! n - integer order of the matrices
599 ! cl - real (l,2:n) lower diagonal matrix elements
600 ! cm - real (l,n) main diagonal matrix elements
601 ! cu - real (l,n-1) upper diagonal matrix elements
602 ! (may be equivalent to au if no longer needed)
603 ! r1 - real (l,n) 1st right-hand-side vector elements
604 ! (may be equivalent to a1 if no longer needed)
605 ! r2 - real (l,n) 2nd right-hand-side vector elements
606 ! (may be equivalent to a2 if no longer needed)
607 !
608 ! output argument list:
609 ! au - real (l,n-1) work array
610 ! a1 - real (l,n) 1st solution vector elements
611 ! a2 - real (l,n) 2nd solution vector elements
612 !
613 ! remarks: this routine can be easily modified to solve a different
614 ! number of right-hand-sides and solutions per matrix besides 2.
615 !
616 !-------------------------------------------------------------------------------
617 implicit none
618 integer :: lons2,l,n,i,k
619 real :: fk
620 real :: cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), &
621 au(l,n-1),a1(l,n),a2(l,n)
622 !-------------------------------------------------------------------------------
623 do i = 1,lons2
624 fk = 1./cm(i,1)
625 au(i,1)= fk*cu(i,1)
626 a1(i,1)= fk*r1(i,1)
627 a2(i,1)= fk*r2(i,1)
628 enddo
629 !
630 do k = 2,n-1
631 do i = 1,lons2
632 fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1))
633 au(i,k) = fk*cu(i,k)
634 a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1))
635 a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1))
636 enddo
637 enddo
638 !
639 do i = 1,lons2
640 fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1))
641 a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1))
642 a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1))
643 enddo
644 do k = n-1,1,-1
645 do i = 1,lons2
646 a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1)
647 a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1)
648 enddo
649 enddo
650 !
651 return
652 end subroutine scv_tri_diagonal_grims
653 !-------------------------------------------------------------------------------
654 !
655 !-------------------------------------------------------------------------------
656 subroutine phys_moist_adiabat_pa(ilev,klev,k1,k2, &
657 prsl,prsik,prslk,tenv,qenv, &
658 klcl,kbot,ktop,tcld,qcld,rd,rv, &
659 ids,ide, jds,jde, kds,kde, &
660 ims,ime, jms,jme, kms,kme, &
661 its,ite, jts,jte, kts,kte)
662 !-------------------------------------------------------------------------------
663 !
664 ! subprogram: phys_moist_adiabat_pa
665 !
666 ! abstract:
667 ! - compute moist adiabatic cloud soundings
668 ! - atmospheric columns of temperature and specific humidity
669 ! are examined by this routine for conditional instability.
670 ! the test parcel is chosen from the layer between layers k1 and k2
671 ! that has the warmest potential wet-bulb temperature.
672 ! excess cloud temperatures and specific humidities are returned
673 ! where the lifted parcel is found to be buoyant.
674 ! fast inlinable functions are invoked to compute
675 ! dewpoint and lifting condensation level temperatures,
676 ! equivalent potential temperature at the lcl, and
677 ! temperature and specific humidity of the ascending parcel.
678 !
679 ! program history log:
680 ! 1983-11-01 phillips
681 ! 1991-05-07 iredell arguments changed, code tidied
682 ! 2000-01-01 song-you hong physcis options
683 ! 2009-10-01 jung-eun kim f90 format with standard physics modules
684 ! 2010-07-01 myung-seo koo dimension allocatable with namelist input
685 !
686 ! usage: call phys_moist_adiabat_pa(ilev,klev,k1,k2, &
687 ! prsl,prslk,prsik,tenv,qenv, &
688 ! klcl,kbot,ktop,tcld,qcld,rd,rv, &
689 ! ids,ide, jds,jde, kds,kde, &
690 ! ims,ime, jms,jme, kms,kme, &
691 ! its,ite, jts,jte, kts,kte)
692 !
693 ! input argument list:
694 ! ilev - integer number of atmospheric columns
695 ! klev - integer number of sigma levels in a column
696 ! k1 - integer lowest level from which a parcel can originate
697 ! k2 - integer highest level from which a parcel can originate
698 ! prsl - real (ilev,klev) pressure values
699 ! prslk,prsik - real (ilev,klev) pressure values to the kappa
700 ! tenv - real (ilev,klev) environment temperatures
701 ! qenv - real (ilev,klev) environment specific humidities
702 !
703 ! output argument list:
704 ! klcl - integer (ilev) level just above lcl (klev+1 if no lcl)
705 ! kbot - integer (ilev) level just above cloud bottom
706 ! ktop - integer (ilev) level just below cloud top
707 ! - note that kbot(i) gt ktop(i) if no cloud.
708 ! tcld - real (ilev,klev) of excess cloud temperatures.
709 ! (parcel t minus environ t, or 0. where no cloud)
710 ! qcld - real (ilev,klev) of excess cloud specific humidities.
711 ! (parcel q minus environ q, or 0. where no cloud)
712 !
713 ! subprograms called:
714 ! ftdp - function to compute dewpoint temperature
715 ! ftlcl - function to compute lcl temperature
716 ! fthe - function to compute equivalent potential temperature
717 ! ftma - function to compute parcel temperature and humidity
718 !
719 ! remarks: all functions are inlined by fpp.
720 ! nonstandard automatic arrays are used.
721 !
722 !-------------------------------------------------------------------------------
723 implicit none
724 !
725 integer,parameter :: nx=151,ny=121
726 integer :: ilev,klev,k1,k2
727 real :: prsl(ilev,klev),prslk(ilev,klev),prsik(ilev,klev)
728 real :: tenv(ilev,klev),qenv(ilev,klev)
729 integer :: klcl(ilev),kbot(ilev),ktop(ilev)
730 real :: tcld(ilev,klev),qcld(ilev,klev)
731 real :: rd,rv
732 integer :: ids,ide, jds,jde, kds,kde, &
733 ims,ime, jms,jme, kms,kme, &
734 its,ite, jts,jte, kts,kte
735 !
736 ! local arrays
737 !
738 real :: slkma(ilev)
739 real :: thema(ilev)
740 real :: pv,tdpd
741 real :: slklcl,thelcl,tlcl
742 real :: xj,yj
743 real :: ftx1,ftx2,ftma1,qx1,qx2,qma,tma,tvcld,tvenv
744 real :: eps,epsm1,ftv
745 integer :: i,k,jx,jy
746 !-------------------------------------------------------------------------------
747 !
748 ! compute parameters
749 !
750 eps=rd/rv
751 epsm1=rd/rv-1.
752 ftv=rv/rd-1.
753 !
754 ! determine warmest potential wet-bulb temperature between k1 and k2.
755 ! compute its lifting condensation level.
756 !
757 do i = 1,ilev
758 slkma(i)=0.
759 thema(i)=0.
760 enddo
761 !
762 do k = k1,k2
763 do i = 1,ilev
764 pv=prsl(i,k)*1.e-3*qenv(i,k)/(eps-epsm1*qenv(i,k))
765 tdpd=tenv(i,k)-ftdp(pv)
766 if(tdpd.gt.0.) then
767 tlcl=ftlcl(tenv(i,k),tdpd)
768 slklcl=prslk(i,k)/prsik(i,1)*tlcl/tenv(i,k)
769 else
770 tlcl=tenv(i,k)
771 slklcl=prslk(i,k)/prsik(i,1)
772 endif
773 thelcl=fthe(tlcl,slklcl*prsik(i,1))
774 if(thelcl.gt.thema(i)) then
775 slkma(i)=slklcl
776 thema(i)=thelcl
777 endif
778 enddo
779 enddo
780 !
781 ! set cloud temperatures and humidities wherever the parcel lifted up
782 ! the moist adiabat is buoyant with respect to the environment.
783 !
784 do i = 1,ilev
785 klcl(i)=klev+1
786 kbot(i)=klev+1
787 ktop(i)=0
788 enddo
789 !
790 do k = kts,klev
791 do i = 1,ilev
792 tcld(i,k)=0.
793 qcld(i,k)=0.
794 enddo
795 enddo
796 !
797 do k = k1,klev
798 do i = 1,ilev
799 if(prslk(i,k)/prsik(i,1).le.slkma(i)) then
800 klcl(i)=min(klcl(i),k)
801 !
802 ! insert ftma tma=ftma(thema(i),prslk(i,k),qma)
803 !
804 xj=min(max(c1xma+c2xma*thema(i),1.),float(nx))
805 yj=min(max(c1yma+c2yma*prslk(i,k),1.),float(ny))
806 jx=min(xj,nx-1.)
807 jy=min(yj,ny-1.)
808 ftx1=tbtma(jx,jy)+(xj-jx)*(tbtma(jx+1,jy)-tbtma(jx,jy))
809 ftx2=tbtma(jx,jy+1)+(xj-jx)*(tbtma(jx+1,jy+1)-tbtma(jx,jy+1))
810 ftma1=ftx1+(yj-jy)*(ftx2-ftx1)
811 qx1=tbqma(jx,jy)+(xj-jx)*(tbqma(jx+1,jy)-tbqma(jx,jy))
812 qx2=tbqma(jx,jy+1)+(xj-jx)*(tbqma(jx+1,jy+1)-tbqma(jx,jy+1))
813 qma=qx1+(yj-jy)*(qx2-qx1)
814 tma=ftma1
815 !
816 tvcld=tma*(1.+ftv*qma)
817 tvenv=tenv(i,k)*(1.+ftv*qenv(i,k))
818 if(tvcld.gt.tvenv) then
819 kbot(i)=min(kbot(i),k)
820 ktop(i)=max(ktop(i),k)
821 tcld(i,k)=tma-tenv(i,k)
822 qcld(i,k)=qma-qenv(i,k)
823 endif
824 endif
825 enddo
826 enddo
827 !
828 return
829 end subroutine phys_moist_adiabat_pa
830 !-------------------------------------------------------------------------------
831 !
832 !-------------------------------------------------------------------------------
833 function ftdp(pv)
834 !-------------------------------------------------------------------------------
835 implicit none
836 !-------------------------------------------------------------------------------
837 real :: ftdp
838 !
839 integer :: jx
840 real :: xj
841 real :: xmax,xmin,xinc
842 real :: pv
843 !
844 xmin= 0.001
845 xmax=10.001
846 xinc=(xmax-xmin)/(nxtdp-1)
847 c1xtdp=1.-xmin/xinc
848 c2xtdp=1./xinc
849 !
850 xj=min(max(c1xtdp+c2xtdp*pv,1.),float(nxtdp))
851 jx=min(xj,nxtdp-1.)
852 ftdp=tbtdp(jx)+(xj-jx)*(tbtdp(jx+1)-tbtdp(jx))
853 !
854 return
855 end function
856 !-------------------------------------------------------------------------------
857 !
858 !-------------------------------------------------------------------------------
859 function ftlcl(t,tdpd)
860 !-------------------------------------------------------------------------------
861 implicit none
862 !-------------------------------------------------------------------------------
863 real :: ftlcl
864 !
865 real,parameter :: clcl1=0.954442e+0, clcl2=0.967772e-3, &
866 clcl3=-0.710321e-3,clcl4=-0.270742e-5
867 real :: t,tdpd
868 !
869 ftlcl=t-tdpd*(clcl1+clcl2*t+tdpd*(clcl3+clcl4*t))
870 !
871 return
872 end function
873 !-------------------------------------------------------------------------------
874 !
875 !-------------------------------------------------------------------------------
876 function fthe(t,pk)
877 !-------------------------------------------------------------------------------
878 implicit none
879 !-------------------------------------------------------------------------------
880 real :: fthe
881 !
882 integer :: jx,jy
883 real :: xmin,xmax,xinc,ymin,ymax,yinc
884 real :: xj,yj
885 real :: ftx1,ftx2
886 real :: t,pk
887 !
888 xmin=ttp-90.
889 xmax=ttp+30.
890 xinc=(xmax-xmin)/(nxthe-1)
891 c1xthe=1.-xmin/xinc
892 c2xthe=1./xinc
893 ymin=0.04**rocp
894 ymax=1.10**rocp
895 yinc=(ymax-ymin)/(nythe-1)
896 c1ythe=1.-ymin/yinc
897 c2ythe=1./yinc
898 !
899 xj=min(c1xthe+c2xthe*t,float(nxthe))
900 yj=min(c1ythe+c2ythe*pk,float(nythe))
901 if(xj.ge.1..and.yj.ge.1.) then
902 jx=min(xj,nxthe-1.)
903 jy=min(yj,nythe-1.)
904 ftx1=tbthe(jx,jy)+(xj-jx)*(tbthe(jx+1,jy)-tbthe(jx,jy))
905 ftx2=tbthe(jx,jy+1)+(xj-jx)*(tbthe(jx+1,jy+1)-tbthe(jx,jy+1))
906 fthe=ftx1+(yj-jy)*(ftx2-ftx1)
907 else
908 fthe=0.
909 endif
910 !
911 return
912 end function
913 !-------------------------------------------------------------------------------
914 !
915 !-------------------------------------------------------------------------------
916 function fpvs_pa(t)
917 !-------------------------------------------------------------------------------
918 implicit none
919 !-------------------------------------------------------------------------------
920 real :: fpvs_pa
921 !
922 integer :: jx
923 real :: xmax,xmin,xinc
924 real :: xj
925 real :: t
926 !
927 xmin=180.0
928 xmax=330.0
929 !
930 xj=min(max(c1xpvs+c2xpvs*t,1.),float(nxsvp))
931 jx=min(xj,nxsvp-1.)
932 fpvs_pa=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
933 fpvs_pa=fpvs_pa * 1.e3
934 !
935 return
936 end function
937 !-------------------------------------------------------------------------------
938 !
939 !-------------------------------------------------------------------------------
940 subroutine grimsinit(rthshten,rqvshten, &
941 restart, &
942 ids,ide, jds,jde, kds,kde, &
943 ims,ime, jms,jme, kms,kme, &
944 its,ite, jts,jte, kts,kte )
945 !-------------------------------------------------------------------------------
946 implicit none
947 !-------------------------------------------------------------------------------
948 logical , intent(in) :: restart
949 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
950 ims, ime, jms, jme, kms, kme, &
951 its, ite, jts, jte, kts, kte
952 real, dimension( ims:ime , kms:kme , jms:jme ) , intent(out) :: &
953 rthshten, &
954 rqvshten
955 integer :: i, j, k, itf, jtf, ktf
956 !
957 jtf=min0(jte,jde-1)
958 ktf=min0(kte,kde-1)
959 itf=min0(ite,ide-1)
960 if(.not.restart)then
961 do j = jts,jtf
962 do k = kts,ktf
963 do i = its,itf
964 rthshten(i,k,j)=0.
965 rqvshten(i,k,j)=0.
966 enddo
967 enddo
968 enddo
969 endif
970 !
971 call funct_dew_point_temp_init
972 call funct_pot_temp_init
973 call funct_moist_adiabat_init
974 call funct_svp_init
975 !
976 end subroutine grimsinit
977 !-------------------------------------------------------------------------------
978 !
979 !-------------------------------------------------------------------------------
980 subroutine funct_dew_point_temp_init
981 !-------------------------------------------------------------------------------
982 implicit none
983 !-------------------------------------------------------------------------------
984 integer :: jx
985 real :: xmax,xmin,xinc,pv,x,t
986 !
987 xmin= 0.001
988 xmax=10.001
989 xinc=(xmax-xmin)/(nxtdp-1)
990 c1xtdp=1.-xmin/xinc
991 c2xtdp=1./xinc
992 t=208.0
993 do jx=1,nxtdp
994 x=xmin+(jx-1)*xinc
995 pv=x
996 t=ftdpxg(t,pv)
997 tbtdp(jx)=t
998 enddo
999 !
1000 return
1001 end subroutine funct_dew_point_temp_init
1002 !-------------------------------------------------------------------------------
1003 !
1004 !-------------------------------------------------------------------------------
1005 subroutine funct_pot_temp_init
1006 !-------------------------------------------------------------------------------
1007 implicit none
1008 !-------------------------------------------------------------------------------
1009 integer :: jx,jy
1010 real :: xmin,xmax,xinc,ymin,ymax,yinc
1011 real :: x,y,t,pk
1012 !
1013 xmin=ttp-90.
1014 xmax=ttp+30.
1015 xinc=(xmax-xmin)/(nxthe-1)
1016 c1xthe=1.-xmin/xinc
1017 c2xthe=1./xinc
1018 ymin=0.04**rocp
1019 ymax=1.10**rocp
1020 yinc=(ymax-ymin)/(nythe-1)
1021 c1ythe=1.-ymin/yinc
1022 c2ythe=1./yinc
1023 !
1024 do jy = 1,nythe
1025 y=ymin+(jy-1)*yinc
1026 pk=y
1027 do jx = 1,nxthe
1028 x=xmin+(jx-1)*xinc
1029 t=x
1030 tbthe(jx,jy)=fthex(t,pk)
1031 enddo
1032 enddo
1033 !
1034 return
1035 end subroutine funct_pot_temp_init
1036 !-------------------------------------------------------------------------------
1037 !
1038 !-------------------------------------------------------------------------------
1039 subroutine funct_moist_adiabat_init
1040 !-------------------------------------------------------------------------------
1041 implicit none
1042 !-------------------------------------------------------------------------------
1043 integer :: jx,jy
1044 real :: xmin,xmax,xinc,ymin,ymax,yinc
1045 real :: y,pk,t,x,the,q
1046 !
1047 xmin=200.
1048 xmax=500.
1049 xinc=(xmax-xmin)/(nxma-1)
1050 c1xma=1.-xmin/xinc
1051 c2xma=1./xinc
1052 ymin=0.01**rocp
1053 ymax=1.10**rocp
1054 yinc=(ymax-ymin)/(nyma-1)
1055 c1yma=1.-ymin/yinc
1056 c2yma=1./yinc
1057 !
1058 do jy = 1,nyma
1059 y=ymin+(jy-1)*yinc
1060 pk=y
1061 t=xmin*y
1062 do jx = 1,nxma
1063 x=xmin+(jx-1)*xinc
1064 the=x
1065 t=ftmaxg(t,the,pk,q)
1066 tbtma(jx,jy)=t
1067 tbqma(jx,jy)=q
1068 enddo
1069 enddo
1070 !
1071 return
1072 end subroutine funct_moist_adiabat_init
1073 !-------------------------------------------------------------------------------
1074 !
1075 !-------------------------------------------------------------------------------
1076 subroutine funct_svp_init
1077 !-------------------------------------------------------------------------------
1078 implicit none
1079 !-------------------------------------------------------------------------------
1080 integer :: jx
1081 real :: xmin,xmax,xinc
1082 real :: t,x
1083 !
1084 xmin=180.0
1085 xmax=330.0
1086 xinc=(xmax-xmin)/(nxsvp-1)
1087 c1xpvs=1.-xmin/xinc
1088 c2xpvs=1./xinc
1089 do jx = 1,nxsvp
1090 x=xmin+(jx-1)*xinc
1091 t=x
1092 tbpvs(jx)=fpvsx(t)
1093 enddo
1094 !
1095 return
1096 end subroutine funct_svp_init
1097 !-------------------------------------------------------------------------------
1098 !
1099 !-------------------------------------------------------------------------------
1100 function ftdpxg(tg,pv)
1101 !-------------------------------------------------------------------------------
1102 implicit none
1103 !-------------------------------------------------------------------------------
1104 real :: ftdpxg
1105 !
1106 real :: tg,pv
1107 real :: t,tr,pvt,el,dpvt,terr
1108 !
1109 t=tg
1110 tr=ttp/t
1111 pvt=psatk*(tr**xa)*exp(xb*(1.-tr))
1112 el=hvap+dldt*(t-ttp)
1113 dpvt=el*pvt/(rv*t**2)
1114 terr=(pvt-pv)/dpvt
1115 t=t-terr
1116 !
1117 do while(abs(terr).gt.terrm)
1118 tr=ttp/t
1119 pvt=psatk*(tr**xa)*exp(xb*(1.-tr))
1120 el=hvap+dldt*(t-ttp)
1121 dpvt=el*pvt/(rv*t**2)
1122 terr=(pvt-pv)/dpvt
1123 t=t-terr
1124 enddo
1125 ftdpxg=t
1126 !
1127 return
1128 end function
1129 !-------------------------------------------------------------------------------
1130 !
1131 !-------------------------------------------------------------------------------
1132 function fthex(t,pk)
1133 !-------------------------------------------------------------------------------
1134 implicit none
1135 !-------------------------------------------------------------------------------
1136 real :: fthex
1137 !
1138 real :: t,pk,p
1139 real :: tr, pv, pd, el, expo
1140 !
1141 p=pk**cpor
1142 tr=ttp/t
1143 pv=psatb*(tr**xa)*exp(xb*(1.-tr))
1144 pd=p-pv
1145 if(pd.gt.0.) then
1146 el=hvap+dldt*(t-ttp)
1147 expo=el*eps*pv/(cp*t*pd)
1148 expo = min(expo,100.0)
1149 fthex=t*pd**(-rocp)*exp(expo)
1150 else
1151 fthex=0.
1152 endif
1153 !
1154 return
1155 end function
1156 !-------------------------------------------------------------------------------
1157 !
1158 !-------------------------------------------------------------------------------
1159 function ftmaxg(tg,the,pk,qma)
1160 !-------------------------------------------------------------------------------
1161 implicit none
1162 !-------------------------------------------------------------------------------
1163 real :: ftmaxg
1164 real :: tg,the,pk,t,p,tr,pv,pd,el,expo,thet,dthet,terr,qma
1165 !
1166 t=tg
1167 p=pk**cpor
1168 tr=ttp/t
1169 pv=psatb*(tr**xa)*exp(xb*(1.-tr))
1170 pd=p-pv
1171 el=hvap+dldt*(t-ttp)
1172 expo=el*eps*pv/(cp*t*pd)
1173 thet=t*pd**(-rocp)*exp(expo)
1174 dthet=thet/t*(1.+expo*(dldt*t/el+el*p/(rv*t*pd)))
1175 terr=(thet-the)/dthet
1176 t=t-terr
1177 !
1178 do while(abs(terr).gt.terrm)
1179 tr=ttp/t
1180 pv=psatb*(tr**xa)*exp(xb*(1.-tr))
1181 pd=p-pv
1182 el=hvap+dldt*(t-ttp)
1183 expo=el*eps*pv/(cp*t*pd)
1184 thet=t*pd**(-rocp)*exp(expo)
1185 dthet=thet/t*(1.+expo*(dldt*t/el+el*p/(rv*t*pd)))
1186 terr=(thet-the)/dthet
1187 t=t-terr
1188 enddo
1189 ftmaxg=t
1190 tr=ttp/t
1191 pv=psatb*(tr**xa)*exp(xb*(1.-tr))
1192 pd=p-pv
1193 qma=eps*pv/(pd+eps*pv)
1194 !
1195 return
1196 end function
1197 !-------------------------------------------------------------------------------
1198 !
1199 !-------------------------------------------------------------------------------
1200 function fpvsx(t)
1201 !-------------------------------------------------------------------------------
1202 implicit none
1203 !-------------------------------------------------------------------------------
1204 real :: fpvsx
1205 !
1206 real :: t
1207 real :: tr
1208 !
1209 tr=ttp/t
1210 if(t.ge.ttp) then
1211 fpvsx=psatk*(tr**xa)*exp(xb*(1.-tr))
1212 else
1213 fpvsx=psatk*(tr**xai)*exp(xbi*(1.-tr))
1214 endif
1215 !
1216 return
1217 end function
1218 !-------------------------------------------------------------------------------
1219 end module module_shcu_grims