wrfv2_fire/phys/module_cu_nsas.F

   1 !
   2 !
   3 !
   4 !
   5 MODULE module_cu_nsas
   6 CONTAINS
   7 !
   8 !-------------------------------------------------------------------------------
   9    subroutine cu_nsas(dt,p3di,p3d,pi3d,qc3d,qi3d,rho3d,itimestep,stepcu,       &
  10                      hbot,htop,cu_act_flag,cudt,curr_secs,adapt_step_flag,     &
  11                      cudtacttime,                                              &
  12                      rthcuten,rqvcuten,rqccuten,rqicuten,                      &
  13                      rucuten,rvcuten,                                          &
  14                      qv3d,t3d,raincv,pratec,xland,dz8w,w,u3d,v3d,              &
  15                      hpbl,hfx,qfx,                                             &
  16                      mp_physics,                                               &
  17                      p_qc,p_qi,p_first_scalar,                                 &
  18                      pgcon,                                                    &
  19                      cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,                      &
  20                      cice,xls,psat,f_qi,f_qc,                                  &
  21                      ids,ide, jds,jde, kds,kde,                                &
  22                      ims,ime, jms,jme, kms,kme,                                &
  23                      its,ite, jts,jte, kts,kte)
  24 !-------------------------------------------------------------------------------
  25   implicit none
  26 !-------------------------------------------------------------------------------
  27 !
  28 !-- dt          time step (s)
  29 !-- p3di        3d pressure (pa) at interface level
  30 !-- p3d         3d pressure (pa)
  31 !-- pi3d        3d exner function (dimensionless)
  32 !-- z           height above sea level (m)
  33 !-- qc3d        cloud water mixing ratio (kg/kg)
  34 !-- qi3d        cloud ice mixing ratio (kg/kg)
  35 !-- qv3d        3d water vapor mixing ratio (kg/kg)
  36 !-- t3d         temperature (k)
  37 !-- raincv      cumulus scheme precipitation (mm)
  38 !-- w           vertical velocity (m/s)
  39 !-- dz8w        dz between full levels (m)
  40 !-- u3d         3d u-velocity interpolated to theta points (m/s)
  41 !-- v3d         3d v-velocity interpolated to theta points (m/s)
  42 !-- ids         start index for i in domain
  43 !-- ide         end index for i in domain
  44 !-- jds         start index for j in domain
  45 !-- jde         end index for j in domain
  46 !-- kds         start index for k in domain
  47 !-- kde         end index for k in domain
  48 !-- ims         start index for i in memory
  49 !-- ime         end index for i in memory
  50 !-- jms         start index for j in memory
  51 !-- jme         end index for j in memory
  52 !-- kms         start index for k in memory
  53 !-- kme         end index for k in memory
  54 !-- its         start index for i in tile
  55 !-- ite         end index for i in tile
  56 !-- jts         start index for j in tile
  57 !-- jte         end index for j in tile
  58 !-- kts         start index for k in tile
  59 !-- kte         end index for k in tile
  60 !-------------------------------------------------------------------------------
  61   integer,  intent(in   )   ::       ids,ide, jds,jde, kds,kde,                &
  62                                      ims,ime, jms,jme, kms,kme,                &
  63                                      its,ite, jts,jte, kts,kte,                &
  64                                      itimestep, stepcu,                        &
  65                                      p_qc,p_qi,p_first_scalar
  66 !
  67   real,     intent(in   )   ::      cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2,       &
  68                                     cice,xls,psat
  69   real,     intent(in   )   ::      dt
  70   real,     optional, intent(in   )   ::      pgcon
  71 !
  72   real,     dimension( ims:ime, kms:kme, jms:jme ),optional                  , &
  73             intent(inout)   ::                                       rthcuten, &
  74                                                                       rucuten, &
  75                                                                       rvcuten, &
  76                                                                      rqccuten, &
  77                                                                      rqicuten, &
  78                                                                      rqvcuten
  79   logical, optional ::                                              F_QC,F_QI
  80   real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
  81             intent(in   )   ::                                           qv3d, &
  82                                                                          qc3d, &
  83                                                                          qi3d, &
  84                                                                         rho3d, &
  85                                                                           p3d, &
  86                                                                          pi3d, &
  87                                                                           t3d
  88   real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
  89             intent(in   )   ::                                           p3di
  90   real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
  91             intent(in   )   ::                                           dz8w, &
  92                                                                             w
  93   real,     dimension( ims:ime, jms:jme )                                    , &
  94             intent(inout) ::                                           raincv, &
  95                                                                        pratec
  96   real,     dimension( ims:ime, jms:jme )                                    , &
  97             intent(out) ::                                               hbot, &
  98                                                                          htop
  99   real,     dimension( ims:ime, jms:jme )                                    , &
 100             intent(in   ) ::                                            xland
 101 !
 102   real,     dimension( ims:ime, kms:kme, jms:jme )                           , &
 103              intent(in   )   ::                                           u3d, &
 104                                                                           v3d
 105   logical,  dimension( ims:ime, jms:jme )                                    , &
 106             intent(inout) ::                                      cu_act_flag
 107   real, intent(   in) ::                                                 cudt
 108   real, intent(   in) ::                                            curr_secs
 109   logical, intent(    in) , optional ::                       adapt_step_flag
 110
 111   real, intent (inout) ::                                         cudtacttime
 112 !
 113   real,     dimension( ims:ime, jms:jme )                                    , &
 114              intent(in   )   ::                                          hpbl, &
 115                                                                           hfx, &
 116                                                                           qfx
 117   integer,   intent(in   )   ::                                    mp_physics
 118   integer :: ncloud
 119 !
 120 !local
 121 !
 122   real,  dimension( its:ite, jts:jte )  ::                            raincv1, &
 123                                                                       raincv2, &
 124                                                                       pratec1, &
 125                                                                       pratec2
 126   real,   dimension( its:ite, kts:kte )  ::                               del, &
 127                                                                         prsll, &
 128                                                                           dot, &
 129                                                                            u1, &
 130                                                                            v1, &
 131                                                                            t1, &
 132                                                                           q1,  &
 133                                                                           qc2, &
 134                                                                           qi2
 135   real,   dimension( its:ite, kts:kte+1 )  ::                           prsii, &
 136                                                                           zii
 137   real,   dimension( its:ite, kts:kte )  ::                               zll
 138   real,   dimension( its:ite)  ::                                         rain
 139   real ::                                                          delt,rdelt
 140   integer, dimension (its:ite)  ::                                       kbot, &
 141                                                                          ktop, &
 142                                                                           kuo
 143   real    :: pgcon_use
 144   logical :: run_param , doing_adapt_dt , decided
 145   integer ::  i,j,k,kp
 146 !
 147 !-------------------------------------------------------------------------------
 148 ! microphysics scheme --> ncloud
 149    if (mp_physics .eq. 0) then
 150      ncloud = 0
 151    elseif ( mp_physics .eq. 1 .or. mp_physics .eq. 3 ) then
 152      ncloud = 1
 153    else
 154      ncloud = 2
 155    endif
 156 !
 157 !-------------------------------------------------------------------------------
 158 !
 159 !*** check to see if this is a convection timestep
 160 !
 161
 162 !  Initialization for adaptive time step.
 163
 164    doing_adapt_dt = .FALSE.
 165    IF ( PRESENT(adapt_step_flag) ) THEN
 166       IF ( adapt_step_flag ) THEN
 167          doing_adapt_dt = .TRUE.
 168          IF ( cudtacttime .EQ. 0. ) THEN
 169             cudtacttime = curr_secs + cudt*60.
 170          END IF
 171       END IF
 172    END IF
 173
 174 !  Do we run through this scheme or not?
 175
 176 !    Test 1:  If this is the initial model time, then yes.
 177 !                ITIMESTEP=1
 178 !    Test 2:  If the user asked for the cumulus to be run every time step, then yes.
 179 !                CUDT=0 or STEPCU=1
 180 !    Test 3:  If not adaptive dt, and this is on the requested cumulus frequency, then yes.
 181 !                MOD(ITIMESTEP,STEPCU)=0
 182 !    Test 4:  If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
 183 !                CURR_SECS >= CUDTACTTIME
 184
 185 !  If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
 186 !  to TRUE.  The decided flag says that one of these tests was able to say "yes", run the scheme.
 187 !  We only proceed to other tests if the previous tests all have left decided as FALSE.
 188
 189 !  If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
 190 !  cumulus run.
 191
 192    decided = .FALSE.
 193    run_param = .FALSE.
 194    IF ( ( .NOT. decided ) .AND. &
 195         ( itimestep .EQ. 1 ) ) THEN
 196       run_param   = .TRUE.
 197       decided     = .TRUE.
 198    END IF
 199
 200    IF ( ( .NOT. decided ) .AND. &
 201         ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
 202       run_param   = .TRUE.
 203       decided     = .TRUE.
 204    END IF
 205
 206    IF ( ( .NOT. decided ) .AND. &
 207         ( .NOT. doing_adapt_dt ) .AND. &
 208         ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
 209       run_param   = .TRUE.
 210       decided     = .TRUE.
 211    END IF
 212
 213    IF ( ( .NOT. decided ) .AND. &
 214         ( doing_adapt_dt ) .AND. &
 215         ( curr_secs .GE. cudtacttime ) ) THEN
 216       run_param   = .TRUE.
 217       decided     = .TRUE.
 218       cudtacttime = curr_secs + cudt*60
 219    END IF
 220
 221 !-------------------------------------------------------------------------------
 222
 223       if(present(pgcon)) then
 224          pgcon_use  = pgcon
 225       else
 226 !        pgcon_use  = 0.7     ! Gregory et al. (1997, QJRMS)
 227          pgcon_use  = 0.55    ! Zhang & Wu (2003,JAS)
 228          ! 0.55 is a physically-based value used by GFS
 229          ! HWRF uses 0.2, for model tuning purposes
 230       endif
 231
 232    if(run_param) then
 233       do j=jts,jte
 234         do i=its,ite
 235           cu_act_flag(i,j)=.TRUE.
 236         enddo
 237       enddo
 238       delt=dt*stepcu
 239       rdelt=1./delt
 240 !
 241    do j = jts,jte  !outer most J_loop
 242       do k = kts,kte
 243         kp = k+1
 244         do i = its,ite
 245           dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
 246           prsll(i,k)=p3d(i,k,j)*0.001
 247           prsii(i,k)=p3di(i,k,j)*0.001
 248         enddo
 249       enddo
 250       do i = its,ite
 251         prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
 252       enddo
 253 !
 254       do i=its,ite
 255         zii(i,1)=0.0
 256       enddo
 257 !
 258       do k=kts,kte
 259         do i=its,ite
 260           zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
 261         enddo
 262       enddo
 263 !
 264       do k=kts,kte
 265         do i=its,ite
 266           zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
 267         enddo
 268       enddo
 269 !
 270       do k=kts,kte
 271         do i=its,ite
 272           del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
 273           u1(i,k)=u3d(i,k,j)
 274           v1(i,k)=v3d(i,k,j)
 275           t1(i,k)=t3d(i,k,j)
 276           q1(i,k)=qv3d(i,k,j)
 277           qi2(i,k) = qi3d(i,k,j)
 278           qc2(i,k) = qc3d(i,k,j)
 279         enddo
 280       enddo
 281 !
 282 ! NCEP SAS
 283       call nsas2d(delt=dt,del=del(its,kts),prsl=prsll(its,kts),                &
 284               prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts),       &
 285               zi=zii(its,kts),ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
 286               q1=q1(its,kts),t1=t1(its,kts),rain=rain(its),                    &
 287               kbot=kbot(its),ktop=ktop(its),                                   &
 288               kuo=kuo(its),                                                    &
 289               lat=j,slimsk=xland(ims,j),dot=dot(its,kts),                      &
 290               u1=u1(its,kts), v1=v1(its,kts),                                  &
 291               cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv,                      &
 292               rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2,                               &
 293               cice=cice,xls=xls,psat=psat,                                     &
 294               pgcon=pgcon_use,                                                 &
 295               ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,               &
 296               ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,               &
 297               its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
 298 !
 299       do i=its,ite
 300         pratec1(i,j)=rain(i)*1000./(stepcu*dt)
 301         raincv1(i,j)=rain(i)*1000./(stepcu)
 302       enddo
 303 !
 304 ! NCEP SCV
 305       call nscv2d(delt=dt,del=del(its,kts),prsl=prsll(its,kts),                &
 306               prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts),       &
 307               zi=zii(its,kts),ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
 308               q1=q1(its,kts),t1=t1(its,kts),rain=rain(its),                    &
 309               kbot=kbot(its),ktop=ktop(its),                                   &
 310               kuo=kuo(its),                                                    &
 311               slimsk=xland(ims,j),dot=dot(its,kts),                            &
 312               u1=u1(its,kts), v1=v1(its,kts),                                  &
 313               cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv,                      &
 314               rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2,                               &
 315               cice=cice,xls=xls,psat=psat,                                     &
 316               hpbl=hpbl(ims,j),hfx=hfx(ims,j),qfx=qfx(ims,j),                  &
 317               pgcon=pgcon_use,                                                 &
 318               ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,               &
 319               ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,               &
 320               its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte   )
 321 !
 322    do i=its,ite
 323      pratec2(i,j)=rain(i)*1000./(stepcu*dt)
 324      raincv2(i,j)=rain(i)*1000./(stepcu)
 325    enddo
 326 !
 327    do i=its,ite
 328      raincv(i,j) = raincv1(i,j) + raincv2(i,j)
 329      pratec(i,j) = pratec1(i,j) + pratec2(i,j)
 330      hbot(i,j) = kbot(i)
 331      htop(i,j) = ktop(i)
 332    enddo
 333 !
 334       IF(PRESENT(rthcuten).AND.PRESENT(rqvcuten)) THEN
 335       do k = kts,kte
 336         do i= its,ite
 337           rthcuten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
 338           rqvcuten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
 339         enddo
 340       enddo
 341       ENDIF
 342 !
 343       IF(PRESENT(rucuten).AND.PRESENT(rvcuten)) THEN
 344       do k = kts,kte
 345         do i= its,ite
 346           rucuten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
 347           rvcuten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
 348         enddo
 349       enddo
 350       ENDIF
 351 !
 352       if(PRESENT( rqicuten )) THEN
 353         IF ( F_QI ) THEN
 354         do k=kts,kte
 355           do i=its,ite
 356             rqicuten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
 357           enddo
 358         enddo
 359         endif
 360       endif
 361       if(PRESENT( rqccuten )) THEN
 362         IF ( F_QC ) THEN
 363         do k=kts,kte
 364           do i=its,ite
 365             rqccuten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
 366           enddo
 367         enddo
 368         endif
 369       endif
 370 !
 371    enddo ! outer most J_loop
 372   endif
 373 !
 374    end subroutine cu_nsas
 375 !
 376 !==============================================================================
 377 ! NCEP SAS (Deep Convection Scheme)
 378 !==============================================================================
 379    subroutine nsas2d(delt,del,prsl,prsi,prslk,zl,zi,                           &
 380             ncloud,                                                            &
 381             qc2,qi2,                                                           &
 382             q1,t1,rain,kbot,ktop,                                              &
 383             kuo,                                                               &
 384             lat,slimsk,dot,u1,v1,cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2,     &
 385             cice,xls,psat,                                                     &
 386             pgcon,                                                             &
 387             ids,ide, jds,jde, kds,kde,                                         &
 388             ims,ime, jms,jme, kms,kme,                                         &
 389             its,ite, jts,jte, kts,kte)
 390 !
 391 !------------------------------------------------------------------------------
 392 !
 393 ! subprogram:    phys_cps_sas      computes convective heating and moistening
 394 !                                                      and momentum transport
 395 !
 396 ! abstract: computes convective heating and moistening using a one
 397 !   cloud type arakawa-schubert convection scheme originally developed
 398 !   by georg grell. the scheme has been revised at ncep since initial
 399 !   implementation in 1993. it includes updraft and downdraft effects.
 400 !   the closure is the cloud work function. both updraft and downdraft
 401 !   are assumed to be saturated and the heating and moistening are
 402 !   accomplished by the compensating environment. convective momentum transport
 403 !   is taken into account. the name comes from "simplified arakawa-schubert
 404 !   convection parameterization".
 405 !
 406 ! developed by hua-lu pan, wan-shu wu, songyou hong, and jongil han
 407 !   implemented into wrf by kyosun sunny lim and songyou hong
 408 !   module with cpp-based options is available in grims
 409 !
 410 ! program history log:
 411 !   92-03-01  hua-lu pan       operational development
 412 !   96-03-01  song-you hong    revised closure, and trigger
 413 !   99-03-01  hua-lu pan       multiple clouds
 414 !   06-03-01  young-hwa byun   closure based on moisture convergence (optional)
 415 !   09-10-01  jung-eun kim     f90 format with standard physics modules
 416 !   10-07-01  jong-il han      revised cloud model,trigger, as in gfs july 2010
 417 !   10-12-01  kyosun sunny lim wrf compatible version
 418 !
 419 !
 420 ! usage:    call phys_cps_sas(delt,delx,del,prsl,prsi,prslk,prsik,zl,zi,       &
 421 !                             q2,q1,t1,u1,v1,rcs,slimsk,dot,cldwrk,rain,       &
 422 !                             jcap,ncloud,lat,kbot,ktop,kuo,                   &
 423 !                             ids,ide, jds,jde, kds,kde,                       &
 424 !                             ims,ime, jms,jme, kms,kme,                       &
 425 !                             its,ite, jts,jte, kts,kte)
 426 !
 427 !   delt     - real model integration time step
 428 !   delx     - real model grid interval
 429 !   del      - real (kms:kme) sigma layer thickness
 430 !   prsl     - real (ims:ime,kms:kme) pressure values
 431 !   prsi     - real (ims:ime,kms:kme) pressure values at interface level
 432 !   prslk    - real (ims:ime,kms:kme) pressure values to the kappa
 433 !   prsik    - real (ims:ime,kms:kme) pressure values to the kappa at interface lev.
 434 !   zl       - real (ims:ime,kms:kme) height above sea level
 435 !   zi       - real (ims:ime,kms:kme) height above sea level at interface level
 436 !   rcs      - real
 437 !   slimsk   - real (ims:ime) land(1),sea(0), ice(2) flag
 438 !   dot      - real (ims:ime,kms:kme) vertical velocity
 439 !   jcap     - integer spectral truncation
 440 !   ncloud   - integer no_cloud(0),no_ice(1),cloud+ice(2)
 441 !   lat      - integer  current latitude index
 442 !
 443 ! output argument list:
 444 !   q2       - real (ims:ime,kms:kme) detrained hydrometeors in kg/kg
 445 !            - in case of the  --> qc2(cloud), qi2(ice)
 446 !   q1       - real (ims:ime,kms:kme) adjusted specific humidity in kg/kg
 447 !   t1       - real (ims:ime,kms:kme) adjusted temperature in kelvin
 448 !   u1       - real (ims:ime,kms:kme) adjusted zonal wind in m/s
 449 !   v1       - real (ims:ime,kms:kme) adjusted meridional wind in m/s
 450 !   cldwrk   - real (ims:ime) cloud work function
 451 !   rain     - real (ims:ime) convective rain in meters
 452 !   kbot     - integer (ims:ime) cloud bottom level
 453 !   ktop     - integer (ims:ime) cloud top level
 454 !   kuo      - integer (ims:ime) bit flag indicating deep convection
 455 !
 456 ! subprograms called:
 457 !   fpvs     - function to compute saturation vapor pressure
 458 !
 459 ! remarks: function fpvs is inlined by fpp.
 460 !          nonstandard automatic arrays are used.
 461 !
 462 ! references :
 463 !   pan and wu    (1995, ncep office note)
 464 !   hong and pan  (1998, mon wea rev)
 465 !   park and hong (2007,jmsj)
 466 !   byun and hong (2007, mon wea rev)
 467 !   han and pan   (2011, wea. forecasting)
 468 !
 469 !------------------------------------------------------------------------------
 470 !------------------------------------------------------------------------------
 471    implicit none
 472 !------------------------------------------------------------------------------
 473 !
 474 ! model tunable parameters
 475 !
 476    real,parameter  ::  alphal = 0.5,    alphas = 0.5
 477    real,parameter  ::  betal  = 0.05,   betas  = 0.05
 478    real,parameter  ::  pdpdwn = 0.0,    pdetrn = 200.0
 479    real,parameter  ::  c0     = 0.002,  c1     = 0.002
 480    real,parameter  ::  xlamdd = 1.0e-4, xlamde = 1.0e-4
 481    real,parameter  ::  clam   = 0.1,    cxlamu = 1.0e-4
 482    real,parameter  ::  aafac  = 0.1
 483    real,parameter  ::  dthk=25.
 484    real,parameter  ::  cincrmax = 180.,cincrmin = 120.
 485    real,parameter  ::  W1l = -8.E-3
 486    real,parameter  ::  W2l = -4.E-2
 487    real,parameter  ::  W3l = -5.E-3
 488    real,parameter  ::  W4l = -5.E-4
 489    real,parameter  ::  W1s = -2.E-4
 490    real,parameter  ::  W2s = -2.E-3
 491    real,parameter  ::  W3s = -1.E-3
 492    real,parameter  ::  W4s = -2.E-5
 493    real,parameter  ::  mbdt = 10., edtmaxl = 0.3, edtmaxs = 0.3
 494    real,parameter  ::  evfacts = 0.3, evfactl = 0.3
 495 !
 496    real,parameter  ::  tf=233.16,tcr=263.16,tcrf=1.0/(tcr-tf)
 497    real,parameter  ::  xk1=2.e-5,xlhor=3.e4,xhver=5000.,theimax=1.
 498    real,parameter  ::  xc1=1.e-7,xc2=1.e4,xc3=3.e3,ecesscr=3.0,edtk1=3.e4
 499 !
 500 !  passing variables
 501 !
 502    real            ::  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
 503    real            ::  pi_,qmin_,t0c_,cice,xlv0,xls,psat
 504    integer         ::  lat,                                                    &
 505                        ncloud,                                                 &
 506                        ids,ide, jds,jde, kds,kde,                              &
 507                        ims,ime, jms,jme, kms,kme,                              &
 508                        its,ite, jts,jte, kts,kte
 509 !
 510    real            ::  delt,rcs
 511    real            ::  del(its:ite,kts:kte),                                   &
 512                        prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme),           &
 513                        prsi(its:ite,kts:kte+1),                                &
 514                        zl(its:ite,kts:kte),zi(its:ite,kts:kte+1),              &
 515                        q1(its:ite,kts:kte),t1(its:ite,kts:kte),                &
 516                        u1(its:ite,kts:kte),v1(its:ite,kts:kte),                &
 517                        dot(its:ite,kts:kte)
 518    real            ::  qi2(its:ite,kts:kte)
 519    real            ::  qc2(its:ite,kts:kte)
 520 !
 521    real            ::  rain(its:ite)
 522    integer         ::  kbot(its:ite),ktop(its:ite),kuo(its:ite)
 523    real            ::  slimsk(ims:ime)
 524    real            ::  pgcon
 525 !
 526 !
 527 !  local variables and arrays
 528 !
 529    integer         ::  i,k,kmax,kbmax,kbm,jmn,indx,indp,kts1,kte1,kmax1,kk
 530    real            ::  p(its:ite,kts:kte),pdot(its:ite),acrtfct(its:ite)
 531    real            ::  uo(its:ite,kts:kte),vo(its:ite,kts:kte)
 532    real            ::  to(its:ite,kts:kte),qo(its:ite,kts:kte)
 533    real            ::  hcko(its:ite,kts:kte)
 534    real            ::  qcko(its:ite,kts:kte),eta(its:ite,kts:kte)
 535    real            ::  etad(its:ite,kts:kte)
 536    real            ::  qrcdo(its:ite,kts:kte)
 537    real            ::  pwo(its:ite,kts:kte),pwdo(its:ite,kts:kte)
 538    real            ::  dtconv(its:ite)
 539    real            ::  deltv(its:ite),acrt(its:ite)
 540    real            ::  qeso(its:ite,kts:kte)
 541    real            ::  tvo(its:ite,kts:kte),dbyo(its:ite,kts:kte)
 542    real            ::  heo(its:ite,kts:kte),heso(its:ite,kts:kte)
 543    real            ::  qrcd(its:ite,kts:kte)
 544    real            ::  dellah(its:ite,kts:kte),dellaq(its:ite,kts:kte)
 545 !
 546    integer         ::  kb(its:ite),kbcon(its:ite)
 547    integer         ::  kbcon1(its:ite)
 548    real            ::  hmax(its:ite),delq(its:ite)
 549    real            ::  hkbo(its:ite),qkbo(its:ite),pbcdif(its:ite)
 550    integer         ::  kbds(its:ite),lmin(its:ite),jmin(its:ite)
 551    integer         ::  ktcon(its:ite)
 552    integer         ::  ktcon1(its:ite)
 553    integer         ::  kbdtr(its:ite),kpbl(its:ite)
 554    integer         ::  klcl(its:ite),ktdown(its:ite)
 555    real            ::  vmax(its:ite)
 556    real            ::  hmin(its:ite),pwavo(its:ite)
 557    real            ::  aa1(its:ite),vshear(its:ite)
 558    real            ::  qevap(its:ite)
 559    real            ::  edt(its:ite)
 560    real            ::  edto(its:ite),pwevo(its:ite)
 561    real            ::  qcond(its:ite)
 562    real            ::  hcdo(its:ite,kts:kte)
 563    real            ::  ddp(its:ite),pp2(its:ite)
 564    real            ::  qcdo(its:ite,kts:kte)
 565    real            ::  adet(its:ite),aatmp(its:ite)
 566    real            ::  xhkb(its:ite),xqkb(its:ite)
 567    real            ::  xpwav(its:ite),xpwev(its:ite),xhcd(its:ite,kts:kte)
 568    real            ::  xaa0(its:ite),f(its:ite),xk(its:ite)
 569    real            ::  xmb(its:ite)
 570    real            ::  edtx(its:ite),xqcd(its:ite,kts:kte)
 571    real            ::  hsbar(its:ite),xmbmax(its:ite)
 572    real            ::  xlamb(its:ite,kts:kte),xlamd(its:ite)
 573    real            ::  excess(its:ite)
 574    real            ::  plcl(its:ite)
 575    real            ::  delhbar(its:ite),delqbar(its:ite),deltbar(its:ite)
 576    real,save       ::  pcrit(15), acritt(15)
 577    real            ::  acrit(15)
 578    real            ::  qcirs(its:ite,kts:kte),qrski(its:ite)
 579    real            ::  dellal(its:ite,kts:kte)
 580    real            ::  rntot(its:ite),delqev(its:ite),delq2(its:ite)
 581 !
 582    real            ::  fent1(its:ite,kts:kte),fent2(its:ite,kts:kte)
 583    real            ::  frh(its:ite,kts:kte)
 584    real            ::  xlamud(its:ite),sumx(its:ite)
 585    real            ::  aa2(its:ite)
 586    real            ::  ucko(its:ite,kts:kte),vcko(its:ite,kts:kte)
 587    real            ::  ucdo(its:ite,kts:kte),vcdo(its:ite,kts:kte)
 588    real            ::  dellau(its:ite,kts:kte),dellav(its:ite,kts:kte)
 589    real            ::  delubar(its:ite),delvbar(its:ite)
 590    real            ::  qlko_ktcon(its:ite)
 591 !
 592    real            ::  alpha,beta,                                             &
 593                        dt2,dtmin,dtmax,dtmaxl,dtmaxs,                          &
 594                        el2orc,eps,fact1,fact2,                                 &
 595                        tem,tem1,cincr
 596    real            ::  dz,dp,es,pprime,qs,                                     &
 597                        dqsdp,desdt,dqsdt,gamma,                                &
 598                        dt,dq,po,thei,delza,dzfac,                              &
 599                        thec,theb,thekb,thekh,theavg,thedif,                    &
 600                        omgkb,omgkbp1,omgdif,omgfac,heom,rh,thermal,chi,        &
 601                        factor,onemf,dz1,qrch,etah,qlk,qc,rfact,shear,          &
 602                        e1,dh,deta,detad,theom,edtmax,dhh,dg,aup,adw,           &
 603                        dv1,dv2,dv3,dv1q,dv2q,dv3q,dvq1,                        &
 604                        dv1u,dv2u,dv3u,dv1v,dv2v,dv3v,                          &
 605                        dellat,xdby,xqrch,xqc,xpw,xpwd,                         &
 606                        w1,w2,w3,w4,qrsk,evef,ptem,ptem1
 607 !
 608    logical         ::  totflg, cnvflg(its:ite),flg(its:ite)
 609 !
 610 !  climatological critical cloud work functions for closure
 611 !
 612    data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,               &
 613               350.,300.,250.,200.,150./
 614    data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216,                 &
 615               .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
 616 !
 617 !-----------------------------------------------------------------------
 618 !
 619 ! define miscellaneous values
 620 !
 621    pi_   = 3.14159
 622    qmin_ = 1.0e-30
 623    t0c_ = 273.15
 624    xlv0 = hvap_
 625    rcs  = 1.
 626    el2orc = hvap_*hvap_/(rv_*cp_)
 627    eps    = rd_/rv_
 628    fact1  = (cvap_-cliq_)/rv_
 629    fact2  = hvap_/rv_-fact1*t0c_
 630    kts1 = kts + 1
 631    kte1 = kte - 1
 632    dt2    = delt
 633    dtmin  = max(dt2,1200.)
 634    dtmax  = max(dt2,3600.)
 635 !
 636 !
 637 !  initialize arrays
 638 !
 639    do i = its,ite
 640      rain(i)     = 0.0
 641      kbot(i)   = kte+1
 642      ktop(i)   = 0
 643      kuo(i)    = 0
 644      cnvflg(i) = .true.
 645      dtconv(i) = 3600.
 646      pdot(i)   = 0.0
 647      edto(i)   = 0.0
 648      edtx(i)   = 0.0
 649      xmbmax(i) = 0.3
 650      excess(i) = 0.0
 651      plcl(i)   = 0.0
 652      kpbl(i)   = 1
 653      aa2(i) = 0.0
 654      qlko_ktcon(i) = 0.0
 655      pbcdif(i)= 0.0
 656      lmin(i) = 1
 657      jmin(i) = 1
 658      edt(i) = 0.0
 659    enddo
 660 !
 661    do k = 1,15
 662      acrit(k) = acritt(k) * (975. - pcrit(k))
 663    enddo
 664 !
 665 ! Define top layer for search of the downdraft originating layer
 666 ! and the maximum thetae for updraft
 667 !
 668    kbmax = kte
 669    kbm   = kte
 670    kmax  = kte
 671    do k = kts,kte
 672      do i = its,ite
 673        if(prsl(i,k).gt.prsi(i,1)*0.45) kbmax = k + 1
 674        if(prsl(i,k).gt.prsi(i,1)*0.70) kbm   = k + 1
 675        if(prsl(i,k).gt.prsi(i,1)*0.04) kmax  = k + 1
 676      enddo
 677    enddo
 678    kmax = min(kmax,kte)
 679    kmax1 = kmax - 1
 680    kbm = min(kbm,kte)
 681 !
 682 ! convert surface pressure to mb from cb
 683 !
 684    do k = kts,kte
 685      do i = its,ite
 686        pwo(i,k)  = 0.0
 687        pwdo(i,k) = 0.0
 688        dellal(i,k) = 0.0
 689        hcko(i,k) = 0.0
 690        qcko(i,k) = 0.0
 691        hcdo(i,k) = 0.0
 692        qcdo(i,k) = 0.0
 693      enddo
 694    enddo
 695 !
 696    do k = kts,kmax
 697      do i = its,ite
 698        p(i,k) = prsl(i,k) * 10.
 699        pwo(i,k) = 0.0
 700        pwdo(i,k) = 0.0
 701        to(i,k) = t1(i,k)
 702        qo(i,k) = q1(i,k)
 703        dbyo(i,k) = 0.0
 704        fent1(i,k) = 1.0
 705        fent2(i,k) = 1.0
 706        frh(i,k) = 0.0
 707        ucko(i,k) = 0.0
 708        vcko(i,k) = 0.0
 709        ucdo(i,k) = 0.0
 710        vcdo(i,k) = 0.0
 711        uo(i,k) = u1(i,k) * rcs
 712        vo(i,k) = v1(i,k) * rcs
 713      enddo
 714    enddo
 715 !
 716 ! column variables
 717 !
 718 !  p is pressure of the layer (mb)
 719 !  t is temperature at t-dt (k)..tn
 720 !  q is mixing ratio at t-dt (kg/kg)..qn
 721 !  to is temperature at t+dt (k)... this is after advection and turbulan
 722 !  qo is mixing ratio at t+dt (kg/kg)..q1
 723 !
 724    do k = kts,kmax
 725      do i = its,ite
 726        qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
 727        qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
 728        qeso(i,k) = max(qeso(i,k),qmin_)
 729        qo(i,k)   = max(qo(i,k), 1.e-10 )
 730        tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
 731      enddo
 732    enddo
 733 !
 734 ! compute moist static energy
 735 !
 736    do k = kts,kmax
 737      do i = its,ite
 738        heo(i,k)  = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qo(i,k)
 739        heso(i,k) = g_ * zl(i,k) + cp_* to(i,k) + hvap_ * qeso(i,k)
 740      enddo
 741    enddo
 742 !
 743 ! Determine level with largest moist static energy
 744 ! This is the level where updraft starts
 745 !
 746    do i = its,ite
 747      hmax(i) = heo(i,1)
 748      kb(i) = 1
 749    enddo
 750 !
 751    do k = kts1,kbm
 752      do i = its,ite
 753        if(heo(i,k).gt.hmax(i)) then
 754          kb(i) = k
 755          hmax(i) = heo(i,k)
 756        endif
 757      enddo
 758    enddo
 759 !
 760    do k = kts,kmax1
 761      do i = its,ite
 762        if(cnvflg(i)) then
 763          dz = .5 * (zl(i,k+1) - zl(i,k))
 764          dp = .5 * (p(i,k+1) - p(i,k))
 765          es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
 766          pprime = p(i,k+1) + (eps-1.) * es
 767          qs = eps * es / pprime
 768          dqsdp = - qs / pprime
 769          desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
 770          dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
 771          gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
 772          dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
 773          dq = dqsdt * dt + dqsdp * dp
 774          to(i,k) = to(i,k+1) + dt
 775          qo(i,k) = qo(i,k+1) + dq
 776          po = .5 * (p(i,k) + p(i,k+1))
 777          qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
 778          qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
 779          qeso(i,k) = max(qeso(i,k),qmin_)
 780          qo(i,k)   = max(qo(i,k), 1.e-10)
 781          frh(i,k)  = 1. - min(qo(i,k)/qeso(i,k), 1.)
 782          heo(i,k)  = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
 783                 cp_ * to(i,k) + hvap_ * qo(i,k)
 784          heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
 785                 cp_ * to(i,k) + hvap_ * qeso(i,k)
 786          uo(i,k)   = .5 * (uo(i,k) + uo(i,k+1))
 787          vo(i,k)   = .5 * (vo(i,k) + vo(i,k+1))
 788        endif
 789      enddo
 790    enddo
 791 !
 792 ! look for convective cloud base as the level of free convection
 793 !
 794    do i = its,ite
 795      if(cnvflg(i)) then
 796        indx = kb(i)
 797        hkbo(i) = heo(i,indx)
 798        qkbo(i) = qo(i,indx)
 799      endif
 800    enddo
 801 !
 802    do i = its,ite
 803      flg(i) = cnvflg(i)
 804      kbcon(i) = kmax
 805    enddo
 806 !
 807      do k = kts,kbmax
 808        do i = its,ite
 809          if(flg(i).and.k.gt.kb(i)) then
 810            hsbar(i) = heso(i,k)
 811            if(hkbo(i).gt.hsbar(i)) then
 812              flg(i) = .false.
 813              kbcon(i) = k
 814            endif
 815          endif
 816        enddo
 817      enddo
 818      do i = its,ite
 819        if(kbcon(i).eq.kmax) cnvflg(i) = .false.
 820      enddo
 821 !
 822      totflg = .true.
 823      do i = its,ite
 824        totflg = totflg .and. (.not. cnvflg(i))
 825      enddo
 826      if(totflg) return
 827 !
 828      do i = its,ite
 829        if(cnvflg(i)) then
 830 !
 831 !  determine critical convective inhibition
 832 !  as a function of vertical velocity at cloud base.
 833 !
 834           pdot(i)  = 10.* dot(i,kbcon(i))
 835           if(slimsk(i).eq.1.) then
 836             w1 = w1l
 837             w2 = w2l
 838             w3 = w3l
 839             w4 = w4l
 840           else
 841             w1 = w1s
 842             w2 = w2s
 843             w3 = w3s
 844             w4 = w4s
 845           endif
 846           if(pdot(i).le.w4) then
 847             tem = (pdot(i) - w4) / (w3 - w4)
 848           elseif(pdot(i).ge.-w4) then
 849             tem = - (pdot(i) + w4) / (w4 - w3)
 850           else
 851             tem = 0.
 852           endif
 853           tem = max(tem,-1.)
 854           tem = min(tem,1.)
 855           tem = 1. - tem
 856           tem1= .5*(cincrmax-cincrmin)
 857           cincr = cincrmax - tem * tem1
 858          pbcdif(i) = -p(i,kbcon(i)) + p(i,kb(i))
 859          if(pbcdif(i).gt.cincr) cnvflg(i) = .false.
 860        endif
 861      enddo
 862 !
 863 !
 864    totflg = .true.
 865    do i = its,ite
 866      totflg = totflg .and. (.not. cnvflg(i))
 867    enddo
 868    if(totflg) return
 869 !
 870    do k = kts,kte1
 871      do i = its,ite
 872        xlamb(i,k) = clam / zi(i,k+1)
 873      enddo
 874    enddo
 875 !
 876 !  assume that updraft entrainment rate above cloud base is
 877 !  same as that at cloud base
 878 !
 879    do k = kts1,kmax1
 880      do i = its,ite
 881        if(cnvflg(i).and.(k.gt.kbcon(i))) then
 882          xlamb(i,k) = xlamb(i,kbcon(i))
 883        endif
 884      enddo
 885    enddo
 886 !
 887 !   assume the detrainment rate for the updrafts to be same as
 888 !   the entrainment rate at cloud base
 889 !
 890    do i = its,ite
 891      if(cnvflg(i)) then
 892        xlamud(i) = xlamb(i,kbcon(i))
 893      endif
 894    enddo
 895 !
 896 !  functions rapidly decreasing with height, mimicking a cloud ensemble
 897 !    (Bechtold et al., 2008)
 898 !
 899    do k = kts1,kmax1
 900      do i = its,ite
 901        if(cnvflg(i).and.(k.gt.kbcon(i))) then
 902          tem = qeso(i,k)/qeso(i,kbcon(i))
 903          fent1(i,k) = tem**2
 904          fent2(i,k) = tem**3
 905        endif
 906      enddo
 907    enddo
 908 !
 909 !  final entrainment rate as the sum of turbulent part and organized entrainment
 910 !    depending on the environmental relative humidity
 911 !    (Bechtold et al., 2008)
 912 !
 913    do k = kts1,kmax1
 914      do i = its,ite
 915        if(cnvflg(i).and.(k.ge.kbcon(i))) then
 916           tem = cxlamu * frh(i,k) * fent2(i,k)
 917           xlamb(i,k) = xlamb(i,k)*fent1(i,k) + tem
 918        endif
 919      enddo
 920    enddo
 921 !
 922 ! determine updraft mass flux
 923 !
 924    do k = kts,kte
 925      do i = its,ite
 926       if(cnvflg(i)) then
 927          eta(i,k) = 1.
 928        endif
 929      enddo
 930    enddo
 931 !
 932    do k = kbmax,kts1,-1
 933      do i = its,ite
 934        if(cnvflg(i).and.k.lt.kbcon(i).and.k.ge.kb(i)) then
 935          dz = zi(i,k+2) - zi(i,k+1)
 936          ptem     = 0.5*(xlamb(i,k)+xlamb(i,k+1))-xlamud(i)
 937          eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
 938        endif
 939      enddo
 940    enddo
 941   do k = kts1,kmax1
 942      do i = its,ite
 943        if(cnvflg(i).and.k.gt.kbcon(i)) then
 944          dz  = zi(i,k+1) - zi(i,k)
 945          ptem = 0.5*(xlamb(i,k)+xlamb(i,k-1))-xlamud(i)
 946          eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
 947        endif
 948      enddo
 949    enddo
 950    do i = its,ite
 951      if(cnvflg(i)) then
 952        dz = zi(i,3) - zi(i,2)
 953        ptem     = 0.5*(xlamb(i,1)+xlamb(i,2))-xlamud(i)
 954        eta(i,1) = eta(i,2) / (1. + ptem * dz)
 955      endif
 956    enddo
 957 !
 958 ! work up updraft cloud properties
 959 !
 960    do i = its,ite
 961      if(cnvflg(i)) then
 962        indx = kb(i)
 963        hcko(i,indx) = hkbo(i)
 964        qcko(i,indx) = qkbo(i)
 965        ucko(i,indx) = uo(i,indx)
 966        vcko(i,indx) = vo(i,indx)
 967        pwavo(i) = 0.
 968      endif
 969    enddo
 970 !
 971 ! cloud property below cloud base is modified by the entrainment proces
 972 !
 973    do k = kts1,kmax1
 974      do i = its,ite
 975        if(cnvflg(i).and.k.gt.kb(i)) then
 976          dz   = zi(i,k+1) - zi(i,k)
 977          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
 978          tem1 = 0.5 * xlamud(i) * dz
 979          factor = 1. + tem - tem1
 980          ptem = 0.5 * tem + pgcon
 981          ptem1= 0.5 * tem - pgcon
 982          hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                          &
 983                      (heo(i,k)+heo(i,k-1)))/factor
 984          ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)                      &
 985                      +ptem1*uo(i,k-1))/factor
 986          vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)                      &
 987                      +ptem1*vo(i,k-1))/factor
 988          dbyo(i,k) = hcko(i,k) - heso(i,k)
 989        endif
 990      enddo
 991    enddo
 992 !
 993 !   taking account into convection inhibition due to existence of
 994 !    dry layers below cloud base
 995 !
 996    do i = its,ite
 997      flg(i) = cnvflg(i)
 998      kbcon1(i) = kmax
 999    enddo
1000 !
1001    do k = kts1,kmax
1002      do i = its,ite
1003        if(flg(i).and.k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
1004          kbcon1(i) = k
1005          flg(i) = .false.
1006        endif
1007      enddo
1008    enddo
1009 !
1010    do i = its,ite
1011      if(cnvflg(i)) then
1012        if(kbcon1(i).eq.kmax) cnvflg(i) = .false.
1013      endif
1014    enddo
1015 !
1016    do i =its,ite
1017      if(cnvflg(i)) then
1018        tem = p(i,kbcon(i)) - p(i,kbcon1(i))
1019        if(tem.gt.dthk) then
1020           cnvflg(i) = .false.
1021        endif
1022      endif
1023    enddo
1024 !
1025    totflg = .true.
1026    do i = its,ite
1027      totflg = totflg .and. (.not. cnvflg(i))
1028    enddo
1029    if(totflg) return
1030 !
1031 !
1032 !  determine cloud top
1033 !
1034 !
1035    do i = its,ite
1036      flg(i) = cnvflg(i)
1037      ktcon(i) = 1
1038    enddo
1039 !
1040 !   check inversion
1041 !
1042    do k = kts1,kmax1
1043      do i = its,ite
1044        if(dbyo(i,k).lt.0..and.flg(i).and.k.gt. kbcon1(i)) then
1045          ktcon(i) = k
1046          flg(i)   = .false.
1047        endif
1048      enddo
1049    enddo
1050 !
1051 !
1052 ! check cloud depth
1053 !
1054    do i = its,ite
1055      if(cnvflg(i).and.(p(i,kbcon(i)) - p(i,ktcon(i))).lt.150.)                 &
1056             cnvflg(i) = .false.
1057    enddo
1058 !
1059    totflg = .true.
1060    do i = its,ite
1061      totflg = totflg .and. (.not. cnvflg(i))
1062    enddo
1063    if(totflg) return
1064 !
1065 !
1066 !  search for downdraft originating level above theta-e minimum
1067 !
1068    do i = its,ite
1069      if(cnvflg(i)) then
1070        hmin(i) = heo(i,kbcon1(i))
1071        lmin(i) = kbmax
1072        jmin(i) = kbmax
1073     endif
1074    enddo
1075 !
1076    do k = kts1,kbmax
1077      do i = its,ite
1078        if(cnvflg(i).and.k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
1079          lmin(i) = k + 1
1080          hmin(i) = heo(i,k)
1081        endif
1082      enddo
1083    enddo
1084 !
1085 ! make sure that jmin is within the cloud
1086 !
1087    do i = its,ite
1088      if(cnvflg(i)) then
1089        jmin(i) = min(lmin(i),ktcon(i)-1)
1090        jmin(i) = max(jmin(i),kbcon1(i)+1)
1091        if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
1092      endif
1093    enddo
1094 !
1095 !  specify upper limit of mass flux at cloud base
1096 !
1097    do i = its,ite
1098      if(cnvflg(i)) then
1099        k = kbcon(i)
1100        dp = 1000. * del(i,k)
1101        xmbmax(i) = dp / (g_ * dt2)
1102      endif
1103    enddo
1104 !
1105 !
1106 ! compute cloud moisture property and precipitation
1107 !
1108    do i = its,ite
1109      aa1(i) = 0.
1110    enddo
1111 !
1112    do k = kts1,kmax
1113      do i = its,ite
1114        if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1115          dz = .5 * (zl(i,k+1) - zl(i,k-1))
1116          dz1 = (zi(i,k+1) - zi(i,k))
1117          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1118          qrch = qeso(i,k)                                                      &
1119               + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1120          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz1
1121          tem1 = 0.5 * xlamud(i) * dz1
1122          factor = 1. + tem - tem1
1123          qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                           &
1124                     (qo(i,k)+qo(i,k-1)))/factor
1125          qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1126 !
1127 ! check if there is excess moisture to release latent heat
1128 !
1129          if(qcirs(i,k).gt.0. .and. k.ge.kbcon(i)) then
1130            etah = .5 * (eta(i,k) + eta(i,k-1))
1131            if(ncloud.gt.0..and.k.gt.jmin(i)) then
1132              dp = 1000. * del(i,k)
1133              qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz1)
1134              dellal(i,k) = etah * c1 * dz1 * qlk * g_ / dp
1135            else
1136              qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz1)
1137            endif
1138            aa1(i) = aa1(i) - dz1 * g_ * qlk
1139            qc = qlk + qrch
1140            pwo(i,k) = etah * c0 * dz1 * qlk
1141            qcko(i,k) = qc
1142            pwavo(i) = pwavo(i) + pwo(i,k)
1143          endif
1144        endif
1145      enddo
1146    enddo
1147 !
1148 ! calculate cloud work function at t+dt
1149 !
1150    do k = kts1,kmax
1151      do i = its,ite
1152        if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon(i)) then
1153          dz1 = zl(i,k+1) - zl(i,k)
1154          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1155          rfact =  1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1156          aa1(i) = aa1(i) +dz1 * (g_ / (cp_ * to(i,k)))                         &
1157                   * dbyo(i,k) / (1. + gamma)* rfact
1158          aa1(i) = aa1(i)+dz1 * g_ * fv_ *                                      &
1159                   max(0.,(qeso(i,k) - qo(i,k)))
1160        endif
1161      enddo
1162    enddo
1163 !
1164    do i = its,ite
1165      if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1166    enddo
1167 !
1168       totflg = .true.
1169       do i=its,ite
1170         totflg = totflg .and. (.not. cnvflg(i))
1171       enddo
1172       if(totflg) return
1173 !
1174 !    estimate the convective overshooting as the level
1175 !    where the [aafac * cloud work function] becomes zero,
1176 !    which is the final cloud top
1177 !
1178       do i = its,ite
1179         if (cnvflg(i)) then
1180           aa2(i) = aafac * aa1(i)
1181         endif
1182       enddo
1183 !
1184       do i = its,ite
1185         flg(i) = cnvflg(i)
1186         ktcon1(i) = kmax1
1187       enddo
1188 !
1189       do k = kts1, kmax
1190         do i = its, ite
1191           if (flg(i)) then
1192             if(k.ge.ktcon(i).and.k.lt.kmax) then
1193               dz1 = zl(i,k+1) - zl(i,k)
1194               gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1195               rfact =  1. + fv_ * cp_ * gamma* to(i,k) / hvap_
1196               aa2(i) = aa2(i) +dz1 * (g_ / (cp_ * to(i,k)))                    &
1197                        * dbyo(i,k) / (1. + gamma)* rfact
1198               if(aa2(i).lt.0.) then
1199                 ktcon1(i) = k
1200                 flg(i) = .false.
1201               endif
1202             endif
1203           endif
1204         enddo
1205       enddo
1206 !
1207 !  compute cloud moisture property, detraining cloud water
1208 !  and precipitation in overshooting layers
1209 !
1210    do k = kts1,kmax
1211      do i = its,ite
1212        if (cnvflg(i)) then
1213          if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1214            dz = (zi(i,k+1) - zi(i,k))
1215            gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1216            qrch = qeso(i,k)+ gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1217            tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1218            tem1 = 0.5 * xlamud(i) * dz
1219            factor = 1. + tem - tem1
1220            qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                         &
1221                       (qo(i,k)+qo(i,k-1)))/factor
1222            qcirs(i,k) = eta(i,k) * qcko(i,k) - eta(i,k) * qrch
1223 !
1224 !  check if there is excess moisture to release latent heat
1225 !
1226            if(qcirs(i,k).gt.0.) then
1227              etah = .5 * (eta(i,k) + eta(i,k-1))
1228              if(ncloud.gt.0.) then
1229                dp = 1000. * del(i,k)
1230                qlk = qcirs(i,k) / (eta(i,k) + etah * (c0 + c1) * dz)
1231                dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
1232              else
1233                qlk = qcirs(i,k) / (eta(i,k) + etah * c0 * dz)
1234              endif
1235              qc = qlk + qrch
1236              pwo(i,k) = etah * c0 * dz * qlk
1237              qcko(i,k) = qc
1238              pwavo(i) = pwavo(i) + pwo(i,k)
1239            endif
1240          endif
1241        endif
1242      enddo
1243    enddo
1244 !
1245 ! exchange ktcon with ktcon1
1246 !
1247    do i = its,ite
1248      if(cnvflg(i)) then
1249        kk = ktcon(i)
1250        ktcon(i) = ktcon1(i)
1251        ktcon1(i) = kk
1252      endif
1253    enddo
1254 !
1255 ! this section is ready for cloud water
1256 !
1257    if (ncloud.gt.0) then
1258 !
1259 !  compute liquid and vapor separation at cloud top
1260 !
1261    do i = its,ite
1262      if(cnvflg(i)) then
1263      k = ktcon(i)-1
1264        gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1265        qrch = qeso(i,k)                                                     &
1266               + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1267        dq = qcko(i,k) - qrch
1268 !
1269 !  check if there is excess moisture to release latent heat
1270 !
1271        if(dq.gt.0.) then
1272          qlko_ktcon(i) = dq
1273          qcko(i,k) = qrch
1274        endif
1275      endif
1276    enddo
1277    endif
1278 !
1279 ! ..... downdraft calculations .....
1280 !
1281 ! determine downdraft strength in terms of wind shear
1282 !
1283    do i = its,ite
1284      if(cnvflg(i)) then
1285        vshear(i) = 0.
1286      endif
1287    enddo
1288 !
1289    do k = kts1,kmax
1290      do i = its,ite
1291        if(k.gt.kb(i).and.k.le.ktcon(i).and.cnvflg(i)) then
1292          shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2                                  &
1293                        + (vo(i,k)-vo(i,k-1)) ** 2)
1294          vshear(i) = vshear(i) + shear
1295        endif
1296      enddo
1297    enddo
1298 !
1299    do i = its,ite
1300      if(cnvflg(i)) then
1301        vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1302        e1 = 1.591-.639*vshear(i)                                               &
1303            +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1304        edt(i)  = 1.-e1
1305        edt(i)  = min(edt(i),.9)
1306        edt(i)  = max(edt(i),.0)
1307        edto(i) = edt(i)
1308        edtx(i) = edt(i)
1309      endif
1310    enddo
1311 !
1312 ! determine detrainment rate between 1 and kbdtr
1313 !
1314    do i = its,ite
1315      if(cnvflg(i)) then
1316        sumx(i) = 0.
1317      endif
1318    enddo
1319 !
1320    do k = kts,kmax1
1321      do i = its,ite
1322        if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
1323          dz = zi(i,k+2) - zi(i,k+1)
1324          sumx(i) = sumx(i) + dz
1325        endif
1326      enddo
1327    enddo
1328 !
1329    do i = its,ite
1330      kbdtr(i) = kbcon(i)
1331      beta = betas
1332      if(slimsk(i).eq.1.) beta = betal
1333      if(cnvflg(i)) then
1334        kbdtr(i) = kbcon(i)
1335        kbdtr(i) = max(kbdtr(i),1)
1336        dz =(sumx(i)+zi(i,2))/float(kbcon(i))
1337        tem = 1./float(kbcon(i))
1338        xlamd(i) = (1.-beta**tem)/dz
1339      endif
1340    enddo
1341 !
1342 ! determine downdraft mass flux
1343 !
1344    do k = kts,kmax
1345      do i = its,ite
1346        if(cnvflg(i)) then
1347          etad(i,k) = 1.
1348        endif
1349        qrcdo(i,k) = 0.
1350        qrcd(i,k) = 0.
1351      enddo
1352    enddo
1353 !
1354    do k = kmax1,kts,-1
1355      do i = its,ite
1356        if(cnvflg(i)) then
1357          if(k.lt.jmin(i).and.k.ge.kbcon(i))then
1358            dz = (zi(i,k+2) - zi(i,k+1))
1359            ptem = xlamdd-xlamde
1360            etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1361          elseif(k.lt.kbcon(i))then
1362            dz = (zi(i,k+2) - zi(i,k+1))
1363            ptem = xlamd(i)+xlamdd-xlamde
1364            etad(i,k) = etad(i,k+1) * (1.-ptem * dz)
1365          endif
1366        endif
1367      enddo
1368    enddo
1369 !
1370 !
1371 ! downdraft moisture properties
1372 !
1373    do i = its,ite
1374      if(cnvflg(i)) then
1375       pwevo(i) = 0.
1376      endif
1377    enddo
1378 !
1379    do i = its,ite
1380      if(cnvflg(i))  then
1381        jmn = jmin(i)
1382        hcdo(i,jmn) = heo(i,jmn)
1383        qcdo(i,jmn) = qo(i,jmn)
1384        qrcdo(i,jmn) = qeso(i,jmn)
1385        ucdo(i,jmn) = uo(i,jmn)
1386        vcdo(i,jmn) = vo(i,jmn)
1387      endif
1388    enddo
1389 !
1390    do k = kmax1,kts,-1
1391      do i = its,ite
1392        if (cnvflg(i) .and. k.lt.jmin(i)) then
1393          dz = zi(i,k+2) - zi(i,k+1)
1394          if(k.ge.kbcon(i)) then
1395            tem  = xlamde * dz
1396            tem1 = 0.5 * xlamdd * dz
1397          else
1398            tem  = xlamde * dz
1399            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1400           endif
1401           factor = 1. + tem - tem1
1402           ptem = 0.5 * tem - pgcon
1403           ptem1= 0.5 * tem + pgcon
1404           hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5*                     &
1405                       (heo(i,k)+heo(i,k+1)))/factor
1406           ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1)               &
1407                      +ptem1*uo(i,k))/factor
1408           vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1)               &
1409                      +ptem1*vo(i,k))/factor
1410           dbyo(i,k) = hcdo(i,k) - heso(i,k)
1411        endif
1412      enddo
1413    enddo
1414 !
1415    do k = kmax1,kts,-1
1416      do i = its,ite
1417        if(cnvflg(i).and.k.lt.jmin(i)) then
1418          dq = qeso(i,k)
1419          dt = to(i,k)
1420          gamma = el2orc * dq / dt**2
1421          qrcdo(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dbyo(i,k)
1422          detad = etad(i,k+1) - etad(i,k)
1423          dz = zi(i,k+2) - zi(i,k+1)
1424          if(k.ge.kbcon(i)) then
1425             tem  = xlamde * dz
1426             tem1 = 0.5 * xlamdd * dz
1427          else
1428             tem  = xlamde * dz
1429             tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1430          endif
1431          factor = 1. + tem - tem1
1432          qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5*                      &
1433                      (qo(i,k)+qo(i,k+1)))/factor
1434          pwdo(i,k) = etad(i,k+1) * qcdo(i,k) -etad(i,k+1) * qrcdo(i,k)
1435          qcdo(i,k) = qrcdo(i,k)
1436          pwevo(i) = pwevo(i) + pwdo(i,k)
1437        endif
1438      enddo
1439    enddo
1440 !
1441 ! final downdraft strength dependent on precip
1442 ! efficiency (edt), normalized condensate (pwav), and
1443 ! evaporate (pwev)
1444 !
1445    do i = its,ite
1446      edtmax = edtmaxl
1447      if(slimsk(i).eq.2.) edtmax = edtmaxs
1448      if(cnvflg(i)) then
1449        if(pwevo(i).lt.0.) then
1450          edto(i) = -edto(i) * pwavo(i) / pwevo(i)
1451          edto(i) = min(edto(i),edtmax)
1452        else
1453          edto(i) = 0.
1454        endif
1455      endif
1456    enddo
1457 !
1458 ! downdraft cloudwork functions
1459 !
1460    do k = kmax1,kts,-1
1461      do i = its,ite
1462        if(cnvflg(i).and.k.lt.jmin(i)) then
1463          gamma = el2orc * qeso(i,k) / to(i,k)**2
1464          dhh=hcdo(i,k)
1465          dt=to(i,k)
1466          dg=gamma
1467          dh=heso(i,k)
1468          dz=-1.*(zl(i,k+1)-zl(i,k))
1469          aa1(i)=aa1(i)+edto(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))             &
1470                 *(1.+fv_*cp_*dg*dt/hvap_)
1471          aa1(i)=aa1(i)+edto(i)*dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1472        endif
1473      enddo
1474    enddo
1475 !
1476    do i = its,ite
1477      if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
1478    enddo
1479 !
1480    totflg = .true.
1481    do i=its,ite
1482      totflg = totflg .and. (.not. cnvflg(i))
1483    enddo
1484    if(totflg) return
1485 !
1486 ! what would the change be, that a cloud with unit mass
1487 ! will do to the environment?
1488 !
1489    do k = kts,kmax
1490      do i = its,ite
1491        if(cnvflg(i)) then
1492          dellah(i,k) = 0.
1493          dellaq(i,k) = 0.
1494          dellau(i,k) = 0.
1495          dellav(i,k) = 0.
1496        endif
1497      enddo
1498    enddo
1499 !
1500    do i = its,ite
1501      if(cnvflg(i)) then
1502        dp = 1000. * del(i,1)
1503        dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1)                          &
1504                    - heo(i,1)) * g_ / dp
1505        dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1)                          &
1506                    - qo(i,1)) * g_ / dp
1507        dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1)                          &
1508                    - uo(i,1)) * g_ / dp
1509        dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1)                          &
1510                    - vo(i,1)) * g_ / dp
1511      endif
1512    enddo
1513 !
1514 ! changed due to subsidence and entrainment
1515 !
1516    do k = kts1,kmax1
1517      do i = its,ite
1518        if(cnvflg(i).and.k.lt.ktcon(i)) then
1519          aup = 1.
1520          if(k.le.kb(i)) aup = 0.
1521          adw = 1.
1522          if(k.gt.jmin(i)) adw = 0.
1523          dv1= heo(i,k)
1524          dv2 = .5 * (heo(i,k) + heo(i,k-1))
1525          dv3= heo(i,k-1)
1526          dv1q= qo(i,k)
1527          dv2q = .5 * (qo(i,k) + qo(i,k-1))
1528          dv3q= qo(i,k-1)
1529          dv1u = uo(i,k)
1530          dv2u = .5 * (uo(i,k) + uo(i,k-1))
1531          dv3u = uo(i,k-1)
1532          dv1v = vo(i,k)
1533          dv2v = .5 * (vo(i,k) + vo(i,k-1))
1534          dv3v = vo(i,k-1)
1535          dp = 1000. * del(i,k)
1536          dz = zi(i,k+1) - zi(i,k)
1537          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1))
1538          tem1 = xlamud(i)
1539          if(k.le.kbcon(i)) then
1540            ptem  = xlamde
1541            ptem1 = xlamd(i)+xlamdd
1542          else
1543            ptem  = xlamde
1544            ptem1 = xlamdd
1545          endif
1546          deta = eta(i,k) - eta(i,k-1)
1547          detad = etad(i,k) - etad(i,k-1)
1548          dellah(i,k) = dellah(i,k) +                                           &
1549              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1               &
1550          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3               &
1551          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2*dz              &
1552          +  aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz                  &
1553          +  adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1))*dz) *g_/dp
1554          dellaq(i,k) = dellaq(i,k) +                                           &
1555              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1q              &
1556          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3q              &
1557          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz             &
1558          +  aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz                  &
1559          +  adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1))*dz) *g_/dp
1560          dellau(i,k) = dellau(i,k) +                                           &
1561              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1u              &
1562          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3u              &
1563          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz             &
1564          +  aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz                  &
1565          +  adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz          &
1566          -  pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u))*g_/dp
1567 !
1568          dellav(i,k) = dellav(i,k) +                                           &
1569              ((aup * eta(i,k) - adw * edto(i) * etad(i,k)) * dv1v              &
1570          - (aup * eta(i,k-1) - adw * edto(i) * etad(i,k-1))* dv3v              &
1571          - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz             &
1572          +  aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz                  &
1573          +  adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz          &
1574          -  pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v))*g_/dp
1575        endif
1576      enddo
1577    enddo
1578 !
1579 ! cloud top
1580 !
1581    do i = its,ite
1582      if(cnvflg(i)) then
1583        indx = ktcon(i)
1584        dp = 1000. * del(i,indx)
1585        dv1 = heo(i,indx-1)
1586        dellah(i,indx) = eta(i,indx-1) *                                        &
1587                         (hcko(i,indx-1) - dv1) * g_ / dp
1588        dvq1 = qo(i,indx-1)
1589        dellaq(i,indx) = eta(i,indx-1) *                                        &
1590                         (qcko(i,indx-1) - dvq1) * g_ / dp
1591        dv1u = uo(i,indx-1)
1592        dellau(i,indx) = eta(i,indx-1) *                                        &
1593                         (ucko(i,indx-1) - dv1u) * g_ / dp
1594        dv1v = vo(i,indx-1)
1595        dellav(i,indx) = eta(i,indx-1) *                                        &
1596                         (vcko(i,indx-1) - dv1v) * g_ / dp
1597 !
1598 !  cloud water
1599 !
1600        dellal(i,indx) = eta(i,indx-1) * qlko_ktcon(i) * g_ / dp
1601      endif
1602    enddo
1603 !
1604 ! final changed variable per unit mass flux
1605 !
1606    do k = kts,kmax
1607      do i = its,ite
1608        if(cnvflg(i).and.k.gt.ktcon(i)) then
1609          qo(i,k) = q1(i,k)
1610          to(i,k) = t1(i,k)
1611        endif
1612        if(cnvflg(i).and.k.le.ktcon(i)) then
1613          qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
1614          dellat  = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1615          to(i,k) = dellat * mbdt + t1(i,k)
1616          qo(i,k) = max(qo(i,k),1.0e-10)
1617        endif
1618      enddo
1619    enddo
1620 !
1621 !------------------------------------------------------------------------------
1622 !
1623 ! the above changed environment is now used to calulate the
1624 ! effect the arbitrary cloud (with unit mass flux)
1625 ! which then is used to calculate the real mass flux,
1626 ! necessary to keep this change in balance with the large-scale
1627 ! destabilization.
1628 !
1629 ! environmental conditions again
1630 !
1631 !------------------------------------------------------------------------------
1632 !
1633    do k = kts,kmax
1634      do i = its,ite
1635        if(cnvflg(i)) then
1636          qeso(i,k)=0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1637          qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1638          qeso(i,k) = max(qeso(i,k),qmin_)
1639          tvo(i,k)  = to(i,k) + fv_ * to(i,k) * max(qo(i,k),qmin_)
1640        endif
1641      enddo
1642    enddo
1643 !
1644    do i = its,ite
1645      if(cnvflg(i)) then
1646        xaa0(i) = 0.
1647        xpwav(i) = 0.
1648      endif
1649    enddo
1650 !
1651 ! moist static energy
1652 !
1653    do k = kts,kmax1
1654      do i = its,ite
1655        if(cnvflg(i)) then
1656          dz = .5 * (zl(i,k+1) - zl(i,k))
1657          dp = .5 * (p(i,k+1) - p(i,k))
1658          es =0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1659          pprime = p(i,k+1) + (eps-1.) * es
1660          qs = eps * es / pprime
1661          dqsdp = - qs / pprime
1662          desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
1663          dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
1664          gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
1665          dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
1666          dq = dqsdt * dt + dqsdp * dp
1667          to(i,k) = to(i,k+1) + dt
1668          qo(i,k) = qo(i,k+1) + dq
1669          po = .5 * (p(i,k) + p(i,k+1))
1670          qeso(i,k) =0.01* fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1671          qeso(i,k) = eps * qeso(i,k) / (po + (eps-1.) * qeso(i,k))
1672          qeso(i,k) = max(qeso(i,k),qmin_)
1673          qo(i,k)   = max(qo(i,k), 1.0e-10)
1674          heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                          &
1675                      cp_ * to(i,k) + hvap_ * qo(i,k)
1676          heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                         &
1677                      cp_ * to(i,k) + hvap_ * qeso(i,k)
1678        endif
1679      enddo
1680    enddo
1681 !
1682    k = kmax
1683    do i = its,ite
1684      if(cnvflg(i)) then
1685        heo(i,k)  = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qo(i,k)
1686        heso(i,k) = g_ * zl(i,k) + cp_ * to(i,k) + hvap_ * qeso(i,k)
1687      endif
1688    enddo
1689 !
1690    do i = its,ite
1691      if(cnvflg(i)) then
1692        xaa0(i) = 0.
1693        xpwav(i) = 0.
1694        indx = kb(i)
1695        xhkb(i) = heo(i,indx)
1696        xqkb(i) = qo(i,indx)
1697        hcko(i,indx) = xhkb(i)
1698        qcko(i,indx) = xqkb(i)
1699      endif
1700    enddo
1701 !
1702 ! ..... static control .....
1703 !
1704 ! moisture and cloud work functions
1705 !
1706    do k = kts1,kmax1
1707      do i = its,ite
1708        if(cnvflg(i).and.k.gt.kb(i).and.k.le.ktcon(i)) then
1709          dz = zi(i,k+1) - zi(i,k)
1710          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1711          tem1 = 0.5 * xlamud(i) * dz
1712          factor = 1. + tem - tem1
1713          hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                         &
1714                     (heo(i,k)+heo(i,k-1)))/factor
1715        endif
1716      enddo
1717    enddo
1718 !
1719    do k = kts1,kmax1
1720      do i = its,ite
1721        if(cnvflg(i).and.k.gt.kb(i).and.k.lt.ktcon(i)) then
1722          dz = zi(i,k+1) - zi(i,k)
1723          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1724          xdby = hcko(i,k) - heso(i,k)
1725          xqrch = qeso(i,k)                                                     &
1726               + gamma * xdby / (hvap_ * (1. + gamma))
1727          tem  = 0.5 * (xlamb(i,k)+xlamb(i,k-1)) * dz
1728          tem1 = 0.5 * xlamud(i) * dz
1729          factor = 1. + tem - tem1
1730          qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*(qo(i,k)+qo(i,k-1)))/factor
1731          dq = eta(i,k) * qcko(i,k) - eta(i,k) * xqrch
1732          if(k.ge.kbcon(i).and.dq.gt.0.) then
1733            etah = .5 * (eta(i,k) + eta(i,k-1))
1734            if(ncloud.gt.0..and.k.gt.jmin(i)) then
1735              qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1736            else
1737              qlk = dq / (eta(i,k) + etah * c0 * dz)
1738            endif
1739            if(k.lt.ktcon1(i)) then
1740              xaa0(i) = xaa0(i) - dz * g_ * qlk
1741            endif
1742            qcko(i,k) = qlk + xqrch
1743            xpw = etah * c0 * dz * qlk
1744            xpwav(i) = xpwav(i) + xpw
1745          endif
1746        endif
1747        if(cnvflg(i).and.k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
1748          dz1 = zl(i,k+1) - zl(i,k)
1749          gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1750          rfact =  1. + fv_ * cp_ * gamma                                       &
1751                   * to(i,k) / hvap_
1752          xdby = hcko(i,k) - heso(i,k)
1753          xaa0(i) = xaa0(i)                                                     &
1754                  + dz1 * (g_ / (cp_ * to(i,k)))                                &
1755                  * xdby / (1. + gamma)                                         &
1756                  * rfact
1757          xaa0(i)=xaa0(i)+                                                      &
1758                   dz1 * g_ * fv_ *                                             &
1759                   max(0.,(qeso(i,k) - qo(i,k)))
1760        endif
1761      enddo
1762    enddo
1763 !
1764 ! ..... downdraft calculations .....
1765 !
1766 !
1767 ! downdraft moisture properties
1768 !
1769    do i = its,ite
1770      xpwev(i) = 0.
1771    enddo
1772    do i = its,ite
1773      if(cnvflg(i)) then
1774        jmn = jmin(i)
1775        xhcd(i,jmn) = heo(i,jmn)
1776        xqcd(i,jmn) = qo(i,jmn)
1777        qrcd(i,jmn) = qeso(i,jmn)
1778      endif
1779    enddo
1780 !
1781    do k = kmax1,kts, -1
1782      do i = its,ite
1783        if(cnvflg(i).and.k.lt.jmin(i)) then
1784          dz = zi(i,k+2) - zi(i,k+1)
1785          if(k.ge.kbcon(i)) then
1786             tem  = xlamde * dz
1787             tem1 = 0.5 * xlamdd * dz
1788          else
1789             tem  = xlamde * dz
1790             tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1791          endif
1792          factor = 1. + tem - tem1
1793          xhcd(i,k) = ((1.-tem1)*xhcd(i,k+1)+tem*0.5*                        &
1794                     (heo(i,k)+heo(i,k+1)))/factor
1795        endif
1796      enddo
1797    enddo
1798 !
1799    do k = kmax1,kts, -1
1800      do i = its,ite
1801        if(cnvflg(i).and.k.lt.jmin(i)) then
1802          dq = qeso(i,k)
1803          dt = to(i,k)
1804          gamma = el2orc * dq / dt**2
1805          dh = xhcd(i,k) - heso(i,k)
1806          qrcd(i,k)=dq+(1./hvap_)*(gamma/(1.+gamma))*dh
1807          dz = zi(i,k+2) - zi(i,k+1)
1808          if(k.ge.kbcon(i)) then
1809            tem  = xlamde * dz
1810            tem1 = 0.5 * xlamdd * dz
1811          else
1812            tem  = xlamde * dz
1813            tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
1814          endif
1815          factor = 1. + tem - tem1
1816          xqcd(i,k) = ((1.-tem1)*xqcd(i,k+1)+tem*0.5*                           &
1817                    (qo(i,k)+qo(i,k+1)))/factor
1818          xpwd     = etad(i,k+1) * (xqcd(i,k) - qrcd(i,k))
1819          xqcd(i,k)= qrcd(i,k)
1820          xpwev(i) = xpwev(i) + xpwd
1821        endif
1822      enddo
1823    enddo
1824 !
1825    do i = its,ite
1826      edtmax = edtmaxl
1827      if(slimsk(i).eq.2.) edtmax = edtmaxs
1828      if(cnvflg(i)) then
1829        if(xpwev(i).ge.0.) then
1830          edtx(i) = 0.
1831        else
1832          edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
1833          edtx(i) = min(edtx(i),edtmax)
1834        endif
1835      endif
1836    enddo
1837 !
1838 ! downdraft cloudwork functions
1839 !
1840    do k = kmax1,kts, -1
1841      do i = its,ite
1842        if(cnvflg(i).and.k.lt.jmin(i)) then
1843          gamma = el2orc * qeso(i,k) / to(i,k)**2
1844          dhh=xhcd(i,k)
1845          dt= to(i,k)
1846          dg= gamma
1847          dh= heso(i,k)
1848          dz=-1.*(zl(i,k+1)-zl(i,k))
1849          xaa0(i)=xaa0(i)+edtx(i)*dz*(g_/(cp_*dt))*((dhh-dh)/(1.+dg))           &
1850                  *(1.+fv_*cp_*dg*dt/hvap_)
1851          xaa0(i)=xaa0(i)+edtx(i)*                                              &
1852                   dz*g_*fv_*max(0.,(qeso(i,k)-qo(i,k)))
1853        endif
1854      enddo
1855    enddo
1856 !
1857 ! calculate critical cloud work function
1858 !
1859    do i = its,ite
1860      acrt(i) = 0.
1861      if(cnvflg(i)) then
1862        if(p(i,ktcon(i)).lt.pcrit(15))then
1863          acrt(i)=acrit(15)*(975.-p(i,ktcon(i)))/(975.-pcrit(15))
1864        else if(p(i,ktcon(i)).gt.pcrit(1))then
1865          acrt(i)=acrit(1)
1866        else
1867          k = int((850. - p(i,ktcon(i)))/50.) + 2
1868          k = min(k,15)
1869          k = max(k,2)
1870          acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))*                               &
1871               (p(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
1872         endif
1873       endif
1874     enddo
1875 !
1876    do i = its,ite
1877      acrtfct(i) = 1.
1878      w1 = w1s
1879      w2 = w2s
1880      w3 = w3s
1881      w4 = w4s
1882      if(slimsk(i).eq.1.) then
1883        w1 = w1l
1884        w2 = w2l
1885        w3 = w3l
1886        w4 = w4l
1887      endif
1888      if(cnvflg(i)) then
1889        if(pdot(i).le.w4) then
1890          acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
1891        elseif(pdot(i).ge.-w4) then
1892        acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
1893        else
1894          acrtfct(i) = 0.
1895        endif
1896        acrtfct(i) = max(acrtfct(i),-1.)
1897        acrtfct(i) = min(acrtfct(i),1.)
1898        acrtfct(i) = 1. - acrtfct(i)
1899        dtconv(i) = dt2 + max((1800. - dt2),0.) * (pdot(i) - w2) / (w1 - w2)
1900        dtconv(i) = max(dtconv(i),dtmin)
1901        dtconv(i) = min(dtconv(i),dtmax)
1902 !
1903      endif
1904    enddo
1905 !
1906 ! large scale forcing
1907 !
1908    do i= its,ite
1909      if(cnvflg(i)) then
1910        f(i) = (aa1(i) - acrt(i) * acrtfct(i)) / dtconv(i)
1911        if(f(i).le.0.) cnvflg(i) = .false.
1912      endif
1913      if(cnvflg(i)) then
1914        xk(i) = (xaa0(i) - aa1(i)) / mbdt
1915        if(xk(i).ge.0.) cnvflg(i) = .false.
1916      endif
1917 !
1918 ! kernel, cloud base mass flux
1919 !
1920      if(cnvflg(i)) then
1921        xmb(i) = -f(i) / xk(i)
1922        xmb(i) = min(xmb(i),xmbmax(i))
1923      endif
1924 !
1925      if(cnvflg(i)) then
1926      endif
1927 !
1928    enddo
1929    totflg = .true.
1930    do i = its,ite
1931      totflg = totflg .and. (.not. cnvflg(i))
1932    enddo
1933    if(totflg) return
1934 !
1935 ! restore t0 and qo to t1 and q1 in case convection stops
1936 !
1937    do k = kts,kmax
1938      do i = its,ite
1939        if (cnvflg(i)) then
1940        to(i,k) = t1(i,k)
1941        qo(i,k) = q1(i,k)
1942        uo(i,k) = u1(i,k)
1943        vo(i,k) = v1(i,k)
1944        qeso(i,k) = 0.01*fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1945        qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.) * qeso(i,k))
1946        qeso(i,k) = max(qeso(i,k),qmin_)
1947        endif
1948      enddo
1949    enddo
1950 !
1951 ! feedback: simply the changes from the cloud with unit mass flux
1952 !           multiplied by  the mass flux necessary to keep the
1953 !           equilibrium with the larger-scale.
1954 !
1955    do i = its,ite
1956      delhbar(i) = 0.
1957      delqbar(i) = 0.
1958      deltbar(i) = 0.
1959      qcond(i) = 0.
1960      qrski(i) = 0.
1961      delubar(i) = 0.
1962      delvbar(i) = 0.
1963    enddo
1964 !
1965    do k = kts,kmax
1966      do i = its,ite
1967        if(cnvflg(i).and.k.le.ktcon(i)) then
1968          aup = 1.
1969          if(k.le.kb(i)) aup = 0.
1970          adw = 1.
1971          if(k.gt.jmin(i)) adw = 0.
1972          dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1973          t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1974          q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1975          tem=1./rcs
1976          u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1977          v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1978          dp = 1000. * del(i,k)
1979          delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
1980          delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
1981          deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
1982          delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
1983          delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
1984        endif
1985      enddo
1986    enddo
1987 !
1988    do i = its,ite
1989      if(cnvflg(i)) then
1990      endif
1991    enddo
1992 !
1993    do k = kts,kmax
1994      do i = its,ite
1995        if (cnvflg(i) .and. k.le.ktcon(i)) then
1996          qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1997          qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
1998          qeso(i,k) = max(qeso(i,k), qmin_ )
1999        endif
2000      enddo
2001    enddo
2002 !
2003    do i = its,ite
2004      rntot(i) = 0.
2005      delqev(i) = 0.
2006      delq2(i) = 0.
2007      flg(i) = cnvflg(i)
2008    enddo
2009 !
2010 !  comptute rainfall
2011 !
2012    do k = kmax,kts,-1
2013      do i = its,ite
2014        if(cnvflg(i).and.k.lt.ktcon(i)) then
2015          aup = 1.
2016          if(k.le.kb(i)) aup = 0.
2017          adw = 1.
2018          if(k.ge.jmin(i)) adw = 0.
2019          rntot(i) = rntot(i)                                                   &
2020                + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
2021                * xmb(i) * .001 * dt2
2022        endif
2023      enddo
2024    enddo
2025 !
2026 !  conversion rainfall (m) and compute the evaporation of falling raindrops
2027 !
2028    do k = kmax,kts,-1
2029      do i = its,ite
2030        delq(i) = 0.0
2031        deltv(i) = 0.0
2032        qevap(i) = 0.0
2033        if(cnvflg(i).and.k.lt.ktcon(i)) then
2034          aup = 1.
2035          if(k.le.kb(i)) aup = 0.
2036          adw = 1.
2037          if(k.ge.jmin(i)) adw = 0.
2038          rain(i) = rain(i)                                                     &
2039                + (aup * pwo(i,k) + adw * edto(i) * pwdo(i,k))                  &
2040                * xmb(i) * .001 * dt2
2041        endif
2042        if(flg(i).and.k.lt.ktcon(i)) then
2043          evef = edt(i) * evfacts
2044          if(slimsk(i).eq.1.) evef = edt(i) * evfactl
2045          qcond(i) = evef * (q1(i,k) - qeso(i,k)) / (1. + el2orc *              &
2046                   qeso(i,k) / t1(i,k)**2)
2047          dp = 1000. * del(i,k)
2048          if(rain(i).gt.0..and.qcond(i).lt.0.) then
2049            qevap(i) = -qcond(i) * (1. - exp(-.32 * sqrt(dt2 * rain(i))))
2050            qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
2051            delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
2052            if (delq2(i).gt.rntot(i)) then
2053              qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
2054              flg(i) = .false.
2055            endif
2056          endif
2057          if(rain(i).gt.0..and.qevap(i).gt.0.) then
2058            q1(i,k) = q1(i,k) + qevap(i)
2059            t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
2060            rain(i) = rain(i) - .001 * qevap(i) * dp / g_
2061            delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
2062            deltv(i) =  - (hvap_/cp_)*qevap(i)/dt2
2063            delq(i) =  + qevap(i)/dt2
2064          endif
2065          dellaq(i,k) = dellaq(i,k) + delq(i)/xmb(i)
2066          delqbar(i)  = delqbar(i) + delq(i)*dp/g_
2067          deltbar(i)  = deltbar(i) + deltv(i)*dp/g_
2068        endif
2069      enddo
2070    enddo
2071 !
2072 !
2073 ! consider the negative rain in the event of rain evaporation and downdrafts
2074 !
2075    do i = its,ite
2076      if(cnvflg(i)) then
2077        if(rain(i).lt.0..and..not.flg(i)) rain(i) = 0.
2078        if(rain(i).le.0.) then
2079          rain(i) = 0.
2080        else
2081          ktop(i) = ktcon(i)
2082          kbot(i) = kbcon(i)
2083          kuo(i) = 1
2084        endif
2085      endif
2086    enddo
2087 !
2088    do k = kts,kmax
2089      do i = its,ite
2090        if(cnvflg(i).and.rain(i).le.0.) then
2091           t1(i,k) = to(i,k)
2092           q1(i,k) = qo(i,k)
2093           u1(i,k) = uo(i,k)
2094           v1(i,k) = vo(i,k)
2095        endif
2096      enddo
2097    enddo
2098 !
2099 !  detrainment of cloud water and ice
2100 !
2101    if (ncloud.gt.0) then
2102      do k = kts,kmax
2103        do i = its,ite
2104          if (cnvflg(i) .and. rain(i).gt.0.) then
2105            if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
2106              tem  = dellal(i,k) * xmb(i) * dt2
2107              tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
2108              if (ncloud.ge.2) then
2109                qi2(i,k) = qi2(i,k) + tem * tem1            ! ice
2110                qc2(i,k) = qc2(i,k) + tem *(1.0-tem1)       ! water
2111              else
2112                qc2(i,k) = qc2(i,k) + tem
2113              endif
2114            endif
2115          endif
2116        enddo
2117      enddo
2118    endif
2119 !
2120    end subroutine nsas2d
2121 !===============================================================================
2122       REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
2123 !-------------------------------------------------------------------------------
2124       IMPLICIT NONE
2125 !-------------------------------------------------------------------------------
2126       REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti,      &
2127            xai,xbi,ttp,tr
2128       INTEGER :: ice
2129 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2130       ttp=t0c+0.01
2131       dldt=cvap-cliq
2132       xa=-dldt/rv
2133       xb=xa+hvap/(rv*ttp)
2134       dldti=cvap-cice
2135       xai=-dldti/rv
2136       xbi=xai+hsub/(rv*ttp)
2137       tr=ttp/t
2138       if(t.lt.ttp.and.ice.eq.1) then
2139         fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2140       else
2141         fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2142       endif
2143 !
2144       if (t.lt.180.) then
2145         tr=ttp/180.
2146         if(t.lt.ttp.and.ice.eq.1) then
2147           fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2148         else
2149           fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2150         endif
2151       endif
2152 !
2153       if (t.ge.330.) then
2154         tr=ttp/330
2155         if(t.lt.ttp.and.ice.eq.1) then
2156           fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
2157         else
2158           fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
2159         endif
2160       endif
2161 !
2162 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2163       END FUNCTION fpvs
2164 !===============================================================================
2165    subroutine nsasinit(rthcuten,rqvcuten,rqccuten,rqicuten,                    &
2166                       rucuten,rvcuten,                                         &
2167                       restart,p_qc,p_qi,p_first_scalar,                        &
2168                       allowed_to_read,                                         &
2169                       ids, ide, jds, jde, kds, kde,                            &
2170                       ims, ime, jms, jme, kms, kme,                            &
2171                       its, ite, jts, jte, kts, kte                  )
2172 !-------------------------------------------------------------------------------
2173    implicit none
2174 !-------------------------------------------------------------------------------
2175    logical , intent(in)           ::  allowed_to_read,restart
2176    integer , intent(in)           ::  ids, ide, jds, jde, kds, kde,            &
2177                                       ims, ime, jms, jme, kms, kme,            &
2178                                       its, ite, jts, jte, kts, kte
2179    integer , intent(in)           ::  p_first_scalar, p_qi, p_qc
2180    real,     dimension( ims:ime , kms:kme , jms:jme ) , intent(out) ::         &
2181                                                               rthcuten,        &
2182                                                               rqvcuten,        &
2183                                                                rucuten,        &
2184                                                                rvcuten,        &
2185                                                               rqccuten,        &
2186                                                               rqicuten
2187    integer :: i, j, k, itf, jtf, ktf
2188    jtf=min0(jte,jde-1)
2189    ktf=min0(kte,kde-1)
2190    itf=min0(ite,ide-1)
2191    if(.not.restart)then
2192      do j=jts,jtf
2193      do k=kts,ktf
2194      do i=its,itf
2195        rthcuten(i,k,j)=0.
2196        rqvcuten(i,k,j)=0.
2197        rucuten(i,k,j)=0.
2198        rvcuten(i,k,j)=0.
2199      enddo
2200      enddo
2201      enddo
2202      if (p_qc .ge. p_first_scalar) then
2203         do j=jts,jtf
2204         do k=kts,ktf
2205         do i=its,itf
2206            rqccuten(i,k,j)=0.
2207         enddo
2208         enddo
2209         enddo
2210      endif
2211      if (p_qi .ge. p_first_scalar) then
2212         do j=jts,jtf
2213         do k=kts,ktf
2214         do i=its,itf
2215            rqicuten(i,k,j)=0.
2216         enddo
2217         enddo
2218         enddo
2219      endif
2220    endif
2221       end subroutine nsasinit
2222 !
2223 !==============================================================================
2224 ! NCEP SCV (Shallow Convection Scheme)
2225 !==============================================================================
2226    subroutine nscv2d(delt,del,prsl,prsi,prslk,zl,zi,                           &
2227                  ncloud,qc2,qi2,q1,t1,rain,kbot,ktop,                          &
2228                  kuo,                                                          &
2229                  slimsk,dot,u1,v1,                                             &
2230                  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2,                     &
2231                  cice,xls,psat,                                                &
2232                  hpbl,hfx,qfx,                                                 &
2233                  pgcon,                                                        &
2234                  ids,ide, jds,jde, kds,kde,                                    &
2235                  ims,ime, jms,jme, kms,kme,                                    &
2236                  its,ite, jts,jte, kts,kte)
2237 !
2238 !-------------------------------------------------------------------------------
2239 !
2240 ! subprogram:    nscv2d           computes shallow-convective heating and moisng
2241 !
2242 ! abstract: computes non-precipitating convective heating and moistening
2243 !   using a one cloud type arakawa-schubert convection scheme as in the ncep
2244 !   sas scheme. the scheme has been operational at ncep gfs model since july 2010
2245 !   the scheme includes updraft and downdraft effects, but the cloud depth is
2246 !   limited less than 150 hpa.
2247 !
2248 ! developed by jong-il han and hua-lu pan
2249 !   implemented into wrf by jiheon jang and songyou hong
2250 !   module with cpp-based options is available in grims
2251 !
2252 ! program history log:
2253 !   10-07-01 jong-il han  initial operational implementation at ncep gfs
2254 !   10-12-01 jihyeon jang implemented into wrf
2255 !
2256 ! subprograms called:
2257 !   fpvs     - function to compute saturation vapor pressure
2258 !
2259 ! references:
2260 !   han and pan (2010, wea. forecasting)
2261 !
2262 !-------------------------------------------------------------------------------
2263    implicit none
2264 !-------------------------------------------------------------------------------
2265 !  in/out variables
2266 !
2267    integer         ::  ids,ide, jds,jde, kds,kde,                              &
2268                        ims,ime, jms,jme, kms,kme,                              &
2269                        its,ite, jts,jte, kts,kte
2270    real            ::  cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
2271    real            ::  pi_,qmin_,t0c_
2272    real            ::  cice,xlv0,xls,psat
2273 !
2274    real            ::  delt
2275    real            ::  del(its:ite,kts:kte),                                   &
2276                        prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme),           &
2277                        prsi(its:ite,kts:kte+1),zl(its:ite,kts:kte)
2278    integer         ::  ncloud
2279    real            ::  slimsk(ims:ime)
2280    real            ::  dot(its:ite,kts:kte)
2281    real            ::  hpbl(ims:ime)
2282    real            ::  rcs
2283    real            ::  hfx(ims:ime),qfx(ims:ime)
2284 !
2285    real            ::  qi2(its:ite,kts:kte),qc2(its:ite,kts:kte)
2286    real            ::  q1(its:ite,kts:kte),                                    &
2287                        t1(its:ite,kts:kte),                                    &
2288                        u1(its:ite,kts:kte),                                    &
2289                        v1(its:ite,kts:kte)
2290    integer         ::  kuo(its:ite)
2291 !
2292    real            ::  rain(its:ite)
2293    integer         ::  kbot(its:ite),ktop(its:ite)
2294    real            ::  pgcon
2295 !
2296 !  local variables and arrays
2297 !
2298    integer         ::  i,j,indx, jmn, k, kk, km1
2299    integer         ::  kpbl(its:ite)
2300 !
2301    real            ::  dellat,                                                 &
2302                        desdt,   deta,    detad,   dg,                          &
2303                        dh,      dhh,     dlnsig,  dp,                          &
2304                        dq,      dqsdp,   dqsdt,   dt,                          &
2305                        dt2,     dtmax,   dtmin,                                &
2306                        dv1h,    dv2h,    dv3h,                                 &
2307                        dv1q,    dv2q,    dv3q,                                 &
2308                        dv1u,    dv2u,    dv3u,                                 &
2309                        dv1v,    dv2v,    dv3v,                                 &
2310                        dz,      dz1,     e1,      clam,                        &
2311                        aafac,                                                  &
2312                        es,      etah,                                          &
2313                        evef,    evfact,  evfactl,                              &
2314                        factor,  fjcap,                                         &
2315                        gamma,   pprime,  betaw,                                &
2316                        qlk,     qrch,    qs,                                   &
2317                        rfact,   shear,   tem1,                                 &
2318                        tem2,    val,     val1,                                 &
2319                        val2,    w1,      w1l,     w1s,                         &
2320                        w2,      w2l,     w2s,     w3,                          &
2321                        w3l,     w3s,     w4,      w4l,                         &
2322                        w4s,     tem,     ptem,    ptem1
2323 !
2324    integer         ::  kb(its:ite), kbcon(its:ite), kbcon1(its:ite),           &
2325                        ktcon(its:ite), ktcon1(its:ite),                        &
2326                        kbm(its:ite), kmax(its:ite)
2327 !
2328    real            ::  aa1(its:ite),                                           &
2329                        delhbar(its:ite), delq(its:ite),                        &
2330                        delq2(its:ite),   delqev(its:ite), rntot(its:ite),      &
2331                        delqbar(its:ite), deltbar(its:ite),                     &
2332                        deltv(its:ite),   edt(its:ite),                         &
2333                        wstar(its:ite),   sflx(its:ite),                        &
2334                        pdot(its:ite),    po(its:ite,kts:kte),                  &
2335                        qcond(its:ite),   qevap(its:ite),  hmax(its:ite),       &
2336                        vshear(its:ite),                                        &
2337                        xlamud(its:ite),  xmb(its:ite),    xmbmax(its:ite)
2338    real            ::  delubar(its:ite), delvbar(its:ite)
2339 !
2340    real            ::  cincr
2341 !
2342    real            ::  thx(its:ite, kts:kte)
2343    real            ::  rhox(its:ite)
2344    real            ::  tvcon
2345 !
2346    real            ::  p(its:ite,kts:kte),       to(its:ite,kts:kte),          &
2347                        qo(its:ite,kts:kte),      qeso(its:ite,kts:kte),        &
2348                        uo(its:ite,kts:kte),      vo(its:ite,kts:kte)
2349 !
2350 !  cloud water
2351 !
2352    real            ::  qlko_ktcon(its:ite),     dellal(its:ite,kts:kte),       &
2353                        dbyo(its:ite,kts:kte),                                  &
2354                        xlamue(its:ite,kts:kte),                                &
2355                        heo(its:ite,kts:kte),    heso(its:ite,kts:kte),         &
2356                        dellah(its:ite,kts:kte), dellaq(its:ite,kts:kte),       &
2357                        dellau(its:ite,kts:kte), dellav(its:ite,kts:kte),       &
2358                        ucko(its:ite,kts:kte),   vcko(its:ite,kts:kte),         &
2359                        hcko(its:ite,kts:kte),   qcko(its:ite,kts:kte),         &
2360                        eta(its:ite,kts:kte),    zi(its:ite,kts:kte+1),         &
2361                        pwo(its:ite,kts:kte)
2362 !
2363    logical         ::  totflg, cnvflg(its:ite), flg(its:ite)
2364 !
2365 !  physical parameters
2366 !
2367    real,parameter  ::  c0=.002,c1=5.e-4
2368    real,parameter  ::  cincrmax=180.,cincrmin=120.,dthk=25.
2369    real            ::  el2orc,fact1,fact2,eps
2370    real,parameter  ::  h1=0.33333333
2371    real,parameter  ::  tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)
2372 !
2373 !-------------------------------------------------------------------------------
2374 !
2375    pi_ = 3.14159
2376    qmin_ = 1.0e-30
2377    t0c_ = 273.15
2378    xlv0 = hvap_
2379       km1 = kte - 1
2380 !
2381 !  compute surface buoyancy flux
2382 !
2383       do k = kts,kte
2384         do i = its,ite
2385           thx(i,k) = t1(i,k)/prslk(i,k)
2386         enddo
2387       enddo
2388 !
2389       do i=its,ite
2390          tvcon = (1.+fv_*q1(i,1))
2391          rhox(i) = prsl(i,1)*1.e3/(rd_*t1(i,1)*tvcon)
2392       enddo
2393 !
2394       do i=its,ite
2395 !        sflx(i) = heat(i)+fv_*t1(i,1)*evap(i)
2396          sflx(i) = hfx(i)/rhox(i)/cp_ + qfx(i)/rhox(i)*fv_*thx(i,1)
2397       enddo
2398 !
2399 !  initialize arrays
2400 !
2401       do i=its,ite
2402         cnvflg(i) = .true.
2403         if(kuo(i).eq.1) cnvflg(i) = .false.
2404         if(sflx(i).le.0.) cnvflg(i) = .false.
2405         if(cnvflg(i)) then
2406           kbot(i)=kte+1
2407           ktop(i)=0
2408         endif
2409         rain(i)=0.
2410         kbcon(i)=kte
2411         ktcon(i)=1
2412         kb(i)=kte
2413         pdot(i) = 0.
2414         qlko_ktcon(i) = 0.
2415         edt(i)  = 0.
2416         aa1(i)  = 0.
2417         vshear(i) = 0.
2418       enddo
2419 !
2420       totflg = .true.
2421       do i=its,ite
2422         totflg = totflg .and. (.not. cnvflg(i))
2423       enddo
2424       if(totflg) return
2425 !
2426       dt2   =  delt
2427       val   =         1200.
2428       dtmin = max(dt2, val )
2429       val   =         3600.
2430       dtmax = max(dt2, val )
2431 !  model tunable parameters are all here
2432       clam    = .3
2433       aafac   = .1
2434       betaw   = .03
2435       evfact  = 0.3
2436       evfactl = 0.3
2437 ! namelist parameter...
2438 !      pgcon   = 0.55    ! Zhang & Wu (2003,JAS)
2439       val     =           1.
2440 !
2441 ! define miscellaneous values
2442 !
2443      el2orc = hvap_*hvap_/(rv_*cp_)
2444      eps    = rd_/rv_
2445      fact1  = (cvap_-cliq_)/rv_
2446      fact2  = hvap_/rv_-fact1*t0c_
2447 !
2448       w1l     = -8.e-3
2449       w2l     = -4.e-2
2450       w3l     = -5.e-3
2451       w4l     = -5.e-4
2452       w1s     = -2.e-4
2453       w2s     = -2.e-3
2454       w3s     = -1.e-3
2455       w4s     = -2.e-5
2456 !
2457 !  define top layer for search of the downdraft originating layer
2458 !  and the maximum thetae for updraft
2459 !
2460       do i=its,ite
2461         kbm(i)   = kte
2462         kmax(i)  = kte
2463       enddo
2464 !
2465       do k = kts, kte
2466         do i=its,ite
2467           if (prsl(i,k).gt.prsi(i,1)*0.70) kbm(i) = k + 1
2468           if (prsl(i,k).gt.prsi(i,1)*0.60) kmax(i) = k + 1
2469         enddo
2470       enddo
2471       do i=its,ite
2472         kbm(i)   = min(kbm(i),kmax(i))
2473       enddo
2474 !
2475 !  hydrostatic height assume zero terr and compute
2476 !  updraft entrainment rate as an inverse function of height
2477 !
2478       do k = kts, km1
2479         do i=its,ite
2480           xlamue(i,k) = clam / zi(i,k)
2481         enddo
2482       enddo
2483       do i=its,ite
2484         xlamue(i,kte) = xlamue(i,km1)
2485       enddo
2486 !
2487 !  pbl height
2488 !
2489       do i=its,ite
2490         flg(i) = cnvflg(i)
2491         kpbl(i)= 1
2492       enddo
2493 !
2494       do k = kts+1, km1
2495         do i=its,ite
2496           if (flg(i).and.zl(i,k).le.hpbl(i)) then
2497             kpbl(i) = k
2498           else
2499             flg(i) = .false.
2500           endif
2501         enddo
2502       enddo
2503 !
2504       do i=its,ite
2505         kpbl(i)= min(kpbl(i),kbm(i))
2506       enddo
2507 !
2508 !   convert surface pressure to mb from cb
2509 !
2510       rcs = 1.
2511       do k = kts, kte
2512         do i =its,ite
2513           if (cnvflg(i) .and. k .le. kmax(i)) then
2514             p(i,k) = prsl(i,k) * 10.0
2515             eta(i,k)  = 1.
2516             hcko(i,k) = 0.
2517             qcko(i,k) = 0.
2518             ucko(i,k) = 0.
2519             vcko(i,k) = 0.
2520             dbyo(i,k) = 0.
2521             pwo(i,k)  = 0.
2522             dellal(i,k) = 0.
2523             to(i,k)   = t1(i,k)
2524             qo(i,k)   = q1(i,k)
2525             uo(i,k)   = u1(i,k) * rcs
2526             vo(i,k)   = v1(i,k) * rcs
2527           endif
2528         enddo
2529       enddo
2530 !
2531 !
2532 !  column variables
2533 !  p is pressure of the layer (mb)
2534 !  t is temperature at t-dt (k)..tn
2535 !  q is mixing ratio at t-dt (kg/kg)..qn
2536 !  to is temperature at t+dt (k)... this is after advection and turbulan
2537 !  qo is mixing ratio at t+dt (kg/kg)..q1
2538 !
2539       do k = kts, kte
2540         do i=its,ite
2541           if (cnvflg(i) .and. k .le. kmax(i)) then
2542             qeso(i,k) = 0.01 * fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2543             qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
2544             val1      =             1.e-8
2545             qeso(i,k) = max(qeso(i,k), val1)
2546             val2      =           1.e-10
2547             qo(i,k)   = max(qo(i,k), val2 )
2548           endif
2549         enddo
2550       enddo
2551 !
2552 !  compute moist static energy
2553 !
2554       do k = kts,kte
2555         do i=its,ite
2556           if (cnvflg(i) .and. k .le. kmax(i)) then
2557             tem       = g_ * zl(i,k) + cp_ * to(i,k)
2558             heo(i,k)  = tem  + hvap_ * qo(i,k)
2559             heso(i,k) = tem  + hvap_ * qeso(i,k)
2560           endif
2561         enddo
2562       enddo
2563 !
2564 !  determine level with largest moist static energy within pbl
2565 !  this is the level where updraft starts
2566 !
2567       do i=its,ite
2568          if (cnvflg(i)) then
2569             hmax(i) = heo(i,1)
2570             kb(i) = 1
2571          endif
2572       enddo
2573 !
2574       do k = kts+1, kte
2575         do i=its,ite
2576           if (cnvflg(i).and.k.le.kpbl(i)) then
2577             if(heo(i,k).gt.hmax(i)) then
2578               kb(i)   = k
2579               hmax(i) = heo(i,k)
2580             endif
2581           endif
2582         enddo
2583       enddo
2584 !
2585       do k = kts, km1
2586         do i=its,ite
2587           if (cnvflg(i) .and. k .le. kmax(i)-1) then
2588             dz      = .5 * (zl(i,k+1) - zl(i,k))
2589             dp      = .5 * (p(i,k+1) - p(i,k))
2590             es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2591             pprime  = p(i,k+1) + (eps-1.) * es
2592             qs      = eps * es / pprime
2593             dqsdp   = - qs / pprime
2594             desdt   = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
2595             dqsdt   = qs * p(i,k+1) * desdt / (es * pprime)
2596             gamma   = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
2597             dt      = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
2598             dq      = dqsdt * dt + dqsdp * dp
2599             to(i,k) = to(i,k+1) + dt
2600             qo(i,k) = qo(i,k+1) + dq
2601             po(i,k) = .5 * (p(i,k) + p(i,k+1))
2602           endif
2603         enddo
2604       enddo
2605 !
2606       do k = kts, km1
2607         do i=its,ite
2608           if (cnvflg(i) .and. k .le. kmax(i)-1) then
2609             qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
2610             qeso(i,k) = eps * qeso(i,k) / (po(i,k) + (eps-1.) * qeso(i,k))
2611             val1      =             1.e-8
2612             qeso(i,k) = max(qeso(i,k), val1)
2613             val2      =           1.e-10
2614             qo(i,k)   = max(qo(i,k), val2 )
2615             heo(i,k)  = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                      &
2616                         cp_ * to(i,k) + hvap_ * qo(i,k)
2617             heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) +                      &
2618                         cp_ * to(i,k) + hvap_ * qeso(i,k)
2619             uo(i,k)   = .5 * (uo(i,k) + uo(i,k+1))
2620             vo(i,k)   = .5 * (vo(i,k) + vo(i,k+1))
2621           endif
2622         enddo
2623       enddo
2624 !
2625 !  look for the level of free convection as cloud base
2626 !
2627       do i=its,ite
2628         flg(i)   = cnvflg(i)
2629         if(flg(i)) kbcon(i) = kmax(i)
2630       enddo
2631 !
2632       do k = kts+1, km1
2633         do i=its,ite
2634           if (flg(i).and.k.lt.kbm(i)) then
2635             if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
2636               kbcon(i) = k
2637               flg(i)   = .false.
2638             endif
2639           endif
2640         enddo
2641       enddo
2642 !
2643       do i=its,ite
2644         if(cnvflg(i)) then
2645           if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
2646         endif
2647       enddo
2648 !
2649       totflg = .true.
2650       do i=its,ite
2651         totflg = totflg .and. (.not. cnvflg(i))
2652       enddo
2653       if(totflg) return
2654 !
2655 !  determine critical convective inhibition
2656 !  as a function of vertical velocity at cloud base.
2657 !
2658       do i=its,ite
2659         if(cnvflg(i)) then
2660           pdot(i)  = 10.* dot(i,kbcon(i))
2661         endif
2662       enddo
2663 !
2664       do i=its,ite
2665         if(cnvflg(i)) then
2666           if(slimsk(i).eq.1.) then
2667             w1 = w1l
2668             w2 = w2l
2669             w3 = w3l
2670             w4 = w4l
2671           else
2672             w1 = w1s
2673             w2 = w2s
2674             w3 = w3s
2675             w4 = w4s
2676           endif
2677           if(pdot(i).le.w4) then
2678             ptem = (pdot(i) - w4) / (w3 - w4)
2679           elseif(pdot(i).ge.-w4) then
2680             ptem = - (pdot(i) + w4) / (w4 - w3)
2681           else
2682             ptem = 0.
2683           endif
2684           val1    =             -1.
2685           ptem = max(ptem,val1)
2686           val2    =             1.
2687           ptem = min(ptem,val2)
2688           ptem = 1. - ptem
2689           ptem1= .5*(cincrmax-cincrmin)
2690           cincr = cincrmax - ptem * ptem1
2691           tem1 = p(i,kb(i)) - p(i,kbcon(i))
2692           if(tem1.gt.cincr) then
2693              cnvflg(i) = .false.
2694           endif
2695         endif
2696       enddo
2697 !
2698       totflg = .true.
2699       do i=its,ite
2700         totflg = totflg .and. (.not. cnvflg(i))
2701       enddo
2702       if(totflg) return
2703 !
2704 !  assume the detrainment rate for the updrafts to be same as
2705 !  the entrainment rate at cloud base
2706 !
2707       do i = its,ite
2708         if(cnvflg(i)) then
2709           xlamud(i) = xlamue(i,kbcon(i))
2710         endif
2711       enddo
2712 !
2713 !  determine updraft mass flux for the subcloud layers
2714 !
2715        do k = km1, kts, -1
2716         do i = its,ite
2717           if (cnvflg(i)) then
2718             if(k.lt.kbcon(i).and.k.ge.kb(i)) then
2719               dz       = zi(i,k+1) - zi(i,k)
2720               ptem     = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
2721               eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
2722             endif
2723           endif
2724         enddo
2725       enddo
2726 !
2727 !  compute mass flux above cloud base
2728 !
2729       do k = kts+1, km1
2730         do i = its,ite
2731          if(cnvflg(i))then
2732            if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
2733               dz       = zi(i,k) - zi(i,k-1)
2734               ptem     = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
2735               eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
2736            endif
2737          endif
2738         enddo
2739       enddo
2740 !
2741 !  compute updraft cloud property
2742 !
2743       do i = its,ite
2744         if(cnvflg(i)) then
2745           indx         = kb(i)
2746           hcko(i,indx) = heo(i,indx)
2747           ucko(i,indx) = uo(i,indx)
2748           vcko(i,indx) = vo(i,indx)
2749         endif
2750       enddo
2751 !
2752       do k = kts+1, km1
2753         do i = its,ite
2754           if (cnvflg(i)) then
2755             if(k.gt.kb(i).and.k.lt.kmax(i)) then
2756               dz   = zi(i,k) - zi(i,k-1)
2757               tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2758               tem1 = 0.5 * xlamud(i) * dz
2759               factor = 1. + tem - tem1
2760               ptem = 0.5 * tem + pgcon
2761               ptem1= 0.5 * tem - pgcon
2762               hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5*                      &
2763                           (heo(i,k)+heo(i,k-1)))/factor
2764               ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k)                  &
2765                           +ptem1*uo(i,k-1))/factor
2766               vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k)                  &
2767                           +ptem1*vo(i,k-1))/factor
2768               dbyo(i,k) = hcko(i,k) - heso(i,k)
2769             endif
2770           endif
2771         enddo
2772       enddo
2773 !
2774 !   taking account into convection inhibition due to existence of
2775 !    dry layers below cloud base
2776 !
2777       do i=its,ite
2778         flg(i) = cnvflg(i)
2779         kbcon1(i) = kmax(i)
2780       enddo
2781 !
2782       do k = kts+1, km1
2783         do i=its,ite
2784           if (flg(i).and.k.lt.kbm(i)) then
2785             if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
2786               kbcon1(i) = k
2787               flg(i)    = .false.
2788             endif
2789           endif
2790         enddo
2791       enddo
2792 !
2793       do i=its,ite
2794         if(cnvflg(i)) then
2795           if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
2796         endif
2797       enddo
2798 !
2799       do i=its,ite
2800         if(cnvflg(i)) then
2801           tem = p(i,kbcon(i)) - p(i,kbcon1(i))
2802           if(tem.gt.dthk) then
2803              cnvflg(i) = .false.
2804           endif
2805         endif
2806       enddo
2807 !
2808       totflg = .true.
2809       do i = its,ite
2810         totflg = totflg .and. (.not. cnvflg(i))
2811       enddo
2812       if(totflg) return
2813 !
2814 !  determine first guess cloud top as the level of zero buoyancy
2815 !    limited to the level of sigma=0.7
2816 !
2817       do i = its,ite
2818         flg(i) = cnvflg(i)
2819         if(flg(i)) ktcon(i) = kbm(i)
2820       enddo
2821 !
2822       do k = kts+1, km1
2823         do i=its,ite
2824           if (flg(i).and.k .lt. kbm(i)) then
2825             if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
2826                ktcon(i) = k
2827                flg(i)   = .false.
2828             endif
2829           endif
2830         enddo
2831       enddo
2832 !
2833 !  specify upper limit of mass flux at cloud base
2834 !
2835       do i = its,ite
2836         if(cnvflg(i)) then
2837           k = kbcon(i)
2838           dp = 1000. * del(i,k)
2839           xmbmax(i) = dp / (g_ * dt2)
2840         endif
2841       enddo
2842 !
2843 !  compute cloud moisture property and precipitation
2844 !
2845       do i = its,ite
2846         if (cnvflg(i)) then
2847           aa1(i) = 0.
2848           qcko(i,kb(i)) = qo(i,kb(i))
2849         endif
2850       enddo
2851 !
2852       do k = kts+1, km1
2853         do i = its,ite
2854           if (cnvflg(i)) then
2855             if(k.gt.kb(i).and.k.lt.ktcon(i)) then
2856               dz    = zi(i,k) - zi(i,k-1)
2857               gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2858               qrch = qeso(i,k)                                                 &
2859                    + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2860               tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2861               tem1 = 0.5 * xlamud(i) * dz
2862               factor = 1. + tem - tem1
2863               qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                      &
2864                           (qo(i,k)+qo(i,k-1)))/factor
2865               dq = eta(i,k) * (qcko(i,k) - qrch)
2866 !
2867 !             rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
2868 !
2869 !  below lfc check if there is excess moisture to release latent heat
2870 !
2871               if(k.ge.kbcon(i).and.dq.gt.0.) then
2872                 etah = .5 * (eta(i,k) + eta(i,k-1))
2873                 if(ncloud.gt.0) then
2874                   dp = 1000. * del(i,k)
2875                   qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2876                   dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2877                 else
2878                   qlk = dq / (eta(i,k) + etah * c0 * dz)
2879                 endif
2880                 aa1(i) = aa1(i) - dz * g_ * qlk
2881                 qcko(i,k)= qlk + qrch
2882                 pwo(i,k) = etah * c0 * dz * qlk
2883               endif
2884             endif
2885           endif
2886         enddo
2887       enddo
2888 !
2889 !  calculate cloud work function
2890 !
2891       do k = kts+1, km1
2892         do i = its,ite
2893           if (cnvflg(i)) then
2894             if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
2895               dz1 = zl(i,k+1) - zl(i,k)
2896               gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2897               rfact =  1. + fv_ * cp_ * gamma                                  &
2898                       * to(i,k) / hvap_
2899               aa1(i) = aa1(i) + dz1 * (g_ / (cp_ * to(i,k)))                   &
2900                       * dbyo(i,k) / (1. + gamma) * rfact
2901               val = 0.
2902               aa1(i)=aa1(i)+ dz1 * g_ * fv_ *                                  &
2903                       max(val,(qeso(i,k) - qo(i,k)))
2904             endif
2905           endif
2906         enddo
2907       enddo
2908 !
2909       do i = its,ite
2910         if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
2911       enddo
2912 !
2913       totflg = .true.
2914       do i=its,ite
2915         totflg = totflg .and. (.not. cnvflg(i))
2916       enddo
2917       if(totflg) return
2918 !
2919 !  estimate the convective overshooting as the level
2920 !    where the [aafac * cloud work function] becomes zero,
2921 !    which is the final cloud top limited to the level of sigma=0.7
2922 !
2923       do i = its,ite
2924         if (cnvflg(i)) then
2925           aa1(i) = aafac * aa1(i)
2926         endif
2927       enddo
2928 !
2929       do i = its,ite
2930         flg(i) = cnvflg(i)
2931         ktcon1(i) = kbm(i)
2932       enddo
2933 !
2934       do k = kts+1, km1
2935         do i = its,ite
2936           if (flg(i)) then
2937             if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
2938               dz1 = zl(i,k+1) - zl(i,k)
2939               gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2940               rfact =  1. + fv_ * cp_ * gamma                                  &
2941                       * to(i,k) / hvap_
2942               aa1(i) = aa1(i) +                                                &
2943                       dz1 * (g_ / (cp_ * to(i,k)))                             &
2944                       * dbyo(i,k) / (1. + gamma) * rfact
2945               if(aa1(i).lt.0.) then
2946                 ktcon1(i) = k
2947                 flg(i) = .false.
2948               endif
2949             endif
2950           endif
2951         enddo
2952       enddo
2953 !
2954 !  compute cloud moisture property, detraining cloud water
2955 !    and precipitation in overshooting layers
2956 !
2957       do k = kts+1, km1
2958         do i = its,ite
2959           if (cnvflg(i)) then
2960             if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
2961               dz    = zi(i,k) - zi(i,k-1)
2962               gamma = el2orc * qeso(i,k) / (to(i,k)**2)
2963               qrch = qeso(i,k)                                                 &
2964                    + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
2965               tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
2966               tem1 = 0.5 * xlamud(i) * dz
2967               factor = 1. + tem - tem1
2968               qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5*                      &
2969                           (qo(i,k)+qo(i,k-1)))/factor
2970               dq = eta(i,k) * (qcko(i,k) - qrch)
2971 !
2972 !  check if there is excess moisture to release latent heat
2973 !
2974               if(dq.gt.0.) then
2975                 etah = .5 * (eta(i,k) + eta(i,k-1))
2976                 if(ncloud.gt.0) then
2977                   dp = 1000. * del(i,k)
2978                   qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
2979                   dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
2980                 else
2981                   qlk = dq / (eta(i,k) + etah * c0 * dz)
2982                 endif
2983                 qcko(i,k) = qlk + qrch
2984                 pwo(i,k) = etah * c0 * dz * qlk
2985               endif
2986             endif
2987           endif
2988         enddo
2989       enddo
2990 !
2991 ! exchange ktcon with ktcon1
2992 !
2993       do i = its,ite
2994         if(cnvflg(i)) then
2995           kk = ktcon(i)
2996           ktcon(i) = ktcon1(i)
2997           ktcon1(i) = kk
2998         endif
2999       enddo
3000 !
3001 !  this section is ready for cloud water
3002 !
3003       if(ncloud.gt.0) then
3004 !
3005 !  compute liquid and vapor separation at cloud top
3006 !
3007       do i = its,ite
3008         if(cnvflg(i)) then
3009           k = ktcon(i) - 1
3010           gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3011           qrch = qeso(i,k)                                                     &
3012                + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
3013           dq = qcko(i,k) - qrch
3014 !
3015 !  check if there is excess moisture to release latent heat
3016 !
3017           if(dq.gt.0.) then
3018             qlko_ktcon(i) = dq
3019             qcko(i,k) = qrch
3020           endif
3021         endif
3022       enddo
3023 !
3024       endif
3025 !
3026 !--- compute precipitation efficiency in terms of windshear
3027 !
3028       do i = its,ite
3029         if(cnvflg(i)) then
3030           vshear(i) = 0.
3031         endif
3032       enddo
3033 !
3034       do k = kts+1,kte
3035         do i = its,ite
3036           if (cnvflg(i)) then
3037             if(k.gt.kb(i).and.k.le.ktcon(i)) then
3038               shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2                             &
3039                         + (vo(i,k)-vo(i,k-1)) ** 2)
3040               vshear(i) = vshear(i) + shear
3041             endif
3042           endif
3043         enddo
3044       enddo
3045 !
3046       do i = its,ite
3047         if(cnvflg(i)) then
3048           vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
3049           e1=1.591-.639*vshear(i)                                              &
3050              +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
3051           edt(i)=1.-e1
3052           val =         .9
3053           edt(i) = min(edt(i),val)
3054           val =         .0
3055           edt(i) = max(edt(i),val)
3056         endif
3057       enddo
3058 !
3059 !--- what would the change be, that a cloud with unit mass
3060 !--- will do to the environment?
3061 !
3062       do k = kts,kte
3063         do i = its,ite
3064           if(cnvflg(i) .and. k .le. kmax(i)) then
3065             dellah(i,k) = 0.
3066             dellaq(i,k) = 0.
3067             dellau(i,k) = 0.
3068             dellav(i,k) = 0.
3069           endif
3070         enddo
3071       enddo
3072 !
3073 !--- changed due to subsidence and entrainment
3074 !
3075       do k = kts+1, km1
3076         do i = its,ite
3077           if (cnvflg(i)) then
3078             if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3079               dp = 1000. * del(i,k)
3080               dz = zi(i,k) - zi(i,k-1)
3081 !
3082               dv1h = heo(i,k)
3083               dv2h = .5 * (heo(i,k) + heo(i,k-1))
3084               dv3h = heo(i,k-1)
3085               dv1q = qo(i,k)
3086               dv2q = .5 * (qo(i,k) + qo(i,k-1))
3087               dv3q = qo(i,k-1)
3088               dv1u = uo(i,k)
3089               dv2u = .5 * (uo(i,k) + uo(i,k-1))
3090               dv3u = uo(i,k-1)
3091               dv1v = vo(i,k)
3092               dv2v = .5 * (vo(i,k) + vo(i,k-1))
3093               dv3v = vo(i,k-1)
3094 !
3095               tem  = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
3096               tem1 = xlamud(i)
3097 !
3098               dellah(i,k) = dellah(i,k) +                                      &
3099           ( eta(i,k)*dv1h - eta(i,k-1)*dv3h                                    &
3100          -  tem*eta(i,k-1)*dv2h*dz                                             &
3101          +  tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz                      &
3102               ) *g_/dp
3103 !
3104               dellaq(i,k) = dellaq(i,k) +                                      &
3105           ( eta(i,k)*dv1q - eta(i,k-1)*dv3q                                    &
3106          -  tem*eta(i,k-1)*dv2q*dz                                             &
3107          +  tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz                      &
3108               ) *g_/dp
3109 !
3110               dellau(i,k) = dellau(i,k) +                                      &
3111           ( eta(i,k)*dv1u - eta(i,k-1)*dv3u                                    &
3112          -  tem*eta(i,k-1)*dv2u*dz                                             &
3113          +  tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz                      &
3114          -  pgcon*eta(i,k-1)*(dv1u-dv3u)                                       &
3115               ) *g_/dp
3116 !
3117               dellav(i,k) = dellav(i,k) +                                      &
3118           ( eta(i,k)*dv1v - eta(i,k-1)*dv3v                                    &
3119          -  tem*eta(i,k-1)*dv2v*dz                                             &
3120          +  tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz                      &
3121          -  pgcon*eta(i,k-1)*(dv1v-dv3v)                                       &
3122               ) *g_/dp
3123 !
3124             endif
3125           endif
3126         enddo
3127       enddo
3128 !
3129 !------- cloud top
3130 !
3131       do i = its,ite
3132         if(cnvflg(i)) then
3133           indx = ktcon(i)
3134           dp = 1000. * del(i,indx)
3135           dv1h = heo(i,indx-1)
3136           dellah(i,indx) = eta(i,indx-1) *                                     &
3137                           (hcko(i,indx-1) - dv1h) * g_ / dp
3138           dv1q = qo(i,indx-1)
3139           dellaq(i,indx) = eta(i,indx-1) *                                     &
3140                           (qcko(i,indx-1) - dv1q) * g_ / dp
3141           dv1u = uo(i,indx-1)
3142           dellau(i,indx) = eta(i,indx-1) *                                     &
3143                           (ucko(i,indx-1) - dv1u) * g_ / dp
3144           dv1v = vo(i,indx-1)
3145           dellav(i,indx) = eta(i,indx-1) *                                     &
3146                           (vcko(i,indx-1) - dv1v) * g_ / dp
3147 !
3148 !  cloud water
3149 !
3150           dellal(i,indx) = eta(i,indx-1) *                                     &
3151                           qlko_ktcon(i) * g_ / dp
3152         endif
3153       enddo
3154 !
3155 !  mass flux at cloud base for shallow convection
3156 !  (Grant, 2001)
3157 !
3158       do i= its,ite
3159         if(cnvflg(i)) then
3160           k = kbcon(i)
3161           ptem = g_*sflx(i)*hpbl(i)/t1(i,1)
3162           wstar(i) = ptem**h1
3163           tem = po(i,k)*100. / (rd_*t1(i,k))
3164           xmb(i) = betaw*tem*wstar(i)
3165           xmb(i) = min(xmb(i),xmbmax(i))
3166         endif
3167       enddo
3168 !
3169       do k = kts,kte
3170         do i = its,ite
3171           if (cnvflg(i) .and. k .le. kmax(i)) then
3172             qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
3173             qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
3174             val     =             1.e-8
3175             qeso(i,k) = max(qeso(i,k), val )
3176           endif
3177         enddo
3178       enddo
3179 !
3180       do i = its,ite
3181         delhbar(i) = 0.
3182         delqbar(i) = 0.
3183         deltbar(i) = 0.
3184         delubar(i) = 0.
3185         delvbar(i) = 0.
3186         qcond(i) = 0.
3187       enddo
3188 !
3189       do k = kts,kte
3190         do i = its,ite
3191           if (cnvflg(i)) then
3192             if(k.gt.kb(i).and.k.le.ktcon(i)) then
3193               dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
3194               t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
3195               q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
3196               tem = 1./rcs
3197               u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
3198               v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
3199               dp = 1000. * del(i,k)
3200               delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
3201               delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
3202               deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
3203               delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
3204               delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
3205             endif
3206           endif
3207         enddo
3208       enddo
3209 !
3210       do k = kts,kte
3211         do i = its,ite
3212           if (cnvflg(i)) then
3213             if(k.gt.kb(i).and.k.le.ktcon(i)) then
3214               qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls &
3215                         ,psat,t0c_)
3216               qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
3217               val     =             1.e-8
3218               qeso(i,k) = max(qeso(i,k), val )
3219             endif
3220           endif
3221         enddo
3222       enddo
3223 !
3224       do i = its,ite
3225         rntot(i) = 0.
3226         delqev(i) = 0.
3227         delq2(i) = 0.
3228         flg(i) = cnvflg(i)
3229       enddo
3230 !
3231       do k = kte, kts, -1
3232         do i = its,ite
3233           if (cnvflg(i)) then
3234             if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3235               rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
3236             endif
3237           endif
3238         enddo
3239       enddo
3240 !
3241 ! evaporating rain
3242 !
3243       do k = kte, kts, -1
3244         do i = its,ite
3245           if (k .le. kmax(i)) then
3246             deltv(i) = 0.
3247             delq(i) = 0.
3248             qevap(i) = 0.
3249             if(cnvflg(i)) then
3250               if(k.lt.ktcon(i).and.k.gt.kb(i)) then
3251                 rain(i) = rain(i) + pwo(i,k) * xmb(i) * .001 * dt2
3252               endif
3253             endif
3254             if(flg(i).and.k.lt.ktcon(i)) then
3255               evef = edt(i) * evfact
3256               if(slimsk(i).eq.1.) evef=edt(i) * evfactl
3257               qcond(i) = evef * (q1(i,k) - qeso(i,k))                          &
3258                        / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
3259               dp = 1000. * del(i,k)
3260               if(rain(i).gt.0..and.qcond(i).lt.0.) then
3261                 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rain(i))))
3262                 qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
3263                 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
3264               endif
3265               if(rain(i).gt.0..and.qcond(i).lt.0..and.                         &
3266                  delq2(i).gt.rntot(i)) then
3267                 qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
3268                 flg(i) = .false.
3269               endif
3270               if(rain(i).gt.0..and.qevap(i).gt.0.) then
3271                 tem  = .001 * dp / g_
3272                 tem1 = qevap(i) * tem
3273                 if(tem1.gt.rain(i)) then
3274                   qevap(i) = rain(i) / tem
3275                   rain(i) = 0.
3276                 else
3277                   rain(i) = rain(i) - tem1
3278                 endif
3279                 q1(i,k) = q1(i,k) + qevap(i)
3280                 t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
3281                 deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
3282                 delq(i) =  + qevap(i)/dt2
3283                 delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
3284               endif
3285               dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
3286               delqbar(i) = delqbar(i) + delq(i)*dp/g_
3287               deltbar(i) = deltbar(i) + deltv(i)*dp/g_
3288             endif
3289           endif
3290         enddo
3291       enddo
3292 !
3293       do i = its,ite
3294         if(cnvflg(i)) then
3295           if(rain(i).lt.0..or..not.flg(i)) rain(i) = 0.
3296           ktop(i) = ktcon(i)
3297           kbot(i) = kbcon(i)
3298           kuo(i) = 0
3299         endif
3300       enddo
3301 !
3302 ! cloud water
3303 !
3304       if (ncloud.gt.0) then
3305 !
3306       do k = kts, km1
3307         do i = its,ite
3308           if (cnvflg(i)) then
3309             if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
3310               tem  = dellal(i,k) * xmb(i) * dt2
3311               tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
3312               if (ncloud.ge.4) then
3313                 qi2(i,k) = qi2(i,k) + tem * tem1            ! ice
3314                 qc2(i,k) = qc2(i,k) + tem *(1.0-tem1)       ! water
3315               else
3316                 qc2(i,k) = qc2(i,k) + tem
3317               endif
3318             endif
3319           endif
3320         enddo
3321       enddo
3322 !
3323       endif
3324 !
3325       end subroutine nscv2d
3326 !-------------------------------------------------------------------------------
3327 !
3328 END MODULE module_cu_nsas