PageRenderTime 63ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_bl_acm.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1210 lines | 722 code | 153 blank | 335 comment | 1 complexity | 596dc859b7d9c664a4454fbd0a5733c4 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. MODULE module_bl_acm
    4. 

  4. 

  5. !  USE module_model_constants
    6. 

  6. 

  7.   REAL, PARAMETER      :: RIC    = 0.25                ! critical Richardson number
    8. 

  8.   REAL, PARAMETER      :: CRANKP = 0.5                 ! CRANK-NIC PARAMETER
    9. 

  9. 

  10. CONTAINS
    11. 

  11. 

  12. !-----------------------------------------------------------------------
    13. 

  13. !-----------------------------------------------------------------------
    14. 

  14.    SUBROUTINE ACMPBL(XTIME,    DTPBL,    ZNW,   SIGMAH,               &
    15. 

  15.                      U3D,      V3D,      PP3D,  DZ8W, TH3D, T3D,      &
    16. 

  16.                      QV3D,     QC3D,     QI3D,  RR3D,                 &
    17. 

  17.                      UST,      HFX,      QFX,   TSK,                  &
    18. 

  18.                      PSFC,     EP1,      G,                           &
    19. 

  19.                      ROVCP,    RD,       CPD,                         &
    20. 

  20.                      PBLH,     KPBL2D,   REGIME,                      &
    21. 

  21.                      GZ1OZ0,   WSPD,     PSIM, MUT,                   &
    22. 

  22.                      RUBLTEN,  RVBLTEN,  RTHBLTEN,                    &
    23. 

  23.                      RQVBLTEN, RQCBLTEN, RQIBLTEN,                    &
    24. 

  24.                      ids,ide, jds,jde, kds,kde,                       &
    25. 

  25.                      ims,ime, jms,jme, kms,kme,                       &
    26. 

  26.                      its,ite, jts,jte, kts,kte)
    27. 

  27. !-----------------------------------------------------------------------
    28. 

  28. !-----------------------------------------------------------------------
    29. 

  29. 

  30. !   THIS MODULE COMPUTES VERTICAL MIXING IN AND ABOVE THE PBL ACCORDING TO 
    31. 

  31. !   THE ASYMMETRICAL CONVECTIVE MODEL, VERSION 2  (ACM2), WHICH IS A COMBINED 
    32. 

  32. !   LOCAL NON-LOCAL CLOSURE SCHEME BASED ON THE ORIGINAL ACM (PLEIM AND CHANG 1992)
    33. 

  33. !
    34. 

  34. !   REFERENCES: 
    35. 

  35. !   Pleim (2007) A combined local and non-local closure model for the atmospheric
    36. 

  36. !                boundary layer.  Part1: Model description and testing.  
    37. 

  37. !                JAMC, 46, 1383-1395
    38. 

  38. !   Pleim (2007) A combined local and non-local closure model for the atmospheric
    39. 

  39. !                boundary layer.  Part2: Application and evaluation in a mesoscale
    40. 

  40. !                meteorology model.  JAMC, 46, 1396-1409
    41. 

  41. !
    42. 

  42. !  REVISION HISTORY:
    43. 

  43. !     AX        3/2005   - developed WRF version based on the MM5 PX LSM
    44. 

  44. !     RG and JP 7/2006   - Finished WRF adaptation
    45. 

  45. !     JP 12/2011 12/2011 - ACM2 modified so it's not dependent on first layer thickness. 
    46. 

  46. !
    47. 

  47. !**********************************************************************
    48. 

  48. !   ARGUMENT LIST:
    49. 

  49. !
    50. 

  50. !... Inputs:
    51. 

  51. !-- XTIME           Time since simulation start (min)
    52. 

  52. !-- DTPBL           PBL time step
    53. 

  53. !-- ZNW             Sigma at full layer
    54. 

  54. !-- SIGMAH          Sigma at half layer
    55. 

  55. !-- U3D             3D u-velocity interpolated to theta points (m/s)
    56. 

  56. !-- V3D             3D v-velocity interpolated to theta points (m/s)
    57. 

  57. !-- PP3D            Pressure at half levels (Pa)
    58. 

  58. !-- DZ8W            dz between full levels (m)
    59. 

  59. !-- TH3D            Potential Temperature (K)
    60. 

  60. !-- T3D             Temperature (K)
    61. 

  61. !-- QV3D            3D water vapor mixing ratio (Kg/Kg)
    62. 

  62. !-- QC3D            3D cloud mixing ratio (Kg/Kg)
    63. 

  63. !-- QI3D            3D ice mixing ratio (Kg/Kg)
    64. 

  64. !-- RR3D            3D dry air density (kg/m^3)
    65. 

  65. !-- UST             Friction Velocity (m/s)
    66. 

  66. !-- HFX		    Upward heat flux at the surface (w/m^2)
    67. 

  67. !-- QFX		    Upward moisture flux at the surface (Kg/m^2/s)
    68. 

  68. !-- TSK             Surface temperature (K)
    69. 

  69. !-- PSFC            Pressure at the surface (Pa)
    70. 

  70. !-- EP1             Constant for virtual temperature (r_v/r_d-1) (dimensionless)
    71. 

  71. !-- G               Gravity (m/s^2)
    72. 

  72. !-- ROVCP           r/cp
    73. 

  73. !-- RD              gas constant for dry air (j/kg/k)
    74. 

  74. !-- CPD             heat capacity at constant pressure for dry air (j/kg/k)
    75. 

  75. !-- GZ1OZ0          log(z/z0) where z0 is roughness length
    76. 

  76. !-- WSPD            wind speed at lowest model level (m/s)
    77. 

  77. !-- PSIM            similarity stability function for momentum
    78. 

  78. !-- MUT             Total Mu : Psfc - Ptop
    79. 

  79.  
    80. 

  80. !-- ids             start index for i in domain
    81. 

  81. !-- ide             end index for i in domain
    82. 

  82. !-- jds             start index for j in domain
    83. 

  83. !-- jde             end index for j in domain
    84. 

  84. !-- kds             start index for k in domain
    85. 

  85. !-- kde             end index for k in domain
    86. 

  86. !-- ims             start index for i in memory
    87. 

  87. !-- ime             end index for i in memory
    88. 

  88. !-- jms             start index for j in memory
    89. 

  89. !-- jme             end index for j in memory
    90. 

  90. !-- kms             start index for k in memory
    91. 

  91. !-- kme             end index for k in memory
    92. 

  92. !-- jts             start index for j in tile
    93. 

  93. !-- jte             end index for j in tile
    94. 

  94. !-- kts             start index for k in tile
    95. 

  95. !-- kte             end index for k in tile
    96. 

  96. !
    97. 

  97. !... Outputs: 
    98. 

  98. !-- PBLH            PBL height (m)
    99. 

  99. !-- KPBL2D          K index for PBL layer
    100. 

  100. !-- REGIME          Flag indicating PBL regime (stable, unstable, etc.)
    101. 

  101. !-- RUBLTEN         U tendency due to PBL parameterization (m/s^2)
    102. 

  102. !-- RVBLTEN         V tendency due to PBL parameterization (m/s^2)
    103. 

  103. !-- RTHBLTEN        Theta tendency due to PBL parameterization (K/s)
    104. 

  104. !-- RQVBLTEN        Qv tendency due to PBL parameterization (kg/kg/s)
    105. 

  105. !-- RQCBLTEN        Qc tendency due to PBL parameterization (kg/kg/s)
    106. 

  106. !-- RQIBLTEN        Qi tendency due to PBL parameterization (kg/kg/s)
    107. 

  107. !-----------------------------------------------------------------------
    108. 

  108. !-----------------------------------------------------------------------
    109. 

  109.      IMPLICIT NONE
    110. 

  110. 

  111. !.......Arguments
    112. 

  112. ! DECLARATIONS - INTEGER
    113. 

  113.     INTEGER,  INTENT(IN   )   ::      ids,ide, jds,jde, kds,kde, &
    114. 

  114.                                       ims,ime, jms,jme, kms,kme, &
    115. 

  115.                                       its,ite, jts,jte, kts,kte, XTIME
    116. 

  116. 

  117. ! DECLARATIONS - REAL
    118. 

  118.     REAL,                                INTENT(IN)  ::  DTPBL, EP1,   &
    119. 

  119.                                                         G, ROVCP, RD, CPD
    120. 

  120. 

  121.     REAL,    DIMENSION( kms:kme ),       INTENT(IN)  :: ZNW, SIGMAH
    122. 

  122. 

  123.     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
    124. 

  124.              INTENT(IN) ::                              U3D, V3D,            &
    125. 

  125.                                                         PP3D, DZ8W, T3D,     &
    126. 

  126.                                                         QV3D, QC3D, QI3D,    &
    127. 

  127.                                                         RR3D, TH3D
    128. 

  128. 

  129.     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: PSIM, GZ1OZ0,     &
    130. 

  130.                                                           HFX, QFX, TSK,    &
    131. 

  131.                                                           PSFC, WSPD, MUT
    132. 

  132. 

  133.     REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::  PBLH, REGIME, UST
    134. 

  134. 

  135.     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ),                         &
    136. 

  136.              INTENT(INOUT)   ::                         RUBLTEN, RVBLTEN,    &
    137. 

  137.                                                         RTHBLTEN, RQVBLTEN,  &
    138. 

  138.                                                         RQCBLTEN, RQIBLTEN
    139. 

  139. 

  140.     INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(OUT  ) ::  KPBL2D
    141. 

  141. 

  142. !... Local Variables
    143. 

  143. 

  144. !... Integer
    145. 

  145.       INTEGER :: J, K
    146. 

  146. !... Real
    147. 

  147.       REAL, DIMENSION( kts:kte ) :: DSIGH, DSIGHI, DSIGFI
    148. 

  148.       REAL, DIMENSION( 0:kte )   :: SIGMAF
    149. 

  149.       REAL  RDT
    150. 

  150.       REAL, PARAMETER :: KARMAN = 0.4
    151. 

  151. !...
    152. 

  152. 

  153.    RDT = 1.0 / DTPBL
    154. 

  154. 

  155.    DO K = 1, kte
    156. 

  156.      SIGMAF(K-1) = ZNW(K)
    157. 

  157.    ENDDO
    158. 

  158.    SIGMAF(kte) = 0.0
    159. 

  159. 

  160.    DO K = 1, kte
    161. 

  161.      DSIGH(K)  = SIGMAF(K) - SIGMAF(K-1)
    162. 

  162.      DSIGHI(K) = 1.0 / DSIGH(K)
    163. 

  163.    ENDDO
    164. 

  164. 

  165.    DO K = kts,kte-1
    166. 

  166.      DSIGFI(K) = 1.0 / (SIGMAH(K+1) - SIGMAH(K))
    167. 

  167.    ENDDO
    168. 

  168. 

  169.    DSIGFI(kte) = DSIGFI(kte-1)
    170. 

  170.    
    171. 

  171.    DO j = jts,jte
    172. 

  172.       CALL ACM2D(j=J,xtime=XTIME, dtpbl=DTPBL, sigmaf=SIGMAF, sigmah=SIGMAH    &
    173. 

  173.               ,dsigfi=DSIGFI,dsighi=DSIGHI,dsigh=DSIGH             &
    174. 

  174.               ,us=u3d(ims,kms,j),vs=v3d(ims,kms,j)                 &
    175. 

  175.               ,theta=th3d(ims,kms,j),tt=t3d(ims,kms,j)             &
    176. 

  176.               ,qvs=qv3d(ims,kms,j),qcs=qc3d(ims,kms,j)             &
    177. 

  177.               ,qis=qi3d(ims,kms,j),dzf=DZ8W(ims,kms,j)             &
    178. 

  178.               ,densx=RR3D(ims,kms,j)                               &
    179. 

  179.               ,utnp=rublten(ims,kms,j),vtnp=rvblten(ims,kms,j)     &
    180. 

  180.               ,ttnp=rthblten(ims,kms,j),qvtnp=rqvblten(ims,kms,j)  &
    181. 

  181.               ,qctnp=rqcblten(ims,kms,j),qitnp=rqiblten(ims,kms,j) &
    182. 

  182.               ,cpd=cpd,g=g,rovcp=rovcp,rd=rd,rdt=rdt               &
    183. 

  183.               ,psfcpa=psfc(ims,j),ust=ust(ims,j)                   &
    184. 

  184.               ,pbl=pblh(ims,j)                                     &
    185. 

  185.               ,regime=regime(ims,j),psim=psim(ims,j)               &
    186. 

  186.               ,hfx=hfx(ims,j),qfx=qfx(ims,j)                       &
    187. 

  187.               ,tg=tsk(ims,j),gz1oz0=gz1oz0(ims,j)                  &
    188. 

  188.               ,wspd=wspd(ims,j) ,klpbl=kpbl2d(ims,j)               &
    189. 

  189.               ,mut=mut(ims,j)                                      &
    190. 

  190.               ,ep1=ep1,karman=karman                               &
    191. 

  191.               ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde   &
    192. 

  192.               ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme   &
    193. 

  193.               ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
    194. 

  194.    ENDDO
    195. 

  195. 

  196.    END SUBROUTINE ACMPBL
    197. 

  197. !-----------------------------------------------------------------------
    198. 

  198. !-----------------------------------------------------------------------
    199. 

  199. 

  200. 

  201. !-----------------------------------------------------------------------
    202. 

  202. !-----------------------------------------------------------------------
    203. 

  203.    SUBROUTINE ACM2D(j,XTIME, DTPBL, sigmaf, sigmah          &
    204. 

  204.               ,dsigfi,dsighi,dsigh                          &
    205. 

  205.               ,us,vs,theta,tt,qvs,qcs,qis                   &
    206. 

  206.               ,dzf,densx,utnp,vtnp,ttnp,qvtnp,qctnp,qitnp   &
    207. 

  207.               ,cpd,g,rovcp,rd,rdt,psfcpa,ust                &
    208. 

  208.               ,pbl,regime,psim                              &
    209. 

  209.               ,hfx,qfx,tg,gz1oz0,wspd ,klpbl                &
    210. 

  210.               ,mut                                          &
    211. 

  211.               ,ep1,karman                                   &
    212. 

  212.               ,ids,ide, jds,jde, kds,kde   &
    213. 

  213.               ,ims,ime, jms,jme, kms,kme   &
    214. 

  214.               ,its,ite, jts,jte, kts,kte   )
    215. 

  215. 

  216. !-----------------------------------------------------------------------
    217. 

  217. !-----------------------------------------------------------------------
    218. 

  218.       IMPLICIT NONE
    219. 

  219. 

  220. !.......Arguments
    221. 

  221. 

  222. !... Real
    223. 

  223.       REAL, DIMENSION( 0:kte ),             INTENT(IN)  :: SIGMAF
    224. 

  224.       REAL, DIMENSION( kms:kme ),           INTENT(IN)  :: SIGMAH
    225. 

  225.       REAL, DIMENSION( kts:kte ),           INTENT(IN)  :: DSIGH, DSIGHI, DSIGFI
    226. 

  226.       REAL ,                                INTENT(IN)  :: DTPBL, G, RD,ep1,karman,CPD,ROVCP,RDT
    227. 

  227.       REAL , DIMENSION( ims:ime ),          INTENT(INOUT)  :: PBL, UST
    228. 

  228.       
    229. 

  229.       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA, TT,   &
    230. 

  230.                                                            QVS, QCS, QIS, DENSX
    231. 

  231.       REAL,  DIMENSION( ims:ime, kms:kme ), intent(in)  :: DZF
    232. 

  232.       REAL,  DIMENSION( ims:ime, kms:kme ), intent(inout)   ::  utnp, &
    233. 

  233. 							        vtnp, &
    234. 

  234. 							        ttnp, &
    235. 

  235. 							        qvtnp, &
    236. 

  236. 							        qctnp, &
    237. 

  237. 							        qitnp
    238. 

  238.       real,     dimension( ims:ime ), intent(in   )   ::   psfcpa
    239. 

  239.       real,     dimension( ims:ime ), intent(in   )   ::   tg
    240. 

  240.       real,     dimension( ims:ime ), intent(inout)   ::   regime
    241. 

  241.       real,     dimension( ims:ime ), intent(in)      ::   wspd, psim, gz1oz0
    242. 

  242.       real,     dimension( ims:ime ), intent(in)      ::   hfx, qfx
    243. 

  243.       real,     dimension( ims:ime ), intent(in)      ::   mut
    244. 

  244. 

  245. !... Integer
    246. 

  246.       INTEGER, DIMENSION( ims:ime ),       INTENT(OUT):: KLPBL
    247. 

  247.       INTEGER,  INTENT(IN)      ::      XTIME
    248. 

  248.       integer,  intent(in   )   ::      ids,ide, jds,jde, kds,kde, &
    249. 

  249.                                         ims,ime, jms,jme, kms,kme, &
    250. 

  250.                                         its,ite, jts,jte, kts,kte, j
    251. 

  251. !--------------------------------------------------------------------
    252. 

  252. !--Local 
    253. 

  253.       INTEGER I, K     
    254. 

  254.       INTEGER :: KPBLHT
    255. 

  255.       INTEGER, DIMENSION( its:ite ) :: KPBLH, NOCONV
    256. 

  256. 

  257. !... Real
    258. 

  258.       REAL    ::  TVCON, WSS, TCONV, TH1, TOG, DTMP, WSSQ
    259. 

  259.       REAL    ::  psix, THV1
    260. 

  260.       REAL, DIMENSION( its:ite )          :: FINT, PSTAR, CPAIR
    261. 

  261.       REAL, DIMENSION( its:ite, kts:kte ) :: THETAV, RIB,             &
    262. 

  262.                                              EDDYZ, UX, VX, THETAX,   &
    263. 

  263.                                              QVX, QCX, QIX, ZA
    264. 

  264.       REAL, DIMENSION( its:ite, 0:kte )   :: ZF
    265. 

  265.       REAL,    DIMENSION( its:ite)           :: WST, TST, QST, USTM, TSTV
    266. 

  266.       REAL,    DIMENSION( its:ite )          :: PBLSIG, MOL
    267. 

  267.       REAL    ::  FINTT, ZMIX, UMIX, VMIX
    268. 

  268.       REAL    ::  TMPFX, TMPVTCON, TMPP, TMPTHS, TMPTH1, TMPVCONV, WS1, DTH
    269. 

  269. 	REAL    ::  A,TST12,RL,ZFUNC
    270. 

  270. !    REAL, PARAMETER :: KARMAN = 0.4
    271. 

  271. 

  272. !... Integer
    273. 

  273.       INTEGER :: KL, jtf, ktf, itf, KMIX, KSRC
    274. 

  274. !...
    275. 

  275.         character*256 :: message
    276. 

  276. !-----initialize vertical tendencies and
    277. 

  277. 

  278.       DO i = its,ite
    279. 

  279.         DO k = kts,kte
    280. 

  280.           utnp(i,k) = 0.
    281. 

  281.           vtnp(i,k) = 0.
    282. 

  282.           ttnp(i,k) = 0.
    283. 

  283.         ENDDO
    284. 

  284.       ENDDO
    285. 

  285. 

  286.       DO k = kts,kte
    287. 

  287.         DO i = its,ite
    288. 

  288.           qvtnp(i,k) = 0.
    289. 

  289.         ENDDO
    290. 

  290.       ENDDO
    291. 

  291. 

  292.       DO k = kts,kte
    293. 

  293.         DO i = its,ite
    294. 

  294.           qctnp(i,k) = 0.
    295. 

  295.           qitnp(i,k) = 0.
    296. 

  296.         ENDDO
    297. 

  297.       ENDDO
    298. 

  298. 

  299. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    300. 

  300. !!!  Compute Micromet Scaling variables, not availiable in WRF for ACM
    301. 

  301.      DO I = its,ite
    302. 

  302.            CPAIR(I)  = CPD * (1.0 + 0.84 * QVS(I,1))                    ! J/(K KG)
    303. 

  303.            TMPFX     = HFX(I)  / (cpair(i) * DENSX(I,1))
    304. 

  304.            TMPVTCON  = 1.0 + EP1 * QVS(I,1)                             ! COnversion factor for virtual temperature
    305. 

  305.            WS1       = SQRT(US(I,1)**2 + VS(I,1)**2)                    ! Level 1 wind speed
    306. 

  306.            TST(I)    = -TMPFX / UST(I)
    307. 

  307.            QST(I)    = -QFX(I) / (UST(I)*DENSX(I,1))
    308. 

  308.            USTM(I)   = UST(I) * WS1 / wspd(i)
    309. 

  309.            THV1      = TMPVTCON * THETA(I,1) 
    310. 

  310.            TSTV(I)   = TST(I)*TMPVTCON + THV1*EP1*QST(I)
    311. 

  311.            IF(ABS(TSTV(I)).LT.1.0E-6) THEN
    312. 

  312.              TSTV(I) = SIGN(1.0E-6,TSTV(I))
    313. 

  313.            ENDIF
    314. 

  314.            MOL(I)    = THV1 * UST(i)**2/(KARMAN*G*TSTV(I))
    315. 

  315.            WST(I)    = UST(I) * (PBL(I)/(KARMAN*ABS(MOL(I)))) ** 0.333333       
    316. 

  316.            PSTAR(I)  =  MUT(I)/1000.                                     ! P* in cb 
    317. 

  317.      ENDDO
    318. 

  318. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    319. 

  319. 

  320. !... Compute PBL height
    321. 

  321. 

  322. !... compute the height of full- and half-sigma level above ground level
    323. 

  323.      DO I = its,ite
    324. 

  324.        ZF(I,0)    = 0.0
    325. 

  325.        KLPBL(I)   = 1
    326. 

  326.      ENDDO
    327. 

  327. 

  328.      DO K = kts, kte
    329. 

  329.        DO I = its,ite
    330. 

  330.          ZF(I,K) = DZF(I,K) + ZF(I,K-1)
    331. 

  331.          ZA(I,K) = 0.5 * (ZF(I,K) + ZF(I,K-1))
    332. 

  332.        ENDDO
    333. 

  333.      ENDDO
    334. 

  334. 

  335.      DO K = kts, kte
    336. 

  336.        DO I = its,ite
    337. 

  337.          TVCON       = 1.0 + EP1 * QVS(I,K)
    338. 

  338.          THETAV(I,K) = THETA(I,K) * TVCON
    339. 

  339.        ENDDO
    340. 

  340.      ENDDO
    341. 

  341. 

  342. 

  343. !...  COMPUTE PBL WHERE RICHARDSON NUMBER = RIC (0.25) HOLTSLAG ET AL 1990  
    344. 

  344.      DO 100 I = its,ite
    345. 

  345.        DO K = 1,kte
    346. 

  346.          KSRC = K
    347. 

  347.          IF (SIGMAF(K).lT.0.9955) GO TO 69
    348. 

  348.        ENDDO
    349. 

  349. 69     CONTINUE
    350. 

  350.        TH1 = 0.0
    351. 

  351.        DO K = 1,KSRC
    352. 

  352.          TH1 = TH1 + THETAV(I,K)  
    353. 

  353.        ENDDO  
    354. 

  354.        TH1 = TH1/KSRC
    355. 

  355.        IF(MOL(I).LT.0.0 .AND. XTIME.GT.1) then
    356. 

  356.          WSS   = (UST(I) ** 3 + 0.6 * WST(I) ** 3) ** 0.33333
    357. 

  357.          TCONV = -8.5 * UST(I) * TSTV(I) / WSS
    358. 

  358.          TH1   = TH1 + TCONV
    359. 

  359.        ENDIF
    360. 

  360. 

  361. 99     KMIX = 1
    362. 

  362.        DO K = 1,kte
    363. 

  363.          DTMP   = THETAV(I,K) - TH1
    364. 

  364.          IF (DTMP.LT.0.0) KMIX = K
    365. 

  365.        ENDDO
    366. 

  366.        IF(KMIX.GT.1) THEN
    367. 

  367.          FINTT = (TH1 - THETAV(I,KMIX)) / (THETAV(I,KMIX+1)               &
    368. 

  368.                - THETAV(I,KMIX))
    369. 

  369.          ZMIX = FINTT * (ZA(I,KMIX+1)-ZA(I,KMIX)) + ZA(I,KMIX)
    370. 

  370.          UMIX = FINTT * (US(I,KMIX+1)-US(I,KMIX)) + US(I,KMIX)
    371. 

  371.          VMIX = FINTT * (VS(I,KMIX+1)-VS(I,KMIX)) + VS(I,KMIX)
    372. 

  372.        ELSE
    373. 

  373.          ZMIX = ZA(I,1)
    374. 

  374.          UMIX = US(I,1)
    375. 

  375.          VMIX = VS(I,1)
    376. 

  376.        ENDIF
    377. 

  377.        DO K = KMIX,kte
    378. 

  378.          DTMP   = THETAV(I,K) - TH1
    379. 

  379.          TOG = 0.5 * (THETAV(I,K) + TH1) / G
    380. 

  380.          WSSQ = (US(I,K)-UMIX)**2                                     &
    381. 

  381.               + (VS(I,K)-VMIX)**2
    382. 

  382.          IF (KMIX == 1) WSSQ = WSSQ + 100.*UST(I)*UST(I) 
    383. 

  383.          WSSQ = MAX( WSSQ, 0.1 )
    384. 

  384.          RIB(I,K) = ABS(ZA(I,K)-ZMIX) * DTMP / (TOG * WSSQ)
    385. 

  385.          IF (RIB(I,K) .GE. RIC) GO TO 201
    386. 

  386.        ENDDO
    387. 

  387. 

  388.        write (message, *)' RIBX never exceeds RIC, RIB(i,kte) = ',rib(i,5),        &
    389. 

  389.                ' THETAV(i,1) = ',thetav(i,1),' MOL=',mol(i),            &
    390. 

  390.                ' TCONV = ',TCONV,' WST = ',WST(I),                      &
    391. 

  391.                ' KMIX = ',kmix,' UST = ',UST(I),                       &
    392. 

  392.                ' TST = ',TST(I),' U,V = ',US(I,1),VS(I,1),              &
    393. 

  393.                ' I,J=',I,J
    394. 

  394.        CALL wrf_error_fatal ( message )
    395. 

  395. 201    CONTINUE
    396. 

  396. 

  397.        KPBLH(I) = K
    398. 

  398. 

  399. 100  CONTINUE
    400. 

  400. 

  401.      DO I = its,ite
    402. 

  402.        IF (KPBLH(I) .NE. 1) THEN
    403. 

  403. !---------INTERPOLATE BETWEEN LEVELS -- jp 7/93
    404. 

  404.          FINT(I) = (RIC - RIB(I,KPBLH(I)-1)) / (RIB(I,KPBLH(I)) -       &
    405. 

  405.                     RIB(I,KPBLH(I)-1))
    406. 

  406.          IF (FINT(I) .GT. 0.5) THEN
    407. 

  407.            KPBLHT  = KPBLH(I)
    408. 

  408.            FINT(I) = FINT(I) - 0.5
    409. 

  409.          ELSE
    410. 

  410.            KPBLHT  = KPBLH(I) - 1
    411. 

  411.            FINT(I) = FINT(I) + 0.5
    412. 

  412.          ENDIF
    413. 

  413.          PBL(I)  = FINT(I) * (ZF(I,KPBLHT) - ZF(I,KPBLHT-1)) +          &
    414. 

  414.                      ZF(I,KPBLHT-1)
    415. 

  415.          KLPBL(I) = KPBLHT
    416. 

  416.          PBLSIG(I)   = FINT(I) * DSIGH(KPBLHT) + SIGMAF(KPBLHT-1)    ! sigma at PBL height
    417. 

  417.        ELSE
    418. 

  418.          KLPBL(I) = 1
    419. 

  419.          PBL(I)    = ZF(I,1)                                                  
    420. 

  420.          PBLSIG(I)   = SIGMAF(1)                                             
    421. 

  421.        ENDIF
    422. 

  422. 

  423.      ENDDO
    424. 

  424. 

  425.      DO I = its,ite       
    426. 

  426.        NOCONV(I) = 0
    427. 

  427.        
    428. 

  428. ! Check for CBL and identify conv. vs. non-conv cells
    429. 

  429.        IF (PBL(I) / MOL(I) .LT. -0.02 .AND. KLPBL(I) .GT. 3        &
    430. 

  430.            .AND. THETAV(I,1) .GT. THETAV(I,2) .AND. XTIME .GT. 1) THEN
    431. 

  431.           NOCONV(I)   = 1
    432. 

  432.           REGIME(I) = 4.0                     ! FREE CONVECTIVE - ACM
    433. 

  433.        ENDIF
    434. 

  434.      ENDDO
    435. 

  435. 

  436. !... Calculate Kz
    437. 

  437.      CALL EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,                &
    438. 

  438.                 US,    VS,  TT,  THETAV, DENSX, PSTAR,              &
    439. 

  439.                 QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,              &
    440. 

  440.                 EDDYZ, its,ite, kts,kte,ims,ime, kms,kme)
    441. 

  441. 

  442.      CALL ACM(DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, J,      &
    443. 

  443.                  KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
    444. 

  444.                  TST, QST,  USTM,   EDDYZ,  DENSX,                  &
    445. 

  445.                  US,    VS,     THETA,  QVS,    QCS,    QIS,        &
    446. 

  446.                  UX,    VX,     THETAX, QVX,    QCX,    QIX,        &
    447. 

  447.                  ids,ide, jds,jde, kds,kde,                         &
    448. 

  448.                  ims,ime, jms,jme, kms,kme,                         &
    449. 

  449.                  its,ite, jts,jte, kts,kte)
    450. 

  450. 

  451. !... Calculate tendency due to PBL parameterization
    452. 

  452. 

  453.      DO K = kts, kte
    454. 

  454.        DO I = its, ite
    455. 

  455.          UTNP(I,K)  = UTNP(I,K) + (UX(I,K) - US(I,K)) * RDT
    456. 

  456.          VTNP(I,K)  = VTNP(I,K) + (VX(I,K) - VS(I,K)) * RDT
    457. 

  457.          TTNP(I,K)  = TTNP(I,K) + (THETAX(I,K) - THETA(I,K)) * RDT
    458. 

  458.          QVTNP(I,K) = QVTNP(I,K) + (QVX(I,K) - QVS(I,K)) * RDT
    459. 

  459.          QCTNP(I,K) = QCTNP(I,K) + (QCX(I,K) - QCS(I,K)) * RDT
    460. 

  460.          QITNP(I,K) = QITNP(I,K) + (QIX(I,K) - QIS(I,K)) * RDT
    461. 

  461.        ENDDO
    462. 

  462.      ENDDO
    463. 

  463. 

  464.    END SUBROUTINE ACM2D
    465. 

  465. !-----------------------------------------------------------------------
    466. 

  466. !-----------------------------------------------------------------------
    467. 

  467. 

  468. !-----------------------------------------------------------------------
    469. 

  469. !-----------------------------------------------------------------------
    470. 

  470.    SUBROUTINE ACMINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
    471. 

  471.                       RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
    472. 

  472.                       restart, allowed_to_read ,                   &
    473. 

  473.                       ids, ide, jds, jde, kds, kde,                &
    474. 

  474.                       ims, ime, jms, jme, kms, kme,                &
    475. 

  475.                       its, ite, jts, jte, kts, kte                 )
    476. 

  476. !-----------------------------------------------------------------------
    477. 

  477. !
    478. 

  478. !    This subroutine is for preparing ACM PBL variables. 
    479. 

  479. !    Called from module_physics_init.F
    480. 

  480. !
    481. 

  481. !  REVISION HISTORY:
    482. 

  482. !     AX     3/2005 - Originally developed
    483. 

  483. !-----------------------------------------------------------------------
    484. 

  484. !   ARGUMENT LIST:
    485. 

  485. !
    486. 

  486. !-----------------------------------------------------------------------
    487. 

  487. !-----------------------------------------------------------------------
    488. 

  488.      IMPLICIT NONE
    489. 

  489. !
    490. 

  490.    LOGICAL , INTENT(IN)          :: restart , allowed_to_read
    491. 

  491. 

  492.    INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
    493. 

  493.                                      ims, ime, jms, jme, kms, kme, &
    494. 

  494.                                      its, ite, jts, jte, kts, kte
    495. 

  495. 

  496.    INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR
    497. 

  497. 

  498. !   REAL , DIMENSION( kms:kme ), INTENT(IN)  :: SHALF
    499. 

  499.    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
    500. 

  500.                                                          RUBLTEN, &
    501. 

  501.                                                          RVBLTEN, &
    502. 

  502.                                                          RTHBLTEN, &
    503. 

  503.                                                          RQVBLTEN, &
    504. 

  504.                                                          RQCBLTEN, & 
    505. 

  505.                                                          RQIBLTEN
    506. 

  506. 

  507. !... Local Variables
    508. 

  508.    INTEGER :: i, j, k, itf, jtf, ktf
    509. 

  509. 

  510. !
    511. 

  511.    jtf=min0(jte,jde-1)
    512. 

  512.    ktf=min0(kte,kde-1)
    513. 

  513.    itf=min0(ite,ide-1)
    514. 

  514. 

  515.    IF(.not.restart)THEN
    516. 

  516.      DO j=jts,jtf
    517. 

  517.      DO k=kts,ktf
    518. 

  518.      DO i=its,itf
    519. 

  519.         RUBLTEN(i,k,j)=0.
    520. 

  520.         RVBLTEN(i,k,j)=0.
    521. 

  521.         RTHBLTEN(i,k,j)=0.
    522. 

  522.         RQVBLTEN(i,k,j)=0.
    523. 

  523.         RQCBLTEN(i,k,j)=0.
    524. 

  524.      ENDDO
    525. 

  525.      ENDDO
    526. 

  526.      ENDDO
    527. 

  527.    ENDIF
    528. 

  528. 

  529.    IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
    530. 

  530.       DO j=jts,jtf
    531. 

  531.       DO k=kts,ktf
    532. 

  532.       DO i=its,itf
    533. 

  533.          RQIBLTEN(i,k,j)=0.
    534. 

  534.       ENDDO
    535. 

  535.       ENDDO
    536. 

  536.       ENDDO
    537. 

  537.    ENDIF
    538. 

  538. 

  539. 

  540.    END SUBROUTINE acminit
    541. 

  541. !-----------------------------------------------------------------------
    542. 

  542. !-----------------------------------------------------------------------
    543. 

  543. 

  544. !-----------------------------------------------------------------------
    545. 

  545. !-------------------------------------------------------------------          
    546. 

  546.    SUBROUTINE EDDYX(DTPBL, ZF,  ZA,     MOL, PBL,  UST,               &
    547. 

  547.                     US,    VS,  TT,  THETAV, DENSX, PSTAR,            &
    548. 

  548.                     QVS,   QCS, QIS, DSIGFI, G, RD, CPAIR,            &
    549. 

  549.                     EDDYZ, its,ite, kts,kte,ims,ime,kms,kme )
    550. 

  550. 

  551. 

  552. !**********************************************************************
    553. 

  553. !   Two methods for computing Kz:
    554. 

  554. !   1.  Boundary scaling similar to Holtslag and Boville (1993)
    555. 

  555. !   2.  Local Kz computed as function of local Richardson # and vertical 
    556. 

  556. !       wind shear, similar to LIU & CARROLL (1996)
    557. 

  557. !
    558. 

  558. !**********************************************************************
    559. 

  559. !
    560. 

  560. !-- DTPBL           time step of the minor loop for the land-surface/pbl model
    561. 

  561. !-- ZF              height of full sigma level
    562. 

  562. !-- ZA              height of half sigma level
    563. 

  563. !-- MOL             Monin-Obukhov length in 1D form
    564. 

  564. !-- PBL             PBL height in 1D form
    565. 

  565. !-- UST             friction velocity U* in 1D form (m/s)
    566. 

  566. !-- US              U wind 
    567. 

  567. !-- VS              V wind
    568. 

  568. !-- TT              temperature
    569. 

  569. !-- THETAV          potential virtual temperature
    570. 

  570. !-- DENSX           dry air density (kg/m^3)
    571. 

  571. !-- PSTAR           P*=Psfc-Ptop
    572. 

  572. !-- QVS             water vapor mixing ratio (Kg/Kg)
    573. 

  573. !-- QCS             cloud mixing ratio (Kg/Kg)
    574. 

  574. !-- QIS             ice mixing ratio (Kg/Kg)
    575. 

  575. !-- DSIGFI          inverse of sigma layer delta
    576. 

  576. !-- G               gravity
    577. 

  577. !-- RD              gas constant for dry air (j/kg/k)
    578. 

  578. !-- CPAIR           specific heat of moist air (M^2 S^-2 K^-1)
    579. 

  579. !-- EDDYZ           eddy diffusivity KZ
    580. 

  580. !-----------------------------------------------------------------------
    581. 

  581. !-----------------------------------------------------------------------
    582. 

  582. 

  583.       IMPLICIT NONE
    584. 

  584. 

  585. !.......Arguments
    586. 

  586.   
    587. 

  587. !... Integer
    588. 

  588.       INTEGER,  INTENT(IN)   ::    its,ite, kts,kte,ims,ime,kms,kme
    589. 

  589. !... Real
    590. 

  590.       REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
    591. 

  591.       REAL ,                                INTENT(IN)  :: DTPBL, G, RD
    592. 

  592.       REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGFI
    593. 

  593.       REAL , DIMENSION( its:ite ),          INTENT(IN)  :: MOL, PSTAR, CPAIR
    594. 

  594. 

  595.       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, TT,   &
    596. 

  596.                                                            QVS, QCS, QIS, DENSX
    597. 

  597.       REAL, DIMENSION( its:ite, kts:kte ), INTENT(IN) :: ZA, THETAV
    598. 

  598.       REAL, DIMENSION( its:ite, 0:kte )  , INTENT(IN) :: ZF
    599. 

  599.       
    600. 

  600.       REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: EDDYZ
    601. 

  601. 

  602. !.......Local variables
    603. 

  603. 

  604. !... Integer
    605. 

  605.       INTEGER  :: ILX, KL, KLM, K, I
    606. 

  606. 

  607. !... Real
    608. 

  608.       REAL     :: ZOVL, PHIH, WT, ZSOL, ZFUNC, DZF, SS, GOTH, EDYZ
    609. 

  609.       REAL     :: RI, QMEAN, TMEAN, XLV, ALPH, CHI, ZK, SQL, DENSF, KZO
    610. 

  610.       REAL     :: FH
    611. 

  611. !... Parameters
    612. 

  612.       REAL, PARAMETER :: RV     = 461.5
    613. 

  613.       REAL, PARAMETER :: RC     = 0.25
    614. 

  614.       REAL, PARAMETER :: RLAM   = 80.0
    615. 

  615.       REAL, PARAMETER :: GAMH   = 16.0 !15.0  !  Holtslag and Boville (1993)
    616. 

  616.       REAL, PARAMETER :: BETAH  = 5.0   !  Holtslag and Boville (1993)
    617. 

  617.       REAL, PARAMETER :: KARMAN = 0.4
    618. 

  618.       REAL, PARAMETER :: EDYZ0  = 0.01  ! New Min Kz
    619. 

  619. !      REAL, PARAMETER :: EDYZ0  = 0.1
    620. 

  620. !--   IMVDIF      imvdif=1 for moist adiabat vertical diffusion
    621. 

  621.       INTEGER, PARAMETER :: imvdif = 1
    622. 

  622. !
    623. 

  623.       ILX = ite 
    624. 

  624.       KL  = kte
    625. 

  625.       KLM = kte - 1
    626. 

  626.       
    627. 

  627.       DO K = kts,KLM
    628. 

  628.         DO I = its,ILX
    629. 

  629.           EDYZ = 0.0
    630. 

  630.           ZOVL = 0.0
    631. 

  631.           DZF  = ZA(I,K+1) - ZA(I,K)
    632. 

  632.           KZO = EDYZ0
    633. 

  633. !--------------------------------------------------------------------------
    634. 

  634.           IF (ZF(I,K) .LT. PBL(I)) THEN
    635. 

  635.             ZOVL = ZF(I,K) / MOL(I)
    636. 

  636.             IF (ZOVL .LT. 0.0) THEN
    637. 

  637.               IF (ZF(I,K) .LT. 0.1 * PBL(I)) THEN
    638. 

  638.                 PHIH = 1.0 / SQRT(1.0 - GAMH * ZOVL)
    639. 

  639.                 WT   = UST(I) / PHIH
    640. 

  640.               ELSE
    641. 

  641.                 ZSOL = 0.1 * PBL(I) / MOL(I)
    642. 

  642.                 PHIH = 1.0 / SQRT(1.0 - GAMH * ZSOL)
    643. 

  643.                 WT   = UST(I) / PHIH
    644. 

  644.               ENDIF
    645. 

  645.             ELSE IF (ZOVL .LT. 1.0) THEN
    646. 

  646.               PHIH = 1.0 + BETAH * ZOVL
    647. 

  647.               WT   = UST(I) / PHIH
    648. 

  648.             ELSE
    649. 

  649.               PHIH = BETAH + ZOVL
    650. 

  650.               WT   = UST(I) / PHIH
    651. 

  651.             ENDIF
    652. 

  652.             ZFUNC      = ZF(I,K) * (1.0 - ZF(I,K) / PBL(I)) ** 2
    653. 

  653.             EDYZ = KARMAN * WT * ZFUNC
    654. 

  654.           ENDIF
    655. 

  655. !--------------------------------------------------------------------------
    656. 

  656.           SS   = ((US(I,K+1) - US(I,K)) ** 2 + (VS(I,K+1) - VS(I,K)) ** 2)   &
    657. 

  657.                   / (DZF * DZF) + 1.0E-9
    658. 

  658.           GOTH = 2.0 * G / (THETAV(I,K+1) + THETAV(I,K))
    659. 

  659.           RI   = GOTH * (THETAV(I,K+1) - THETAV(I,K)) / (DZF * SS)
    660. 

  660. !--------------------------------------------------------------------------
    661. 

  661. !         Adjustment to vert diff in Moist air
    662. 

  662.           IF(imvdif.eq.1)then
    663. 

  663.             IF ((QCS(I,K)+QIS(I,K)) .GT. 0.01E-3 .OR. (QCS(I,K+1)+             &
    664. 

  664.                  QIS(I,K+1)) .GT. 0.01E-3) THEN
    665. 

  665.               QMEAN = 0.5 * (QVS(I,K) + QVS(I,K+1))
    666. 

  666.               TMEAN = 0.5 * (TT(I,K) + TT(I,K+1))
    667. 

  667.               XLV   = (2.501 - 0.00237 * (TMEAN - 273.15)) * 1.E6
    668. 

  668.               ALPH  =  XLV * QMEAN / RD / TMEAN
    669. 

  669.               CHI   =  XLV * XLV * QMEAN / CPAIR(I) / RV / TMEAN / TMEAN
    670. 

  670.               RI    = (1.0 + ALPH) * (RI -G * G / SS / TMEAN / CPAIR(I) *       &
    671. 

  671.                       ((CHI - ALPH) / (1.0 + CHI)))
    672. 

  672.             ENDIF
    673. 

  673.           ENDIF
    674. 

  674. !--------------------------------------------------------------------------
    675. 

  675.             
    676. 

  676. 	        ZK  = 0.4 * ZF(I,K)
    677. 

  677.           SQL = (ZK * RLAM / (RLAM + ZK)) ** 2
    678. 

  678.             
    679. 

  679.           IF (RI .GE. 0.0) THEN
    680. 

  680. 	          IF (ZF(I,K).LT.PBL(I).AND.ZOVL.GT.0.0) THEN
    681. 

  681. 	            FH = MAX((1.-ZF(I,K)/PBL(I))**2,0.01) * PHIH **(-2)
    682. 

  682.                    SQL = ZK ** 2
    683. 

  683. 	          ELSE
    684. 

  684. 	            FH = (MAX(1.-RI/RC,0.01))**2
    685. 

  685. 	          ENDIF
    686. 

  686.             EDDYZ(I,K) = KZO + SQRT(SS) * FH * SQL
    687. 

  687.           ELSE
    688. 

  688.             EDDYZ(I,K) = KZO + SQRT(SS * (1.0 - 25.0 * RI)) * SQL
    689. 

  689.           ENDIF
    690. 

  690. 	  
    691. 

  691.           IF(EDYZ.GT.EDDYZ(I,K)) THEN
    692. 

  692.             EDDYZ(I,K) = EDYZ
    693. 

  693.           ENDIF
    694. 

  694. 

  695.           EDDYZ(I,K) = MIN(1000.0,EDDYZ(I,K))
    696. 

  696.           EDDYZ(I,K) = MAX(KZO,EDDYZ(I,K))
    697. 

  697. 

  698.           DENSF     = 0.5 * (DENSX(I,K+1) + DENSX(I,K))
    699. 

  699. 

  700.           EDDYZ(I,K) = EDDYZ(I,K) * (DENSF * G / PSTAR(I)) ** 2 *       &
    701. 

  701.                        DTPBL * DSIGFI(K)*1E-6
    702. 

  702. 

  703.         ENDDO             ! for I loop
    704. 

  704.       ENDDO               ! for k loop
    705. 

  705. !
    706. 

  706.       DO I = its,ILX
    707. 

  707.         EDDYZ(I,KL) = 0.0 ! EDDYZ(I,KLM) -- changed jp 3/08
    708. 

  708.       ENDDO
    709. 

  709. 

  710.    END SUBROUTINE EDDYX
    711. 

  711. !-----------------------------------------------------------------------
    712. 

  712. !-----------------------------------------------------------------------
    713. 

  713. 

  714. !-----------------------------------------------------------------------
    715. 

  715. !-------------------------------------------------------------------          
    716. 

  716.    SUBROUTINE ACM (DTPBL, PSTAR,  NOCONV, SIGMAF, DSIGH, DSIGHI, JX, &
    717. 

  717.                    KLPBL, PBL,   PBLSIG, MOL,  UST,                  &
    718. 

  718.                    TST, QST,  USTM,   EDDYZ,  DENSX,               &
    719. 

  719.                    US,    VS,     THETA,  QVS,    QCS,    QIS,     &
    720. 

  720.                    UX,    VX,     THETAX, QVX,    QCX,    QIX,     &
    721. 

  721.                    ids,ide, jds,jde, kds,kde,                      &
    722. 

  722.                    ims,ime, jms,jme, kms,kme,                      &
    723. 

  723.                    its,ite, jts,jte, kts,kte)
    724. 

  724. !**********************************************************************
    725. 

  725. !   PBL model called the Asymmetric Convective Model, Version 2 (ACM2) 
    726. 

  726. !   -- See top of module for summary and references
    727. 

  727. !
    728. 

  728. !---- REVISION HISTORY:
    729. 

  729. !   AX     3/2005 - developed WRF version based on ACM2 in the MM5 PX LSM
    730. 

  730. !   JP and RG 8/2006 - updates
    731. 

  731. !
    732. 

  732. !**********************************************************************
    733. 

  733. !  ARGUMENTS:
    734. 

  734. !-- DTPBL           PBL time step
    735. 

  735. !-- PSTAR           Psurf - Ptop in cb
    736. 

  736. !-- NOCONV          If free convection =0, no; =1, yes
    737. 

  737. !-- SIGMAF          Sigma for full layer
    738. 

  738. !-- DSIGH           Sigma thickness
    739. 

  739. !-- DSIGHI          Inverse of sigma thickness
    740. 

  740. !-- JX              N-S index
    741. 

  741. !-- KLPBL           PBL level at K index
    742. 

  742. !-- PBL             PBL height in m
    743. 

  743. !-- PBLSIG          Sigma level for PBL 
    744. 

  744. !-- MOL             Monin-Obukhov length in 1D form
    745. 

  745. !-- UST             U* in 1D form
    746. 

  746. !-- TST             Theta* in 1D form
    747. 

  747. !-- QST             Q* in 1D form
    748. 

  748. !-- USTM            U* for computation of momemtum flux 
    749. 

  749. !-- EDDYZ           eddy diffusivity KZ
    750. 

  750. !-- DENSX           dry air density (kg/m^3)
    751. 

  751. !-- US              U wind 
    752. 

  752. !-- VS              V wind
    753. 

  753. !-- THETA           potential temperature
    754. 

  754. !-- QVS             water vapor mixing ratio (Kg/Kg)
    755. 

  755. !-- QCS             cloud mixing ratio (Kg/Kg)
    756. 

  756. !-- QIS             ice mixing ratio (Kg/Kg)
    757. 

  757. !-- UX              new U wind 
    758. 

  758. !-- VX              new V wind
    759. 

  759. !-- THETAX          new potential temperature
    760. 

  760. !-- QVX             new water vapor mixing ratio (Kg/Kg)
    761. 

  761. !-- QCX             new cloud mixing ratio (Kg/Kg)
    762. 

  762. !-- QIX             new ice mixing ratio (Kg/Kg)
    763. 

  763. !-----------------------------------------------------------------------
    764. 

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

  765. 

  766.       IMPLICIT NONE
    767. 

  767. 

  768. !.......Arguments
    769. 

  769. 

  770. !... Integer
    771. 

  771.       INTEGER,  INTENT(IN)      ::      ids,ide, jds,jde, kds,kde, &
    772. 

  772.                                         ims,ime, jms,jme, kms,kme, &
    773. 

  773.                                         its,ite, jts,jte, kts,kte, JX
    774. 

  774.       INTEGER,  DIMENSION( its:ite ), INTENT(IN)  :: NOCONV
    775. 

  775.       INTEGER,  DIMENSION( ims:ime ), INTENT(IN)  :: KLPBL
    776. 

  776. 

  777. !... Real
    778. 

  778.       REAL , DIMENSION( ims:ime ),          INTENT(IN)  :: PBL, UST
    779. 

  779.       REAL ,                                INTENT(IN)  :: DTPBL
    780. 

  780.       REAL , DIMENSION( its:ite ),          INTENT(IN)  :: PSTAR, PBLSIG,  &
    781. 

  781.                                                            MOL, TST, &
    782. 

  782.                                                            QST, USTM
    783. 

  783.       REAL , DIMENSION( kts:kte ),          INTENT(IN)  :: DSIGHI, DSIGH
    784. 

  784.       REAL , DIMENSION( 0:kte ),            INTENT(IN)  :: SIGMAF
    785. 

  785.       REAL , DIMENSION( its:ite, kts:kte ), INTENT(INOUT)  :: EDDYZ
    786. 

  786.       REAL , DIMENSION( ims:ime, kms:kme ), INTENT(IN)  :: US,VS, THETA,   &
    787. 

  787.                                                            QVS, QCS, QIS, DENSX
    788. 

  788.       REAL , DIMENSION( its:ite, kts:kte ), INTENT(OUT) :: UX, VX, THETAX,      &
    789. 

  789.                                                            QVX, QCX, QIX
    790. 

  790. !.......Local variables
    791. 

  791. 

  792. !... Parameters
    793. 

  793.       INTEGER, PARAMETER :: NSP   = 6
    794. 

  794. !
    795. 

  795. !......ACM2 Parameters
    796. 

  796. !     INTEGER, PARAMETER :: IFACM = 0
    797. 

  797. !
    798. 

  798.       REAL,    PARAMETER :: G1000 = 9.8 * 1.0E-3
    799. 

  799.       REAL,    PARAMETER :: XX    = 0.5          ! FACTOR APPLIED TO CONV MIXING TIME STEP
    800. 

  800.       REAL,    PARAMETER :: KARMAN = 0.4
    801. 

  801. 

  802. !... Integer
    803. 

  803.       INTEGER :: ILX, KL, KLM, I, K, NSPX, NLP, NL, JJ, L
    804. 

  804.       INTEGER :: KCBLMX
    805. 

  805.       INTEGER, DIMENSION( its:ite ) :: KCBL
    806. 

  806. 

  807. !... Real
    808. 

  808.       REAL                               :: G1000I, MBMAX, HOVL, MEDDY, MBAR
    809. 

  809.       REAL                               :: EKZ, RZ, FM, WSPD, DTS, DTRAT, F1
    810. 

  810.       REAL, DIMENSION( its:ite )         :: PSTARI, FSACM, DTLIM
    811. 

  811.       REAL, DIMENSION( kts:kte, its:ite) :: MBARKS, MDWN
    812. 

  812.       REAL, DIMENSION( 1:NSP, its:ite )  :: FS, BCBOTN
    813. 

  813.       REAL, DIMENSION( kts:kte )         :: XPLUS, XMINUS
    814. 

  814.       REAL  DELC
    815. 

  815.       REAL, DIMENSION( 1:NSP,its:ite,kts:kte  ) :: VCI
    816. 

  816. 

  817.       REAL, DIMENSION( kts:kte )               :: AI, BI, CI, EI !, Y
    818. 

  818.       REAL, DIMENSION( 1:NSP,kts:kte )         :: DI, UI    
    819. 

  819. !
    820. 

  820. !--Start Exicutable ----
    821. 

  821. 

  822.       ILX = ite
    823. 

  823.       KL  = kte
    824. 

  824.       KLM = kte - 1
    825. 

  825. 

  826.       G1000I = 1.0 / G1000
    827. 

  827.       KCBLMX = 0
    828. 

  828.       MBMAX  = 0.0
    829. 

  829. 

  830. !---COMPUTE ACM MIXING RATE
    831. 

  831.       DO I = its, ILX
    832. 

  832.         DTLIM(I)  = DTPBL
    833. 

  833.         PSTARI(I) = 1.0 / PSTAR(I)
    834. 

  834.         KCBL(I)   = 1
    835. 

  835.         FSACM(I)  = 0.0
    836. 

  836. 

  837.         IF (NOCONV(I) .EQ. 1) THEN
    838. 

  838.           KCBL(I) = KLPBL(I)
    839. 

  839. 

  840. !-------MBARKS IS UPWARD MIXING RATE; MDWN IS DOWNWARD MIXING RATE
    841. 

  841. !--New couple ACM & EDDY-------------------------------------------------------------
    842. 

  842.           HOVL     = -PBL(I) / MOL(I)
    843. 

  843.           FSACM(I) = 1./(1.+((KARMAN/(HOVL))**0.3333)/(0.72*KARMAN))
    844. 

  844.           MEDDY    = EDDYZ(I,1) / (DTPBL * (PBLSIG(I) - SIGMAF(1)))
    845. 

  845.           MBAR     = MEDDY * FSACM(I)
    846. 

  846.           DO K = kts,KCBL(I)-1
    847. 

  847.             EDDYZ(I,K) = EDDYZ(I,K) * (1.0 - FSACM(I))
    848. 

  848.           ENDDO
    849. 

  849. 

  850.           MBMAX = AMAX1(MBMAX,MBAR)
    851. 

  851.           DO K = kts+1,KCBL(I)
    852. 

  852.             MBARKS(K,I) = MBAR
    853. 

  853.             MDWN(K,I)   = MBAR * (PBLSIG(I) - SIGMAF(K-1)) * DSIGHI(K)
    854. 

  854.           ENDDO
    855. 

  855.           MBARKS(1,I) = MBAR
    856. 

  856.           MBARKS(KCBL(I),I) = MDWN(KCBL(I),I)
    857. 

  857.           MDWN(KCBL(I)+1,I) = 0.0
    858. 

  858.         ENDIF
    859. 

  859.       ENDDO                              ! end of I loop
    860. 

  860. 

  861.       DO K = kts,KLM
    862. 

  862.         DO I = its,ILX
    863. 

  863.           EKZ   = EDDYZ(I,K) / DTPBL * DSIGHI(K)
    864. 

  864.           DTLIM(I) = AMIN1(0.75 / EKZ,DTLIM(I))
    865. 

  865.         ENDDO
    866. 

  866.       ENDDO
    867. 

  867.        
    868. 

  868.       DO I = its,ILX 
    869. 

  869.         IF (NOCONV(I) .EQ. 1) THEN
    870. 

  870.           KCBLMX = AMAX0(KLPBL(I),KCBLMX)
    871. 

  871.           RZ     = (SIGMAF(KCBL(I)) - SIGMAF(1)) * DSIGHI(1)
    872. 

  872.           DTLIM(I)  = AMIN1(XX / (MBARKS(1,I) * RZ),DTLIM(I))
    873. 

  873.         ENDIF
    874. 

  874.       ENDDO
    875. 

  875. 

  876.       DO K = kts,KL
    877. 

  877.         DO I = its,ILX
    878. 

  878.           VCI(1,I,K) = THETA(I,K)
    879. 

  879.           VCI(2,I,K) = QVS(I,K)
    880. 

  880.           VCI(3,I,K) = US(I,K)
    881. 

  881.           VCI(4,I,K) = VS(I,K)
    882. 

  882.           ! -- Also mix cloud water and ice IF necessary
    883. 

  883.           ! IF (IMOISTX.NE.1.AND.IMOISTX.NE.3) THEN  !!! Check other PBL models
    884. 

  884.           VCI(5,I,K) = QCS(I,K)
    885. 

  885.           VCI(6,I,K) = QIS(I,K)
    886. 

  886.         ENDDO
    887. 

  887.       ENDDO
    888. 

  888. 

  889.       NSPX=6
    890. 

  890. 

  891.       DO I = its,ILX
    892. 

  892.         FS(1,I) = -UST(I) * TST(I) * DENSX(I,1) * PSTARI(I)
    893. 

  893.         FS(2,I) = -UST(I) * QST(I) * DENSX(I,1) * PSTARI(I)
    894. 

  894.         FM      = -USTM(I) * USTM(I) * DENSX(I,1) * PSTARI(I)
    895. 

  895.         WSPD    = SQRT(US(I,1) * US(I,1) + VS(I,1) * VS(I,1)) + 1.E-9
    896. 

  896.         FS(3,I) = FM * US(I,1) / WSPD
    897. 

  897.         FS(4,I) = FM * VS(I,1) / WSPD
    898. 

  898.         FS(5,I) = 0.0
    899. 

  899.         FS(6,I) = 0.0                      ! SURFACE FLUXES OF CLOUD WATER AND ICE = 0
    900. 

  900.       ENDDO
    901. 

  901. !
    902. 

  902. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    903. 

  903.       DO I = its,ILX      
    904. 

  904. 

  905.         NLP   = INT(DTPBL / DTLIM(I) + 1.0)
    906. 

  906.         DTS   = (DTPBL / NLP)
    907. 

  907.         DTRAT = DTS / DTPBL
    908. 

  908.         DO NL = 1,NLP           ! LOOP OVER SUB TIME LOOP              
    909. 

  909. 

  910. !-- COMPUTE ARRAY ELEMENTS THAT ARE INDEPENDANT OF SPECIES
    911. 

  911. 

  912.           DO K = kts,KL
    913. 

  913.             AI(K) = 0.0
    914. 

  914.             BI(K) = 0.0
    915. 

  915.             CI(K) = 0.0
    916. 

  916.             EI(K) = 0.0
    917. 

  917.           ENDDO
    918. 

  918. 

  919.           DO K = 2, KCBL(I)
    920. 

  920.             EI(K-1) = -CRANKP * MDWN(K,I) * DTS * DSIGH(K) * DSIGHI(K-1)
    921. 

  921.             BI(K)   = 1.0 + CRANKP * MDWN(K,I) * DTS
    922. 

  922.             AI(K)   = -CRANKP * MBARKS(K,I) * DTS
    923. 

  923.           ENDDO
    924. 

  924. 

  925.           EI(1) = EI(1) -EDDYZ(I,1) * CRANKP * DSIGHI(1 )* DTRAT
    926. 

  926.           AI(2) = AI(2) -EDDYZ(I,1) * CRANKP * DSIGHI(2) * DTRAT
    927. 

  927. 

  928.           DO K =  KCBL(I)+1, KL
    929. 

  929.             BI(K) = 1.0
    930. 

  930.           ENDDO
    931. 

  931. 

  932.           DO K = 2,KL
    933. 

  933.             XPLUS(K)  = EDDYZ(I,K) * DSIGHI(K) * DTRAT
    934. 

  934.             XMINUS(K) = EDDYZ(I,K-1) * DSIGHI(K) * DTRAT
    935. 

  935.             CI(K)     = - XMINUS(K) * CRANKP
    936. 

  936.             EI(K)     = EI(K) - XPLUS(K) * CRANKP
    937. 

  937.             BI(K)     = BI(K) + XPLUS(K) * CRANKP + XMINUS(K) * CRANKP
    938. 

  938.           ENDDO
    939. 

  939. 

  940.           IF (NOCONV(I) .EQ. 1) THEN
    941. 

  941.             BI(1) = 1.0 + CRANKP * MBARKS(1,I) * (PBLSIG(I) - SIGMAF(1)) *    &
    942. 

  942.                     DTS * DSIGHI(1) + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
    943. 

  943.           ELSE
    944. 

  944.             BI(1) = 1.0  + EDDYZ(I,1) * DSIGHI(1) * CRANKP * DTRAT
    945. 

  945.           ENDIF
    946. 

  946. 

  947. 

  948.           DO K = 1,KL
    949. 

  949.             DO L = 1,NSPX                    
    950. 

  950.               DI(L,K) = 0.0
    951. 

  951.             ENDDO
    952. 

  952.           ENDDO
    953. 

  953. !
    954. 

  954. !**   COMPUTE TENDENCY OF CBL CONCENTRATIONS - SEMI-IMPLICIT SOLUTION
    955. 

  955.           DO K = 2,KCBL(I)
    956. 

  956.             DO L = 1,NSPX                    
    957. 

  957.               DELC = DTS * (MBARKS(K,I) * VCI(L,I,1) - MDWN(K,I) *          &
    958. 

  958.                  VCI(L,I,K) + DSIGH(K+1) * DSIGHI(K) *                  &
    959. 

  959.                         MDWN(K+1,I) * VCI(L,I,K+1))
    960. 

  960.               DI(L,K)   = VCI(L,I,K) + (1.0 - CRANKP) * DELC
    961. 

  961.             ENDDO
    962. 

  962.           ENDDO
    963. 

  963. 

  964.           DO K = KCBL(I)+1, KL
    965. 

  965.             DO L = 1,NSPX                    
    966. 

  966.               DI(L,K) = VCI(L,I,K)
    967. 

  967.             ENDDO
    968. 

  968.           ENDDO
    969. 

  969. 

  970.           DO K = 2,KL
    971. 

  971.             IF (K .EQ. KL) THEN
    972. 

  972.               DO L = 1,NSPX                    
    973. 

  973.                 DI(L,K) = DI(L,K)  - (1.0 - CRANKP) * XMINUS(K) *                  &
    974. 

  974.                           (VCI(L,I,K) - VCI(L,I,K-1))
    975. 

  975.               ENDDO
    976. 

  976.             ELSE
    977. 

  977.               DO L = 1,NSPX                    
    978. 

  978.                 DI(L,K) = DI(L,K) + (1.0 - CRANKP) * XPLUS(K) *                   &
    979. 

  979.                           (VCI(L,I,K+1) - VCI(L,I,K))  -                         &
    980. 

  980.                           (1.0 - CRANKP) * XMINUS(K) *                           &
    981. 

  981.                           (VCI(L,I,K) - VCI(L,I,K-1))
    982. 

  982.               ENDDO
    983. 

  983.             ENDIF
    984. 

  984.           ENDDO
    985. 

  985. 

  986.           IF (NOCONV(I) .EQ. 1) THEN
    987. 

  987.             DO L = 1,NSPX                    
    988. 

  988.               F1    = -G1000I * (MBARKS(1,I) *                                &
    989. 

  989.                       (PBLSIG(I) - SIGMAF(1)) * VCI(L,I,1) -                  &
    990. 

  990.                       MDWN(2,I) * VCI(L,I,2) * DSIGH(2))
    991. 

  991. 

  992.               DI(L,1) = VCI(L,I,1) - G1000 * (FS(L,I) - (1.0 - CRANKP)        &
    993. 

  993.                         * F1) * DSIGHI(1) * DTS
    994. 

  994.             ENDDO
    995. 

  995.           ELSE
    996. 

  996.             DO L = 1,NSPX                    
    997. 

  997.               DI(L,1) = VCI(L,I,1) - G1000 * FS(L,I) * DSIGHI(1) * DTS
    998. 

  998.             ENDDO
    999. 

  999.           ENDIF
    1000. 

  1000.           DO L = 1,NSPX                    
    1001. 

  1001.             DI(L,1) = DI(L,1) + (1.0 - CRANKP) * EDDYZ(I,1) * DSIGHI(1)      &
    1002. 

  1002.                      * DTRAT * (VCI(L,I,2) - VCI(L,I,1))
    1003. 

  1003.           ENDDO
    1004. 

  1004.           IF ( NOCONV(I) .EQ. 1 ) THEN
    1005. 

  1005.             CALL MATRIX (AI, BI, CI, DI, EI, UI, KL, NSPX)
    1006. 

  1006.           ELSE
    1007. 

  1007.             CALL TRI (CI, BI, EI, DI, UI, KL, NSPX)
    1008. 

  1008.           END IF
    1009. 

  1009. !
    1010. 

  1010. !-- COMPUTE NEW THETAV AND Q
    1011. 

  1011.           DO K = 1,KL
    1012. 

  1012.             DO L = 1,NSPX                    
    1013. 

  1013.               VCI(L,I,K) = UI(L,K)
    1014. 

  1014.             ENDDO
    1015. 

  1015.           ENDDO
    1016. 

  1016. 

  1017.         ENDDO                   ! END I LOOP
    1018. 

  1018.       ENDDO                     ! END SUB TIME LOOP
    1019. 

  1019. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1020. 

  1020. 

  1021. !
    1022. 

  1022.       DO K = kts,KL
    1023. 

  1023.         DO I = its,ILX
    1024. 

  1024.           THETAX(I,K) = VCI(1,I,K)
    1025. 

  1025.           QVX(I,K)    = VCI(2,I,K)
    1026. 

  1026.           UX(I,K)     = VCI(3,I,K)
    1027. 

  1027.           VX(I,K)     = VCI(4,I,K)
    1028. 

  1028.         ENDDO
    1029. 

  1029.       ENDDO
    1030. 

  1030. 

  1031.       DO K = kts,KL
    1032. 

  1032.         DO I = its,ILX
    1033. 

  1033.           QCX(I,K) = VCI(5,I,K)
    1034. 

  1034.           QIX(I,K) = VCI(6,I,K)
    1035. 

  1035.         ENDDO
    1036. 

  1036.       ENDDO
    1037. 

  1037. 

  1038.    END SUBROUTINE ACM
    1039. 

  1039. !-----------------------------------------------------------------------
    1040. 

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

  1041. 

  1042. !-----------------------------------------------------------------------
    1043. 

  1043. !-----------------------------------------------------------------------
    1044. 

  1044.    SUBROUTINE MATRIX(A,B,C,D,E,X,KL,NSP)
    1045. 

  1045.    
    1046. 

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

  1047. !-----------------------------------------------------------------------
    1048. 

  1048.    IMPLICIT NONE
    1049. 

  1049. !
    1050. 

  1050. !-- Bordered band diagonal matrix solver for ACM2
    1051. 

  1051. 

  1052. !-- ACM2 Matrix is in this form:
    1053. 

  1053. !   B1 E1
    1054. 

  1054. !   A2 B2 E2
    1055. 

  1055. !   A3 C3 B3 E3
    1056. 

  1056. !   A4    C4 B4 E4
    1057. 

  1057. !   A5       C5 B5 E5
    1058. 

  1058. !   A6          C6 B6
    1059. 

  1059. 

  1060. !--Upper Matrix is
    1061. 

  1061. !  U11 U12
    1062. 

  1062. !      U22 U23
    1063. 

  1063. !          U33 U34
    1064. 

  1064. !              U44 U45
    1065. 

  1065. !                  U55 U56
    1066. 

  1066. !                      U66
    1067. 

  1067. 

  1068. !--Lower Matrix is:
    1069. 

  1069. !  1
    1070. 

  1070. ! L21  1
    1071. 

  1071. ! L31 L32  1
    1072. 

  1072. ! L41 L42 L43  1
    1073. 

  1073. ! L51 L52 L53 L54  1
    1074. 

  1074. ! L61 L62 L63 L64 L65 1
    1075. 

  1075. !---------------------------------------------------------
    1076. 

  1076. !...Arguments
    1077. 

  1077.       INTEGER, INTENT(IN)   :: KL
    1078. 

  1078.       INTEGER, INTENT(IN)   :: NSP
    1079. 

  1079.       REAL A(KL),B(KL),E(KL)
    1080. 

  1080.       REAL C(KL),D(NSP,KL),X(NSP,KL)
    1081. 

  1081. 

  1082. !...Locals
    1083. 

  1083.       REAL Y(NSP,KL),AIJ,SUM
    1084. 

  1084.       REAL L(KL,KL),UII(KL),UIIP1(KL),RUII(KL)
    1085. 

  1085.       INTEGER I,J,V
    1086. 

  1086. 

  1087. !-- Define Upper and Lower matrices
    1088. 

  1088.       L(1,1) = 1.
    1089. 

  1089.       UII(1) = B(1)
    1090. 

  1090.       RUII(1) = 1./UII(1)
    1091. 

  1091.       DO I = 2, KL
    1092. 

  1092. 	      L(I,I) = 1.
    1093. 

  1093. 	      L(I,1) = A(I)/B(1)
    1094. 

  1094.         UIIP1(I-1)=E(I-1)
    1095. 

  1095. 	      IF(I.GE.3) THEN
    1096. 

  1096. 	        DO J = 2,I-1
    1097. 

  1097. 	          IF(I.EQ.J+1) THEN
    1098. 

  1098. 	            AIJ = C(I)
    1099. 

  1099. 	          ELSE
    1100. 

  1100. 	            AIJ = 0.
    1101. 

  1101. 	          ENDIF
    1102. 

  1102. 	          L(I,J) = (AIJ-L(I,J-1)*E(J-1))/      &
    1103. 

  1103.                       (B(J)-L(J,J-1)*E(J-1))
    1104. 

  1104. 	        ENDDO
    1105. 

  1105. 	      ENDIF
    1106. 

  1106.       ENDDO
    1107. 

  1107. 	  
    1108. 

  1108.       DO I = 2,KL
    1109. 

  1109.         UII(I) = B(I)-L(I,I-1)*E(I-1)
    1110. 

  1110.         RUII(I) = 1./UII(I)
    1111. 

  1111.       ENDDO
    1112. 

  1112.   
    1113. 

  1113. !-- Forward sub for Ly=d
    1114. 

  1114.       DO V= 1, NSP
    1115. 

  1115.         Y(V,1) = D(V,1)
    1116. 

  1116.         DO I=2,KL
    1117. 

  1117. 	        SUM = D(V,I)
    1118. 

  1118. 	        DO J=1,I-1
    1119. 

  1119. 	          SUM = SUM - L(I,J)*Y(V,J)
    1120. 

  1120. 	        ENDDO
    1121. 

  1121. 	        Y(V,I) = SUM
    1122. 

  1122.         ENDDO
    1123. 

  1123.       ENDDO
    1124. 

  1124. 

  1125. !-- Back sub for Ux=y
    1126. 

  1126. 

  1127.       DO V= 1, NSP
    1128. 

  1128.         X(V,KL) = Y(V,KL)*RUII(KL)
    1129. 

  1129.       ENDDO
    1130. 

  1130.       DO I = KL-1,1,-1
    1131. 

  1131.         DO V= 1, NSP
    1132. 

  1132.          X(V,I) = (Y(V,I)-UIIP1(I)*X(V,I+1))*RUII(I)
    1133. 

  1133.         ENDDO
    1134. 

  1134.       ENDDO
    1135. 

  1135. 

  1136.    END SUBROUTINE MATRIX
    1137. 

  1137. 

  1138. 

  1139. !-----------------------------------------------------------------------
    1140. 

  1140. !-----------------------------------------------------------------------
    1141. 

  1141.       SUBROUTINE TRI ( L, D, U, B, X,KL,NSP)
    1142. 

  1142. !-----------------------------------------------------------------------
    1143. 

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

  1144. 

  1145. !  FUNCTION:
    1146. 

  1146. !    Solves tridiagonal system by Thomas algorithm. 
    1147. 

  1147. !   The associated tri-diagonal system is stored in 3 arrays
    1148. 

  1148. !   D : diagonal
    1149. 

  1149. !   L : sub-diagonal
    1150. 

  1150. !   U : super-diagonal
    1151. 

  1151. !   B : right hand side function
    1152. 

  1152. !   X : return solution from tridiagonal solver
    1153. 

  1153. 

  1154. !     [ D(1) U(1) 0    0    0 ...       0     ]
    1155. 

  1155. !     [ L(2) D(2) U(2) 0    0 ...       .     ]
    1156. 

  1156. !     [ 0    L(3) D(3) U(3) 0 ...       .     ]
    1157. 

  1157. !     [ .       .     .     .           .     ] X(i) = B(i)
    1158. 

  1158. !     [ .             .     .     .     0     ]
    1159. 

  1159. !     [ .                   .     .     .     ]
    1160. 

  1160. !     [ 0                           L(n) D(n) ]
    1161. 

  1161. 

  1162. !-----------------------------------------------------------------------
    1163. 

  1163. 

  1164.       IMPLICIT NONE
    1165. 

  1165. 

  1166. ! Arguments:
    1167. 

  1167. 

  1168.       INTEGER, INTENT(IN)   :: KL
    1169. 

  1169.       INTEGER, INTENT(IN)   :: NSP
    1170. 

  1170. 

  1171.       REAL        L( KL )               ! subdiagonal
    1172. 

  1172.       REAL        D(KL)   ! diagonal
    1173. 

  1173.       REAL        U( KL )               ! superdiagonal
    1174. 

  1174.       REAL        B(NSP,KL )   ! R.H. side
    1175. 

  1175.       REAL        X( NSP,KL )   ! solution
    1176. 

  1176. 

  1177. ! Local Variables:
    1178. 

  1178. 

  1179.       REAL        GAM( KL )
    1180. 

  1180.       REAL        BET
    1181. 

  1181.       INTEGER     V, K
    1182. 

  1182. 

  1183. ! Decomposition and forward substitution:
    1184. 

  1184.       BET = 1.0 / D( 1 )
    1185. 

  1185.       DO V = 1, NSP
    1186. 

  1186.          X( V,1 ) = BET * B(V,1 )
    1187. 

  1187.       ENDDO
    1188. 

  1188. 

  1189.       DO K = 2, KL
    1190. 

  1190.         GAM(K ) = BET * U( K-1 )
    1191. 

  1191.         BET = 1.0 / ( D( K ) - L( K ) * GAM( K ) )
    1192. 

  1192. 	      DO V = 1, NSP
    1193. 

  1193.            X( V, K ) = BET * ( B( V,K ) - L( K ) * X( V,K-1 ) )
    1194. 

  1194. 	      ENDDO
    1195. 

  1195.       ENDDO
    1196. 

  1196. 

  1197. ! Back-substitution:
    1198. 

  1198. 

  1199.       DO K = KL - 1, 1, -1
    1200. 

  1200.         DO V = 1, NSP
    1201. 

  1201.           X( V,K ) = X( V,K ) - GAM( K+1 ) * X( V,K+1 )
    1202. 

  1202.         ENDDO
    1203. 

  1203.       ENDDO
    1204. 

  1204.      
    1205. 

  1205.   END SUBROUTINE TRI
    1206. 

  1206. !-----------------------------------------------------------------------
    1207. 

  1207. !-----------------------------------------------------------------------
    1208. 

  1208. 

  1209. END MODULE module_bl_acm
    1210. 

  1210.                         
    1211. 

  1211.