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