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