phys/module_cu_gf_sh.F

   1 ! module cup_gf_sh will call shallow convection as described in Grell and
   2 ! Freitas (2016). Input variables are:
   3 !    zo               Height at model levels
   4 !    t,tn             Temperature without and with forcing at model levels
   5 !    q,qo             mixing ratio without and with forcing at model levels
   6 !    po               pressure at model levels (mb)
   7 !    psur             surface pressure (mb)
   8 !    z1               surface height
   9 !    dhdt             forcing for boundary layer equilibrium
  10 !    hfx,qfx          in w/m2 (positive, if upward from sfc)
  11 !    kpbl             level of boundaty layer height
  12 !    xland            land mask (1. for land)
  13 !    ichoice          which closure to choose
  14 !                     1: old g
  15 !                     2: zws
  16 !                     3: dhdt
  17 !                     0: average
  18 !    tcrit            parameter for water/ice conversion (258)
  19 !
  20 !!!!!!!!!!!! Variables that are diagnostic
  21 !
  22 !    zuo               normalized mass flux profile
  23 !    xmb_out           base mass flux
  24 !    kbcon             convective cloud base
  25 !    ktop              cloud top
  26 !    k22               level of updraft originating air
  27 !    ierr              error flag
  28 !    ierrc             error description
  29 !
  30 !!!!!!!!!!!! Variables that are on output
  31 !    outt               temperature tendency (K/s)
  32 !    outq               mixing ratio tendency (kg/kg/s)
  33 !    outqc              cloud water/ice tendency (kg/kg/s)
  34 !    pre                precip rate (mm/s)
  35 !    cupclw             incloud mixing ratio of cloudwater/ice (for radiation)
  36 !                       this needs heavy tuning factors, since cloud fraction is
  37 !                       not included (kg/kg)
  38 !    cnvwt              required for GFS physics
  39 !
  40 !    itf,ktf,its,ite, kts,kte are dimensions
  41 !    ztexec,zqexec    excess temperature and moisture for updraft
  42 MODULE module_cu_gf_sh
  43
  44 #if ( WRF_CHEM == 1 )
  45      USE module_cu_gf_ctrans,only: ctrans_gf
  46 #endif
  47
  48     real, parameter:: c1_shal=0.! .0005
  49     real, parameter:: g  =9.81
  50     real, parameter:: cp =1004.
  51     real, parameter:: xlv=2.5e6
  52     real, parameter:: r_v=461.
  53     real, parameter:: c0_shal=.001
  54     real, parameter:: fluxtune=1.5
  55
  56
  57 contains
  58   SUBROUTINE CUP_gf_sh (                                              &
  59 ! input variables, must be supplied
  60                          zo,T,Q,Z1,TN,QO,PO,PSUR,dhdt,kpbl,rho,     &
  61                          hfx,qfx,xland,ichoice,tcrit,dtime, &
  62 ! input variables. Ierr should be initialized to zero or larger than zero for
  63 ! turning off shallow convection for grid points
  64                          zuo,xmb_out,kbcon,ktop,k22,ierr,ierrc,    &
  65 ! output tendencies
  66                          OUTT,OUTQ,OUTQC,cnvwt,pre,cupclw,             &
  67 #if ( WRF_CHEM == 1 )
  68                          num_chem,chem2d,outchemt,             &
  69                          num_tracer,tracer2d,outtracert,       &
  70                          numgas,chemopt,traceropt,             &
  71                          conv_tr_wetscav,conv_tr_aqchem,       &
  72                          chem_conv_tr,                         &
  73 #endif
  74 ! dimesnional variables
  75                          itf,ktf,its,ite, kts,kte,ipr)
  76 !
  77 ! this module needs some subroutines from gf_deep
  78 !
  79   use module_cu_gf_deep,only:cup_env,cup_env_clev,get_cloud_bc,cup_minimi,  &
  80                       get_inversion_layers,rates_up_pdf,get_cloud_bc,     &
  81                       cup_up_aa0,cup_kbcon,get_lateral_massflux
  82      implicit none
  83      integer                                                           &
  84         ,intent (in   )                   ::                           &
  85         itf,ktf,        &
  86         its,ite, kts,kte,ipr
  87      logical :: MAKE_CALC_FOR_XK = .true.
  88      integer, intent (in   )              ::                           &
  89         ichoice
  90   !
  91   !
  92   !
  93   ! outtem = output temp tendency (per s)
  94   ! outq   = output q tendency (per s)
  95   ! outqc  = output qc tendency (per s)
  96   ! pre    = output precip
  97      real,    dimension (its:ite,kts:kte)                              &
  98         ,intent (inout  )                   ::                           &
  99         cnvwt,OUTT,OUTQ,OUTQC,cupclw,zuo
 100      real,    dimension (its:ite)                                      &
 101         ,intent (out  )                   ::                           &
 102         xmb_out
 103      integer,    dimension (its:ite)                                   &
 104         ,intent (inout  )                   ::                           &
 105         ierr
 106      integer,    dimension (its:ite)                                   &
 107         ,intent (out  )                   ::                           &
 108         kbcon,ktop,k22
 109      integer,    dimension (its:ite)                                   &
 110         ,intent (in  )                   ::                           &
 111         kpbl
 112   !
 113   ! basic environmental input includes a flag (ierr) to turn off
 114   ! convection for this call only and at that particular gridpoint
 115   !
 116      real,    dimension (its:ite,kts:kte)                              &
 117         ,intent (in   )                   ::                           &
 118         T,PO,tn,dhdt,rho
 119      real,    dimension (its:ite,kts:kte)                              &
 120         ,intent (inout)                   ::                           &
 121          Q,QO
 122      real, dimension (its:ite)                                         &
 123         ,intent (in   )                   ::                           &
 124         xland,Z1,PSUR,hfx,qfx
 125
 126        real                                                            &
 127         ,intent (in   )                   ::                           &
 128         dtime,tcrit
 129
 130 #if ( WRF_CHEM == 1 )
 131        INTEGER,INTENT(IN   ) ::                                        &
 132               num_chem,num_tracer,numgas,chemopt,traceropt,            &
 133               conv_tr_wetscav,conv_tr_aqchem,chem_conv_tr
 134        REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN)::     &
 135                         tracer2d
 136        REAL,DIMENSION(its:ite , kts:kte , num_chem),INTENT(IN)::     &
 137                         chem2d
 138        REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT)::  &
 139                         outtracert
 140        REAL,DIMENSION(its:ite , kts:kte , num_tracer),INTENT(INOUT)::  &
 141                         outchemt
 142        INTEGER :: nv
 143        real,dimension(its:ite,kts:kte) ::  tempco
 144        real,dimension(its:ite,kts:kte) ::                           &
 145                      zdo,clw_all,pwdo,dd_massentro,dd_massdetro
 146        real,dimension(its:ite)::                               &
 147                      pwevo,pwavo,edto
 148        integer,dimension(its:ite)::                            &
 149                      jmin
 150 #endif
 151
 152   !
 153   !***************** the following are your basic environmental
 154   !                  variables. They carry a "_cup" if they are
 155   !                  on model cloud levels (staggered). They carry
 156   !                  an "o"-ending (z becomes zo), if they are the forced
 157   !                  variables.
 158   !
 159   ! z           = heights of model levels
 160   ! q           = environmental mixing ratio
 161   ! qes         = environmental saturation mixing ratio
 162   ! t           = environmental temp
 163   ! p           = environmental pressure
 164   ! he          = environmental moist static energy
 165   ! hes         = environmental saturation moist static energy
 166   ! z_cup       = heights of model cloud levels
 167   ! q_cup       = environmental q on model cloud levels
 168   ! qes_cup     = saturation q on model cloud levels
 169   ! t_cup       = temperature (Kelvin) on model cloud levels
 170   ! p_cup       = environmental pressure
 171   ! he_cup = moist static energy on model cloud levels
 172   ! hes_cup = saturation moist static energy on model cloud levels
 173   ! gamma_cup = gamma on model cloud levels
 174   ! dby = buoancy term
 175   ! entr = entrainment rate
 176   ! bu = buoancy term
 177   ! gamma_cup = gamma on model cloud levels
 178   ! qrch = saturation q in cloud
 179   ! pwev = total normalized integrated evaoprate (I2)
 180   ! z1 = terrain elevation
 181   ! psur        = surface pressure
 182   ! zu      = updraft normalized mass flux
 183   ! kbcon       = LFC of parcel from k22
 184   ! k22         = updraft originating level
 185   ! ichoice       = flag if only want one closure (usually set to zero!)
 186   ! dby = buoancy term
 187   ! ktop = cloud top (output)
 188   ! xmb    = total base mass flux
 189   ! hc = cloud moist static energy
 190   ! hkb = moist static energy at originating level
 191
 192      real,    dimension (its:ite,kts:kte) ::                           &
 193         entr_rate_2d,he,hes,qes,z,                      &
 194         heo,heso,qeso,zo,                                              &
 195         xhe,xhes,xqes,xz,xt,xq,                                        &
 196         qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup,      &
 197         qeso_cup,qo_cup,heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,     &
 198         tn_cup,                                                        &
 199         xqes_cup,xq_cup,xhe_cup,xhes_cup,xz_cup,     &
 200         xt_cup,dby,hc,zu,   &
 201         dbyo,qco,pwo,hco,qrco,     &
 202         dbyt,xdby,xhc,xzu,            &
 203
 204   ! cd  = detrainment function for updraft
 205   ! dellat = change of temperature per unit mass flux of cloud ensemble
 206   ! dellaq = change of q per unit mass flux of cloud ensemble
 207   ! dellaqc = change of qc per unit mass flux of cloud ensemble
 208
 209         cd,DELLAH,DELLAQ,DELLAT,DELLAQC
 210
 211   ! aa0 cloud work function for downdraft
 212   ! aa0     = cloud work function without forcing effects
 213   ! aa1     = cloud work function with forcing effects
 214   ! xaa0    = cloud work function with cloud effects (ensemble dependent)
 215
 216      real,    dimension (its:ite) ::                                   &
 217        zws,ztexec,zqexec,pre,AA1,AA0,XAA0,HKB,                          &
 218        flux_tun,HKBO,XHKB,                                    &
 219        rand_vmas,xmbmax,XMB,                         &
 220        cap_max,entr_rate,                                    &
 221        cap_max_increment
 222      integer,    dimension (its:ite) ::                                &
 223        kstabi,xland1,KBMAX,ktopx
 224
 225      integer                              ::                           &
 226        I,K,ki
 227      real                                 ::                           &
 228       dz,mbdt,zkbmax,      &
 229       cap_maxs,trash,trash2,frh
 230
 231       real buo_flux,pgeoh,dp,entup,detup,totmas
 232
 233      real xff_shal(3),blqe,xkshal
 234      character*50 :: ierrc(its:ite)
 235      real,    dimension (its:ite,kts:kte) ::                           &
 236        up_massentr,up_massdetr,up_massentro,up_massdetro
 237      real :: C_up,x_add,qaver
 238      real,    dimension (its:ite,kts:kte) :: dtempdz
 239      integer, dimension (its:ite,kts:kte) ::  k_inv_layers
 240      integer, dimension (its:ite) ::  start_level
 241      start_level(:)=0
 242      rand_vmas(:)=0.
 243      flux_tun=fluxtune
 244       do i=its,itf
 245         xland1(i)=int(xland(i)+.001) ! 1.
 246         ktopx(i)=0
 247         if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
 248             xland1(i)=0
 249 !            ierr(i)=100
 250         endif
 251         pre(i)=0.
 252         xmb_out(i)=0.
 253         cap_max_increment(i)=25.
 254         ierrc(i)=" "
 255         entr_rate(i) = 9.e-5 ! 1.75e-3 ! 1.2e-3 ! .2/50.
 256       enddo
 257 !
 258 !--- initial entrainment rate (these may be changed later on in the
 259 !--- program
 260 !
 261
 262 !
 263 !--- initial detrainmentrates
 264 !
 265       do k=kts,ktf
 266       do i=its,itf
 267         up_massentro(i,k)=0.
 268         up_massdetro(i,k)=0.
 269         z(i,k)=zo(i,k)
 270         xz(i,k)=zo(i,k)
 271         qrco(i,k)=0.
 272         pwo(i,k)=0.
 273         cd(i,k)=1.*entr_rate(i)
 274         dellaqc(i,k)=0.
 275         cupclw(i,k)=0.
 276       enddo
 277       enddo
 278 !
 279 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
 280 !
 281 !--- minimum depth (m), clouds must have
 282 !
 283 !
 284 !--- maximum depth (mb) of capping
 285 !--- inversion (larger cap = no convection)
 286 !
 287       cap_maxs=125.
 288       DO i=its,itf
 289         kbmax(i)=1
 290         aa0(i)=0.
 291         aa1(i)=0.
 292       enddo
 293       do i=its,itf
 294           cap_max(i)=cap_maxs
 295           ztexec(i)  = 0.
 296           zqexec(i)  = 0.
 297           zws(i)     = 0.
 298       enddo
 299       do i=its,itf
 300          !- buoyancy flux (H+LE)
 301          buo_flux= (hfx(i)/cp+0.608*t(i,1)*qfx(i)/xlv)/rho(i,1)
 302          pgeoh = zo(i,2)*g
 303          !-convective-scale velocity w*
 304          zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,2)*g/t(i,1))
 305          if(zws(i) > TINY(pgeoh)) then
 306           !-convective-scale velocity w*
 307           zws(i) = 1.2*zws(i)**.3333
 308           !- temperature excess
 309           ztexec(i)     = MAX(flux_tun(i)*hfx(i)/(rho(i,1)*zws(i)*cp),0.0)
 310           !- moisture  excess
 311           zqexec(i)     = MAX(flux_tun(i)*qfx(i)/xlv/(rho(i,1)*zws(i)),0.)
 312          endif
 313        !- zws for shallow convection closure (Grant 2001)
 314        !- height of the pbl
 315        zws(i) = max(0.,flux_tun(i)*0.41*buo_flux*zo(i,kpbl(i))*g/t(i,kpbl(i)))
 316        zws(i) = 1.2*zws(i)**.3333
 317        zws(i) = zws(i)*rho(i,kpbl(i)) !check if zrho is correct
 318
 319       enddo
 320
 321 !
 322 !--- max height(m) above ground where updraft air can originate
 323 !
 324       zkbmax=3000.
 325 !
 326 !--- calculate moist static energy, heights, qes
 327 !
 328       call cup_env(z,qes,he,hes,t,q,po,z1, &
 329            psur,ierr,tcrit,-1,   &
 330            itf,ktf, &
 331            its,ite, kts,kte)
 332       call cup_env(zo,qeso,heo,heso,tn,qo,po,z1, &
 333            psur,ierr,tcrit,-1,   &
 334            itf,ktf, &
 335            its,ite, kts,kte)
 336
 337 !
 338 !--- environmental values on cloud levels
 339 !
 340       call cup_env_clev(t,qes,q,he,hes,z,po,qes_cup,q_cup,he_cup, &
 341            hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
 342            ierr,z1,          &
 343            itf,ktf, &
 344            its,ite, kts,kte)
 345       call cup_env_clev(tn,qeso,qo,heo,heso,zo,po,qeso_cup,qo_cup, &
 346            heo_cup,heso_cup,zo_cup,po_cup,gammao_cup,tn_cup,psur,  &
 347            ierr,z1,          &
 348            itf,ktf, &
 349            its,ite, kts,kte)
 350       do i=its,itf
 351         if(ierr(i).eq.0)then
 352 !
 353       do k=kts,ktf
 354         if(zo_cup(i,k).gt.zkbmax+z1(i))then
 355           kbmax(i)=k
 356           go to 25
 357         endif
 358       enddo
 359  25   continue
 360 !
 361       kbmax(i)=min(kbmax(i),ktf/2)
 362       endif
 363       enddo
 364
 365 !
 366 !
 367 !
 368 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
 369 !
 370        DO 36 i=its,itf
 371          if(kpbl(i).gt.3)cap_max(i)=po_cup(i,kpbl(i))
 372          IF(ierr(I) == 0)THEN
 373           k22(i)=maxloc(HEO_CUP(i,2:kbmax(i)),1)
 374           k22(i)=max(2,k22(i))
 375           IF(K22(I).GT.KBMAX(i))then
 376            ierr(i)=2
 377            ierrc(i)="could not find k22"
 378            ktop(i)=0
 379            k22(i)=0
 380            kbcon(i)=0
 381          endif
 382          endif
 383  36   CONTINUE
 384 !
 385 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE  - KBCON
 386 !
 387       do i=its,itf
 388        if(ierr(I).eq.0)then
 389              x_add = xlv*zqexec(i)+cp*ztexec(i)
 390              call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
 391              call get_cloud_bc(kte,heo_cup(i,1:kte),hkbo(i),k22(i),x_add)
 392        endif ! ierr
 393       enddo
 394
 395 !JOE-Georg and Saulo's new idea:
 396       do i=its,itf
 397       do k=kts,ktf
 398           dbyo(i,k)= 0. !hkbo(i)-heso_cup(i,k)
 399       enddo
 400       enddo
 401
 402       call cup_kbcon(ierrc,cap_max_increment,5,k22,kbcon,heo_cup,heso_cup, &
 403            hkbo,ierr,kbmax,po_cup,cap_max, &
 404            ztexec,zqexec, &
 405            0,itf,ktf, &
 406            its,ite, kts,kte, &
 407            z_cup,entr_rate,heo,0)
 408 !--- get inversion layers for cloud tops
 409       call cup_minimi(HEso_cup,Kbcon,kbmax,kstabi,ierr,  &
 410            itf,ktf, &
 411            its,ite, kts,kte)
 412 !
 413       call get_inversion_layers(ierr,p_cup,t_cup,z_cup,q_cup,qes_cup,k_inv_layers,&
 414                            kbcon,kstabi,dtempdz,itf,ktf,its,ite, kts,kte)
 415 !
 416 !
 417       DO i=its,itf
 418          entr_rate_2d(i,:)=entr_rate(i)
 419          IF(ierr(I) == 0)THEN
 420             start_level(i)=k22(i)
 421             x_add = xlv*zqexec(i)+cp*ztexec(i)
 422             call get_cloud_bc(kte,he_cup (i,1:kte),hkb (i),k22(i),x_add)
 423             if(kbcon(i).gt.ktf-4)then
 424                 ierr(i)=231
 425             endif
 426             do k=kts,ktf
 427                frh = 2.*min(qo_cup(i,k)/qeso_cup(i,k),1.)
 428                entr_rate_2d(i,k)=entr_rate(i)*(2.3-frh)
 429                cd(i,k)=entr_rate_2d(i,k)
 430             enddo
 431 !
 432 ! first estimate for shallow convection
 433 !
 434             ktop(i)=1
 435 !            if(k_inv_layers(i,1).gt.0)then
 436 !!               ktop(i)=min(k_inv_layers(i,1),k_inv_layers(i,2))
 437             if(k_inv_layers(i,1).gt.0 .and.   &
 438                (po_cup(i,kbcon(i))-po_cup(i,k_inv_layers(i,1))).lt.200.)then
 439                ktop(i)=k_inv_layers(i,1)
 440             else
 441                do k=kbcon(i)+1,ktf
 442                   if((po_cup(i,kbcon(i))-po_cup(i,k)).gt.200.)then
 443                     ktop(i)=k
 444                     exit
 445                   endif
 446                enddo
 447             endif
 448          endif
 449       enddo
 450 ! get normalized mass flux profile
 451       call rates_up_pdf(rand_vmas,ipr,'shallow',ktop,ierr,po_cup,entr_rate_2d,hkbo,heo,heso_cup,zo_cup, &
 452            xland1,kstabi,k22,kbcon,its,ite,itf,kts,kte,ktf,zuo,kpbl,ktopx,ktopx,kbcon)
 453       do i=its,itf
 454         if(ierr(i).eq.0)then
 455 !           do k=maxloc(zuo(i,:),1),1,-1 ! ktop(i)-1,1,-1
 456 !             if(zuo(i,k).lt.1.e-6)then
 457 !               k22(i)=k+1
 458 !               start_level(i)=k22(i)
 459 !               exit
 460 !             endif
 461 !           enddo
 462            if(k22(i).gt.1)then
 463              do k=1,k22(i)-1
 464               zuo(i,k)=0.
 465               zu (i,k)=0.
 466               xzu(i,k)=0.
 467              enddo
 468            endif
 469            do k=maxloc(zuo(i,:),1),ktop(i)
 470              if(zuo(i,k).lt.1.e-6)then
 471                ktop(i)=k-1
 472                exit
 473              endif
 474            enddo
 475            do k=k22(i),ktop(i)
 476              xzu(i,k)= zuo(i,k)
 477               zu(i,k)= zuo(i,k)
 478            enddo
 479            do k=ktop(i)+1,ktf
 480              zuo(i,k)=0.
 481              zu (i,k)=0.
 482              xzu(i,k)=0.
 483            enddo
 484            k22(i)=max(2,k22(i))
 485         endif
 486       enddo
 487 !
 488 ! calculate mass entrainment and detrainment
 489 !
 490       CALL get_lateral_massflux(itf,ktf, its,ite, kts,kte &
 491                                 ,ierr,ktop,zo_cup,zuo,cd,entr_rate_2d        &
 492                                 ,up_massentro, up_massdetro ,up_massentr, up_massdetr &
 493                                 ,'shallow',kbcon,k22)
 494
 495       do k=kts,ktf
 496       do i=its,itf
 497          hc(i,k)=0.
 498          qco(i,k)=0.
 499          qrco(i,k)=0.
 500          DBY(I,K)=0.
 501          hco(i,k)=0.
 502          DBYo(I,K)=0.
 503       enddo
 504       enddo
 505       do i=its,itf
 506        IF(ierr(I) /= 0) cycle
 507          do k=1,start_level(i)-1
 508             hc(i,k)=he_cup(i,k)
 509             hco(i,k)=heo_cup(i,k)
 510          enddo
 511          k=start_level(i)
 512          hc(i,k)=hkb(i)
 513          hco(i,k)=hkbo(i)
 514       enddo
 515 !
 516 !
 517       do 42 i=its,itf
 518         dbyt(i,:)=0.
 519         IF(ierr(I) /= 0) cycle
 520          do k=start_level(i)+1,ktop(i)
 521           hc(i,k)=(hc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)*hc(i,k-1)+ &
 522                          up_massentr(i,k-1)*he(i,k-1))   /            &
 523                          (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
 524           dby(i,k)=max(0.,hc(i,k)-hes_cup(i,k))
 525           hco(i,k)=(hco(i,k-1)*zuo(i,k-1)-.5*up_massdetro(i,k-1)*hco(i,k-1)+ &
 526                          up_massentro(i,k-1)*heo(i,k-1))   /            &
 527                          (zuo(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
 528           dbyo(i,k)=hco(i,k)-heso_cup(i,k)
 529           DZ=Zo_cup(i,K+1)-Zo_cup(i,K)
 530           dbyt(i,k)=dbyt(i,k-1)+dbyo(i,k)*dz
 531          enddo
 532        ki=maxloc(dbyt(i,:),1)
 533        if(ktop(i).gt.ki+1)then
 534          ktop(i)=ki+1
 535          zuo(i,ktop(i)+1:ktf)=0.
 536          zu(i,ktop(i)+1:ktf)=0.
 537          cd(i,ktop(i)+1:ktf)=0.
 538          up_massdetro(i,ktop(i))=zuo(i,ktop(i))
 539 !         up_massentro(i,ktop(i))=0.
 540          up_massentro(i,ktop(i):ktf)=0.
 541          up_massdetro(i,ktop(i)+1:ktf)=0.
 542          entr_rate_2d(i,ktop(i)+1:ktf)=0.
 543
 544 !         ierr(i)=423
 545        endif
 546
 547          if(ktop(i).lt.kbcon(i)+1)then
 548             ierr(i)=5
 549             ierrc(i)='ktop is less than kbcon+1'
 550              go to 42
 551          endif
 552          if(ktop(i).gt.ktf-2)then
 553              ierr(i)=5
 554              ierrc(i)="ktop is larger than ktf-2"
 555              go to 42
 556          endif
 557 !
 558          call get_cloud_bc(kte,qo_cup (i,1:kte),qaver,k22(i))
 559          qaver = qaver + zqexec(i)
 560          do k=1,start_level(i)-1
 561            qco (i,k)= qo_cup(i,k)
 562          enddo
 563          k=start_level(i)
 564          qco (i,k)= qaver
 565 !
 566          do k=start_level(i)+1,ktop(i)
 567           trash=QESo_cup(I,K)+(1./XLV)*(GAMMAo_cup(i,k) &
 568                 /(1.+GAMMAo_cup(i,k)))*DBYo(I,K)
 569           !- total water liq+vapour
 570           trash2  = qco(i,k-1) ! +qrco(i,k-1)
 571           qco (i,k)=   (trash2* ( zuo(i,k-1)-0.5*up_massdetr(i,k-1)) + &
 572                        up_massentr(i,k-1)*qo(i,k-1))   /            &
 573                        (zuo(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1))
 574
 575           if(qco(i,k)>=trash ) then
 576               DZ=Z_cup(i,K)-Z_cup(i,K-1)
 577               ! cloud liquid water
 578               qrco(i,k)= (qco(i,k)-trash)/(1.+(c0_shal+c1_shal)*dz)
 579 !              qrco(i,k)= (qco(i,k)-trash)/(1.+c0_shal*dz)
 580               pwo(i,k)=c0_shal*dz*qrco(i,k)*zuo(i,k)
 581               ! cloud water vapor
 582               qco (i,k)= trash+qrco(i,k)
 583
 584           else
 585               qrco(i,k)= 0.0
 586           endif
 587           cupclw(i,k)=qrco(i,k)
 588          enddo
 589          trash=0.
 590          trash2=0.
 591          do k=k22(i)+1,ktop(i)
 592           dp=100.*(po_cup(i,k)-po_cup(i,k+1))
 593           cnvwt(i,k)=zuo(i,k)*cupclw(i,k)*g/dp
 594           trash2=trash2+entr_rate_2d(i,k)
 595           qco(i,k)=qco(i,k)-qrco(i,k)
 596          enddo
 597          do k=k22(i)+1,max(kbcon(i),k22(i)+1)
 598           trash=trash+entr_rate_2d(i,k)
 599          enddo
 600          do k=ktop(i)+1,ktf-1
 601            hc  (i,k)=hes_cup (i,k)
 602            hco (i,k)=heso_cup(i,k)
 603            qco (i,k)=qeso_cup(i,k)
 604            qrco(i,k)=0.
 605            dby (i,k)=0.
 606            dbyo(i,k)=0.
 607            zu  (i,k)=0.
 608            xzu (i,k)=0.
 609            zuo (i,k)=0.
 610          enddo
 611  42 continue
 612 !
 613 !--- calculate workfunctions for updrafts
 614 !
 615       IF(MAKE_CALC_FOR_XK) THEN
 616         call cup_up_aa0(aa0,z,zu,dby,GAMMA_CUP,t_cup, &
 617             kbcon,ktop,ierr,           &
 618             itf,ktf, its,ite, kts,kte)
 619         call cup_up_aa0(aa1,zo,zuo,dbyo,GAMMAo_CUP,tn_cup, &
 620             kbcon,ktop,ierr,           &
 621             itf,ktf, its,ite, kts,kte)
 622         do i=its,itf
 623           if(ierr(i) == 0)then
 624            if(aa1(i) <= 0.)then
 625                ierr(i)=17
 626                ierrc(i)="cloud work function zero"
 627            endif
 628          endif
 629        enddo
 630       ENDIF
 631 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 632 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 633
 634 !
 635 !--- change per unit mass that a model cloud would modify the environment
 636 !
 637 !--- 1. in bottom layer
 638 !
 639       do k=kts,kte
 640        do i=its,itf
 641         dellah(i,k)=0.
 642         dellaq(i,k)=0.
 643        enddo
 644       enddo
 645 !
 646 !----------------------------------------------  cloud level ktop
 647 !
 648 !- - - - - - - - - - - - - - - - - - - - - - - - model level ktop-1
 649 !      .               .                 .
 650 !      .               .                 .
 651 !      .               .                 .
 652 !      .               .                 .
 653 !      .               .                 .
 654 !      .               .                 .
 655 !
 656 !----------------------------------------------  cloud level k+2
 657 !
 658 !- - - - - - - - - - - - - - - - - - - - - - - - model level k+1
 659 !
 660 !----------------------------------------------  cloud level k+1
 661 !
 662 !- - - - - - - - - - - - - - - - - - - - - - - - model level k
 663 !
 664 !----------------------------------------------  cloud level k
 665 !
 666 !      .               .                 .
 667 !      .               .                 .
 668 !      .               .                 .
 669 !      .               .                 .
 670 !      .               .                 .
 671 !      .               .                 .
 672 !      .               .                 .
 673 !      .               .                 .
 674 !      .               .                 .
 675 !      .               .                 .
 676 !
 677 !----------------------------------------------  cloud level 3
 678 !
 679 !- - - - - - - - - - - - - - - - - - - - - - - - model level 2
 680 !
 681 !----------------------------------------------  cloud level 2
 682 !
 683 !- - - - - - - - - - - - - - - - - - - - - - - - model level 1
 684       trash2=0.
 685       do i=its,itf
 686         if(ierr(i).eq.0)then
 687          do k=k22(i),ktop(i)
 688             ! entrainment/detrainment for updraft
 689             entup=up_massentro(i,k)
 690             detup=up_massdetro(i,k)
 691             totmas=detup-entup+zuo(i,k+1)-zuo(i,k)
 692             if(abs(totmas).gt.1.e-6)then
 693                write(0,*)'*********************',i,k,totmas
 694                write(0,*)k22(i),kbcon(i),ktop(i)
 695             endif
 696             dp=100.*(po_cup(i,k)-po_cup(i,k+1))
 697             dellah(i,k) =-(zuo(i,k+1)*(hco(i,k+1)-heo_cup(i,k+1) )-     &
 698                            zuo(i,k  )*(hco(i,k  )-heo_cup(i,k  ) ))*g/dp
 699
 700             !-- take out cloud liquid water for detrainment
 701             dz=zo_cup(i,k+1)-zo_cup(i,k)
 702             if(k.lt.ktop(i))then
 703              dellaqc(i,k)= zuo(i,k)*c1_shal*qrco(i,k)*dz/dp*g !  detup*0.5*(qrco(i,k+1)+qrco(i,k)) *g/dp
 704             else
 705              dellaqc(i,k)=   detup*qrco(i,k) *g/dp
 706             endif
 707
 708             !-- condensation source term = detrained + flux divergence of
 709             !-- cloud liquid water (qrco)
 710             C_up = dellaqc(i,k)+(zuo(i,k+1)* qrco(i,k+1) -       &
 711                                   zuo(i,k  )* qrco(i,k  )  )*g/dp
 712 !            C_up = dellaqc(i,k)
 713             !-- water vapor budget (flux divergence of Q_up-Q_env - condensation
 714             !term)
 715             dellaq(i,k) =-(zuo(i,k+1)*(qco(i,k+1)-qo_cup(i,k+1) ) -      &
 716                            zuo(i,k  )*(qco(i,k  )-qo_cup(i,k  ) ) )*g/dp &
 717                            - C_up - 0.5*(pwo (i,k)+pwo (i,k+1))*g/dp
 718           enddo
 719         endif
 720       enddo
 721
 722 !
 723 !--- using dellas, calculate changed environmental profiles
 724 !
 725       mbdt=.5 !3.e-4
 726
 727       do k=kts,ktf
 728        do i=its,itf
 729          dellat(i,k)=0.
 730          if(ierr(i)/=0)cycle
 731          xhe(i,k)=dellah(i,k)*mbdt+heo(i,k)
 732          xq (i,k)=max(1.e-16,(dellaq(i,k)+dellaqc(i,k))*mbdt+qo(i,k))
 733          dellat(i,k)=(1./cp)*(dellah(i,k)-xlv*(dellaq(i,k)))
 734          xt (i,k)= (-dellaqc(i,k)*xlv/cp+dellat(i,k))*mbdt+tn(i,k)
 735          xt (i,k)=  max(190.,xt(i,k))
 736
 737        enddo
 738       enddo
 739       do i=its,itf
 740        if(ierr(i).eq.0)then
 741 !        xhkb(i)=hkbo(i)+(dellah(i,k22(i)))*mbdt
 742         xhe(i,ktf)=heo(i,ktf)
 743         xq(i,ktf)=qo(i,ktf)
 744         xt(i,ktf)=tn(i,ktf)
 745        endif
 746       enddo
 747 !
 748 !
 749      IF(MAKE_CALC_FOR_XK) THEN
 750 !
 751 !--- calculate moist static energy, heights, qes
 752 !
 753       call cup_env(xz,xqes,xhe,xhes,xt,xq,po,z1, &
 754            psur,ierr,tcrit,-1,   &
 755            itf,ktf, &
 756            its,ite, kts,kte)
 757 !
 758 !--- environmental values on cloud levels
 759 !
 760       call cup_env_clev(xt,xqes,xq,xhe,xhes,xz,po,xqes_cup,xq_cup, &
 761            xhe_cup,xhes_cup,xz_cup,po_cup,gamma_cup,xt_cup,psur,   &
 762            ierr,z1,          &
 763            itf,ktf, &
 764            its,ite, kts,kte)
 765 !
 766 !
 767 !**************************** static control
 768       do k=kts,ktf
 769       do i=its,itf
 770          xhc(i,k)=0.
 771          xDBY(I,K)=0.
 772       enddo
 773       enddo
 774       do i=its,itf
 775         if(ierr(i).eq.0)then
 776          x_add = xlv*zqexec(i)+cp*ztexec(i)
 777          call get_cloud_bc(kte,xhe_cup (i,1:kte),xhkb (i),k22(i),x_add)
 778          do k=1,start_level(i)-1
 779             xhc(i,k)=xhe_cup(i,k)
 780          enddo
 781          k=start_level(i)
 782          xhc(i,k)=xhkb(i)
 783         endif !ierr
 784       enddo
 785 !
 786 !
 787       do i=its,itf
 788        if(ierr(i).eq.0)then
 789         xzu(i,1:ktf)=zuo(i,1:ktf)
 790         do k=start_level(i)+1,ktop(i)
 791          xhc(i,k)=(xhc(i,k-1)*xzu(i,k-1)-.5*up_massdetro(i,k-1)*xhc(i,k-1)+ &
 792                           up_massentro(i,k-1)*xhe(i,k-1))   /            &
 793                           (xzu(i,k-1)-.5*up_massdetro(i,k-1)+up_massentro(i,k-1))
 794          xdby(i,k)=xhc(i,k)-xhes_cup(i,k)
 795         enddo
 796         do k=ktop(i)+1,ktf
 797            xHC (i,K)=xhes_cup(i,k)
 798            xDBY(I,K)=0.
 799            xzu (i,k)=0.
 800         enddo
 801        endif
 802       enddo
 803
 804 !
 805 !--- workfunctions for updraft
 806 !
 807       call cup_up_aa0(xaa0,xz,xzu,xdby,GAMMA_CUP,xt_cup, &
 808            kbcon,ktop,ierr,           &
 809            itf,ktf, &
 810            its,ite, kts,kte)
 811 !
 812      ENDIF
 813 !
 814 !
 815 ! now for shallow forcing
 816 !
 817        do i=its,itf
 818         xmb(i)=0.
 819         xff_shal(1:3)=0.
 820         if(ierr(i).eq.0)then
 821           xmbmax(i)=1.0
 822 !         xmbmax(i)=100.*(p(i,kbcon(i))-p(i,kbcon(i)+1))/(g*dtime)
 823 !
 824 !-stabilization closure
 825           xkshal=(xaa0(i)-aa1(i))/mbdt
 826              if(xkshal.le.0.and.xkshal.gt.-.01*mbdt) &
 827                            xkshal=-.01*mbdt
 828              if(xkshal.gt.0.and.xkshal.lt.1.e-2) &
 829                            xkshal=1.e-2
 830
 831           xff_shal(1)=max(0.,-(aa1(i)-aa0(i))/(xkshal*dtime))
 832 !
 833 !- closure from Grant (2001)
 834           xff_shal(2)=.03*zws(i)
 835 !- boundary layer QE closure
 836           blqe=0.
 837           trash=0.
 838           do k=1,kpbl(i)
 839                 blqe=blqe+100.*dhdt(i,k)*(po_cup(i,k)-po_cup(i,k+1))/g
 840           enddo
 841           trash=max((hc(i,kbcon(i))-he_cup(i,kbcon(i))),1.e1)
 842           xff_shal(3)=max(0.,blqe/trash)
 843           xff_shal(3)=min(xmbmax(i),xff_shal(3))
 844 !- average
 845           xmb(i)=(xff_shal(1)+xff_shal(2)+xff_shal(3))/3.
 846           xmb(i)=min(xmbmax(i),xmb(i))
 847           if(ichoice > 0)xmb(i)=min(xmbmax(i),xff_shal(ichoice))
 848           if(xmb(i) <= 0.)then
 849              ierr(i)=21
 850              ierrc(i)="21"
 851           endif
 852         endif
 853         if(ierr(i).ne.0)then
 854            k22  (i)=0
 855            kbcon(i)=0
 856            ktop (i)=0
 857            xmb  (i)=0.
 858            outt (i,:)=0.
 859            outq (i,:)=0.
 860            outqc(i,:)=0.
 861         else if(ierr(i).eq.0)then
 862           xmb_out(i)=xmb(i)
 863 !
 864 ! final tendencies
 865 !
 866           pre(i)=0.
 867           do k=2,ktop(i)
 868            outt (i,k)= dellat (i,k)*xmb(i)
 869            outq (i,k)= dellaq (i,k)*xmb(i)
 870            outqc(i,k)= dellaqc(i,k)*xmb(i)
 871            pre  (i)  = pre(i)+pwo(i,k)*xmb(i)
 872           enddo
 873         endif
 874        enddo
 875
 876 #if ( WRF_CHEM == 1 )
 877 !--- calculate in-cloud/updraft air temperature
 878     do i=its,itf
 879       if (ierr(i)==0) then
 880         do k=kts,ktf
 881           tempco(i,k)=(1./cp)*(hco(i,k)-g*zo_cup(i,k)-xlv*qco(i,k))
 882         enddo !k
 883       else
 884         do k=kts,ktf
 885           tempco(i,k)=tn_cup(i,k)
 886         enddo !k
 887       endif !ierr
 888     enddo !i
 889     do i=its,itf
 890       pwevo(i)=0.
 891       pwavo(i)=0.
 892       jmin(i)=0
 893       edto(i)=0.
 894       do k=kts,kte
 895         zdo(i,k)=0.
 896         pwdo(i,k)=0.
 897         dd_massentro(i,k)=0.
 898         dd_massdetro(i,k)=0.
 899         clw_all(i,k)=0.
 900       enddo
 901     enddo
 902
 903     if ((chem_conv_tr>0).and.(chemopt>0)) then
 904       call ctrans_gf(numgas,num_chem,chem2d,chemopt,0  &
 905                      ,outchemt,conv_tr_wetscav,conv_tr_aqchem &
 906                      ,po,po_cup,zo_cup &
 907                      ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
 908                      ,up_massentro,up_massdetro &
 909                      ,dd_massentro,dd_massdetro &
 910                      ,tempco,clw_all  &
 911                      ,ktop,k22,kbcon,jmin  &
 912                      ,xmb,ierr,edto  &
 913                      ,itf,ktf,its,ite,kts,kte  &
 914                      ,1)
 915     endif
 916     if ((chem_conv_tr>0).and.(traceropt>0)) then
 917       call ctrans_gf(0,num_tracer,tracer2d,0,traceropt  &
 918                      ,outtracert,0,0        &
 919                      ,po,po_cup,zo_cup &
 920                      ,zuo,zdo,pwo,pwdo,pwevo,pwavo &
 921                      ,up_massentro,up_massdetro &
 922                      ,dd_massentro,dd_massdetro &
 923                      ,tempco,clw_all  &
 924                      ,ktop,k22,kbcon,jmin  &
 925                      ,xmb,ierr,edto  &
 926                      ,itf,ktf,its,ite,kts,kte  &
 927                      ,1)
 928     endif
 929 #endif
 930 !
 931 ! done shallow
 932 !--------------------------done------------------------------
 933 !
 934
 935    END SUBROUTINE CUP_gf_sh
 936 END MODULE module_cu_gf_sh