module_bl_mrf.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_bl_mrf
4
5 CONTAINS
6
7 !-------------------------------------------------------------------
8 SUBROUTINE MRF(U3D,V3D,TH3D,T3D,QV3D,QC3D,P3D,PI3D, &
9 RUBLTEN,RVBLTEN,RTHBLTEN, &
10 RQVBLTEN,RQCBLTEN, &
11 CP,G,ROVCP,R,ROVG, &
12 dz8w,z,XLV,RV,PSFC, &
13 p1000mb, &
14 ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
15 XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
16 DT,DTMIN,KPBL2D, &
17 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
18 flag_qi, &
19 ids,ide, jds,jde, kds,kde, &
20 ims,ime, jms,jme, kms,kme, &
21 its,ite, jts,jte, kts,kte, &
22 ! Optional
23 QI3D,RQIBLTEN, &
24 regime )
25 !-------------------------------------------------------------------
26 IMPLICIT NONE
27 !-------------------------------------------------------------------
28 !-- U3D 3D u-velocity interpolated to theta points (m/s)
29 !-- V3D 3D v-velocity interpolated to theta points (m/s)
30 !-- TH3D 3D potential temperature (K)
31 !-- T3D temperature (K)
32 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
33 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
34 !-- QI3D 3D ice mixing ratio (Kg/Kg)
35 !-- P3D 3D pressure (Pa)
36 !-- PI3D 3D exner function (dimensionless)
37 !-- rr3D 3D dry air density (kg/m^3)
38 !-- RUBLTEN U tendency due to
39 ! PBL parameterization (m/s^2)
40 !-- RVBLTEN V tendency due to
41 ! PBL parameterization (m/s^2)
42 !-- RTHBLTEN Theta tendency due to
43 ! PBL parameterization (K/s)
44 !-- RQVBLTEN Qv tendency due to
45 ! PBL parameterization (kg/kg/s)
46 !-- RQCBLTEN Qc tendency due to
47 ! PBL parameterization (kg/kg/s)
48 !-- RQIBLTEN Qi tendency due to
49 ! PBL parameterization (kg/kg/s)
50 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
51 !-- G acceleration due to gravity (m/s^2)
52 !-- ROVCP R/CP
53 !-- R gas constant for dry air (J/kg/K)
54 !-- ROVG R/G
55 !-- dz8w dz between full levels (m)
56 !-- z height above sea level (m)
57 !-- XLV latent heat of vaporization (J/kg)
58 !-- RV gas constant for water vapor (J/kg/K)
59 !-- PSFC pressure at the surface (Pa)
60 !-- ZNT roughness length (m)
61 !-- UST u* in similarity theory (m/s)
62 !-- ZOL z/L height over Monin-Obukhov length
63 !-- HOL PBL height over Monin-Obukhov length
64 !-- PBL PBL height (m)
65 !-- REGIME flag indicating PBL regime (stable, unstable, etc.)
66 !-- PSIM similarity stability function for momentum
67 !-- PSIH similarity stability function for heat
68 !-- XLAND land mask (1 for land, 2 for water)
69 !-- HFX upward heat flux at the surface (W/m^2)
70 !-- QFX upward moisture flux at the surface (kg/m^2/s)
71 !-- TSK surface temperature (K)
72 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
73 !-- WSPD wind speed at lowest model level (m/s)
74 !-- BR bulk Richardson number in surface layer
75 !-- DT time step (s)
76 !-- DTMIN time step (minute)
77 !-- rvovrd R_v divided by R_d (dimensionless)
78 !-- SVP1 constant for saturation vapor pressure (kPa)
79 !-- SVP2 constant for saturation vapor pressure (dimensionless)
80 !-- SVP3 constant for saturation vapor pressure (K)
81 !-- SVPT0 constant for saturation vapor pressure (K)
82 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
83 !-- EP2 constant for specific humidity calculation
84 !-- KARMAN Von Karman constant
85 !-- EOMEG angular velocity of earth's rotation (rad/s)
86 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
87 !-- ids start index for i in domain
88 !-- ide end index for i in domain
89 !-- jds start index for j in domain
90 !-- jde end index for j in domain
91 !-- kds start index for k in domain
92 !-- kde end index for k in domain
93 !-- ims start index for i in memory
94 !-- ime end index for i in memory
95 !-- jms start index for j in memory
96 !-- jme end index for j in memory
97 !-- kms start index for k in memory
98 !-- kme end index for k in memory
99 !-- its start index for i in tile
100 !-- ite end index for i in tile
101 !-- jts start index for j in tile
102 !-- jte end index for j in tile
103 !-- kts start index for k in tile
104 !-- kte end index for k in tile
105 !-------------------------------------------------------------------
106
107 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
108 ims,ime, jms,jme, kms,kme, &
109 its,ite, jts,jte, kts,kte
110
111 !
112 REAL, INTENT(IN ) :: P1000mb
113 REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
114 ROVG,R,XLV,RV
115
116 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
117 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
118
119 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
120 INTENT(IN ) :: QV3D, &
121 QC3D, &
122 P3D, &
123 PI3D, &
124 TH3D, &
125 T3D, &
126 dz8w, &
127 z
128 !
129 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
130 INTENT(INOUT) :: RUBLTEN, &
131 RVBLTEN, &
132 RTHBLTEN, &
133 RQVBLTEN, &
134 RQCBLTEN
135
136 REAL, DIMENSION( ims:ime, jms:jme ) , &
137 INTENT(IN ) :: XLAND, &
138 HFX, &
139 QFX
140
141 REAL, DIMENSION( ims:ime, jms:jme ) , &
142 INTENT(INOUT) :: HOL, &
143 UST, &
144 PBL, &
145 ZNT
146
147 LOGICAL, INTENT(IN) :: FLAG_QI
148 !
149 !m The following 5 variables are changed to memory size from tile size--
150 !
151 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
152 PSIM, &
153 PSIH
154
155 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: &
156 WSPD
157
158 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: &
159 GZ1OZ0, &
160 BR
161
162 REAL, DIMENSION( ims:ime, jms:jme ) , &
163 INTENT(IN ) :: PSFC
164
165 REAL, DIMENSION( ims:ime, jms:jme ) , &
166 INTENT(IN ) :: TSK
167
168 REAL, DIMENSION( ims:ime, jms:jme ) , &
169 INTENT(INOUT) :: ZOL
170
171 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
172 INTENT(OUT ) :: KPBL2D
173
174 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
175 INTENT(IN ) :: U3D, &
176 V3D
177 !
178 ! Optional
179 !
180 REAL, DIMENSION( ims:ime, jms:jme ) , &
181 OPTIONAL , &
182 INTENT(INOUT) :: REGIME
183
184 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
185 INTENT(INOUT) :: RQIBLTEN
186
187 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
188 OPTIONAL , &
189 INTENT(IN ) :: QI3D
190
191 ! LOCAL VARS
192 REAL, DIMENSION( its:ite, kts:kte ) :: dz8w2d, &
193 z2d
194
195
196 INTEGER :: I,J,K,NK
197
198 !
199 DO J=jts,jte
200 DO k=kts,kte
201 NK=kme-k
202 DO i=its,ite
203 dz8w2d(I,K) = dz8w(i,NK,j)
204 z2d(I,K) = z(i,NK,j)
205 ENDDO
206 ENDDO
207
208
209 CALL MRF2D(J,U3D(ims,kms,j),V3D(ims,kms,j),T3D(ims,kms,j), &
210 QV3D(ims,kms,j),QC3D(ims,kms,j), &
211 P3D(ims,kms,j),RUBLTEN(ims,kms,j),RVBLTEN(ims,kms,j),&
212 RTHBLTEN(ims,kms,j),RQVBLTEN(ims,kms,j), &
213 RQCBLTEN(ims,kms,j), &
214 p1000mb, &
215 CP,G,ROVCP,R,ROVG, &
216 dz8w2d,z2d,XLV,Rv, &
217 PSFC(ims,j),ZNT(ims,j), &
218 UST(ims,j),ZOL(ims,j), &
219 HOL(ims,j),PBL(ims,j),PSIM(ims,j), &
220 PSIH(ims,j),XLAND(ims,j),HFX(ims,j),QFX(ims,j), &
221 TSK(ims,j),GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j), &
222 DT,DTMIN,KPBL2D(ims,j), &
223 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT, &
224 flag_qi, &
225 ids,ide, jds,jde, kds,kde, &
226 ims,ime, jms,jme, kms,kme, &
227 its,ite, jts,jte, kts,kte, &
228 !optional
229 QI2DTEN=RQIBLTEN(ims,kms,j), &
230 REGIME=REGIME(ims,j),QI2D=QI3D(ims,kms,j) )
231
232
233 DO k=kts,kte
234 DO i=its,ite
235 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/PI3D(I,K,J)
236 ENDDO
237 ENDDO
238
239 ENDDO
240
241 END SUBROUTINE MRF
242
243 !-------------------------------------------------------------------
244 SUBROUTINE MRF2D(J,U2D,V2D,T2D,QV2D,QC2D, P2D, &
245 U2DTEN,V2DTEN,T2DTEN, &
246 QV2DTEN,QC2DTEN, &
247 p1000mb, &
248 CP,G,ROVCP,R,ROVG, &
249 dz8w2d,z2d,XLV,RV,PSFCPA, &
250 ZNT,UST,ZOL,HOL,PBL,PSIM,PSIH, &
251 XLAND,HFX,QFX,TSK,GZ1OZ0,WSPD,BR, &
252 DT,DTMIN,KPBL1D, &
253 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,&
254 flag_qi, &
255 ids,ide, jds,jde, kds,kde, &
256 ims,ime, jms,jme, kms,kme, &
257 its,ite, jts,jte, kts,kte, &
258 ! optional
259 regime, qi2d, QI2DTEN )
260 !-------------------------------------------------------------------
261 IMPLICIT NONE
262 !-------------------------------------------------------------------
263 ! BASED ON THE "COUNTERGRADIENT" TRANSPORT TERM OF TROEN
264 ! AND MAHRT (1986) FOR THE UNSTABLE PBL.
265 ! THIS ROUTINE USES AN IMPLICIT APPROACH FOR VERTICAL FLUX
266 ! DIVERGENCE AND DOES NOT REQUIRE "MITER" TIMESTEPS.
267 ! IT INCLUDES VERTICAL DIFFUSION IN THE STABLE ATMOSPHERE
268 ! AND MOIST VERTICAL DIFFUSION IN CLOUDS.
269 ! SURFACE FLUXES CALCULATED AS IN HIRPBL.
270 ! 5-LAYER SOIL MODEL OPTION REQUIRED IN SLAB DUE TO LONG TIMESTEP
271 !
272 ! CODED BY SONG-YOU HONG (NCEP), IMPLEMENTED BY JIMY DUDHIA (NCAR)
273 ! FALL 1996
274 !
275 ! REFERENCES:
276 !
277 ! HONG AND PAN (1996), MON. WEA. REV.
278 ! TROEN AND MAHRT (1986), BOUNDARY LAYER MET.
279 !
280 ! CHANGES:
281 ! INCREASE RLAM FROM 30 TO 150, AND CHANGE FREE ATMOSPHERE
282 ! STABILITY FUNCTION TO INCREASE VERTICAL DIFFUSION
283 ! (HONG, JUNE 1997)
284 !
285 ! PUT LOWER LIMIT ON PSI FOR STABLE CONDITIONS. THIS WILL
286 ! PREVENT FLUXES FROM BECOMING TOO SMALL (DUDHIA, OCTOBER 1997)
287 !
288 ! CORRECTION TO REGIME CALCULATION. THIS WILL ALLOW POINTS IN
289 ! REGIME 4 MUCH MORE FREQUENTLY GIVING LARGER SURFACE FLUXES
290 ! REGIME 3 NO LONGER USES HOL < 1.5 OR THVX LAPSE-RATE CHECK
291 ! IN MRF SCHEME. THIS WILL MAKE REGIME 3 MUCH LESS FREQUENT.
292 !
293 ! ADD SURFACE PRESSURE, PS(I), ARRAY FOR EFFICIENCY
294 !
295 ! FIX FOR PROBLEM WITH THIN LAYERS AND HIGH ROUGHNESS
296 !
297 ! CHARNOCK CONSTANT NOW COMES FROM NAMELIST (DEFAULT SAME)
298 !
299 !-------------------------------------------------------------------
300
301 REAL RLAM,PRMIN,PRMAX,XKZMIN,XKZMAX,RIMIN,BRCR, &
302 CFAC,PFAC,SFCFRAC,CKZ,ZFMIN,APHI5,APHI16,GAMCRT, &
303 GAMCRQ,XKA,PRT
304
305 PARAMETER (RLAM=150.,PRMIN=0.5,PRMAX=4.)
306 PARAMETER (XKZMIN=0.01,XKZMAX=1000.,RIMIN=-100.)
307 PARAMETER (BRCR=0.5,CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
308 PARAMETER (CKZ=0.001,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
309 PARAMETER (GAMCRT=3.,GAMCRQ=2.E-3)
310 PARAMETER (XKA=2.4E-5)
311 PARAMETER (PRT=1.)
312 !
313 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
314 ims,ime, jms,jme, kms,kme, &
315 its,ite, jts,jte, kts,kte, &
316 J
317 !
318 LOGICAL, INTENT(IN) :: FLAG_QI
319 !
320 REAL, INTENT(IN ) :: P1000mb
321 REAL, INTENT(IN ) :: DT,DTMIN,CP,G,ROVCP, &
322 ROVG,R,XLV,RV
323
324 REAL, INTENT(IN ) :: SVP1,SVP2,SVP3,SVPT0
325 REAL, INTENT(IN ) :: EP1,EP2,KARMAN,EOMEG,STBOLT
326
327 REAL, DIMENSION( ims:ime, kms:kme ) , &
328 INTENT(IN ) :: QV2D, &
329 QC2D, &
330 P2D, &
331 T2D
332 !
333 REAL, DIMENSION( ims:ime, kms:kme ) , &
334 INTENT(INOUT) :: U2DTEN, &
335 V2DTEN, &
336 T2DTEN, &
337 QV2DTEN, &
338 QC2DTEN
339
340 REAL, DIMENSION( ims:ime ) , &
341 INTENT(INOUT) :: HOL, &
342 UST, &
343 PBL, &
344 ZNT
345
346 REAL, DIMENSION( ims:ime ) , &
347 INTENT(IN ) :: XLAND, &
348 HFX, &
349 QFX
350 !
351 !m The following 5 are changed to memory size---
352 !
353 REAL, DIMENSION( ims:ime ), INTENT(IN ) :: PSIM, &
354 PSIH
355
356 REAL, DIMENSION( ims:ime ), INTENT(INOUT) :: WSPD
357
358 REAL, DIMENSION( ims:ime ), INTENT(IN ) :: GZ1OZ0, &
359 BR
360
361 REAL, DIMENSION( ims:ime ) , &
362 INTENT(IN ) :: PSFCPA
363
364 REAL, DIMENSION( ims:ime ) , &
365 INTENT(IN ) :: TSK
366
367 REAL, DIMENSION( ims:ime ) , &
368 INTENT(INOUT) :: ZOL
369
370 INTEGER, DIMENSION( ims:ime ) , &
371 INTENT(OUT ) :: KPBL1D
372
373 REAL, DIMENSION( ims:ime, kms:kme ) , &
374 INTENT(IN ) :: U2D, &
375 V2D
376
377 ! MODULE-LOCAL VARIABLES (DEFINED IN SUBROUTINE MRF)
378 !
379 REAL, DIMENSION( its:ite, kts:kte ) , &
380 INTENT(IN) :: dz8w2d, &
381 z2d
382 !
383 !
384 ! Optional
385 !
386 REAL, DIMENSION( ims:ime ) , &
387 OPTIONAL , &
388 INTENT(INOUT) :: REGIME
389
390 REAL, DIMENSION( ims:ime, kms:kme ) , &
391 OPTIONAL , &
392 INTENT(IN ) :: QI2D
393
394 REAL, DIMENSION( ims:ime, kms:kme ) , &
395 OPTIONAL , &
396 INTENT(INOUT) :: QI2DTEN
397
398 ! LOCAL VARS
399
400 REAL, DIMENSION( its:ite, kts:kte+1 ) :: ZQ
401
402 REAL, DIMENSION( its:ite, kts:kte ) :: &
403 UX,VX,QX, &
404 QCX,THX,THVX, &
405 DZQ,DZA, &
406 TTNP,QTNP, &
407 QCTNP,ZA, &
408 UXS,VXS, &
409 THXS,QXS, &
410 QCXS,QIX, &
411 QITNP,QIXS, &
412 UTNP,VTNP
413 !
414 REAL, DIMENSION( its:ite ) :: QIXSV,RHOX, &
415 WSPD1,GOVRTH, &
416 PBL0,THXSV, &
417 UXSV,VXSV, &
418 QXSV,QCXSV, &
419 QGH,TGDSA,PS
420
421 INTEGER :: ILXM,JLXM,KL, &
422 KLM,KLP1,KLPBL
423 !
424 INTEGER, DIMENSION( its:ite ) :: KPBL,KPBL0
425 !
426 REAL, DIMENSION( its:ite, kts:kte ) :: SCR3,SCR4
427 !
428 REAL, DIMENSION( its:ite ) :: DUM1, &
429 XKZMKL
430 !
431 REAL, DIMENSION( its:ite ) :: ZL1,THERMAL, &
432 WSCALE,HGAMT, &
433 HGAMQ,BRDN, &
434 BRUP,PHIM, &
435 PHIH,CPM, &
436 DUSFC,DVSFC, &
437 DTSFC,DQSFC
438
439 !
440 REAL, DIMENSION( its:ite, kts:kte ) :: XKZM,XKZH, &
441 A1,A2, &
442 AD,AU, &
443 TX
444 !
445 REAL, DIMENSION( its:ite, kts:kte ) :: AL
446 !
447 LOGICAL, DIMENSION( its:ite ) :: PBLFLG, &
448 SFCFLG, &
449 STABLE
450 !
451 REAL, DIMENSION( its:ite ) :: THGB
452
453 INTEGER :: N,I,K,KK,L,NZOL,IMVDIF
454
455 INTEGER :: JBGN,JEND,IBGN,IEND,NK
456
457 REAL :: ZOLN,X,Y,CONT,CONQ,CONW,PL,THCON,TVCON,E1,DTSTEP
458 REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL
459 REAL :: DTTHX,PSIX,DTG,PSIQ,USTM
460 REAL :: DT4,RDT,SPDK2,FM,FH,HOL1,GAMFAC,VPERT,PRNUM
461 REAL :: ZFAC,XKZO,SS,RI,QMEAN,TMEAN,ALPH,CHI,ZK,RL2,DK,SRI
462 REAL :: BRINT,DTODSD,DSIG,RDZ,DSDZT,DSDZQ,DSDZ2,TTEND,QTEND
463 REAL :: UTEND,VTEND,QCTEND,QITEND,TGC,DTODSU
464
465 !----------------------------------------------------------------------
466
467 KLPBL=1
468 KL=kte
469 ILXM=ite-1
470 JLXM=jte-1
471 KLM=kte-1
472 KLP1=kte+1
473 !
474 CONT=1000.*CP/G
475 CONQ=1000.*XLV/G
476 CONW=1000./G
477
478 !-- IMVDIF imvdif=1 for moist adiabat vertical diffusion
479
480 IMVDIF=1
481
482 ! DO i=its,ite
483 !!PS PSFC cmb
484 ! PSFC(I)=PSFCPA(I)/1000.
485 ! ENDDO
486
487
488 !
489 !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
490 !
491 DO 5 I=its,ite
492 TGDSA(I)=TSK(I)
493 ! PS PSFC cmb
494 PS(I)=PSFCPA(I)/1000.
495 ! THGB(I)=TSK(I)*(100./PS(I))**ROVCP
496 THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
497 5 CONTINUE
498 !
499 !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
500 ! T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
501 !
502 ! *** NOTE ***
503 ! THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
504 ! SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
505 ! TENDENCIES.
506 !
507 DO 24 K=kts,kte
508 NK=kme-K
509 DO 24 I=its,ite
510 UX(I,K)=U2D(I,NK)
511 VX(I,K)=V2D(I,NK)
512 24 CONTINUE
513 !
514 !.....SCR3(I,K) STORE TEMPERATURE,
515 ! SCR4(I,K) STORE VIRTUAL TEMPERATURE.
516 !
517 DO 30 K=kts,kte
518 NK=kme-K
519 DO 30 I=its,ite
520 ! PL cmb
521 PL=P2D(I,NK)/1000.
522 SCR3(I,K)=T2D(I,NK)
523 ! THCON=(100./PL)**ROVCP
524 THCON=(P1000mb/(PL*1000.))**ROVCP
525 THX(I,K)=SCR3(I,K)*THCON
526 TX(I,K)=SCR3(I,K)
527 SCR4(I,K)=SCR3(I,K)
528 THVX(I,K)=THX(I,K)
529 QX(I,K)=0.
530 30 CONTINUE
531 !
532 DO I=its,ite
533 QGH(i)=0.
534 CPM(i)=CP
535 ENDDO
536 !
537 ! IF(IDRY.EQ.1)GOTO 80
538 DO 50 K=kts,kte
539 NK=kme-K
540 DO 50 I=its,ite
541 QX(I,K)=QV2D(I,NK)
542 TVCON=(1.+EP1*QX(I,K))
543 THVX(I,K)=THX(I,K)*TVCON
544 SCR4(I,K)=SCR3(I,K)*TVCON
545 50 CONTINUE
546 !
547 DO 60 I=its,ite
548 E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
549 QGH(I)=EP2*E1/(PS(I)-E1)
550 CPM(I)=CP*(1.+0.8*QX(I,KL))
551 60 CONTINUE
552 !
553 ! IF(IMOIST.EQ.1)GOTO 80
554 DO 70 K=kts,kte
555 NK=kme-K
556 DO 70 I=its,ite
557 QCX(I,K)=QC2D(I,NK)
558 70 CONTINUE
559
560 IF (flag_QI .AND. PRESENT( QI2D ) ) THEN
561 DO K=kts,kte
562 NK=kme-K
563 DO I=its,ite
564 QIX(I,K)=QI2D(I,NK)
565 ENDDO
566 ENDDO
567 ELSE
568 DO K=kts,kte
569 NK=kme-K
570 DO I=its,ite
571 QIX(I,K)=0.
572 ENDDO
573 ENDDO
574 ENDIF
575
576 80 CONTINUE
577
578 !
579 !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
580 ! LEVEL, AND THE LAYER THICKNESSES.
581 !
582 DO 90 I=its,ite
583 ZQ(I,KLP1)=0.
584 RHOX(I)=PS(I)*1000./(R*SCR4(I,KL))
585 90 CONTINUE
586 !
587 DO 110 KK=kts,kte
588 K=kme-KK
589 DO 100 I=its,ite
590 DUM1(I)=ZQ(I,K+1)
591 100 CONTINUE
592 !
593 DO 110 I=its,ite
594 ZQ(I,K)=dz8w2d(I,K)+DUM1(I)
595 110 CONTINUE
596 !
597 DO 120 K=kts,kte
598 DO 120 I=its,ite
599 ZA(I,K)=0.5*(ZQ(I,K)+ZQ(I,K+1))
600 DZQ(I,K)=ZQ(I,K)-ZQ(I,K+1)
601 120 CONTINUE
602 !
603 DO 130 K=kts,kte-1
604 DO 130 I=its,ite
605 DZA(I,K)=ZA(I,K)-ZA(I,K+1)
606 130 CONTINUE
607
608 DTSTEP=DT
609 !
610 DO 160 I=its,ite
611 GOVRTH(I)=G/THX(I,KL)
612 160 CONTINUE
613 !
614 !-----INITIALIZE VERTICAL TENDENCIES AND
615 !
616 DO I=its,ite
617 DO K=kts,kte
618 UTNP(i,k)=0.
619 VTNP(i,k)=0.
620 TTNP(i,k)=0.
621 ENDDO
622 ENDDO
623 !
624 ! IF(IDRY.EQ.1)GOTO 250
625 DO 230 K=kts,kte
626 DO 230 I=its,ite
627 QTNP(I,K)=0.
628 230 CONTINUE
629 !
630 ! IF(IMOIST.EQ.1)GOTO 250
631 DO 240 K=kts,kte
632 DO 240 I=its,ite
633 QCTNP(I,K)=0.
634 QITNP(I,K)=0.
635 240 CONTINUE
636
637 250 CONTINUE
638 !
639 !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
640 ! AKB(1976), EQ(12).
641
642 ! DO 260 I=its,ite
643 ! GZ1OZ0(I)=ALOG(ZA(I,KL)/ZNT(I))
644 ! IF((XLAND(I)-1.5).GE.0)THEN
645 ! ZL=ZNT(I)
646 ! ELSE
647 ! ZL=0.01
648 ! ENDIF
649 ! WSPD(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))
650 ! TSKV=THGB(I)*(1.+EP1*QGH(I)*MAVAIL(I))
651 ! DTHVDZ=(THVX(I,KL)-TSKV)
652 ! IF(-DTHVDZ.GE.0)THEN
653 ! DTHVM=-DTHVDZ
654 ! ELSE
655 ! DTHVM=0.
656 ! ENDIF
657 ! VCONV=VCONVC*SQRT(DTHVM)
658 ! WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV)
659 ! WSPD(I)=AMAX1(WSPD(I),1.)
660 ! BR(I)=GOVRTH(I)*ZA(I,KL)*DTHVDZ/(WSPD(I)*WSPD(I))
661 ! 260 CONTINUE
662
663 !!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
664 !!
665 !! THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
666 !! AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
667 !!
668 !! CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
669 !!
670 !! 1. BR .GE. 0.2;
671 !! REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
672 !!
673 !! 2. BR .LT. 0.2 .AND. BR .GT. 0.0;
674 !! REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
675 !! (REGIME=2),
676 !!
677 !! 3. BR .EQ. 0.0
678 !! REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
679 !!
680 !! 4. BR .LT. 0.0
681 !! REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
682 !!
683 !!-----
684 !
685 ! DO 320 I=its,ite
686 !!----
687 !!-- REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
688 !!-- IF(BR(I).LT.0..AND.HOL(I).GT.1.5)GOTO 310
689 !
690 ! IF(BR(I).LT.0.)GOTO 310
691 !!
692 !!-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
693 !!
694 ! IF(BR(I).LT.0.2)GOTO 270
695 ! REGIME(I)=1.
696 ! PSIM(I)=-10.*GZ1OZ0(I)
697 !! LOWER LIMIT ON PSI IN STABLE CONDITIONS
698 ! PSIM(I)=AMAX1(PSIM(I),-10.)
699 ! PSIH(I)=PSIM(I)
700 ! HOL(I)=0.0
701 ! PBL(I)=0.0
702 ! GOTO 320
703 !!
704 !!-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
705 !!
706 ! 270 IF(BR(I).EQ.0.0)GOTO 280
707 ! REGIME(I)=2.
708 ! PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))
709 !! LOWER LIMIT ON PSI IN STABLE CONDITIONS
710 ! PSIM(I)=AMAX1(PSIM(I),-10.)
711 !!.....AKB(1976), EQ(16).
712 ! PSIH(I)=PSIM(I)
713 ! HOL(I)=0.0
714 ! PBL(I)=0.0
715 ! GOTO 320
716 !!
717 !!-----CLASS 3; FORCED CONVECTION:
718 !!
719 ! 280 REGIME(I)=3.
720 ! PSIM(I)=0.0
721 ! PSIH(I)=PSIM(I)
722 !
723 !! special use kte instead of kme
724 !
725 ! DO 290 KK=kts,kte-1
726 ! K=kte-KK
727 ! IF(THVX(I,K).GT.THVX(I,KL))GOTO 300
728 ! 290 CONTINUE
729 ! STOP 290
730 ! 300 PBL(I)=ZQ(I,K+1)
731 ! IF(UST(I).LT.0.01)THEN
732 ! ZOL(I)=BR(I)*GZ1OZ0(I)
733 ! ELSE
734 ! ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I))
735 ! ENDIF
736 ! HOL(I)=-ZOL(I)*PBL(I)/ZA(I,KL)
737 ! GOTO 320
738 !
739 !!-----CLASS 4; FREE CONVECTION:
740 !
741 !! 310 IF(THVX(I,KLM).GT.THVX(I,KL))GOTO 280
742 !
743 ! 310 CONTINUE
744 ! REGIME(I)=4.
745 ! IF(UST(I).LT.0.01)THEN
746 ! ZOL(I)=BR(I)*GZ1OZ0(I)
747 ! ELSE
748 ! ZOL(I)=KARMAN*GOVRTH(I)*ZA(I,KL)*MOL(I,J)/(UST(I)*UST(I))
749 ! ENDIF
750 ! ZOL(I)=AMIN1(ZOL(I),0.)
751 ! ZOL(I)=AMAX1(ZOL(I),-9.9999)
752 ! NZOL=INT(-ZOL(I)*100.)
753 ! RZOL=-ZOL(I)*100.-NZOL
754 ! PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
755 ! PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
756 !!---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
757 !!--- THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
758 ! PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
759 ! PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
760 ! 320 CONTINUE
761
762 !-----COMPUTE THE FRICTIONAL VELOCITY:
763 ! ZA(1982) EQS(2.60),(2.61).
764
765 DO 330 I=its,ite
766 DTG=THX(I,KL)-THGB(I)
767 PSIX=GZ1OZ0(I)-PSIM(I)
768 IF((XLAND(I)-1.5).GE.0)THEN
769 ZL=ZNT(I)
770 ELSE
771 ZL=0.01
772 ENDIF
773 PSIQ=ALOG(KARMAN*UST(I)*ZA(I,KL)/XKA+ZA(I,KL)/ZL)-PSIH(I)
774 UST(I)=KARMAN*WSPD(I)/PSIX
775 !
776 USTM=AMAX1(UST(I),0.1)
777 IF((XLAND(I)-1.5).GE.0)THEN
778 UST(I)=UST(I)
779 ELSE
780 UST(I)=USTM
781 ENDIF
782 ! MOL(I,J)=KARMAN*DTG/(GZ1OZ0(I)-PSIH(I))/PRT
783 330 CONTINUE
784 !
785 DO 420 I=its,ite
786 WSPD1(I)=SQRT(UX(I,KL)*UX(I,KL)+VX(I,KL)*VX(I,KL))+1.E-9
787 420 CONTINUE
788 !
789 !---- COMPUTE VERTICAL DIFFUSION
790 !
791 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
792 ! COMPUTE PRELIMINARY VARIABLES
793 !
794 !
795 DT4=2.*DTSTEP
796 RDT=1./DT4
797 !
798 DO I=its,ite
799 HGAMT(I)=0.
800 HGAMQ(I)=0.
801 WSCALE(I)=0.
802 KPBL(I)=KL
803 PBL(I)=ZQ(I,KL)
804 KPBL0(I)=KL
805 PBL0(I)=ZQ(I,KL)
806 PBLFLG(I)=.TRUE.
807 SFCFLG(I)=.TRUE.
808 IF(BR(I).GT.0.0)SFCFLG(I)=.FALSE.
809 ZL1(I)=ZA(I,KL)
810 THERMAL(I)=THVX(I,KL)
811 ENDDO
812
813 ! COMPUTE THE FIRST GUESS OF PBL HEIGHT
814
815 DO I=its,ite
816 STABLE(I)=.FALSE.
817 BRUP(I)=BR(I)
818 ENDDO
819 DO K=KLM,KLPBL,-1
820 DO I=its,ite
821 IF(.NOT.STABLE(I))THEN
822 BRDN(I)=BRUP(I)
823 SPDK2=MAX(UX(I,K)**2+VX(I,K)**2,1.)
824 BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
825 KPBL(I)=K
826 STABLE(I)=BRUP(I).GT.BRCR
827 ENDIF
828 ENDDO
829 ENDDO
830 !
831 DO I=its,ite
832 K=KPBL(I)
833 IF(BRDN(I).GE.BRCR)THEN
834 BRINT=0.
835 ELSEIF(BRUP(I).LE.BRCR)THEN
836 BRINT=1.
837 ELSE
838 BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
839 ENDIF
840 PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
841 IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
842 ENDDO
843 !
844 DO I=its,ite
845 FM=GZ1OZ0(I)-PSIM(I)
846 FH=GZ1OZ0(I)-PSIH(I)
847 HOL(I)=MAX(BR(I)*FM*FM/FH,RIMIN)
848 IF(SFCFLG(I))THEN
849 HOL(I)=MIN(HOL(I),-ZFMIN)
850 ELSE
851 HOL(I)=MAX(HOL(I),ZFMIN)
852 ENDIF
853 !
854 HOL1=HOL(I)*PBL(I)/ZL1(I)*SFCFRAC
855 HOL(I)=-HOL(I)*PBL(I)/ZL1(I)
856 IF(SFCFLG(I))THEN
857 PHIM(I)=(1.-APHI16*HOL1)**(-1./4.)
858 PHIH(I)=(1.-APHI16*HOL1)**(-1./2.)
859 ELSE
860 PHIM(I)=(1.+APHI5*HOL1)
861 PHIH(I)=PHIM(I)
862 ENDIF
863 WSCALE(I)=UST(I)/PHIM(I)
864 WSCALE(I)=MIN(WSCALE(I),UST(I)*APHI16)
865 WSCALE(I)=MAX(WSCALE(I),UST(I)/APHI5)
866 ENDDO
867
868 ! COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
869 ! UNDER UNSTABLE CONDITIONS
870
871 DO I=its,ite
872 IF(SFCFLG(I))THEN
873 GAMFAC=CFAC/RHOX(I)/WSCALE(I)
874 HGAMT(I)=MIN(GAMFAC*HFX(I)/CPM(I),GAMCRT)
875 HGAMQ(I)=MIN(GAMFAC*QFX(I),GAMCRQ)
876 IF((XLAND(I)-1.5).GE.0)HGAMQ(I)=0.
877 VPERT=HGAMT(I)+EP1*THX(I,KL)*HGAMQ(I)
878 VPERT=MIN(VPERT,GAMCRT)
879 THERMAL(I)=THERMAL(I)+MAX(VPERT,0.)
880 HGAMT(I)=MAX(HGAMT(I),0.0)
881 HGAMQ(I)=MAX(HGAMQ(I),0.0)
882 ELSE
883 PBLFLG(I)=.FALSE.
884 ENDIF
885 ENDDO
886 !
887 DO I=its,ite
888 IF(PBLFLG(I))THEN
889 KPBL(I)=KL
890 PBL(I)=ZQ(I,KL)
891 ENDIF
892 ENDDO
893 !
894 ! ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
895 !
896 DO I=its,ite
897 IF(PBLFLG(I))THEN
898 STABLE(I)=.FALSE.
899 BRUP(I)=BR(I)
900 ENDIF
901 ENDDO
902 DO K=KLM,KLPBL,-1
903 DO I=its,ite
904 IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
905 BRDN(I)=BRUP(I)
906 SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)
907 BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
908 KPBL(I)=K
909 STABLE(I)=BRUP(I).GT.BRCR
910 ENDIF
911 ENDDO
912 ENDDO
913 !
914 DO I=its,ite
915 IF(PBLFLG(I))THEN
916 K=KPBL(I)
917 IF(BRDN(I).GE.BRCR)THEN
918 BRINT=0.
919 ELSEIF(BRUP(I).LE.BRCR)THEN
920 BRINT=1.
921 ELSE
922 BRINT=(BRCR-BRDN(I))/(BRUP(I)-BRDN(I))
923 ENDIF
924 PBL(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
925 IF(PBL(I).LT.ZQ(I,KPBL(I)+1))KPBL(I)=KPBL(I)+1
926 IF(KPBL(I).LE.1)PBLFLG(I)=.FALSE.
927 ENDIF
928 ENDDO
929 !
930 ! DIAGNOSTIC PBL HEIGHT WITH BRCR EFFECTIVELY ZERO (PBL0)
931 !
932 DO I=its,ite
933 IF(PBLFLG(I))THEN
934 STABLE(I)=.FALSE.
935 BRUP(I)=BR(I)
936 ENDIF
937 ENDDO
938 DO K=KLM,KLPBL,-1
939 DO I=its,ite
940 IF(.NOT.STABLE(I).AND.PBLFLG(I))THEN
941 BRDN(I)=BRUP(I)
942 SPDK2=MAX((UX(I,K)**2+VX(I,K)**2),1.)
943 BRUP(I)=(THVX(I,K)-THERMAL(I))*(G*ZA(I,K)/THVX(I,KL))/SPDK2
944 KPBL0(I)=K
945 STABLE(I)=BRUP(I).GT.0.0
946 ENDIF
947
948 ENDDO
949 ENDDO
950 !
951 DO I=its,ite
952 IF(PBLFLG(I))THEN
953 K=KPBL0(I)
954 IF(BRDN(I).GE.0.0)THEN
955 BRINT=0.
956 ELSEIF(BRUP(I).LE.0.0)THEN
957 BRINT=1.
958 ELSE
959 BRINT=(0.0-BRDN(I))/(BRUP(I)-BRDN(I))
960 ENDIF
961 PBL0(I)=ZA(I,K+1)+BRINT*(ZA(I,K)-ZA(I,K+1))
962 IF(PBL0(I).LT.ZQ(I,KPBL0(I)+1))KPBL0(I)=KPBL0(I)+1
963 IF(KPBL0(I).LE.1)PBLFLG(I)=.FALSE.
964 ENDIF
965 ENDDO
966
967 !
968 ! COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
969 !
970 DO K=kte,KLPBL,-1
971 DO I=its,ite
972 IF(KPBL(I).LT.K)THEN
973 PRNUM=(PHIH(I)/PHIM(I)+CFAC*KARMAN*SFCFRAC)
974 PRNUM=MIN(PRNUM,PRMAX)
975 PRNUM=MAX(PRNUM,PRMIN)
976 ZFAC=MAX((1.-(ZQ(I,K)-ZL1(I))/(PBL(I)-ZL1(I))),ZFMIN)
977 XKZO=CKZ*DZA(I,K-1)
978 XKZM(I,K)=XKZO+WSCALE(I)*KARMAN*ZQ(I,K)*ZFAC**PFAC
979 XKZH(I,K)=XKZM(I,K)/PRNUM
980 XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
981 XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
982 XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
983 XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
984 ENDIF
985 ENDDO
986 ENDDO
987 !
988 ! COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
989 !
990 DO K=kts+1,kte
991 DO I=its,ite
992 XKZO=CKZ*DZA(I,K-1)
993 IF(K.LE.KPBL(I))THEN
994 SS=((UX(I,K-1)-UX(I,K))*(UX(I,K-1)-UX(I,K))+(VX(I,K-1)- &
995 VX(I,K))*(VX(I,K-1)-VX(I,K)))/(DZA(I,K-1)*DZA(I,K-1))+ &
996 1.E-9
997 RI=GOVRTH(I)*(THVX(I,K-1)-THVX(I,K))/(SS*DZA(I,K-1))
998 IF(IMVDIF.EQ.1)THEN
999 IF((QCX(I,K)+QIX(I,K)).GT.0.01E-3.AND.(QCX(I,K-1)+ &
1000 QIX(I,K-1)).GT.0.01E-3)THEN
1001 ! IN CLOUD
1002 QMEAN=0.5*(QX(I,K)+QX(I,K-1))
1003 TMEAN=0.5*(SCR3(I,K)+SCR3(I,K-1))
1004 ALPH=XLV*QMEAN/R/TMEAN
1005 CHI=XLV*XLV*QMEAN/CP/RV/TMEAN/TMEAN
1006 RI=(1.+ALPH)*(RI-G*G/SS/TMEAN/CP*((CHI-ALPH)/(1.+CHI)))
1007 ENDIF
1008 ENDIF
1009 ZK=KARMAN*ZQ(I,K)
1010 RL2=(ZK*RLAM/(RLAM+ZK))**2
1011 DK=RL2*SQRT(SS)
1012 IF(RI.LT.0.)THEN
1013 ! UNSTABLE REGIME
1014 SRI=SQRT(-RI)
1015 XKZM(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.746*SRI))
1016 XKZH(I,K)=XKZO+DK*(1+8.*(-RI)/(1+1.286*SRI))
1017 ELSE
1018 ! STABLE REGIME
1019 XKZH(I,K)=XKZO+DK/(1+5.*RI)**2
1020 PRNUM=1.0+2.1*RI
1021 PRNUM=MIN(PRNUM,PRMAX)
1022 XKZM(I,K)=(XKZH(I,K)-XKZO)*PRNUM+XKZO
1023 ENDIF
1024 !
1025 XKZM(I,K)=MIN(XKZM(I,K),XKZMAX)
1026 XKZM(I,K)=MAX(XKZM(I,K),XKZMIN)
1027 XKZH(I,K)=MIN(XKZH(I,K),XKZMAX)
1028 XKZH(I,K)=MAX(XKZH(I,K),XKZMIN)
1029 ENDIF
1030 !
1031 ENDDO
1032 ENDDO
1033
1034 ! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
1035
1036 DO I=its,ite
1037 DO K=kts,kte
1038 AU(i,k)=0.
1039 AL(i,k)=0.
1040 AD(i,k)=0.
1041 A1(i,k)=0.
1042 A2(i,k)=0.
1043 ENDDO
1044 ENDDO
1045
1046 DO I=its,ite
1047 AD(I,1)=1.
1048 A1(I,1)=SCR3(I,KL)+HFX(I)/(RHOX(I)*CPM(I))/ZQ(I,KL)*DT4
1049 A2(I,1)=QX(I,KL)+QFX(I)/(RHOX(I))/ZQ(I,KL)*DT4
1050 ENDDO
1051 !
1052 DO K=kte,kts+1,-1
1053 KK=kme-K
1054 DO I=its,ite
1055 DTODSD=DT4/dz8w2d(I,K)
1056 DTODSU=DT4/dz8w2d(I,K-1)
1057 DSIG=z2d(I,K)-z2d(I,K-1)
1058 DSIG=-DSIG
1059 RDZ=1./DZA(I,K-1)
1060 IF(PBLFLG(I).AND.KPBL(I).LT.K)THEN
1061 DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP-HGAMT(I)/PBL(I))
1062 DSDZQ=DSIG*XKZH(I,K)*RDZ*(-HGAMQ(I)/PBL(I))
1063 A2(I,KK)=A2(I,KK)+DTODSD*DSDZQ
1064 A2(I,KK+1)=QX(I,K-1)-DTODSU*DSDZQ
1065 ELSE
1066 DSDZT=DSIG*XKZH(I,K)*RDZ*(G/CP)
1067 A2(I,KK+1)=QX(I,K-1)
1068 ENDIF
1069 DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ
1070 AU(I,KK)=-DTODSD*DSDZ2
1071 AL(I,KK)=-DTODSU*DSDZ2
1072 AD(I,KK)=AD(I,KK)-AU(I,KK)
1073 AD(I,KK+1)=1.-AL(I,KK)
1074 A1(I,KK)=A1(I,KK)+DTODSD*DSDZT
1075 A1(I,KK+1)=SCR3(I,K-1)-DTODSU*DSDZT
1076 ENDDO
1077 ENDDO
1078
1079 ! SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
1080
1081 CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1082 its,ite,kts,kte )
1083
1084 ! RECOVER TENDENCIES OF HEAT AND MOISTURE
1085
1086 DO K=kte,kts,-1
1087 KK=kme-K
1088 DO I=its,ite
1089 TTEND=(A1(I,KK)-SCR3(I,K))*RDT
1090 QTEND=(A2(I,KK)-QX(I,K))*RDT
1091 TTNP(I,K)=TTNP(I,K)+TTEND
1092 QTNP(I,K)=QTNP(I,K)+QTEND
1093 ENDDO
1094 ENDDO
1095
1096 ! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
1097
1098 DO I=its,ite
1099 DO K=kts,kte
1100 AU(i,k)=0.
1101 AL(i,k)=0.
1102 AD(i,k)=0.
1103 A1(i,k)=0.
1104 A2(i,k)=0.
1105 ENDDO
1106 ENDDO
1107
1108 DO I=its,ite
1109 AD(I,1)=1.
1110 A1(I,1)=UX(I,KL)-UX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL) &
1111 *DT4*(WSPD1(I)/WSPD(I))**2
1112 A2(I,1)=VX(I,KL)-VX(I,KL)/WSPD1(I)*UST(I)*UST(I)/ZQ(I,KL) &
1113 *DT4*(WSPD1(I)/WSPD(I))**2
1114 ENDDO
1115 !
1116 DO K=kte,kts+1,-1
1117 KK=kme-K
1118 DO I=its,ite
1119 DTODSD=DT4/dz8w2d(I,K)
1120 DTODSU=DT4/dz8w2d(I,K-1)
1121 DSIG=z2d(I,K)-z2d(I,K-1)
1122 DSIG=-DSIG
1123 RDZ=1./DZA(I,K-1)
1124 DSDZ2=DSIG*XKZM(I,K)*RDZ*RDZ
1125 AU(I,KK)=-DTODSD*DSDZ2
1126 AL(I,KK)=-DTODSU*DSDZ2
1127 AD(I,KK)=AD(I,KK)-AU(I,KK)
1128 AD(I,KK+1)=1.-AL(I,KK)
1129 A1(I,KK+1)=UX(I,K-1)
1130 A2(I,KK+1)=VX(I,K-1)
1131 ENDDO
1132 ENDDO
1133
1134 ! SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
1135
1136 CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1137 its,ite,kts,kte )
1138
1139 ! RECOVER TENDENCIES OF MOMENTUM
1140
1141 DO K=kte,kts,-1
1142 KK=kme-K
1143 DO I=its,ite
1144 UTEND=(A1(I,KK)-UX(I,K))*RDT
1145 VTEND=(A2(I,KK)-VX(I,K))*RDT
1146 UTNP(I,K)=UTNP(I,K)+UTEND
1147 VTNP(I,K)=VTNP(I,K)+VTEND
1148 ENDDO
1149 ENDDO
1150
1151 ! COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR CLOUD
1152
1153 DO I=its,ite
1154 DO K=kts,kte
1155 AU(i,k)=0.
1156 AL(i,k)=0.
1157 AD(i,k)=0.
1158 A1(i,k)=0.
1159 A2(i,k)=0.
1160 ENDDO
1161 ENDDO
1162
1163 ! IF(IMOIST.EQ.1)GOTO 690
1164 DO I=its,ite
1165 AD(I,1)=1.
1166 A1(I,1)=QCX(I,KL)
1167 A2(I,1)=QIX(I,KL)
1168 ENDDO
1169 !
1170 DO K=kte,kts+1,-1
1171 KK=kme-K
1172 DO I=its,ite
1173 DTODSD=DT4/dz8w2d(I,K)
1174 DTODSU=DT4/dz8w2d(I,K-1)
1175 DSIG=z2d(I,K)-z2d(I,K-1)
1176 DSIG=-DSIG
1177 RDZ=1./DZA(I,K-1)
1178 A1(I,KK+1)=QCX(I,K-1)
1179 A2(I,KK+1)=QIX(I,K-1)
1180 DSDZ2=DSIG*XKZH(I,K)*RDZ*RDZ
1181 AU(I,KK)=-DTODSD*DSDZ2
1182 AL(I,KK)=-DTODSU*DSDZ2
1183 AD(I,KK)=AD(I,KK)-AU(I,KK)
1184 AD(I,KK+1)=1.-AL(I,KK)
1185 ENDDO
1186 ENDDO
1187
1188 ! SOLVE TRIDIAGONAL PROBLEM FOR CLOUD
1189
1190 CALL TRIDI2(AL,AD,AU,A1,A2,AU,A1,A2, &
1191 its,ite,kts,kte )
1192 !
1193 DO K=kte,kts,-1
1194 KK=kme-K
1195 DO I=its,ite
1196 QCTEND=(A1(I,KK)-QCX(I,K))*RDT
1197 QITEND=(A2(I,KK)-QIX(I,K))*RDT
1198 QCTNP(I,K)=QCTNP(I,K)+QCTEND
1199 QITNP(I,K)=QITNP(I,K)+QITEND
1200 ENDDO
1201 ENDDO
1202 !
1203 !---- END OF VERTICAL DIFFUSION
1204 !
1205 690 CONTINUE
1206 !
1207 !-----CALCULATION OF NEW VALUES DUE TO VERTICAL EXCHANGE PROCESSES IS
1208 ! COMPLETED. THE FINAL STEP IS TO ADD THE TENDENCIES CALCULATED
1209 ! IN HIRPBL TO THOSE OF MM4.
1210
1211 DO 820 K=kts,kte
1212 NK=kme-K
1213 DO 820 I=its,ite
1214 U2DTEN(I,NK)=UTNP(I,K)
1215 V2DTEN(I,NK)=VTNP(I,K)
1216 820 CONTINUE
1217 !
1218 ! IF(J.EQ.1.AND.IN.GT.1)GOTO 860
1219 !SUE JBGN=3
1220 !SUE JEND=JLXM-1
1221
1222 ! change when nest
1223 !SUE JBGN=2
1224 !SUE JEND=JLXM
1225
1226 JBGN=jts
1227 JEND=jte
1228 IBGN=its
1229 IEND=ite
1230
1231 ! IF(J.LT.JBGN.OR.J.GT.JEND)GOTO 860
1232 !SUE IBGN=3
1233 !SUE IEND=ILXM-1
1234
1235 ! change when nest
1236 !SUE IBGN=2
1237 !SUE IEND=ILXM
1238
1239 DO 830 K=kts,kte
1240 NK=kme-K
1241 DO 830 I=IBGN,IEND
1242 T2DTEN(I,NK)=TTNP(I,K)
1243 830 CONTINUE
1244 !
1245 ! IF(IDRY.EQ.1)GOTO 860
1246 DO 840 K=kts,kte
1247 NK=kme-K
1248 DO 840 I=IBGN,IEND
1249 QV2DTEN(I,NK)=QTNP(I,K)
1250 840 CONTINUE
1251
1252 ! IF(IMOIST.EQ.1)GOTO 860
1253 DO 850 K=kts,kte
1254 NK=kme-K
1255 DO 850 I=IBGN,IEND
1256 QC2DTEN(I,NK)=QCTNP(I,K)
1257 850 CONTINUE
1258
1259 IF(flag_QI .AND. PRESENT( QI2DTEN ) ) THEN
1260 DO K=kts,kte
1261 NK=kme-K
1262 DO I=IBGN,IEND
1263 QI2DTEN(I,NK)=QITNP(I,K)
1264 ENDDO
1265 ENDDO
1266 ENDIF
1267
1268 860 CONTINUE
1269 !
1270 !-----APPLY ASSELIN FILTER TO TGD FOR LARGE TIME STEP:
1271 !
1272 ! DO 885 I=its,ite
1273 ! TSK(I)=TSK(I)*(PS(I)/100.)**ROVCP
1274 ! 885 CONTINUE
1275 !
1276 940 CONTINUE
1277 !
1278 ! KPBL IS NEEDED FOR THE FDDA, AND SINCE THERE IS NO LONGER JUST ONE
1279 ! LARGE "J LOOP" IT MUST BE STORED AS (I,J)...
1280 !
1281 ! USE NEW DIAGNOSED PBL DEPTH (CALCULATED WITH brcr=0.0)
1282 ! PBL IS USED FOR OUTPUT AND NEXT-TIME-STEP BELJAARS CALC IN SFCLAY
1283 DO 950 I=its,ite
1284 KPBL1D(I)=KPBL0(I)
1285 PBL(I)=PBL0(I)
1286 950 CONTINUE
1287
1288 END SUBROUTINE MRF2D
1289
1290 !================================================================
1291 SUBROUTINE TRIDI2(CL,CM,CU,R1,R2,AU,A1,A2, &
1292 its,ite,kts,kte )
1293 !----------------------------------------------------------------
1294 IMPLICIT NONE
1295 !----------------------------------------------------------------
1296
1297 INTEGER, INTENT(IN ) :: its,ite, kts,kte
1298
1299 REAL, DIMENSION( its:ite, kts+1:kte+1 ) , &
1300 INTENT(IN ) :: CL
1301
1302 REAL, DIMENSION( its:ite, kts:kte ) , &
1303 INTENT(IN ) :: CM, &
1304 R1, &
1305 R2
1306 REAL, DIMENSION( its:ite, kts:kte ) , &
1307 INTENT(INOUT) :: AU, &
1308 CU, &
1309 A1, &
1310 A2
1311
1312 REAL :: FK
1313 INTEGER :: I,K,L,N
1314
1315 !----------------------------------------------------------------
1316
1317 L=ite
1318 N=kte
1319
1320 DO I=its,L
1321 FK=1./CM(I,1)
1322 AU(I,1)=FK*CU(I,1)
1323 A1(I,1)=FK*R1(I,1)
1324 A2(I,1)=FK*R2(I,1)
1325 ENDDO
1326 DO K=2,N-1
1327 DO I=its,L
1328 FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1329 AU(I,K)=FK*CU(I,K)
1330 A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
1331 A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
1332 ENDDO
1333 ENDDO
1334 DO I=its,L
1335 FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1336 A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
1337 A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
1338
1339 ENDDO
1340 DO K=N-1,kts,-1
1341 DO I=its,L
1342 A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
1343 A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
1344 ENDDO
1345 ENDDO
1346
1347 END SUBROUTINE TRIDI2
1348
1349 !===================================================================
1350 SUBROUTINE mrfinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1351 RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR, &
1352 restart, allowed_to_read , &
1353 ids, ide, jds, jde, kds, kde, &
1354 ims, ime, jms, jme, kms, kme, &
1355 its, ite, jts, jte, kts, kte )
1356 !-------------------------------------------------------------------
1357 IMPLICIT NONE
1358 !-------------------------------------------------------------------
1359 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1360 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1361 ims, ime, jms, jme, kms, kme, &
1362 its, ite, jts, jte, kts, kte
1363 INTEGER , INTENT(IN) :: P_QI,P_FIRST_SCALAR
1364
1365 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
1366 RUBLTEN, &
1367 RVBLTEN, &
1368 RTHBLTEN, &
1369 RQVBLTEN, &
1370 RQCBLTEN, &
1371 RQIBLTEN
1372 INTEGER :: i, j, k, itf, jtf, ktf
1373
1374 jtf=min0(jte,jde-1)
1375 ktf=min0(kte,kde-1)
1376 itf=min0(ite,ide-1)
1377
1378 IF(.not.restart)THEN
1379 DO j=jts,jtf
1380 DO k=kts,ktf
1381 DO i=its,itf
1382 RUBLTEN(i,k,j)=0.
1383 RVBLTEN(i,k,j)=0.
1384 RTHBLTEN(i,k,j)=0.
1385 RQVBLTEN(i,k,j)=0.
1386 RQCBLTEN(i,k,j)=0.
1387 ENDDO
1388 ENDDO
1389 ENDDO
1390 ENDIF
1391
1392 IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
1393 DO j=jts,jtf
1394 DO k=kts,ktf
1395 DO i=its,itf
1396 RQIBLTEN(i,k,j)=0.
1397 ENDDO
1398 ENDDO
1399 ENDDO
1400 ENDIF
1401
1402 END SUBROUTINE mrfinit
1403
1404 !-------------------------------------------------------------------
1405
1406 END MODULE module_bl_mrf
1407