source: lmdz_wrf/WRFV3/phys/module_bl_myjurb.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

  • Property svn:executable set to *
File size: 67.3 KB
Line 
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
158
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
196
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
216
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!
257     
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)
338
339            DZ_KFLIP(K)=DZ(I,KFLIP,J)
340 
341
342
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
352
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!
469
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)))
526
527          ENDDO
528!            SFK(1)=1.
529          if(flag_bep)then
530           ! Multi-layer UCM
531           
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!
599         
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!
627         
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!
694
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!
703
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)
709
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
734
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!
1560
1561   
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
1652
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
1737     
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!
1750     
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
1768
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!-----------------------------------------------------------------------
1778
1779      JTF=MIN0(JTE,JDE-1)
1780      KTF=MIN0(KTE,KDE-1)
1781      ITF=MIN0(ITE,IDE-1)
1782
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
1797
1798      END SUBROUTINE MYJURBINIT
1799!-----------------------------------------------------------------------
1800!
1801      END MODULE MODULE_BL_MYJURB
1802!
1803!-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.