wrfv2_fire/phys/module_cu_nsas.F

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