phys/module_cu_tiedtke.F

   1 !-----------------------------------------------------------------------
   2 !
   3 !wrf:model_layer:physics
   4 !
   5 !####################tiedtke scheme#########################
   6 !   taken from the IPRC IRAM - Yuqing Wang, university of hawaii
   7 !   added by Chunxi Zhang and Yuqing Wang to wrf3.2, may, 2010
   8 !   refenrence: Tiedtke (1989, mwr, 117, 1779-1800)
   9 !               Nordeng, t.e., (1995), cape closure and organized entrainment/detrainment
  10 !               Yuqing Wang et al. (2003,j. climate, 16, 1721-1738) for improvements
  11 !                                                  for cloud top detrainment
  12 !                       (2004, mon. wea. rev., 132, 274-296), improvements for pbl clouds
  13 !                       (2007,mon. wea. rev., 135, 567-585), diurnal cycle of precipitation
  14 !###########################################################
  15 module module_cu_tiedtke
  16 !
  17      use module_model_constants, only:rd=>r_d, rv=>r_v, &
  18    &       cpd=>cp, alv=>xlv, als=>xls, alf=>xlf, g
  19
  20      implicit none
  21
  22      real :: rcpd,vtmpc1,t000,hgfr,rhoh2o,tmelt, &
  23              c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg
  24
  25      real :: entrpen,entrscv,entrmid,entrdd,cmfctop,rhm,rhc,    &
  26              cmfcmax,cmfcmin,cmfdeps,cprcon,crirh,zbuo0,  &
  27              fdbk,ztau
  28
  29      real :: cevapcu1, cevapcu2, zdnoprc
  30
  31      parameter(                     &
  32       rcpd=1.0/cpd,                 &
  33       rhoh2o=1.0e03,                &
  34       tmelt=273.16,                 &
  35       t000= 273.15,                 &
  36       hgfr= 233.15,                 &
  37       zrg=1.0/g,                    &
  38       c1es=610.78,                  &
  39       c2es=c1es*rd/rv,              &
  40       c3les=17.269,                 &
  41       c4les=35.86,                  &
  42       c5les=c3les*(tmelt-c4les),    &
  43       c3ies=21.875,                 &
  44       c4ies=7.66,                   &
  45       c5ies=c3ies*(tmelt-c4ies),    &
  46       vtmpc1=rv/rd-1.0,             &
  47       cevapcu1=1.93e-6*261.0*0.5/g, &
  48       cevapcu2=1.e3/(38.3*0.293) )
  49
  50
  51 !                specify parameters for massflux-scheme
  52 !                  --------------------------------------
  53 !                   these are tunable parameters
  54 !
  55 !     entrpen: average entrainment rate for penetrative convection
  56 !     -------
  57 !
  58       parameter(entrpen=1.0e-4)
  59 !
  60 !     entrscv: average entrainment rate for shallow convection
  61 !     -------
  62 !
  63       parameter(entrscv=1.2e-3)
  64 !
  65 !     entrmid: average entrainment rate for midlevel convection
  66 !     -------
  67 !
  68       parameter(entrmid=1.0e-4)
  69 !
  70 !     entrdd: average entrainment rate for downdrafts
  71 !     ------
  72 !
  73       parameter(entrdd =2.0e-4)
  74 !
  75 !     cmfctop:   relative cloud massflux at level above nonbuoyancy level
  76 !     -------
  77 !
  78       parameter(cmfctop=0.30)
  79 !
  80 !     cmfcmax:   maximum massflux value allowed for updrafts etc
  81 !     -------
  82 !
  83       parameter(cmfcmax=1.0)
  84 !
  85 !     cmfcmin:   minimum massflux value (for safety)
  86 !     -------
  87 !
  88       parameter(cmfcmin=1.e-10)
  89 !
  90 !     cmfdeps:   fractional massflux for downdrafts at lfs
  91 !     -------
  92 !
  93       parameter(cmfdeps=0.30)
  94 !
  95 !     cprcon:  coefficients for determining conversion from cloud water
  96 !
  97       parameter(cprcon = 1.1e-3/g)
  98 !
  99 !     zdnoprc: the pressure depth below which no precipitation
 100 !
 101       parameter(zdnoprc = 1.5e4)
 102 !--------------------
 103       parameter(rhc=0.80,rhm=1.0,zbuo0=0.50)
 104 !--------------------
 105       parameter(crirh=0.70,fdbk = 1.0,ztau = 2400.0)
 106 !--------------------
 107       logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv
 108       parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.)
 109 !--------------------
 110 !#################### end of variables definition##########################
 111 !-----------------------------------------------------------------------
 112 !
 113 contains
 114 !-----------------------------------------------------------------------
 115       subroutine cu_tiedtke(                                    &
 116                  dt,itimestep,stepcu                            &
 117                 ,raincv,pratec,qfx,znu                          &
 118                 ,u3d,v3d,w,t3d,qv3d,qc3d,qi3d,pi3d,rho3d        &
 119                 ,qvften,qvpblten                                &
 120                 ,dz8w,pcps,p8w,xland,cu_act_flag                &
 121                 ,ids,ide, jds,jde, kds,kde                      &
 122                 ,ims,ime, jms,jme, kms,kme                      &
 123                 ,its,ite, jts,jte, kts,kte                      &
 124                 ,rthcuten,rqvcuten,rqccuten,rqicuten            &
 125                 ,rucuten, rvcuten                               &
 126                 ,f_qv    ,f_qc    ,f_qr    ,f_qi    ,f_qs       &
 127                                                                 )
 128
 129 !-------------------------------------------------------------------
 130       implicit none
 131 !-------------------------------------------------------------------
 132 !-- u3d         3d u-velocity interpolated to theta points (m/s)
 133 !-- v3d         3d v-velocity interpolated to theta points (m/s)
 134 !-- th3d        3d potential temperature (k)
 135 !-- t3d         temperature (k)
 136 !-- qv3d        3d water vapor mixing ratio (kg/kg)
 137 !-- qc3d        3d cloud mixing ratio (kg/kg)
 138 !-- qi3d        3d ice mixing ratio (kg/kg)
 139 !-- rho3d       3d air density (kg/m^3)
 140 !-- p8w         3d hydrostatic pressure at full levels (pa)
 141 !-- pcps        3d hydrostatic pressure at half levels (pa)
 142 !-- pi3d        3d exner function (dimensionless)
 143 !-- rthcuten      theta tendency due to
 144 !                 cumulus scheme precipitation (k/s)
 145 !-- rucuten       u wind tendency due to
 146 !                 cumulus scheme precipitation (k/s)
 147 !-- rvcuten       v wind tendency due to
 148 !                 cumulus scheme precipitation (k/s)
 149 !-- rqvcuten      qv tendency due to
 150 !                 cumulus scheme precipitation (kg/kg/s)
 151 !-- rqrcuten      qr tendency due to
 152 !                 cumulus scheme precipitation (kg/kg/s)
 153 !-- rqccuten      qc tendency due to
 154 !                 cumulus scheme precipitation (kg/kg/s)
 155 !-- rqscuten      qs tendency due to
 156 !                 cumulus scheme precipitation (kg/kg/s)
 157 !-- rqicuten      qi tendency due to
 158 !                 cumulus scheme precipitation (kg/kg/s)
 159 !-- rainc         accumulated total cumulus scheme precipitation (mm)
 160 !-- raincv        cumulus scheme precipitation (mm)
 161 !-- pratec        precipitiation rate from cumulus scheme (mm/s)
 162 !-- dz8w        dz between full levels (m)
 163 !-- qfx         upward moisture flux at the surface (kg/m^2/s)
 164 !-- dt          time step (s)
 165 !-- ids         start index for i in domain
 166 !-- ide         end index for i in domain
 167 !-- jds         start index for j in domain
 168 !-- jde         end index for j in domain
 169 !-- kds         start index for k in domain
 170 !-- kde         end index for k in domain
 171 !-- ims         start index for i in memory
 172 !-- ime         end index for i in memory
 173 !-- jms         start index for j in memory
 174 !-- jme         end index for j in memory
 175 !-- kms         start index for k in memory
 176 !-- kme         end index for k in memory
 177 !-- its         start index for i in tile
 178 !-- ite         end index for i in tile
 179 !-- jts         start index for j in tile
 180 !-- jte         end index for j in tile
 181 !-- kts         start index for k in tile
 182 !-- kte         end index for k in tile
 183 !-------------------------------------------------------------------
 184       integer, intent(in) ::            ids,ide, jds,jde, kds,kde,      &
 185                                         ims,ime, jms,jme, kms,kme,      &
 186                                         its,ite, jts,jte, kts,kte,      &
 187                                         itimestep,                      &
 188                                         stepcu
 189
 190       real,    intent(in) ::                                            &
 191                                         dt
 192
 193
 194       real,    dimension(ims:ime, jms:jme), intent(in) ::               &
 195                                         xland
 196
 197       real,    dimension(ims:ime, jms:jme), intent(inout) ::            &
 198                                         raincv, pratec
 199
 200       logical, dimension(ims:ime,jms:jme), intent(inout) ::             &
 201                                         cu_act_flag
 202
 203
 204       real,    dimension(ims:ime, kms:kme, jms:jme), intent(in) ::      &
 205                                         dz8w,                           &
 206                                         p8w,                            &
 207                                         pcps,                           &
 208                                         pi3d,                           &
 209                                         qc3d,                           &
 210                                         qvften,                         &
 211                                         qvpblten,                       &
 212                                         qi3d,                           &
 213                                         qv3d,                           &
 214                                         rho3d,                          &
 215                                         t3d,                            &
 216                                         u3d,                            &
 217                                         v3d,                            &
 218                                         w
 219
 220 !--------------------------- optional vars ----------------------------
 221
 222       real, dimension(ims:ime, kms:kme, jms:jme),                       &
 223                optional, intent(inout) ::                               &
 224                                         rqccuten,                       &
 225                                         rqicuten,                       &
 226                                         rqvcuten,                       &
 227                                         rthcuten,                       &
 228                                         rucuten,                        &
 229                                         rvcuten
 230
 231 !
 232 ! flags relating to the optional tendency arrays declared above
 233 ! models that carry the optional tendencies will provdide the
 234 ! optional arguments at compile time; these flags all the model
 235 ! to determine at run-time whether a particular tracer is in
 236 ! use or not.
 237 !
 238      logical, optional ::                                    &
 239                                                    f_qv      &
 240                                                   ,f_qc      &
 241                                                   ,f_qr      &
 242                                                   ,f_qi      &
 243                                                   ,f_qs
 244
 245 !--------------------------- local vars ------------------------------
 246
 247       real,    dimension(ims:ime, jms:jme) ::                           &
 248                                         qfx
 249
 250       real      ::                                      &
 251                                         delt,                           &
 252                                         rdelt
 253
 254       real     , dimension(its:ite) ::                  &
 255                                         rcs,                            &
 256                                         rn,                             &
 257                                         evap
 258       integer  , dimension(its:ite) ::  slimsk
 259
 260
 261       real     , dimension(its:ite, kts:kte+1) ::       &
 262                                         prsi
 263
 264       real     , dimension(its:ite, kts:kte) ::         &
 265                                         del,                            &
 266                                         dot,                            &
 267                                         phil,                           &
 268                                         prsl,                           &
 269                                         q1,                             &
 270                                         q2,                             &
 271                                         q3,                             &
 272                                         q1b,                            &
 273                                         q1bl,                           &
 274                                         q11,                            &
 275                                         q12,                            &
 276                                         t1,                             &
 277                                         u1,                             &
 278                                         v1,                             &
 279                                         zi,                             &
 280                                         zl,                             &
 281                                         omg,                            &
 282                                         ght
 283
 284       integer, dimension(its:ite) ::                                    &
 285                                         kbot,                           &
 286                                         ktop
 287
 288       integer ::                                                        &
 289                                         i,                              &
 290                                         im,                             &
 291                                         j,                              &
 292                                         k,                              &
 293                                         km,                             &
 294                                         kp,                             &
 295                                         kx
 296
 297
 298       logical :: run_param , doing_adapt_dt , decided
 299
 300 !-------other local variables----
 301       integer,dimension( its:ite ) :: ktype
 302       real, dimension( kts:kte )   :: sig1      ! half sigma levels
 303       real, dimension( kms:kme )   :: znu
 304       integer                      :: zz
 305 !-----------------------------------------------------------------------
 306
 307       do j=jts,jte
 308          do i=its,ite
 309             cu_act_flag(i,j)=.true.
 310          enddo
 311       enddo
 312
 313       im=ite-its+1
 314       kx=kte-kts+1
 315       delt=dt*stepcu
 316       rdelt=1./delt
 317
 318 !-------------  j loop (outer) --------------------------------------------------
 319
 320    do j=jts,jte
 321
 322 ! --------------- compute zi and zl -----------------------------------------
 323       do i=its,ite
 324         zi(i,kts)=0.0
 325       enddo
 326
 327       do k=kts+1,kte
 328         km=k-1
 329         do i=its,ite
 330           zi(i,k)=zi(i,km)+dz8w(i,km,j)
 331         enddo
 332       enddo
 333
 334       do k=kts+1,kte
 335         km=k-1
 336         do i=its,ite
 337           zl(i,km)=(zi(i,k)+zi(i,km))*0.5
 338         enddo
 339       enddo
 340
 341       do i=its,ite
 342         zl(i,kte)=2.*zi(i,kte)-zl(i,kte-1)
 343       enddo
 344
 345 ! --------------- end compute zi and zl -------------------------------------
 346       do i=its,ite
 347         slimsk(i)=int(abs(xland(i,j)-2.))
 348       enddo
 349
 350       do k=kts,kte
 351         kp=k+1
 352         do i=its,ite
 353           dot(i,k)=-0.5*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
 354         enddo
 355       enddo
 356
 357       do k=kts,kte
 358         zz = kte+1-k
 359         do i=its,ite
 360           u1(i,zz)=u3d(i,k,j)
 361           v1(i,zz)=v3d(i,k,j)
 362           t1(i,zz)=t3d(i,k,j)
 363           q1(i,zz)= qv3d(i,k,j)
 364           if(itimestep == 1) then
 365              q1b(i,zz)=0.
 366              q1bl(i,zz)=0.
 367           else
 368              q1b(i,zz)=qvften(i,k,j)
 369              q1bl(i,zz)=qvpblten(i,k,j)
 370           endif
 371           q2(i,zz)=qc3d(i,k,j)
 372           q3(i,zz)=qi3d(i,k,j)
 373           omg(i,zz)=dot(i,k)
 374           ght(i,zz)=zl(i,k)
 375           prsl(i,zz) = pcps(i,k,j)
 376         enddo
 377       enddo
 378
 379       do k=kts,kte+1
 380         zz = kte+2-k
 381         do i=its,ite
 382           prsi(i,zz) = p8w(i,k,j)
 383         enddo
 384       enddo
 385
 386       do k=kts,kte
 387          zz = kte+1-k
 388          sig1(zz) = znu(k)
 389       enddo
 390
 391 !###############before call tiecnv, we need evap########################
 392 !       evap is the vapor flux at the surface
 393 !########################################################################
 394 !
 395       do i=its,ite
 396         evap(i) = qfx(i,j)
 397       enddo
 398 !########################################################################
 399       call tiecnv(u1,v1,t1,q1,q2,q3,q1b,q1bl,ght,omg,prsl,prsi,evap,             &
 400                   rn,slimsk,ktype,im,kx,kx+1,sig1,delt)
 401
 402       do i=its,ite
 403          raincv(i,j)=rn(i)/stepcu
 404          pratec(i,j)=rn(i)/(stepcu * dt)
 405       enddo
 406
 407       do k=kts,kte
 408         zz = kte+1-k
 409         do i=its,ite
 410           rthcuten(i,k,j)=(t1(i,zz)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
 411           rqvcuten(i,k,j)=(q1(i,zz)-qv3d(i,k,j))*rdelt
 412           rucuten(i,k,j) =(u1(i,zz)-u3d(i,k,j))*rdelt
 413           rvcuten(i,k,j) =(v1(i,zz)-v3d(i,k,j))*rdelt
 414         enddo
 415       enddo
 416
 417       if(present(rqccuten))then
 418         if ( f_qc ) then
 419           do k=kts,kte
 420             zz = kte+1-k
 421             do i=its,ite
 422               rqccuten(i,k,j)=(q2(i,zz)-qc3d(i,k,j))*rdelt
 423             enddo
 424           enddo
 425         endif
 426       endif
 427
 428       if(present(rqicuten))then
 429         if ( f_qi ) then
 430           do k=kts,kte
 431             zz = kte+1-k
 432             do i=its,ite
 433               rqicuten(i,k,j)=(q3(i,zz)-qi3d(i,k,j))*rdelt
 434             enddo
 435           enddo
 436         endif
 437       endif
 438
 439
 440    enddo
 441
 442    end subroutine cu_tiedtke
 443
 444 !====================================================================
 445    subroutine tiedtkeinit(rthcuten,rqvcuten,rqccuten,rqicuten,          &
 446                      rucuten,rvcuten,                                   &
 447                      restart,p_qc,p_qi,p_first_scalar,                  &
 448                      allowed_to_read,                                   &
 449                      ids, ide, jds, jde, kds, kde,                      &
 450                      ims, ime, jms, jme, kms, kme,                      &
 451                      its, ite, jts, jte, kts, kte)
 452 !--------------------------------------------------------------------
 453    implicit none
 454 !--------------------------------------------------------------------
 455    logical , intent(in)           ::  allowed_to_read,restart
 456    integer , intent(in)           ::  ids, ide, jds, jde, kds, kde, &
 457                                       ims, ime, jms, jme, kms, kme, &
 458                                       its, ite, jts, jte, kts, kte
 459    integer , intent(in)           ::  p_first_scalar, p_qi, p_qc
 460
 461    real,     dimension( ims:ime , kms:kme , jms:jme ) , intent(out) ::  &
 462                                                               rthcuten, &
 463                                                               rqvcuten, &
 464                                                               rqccuten, &
 465                                                               rqicuten, &
 466                                                               rucuten,rvcuten
 467
 468    integer :: i, j, k, itf, jtf, ktf
 469
 470    jtf=min0(jte,jde-1)
 471    ktf=min0(kte,kde-1)
 472    itf=min0(ite,ide-1)
 473
 474    if(.not.restart)then
 475      do j=jts,jtf
 476      do k=kts,ktf
 477      do i=its,itf
 478        rthcuten(i,k,j)=0.
 479        rqvcuten(i,k,j)=0.
 480        rucuten(i,k,j)=0.
 481        rvcuten(i,k,j)=0.
 482      enddo
 483      enddo
 484      enddo
 485
 486      if (p_qc .ge. p_first_scalar) then
 487         do j=jts,jtf
 488         do k=kts,ktf
 489         do i=its,itf
 490            rqccuten(i,k,j)=0.
 491         enddo
 492         enddo
 493         enddo
 494      endif
 495
 496      if (p_qi .ge. p_first_scalar) then
 497         do j=jts,jtf
 498         do k=kts,ktf
 499         do i=its,itf
 500            rqicuten(i,k,j)=0.
 501         enddo
 502         enddo
 503         enddo
 504      endif
 505    endif
 506
 507       end subroutine tiedtkeinit
 508
 509 ! ------------------------------------------------------------------------
 510
 511 !------------this is the combined version for tiedtke---------------
 512 !----------------------------------------------------------------
 513 !  in this module only the mass flux convection scheme of the ecmwf is included
 514 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 515 !#############################################################
 516 !
 517 !             level 1 subroutines
 518 !
 519 !#############################################################
 520 !********************************************************
 521 !        subroutine tiecnv
 522 !********************************************************
 523       subroutine tiecnv(pu,pv,pt,pqv,pqc,pqi,pqvf,pqvbl,poz,pomg,  &
 524                pap,paph,evap,zprecc,lndj,ktype,lq,km,km1,sig1,dt)
 525 !-----------------------------------------------------------------
 526 !  this is the interface between the meso-scale model and the mass
 527 !  flux convection module
 528 !-----------------------------------------------------------------
 529       implicit none
 530
 531       real pu(lq,km),pv(lq,km),pt(lq,km),pqv(lq,km),pqvf(lq,km)
 532       real poz(lq,km),pomg(lq,km),evap(lq),zprecc(lq),pqvbl(lq,km)
 533
 534       real pum1(lq,km),    pvm1(lq,km),                             &
 535           ptte(lq,km),    pqte(lq,km),  pvom(lq,km),  pvol(lq,km),  &
 536           pverv(lq,km),   pgeo(lq,km),  pap(lq,km),   paph(lq,km1)
 537       real pqhfl(lq),      zqq(lq,km),   paprc(lq),    paprs(lq),   &
 538           prsfc(lq),      pssfc(lq),    paprsm(lq),   pcte(lq,km)
 539       real ztp1(lq,km),    zqp1(lq,km),  ztu(lq,km),   zqu(lq,km),  &
 540           zlu(lq,km),     zlude(lq,km), zmfu(lq,km),  zmfd(lq,km),  &
 541           zqsat(lq,km),   pqc(lq,km),   pqi(lq,km),   zrain(lq)
 542
 543       real sig(km1),sig1(km)
 544       integer icbot(lq),   ictop(lq),     ktype(lq),   lndj(lq)
 545       real  dt
 546       logical locum(lq)
 547
 548       real psheat,psrain,psevap,psmelt,psdiss,tt
 549       real ztmst,ztpp1,fliq,fice,ztc,zalf
 550       integer i,j,k,lq,lp,km,km1
 551 !     real tlucua
 552 !     external tlucua
 553
 554       ztmst=dt
 555 !  masv flux diagnostics.
 556
 557       psheat=0.0
 558       psrain=0.0
 559       psevap=0.0
 560       psmelt=0.0
 561       psdiss=0.0
 562       do  j=1,lq
 563         zrain(j)=0.0
 564         locum(j)=.false.
 565         prsfc(j)=0.0
 566         pssfc(j)=0.0
 567         paprc(j)=0.0
 568         paprs(j)=0.0
 569         paprsm(j)=0.0
 570         pqhfl(j)=evap(j)
 571      end do
 572
 573 !     convert model variables for mflux scheme
 574
 575       do  k=1,km
 576         do  j=1,lq
 577           ptte(j,k)=0.0
 578           pcte(j,k)=0.0
 579           pvom(j,k)=0.0
 580           pvol(j,k)=0.0
 581           ztp1(j,k)=pt(j,k)
 582           zqp1(j,k)=pqv(j,k)/(1.0+pqv(j,k))
 583           pum1(j,k)=pu(j,k)
 584           pvm1(j,k)=pv(j,k)
 585           pverv(j,k)=pomg(j,k)
 586           pgeo(j,k)=g*poz(j,k)
 587           tt=ztp1(j,k)
 588           zqsat(j,k)=tlucua(tt)/pap(j,k)
 589           zqsat(j,k)=min(0.5,zqsat(j,k))
 590           zqsat(j,k)=zqsat(j,k)/(1.-vtmpc1*zqsat(j,k))
 591           pqte(j,k)=pqvf(j,k)+pqvbl(j,k)
 592           zqq(j,k)=pqte(j,k)
 593         end do
 594       end do
 595 !
 596 !-----------------------------------------------------------------------
 597 !*    2.     call 'cumastr'(master-routine for cumulus parameterization)
 598 !
 599       call cumastr_new &
 600          (lq,       km,       km1,      km-1,    ztp1,   &
 601           zqp1,     pum1,     pvm1,     pverv,   zqsat,  &
 602           pqhfl,    ztmst,    pap,      paph,    pgeo,   &
 603           ptte,     pqte,     pvom,     pvol,    prsfc,  &
 604           pssfc,    paprc,    paprsm,   paprs,   locum,  &
 605           ktype,    icbot,    ictop,    ztu,     zqu,    &
 606           zlu,      zlude,    zmfu,     zmfd,    zrain,  &
 607           psrain,   psevap,   psheat,   psdiss,  psmelt, &
 608           pcte,     sig1,     lndj)
 609 !
 610 !     to include the cloud water and cloud ice detrained from convection
 611 !
 612       if(fdbk.ge.1.0e-9) then
 613       do k=1,km
 614       do j=1,lq
 615       if(pcte(j,k).gt.0.0) then
 616         ztpp1=pt(j,k)+ptte(j,k)*ztmst
 617         if(ztpp1.ge.t000) then
 618            fliq=1.0
 619            zalf=0.0
 620         else if(ztpp1.le.hgfr) then
 621            fliq=0.0
 622            zalf=alf
 623         else
 624            ztc=ztpp1-t000
 625            fliq=0.0059+0.9941*exp(-0.003102*ztc*ztc)
 626            zalf=alf
 627         endif
 628         fice=1.0-fliq
 629         pqc(j,k)=pqc(j,k)+fliq*pcte(j,k)*ztmst
 630         pqi(j,k)=pqi(j,k)+fice*pcte(j,k)*ztmst
 631         ptte(j,k)=ptte(j,k)-zalf*rcpd*fliq*pcte(j,k)
 632       endif
 633       end do
 634       end do
 635       endif
 636 !
 637       do k=1,km
 638         do j=1,lq
 639           pt(j,k)=ztp1(j,k)+ptte(j,k)*ztmst
 640           zqp1(j,k)=zqp1(j,k)+(pqte(j,k)-zqq(j,k))*ztmst
 641           pqv(j,k)=zqp1(j,k)/(1.0-zqp1(j,k))
 642         end do
 643       end do
 644       do j=1,lq
 645         zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst)
 646       end do
 647       if (lmfdudv) then
 648         do  k=1,km
 649           do j=1,lq
 650             pu(j,k)=pu(j,k)+pvom(j,k)*ztmst
 651             pv(j,k)=pv(j,k)+pvol(j,k)*ztmst
 652           end do
 653         end do
 654       endif
 655 !
 656       return
 657       end subroutine tiecnv
 658
 659 !#############################################################
 660 !
 661 !             level 2 subroutines
 662 !
 663 !#############################################################
 664 !***********************************************************
 665 !           subroutine cumastr_new
 666 !***********************************************************
 667       subroutine cumastr_new                             &
 668          (klon,     klev,     klevp1,   klevm1,   pten,  &
 669           pqen,     puen,     pven,     pverv,    pqsen, &
 670           pqhfl,    ztmst,    pap,      paph,     pgeo,  &
 671           ptte,     pqte,     pvom,     pvol,     prsfc, &
 672           pssfc,    paprc,    paprsm,   paprs,    ldcum, &
 673           ktype,    kcbot,    kctop,    ptu,      pqu,   &
 674           plu,      plude,    pmfu,     pmfd,     prain, &
 675           psrain,   psevap,   psheat,   psdiss,   psmelt,&
 676           pcte,     sig1,     lndj)
 677 !
 678 !***cumastr*  master routine for cumulus massflux-scheme
 679 !     m.tiedtke      e.c.m.w.f.     1986/1987/1989
 680 !***purpose
 681 !   -------
 682 !          this routine computes the physical tendencies of the
 683 !     prognostic variables t,q,u and v due to convective processes.
 684 !     processes considered are: convective fluxes, formation of
 685 !     precipitation, evaporation of falling rain below cloud base,
 686 !     saturated cumulus downdrafts.
 687 !***interface.
 688 !   ----------
 689 !          *cumastr* is called from *mssflx*
 690 !     the routine takes its input from the long-term storage
 691 !     t,q,u,v,phi and p and moisture tendencies.
 692 !     it returns its output to the same space
 693 !      1.modified tendencies of model variables
 694 !      2.rates of convective precipitation
 695 !        (used in subroutine surf)
 696 !      3.cloud base, cloud top and precip for radiation
 697 !        (used in subroutine cloud)
 698 !***method
 699 !   ------
 700 !     parameterization is done using a massflux-scheme.
 701 !        (1) define constants and parameters
 702 !        (2) specify values (t,q,qs...) at half levels and
 703 !            initialize updraft- and downdraft-values in 'cuini'
 704 !        (3) calculate cloud base in 'cubase'
 705 !            and specify cloud base massflux from pbl moisture budget
 706 !        (4) do cloud ascent in 'cuasc' in absence of downdrafts
 707 !        (5) do downdraft calculations:
 708 !              (a) determine values at lfs in 'cudlfs'
 709 !              (b) determine moist descent in 'cuddraf'
 710 !              (c) recalculate cloud base massflux considering the
 711 !                  effect of cu-downdrafts
 712 !        (6) do final cloud ascent in 'cuasc'
 713 !        (7) do final adjusments to convective fluxes in 'cuflx',
 714 !            do evaporation in subcloud layer
 715 !        (8) calculate increments of t and q in 'cudtdq'
 716 !        (9) calculate increments of u and v in 'cududv'
 717 !***externals.
 718 !   ----------
 719 !       cuini:  initializes values at vertical grid used in cu-parametr.
 720 !       cubase: cloud base calculation for penetr.and shallow convection
 721 !       cuasc:  cloud ascent for entraining plume
 722 !       cudlfs: determines values at lfs for downdrafts
 723 !       cuddraf:does moist descent for cumulus downdrafts
 724 !       cuflx:  final adjustments to convective fluxes (also in pbl)
 725 !       cudqdt: updates tendencies for t and q
 726 !       cududv: updates tendencies for u and v
 727 !***switches.
 728 !   --------
 729 !          lmfpen=.t.   penetrative convection is switched on
 730 !          lmfscv=.t.   shallow convection is switched on
 731 !          lmfmid=.t.   midlevel convection is switched on
 732 !          lmfdd=.t.    cumulus downdrafts switched on
 733 !          lmfdudv=.t.  cumulus friction switched on
 734 !***
 735 !     model parameters (defined in subroutine cuparam)
 736 !     ------------------------------------------------
 737 !     entrpen    entrainment rate for penetrative convection
 738 !     entrscv    entrainment rate for shallow convection
 739 !     entrmid    entrainment rate for midlevel convection
 740 !     entrdd     entrainment rate for cumulus downdrafts
 741 !     cmfctop    relative cloud massflux at level above nonbuoyancy
 742 !                level
 743 !     cmfcmax    maximum massflux value allowed for
 744 !     cmfcmin    minimum massflux value (for safety)
 745 !     cmfdeps    fractional massflux for downdrafts at lfs
 746 !     cprcon     coefficient for conversion from cloud water to rain
 747 !***reference.
 748 !   ----------
 749 !          paper on massflux scheme (tiedtke,1989)
 750 !-----------------------------------------------------------------
 751 !-------------------------------------------------------------------
 752       implicit none
 753 !-------------------------------------------------------------------
 754       integer   klon, klev, klevp1
 755       integer   klevm1
 756       real      ztmst
 757       real      psrain, psevap, psheat, psdiss, psmelt, zcons2
 758       integer   jk,jl,ikb
 759       real      zqumqe, zdqmin, zmfmax, zalvdcp, zqalv
 760       real      zhsat, zgam, zzz, zhhat, zbi, zro, zdz, zdhdz, zdepth
 761       real      zfac, zrh, zpbmpt, dept, zht, zeps
 762       integer   icum, itopm2
 763       real     pten(klon,klev),        pqen(klon,klev), &
 764               puen(klon,klev),        pven(klon,klev),  &
 765               ptte(klon,klev),        pqte(klon,klev),  &
 766               pvom(klon,klev),        pvol(klon,klev),  &
 767               pqsen(klon,klev),       pgeo(klon,klev),  &
 768               pap(klon,klev),         paph(klon,klevp1),&
 769               pverv(klon,klev),       pqhfl(klon)
 770       real     ptu(klon,klev),         pqu(klon,klev),  &
 771               plu(klon,klev),         plude(klon,klev), &
 772               pmfu(klon,klev),        pmfd(klon,klev),  &
 773               paprc(klon),            paprs(klon),      &
 774               paprsm(klon),           prain(klon),      &
 775               prsfc(klon),            pssfc(klon)
 776       real     ztenh(klon,klev),       zqenh(klon,klev),&
 777               zgeoh(klon,klev),       zqsenh(klon,klev),&
 778               ztd(klon,klev),         zqd(klon,klev),   &
 779               zmfus(klon,klev),       zmfds(klon,klev), &
 780               zmfuq(klon,klev),       zmfdq(klon,klev), &
 781               zdmfup(klon,klev),      zdmfdp(klon,klev),&
 782               zmful(klon,klev),       zrfl(klon),       &
 783               zuu(klon,klev),         zvu(klon,klev),   &
 784               zud(klon,klev),         zvd(klon,klev)
 785       real     zentr(klon),            zhcbase(klon),   &
 786               zmfub(klon),            zmfub1(klon),     &
 787               zdqpbl(klon),           zdqcv(klon)
 788       real     zsfl(klon),             zdpmel(klon,klev), &
 789               pcte(klon,klev),        zcape(klon),        &
 790               zheat(klon),            zhhatt(klon,klev),  &
 791               zhmin(klon),            zrelh(klon)
 792       real     sig1(klev)
 793       integer  ilab(klon,klev),        idtop(klon),   &
 794               ictop0(klon),           ilwmin(klon)
 795       integer  kcbot(klon),            kctop(klon),   &
 796               ktype(klon),            ihmin(klon),    &
 797               ktop0,                  lndj(klon)
 798       logical  ldcum(klon)
 799       logical  loddraf(klon),          llo1
 800 !-------------------------------------------
 801 !     1.    specify constants and parameters
 802 !-------------------------------------------
 803       zcons2=1./(g*ztmst)
 804 !--------------------------------------------------------------
 805 !*    2.    initialize values at vertical grid points in 'cuini'
 806 !--------------------------------------------------------------
 807       call cuini &
 808          (klon,     klev,     klevp1,   klevm1,   pten,  &
 809           pqen,     pqsen,    puen,     pven,     pverv, &
 810           pgeo,     paph,     zgeoh,    ztenh,    zqenh,  &
 811           zqsenh,   ilwmin,   ptu,      pqu,      ztd,   &
 812           zqd,      zuu,      zvu,      zud,      zvd,   &
 813           pmfu,     pmfd,     zmfus,    zmfds,    zmfuq, &
 814           zmfdq,    zdmfup,   zdmfdp,   zdpmel,   plu,  &
 815           plude,    ilab)
 816 !----------------------------------
 817 !*    3.0   cloud base calculations
 818 !----------------------------------
 819 !*         (a) determine cloud base values in 'cubase'
 820 !          -------------------------------------------
 821       call cubase &
 822          (klon,     klev,     klevp1,   klevm1,   ztenh, &
 823           zqenh,    zgeoh,    paph,     ptu,      pqu,   &
 824           plu,      puen,     pven,     zuu,      zvu,   &
 825           ldcum,    kcbot,    ilab)
 826 !*          (b) determine total moisture convergence and
 827 !*              then decide on type of cumulus convection
 828 !               -----------------------------------------
 829        jk=1
 830        do jl=1,klon
 831        zdqcv(jl) =pqte(jl,jk)*(paph(jl,jk+1)-paph(jl,jk))
 832        zdqpbl(jl)=0.0
 833        idtop(jl)=0
 834        end do
 835
 836        do jk=2,klev
 837        do jl=1,klon
 838        zdqcv(jl)=zdqcv(jl)+pqte(jl,jk)*(paph(jl,jk+1)-paph(jl,jk))
 839        if(jk.ge.kcbot(jl)) zdqpbl(jl)=zdqpbl(jl)+pqte(jl,jk)  &
 840                                     *(paph(jl,jk+1)-paph(jl,jk))
 841        end do
 842        end do
 843
 844       do jl=1,klon
 845          ktype(jl)=0
 846       if(zdqcv(jl).gt.max(0.,1.1*pqhfl(jl)*g)) then
 847          ktype(jl)=1
 848       else
 849          ktype(jl)=2
 850       endif
 851
 852 !*         (c) determine moisture supply for boundary layer
 853 !*             and determine cloud base massflux ignoring
 854 !*             the effects of downdrafts at this stage
 855 !              ------------------------------------------
 856       ikb=kcbot(jl)
 857       zqumqe=pqu(jl,ikb)+plu(jl,ikb)-zqenh(jl,ikb)
 858       zdqmin=max(0.01*zqenh(jl,ikb),1.e-10)
 859       if(zdqpbl(jl).gt.0..and.zqumqe.gt.zdqmin.and.ldcum(jl)) then
 860          zmfub(jl)=zdqpbl(jl)/(g*max(zqumqe,zdqmin))
 861       else
 862          zmfub(jl)=0.01
 863          ldcum(jl)=.false.
 864       endif
 865       zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
 866       zmfub(jl)=min(zmfub(jl),zmfmax)
 867 !------------------------------------------------------
 868 !*    4.0   determine cloud ascent for entraining plume
 869 !------------------------------------------------------
 870 !*         (a) estimate cloud height for entrainment/detrainment
 871 !*             calculations in cuasc (max.possible cloud height
 872 !*             for non-entraining plume, following a.-s.,1974)
 873 ! -------------------------------------------------------------
 874       ikb=kcbot(jl)
 875       zhcbase(jl)=cpd*ptu(jl,ikb)+zgeoh(jl,ikb)+alv*pqu(jl,ikb)
 876       ictop0(jl)=kcbot(jl)-1
 877       end do
 878
 879       zalvdcp=alv/cpd
 880       zqalv=1./alv
 881       do jk=klevm1,3,-1
 882       do jl=1,klon
 883       zhsat=cpd*ztenh(jl,jk)+zgeoh(jl,jk)+alv*zqsenh(jl,jk)
 884       zgam=c5les*zalvdcp*zqsenh(jl,jk)/  &
 885           ((1.-vtmpc1*zqsenh(jl,jk))*(ztenh(jl,jk)-c4les)**2)
 886       zzz=cpd*ztenh(jl,jk)*0.608
 887       zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz*zqalv)* &
 888                  max(zqsenh(jl,jk)-zqenh(jl,jk),0.)
 889       zhhatt(jl,jk)=zhhat
 890       if(jk.lt.ictop0(jl).and.zhcbase(jl).gt.zhhat) ictop0(jl)=jk
 891       end do
 892       end do
 893
 894       do jl=1,klon
 895       jk=kcbot(jl)
 896       zhsat=cpd*ztenh(jl,jk)+zgeoh(jl,jk)+alv*zqsenh(jl,jk)
 897       zgam=c5les*zalvdcp*zqsenh(jl,jk)/   &
 898           ((1.-vtmpc1*zqsenh(jl,jk))*(ztenh(jl,jk)-c4les)**2)
 899       zzz=cpd*ztenh(jl,jk)*0.608
 900       zhhat=zhsat-(zzz+zgam*zzz)/(1.+zgam*zzz*zqalv)* &
 901                  max(zqsenh(jl,jk)-zqenh(jl,jk),0.)
 902       zhhatt(jl,jk)=zhhat
 903       end do
 904 !
 905 ! find lowest possible org. detrainment level
 906 !
 907       do jl = 1, klon
 908          zhmin(jl) = 0.
 909          if( ldcum(jl).and.ktype(jl).eq.1 ) then
 910             ihmin(jl) = kcbot(jl)
 911          else
 912             ihmin(jl) = -1
 913          end if
 914       end do
 915 !
 916       zbi = 1./(25.*g)
 917       do jk = klev, 1, -1
 918       do jl = 1, klon
 919       llo1 = ldcum(jl).and.ktype(jl).eq.1.and.ihmin(jl).eq.kcbot(jl)
 920       if (llo1.and.jk.lt.kcbot(jl).and.jk.ge.ictop0(jl)) then
 921         ikb = kcbot(jl)
 922         zro = rd*ztenh(jl,jk)/(g*paph(jl,jk))
 923         zdz = (paph(jl,jk)-paph(jl,jk-1))*zro
 924         zdhdz=(cpd*(pten(jl,jk-1)-pten(jl,jk))+alv*(pqen(jl,jk-1)-   &
 925           pqen(jl,jk))+(pgeo(jl,jk-1)-pgeo(jl,jk)))*g/(pgeo(jl,      &
 926           jk-1)-pgeo(jl,jk))
 927         zdepth = zgeoh(jl,jk) - zgeoh(jl,ikb)
 928         zfac = sqrt(1.+zdepth*zbi)
 929         zhmin(jl) = zhmin(jl) + zdhdz*zfac*zdz
 930         zrh = -alv*(zqsenh(jl,jk)-zqenh(jl,jk))*zfac
 931         if (zhmin(jl).gt.zrh) ihmin(jl) = jk
 932       end if
 933       end do
 934       end do
 935
 936       do jl = 1, klon
 937       if (ldcum(jl).and.ktype(jl).eq.1) then
 938         if (ihmin(jl).lt.ictop0(jl)) ihmin(jl) = ictop0(jl)
 939       end if
 940       if(ktype(jl).eq.1) then
 941         zentr(jl)=entrpen
 942       else
 943         zentr(jl)=entrscv
 944       endif
 945       if(lndj(jl).eq.1) zentr(jl)=zentr(jl)*1.1
 946       end do
 947 !*         (b) do ascent in 'cuasc'in absence of downdrafts
 948 !----------------------------------------------------------
 949       call cuasc_new &
 950          (klon,     klev,     klevp1,   klevm1,   ztenh,   &
 951           zqenh,    puen,     pven,     pten,     pqen,    &
 952           pqsen,    pgeo,     zgeoh,    pap,      paph,    &
 953           pqte,     pverv,    ilwmin,   ldcum,    zhcbase, &
 954           ktype,    ilab,     ptu,      pqu,      plu,     &
 955           zuu,      zvu,      pmfu,     zmfub,    zentr,   &
 956           zmfus,    zmfuq,    zmful,    plude,    zdmfup,  &
 957           kcbot,    kctop,    ictop0,   icum,     ztmst,   &
 958           ihmin,    zhhatt,   zqsenh)
 959       if(icum.eq.0) return
 960 !*     (c) check cloud depth and change entrainment rate accordingly
 961 !          calculate precipitation rate (for downdraft calculation)
 962 !------------------------------------------------------------------
 963       do jl=1,klon
 964       zpbmpt=paph(jl,kcbot(jl))-paph(jl,kctop(jl))
 965       if(ldcum(jl)) ictop0(jl)=kctop(jl)
 966       if(ldcum(jl).and.ktype(jl).eq.1.and.zpbmpt.lt.zdnoprc) ktype(jl)=2
 967       if(ktype(jl).eq.2) then
 968         zentr(jl)=entrscv
 969         if(lndj(jl).eq.1) zentr(jl)=zentr(jl)*1.1
 970       endif
 971       zrfl(jl)=zdmfup(jl,1)
 972       end do
 973
 974       do jk=2,klev
 975       do jl=1,klon
 976           zrfl(jl)=zrfl(jl)+zdmfup(jl,jk)
 977       end do
 978       end do
 979 !-----------------------------------------
 980 !*    5.0   cumulus downdraft calculations
 981 !-----------------------------------------
 982       if(lmfdd) then
 983 !*      (a) determine lfs in 'cudlfs'
 984 !--------------------------------------
 985          call cudlfs &
 986          (klon,     klev,     klevp1,   ztenh,    zqenh,  &
 987           puen,     pven,     zgeoh,    paph,     ptu,    &
 988           pqu,      zuu,      zvu,      ldcum,    kcbot,  &
 989           kctop,    zmfub,    zrfl,     ztd,      zqd,    &
 990           zud,      zvd,      pmfd,     zmfds,    zmfdq,  &
 991           zdmfdp,   idtop,    loddraf)
 992 !*     (b)  determine downdraft t,q and fluxes in 'cuddraf'
 993 !------------------------------------------------------------
 994          call cuddraf &
 995          (klon,     klev,     klevp1,   ztenh,    zqenh,  &
 996           puen,     pven,     zgeoh,    paph,     zrfl,   &
 997           loddraf,  ztd,      zqd,      zud,      zvd,    &
 998           pmfd,     zmfds,    zmfdq,    zdmfdp)
 999 !*     (c)  recalculate convective fluxes due to effect of
1000 !           downdrafts on boundary layer moisture budget
1001 !-----------------------------------------------------------
1002       end if
1003 !
1004 !-- 5.1 recalculate cloud base massflux from a cape closure
1005 !       for deep convection (ktype=1) and by pbl equilibrium
1006 !       taking downdrafts into account for shallow convection
1007 !       (ktype=2)
1008 !       implemented by y. wang based on echam4 in nov. 2001.
1009 !
1010       do jl=1,klon
1011         zheat(jl)=0.0
1012         zcape(jl)=0.0
1013         zrelh(jl)=0.0
1014         zmfub1(jl)=zmfub(jl)
1015       end do
1016 !
1017       do jl=1,klon
1018       if(ldcum(jl).and.ktype(jl).eq.1) then
1019        ktop0=max(12,kctop(jl))
1020        ikb = kcbot(jl)
1021        do jk=2,klev
1022        if(jk.le.kcbot(jl).and.jk.gt.kctop(jl)) then
1023          zro=paph(jl,jk)/(rd*ztenh(jl,jk))
1024          zdz=(paph(jl,jk)-paph(jl,jk-1))/(g*zro)
1025          zheat(jl)=zheat(jl)+((pten(jl,jk-1)-pten(jl,jk)   &
1026            +g*zdz/cpd)/ztenh(jl,jk)+0.608*(pqen(jl,jk-1)-  &
1027            pqen(jl,jk)))*(pmfu(jl,jk)+pmfd(jl,jk))*g/zro
1028          zcape(jl)=zcape(jl)+g*((ptu(jl,jk)*(1.+.608*pqu(jl,jk) &
1029            -plu(jl,jk)))/(ztenh(jl,jk)*(1.+.608*zqenh(jl,jk))) &
1030            -1.0)*zdz
1031        endif
1032        if(jk.le.kcbot(jl).and.jk.gt.ktop0) then
1033          dept=(paph(jl,jk+1)-paph(jl,jk))/(paph(jl,ikb+1)-  &
1034             paph(jl,ktop0+1))
1035          zrelh(jl)=zrelh(jl)+dept*pqen(jl,jk)/pqsen(jl,jk)
1036        endif
1037        enddo
1038 !
1039        if(zrelh(jl).ge.crirh) then
1040          zht=max(0.0,(zcape(jl)-0.0))/(ztau*zheat(jl))
1041          zmfub1(jl)=max(zmfub(jl)*zht,0.01)
1042          zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
1043          zmfub1(jl)=min(zmfub1(jl),zmfmax)
1044        else
1045          zmfub1(jl)=0.01
1046          zmfub(jl)=0.01
1047          ldcum(jl)=.false.
1048         endif
1049        endif
1050        end do
1051 !
1052 !*  5.2   recalculate convective fluxes due to effect of
1053 !         downdrafts on boundary layer moisture budget
1054 !--------------------------------------------------------
1055        do jl=1,klon
1056         if(ktype(jl).ne.1) then
1057            ikb=kcbot(jl)
1058            if(pmfd(jl,ikb).lt.0.0.and.loddraf(jl)) then
1059               zeps=cmfdeps
1060            else
1061               zeps=0.
1062            endif
1063            zqumqe=pqu(jl,ikb)+plu(jl,ikb)-          &
1064                  zeps*zqd(jl,ikb)-(1.-zeps)*zqenh(jl,ikb)
1065            zdqmin=max(0.01*zqenh(jl,ikb),1.e-10)
1066            zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
1067            if(zdqpbl(jl).gt.0..and.zqumqe.gt.zdqmin.and.ldcum(jl) &
1068              .and.zmfub(jl).lt.zmfmax) then
1069               zmfub1(jl)=zdqpbl(jl)/(g*max(zqumqe,zdqmin))
1070            else
1071               zmfub1(jl)=zmfub(jl)
1072            endif
1073            llo1=(ktype(jl).eq.2).and.abs(zmfub1(jl)  &
1074                 -zmfub(jl)).lt.0.2*zmfub(jl)
1075            if(.not.llo1) zmfub1(jl)=zmfub(jl)
1076            zmfub1(jl)=min(zmfub1(jl),zmfmax)
1077         end if
1078         end do
1079
1080         do jk=1,klev
1081         do jl=1,klon
1082         if(ldcum(jl)) then
1083            zfac=zmfub1(jl)/max(zmfub(jl),1.e-10)
1084            pmfd(jl,jk)=pmfd(jl,jk)*zfac
1085            zmfds(jl,jk)=zmfds(jl,jk)*zfac
1086            zmfdq(jl,jk)=zmfdq(jl,jk)*zfac
1087            zdmfdp(jl,jk)=zdmfdp(jl,jk)*zfac
1088         else
1089            pmfd(jl,jk)=0.0
1090            zmfds(jl,jk)=0.0
1091            zmfdq(jl,jk)=0.0
1092            zdmfdp(jl,jk)=0.0
1093         endif
1094         end do
1095         end do
1096
1097         do jl=1,klon
1098            if(ldcum(jl)) then
1099               zmfub(jl)=zmfub1(jl)
1100            else
1101               zmfub(jl)=0.0
1102            endif
1103         end do
1104 !
1105 !---------------------------------------------------------------
1106 !*    6.0      determine final cloud ascent for entraining plume
1107 !*             for penetrative convection (type=1),
1108 !*             for shallow to medium convection (type=2)
1109 !*             and for mid-level convection (type=3).
1110 !---------------------------------------------------------------
1111       call cuasc_new &
1112          (klon,     klev,     klevp1,   klevm1,   ztenh,  &
1113           zqenh,    puen,     pven,     pten,     pqen,   &
1114           pqsen,    pgeo,     zgeoh,    pap,      paph,   &
1115           pqte,     pverv,    ilwmin,   ldcum,    zhcbase,&
1116           ktype,    ilab,     ptu,      pqu,      plu,    &
1117           zuu,      zvu,      pmfu,     zmfub,    zentr,  &
1118           zmfus,    zmfuq,    zmful,    plude,    zdmfup, &
1119           kcbot,    kctop,    ictop0,   icum,     ztmst,  &
1120           ihmin,    zhhatt,   zqsenh)
1121 !----------------------------------------------------------
1122 !*    7.0      determine final convective fluxes in 'cuflx'
1123 !----------------------------------------------------------
1124       call cuflx &
1125          (klon,     klev,     klevp1,   pqen,     pqsen,  &
1126           ztenh,    zqenh,    paph,     zgeoh,    kcbot,  &
1127           kctop,    idtop,    ktype,    loddraf,  ldcum,  &
1128           pmfu,     pmfd,     zmfus,    zmfds,    zmfuq,  &
1129           zmfdq,    zmful,    plude,    zdmfup,   zdmfdp, &
1130           zrfl,     prain,    pten,     zsfl,     zdpmel, &
1131           itopm2,   ztmst,    sig1)
1132 !----------------------------------------------------------------
1133 !*    8.0      update tendencies for t and q in subroutine cudtdq
1134 !----------------------------------------------------------------
1135       call cudtdq                                          &
1136          (klon,     klev,     klevp1,   itopm2,   paph,    &
1137           ldcum,    pten,     ptte,     pqte,     zmfus,   &
1138           zmfds,    zmfuq,    zmfdq,    zmful,    zdmfup,  &
1139           zdmfdp,   ztmst,    zdpmel,   prain,    zrfl,    &
1140           zsfl,     psrain,   psevap,   psheat,   psmelt,  &
1141           prsfc,    pssfc,    paprc,    paprsm,   paprs,   &
1142           pqen,     pqsen,    plude,    pcte)
1143 !----------------------------------------------------------------
1144 !*    9.0      update tendencies for u and u in subroutine cududv
1145 !----------------------------------------------------------------
1146       if(lmfdudv) then
1147       call cududv  &
1148          (klon,     klev,     klevp1,   itopm2,   ktype,   &
1149           kcbot,    paph,     ldcum,    puen,     pven,    &
1150           pvom,     pvol,     zuu,      zud,      zvu,     &
1151           zvd,      pmfu,     pmfd,     psdiss)
1152       end if
1153       return
1154       end subroutine cumastr_new
1155 !
1156
1157 !#############################################################
1158 !
1159 !             level 3 subroutines
1160 !
1161 !#############################################################
1162 !**********************************************
1163 !       subroutine cuini
1164 !**********************************************
1165 !
1166       subroutine cuini                                    &
1167          (klon,     klev,     klevp1,   klevm1,   pten,   &
1168           pqen,     pqsen,    puen,     pven,     pverv,  &
1169           pgeo,     paph,     pgeoh,    ptenh,    pqenh,  &
1170           pqsenh,   klwmin,   ptu,      pqu,      ptd,    &
1171           pqd,      puu,      pvu,      pud,      pvd,    &
1172           pmfu,     pmfd,     pmfus,    pmfds,    pmfuq,  &
1173           pmfdq,    pdmfup,   pdmfdp,   pdpmel,   plu,    &
1174           plude,    klab)
1175 !      m.tiedtke         e.c.m.w.f.     12/89
1176 !***purpose
1177 !   -------
1178 !          this routine interpolates large-scale fields of t,q etc.
1179 !          to half levels (i.e. grid for massflux scheme),
1180 !          and initializes values for updrafts and downdrafts
1181 !***interface
1182 !   ---------
1183 !          this routine is called from *cumastr*.
1184 !***method.
1185 !  --------
1186 !          for extrapolation to half levels see tiedtke(1989)
1187 !***externals
1188 !   ---------
1189 !          *cuadjtq* to specify qs at half levels
1190 ! ----------------------------------------------------------------
1191 !-------------------------------------------------------------------
1192       implicit none
1193 !-------------------------------------------------------------------
1194       integer   klon, klev, klevp1
1195       integer   klevm1
1196       integer   jk,jl,ik, icall
1197       real      zdp, zzs
1198       real     pten(klon,klev),        pqen(klon,klev),    &
1199               puen(klon,klev),        pven(klon,klev),     &
1200               pqsen(klon,klev),       pverv(klon,klev),    &
1201               pgeo(klon,klev),        pgeoh(klon,klev),    &
1202               paph(klon,klevp1),      ptenh(klon,klev),    &
1203               pqenh(klon,klev),       pqsenh(klon,klev)
1204       real     ptu(klon,klev),         pqu(klon,klev),     &
1205               ptd(klon,klev),         pqd(klon,klev),      &
1206               puu(klon,klev),         pud(klon,klev),      &
1207               pvu(klon,klev),         pvd(klon,klev),      &
1208               pmfu(klon,klev),        pmfd(klon,klev),     &
1209               pmfus(klon,klev),       pmfds(klon,klev),    &
1210               pmfuq(klon,klev),       pmfdq(klon,klev),    &
1211               pdmfup(klon,klev),      pdmfdp(klon,klev),   &
1212               plu(klon,klev),         plude(klon,klev)
1213       real     zwmax(klon),            zph(klon),          &
1214               pdpmel(klon,klev)
1215       integer  klab(klon,klev),        klwmin(klon)
1216       logical  loflag(klon)
1217 !------------------------------------------------------------
1218 !*    1.       specify large scale parameters at half levels
1219 !*             adjust temperature fields if staticly unstable
1220 !*             find level of maximum vertical velocity
1221 ! -----------------------------------------------------------
1222       zdp=0.5
1223       do jk=2,klev
1224       do jl=1,klon
1225       pgeoh(jl,jk)=pgeo(jl,jk)+(pgeo(jl,jk-1)-pgeo(jl,jk))*zdp
1226       ptenh(jl,jk)=(max(cpd*pten(jl,jk-1)+pgeo(jl,jk-1),   &
1227                   cpd*pten(jl,jk)+pgeo(jl,jk))-pgeoh(jl,jk))*rcpd
1228       pqsenh(jl,jk)=pqsen(jl,jk-1)
1229       zph(jl)=paph(jl,jk)
1230       loflag(jl)=.true.
1231       end do
1232
1233       ik=jk
1234       icall=0
1235       call cuadjtq(klon,klev,ik,zph,ptenh,pqsenh,loflag,icall)
1236       do jl=1,klon
1237       pqenh(jl,jk)=min(pqen(jl,jk-1),pqsen(jl,jk-1))    &
1238                  +(pqsenh(jl,jk)-pqsen(jl,jk-1))
1239       pqenh(jl,jk)=max(pqenh(jl,jk),0.)
1240       end do
1241       end do
1242
1243       do jl=1,klon
1244       ptenh(jl,klev)=(cpd*pten(jl,klev)+pgeo(jl,klev)-   &
1245                      pgeoh(jl,klev))*rcpd
1246       pqenh(jl,klev)=pqen(jl,klev)
1247       ptenh(jl,1)=pten(jl,1)
1248       pqenh(jl,1)=pqen(jl,1)
1249       pgeoh(jl,1)=pgeo(jl,1)
1250       klwmin(jl)=klev
1251       zwmax(jl)=0.
1252       end do
1253
1254       do jk=klevm1,2,-1
1255       do jl=1,klon
1256       zzs=max(cpd*ptenh(jl,jk)+pgeoh(jl,jk),   &
1257              cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))
1258       ptenh(jl,jk)=(zzs-pgeoh(jl,jk))*rcpd
1259       end do
1260       end do
1261
1262       do jk=klev,3,-1
1263       do jl=1,klon
1264       if(pverv(jl,jk).lt.zwmax(jl)) then
1265          zwmax(jl)=pverv(jl,jk)
1266          klwmin(jl)=jk
1267       end if
1268       end do
1269       end do
1270 !-----------------------------------------------------------
1271 !*    2.0      initialize values for updrafts and downdrafts
1272 !-----------------------------------------------------------
1273       do jk=1,klev
1274       ik=jk-1
1275       if(jk.eq.1) ik=1
1276       do jl=1,klon
1277       ptu(jl,jk)=ptenh(jl,jk)
1278       ptd(jl,jk)=ptenh(jl,jk)
1279       pqu(jl,jk)=pqenh(jl,jk)
1280       pqd(jl,jk)=pqenh(jl,jk)
1281       plu(jl,jk)=0.
1282       puu(jl,jk)=puen(jl,ik)
1283       pud(jl,jk)=puen(jl,ik)
1284       pvu(jl,jk)=pven(jl,ik)
1285       pvd(jl,jk)=pven(jl,ik)
1286       pmfu(jl,jk)=0.
1287       pmfd(jl,jk)=0.
1288       pmfus(jl,jk)=0.
1289       pmfds(jl,jk)=0.
1290       pmfuq(jl,jk)=0.
1291       pmfdq(jl,jk)=0.
1292       pdmfup(jl,jk)=0.
1293       pdmfdp(jl,jk)=0.
1294       pdpmel(jl,jk)=0.
1295       plude(jl,jk)=0.
1296       klab(jl,jk)=0
1297       end do
1298       end do
1299       return
1300       end subroutine cuini
1301
1302 !**********************************************
1303 !       subroutine cubase
1304 !**********************************************
1305       subroutine cubase &
1306          (klon,     klev,     klevp1,   klevm1,   ptenh, &
1307           pqenh,    pgeoh,    paph,     ptu,      pqu,   &
1308           plu,      puen,     pven,     puu,      pvu,   &
1309           ldcum,    kcbot,    klab)
1310 !      this routine calculates cloud base values (t and q)
1311 !      for cumulus parameterization
1312 !      m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
1313 !***purpose.
1314 !   --------
1315 !          to produce cloud base values for cu-parametrization
1316 !***interface
1317 !   ---------
1318 !          this routine is called from *cumastr*.
1319 !          input are environm. values of t,q,p,phi at half levels.
1320 !          it returns cloud base values and flags as follows;
1321 !                 klab=1 for subcloud levels
1322 !                 klab=2 for condensation level
1323 !***method.
1324 !  --------
1325 !          lift surface air dry-adiabatically to cloud base
1326 !          (non entraining plume,i.e.constant massflux)
1327 !***externals
1328 !   ---------
1329 !          *cuadjtq* for adjusting t and q due to condensation in ascent
1330 ! ----------------------------------------------------------------
1331 !-------------------------------------------------------------------
1332       implicit none
1333 !-------------------------------------------------------------------
1334       integer   klon, klev, klevp1
1335       integer   klevm1
1336       integer   jl,jk,is,ik,icall,ikb
1337       real      zbuo,zz
1338       real     ptenh(klon,klev),       pqenh(klon,klev),  &
1339               pgeoh(klon,klev),       paph(klon,klevp1)
1340       real     ptu(klon,klev),         pqu(klon,klev),   &
1341               plu(klon,klev)
1342       real     puen(klon,klev),        pven(klon,klev),  &
1343               puu(klon,klev),         pvu(klon,klev)
1344       real     zqold(klon,klev),       zph(klon)
1345       integer  klab(klon,klev),        kcbot(klon)
1346       logical  ldcum(klon),            loflag(klon)
1347       logical  ldbase(klon)
1348       logical  llo1
1349 !***input variables:
1350 !       ptenh [ztenh] - environment temperature on half levels. (cuini)
1351 !       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
1352 !       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
1353 !       paph - pressure of half levels. (mssflx)
1354 !***variables modified by cubase:
1355 !       ldcum - logical denoting profiles. (cubase)
1356 !       ktype - convection type - 1: penetrative  (cumastr)
1357 !                                 2: stratocumulus (cumastr)
1358 !                                 3: mid-level  (cuasc)
1359 !       ptu - cloud temperature.
1360 !       pqu - cloud specific humidity.
1361 !       plu - cloud liquid water (moisture condensed out)
1362 !       kcbot - cloud base level. (cubase)
1363 !       klab [ilab] - level label - 1: sub-cloud layer (cubase)
1364 !------------------------------------------------
1365 !     1.       initialize values at lifting level
1366 !------------------------------------------------
1367       do jl=1,klon
1368         klab(jl,klev)=1
1369         kcbot(jl)=klevm1
1370         ldcum(jl)=.false.
1371         ldbase(jl)=.false.
1372         puu(jl,klev)=puen(jl,klev)*(paph(jl,klevp1)-paph(jl,klev))
1373         pvu(jl,klev)=pven(jl,klev)*(paph(jl,klevp1)-paph(jl,klev))
1374       end do
1375 !-------------------------------------------------------
1376 !     2.0      do ascent in subcloud layer,
1377 !              check for existence of condensation level,
1378 !              adjust t,q and l accordingly in *cuadjtq*,
1379 !              check for buoyancy and set flags
1380 !-------------------------------------------------------
1381       do jk=1,klev
1382       do jl=1,klon
1383         zqold(jl,jk)=0.0
1384       end do
1385       end do
1386
1387       do jk=klevm1,2,-1
1388         is=0
1389         do jl=1,klon
1390           if(klab(jl,jk+1).eq.1 .or.(ldcum(jl).and.kcbot(jl).eq.jk+1)) then
1391              is=is+1
1392              loflag(jl)=.true.
1393           else
1394              loflag(jl)=.false.
1395           endif
1396           zph(jl)=paph(jl,jk)
1397         end do
1398         if(is.eq.0) cycle
1399
1400 !             calculate averages of u and v for subcloud area,
1401 !             the values will be used to define cloud base values.
1402
1403         if(lmfdudv) then
1404           do jl=1,klon
1405             if(.not.ldbase(jl)) then
1406               puu(jl,klev)=puu(jl,klev)+ &
1407                   puen(jl,jk)*(paph(jl,jk+1)-paph(jl,jk))
1408               pvu(jl,klev)=pvu(jl,klev)+ &
1409                   pven(jl,jk)*(paph(jl,jk+1)-paph(jl,jk))
1410             endif
1411           enddo
1412         endif
1413
1414         do jl=1,klon
1415           if(loflag(jl)) then
1416              pqu(jl,jk)=pqu(jl,jk+1)
1417              ptu(jl,jk)=(cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1)  &
1418                        -pgeoh(jl,jk))*rcpd
1419              zqold(jl,jk)=pqu(jl,jk)
1420           end if
1421         end do
1422
1423         ik=jk
1424         icall=1
1425         call cuadjtq(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1426
1427         do jl=1,klon
1428           if(loflag(jl)) then
1429            if(pqu(jl,jk).eq.zqold(jl,jk)) then
1430              zbuo=ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk))-      &
1431                  ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))+zbuo0
1432              if(zbuo.gt.0.) klab(jl,jk)=1
1433            else
1434              klab(jl,jk)=2
1435              plu(jl,jk)=plu(jl,jk)+zqold(jl,jk)-pqu(jl,jk)
1436              zbuo=ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk))-      &
1437                  ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))+zbuo0
1438              llo1=zbuo.gt.0..and.klab(jl,jk+1).eq.1
1439              if(llo1) then
1440                 kcbot(jl)=jk
1441                 ldcum(jl)=.true.
1442                 ldbase(jl)=.true.
1443              end if
1444            end if
1445           end if
1446         end do
1447       end do
1448
1449       if(lmfdudv) then
1450          do jl=1,klon
1451          if(ldcum(jl)) then
1452             ikb=kcbot(jl)
1453             zz=1./(paph(jl,klevp1)-paph(jl,ikb))
1454             puu(jl,klev)=puu(jl,klev)*zz
1455             pvu(jl,klev)=pvu(jl,klev)*zz
1456          else
1457             puu(jl,klev)=puen(jl,klevm1)
1458             pvu(jl,klev)=pven(jl,klevm1)
1459          end if
1460          end do
1461       end if
1462       return
1463       end subroutine cubase
1464 !
1465 !**********************************************
1466 !       subroutine cuasc_new
1467 !**********************************************
1468       subroutine cuasc_new &
1469          (klon,     klev,     klevp1,   klevm1,   ptenh,  &
1470           pqenh,    puen,     pven,     pten,     pqen,   &
1471           pqsen,    pgeo,     pgeoh,    pap,      paph,   &
1472           pqte,     pverv,    klwmin,   ldcum,    phcbase,&
1473           ktype,    klab,     ptu,      pqu,      plu,    &
1474           puu,      pvu,      pmfu,     pmfub,    pentr,  &
1475           pmfus,    pmfuq,    pmful,    plude,    pdmfup, &
1476           kcbot,    kctop,    kctop0,   kcum,     ztmst,  &
1477           khmin,    phhatt,   pqsenh)
1478 !     this routine does the calculations for cloud ascents
1479 !     for cumulus parameterization
1480 !     m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
1481 !     y.wang            iprc           11/01 modif.
1482 !***purpose.
1483 !   --------
1484 !          to produce cloud ascents for cu-parametrization
1485 !          (vertical profiles of t,q,l,u and v and corresponding
1486 !           fluxes as well as precipitation rates)
1487 !***interface
1488 !   ---------
1489 !          this routine is called from *cumastr*.
1490 !***method.
1491 !  --------
1492 !          lift surface air dry-adiabatically to cloud base
1493 !          and then calculate moist ascent for
1494 !          entraining/detraining plume.
1495 !          entrainment and detrainment rates differ for
1496 !          shallow and deep cumulus convection.
1497 !          in case there is no penetrative or shallow convection
1498 !          check for possibility of mid level convection
1499 !          (cloud base values calculated in *cubasmc*)
1500 !***externals
1501 !   ---------
1502 !          *cuadjtq* adjust t and q due to condensation in ascent
1503 !          *cuentr_new*  calculate entrainment/detrainment rates
1504 !          *cubasmc* calculate cloud base values for midlevel convection
1505 !***reference
1506 !   ---------
1507 !          (tiedtke,1989)
1508 !***input variables:
1509 !       ptenh [ztenh] - environ temperature on half levels. (cuini)
1510 !       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
1511 !       puen - environment wind u-component. (mssflx)
1512 !       pven - environment wind v-component. (mssflx)
1513 !       pten - environment temperature. (mssflx)
1514 !       pqen - environment specific humidity. (mssflx)
1515 !       pqsen - environment saturation specific humidity. (mssflx)
1516 !       pgeo - geopotential. (mssflx)
1517 !       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
1518 !       pap - pressure in pa.  (mssflx)
1519 !       paph - pressure of half levels. (mssflx)
1520 !       pqte - moisture convergence (delta q/delta t). (mssflx)
1521 !       pverv - large scale vertical velocity (omega). (mssflx)
1522 !       klwmin [ilwmin] - level of minimum omega. (cuini)
1523 !       klab [ilab] - level label - 1: sub-cloud layer.
1524 !                                   2: condensation level (cloud base)
1525 !       pmfub [zmfub] - updraft mass flux at cloud base. (cumastr)
1526 !***variables modified by cuasc:
1527 !       ldcum - logical denoting profiles. (cubase)
1528 !       ktype - convection type - 1: penetrative  (cumastr)
1529 !                                 2: stratocumulus (cumastr)
1530 !                                 3: mid-level  (cuasc)
1531 !       ptu - cloud temperature.
1532 !       pqu - cloud specific humidity.
1533 !       plu - cloud liquid water (moisture condensed out)
1534 !       puu [zuu] - cloud momentum u-component.
1535 !       pvu [zvu] - cloud momentum v-component.
1536 !       pmfu - updraft mass flux.
1537 !       pentr [zentr] - entrainment rate. (cumastr ) (cubasmc)
1538 !       pmfus [zmfus] - updraft flux of dry static energy. (cubasmc)
1539 !       pmfuq [zmfuq] - updraft flux of specific humidity.
1540 !       pmful [zmful] - updraft flux of cloud liquid water.
1541 !       plude - liquid water returned to environment by detrainment.
1542 !       pdmfup [zmfup] -
1543 !       kcbot - cloud base level. (cubase)
1544 !       kctop -
1545 !       kctop0 [ictop0] - estimate of cloud top. (cumastr)
1546 !       kcum [icum] -
1547 !-------------------------------------------------------------------
1548       implicit none
1549 !-------------------------------------------------------------------
1550       integer   klon, klev, klevp1
1551       integer   klevm1,kcum
1552       real      ztmst,zcons2,zdz,zdrodz
1553       integer   jl,jk,ikb,ik,is,ikt,icall
1554       real      zmfmax,zfac,zmftest,zdprho,zmse,znevn,zodmax
1555       real      zqeen,zseen,zscde,zga,zdt,zscod
1556       real      zqude,zqcod, zmfusk, zmfuqk,zmfulk
1557       real      zbuo, zprcon, zlnew, zz, zdmfeu, zdmfdu
1558       real      zbuoyz,zzdmf
1559       real     ptenh(klon,klev),       pqenh(klon,klev), &
1560               puen(klon,klev),        pven(klon,klev),   &
1561               pten(klon,klev),        pqen(klon,klev),   &
1562               pgeo(klon,klev),        pgeoh(klon,klev),  &
1563               pap(klon,klev),         paph(klon,klevp1), &
1564               pqsen(klon,klev),       pqte(klon,klev),   &
1565               pverv(klon,klev),       pqsenh(klon,klev)
1566       real     ptu(klon,klev),         pqu(klon,klev),   &
1567               puu(klon,klev),         pvu(klon,klev),    &
1568               pmfu(klon,klev),        zph(klon),         &
1569               pmfub(klon),            pentr(klon),       &
1570               pmfus(klon,klev),       pmfuq(klon,klev),  &
1571               plu(klon,klev),         plude(klon,klev),  &
1572               pmful(klon,klev),       pdmfup(klon,klev)
1573       real     zdmfen(klon),           zdmfde(klon),     &
1574               zmfuu(klon),            zmfuv(klon),       &
1575               zpbase(klon),           zqold(klon),       &
1576               phhatt(klon,klev),      zodetr(klon,klev), &
1577               zoentr(klon,klev),      zbuoy(klon)
1578       real     phcbase(klon)
1579       integer  klwmin(klon),           ktype(klon),      &
1580               klab(klon,klev),        kcbot(klon),       &
1581               kctop(klon),            kctop0(klon),      &
1582               khmin(klon)
1583       logical  ldcum(klon),            loflag(klon)
1584 !--------------------------------
1585 !*    1.       specify parameters
1586 !--------------------------------
1587       zcons2=1./(g*ztmst)
1588 !---------------------------------
1589 !     2.        set default values
1590 !---------------------------------
1591       do jl=1,klon
1592         zmfuu(jl)=0.
1593         zmfuv(jl)=0.
1594         zbuoy(jl)=0.
1595         if(.not.ldcum(jl)) ktype(jl)=0
1596       end do
1597
1598       do jk=1,klev
1599       do jl=1,klon
1600           plu(jl,jk)=0.
1601           pmfu(jl,jk)=0.
1602           pmfus(jl,jk)=0.
1603           pmfuq(jl,jk)=0.
1604           pmful(jl,jk)=0.
1605           plude(jl,jk)=0.
1606           pdmfup(jl,jk)=0.
1607           zoentr(jl,jk)=0.
1608           zodetr(jl,jk)=0.
1609           if(.not.ldcum(jl).or.ktype(jl).eq.3) klab(jl,jk)=0
1610           if(.not.ldcum(jl).and.paph(jl,jk).lt.4.e4) kctop0(jl)=jk
1611       end do
1612       end do
1613 !------------------------------------------------
1614 !     3.0      initialize values at lifting level
1615 !------------------------------------------------
1616       do jl=1,klon
1617         kctop(jl)=klevm1
1618         if(.not.ldcum(jl)) then
1619            kcbot(jl)=klevm1
1620            pmfub(jl)=0.
1621            pqu(jl,klev)=0.
1622         end if
1623         pmfu(jl,klev)=pmfub(jl)
1624         pmfus(jl,klev)=pmfub(jl)*(cpd*ptu(jl,klev)+pgeoh(jl,klev))
1625         pmfuq(jl,klev)=pmfub(jl)*pqu(jl,klev)
1626         if(lmfdudv) then
1627            zmfuu(jl)=pmfub(jl)*puu(jl,klev)
1628            zmfuv(jl)=pmfub(jl)*pvu(jl,klev)
1629         end if
1630       end do
1631 !
1632 !-- 3.1 find organized entrainment at cloud base
1633 !
1634       do jl=1,klon
1635       ldcum(jl)=.false.
1636       if (ktype(jl).eq.1) then
1637       ikb = kcbot(jl)
1638       zbuoy(jl)=g*((ptu(jl,ikb)-ptenh(jl,ikb))/ptenh(jl,ikb)+ &
1639                0.608*(pqu(jl,ikb)-pqenh(jl,ikb)))
1640        if (zbuoy(jl).gt.0.) then
1641         zdz = (pgeo(jl,ikb-1)-pgeo(jl,ikb))*zrg
1642         zdrodz = -log(pten(jl,ikb-1)/pten(jl,ikb))/zdz -  &
1643                  g/(rd*ptenh(jl,ikb))
1644         zoentr(jl,ikb-1)=zbuoy(jl)*0.5/(1.+zbuoy(jl)*zdz) &
1645                 +zdrodz
1646         zoentr(jl,ikb-1) = min(zoentr(jl,ikb-1),1.e-3)
1647         zoentr(jl,ikb-1) = max(zoentr(jl,ikb-1),0.)
1648        end if
1649       end if
1650       end do
1651 !
1652 !-----------------------------------------------------------------
1653 !     4.       do ascent: subcloud layer (klab=1) ,clouds (klab=2)
1654 !              by doing first dry-adiabatic ascent and then
1655 !              by adjusting t,q and l accordingly in *cuadjtq*,
1656 !              then check for buoyancy and set flags accordingly
1657 !-----------------------------------------------------------------
1658       do jk=klevm1,2,-1
1659 !                  specify cloud base values for midlevel convection
1660 !                  in *cubasmc* in case there is not already convection
1661 ! ---------------------------------------------------------------------
1662       ik=jk
1663       if(lmfmid.and.ik.lt.klevm1.and.ik.gt.klev-13) then
1664       call cubasmc  &
1665          (klon,     klev,     klevm1,   ik,      pten,  &
1666           pqen,     pqsen,    puen,     pven,    pverv, &
1667           pgeo,     pgeoh,    ldcum,    ktype,   klab,  &
1668           pmfu,     pmfub,    pentr,    kcbot,   ptu,   &
1669           pqu,      plu,      puu,     pvu,      pmfus, &
1670           pmfuq,    pmful,    pdmfup,  zmfuu,    zmfuv)
1671       endif
1672       is=0
1673       do jl=1,klon
1674         zqold(jl)=0.0
1675         is=is+klab(jl,jk+1)
1676         if(klab(jl,jk+1).eq.0) klab(jl,jk)=0
1677         loflag(jl)=klab(jl,jk+1).gt.0
1678         zph(jl)=paph(jl,jk)
1679         if(ktype(jl).eq.3.and.jk.eq.kcbot(jl)) then
1680            zmfmax=(paph(jl,jk)-paph(jl,jk-1))*zcons2
1681            if(pmfub(jl).gt.zmfmax) then
1682               zfac=zmfmax/pmfub(jl)
1683               pmfu(jl,jk+1)=pmfu(jl,jk+1)*zfac
1684               pmfus(jl,jk+1)=pmfus(jl,jk+1)*zfac
1685               pmfuq(jl,jk+1)=pmfuq(jl,jk+1)*zfac
1686               zmfuu(jl)=zmfuu(jl)*zfac
1687               zmfuv(jl)=zmfuv(jl)*zfac
1688               pmfub(jl)=zmfmax
1689            end if
1690         end if
1691       end do
1692
1693       if(is.eq.0) cycle
1694 !
1695 !*     specify entrainment rates in *cuentr_new*
1696 ! -------------------------------------
1697       ik=jk
1698       call cuentr_new &
1699          (klon,     klev,     klevp1,   ik,       ptenh,&
1700           paph,     pap,      pgeoh,    klwmin,   ldcum,&
1701           ktype,    kcbot,    kctop0,   zpbase,   pmfu, &
1702           pentr,    zdmfen,   zdmfde,   zodetr,   khmin)
1703 !
1704 !      do adiabatic ascent for entraining/detraining plume
1705 ! -------------------------------------------------------
1706 ! do adiabatic ascent for entraining/detraining plume
1707 ! the cloud ensemble entrains environmental values
1708 ! in turbulent detrainment cloud ensemble values are detrained
1709 ! in organized detrainment the dry static energy and
1710 ! moisture that are neutral compared to the
1711 ! environmental air are detrained
1712 !
1713       do jl=1,klon
1714       if(loflag(jl)) then
1715         if(jk.lt.kcbot(jl)) then
1716          zmftest=pmfu(jl,jk+1)+zdmfen(jl)-zdmfde(jl)
1717          zmfmax=min(zmftest,(paph(jl,jk)-paph(jl,jk-1))*zcons2)
1718          zdmfen(jl)=max(zdmfen(jl)-max(zmftest-zmfmax,0.),0.)
1719         end if
1720         zdmfde(jl)=min(zdmfde(jl),0.75*pmfu(jl,jk+1))
1721         pmfu(jl,jk)=pmfu(jl,jk+1)+zdmfen(jl)-zdmfde(jl)
1722         if (jk.lt.kcbot(jl)) then
1723           zdprho = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1724           zoentr(jl,jk) = zoentr(jl,jk)*zdprho*pmfu(jl,jk+1)
1725           zmftest = pmfu(jl,jk) + zoentr(jl,jk)-zodetr(jl,jk)
1726           zmfmax = min(zmftest,(paph(jl,jk)-paph(jl,jk-1))*zcons2)
1727           zoentr(jl,jk) = max(zoentr(jl,jk)-max(zmftest-zmfmax,0.),0.)
1728         end if
1729 !
1730 ! limit organized detrainment to not allowing for too deep clouds
1731 !
1732         if (ktype(jl).eq.1.and.jk.lt.kcbot(jl).and.jk.le.khmin(jl)) then
1733           zmse = cpd*ptu(jl,jk+1) + alv*pqu(jl,jk+1) + pgeoh(jl,jk+1)
1734           ikt = kctop0(jl)
1735           znevn=(pgeoh(jl,ikt)-pgeoh(jl,jk+1))*(zmse-phhatt(jl,  &
1736                jk+1))*zrg
1737           if (znevn.le.0.) znevn = 1.
1738           zdprho = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1739           zodmax = ((phcbase(jl)-zmse)/znevn)*zdprho*pmfu(jl,jk+1)
1740           zodmax = max(zodmax,0.)
1741           zodetr(jl,jk) = min(zodetr(jl,jk),zodmax)
1742         end if
1743         zodetr(jl,jk) = min(zodetr(jl,jk),0.75*pmfu(jl,jk))
1744         pmfu(jl,jk) = pmfu(jl,jk) + zoentr(jl,jk) - zodetr(jl,jk)
1745         zqeen=pqenh(jl,jk+1)*zdmfen(jl)
1746         zqeen=zqeen + pqenh(jl,jk+1)*zoentr(jl,jk)
1747         zseen=(cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*zdmfen(jl)
1748         zseen=zseen+(cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*  &
1749              zoentr(jl,jk)
1750         zscde=(cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1))*zdmfde(jl)
1751 ! find moist static energy that give nonbuoyant air
1752         zga = alv*pqsenh(jl,jk+1)/(rv*(ptenh(jl,jk+1)**2))
1753         zdt = (plu(jl,jk+1)-0.608*(pqsenh(jl,jk+1)-pqenh(jl, &
1754                jk+1)))/(1./ptenh(jl,jk+1)+0.608*zga)
1755         zscod = cpd*ptenh(jl,jk+1) + pgeoh(jl,jk+1) + cpd*zdt
1756         zscde = zscde + zodetr(jl,jk)*zscod
1757         zqude = pqu(jl,jk+1)*zdmfde(jl)
1758         zqcod = pqsenh(jl,jk+1) + zga*zdt
1759         zqude = zqude + zodetr(jl,jk)*zqcod
1760         plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
1761         plude(jl,jk) = plude(jl,jk)+plu(jl,jk+1)*zodetr(jl,jk)
1762         zmfusk = pmfus(jl,jk+1) + zseen - zscde
1763         zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude
1764         zmfulk = pmful(jl,jk+1) - plude(jl,jk)
1765         plu(jl,jk) = zmfulk*(1./max(cmfcmin,pmfu(jl,jk)))
1766         pqu(jl,jk) = zmfuqk*(1./max(cmfcmin,pmfu(jl,jk)))
1767         ptu(jl,jk)=(zmfusk*(1./max(cmfcmin,pmfu(jl,jk)))-  &
1768             pgeoh(jl,jk))*rcpd
1769         ptu(jl,jk) = max(100.,ptu(jl,jk))
1770         ptu(jl,jk) = min(400.,ptu(jl,jk))
1771         zqold(jl) = pqu(jl,jk)
1772       end if
1773       end do
1774 !*             do corrections for moist ascent
1775 !*             by adjusting t,q and l in *cuadjtq*
1776 !------------------------------------------------
1777       ik=jk
1778       icall=1
1779 !
1780       call cuadjtq(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1781 !
1782       do jl=1,klon
1783       if(loflag(jl).and.pqu(jl,jk).ne.zqold(jl)) then
1784          klab(jl,jk)=2
1785          plu(jl,jk)=plu(jl,jk)+zqold(jl)-pqu(jl,jk)
1786          zbuo=ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))-  &
1787         ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
1788          if(klab(jl,jk+1).eq.1) zbuo=zbuo+zbuo0
1789          if(zbuo.gt.0..and.pmfu(jl,jk).gt.0.01*pmfub(jl).and. &
1790                             jk.ge.kctop0(jl)) then
1791             kctop(jl)=jk
1792             ldcum(jl)=.true.
1793             if(zpbase(jl)-paph(jl,jk).ge.zdnoprc) then
1794                zprcon=cprcon
1795             else
1796                zprcon=0.
1797             endif
1798             zlnew=plu(jl,jk)/(1.+zprcon*(pgeoh(jl,jk)-pgeoh(jl,jk+1)))
1799             pdmfup(jl,jk)=max(0.,(plu(jl,jk)-zlnew)*pmfu(jl,jk))
1800             plu(jl,jk)=zlnew
1801          else
1802             klab(jl,jk)=0
1803             pmfu(jl,jk)=0.
1804          end if
1805       end if
1806       if(loflag(jl)) then
1807          pmful(jl,jk)=plu(jl,jk)*pmfu(jl,jk)
1808          pmfus(jl,jk)=(cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk)
1809          pmfuq(jl,jk)=pqu(jl,jk)*pmfu(jl,jk)
1810       end if
1811      end do
1812 !
1813       if(lmfdudv) then
1814 !
1815         do jl=1,klon
1816         zdmfen(jl) = zdmfen(jl) + zoentr(jl,jk)
1817         zdmfde(jl) = zdmfde(jl) + zodetr(jl,jk)
1818            if(loflag(jl)) then
1819               if(ktype(jl).eq.1.or.ktype(jl).eq.3) then
1820                  if(zdmfen(jl).le.1.e-20) then
1821                     zz=3.
1822                  else
1823                     zz=2.
1824                  endif
1825               else
1826                  if(zdmfen(jl).le.1.0e-20) then
1827                     zz=1.
1828                  else
1829                     zz=0.
1830                  endif
1831               end if
1832               zdmfeu=zdmfen(jl)+zz*zdmfde(jl)
1833               zdmfdu=zdmfde(jl)+zz*zdmfde(jl)
1834               zdmfdu=min(zdmfdu,0.75*pmfu(jl,jk+1))
1835               zmfuu(jl)=zmfuu(jl)+                              &
1836                        zdmfeu*puen(jl,jk)-zdmfdu*puu(jl,jk+1)
1837               zmfuv(jl)=zmfuv(jl)+                              &
1838                        zdmfeu*pven(jl,jk)-zdmfdu*pvu(jl,jk+1)
1839               if(pmfu(jl,jk).gt.0.) then
1840                  puu(jl,jk)=zmfuu(jl)*(1./pmfu(jl,jk))
1841                  pvu(jl,jk)=zmfuv(jl)*(1./pmfu(jl,jk))
1842               end if
1843            end if
1844         end do
1845 !
1846         end if
1847 !
1848 ! compute organized entrainment
1849 ! for use at next level
1850 !
1851       do jl = 1, klon
1852        if (loflag(jl).and.ktype(jl).eq.1) then
1853         zbuoyz=g*((ptu(jl,jk)-ptenh(jl,jk))/ptenh(jl,jk)+  &
1854               0.608*(pqu(jl,jk)-pqenh(jl,jk))-plu(jl,jk))
1855         zbuoyz = max(zbuoyz,0.0)
1856         zdz = (pgeo(jl,jk-1)-pgeo(jl,jk))*zrg
1857         zdrodz = -log(pten(jl,jk-1)/pten(jl,jk))/zdz -  &
1858                  g/(rd*ptenh(jl,jk))
1859         zbuoy(jl) = zbuoy(jl) + zbuoyz*zdz
1860         zoentr(jl,jk-1) = zbuoyz*0.5/(1.+zbuoy(jl))+zdrodz
1861         zoentr(jl,jk-1) = min(zoentr(jl,jk-1),1.e-3)
1862         zoentr(jl,jk-1) = max(zoentr(jl,jk-1),0.)
1863        end if
1864       end do
1865 !
1866     end do ! end outer cycle
1867 ! -----------------------------------------------------------------
1868 !     5.       determine convective fluxes above non-buoyancy level
1869 ! -----------------------------------------------------------------
1870 !                  (note: cloud variables like t,q and l are not
1871 !                         affected by detrainment and are already known
1872 !                         from previous calculations above)
1873       do jl=1,klon
1874       if(kctop(jl).eq.klevm1) ldcum(jl)=.false.
1875       kcbot(jl)=max(kcbot(jl),kctop(jl))
1876       end do
1877
1878       is=0
1879       do jl=1,klon
1880       if(ldcum(jl)) then
1881          is=is+1
1882       endif
1883       end do
1884       kcum=is
1885       if(is.eq.0) return
1886       do jl=1,klon
1887       if(ldcum(jl)) then
1888          jk=kctop(jl)-1
1889          zzdmf=cmfctop
1890          zdmfde(jl)=(1.-zzdmf)*pmfu(jl,jk+1)
1891          plude(jl,jk)=zdmfde(jl)*plu(jl,jk+1)
1892          pmfu(jl,jk)=pmfu(jl,jk+1)-zdmfde(jl)
1893          pmfus(jl,jk)=(cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk)
1894          pmfuq(jl,jk)=pqu(jl,jk)*pmfu(jl,jk)
1895          pmful(jl,jk)=plu(jl,jk)*pmfu(jl,jk)
1896          plude(jl,jk-1)=pmful(jl,jk)
1897          pdmfup(jl,jk)=0.
1898       end if
1899       end do
1900
1901       if(lmfdudv) then
1902          do jl=1,klon
1903          if(ldcum(jl)) then
1904             jk=kctop(jl)-1
1905             puu(jl,jk)=puu(jl,jk+1)
1906             pvu(jl,jk)=pvu(jl,jk+1)
1907          end if
1908          end do
1909       end if
1910       return
1911       end subroutine cuasc_new
1912 !
1913
1914 !**********************************************
1915 !       subroutine cudlfs
1916 !**********************************************
1917       subroutine cudlfs &
1918          (klon,     klev,     klevp1,   ptenh,    pqenh,  &
1919           puen,     pven,     pgeoh,    paph,     ptu,    &
1920           pqu,      puu,      pvu,      ldcum,    kcbot,  &
1921           kctop,    pmfub,    prfl,     ptd,      pqd,    &
1922           pud,      pvd,      pmfd,     pmfds,    pmfdq,  &
1923           pdmfdp,   kdtop,    lddraf)
1924 !      this routine calculates level of free sinking for
1925 !      cumulus downdrafts and specifies t,q,u and v values
1926 !      m.tiedtke         e.c.m.w.f.    12/86 modif.  12/89
1927 !***purpose.
1928 !   --------
1929 !          to produce lfs-values for cumulus downdrafts
1930 !          for massflux cumulus parameterization
1931 !***interface
1932 !   ---------
1933 !          this routine is called from *cumastr*.
1934 !          input are environmental values of t,q,u,v,p,phi
1935 !          and updraft values t,q,u and v and also
1936 !          cloud base massflux and cu-precipitation rate.
1937 !          it returns t,q,u and v values and massflux at lfs.
1938 !***method.
1939 !  --------
1940 !          check for negative buoyancy of air of equal parts of
1941 !          moist environmental air and cloud air.
1942 !***externals
1943 !   ---------
1944 !          *cuadjtq* for calculating wet bulb t and q at lfs
1945 ! ----------------------------------------------------------------
1946 !-------------------------------------------------------------------
1947       implicit none
1948 !-------------------------------------------------------------------
1949       integer   klon, klev, klevp1
1950       integer   jl,ke,jk,is,ik,icall
1951       real      zttest, zqtest, zbuo, zmftop
1952       real     ptenh(klon,klev),       pqenh(klon,klev),   &
1953               puen(klon,klev),        pven(klon,klev),     &
1954               pgeoh(klon,klev),       paph(klon,klevp1),   &
1955               ptu(klon,klev),         pqu(klon,klev),      &
1956               puu(klon,klev),         pvu(klon,klev),      &
1957               pmfub(klon),            prfl(klon)
1958       real     ptd(klon,klev),         pqd(klon,klev),     &
1959               pud(klon,klev),         pvd(klon,klev),      &
1960               pmfd(klon,klev),        pmfds(klon,klev),    &
1961               pmfdq(klon,klev),       pdmfdp(klon,klev)
1962       real     ztenwb(klon,klev),      zqenwb(klon,klev),  &
1963               zcond(klon),            zph(klon)
1964       integer  kcbot(klon),            kctop(klon),        &
1965               kdtop(klon)
1966       logical  ldcum(klon),            llo2(klon),         &
1967               lddraf(klon)
1968 !-----------------------------------------------
1969 !     1.       set default values for downdrafts
1970 !-----------------------------------------------
1971       do jl=1,klon
1972       lddraf(jl)=.false.
1973       kdtop(jl)=klevp1
1974       end do
1975       if(.not.lmfdd) return
1976 !------------------------------------------------------------
1977 !     2.       determine level of free sinking by
1978 !              doing a scan from top to base of cumulus clouds
1979 !              for every point and proceed as follows:
1980 !                (1) detemine wet bulb environmental t and q
1981 !                (2) do mixing with cumulus cloud air
1982 !                (3) check for negative buoyancy
1983 !              the assumption is that air of downdrafts is mixture
1984 !              of 50% cloud air + 50% environmental air at wet bulb
1985 !              temperature (i.e. which became saturated due to
1986 !              evaporation of rain and cloud water)
1987 !------------------------------------------------------------------
1988       ke=klev-3
1989       do jk=3,ke
1990 !   2.1      calculate wet-bulb temperature and moisture
1991 !            for environmental air in *cuadjtq*
1992 ! -----------------------------------------------------
1993       is=0
1994       do jl=1,klon
1995       ztenwb(jl,jk)=ptenh(jl,jk)
1996       zqenwb(jl,jk)=pqenh(jl,jk)
1997       zph(jl)=paph(jl,jk)
1998       llo2(jl)=ldcum(jl).and.prfl(jl).gt.0..and..not.lddraf(jl).and. &
1999               (jk.lt.kcbot(jl).and.jk.gt.kctop(jl))
2000       if(llo2(jl))then
2001          is=is+1
2002       endif
2003       end do
2004
2005       if(is.eq.0) cycle
2006       ik=jk
2007       icall=2
2008       call cuadjtq(klon,klev,ik,zph,ztenwb,zqenwb,llo2,icall)
2009 !   2.2      do mixing of cumulus and environmental air
2010 !            and check for negative buoyancy.
2011 !            then set values for downdraft at lfs.
2012 ! -----------------------------------------------------
2013       do jl=1,klon
2014       if(llo2(jl)) then
2015          zttest=0.5*(ptu(jl,jk)+ztenwb(jl,jk))
2016          zqtest=0.5*(pqu(jl,jk)+zqenwb(jl,jk))
2017          zbuo=zttest*(1.+vtmpc1*zqtest)-  &
2018              ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2019          zcond(jl)=pqenh(jl,jk)-zqenwb(jl,jk)
2020          zmftop=-cmfdeps*pmfub(jl)
2021          if(zbuo.lt.0..and.prfl(jl).gt.10.*zmftop*zcond(jl)) then
2022             kdtop(jl)=jk
2023             lddraf(jl)=.true.
2024             ptd(jl,jk)=zttest
2025             pqd(jl,jk)=zqtest
2026             pmfd(jl,jk)=zmftop
2027             pmfds(jl,jk)=pmfd(jl,jk)*(cpd*ptd(jl,jk)+pgeoh(jl,jk))
2028             pmfdq(jl,jk)=pmfd(jl,jk)*pqd(jl,jk)
2029             pdmfdp(jl,jk-1)=-0.5*pmfd(jl,jk)*zcond(jl)
2030             prfl(jl)=prfl(jl)+pdmfdp(jl,jk-1)
2031          end if
2032       end if
2033       end do
2034
2035       if(lmfdudv) then
2036           do jl=1,klon
2037           if(pmfd(jl,jk).lt.0.) then
2038              pud(jl,jk)=0.5*(puu(jl,jk)+puen(jl,jk-1))
2039              pvd(jl,jk)=0.5*(pvu(jl,jk)+pven(jl,jk-1))
2040           end if
2041           end do
2042       end if
2043
2044       end do
2045       return
2046       end subroutine cudlfs
2047 !
2048
2049 !**********************************************
2050 !       subroutine cuddraf
2051 !**********************************************
2052       subroutine cuddraf &
2053          (klon,     klev,     klevp1,   ptenh,    pqenh, &
2054           puen,     pven,     pgeoh,    paph,     prfl,  &
2055           lddraf,   ptd,      pqd,      pud,      pvd,   &
2056           pmfd,     pmfds,    pmfdq,    pdmfdp)
2057 !     this routine calculates cumulus downdraft descent
2058 !     m.tiedtke         e.c.m.w.f.    12/86 modif.  12/89
2059 !***purpose.
2060 !   --------
2061 !          to produce the vertical profiles for cumulus downdrafts
2062 !          (i.e. t,q,u and v and fluxes)
2063 !***interface
2064 !   ---------
2065 !          this routine is called from *cumastr*.
2066 !          input is t,q,p,phi,u,v at half levels.
2067 !          it returns fluxes of s,q and evaporation rate
2068 !          and u,v at levels where downdraft occurs
2069 !***method.
2070 !  --------
2071 !          calculate moist descent for entraining/detraining plume by
2072 !          a) moving air dry-adiabatically to next level below and
2073 !          b) correcting for evaporation to obtain saturated state.
2074 !***externals
2075 !   ---------
2076 !          *cuadjtq* for adjusting t and q due to evaporation in
2077 !          saturated descent
2078 !***reference
2079 !   ---------
2080 !          (tiedtke,1989)
2081 ! ----------------------------------------------------------------
2082 !-------------------------------------------------------------------
2083       implicit none
2084 !-------------------------------------------------------------------
2085       integer   klon, klev, klevp1
2086       integer   jk,is,jl,itopde, ik, icall
2087       real      zentr,zseen, zqeen, zsdde, zqdde,zmfdsk, zmfdqk
2088       real      zbuo, zdmfdp, zmfduk, zmfdvk
2089       real     ptenh(klon,klev),       pqenh(klon,klev),  &
2090               puen(klon,klev),        pven(klon,klev),    &
2091               pgeoh(klon,klev),       paph(klon,klevp1)
2092       real     ptd(klon,klev),         pqd(klon,klev),    &
2093               pud(klon,klev),         pvd(klon,klev),     &
2094               pmfd(klon,klev),        pmfds(klon,klev),   &
2095               pmfdq(klon,klev),       pdmfdp(klon,klev),  &
2096               prfl(klon)
2097       real     zdmfen(klon),           zdmfde(klon),      &
2098               zcond(klon),            zph(klon)
2099       logical  lddraf(klon),           llo2(klon)
2100 !--------------------------------------------------------------
2101 !     1.       calculate moist descent for cumulus downdraft by
2102 !                (a) calculating entrainment rates, assuming
2103 !                     linear decrease of massflux in pbl
2104 !                 (b) doing moist descent - evaporative cooling
2105 !                     and moistening is calculated in *cuadjtq*
2106 !                 (c) checking for negative buoyancy and
2107 !                     specifying final t,q,u,v and downward fluxes
2108 ! ----------------------------------------------------------------
2109       do jk=3,klev
2110       is=0
2111       do jl=1,klon
2112       zph(jl)=paph(jl,jk)
2113       llo2(jl)=lddraf(jl).and.pmfd(jl,jk-1).lt.0.
2114       if(llo2(jl)) then
2115          is=is+1
2116       endif
2117       end do
2118
2119       if(is.eq.0) cycle
2120       do jl=1,klon
2121       if(llo2(jl)) then
2122          zentr=entrdd*pmfd(jl,jk-1)*rd*ptenh(jl,jk-1)/   &
2123               (g*paph(jl,jk-1))*(paph(jl,jk)-paph(jl,jk-1))
2124          zdmfen(jl)=zentr
2125          zdmfde(jl)=zentr
2126       end if
2127       end do
2128
2129       itopde=klev-2
2130       if(jk.gt.itopde) then
2131          do jl=1,klon
2132          if(llo2(jl)) then
2133             zdmfen(jl)=0.
2134             zdmfde(jl)=pmfd(jl,itopde)*      &
2135             (paph(jl,jk)-paph(jl,jk-1))/     &
2136             (paph(jl,klevp1)-paph(jl,itopde))
2137          end if
2138          end do
2139       end if
2140
2141       do jl=1,klon
2142          if(llo2(jl)) then
2143             pmfd(jl,jk)=pmfd(jl,jk-1)+zdmfen(jl)-zdmfde(jl)
2144             zseen=(cpd*ptenh(jl,jk-1)+pgeoh(jl,jk-1))*zdmfen(jl)
2145             zqeen=pqenh(jl,jk-1)*zdmfen(jl)
2146             zsdde=(cpd*ptd(jl,jk-1)+pgeoh(jl,jk-1))*zdmfde(jl)
2147             zqdde=pqd(jl,jk-1)*zdmfde(jl)
2148             zmfdsk=pmfds(jl,jk-1)+zseen-zsdde
2149             zmfdqk=pmfdq(jl,jk-1)+zqeen-zqdde
2150             pqd(jl,jk)=zmfdqk*(1./min(-cmfcmin,pmfd(jl,jk)))
2151             ptd(jl,jk)=(zmfdsk*(1./min(-cmfcmin,pmfd(jl,jk)))- &
2152                        pgeoh(jl,jk))*rcpd
2153             ptd(jl,jk)=min(400.,ptd(jl,jk))
2154             ptd(jl,jk)=max(100.,ptd(jl,jk))
2155             zcond(jl)=pqd(jl,jk)
2156          end if
2157       end do
2158
2159       ik=jk
2160       icall=2
2161       call cuadjtq(klon,klev,ik,zph,ptd,pqd,llo2,icall)
2162       do jl=1,klon
2163          if(llo2(jl)) then
2164             zcond(jl)=zcond(jl)-pqd(jl,jk)
2165             zbuo=ptd(jl,jk)*(1.+vtmpc1*pqd(jl,jk))- &
2166            ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2167             if(zbuo.ge.0..or.prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then
2168                pmfd(jl,jk)=0.
2169             endif
2170             pmfds(jl,jk)=(cpd*ptd(jl,jk)+pgeoh(jl,jk))*pmfd(jl,jk)
2171             pmfdq(jl,jk)=pqd(jl,jk)*pmfd(jl,jk)
2172             zdmfdp=-pmfd(jl,jk)*zcond(jl)
2173             pdmfdp(jl,jk-1)=zdmfdp
2174             prfl(jl)=prfl(jl)+zdmfdp
2175          end if
2176       end do
2177
2178       if(lmfdudv) then
2179           do jl=1,klon
2180              if(llo2(jl).and.pmfd(jl,jk).lt.0.) then
2181                 zmfduk=pmfd(jl,jk-1)*pud(jl,jk-1)+   &
2182                zdmfen(jl)*puen(jl,jk-1)-zdmfde(jl)*pud(jl,jk-1)
2183                 zmfdvk=pmfd(jl,jk-1)*pvd(jl,jk-1)+   &
2184                zdmfen(jl)*pven(jl,jk-1)-zdmfde(jl)*pvd(jl,jk-1)
2185                 pud(jl,jk)=zmfduk*(1./min(-cmfcmin,pmfd(jl,jk)))
2186                 pvd(jl,jk)=zmfdvk*(1./min(-cmfcmin,pmfd(jl,jk)))
2187              end if
2188           end do
2189        end if
2190
2191       end do
2192       return
2193       end subroutine cuddraf
2194 !
2195
2196 !**********************************************
2197 !       subroutine cuflx
2198 !**********************************************
2199       subroutine cuflx &
2200          (klon,     klev,     klevp1,   pqen,    pqsen,     &
2201           ptenh,    pqenh,    paph,     pgeoh,   kcbot,    &
2202           kctop,    kdtop,    ktype,    lddraf,  ldcum,  &
2203           pmfu,     pmfd,     pmfus,    pmfds,   pmfuq,  &
2204           pmfdq,    pmful,    plude,    pdmfup,  pdmfdp, &
2205           prfl,     prain,    pten,     psfl,    pdpmel, &
2206           ktopm2,   ztmst,    sig1)
2207 !      m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
2208 !***purpose
2209 !   -------
2210 !          this routine does the final calculation of convective
2211 !          fluxes in the cloud layer and in the subcloud layer
2212 !***interface
2213 !   ---------
2214 !          this routine is called from *cumastr*.
2215 !***externals
2216 !   ---------
2217 !          none
2218 ! ----------------------------------------------------------------
2219 !-------------------------------------------------------------------
2220       implicit none
2221 !-------------------------------------------------------------------
2222       integer   klon, klev, klevp1
2223       integer   ktopm2, itop, jl, jk, ikb
2224       real      ztmst, zcons1, zcons2, zcucov, ztmelp2
2225       real      zzp, zfac, zsnmlt, zrfl, cevapcu, zrnew
2226       real      zrmin, zrfln, zdrfl, zdpevap
2227       real     pqen(klon,klev),        pqsen(klon,klev),  &
2228               ptenh(klon,klev),       pqenh(klon,klev),   &
2229               paph(klon,klevp1),      pgeoh(klon,klev)
2230       real     pmfu(klon,klev),        pmfd(klon,klev),   &
2231               pmfus(klon,klev),       pmfds(klon,klev),   &
2232               pmfuq(klon,klev),       pmfdq(klon,klev),   &
2233               pdmfup(klon,klev),      pdmfdp(klon,klev),  &
2234               pmful(klon,klev),       plude(klon,klev),   &
2235               prfl(klon),             prain(klon)
2236       real     pten(klon,klev),        pdpmel(klon,klev), &
2237               psfl(klon),             zpsubcl(klon)
2238       real     sig1(klev)
2239       integer  kcbot(klon),            kctop(klon),     &
2240               kdtop(klon),            ktype(klon)
2241       logical  lddraf(klon),           ldcum(klon)
2242 !*       specify constants
2243       zcons1=cpd/(alf*g*ztmst)
2244       zcons2=1./(g*ztmst)
2245       zcucov=0.05
2246       ztmelp2=tmelt+2.
2247 !*  1.0      determine final convective fluxes
2248 !---------------------------------------------
2249       itop=klev
2250       do jl=1,klon
2251       prfl(jl)=0.
2252       psfl(jl)=0.
2253       prain(jl)=0.
2254 !     switch off shallow convection
2255       if(.not.lmfscv.and.ktype(jl).eq.2)then
2256         ldcum(jl)=.false.
2257         lddraf(jl)=.false.
2258       endif
2259       itop=min(itop,kctop(jl))
2260       if(.not.ldcum(jl).or.kdtop(jl).lt.kctop(jl)) lddraf(jl)=.false.
2261       if(.not.ldcum(jl)) ktype(jl)=0
2262       end do
2263
2264       ktopm2=itop-2
2265       do jk=ktopm2,klev
2266       do jl=1,klon
2267       if(ldcum(jl).and.jk.ge.kctop(jl)-1) then
2268          pmfus(jl,jk)=pmfus(jl,jk)-pmfu(jl,jk)*  &
2269                      (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
2270          pmfuq(jl,jk)=pmfuq(jl,jk)-pmfu(jl,jk)*pqenh(jl,jk)
2271          if(lddraf(jl).and.jk.ge.kdtop(jl)) then
2272             pmfds(jl,jk)=pmfds(jl,jk)-pmfd(jl,jk)*  &
2273                         (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
2274             pmfdq(jl,jk)=pmfdq(jl,jk)-pmfd(jl,jk)*pqenh(jl,jk)
2275          else
2276             pmfd(jl,jk)=0.
2277             pmfds(jl,jk)=0.
2278             pmfdq(jl,jk)=0.
2279             pdmfdp(jl,jk-1)=0.
2280          end if
2281       else
2282          pmfu(jl,jk)=0.
2283          pmfd(jl,jk)=0.
2284          pmfus(jl,jk)=0.
2285          pmfds(jl,jk)=0.
2286          pmfuq(jl,jk)=0.
2287          pmfdq(jl,jk)=0.
2288          pmful(jl,jk)=0.
2289          pdmfup(jl,jk-1)=0.
2290          pdmfdp(jl,jk-1)=0.
2291          plude(jl,jk-1)=0.
2292       end if
2293       end do
2294       end do
2295
2296       do jk=ktopm2,klev
2297       do jl=1,klon
2298       if(ldcum(jl).and.jk.gt.kcbot(jl)) then
2299          ikb=kcbot(jl)
2300          zzp=((paph(jl,klevp1)-paph(jl,jk))/  &
2301              (paph(jl,klevp1)-paph(jl,ikb)))
2302          if(ktype(jl).eq.3) then
2303             zzp=zzp**2
2304          endif
2305          pmfu(jl,jk)=pmfu(jl,ikb)*zzp
2306          pmfus(jl,jk)=pmfus(jl,ikb)*zzp
2307          pmfuq(jl,jk)=pmfuq(jl,ikb)*zzp
2308          pmful(jl,jk)=pmful(jl,ikb)*zzp
2309       end if
2310 !*    2.        calculate rain/snow fall rates
2311 !*              calculate melting of snow
2312 !*              calculate evaporation of precip
2313 !----------------------------------------------
2314       if(ldcum(jl)) then
2315          prain(jl)=prain(jl)+pdmfup(jl,jk)
2316          if(pten(jl,jk).gt.tmelt) then
2317             prfl(jl)=prfl(jl)+pdmfup(jl,jk)+pdmfdp(jl,jk)
2318             if(psfl(jl).gt.0..and.pten(jl,jk).gt.ztmelp2) then
2319                zfac=zcons1*(paph(jl,jk+1)-paph(jl,jk))
2320                zsnmlt=min(psfl(jl),zfac*(pten(jl,jk)-ztmelp2))
2321                pdpmel(jl,jk)=zsnmlt
2322                psfl(jl)=psfl(jl)-zsnmlt
2323                prfl(jl)=prfl(jl)+zsnmlt
2324             end if
2325          else
2326             psfl(jl)=psfl(jl)+pdmfup(jl,jk)+pdmfdp(jl,jk)
2327          end if
2328       end if
2329       end do
2330       end do
2331
2332       do jl=1,klon
2333         prfl(jl)=max(prfl(jl),0.)
2334         psfl(jl)=max(psfl(jl),0.)
2335         zpsubcl(jl)=prfl(jl)+psfl(jl)
2336       end do
2337
2338       do jk=ktopm2,klev
2339       do jl=1,klon
2340       if(ldcum(jl).and.jk.ge.kcbot(jl).and. &
2341              zpsubcl(jl).gt.1.e-20) then
2342           zrfl=zpsubcl(jl)
2343           cevapcu=cevapcu1*sqrt(cevapcu2*sqrt(sig1(jk)))
2344           zrnew=(max(0.,sqrt(zrfl/zcucov)-   &
2345                   cevapcu*(paph(jl,jk+1)-paph(jl,jk))* &
2346                 max(0.,pqsen(jl,jk)-pqen(jl,jk))))**2*zcucov
2347           zrmin=zrfl-zcucov*max(0.,0.8*pqsen(jl,jk)-pqen(jl,jk)) &
2348                *zcons2*(paph(jl,jk+1)-paph(jl,jk))
2349           zrnew=max(zrnew,zrmin)
2350           zrfln=max(zrnew,0.)
2351           zdrfl=min(0.,zrfln-zrfl)
2352           pdmfup(jl,jk)=pdmfup(jl,jk)+zdrfl
2353           zpsubcl(jl)=zrfln
2354       end if
2355       end do
2356       end do
2357
2358       do jl=1,klon
2359         zdpevap=zpsubcl(jl)-(prfl(jl)+psfl(jl))
2360         prfl(jl)=prfl(jl)+zdpevap*prfl(jl)*  &
2361                   (1./max(1.e-20,prfl(jl)+psfl(jl)))
2362         psfl(jl)=psfl(jl)+zdpevap*psfl(jl)*  &
2363                   (1./max(1.e-20,prfl(jl)+psfl(jl)))
2364       end do
2365
2366       return
2367       end subroutine cuflx
2368 !
2369
2370 !**********************************************
2371 !       subroutine cudtdq
2372 !**********************************************
2373       subroutine cudtdq &
2374          (klon,     klev,     klevp1,   ktopm2,   paph,   &
2375           ldcum,    pten,     ptte,     pqte,     pmfus,  &
2376           pmfds,    pmfuq,    pmfdq,    pmful,    pdmfup, &
2377           pdmfdp,   ztmst,    pdpmel,   prain,    prfl,   &
2378           psfl,     psrain,   psevap,   psheat,   psmelt, &
2379           prsfc,    pssfc,    paprc,    paprsm,   paprs,  &
2380           pqen,     pqsen,    plude,    pcte)
2381 !**** *cudtdq* - updates t and q tendencies, precipitation rates
2382 !                does global diagnostics
2383 !      m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
2384 !***interface.
2385 !   ----------
2386 !          *cudtdq* is called from *cumastr*
2387 ! ----------------------------------------------------------------
2388 !-------------------------------------------------------------------
2389       implicit none
2390 !-------------------------------------------------------------------
2391       integer   klon, klev, klevp1
2392       integer   ktopm2,jl, jk
2393       real      ztmst, psrain, psevap, psheat, psmelt, zdiagt, zdiagw
2394       real      zalv, rhk, rhcoe, pldfd, zdtdt, zdqdt
2395       real     ptte(klon,klev),        pqte(klon,klev),  &
2396               pten(klon,klev),        plude(klon,klev),  &
2397               pgeo(klon,klev),        paph(klon,klevp1), &
2398               paprc(klon),            paprs(klon),       &
2399               paprsm(klon),           pcte(klon,klev),   &
2400               prsfc(klon),            pssfc(klon)
2401       real     pmfus(klon,klev),       pmfds(klon,klev), &
2402               pmfuq(klon,klev),       pmfdq(klon,klev), &
2403               pmful(klon,klev),       pqsen(klon,klev), &
2404               pdmfup(klon,klev),      pdmfdp(klon,klev),&
2405               prfl(klon),             prain(klon),      &
2406               pqen(klon,klev)
2407       real     pdpmel(klon,klev),      psfl(klon)
2408       real     zsheat(klon),           zmelt(klon)
2409       logical  ldcum(klon)
2410 !--------------------------------
2411 !*    1.0      specify parameters
2412 !--------------------------------
2413       zdiagt=ztmst
2414       zdiagw=zdiagt/rhoh2o
2415 !--------------------------------------------------
2416 !*    2.0      incrementation of t and q tendencies
2417 !--------------------------------------------------
2418       do jl=1,klon
2419         zmelt(jl)=0.
2420         zsheat(jl)=0.
2421       end do
2422
2423       do jk=ktopm2,klev
2424       if(jk.lt.klev) then
2425          do jl=1,klon
2426          if(ldcum(jl)) then
2427             if(pten(jl,jk).gt.tmelt) then
2428                zalv=alv
2429             else
2430                zalv=als
2431             endif
2432             rhk=min(1.0,pqen(jl,jk)/pqsen(jl,jk))
2433             rhcoe=max(0.0,(rhk-rhc)/(rhm-rhc))
2434             pldfd=max(0.0,rhcoe*fdbk*plude(jl,jk))
2435             zdtdt=(g/(paph(jl,jk+1)-paph(jl,jk)))*rcpd*      &
2436               (pmfus(jl,jk+1)-pmfus(jl,jk)+                  &
2437               pmfds(jl,jk+1)-pmfds(jl,jk)-alf*pdpmel(jl,jk)  &
2438               -zalv*(pmful(jl,jk+1)-pmful(jl,jk)-pldfd-      &
2439               (pdmfup(jl,jk)+pdmfdp(jl,jk))))
2440             ptte(jl,jk)=ptte(jl,jk)+zdtdt
2441             zdqdt=(g/(paph(jl,jk+1)-paph(jl,jk)))*&
2442               (pmfuq(jl,jk+1)-pmfuq(jl,jk)+       &
2443               pmfdq(jl,jk+1)-pmfdq(jl,jk)+        &
2444               pmful(jl,jk+1)-pmful(jl,jk)-pldfd-  &
2445               (pdmfup(jl,jk)+pdmfdp(jl,jk)))
2446             pqte(jl,jk)=pqte(jl,jk)+zdqdt
2447             pcte(jl,jk)=(g/(paph(jl,jk+1)-paph(jl,jk)))*pldfd
2448             zsheat(jl)=zsheat(jl)+zalv*(pdmfup(jl,jk)+pdmfdp(jl,jk))
2449             zmelt(jl)=zmelt(jl)+pdpmel(jl,jk)
2450          end if
2451          end do
2452       else
2453          do jl=1,klon
2454          if(ldcum(jl)) then
2455             if(pten(jl,jk).gt.tmelt) then
2456                zalv=alv
2457             else
2458                zalv=als
2459             endif
2460             rhk=min(1.0,pqen(jl,jk)/pqsen(jl,jk))
2461             rhcoe=max(0.0,(rhk-rhc)/(rhm-rhc))
2462             pldfd=max(0.0,rhcoe*fdbk*plude(jl,jk))
2463             zdtdt=-(g/(paph(jl,jk+1)-paph(jl,jk)))*rcpd*           &
2464                 (pmfus(jl,jk)+pmfds(jl,jk)+alf*pdpmel(jl,jk)-zalv* &
2465                 (pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)+pldfd))
2466             ptte(jl,jk)=ptte(jl,jk)+zdtdt
2467             zdqdt=-(g/(paph(jl,jk+1)-paph(jl,jk)))*                &
2468                      (pmfuq(jl,jk)+pmfdq(jl,jk)+pldfd+             &
2469                      (pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)))
2470             pqte(jl,jk)=pqte(jl,jk)+zdqdt
2471             pcte(jl,jk)=(g/(paph(jl,jk+1)-paph(jl,jk)))*pldfd
2472             zsheat(jl)=zsheat(jl)+zalv*(pdmfup(jl,jk)+pdmfdp(jl,jk))
2473             zmelt(jl)=zmelt(jl)+pdpmel(jl,jk)
2474          end if
2475          end do
2476       end if
2477       end do
2478 !---------------------------------------------------------
2479 !      3.      update surface fields and do global budgets
2480 !---------------------------------------------------------
2481       do jl=1,klon
2482       prsfc(jl)=prfl(jl)
2483       pssfc(jl)=psfl(jl)
2484       paprc(jl)=paprc(jl)+zdiagw*(prfl(jl)+psfl(jl))
2485       paprs(jl)=paprsm(jl)+zdiagw*psfl(jl)
2486       psheat=psheat+zsheat(jl)
2487       psrain=psrain+prain(jl)
2488       psevap=psevap-(prfl(jl)+psfl(jl))
2489       psmelt=psmelt+zmelt(jl)
2490       end do
2491       psevap=psevap+psrain
2492       return
2493       end subroutine cudtdq
2494
2495 !
2496 !**********************************************
2497 !       subroutine cududv
2498 !**********************************************
2499       subroutine cududv &
2500          (klon,     klev,     klevp1,   ktopm2,   ktype,  &
2501           kcbot,    paph,     ldcum,    puen,     pven,   &
2502           pvom,     pvol,     puu,      pud,      pvu,    &
2503           pvd,      pmfu,     pmfd,     psdiss)
2504 !**** *cududv* - updates u and v tendencies,
2505 !                does global diagnostic of dissipation
2506 !      m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
2507 !***interface.
2508 !   ----------
2509 !          *cududv* is called from *cumastr*
2510 ! ----------------------------------------------------------------
2511 !-------------------------------------------------------------------
2512       implicit none
2513 !-------------------------------------------------------------------
2514       integer   klon, klev, klevp1
2515       integer   ktopm2, jk, ik, jl, ikb
2516       real      psdiss,zzp, zdudt ,zdvdt, zsum
2517       real     puen(klon,klev),        pven(klon,klev),   &
2518               pvol(klon,klev),        pvom(klon,klev),    &
2519               paph(klon,klevp1)
2520       real     puu(klon,klev),         pud(klon,klev),    &
2521               pvu(klon,klev),         pvd(klon,klev),     &
2522               pmfu(klon,klev),        pmfd(klon,klev)
2523       real     zmfuu(klon,klev),       zmfdu(klon,klev),  &
2524               zmfuv(klon,klev),       zmfdv(klon,klev),   &
2525               zdiss(klon)
2526       integer  ktype(klon),            kcbot(klon)
2527       logical  ldcum(klon)
2528 !------------------------------------------------------------
2529 !*    1.0      calculate fluxes and update u and v tendencies
2530 ! -----------------------------------------------------------
2531       do jk=ktopm2,klev
2532       ik=jk-1
2533       do jl=1,klon
2534       if(ldcum(jl)) then
2535         zmfuu(jl,jk)=pmfu(jl,jk)*(puu(jl,jk)-puen(jl,ik))
2536         zmfuv(jl,jk)=pmfu(jl,jk)*(pvu(jl,jk)-pven(jl,ik))
2537         zmfdu(jl,jk)=pmfd(jl,jk)*(pud(jl,jk)-puen(jl,ik))
2538         zmfdv(jl,jk)=pmfd(jl,jk)*(pvd(jl,jk)-pven(jl,ik))
2539       end if
2540       end do
2541       end do
2542
2543       do jk=ktopm2,klev
2544       do jl=1,klon
2545       if(ldcum(jl).and.jk.gt.kcbot(jl)) then
2546          ikb=kcbot(jl)
2547          zzp=((paph(jl,klevp1)-paph(jl,jk))/  &
2548              (paph(jl,klevp1)-paph(jl,ikb)))
2549          if(ktype(jl).eq.3) then
2550             zzp=zzp**2
2551          endif
2552          zmfuu(jl,jk)=zmfuu(jl,ikb)*zzp
2553          zmfuv(jl,jk)=zmfuv(jl,ikb)*zzp
2554          zmfdu(jl,jk)=zmfdu(jl,ikb)*zzp
2555          zmfdv(jl,jk)=zmfdv(jl,ikb)*zzp
2556       end if
2557       end do
2558       end do
2559
2560       do jl=1,klon
2561       zdiss(jl)=0.
2562       end do
2563
2564       do jk=ktopm2,klev
2565       if(jk.lt.klev) then
2566          do jl=1,klon
2567             if(ldcum(jl)) then
2568                zdudt=(g/(paph(jl,jk+1)-paph(jl,jk)))* &
2569                     (zmfuu(jl,jk+1)-zmfuu(jl,jk)+     &
2570                      zmfdu(jl,jk+1)-zmfdu(jl,jk))
2571                zdvdt=(g/(paph(jl,jk+1)-paph(jl,jk)))* &
2572                     (zmfuv(jl,jk+1)-zmfuv(jl,jk)+     &
2573                      zmfdv(jl,jk+1)-zmfdv(jl,jk))
2574                zdiss(jl)=zdiss(jl)+        &
2575                         puen(jl,jk)*(zmfuu(jl,jk+1)-zmfuu(jl,jk)+   &
2576                                      zmfdu(jl,jk+1)-zmfdu(jl,jk))+  &
2577                         pven(jl,jk)*(zmfuv(jl,jk+1)-zmfuv(jl,jk)+   &
2578                                      zmfdv(jl,jk+1)-zmfdv(jl,jk))
2579                pvom(jl,jk)=pvom(jl,jk)+zdudt
2580                pvol(jl,jk)=pvol(jl,jk)+zdvdt
2581             end if
2582          end do
2583       else
2584          do jl=1,klon
2585             if(ldcum(jl)) then
2586                zdudt=-(g/(paph(jl,jk+1)-paph(jl,jk)))* &
2587                         (zmfuu(jl,jk)+zmfdu(jl,jk))
2588                zdvdt=-(g/(paph(jl,jk+1)-paph(jl,jk)))* &
2589                         (zmfuv(jl,jk)+zmfdv(jl,jk))
2590                zdiss(jl)=zdiss(jl)-        &
2591                          (puen(jl,jk)*(zmfuu(jl,jk)+zmfdu(jl,jk))+ &
2592                           pven(jl,jk)*(zmfuv(jl,jk)+zmfdv(jl,jk)))
2593                pvom(jl,jk)=pvom(jl,jk)+zdudt
2594                pvol(jl,jk)=pvol(jl,jk)+zdvdt
2595             end if
2596          end do
2597        end if
2598       end do
2599       zsum=ssum(klon,zdiss(1),1)
2600       psdiss=psdiss+zsum
2601       return
2602       end subroutine cududv
2603 !
2604
2605 !#################################################################
2606 !
2607 !                 level 4 subroutines
2608 !
2609 !#################################################################
2610 !**************************************************************
2611 !             subroutine cubasmc
2612 !**************************************************************
2613       subroutine cubasmc   &
2614          (klon,     klev,     klevm1,  kk,     pten,  &
2615           pqen,     pqsen,    puen,    pven,   pverv, &
2616           pgeo,     pgeoh,    ldcum,   ktype,  klab,  &
2617           pmfu,     pmfub,    pentr,   kcbot,  ptu,   &
2618           pqu,      plu,      puu,     pvu,    pmfus, &
2619           pmfuq,    pmful,    pdmfup,  pmfuu,  pmfuv)
2620 !      m.tiedtke         e.c.m.w.f.     12/89
2621 !***purpose.
2622 !   --------
2623 !          this routine calculates cloud base values
2624 !          for midlevel convection
2625 !***interface
2626 !   ---------
2627 !          this routine is called from *cuasc*.
2628 !          input are environmental values t,q etc
2629 !          it returns cloudbase values for midlevel convection
2630 !***method.
2631 !   -------
2632 !          s. tiedtke (1989)
2633 !***externals
2634 !   ---------
2635 !          none
2636 ! ----------------------------------------------------------------
2637 !-------------------------------------------------------------------
2638       implicit none
2639 !-------------------------------------------------------------------
2640       integer   klon, klev, klevp1
2641       integer   klevm1,kk, jl
2642       real      zzzmb
2643       real     pten(klon,klev),        pqen(klon,klev),  &
2644               puen(klon,klev),        pven(klon,klev),   &
2645               pqsen(klon,klev),       pverv(klon,klev),  &
2646               pgeo(klon,klev),        pgeoh(klon,klev)
2647       real     ptu(klon,klev),         pqu(klon,klev),   &
2648               puu(klon,klev),         pvu(klon,klev),    &
2649               plu(klon,klev),         pmfu(klon,klev),   &
2650               pmfub(klon),            pentr(klon),       &
2651               pmfus(klon,klev),       pmfuq(klon,klev),  &
2652               pmful(klon,klev),       pdmfup(klon,klev), &
2653               pmfuu(klon),            pmfuv(klon)
2654       integer  ktype(klon),            kcbot(klon),      &
2655               klab(klon,klev)
2656       logical  ldcum(klon)
2657 !--------------------------------------------------------
2658 !*    1.      calculate entrainment and detrainment rates
2659 ! -------------------------------------------------------
2660          do jl=1,klon
2661           if( .not. ldcum(jl).and.klab(jl,kk+1).eq.0.0.and.  &
2662              pqen(jl,kk).gt.0.80*pqsen(jl,kk)) then
2663             ptu(jl,kk+1)=(cpd*pten(jl,kk)+pgeo(jl,kk)-pgeoh(jl,kk+1)) &
2664                                *rcpd
2665             pqu(jl,kk+1)=pqen(jl,kk)
2666             plu(jl,kk+1)=0.
2667             zzzmb=max(cmfcmin,-pverv(jl,kk)/g)
2668             zzzmb=min(zzzmb,cmfcmax)
2669             pmfub(jl)=zzzmb
2670             pmfu(jl,kk+1)=pmfub(jl)
2671             pmfus(jl,kk+1)=pmfub(jl)*(cpd*ptu(jl,kk+1)+pgeoh(jl,kk+1))
2672             pmfuq(jl,kk+1)=pmfub(jl)*pqu(jl,kk+1)
2673             pmful(jl,kk+1)=0.
2674             pdmfup(jl,kk+1)=0.
2675             kcbot(jl)=kk
2676             klab(jl,kk+1)=1
2677             ktype(jl)=3
2678             pentr(jl)=entrmid
2679                if(lmfdudv) then
2680                   puu(jl,kk+1)=puen(jl,kk)
2681                   pvu(jl,kk+1)=pven(jl,kk)
2682                   pmfuu(jl)=pmfub(jl)*puu(jl,kk+1)
2683                   pmfuv(jl)=pmfub(jl)*pvu(jl,kk+1)
2684                end if
2685          end if
2686         end do
2687       return
2688       end subroutine cubasmc
2689
2690 !
2691 !**************************************************************
2692 !             subroutine cuadjtq
2693 !**************************************************************
2694       subroutine cuadjtq(klon,klev,kk,pp,pt,pq,ldflag,kcall)
2695 !      m.tiedtke         e.c.m.w.f.     12/89
2696 !      d.salmond         cray(uk))      12/8/91
2697 !***purpose.
2698 !   --------
2699 !          to produce t,q and l values for cloud ascent
2700 !***interface
2701 !   ---------
2702 !          this routine is called from subroutines:
2703 !              *cubase*   (t and q at condenstion level)
2704 !              *cuasc*    (t and q at cloud levels)
2705 !              *cuini*    (environmental t and qs values at half levels)
2706 !          input are unadjusted t and q values,
2707 !          it returns adjusted values of t and q
2708 !          note: input parameter kcall defines calculation as
2709 !               kcall=0    env. t and qs in*cuini*
2710 !               kcall=1  condensation in updrafts  (e.g.  cubase, cuasc)
2711 !               kcall=2  evaporation in downdrafts (e.g.  cudlfs,cuddraf
2712 !***externals
2713 !   ---------
2714 !          3 lookup tables ( tlucua, tlucub, tlucuc )
2715 !          for condensation calculations.
2716 !          the tables are initialised in *setphys*.
2717 ! ----------------------------------------------------------------
2718 !-------------------------------------------------------------------
2719       implicit none
2720 !-------------------------------------------------------------------
2721       integer   klon, klev
2722       integer   kk, kcall, isum, jl
2723       real      zqsat, zcor, zcond1, tt
2724       real     pt(klon,klev),          pq(klon,klev),  &
2725               zcond(klon),            zqp(klon),       &
2726               pp(klon)
2727       logical  ldflag(klon)
2728 !------------------------------------------------------------------
2729 !     2.      calculate condensation and adjust t and q accordingly
2730 !------------------------------------------------------------------
2731       if (kcall.eq.1 ) then
2732          isum=0
2733          do jl=1,klon
2734          zcond(jl)=0.
2735          if(ldflag(jl)) then
2736             zqp(jl)=1./pp(jl)
2737             tt=pt(jl,kk)
2738             zqsat=tlucua(tt)*zqp(jl)
2739             zqsat=min(0.5,zqsat)
2740             zcor=1./(1.-vtmpc1*zqsat)
2741             zqsat=zqsat*zcor
2742             zcond(jl)=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2743             zcond(jl)=max(zcond(jl),0.)
2744             pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond(jl)
2745             pq(jl,kk)=pq(jl,kk)-zcond(jl)
2746             if(zcond(jl).ne.0.0) isum=isum+1
2747          end if
2748          end do
2749
2750          if(isum.eq.0) return
2751          do jl=1,klon
2752          if(ldflag(jl).and.zcond(jl).ne.0.) then
2753             tt=pt(jl,kk)
2754             zqsat=tlucua(tt)*zqp(jl)
2755             zqsat=min(0.5,zqsat)
2756             zcor=1./(1.-vtmpc1*zqsat)
2757             zqsat=zqsat*zcor
2758             zcond1=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2759             pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond1
2760             pq(jl,kk)=pq(jl,kk)-zcond1
2761          end if
2762          end do
2763       end if
2764       if(kcall.eq.2) then
2765          isum=0
2766          do jl=1,klon
2767          zcond(jl)=0.
2768          if(ldflag(jl)) then
2769             tt=pt(jl,kk)
2770             zqp(jl)=1./pp(jl)
2771             zqsat=tlucua(tt)*zqp(jl)
2772             zqsat=min(0.5,zqsat)
2773             zcor=1./(1.-vtmpc1*zqsat)
2774             zqsat=zqsat*zcor
2775             zcond(jl)=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2776             zcond(jl)=min(zcond(jl),0.)
2777             pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond(jl)
2778             pq(jl,kk)=pq(jl,kk)-zcond(jl)
2779             if(zcond(jl).ne.0.0) isum=isum+1
2780          end if
2781          end do
2782
2783          if(isum.eq.0) return
2784          do jl=1,klon
2785          if(ldflag(jl).and.zcond(jl).ne.0.) then
2786             tt=pt(jl,kk)
2787             zqsat=tlucua(tt)*zqp(jl)
2788             zqsat=min(0.5,zqsat)
2789             zcor=1./(1.-vtmpc1*zqsat)
2790             zqsat=zqsat*zcor
2791             zcond1=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2792             pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond1
2793             pq(jl,kk)=pq(jl,kk)-zcond1
2794          end if
2795          end do
2796       end if
2797       if(kcall.eq.0) then
2798          isum=0
2799          do jl=1,klon
2800            tt=pt(jl,kk)
2801            zqp(jl)=1./pp(jl)
2802            zqsat=tlucua(tt)*zqp(jl)
2803            zqsat=min(0.5,zqsat)
2804            zcor=1./(1.-vtmpc1*zqsat)
2805            zqsat=zqsat*zcor
2806            zcond(jl)=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2807            pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond(jl)
2808            pq(jl,kk)=pq(jl,kk)-zcond(jl)
2809            if(zcond(jl).ne.0.0) isum=isum+1
2810          end do
2811
2812          if(isum.eq.0) return
2813          do jl=1,klon
2814            tt=pt(jl,kk)
2815            zqsat=tlucua(tt)*zqp(jl)
2816            zqsat=min(0.5,zqsat)
2817            zcor=1./(1.-vtmpc1*zqsat)
2818            zqsat=zqsat*zcor
2819            zcond1=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2820            pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond1
2821            pq(jl,kk)=pq(jl,kk)-zcond1
2822          end do
2823       end if
2824       if(kcall.eq.4) then
2825          do jl=1,klon
2826            tt=pt(jl,kk)
2827            zqp(jl)=1./pp(jl)
2828            zqsat=tlucua(tt)*zqp(jl)
2829            zqsat=min(0.5,zqsat)
2830            zcor=1./(1.-vtmpc1*zqsat)
2831            zqsat=zqsat*zcor
2832            zcond(jl)=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2833            pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond(jl)
2834            pq(jl,kk)=pq(jl,kk)-zcond(jl)
2835          end do
2836
2837          do jl=1,klon
2838            tt=pt(jl,kk)
2839            zqsat=tlucua(tt)*zqp(jl)
2840            zqsat=min(0.5,zqsat)
2841            zcor=1./(1.-vtmpc1*zqsat)
2842            zqsat=zqsat*zcor
2843            zcond1=(pq(jl,kk)-zqsat)/(1.+zqsat*zcor*tlucub(tt))
2844            pt(jl,kk)=pt(jl,kk)+tlucuc(tt)*zcond1
2845            pq(jl,kk)=pq(jl,kk)-zcond1
2846          end do
2847       end if
2848       return
2849       end subroutine cuadjtq
2850
2851 !
2852 !**********************************************************
2853 !        subroutine cuentr_new
2854 !**********************************************************
2855       subroutine cuentr_new                              &
2856          (klon,     klev,     klevp1,   kk,       ptenh, &
2857           paph,     pap,      pgeoh,    klwmin,   ldcum, &
2858           ktype,    kcbot,    kctop0,   zpbase,   pmfu,  &
2859           pentr,    zdmfen,   zdmfde,   zodetr,   khmin)
2860 !      m.tiedtke         e.c.m.w.f.     12/89
2861 !      y.wang            iprc           11/01
2862 !***purpose.
2863 !   --------
2864 !          this routine calculates entrainment/detrainment rates
2865 !          for updrafts in cumulus parameterization
2866 !***interface
2867 !   ---------
2868 !          this routine is called from *cuasc*.
2869 !          input are environmental values t,q etc
2870 !          and updraft values t,q etc
2871 !          it returns entrainment/detrainment rates
2872 !***method.
2873 !  --------
2874 !          s. tiedtke (1989), nordeng(1996)
2875 !***externals
2876 !   ---------
2877 !          none
2878 ! ----------------------------------------------------------------
2879 !-------------------------------------------------------------------
2880       implicit none
2881 !-------------------------------------------------------------------
2882       integer   klon, klev, klevp1
2883       integer   kk, jl, iklwmin,ikb, ikt, ikh
2884       real      zrrho, zdprho, zpmid, zentr, zzmzk, ztmzk, arg, zorgde
2885       real     ptenh(klon,klev),                           &
2886               pap(klon,klev),         paph(klon,klevp1),   &
2887               pmfu(klon,klev),        pgeoh(klon,klev),    &
2888               pentr(klon),            zpbase(klon),        &
2889               zdmfen(klon),           zdmfde(klon),        &
2890               zodetr(klon,klev)
2891       integer  klwmin(klon),           ktype(klon),        &
2892               kcbot(klon),            kctop0(klon),        &
2893               khmin(klon)
2894       logical  ldcum(klon),llo1,llo2
2895 !---------------------------------------------------------
2896 !*    1.       calculate entrainment and detrainment rates
2897 !---------------------------------------------------------
2898 !*    1.1      specify entrainment rates for shallow clouds
2899 !----------------------------------------------------------
2900 !*    1.2      specify entrainment rates for deep clouds
2901 !-------------------------------------------------------
2902       do jl = 1, klon
2903         zpbase(jl) = paph(jl,kcbot(jl))
2904         zrrho = (rd*ptenh(jl,kk+1))/paph(jl,kk+1)
2905         zdprho = (paph(jl,kk+1)-paph(jl,kk))*zrg
2906         zpmid = 0.5*(zpbase(jl)+paph(jl,kctop0(jl)))
2907         zentr = pentr(jl)*pmfu(jl,kk+1)*zdprho*zrrho
2908         llo1 = kk.lt.kcbot(jl).and.ldcum(jl)
2909         if(llo1) then
2910            zdmfde(jl) = zentr
2911         else
2912            zdmfde(jl) = 0.0
2913         endif
2914         llo2 = llo1.and.ktype(jl).eq.2.and.((zpbase(jl)-paph(jl,kk)) &
2915              .lt.zdnoprc.or.paph(jl,kk).gt.zpmid)
2916         if(llo2) then
2917             zdmfen(jl) = zentr
2918         else
2919             zdmfen(jl) = 0.0
2920         endif
2921         iklwmin = max(klwmin(jl),kctop0(jl)+2)
2922         llo2 = llo1.and.ktype(jl).eq.3.and.(kk.ge.iklwmin.or.pap(jl,kk) &
2923              .gt.zpmid)
2924         if (llo2) zdmfen(jl) = zentr
2925         llo2 = llo1.and.ktype(jl).eq.1
2926 ! turbulent entrainment
2927         if (llo2) zdmfen(jl) = zentr
2928 ! organized detrainment, detrainment starts at khmin
2929         ikb = kcbot(jl)
2930         zodetr(jl,kk) = 0.
2931         if (llo2.and.kk.le.khmin(jl).and.kk.ge.kctop0(jl)) then
2932           ikt = kctop0(jl)
2933           ikh = khmin(jl)
2934           if (ikh.gt.ikt) then
2935             zzmzk = -(pgeoh(jl,ikh)-pgeoh(jl,kk))*zrg
2936             ztmzk = -(pgeoh(jl,ikh)-pgeoh(jl,ikt))*zrg
2937             arg = 3.1415*(zzmzk/ztmzk)*0.5
2938 #ifndef AARCH64_X86_CORRECTNESS_FIX
2939             zorgde = tan(arg)*3.1415*0.5/ztmzk
2940 #else
2941             zorgde = real(tan(dble(arg))*3.1415*0.5/dble(ztmzk))
2942 #endif
2943             zdprho = (paph(jl,kk+1)-paph(jl,kk))*(zrg*zrrho)
2944             zodetr(jl,kk) = min(zorgde,1.e-3)*pmfu(jl,kk+1)*zdprho
2945           end if
2946         end if
2947       enddo
2948 !
2949       return
2950       end subroutine cuentr_new
2951 !
2952
2953 !**********************************************************
2954 !        function ssum, tlucua, tlucub, tlucuc
2955 !**********************************************************
2956       real function ssum ( n, x, ix )
2957 !
2958 ! computes ssum = sum of [x(i)]
2959 !     for n elements of x with skip increment ix for vector x
2960 !
2961       implicit none
2962       real x(*)
2963       real zsum
2964       integer n, ix, jx, jl
2965 !
2966       jx = 1
2967       zsum = 0.0
2968       do jl = 1, n
2969         zsum = zsum + x(jx)
2970         jx = jx + ix
2971       enddo
2972 !
2973       ssum=zsum
2974 !
2975       return
2976       end function ssum
2977
2978       real function tlucua(tt)
2979 !
2980 !  set up lookup tables for cloud ascent calculations.
2981 !
2982       implicit none
2983       real zcvm3,zcvm4,tt !,tlucua
2984 !
2985       if(tt-tmelt.gt.0.) then
2986          zcvm3=c3les
2987          zcvm4=c4les
2988       else
2989          zcvm3=c3ies
2990          zcvm4=c4ies
2991       end if
2992       tlucua=c2es*exp(zcvm3*(tt-tmelt)*(1./(tt-zcvm4)))
2993 !
2994       return
2995       end function tlucua
2996 !
2997       real function tlucub(tt)
2998 !
2999 !  set up lookup tables for cloud ascent calculations.
3000 !
3001       implicit none
3002       real z5alvcp,z5alscp,zcvm4,zcvm5,tt !,tlucub
3003 !
3004       z5alvcp=c5les*alv/cpd
3005       z5alscp=c5ies*als/cpd
3006       if(tt-tmelt.gt.0.) then
3007          zcvm4=c4les
3008          zcvm5=z5alvcp
3009       else
3010          zcvm4=c4ies
3011          zcvm5=z5alscp
3012       end if
3013       tlucub=zcvm5*(1./(tt-zcvm4))**2
3014 !
3015       return
3016       end function tlucub
3017 !
3018       real function tlucuc(tt)
3019 !
3020 !  set up lookup tables for cloud ascent calculations.
3021 !
3022       implicit none
3023       real zalvdcp,zalsdcp,tt,zldcp !,tlucuc
3024 !
3025       zalvdcp=alv/cpd
3026       zalsdcp=als/cpd
3027       if(tt-tmelt.gt.0.) then
3028          zldcp=zalvdcp
3029       else
3030          zldcp=zalsdcp
3031       end if
3032       tlucuc=zldcp
3033 !
3034       return
3035       end function tlucuc
3036 !
3037
3038 end module module_cu_tiedtke