| blob | 4768192f6fe3d3c757e430d9419346b404d22242 |
1 !-----------------------------------------------------------------------
2 !
3 MODULE MODULE_BL_QNSEPBL
4 !
5 !-----------------------------------------------------------------------
6 !
7 USE MODULE_MODEL_CONSTANTS
8 !
9 !-----------------------------------------------------------------------
10 !
11 ! REFERENCES: Janjic (2002), NCEP Office Note 437
12 ! Mellor and Yamada (1982), Rev. Geophys. Space Phys.
13 !
14 ! ABSTRACT:
15 ! MYJ UPDATES THE TURBULENT KINETIC ENERGY WITH THE PRODUCTION/
16 ! DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
17 ! (USING AN IMPLICIT FORMULATION) FROM MELLOR-YAMADA
18 ! LEVEL 2.5 AS EXTENDED BY JANJIC. EXCHANGE COEFFICIENTS FOR
19 ! THE SURFACE AND FOR ALL LAYER INTERFACES ARE COMPUTED FROM
20 ! MONIN-OBUKHOV THEORY.
21 ! THE TURBULENT VERTICAL EXCHANGE IS THEN EXECUTED.
22 !
23 !-----------------------------------------------------------------------
24 !
25 INTEGER :: ITRMX=5 ! Iteration count for mixing length computation
26 !
27 ! REAL,PARAMETER :: G=9.81,PI=3.1415926,R_D=287.04,R_V=461.6 &
28 ! & ,VKARMAN=0.4
29 REAL,PARAMETER :: PI=3.1415926,VKARMAN=0.4
30 !
31 !-----------------------------------------------------------------------
32 !*** QNSE MODEL CONSTANTS
33 !-----------------------------------------------------------------------
34 !
35 REAL,PARAMETER :: EPSQ2L=0.01
36 REAL,PARAMETER :: C0=0.55,CEPS=C0**3,BLCKDR=0.0063,CN=0.75 &
37 & ,AM1=8.0,AM2=2.3,AM3=35.0,AH1=1.4,AH2=-0.01 &
38 & ,AH3=1.29,AH4=2.44,AH5=19.8 &
39 & ,ARIMIN=0.127,BM1=2.88,BM2=16.0,BH1=3.6,BH2=16.0 &
40 & ,BH3=720.0,EPSKM=1.E-3
41 !
42 REAL,PARAMETER :: CAPA=R_D/CP
43 REAL,PARAMETER :: RLIVWV=XLS/XLV,ELOCP=2.72E6/CP
44 REAL,PARAMETER :: EPS1=1.E-12,EPS2=0.
45 REAL,PARAMETER :: EPSL=0.32,EPSRU=1.E-7,EPSRS=1.E-7 &
46 & ,EPSTRB=1.E-24
47 REAL,PARAMETER :: EPSA=1.E-8,EPSIT=1.E-4,EPSU2=1.E-4,EPSUST=0.07 &
48 & ,FH=1.01
49 REAL,PARAMETER :: ALPH=0.30,BETA=1./273.,EL0MAX=1000.,EL0MIN=1. &
50 & ,ELFC=0.23*0.5,GAM1=0.2222222222222222222 &
51 & ,PRT=1.
52 REAL,PARAMETER :: A1=0.659888514560862645 &
53 & ,A2x=0.6574209922667784586 &
54 & ,B1=11.87799326209552761 &
55 & ,B2=7.226971804046074028 &
56 & ,C1=0.000830955950095854396
57 REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
58 REAL,PARAMETER :: ELZ0=0.,ESQ=5.0,EXCM=0.001 &
59 & ,FHNEU=0.8,GLKBR=10.,GLKBS=30. &
60 & ,QVISC=2.1E-5,RFC=0.191,RIC=0.505,SMALL=0.35 &
61 & ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5 &
62 & ,USTC=0.7,USTR=0.225,VISC=1.5E-5 &
63 & ,WOLD=0.15,WWST=1.2,ZTMAX=1.,ZTFC=1.,ZTMIN=-5.
64 !
65 REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
66 !
67 REAL,PARAMETER :: BTG=BETA*G,CZIV=SMALL*GLKBS &
68 ! & ,EP_1=R_V/R_D-1.,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
69 & ,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS &
70 & ,RB1=1./B1,RTVISC=1./TVISC,RVISC=1./VISC &
71 & ,ZQRZT=SQSC/SQPR
72 !
73 REAL,PARAMETER :: ADNH= 9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG &
74 & ,ADNM=18.*A1*A1*A2x*(B2-3.*A2x)*BTG &
75 & ,ANMH=-9.*A1*A2x*A2x*BTG*BTG &
76 & ,ANMM=-3.*A1*A2x*(3.*A2x+3.*B2*C1+18.*A1*C1-B2) &
77 & *BTG &
78 & ,BDNH= 3.*A2x*(7.*A1+B2)*BTG &
79 & ,BDNM= 6.*A1*A1 &
80 & ,BEQH= A2x*B1*BTG+3.*A2x*(7.*A1+B2)*BTG &
81 & ,BEQM=-A1*B1*(1.-3.*C1)+6.*A1*A1 &
82 & ,BNMH=-A2x*BTG &
83 & ,BNMM=A1*(1.-3.*C1) &
84 & ,BSHH=9.*A1*A2x*A2x*BTG &
85 & ,BSHM=18.*A1*A1*A2x*C1 &
86 & ,BSMH=-3.*A1*A2x*(3.*A2x+3.*B2*C1+12.*A1*C1-B2) &
87 & *BTG &
88 & ,CESH=A2x &
89 & ,CESM=A1*(1.-3.*C1) &
90 & ,CNV=EP_1*G/BTG &
91 & ,ELFCS=VKARMAN*BTG &
92 & ,FZQ1=RTVISC*QVISC*ZQRZT &
93 & ,FZQ2=RTVISC*QVISC*ZQRZT &
94 & ,FZT1=RVISC *TVISC*SQPR &
95 & ,FZT2=CZIV*GRRS*TVISC*SQPR &
96 & ,FZU1=CZIV*VISC &
97 & ,PIHF=0.5*PI &
98 & ,RFAC=RIC/(FHNEU*RFC*RFC) &
99 & ,RQVISC=1./QVISC &
100 & ,RRIC=1./RIC &
101 & ,USTFC=0.018/G &
102 & ,WNEW=1.-WOLD &
103 & ,WWST2=WWST*WWST
104 !
105 !-----------------------------------------------------------------------
106 !*** FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)**2
107 !-----------------------------------------------------------------------
108 !
109 REAL,PARAMETER :: AEQH=9.*A1*A2x*A2x*B1*BTG*BTG &
110 & +9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG &
111 & ,AEQM=3.*A1*A2x*B1*(3.*A2x+3.*B2*C1+18.*A1*C1-B2)&
112 & *BTG+18.*A1*A1*A2x*(B2-3.*A2x)*BTG
113 !
114 !-----------------------------------------------------------------------
115 !*** FORBIDDEN TURBULENCE AREA
116 !-----------------------------------------------------------------------
117 !
118 REAL,PARAMETER :: REQU=-AEQH/AEQM &
119 & ,EPSGH=1.E-9,EPSGM=REQU*EPSGH
120 !
121 !-----------------------------------------------------------------------
122 !*** NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
123 !-----------------------------------------------------------------------
124 !
125 REAL,PARAMETER :: UBRYL=(18.*REQU*A1*A1*A2x*B2*C1*BTG &
126 & +9.*A1*A2x*A2x*B2*BTG*BTG) &
127 & /(REQU*ADNM+ADNH) &
128 & ,UBRY=(1.+EPSRS)*UBRYL,UBRY3=3.*UBRY
129 !
130 REAL,PARAMETER :: AUBH=27.*A1*A2x*A2x*B2*BTG*BTG-ADNH*UBRY3 &
131 & ,AUBM=54.*A1*A1*A2x*B2*C1*BTG -ADNM*UBRY3 &
132 & ,BUBH=(9.*A1*A2x+3.*A2x*B2)*BTG-BDNH*UBRY3 &
133 & ,BUBM=18.*A1*A1*C1 -BDNM*UBRY3 &
134 & ,CUBR=1. - UBRY3 &
135 & ,RCUBR=1./CUBR
136 !
137 !-----------------------------------------------------------------------
138 !
139 CONTAINS
140 !
141 !----------------------------------------------------------------------
142 SUBROUTINE QNSEPBL(DT,STEPBL,HT,DZ &
143 & ,PMID,PINT,TH,T,EXNER,QV,CWM,U,V,RHO &
144 & ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0,CORF &
145 & ,LOWLYR,XLAND,SICE,SNOW &
146 & ,TKE,EXCH_H,EXCH_M,USTAR,ZNT,EL_MYJ,PBLH,KPBL,CT &
147 & ,AKHS,AKMS,ELFLX,HFX & !JP HFX
148 & ,LM_BL89,WTHV_MF &
149 & ,RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN &
150 & ,IDS,IDE,JDS,JDE,KDS,KDE &
151 & ,IMS,IME,JMS,JME,KMS,KME &
152 & ,ITS,ITE,JTS,JTE,KTS,KTE)
153 !----------------------------------------------------------------------
154 !
155 IMPLICIT NONE
156 !
157 !----------------------------------------------------------------------
158 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
159 & ,IMS,IME,JMS,JME,KMS,KME &
160 & ,ITS,ITE,JTS,JTE,KTS,KTE
161 !
162 INTEGER,INTENT(IN) :: STEPBL
164 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
165 !
166 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: KPBL
167 !
168 REAL,INTENT(IN) :: DT
169 !
170 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,SICE,SNOW &
171 & ,TSK,XLAND,HFX
172 !
173 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DZ &
174 & ,EXNER &
175 & ,PMID,PINT &
176 & ,QV,RHO &
177 & ,T,TH,U,V
179 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN),OPTIONAL :: &
180 & WTHV_MF &
181 & ,LM_BL89
182 !
183 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PBLH
184 !
185 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS
186 !
187 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
188 & ,INTENT(OUT) :: EL_MYJ
189 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
190 & ,INTENT(INOUT) :: RQCBLTEN,RQVBLTEN &
191 & ,RTHBLTEN &
192 & ,RUBLTEN,RVBLTEN
193 !
194 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CT,QSFC,QZ0 &
195 & ,THZ0,USTAR &
196 & ,UZ0,VZ0,ZNT
197 !
198 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
199 & ,INTENT(INOUT) :: EXCH_H,EXCH_M,TKE
200 !
201 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHKLOWQ,ELFLX,CORF
202 !
203 !----------------------------------------------------------------------
204 !***
205 !*** LOCAL VARIABLES
206 !***
207 INTEGER :: I,J,K,KFLIP,LLOW,LMH,LMXL
208 !
209 INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: LPBL
210 !
211 REAL :: AKHS_DENS,AKMS_DENS,APEX,DCDT,DELTAZ,DQDT,DTDIF,DTDT &
212 & ,DTTURBL,DUDT,DVDT,EXNSFC,PSFC,PTOP,QFC1,QLOW,QOLD &
213 & ,RATIOMX,RDTTURBL,RG,RWMSK,SEAMASK,THNEW,THOLD,TX &
214 & ,ULOW,VLOW,WMSK,RLOW
215 !
216 REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,TK,UK,VK,&
217 WTHV_MFK,LM_BL89K
218 !
219 REAL,DIMENSION(KTS:KTE-1) :: AKHK,AKMK,EL,RI,GH,S2
220 !
221 REAL,DIMENSION(KTS:KTE+1) :: ZHK
222 !
223 REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
224 !
225 REAL,DIMENSION(KTS:KTE,ITS:ITE) :: RHOK
226 !
227 REAL,DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: APE,THE
228 !
229 REAL,DIMENSION(ITS:ITE,KTS:KTE-1,JTS:JTE) :: AKH,AKM
230 !
231 REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
232 !
233 !*** Begin debugging
234 REAL :: ZSL_DIAG
235 INTEGER :: IMD,JMD,PRINT_DIAG
236 !*** End debugging
237 !
238 !----------------------------------------------------------------------
239 !**********************************************************************
240 !----------------------------------------------------------------------
241 !
242 !*** Begin debugging
243 IMD=(IMS+IME)/2
244 JMD=(JMS+JME)/2
245 !*** End debugging
246 !
247 !*** MAKE PREPARATIONS
248 !
249 !----------------------------------------------------------------------
250 DTTURBL=DT*STEPBL
251 RDTTURBL=1./DTTURBL
252 DTDIF=DTTURBL
253 RG=1./G
254 !
255 DO J=JTS,JTE
256 DO K=KTS,KTE-1
257 DO I=ITS,ITE
258 AKM(I,K,J)=0.
259 ENDDO
260 ENDDO
261 ENDDO
262 !
263 DO J=JTS,JTE
264 DO K=KTS,KTE+1
265 DO I=ITS,ITE
266 ZINT(I,K,J)=0.
267 ENDDO
268 ENDDO
269 ENDDO
270 !
271 DO J=JTS,JTE
272 DO I=ITS,ITE
273 ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer
274 !
275 !!!!!!!!!
276 !!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
277 !!!!!!!!!
278 !!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer
279 !!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer
280 !
281 ENDDO
282 ENDDO
283 !
284 DO J=JTS,JTE
285 DO K=KTE,KTS,-1
286 KFLIP=KTE+1-K
287 DO I=ITS,ITE
288 ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
289 APEX=1./EXNER(I,K,J)
290 APE(I,K,J)=APEX
291 TX=T(I,K,J)
292 THE(I,K,J)=(CWM(I,K,J)*(-ELOCP/TX)+1.)*TH(I,K,J)
293 ENDDO
294 ENDDO
295 ENDDO
296 !
297 EL_MYJ(its:ite,:,jts:jte) = 0.
298 !
299 !----------------------------------------------------------------------
300 setup_integration: DO J=JTS,JTE
301 !----------------------------------------------------------------------
302 !
303 DO I=ITS,ITE
304 !
305 !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
306 !
307 LMH=KTE-LOWLYR(I,J)+1
308 !
309 PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
310 PSFC=PINT(I,LOWLYR(I,J),J)
311 !
312 !*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
313 !
314 SEAMASK=XLAND(I,J)-1.
315 !
316 !*** FILL 1-D VERTICAL ARRAYS
317 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
318 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
319 !
320 DO K=KTE,KTS,-1
321 KFLIP=KTE+1-K
322 TK(K)=T(I,KFLIP,J)
323 THEK(K)=THE(I,KFLIP,J)
324 RATIOMX=QV(I,KFLIP,J)
325 QK(K)=RATIOMX/(1.+RATIOMX)
326 CWMK(K)=CWM(I,KFLIP,J)
327 PK(K)=PMID(I,KFLIP,J)
328 UK(K)=U(I,KFLIP,J)
329 VK(K)=V(I,KFLIP,J)
330 !
331 !*** TKE=0.5*(q**2) ==> q**2=2.*TKE
332 !
333 Q2K(K)=2.*TKE(I,KFLIP,J)
335 !**** buoyancy flux of mass flux shallow convection scheme
336 !
337 IF ( PRESENT (WTHV_MF) .AND. PRESENT (LM_BL89) ) THEN
338 WTHV_MFK(K)=WTHV_MF(I,KFLIP,J)
339 LM_BL89K(K)=LM_BL89(I,KFLIP,J)
340 ELSE
341 WTHV_MFK(K)=0.
342 LM_BL89K(K)=0.
343 ENDIF
345 !*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
346 !
347 ZHK(K)=ZINT(I,K,J)
348 !
349 ENDDO
350 ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer
351 !
352 !*** Begin debugging
353 ! IF(I==IMD.AND.J==JMD)THEN
354 ! PRINT_DIAG=1
355 ! ELSE
356 ! PRINT_DIAG=0
357 ! ENDIF
358 ! IF(I==227.AND.J==363)PRINT_DIAG=2
359 !*** End debugging
360 !
361 !----------------------------------------------------------------------
362 !***
363 !*** FIND THE MIXING LENGTH
364 !***
365 CALL MIXLEN(LMH,UK,VK,TK,THEK,QK,CWMK &
366 & ,Q2K,ZHK,USTAR(I,J),CORF(I,J),S2,GH,RI,EL &
367 & ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J),LM_BL89K &
368 & ,IDS,IDE,JDS,JDE,KDS,KDE &
369 & ,IMS,IME,JMS,JME,KMS,KME &
370 & ,ITS,ITE,JTS,JTE,KTS,KTE)
371 !
372 !----------------------------------------------------------------------
373 !*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
374 !----------------------------------------------------------------------
375 !
376 KPBL(I,J)=KTE-LPBL(I,J)+1
377 !
378 !----------------------------------------------------------------------
379 !***
380 !*** FIND THE QNSE EXCHANGE COEFFICIENTS
381 !***
382 CALL DIFCOF(LMH,EL,RI,Q2K,ZHK,AKMK,AKHK &
383 & ,IDS,IDE,JDS,JDE,KDS,KDE &
384 & ,IMS,IME,JMS,JME,KMS,KME &
385 & ,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
386 !
387 !----------------------------------------------------------------------
388 !*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
389 !*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. COUNTING
390 !*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
391 !*** ARE DEFINED ON THE TOPS OF THE LAYERS KTS TO KTE-1.
392 !
393 DO K=KTS,KTE-1
394 KFLIP=KTE-K
395 AKH(I,K,J)=AKHK(K)
396 AKM(I,K,J)=AKMK(K)
397 DELTAZ=0.5*(ZHK(KFLIP)-ZHK(KFLIP+2))
398 EXCH_H(I,K+1,J)=AKHK(KFLIP)*DELTAZ
399 EXCH_M(I,K+1,J)=AKMK(KFLIP)*DELTAZ
400 ENDDO
401 !
402 !----------------------------------------------------------------------
403 !***
404 !*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
405 !*** THE TURBULENT KINETIC ENERGY
406 !***
407 !
408 CALL PRODQ2(LMH,DTTURBL,USTAR(I,J),S2,RI,Q2K,EL,ZHK,AKMK,AKHK &
409 & ,WTHV_MFK,IDS,IDE,JDS,JDE,KDS,KDE &
410 & ,IMS,IME,JMS,JME,KMS,KME &
411 & ,ITS,ITE,JTS,JTE,KTS,KTE)
412 !
413 !----------------------------------------------------------------------
414 !***
415 !*** CARRY OUT THE VERTICAL DIFFUSION OF
416 !*** TURBULENT KINETIC ENERGY
417 !***
418 !
419 CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK &
420 & ,IDS,IDE,JDS,JDE,KDS,KDE &
421 & ,IMS,IME,JMS,JME,KMS,KME &
422 & ,ITS,ITE,JTS,JTE,KTS,KTE)
423 !
424 !*** SAVE THE NEW TKE AND MIXING LENGTH.
425 !
426 DO K=KTS,KTE
427 KFLIP=KTE+1-K
428 Q2K(KFLIP)=AMAX1(Q2K(KFLIP),EPSQ2L)
429 TKE(I,K,J)=0.5*Q2K(KFLIP)
430 IF(KFLIP/=KTE)EL_MYJ(I,K,J)=EL(KFLIP) ! EL IS NOT DEFINED AT KTE
431 ENDDO
432 !
433 ENDDO
434 !
435 !----------------------------------------------------------------------
436 ENDDO setup_integration
437 !----------------------------------------------------------------------
438 !
439 !*** CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
440 !
442 DO J=JTS,JTE
443 DO I=ITS,ITE
444 PSFC=PINT(I,LOWLYR(I,J),J)
445 RLOW=PSFC/(R_D*T(I,1,J))
447 THSK(I,J)=THZ0(I,J)
448 THSK(I,J)=HFX(I,J)/(AKHS(I,J)*RLOW*CP)+TH(I,1,J)
450 ENDDO
451 ENDDO
454 !
455 !----------------------------------------------------------------------
456 !
457 !----------------------------------------------------------------------
458 main_integration: DO J=JTS,JTE
459 !----------------------------------------------------------------------
460 !
461 DO I=ITS,ITE
462 !
463 !*** FILL 1-D VERTICAL ARRAYS
464 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
465 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
466 !
467 DO K=KTE,KTS,-1
468 KFLIP=KTE+1-K
469 THEK(K)=THE(I,KFLIP,J)
470 RATIOMX=QV(I,KFLIP,J)
471 QK(K)=RATIOMX/(1.+RATIOMX)
472 CWMK(K)=CWM(I,KFLIP,J)
473 ZHK(K)=ZINT(I,K,J)
474 RHOK(K,I)=PMID(I,KFLIP,J)/(R_D*T(I,KFLIP,J)* &
475 & (1.+P608*QK(K)-CWMK(K)))
476 ENDDO
477 !
478 !*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
479 !*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. THESE COEFFICIENTS
480 !*** ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
481 !
482 DO K=KTS,KTE-1
483 AKHK(K)=AKH(I,K,J)*0.5*(RHOK(K,I)+RHOK(K+1,I))
484 ENDDO
485 !
486 ZHK(KTE+1)=ZINT(I,KTE+1,J)
487 !
488 SEAMASK=XLAND(I,J)-1.
489 THZ0(I,J)=(1.-SEAMASK)*THSK(I,J)+SEAMASK*THZ0(I,J)
490 !
491 LLOW=LOWLYR(I,J)
492 AKHS_DENS=AKHS(I,J)*RHOK(KTE+1-LLOW,I)
493 !
494 IF(SEAMASK<0.5)THEN
495 QFC1=XLV*CHKLOWQ(I,J)*AKHS_DENS
496 !
497 IF(SNOW(I,J)>0..OR.SICE(I,J)>0.5)THEN
498 QFC1=QFC1*RLIVWV
499 ENDIF
500 !
501 IF(QFC1>0.)THEN
502 QLOW=QK(KTE+1-LLOW)
503 QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
504 ENDIF
505 !
506 ELSE
507 PSFC=PINT(I,LOWLYR(I,J),J)
508 EXNSFC=(1.E5/PSFC)**CAPA
509 QSFC(I,J)=PQ0SEA/PSFC &
510 & *EXP(A2*(THSK(I,J)-A3*EXNSFC)/(THSK(I,J)-A4*EXNSFC))
511 ENDIF
512 !
513 QZ0 (I,J)=(1.-SEAMASK)*QSFC(I,J)+SEAMASK*QZ0 (I,J)
514 !
515 !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
516 !
517 LMH=KTE-LOWLYR(I,J)+1
518 !
519 !----------------------------------------------------------------------
520 !*** CARRY OUT THE VERTICAL DIFFUSION OF
521 !*** TEMPERATURE AND WATER VAPOR
522 !----------------------------------------------------------------------
523 !
524 CALL VDIFH(DTDIF,LMH,THZ0(I,J),QZ0(I,J) &
525 & ,AKHS_DENS,CHKLOWQ(I,J),CT(I,J) &
526 & ,THEK,QK,CWMK,AKHK,ZHK,RHOK(KTS,I) &
527 & ,IDS,IDE,JDS,JDE,KDS,KDE &
528 & ,IMS,IME,JMS,JME,KMS,KME &
529 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
530 !----------------------------------------------------------------------
531 !***
532 !*** COMPUTE PRIMARY VARIABLE TENDENCIES
533 !***
534 DO K=KTS,KTE
535 KFLIP=KTE+1-K
536 THOLD=TH(I,K,J)
537 THNEW=THEK(KFLIP)+CWMK(KFLIP)*ELOCP*APE(I,K,J)
538 DTDT=(THNEW-THOLD)*RDTTURBL
539 QOLD=QV(I,K,J)/(1.+QV(I,K,J))
540 DQDT=(QK(KFLIP)-QOLD)*RDTTURBL
541 DCDT=(CWMK(KFLIP)-CWM(I,K,J))*RDTTURBL
542 !
543 RTHBLTEN(I,K,J)=DTDT + RTHBLTEN(I,K,J)
544 RQVBLTEN(I,K,J)=DQDT/(1.-QK(KFLIP))**2 + RQVBLTEN(I,K,J)
545 RQCBLTEN(I,K,J)=DCDT + RQCBLTEN(I,K,J)
546 ENDDO
547 !
548 !*** Begin debugging
549 ! IF(I==IMD.AND.J==JMD)THEN
550 ! PRINT_DIAG=0
551 ! ELSE
552 ! PRINT_DIAG=0
553 ! ENDIF
554 ! IF(I==227.AND.J==363)PRINT_DIAG=0
555 !*** End debugging
556 !
557 PSFC=.01*PINT(I,LOWLYR(I,J),J)
558 ZSL_DIAG=0.5*DZ(I,1,J)
559 !
560 !*** Begin debugging
561 ! IF(PRINT_DIAG==1)THEN
562 !
563 ! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
564 ! '{turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
565 ! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
566 ! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
567 ! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
568 ! '{turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
569 ! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
570 ! , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
571 ! write(6,"(a)") &
572 ! '{turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, S2, EL, Q**2, Akh, EXCH_h, Dz, Dp'
573 ! do k=kts,kte/2
574 ! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
575 ! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
576 ! '{turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
577 ! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), S2(KFLIP) &
578 ! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
579 ! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
580 ! enddo
581 !
582 ! ELSEIF(PRINT_DIAG==2)THEN
583 !
584 ! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
585 ! '}turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
586 ! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
587 ! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
588 ! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
589 ! '}turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
590 ! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
591 ! , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
592 ! write(6,"(a)") &
593 ! '}turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, S2, EL, Q**2, Akh, EXCH_h, Dz, Dp'
594 ! do k=kts,kte/2
595 ! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
596 ! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
597 ! '}turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
598 ! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), S2(KFLIP) &
599 ! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
600 ! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
601 ! enddo
602 ! ENDIF
603 !*** End debugging
604 !
605 !----------------------------------------------------------------------
606 ENDDO
607 !----------------------------------------------------------------------
608 DO I=ITS,ITE
609 !
610 !*** FILL 1-D VERTICAL ARRAYS
611 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
612 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
613 !
614 DO K=KTS,KTE-1
615 AKMK(K)=AKM(I,K,J)
616 AKMK(K)=AKMK(K)*(RHOK(K,I)+RHOK(K+1,I))*0.5
617 ENDDO
618 !
619 LLOW=LOWLYR(I,J)
620 AKMS_DENS=AKMS(I,J)*RHOK(KTE+1-LLOW,I)
621 !
622 DO K=KTE,KTS,-1
623 KFLIP=KTE+1-K
624 UK(K)=U(I,KFLIP,J)
625 VK(K)=V(I,KFLIP,J)
626 ZHK(K)=ZINT(I,K,J)
627 ENDDO
628 ZHK(KTE+1)=ZINT(I,KTE+1,J)
629 !
630 !----------------------------------------------------------------------
631 !*** CARRY OUT THE VERTICAL DIFFUSION OF
632 !*** VELOCITY COMPONENTS
633 !----------------------------------------------------------------------
634 !
635 CALL VDIFV(LMH,DTDIF,UZ0(I,J),VZ0(I,J) &
636 & ,AKMS_DENS,UK,VK,AKMK,ZHK,RHOK(KTS,I) &
637 & ,IDS,IDE,JDS,JDE,KDS,KDE &
638 & ,IMS,IME,JMS,JME,KMS,KME &
639 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
640 !
641 !----------------------------------------------------------------------
642 !***
643 !*** COMPUTE PRIMARY VARIABLE TENDENCIES
644 !***
645 DO K=KTS,KTE
646 KFLIP=KTE+1-K
647 DUDT=(UK(KFLIP)-U(I,K,J))*RDTTURBL
648 DVDT=(VK(KFLIP)-V(I,K,J))*RDTTURBL
649 RUBLTEN(I,K,J)=DUDT + RUBLTEN(I,K,J)
650 RVBLTEN(I,K,J)=DVDT + RVBLTEN(I,K,J)
651 ENDDO
652 !
653 ENDDO
654 !----------------------------------------------------------------------
655 !
656 ENDDO main_integration
657 !
658 !----------------------------------------------------------------------
659 !
660 END SUBROUTINE QNSEPBL
661 !
662 !----------------------------------------------------------------------
663 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
664 !----------------------------------------------------------------------
665 SUBROUTINE MIXLEN &
666 !----------------------------------------------------------------------
667 ! ******************************************************************
668 ! * *
669 ! * LEVEL 2.5 MIXING LENGTH *
670 ! * *
671 ! ******************************************************************
672 !
673 &(LMH,U,V,T,THE,Q,CWM,Q2,Z,USTAR,CORF &
674 &,S2,GH,RI,EL,PBLH,LPBL,LMXL,CT,LM_BL89 &
675 &,IDS,IDE,JDS,JDE,KDS,KDE &
676 &,IMS,IME,JMS,JME,KMS,KME &
677 &,ITS,ITE,JTS,JTE,KTS,KTE)
678 !----------------------------------------------------------------------
679 !
680 IMPLICIT NONE
681 !
682 !----------------------------------------------------------------------
683 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
684 & ,IMS,IME,JMS,JME,KMS,KME &
685 & ,ITS,ITE,JTS,JTE,KTS,KTE
686 !
687 INTEGER,INTENT(IN) :: LMH
688 !
689 INTEGER,INTENT(OUT) :: LMXL,LPBL
690 !
691 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: CWM,Q,Q2,T,THE,U,V,LM_BL89
692 !
693 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
694 !
695 REAL,INTENT(OUT) :: PBLH
696 !
697 REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: EL,RI,GH,S2
698 !
699 REAL,INTENT(INOUT) :: CT
700 !
701 REAL,INTENT(IN) :: CORF,USTAR
702 !----------------------------------------------------------------------
703 !***
704 !*** LOCAL VARIABLES
705 !***
706 INTEGER :: K,LPBLM
707 !
708 REAL :: A,ADEN,B,BDEN,AUBR,BUBR,BLMX,EL0,ELOQ2X,GHL,S2L &
709 & ,QOL2ST,QOL2UN,QDZL,RDZ,SQ,SREL,SZQ,TEM,THM,VKRMZ,RLAMBDA &
710 & ,RLB,RLN,F
711 !
712 REAL,DIMENSION(KTS:KTE) :: Q1,EN2
713 !
714 REAL,DIMENSION(KTS:KTE-1) :: DTH,ELM,REL
715 !
716 !----------------------------------------------------------------------
717 !**********************************************************************
718 !--------------FIND THE HEIGHT OF THE PBL-------------------------------
719 LPBL=LMH
720 !
721 DO K=LMH-1,1,-1
722 IF(Q2(K)<=EPSQ2L*FH)THEN
723 LPBL=K
724 GO TO 110
725 ENDIF
726 ENDDO
727 !
728 LPBL=1
729 !
730 !--------------THE HEIGHT OF THE PBL------------------------------------
731 !
732 110 PBLH=Z(LPBL)-Z(LMH+1)
733 !
734 !-----------------------------------------------------------------------
735 DO K=KTS,LMH
736 Q1(K)=0.
737 ENDDO
738 !
739 DO K=1,LMH-1
740 DTH(K)=THE(K)-THE(K+1)
741 ENDDO
742 !
743 DO K=LMH-2,1,-1
744 IF(DTH(K)>0..AND.DTH(K+1)<=0.)THEN
745 DTH(K)=DTH(K)+CT
746 EXIT
747 ENDIF
748 ENDDO
749 !
750 CT=0.
751 !----------------------------------------------------------------------
752 !*** COMPUTE LOCAL GRADIENT RICHARDSON NUMBER
753 !----------------------------------------------------------------------
754 DO K=KTS,LMH-1
755 RDZ=2./(Z(K)-Z(K+2))
756 S2L=((U(K)-U(K+1))**2+(V(K)-V(K+1))**2)*RDZ*RDZ ! S**2
757 S2L=MAX(S2L,EPSGM)
758 S2(K)=S2L
759 !
760 TEM=(T(K)+T(K+1))*0.5
761 THM=(THE(K)+THE(K+1))*0.5
762 !
763 A=THM*P608
764 B=(ELOCP/TEM-1.-P608)*THM
765 !
766 GHL=(DTH(K)*((Q(K)+Q(K+1)+CWM(K)+CWM(K+1))*(0.5*P608)+1.) &
767 & +(Q(K)-Q(K+1)+CWM(K)-CWM(K+1))*A &
768 & +(CWM(K)-CWM(K+1))*B)*RDZ ! dTheta/dz
769 !
770 IF(ABS(GHL)<=EPSGH)GHL=EPSGH
771 !
772 EN2(K)=GHL*G/THM ! N**2
773 !
774 GH(K)=GHL
775 RI(K)=EN2(K)/S2L
776 ENDDO
777 !
778 !----------------------------------------------------------------------
779 !*** FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
780 !----------------------------------------------------------------------
781 !
782 LMXL=LMH
783 !
784 DO K=KTS,LMH-1
785 S2L=S2(K)
786 GHL=GH(K)
787 !
788 IF(GHL>=EPSGH)THEN
789 IF(S2L/GHL<=REQU)THEN
790 ELM(K)=EPSL
791 LMXL=K
792 ELSE
793 AUBR=(AUBM*S2L+AUBH*GHL)*GHL
794 BUBR= BUBM*S2L+BUBH*GHL
795 QOL2ST=(-0.5*BUBR+SQRT(BUBR*BUBR*0.25-AUBR*CUBR))*RCUBR
796 ELOQ2X=1./QOL2ST
797 ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
798 ENDIF
799 ELSE
800 ADEN=(ADNM*S2L+ADNH*GHL)*GHL
801 BDEN= BDNM*S2L+BDNH*GHL
802 QOL2UN=-0.5*BDEN+SQRT(BDEN*BDEN*0.25-ADEN)
803 ELOQ2X=1./(QOL2UN+EPSRU) ! repsr1/qol2un
804 ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
805 ENDIF
806 ENDDO
807 !
808 IF(ELM(LMH-1)==EPSL)LMXL=LMH
809 !
810 !----------------------------------------------------------------------
811 !*** THE HEIGHT OF THE MIXED LAYER
812 !----------------------------------------------------------------------
813 !
814 BLMX=Z(LMXL)-Z(LMH+1)
815 !
816 !----------------------------------------------------------------------
817 DO K=LPBL,LMH
818 Q1(K)=SQRT(Q2(K))
819 ENDDO
820 !----------------------------------------------------------------------
821 SZQ=0.
822 SQ =0.
823 !
824 DO K=KTS,LMH-1
825 QDZL=(Q1(K)+Q1(K+1))*(Z(K+1)-Z(K+2))
826 SZQ=(Z(K+1)+Z(K+2)-Z(LMH+1)-Z(LMH+1))*QDZL+SZQ
827 SQ=QDZL+SQ
828 ENDDO
829 !
830 !----------------------------------------------------------------------
831 !*** COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
832 !----------------------------------------------------------------------
833 !
834 EL0=MIN(ALPH*SZQ*0.5/SQ,EL0MAX)
835 EL0=MAX(EL0 ,EL0MIN)
836 !
837 !----------------------------------------------------------------------
838 !*** ABOVE THE PBL TOP
839 !----------------------------------------------------------------------
840 !
841 LPBLM=MAX(LPBL-1,1)
842 !
843 DO K=KTS,LPBLM
844 EL(K)=MIN((Z(K)-Z(K+2))*ELFC,ELM(K))
845 REL(K)=EL(K)/ELM(K)
846 ENDDO
847 !
848 !----------------------------------------------------------------------
849 !*** INSIDE THE PBL
850 !----------------------------------------------------------------------
851 !
852 IF(LPBL<LMH)THEN
853 DO K=LPBL,LMH-1
854 VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
855 EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))
856 REL(K)=EL(K)/ELM(K)
857 ENDDO
858 ENDIF
859 !
860 DO K=LPBL+1,LMH-2
861 SREL=MIN(((REL(K-1)+REL(K+1))*0.5+REL(K))*0.5,REL(K))
862 EL(K)=MAX(SREL*ELM(K),EPSL)
863 ENDDO
864 !
865 !----------------------------------------------------------------------
866 !*** MIXING LENGTH FOR THE QNSE MODEL IN STABLE CASE
867 !----------------------------------------------------------------------
868 !
869 F=MAX(CORF,EPS1)
870 RLAMBDA=F/(BLCKDR*USTAR)
871 DO K=KTS,LMH-1
872 IF(EN2(K)>=0.0)THEN ! Stable case
873 VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
874 RLB=RLAMBDA+1./VKRMZ
875 RLN=SQRT(2.*EN2(K)/Q2(K))/CN
876 ! EL(K)=MIN(1./(RLB+RLN),ELM(K))
877 EL(K)=1./(RLB+RLN)
878 ENDIF
879 ENDDO
880 !
881 !----------------------------------------------------------------------
882 END SUBROUTINE MIXLEN
883 !----------------------------------------------------------------------
884 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
885 !----------------------------------------------------------------------
886 SUBROUTINE PRODQ2 &
887 !----------------------------------------------------------------------
888 ! ******************************************************************
889 ! * *
890 ! * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
891 ! * *
892 ! ******************************************************************
893 !
894 &(LMH,DTTURBL,USTAR,S2,RI,Q2,EL,Z,AKM,AKH,WTHV_MF &
895 &,IDS,IDE,JDS,JDE,KDS,KDE &
896 &,IMS,IME,JMS,JME,KMS,KME &
897 &,ITS,ITE,JTS,JTE,KTS,KTE)
898 !----------------------------------------------------------------------
899 !
900 IMPLICIT NONE
901 !
902 !----------------------------------------------------------------------
903 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
904 & ,IMS,IME,JMS,JME,KMS,KME &
905 & ,ITS,ITE,JTS,JTE,KTS,KTE
906 !
907 INTEGER,INTENT(IN) :: LMH
908 !
909 REAL,INTENT(IN) :: DTTURBL,USTAR
910 !
911 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: S2,RI,AKM,AKH,EL
912 !
913 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
914 !
916 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: WTHV_MF
918 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
920 !----------------------------------------------------------------------
921 !***
922 !*** LOCAL VARIABLES
923 !***
924 INTEGER :: K
925 !
926 REAL :: S2L,Q2L,DELTAZ,AKML,AKHL,EN2,PR,BPR,DIS,RC02
927 !
928 !----------------------------------------------------------------------
929 !**********************************************************************
930 !----------------------------------------------------------------------
931 !
932 RC02=2.0/(C0*C0)
933 main_integration: DO K=1,LMH-1
934 S2L=S2(K)
935 Q2L=Q2(K)
936 DELTAZ=0.5*(Z(K)-Z(K+2))
937 AKML=AKM(K)*DELTAZ
938 AKHL=AKH(K)*DELTAZ
939 EN2=RI(K)*S2L !N**2
940 !
941 !*** TURBULENCE PRODUCTION TERM
942 !
943 PR=AKML*S2L
944 !
945 !*** BUOYANCY PRODUCTION
946 !
947 BPR=AKHL*EN2
948 ! BPR=BPR-WTHV_MF(K)
949 !*** DISSIPATION
950 !
951 DIS=CEPS*(0.5*Q2L)**1.5/EL(K)
952 !
953 Q2L=Q2L+2.0*(PR-BPR-DIS)*DTTURBL
954 Q2(K)=AMAX1(Q2L,EPSQ2L)
955 !----------------------------------------------------------------------
956 !*** END OF PRODUCTION/DISSIPATION LOOP
957 !----------------------------------------------------------------------
958 !
959 ENDDO main_integration
960 !
961 !----------------------------------------------------------------------
962 !*** LOWER BOUNDARY CONDITION FOR Q2
963 !----------------------------------------------------------------------
964 !
965 Q2(LMH)=AMAX1(RC02*USTAR*USTAR,EPSQ2L)
966 !----------------------------------------------------------------------
967 !
968 END SUBROUTINE PRODQ2
969 !
970 !----------------------------------------------------------------------
971 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
972 !----------------------------------------------------------------------
973 SUBROUTINE DIFCOF &
974 ! ******************************************************************
975 ! * *
976 ! * DIFFUSION COEFFICIENTS KM, KH BASED ON THE QNSE THEORY *
977 ! * *
978 ! ******************************************************************
979 &(LMH,EL,RI,Q2,Z,AKM,AKH &
980 &,IDS,IDE,JDS,JDE,KDS,KDE &
981 &,IMS,IME,JMS,JME,KMS,KME &
982 &,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
983 !----------------------------------------------------------------------
984 !
985 IMPLICIT NONE
986 !
987 !----------------------------------------------------------------------
988 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
989 & ,IMS,IME,JMS,JME,KMS,KME &
990 & ,ITS,ITE,JTS,JTE,KTS,KTE
991 !
992 INTEGER,INTENT(IN) :: LMH
993 !
994 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: Q2
995 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL,RI
996 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
997 !
998 REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: AKH,AKM
999 !----------------------------------------------------------------------
1000 !***
1001 !*** LOCAL VARIABLES
1002 !***
1003 INTEGER :: K
1004 !
1005 REAL :: AK0,ALPHAM,ALPHAH,RIL,RIL2,ARIL,ARIL2,ARIL4,ELL,Q1L,RDZ &
1006 & ,AK0DZ,AKMIN
1007 !
1008 !*** Begin debugging
1009 INTEGER,INTENT(IN) :: PRINT_DIAG
1010 ! REAL :: D2Tmin
1011 !*** End debugging
1012 !
1013 !----------------------------------------------------------------------
1014 !**********************************************************************
1015 !----------------------------------------------------------------------
1016 !
1017 DO K=1,LMH-1
1018 ELL=EL(K)
1019 Q1L=SQRT(0.5*Q2(K)) !Note that Q1L is SQRT(TKE)
1020 RIL=RI(K)
1021 AK0=C0*ELL*Q1L !KM in neutral case
1022 !
1023 !----------------------------------------------------------------------
1024 !*** STABILITY FUNCTIONS ALPHAM AND ALPHAH
1025 !----------------------------------------------------------------------
1026 !
1027 !!! UNSTABLE CASE
1028 !
1029 IF(RIL<=0) THEN
1030 ARIL=MIN(ABS(RIL),2.*ARIMIN)
1031 ARIL2=ARIL*ARIL
1032 ARIL4=ARIL2*ARIL2
1033 ALPHAM=1.0+BM1*ARIL+BM2*ARIL2
1034 ALPHAH=AH1+BH1*ARIL+BH2*ARIL2+BH3*ARIL4
1035 !
1036 !!! STABLE CASE
1037 !
1038 ELSE
1039 RIL2=RIL*RIL
1040 ALPHAM=(1.0+AM1*RIL2)/(1.0+AM2*RIL+AM3*RIL2)
1041 ALPHAH=(AH1+AH2*RIL+AH3*RIL2)/(1.0+AH4*RIL+AH5*RIL2)
1042 ENDIF
1043 !
1044 !-----------------------------------------------------------------------
1045 !*** END OF STABILITY FUNCTIONS COMPUTATIONS
1046 !-----------------------------------------------------------------------
1047 !
1048 !----------------------------------------------------------------------
1049 !*** DIFFUSION COEFFICIENTS
1050 !----------------------------------------------------------------------
1051 !
1052 RDZ=2./(Z(K)-Z(K+2))
1053 AK0DZ=AK0*RDZ
1054 AKMIN=EPSKM*RDZ
1055 AKM(K)=MAX(AK0DZ*ALPHAM,AKMIN)
1056 AKH(K)=MAX(AK0DZ*ALPHAH,AKMIN)
1057 !----------------------------------------------------------------------
1058 ENDDO
1059 !----------------------------------------------------------------------
1060 !
1061 END SUBROUTINE DIFCOF
1062 !
1063 !----------------------------------------------------------------------
1064 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1065 !----------------------------------------------------------------------
1066 SUBROUTINE VDIFQ &
1067 ! ******************************************************************
1068 ! * *
1069 ! * VERTICAL DIFFUSION OF Q2 (TKE) *
1070 ! * *
1071 ! ******************************************************************
1072 &(LMH,DTDIF,Q2,EL,Z &
1073 &,IDS,IDE,JDS,JDE,KDS,KDE &
1074 &,IMS,IME,JMS,JME,KMS,KME &
1075 &,ITS,ITE,JTS,JTE,KTS,KTE)
1076 !----------------------------------------------------------------------
1077 !
1078 IMPLICIT NONE
1079 !
1080 !----------------------------------------------------------------------
1081 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1082 & ,IMS,IME,JMS,JME,KMS,KME &
1083 & ,ITS,ITE,JTS,JTE,KTS,KTE
1084 !
1085 INTEGER,INTENT(IN) :: LMH
1086 !
1087 REAL,INTENT(IN) :: DTDIF
1088 !
1089 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL
1090 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1091 !
1092 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
1093 !----------------------------------------------------------------------
1094 !***
1095 !*** LOCAL VARIABLES
1096 !***
1097 INTEGER :: K
1098 !
1099 REAL :: ADEN,AKQS,BDEN,BESH,BESM,CDEN,CF,DTOZS,ELL,ELOQ2,ELOQ4 &
1100 & ,ELQDZ,ESH,ESM,ESQHF,GHL,GML,Q1L,RDEN,RDZ
1101 !
1102 REAL,DIMENSION(KTS:KTE-2) :: AKQ,CM,CR,DTOZ,RSQ2
1103 !----------------------------------------------------------------------
1104 !**********************************************************************
1105 !----------------------------------------------------------------------
1106 !***
1107 !*** VERTICAL TURBULENT DIFFUSION
1108 !***
1109 !----------------------------------------------------------------------
1110 ESQHF=0.5*ESQ
1111 !
1112 DO K=KTS,LMH-2
1113 DTOZ(K)=(DTDIF+DTDIF)/(Z(K)-Z(K+2))
1114 AKQ(K)=SQRT((Q2(K)+Q2(K+1))*0.5)*(EL(K)+EL(K+1))*ESQHF &
1115 & /(Z(K+1)-Z(K+2))
1116 CR(K)=-DTOZ(K)*AKQ(K)
1117 ENDDO
1118 !
1119 CM(1)=DTOZ(1)*AKQ(1)+1.
1120 RSQ2(1)=Q2(1)
1121 !
1122 DO K=KTS+1,LMH-2
1123 CF=-DTOZ(K)*AKQ(K-1)/CM(K-1)
1124 CM(K)=-CR(K-1)*CF+(AKQ(K-1)+AKQ(K))*DTOZ(K)+1.
1125 RSQ2(K)=-RSQ2(K-1)*CF+Q2(K)
1126 ENDDO
1127 !
1128 DTOZS=(DTDIF+DTDIF)/(Z(LMH-1)-Z(LMH+1))
1129 AKQS=SQRT((Q2(LMH-1)+Q2(LMH))*0.5)*(EL(LMH-1)+ELZ0)*ESQHF &
1130 & /(Z(LMH)-Z(LMH+1))
1131 !
1132 CF=-DTOZS*AKQ(LMH-2)/CM(LMH-2)
1133 !
1134 Q2(LMH-1)=(DTOZS*AKQS*Q2(LMH)-RSQ2(LMH-2)*CF+Q2(LMH-1)) &
1135 & /((AKQ(LMH-2)+AKQS)*DTOZS-CR(LMH-2)*CF+1.)
1136 !
1137 DO K=LMH-2,KTS,-1
1138 Q2(K)=(-CR(K)*Q2(K+1)+RSQ2(K))/CM(K)
1139 ENDDO
1140 !----------------------------------------------------------------------
1141 !
1142 END SUBROUTINE VDIFQ
1143 !
1144 !----------------------------------------------------------------------
1145 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1146 !---------------------------------------------------------------------
1147 SUBROUTINE VDIFH(DTDIF,LMH,THZ0,QZ0,RKHS,CHKLOWQ,CT &
1148 & ,THE,Q,CWM,RKH,Z,RHO &
1149 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1150 & ,IMS,IME,JMS,JME,KMS,KME &
1151 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1152 ! ***************************************************************
1153 ! * *
1154 ! * VERTICAL DIFFUSION OF MASS VARIABLES *
1155 ! * *
1156 ! ***************************************************************
1157 !---------------------------------------------------------------------
1158 !
1159 IMPLICIT NONE
1160 !
1161 !---------------------------------------------------------------------
1162 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1163 & ,IMS,IME,JMS,JME,KMS,KME &
1164 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1165 !
1166 INTEGER,INTENT(IN) :: LMH
1167 !
1168 REAL,INTENT(IN) :: CHKLOWQ,CT,DTDIF,QZ0,RKHS,THZ0
1169 !
1170 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKH
1171 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1172 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1173 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: CWM,Q,THE
1174 !
1175 !----------------------------------------------------------------------
1176 !***
1177 !*** LOCAL VARIABLES
1178 !***
1179 INTEGER :: K
1180 !
1181 REAL :: CF,CMB,CMCB,CMQB,CMTB,CTHF,DTOZL,DTOZS &
1182 & ,RCML,RKHH,RKQS,RSCB,RSQB,RSTB
1183 !
1184 REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RKCT,RSC,RSQ,RST
1185 !
1186 !----------------------------------------------------------------------
1187 !**********************************************************************
1188 !----------------------------------------------------------------------
1189 CTHF=0.5*CT
1190 !
1191 DO K=KTS,LMH-1
1192 DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1193 CR(K)=-DTOZ(K)*RKH(K)
1194 RKCT(K)=RKH(K)*(Z(K)-Z(K+2))*CTHF
1195 ENDDO
1196 !
1197 CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHO(KTS)
1198 !----------------------------------------------------------------------
1199 RST(KTS)=-RKCT(KTS)*DTOZ(KTS) &
1200 & +THE(KTS)*RHO(KTS)
1201 RSQ(KTS)=Q(KTS) *RHO(KTS)
1202 RSC(KTS)=CWM(KTS)*RHO(KTS)
1203 !----------------------------------------------------------------------
1204 DO K=KTS+1,LMH-1
1205 DTOZL=DTOZ(K)
1206 CF=-DTOZL*RKH(K-1)/CM(K-1)
1207 CM(K)=-CR(K-1)*CF+(RKH(K-1)+RKH(K))*DTOZL+RHO(K)
1208 RST(K)=-RST(K-1)*CF+(RKCT(K-1)-RKCT(K))*DTOZL+THE(K)*RHO(K)
1209 RSQ(K)=-RSQ(K-1)*CF+Q(K) *RHO(K)
1210 RSC(K)=-RSC(K-1)*CF+CWM(K)*RHO(K)
1211 ENDDO
1212 !
1213 DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1214 RKHH=RKH(LMH-1)
1215 !
1216 CF=-DTOZS*RKHH/CM(LMH-1)
1217 RKQS=RKHS*CHKLOWQ
1218 !
1219 CMB=CR(LMH-1)*CF
1220 CMTB=-CMB+(RKHH+RKHS)*DTOZS+RHO(LMH)
1221 CMQB=-CMB+(RKHH+RKQS)*DTOZS+RHO(LMH)
1222 CMCB=-CMB+(RKHH )*DTOZS+RHO(LMH)
1223 !
1224 RSTB=-RST(LMH-1)*CF+RKCT(LMH-1)*DTOZS+THE(LMH)*RHO(LMH)
1225 RSQB=-RSQ(LMH-1)*CF+Q(LMH) *RHO(LMH)
1226 RSCB=-RSC(LMH-1)*CF+CWM(LMH)*RHO(LMH)
1227 !----------------------------------------------------------------------
1228 THE(LMH)=(DTOZS*RKHS*THZ0+RSTB)/CMTB
1229 Q(LMH) =(DTOZS*RKQS*QZ0 +RSQB)/CMQB
1230 CWM(LMH)=( RSCB)/CMCB
1231 !----------------------------------------------------------------------
1232 DO K=LMH-1,KTS,-1
1233 RCML=1./CM(K)
1234 THE(K)=(-CR(K)*THE(K+1)+RST(K))*RCML
1235 Q(K) =(-CR(K)* Q(K+1)+RSQ(K))*RCML
1236 CWM(K)=(-CR(K)*CWM(K+1)+RSC(K))*RCML
1237 ENDDO
1238 !----------------------------------------------------------------------
1239 !
1240 END SUBROUTINE VDIFH
1241 !
1242 !---------------------------------------------------------------------
1243 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1244 !---------------------------------------------------------------------
1245 SUBROUTINE VDIFV(LMH,DTDIF,UZ0,VZ0,RKMS,U,V,RKM,Z,RHO &
1246 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1247 & ,IMS,IME,JMS,JME,KMS,KME &
1248 ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1249 ! ***************************************************************
1250 ! * *
1251 ! * VERTICAL DIFFUSION OF VELOCITY COMPONENTS *
1252 ! * *
1253 ! ***************************************************************
1254 !---------------------------------------------------------------------
1255 !
1256 IMPLICIT NONE
1257 !
1258 !---------------------------------------------------------------------
1259 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1260 & ,IMS,IME,JMS,JME,KMS,KME &
1261 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1262 !
1263 INTEGER,INTENT(IN) :: LMH
1264 !
1265 REAL,INTENT(IN) :: RKMS,DTDIF,UZ0,VZ0
1266 !
1267 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKM
1268 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1269 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1270 !
1271 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: U,V
1272 !----------------------------------------------------------------------
1273 !***
1274 !*** LOCAL VARIABLES
1275 !***
1276 INTEGER :: K
1277 !
1278 REAL :: CF,DTOZAK,DTOZL,DTOZS,RCML,RCMVB,RHOK,RKMH
1279 !
1280 REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RSU,RSV
1281 !----------------------------------------------------------------------
1282 !**********************************************************************
1283 !----------------------------------------------------------------------
1284 DO K=1,LMH-1
1285 DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1286 CR(K)=-DTOZ(K)*RKM(K)
1287 ENDDO
1288 !
1289 RHOK=RHO(1)
1290 CM(1)=DTOZ(1)*RKM(1)+RHOK
1291 RSU(1)=U(1)*RHOK
1292 RSV(1)=V(1)*RHOK
1293 !----------------------------------------------------------------------
1294 DO K=2,LMH-1
1295 DTOZL=DTOZ(K)
1296 CF=-DTOZL*RKM(K-1)/CM(K-1)
1297 RHOK=RHO(K)
1298 CM(K)=-CR(K-1)*CF+(RKM(K-1)+RKM(K))*DTOZL+RHOK
1299 RSU(K)=-RSU(K-1)*CF+U(K)*RHOK
1300 RSV(K)=-RSV(K-1)*CF+V(K)*RHOK
1301 ENDDO
1302 !----------------------------------------------------------------------
1303 DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1304 RKMH=RKM(LMH-1)
1305 !
1306 CF=-DTOZS*RKMH/CM(LMH-1)
1307 RHOK=RHO(LMH)
1308 RCMVB=1./((RKMH+RKMS)*DTOZS-CR(LMH-1)*CF+RHOK)
1309 DTOZAK=DTOZS*RKMS
1310 !----------------------------------------------------------------------
1311 U(LMH)=(DTOZAK*UZ0-RSU(LMH-1)*CF+U(LMH)*RHOK)*RCMVB
1312 V(LMH)=(DTOZAK*VZ0-RSV(LMH-1)*CF+V(LMH)*RHOK)*RCMVB
1313 !----------------------------------------------------------------------
1314 DO K=LMH-1,1,-1
1315 RCML=1./CM(K)
1316 U(K)=(-CR(K)*U(K+1)+RSU(K))*RCML
1317 V(K)=(-CR(K)*V(K+1)+RSV(K))*RCML
1318 ENDDO
1319 !----------------------------------------------------------------------
1320 !
1321 END SUBROUTINE VDIFV
1322 !
1323 !-----------------------------------------------------------------------
1324 !
1325 !=======================================================================
1326 SUBROUTINE QNSEPBLINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1327 & TKE,EXCH_H,RESTART,ALLOWED_TO_READ, &
1328 & IDS,IDE,JDS,JDE,KDS,KDE, &
1329 & IMS,IME,JMS,JME,KMS,KME, &
1330 & ITS,ITE,JTS,JTE,KTS,KTE )
1331 !-----------------------------------------------------------------------
1332 IMPLICIT NONE
1333 !-----------------------------------------------------------------------
1334 LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1335 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1336 & IMS,IME,JMS,JME,KMS,KME, &
1337 & ITS,ITE,JTS,JTE,KTS,KTE
1339 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
1340 & RUBLTEN, &
1341 & RVBLTEN, &
1342 & RTHBLTEN, &
1343 & RQVBLTEN, &
1344 & TKE
1345 INTEGER :: I,J,K,ITF,JTF,KTF
1346 !-----------------------------------------------------------------------
1347 !-----------------------------------------------------------------------
1349 JTF=MIN0(JTE,JDE-1)
1350 KTF=MIN0(KTE,KDE-1)
1351 ITF=MIN0(ITE,IDE-1)
1353 IF(.NOT.RESTART)THEN
1354 DO J=JTS,JTF
1355 DO K=KTS,KTF
1356 DO I=ITS,ITF
1357 TKE(I,K,J)=EPSQ2L
1358 RUBLTEN(I,K,J)=0.
1359 RVBLTEN(I,K,J)=0.
1360 RTHBLTEN(I,K,J)=0.
1361 RQVBLTEN(I,K,J)=0.
1362 EXCH_H(I,K,J)=0.
1363 ENDDO
1364 ENDDO
1365 ENDDO
1366 ENDIF
1368 END SUBROUTINE QNSEPBLINIT
1369 !-----------------------------------------------------------------------
1370 !
1371 END MODULE MODULE_BL_QNSEPBL
1372 !
1373 !-----------------------------------------------------------------------