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