phys/module_cu_osas.F

   1 !!
   2 MODULE module_cu_osas
   3
   4 CONTAINS
   5
   6 !-----------------------------------------------------------------
   7       SUBROUTINE CU_OSAS(DT,ITIMESTEP,STEPCU,                        &
   8                  RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN,               &
   9                  RUCUTEN,RVCUTEN,                                   & ! gopal's doing for SAS
  10                  RAINCV,PRATEC,HTOP,HBOT,                           &
  11                  U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D,           &
  12                  DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG,                   &
  13                  P_QC,                                              &
  14                  STORE_RAND,MOMMIX, & ! gopal's doing
  15                  P_QI,P_FIRST_SCALAR,                               &
  16                  ids,ide, jds,jde, kds,kde,                         &
  17                  ims,ime, jms,jme, kms,kme,                         &
  18                  its,ite, jts,jte, kts,kte                          )
  19
  20 !-------------------------------------------------------------------
  21       USE MODULE_GFS_MACHINE , ONLY : kind_phys
  22       USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
  23       USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP  &
  24      &,             RV => con_RV, FV => con_fvirt, T0C => con_T0C       &
  25      &,             CVAP => con_CVAP, CLIQ => con_CLIQ                  &
  26      &,             EPS => con_eps, EPSM1 => con_epsm1                  &
  27      &,             ROVCP => con_rocp, RD => con_rd
  28 !-------------------------------------------------------------------
  29       IMPLICIT NONE
  30 !-------------------------------------------------------------------
  31 !-- U3D         3D u-velocity interpolated to theta points (m/s)
  32 !-- V3D         3D v-velocity interpolated to theta points (m/s)
  33 !-- TH3D        3D potential temperature (K)
  34 !-- T3D         temperature (K)
  35 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  36 !-- QC3D        3D cloud mixing ratio (Kg/Kg)
  37 !-- QI3D        3D ice mixing ratio (Kg/Kg)
  38 !-- P8w         3D pressure at full levels (Pa)
  39 !-- Pcps        3D pressure (Pa)
  40 !-- PI3D        3D exner function (dimensionless)
  41 !-- rr3D        3D dry air density (kg/m^3)
  42 !-- RUBLTEN     U tendency due to
  43 !               PBL parameterization (m/s^2)
  44 !-- RVBLTEN     V tendency due to
  45 !               PBL parameterization (m/s^2)
  46 !-- RTHBLTEN    Theta tendency due to
  47 !               PBL parameterization (K/s)
  48 !-- RQVBLTEN    Qv tendency due to
  49 !               PBL parameterization (kg/kg/s)
  50 !-- RQCBLTEN    Qc tendency due to
  51 !               PBL parameterization (kg/kg/s)
  52 !-- RQIBLTEN    Qi tendency due to
  53 !               PBL parameterization (kg/kg/s)
  54 !
  55 !-- MOMMIX      MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
  56 !-- RUCUTEN     U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
  57 !-- RVCUTEN     V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
  58 !
  59 !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
  60 !-- GRAV        acceleration due to gravity (m/s^2)
  61 !-- ROVCP       R/CP
  62 !-- RD          gas constant for dry air (J/kg/K)
  63 !-- ROVG        R/G
  64 !-- P_QI        species index for cloud ice
  65 !-- dz8w        dz between full levels (m)
  66 !-- z           height above sea level (m)
  67 !-- PSFC        pressure at the surface (Pa)
  68 !-- UST         u* in similarity theory (m/s)
  69 !-- PBL         PBL height (m)
  70 !-- PSIM        similarity stability function for momentum
  71 !-- PSIH        similarity stability function for heat
  72 !-- HFX         upward heat flux at the surface (W/m^2)
  73 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
  74 !-- TSK         surface temperature (K)
  75 !-- GZ1OZ0      log(z/z0) where z0 is roughness length
  76 !-- WSPD        wind speed at lowest model level (m/s)
  77 !-- BR          bulk Richardson number in surface layer
  78 !-- DT          time step (s)
  79 !-- rvovrd      R_v divided by R_d (dimensionless)
  80 !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  81 !-- KARMAN      Von Karman constant
  82 !-- ids         start index for i in domain
  83 !-- ide         end index for i in domain
  84 !-- jds         start index for j in domain
  85 !-- jde         end index for j in domain
  86 !-- kds         start index for k in domain
  87 !-- kde         end index for k in domain
  88 !-- ims         start index for i in memory
  89 !-- ime         end index for i in memory
  90 !-- jms         start index for j in memory
  91 !-- jme         end index for j in memory
  92 !-- kms         start index for k in memory
  93 !-- kme         end index for k in memory
  94 !-- its         start index for i in tile
  95 !-- ite         end index for i in tile
  96 !-- jts         start index for j in tile
  97 !-- jte         end index for j in tile
  98 !-- kts         start index for k in tile
  99 !-- kte         end index for k in tile
 100 !-------------------------------------------------------------------
 101
 102       INTEGER ::                        ICLDCK
 103
 104       INTEGER, INTENT(IN) ::            ids,ide, jds,jde, kds,kde,      &
 105                                         ims,ime, jms,jme, kms,kme,      &
 106                                         its,ite, jts,jte, kts,kte,      &
 107                                         ITIMESTEP,                      &     !NSTD
 108                                         P_FIRST_SCALAR,                 &
 109                                         P_QC,                           &
 110                                         P_QI,                           &
 111                                         STEPCU
 112
 113       REAL,    INTENT(IN) ::                                            &
 114                                         DT
 115
 116
 117       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::      &
 118                                         RQCCUTEN,                       &
 119                                         RQICUTEN,                       &
 120                                         RQVCUTEN,                       &
 121                                         RTHCUTEN
 122       REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) ::      &
 123                                         RUCUTEN,                        &  ! gopal's doing for SAS
 124                                         RVCUTEN                            ! gopal's doing for SAS
 125       REAL, OPTIONAL,   INTENT(IN) ::    MOMMIX
 126       REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL,                   &
 127                          INTENT(IN) :: STORE_RAND
 128
 129       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
 130                                         XLAND
 131
 132       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::            &
 133                                         RAINCV, PRATEC
 134
 135       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::              &
 136                                         HBOT,                           &
 137                                         HTOP
 138
 139       LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) ::             &
 140                                         CU_ACT_FLAG
 141
 142
 143       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) ::      &
 144                                         DZ8W,                           &
 145                                         P8w,                            &
 146                                         Pcps,                           &
 147                                         PI3D,                           &
 148                                         QC3D,                           &
 149                                         QI3D,                           &
 150                                         QV3D,                           &
 151                                         RHO3D,                          &
 152                                         T3D,                            &
 153                                         U3D,                            &
 154                                         V3D,                            &
 155                                         W
 156
 157 !--------------------------- LOCAL VARS ------------------------------
 158
 159       REAL,    DIMENSION(ims:ime, jms:jme) ::                           &
 160                                         PSFC
 161
 162
 163       REAL     (kind=kind_phys) ::                                      &
 164                                         DELT,                           &
 165                                         DPSHC,                          &
 166                                         RDELT,                          &
 167                                         RSEED
 168
 169       REAL     (kind=kind_phys), DIMENSION(its:ite) ::                  &
 170                                         CLDWRK,                         &
 171                                         PS,                             &
 172                                         RCS,                            &
 173                                         RN,                             &
 174                                         SLIMSK,                         &
 175                                         XKT2
 176
 177       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) ::       &
 178                                         PRSI
 179
 180       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte) ::         &
 181                                         DEL,                            &
 182                                         DOT,                            &
 183                                         PHIL,                           &
 184                                         PRSL,                           &
 185                                         PRSLK,                          &
 186                                         Q1,                             &
 187                                         T1,                             &
 188                                         U1,                             &
 189                                         V1,                             &
 190                                         ZI,                             &
 191                                         ZL
 192
 193       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) ::      &
 194                                         QL
 195
 196       INTEGER, DIMENSION(its:ite) ::                                    &
 197                                         KBOT,                           &
 198                                         KTOP,                           &
 199                                         KUO
 200
 201       INTEGER ::                                                        &
 202                                         I,                              &
 203                                         IGPVS,                          &
 204                                         IM,                             &
 205                                         J,                              &
 206                                         JCAP,                           &
 207                                         K,                              &
 208                                         KM,                             &
 209                                         KP,                             &
 210                                         KX,                             &
 211                                         NCLOUD
 212
 213       DATA IGPVS/0/
 214
 215 !-----------------------------------------------------------------------
 216 !
 217
 218       DO J=JTS,JTE
 219          DO I=ITS,ITE
 220             CU_ACT_FLAG(I,J)=.TRUE.
 221          ENDDO
 222       ENDDO
 223
 224       IM=ITE-ITS+1
 225       KX=KTE-KTS+1
 226       JCAP=126
 227       DPSHC=30_kind_phys
 228       DELT=DT*STEPCU
 229       RDELT=1./DELT
 230       NCLOUD=1
 231
 232
 233    DO J=jms,jme
 234      DO I=ims,ime
 235        PSFC(i,j)=P8w(i,kms,j)
 236      ENDDO
 237    ENDDO
 238
 239    if(igpvs.eq.0) CALL GFUNCPHYS
 240    igpvs=1
 241
 242 !-------------  J LOOP (OUTER) --------------------------------------------------
 243
 244    DO J=jts,jte
 245
 246 ! --------------- compute zi and zl -----------------------------------------
 247       DO i=its,ite
 248         ZI(I,KTS)=0.0
 249       ENDDO
 250
 251       DO k=kts+1,kte
 252         KM=K-1
 253         DO i=its,ite
 254           ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
 255         ENDDO
 256       ENDDO
 257
 258       DO k=kts+1,kte
 259         KM=K-1
 260         DO i=its,ite
 261           ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
 262         ENDDO
 263       ENDDO
 264
 265       DO i=its,ite
 266         ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
 267       ENDDO
 268
 269 ! --------------- end compute zi and zl -------------------------------------
 270
 271 !    Based on some important findings from Morris Bender, XKT2 was defined in
 272 !    terms of random number instead of random number based cloud tops
 273 !    Also, these random numbers are stored and are changed only once in
 274 !    approximately 5 minutes interval now. This is gopal's doing for HWRF.
 275
 276 !     call random_number(XKT2)
 277
 278 #if (EM_CORE == 1)
 279 !    XKT2 was defined in terms of random number instead of random number based cloud tops
 280 !    ZCX
 281      call init_random_seed()
 282      call random_number(XKT2)
 283 #ifdef REGTEST
 284 ! for regtest only
 285      xkt2 = 0.1
 286 #endif
 287 #endif
 288 !
 289 #if (NMM_CORE == 1)
 290       DO i=its,ite
 291          XKT2(i) = STORE_RAND(i,j)
 292       ENDDO
 293 #endif
 294
 295       DO i=its,ite
 296         PS(i)=PSFC(i,j)*.001
 297         RCS(i)=1.
 298         SLIMSK(i)=ABS(XLAND(i,j)-2.)
 299       ENDDO
 300
 301       DO i=its,ite
 302         PRSI(i,kts)=PS(i)
 303       ENDDO
 304
 305       DO k=kts,kte
 306         kp=k+1
 307         DO i=its,ite
 308           PRSL(I,K)=Pcps(i,k,j)*.001
 309           PHIL(I,K)=ZL(I,K)*GRAV
 310           DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
 311         ENDDO
 312       ENDDO
 313
 314       DO k=kts,kte
 315         DO i=its,ite
 316           DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
 317           U1(i,k)=U3D(i,k,j)
 318           V1(i,k)=V3D(i,k,j)
 319           Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
 320           T1(i,k)=T3D(i,k,j)
 321           QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
 322           QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
 323           PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
 324         ENDDO
 325       ENDDO
 326
 327       DO k=kts+1,kte+1
 328         km=k-1
 329         DO i=its,ite
 330           PRSI(i,k)=PRSI(i,km)-del(i,km)
 331         ENDDO
 332       ENDDO
 333
 334
 335       CALL OSASCNV(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL,                  &
 336                   QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,                    &
 337                   KTOP,KUO,SLIMSK,DOT,XKT2,NCLOUD)
 338
 339 !!!   make more like GFDL ... eliminate shallow convection.....
 340 !!!   CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
 341 #if (EM_CORE == 1)
 342       CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
 343 #endif
 344
 345       DO I=ITS,ITE
 346         RAINCV(I,J)=RN(I)*1000./STEPCU
 347         PRATEC(I,J)=RN(I)*1000./(STEPCU * DT)
 348         HBOT(I,J)=KBOT(I)
 349         HTOP(I,J)=KTOP(I)
 350       ENDDO
 351
 352       DO K=KTS,KTE
 353         DO I=ITS,ITE
 354           RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
 355           RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
 356         ENDDO
 357       ENDDO
 358
 359 !===============================================================================
 360 !     ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
 361 !     MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
 362 !     divergence damping term, a reducion factor for cumulum mixing may be
 363 !     required otherwise storms were too weak.
 364 !===============================================================================
 365 !
 366 #if (NMM_CORE == 1)
 367       DO K=KTS,KTE
 368         DO I=ITS,ITE
 369          RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
 370          RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
 371         ENDDO
 372       ENDDO
 373 #endif
 374
 375
 376       IF(P_QC .ge. P_FIRST_SCALAR)THEN
 377         DO K=KTS,KTE
 378           DO I=ITS,ITE
 379             RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
 380           ENDDO
 381         ENDDO
 382       ENDIF
 383
 384       IF(P_QI .ge. P_FIRST_SCALAR)THEN
 385         DO K=KTS,KTE
 386           DO I=ITS,ITE
 387             RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
 388           ENDDO
 389         ENDDO
 390       ENDIF
 391
 392    ENDDO    ! Outer most J loop
 393
 394    END SUBROUTINE CU_OSAS
 395
 396 !====================================================================
 397    SUBROUTINE osasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN,          &
 398                       RUCUTEN,RVCUTEN,                              &   ! gopal's doing for SAS
 399                       RESTART,P_QC,P_QI,P_FIRST_SCALAR,             &
 400                       allowed_to_read,                              &
 401                       ids, ide, jds, jde, kds, kde,                 &
 402                       ims, ime, jms, jme, kms, kme,                 &
 403                       its, ite, jts, jte, kts, kte                  )
 404 !--------------------------------------------------------------------
 405    IMPLICIT NONE
 406 !--------------------------------------------------------------------
 407    LOGICAL , INTENT(IN)           ::  allowed_to_read,restart
 408    INTEGER , INTENT(IN)           ::  ids, ide, jds, jde, kds, kde, &
 409                                       ims, ime, jms, jme, kms, kme, &
 410                                       its, ite, jts, jte, kts, kte
 411    INTEGER , INTENT(IN)           ::  P_FIRST_SCALAR, P_QI, P_QC
 412
 413    REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::  &
 414                                                               RTHCUTEN, &
 415                                                               RQVCUTEN, &
 416                                                               RQCCUTEN, &
 417                                                               RQICUTEN
 418    REAL,     DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) ::  &
 419                                                               RUCUTEN,  & ! gopal's doing for SAS
 420                                                               RVCUTEN
 421
 422    INTEGER :: i, j, k, itf, jtf, ktf
 423
 424    jtf=min0(jte,jde-1)
 425    ktf=min0(kte,kde-1)
 426    itf=min0(ite,ide-1)
 427
 428 #if ( HWRF == 1 )
 429 !zhang's doing
 430    IF(.not.restart .or. .not.allowed_to_read)THEN
 431 !end of zhang's doing
 432 #else
 433    IF(.not.restart)THEN
 434 #endif
 435      DO j=jts,jtf
 436      DO k=kts,ktf
 437      DO i=its,itf
 438        RTHCUTEN(i,k,j)=0.
 439        RQVCUTEN(i,k,j)=0.
 440        RUCUTEN(i,j,k)=0.   ! gopal's doing for SAS
 441        RVCUTEN(i,j,k)=0.    ! gopal's doing for SAS
 442      ENDDO
 443      ENDDO
 444      ENDDO
 445
 446      IF (P_QC .ge. P_FIRST_SCALAR) THEN
 447         DO j=jts,jtf
 448         DO k=kts,ktf
 449         DO i=its,itf
 450            RQCCUTEN(i,k,j)=0.
 451         ENDDO
 452         ENDDO
 453         ENDDO
 454      ENDIF
 455
 456      IF (P_QI .ge. P_FIRST_SCALAR) THEN
 457         DO j=jts,jtf
 458         DO k=kts,ktf
 459         DO i=its,itf
 460            RQICUTEN(i,k,j)=0.
 461         ENDDO
 462         ENDDO
 463         ENDDO
 464      ENDIF
 465    ENDIF
 466
 467       END SUBROUTINE osasinit
 468
 469 ! ------------------------------------------------------------------------
 470
 471       SUBROUTINE OSASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL,         &
 472      &       Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,            &
 473      &       DOT,XKT2,ncloud)
 474 !  for cloud water version
 475 !     parameter(ncloud=0)
 476 !     SUBROUTINE OSASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
 477 !    &       Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
 478 !    &       DOT,xkt2,ncloud)
 479 !
 480       USE MODULE_GFS_MACHINE , ONLY : kind_phys
 481       USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
 482       USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
 483      &,             RV => con_RV, FV => con_fvirt, T0C => con_T0C       &
 484      &,             CVAP => con_CVAP, CLIQ => con_CLIQ                  &
 485      &,             EPS => con_eps, EPSM1 => con_epsm1
 486
 487       implicit none
 488 !
 489 !     include 'constant.h'
 490 !
 491       integer            IM, IX,  KM, JCAP, ncloud,                     &
 492      &                   KBOT(IM), KTOP(IM), KUO(IM), J
 493       real(kind=kind_phys) DELT
 494       real(kind=kind_phys) PS(IM),      DEL(IX,KM),  PRSL(IX,KM),       &
 495 !     real(kind=kind_phys)              DEL(IX,KM),  PRSL(IX,KM),
 496      &                     QL(IX,KM,2), Q1(IX,KM),   T1(IX,KM),         &
 497      &                     U1(IX,KM),   V1(IX,KM),   RCS(IM),           &
 498      &                     CLDWRK(IM),  RN(IM),      SLIMSK(IM),        &
 499      &                     DOT(IX,KM),  XKT2(IM),    PHIL(IX,KM)
 500 !
 501       integer              I, INDX, jmn, k, knumb, latd, lond, km1
 502 !
 503       real(kind=kind_phys) adw,     alpha,   alphal,  alphas,           &
 504      &                     aup,     beta,    betal,   betas,            &
 505      &                     c0,      cpoel,   dellat,  delta,            &
 506      &                     desdt,   deta,    detad,   dg,               &
 507      &                     dh,      dhh,     dlnsig,  dp,               &
 508      &                     dq,      dqsdp,   dqsdt,   dt,               &
 509      &                     dt2,     dtmax,   dtmin,   dv1,              &
 510      &                     dv1q,    dv2,     dv2q,    dv1u,             &
 511      &                     dv1v,    dv2u,    dv2v,    dv3u,             &
 512      &                     dv3v,    dv3,     dv3q,    dvq1,             &
 513      &                     dz,      dz1,     e1,      edtmax,           &
 514      &                     edtmaxl, edtmaxs, el2orc,  elocp,            &
 515      &                     es,      etah,                               &
 516      &                     evef,    evfact,  evfactl, fact1,            &
 517      &                     fact2,   factor,  fjcap,   fkm,              &
 518      &                     fuv,     g,       gamma,   onemf,            &
 519      &                     onemfu,  pdetrn,  pdpdwn,  pprime,           &
 520      &                     qc,      qlk,     qrch,    qs,               &
 521      &                     rain,    rfact,   shear,   tem1,             &
 522      &                     tem2,    terr,    val,     val1,             &
 523      &                     val2,    w1,      w1l,     w1s,              &
 524      &                     w2,      w2l,     w2s,     w3,               &
 525      &                     w3l,     w3s,     w4,      w4l,              &
 526      &                     w4s,     xdby,    xpw,     xpwd,             &
 527      &                     xqc,     xqrch,   xlambu,  mbdt,             &
 528      &                     tem
 529 !
 530 !
 531       integer              JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM),      &
 532      &                     KT2(IM),  KTCON(IM), LMIN(IM),               &
 533      &                     kbm(IM),  kbmax(IM), kmax(IM)
 534 !
 535       real(kind=kind_phys) AA1(IM),     ACRT(IM),   ACRTFCT(IM),        &
 536      &                     DELHBAR(IM), DELQ(IM),   DELQ2(IM),          &
 537      &                     DELQBAR(IM), DELQEV(IM), DELTBAR(IM),        &
 538      &                     DELTV(IM),   DTCONV(IM), EDT(IM),            &
 539      &                     EDTO(IM),    EDTX(IM),   FLD(IM),            &
 540      &                     HCDO(IM),    HKBO(IM),   HMAX(IM),           &
 541      &                     HMIN(IM),    HSBAR(IM),  UCDO(IM),           &
 542      &                     UKBO(IM),    VCDO(IM),   VKBO(IM),           &
 543      &                     PBCDIF(IM),  PDOT(IM),   PO(IM,KM),          &
 544      &                                  PWAVO(IM),  PWEVO(IM),          &
 545 !    &                     PSFC(IM),    PWAVO(IM),  PWEVO(IM),          &
 546      &                     QCDO(IM),    QCOND(IM),  QEVAP(IM),          &
 547      &                     QKBO(IM),    RNTOT(IM),  VSHEAR(IM),         &
 548      &                     XAA0(IM),    XHCD(IM),   XHKB(IM),           &
 549      &                     XK(IM),      XLAMB(IM),  XLAMD(IM),          &
 550      &                     XMB(IM),     XMBMAX(IM), XPWAV(IM),          &
 551      &                     XPWEV(IM),   XQCD(IM),   XQKB(IM)
 552 !
 553 !  PHYSICAL PARAMETERS
 554       PARAMETER(G=grav)
 555       PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP,                            &
 556      &          EL2ORC=HVAP*HVAP/(RV*CP))
 557       PARAMETER(TERR=0.,C0=.002,DELTA=fv)
 558       PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
 559 !  LOCAL VARIABLES AND ARRAYS
 560       real(kind=kind_phys) PFLD(IM,KM),    TO(IM,KM),     QO(IM,KM),    &
 561      &                     UO(IM,KM),      VO(IM,KM),     QESO(IM,KM)
 562 !  cloud water
 563       real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM),    TVO(IM,KM),   &
 564      &                     DBYO(IM,KM),    ZO(IM,KM),     SUMZ(IM,KM),  &
 565      &                     SUMH(IM,KM),    HEO(IM,KM),    HESO(IM,KM),  &
 566      &                     QRCD(IM,KM),    DELLAH(IM,KM), DELLAQ(IM,KM),&
 567      &                     DELLAU(IM,KM),  DELLAV(IM,KM), HCKO(IM,KM),  &
 568      &                     UCKO(IM,KM),    VCKO(IM,KM),   QCKO(IM,KM),  &
 569      &                     ETA(IM,KM),     ETAU(IM,KM),   ETAD(IM,KM),  &
 570      &                     QRCDO(IM,KM),   PWO(IM,KM),    PWDO(IM,KM),  &
 571      &                     RHBAR(IM),      TX1(IM)
 572 !
 573       LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
 574 !
 575       real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
 576 !     SAVE PCRIT, ACRITT
 577       DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,     &
 578      &           350.,300.,250.,200.,150./
 579       DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,       &
 580      &           .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
 581 !  GDAS DERIVED ACRIT
 582 !     DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688,              &
 583 !    &            .743,.813,.886,.947,1.138,1.377,1.896/
 584 !
 585       real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
 586       parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
 587       parameter (RZERO=0.0,RONE=1.0)
 588 !-----------------------------------------------------------------------
 589 !
 590       km1 = km - 1
 591 !  INITIALIZE ARRAYS
 592 !
 593       DO I=1,IM
 594         RN(I)=0.
 595         KBOT(I)=KM+1
 596         KTOP(I)=0
 597         KUO(I)=0
 598         CNVFLG(I) = .TRUE.
 599         DTCONV(I) = 3600.
 600         CLDWRK(I) = 0.
 601         PDOT(I) = 0.
 602         KT2(I) = 0
 603         QLKO_KTCON(I) = 0.
 604         DELLAL(I) = 0.
 605       ENDDO
 606 !!
 607       DO K = 1, 15
 608         ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
 609       ENDDO
 610       DT2 = DELT
 611 !cmr  dtmin = max(dt2,1200.)
 612       val   =         1200.
 613       dtmin = max(dt2, val )
 614 !cmr  dtmax = max(dt2,3600.)
 615       val   =         3600.
 616       dtmax = max(dt2, val )
 617 !  MODEL TUNABLE PARAMETERS ARE ALL HERE
 618       MBDT    = 10.
 619       EDTMAXl = .3
 620       EDTMAXs = .3
 621       ALPHAl  = .5
 622       ALPHAs  = .5
 623       BETAl   = .15
 624       betas   = .15
 625       BETAl   = .05
 626       betas   = .05
 627 !     change for hurricane model
 628         BETAl = .5
 629         betas = .5
 630 !     EVEF    = 0.07
 631       evfact  = 0.3
 632       evfactl = 0.3
 633 !     change for hurricane model
 634          evfact = 0.6
 635          evfactl = .6
 636 #if ( EM_CORE == 1 )
 637 !  HAWAII TEST - ZCX
 638       ALPHAl  = .5
 639       ALPHAs  = .75
 640       BETAl   = .05
 641       betas   = .05
 642       evfact  = 0.5
 643       evfactl = 0.5
 644 #endif
 645       PDPDWN  = 0.
 646       PDETRN  = 200.
 647       xlambu  = 1.e-4
 648       fjcap   = (float(jcap) / 126.) ** 2
 649 !cmr  fjcap   = max(fjcap,1.)
 650       val     =           1.
 651       fjcap   = max(fjcap,val)
 652       fkm     = (float(km) / 28.) ** 2
 653 !cmr  fkm     = max(fkm,1.)
 654       fkm     = max(fkm,val)
 655       W1l     = -8.E-3
 656       W2l     = -4.E-2
 657       W3l     = -5.E-3
 658       W4l     = -5.E-4
 659       W1s     = -2.E-4
 660       W2s     = -2.E-3
 661       W3s     = -1.E-3
 662       W4s     = -2.E-5
 663 !CCCC IF(IM.EQ.384) THEN
 664         LATD  = 92
 665         lond  = 189
 666 !CCCC ELSEIF(IM.EQ.768) THEN
 667 !CCCC   LATD = 80
 668 !CCCC ELSE
 669 !CCCC   LATD = 0
 670 !CCCC ENDIF
 671 !
 672 !  DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
 673 !  AND THE MAXIMUM THETAE FOR UPDRAFT
 674 !
 675       DO I=1,IM
 676         KBMAX(I) = KM
 677         KBM(I)   = KM
 678         KMAX(I)  = KM
 679         TX1(I)   = 1.0 / PS(I)
 680       ENDDO
 681 !
 682       DO K = 1, KM
 683         DO I=1,IM
 684           IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
 685           IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I)   = K + 1
 686           IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I)  = MIN(KM,K + 1)
 687         ENDDO
 688       ENDDO
 689       DO I=1,IM
 690         KBMAX(I) = MIN(KBMAX(I),KMAX(I))
 691         KBM(I)   = MIN(KBM(I),KMAX(I))
 692       ENDDO
 693 !
 694 !   CONVERT SURFACE PRESSURE TO MB FROM CB
 695 !
 696 !!
 697       DO K = 1, KM
 698         DO I=1,IM
 699           if (K .le. kmax(i)) then
 700             PFLD(I,k) = PRSL(I,K) * 10.0
 701             PWO(I,k)  = 0.
 702             PWDO(I,k) = 0.
 703             TO(I,k)   = T1(I,k)
 704             QO(I,k)   = Q1(I,k)
 705             UO(I,k)   = U1(I,k)
 706             VO(I,k)   = V1(I,k)
 707             DBYO(I,k) = 0.
 708             SUMZ(I,k) = 0.
 709             SUMH(I,k) = 0.
 710           endif
 711         ENDDO
 712       ENDDO
 713
 714 !
 715 !  COLUMN VARIABLES
 716 !  P IS PRESSURE OF THE LAYER (MB)
 717 !  T IS TEMPERATURE AT T-DT (K)..TN
 718 !  Q IS MIXING RATIO AT T-DT (KG/KG)..QN
 719 !  TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
 720 !  QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
 721 !
 722       DO K = 1, KM
 723         DO I=1,IM
 724           if (k .le. kmax(i)) then
 725          !jfe        QESO(I,k) = 10. * FPVS(T1(I,k))
 726          !
 727             QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
 728          !
 729             QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
 730          !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
 731             val1      =             1.E-8
 732             QESO(I,k) = MAX(QESO(I,k), val1)
 733          !cmr        QO(I,k)   = max(QO(I,k),1.e-10)
 734             val2      =           1.e-10
 735             QO(I,k)   = max(QO(I,k), val2 )
 736          !           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
 737             TVO(I,k)  = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
 738           endif
 739         ENDDO
 740       ENDDO
 741
 742 !
 743 !  HYDROSTATIC HEIGHT ASSUME ZERO TERR
 744 !
 745       DO K = 1, KM
 746         DO I=1,IM
 747           ZO(I,k) = PHIL(I,k) / G
 748         ENDDO
 749       ENDDO
 750 !  COMPUTE MOIST STATIC ENERGY
 751       DO K = 1, KM
 752         DO I=1,IM
 753           if (K .le. kmax(i)) then
 754 !           tem       = G * ZO(I,k) + CP * TO(I,k)
 755             tem       = PHIL(I,k) + CP * TO(I,k)
 756             HEO(I,k)  = tem  + HVAP * QO(I,k)
 757             HESO(I,k) = tem  + HVAP * QESO(I,k)
 758 !           HEO(I,k)  = MIN(HEO(I,k),HESO(I,k))
 759           endif
 760         ENDDO
 761       ENDDO
 762 !
 763 !  DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
 764 !  THIS IS THE LEVEL WHERE UPDRAFT STARTS
 765 !
 766       DO I=1,IM
 767         HMAX(I) = HEO(I,1)
 768         KB(I) = 1
 769       ENDDO
 770 !!
 771       DO K = 2, KM
 772         DO I=1,IM
 773           if (k .le. kbm(i)) then
 774             IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
 775               KB(I)   = K
 776               HMAX(I) = HEO(I,k)
 777             ENDIF
 778           endif
 779         ENDDO
 780       ENDDO
 781 !     DO K = 1, KMAX - 1
 782 !         TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
 783 !         QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
 784 !         QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
 785 !         HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
 786 !         HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
 787 !     ENDDO
 788       DO K = 1, KM1
 789         DO I=1,IM
 790           if (k .le. kmax(i)-1) then
 791             DZ      = .5 * (ZO(I,k+1) - ZO(I,k))
 792             DP      = .5 * (PFLD(I,k+1) - PFLD(I,k))
 793 !jfe        ES      = 10. * FPVS(TO(I,k+1))
 794 !
 795             ES      = 0.01 * fpvs(TO(I,K+1))      ! fpvs is in Pa
 796 !
 797             PPRIME  = PFLD(I,k+1) + EPSM1 * ES
 798             QS      = EPS * ES / PPRIME
 799             DQSDP   = - QS / PPRIME
 800             DESDT   = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
 801             DQSDT   = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
 802             GAMMA   = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
 803             DT      = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
 804             DQ      = DQSDT * DT + DQSDP * DP
 805             TO(I,k) = TO(I,k+1) + DT
 806             QO(I,k) = QO(I,k+1) + DQ
 807             PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
 808           endif
 809         ENDDO
 810       ENDDO
 811 !
 812       DO K = 1, KM1
 813         DO I=1,IM
 814           if (k .le. kmax(I)-1) then
 815 !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
 816 !
 817             QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
 818 !
 819             QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
 820 !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
 821             val1      =             1.E-8
 822             QESO(I,k) = MAX(QESO(I,k), val1)
 823 !cmr        QO(I,k)   = max(QO(I,k),1.e-10)
 824             val2      =           1.e-10
 825             QO(I,k)   = max(QO(I,k), val2 )
 826 !           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
 827             HEO(I,k)  = .5 * G * (ZO(I,k) + ZO(I,k+1)) +                &
 828      &                  CP * TO(I,k) + HVAP * QO(I,k)
 829             HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) +                &
 830      &                  CP * TO(I,k) + HVAP * QESO(I,k)
 831             UO(I,k)   = .5 * (UO(I,k) + UO(I,k+1))
 832             VO(I,k)   = .5 * (VO(I,k) + VO(I,k+1))
 833           endif
 834         ENDDO
 835       ENDDO
 836 !     k = kmax
 837 !       HEO(I,k) = HEO(I,k)
 838 !       hesol(k) = HESO(I,k)
 839 !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
 840 !        PRINT *, '   HEO ='
 841 !        PRINT 6001, (HEO(I,K),K=1,KMAX)
 842 !        PRINT *, '   HESO ='
 843 !        PRINT 6001, (HESO(I,K),K=1,KMAX)
 844 !        PRINT *, '   TO ='
 845 !        PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
 846 !        PRINT *, '   QO ='
 847 !        PRINT 6003, (QO(I,K),K=1,KMAX)
 848 !        PRINT *, '   QSO ='
 849 !        PRINT 6003, (QESO(I,K),K=1,KMAX)
 850 !      ENDIF
 851 !
 852 !  LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
 853 !
 854       DO I=1,IM
 855         IF(CNVFLG(I)) THEN
 856           INDX    = KB(I)
 857           HKBO(I) = HEO(I,INDX)
 858           QKBO(I) = QO(I,INDX)
 859           UKBO(I) = UO(I,INDX)
 860           VKBO(I) = VO(I,INDX)
 861         ENDIF
 862         FLG(I)    = CNVFLG(I)
 863         KBCON(I)  = KMAX(I)
 864       ENDDO
 865 !!
 866       DO K = 1, KM
 867         DO I=1,IM
 868           if (k .le. kbmax(i)) then
 869             IF(FLG(I).AND.K.GT.KB(I)) THEN
 870               HSBAR(I)   = HESO(I,k)
 871               IF(HKBO(I).GT.HSBAR(I)) THEN
 872                 FLG(I)   = .FALSE.
 873                 KBCON(I) = K
 874               ENDIF
 875             ENDIF
 876           endif
 877         ENDDO
 878       ENDDO
 879       DO I=1,IM
 880         IF(CNVFLG(I)) THEN
 881           PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
 882           PDOT(I)   = 10.* DOT(I,KBCON(I))
 883           IF(PBCDIF(I).GT.150.)    CNVFLG(I) = .FALSE.
 884           IF(KBCON(I).EQ.KMAX(I))  CNVFLG(I) = .FALSE.
 885         ENDIF
 886       ENDDO
 887 !!
 888       TOTFLG = .TRUE.
 889       DO I=1,IM
 890         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
 891       ENDDO
 892       IF(TOTFLG) RETURN
 893 !  FOUND LFC, CAN DEFINE REST OF VARIABLES
 894  6001 FORMAT(2X,-2P10F12.2)
 895  6002 FORMAT(2X,10F12.2)
 896  6003 FORMAT(2X,3P10F12.2)
 897
 898 !
 899 !  DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
 900 !
 901       DO I = 1, IM
 902         alpha = alphas
 903         if(SLIMSK(I).eq.1.) alpha = alphal
 904         IF(CNVFLG(I)) THEN
 905           IF(KB(I).EQ.1) THEN
 906             DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
 907           ELSE
 908             DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1))               &
 909      &         - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
 910           ENDIF
 911           IF(KBCON(I).NE.KB(I)) THEN
 912 !cmr        XLAMB(I) = -ALOG(ALPHA) / DZ
 913             XLAMB(I) = - LOG(ALPHA) / DZ
 914           ELSE
 915             XLAMB(I) = 0.
 916           ENDIF
 917         ENDIF
 918       ENDDO
 919 !  DETERMINE UPDRAFT MASS FLUX
 920       DO K = 1, KM
 921         DO I = 1, IM
 922           if (k .le. kmax(i) .and. CNVFLG(I)) then
 923             ETA(I,k)  = 1.
 924             ETAU(I,k) = 1.
 925           ENDIF
 926         ENDDO
 927       ENDDO
 928       DO K = KM1, 2, -1
 929         DO I = 1, IM
 930           if (k .le. kbmax(i)) then
 931             IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
 932               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
 933               ETA(I,k)  = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
 934               ETAU(I,k) = ETA(I,k)
 935             ENDIF
 936           endif
 937         ENDDO
 938       ENDDO
 939       DO I = 1, IM
 940         IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
 941           DZ = .5 * (ZO(I,2) - ZO(I,1))
 942           ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
 943           ETAU(I,1) = ETA(I,1)
 944         ENDIF
 945       ENDDO
 946 !
 947 !  WORK UP UPDRAFT CLOUD PROPERTIES
 948 !
 949       DO I = 1, IM
 950         IF(CNVFLG(I)) THEN
 951           INDX         = KB(I)
 952           HCKO(I,INDX) = HKBO(I)
 953           QCKO(I,INDX) = QKBO(I)
 954           UCKO(I,INDX) = UKBO(I)
 955           VCKO(I,INDX) = VKBO(I)
 956           PWAVO(I)     = 0.
 957         ENDIF
 958       ENDDO
 959 !
 960 !  CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
 961 !
 962       DO K = 2, KM1
 963         DO I = 1, IM
 964           if (k .le. kmax(i)-1) then
 965             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
 966               FACTOR = ETA(I,k-1) / ETA(I,k)
 967               ONEMF = 1. - FACTOR
 968               HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *                &
 969      &                    .5 * (HEO(I,k) + HEO(I,k+1))
 970               UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF *                &
 971      &                    .5 * (UO(I,k) + UO(I,k+1))
 972               VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF *                &
 973      &                    .5 * (VO(I,k) + VO(I,k+1))
 974               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
 975             ENDIF
 976             IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
 977               HCKO(I,k) = HCKO(I,k-1)
 978               UCKO(I,k) = UCKO(I,k-1)
 979               VCKO(I,k) = VCKO(I,k-1)
 980               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
 981             ENDIF
 982           endif
 983         ENDDO
 984       ENDDO
 985 !  DETERMINE CLOUD TOP
 986       DO I = 1, IM
 987         FLG(I) = CNVFLG(I)
 988         KTCON(I) = 1
 989       ENDDO
 990 !     DO K = 2, KMAX
 991 !       KK = KMAX - K + 1
 992 !         IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
 993 !           KTCON(I) = KK + 1
 994 !           FLG(I) = .FALSE.
 995 !         ENDIF
 996 !     ENDDO
 997       DO K = 2, KM
 998         DO I = 1, IM
 999           if (k .le. kmax(i)) then
1000             IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
1001               KTCON(I) = K
1002               FLG(I) = .FALSE.
1003             ENDIF
1004           endif
1005         ENDDO
1006       ENDDO
1007       DO I = 1, IM
1008         IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
1009      &  CNVFLG(I) = .FALSE.
1010       ENDDO
1011       TOTFLG = .TRUE.
1012       DO I = 1, IM
1013         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1014       ENDDO
1015       IF(TOTFLG) RETURN
1016 !
1017 !  SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
1018 !
1019       DO I = 1, IM
1020         HMIN(I) = HEO(I,KBCON(I))
1021         LMIN(I) = KBMAX(I)
1022         JMIN(I) = KBMAX(I)
1023       ENDDO
1024       DO I = 1, IM
1025         DO K = KBCON(I), KBMAX(I)
1026           IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
1027             LMIN(I) = K + 1
1028             HMIN(I) = HEO(I,k)
1029           ENDIF
1030         ENDDO
1031       ENDDO
1032 !
1033 !  Make sure that JMIN(I) is within the cloud
1034 !
1035       DO I = 1, IM
1036         IF(CNVFLG(I)) THEN
1037           JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
1038           XMBMAX(I) = .1
1039           JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
1040         ENDIF
1041       ENDDO
1042 !
1043 !  ENTRAINING CLOUD
1044 !
1045       do k = 2, km1
1046         DO I = 1, IM
1047           if (k .le. kmax(i)-1) then
1048             if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
1049               SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
1050               SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))    &
1051      &                  * HEO(I,k)
1052             ENDIF
1053           endif
1054         enddo
1055       enddo
1056 !!
1057       DO I = 1, IM
1058         IF(CNVFLG(I)) THEN
1059 !         call random_number(XKT2)
1060 !         call srand(fhour)
1061 !         XKT2(I) = rand()
1062           KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
1063 !         KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1064 !         KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1065           tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1066           tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1067           if (abs(tem2) .gt. 0.000001) THEN
1068             XLAMB(I) = tem1 / tem2
1069           else
1070             CNVFLG(I) = .false.
1071           ENDIF
1072 !         XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1073 !    &          / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1074           XLAMB(I) = max(XLAMB(I),RZERO)
1075           XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
1076         ENDIF
1077       ENDDO
1078 !!
1079       DO I = 1, IM
1080        DWNFLG(I)  = CNVFLG(I)
1081        DWNFLG2(I) = CNVFLG(I)
1082        IF(CNVFLG(I)) THEN
1083         if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
1084       if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
1085      &  DWNFLG(I) = .false.
1086         do k = JMIN(I), KT2(I)
1087           if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
1088         enddo
1089 !       IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
1090 !    &     DWNFLG(I)=.FALSE.
1091         IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
1092      &     DWNFLG2(I)=.FALSE.
1093        ENDIF
1094       ENDDO
1095 !!
1096       DO K = 2, KM1
1097         DO I = 1, IM
1098           if (k .le. kmax(i)-1) then
1099             IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1100               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
1101 !             ETA(I,k)  = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
1102 !  to simplify matter, we will take the linear approach here
1103 !
1104               ETA(I,k)  = ETA(I,k-1) * (1. + XLAMB(I) * dz)
1105               ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
1106             ENDIF
1107           endif
1108         ENDDO
1109       ENDDO
1110 !!
1111       DO K = 2, KM1
1112         DO I = 1, IM
1113           if (k .le. kmax(i)-1) then
1114 !           IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1115             IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
1116               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
1117               ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
1118             ENDIF
1119           endif
1120         ENDDO
1121       ENDDO
1122 !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1123 !        PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
1124 !        PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
1125 !      ENDIF
1126 !     IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
1127 !       print *, ' xlamb =', xlamb
1128 !       print *, ' eta =', (eta(k),k=1,KT2(I))
1129 !       print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
1130 !       print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
1131 !       print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
1132 !       print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
1133 !     ENDIF
1134       DO I = 1, IM
1135         if(DWNFLG(I)) THEN
1136           KTCON(I) = KT2(I)
1137         ENDIF
1138       ENDDO
1139 !
1140 !  CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
1141 !
1142       DO K = 2, KM1
1143         DO I = 1, IM
1144           if (k .le. kmax(i)-1) then
1145 !jfe
1146             IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1147 !jfe      IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1148               FACTOR    = ETA(I,k-1) / ETA(I,k)
1149               ONEMF     = 1. - FACTOR
1150               fuv       = ETAU(I,k-1) / ETAU(I,k)
1151               onemfu    = 1. - fuv
1152               HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *                &
1153      &                    .5 * (HEO(I,k) + HEO(I,k+1))
1154               UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu *                  &
1155      &                    .5 * (UO(I,k) + UO(I,k+1))
1156               VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu *                  &
1157      &                    .5 * (VO(I,k) + VO(I,k+1))
1158               DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1159             ENDIF
1160           endif
1161         ENDDO
1162       ENDDO
1163 !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1164 !        PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
1165 !        PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
1166 !      ENDIF
1167       DO I = 1, IM
1168         if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I))            &
1169      &     THEN
1170           CNVFLG(I) = .false.
1171           DWNFLG(I) = .false.
1172           DWNFLG2(I) = .false.
1173         ENDIF
1174       ENDDO
1175 !!
1176       TOTFLG = .TRUE.
1177       DO I = 1, IM
1178         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1179       ENDDO
1180       IF(TOTFLG) RETURN
1181 !!
1182 !
1183 !  COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
1184 !
1185       DO I = 1, IM
1186           AA1(I) = 0.
1187           RHBAR(I) = 0.
1188       ENDDO
1189       DO K = 1, KM
1190         DO I = 1, IM
1191           if (k .le. kmax(i)) then
1192             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1193               DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1194               DZ1 = (ZO(I,k) - ZO(I,k-1))
1195               GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1196               QRCH = QESO(I,k)                                          &
1197      &             + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
1198               FACTOR = ETA(I,k-1) / ETA(I,k)
1199               ONEMF = 1. - FACTOR
1200               QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF *                &
1201      &                    .5 * (QO(I,k) + QO(I,k+1))
1202               DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
1203               RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
1204 !
1205 !  BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1206 !
1207               IF(DQ.GT.0.) THEN
1208                 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1209                 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1210                 AA1(I) = AA1(I) - DZ1 * G * QLK
1211                 QC = QLK + QRCH
1212                 PWO(I,k) = ETAH * C0 * DZ * QLK
1213                 QCKO(I,k) = QC
1214                 PWAVO(I) = PWAVO(I) + PWO(I,k)
1215               ENDIF
1216             ENDIF
1217           endif
1218         ENDDO
1219       ENDDO
1220       DO I = 1, IM
1221         RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
1222       ENDDO
1223 !
1224 !  this section is ready for cloud water
1225 !
1226       if(ncloud.gt.0) THEN
1227 !
1228 !  compute liquid and vapor separation at cloud top
1229 !
1230       DO I = 1, IM
1231         k = KTCON(I)
1232         IF(CNVFLG(I)) THEN
1233           GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
1234           QRCH = QESO(I,K)                                              &
1235      &         + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
1236           DQ = QCKO(I,K-1) - QRCH
1237 !
1238 !  CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1239 !
1240           IF(DQ.GT.0.) THEN
1241             QLKO_KTCON(I) = dq
1242             QCKO(I,K-1) = QRCH
1243           ENDIF
1244         ENDIF
1245       ENDDO
1246       ENDIF
1247 !
1248 !  CALCULATE CLOUD WORK FUNCTION AT T+DT
1249 !
1250       DO K = 1, KM
1251         DO I = 1, IM
1252           if (k .le. kmax(i)) then
1253             IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1254               DZ1 = ZO(I,k) - ZO(I,k-1)
1255               GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1256               RFACT =  1. + DELTA * CP * GAMMA                          &
1257      &                 * TO(I,k-1) / HVAP
1258               AA1(I) = AA1(I) +                                         &
1259      &                 DZ1 * (G / (CP * TO(I,k-1)))                     &
1260      &                 * DBYO(I,k-1) / (1. + GAMMA)                     &
1261      &                 * RFACT
1262               val = 0.
1263               AA1(I)=AA1(I)+                                            &
1264      &                 DZ1 * G * DELTA *                                &
1265 !cmr &                 MAX( 0.,(QESO(I,k-1) - QO(I,k-1)))               &
1266      &                 MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1267             ENDIF
1268           endif
1269         ENDDO
1270       ENDDO
1271       DO I = 1, IM
1272         IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I)  = .FALSE.
1273         IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1274         IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I)  = .FALSE.
1275       ENDDO
1276 !!
1277       TOTFLG = .TRUE.
1278       DO I = 1, IM
1279         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1280       ENDDO
1281       IF(TOTFLG) RETURN
1282 !!
1283 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1284 !cccc   PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
1285 !cccc ENDIF
1286 !
1287 !------- DOWNDRAFT CALCULATIONS
1288 !
1289 !
1290 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1291 !
1292       DO I = 1, IM
1293         IF(CNVFLG(I)) THEN
1294           VSHEAR(I) = 0.
1295         ENDIF
1296       ENDDO
1297       DO K = 1, KM
1298         DO I = 1, IM
1299           if (k .le. kmax(i)) then
1300             IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
1301               shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2              &
1302      &                          + (VO(I,k+1)-VO(I,k)) ** 2)
1303               VSHEAR(I) = VSHEAR(I) + SHEAR
1304             ENDIF
1305           endif
1306         ENDDO
1307       ENDDO
1308       DO I = 1, IM
1309         EDT(I) = 0.
1310         IF(CNVFLG(I)) THEN
1311           KNUMB = KTCON(I) - KB(I) + 1
1312           KNUMB = MAX(KNUMB,1)
1313           VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
1314           E1=1.591-.639*VSHEAR(I)                                       &
1315      &       +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
1316           EDT(I)=1.-E1
1317 !cmr      EDT(I) = MIN(EDT(I),.9)
1318           val =         .9
1319           EDT(I) = MIN(EDT(I),val)
1320 !cmr      EDT(I) = MAX(EDT(I),.0)
1321           val =         .0
1322           EDT(I) = MAX(EDT(I),val)
1323           EDTO(I)=EDT(I)
1324           EDTX(I)=EDT(I)
1325         ENDIF
1326       ENDDO
1327 !  DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
1328       DO I = 1, IM
1329         KBDTR(I) = KBCON(I)
1330         beta = betas
1331         if(SLIMSK(I).eq.1.) beta = betal
1332         IF(CNVFLG(I)) THEN
1333           KBDTR(I) = KBCON(I)
1334           KBDTR(I) = MAX(KBDTR(I),1)
1335           XLAMD(I) = 0.
1336           IF(KBDTR(I).GT.1) THEN
1337             DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1)            &
1338      &         - ZO(I,1)
1339             XLAMD(I) =  LOG(BETA) / DZ
1340           ENDIF
1341         ENDIF
1342       ENDDO
1343 !  DETERMINE DOWNDRAFT MASS FLUX
1344       DO K = 1, KM
1345         DO I = 1, IM
1346           IF(k .le. kmax(i)) then
1347             IF(CNVFLG(I)) THEN
1348               ETAD(I,k) = 1.
1349             ENDIF
1350             QRCDO(I,k) = 0.
1351           endif
1352         ENDDO
1353       ENDDO
1354       DO K = KM1, 2, -1
1355         DO I = 1, IM
1356           if (k .le. kbmax(i)) then
1357             IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
1358               DZ        = .5 * (ZO(I,k+1) - ZO(I,k-1))
1359               ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1360             ENDIF
1361           endif
1362         ENDDO
1363       ENDDO
1364       K = 1
1365       DO I = 1, IM
1366         IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
1367           DZ = .5 * (ZO(I,2) - ZO(I,1))
1368           ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1369         ENDIF
1370       ENDDO
1371 !
1372 !--- DOWNDRAFT MOISTURE PROPERTIES
1373 !
1374       DO I = 1, IM
1375         PWEVO(I) = 0.
1376         FLG(I) = CNVFLG(I)
1377       ENDDO
1378       DO I = 1, IM
1379         IF(CNVFLG(I)) THEN
1380           JMN = JMIN(I)
1381           HCDO(I) = HEO(I,JMN)
1382           QCDO(I) = QO(I,JMN)
1383           QRCDO(I,JMN) = QESO(I,JMN)
1384           UCDO(I) = UO(I,JMN)
1385           VCDO(I) = VO(I,JMN)
1386         ENDIF
1387       ENDDO
1388       DO K = KM1, 1, -1
1389         DO I = 1, IM
1390           if (k .le. kmax(i)-1) then
1391             IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
1392               DQ = QESO(I,k)
1393               DT = TO(I,k)
1394               GAMMA      = EL2ORC * DQ / DT**2
1395               DH         = HCDO(I) - HESO(I,k)
1396               QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1397               DETAD      = ETAD(I,k+1) - ETAD(I,k)
1398               PWDO(I,k)  = ETAD(I,k+1) * QCDO(I) -                      &
1399      &                     ETAD(I,k) * QRCDO(I,k)
1400               PWDO(I,k)  = PWDO(I,k) - DETAD *                          &
1401      &                    .5 * (QRCDO(I,k) + QRCDO(I,k+1))
1402               QCDO(I)    = QRCDO(I,k)
1403               PWEVO(I)   = PWEVO(I) + PWDO(I,k)
1404             ENDIF
1405           endif
1406         ENDDO
1407       ENDDO
1408 !     IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
1409 !       PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
1410 !     ENDIF
1411 !
1412 !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
1413 !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
1414 !--- EVAPORATE (PWEV)
1415 !
1416       DO I = 1, IM
1417         edtmax = edtmaxl
1418         if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1419         IF(DWNFLG2(I)) THEN
1420           IF(PWEVO(I).LT.0.) THEN
1421             EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
1422             EDTO(I) = MIN(EDTO(I),EDTMAX)
1423           ELSE
1424             EDTO(I) = 0.
1425           ENDIF
1426         ELSE
1427           EDTO(I) = 0.
1428         ENDIF
1429       ENDDO
1430 !
1431 !
1432 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1433 !
1434 !
1435       DO K = KM1, 1, -1
1436         DO I = 1, IM
1437           if (k .le. kmax(i)-1) then
1438             IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1439               GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1440               DHH=HCDO(I)
1441               DT=TO(I,k+1)
1442               DG=GAMMA
1443               DH=HESO(I,k+1)
1444               DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1445               AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG))   &
1446      &               *(1.+DELTA*CP*DG*DT/HVAP)
1447               val=0.
1448               AA1(I)=AA1(I)+EDTO(I)*                                    &
1449 !cmr &        DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1)))               &
1450      &        DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1451             ENDIF
1452           endif
1453         ENDDO
1454       ENDDO
1455 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1456 !cccc   PRINT *, '  AA1(I) AFTER DWNDRFT =', AA1(I)
1457 !cccc ENDIF
1458       DO I = 1, IM
1459         IF(AA1(I).LE.0.) CNVFLG(I)  = .FALSE.
1460         IF(AA1(I).LE.0.) DWNFLG(I)  = .FALSE.
1461         IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1462       ENDDO
1463 !!
1464       TOTFLG = .TRUE.
1465       DO I = 1, IM
1466         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1467       ENDDO
1468       IF(TOTFLG) RETURN
1469 !!
1470 !
1471 !
1472 !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
1473 !--- WILL DO TO THE ENVIRONMENT?
1474 !
1475       DO K = 1, KM
1476         DO I = 1, IM
1477           IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1478             DELLAH(I,k) = 0.
1479             DELLAQ(I,k) = 0.
1480             DELLAU(I,k) = 0.
1481             DELLAV(I,k) = 0.
1482           ENDIF
1483         ENDDO
1484       ENDDO
1485       DO I = 1, IM
1486         IF(CNVFLG(I)) THEN
1487           DP = 1000. * DEL(I,1)
1488           DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I)                  &
1489      &                - HEO(I,1)) * G / DP
1490           DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I)                  &
1491      &                - QO(I,1)) * G / DP
1492           DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I)                  &
1493      &                - UO(I,1)) * G / DP
1494           DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I)                  &
1495      &                - VO(I,1)) * G / DP
1496         ENDIF
1497       ENDDO
1498 !
1499 !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1500 !
1501       DO K = 2, KM1
1502         DO I = 1, IM
1503           if (k .le. kmax(i)-1) then
1504             IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1505               AUP = 1.
1506               IF(K.LE.KB(I)) AUP = 0.
1507               ADW = 1.
1508               IF(K.GT.JMIN(I)) ADW = 0.
1509               DV1= HEO(I,k)
1510               DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1511               DV3= HEO(I,k-1)
1512               DV1Q= QO(I,k)
1513               DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1514               DV3Q= QO(I,k-1)
1515               DV1U= UO(I,k)
1516               DV2U = .5 * (UO(I,k) + UO(I,k+1))
1517               DV3U= UO(I,k-1)
1518               DV1V= VO(I,k)
1519               DV2V = .5 * (VO(I,k) + VO(I,k+1))
1520               DV3V= VO(I,k-1)
1521               DP = 1000. * DEL(I,K)
1522               DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1523               DETA = ETA(I,k) - ETA(I,k-1)
1524               DETAD = ETAD(I,k) - ETAD(I,k-1)
1525               DELLAH(I,k) = DELLAH(I,k) +                               &
1526      &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1   &
1527      &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3   &
1528      &                    - AUP * DETA * DV2                            &
1529      &                    + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1530               DELLAQ(I,k) = DELLAQ(I,k) +                               &
1531      &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q  &
1532      &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q  &
1533      &                    - AUP * DETA * DV2Q                           &
1534      &       +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1535               DELLAU(I,k) = DELLAU(I,k) +                               &
1536      &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U  &
1537      &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U  &
1538      &                     - AUP * DETA * DV2U                          &
1539      &                    + ADW * EDTO(I) * DETAD * UCDO(I)             &
1540      &                    ) * G / DP
1541               DELLAV(I,k) = DELLAV(I,k) +                               &
1542      &            ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V  &
1543      &        - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V  &
1544      &                     - AUP * DETA * DV2V                          &
1545      &                    + ADW * EDTO(I) * DETAD * VCDO(I)             &
1546      &                    ) * G / DP
1547             ENDIF
1548           endif
1549         ENDDO
1550       ENDDO
1551 !
1552 !------- CLOUD TOP
1553 !
1554       DO I = 1, IM
1555         IF(CNVFLG(I)) THEN
1556           INDX = KTCON(I)
1557           DP = 1000. * DEL(I,INDX)
1558           DV1 = HEO(I,INDX-1)
1559           DELLAH(I,INDX) = ETA(I,INDX-1) *                              &
1560      &                     (HCKO(I,INDX-1) - DV1) * G / DP
1561           DVQ1 = QO(I,INDX-1)
1562           DELLAQ(I,INDX) = ETA(I,INDX-1) *                              &
1563      &                     (QCKO(I,INDX-1) - DVQ1) * G / DP
1564           DV1U = UO(I,INDX-1)
1565           DELLAU(I,INDX) = ETA(I,INDX-1) *                              &
1566      &                     (UCKO(I,INDX-1) - DV1U) * G / DP
1567           DV1V = VO(I,INDX-1)
1568           DELLAV(I,INDX) = ETA(I,INDX-1) *                              &
1569      &                     (VCKO(I,INDX-1) - DV1V) * G / DP
1570 !
1571 !  cloud water
1572 !
1573           DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1574         ENDIF
1575       ENDDO
1576 !
1577 !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1578 !
1579       DO K = 1, KM
1580         DO I = 1, IM
1581           if (k .le. kmax(i)) then
1582             IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1583               QO(I,k) = Q1(I,k)
1584               TO(I,k) = T1(I,k)
1585               UO(I,k) = U1(I,k)
1586               VO(I,k) = V1(I,k)
1587             ENDIF
1588             IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1589               QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1590               DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1591               TO(I,k) = DELLAT * MBDT + T1(I,k)
1592 !cmr          QO(I,k) = max(QO(I,k),1.e-10)
1593               val   =           1.e-10
1594               QO(I,k) = max(QO(I,k), val  )
1595             ENDIF
1596           endif
1597         ENDDO
1598       ENDDO
1599 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1600 !
1601 !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1602 !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1603 !--- WOULD HAVE ON THE STABILITY,
1604 !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1605 !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1606 !--- DESTABILIZATION.
1607 !
1608 !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1609 !
1610       DO K = 1, KM
1611         DO I = 1, IM
1612           IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1613 !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
1614 !
1615             QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
1616 !
1617             QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1618 !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
1619             val       =             1.E-8
1620             QESO(I,k) = MAX(QESO(I,k), val )
1621             TVO(I,k)  = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1622           ENDIF
1623         ENDDO
1624       ENDDO
1625       DO I = 1, IM
1626         IF(CNVFLG(I)) THEN
1627           XAA0(I) = 0.
1628           XPWAV(I) = 0.
1629         ENDIF
1630       ENDDO
1631 !
1632 !  HYDROSTATIC HEIGHT ASSUME ZERO TERR
1633 !
1634 !     DO I = 1, IM
1635 !       IF(CNVFLG(I)) THEN
1636 !         DLNSIG =  LOG(PRSL(I,1)/PS(I))
1637 !         ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1638 !       ENDIF
1639 !     ENDDO
1640 !     DO K = 2, KM
1641 !       DO I = 1, IM
1642 !         IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1643 !           DLNSIG =  LOG(PRSL(I,K) / PRSL(I,K-1))
1644 !           ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1645 !    &             * .5 * (TVO(I,k) + TVO(I,k-1))
1646 !         ENDIF
1647 !       ENDDO
1648 !     ENDDO
1649 !
1650 !--- MOIST STATIC ENERGY
1651 !
1652       DO K = 1, KM1
1653         DO I = 1, IM
1654           IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1655             DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1656             DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1657 !jfe        ES = 10. * FPVS(TO(I,k+1))
1658 !
1659             ES = 0.01 * fpvs(TO(I,K+1))      ! fpvs is in Pa
1660 !
1661             PPRIME = PFLD(I,k+1) + EPSM1 * ES
1662             QS = EPS * ES / PPRIME
1663             DQSDP = - QS / PPRIME
1664             DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1665             DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1666             GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1667             DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1668             DQ = DQSDT * DT + DQSDP * DP
1669             TO(I,k) = TO(I,k+1) + DT
1670             QO(I,k) = QO(I,k+1) + DQ
1671             PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1672           ENDIF
1673         ENDDO
1674       ENDDO
1675       DO K = 1, KM1
1676         DO I = 1, IM
1677           IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1678 !jfe        QESO(I,k) = 10. * FPVS(TO(I,k))
1679 !
1680             QESO(I,k) = 0.01 * fpvs(TO(I,K))      ! fpvs is in Pa
1681 !
1682             QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1683 !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
1684             val1      =             1.E-8
1685             QESO(I,k) = MAX(QESO(I,k), val1)
1686 !cmr        QO(I,k)   = max(QO(I,k),1.e-10)
1687             val2      =           1.e-10
1688             QO(I,k)   = max(QO(I,k), val2 )
1689 !           QO(I,k)   = MIN(QO(I,k),QESO(I,k))
1690             HEO(I,k)   = .5 * G * (ZO(I,k) + ZO(I,k+1)) +               &
1691      &                    CP * TO(I,k) + HVAP * QO(I,k)
1692             HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) +                &
1693      &                  CP * TO(I,k) + HVAP * QESO(I,k)
1694           ENDIF
1695         ENDDO
1696       ENDDO
1697       DO I = 1, IM
1698         k = kmax(i)
1699         IF(CNVFLG(I)) THEN
1700           HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1701           HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1702 !         HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1703         ENDIF
1704       ENDDO
1705       DO I = 1, IM
1706         IF(CNVFLG(I)) THEN
1707           INDX = KB(I)
1708           XHKB(I) = HEO(I,INDX)
1709           XQKB(I) = QO(I,INDX)
1710           HCKO(I,INDX) = XHKB(I)
1711           QCKO(I,INDX) = XQKB(I)
1712         ENDIF
1713       ENDDO
1714 !
1715 !
1716 !**************************** STATIC CONTROL
1717 !
1718 !
1719 !------- MOISTURE AND CLOUD WORK FUNCTIONS
1720 !
1721       DO K = 2, KM1
1722         DO I = 1, IM
1723           if (k .le. kmax(i)-1) then
1724 !           IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1725             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1726               FACTOR = ETA(I,k-1) / ETA(I,k)
1727               ONEMF = 1. - FACTOR
1728               HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF *                &
1729      &                    .5 * (HEO(I,k) + HEO(I,k+1))
1730             ENDIF
1731 !           IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1732 !             HEO(I,k) = HEO(I,k-1)
1733 !           ENDIF
1734           endif
1735         ENDDO
1736       ENDDO
1737       DO K = 2, KM1
1738         DO I = 1, IM
1739           if (k .le. kmax(i)-1) then
1740             IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1741               DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1742               GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1743               XDBY = HCKO(I,k) - HESO(I,k)
1744 !cmr          XDBY = MAX(XDBY,0.)
1745               val  =          0.
1746               XDBY = MAX(XDBY,val)
1747               XQRCH = QESO(I,k)                                         &
1748      &              + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1749               FACTOR = ETA(I,k-1) / ETA(I,k)
1750               ONEMF = 1. - FACTOR
1751               QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF *                &
1752      &                    .5 * (QO(I,k) + QO(I,k+1))
1753               DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1754               IF(DQ.GT.0.) THEN
1755                 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1756                 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1757                 XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1758                 XQC = QLK + XQRCH
1759                 XPW = ETAH * C0 * DZ * QLK
1760                 QCKO(I,k) = XQC
1761                 XPWAV(I) = XPWAV(I) + XPW
1762               ENDIF
1763             ENDIF
1764 !           IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1765             IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1766               DZ1 = ZO(I,k) - ZO(I,k-1)
1767               GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1768               RFACT =  1. + DELTA * CP * GAMMA                          &
1769      &                 * TO(I,k-1) / HVAP
1770               XDBY = HCKO(I,k-1) - HESO(I,k-1)
1771               XAA0(I) = XAA0(I)                                         &
1772      &                + DZ1 * (G / (CP * TO(I,k-1)))                    &
1773      &                * XDBY / (1. + GAMMA)                             &
1774      &                * RFACT
1775               val=0.
1776               XAA0(I)=XAA0(I)+                                          &
1777      &                 DZ1 * G * DELTA *                                &
1778 !cmr &                 MAX( 0.,(QESO(I,k-1) - QO(I,k-1)))               &
1779      &                 MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1780             ENDIF
1781           endif
1782         ENDDO
1783       ENDDO
1784 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1785 !cccc   PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1786 !cccc ENDIF
1787 !
1788 !------- DOWNDRAFT CALCULATIONS
1789 !
1790 !
1791 !--- DOWNDRAFT MOISTURE PROPERTIES
1792 !
1793       DO I = 1, IM
1794         XPWEV(I) = 0.
1795       ENDDO
1796       DO I = 1, IM
1797         IF(DWNFLG2(I)) THEN
1798           JMN = JMIN(I)
1799           XHCD(I) = HEO(I,JMN)
1800           XQCD(I) = QO(I,JMN)
1801           QRCD(I,JMN) = QESO(I,JMN)
1802         ENDIF
1803       ENDDO
1804       DO K = KM1, 1, -1
1805         DO I = 1, IM
1806           if (k .le. kmax(i)-1) then
1807             IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1808               DQ = QESO(I,k)
1809               DT = TO(I,k)
1810               GAMMA    = EL2ORC * DQ / DT**2
1811               DH       = XHCD(I) - HESO(I,k)
1812               QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1813               DETAD    = ETAD(I,k+1) - ETAD(I,k)
1814               XPWD     = ETAD(I,k+1) * QRCD(I,k+1) -                    &
1815      &                   ETAD(I,k) * QRCD(I,k)
1816               XPWD     = XPWD - DETAD *                                 &
1817      &                 .5 * (QRCD(I,k) + QRCD(I,k+1))
1818               XPWEV(I) = XPWEV(I) + XPWD
1819             ENDIF
1820           endif
1821         ENDDO
1822       ENDDO
1823 !
1824       DO I = 1, IM
1825         edtmax = edtmaxl
1826         if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1827         IF(DWNFLG2(I)) THEN
1828           IF(XPWEV(I).GE.0.) THEN
1829             EDTX(I) = 0.
1830           ELSE
1831             EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1832             EDTX(I) = MIN(EDTX(I),EDTMAX)
1833           ENDIF
1834         ELSE
1835           EDTX(I) = 0.
1836         ENDIF
1837       ENDDO
1838 !
1839 !
1840 !
1841 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1842 !
1843 !
1844       DO K = KM1, 1, -1
1845         DO I = 1, IM
1846           if (k .le. kmax(i)-1) then
1847             IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1848               GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1849               DHH=XHCD(I)
1850               DT= TO(I,k+1)
1851               DG= GAMMA
1852               DH= HESO(I,k+1)
1853               DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1854               XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1855      &                *(1.+DELTA*CP*DG*DT/HVAP)
1856               val=0.
1857               XAA0(I)=XAA0(I)+EDTX(I)*                                  &
1858 !cmr &        DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1)))               &
1859      &        DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1860             ENDIF
1861           endif
1862         ENDDO
1863       ENDDO
1864 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1865 !cccc   PRINT *, '  XAA AFTER DWNDRFT =', XAA0(I)
1866 !cccc ENDIF
1867 !
1868 !  CALCULATE CRITICAL CLOUD WORK FUNCTION
1869 !
1870       DO I = 1, IM
1871         ACRT(I) = 0.
1872         IF(CNVFLG(I)) THEN
1873 !       IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
1874           IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
1875             ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I)))                   &
1876      &              /(975.-PCRIT(15))
1877           ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
1878             ACRT(I)=ACRIT(1)
1879           ELSE
1880 !cmr        K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
1881             K =  int((850. - PFLD(I,KTCON(I)))/50.) + 2
1882             K = MIN(K,15)
1883             K = MAX(K,2)
1884             ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))*                     &
1885      &           (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
1886            ENDIF
1887 !        ELSE
1888 !          ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
1889          ENDIF
1890       ENDDO
1891       DO I = 1, IM
1892         ACRTFCT(I) = 1.
1893         IF(CNVFLG(I)) THEN
1894           if(SLIMSK(I).eq.1.) THEN
1895             w1 = w1l
1896             w2 = w2l
1897             w3 = w3l
1898             w4 = w4l
1899           else
1900             w1 = w1s
1901             w2 = w2s
1902             w3 = w3s
1903             w4 = w4s
1904           ENDIF
1905 !C       IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
1906 !         ACRTFCT(I) = PDOT(I) / W3
1907 !
1908 !  modify critical cloud workfunction by cloud base vertical velocity
1909 !
1910           IF(PDOT(I).LE.W4) THEN
1911             ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
1912           ELSEIF(PDOT(I).GE.-W4) THEN
1913             ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
1914           ELSE
1915             ACRTFCT(I) = 0.
1916           ENDIF
1917 !cmr      ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
1918           val1    =             -1.
1919           ACRTFCT(I) = MAX(ACRTFCT(I),val1)
1920 !cmr      ACRTFCT(I) = MIN(ACRTFCT(I),1.)
1921           val2    =             1.
1922           ACRTFCT(I) = MIN(ACRTFCT(I),val2)
1923           ACRTFCT(I) = 1. - ACRTFCT(I)
1924 !
1925 !  modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
1926 !
1927 !         if(RHBAR(I).ge..8) THEN
1928 !           ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
1929 !         ENDIF
1930 !
1931 !  modify adjustment time scale by cloud base vertical velocity
1932 !
1933           DTCONV(I) = DT2 + max((1800. - DT2),RZERO) *                  &
1934      &                (PDOT(I) - W2) / (W1 - W2)
1935 !         DTCONV(I) = MAX(DTCONV(I), DT2)
1936 !         DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
1937           DTCONV(I) = max(DTCONV(I),dtmin)
1938           DTCONV(I) = min(DTCONV(I),dtmax)
1939
1940         ENDIF
1941       ENDDO
1942 !
1943 !--- LARGE SCALE FORCING
1944 !
1945       DO I= 1, IM
1946         FLG(I) = CNVFLG(I)
1947         IF(CNVFLG(I)) THEN
1948 !         F = AA1(I) / DTCONV(I)
1949           FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
1950           IF(FLD(I).LE.0.) FLG(I) = .FALSE.
1951         ENDIF
1952         CNVFLG(I) = FLG(I)
1953         IF(CNVFLG(I)) THEN
1954 !         XAA0(I) = MAX(XAA0(I),0.)
1955           XK(I) = (XAA0(I) - AA1(I)) / MBDT
1956           IF(XK(I).GE.0.) FLG(I) = .FALSE.
1957         ENDIF
1958 !
1959 !--- KERNEL, CLOUD BASE MASS FLUX
1960 !
1961         CNVFLG(I) = FLG(I)
1962         IF(CNVFLG(I)) THEN
1963           XMB(I) = -FLD(I) / XK(I)
1964           XMB(I) = MIN(XMB(I),XMBMAX(I))
1965         ENDIF
1966       ENDDO
1967 !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1968 !        print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
1969 !        PRINT *, '  A1, XA =', AA1(I), XAA0(I)
1970 !        PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
1971 !      ENDIF
1972       TOTFLG = .TRUE.
1973       DO I = 1, IM
1974         TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1975       ENDDO
1976       IF(TOTFLG) RETURN
1977 !
1978 !  restore t0 and QO to t1 and q1 in case convection stops
1979 !
1980       do k = 1, km
1981         DO I = 1, IM
1982           if (k .le. kmax(i)) then
1983             TO(I,k) = T1(I,k)
1984             QO(I,k) = Q1(I,k)
1985 !jfe        QESO(I,k) = 10. * FPVS(T1(I,k))
1986 !
1987             QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
1988 !
1989             QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
1990 !cmr        QESO(I,k) = MAX(QESO(I,k),1.E-8)
1991             val     =             1.E-8
1992             QESO(I,k) = MAX(QESO(I,k), val )
1993           endif
1994         enddo
1995       enddo
1996 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1997 !
1998 !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
1999 !---           MULTIPLIED BY  THE MASS FLUX NECESSARY TO KEEP THE
2000 !---           EQUILIBRIUM WITH THE LARGER-SCALE.
2001 !
2002       DO I = 1, IM
2003         DELHBAR(I) = 0.
2004         DELQBAR(I) = 0.
2005         DELTBAR(I) = 0.
2006         QCOND(I) = 0.
2007       ENDDO
2008       DO K = 1, KM
2009         DO I = 1, IM
2010           if (k .le. kmax(i)) then
2011             IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2012               AUP = 1.
2013               IF(K.Le.KB(I)) AUP = 0.
2014               ADW = 1.
2015               IF(K.GT.JMIN(I)) ADW = 0.
2016               DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
2017               T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
2018               Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
2019               U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
2020               V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
2021               DP = 1000. * DEL(I,K)
2022               DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
2023               DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
2024               DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
2025             ENDIF
2026           endif
2027         ENDDO
2028       ENDDO
2029       DO K = 1, KM
2030         DO I = 1, IM
2031           if (k .le. kmax(i)) then
2032             IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2033 !jfe          QESO(I,k) = 10. * FPVS(T1(I,k))
2034 !
2035               QESO(I,k) = 0.01 * fpvs(T1(I,K))      ! fpvs is in Pa
2036 !
2037               QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
2038 !cmr          QESO(I,k) = MAX(QESO(I,k),1.E-8)
2039               val     =             1.E-8
2040               QESO(I,k) = MAX(QESO(I,k), val )
2041 !
2042 !  cloud water
2043 !
2044               if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
2045                 tem  = DELLAL(I) * XMB(I) * dt2
2046                 tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
2047                 if (QL(I,k,2) .gt. -999.0) then
2048                   QL(I,k,1) = QL(I,k,1) + tem * tem1            ! Ice
2049                   QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1)       ! Water
2050                 else
2051                   tem2      = QL(I,k,1) + tem
2052                   QL(I,k,1) = tem2 * tem1                       ! Ice
2053                   QL(I,k,2) = tem2 - QL(I,k,1)                  ! Water
2054                 endif
2055 !               QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
2056                 dp = 1000. * del(i,k)
2057                 DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
2058               ENDIF
2059             ENDIF
2060           endif
2061         ENDDO
2062       ENDDO
2063 !     IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2064 !       PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2065 !       PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2066 !       PRINT *, '   DELLBAR ='
2067 !       PRINT 6003,  HVAP*DELLbar
2068 !       PRINT *, '   DELLAQ ='
2069 !       PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2070 !       PRINT *, '   DELLAT ='
2071 !       PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I),         &
2072 !    &               K=1,KMAX)
2073 !     ENDIF
2074       DO I = 1, IM
2075         RNTOT(I) = 0.
2076         DELQEV(I) = 0.
2077         DELQ2(I) = 0.
2078         FLG(I) = CNVFLG(I)
2079       ENDDO
2080       DO K = KM, 1, -1
2081         DO I = 1, IM
2082           if (k .le. kmax(i)) then
2083             IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2084               AUP = 1.
2085               IF(K.Le.KB(I)) AUP = 0.
2086               ADW = 1.
2087               IF(K.GT.JMIN(I)) ADW = 0.
2088               rain =  AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2089               RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
2090             ENDIF
2091           endif
2092         ENDDO
2093       ENDDO
2094       DO K = KM, 1, -1
2095         DO I = 1, IM
2096           if (k .le. kmax(i)) then
2097             DELTV(I) = 0.
2098             DELQ(I) = 0.
2099             QEVAP(I) = 0.
2100             IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2101               AUP = 1.
2102               IF(K.Le.KB(I)) AUP = 0.
2103               ADW = 1.
2104               IF(K.GT.JMIN(I)) ADW = 0.
2105               rain =  AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2106               RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
2107             ENDIF
2108             IF(FLG(I).AND.K.LE.KTCON(I)) THEN
2109               evef = EDT(I) * evfact
2110               if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
2111 !             if(SLIMSK(I).eq.1.) evef=.07
2112 !             if(SLIMSK(I).ne.1.) evef = 0.
2113               QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k))                   &
2114      &                 / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
2115               DP = 1000. * DEL(I,K)
2116               IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
2117                 QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
2118                 QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
2119                 DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
2120               ENDIF
2121               if(RN(I).gt.0..and.QCOND(I).LT.0..and.                    &
2122      &           DELQ2(I).gt.RNTOT(I)) THEN
2123                 QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
2124                 FLG(I) = .false.
2125               ENDIF
2126               IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
2127                 Q1(I,k) = Q1(I,k) + QEVAP(I)
2128                 T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
2129                 RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
2130                 DELTV(I) = - ELOCP*QEVAP(I)/DT2
2131                 DELQ(I) =  + QEVAP(I)/DT2
2132                 DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
2133               ENDIF
2134               DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
2135               DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
2136               DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
2137             ENDIF
2138           endif
2139         ENDDO
2140       ENDDO
2141 !      IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2142 !        PRINT *, '   DELLAH ='
2143 !        PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
2144 !        PRINT *, '   DELLAQ ='
2145 !        PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2146 !        PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2147 !        PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2148 !        PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
2149 !CCCC   PRINT *, '   DELLBAR ='
2150 !CCCC   PRINT *,  HVAP*DELLbar
2151 !      ENDIF
2152 !
2153 !  PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
2154 !  IN UNIT OF M INSTEAD OF KG
2155 !
2156       DO I = 1, IM
2157         IF(CNVFLG(I)) THEN
2158 !
2159 !  IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
2160 !    MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
2161 !    HEATING AND THE MOISTENING
2162 !
2163           if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
2164           IF(RN(I).LE.0.) THEN
2165             RN(I) = 0.
2166           ELSE
2167             KTOP(I) = KTCON(I)
2168             KBOT(I) = KBCON(I)
2169             KUO(I) = 1
2170             CLDWRK(I) = AA1(I)
2171           ENDIF
2172         ENDIF
2173       ENDDO
2174       DO K = 1, KM
2175         DO I = 1, IM
2176           if (k .le. kmax(i)) then
2177             IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
2178               T1(I,k) = TO(I,k)
2179               Q1(I,k) = QO(I,k)
2180             ENDIF
2181           endif
2182         ENDDO
2183       ENDDO
2184 !!
2185       RETURN
2186    END SUBROUTINE OSASCNV
2187
2188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2189
2190       SUBROUTINE SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
2191 !
2192       USE MODULE_GFS_MACHINE , ONLY : kind_phys
2193       USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
2194      &,             RD => con_RD
2195
2196       implicit none
2197 !
2198 !     include 'constant.h'
2199 !
2200       integer              IM, IX, KM, KUO(IM)
2201       real(kind=kind_phys) DEL(IX,KM),   PRSI(IX,KM+1), PRSL(IX,KM),    &
2202      &                     PRSLK(IX,KM),                                &
2203      &                     Q(IX,KM),     T(IX,KM),      DT, DPSHC
2204 !
2205 !     Locals
2206 !
2207       real(kind=kind_phys) ck,    cpdt,   dmse,   dsdz1, dsdz2,         &
2208      &                     dsig,  dtodsl, dtodsu, eldq,  g,             &
2209      &                     gocp,  rtdls
2210 !
2211       integer              k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
2212       integer              INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk  &
2213      &,                    KTOPM(IM)
2214 !!
2215 !  PHYSICAL PARAMETERS
2216       PARAMETER(G=GRAV, GOCP=G/CP)
2217 !  BOUNDS OF PARCEL ORIGIN
2218       PARAMETER(KLIFTL=2,KLIFTU=2)
2219       LOGICAL   LSHC(IM)
2220       real(kind=kind_phys) Q2(IM*KM),     T2(IM*KM),                    &
2221      &                     PRSL2(IM*KM),  PRSLK2(IM*KM),                &
2222      &                     AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
2223 !-----------------------------------------------------------------------
2224 !  COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
2225 !  AND MOIST STATIC INSTABILITY.
2226       DO I=1,IM
2227         LSHC(I)=.FALSE.
2228       ENDDO
2229       DO K=1,KM-1
2230         DO I=1,IM
2231           IF(KUO(I).EQ.0) THEN
2232             ELDQ    = HVAP*(Q(I,K)-Q(I,K+1))
2233             CPDT    = CP*(T(I,K)-T(I,K+1))
2234             RTDLS   = (PRSL(I,K)-PRSL(I,K+1)) /                         &
2235      &                 PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
2236             DMSE    = ELDQ+CPDT-RTDLS
2237             LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2238           ENDIF
2239         ENDDO
2240       ENDDO
2241       N2 = 0
2242       DO I=1,IM
2243         IF(LSHC(I)) THEN
2244           N2         = N2 + 1
2245           INDEX2(N2) = I
2246         ENDIF
2247       ENDDO
2248       IF(N2.EQ.0) RETURN
2249       DO K=1,KM
2250         KK = (K-1)*N2
2251         DO I=1,N2
2252           IK         = KK + I
2253           ii         = index2(i)
2254           Q2(IK)     = Q(II,K)
2255           T2(IK)     = T(II,K)
2256           PRSL2(IK)  = PRSL(II,K)
2257           PRSLK2(IK) = PRSLK(II,K)
2258         ENDDO
2259       ENDDO
2260       do i=1,N2
2261         ktopm(i) = KM
2262       enddo
2263       do k=2,KM
2264         do i=1,N2
2265           ii = index2(i)
2266           if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
2267         enddo
2268       enddo
2269
2270 !-----------------------------------------------------------------------
2271 !  COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
2272 !  CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
2273       CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2,           &
2274      &            KLCL,KBOT,KTOP,AL,AU)
2275       DO I=1,N2
2276         KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
2277         KTOP(I) = min(KTOP(I)+1, ktopm(i))
2278         LSHC(I) = .FALSE.
2279       ENDDO
2280       DO K=1,KM-1
2281         KK = (K-1)*N2
2282         DO I=1,N2
2283           IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
2284             IK      = KK + I
2285             IKU     = IK + N2
2286             ELDQ    = HVAP * (Q2(IK)-Q2(IKU))
2287             CPDT    = CP   * (T2(IK)-T2(IKU))
2288             RTDLS   = (PRSL2(IK)-PRSL2(IKU)) /                          &
2289      &                 PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
2290             DMSE    = ELDQ + CPDT - RTDLS
2291             LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2292             AU(IK)  = G/RTDLS
2293           ENDIF
2294         ENDDO
2295       ENDDO
2296       K1=KM+1
2297       K2=0
2298       DO I=1,N2
2299         IF(.NOT.LSHC(I)) THEN
2300           KBOT(I) = KM+1
2301           KTOP(I) = 0
2302         ENDIF
2303         K1 = MIN(K1,KBOT(I))
2304         K2 = MAX(K2,KTOP(I))
2305       ENDDO
2306       KT = K2-K1+1
2307       IF(KT.LT.2) RETURN
2308 !-----------------------------------------------------------------------
2309 !  SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
2310 !  COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
2311 !  EXPAND FINAL FIELDS.
2312       KK = (K1-1) * N2
2313       DO I=1,N2
2314         IK     = KK + I
2315         AD(IK) = 1.
2316       ENDDO
2317 !
2318 !     DTODSU=DT/DEL(K1)
2319       DO K=K1,K2-1
2320 !       DTODSL=DTODSU
2321 !       DTODSU=   DT/DEL(K+1)
2322 !       DSIG=SL(K)-SL(K+1)
2323         KK = (K-1) * N2
2324         DO I=1,N2
2325           ii     = index2(i)
2326           DTODSL = DT/DEL(II,K)
2327           DTODSU = DT/DEL(II,K+1)
2328           DSIG   = PRSL(II,K) - PRSL(II,K+1)
2329           IK     = KK + I
2330           IKU    = IK + N2
2331           IF(K.EQ.KBOT(I)) THEN
2332             CK=1.5
2333           ELSEIF(K.EQ.KTOP(I)-1) THEN
2334             CK=1.
2335           ELSEIF(K.EQ.KTOP(I)-2) THEN
2336             CK=3.
2337           ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
2338             CK=5.
2339           ELSE
2340             CK=0.
2341           ENDIF
2342           DSDZ1   = CK*DSIG*AU(IK)*GOCP
2343           DSDZ2   = CK*DSIG*AU(IK)*AU(IK)
2344           AU(IK)  = -DTODSL*DSDZ2
2345           AL(IK)  = -DTODSU*DSDZ2
2346           AD(IK)  = AD(IK)-AU(IK)
2347           AD(IKU) = 1.-AL(IK)
2348           T2(IK)  = T2(IK)+DTODSL*DSDZ1
2349           T2(IKU) = T2(IKU)-DTODSU*DSDZ1
2350         ENDDO
2351       ENDDO
2352       IK1=(K1-1)*N2+1
2353       CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1),      &
2354      &                                  AU(IK1),Q2(IK1),T2(IK1))
2355       DO K=K1,K2
2356         KK = (K-1)*N2
2357         DO I=1,N2
2358           IK = KK + I
2359           Q(INDEX2(I),K) = Q2(IK)
2360           T(INDEX2(I),K) = T2(IK)
2361         ENDDO
2362       ENDDO
2363 !-----------------------------------------------------------------------
2364       RETURN
2365       END SUBROUTINE SHALCV
2366 !-----------------------------------------------------------------------
2367       SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
2368 !yt      INCLUDE DBTRIDI2;
2369 !!
2370       USE MODULE_GFS_MACHINE , ONLY : kind_phys
2371       implicit none
2372       integer             k,n,l,i
2373       real(kind=kind_phys) fk
2374 !!
2375       real(kind=kind_phys)                                              &
2376      &          CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N),            &
2377      &          AU(L,N-1),A1(L,N),A2(L,N)
2378 !-----------------------------------------------------------------------
2379       DO I=1,L
2380         FK=1./CM(I,1)
2381         AU(I,1)=FK*CU(I,1)
2382         A1(I,1)=FK*R1(I,1)
2383         A2(I,1)=FK*R2(I,1)
2384       ENDDO
2385       DO K=2,N-1
2386         DO I=1,L
2387           FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
2388           AU(I,K)=FK*CU(I,K)
2389           A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
2390           A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
2391         ENDDO
2392       ENDDO
2393       DO I=1,L
2394         FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
2395         A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
2396         A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
2397       ENDDO
2398       DO K=N-1,1,-1
2399         DO I=1,L
2400           A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
2401           A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
2402         ENDDO
2403       ENDDO
2404 !-----------------------------------------------------------------------
2405       RETURN
2406       END SUBROUTINE TRIDI2T3
2407 !-----------------------------------------------------------------------
2408
2409       SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV,             &
2410      &                  KLCL,KBOT,KTOP,TCLD,QCLD)
2411 !yt      INCLUDE DBMSTADB;
2412 !!
2413       USE MODULE_GFS_MACHINE, ONLY : kind_phys
2414       USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
2415       USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
2416
2417       implicit none
2418 !!
2419 !     include 'constant.h'
2420 !!
2421       integer              k,k1,k2,km,i,im
2422       real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
2423       real(kind=kind_phys) tma,tvcld,tvenv
2424 !!
2425       real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM),      &
2426      &                     QENV(IM,KM), TCLD(IM,KM),  QCLD(IM,KM)
2427       INTEGER              KLCL(IM),    KBOT(IM),      KTOP(IM)
2428 !  LOCAL ARRAYS
2429       real(kind=kind_phys) SLKMA(IM), THEMA(IM)
2430 !-----------------------------------------------------------------------
2431 !  DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
2432 !  COMPUTE ITS LIFTING CONDENSATION LEVEL.
2433 !
2434       DO I=1,IM
2435         SLKMA(I) = 0.
2436         THEMA(I) = 0.
2437       ENDDO
2438       DO K=K1,K2
2439         DO I=1,IM
2440           PV   = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
2441           TDPD = TENV(I,K)-FTDP(PV)
2442           IF(TDPD.GT.0.) THEN
2443             TLCL   = FTLCL(TENV(I,K),TDPD)
2444             SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
2445           ELSE
2446             TLCL   = TENV(I,K)
2447             SLKLCL = PRSLK(I,K)
2448           ENDIF
2449           THELCL=FTHE(TLCL,SLKLCL)
2450           IF(THELCL.GT.THEMA(I)) THEN
2451             SLKMA(I) = SLKLCL
2452             THEMA(I) = THELCL
2453           ENDIF
2454         ENDDO
2455       ENDDO
2456 !-----------------------------------------------------------------------
2457 !  SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
2458 !  THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
2459       DO I=1,IM
2460         KLCL(I)=KM+1
2461         KBOT(I)=KM+1
2462         KTOP(I)=0
2463       ENDDO
2464       DO K=1,KM
2465         DO I=1,IM
2466           TCLD(I,K)=0.
2467           QCLD(I,K)=0.
2468         ENDDO
2469       ENDDO
2470       DO K=K1,KM
2471         DO I=1,IM
2472           IF(PRSLK(I,K).LE.SLKMA(I)) THEN
2473             KLCL(I)=MIN(KLCL(I),K)
2474             CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
2475 !           TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
2476             TVCLD=TMA*(1.+FV*QMA)
2477             TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
2478             IF(TVCLD.GT.TVENV) THEN
2479               KBOT(I)=MIN(KBOT(I),K)
2480               KTOP(I)=MAX(KTOP(I),K)
2481               TCLD(I,K)=TMA-TENV(I,K)
2482               QCLD(I,K)=QMA-QENV(I,K)
2483             ENDIF
2484           ENDIF
2485         ENDDO
2486       ENDDO
2487 !-----------------------------------------------------------------------
2488       RETURN
2489       END SUBROUTINE MSTADBT3
2490
2491 #if (EM_CORE == 1)
2492 !   random seeds - ZCX
2493       SUBROUTINE init_random_seed()
2494             INTEGER :: i, n, clock
2495             INTEGER, DIMENSION(:), ALLOCATABLE :: seed
2496
2497             CALL RANDOM_SEED(size = n)
2498             ALLOCATE(seed(n))
2499
2500             CALL SYSTEM_CLOCK(COUNT=clock)
2501
2502             seed = clock + 37 * (/ (i - 1, i = 1, n) /)
2503             CALL RANDOM_SEED(PUT = seed)
2504
2505             DEALLOCATE(seed)
2506       END SUBROUTINE
2507 #endif
2508       END MODULE module_cu_osas