| blob | de60d49311bcab80f375450b5b135698648f73a2 |
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,QCW,QCI,QCS,QCR,QCG & !BSF
133 & ,U,V,RHO,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 & !BSF
134 & ,LOWLYR,XLAND,SICE,SNOW &
135 & ,TKE_MYJ,EXCH_H,USTAR,ZNT,EL_MYJ,PBLH,KPBL,CT &
136 & ,AKHS,AKMS,ELFLX,MIXHT & !PLee (3/07)
137 & ,RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN &
138 & ,RQIBLTEN,RQSBLTEN,RQRBLTEN,RQGBLTEN & !BSF
139 & ,IDS,IDE,JDS,JDE,KDS,KDE &
140 & ,IMS,IME,JMS,JME,KMS,KME &
141 & ,ITS,ITE,JTS,JTE,KTS,KTE)
142 !----------------------------------------------------------------------
143 !
144 IMPLICIT NONE
145 !
146 !----------------------------------------------------------------------
147 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
148 & ,IMS,IME,JMS,JME,KMS,KME &
149 & ,ITS,ITE,JTS,JTE,KTS,KTE
150 !
151 INTEGER,INTENT(IN) :: STEPBL
153 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
154 !
155 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: KPBL
156 !
157 REAL,INTENT(IN) :: DT
158 !
159 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,SICE,SNOW &
160 & ,TSK,XLAND
161 !
162 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ,EXNER & !BSF
163 & ,PMID,PINT &
164 & ,RHO & !BSF
165 & ,T,TH,U,V
166 !
167 REAL,OPTIONAL,INTENT(IN),DIMENSION(IMS:IME,KMS:KME,JMS:JME):: & !BSF
168 QV,QCW,QCI,QCS,QCR,QCG !BSF
169 !
170 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PBLH,MIXHT !PLee (3/07)
171 !
172 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS
173 !
174 REAL, INTENT(OUT), DIMENSION(IMS:IME,KMS:KME,JMS:JME) :: & !BSF
175 & EL_MYJ,RTHBLTEN & !BSF
176 & ,RUBLTEN,RVBLTEN
177 !
178 REAL,OPTIONAL,INTENT(OUT),DIMENSION(IMS:IME,KMS:KME,JMS:JME):: & !BSF
179 & RQCBLTEN,RQVBLTEN & !BSF
180 & ,RQIBLTEN,RQSBLTEN & !BSF
181 & ,RQRBLTEN,RQGBLTEN !BSF
182 !
183 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CT,QSFC,QZ0 &
184 & ,THZ0,USTAR &
185 & ,UZ0,VZ0,ZNT
186 !
187 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
188 & ,INTENT(INOUT) :: EXCH_H,TKE_MYJ
189 !
190 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHKLOWQ,ELFLX
191 !
192 !----------------------------------------------------------------------
193 !***
194 !*** LOCAL VARIABLES
195 !***
196 !
197 !-- NSPECE is the number of mass SPECies to be vertically mixed
198 ! (NSPECS - starting number, NSPECE - ending number)
199 ! => Up to 7 species are possible (T,Q,QCW,QCI,QCS,QCR,QCG)
200 !
201 !-- Logic behind the order of QCW,QCI,QCS,QCR,QCG is based on
202 ! increasing terminal fall speed of hydrometeors from cloud
203 ! droplets with very low fall speeds (QCW) to fast falling
204 ! graupel (QCG)
205 !
206 INTEGER, PARAMETER :: NSPECS=1, NSPECE=7 !BSF
207 !
208 !-- flQI,flQS,flQR,flQG controls what hydrometeors are present
209 ! that contribute to total condensate, CWM. These flags may also be
210 ! used in future high resolution runs where all condensate fields
211 ! could be vertically mixed (to be determined). It is assumed
212 ! that cloud water is always being passed into the scheme.
213 !
214 LOGICAL :: flQI,flQS,flQR,flQG !BSF
215 !
216 INTEGER :: I,J,K,KFLIP,LLOW,LMH,LMXL, NUMS,NSPEC,NSP !BSF
217 !
218 INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: LPBL
219 !
220 REAL :: AKHS_DENS,AKMS_DENS,APEX,DELTAZ,DQDT,DTDIF,DTDT & !BSF
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(NSPECS:NSPECE) :: CLOW,CTS,SZ0 !BSF
226 !
227 REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,TK,UK,VK & !BSF
228 & ,QCWK,QCIK,QCSK,QCRK,QCGK !BSF
229 !
230 REAL, DIMENSION(NSPECS:NSPECE,KTS:KTE) :: SPECIES
231 !
232 REAL,DIMENSION(KTS:KTE-1) :: AKHK,AKMK,EL,GH,GM
233 !
234 REAL,DIMENSION(KTS:KTE+1) :: ZHK
235 !
236 REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
237 !
238 REAL,DIMENSION(KTS:KTE,ITS:ITE) :: RHOK
239 !
240 REAL,DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: APE,THE, CWM !BSF
241 !
242 REAL,DIMENSION(ITS:ITE,KTS:KTE-1,JTS:JTE) :: AKH,AKM
243 !
244 REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
245 !
246 !*** Begin debugging
247 REAL :: ZSL_DIAG
248 INTEGER :: IMD,JMD,PRINT_DIAG
249 !*** End debugging
250 !
251 !----------------------------------------------------------------------
252 !**********************************************************************
253 !----------------------------------------------------------------------
254 !
255 !*** Begin debugging
256 IMD=(IMS+IME)/2
257 JMD=(JMS+JME)/2
258 !*** End debugging
259 !
260 !*** MAKE PREPARATIONS
261 !
262 !----------------------------------------------------------------------
263 !BSF Begin:
264 NSPEC=3 !-- T, Q, Cloud water
265 flQI=.FALSE.
266 flQS=.FALSE.
267 flQR=.FALSE.
268 flQG=.FALSE.
269 !
270 IF (PRESENT(QCI)) THEN
271 flQI=.TRUE.
272 NSPEC=NSPEC+1
273 ENDIF
274 !BSF IF (PRESENT(QCS)) THEN
275 !BSF flQS=.TRUE.
276 !BSF NSPEC=NSPEC+1
277 !BSF ENDIF
278 !BSF IF (PRESENT(QCR)) THEN
279 !BSF flQR=.TRUE.
280 !BSF NSPEC=NSPEC+1
281 !BSF ENDIF
282 !BSF IF (PRESENT(QCG)) THEN
283 !BSF flQG=.TRUE.
284 !BSF NSPEC=NSPEC+1
285 !BSF ENDIF
286 !
287 !BSF End:
289 DTTURBL=DT*STEPBL
290 RDTTURBL=1./DTTURBL
291 DTDIF=DTTURBL
292 RG=1./G
293 !
294 DO J=JTS,JTE
295 DO K=KTS,KTE-1
296 DO I=ITS,ITE
297 AKM(I,K,J)=0.
298 ENDDO
299 ENDDO
300 ENDDO
301 !
302 DO J=JTS,JTE
303 DO K=KTS,KTE+1
304 DO I=ITS,ITE
305 ZINT(I,K,J)=0.
306 ENDDO
307 ENDDO
308 ENDDO
309 !
310 DO J=JTS,JTE
311 DO I=ITS,ITE
312 ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer
313 !
314 !!!!!!!!!
315 !!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
316 !!!!!!!!!
317 !!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer
318 !!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer
319 !
320 ENDDO
321 ENDDO
322 !
323 DO J=JTS,JTE
324 DO K=KTE,KTS,-1
325 KFLIP=KTE+1-K
326 DO I=ITS,ITE
327 ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
328 APEX=1./EXNER(I,K,J)
329 APE(I,K,J)=APEX
330 TX=T(I,K,J)
331 CWM(I,K,J)=QCW(I,K,J) !BSF
332 IF (flQI) CWM(I,K,J)=CWM(I,K,J)+QCI(I,K,J) !BSF
333 !BSF IF (flQS) CWM(I,K,J)=CWM(I,K,J)+QCS(I,K,J)
334 !BSF IF (flQR) CWM(I,K,J)=CWM(I,K,J)+QCR(I,K,J)
335 !BSF IF (flQG) CWM(I,K,J)=CWM(I,K,J)+QCG(I,K,J)
336 THE(I,K,J)=(CWM(I,K,J)*(-ELOCP/TX)+1.)*TH(I,K,J)
337 ENDDO
338 ENDDO
339 ENDDO
340 !
341 EL_MYJ(its:ite,:,jts:jte) = 0.
342 !
343 !----------------------------------------------------------------------
344 setup_integration: DO J=JTS,JTE
345 !----------------------------------------------------------------------
346 !
347 DO I=ITS,ITE
348 !
349 !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
350 !
351 LMH=KTE-LOWLYR(I,J)+1
352 !
353 PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
354 PSFC=PINT(I,LOWLYR(I,J),J)
355 !
356 !*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
357 !
358 SEAMASK=XLAND(I,J)-1.
359 !
360 !*** FILL 1-D VERTICAL ARRAYS
361 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
362 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
363 !
364 DO K=KTE,KTS,-1
365 KFLIP=KTE+1-K
366 TK(K)=T(I,KFLIP,J)
367 THEK(K)=THE(I,KFLIP,J)
368 RATIOMX=QV(I,KFLIP,J)
369 QK(K)=RATIOMX/(1.+RATIOMX)
370 CWMK(K)=CWM(I,KFLIP,J)
371 PK(K)=PMID(I,KFLIP,J)
372 UK(K)=U(I,KFLIP,J)
373 VK(K)=V(I,KFLIP,J)
374 !
375 !*** TKE=0.5*(q**2) ==> q**2=2.*TKE
376 !
377 Q2K(K)=2.*TKE_MYJ(I,KFLIP,J)
378 !
379 !*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
380 !
381 ZHK(K)=ZINT(I,K,J)
382 !
383 ENDDO
384 ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer
385 !
386 !*** Begin debugging
387 ! IF(I==IMD.AND.J==JMD)THEN
388 ! PRINT_DIAG=1
389 ! ELSE
390 ! PRINT_DIAG=0
391 ! ENDIF
392 ! IF(I==227.AND.J==363)PRINT_DIAG=2
393 !*** End debugging
394 !
395 !----------------------------------------------------------------------
396 !***
397 !*** FIND THE MIXING LENGTH
398 !***
399 CALL MIXLEN(LMH,UK,VK,TK,THEK,QK,CWMK &
400 & ,Q2K,ZHK,GM,GH,EL &
401 & ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J),MIXHT(I,J) & !PLee (3/07)
402 & ,IDS,IDE,JDS,JDE,KDS,KDE &
403 & ,IMS,IME,JMS,JME,KMS,KME &
404 & ,ITS,ITE,JTS,JTE,KTS,KTE)
405 !
406 !----------------------------------------------------------------------
407 !***
408 !*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
409 !*** THE TURBULENT KINETIC ENERGY
410 !***
411 !
412 CALL PRODQ2(LMH,DTTURBL,USTAR(I,J),GM,GH,EL,Q2K &
413 & ,IDS,IDE,JDS,JDE,KDS,KDE &
414 & ,IMS,IME,JMS,JME,KMS,KME &
415 & ,ITS,ITE,JTS,JTE,KTS,KTE)
416 !
417 !----------------------------------------------------------------------
418 !*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
419 !----------------------------------------------------------------------
420 !
421 KPBL(I,J)=KTE-LPBL(I,J)+1
422 !
423 !----------------------------------------------------------------------
424 !***
425 !*** FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
426 !***
427 CALL DIFCOF(LMH,LMXL,GM,GH,EL,TK,Q2K,ZHK,AKMK,AKHK &
428 & ,IDS,IDE,JDS,JDE,KDS,KDE &
429 & ,IMS,IME,JMS,JME,KMS,KME &
430 & ,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
431 !
432 !*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
433 !*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1.
434 !*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS AS
435 !*** EXCH_H ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1
436 !*** THUS EXCH_H INDICES INCREASE UPWARD WITH K=KTS AT THE GROUND.
437 !
438 DO K=KTS,KTE-1
439 KFLIP=KTE-K
440 AKH(I,K,J)=AKHK(K)
441 AKM(I,K,J)=AKMK(K)
442 DELTAZ=0.5*(ZHK(KFLIP)-ZHK(KFLIP+2))
443 EXCH_H(I,K,J)=AKHK(KFLIP)*DELTAZ
444 ENDDO
445 !
446 !----------------------------------------------------------------------
447 !***
448 !*** CARRY OUT THE VERTICAL DIFFUSION OF
449 !*** TURBULENT KINETIC ENERGY
450 !***
451 !
452 CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK &
453 & ,IDS,IDE,JDS,JDE,KDS,KDE &
454 & ,IMS,IME,JMS,JME,KMS,KME &
455 & ,ITS,ITE,JTS,JTE,KTS,KTE)
456 !
457 !*** SAVE THE NEW TKE AND MIXING LENGTH.
458 !
459 DO K=KTS,KTE
460 KFLIP=KTE+1-K
461 Q2K(KFLIP)=AMAX1(Q2K(KFLIP),EPSQ2)
462 TKE_MYJ(I,K,J)=0.5*Q2K(KFLIP)
463 IF(KFLIP<KTE)EL_MYJ(I,K,J)=EL(KFLIP) ! EL IS NOT DEFINED AT KTE (ground surface)
464 ENDDO
465 !
466 ENDDO
467 !
468 !----------------------------------------------------------------------
469 ENDDO setup_integration
470 !----------------------------------------------------------------------
471 !
472 !*** CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
473 !
474 DO J=JTS,JTE
475 DO I=ITS,ITE
476 PSFC=PINT(I,LOWLYR(I,J),J)
477 THSK(I,J)=TSK(I,J)*(1.E5/PSFC)**CAPA
478 ENDDO
479 ENDDO
480 !
481 !----------------------------------------------------------------------
482 !
483 !----------------------------------------------------------------------
484 main_integration: DO J=JTS,JTE
485 !----------------------------------------------------------------------
486 !
487 DO I=ITS,ITE
488 !
489 !*** FILL 1-D VERTICAL ARRAYS
490 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
491 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
492 !
493 DO K=KTE,KTS,-1
494 KFLIP=KTE+1-K
495 THEK(K)=THE(I,KFLIP,J)
496 RATIOMX=QV(I,KFLIP,J)
497 QK(K)=RATIOMX/(1.+RATIOMX)
498 CWMK(K)=CWM(I,KFLIP,J)
499 QCWK(K)=QCW(I,KFLIP,J)
500 !BSF Begin:
501 SPECIES(1,K)=THEK(K)
502 SPECIES(2,K)=QK(K)
503 SPECIES(3,K)=QCWK(K)
504 NSP=3
505 IF (flQI) THEN
506 NSP=NSP+1
507 QCIK(K)=QCI(I,KFLIP,J)
508 SPECIES(NSP,K)=QCIK(K)
509 ENDIF
510 !BSF IF (flQS) THEN
511 !BSF NSP=NSP+1
512 !BSF QCSK(K)=QCS(I,KFLIP,J)
513 !BSF SPECIES(NSP,K)=QCSK(K)
514 !BSF ENDIF
515 !BSF IF (flQR) THEN
516 !BSF NSP=NSP+1
517 !BSF QCRK(K)=QCR(I,KFLIP,J)
518 !BSF SPECIES(NSP,K)=QCRK(K)
519 !BSF ENDIF
520 !BSF IF (flQG) THEN
521 !BSF NSP=NSP+1
522 !BSF QCGK(K)=QCG(I,KFLIP,J)
523 !BSF SPECIES(NSP,K)=QCGK(K)
524 !BSF ENDIF
525 !BSF End:
526 ZHK(K)=ZINT(I,K,J)
527 RHOK(K,I)=PMID(I,KFLIP,J)/(R_D*T(I,KFLIP,J)* &
528 & (1.+P608*QK(K)-CWMK(K)))
529 ENDDO
530 !
531 !*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
532 !*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1. THESE COEFFICIENTS
533 !*** ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
534 !
535 DO K=KTS,KTE-1
536 AKHK(K)=AKH(I,K,J)*0.5*(RHOK(K,I)+RHOK(K+1,I))
537 ENDDO
538 !
539 ZHK(KTE+1)=ZINT(I,KTE+1,J)
540 !
541 SEAMASK=XLAND(I,J)-1.
542 THZ0(I,J)=(1.-SEAMASK)*THSK(I,J)+SEAMASK*THZ0(I,J)
543 !
544 LLOW=LOWLYR(I,J)
545 AKHS_DENS=AKHS(I,J)*RHOK(KTE+1-LLOW,I)
546 !
547 IF(SEAMASK<0.5)THEN
548 QFC1=XLV*CHKLOWQ(I,J)*AKHS_DENS
549 !
550 IF(SNOW(I,J)>0..OR.SICE(I,J)>0.5)THEN
551 QFC1=QFC1*RLIVWV
552 ENDIF
553 !
554 IF(QFC1>0.)THEN
555 QLOW=QK(KTE+1-LLOW)
556 QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
557 ENDIF
558 !
559 ELSE
560 PSFC=PINT(I,LOWLYR(I,J),J)
561 EXNSFC=(1.E5/PSFC)**CAPA
562 QSFC(I,J)=PQ0SEA/PSFC &
563 & *EXP(A2*(THSK(I,J)-A3*EXNSFC)/(THSK(I,J)-A4*EXNSFC))
564 ENDIF
565 !
566 QZ0 (I,J)=(1.-SEAMASK)*QSFC(I,J)+SEAMASK*QZ0 (I,J)
567 !
568 ! BSF: Begin
569 SZ0(1)=THZ0(I,J)
570 SZ0(2)=QZ0 (I,J)
571 DO NUMS=3,NSPECE
572 SZ0(NUMS)=0.
573 ENDDO
574 !
575 CLOW(1)=1.
576 CLOW(2)=CHKLOWQ(I,J)
577 DO NUMS=3,NSPECE
578 CLOW(NUMS)=0.
579 ENDDO
580 !
581 CTS(1)=CT(I,J)
582 DO NUMS=2,NSPECE
583 CTS(NUMS)=0.
584 ENDDO
585 ! BSF: End
586 !
587 !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
588 !
589 LMH=KTE-LOWLYR(I,J)+1
590 !
591 !----------------------------------------------------------------------
592 !*** CARRY OUT THE VERTICAL DIFFUSION OF
593 !*** TEMPERATURE AND WATER VAPOR
594 !----------------------------------------------------------------------
595 !
596 CALL VDIFH(DTDIF,LMH,NSPECS,NSPECE,LPBL(I,J),SZ0 & !BSF
597 & ,AKHS_DENS,CLOW,CTS & !BSF
598 & ,SPECIES,NSPEC,AKHK,ZHK,RHOK(KTS,I) & !BSF
599 & ,IDS,IDE,JDS,JDE,KDS,KDE &
600 & ,IMS,IME,JMS,JME,KMS,KME &
601 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
602 !----------------------------------------------------------------------
603 !***
604 !*** COMPUTE PRIMARY VARIABLE TENDENCIES
605 !***
606 !BSF Begin:
607 DO K=KTS,KTE
608 THEK(K)=SPECIES(1,K)
609 QK (K)=SPECIES(2,K)
610 QCWK(K)=SPECIES(3,K)
611 CWMK(K)=QCWK(K)
612 NSP=3
613 IF (flQI) THEN
614 NSP=NSP+1
615 QCIK(K)=SPECIES(NSP,K)
616 CWMK(K)=CWMK(K)+QCIK(K)
617 ENDIF
618 !BSF IF (flQS) THEN
619 !BSF NSP=NSP+1
620 !BSF QCSK(K)=SPECIES(NSP,K)
621 !BSF CWMK(K)=CWMK(K)+QCSK(K)
622 !BSF ENDIF
623 !BSF IF (flQR) THEN
624 !BSF NSP=NSP+1
625 !BSF QCRK(K)=SPECIES(NSP,K)
626 !BSF CWMK(K)=CWMK(K)+QCRK(K)
627 !BSF ENDIF
628 !BSF IF (flQG) THEN
629 !BSF NSP=NSP+1
630 !BSF QCGK(K)=SPECIES(NSP,K)
631 !BSF CWMK(K)=CWMK(K)+QCGK(K)
632 !BSF ENDIF
633 ENDDO
634 !BSF End:
635 !
636 DO K=KTS,KTE
637 KFLIP=KTE+1-K
638 THOLD=TH(I,K,J)
639 THNEW=THEK(KFLIP)+CWMK(KFLIP)*ELOCP*APE(I,K,J)
640 DTDT=(THNEW-THOLD)*RDTTURBL
641 QOLD=QV(I,K,J)/(1.+QV(I,K,J))
642 DQDT=(QK(KFLIP)-QOLD)*RDTTURBL
643 RTHBLTEN(I,K,J)=DTDT
644 RQVBLTEN(I,K,J)=DQDT/(1.-QK(KFLIP))**2
645 !BSF Begin:
646 RQCBLTEN(I,K,J)=(QCWK(KFLIP)-QCW(I,K,J))*RDTTURBL
647 IF (flQI) RQIBLTEN(I,K,J)=(QCIK(KFLIP)-QCI(I,K,J))*RDTTURBL
648 !BSF IF (flQS) RQSBLTEN(I,K,J)=(QCSK(KFLIP)-QCS(I,K,J))*RDTTURBL
649 !BSF IF (flQR) RQRBLTEN(I,K,J)=(QCRK(KFLIP)-QCR(I,K,J))*RDTTURBL
650 !BSF IF (flQG) RQGBLTEN(I,K,J)=(QCGK(KFLIP)-QCG(I,K,J))*RDTTURBL
651 !BSF End:
652 ENDDO
653 !
654 !*** Begin debugging
655 ! IF(I==IMD.AND.J==JMD)THEN
656 ! PRINT_DIAG=0
657 ! ELSE
658 ! PRINT_DIAG=0
659 ! ENDIF
660 ! IF(I==227.AND.J==363)PRINT_DIAG=0
661 !*** End debugging
662 !
663 PSFC=.01*PINT(I,LOWLYR(I,J),J)
664 ZSL_DIAG=0.5*DZ(I,1,J)
665 !
666 !*** Begin debugging
667 ! IF(PRINT_DIAG==1)THEN
668 !
669 ! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
670 ! '{turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
671 ! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
672 ! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
673 ! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
674 ! '{turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
675 ! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
676 ! , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
677 ! write(6,"(a)") &
678 ! '{turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
679 ! do k=kts,kte/2
680 ! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
681 ! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
682 ! '{turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
683 ! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
684 ! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
685 ! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
686 ! enddo
687 !
688 ! ELSEIF(PRINT_DIAG==2)THEN
689 !
690 ! write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
691 ! '}turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
692 ! , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag &
693 ! , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
694 ! write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
695 ! '}turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
696 ! , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
697 ! , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
698 ! write(6,"(a)") &
699 ! '}turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
700 ! do k=kts,kte/2
701 ! KFLIP=KTE-K !-- Includes the KFLIP-1 in earlier versions
702 ! write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
703 ! '}turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
704 ! , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
705 ! , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
706 ! , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
707 ! enddo
708 ! ENDIF
709 !*** End debugging
710 !
711 !----------------------------------------------------------------------
712 ENDDO
713 !----------------------------------------------------------------------
714 DO I=ITS,ITE
715 !
716 !*** FILL 1-D VERTICAL ARRAYS
717 !*** AND FLIP DIRECTION SINCE MYJ SCHEME
718 !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
719 !
720 DO K=KTS,KTE-1
721 AKMK(K)=AKM(I,K,J)
722 AKMK(K)=AKMK(K)*(RHOK(K,I)+RHOK(K+1,I))*0.5
723 ENDDO
724 !
725 LLOW=LOWLYR(I,J)
726 AKMS_DENS=AKMS(I,J)*RHOK(KTE+1-LLOW,I)
727 !
728 DO K=KTE,KTS,-1
729 KFLIP=KTE+1-K
730 UK(K)=U(I,KFLIP,J)
731 VK(K)=V(I,KFLIP,J)
732 ZHK(K)=ZINT(I,K,J)
733 ENDDO
734 ZHK(KTE+1)=ZINT(I,KTE+1,J)
735 !
736 !----------------------------------------------------------------------
737 !*** CARRY OUT THE VERTICAL DIFFUSION OF
738 !*** VELOCITY COMPONENTS
739 !----------------------------------------------------------------------
740 !
741 CALL VDIFV(LMH,DTDIF,UZ0(I,J),VZ0(I,J) &
742 & ,AKMS_DENS,UK,VK,AKMK,ZHK,RHOK(KTS,I) &
743 & ,IDS,IDE,JDS,JDE,KDS,KDE &
744 & ,IMS,IME,JMS,JME,KMS,KME &
745 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
746 !
747 !----------------------------------------------------------------------
748 !***
749 !*** COMPUTE PRIMARY VARIABLE TENDENCIES
750 !***
751 DO K=KTS,KTE
752 KFLIP=KTE+1-K
753 DUDT=(UK(KFLIP)-U(I,K,J))*RDTTURBL
754 DVDT=(VK(KFLIP)-V(I,K,J))*RDTTURBL
755 RUBLTEN(I,K,J)=DUDT
756 RVBLTEN(I,K,J)=DVDT
757 ENDDO
758 !
759 ENDDO
760 !----------------------------------------------------------------------
761 !
762 ENDDO main_integration
763 !
764 !----------------------------------------------------------------------
765 !
766 END SUBROUTINE MYJPBL
767 !
768 !----------------------------------------------------------------------
769 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
770 !----------------------------------------------------------------------
771 SUBROUTINE MIXLEN &
772 !----------------------------------------------------------------------
773 ! ******************************************************************
774 ! * *
775 ! * LEVEL 2.5 MIXING LENGTH *
776 ! * *
777 ! ******************************************************************
778 !
779 &(LMH,U,V,T,THE,Q,CWM,Q2,Z,GM,GH,EL,PBLH,LPBL,LMXL,CT,MIXHT & !PLee(3/07)
780 &,IDS,IDE,JDS,JDE,KDS,KDE &
781 &,IMS,IME,JMS,JME,KMS,KME &
782 &,ITS,ITE,JTS,JTE,KTS,KTE)
783 !----------------------------------------------------------------------
784 !
785 IMPLICIT NONE
786 !
787 !----------------------------------------------------------------------
788 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
789 & ,IMS,IME,JMS,JME,KMS,KME &
790 & ,ITS,ITE,JTS,JTE,KTS,KTE
791 !
792 INTEGER,INTENT(IN) :: LMH
793 !
794 INTEGER,INTENT(OUT) :: LMXL,LPBL
795 !
796 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: CWM,Q,Q2,T,THE,U,V
797 !
798 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
799 !
800 REAL,INTENT(OUT) :: PBLH,MIXHT !PLee (3/07)
801 !
802 REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: EL,GH,GM
803 !
804 REAL,INTENT(INOUT) :: CT
805 !----------------------------------------------------------------------
806 !***
807 !*** LOCAL VARIABLES
808 !***
809 INTEGER :: K,LPBLM
810 !
811 REAL :: A,ADEN,B,BDEN,AUBR,BUBR,BLMX,EL0,ELOQ2X,GHL,GML &
812 & ,QOL2ST,QOL2UN,QDZL,RDZ,SQ,SREL,SZQ,TEM,THM,VKRMZ
813 !
814 REAL,DIMENSION(KTS:KTE) :: Q1
815 !
816 REAL,DIMENSION(KTS:KTE-1) :: DTH,ELM,REL
817 !
818 !----------------------------------------------------------------------
819 !**********************************************************************
820 !--------------FIND THE HEIGHT OF THE PBL-------------------------------
821 LPBL=LMH
822 !
823 DO K=LMH-1,1,-1
824 IF(Q2(K)<=EPSQ2*FH)THEN
825 LPBL=K
826 GO TO 110
827 ENDIF
828 ENDDO
829 !
830 LPBL=1
831 !
832 !--------------THE HEIGHT OF THE PBL------------------------------------
833 !
834 110 PBLH=Z(LPBL+1)-Z(LMH+1)
835 !
836 !-----------------------------------------------------------------------
837 DO K=KTS,LMH
838 Q1(K)=0.
839 ENDDO
840 !
841 DO K=1,LMH-1
842 DTH(K)=THE(K)-THE(K+1)
843 ENDDO
844 !
845 DO K=LMH-2,1,-1
846 IF(DTH(K)>0..AND.DTH(K+1)<=0.)THEN
847 DTH(K)=DTH(K)+CT
848 EXIT
849 ENDIF
850 ENDDO
851 !
852 CT=0.
853 !----------------------------------------------------------------------
854 DO K=KTS,LMH-1
855 RDZ=2./(Z(K)-Z(K+2))
856 GML=((U(K)-U(K+1))**2+(V(K)-V(K+1))**2)*RDZ*RDZ
857 GM(K)=MAX(GML,EPSGM)
858 !
859 TEM=(T(K)+T(K+1))*0.5
860 THM=(THE(K)+THE(K+1))*0.5
861 !
862 A=THM*P608
863 B=(ELOCP/TEM-1.-P608)*THM
864 !
865 GHL=(DTH(K)*((Q(K)+Q(K+1)+CWM(K)+CWM(K+1))*(0.5*P608)+1.) &
866 & +(Q(K)-Q(K+1)+CWM(K)-CWM(K+1))*A &
867 & +(CWM(K)-CWM(K+1))*B)*RDZ
868 !
869 IF(ABS(GHL)<=EPSGH)GHL=EPSGH
870 !
871 GH(K)=GHL
872 ENDDO
873 !
874 !----------------------------------------------------------------------
875 !*** FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
876 !----------------------------------------------------------------------
877 !
878 LMXL=LMH
879 !
880 DO K=KTS,LMH-1
881 GML=GM(K)
882 GHL=GH(K)
883 !
884 IF(GHL>=EPSGH)THEN
885 IF(GML/GHL<=REQU)THEN
886 ELM(K)=EPSL
887 LMXL=K
888 ELSE
889 AUBR=(AUBM*GML+AUBH*GHL)*GHL
890 BUBR= BUBM*GML+BUBH*GHL
891 QOL2ST=(-0.5*BUBR+SQRT(BUBR*BUBR*0.25-AUBR*CUBR))*RCUBR
892 ELOQ2X=1./QOL2ST
893 ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
894 ENDIF
895 ELSE
896 ADEN=(ADNM*GML+ADNH*GHL)*GHL
897 BDEN= BDNM*GML+BDNH*GHL
898 QOL2UN=-0.5*BDEN+SQRT(BDEN*BDEN*0.25-ADEN)
899 ELOQ2X=1./(QOL2UN+EPSRU) ! repsr1/qol2un
900 ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
901 ENDIF
902 ENDDO
903 !
904 IF(ELM(LMH-1)==EPSL)LMXL=LMH
905 !
906 !----------------------------------------------------------------------
907 !*** THE HEIGHT OF THE MIXED LAYER
908 !----------------------------------------------------------------------
909 !
910 BLMX=Z(LMXL)-Z(LMH+1)
911 MIXHT=BLMX !PLee (3/07)
912 !
913 !----------------------------------------------------------------------
914 DO K=LPBL,LMH
915 Q1(K)=SQRT(Q2(K))
916 ENDDO
917 !----------------------------------------------------------------------
918 SZQ=0.
919 SQ =0.
920 !
921 DO K=KTS,LMH-1
922 QDZL=(Q1(K)+Q1(K+1))*(Z(K+1)-Z(K+2))
923 SZQ=(Z(K+1)+Z(K+2)-Z(LMH+1)-Z(LMH+1))*QDZL+SZQ
924 SQ=QDZL+SQ
925 ENDDO
926 !
927 !----------------------------------------------------------------------
928 !*** COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
929 !----------------------------------------------------------------------
930 !
931 EL0=MIN(ALPH*SZQ*0.5/SQ,EL0MAX)
932 EL0=MAX(EL0 ,EL0MIN)
933 !
934 !----------------------------------------------------------------------
935 !*** ABOVE THE PBL TOP
936 !----------------------------------------------------------------------
937 !
938 LPBLM=MAX(LPBL-1,1)
939 !
940 DO K=KTS,LPBLM
941 EL(K)=MIN((Z(K)-Z(K+2))*ELFC,ELM(K))
942 REL(K)=EL(K)/ELM(K)
943 ENDDO
944 !
945 !----------------------------------------------------------------------
946 !*** INSIDE THE PBL
947 !----------------------------------------------------------------------
948 !
949 IF(LPBL<LMH)THEN
950 DO K=LPBL,LMH-1
951 VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
952 EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))
953 REL(K)=EL(K)/ELM(K)
954 ENDDO
955 ENDIF
956 !
957 DO K=LPBL+1,LMH-2
958 SREL=MIN(((REL(K-1)+REL(K+1))*0.5+REL(K))*0.5,REL(K))
959 EL(K)=MAX(SREL*ELM(K),EPSL)
960 ENDDO
961 !
962 !----------------------------------------------------------------------
963 END SUBROUTINE MIXLEN
964 !----------------------------------------------------------------------
965 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
966 !----------------------------------------------------------------------
967 SUBROUTINE PRODQ2 &
968 !----------------------------------------------------------------------
969 ! ******************************************************************
970 ! * *
971 ! * LEVEL 2.5 Q2 PRODUCTION/DISSIPATION *
972 ! * *
973 ! ******************************************************************
974 !
975 &(LMH,DTTURBL,USTAR,GM,GH,EL,Q2 &
976 &,IDS,IDE,JDS,JDE,KDS,KDE &
977 &,IMS,IME,JMS,JME,KMS,KME &
978 &,ITS,ITE,JTS,JTE,KTS,KTE)
979 !----------------------------------------------------------------------
980 !
981 IMPLICIT NONE
982 !
983 !----------------------------------------------------------------------
984 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
985 & ,IMS,IME,JMS,JME,KMS,KME &
986 & ,ITS,ITE,JTS,JTE,KTS,KTE
987 !
988 INTEGER,INTENT(IN) :: LMH
989 !
990 REAL,INTENT(IN) :: DTTURBL,USTAR
991 !
992 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: GH,GM
993 REAL,DIMENSION(KTS:KTE-1),INTENT(INOUT) :: EL
994 !
995 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
996 !----------------------------------------------------------------------
997 !***
998 !*** LOCAL VARIABLES
999 !***
1000 INTEGER :: K
1001 !
1002 REAL :: ADEN,AEQU,ANUM,ARHS,BDEN,BEQU,BNUM,BRHS,CDEN,CRHS &
1003 & ,DLOQ1,ELOQ11,ELOQ12,ELOQ13,ELOQ21,ELOQ22,ELOQ31,ELOQ32 &
1004 & ,ELOQ41,ELOQ42,ELOQ51,ELOQ52,ELOQN,EQOL2,GHL,GML &
1005 & ,RDEN1,RDEN2,RHS2,RHSP1,RHSP2,RHST2
1006 !
1007 !----------------------------------------------------------------------
1008 !**********************************************************************
1009 !----------------------------------------------------------------------
1010 !
1011 main_integration: DO K=1,LMH-1
1012 GML=GM(K)
1013 GHL=GH(K)
1014 !
1015 !----------------------------------------------------------------------
1016 !*** COEFFICIENTS OF THE EQUILIBRIUM EQUATION
1017 !----------------------------------------------------------------------
1018 !
1019 AEQU=(AEQM*GML+AEQH*GHL)*GHL
1020 BEQU= BEQM*GML+BEQH*GHL
1021 !
1022 !----------------------------------------------------------------------
1023 !*** EQUILIBRIUM SOLUTION FOR L/Q
1024 !----------------------------------------------------------------------
1025 !
1026 EQOL2=-0.5*BEQU+SQRT(BEQU*BEQU*0.25-AEQU)
1027 !
1028 !----------------------------------------------------------------------
1029 !*** IS THERE PRODUCTION/DISSIPATION ?
1030 !----------------------------------------------------------------------
1031 !
1032 IF((GML+GHL*GHL<=EPSTRB) &
1033 & .OR.(GHL>=EPSGH.AND.GML/GHL<=REQU) &
1034 & .OR.(EQOL2<=EPS2))THEN
1035 !
1036 !----------------------------------------------------------------------
1037 !*** NO TURBULENCE
1038 !----------------------------------------------------------------------
1039 !
1040 Q2(K)=EPSQ2
1041 EL(K)=EPSL
1042 !----------------------------------------------------------------------
1043 !
1044 ELSE
1045 !
1046 !----------------------------------------------------------------------
1047 !*** TURBULENCE
1048 !----------------------------------------------------------------------
1049 !----------------------------------------------------------------------
1050 !*** COEFFICIENTS OF THE TERMS IN THE NUMERATOR
1051 !----------------------------------------------------------------------
1052 !
1053 ANUM=(ANMM*GML+ANMH*GHL)*GHL
1054 BNUM= BNMM*GML+BNMH*GHL
1055 !
1056 !----------------------------------------------------------------------
1057 !*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
1058 !----------------------------------------------------------------------
1059 !
1060 ADEN=(ADNM*GML+ADNH*GHL)*GHL
1061 BDEN= BDNM*GML+BDNH*GHL
1062 CDEN= 1.
1063 !
1064 !----------------------------------------------------------------------
1065 !*** COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.
1066 !----------------------------------------------------------------------
1067 !
1068 ARHS=-(ANUM*BDEN-BNUM*ADEN)*2.
1069 BRHS=- ANUM*4.
1070 CRHS=- BNUM*2.
1071 !
1072 !----------------------------------------------------------------------
1073 !*** INITIAL VALUE OF L/Q
1074 !----------------------------------------------------------------------
1075 !
1076 DLOQ1=EL(K)/SQRT(Q2(K))
1077 !
1078 !----------------------------------------------------------------------
1079 !*** FIRST ITERATION FOR L/Q, RHS=0
1080 !----------------------------------------------------------------------
1081 !
1082 ELOQ21=1./EQOL2
1083 ELOQ11=SQRT(ELOQ21)
1084 ELOQ31=ELOQ21*ELOQ11
1085 ELOQ41=ELOQ21*ELOQ21
1086 ELOQ51=ELOQ21*ELOQ31
1087 !
1088 !----------------------------------------------------------------------
1089 !*** 1./DENOMINATOR
1090 !----------------------------------------------------------------------
1091 !
1092 RDEN1=1./(ADEN*ELOQ41+BDEN*ELOQ21+CDEN)
1093 !
1094 !----------------------------------------------------------------------
1095 !*** D(RHS)/D(L/Q)
1096 !----------------------------------------------------------------------
1097 !
1098 RHSP1=(ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1
1099 !
1100 !----------------------------------------------------------------------
1101 !*** FIRST-GUESS SOLUTION
1102 !----------------------------------------------------------------------
1103 !
1104 ELOQ12=ELOQ11+(DLOQ1-ELOQ11)*EXP(RHSP1*DTTURBL)
1105 ELOQ12=MAX(ELOQ12,EPS1)
1106 !
1107 !----------------------------------------------------------------------
1108 !*** SECOND ITERATION FOR L/Q
1109 !----------------------------------------------------------------------
1110 !
1111 ELOQ22=ELOQ12*ELOQ12
1112 ELOQ32=ELOQ22*ELOQ12
1113 ELOQ42=ELOQ22*ELOQ22
1114 ELOQ52=ELOQ22*ELOQ32
1115 !
1116 !----------------------------------------------------------------------
1117 !*** 1./DENOMINATOR
1118 !----------------------------------------------------------------------
1119 !
1120 RDEN2=1./(ADEN*ELOQ42+BDEN*ELOQ22+CDEN)
1121 RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1
1122 RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2
1123 RHST2=RHS2/RHSP2
1124 !
1125 !----------------------------------------------------------------------
1126 !*** CORRECTED SOLUTION
1127 !----------------------------------------------------------------------
1128 !
1129 ELOQ13=ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*EXP(RHSP2*DTTURBL)
1130 ELOQ13=AMAX1(ELOQ13,EPS1)
1131 !
1132 !----------------------------------------------------------------------
1133 !*** TWO ITERATIONS IS ENOUGH IN MOST CASES ...
1134 !----------------------------------------------------------------------
1135 !
1136 ELOQN=ELOQ13
1137 !
1138 IF(ELOQN>EPS1)THEN
1139 Q2(K)=EL(K)*EL(K)/(ELOQN*ELOQN)
1140 Q2(K)=AMAX1(Q2(K),EPSQ2)
1141 IF(Q2(K)==EPSQ2)THEN
1142 EL(K)=EPSL
1143 ENDIF
1144 ELSE
1145 Q2(K)=EPSQ2
1146 EL(K)=EPSL
1147 ENDIF
1148 !
1149 !----------------------------------------------------------------------
1150 !*** END OF TURBULENT BRANCH
1151 !----------------------------------------------------------------------
1152 !
1153 ENDIF
1154 !----------------------------------------------------------------------
1155 !*** END OF PRODUCTION/DISSIPATION LOOP
1156 !----------------------------------------------------------------------
1157 !
1158 ENDDO main_integration
1159 !
1160 !----------------------------------------------------------------------
1161 !*** LOWER BOUNDARY CONDITION FOR Q2
1162 !----------------------------------------------------------------------
1163 !
1164 Q2(LMH)=AMAX1(B1**(2./3.)*USTAR*USTAR,EPSQ2)
1165 !----------------------------------------------------------------------
1166 !
1167 END SUBROUTINE PRODQ2
1168 !
1169 !----------------------------------------------------------------------
1170 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1171 !----------------------------------------------------------------------
1172 SUBROUTINE DIFCOF &
1173 ! ******************************************************************
1174 ! * *
1175 ! * LEVEL 2.5 DIFFUSION COEFFICIENTS *
1176 ! * *
1177 ! ******************************************************************
1178 &(LMH,LMXL,GM,GH,EL,T,Q2,Z,AKM,AKH &
1179 &,IDS,IDE,JDS,JDE,KDS,KDE &
1180 &,IMS,IME,JMS,JME,KMS,KME &
1181 &,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG) ! debug
1182 !----------------------------------------------------------------------
1183 !
1184 IMPLICIT NONE
1185 !
1186 !----------------------------------------------------------------------
1187 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1188 & ,IMS,IME,JMS,JME,KMS,KME &
1189 & ,ITS,ITE,JTS,JTE,KTS,KTE
1190 !
1191 INTEGER,INTENT(IN) :: LMH,LMXL
1192 !
1193 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: Q2,T
1194 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL,GH,GM
1195 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1196 !
1197 REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: AKH,AKM
1198 !----------------------------------------------------------------------
1199 !***
1200 !*** LOCAL VARIABLES
1201 !***
1202 INTEGER :: K,KINV
1203 !
1204 REAL :: ADEN,AKMIN,BDEN,BESH,BESM,CDEN,D2T,ELL,ELOQ2,ELOQ4,ELQDZ &
1205 & ,ESH,ESM,GHL,GML,Q1L,RDEN,RDZ
1206 !
1207 !*** Begin debugging
1208 INTEGER,INTENT(IN) :: PRINT_DIAG
1209 ! REAL :: D2Tmin
1210 !*** End debugging
1211 !
1212 !----------------------------------------------------------------------
1213 !**********************************************************************
1214 !----------------------------------------------------------------------
1215 !
1216 DO K=1,LMH-1
1217 ELL=EL(K)
1218 !
1219 ELOQ2=ELL*ELL/Q2(K)
1220 ELOQ4=ELOQ2*ELOQ2
1221 !
1222 GML=GM(K)
1223 GHL=GH(K)
1224 !
1225 !----------------------------------------------------------------------
1226 !*** COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
1227 !----------------------------------------------------------------------
1228 !
1229 ADEN=(ADNM*GML+ADNH*GHL)*GHL
1230 BDEN= BDNM*GML+BDNH*GHL
1231 CDEN= 1.
1232 !
1233 !----------------------------------------------------------------------
1234 !*** COEFFICIENTS FOR THE SM DETERMINANT
1235 !----------------------------------------------------------------------
1236 !
1237 BESM=BSMH*GHL
1238 !
1239 !----------------------------------------------------------------------
1240 !*** COEFFICIENTS FOR THE SH DETERMINANT
1241 !----------------------------------------------------------------------
1242 !
1243 BESH=BSHM*GML+BSHH*GHL
1244 !
1245 !----------------------------------------------------------------------
1246 !*** 1./DENOMINATOR
1247 !----------------------------------------------------------------------
1248 !
1249 RDEN=1./(ADEN*ELOQ4+BDEN*ELOQ2+CDEN)
1250 !
1251 !----------------------------------------------------------------------
1252 !*** SM AND SH
1253 !----------------------------------------------------------------------
1254 !
1255 ESM=(BESM*ELOQ2+CESM)*RDEN
1256 ESH=(BESH*ELOQ2+CESH)*RDEN
1257 !
1258 !----------------------------------------------------------------------
1259 !*** DIFFUSION COEFFICIENTS
1260 !----------------------------------------------------------------------
1261 !
1262 RDZ=2./(Z(K)-Z(K+2))
1263 Q1L=SQRT(Q2(K))
1264 ELQDZ=ELL*Q1L*RDZ
1265 AKM(K)=ELQDZ*ESM
1266 AKH(K)=ELQDZ*ESH
1267 !----------------------------------------------------------------------
1268 ENDDO
1269 !----------------------------------------------------------------------
1270 !
1271 !----------------------------------------------------------------------
1272 !*** INVERSIONS
1273 !----------------------------------------------------------------------
1274 !
1275 ! IF(LMXL==LMH)THEN
1276 ! KINV=LMH
1277 ! D2Tmin=0.
1278 !
1279 ! DO K=LMH/2,LMH-1
1280 ! D2T=T(K-1)-2.*T(K)+T(K+1)
1281 ! IF(D2T<D2Tmin)THEN
1282 ! D2Tmin=D2T
1283 ! IF(D2T<0)KINV=K
1284 ! ENDIF
1285 ! ENDDO
1286 !
1287 ! IF(KINV<LMH)THEN
1288 ! DO K=KINV-1,LMH-1
1289 ! RDZ=2./(Z(K)-Z(K+2))
1290 ! AKMIN=0.5*RDZ
1291 ! AKM(K)=MAX(AKM(K),AKMIN)
1292 ! AKH(K)=MAX(AKH(K),AKMIN)
1293 ! ENDDO
1294 !
1295 !*** Begin debugging
1296 ! IF(PRINT_DIAG>0)THEN
1297 ! write(6,"(a,3i3)") '{turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
1298 ! write(6,"(a,3i3)") '}turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
1299 ! IF(PRINT_DIAG==1)THEN
1300 ! write(6,"(a)") &
1301 ! '{turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
1302 ! ELSE
1303 ! write(6,"(a)") &
1304 ! '}turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
1305 ! ENDIF
1306 ! DO K=LMH-1,KINV-1,-1
1307 ! D2T=T(K-1)-2.*T(K)+T(K+1)
1308 ! RDZ=2./(Z(K)-Z(K+2))
1309 ! AKMIN=0.5*RDZ
1310 ! IF(PRINT_DIAG==1)THEN
1311 ! write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '{turb3 ' &
1312 ! ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
1313 ! ELSE
1314 ! write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '}turb3 ' &
1315 ! ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
1316 ! ENDIF
1317 ! ENDDO
1318 ! ENDIF !- IF (print_diag > 0) THEN
1319 ! ENDIF !- IF(KINV<LMH)THEN
1320 !*** End debugging
1321 !
1322 ! ENDIF
1323 !----------------------------------------------------------------------
1324 !
1325 END SUBROUTINE DIFCOF
1326 !
1327 !----------------------------------------------------------------------
1328 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1329 !----------------------------------------------------------------------
1330 SUBROUTINE VDIFQ &
1331 ! ******************************************************************
1332 ! * *
1333 ! * VERTICAL DIFFUSION OF Q2 (TKE) *
1334 ! * *
1335 ! ******************************************************************
1336 &(LMH,DTDIF,Q2,EL,Z &
1337 &,IDS,IDE,JDS,JDE,KDS,KDE &
1338 &,IMS,IME,JMS,JME,KMS,KME &
1339 &,ITS,ITE,JTS,JTE,KTS,KTE)
1340 !----------------------------------------------------------------------
1341 !
1342 IMPLICIT NONE
1343 !
1344 !----------------------------------------------------------------------
1345 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1346 & ,IMS,IME,JMS,JME,KMS,KME &
1347 & ,ITS,ITE,JTS,JTE,KTS,KTE
1348 !
1349 INTEGER,INTENT(IN) :: LMH
1350 !
1351 REAL,INTENT(IN) :: DTDIF
1352 !
1353 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL
1354 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1355 !
1356 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
1357 !----------------------------------------------------------------------
1358 !***
1359 !*** LOCAL VARIABLES
1360 !***
1361 INTEGER :: K
1362 !
1363 REAL :: ADEN,AKQS,BDEN,BESH,BESM,CDEN,CF,DTOZS,ELL,ELOQ2,ELOQ4 &
1364 & ,ELQDZ,ESH,ESM,ESQHF,GHL,GML,Q1L,RDEN,RDZ
1365 !
1366 REAL,DIMENSION(KTS:KTE-2) :: AKQ,CM,CR,DTOZ,RSQ2
1367 !----------------------------------------------------------------------
1368 !**********************************************************************
1369 !----------------------------------------------------------------------
1370 !***
1371 !*** VERTICAL TURBULENT DIFFUSION
1372 !***
1373 !----------------------------------------------------------------------
1374 ESQHF=0.5*ESQ
1375 !
1376 DO K=KTS,LMH-2
1377 DTOZ(K)=(DTDIF+DTDIF)/(Z(K)-Z(K+2))
1378 AKQ(K)=SQRT((Q2(K)+Q2(K+1))*0.5)*(EL(K)+EL(K+1))*ESQHF &
1379 & /(Z(K+1)-Z(K+2))
1380 CR(K)=-DTOZ(K)*AKQ(K)
1381 ENDDO
1382 !
1383 CM(1)=DTOZ(1)*AKQ(1)+1.
1384 RSQ2(1)=Q2(1)
1385 !
1386 DO K=KTS+1,LMH-2
1387 CF=-DTOZ(K)*AKQ(K-1)/CM(K-1)
1388 CM(K)=-CR(K-1)*CF+(AKQ(K-1)+AKQ(K))*DTOZ(K)+1.
1389 RSQ2(K)=-RSQ2(K-1)*CF+Q2(K)
1390 ENDDO
1391 !
1392 DTOZS=(DTDIF+DTDIF)/(Z(LMH-1)-Z(LMH+1))
1393 AKQS=SQRT((Q2(LMH-1)+Q2(LMH))*0.5)*(EL(LMH-1)+ELZ0)*ESQHF &
1394 & /(Z(LMH)-Z(LMH+1))
1395 !
1396 CF=-DTOZS*AKQ(LMH-2)/CM(LMH-2)
1397 !
1398 Q2(LMH-1)=(DTOZS*AKQS*Q2(LMH)-RSQ2(LMH-2)*CF+Q2(LMH-1)) &
1399 & /((AKQ(LMH-2)+AKQS)*DTOZS-CR(LMH-2)*CF+1.)
1400 !
1401 DO K=LMH-2,KTS,-1
1402 Q2(K)=(-CR(K)*Q2(K+1)+RSQ2(K))/CM(K)
1403 ENDDO
1404 !----------------------------------------------------------------------
1405 !
1406 END SUBROUTINE VDIFQ
1407 !
1408 !----------------------------------------------------------------------
1409 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1410 !---------------------------------------------------------------------
1411 SUBROUTINE VDIFH(DTDIF,LMH,MSS,MSE,LPBL,SZ0 &
1412 & ,RKHS,CLOW,CTS &
1413 & ,SPECIES,NSPEC,RKH,ZHK,RHO &
1414 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1415 & ,IMS,IME,JMS,JME,KMS,KME &
1416 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1417 ! ***************************************************************
1418 ! * *
1419 ! * VERTICAL DIFFUSION OF MASS VARIABLES *
1420 ! * *
1421 ! * Z. JANJIC MODIFIED TO WORK WITH ARBITRARY NUMBER OF SPECIES *
1422 ! * ON 05/17/2008 *
1423 ! * *
1424 ! ***************************************************************
1425 !---------------------------------------------------------------------
1426 !
1427 IMPLICIT NONE
1428 !
1429 !---------------------------------------------------------------------
1430 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1431 & ,IMS,IME,JMS,JME,KMS,KME &
1432 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1433 !
1434 INTEGER,INTENT(IN):: LMH,LPBL,MSS,MSE,NSPEC !BSF
1435 !
1436 ! REAL,INTENT(IN) :: CHKLOWQ,CT,DTDIF,QZ0,RKHS,THZ0
1437 !
1438 REAL,INTENT(IN):: DTDIF,RKHS
1439 !
1440 REAL,DIMENSION(MSS:MSE),INTENT(IN):: CLOW,CTS,SZ0
1441 !
1442 REAL,DIMENSION(KTS:KTE-1),INTENT(IN):: RKH
1443 !
1444 REAL,DIMENSION(KTS:KTE ),INTENT(IN):: RHO
1445 !
1446 REAL,DIMENSION(KTS:KTE+1),INTENT(IN):: ZHK
1447 !
1448 REAL,DIMENSION(MSS:MSE,KTS:KTE),INTENT(INOUT):: SPECIES
1450 !----------------------------------------------------------------------
1451 !***
1452 !*** LOCAL VARIABLES
1453 !***
1454 INTEGER:: K,M
1455 !
1456 REAL:: CF,CMB,CMSB,DTOZL,DTOZS,RCML,RHOK,RKHH,RKHZ,RKSS,RSSB
1457 !
1458 REAL,DIMENSION(KTS:KTE-1):: CM,CR,DTOZ
1459 !
1460 REAL,DIMENSION(MSS:MSE,KTS:KTE-1):: RKCT,RSS
1461 !
1462 !----------------------------------------------------------------------
1463 !**********************************************************************
1464 !---PREPARATIONS-------------------------------------------------------
1465 DO K=KTS,LMH-1
1466 DTOZ(K)=DTDIF/(ZHK(K)-ZHK(K+1))
1467 CR(K)=-DTOZ(K)*RKH(K)
1468 IF(K.LT.LPBL) THEN
1469 DO M=MSS,NSPEC
1470 RKCT(M,K)=0.
1471 ENDDO
1472 ELSE
1473 RKHZ=RKH(K)*(ZHK(K)-ZHK(K+2))
1474 DO M=MSS,NSPEC
1475 RKCT(M,K)=RKHZ*CTS(M)*0.5
1476 ENDDO
1477 ENDIF
1478 ENDDO
1479 !---TOP LEVEL----------------------------------------------------------
1480 RHOK=RHO(KTS)
1481 CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHOK
1482 DO M=MSS,NSPEC
1483 RSS(M,KTS)=-RKCT(M,KTS)*DTOZ(KTS)+SPECIES(M,KTS)*RHOK
1484 ENDDO
1485 !---INTERMEDIATE LEVELS------------------------------------------------
1486 DO K=KTS+1,LMH-1
1487 DTOZL=DTOZ(K)
1488 CF=-DTOZL*RKH(K-1)/CM(K-1)
1489 RHOK=RHO(K)
1490 CM(K)=-CR(K-1)*CF+(RKH(K-1)+RKH(K))*DTOZL+RHOK
1491 DO M=MSS,NSPEC
1492 RSS(M,K)=-RSS(M,K-1)*CF &
1493 +(RKCT(M,K-1)-RKCT(M,K))*DTOZL+SPECIES(M,K)*RHOK
1494 ENDDO
1495 ENDDO
1496 !---BOTTOM LEVEL-------------------------------------------------------
1497 DTOZS=DTDIF/(ZHK(LMH)-ZHK(LMH+1))
1498 RKHH=RKH(LMH-1)
1499 !
1500 CF=-DTOZS*RKHH/CM(LMH-1)
1501 CMB=CR(LMH-1)*CF
1502 RHOK=RHO(LMH)
1503 !
1504 DO M=MSS,NSPEC
1505 RKSS=RKHS*CLOW(M)
1506 CMSB=-CMB+(RKHH+RKSS)*DTOZS+RHOK
1507 RSSB=-RSS(M,LMH-1)*CF+RKCT(M,LMH-1)*DTOZS+SPECIES(M,LMH)*RHOK
1508 SPECIES(M,LMH)=(DTOZS*RKSS*SZ0(M)+RSSB)/CMSB
1509 ENDDO
1510 !---BACKSUBSTITUTION---------------------------------------------------
1511 DO K=LMH-1,KTS,-1
1512 RCML=1./CM(K)
1513 DO M=MSS,NSPEC
1514 SPECIES(M,K)=(-CR(K)*SPECIES(M,K+1)+RSS(M,K))*RCML
1515 ENDDO
1516 ENDDO
1517 !----------------------------------------------------------------------
1518 !
1519 END SUBROUTINE VDIFH
1520 !
1521 !----------------------------------------------------------------------
1522 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1523 !---------------------------------------------------------------------
1524 SUBROUTINE VDIFX(DTDIF,LMH,RKHS,CT &
1525 ,DUST1,DUST2,RKH,Z,RHO &
1526 ,IDS,IDE,JDS,JDE,KDS,KDE &
1527 ,IMS,IME,JMS,JME,KMS,KME &
1528 ,ITS,ITE,JTS,JTE,KTS,KTE)
1529 ! ***************************************************************
1530 ! * *
1531 ! * VERTICAL DIFFUSION OF MASS VARIABLES *
1532 ! * *
1533 ! ***************************************************************
1534 !---------------------------------------------------------------------
1535 !
1536 IMPLICIT NONE
1537 !
1538 !---------------------------------------------------------------------
1539 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1540 ,IMS,IME,JMS,JME,KMS,KME &
1541 ,ITS,ITE,JTS,JTE,KTS,KTE
1542 !
1543 INTEGER,INTENT(IN) :: LMH
1544 !
1545 REAL,INTENT(IN) :: CT,DTDIF,RKHS
1546 !
1547 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKH
1548 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1549 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1550 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: DUST1,DUST2
1551 !
1552 !----------------------------------------------------------------------
1553 !***
1554 !*** LOCAL VARIABLES
1555 !***
1556 INTEGER :: K
1557 !
1558 REAL :: CF,CMB,CMDB,CTHF,DTOZL,DTOZS &
1559 ,RCML,RKHH,RSD1B,RSD2B
1560 !
1561 REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RKCT,RSD1,RSD2
1562 !
1563 !----------------------------------------------------------------------
1564 !**********************************************************************
1565 !----------------------------------------------------------------------
1566 CTHF=0.5*CT
1567 !
1568 DO K=KTS,LMH-1
1569 DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1570 CR(K)=-DTOZ(K)*RKH(K)
1571 RKCT(K)=RKH(K)*(Z(K)-Z(K+2))*CTHF
1572 ENDDO
1573 !
1574 CM(KTS)=DTOZ(KTS)*RKH(KTS)+RHO(KTS)
1575 !----------------------------------------------------------------------
1576 RSD1(KTS)=DUST1(KTS)*RHO(KTS)
1577 RSD2(KTS)=DUST2(KTS)*RHO(KTS)
1578 !----------------------------------------------------------------------
1579 DO K=KTS+1,LMH-1
1580 DTOZL=DTOZ(K)
1581 CF=-DTOZL*RKH(K-1)/CM(K-1)
1582 CM(K)=-CR(K-1)*CF+(RKH(K-1)+RKH(K))*DTOZL+RHO(K)
1583 RSD1(K)=-RSD1(K-1)*CF+DUST1(K)*RHO(K)
1584 RSD2(K)=-RSD2(K-1)*CF+DUST2(K)*RHO(K)
1585 ENDDO
1586 !
1587 DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1588 RKHH=RKH(LMH-1)
1589 !
1590 CF=-DTOZS*RKHH/CM(LMH-1)
1591 !
1592 CMB=CR(LMH-1)*CF
1593 CMDB=-CMB+RKHH*DTOZS+RHO(LMH)
1594 !
1595 RSD1B=-RSD1(LMH-1)*CF+DUST1(LMH)*RHO(LMH)
1596 RSD2B=-RSD2(LMH-1)*CF+DUST2(LMH)*RHO(LMH)
1597 !----------------------------------------------------------------------
1598 DUST1(LMH)=RSD1B/CMDB
1599 DUST2(LMH)=RSD2B/CMDB
1600 !----------------------------------------------------------------------
1601 DO K=LMH-1,KTS,-1
1602 RCML=1./CM(K)
1603 DUST1(K)=(-CR(K)*DUST1(K+1)+RSD1(K))*RCML
1604 DUST2(K)=(-CR(K)*DUST2(K+1)+RSD2(K))*RCML
1605 ENDDO
1606 !----------------------------------------------------------------------
1607 !
1608 END SUBROUTINE VDIFX
1610 !---------------------------------------------------------------------
1611 !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
1612 !---------------------------------------------------------------------
1613 SUBROUTINE VDIFV(LMH,DTDIF,UZ0,VZ0,RKMS,U,V,RKM,Z,RHO &
1614 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1615 & ,IMS,IME,JMS,JME,KMS,KME &
1616 ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
1617 ! ***************************************************************
1618 ! * *
1619 ! * VERTICAL DIFFUSION OF VELOCITY COMPONENTS *
1620 ! * *
1621 ! ***************************************************************
1622 !---------------------------------------------------------------------
1623 !
1624 IMPLICIT NONE
1625 !
1626 !---------------------------------------------------------------------
1627 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1628 & ,IMS,IME,JMS,JME,KMS,KME &
1629 & ,ITS,ITE,JTS,JTE,KTS,KTE,I,J
1630 !
1631 INTEGER,INTENT(IN) :: LMH
1632 !
1633 REAL,INTENT(IN) :: RKMS,DTDIF,UZ0,VZ0
1634 !
1635 REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: RKM
1636 REAL,DIMENSION(KTS:KTE),INTENT(IN) :: RHO
1637 REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
1638 !
1639 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: U,V
1640 !----------------------------------------------------------------------
1641 !***
1642 !*** LOCAL VARIABLES
1643 !***
1644 INTEGER :: K
1645 !
1646 REAL :: CF,DTOZAK,DTOZL,DTOZS,RCML,RCMVB,RHOK,RKMH
1647 !
1648 REAL,DIMENSION(KTS:KTE-1) :: CM,CR,DTOZ,RSU,RSV
1649 !----------------------------------------------------------------------
1650 !**********************************************************************
1651 !----------------------------------------------------------------------
1652 DO K=1,LMH-1
1653 DTOZ(K)=DTDIF/(Z(K)-Z(K+1))
1654 CR(K)=-DTOZ(K)*RKM(K)
1655 ENDDO
1656 !
1657 RHOK=RHO(1)
1658 CM(1)=DTOZ(1)*RKM(1)+RHOK
1659 RSU(1)=U(1)*RHOK
1660 RSV(1)=V(1)*RHOK
1661 !----------------------------------------------------------------------
1662 DO K=2,LMH-1
1663 DTOZL=DTOZ(K)
1664 CF=-DTOZL*RKM(K-1)/CM(K-1)
1665 RHOK=RHO(K)
1666 CM(K)=-CR(K-1)*CF+(RKM(K-1)+RKM(K))*DTOZL+RHOK
1667 RSU(K)=-RSU(K-1)*CF+U(K)*RHOK
1668 RSV(K)=-RSV(K-1)*CF+V(K)*RHOK
1669 ENDDO
1670 !----------------------------------------------------------------------
1671 DTOZS=DTDIF/(Z(LMH)-Z(LMH+1))
1672 RKMH=RKM(LMH-1)
1673 !
1674 CF=-DTOZS*RKMH/CM(LMH-1)
1675 RHOK=RHO(LMH)
1676 RCMVB=1./((RKMH+RKMS)*DTOZS-CR(LMH-1)*CF+RHOK)
1677 DTOZAK=DTOZS*RKMS
1678 !----------------------------------------------------------------------
1679 U(LMH)=(DTOZAK*UZ0-RSU(LMH-1)*CF+U(LMH)*RHOK)*RCMVB
1680 V(LMH)=(DTOZAK*VZ0-RSV(LMH-1)*CF+V(LMH)*RHOK)*RCMVB
1681 !----------------------------------------------------------------------
1682 DO K=LMH-1,1,-1
1683 RCML=1./CM(K)
1684 U(K)=(-CR(K)*U(K+1)+RSU(K))*RCML
1685 V(K)=(-CR(K)*V(K+1)+RSV(K))*RCML
1686 ENDDO
1687 !----------------------------------------------------------------------
1688 !
1689 END SUBROUTINE VDIFV
1690 !
1691 !-----------------------------------------------------------------------
1692 !
1693 !=======================================================================
1694 SUBROUTINE MYJPBLINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN, &
1695 & TKE_MYJ,EXCH_H,RESTART,ALLOWED_TO_READ, &
1696 & IDS,IDE,JDS,JDE,KDS,KDE, &
1697 & IMS,IME,JMS,JME,KMS,KME, &
1698 & ITS,ITE,JTS,JTE,KTS,KTE )
1699 !-----------------------------------------------------------------------
1700 IMPLICIT NONE
1701 !-----------------------------------------------------------------------
1702 LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
1703 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
1704 & IMS,IME,JMS,JME,KMS,KME, &
1705 & ITS,ITE,JTS,JTE,KTS,KTE
1707 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
1708 & RUBLTEN, &
1709 & RVBLTEN, &
1710 & RTHBLTEN, &
1711 & RQVBLTEN, &
1712 & TKE_MYJ
1713 INTEGER :: I,J,K,ITF,JTF,KTF
1714 !-----------------------------------------------------------------------
1715 !-----------------------------------------------------------------------
1717 JTF=MIN0(JTE,JDE-1)
1718 KTF=MIN0(KTE,KDE-1)
1719 ITF=MIN0(ITE,IDE-1)
1721 IF(.NOT.RESTART)THEN
1722 DO J=JTS,JTF
1723 DO K=KTS,KTF
1724 DO I=ITS,ITF
1725 TKE_MYJ(I,K,J)=EPSQ2
1726 RUBLTEN(I,K,J)=0.
1727 RVBLTEN(I,K,J)=0.
1728 RTHBLTEN(I,K,J)=0.
1729 RQVBLTEN(I,K,J)=0.
1730 EXCH_H(I,K,J)=0.
1731 ENDDO
1732 ENDDO
1733 ENDDO
1734 ENDIF
1736 END SUBROUTINE MYJPBLINIT
1737 !-----------------------------------------------------------------------
1738 !
1739 END MODULE MODULE_BL_MYJPBL
1740 !
1741 !-----------------------------------------------------------------------