PageRenderTime 70ms CodeModel.GetById 27ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_bl_myjurb.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1803 lines | 874 code | 27 blank | 902 comment | 4 complexity | 6ef59fc7c9dc406dd4af0a2324195dc7 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. !
    2. 

  2.       MODULE MODULE_BL_MYJURB
    3. 

  3. !
    4. 

  4. !-----------------------------------------------------------------------
    5. 

  5. !
    6. 

  6.       USE MODULE_MODEL_CONSTANTS
    7. 

  7. !
    8. 

  8. !-----------------------------------------------------------------------
    9. 

  9. !
    10. 

  10. ! REFERENCES:  Janjic (2002), NCEP Office Note 437
    11. 

  11. !              Mellor and Yamada (1982), Rev. Geophys. Space Phys.
    12. 

  12. !
    13. 

  13. ! ABSTRACT:
    14. 

  14. !     MYJ UPDATES THE TURBULENT KINETIC ENERGY WITH THE PRODUCTION/
    15. 

  15. !     DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
    16. 

  16. !     (USING AN IMPLICIT FORMULATION) FROM MELLOR-YAMADA
    17. 

  17. !     LEVEL 2.5 AS EXTENDED BY JANJIC.  EXCHANGE COEFFICIENTS FOR
    18. 

  18. !     THE SURFACE AND FOR ALL LAYER INTERFACES ARE COMPUTED FROM
    19. 

  19. !     MONIN-OBUKHOV THEORY.
    20. 

  20. !     THE TURBULENT VERTICAL EXCHANGE IS THEN EXECUTED.
    21. 

  21. !
    22. 

  22. !-----------------------------------------------------------------------
    23. 

  23. !
    24. 

  24.       INTEGER :: ITRMX=5 ! Iteration count for mixing length computation
    25. 

  25. !
    26. 

  26. !     REAL,PARAMETER :: G=9.81,PI=3.1415926,R_D=287.04,R_V=461.6        &
    27. 

  27. !    &                 ,VKARMAN=0.4
    28. 

  28.       REAL,PARAMETER :: PI=3.1415926,VKARMAN=0.4
    29. 

  29. !     REAL,PARAMETER :: CP=7.*R_D/2.
    30. 

  30.       REAL,PARAMETER :: CAPA=R_D/CP
    31. 

  31.       REAL,PARAMETER :: RLIVWV=XLS/XLV,ELOCP=2.72E6/CP
    32. 

  32.       REAL,PARAMETER :: EPS1=1.E-12,EPS2=0.
    33. 

  33.       REAL,PARAMETER :: EPSL=0.32,EPSRU=1.E-7,EPSRS=1.E-7               &
    34. 

  34.      &                 ,EPSTRB=1.E-24
    35. 

  35.       REAL,PARAMETER :: EPSA=1.E-8,EPSIT=1.E-4,EPSU2=1.E-4,EPSUST=0.07  &
    36. 

  36.      &                 ,FH=1.01 
    37. 

  37.       REAL,PARAMETER :: ALPH=0.30,BETA=1./273.,EL0MAX=1000.,EL0MIN=1.   &
    38. 

  38.      &                 ,ELFC=0.23*0.5,GAM1=0.2222222222222222222        &
    39. 

  39.      &                 ,PRT=1.
    40. 

  40.       REAL,PARAMETER :: A1=0.659888514560862645                         &
    41. 

  41.      &                 ,A2x=0.6574209922667784586                       &
    42. 

  42.      &                 ,B1=11.87799326209552761                         &
    43. 

  43.      &                 ,B2=7.226971804046074028                         &
    44. 

  44.      &                 ,C1=0.000830955950095854396
    45. 

  45.       REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86
    46. 

  46.       REAL,PARAMETER :: ELZ0=0.,ESQ=5.0,EXCM=0.001                      &
    47. 

  47.      &                 ,FHNEU=0.8,GLKBR=10.,GLKBS=30.                   &
    48. 

  48.      &                 ,QVISC=2.1E-5,RFC=0.191,RIC=0.505,SMALL=0.35     &
    49. 

  49.      &                 ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5   &
    50. 

  50.      &                 ,USTC=0.7,USTR=0.225,VISC=1.5E-5                 &
    51. 

  51.      &                 ,WOLD=0.15,WWST=1.2,ZTMAX=1.,ZTFC=1.,ZTMIN=-5.
    52. 

  52. !
    53. 

  53.       REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC
    54. 

  54. !
    55. 

  55.       REAL,PARAMETER :: BTG=BETA*G,CZIV=SMALL*GLKBS                     &
    56. 

  56. !    &                 ,EP_1=R_V/R_D-1.,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS  &
    57. 

  57.      &                 ,ESQHF=0.5*5.0,GRRS=GLKBR/GLKBS                  &
    58. 

  58.      &                 ,RB1=1./B1,RTVISC=1./TVISC,RVISC=1./VISC         &
    59. 

  59.      &                 ,ZQRZT=SQSC/SQPR
    60. 

  60. !
    61. 

  61.       REAL,PARAMETER :: ADNH= 9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG      &                  
    62. 

  62.      &                 ,ADNM=18.*A1*A1*A2x*(B2-3.*A2x)*BTG              & 
    63. 

  63.      &                 ,ANMH=-9.*A1*A2x*A2x*BTG*BTG                     &
    64. 

  64.      &                 ,ANMM=-3.*A1*A2x*(3.*A2x+3.*B2*C1+18.*A1*C1-B2)  &
    65. 

  65.      &                                *BTG                              &   
    66. 

  66.      &                 ,BDNH= 3.*A2x*(7.*A1+B2)*BTG                     &
    67. 

  67.      &                 ,BDNM= 6.*A1*A1                                  &
    68. 

  68.      &                 ,BEQH= A2x*B1*BTG+3.*A2x*(7.*A1+B2)*BTG          &
    69. 

  69.      &                 ,BEQM=-A1*B1*(1.-3.*C1)+6.*A1*A1                 &
    70. 

  70.      &                 ,BNMH=-A2x*BTG                                   &     
    71. 

  71.      &                 ,BNMM=A1*(1.-3.*C1)                              &
    72. 

  72.      &                 ,BSHH=9.*A1*A2x*A2x*BTG                          &
    73. 

  73.      &                 ,BSHM=18.*A1*A1*A2x*C1                           &
    74. 

  74.      &                 ,BSMH=-3.*A1*A2x*(3.*A2x+3.*B2*C1+12.*A1*C1-B2)  &
    75. 

  75.      &                                *BTG                              &
    76. 

  76.      &                 ,CESH=A2x                                        &
    77. 

  77.      &                 ,CESM=A1*(1.-3.*C1)                              &
    78. 

  78.      &                 ,CNV=EP_1*G/BTG                                  &
    79. 

  79.      &                 ,ELFCS=VKARMAN*BTG                               &
    80. 

  80.      &                 ,FZQ1=RTVISC*QVISC*ZQRZT                         &
    81. 

  81.      &                 ,FZQ2=RTVISC*QVISC*ZQRZT                         &
    82. 

  82.      &                 ,FZT1=RVISC *TVISC*SQPR                          &
    83. 

  83.      &                 ,FZT2=CZIV*GRRS*TVISC*SQPR                       &
    84. 

  84.      &                 ,FZU1=CZIV*VISC                                  &
    85. 

  85.      &                 ,PIHF=0.5*PI                                     &
    86. 

  86.      &                 ,RFAC=RIC/(FHNEU*RFC*RFC)                        &
    87. 

  87.      &                 ,RQVISC=1./QVISC                                 &
    88. 

  88.      &                 ,RRIC=1./RIC                                     &
    89. 

  89.      &                 ,USTFC=0.018/G                                   &
    90. 

  90.      &                 ,WNEW=1.-WOLD                                    &
    91. 

  91.      &                 ,WWST2=WWST*WWST
    92. 

  92. !
    93. 

  93. !-----------------------------------------------------------------------
    94. 

  94. !***  FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)**2
    95. 

  95. !-----------------------------------------------------------------------
    96. 

  96. !
    97. 

  97.       REAL,PARAMETER :: AEQH=9.*A1*A2x*A2x*B1*BTG*BTG                   &
    98. 

  98.      &                      +9.*A1*A2x*A2x*(12.*A1+3.*B2)*BTG*BTG       &
    99. 

  99.      &                 ,AEQM=3.*A1*A2x*B1*(3.*A2x+3.*B2*C1+18.*A1*C1-B2)&
    100. 

  100.      &                      *BTG+18.*A1*A1*A2x*(B2-3.*A2x)*BTG
    101. 

  101. !
    102. 

  102. !-----------------------------------------------------------------------
    103. 

  103. !***  FORBIDDEN TURBULENCE AREA
    104. 

  104. !-----------------------------------------------------------------------
    105. 

  105. !
    106. 

  106.       REAL,PARAMETER :: REQU=-AEQH/AEQM                                 &
    107. 

  107.      &                 ,EPSGH=1.E-9,EPSGM=REQU*EPSGH
    108. 

  108. !
    109. 

  109. !-----------------------------------------------------------------------
    110. 

  110. !***  NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
    111. 

  111. !-----------------------------------------------------------------------
    112. 

  112. ! 
    113. 

  113.       REAL,PARAMETER :: UBRYL=(18.*REQU*A1*A1*A2x*B2*C1*BTG             &
    114. 

  114.      &                         +9.*A1*A2x*A2x*B2*BTG*BTG)               &
    115. 

  115.      &                        /(REQU*ADNM+ADNH)                         &
    116. 

  116.      &                 ,UBRY=(1.+EPSRS)*UBRYL,UBRY3=3.*UBRY
    117. 

  117. !
    118. 

  118.       REAL,PARAMETER :: AUBH=27.*A1*A2x*A2x*B2*BTG*BTG-ADNH*UBRY3       &
    119. 

  119.      &                 ,AUBM=54.*A1*A1*A2x*B2*C1*BTG -ADNM*UBRY3        &
    120. 

  120.      &                 ,BUBH=(9.*A1*A2x+3.*A2x*B2)*BTG-BDNH*UBRY3       &
    121. 

  121.      &                 ,BUBM=18.*A1*A1*C1           -BDNM*UBRY3         &
    122. 

  122.      &                 ,CUBR=1.                     -     UBRY3         &
    123. 

  123.      &                 ,RCUBR=1./CUBR
    124. 

  124. !
    125. 

  125. !-----------------------------------------------------------------------
    126. 

  126. !
    127. 

  127.       CONTAINS
    128. 

  128. !
    129. 

  129. !----------------------------------------------------------------------
    130. 

  130.       SUBROUTINE MYJURB(IDIFF,FLAG_BEP,DT,STEPBL,HT,DZ                                &
    131. 

  131.      &                 ,PMID,PINT,TH,T,EXNER,QV,CWM,U,V,RHO            &
    132. 

  132.      &                 ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0              &
    133. 

  133.      &                 ,LOWLYR,XLAND,SICE,SNOW                         &
    134. 

  134.      &                 ,TKE_MYJ,EXCH_H,EXCH_M,USTAR,ZNT,EL_MYJ,PBLH,KPBL,CT   &
    135. 

  135.      &                 ,AKHS,AKMS,ELFLX                                &
    136. 

  136.      &                 ,RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN     &
    137. 

  137. ! Variables added for multi-layer UCM
    138. 

  138.      &                 ,FRC_URB2D                                      & 
    139. 

  139.      &                 ,A_U_BEP,A_V_BEP,A_T_BEP,A_Q_BEP                &
    140. 

  140.      &                 ,A_E_BEP,B_U_BEP,B_V_BEP                        &
    141. 

  141.      &                 ,B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP                &
    142. 

  142.      &                 ,DL_U_BEP,SF_BEP,VL_BEP                         &
    143. 

  143. ! Multi-layer UCM end
    144. 

  144.      &                 ,IDS,IDE,JDS,JDE,KDS,KDE                        &
    145. 

  145.      &                 ,IMS,IME,JMS,JME,KMS,KME                        &
    146. 

  146.      &                 ,ITS,ITE,JTS,JTE,KTS,KTE)
    147. 

  147. !----------------------------------------------------------------------
    148. 

  148. !
    149. 

  149.       IMPLICIT NONE
    150. 

  150. !
    151. 

  151. !----------------------------------------------------------------------
    152. 

  152.       INTEGER,INTENT(IN) :: IDIFF
    153. 

  153.       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
    154. 

  154.      &                     ,IMS,IME,JMS,JME,KMS,KME                    &
    155. 

  155.      &                     ,ITS,ITE,JTS,JTE,KTS,KTE
    156. 

  156. !
    157. 

  157.       INTEGER,INTENT(IN) :: STEPBL
    158. 

  158. 

  159.       INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR
    160. 

  160. !
    161. 

  161.       INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: KPBL
    162. 

  162. !
    163. 

  163.       REAL,INTENT(IN) :: DT
    164. 

  164. !
    165. 

  165.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,SICE,SNOW       &
    166. 

  166.      &                                             ,TSK,XLAND
    167. 

  167. !
    168. 

  168.       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: CWM,DZ     &
    169. 

  169.      &                                                     ,EXNER      &
    170. 

  170.      &                                                     ,PMID,PINT  &
    171. 

  171.      &                                                     ,QV,RHO     &
    172. 

  172.      &                                                     ,T,TH,U,V   
    173. 

  173. !
    174. 

  174.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: PBLH
    175. 

  175. !
    176. 

  176.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS
    177. 

  177. !
    178. 

  178.       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME)                          &
    179. 

  179.      &    ,INTENT(OUT) ::                      EL_MYJ                  &
    180. 

  180.      &                                        ,RQCBLTEN,RQVBLTEN       &
    181. 

  181.      &                                        ,RTHBLTEN                &
    182. 

  182.      &                                        ,RUBLTEN,RVBLTEN        
    183. 

  183. ! Variables added for multi-layer UCM
    184. 

  184.       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: A_U_BEP    &
    185. 

  185.      &                                        ,A_V_BEP,A_T_BEP         &
    186. 

  186.      &                                        ,A_E_BEP,B_U_BEP         &
    187. 

  187.      &                                        ,A_Q_BEP,B_Q_BEP         &
    188. 

  188.      &                                        ,B_V_BEP,B_T_BEP         &
    189. 

  189.      &                                        ,B_E_BEP,DLG_BEP         &
    190. 

  190.      &                                        ,DL_U_BEP                &
    191. 

  191.      &                                        ,VL_BEP,SF_BEP
    192. 

  192.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: FRC_URB2D
    193. 

  193.       REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME)                          &
    194. 

  194.      &    ,INTENT(INOUT) ::                    EXCH_H,EXCH_M,TKE_MYJ
    195. 

  195. ! UCM end
    196. 

  196. 

  197.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CT,QSFC,QZ0     &
    198. 

  198.      &                                                ,THZ0,USTAR      &
    199. 

  199.      &                                                ,UZ0,VZ0,ZNT
    200. 

  200. !
    201. 

  201.       REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHKLOWQ,ELFLX
    202. 

  202.       LOGICAL,    INTENT(IN   )    ::     FLAG_BEP
    203. 

  203. !
    204. 

  204. !----------------------------------------------------------------------
    205. 

  205. !***
    206. 

  206. !***  LOCAL VARIABLES
    207. 

  207. !***
    208. 

  208.       INTEGER :: I,J,K,KFLIP,LLOW,LMH,LMXL
    209. 

  209. ! Variables added for multi-layer UCM
    210. 

  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. 

  211.       REAL,DIMENSION(KTS:KTE) :: A_Q,B_Q,VLK,DL_U,B_E_FACE,DZ_KFLIP   !Vertically flipped fields
    212. 

  212.       REAL,DIMENSION(KTS:KTE) :: SFK                             !Vertically flipped SF,VL
    213. 

  213.       INTEGER :: IURB
    214. 

  214.       REAL :: SLOPE,ICEPT
    215. 

  215. ! UCM end
    216. 

  216. 

  217. !
    218. 

  218.       INTEGER,DIMENSION(ITS:ITE,JTS:JTE) :: LPBL
    219. 

  219. !
    220. 

  220.       REAL :: AKHS_DENS,AKMS_DENS,APEX,DCDT,DELTAZ,DQDT,DTDIF,DTDT     &
    221. 

  221.      &       ,DTTURBL,DUDT,DVDT,EXNSFC,PSFC,PTOP,QFC1,QLOW,QOLD        &
    222. 

  222.      &       ,RATIOMX,RDTTURBL,RG,RWMSK,SEAMASK,THNEW,THOLD,TX         &
    223. 

  223.      &       ,ULOW,VLOW,WMSK
    224. 

  224. !
    225. 

  225.       REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,TK,UK,VK
    226. 

  226. !      
    227. 

  227.       REAL,DIMENSION(KTS:KTE-1) :: AKHK,AKMK,EL,GH,GM
    228. 

  228. !
    229. 

  229.       REAL,DIMENSION(KTS:KTE+1) :: ZHK
    230. 

  230. !
    231. 

  231.       REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK
    232. 

  232. !
    233. 

  233.       REAL,DIMENSION(KTS:KTE,ITS:ITE) :: RHOK
    234. 

  234. !
    235. 

  235.       REAL,DIMENSION(ITS:ITE,KTS:KTE,JTS:JTE) :: APE,THE
    236. 

  236. !
    237. 

  237.       REAL,DIMENSION(ITS:ITE,KTS:KTE-1,JTS:JTE) :: AKH,AKM
    238. 

  238. !
    239. 

  239.       REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT
    240. 

  240. !
    241. 

  241. !***  Begin debugging
    242. 

  242.       REAL :: ZSL_DIAG
    243. 

  243.       INTEGER :: IMD,JMD,PRINT_DIAG
    244. 

  244. !***  End debugging
    245. 

  245. !
    246. 

  246. !----------------------------------------------------------------------
    247. 

  247. !**********************************************************************
    248. 

  248. !----------------------------------------------------------------------
    249. 

  249. !
    250. 

  250. !***  Begin debugging
    251. 

  251.       IMD=(IMS+IME)/2
    252. 

  252.       JMD=(JMS+JME)/2
    253. 

  253. !***  End debugging
    254. 

  254. !
    255. 

  255. !***  MAKE PREPARATIONS
    256. 

  256. !
    257. 

  257.      
    258. 

  258. !----------------------------------------------------------------------
    259. 

  259.       DTTURBL=DT*STEPBL
    260. 

  260.       RDTTURBL=1./DTTURBL
    261. 

  261.       DTDIF=DTTURBL
    262. 

  262.       RG=1./G
    263. 

  263. !
    264. 

  264.       DO J=JTS,JTE
    265. 

  265.       DO K=KTS,KTE-1
    266. 

  266.       DO I=ITS,ITE
    267. 

  267.         AKM(I,K,J)=0.
    268. 

  268.       ENDDO
    269. 

  269.       ENDDO
    270. 

  270.       ENDDO
    271. 

  271. !
    272. 

  272.       DO J=JTS,JTE
    273. 

  273.       DO K=KTS,KTE+1
    274. 

  274.       DO I=ITS,ITE
    275. 

  275.         ZINT(I,K,J)=0.
    276. 

  276.       ENDDO
    277. 

  277.       ENDDO
    278. 

  278.       ENDDO
    279. 

  279. !
    280. 

  280.       DO J=JTS,JTE
    281. 

  281.       DO I=ITS,ITE
    282. 

  282.         ZINT(I,KTE+1,J)=HT(I,J)     ! Z at bottom of lowest sigma layer
    283. 

  283. !
    284. 

  284. !!!!!!!!!
    285. 

  285. !!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES
    286. 

  286. !!!!!!!!!
    287. 

  287. !!!!!!  ZINT(I,KTE+1,J)=1.E-4       ! Z of bottom of lowest eta layer
    288. 

  288. !!!!!!  ZHK(KTE+1)=1.E-4            ! Z of bottom of lowest eta layer
    289. 

  289. !
    290. 

  290.       ENDDO
    291. 

  291.       ENDDO
    292. 

  292. !
    293. 

  293.       DO J=JTS,JTE
    294. 

  294.       DO K=KTE,KTS,-1
    295. 

  295.         KFLIP=KTE+1-K
    296. 

  296.         DO I=ITS,ITE
    297. 

  297.           ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J)
    298. 

  298.           APEX=1./EXNER(I,K,J)
    299. 

  299.           APE(I,K,J)=APEX
    300. 

  300.           TX=T(I,K,J)
    301. 

  301.           THE(I,K,J)=(CWM(I,K,J)*(-ELOCP/TX)+1.)*TH(I,K,J)
    302. 

  302.         ENDDO
    303. 

  303.       ENDDO
    304. 

  304.       ENDDO
    305. 

  305.       EL_MYJ(its:ite,:,jts:jte) = 0.
    306. 

  306. !
    307. 

  307. !----------------------------------------------------------------------
    308. 

  308.       setup_integration:  DO J=JTS,JTE
    309. 

  309. !----------------------------------------------------------------------
    310. 

  310. !
    311. 

  311.         DO I=ITS,ITE
    312. 

  312. !
    313. 

  313. !***  LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
    314. 

  314. !
    315. 

  315.           LMH=KTE-LOWLYR(I,J)+1
    316. 

  316. !
    317. 

  317.           PTOP=PINT(I,KTE+1,J)      ! KTE+1=KME
    318. 

  318.           PSFC=PINT(I,LOWLYR(I,J),J)
    319. 

  319. !
    320. 

  320. !***  CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
    321. 

  321. !
    322. 

  322.           SEAMASK=XLAND(I,J)-1.
    323. 

  323. !
    324. 

  324. !***  FILL 1-D VERTICAL ARRAYS
    325. 

  325. !***  AND FLIP DIRECTION SINCE MYJ SCHEME
    326. 

  326. !***  COUNTS DOWNWARD FROM THE DOMAIN'S TOP
    327. 

  327. !
    328. 

  328.           DO K=KTE,KTS,-1
    329. 

  329.             KFLIP=KTE+1-K
    330. 

  330.             TK(K)=T(I,KFLIP,J)
    331. 

  331.             THEK(K)=THE(I,KFLIP,J)
    332. 

  332.             RATIOMX=QV(I,KFLIP,J)
    333. 

  333.             QK(K)=RATIOMX/(1.+RATIOMX)
    334. 

  334.             CWMK(K)=CWM(I,KFLIP,J)
    335. 

  335.             PK(K)=PMID(I,KFLIP,J)
    336. 

  336.             UK(K)=U(I,KFLIP,J)
    337. 

  337.             VK(K)=V(I,KFLIP,J)
    338. 

  338. 

  339.             DZ_KFLIP(K)=DZ(I,KFLIP,J)
    340. 

  340.  
    341. 

  341. 

  342. 

  343. !
    344. 

  344. !***  TKE=0.5*(q**2) ==> q**2=2.*TKE
    345. 

  345. !
    346. 

  346.             Q2K(K)=2.*TKE_MYJ(I,KFLIP,J)
    347. 

  347. !
    348. 

  348. !***  COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
    349. 

  349. !
    350. 

  350.             ZHK(K)=ZINT(I,K,J)
    351. 

  351.           ENDDO
    352. 

  352. 

  353.           ZHK(KTE+1)=HT(I,J)          ! Z at bottom of lowest sigma layer
    354. 

  354. ! Multi-layer UCM
    355. 

  355.           if(flag_bep)then
    356. 

  356.            DO K=KTE,KTS,-1
    357. 

  357.             KFLIP=KTE+1-K
    358. 

  358.             B_E(K)=B_E_BEP(I,KFLIP,J)
    359. 

  359.             A_E(K)=A_E_BEP(I,KFLIP,J)
    360. 

  360.             SFK(K)=SF_BEP(I,KFLIP,J)
    361. 

  361.             VLK(K)=VL_BEP(I,KFLIP,J)
    362. 

  362.             B_E_FACE(K)=0.
    363. 

  363.             DL_U(K)=DL_U_BEP(I,KFLIP,J)
    364. 

  364.            enddo
    365. 

  365.           else
    366. 

  366.            DO K=KTE,KTS,-1
    367. 

  367.             B_E(K)=0.
    368. 

  368.             A_E(K)=0.
    369. 

  369.             SFK(K)=1.
    370. 

  370.             VLK(K)=1.
    371. 

  371.             B_E_FACE(K)=0.
    372. 

  372.             DL_U(K)=0.
    373. 

  373.            enddo
    374. 

  374.           endif
    375. 

  375. !UCM end
    376. 

  376. ! Multi-layer UCM - caculate B_E at grid faces
    377. 

  377. !          SFK(1)=1
    378. 

  378.          IF (FLAG_BEP) THEN
    379. 

  379.           DO K=KTE,KTS,-1
    380. 

  380.              IF (K==KTE) THEN
    381. 

  381.                 B_E_FACE(K)=0.0
    382. 

  382.              ELSEIF ( K == KTS ) THEN
    383. 

  383.                 B_E_FACE(K) = 0.0
    384. 

  384.              ELSE
    385. 

  385.                 SLOPE=(B_E(K-1)-B_E(K))/0.5/(DZ_KFLIP(K-1)+DZ_KFLIP(K))
    386. 

  386.                 ICEPT=B_E(K)-SLOPE*(ZHK(K)-0.5*DZ_KFLIP(K))
    387. 

  387.                 B_E_FACE(K)=2.*(SLOPE*ZHK(K)+ICEPT)                         !*2 accounts for Q2K=2.*TKE_MYJ
    388. 

  388.              ENDIF
    389. 

  389.           ENDDO
    390. 

  390.          ENDIF
    391. 

  391. !
    392. 

  392. !***  Begin debugging
    393. 

  393. !         IF(I==IMD.AND.J==JMD)THEN
    394. 

  394. !           PRINT_DIAG=1
    395. 

  395. !         ELSE
    396. 

  396. !           PRINT_DIAG=0
    397. 

  397. !         ENDIF
    398. 

  398. !         IF(I==227.AND.J==363)PRINT_DIAG=2
    399. 

  399. !***  End debugging
    400. 

  400. !
    401. 

  401. !----------------------------------------------------------------------
    402. 

  402. !***
    403. 

  403. !***  FIND THE MIXING LENGTH
    404. 

  404. !***
    405. 

  405.           CALL MIXLEN(LMH,UK,VK,TK,THEK,QK,CWMK                    &
    406. 

  406.      &               ,Q2K,ZHK,GM,GH,EL                                 &
    407. 

  407.      &               ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J)                 &
    408. 

  408. ! Multi-layer UCM
    409. 

  409.                      ,DL_U                                             &
    410. 

  410. ! UCM end
    411. 

  411.      &               ,IDS,IDE,JDS,JDE,KDS,KDE                          &
    412. 

  412.      &               ,IMS,IME,JMS,JME,KMS,KME                          &
    413. 

  413.      &               ,ITS,ITE,JTS,JTE,KTS,KTE)
    414. 

  414. !
    415. 

  415. !----------------------------------------------------------------------
    416. 

  416. !***
    417. 

  417. !***  SOLVE FOR THE PRODUCTION/DISSIPATION OF
    418. 

  418. !***  THE TURBULENT KINETIC ENERGY
    419. 

  419. !***
    420. 

  420. !
    421. 

  421.           CALL PRODQ2(LMH,DTTURBL,USTAR(I,J),GM,GH,EL,Q2K              &
    422. 

  422.      &               ,IDS,IDE,JDS,JDE,KDS,KDE                          &
    423. 

  423.      &               ,IMS,IME,JMS,JME,KMS,KME                          &
    424. 

  424.      &               ,ITS,ITE,JTS,JTE,KTS,KTE)
    425. 

  425. !
    426. 

  426. ! Multi-layer UCM
    427. 

  427.          IF (FLAG_BEP) THEN
    428. 

  428.             Q2K(LMH)=Q2K(LMH)*(1-FRC_URB2D(I,J))+(B_E_FACE(LMH)*       &
    429. 

  429.                  &         0.5*EL(LMH-1)/11.788/2.)**(2./3.)
    430. 

  430.          ENDIF
    431. 

  431. !UCM end
    432. 

  432. !
    433. 

  433. !----------------------------------------------------------------------
    434. 

  434. !*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
    435. 

  435. !----------------------------------------------------------------------
    436. 

  436. !
    437. 

  437.           KPBL(I,J)=KTE-LPBL(I,J)+1
    438. 

  438. !
    439. 

  439. !----------------------------------------------------------------------
    440. 

  440. !***
    441. 

  441. !***  FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
    442. 

  442. !***
    443. 

  443.           CALL DIFCOF(LMH,LMXL,GM,GH,EL,TK,Q2K,ZHK,AKMK,AKHK      &
    444. 

  444.      &               ,IDS,IDE,JDS,JDE,KDS,KDE                          &
    445. 

  445.      &               ,IMS,IME,JMS,JME,KMS,KME                          &
    446. 

  446.      &               ,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG)   ! debug
    447. 

  447. !
    448. 

  448. !***  COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH 
    449. 

  449. !***  ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1.  COUNTING 
    450. 

  450. !***  COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
    451. 

  451. !***  ARE DEFINED ON THE TOPS OF THE LAYERS KTS TO KTE-1.
    452. 

  452. !
    453. 

  453.           DO K=KTS,KTE-1
    454. 

  454.             KFLIP=KTE-K
    455. 

  455.             AKH(I,K,J)=AKHK(K)
    456. 

  456.             AKM(I,K,J)=AKMK(K)
    457. 

  457.             DELTAZ=0.5*(ZHK(KFLIP)-ZHK(KFLIP+2))
    458. 

  458. ! Multi-layer UCM
    459. 

  459.             EXCH_H(I,K+1,J)=AKHK(KFLIP)*DELTAZ
    460. 

  460.             EXCH_M(I,K+1,J)=AKMK(KFLIP)*DELTAZ
    461. 

  461. ! UCM end
    462. 

  462.           ENDDO
    463. 

  463. !----------------------------------------------------------------------
    464. 

  464. !***
    465. 

  465. !***  CARRY OUT THE VERTICAL DIFFUSION OF
    466. 

  466. !***  TURBULENT KINETIC ENERGY
    467. 

  467. !***
    468. 

  468. !
    469. 

  469. 

  470.           CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK                             &
    471. 

  471. ! Multi-layer UCM
    472. 

  472.      &              ,A_E,B_E_FACE,SFK,VLK                              &
    473. 

  473. ! UCM end
    474. 

  474.      &              ,IDS,IDE,JDS,JDE,KDS,KDE                           &
    475. 

  475.      &              ,IMS,IME,JMS,JME,KMS,KME                           &
    476. 

  476.      &              ,ITS,ITE,JTS,JTE,KTS,KTE)
    477. 

  477. !
    478. 

  478. !***  SAVE THE NEW TKE AND MIXING LENGTH.
    479. 

  479. !
    480. 

  480.           DO K=KTS,KTE
    481. 

  481.             KFLIP=KTE+1-K
    482. 

  482.             Q2K(KFLIP)=AMAX1(Q2K(KFLIP),EPSQ2)
    483. 

  483.             TKE_MYJ(I,K,J)=0.5*Q2K(KFLIP)
    484. 

  484.             IF(KFLIP/=KTE)EL_MYJ(I,K,J)=EL(KFLIP)   ! EL IS NOT DEFINED AT KTE (ground surface)
    485. 

  485.           ENDDO
    486. 

  486. !
    487. 

  487.         ENDDO
    488. 

  488. !
    489. 

  489. !----------------------------------------------------------------------
    490. 

  490.       ENDDO setup_integration
    491. 

  491. !----------------------------------------------------------------------
    492. 

  492. ! Multi-layer UCM - if idiff=1 then integration in PBL driver
    493. 

  493.       IF(IDIFF.NE.1)then 
    494. 

  494. ! UCM end
    495. 

  495. !
    496. 

  496. !***  CONVERT SURFACE SENSIBLE TEMPERATURE TO POTENTIAL TEMPERATURE.
    497. 

  497. !
    498. 

  498.       DO J=JTS,JTE
    499. 

  499.       DO I=ITS,ITE
    500. 

  500.         PSFC=PINT(I,LOWLYR(I,J),J)
    501. 

  501.         THSK(I,J)=TSK(I,J)*(1.E5/PSFC)**CAPA
    502. 

  502.       ENDDO
    503. 

  503.       ENDDO
    504. 

  504. !
    505. 

  505. !----------------------------------------------------------------------
    506. 

  506. !
    507. 

  507. !----------------------------------------------------------------------
    508. 

  508.       main_integration:  DO J=JTS,JTE
    509. 

  509. !----------------------------------------------------------------------
    510. 

  510. !
    511. 

  511.         DO I=ITS,ITE
    512. 

  512. !
    513. 

  513. !***  FILL 1-D VERTICAL ARRAYS
    514. 

  514. !***  AND FLIP DIRECTION SINCE MYJ SCHEME
    515. 

  515. !***  COUNTS DOWNWARD FROM THE DOMAIN'S TOP
    516. 

  516. !
    517. 

  517.           DO K=KTE,KTS,-1
    518. 

  518.             KFLIP=KTE+1-K
    519. 

  519.             THEK(K)=THE(I,KFLIP,J)
    520. 

  520.             RATIOMX=QV(I,KFLIP,J)                            !mixing ratio
    521. 

  521.             QK(K)=RATIOMX/(1.+RATIOMX)                       !specific humidity
    522. 

  522.             CWMK(K)=CWM(I,KFLIP,J)
    523. 

  523.             ZHK(K)=ZINT(I,K,J)
    524. 

  524.             RHOK(K,I)=PMID(I,KFLIP,J)/(R_D*T(I,KFLIP,J)*               &
    525. 

  525.      &                                (1.+P608*QK(K)-CWMK(K)))
    526. 

  526. 

  527.           ENDDO
    528. 

  528. !            SFK(1)=1.
    529. 

  529.           if(flag_bep)then
    530. 

  530.            ! Multi-layer UCM
    531. 

  531.            
    532. 

  532.             DO K=KTE,KTS,-1
    533. 

  533.              KFLIP=KTE+1-K
    534. 

  534.              A_Q(K)=A_Q_BEP(I,KFLIP,J)
    535. 

  535.              B_Q(K)=B_Q_BEP(I,KFLIP,J)
    536. 

  536.              A_T(K)=A_T_BEP(I,KFLIP,J)
    537. 

  537.              B_T(K)=B_T_BEP(I,KFLIP,J)
    538. 

  538.              SFK(K)=SF_BEP(I,KFLIP,J)
    539. 

  539.              VLK(K)=VL_BEP(I,KFLIP,J)             
    540. 

  540.             ENDDO
    541. 

  541. ! UCM end
    542. 

  542.           else
    543. 

  543.            DO K=KTE,KTS,-1
    544. 

  544.              A_Q(K)=0.
    545. 

  545.              B_Q(K)=0.
    546. 

  546.              A_T(K)=0.
    547. 

  547.              B_T(K)=0.
    548. 

  548.              SFK(K)=1.
    549. 

  549.              VLK(K)=1.
    550. 

  550.             ENDDO
    551. 

  551.           endif           
    552. 

  552. !
    553. 

  553. !***  COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
    554. 

  554. !***  ARE DEFINED ON THE BOTTOMS OF THE LAYERS KTS TO KTE-1.  THESE COEFFICIENTS
    555. 

  555. !***  ARE ALSO MULTIPLIED BY THE DENSITY AT THE BOTTOM INTERFACE LEVEL.
    556. 

  556. !
    557. 

  557.           DO K=KTS,KTE-1
    558. 

  558.             AKHK(K)=AKH(I,K,J)*0.5*(RHOK(K,I)+RHOK(K+1,I))
    559. 

  559.           ENDDO
    560. 

  560. !
    561. 

  561.           ZHK(KTE+1)=ZINT(I,KTE+1,J)
    562. 

  562. !
    563. 

  563.           SEAMASK=XLAND(I,J)-1.
    564. 

  564.           THZ0(I,J)=(1.-SEAMASK)*THSK(I,J)+SEAMASK*THZ0(I,J)
    565. 

  565. !
    566. 

  566.           LLOW=LOWLYR(I,J)
    567. 

  567.           AKHS_DENS=AKHS(I,J)*RHOK(KTE+1-LLOW,I)
    568. 

  568. !
    569. 

  569.           IF(SEAMASK<0.5)THEN
    570. 

  570.             QFC1=XLV*CHKLOWQ(I,J)*AKHS_DENS
    571. 

  571. !
    572. 

  572.             IF(SNOW(I,J)>0..OR.SICE(I,J)>0.5)THEN
    573. 

  573.               QFC1=QFC1*RLIVWV
    574. 

  574.             ENDIF
    575. 

  575. !
    576. 

  576.             IF(QFC1>0.)THEN
    577. 

  577.               QLOW=QK(KTE+1-LLOW)
    578. 

  578.               QSFC(I,J)=QLOW+ELFLX(I,J)/QFC1
    579. 

  579.             ENDIF
    580. 

  580. !
    581. 

  581.           ELSE
    582. 

  582.             PSFC=PINT(I,LOWLYR(I,J),J)
    583. 

  583.             EXNSFC=(1.E5/PSFC)**CAPA
    584. 

  584.             QSFC(I,J)=PQ0SEA/PSFC                                      &
    585. 

  585.      &         *EXP(A2*(THSK(I,J)-A3*EXNSFC)/(THSK(I,J)-A4*EXNSFC))
    586. 

  586.           ENDIF
    587. 

  587. !
    588. 

  588.           QZ0 (I,J)=(1.-SEAMASK)*QSFC(I,J)+SEAMASK*QZ0 (I,J)
    589. 

  589. !
    590. 

  590. !***  LOWEST LAYER ABOVE GROUND MUST BE FLIPPED
    591. 

  591. !
    592. 

  592.           LMH=KTE-LOWLYR(I,J)+1
    593. 

  593. !
    594. 

  594. !----------------------------------------------------------------------
    595. 

  595. !***  CARRY OUT THE VERTICAL DIFFUSION OF
    596. 

  596. !***  TEMPERATURE AND WATER VAPOR
    597. 

  597. !----------------------------------------------------------------------
    598. 

  598. !
    599. 

  599.          
    600. 

  600.           CALL VDIFH(DTDIF,LMH,THZ0(I,J),QZ0(I,J)                      &
    601. 

  601.      &              ,AKHS_DENS,CHKLOWQ(I,J),CT(I,J)                    &
    602. 

  602.      &              ,THEK,QK,CWMK,AKHK,ZHK,RHOK(KTS,I)                 &
    603. 

  603. ! Multi-layer UCM
    604. 

  604.      &              ,A_T,B_T,A_Q,B_Q,SFK,VLK                           &
    605. 

  605. ! UCM end
    606. 

  606.      &              ,IDS,IDE,JDS,JDE,KDS,KDE                           &
    607. 

  607.      &              ,IMS,IME,JMS,JME,KMS,KME                           &
    608. 

  608.      &              ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
    609. 

  609. !----------------------------------------------------------------------
    610. 

  610. !***
    611. 

  611. !***  COMPUTE PRIMARY VARIABLE TENDENCIES
    612. 

  612. !***
    613. 

  613.           DO K=KTS,KTE
    614. 

  614.             KFLIP=KTE+1-K
    615. 

  615.             THOLD=TH(I,K,J)
    616. 

  616.             THNEW=THEK(KFLIP)+CWMK(KFLIP)*ELOCP*APE(I,K,J)
    617. 

  617.             DTDT=(THNEW-THOLD)*RDTTURBL
    618. 

  618.             QOLD=QV(I,K,J)/(1.+QV(I,K,J))
    619. 

  619.             DQDT=(QK(KFLIP)-QOLD)*RDTTURBL
    620. 

  620.             DCDT=(CWMK(KFLIP)-CWM(I,K,J))*RDTTURBL
    621. 

  621. !
    622. 

  622.             RTHBLTEN(I,K,J)=DTDT
    623. 

  623.             RQVBLTEN(I,K,J)=DQDT/(1.-QK(KFLIP))**2
    624. 

  624.             RQCBLTEN(I,K,J)=DCDT
    625. 

  625.           ENDDO
    626. 

  626. !
    627. 

  627.           
    628. 

  628. !*** Begin debugging
    629. 

  629. !         IF(I==IMD.AND.J==JMD)THEN
    630. 

  630. !           PRINT_DIAG=0
    631. 

  631. !         ELSE
    632. 

  632. !           PRINT_DIAG=0
    633. 

  633. !         ENDIF
    634. 

  634. !         IF(I==227.AND.J==363)PRINT_DIAG=0
    635. 

  635. !*** End debugging
    636. 

  636. !
    637. 

  637.         PSFC=.01*PINT(I,LOWLYR(I,J),J)
    638. 

  638.         ZSL_DIAG=0.5*DZ(I,1,J)
    639. 

  639. !
    640. 

  640. !*** Begin debugging
    641. 

  641. !         IF(PRINT_DIAG==1)THEN
    642. 

  642. !
    643. 

  643. !           write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
    644. 

  644. !           '{turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
    645. 

  645. !           , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag  &
    646. 

  646. !           , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
    647. 

  647. !           write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
    648. 

  648. !           '{turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
    649. 

  649. !           , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
    650. 

  650. !           , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
    651. 

  651. !           write(6,"(a)") &
    652. 

  652. !           '{turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
    653. 

  653. !           do k=kts,kte/2
    654. 

  654. !             KFLIP=KTE-K   !-- Includes the KFLIP-1 in earlier versions
    655. 

  655. !             write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
    656. 

  656. !            '{turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
    657. 

  657. !            , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
    658. 

  658. !            , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
    659. 

  659. !            , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
    660. 

  660. !           enddo
    661. 

  661. !
    662. 

  662. !         ELSEIF(PRINT_DIAG==2)THEN
    663. 

  663. !
    664. 

  664. !           write(6,"(a, 2i5, 2i3, 2f8.2, f6.2, 2f8.2)") &
    665. 

  665. !           '}turb4 i,j, Kpbl, Kmxl, Psfc, Zsfc, Zsl, Zpbl, Zmxl = ' &
    666. 

  666. !           , i, j, KPBL(i,j), KTE-LMXL+1, PSFC, ZHK(LMH+1), ZSL_diag  &
    667. 

  667. !           , PBLH(i,j), ZHK(LMXL)-ZHK(LMH+1)
    668. 

  668. !           write(6,"(a, 2f7.2, f7.3, 3e11.4)") &
    669. 

  669. !           '}turb4 tsk, thsk, qz0, q**2_0, akhs, exch_0 = ' &
    670. 

  670. !           , tsk(i,j)-273.15, thsk(i,j), 1000.*qz0(i,j) &
    671. 

  671. !           , 2.*tke_myj(i,1,j), akhs(i,j), akhs(i,j)*ZSL_diag
    672. 

  672. !           write(6,"(a)") &
    673. 

  673. !           '}turb5 k, Pmid, Pint_1, Tc, TH, DTH, GH, GM, EL, Q**2, Akh, EXCH_h, Dz, Dp'
    674. 

  674. !           do k=kts,kte/2
    675. 

  675. !             KFLIP=KTE-K   !-- Includes the KFLIP-1 in earlier versions
    676. 

  676. !             write(6,"(a,i3, 2f8.2, 2f8.3, 3e12.4, 4e11.4, f7.2, f6.2)") &
    677. 

  677. !            '}turb5 ', k, .01*pmid(i,k,j),.01*pint(i,k,j), T(i,k,j)-273.15 &
    678. 

  678. !            , th(i,k,j), DTTURBL*rthblten(i,k,j), GH(KFLIP), GM(KFLIP) &
    679. 

  679. !            , el_myj(i,KFLIP,j), 2.*tke_myj(i,k+1,j), akh(i,KFLIP,j) &
    680. 

  680. !            , exch_h(i,k,j), dz(i,k,j), .01*(pint(i,k,j)-pint(i,k+1,j))
    681. 

  681. !           enddo
    682. 

  682. !         ENDIF
    683. 

  683. !*** End debugging
    684. 

  684. !
    685. 

  685. !----------------------------------------------------------------------
    686. 

  686.         ENDDO
    687. 

  687. !----------------------------------------------------------------------
    688. 

  688.         DO I=ITS,ITE
    689. 

  689. !
    690. 

  690. !***  FILL 1-D VERTICAL ARRAYS
    691. 

  691. !***  AND FLIP DIRECTION SINCE MYJ SCHEME
    692. 

  692. !***  COUNTS DOWNWARD FROM THE DOMAIN'S TOP
    693. 

  693. !
    694. 

  694. 

  695.           DO K=KTS,KTE-1
    696. 

  696.             AKMK(K)=AKM(I,K,J)
    697. 

  697.             AKMK(K)=AKMK(K)*(RHOK(K,I)+RHOK(K+1,I))*0.5
    698. 

  698.           ENDDO
    699. 

  699. !
    700. 

  700.           LLOW=LOWLYR(I,J)
    701. 

  701.           AKMS_DENS=AKMS(I,J)*RHOK(KTE+1-LLOW,I)
    702. 

  702. !
    703. 

  703. 

  704.           DO K=KTE,KTS,-1
    705. 

  705.             KFLIP=KTE+1-K
    706. 

  706.             UK(K)=U(I,KFLIP,J)
    707. 

  707.             VK(K)=V(I,KFLIP,J)
    708. 

  708.             ZHK(K)=ZINT(I,K,J)
    709. 

  709. 

  710.           ENDDO            
    711. 

  711. !         SFK(1)=1.
    712. 

  712.           if(flag_bep)then
    713. 

  713. ! Multi-layer UCM           
    714. 

  714.            DO K=KTE,KTS,-1
    715. 

  715.             KFLIP=KTE+1-K            
    716. 

  716.             A_U(K)=A_U_BEP(I,KFLIP,J)
    717. 

  717.             A_V(K)=A_V_BEP(I,KFLIP,J)
    718. 

  718.             B_U(K)=B_U_BEP(I,KFLIP,J)
    719. 

  719.             B_V(K)=B_V_BEP(I,KFLIP,J)
    720. 

  720.             SFK(K)=SF_BEP(I,KFLIP,J)
    721. 

  721.             VLK(K)=VL_BEP(I,KFLIP,J)
    722. 

  722.            ENDDO
    723. 

  723. ! UCM end
    724. 

  724.          else
    725. 

  725.            DO K=KTE,KTS,-1            
    726. 

  726.             A_U(K)=0.
    727. 

  727.             A_V(K)=0.
    728. 

  728.             B_U(K)=0.
    729. 

  729.             B_V(K)=0.
    730. 

  730.             SFK(K)=1.
    731. 

  731.             VLK(K)=1.
    732. 

  732.            ENDDO
    733. 

  733.           endif
    734. 

  734. 

  735.           ZHK(KTE+1)=ZINT(I,KTE+1,J)
    736. 

  736. !
    737. 

  737. !----------------------------------------------------------------------
    738. 

  738. !***  CARRY OUT THE VERTICAL DIFFUSION OF
    739. 

  739. !***  VELOCITY COMPONENTS
    740. 

  740. !----------------------------------------------------------------------
    741. 

  741. !
    742. 

  742.           CALL VDIFV(LMH,DTDIF,UZ0(I,J),VZ0(I,J)                   &
    743. 

  743.      &              ,AKMS_DENS,UK,VK,AKMK,ZHK,RHOK(KTS,I)              &
    744. 

  744. ! Multi-layer UCM
    745. 

  745.      &              ,A_U,A_V,B_U,B_V,SFK,VLK                           &
    746. 

  746. ! UCM end
    747. 

  747.      &              ,IDS,IDE,JDS,JDE,KDS,KDE                           &
    748. 

  748.      &              ,IMS,IME,JMS,JME,KMS,KME                           &
    749. 

  749.      &              ,ITS,ITE,JTS,JTE,KTS,KTE,I,J)
    750. 

  750. !
    751. 

  751. !----------------------------------------------------------------------
    752. 

  752. !***
    753. 

  753. !***  COMPUTE PRIMARY VARIABLE TENDENCIES
    754. 

  754. !***
    755. 

  755.           DO K=KTS,KTE
    756. 

  756.             KFLIP=KTE+1-K
    757. 

  757.             DUDT=(UK(KFLIP)-U(I,K,J))*RDTTURBL
    758. 

  758.             DVDT=(VK(KFLIP)-V(I,K,J))*RDTTURBL
    759. 

  759.             RUBLTEN(I,K,J)=DUDT
    760. 

  760.             RVBLTEN(I,K,J)=DVDT
    761. 

  761.           ENDDO
    762. 

  762. !
    763. 

  763.         ENDDO
    764. 

  764. !----------------------------------------------------------------------
    765. 

  765. !
    766. 

  766.       ENDDO main_integration
    767. 

  767. !
    768. 

  768.       ENDIF
    769. 

  769. !----------------------------------------------------------------------
    770. 

  770. !
    771. 

  771.       END SUBROUTINE MYJURB
    772. 

  772. !
    773. 

  773. !----------------------------------------------------------------------
    774. 

  774. !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    775. 

  775. !----------------------------------------------------------------------
    776. 

  776.                           SUBROUTINE MIXLEN                           &
    777. 

  777. !----------------------------------------------------------------------
    778. 

  778. !   ******************************************************************
    779. 

  779. !   *                                                                *
    780. 

  780. !   *                   LEVEL 2.5 MIXING LENGTH                      *
    781. 

  781. !   *                                                                *
    782. 

  782. !   ******************************************************************
    783. 

  783. !
    784. 

  784.      &(LMH,U,V,T,THE,Q,CWM,Q2,Z,GM,GH,EL,PBLH,LPBL,LMXL,CT             &
    785. 

  785. ! Multi-layer UCM
    786. 

  786.      &,DL_U                                                            &
    787. 

  787. ! UCM end
    788. 

  788.      &,IDS,IDE,JDS,JDE,KDS,KDE                                         &
    789. 

  789.      &,IMS,IME,JMS,JME,KMS,KME                                         &
    790. 

  790.      &,ITS,ITE,JTS,JTE,KTS,KTE)
    791. 

  791. !----------------------------------------------------------------------
    792. 

  792. !
    793. 

  793.       IMPLICIT NONE
    794. 

  794. !
    795. 

  795. !----------------------------------------------------------------------
    796. 

  796.       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
    797. 

  797.      &                     ,IMS,IME,JMS,JME,KMS,KME                    &
    798. 

  798.      &                     ,ITS,ITE,JTS,JTE,KTS,KTE
    799. 

  799. !
    800. 

  800.       INTEGER,INTENT(IN) :: LMH
    801. 

  801. !
    802. 

  802.       INTEGER,INTENT(OUT) :: LMXL,LPBL
    803. 

  803. !
    804. 

  804.       REAL,DIMENSION(KTS:KTE),INTENT(IN) :: CWM,Q,Q2,T,THE,U,V
    805. 

  805. ! Multi-layer UCM
    806. 

  806.       REAL,DIMENSION(KTS:KTE),INTENT(IN) :: DL_U
    807. 

  807. ! UCM end
    808. 

  808. !
    809. 

  809.       REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
    810. 

  810. !
    811. 

  811.       REAL,INTENT(OUT) :: PBLH
    812. 

  812. !
    813. 

  813.       REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: EL,GH,GM
    814. 

  814. !
    815. 

  815.       REAL,INTENT(INOUT) :: CT
    816. 

  816. !----------------------------------------------------------------------
    817. 

  817. !***
    818. 

  818. !***  LOCAL VARIABLES
    819. 

  819. !***
    820. 

  820.       INTEGER :: K,LPBLM
    821. 

  821. !
    822. 

  822.       REAL :: A,ADEN,B,BDEN,AUBR,BUBR,BLMX,EL0,ELOQ2X,GHL,GML           &
    823. 

  823.      &       ,QOL2ST,QOL2UN,QDZL,RDZ,SQ,SREL,SZQ,TEM,THM,VKRMZ
    824. 

  824. !
    825. 

  825.       REAL,DIMENSION(KTS:KTE) :: Q1
    826. 

  826. !
    827. 

  827.       REAL,DIMENSION(KTS:KTE-1) :: DTH,ELM,REL
    828. 

  828. !
    829. 

  829. !----------------------------------------------------------------------
    830. 

  830. !**********************************************************************
    831. 

  831. !--------------FIND THE HEIGHT OF THE PBL-------------------------------
    832. 

  832.       LPBL=LMH
    833. 

  833. !
    834. 

  834.       DO K=LMH-1,1,-1
    835. 

  835.         IF(Q2(K)<=EPSQ2*FH)THEN
    836. 

  836.           LPBL=K
    837. 

  837.           GO TO 110
    838. 

  838.         ENDIF
    839. 

  839.       ENDDO
    840. 

  840. !
    841. 

  841.       LPBL=1
    842. 

  842. !
    843. 

  843. !--------------THE HEIGHT OF THE PBL------------------------------------
    844. 

  844. !
    845. 

  845.  110  PBLH=Z(LPBL+1)-Z(LMH+1)
    846. 

  846. !
    847. 

  847. !-----------------------------------------------------------------------
    848. 

  848.       DO K=KTS,LMH
    849. 

  849.         Q1(K)=0.
    850. 

  850.       ENDDO
    851. 

  851. !
    852. 

  852.       DO K=1,LMH-1
    853. 

  853.         DTH(K)=THE(K)-THE(K+1)
    854. 

  854.       ENDDO
    855. 

  855. !
    856. 

  856.       DO K=LMH-2,1,-1
    857. 

  857.         IF(DTH(K)>0..AND.DTH(K+1)<=0.)THEN
    858. 

  858.           DTH(K)=DTH(K)+CT
    859. 

  859.           EXIT
    860. 

  860.         ENDIF
    861. 

  861.       ENDDO
    862. 

  862. !
    863. 

  863.       CT=0.
    864. 

  864. !----------------------------------------------------------------------
    865. 

  865.       DO K=KTS,LMH-1
    866. 

  866.         RDZ=2./(Z(K)-Z(K+2))
    867. 

  867.         GML=((U(K)-U(K+1))**2+(V(K)-V(K+1))**2)*RDZ*RDZ
    868. 

  868.         GM(K)=MAX(GML,EPSGM)
    869. 

  869. !
    870. 

  870.         TEM=(T(K)+T(K+1))*0.5
    871. 

  871.         THM=(THE(K)+THE(K+1))*0.5
    872. 

  872. !
    873. 

  873.         A=THM*P608
    874. 

  874.         B=(ELOCP/TEM-1.-P608)*THM
    875. 

  875. !
    876. 

  876.         GHL=(DTH(K)*((Q(K)+Q(K+1)+CWM(K)+CWM(K+1))*(0.5*P608)+1.)      &
    877. 

  877.      &     +(Q(K)-Q(K+1)+CWM(K)-CWM(K+1))*A                            &
    878. 

  878.      &     +(CWM(K)-CWM(K+1))*B)*RDZ
    879. 

  879. !
    880. 

  880.         IF(ABS(GHL)<=EPSGH)GHL=EPSGH
    881. 

  881.         GH(K)=GHL
    882. 

  882.       ENDDO
    883. 

  883. !
    884. 

  884. !----------------------------------------------------------------------
    885. 

  885. !***  FIND MAXIMUM MIXING LENGTHS AND THE LEVEL OF THE PBL TOP
    886. 

  886. !----------------------------------------------------------------------
    887. 

  887. !
    888. 

  888.       LMXL=LMH
    889. 

  889. !
    890. 

  890.       DO K=KTS,LMH-1
    891. 

  891.         GML=GM(K)
    892. 

  892.         GHL=GH(K)
    893. 

  893. !
    894. 

  894.         IF(GHL>=EPSGH)THEN
    895. 

  895.           IF(GML/GHL<=REQU)THEN
    896. 

  896.             ELM(K)=EPSL
    897. 

  897.             LMXL=K
    898. 

  898.           ELSE
    899. 

  899.             AUBR=(AUBM*GML+AUBH*GHL)*GHL
    900. 

  900.             BUBR= BUBM*GML+BUBH*GHL
    901. 

  901.             QOL2ST=(-0.5*BUBR+SQRT(BUBR*BUBR*0.25-AUBR*CUBR))*RCUBR
    902. 

  902.             ELOQ2X=1./QOL2ST
    903. 

  903.             ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
    904. 

  904.           ENDIF
    905. 

  905.         ELSE
    906. 

  906.           ADEN=(ADNM*GML+ADNH*GHL)*GHL
    907. 

  907.           BDEN= BDNM*GML+BDNH*GHL
    908. 

  908.           QOL2UN=-0.5*BDEN+SQRT(BDEN*BDEN*0.25-ADEN)
    909. 

  909.           ELOQ2X=1./(QOL2UN+EPSRU)       ! repsr1/qol2un
    910. 

  910.           ELM(K)=MAX(SQRT(ELOQ2X*Q2(K)),EPSL)
    911. 

  911.         ENDIF
    912. 

  912.       ENDDO
    913. 

  913. !
    914. 

  914.       IF(ELM(LMH-1)==EPSL)LMXL=LMH
    915. 

  915. !
    916. 

  916. !----------------------------------------------------------------------
    917. 

  917. !***  THE HEIGHT OF THE MIXED LAYER
    918. 

  918. !----------------------------------------------------------------------
    919. 

  919. !
    920. 

  920.       BLMX=Z(LMXL)-Z(LMH+1)
    921. 

  921. !
    922. 

  922. !----------------------------------------------------------------------
    923. 

  923.       DO K=LPBL,LMH
    924. 

  924.         Q1(K)=SQRT(Q2(K))
    925. 

  925.       ENDDO
    926. 

  926. !----------------------------------------------------------------------
    927. 

  927.       SZQ=0.
    928. 

  928.       SQ =0.
    929. 

  929. !
    930. 

  930.       DO K=KTS,LMH-1
    931. 

  931.         QDZL=(Q1(K)+Q1(K+1))*(Z(K+1)-Z(K+2))
    932. 

  932.         SZQ=(Z(K+1)+Z(K+2)-Z(LMH+1)-Z(LMH+1))*QDZL+SZQ
    933. 

  933.         SQ=QDZL+SQ
    934. 

  934.       ENDDO
    935. 

  935. !
    936. 

  936. !----------------------------------------------------------------------
    937. 

  937. !***  COMPUTATION OF ASYMPTOTIC L IN BLACKADAR FORMULA
    938. 

  938. !----------------------------------------------------------------------
    939. 

  939. !
    940. 

  940.       EL0=MIN(ALPH*SZQ*0.5/SQ,EL0MAX)
    941. 

  941.       EL0=MAX(EL0            ,EL0MIN)
    942. 

  942. !
    943. 

  943. !----------------------------------------------------------------------
    944. 

  944. !***  ABOVE THE PBL TOP
    945. 

  945. !----------------------------------------------------------------------
    946. 

  946. !
    947. 

  947.       LPBLM=MAX(LPBL-1,1)
    948. 

  948. !
    949. 

  949.       DO K=KTS,LPBLM
    950. 

  950.         EL(K)=MIN((Z(K)-Z(K+2))*ELFC,ELM(K))
    951. 

  951.         REL(K)=EL(K)/ELM(K)
    952. 

  952.       ENDDO
    953. 

  953. !
    954. 

  954. !----------------------------------------------------------------------
    955. 

  955. !***  INSIDE THE PBL
    956. 

  956. !----------------------------------------------------------------------
    957. 

  957. !
    958. 

  958.       IF(LPBL<LMH)THEN
    959. 

  959.         DO K=LPBL,LMH-1
    960. 

  960.           VKRMZ=(Z(K+1)-Z(LMH+1))*VKARMAN
    961. 

  961. ! Multi-layer UCM
    962. 

  962.           IF (DL_U(K).GT.0.0) THEN
    963. 

  963.                EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))
    964. 

  964.                EL(K)=1./(1./EL(K)+1./DL_U(K))
    965. 

  965.             ELSE
    966. 

  966.                EL(K)=MIN(VKRMZ/(VKRMZ/EL0+1.),ELM(K))       
    967. 

  967.           ENDIF
    968. 

  968. ! UCM end
    969. 

  969.           REL(K)=EL(K)/ELM(K)
    970. 

  970.         ENDDO
    971. 

  971.       ENDIF
    972. 

  972. !
    973. 

  973.       DO K=LPBL+1,LMH-2
    974. 

  974.         SREL=MIN(((REL(K-1)+REL(K+1))*0.5+REL(K))*0.5,REL(K))
    975. 

  975.         EL(K)=MAX(SREL*ELM(K),EPSL)
    976. 

  976.       ENDDO
    977. 

  977. !
    978. 

  978. !----------------------------------------------------------------------
    979. 

  979.       END SUBROUTINE MIXLEN
    980. 

  980. !----------------------------------------------------------------------
    981. 

  981. !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    982. 

  982. !----------------------------------------------------------------------
    983. 

  983.                           SUBROUTINE PRODQ2                            &
    984. 

  984. !----------------------------------------------------------------------
    985. 

  985. !   ******************************************************************
    986. 

  986. !   *                                                                *
    987. 

  987. !   *            LEVEL 2.5 Q2 PRODUCTION/DISSIPATION                 *
    988. 

  988. !   *                                                                *
    989. 

  989. !   ******************************************************************
    990. 

  990. !
    991. 

  991.      &(LMH,DTTURBL,USTAR,GM,GH,EL,Q2                                   &
    992. 

  992.      &,IDS,IDE,JDS,JDE,KDS,KDE                                         &
    993. 

  993.      &,IMS,IME,JMS,JME,KMS,KME                                         &
    994. 

  994.      &,ITS,ITE,JTS,JTE,KTS,KTE)
    995. 

  995. !----------------------------------------------------------------------
    996. 

  996. !
    997. 

  997.       IMPLICIT NONE
    998. 

  998. !
    999. 

  999. !----------------------------------------------------------------------
    1000. 

  1000.       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
    1001. 

  1001.      &                     ,IMS,IME,JMS,JME,KMS,KME                    &
    1002. 

  1002.      &                     ,ITS,ITE,JTS,JTE,KTS,KTE
    1003. 

  1003. !
    1004. 

  1004.       INTEGER,INTENT(IN) :: LMH
    1005. 

  1005. !
    1006. 

  1006.       REAL,INTENT(IN) :: DTTURBL,USTAR
    1007. 

  1007. !
    1008. 

  1008.       REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: GH,GM
    1009. 

  1009.       REAL,DIMENSION(KTS:KTE-1),INTENT(INOUT) :: EL
    1010. 

  1010. !
    1011. 

  1011.       REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: Q2
    1012. 

  1012. !----------------------------------------------------------------------
    1013. 

  1013. !***
    1014. 

  1014. !***  LOCAL VARIABLES
    1015. 

  1015. !***
    1016. 

  1016.       INTEGER :: K
    1017. 

  1017. !
    1018. 

  1018.       REAL :: ADEN,AEQU,ANUM,ARHS,BDEN,BEQU,BNUM,BRHS,CDEN,CRHS        &
    1019. 

  1019.      &       ,DLOQ1,ELOQ11,ELOQ12,ELOQ13,ELOQ21,ELOQ22,ELOQ31,ELOQ32   &
    1020. 

  1020.      &       ,ELOQ41,ELOQ42,ELOQ51,ELOQ52,ELOQN,EQOL2,GHL,GML          &
    1021. 

  1021.      &       ,RDEN1,RDEN2,RHS2,RHSP1,RHSP2,RHST2
    1022. 

  1022. !
    1023. 

  1023. !----------------------------------------------------------------------
    1024. 

  1024. !**********************************************************************
    1025. 

  1025. !----------------------------------------------------------------------
    1026. 

  1026. !
    1027. 

  1027.       main_integration: DO K=1,LMH-1
    1028. 

  1028.         GML=GM(K)
    1029. 

  1029.         GHL=GH(K)
    1030. 

  1030. !
    1031. 

  1031. !----------------------------------------------------------------------
    1032. 

  1032. !***  COEFFICIENTS OF THE EQUILIBRIUM EQUATION
    1033. 

  1033. !----------------------------------------------------------------------
    1034. 

  1034. !
    1035. 

  1035.         AEQU=(AEQM*GML+AEQH*GHL)*GHL
    1036. 

  1036.         BEQU= BEQM*GML+BEQH*GHL
    1037. 

  1037. !
    1038. 

  1038. !----------------------------------------------------------------------
    1039. 

  1039. !***  EQUILIBRIUM SOLUTION FOR L/Q
    1040. 

  1040. !----------------------------------------------------------------------
    1041. 

  1041. !
    1042. 

  1042.         EQOL2=-0.5*BEQU+SQRT(BEQU*BEQU*0.25-AEQU)
    1043. 

  1043. !
    1044. 

  1044. !----------------------------------------------------------------------
    1045. 

  1045. !***  IS THERE PRODUCTION/DISSIPATION ?
    1046. 

  1046. !----------------------------------------------------------------------
    1047. 

  1047. !
    1048. 

  1048.         IF((GML+GHL*GHL<=EPSTRB)                                       &
    1049. 

  1049.      &   .OR.(GHL>=EPSGH.AND.GML/GHL<=REQU)                            &
    1050. 

  1050.      &   .OR.(EQOL2<=EPS2))THEN
    1051. 

  1051. !
    1052. 

  1052. !----------------------------------------------------------------------
    1053. 

  1053. !***  NO TURBULENCE
    1054. 

  1054. !----------------------------------------------------------------------
    1055. 

  1055. !
    1056. 

  1056.           Q2(K)=EPSQ2
    1057. 

  1057.           EL(K)=EPSL
    1058. 

  1058. !----------------------------------------------------------------------
    1059. 

  1059. !
    1060. 

  1060.         ELSE
    1061. 

  1061. !
    1062. 

  1062. !----------------------------------------------------------------------
    1063. 

  1063. !***  TURBULENCE
    1064. 

  1064. !----------------------------------------------------------------------
    1065. 

  1065. !----------------------------------------------------------------------
    1066. 

  1066. !***  COEFFICIENTS OF THE TERMS IN THE NUMERATOR
    1067. 

  1067. !----------------------------------------------------------------------
    1068. 

  1068. !
    1069. 

  1069.           ANUM=(ANMM*GML+ANMH*GHL)*GHL
    1070. 

  1070.           BNUM= BNMM*GML+BNMH*GHL
    1071. 

  1071. !
    1072. 

  1072. !----------------------------------------------------------------------
    1073. 

  1073. !***  COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
    1074. 

  1074. !----------------------------------------------------------------------
    1075. 

  1075. !
    1076. 

  1076.           ADEN=(ADNM*GML+ADNH*GHL)*GHL
    1077. 

  1077.           BDEN= BDNM*GML+BDNH*GHL
    1078. 

  1078.           CDEN= 1.
    1079. 

  1079. !
    1080. 

  1080. !----------------------------------------------------------------------
    1081. 

  1081. !***  COEFFICIENTS OF THE NUMERATOR OF THE LINEARIZED EQ.
    1082. 

  1082. !----------------------------------------------------------------------
    1083. 

  1083. !
    1084. 

  1084.           ARHS=-(ANUM*BDEN-BNUM*ADEN)*2.
    1085. 

  1085.           BRHS=- ANUM*4.
    1086. 

  1086.           CRHS=- BNUM*2.
    1087. 

  1087. !
    1088. 

  1088. !----------------------------------------------------------------------
    1089. 

  1089. !***  INITIAL VALUE OF L/Q
    1090. 

  1090. !----------------------------------------------------------------------
    1091. 

  1091. !
    1092. 

  1092.           DLOQ1=EL(K)/SQRT(Q2(K))
    1093. 

  1093. !
    1094. 

  1094. !----------------------------------------------------------------------
    1095. 

  1095. !***  FIRST ITERATION FOR L/Q, RHS=0
    1096. 

  1096. !----------------------------------------------------------------------
    1097. 

  1097. !
    1098. 

  1098.           ELOQ21=1./EQOL2
    1099. 

  1099.           ELOQ11=SQRT(ELOQ21)
    1100. 

  1100.           ELOQ31=ELOQ21*ELOQ11
    1101. 

  1101.           ELOQ41=ELOQ21*ELOQ21
    1102. 

  1102.           ELOQ51=ELOQ21*ELOQ31
    1103. 

  1103. !
    1104. 

  1104. !----------------------------------------------------------------------
    1105. 

  1105. !***  1./DENOMINATOR
    1106. 

  1106. !----------------------------------------------------------------------
    1107. 

  1107. !
    1108. 

  1108.           RDEN1=1./(ADEN*ELOQ41+BDEN*ELOQ21+CDEN)
    1109. 

  1109. !
    1110. 

  1110. !----------------------------------------------------------------------
    1111. 

  1111. !***  D(RHS)/D(L/Q)
    1112. 

  1112. !----------------------------------------------------------------------
    1113. 

  1113. !
    1114. 

  1114.           RHSP1=(ARHS*ELOQ51+BRHS*ELOQ31+CRHS*ELOQ11)*RDEN1*RDEN1
    1115. 

  1115. !
    1116. 

  1116. !----------------------------------------------------------------------
    1117. 

  1117. !***  FIRST-GUESS SOLUTION
    1118. 

  1118. !----------------------------------------------------------------------
    1119. 

  1119. !
    1120. 

  1120.           ELOQ12=ELOQ11+(DLOQ1-ELOQ11)*EXP(RHSP1*DTTURBL)
    1121. 

  1121.           ELOQ12=MAX(ELOQ12,EPS1)
    1122. 

  1122. !
    1123. 

  1123. !----------------------------------------------------------------------
    1124. 

  1124. !***  SECOND ITERATION FOR L/Q
    1125. 

  1125. !----------------------------------------------------------------------
    1126. 

  1126. !
    1127. 

  1127.           ELOQ22=ELOQ12*ELOQ12
    1128. 

  1128.           ELOQ32=ELOQ22*ELOQ12
    1129. 

  1129.           ELOQ42=ELOQ22*ELOQ22
    1130. 

  1130.           ELOQ52=ELOQ22*ELOQ32
    1131. 

  1131. !
    1132. 

  1132. !----------------------------------------------------------------------
    1133. 

  1133. !***  1./DENOMINATOR
    1134. 

  1134. !----------------------------------------------------------------------
    1135. 

  1135. !
    1136. 

  1136.           RDEN2=1./(ADEN*ELOQ42+BDEN*ELOQ22+CDEN)
    1137. 

  1137.           RHS2 =-(ANUM*ELOQ42+BNUM*ELOQ22)*RDEN2+RB1
    1138. 

  1138.           RHSP2= (ARHS*ELOQ52+BRHS*ELOQ32+CRHS*ELOQ12)*RDEN2*RDEN2
    1139. 

  1139.           RHST2=RHS2/RHSP2
    1140. 

  1140. !
    1141. 

  1141. !----------------------------------------------------------------------
    1142. 

  1142. !***  CORRECTED SOLUTION
    1143. 

  1143. !----------------------------------------------------------------------
    1144. 

  1144. !
    1145. 

  1145.           ELOQ13=ELOQ12-RHST2+(RHST2+DLOQ1-ELOQ12)*EXP(RHSP2*DTTURBL)
    1146. 

  1146.           ELOQ13=AMAX1(ELOQ13,EPS1)
    1147. 

  1147. !
    1148. 

  1148. !----------------------------------------------------------------------
    1149. 

  1149. !***  TWO ITERATIONS IS ENOUGH IN MOST CASES ...
    1150. 

  1150. !----------------------------------------------------------------------
    1151. 

  1151. !
    1152. 

  1152.           ELOQN=ELOQ13
    1153. 

  1153. !
    1154. 

  1154.           IF(ELOQN>EPS1)THEN
    1155. 

  1155.             Q2(K)=EL(K)*EL(K)/(ELOQN*ELOQN)
    1156. 

  1156.             Q2(K)=AMAX1(Q2(K),EPSQ2)
    1157. 

  1157.             IF(Q2(K)==EPSQ2)THEN
    1158. 

  1158.               EL(K)=EPSL
    1159. 

  1159.             ENDIF
    1160. 

  1160.           ELSE
    1161. 

  1161.             Q2(K)=EPSQ2
    1162. 

  1162.             EL(K)=EPSL
    1163. 

  1163.           ENDIF
    1164. 

  1164. !
    1165. 

  1165. !----------------------------------------------------------------------
    1166. 

  1166. !***  END OF TURBULENT BRANCH
    1167. 

  1167. !----------------------------------------------------------------------
    1168. 

  1168. !
    1169. 

  1169.         ENDIF
    1170. 

  1170. !----------------------------------------------------------------------
    1171. 

  1171. !***  END OF PRODUCTION/DISSIPATION LOOP
    1172. 

  1172. !----------------------------------------------------------------------
    1173. 

  1173. !
    1174. 

  1174.       ENDDO main_integration
    1175. 

  1175. !
    1176. 

  1176. !----------------------------------------------------------------------
    1177. 

  1177. !***  LOWER BOUNDARY CONDITION FOR Q2
    1178. 

  1178. !----------------------------------------------------------------------
    1179. 

  1179. !
    1180. 

  1180.       Q2(LMH)=AMAX1(B1**(2./3.)*USTAR*USTAR,EPSQ2)
    1181. 

  1181. !----------------------------------------------------------------------
    1182. 

  1182. !
    1183. 

  1183.       END SUBROUTINE PRODQ2
    1184. 

  1184. !
    1185. 

  1185. !----------------------------------------------------------------------
    1186. 

  1186. !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    1187. 

  1187. !----------------------------------------------------------------------
    1188. 

  1188.                            SUBROUTINE DIFCOF                           &
    1189. 

  1189. !   ******************************************************************
    1190. 

  1190. !   *                                                                *
    1191. 

  1191. !   *                LEVEL 2.5 DIFFUSION COEFFICIENTS                *
    1192. 

  1192. !   *                                                                *
    1193. 

  1193. !   ******************************************************************
    1194. 

  1194.      &(LMH,LMXL,GM,GH,EL,T,Q2,Z,AKM,AKH                                &
    1195. 

  1195.      &,IDS,IDE,JDS,JDE,KDS,KDE                                         &
    1196. 

  1196.      &,IMS,IME,JMS,JME,KMS,KME                                         &
    1197. 

  1197.      &,ITS,ITE,JTS,JTE,KTS,KTE,PRINT_DIAG)   ! debug
    1198. 

  1198. !----------------------------------------------------------------------
    1199. 

  1199. !
    1200. 

  1200.       IMPLICIT NONE
    1201. 

  1201. !
    1202. 

  1202. !----------------------------------------------------------------------
    1203. 

  1203.       INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                    &
    1204. 

  1204.      &                     ,IMS,IME,JMS,JME,KMS,KME                    &
    1205. 

  1205.      &                     ,ITS,ITE,JTS,JTE,KTS,KTE
    1206. 

  1206. !
    1207. 

  1207.       INTEGER,INTENT(IN) :: LMH,LMXL
    1208. 

  1208. !
    1209. 

  1209.       REAL,DIMENSION(KTS:KTE),INTENT(IN) :: Q2,T
    1210. 

  1210.       REAL,DIMENSION(KTS:KTE-1),INTENT(IN) :: EL,GH,GM
    1211. 

  1211.       REAL,DIMENSION(KTS:KTE+1),INTENT(IN) :: Z
    1212. 

  1212. !
    1213. 

  1213.       REAL,DIMENSION(KTS:KTE-1),INTENT(OUT) :: AKH,AKM
    1214. 

  1214. !----------------------------------------------------------------------
    1215. 

  1215. !***
    1216. 

  1216. !***  LOCAL VARIABLES
    1217. 

  1217. !***
    1218. 

  1218.       INTEGER :: K,KINV
    1219. 

  1219. !
    1220. 

  1220.       REAL :: ADEN,AKMIN,BDEN,BESH,BESM,CDEN,D2T,ELL,ELOQ2,ELOQ4,ELQDZ &
    1221. 

  1221.      &       ,ESH,ESM,GHL,GML,Q1L,RDEN,RDZ
    1222. 

  1222. !
    1223. 

  1223. !*** Begin debugging
    1224. 

  1224.       INTEGER,INTENT(IN) :: PRINT_DIAG
    1225. 

  1225. !     REAL :: D2Tmin
    1226. 

  1226. !*** End debugging
    1227. 

  1227. !
    1228. 

  1228. !----------------------------------------------------------------------
    1229. 

  1229. !**********************************************************************
    1230. 

  1230. !----------------------------------------------------------------------
    1231. 

  1231. !
    1232. 

  1232.       DO K=1,LMH-1
    1233. 

  1233.         ELL=EL(K)
    1234. 

  1234. !
    1235. 

  1235.         ELOQ2=ELL*ELL/Q2(K)
    1236. 

  1236.         ELOQ4=ELOQ2*ELOQ2
    1237. 

  1237. !
    1238. 

  1238.         GML=GM(K)
    1239. 

  1239.         GHL=GH(K)
    1240. 

  1240. !
    1241. 

  1241. !----------------------------------------------------------------------
    1242. 

  1242. !***  COEFFICIENTS OF THE TERMS IN THE DENOMINATOR
    1243. 

  1243. !----------------------------------------------------------------------
    1244. 

  1244. !
    1245. 

  1245.         ADEN=(ADNM*GML+ADNH*GHL)*GHL
    1246. 

  1246.         BDEN= BDNM*GML+BDNH*GHL
    1247. 

  1247.         CDEN= 1.
    1248. 

  1248. !
    1249. 

  1249. !----------------------------------------------------------------------
    1250. 

  1250. !***  COEFFICIENTS FOR THE SM DETERMINANT
    1251. 

  1251. !----------------------------------------------------------------------
    1252. 

  1252. !
    1253. 

  1253.         BESM=BSMH*GHL
    1254. 

  1254. !
    1255. 

  1255. !----------------------------------------------------------------------
    1256. 

  1256. !***  COEFFICIENTS FOR THE SH DETERMINANT
    1257. 

  1257. !----------------------------------------------------------------------
    1258. 

  1258. !
    1259. 

  1259.         BESH=BSHM*GML+BSHH*GHL
    1260. 

  1260. !
    1261. 

  1261. !----------------------------------------------------------------------
    1262. 

  1262. !***  1./DENOMINATOR
    1263. 

  1263. !----------------------------------------------------------------------
    1264. 

  1264. !
    1265. 

  1265.         RDEN=1./(ADEN*ELOQ4+BDEN*ELOQ2+CDEN)
    1266. 

  1266. !
    1267. 

  1267. !----------------------------------------------------------------------
    1268. 

  1268. !***  SM AND SH
    1269. 

  1269. !----------------------------------------------------------------------
    1270. 

  1270. !
    1271. 

  1271.         ESM=(BESM*ELOQ2+CESM)*RDEN
    1272. 

  1272.         ESH=(BESH*ELOQ2+CESH)*RDEN
    1273. 

  1273. !
    1274. 

  1274. !----------------------------------------------------------------------
    1275. 

  1275. !***  DIFFUSION COEFFICIENTS
    1276. 

  1276. !----------------------------------------------------------------------
    1277. 

  1277. !
    1278. 

  1278.         RDZ=2./(Z(K)-Z(K+2))
    1279. 

  1279.         Q1L=SQRT(Q2(K))
    1280. 

  1280.         ELQDZ=ELL*Q1L*RDZ
    1281. 

  1281.         AKM(K)=ELQDZ*ESM
    1282. 

  1282.         AKH(K)=ELQDZ*ESH
    1283. 

  1283. !----------------------------------------------------------------------
    1284. 

  1284.       ENDDO
    1285. 

  1285. !----------------------------------------------------------------------
    1286. 

  1286. !
    1287. 

  1287. !----------------------------------------------------------------------
    1288. 

  1288. !***  INVERSIONS
    1289. 

  1289. !----------------------------------------------------------------------
    1290. 

  1290. !
    1291. 

  1291. !     IF(LMXL==LMH)THEN
    1292. 

  1292. !       KINV=LMH
    1293. 

  1293. !       D2Tmin=0.
    1294. 

  1294. !
    1295. 

  1295. !       DO K=LMH/2,LMH-1
    1296. 

  1296. !         D2T=T(K-1)-2.*T(K)+T(K+1)
    1297. 

  1297. !         IF(D2T<D2Tmin)THEN
    1298. 

  1298. !           D2Tmin=D2T
    1299. 

  1299. !           IF(D2T<0)KINV=K
    1300. 

  1300. !         ENDIF
    1301. 

  1301. !       ENDDO
    1302. 

  1302. !
    1303. 

  1303. !       IF(KINV<LMH)THEN
    1304. 

  1304. !         DO K=KINV-1,LMH-1
    1305. 

  1305. !           RDZ=2./(Z(K)-Z(K+2))
    1306. 

  1306. !           AKMIN=0.5*RDZ
    1307. 

  1307. !           AKM(K)=MAX(AKM(K),AKMIN)
    1308. 

  1308. !           AKH(K)=MAX(AKH(K),AKMIN)
    1309. 

  1309. !         ENDDO
    1310. 

  1310. !
    1311. 

  1311. !*** Begin debugging
    1312. 

  1312. !         IF(PRINT_DIAG>0)THEN
    1313. 

  1313. !           write(6,"(a,3i3)") '{turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
    1314. 

  1314. !           write(6,"(a,3i3)") '}turb1 lmxl,lmh,kinv=',lmxl,lmh,kinv
    1315. 

  1315. !           IF(PRINT_DIAG==1)THEN
    1316. 

  1316. !             write(6,"(a)") &
    1317. 

  1317. !               '{turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
    1318. 

  1318. !           ELSE
    1319. 

  1319. !             write(6,"(a)") &
    1320. 

  1320. !               '}turb3 k, t, d2t, rdz, z(k), z(k+2), akmin, akh '
    1321. 

  1321. !           ENDIF
    1322. 

  1322. !           DO K=LMH-1,KINV-1,-1
    1323. 

  1323. !             D2T=T(K-1)-2.*T(K)+T(K+1)
    1324. 

  1324. !             RDZ=2./(Z(K)-Z(K+2))
    1325. 

  1325. !             AKMIN=0.5*RDZ
    1326. 

  1326. !             IF(PRINT_DIAG==1)THEN
    1327. 

  1327. !               write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '{turb3 ' &
    1328. 

  1328. !               ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
    1329. 

  1329. !             ELSE
    1330. 

  1330. !               write(6,"(a,i3,f8.3,2e12.5,2f9.2,2e12.5)") '}turb3 ' &
    1331. 

  1331. !               ,k,t(k)-273.15,d2t,rdz,z(k),z(k+2),akmin,akh(k)
    1332. 

  1332. !             ENDIF
    1333. 

  1333. !           ENDDO
    1334. 

  1334. !         ENDIF     !- IF (print_diag > 0) THEN
    1335. 

  1335. !       ENDIF       !- IF(KINV<LMH)THEN
    1336. 

  1336. !*** End debugging
    1337. 

  1337. !
    1338. 

  1338. !     ENDIF
    1339. 

  1339. !----------------------------------------------------------------------
    1340. 

  1340. !
    1341. 

  1341.       END SUBROUTINE DIFCOF
    1342. 

  1342. !
    1343. 

  1343. !----------------------------------------------------------------------
    1344. 

  1344. !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    1345. 

  1345. !----------------------------------------------------------------------
    1346. 

  1346.                            SUBROUTINE VDIFQ                            &
    1347. 

  1347. !   ******************************************************************
    1348. 

  1348. !   *                                                                *
    1349. 

  1349. !   *               VERTICAL DIFFUSION OF Q2 (TKE)                   *
    1350. 

  1350. !   *                                                                *
    1351. 

  1351. !   *********************************************…
    1352.

Large files files are truncated, but you can click here to view the full file