phys/module_cu_ntiedtke.F

   1 !-----------------------------------------------------------------------
   2 !
   3 !wrf:model_layer:physics
   4 !
   5 !####################tiedtke scheme#########################
   6 !    m.tiedtke      e.c.m.w.f.      1989
   7 !    j.morcrette                    1992
   8 !--------------------------------------------
   9 !    modifications
  10 !    C. zhang & Yuqing Wang         2011-2017
  11 !
  12 !   modified from IPRC IRAM - yuqing wang, university of hawaii
  13 !               & ICTP REGCM4.4
  14 !
  15 !   The current version is stable.  There are many updates to the old Tiedtke scheme (cu_physics=6)
  16 !   update notes:
  17 !        the new Tiedtke scheme is similar to the Tiedtke scheme used in REGCM4 and ECMWF cy40r1.
  18 !        the major differences to the old Tiedtke (cu_physics=6) scheme are,
  19 !        (a) New trigger functions for deep and shallow convections (Jakob and Siebesma 2003;
  20 !            Bechtold et al. 2004, 2008, 2014).
  21 !        (b) Non-equilibrium situations are considered in the closure for deep convection
  22 !            (Bechtold et al. 2014).
  23 !        (c) New convection time scale for the deep convection closure (Bechtold et al. 2008).
  24 !        (d) New entrainment and detrainment rates for all convection types (Bechtold et al. 2008).
  25 !        (e) New formula for the conversion from cloud water/ice to rain/snow (Sundqvist 1978).
  26 !        (f) Different way to include cloud scale pressure gradients (Gregory et al. 1997;
  27 !            Wu and Yanai 1994)
  28 !
  29 !   other refenrence: tiedtke (1989, mwr, 117, 1779-1800)
  30 !                     IFS documentation - cy33r1, cy37r2, cy38r1, cy40r1
  31 !
  32 !===========================================================
  33 !  Note for climate simulation of Tropical Cyclones
  34 !  This version of Tiedtke scheme was tested with YSU PBL scheme, RRTMG radation
  35 !  schemes, and WSM6 microphysics schemes, at horizontal resolution around 20 km
  36 !  Set: momtrans = 2.
  37 !       pgcoef   = 0.7 to 1.0 is good depends on the basin
  38 !       nonequil = .false.
  39 !===========================================================
  40 !  Note for the diurnal simulation of precipitaton
  41 !  When nonequil = .true., the CAPE is relaxed toward to a value from PBL
  42 !  It can improve the diurnal precipitation over land.
  43 !===========================================================
  44 !###########################################################
  45
  46 module module_cu_ntiedtke
  47
  48 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  49 #if defined(mpas)
  50      use mpas_atmphys_constants, only: rd=>R_d, rv=>R_v, &
  51    &       cpd=>cp, alv=>xlv, als=>xls, alf=>xlf, g=>gravity
  52 #else
  53      use module_model_constants, only:rd=>r_d, rv=>r_v, &
  54    &       cpd=>cp, alv=>xlv, als=>xls, alf=>xlf, g
  55 #endif
  56
  57      implicit none
  58      real,private :: t13,rcpd,vtmpc1,tmelt,            &
  59              c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg
  60
  61      real,private :: r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp,rtwat,rtber,rtice
  62      real,private :: entrdd,cmfcmax,cmfcmin,cmfdeps,zdnoprc,cprcon,pgcoef
  63      integer,private :: momtrans
  64
  65      parameter(         &
  66       t13=1.0/3.0,      &
  67       rcpd=1.0/cpd,     &
  68       tmelt=273.16,     &
  69       zrg=1.0/g,        &
  70       c1es=610.78,      &
  71       c2es=c1es*rd/rv,  &
  72       c3les=17.2693882, &
  73       c3ies=21.875,     &
  74       c4les=35.86,      &
  75       c4ies=7.66,       &
  76       c5les=c3les*(tmelt-c4les),  &
  77       c5ies=c3ies*(tmelt-c4ies),  &
  78       r5alvcp=c5les*alv*rcpd,     &
  79       r5alscp=c5ies*als*rcpd,     &
  80       ralvdcp=alv*rcpd,           &
  81       ralsdcp=als*rcpd,           &
  82       ralfdcp=alf*rcpd,           &
  83       rtwat=tmelt,                &
  84       rtber=tmelt-5.,             &
  85       rtice=tmelt-23.,            &
  86       vtmpc1=rv/rd-1.0 )
  87 !
  88 !     entrdd: average entrainment & detrainment rate for downdrafts
  89 !     ------
  90 !
  91       parameter(entrdd = 2.0e-4)
  92 !
  93 !     cmfcmax:   maximum massflux value allowed for updrafts etc
  94 !     -------
  95 !
  96       parameter(cmfcmax = 1.0)
  97 !
  98 !     cmfcmin:   minimum massflux value (for safety)
  99 !     -------
 100 !
 101       parameter(cmfcmin = 1.e-10)
 102 !
 103 !     cmfdeps:   fractional massflux for downdrafts at lfs
 104 !     -------
 105 !
 106       parameter(cmfdeps = 0.30)
 107
 108 !     zdnoprc:   deep cloud is thicker than this height (Unit:Pa)
 109 !
 110       parameter(zdnoprc = 2.0e4)
 111 !     -------
 112 !
 113 !     cprcon:    coefficient from cloud water to rain water
 114 !
 115       parameter(cprcon = 1.4e-3)
 116 !     -------
 117 !
 118 !     momtrans: momentum transport method
 119 !     ( 1 = IFS40r1 method; 2 = new method )
 120 !
 121       parameter(momtrans = 2 )
 122 !     -------
 123 !
 124 !     coefficient for pressure gradient intensity
 125 !     (0.7 - 1.0 is recommended in this vesion of Tiedtke scheme)
 126       parameter(pgcoef=0.7)
 127 !     -------
 128 !
 129       logical :: nonequil
 130 !     nonequil: representing equilibrium and nonequilibrium convection
 131 !     ( .false. [equilibrium: removing all CAPE]; .true. [nonequilibrium: relaxing CAPE toward CAPE from PBL].
 132 !       Ref. Bechtold et al. 2014 JAS )
 133 !
 134       parameter(nonequil = .true. )
 135 !
 136 !--------------------
 137 !     switches for deep, mid, shallow convections, downdraft, and momentum transport
 138 !     ------------------
 139       logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv
 140       parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.)
 141 !--------------------
 142 !#################### end of variables definition##########################
 143 !-----------------------------------------------------------------------
 144 !
 145 contains
 146 !-----------------------------------------------------------------------
 147      subroutine cu_ntiedtke(                                    &
 148                  dt,itimestep,stepcu                            &
 149                 ,raincv,pratec,qfx,hfx                          &
 150                 ,u3d,v3d,w,t3d,qv3d,qc3d,qi3d,pi3d,rho3d        &
 151                 ,qvften,thften                                  &
 152                 ,dz8w,pcps,p8w,xland,cu_act_flag,dx             &
 153                 ,ids,ide, jds,jde, kds,kde                      &
 154                 ,ims,ime, jms,jme, kms,kme                      &
 155                 ,its,ite, jts,jte, kts,kte                      &
 156                 ,rthcuten,rqvcuten,rqccuten,rqicuten            &
 157                 ,rucuten, rvcuten                               &
 158                 ,f_qv    ,f_qc    ,f_qr    ,f_qi    ,f_qs       &
 159                                                                 )
 160 !-------------------------------------------------------------------
 161       implicit none
 162 !-------------------------------------------------------------------
 163 !-- u3d         3d u-velocity interpolated to theta points (m/s)
 164 !-- v3d         3d v-velocity interpolated to theta points (m/s)
 165 !-- th3d        3d potential temperature (k)
 166 !-- t3d         temperature (k)
 167 !-- qv3d        3d water vapor mixing ratio (kg/kg)
 168 !-- qc3d        3d cloud mixing ratio (kg/kg)
 169 !-- qi3d        3d ice mixing ratio (kg/kg)
 170 !-- rho3d       3d air density (kg/m^3)
 171 !-- p8w         3d hydrostatic pressure at full levels (pa)
 172 !-- pcps        3d hydrostatic pressure at half levels (pa)
 173 !-- pi3d        3d exner function (dimensionless)
 174 !-- qvften      3d total advective + PBL moisture tendency (kg kg-1 s-1)
 175 !-- thften      3d total advective + PBL + radiative temperature tendency (k s-1)
 176 !-- rthcuten      theta tendency due to
 177 !                 cumulus scheme precipitation (k/s)
 178 !-- rucuten       u wind tendency due to
 179 !                 cumulus scheme precipitation (k/s)
 180 !-- rvcuten       v wind tendency due to
 181 !                 cumulus scheme precipitation (k/s)
 182 !-- rqvcuten      qv tendency due to
 183 !                 cumulus scheme precipitation (kg/kg/s)
 184 !-- rqccuten      qc tendency due to
 185 !                 cumulus scheme precipitation (kg/kg/s)
 186 !-- rqicuten      qi tendency due to
 187 !                 cumulus scheme precipitation (kg/kg/s)
 188 !-- rainc         accumulated total cumulus scheme precipitation (mm)
 189 !-- raincv        cumulus scheme precipitation (mm)
 190 !-- pratec        precipitiation rate from cumulus scheme (mm/s)
 191 !-- dz8w        dz between full levels (m)
 192 !-- qfx         upward moisture flux at the surface (kg/m^2/s)
 193 !-- hfx         upward heat flux at the surface (w/m^2)
 194 !-- dt          time step (s)
 195 !-- ids         start index for i in domain
 196 !-- ide         end index for i in domain
 197 !-- jds         start index for j in domain
 198 !-- jde         end index for j in domain
 199 !-- kds         start index for k in domain
 200 !-- kde         end index for k in domain
 201 !-- ims         start index for i in memory
 202 !-- ime         end index for i in memory
 203 !-- jms         start index for j in memory
 204 !-- jme         end index for j in memory
 205 !-- kms         start index for k in memory
 206 !-- kme         end index for k in memory
 207 !-- its         start index for i in tile
 208 !-- ite         end index for i in tile
 209 !-- jts         start index for j in tile
 210 !-- jte         end index for j in tile
 211 !-- kts         start index for k in tile
 212 !-- kte         end index for k in tile
 213 !-------------------------------------------------------------------
 214       integer, intent(in) ::            ids,ide, jds,jde, kds,kde,      &
 215                                         ims,ime, jms,jme, kms,kme,      &
 216                                         its,ite, jts,jte, kts,kte,      &
 217                                         itimestep,                      &
 218                                         stepcu
 219
 220       real,    intent(in) ::                                            &
 221                                         dt
 222       real,    dimension(ims:ime, jms:jme), intent(in) ::               &
 223                                         dx
 224
 225       real,    dimension(ims:ime, jms:jme), intent(in) ::               &
 226                                         xland
 227
 228       real,    dimension(ims:ime, jms:jme), intent(inout) ::            &
 229                                         raincv, pratec
 230
 231       logical, dimension(ims:ime,jms:jme), intent(inout) ::             &
 232                                         cu_act_flag
 233
 234       real,    dimension(ims:ime, kms:kme, jms:jme), intent(in) ::      &
 235                                         dz8w,                           &
 236                                         pcps,                           &
 237                                         p8w,                            &
 238                                         pi3d,                           &
 239                                         qc3d,                           &
 240                                         qvften,                         &
 241                                         thften,                         &
 242                                         qi3d,                           &
 243                                         qv3d,                           &
 244                                         rho3d,                          &
 245                                         t3d,                            &
 246                                         u3d,                            &
 247                                         v3d,                            &
 248                                         w
 249       real,    dimension(ims:ime, jms:jme) ::                           &
 250                                         qfx,                            &
 251                                         hfx
 252
 253 !--------------------------- optional vars ----------------------------
 254
 255       real, dimension(ims:ime, kms:kme, jms:jme),                       &
 256                optional, intent(inout) ::                               &
 257                                         rqccuten,                       &
 258                                         rqicuten,                       &
 259                                         rqvcuten,                       &
 260                                         rthcuten,                       &
 261                                         rucuten,                        &
 262                                         rvcuten
 263
 264 !
 265 ! flags relating to the optional tendency arrays declared above
 266 ! models that carry the optional tendencies will provdide the
 267 ! optional arguments at compile time; these flags all the model
 268 ! to determine at run-time whether a particular tracer is in
 269 ! use or not.
 270 !
 271      logical, optional ::                                      &
 272                                                    f_qv      &
 273                                                   ,f_qc      &
 274                                                   ,f_qr      &
 275                                                   ,f_qi      &
 276                                                   ,f_qs
 277
 278 !--------------------------- local vars ------------------------------
 279       real      ::                                      &
 280                                         delt,                           &
 281                                         rdelt
 282
 283       real     , dimension(its:ite) ::                  &
 284                                         rcs,                            &
 285                                         rn,                             &
 286                                         evap,                           &
 287                                         heatflux,                       &
 288                                         dx2d
 289
 290       integer  , dimension(its:ite) ::  slimsk
 291
 292
 293       real     , dimension(its:ite, kts:kte+1) ::         &
 294                                         prsi,           &
 295                                         ghti,           &
 296                                           zi
 297
 298       real     , dimension(its:ite, kts:kte) ::         &
 299                                         dot,                            &
 300                                         prsl,                           &
 301                                         q1,                             &
 302                                         q2,                             &
 303                                         q3,                             &
 304                                         q1b,                            &
 305                                         t1b,                            &
 306                                         q11,                            &
 307                                         q12,                            &
 308                                         t1,                             &
 309                                         u1,                             &
 310                                         v1,                             &
 311                                         zl,                             &
 312                                         omg,                            &
 313                                         ghtl
 314
 315       integer, dimension(its:ite) ::                                    &
 316                                         kbot,                           &
 317                                         ktop
 318
 319       integer ::                                                        &
 320                                         i,                              &
 321                                         im,                             &
 322                                         j,                              &
 323                                         k,                              &
 324                                         km,                             &
 325                                         kp,                             &
 326                                         kx,                             &
 327                                         kx1
 328
 329 !-------other local variables----
 330       integer                      :: zz, pp
 331 !-----------------------------------------------------------------------
 332 !
 333 !
 334 !***  check to see if this is a convection timestep
 335 !
 336
 337 !-----------------------------------------------------------------------
 338       do j=jts,jte
 339          do i=its,ite
 340             cu_act_flag(i,j)=.true.
 341          enddo
 342       enddo
 343
 344       im=ite-its+1
 345       kx=kte-kts+1
 346       kx1=kx+1
 347       delt=dt*stepcu
 348       rdelt=1./delt
 349
 350 !-------------  j loop (outer) --------------------------------------------------
 351
 352    do j=jts,jte
 353
 354 ! --------------- compute zi and zl -----------------------------------------
 355       do i=its,ite
 356         zi(i,kts)=0.0
 357       enddo
 358 !
 359       do k=kts,kte
 360         do i=its,ite
 361           zi(i,k+1)=zi(i,k)+dz8w(i,k,j)
 362         enddo
 363       enddo
 364 !
 365       do k=kts,kte
 366         do i=its,ite
 367           zl(i,k)=0.5*(zi(i,k)+zi(i,k+1))
 368         enddo
 369       enddo
 370
 371 ! --------------- end compute zi and zl -------------------------------------
 372       do i=its,ite
 373         slimsk(i)=int(abs(xland(i,j)-2.))
 374       enddo
 375
 376       do i=its,ite
 377          dx2d(i) = dx(i,j)
 378       enddo
 379
 380       do k=kts,kte
 381         kp=k+1
 382         do i=its,ite
 383           dot(i,k)=-0.5*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
 384         enddo
 385       enddo
 386
 387       pp = 0
 388       do k=kts,kte
 389         zz = kte-pp
 390         do i=its,ite
 391           u1(i,zz)=u3d(i,k,j)
 392           v1(i,zz)=v3d(i,k,j)
 393           t1(i,zz)=t3d(i,k,j)
 394           q1(i,zz)=qv3d(i,k,j)
 395           if(itimestep == 1) then
 396              q1b(i,zz)=0.
 397              t1b(i,zz)=0.
 398           else
 399              q1b(i,zz)=qvften(i,k,j)
 400              t1b(i,zz)=thften(i,k,j)
 401           endif
 402           q2(i,zz)=qc3d(i,k,j)
 403           q3(i,zz)=qi3d(i,k,j)
 404           omg(i,zz)=dot(i,k)
 405           ghtl(i,zz)=zl(i,k)
 406           prsl(i,zz) = pcps(i,k,j)
 407         enddo
 408         pp = pp + 1
 409       enddo
 410
 411       pp = 0
 412       do k=kts,kte+1
 413         zz = kte+1-pp
 414         do i=its,ite
 415           ghti(i,zz) = zi(i,k)
 416           prsi(i,zz) = p8w(i,k,j)
 417         enddo
 418         pp = pp + 1
 419       enddo
 420 !
 421       do i=its,ite
 422         evap(i)    = qfx(i,j)
 423         heatflux(i)= hfx(i,j)
 424       enddo
 425 !
 426 !########################################################################
 427       call tiecnvn(u1,v1,t1,q1,q2,q3,q1b,t1b,ghtl,ghti,omg,prsl,prsi,evap,heatflux,  &
 428                   rn,slimsk,im,kx,kx1,delt,dx2d)
 429
 430       do i=its,ite
 431          raincv(i,j)=rn(i)/stepcu
 432          pratec(i,j)=rn(i)/(stepcu * dt)
 433       enddo
 434
 435       pp = 0
 436       do k=kts,kte
 437         zz = kte-pp
 438         do i=its,ite
 439           rthcuten(i,k,j)=(t1(i,zz)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
 440           rqvcuten(i,k,j)=(q1(i,zz)-qv3d(i,k,j))*rdelt
 441           rucuten(i,k,j) =(u1(i,zz)-u3d(i,k,j))*rdelt
 442           rvcuten(i,k,j) =(v1(i,zz)-v3d(i,k,j))*rdelt
 443         enddo
 444         pp = pp + 1
 445       enddo
 446
 447       if(present(rqccuten))then
 448         if ( f_qc ) then
 449           pp = 0
 450           do k=kts,kte
 451             zz = kte-pp
 452             do i=its,ite
 453               rqccuten(i,k,j)=(q2(i,zz)-qc3d(i,k,j))*rdelt
 454             enddo
 455             pp = pp + 1
 456           enddo
 457         endif
 458       endif
 459
 460       if(present(rqicuten))then
 461         if ( f_qi ) then
 462           pp = 0
 463           do k=kts,kte
 464             zz = kte-pp
 465             do i=its,ite
 466               rqicuten(i,k,j)=(q3(i,zz)-qi3d(i,k,j))*rdelt
 467             enddo
 468             pp = pp + 1
 469           enddo
 470         endif
 471       endif
 472
 473
 474    enddo
 475
 476    end subroutine cu_ntiedtke
 477
 478 !====================================================================
 479    subroutine ntiedtkeinit(rthcuten,rqvcuten,rqccuten,rqicuten,         &
 480                      rucuten,rvcuten,rthften,rqvften,                   &
 481                      restart,p_qc,p_qi,p_first_scalar,                  &
 482                      allowed_to_read,                                   &
 483                      ids, ide, jds, jde, kds, kde,                      &
 484                      ims, ime, jms, jme, kms, kme,                      &
 485                      its, ite, jts, jte, kts, kte)
 486 !--------------------------------------------------------------------
 487    implicit none
 488 !--------------------------------------------------------------------
 489    logical , intent(in)           ::  allowed_to_read,restart
 490    integer , intent(in)           ::  ids, ide, jds, jde, kds, kde, &
 491                                       ims, ime, jms, jme, kms, kme, &
 492                                       its, ite, jts, jte, kts, kte
 493    integer , intent(in)           ::  p_first_scalar, p_qi, p_qc
 494
 495    real,     dimension( ims:ime , kms:kme , jms:jme ) , intent(out) ::  &
 496                                                               rthcuten, &
 497                                                               rqvcuten, &
 498                                                               rqccuten, &
 499                                                               rqicuten, &
 500                                                               rucuten,rvcuten,&
 501                                                               rthften,rqvften
 502
 503    integer :: i, j, k, itf, jtf, ktf
 504
 505    jtf=min0(jte,jde-1)
 506    ktf=min0(kte,kde-1)
 507    itf=min0(ite,ide-1)
 508
 509    if(.not.restart)then
 510      do j=jts,jtf
 511      do k=kts,ktf
 512      do i=its,itf
 513        rthcuten(i,k,j)=0.
 514        rqvcuten(i,k,j)=0.
 515        rucuten(i,k,j)=0.
 516        rvcuten(i,k,j)=0.
 517      enddo
 518      enddo
 519      enddo
 520
 521      DO j=jts,jtf
 522      DO k=kts,ktf
 523      DO i=its,itf
 524         rthften(i,k,j)=0.
 525         rqvften(i,k,j)=0.
 526      ENDDO
 527      ENDDO
 528      ENDDO
 529
 530      if (p_qc .ge. p_first_scalar) then
 531         do j=jts,jtf
 532         do k=kts,ktf
 533         do i=its,itf
 534            rqccuten(i,k,j)=0.
 535         enddo
 536         enddo
 537         enddo
 538      endif
 539
 540      if (p_qi .ge. p_first_scalar) then
 541         do j=jts,jtf
 542         do k=kts,ktf
 543         do i=its,itf
 544            rqicuten(i,k,j)=0.
 545         enddo
 546         enddo
 547         enddo
 548      endif
 549    endif
 550
 551    end subroutine ntiedtkeinit
 552
 553 !-----------------------------------------------------------------
 554 !          level 1 subroutine 'tiecnvn'
 555 !-----------------------------------------------------------------
 556       subroutine tiecnvn(pu,pv,pt,pqv,pqc,pqi,pqvf,ptf,poz,pzz,pomg, &
 557      &         pap,paph,evap,hfx,zprecc,lndj,lq,km,km1,dt,dx)
 558 !-----------------------------------------------------------------
 559 !  this is the interface between the model and the mass
 560 !  flux convection module
 561 !-----------------------------------------------------------------
 562       implicit none
 563 !
 564       real pu(lq,km),   pv(lq,km),   pt(lq,km),   pqv(lq,km)
 565       real poz(lq,km),  pomg(lq,km), evap(lq),    zprecc(lq)
 566       real pzz(lq,km1)
 567
 568       real pum1(lq,km),    pvm1(lq,km),  ztt(lq,km),                &
 569      &     ptte(lq,km),    pqte(lq,km),  pvom(lq,km),  pvol(lq,km), &
 570      &     pverv(lq,km),   pgeo(lq,km),  pap(lq,km),   paph(lq,km1)
 571       real pqhfl(lq),      zqq(lq,km),    &
 572      &     prsfc(lq),      pssfc(lq),    pcte(lq,km), &
 573      &     phhfl(lq),      hfx(lq),      pgeoh(lq,km1)
 574       real ztp1(lq,km),    zqp1(lq,km),  ztu(lq,km),   zqu(lq,km),  &
 575      &     zlu(lq,km),     zlude(lq,km), zmfu(lq,km),  zmfd(lq,km), &
 576      &     zqsat(lq,km),   pqc(lq,km),   pqi(lq,km),   zrain(lq)
 577       real pqvf(lq,km),    ptf(lq,km)
 578       real dx(lq)
 579
 580       integer icbot(lq),   ictop(lq),     ktype(lq),   lndj(lq)
 581       logical locum(lq)
 582 !
 583       real ztmst,fliq,fice,ztc,zalf,tt
 584       integer i,j,k,lq,km,km1
 585       real dt,ztpp1
 586       real zew,zqs,zcor
 587       real scale_fac(lq), scale_fac2(lq), dxref
 588 !
 589 !  set scale-dependency factor when dx is < 15 km
 590 !
 591       dxref = 15000.
 592       do j=1,lq
 593       if (dx(j).lt.dxref) then
 594           scale_fac(j) = (1.06133+log(dxref/dx(j)))**3
 595           scale_fac2(j) = scale_fac(j)**0.5
 596       else
 597           scale_fac(j) = 1.+1.33e-5*dx(j)
 598           scale_fac2(j) = 1.
 599       end if
 600       end do
 601 !
 602       ztmst=dt
 603 !
 604 !  masv flux diagnostics.
 605 !
 606       do j=1,lq
 607         zrain(j)=0.0
 608         locum(j)=.false.
 609         prsfc(j)=0.0
 610         pssfc(j)=0.0
 611         pqhfl(j)=evap(j)
 612         phhfl(j)=hfx(j)
 613         pgeoh(j,km1)=g*pzz(j,km1)
 614       end do
 615 !
 616 !     convert model variables for mflux scheme
 617 !
 618       do k=1,km
 619         do j=1,lq
 620           pcte(j,k)=0.0
 621           pvom(j,k)=0.0
 622           pvol(j,k)=0.0
 623           ztp1(j,k)=pt(j,k)
 624           zqp1(j,k)=pqv(j,k)/(1.0+pqv(j,k))
 625           pum1(j,k)=pu(j,k)
 626           pvm1(j,k)=pv(j,k)
 627           pverv(j,k)=pomg(j,k)
 628           pgeo(j,k)=g*poz(j,k)
 629           pgeoh(j,k)=g*pzz(j,k)
 630           tt=ztp1(j,k)
 631           zew  = foeewm(tt)
 632           zqs  = zew/pap(j,k)
 633           zqs  = min(0.5,zqs)
 634           zcor = 1./(1.-vtmpc1*zqs)
 635           zqsat(j,k)=zqs*zcor
 636           pqte(j,k)=pqvf(j,k)
 637           zqq(j,k) =pqte(j,k)
 638           ptte(j,k)=ptf(j,k)
 639           ztt(j,k) =ptte(j,k)
 640         end do
 641       end do
 642 !
 643 !-----------------------------------------------------------------------
 644 !*    2.     call 'cumastrn'(master-routine for cumulus parameterization)
 645 !
 646       call cumastrn        &
 647      &    (lq,       km,       km1,      km-1,    ztp1, &
 648      &     zqp1,     pum1,     pvm1,     pverv,   zqsat,&
 649      &     pqhfl,    ztmst,    pap,      paph,    pgeo, &
 650      &     ptte,     pqte,     pvom,     pvol,    prsfc,&
 651      &     pssfc,    locum,                             &
 652      &     ktype,    icbot,    ictop,    ztu,     zqu,  &
 653      &     zlu,      zlude,    zmfu,     zmfd,    zrain,&
 654      &     pcte,     phhfl,    lndj,     pgeoh,   dx,   &
 655      &     scale_fac, scale_fac2)
 656 !
 657 !     to include the cloud water and cloud ice detrained from convection
 658 !
 659       do k=1,km
 660       do j=1,lq
 661       if(pcte(j,k).gt.0.) then
 662         fliq=foealfa(ztp1(j,k))
 663         fice=1.0-fliq
 664         pqc(j,k)=pqc(j,k)+fliq*pcte(j,k)*ztmst
 665         pqi(j,k)=pqi(j,k)+fice*pcte(j,k)*ztmst
 666       endif
 667       end do
 668       end do
 669 !
 670       do k=1,km
 671         do j=1,lq
 672           pt(j,k)=  ztp1(j,k)+(ptte(j,k)-ztt(j,k))*ztmst
 673           zqp1(j,k)=zqp1(j,k)+(pqte(j,k)-zqq(j,k))*ztmst
 674           pqv(j,k)=zqp1(j,k)/(1.0-zqp1(j,k))
 675         end do
 676       end do
 677
 678       do j=1,lq
 679         zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst)
 680       end do
 681
 682       if (lmfdudv) then
 683         do k=1,km
 684           do j=1,lq
 685             pu(j,k)=pu(j,k)+pvom(j,k)*ztmst
 686             pv(j,k)=pv(j,k)+pvol(j,k)*ztmst
 687           end do
 688         end do
 689       endif
 690 !
 691       return
 692       end subroutine tiecnvn
 693
 694 !#############################################################
 695 !
 696 !             level 2 subroutines
 697 !
 698 !#############################################################
 699 !***********************************************************
 700 !           subroutine cumastrn
 701 !***********************************************************
 702       subroutine cumastrn  &
 703      &    (klon,     klev,     klevp1,   klevm1,   pten, &
 704      &     pqen,     puen,     pven,     pverv,    pqsen,&
 705      &     pqhfl,    ztmst,    pap,      paph,     pgeo, &
 706      &     ptte,     pqte,     pvom,     pvol,     prsfc,&
 707      &     pssfc,    ldcum,                              &
 708      &     ktype,    kcbot,    kctop,    ptu,      pqu,&
 709      &     plu,      plude,    pmfu,     pmfd,     prain,&
 710      &     pcte,     phhfl,    lndj,     zgeoh,    dx,   &
 711      &     scale_fac,  scale_fac2)
 712       implicit none
 713 !
 714 !***cumastrn*  master routine for cumulus massflux-scheme
 715 !     m.tiedtke      e.c.m.w.f.      1986/1987/1989
 716 !     modifications
 717 !     y.wang         i.p.r.c         2001
 718 !     c.zhang                        2012
 719 !***purpose
 720 !   -------
 721 !          this routine computes the physical tendencies of the
 722 !     prognostic variables t,q,u and v due to convective processes.
 723 !     processes considered are: convective fluxes, formation of
 724 !     precipitation, evaporation of falling rain below cloud base,
 725 !     saturated cumulus downdrafts.
 726 !***method
 727 !   ------
 728 !     parameterization is done using a massflux-scheme.
 729 !        (1) define constants and parameters
 730 !        (2) specify values (t,q,qs...) at half levels and
 731 !            initialize updraft- and downdraft-values in 'cuinin'
 732 !        (3) calculate cloud base in 'cutypen', calculate cloud types in cutypen,
 733 !            and specify cloud base massflux
 734 !        (4) do cloud ascent in 'cuascn' in absence of downdrafts
 735 !        (5) do downdraft calculations:
 736 !              (a) determine values at lfs in 'cudlfsn'
 737 !              (b) determine moist descent in 'cuddrafn'
 738 !              (c) recalculate cloud base massflux considering the
 739 !                  effect of cu-downdrafts
 740 !        (6) do final adjusments to convective fluxes in 'cuflxn',
 741 !            do evaporation in subcloud layer
 742 !        (7) calculate increments of t and q in 'cudtdqn'
 743 !        (8) calculate increments of u and v in 'cududvn'
 744 !***externals.
 745 !   ----------
 746 !       cuinin:  initializes values at vertical grid used in cu-parametr.
 747 !       cutypen: cloud bypes, 1: deep cumulus 2: shallow cumulus
 748 !       cuascn:  cloud ascent for entraining plume
 749 !       cudlfsn: determines values at lfs for downdrafts
 750 !       cuddrafn:does moist descent for cumulus downdrafts
 751 !       cuflxn:  final adjustments to convective fluxes (also in pbl)
 752 !       cudqdtn: updates tendencies for t and q
 753 !       cududvn: updates tendencies for u and v
 754 !***switches.
 755 !   --------
 756 !          lmfmid=.t.   midlevel convection is switched on
 757 !          lmfdd=.t.    cumulus downdrafts switched on
 758 !          lmfdudv=.t.  cumulus friction switched on
 759 !***
 760 !     model parameters (defined in subroutine cuparam)
 761 !     ------------------------------------------------
 762 !     entrdd     entrainment rate for cumulus downdrafts
 763 !     cmfcmax    maximum massflux value allowed for
 764 !     cmfcmin    minimum massflux value (for safety)
 765 !     cmfdeps    fractional massflux for downdrafts at lfs
 766 !     cprcon     coefficient for conversion from cloud water to rain
 767 !***reference.
 768 !   ----------
 769 !          paper on massflux scheme (tiedtke,1989)
 770 !-----------------------------------------------------------------
 771       integer  klev,klon,klevp1,klevm1
 772       real     pten(klon,klev),        pqen(klon,klev),&
 773      &         puen(klon,klev),        pven(klon,klev),&
 774      &         ptte(klon,klev),        pqte(klon,klev),&
 775      &         pvom(klon,klev),        pvol(klon,klev),&
 776      &         pqsen(klon,klev),       pgeo(klon,klev),&
 777      &         pap(klon,klev),         paph(klon,klevp1),&
 778      &         pverv(klon,klev),       pqhfl(klon),&
 779      &         phhfl(klon)
 780       real     ptu(klon,klev),         pqu(klon,klev),&
 781      &         plu(klon,klev),         plude(klon,klev),&
 782      &         pmfu(klon,klev),        pmfd(klon,klev),&
 783      &         prain(klon),&
 784      &         prsfc(klon),            pssfc(klon)
 785       real     ztenh(klon,klev),       zqenh(klon,klev),&
 786      &         zgeoh(klon,klevp1),     zqsenh(klon,klev),&
 787      &         ztd(klon,klev),         zqd(klon,klev),&
 788      &         zmfus(klon,klev),       zmfds(klon,klev),&
 789      &         zmfuq(klon,klev),       zmfdq(klon,klev),&
 790      &         zdmfup(klon,klev),      zdmfdp(klon,klev),&
 791      &         zmful(klon,klev),       zrfl(klon),&
 792      &         zuu(klon,klev),         zvu(klon,klev),&
 793      &         zud(klon,klev),         zvd(klon,klev),&
 794      &         zlglac(klon,klev)
 795       real     pmflxr(klon,klevp1),    pmflxs(klon,klevp1)
 796       real     zhcbase(klon),&
 797      &         zmfub(klon),            zmfub1(klon),&
 798      &         zdhpbl(klon)
 799       real     zsfl(klon),             zdpmel(klon,klev),&
 800      &         pcte(klon,klev),        zcape(klon),&
 801      &         zcape1(klon),           zcape2(klon),&
 802      &         ztauc(klon),            ztaubl(klon),&
 803      &         zheat(klon)
 804       real     wup(klon),              zdqcv(klon)
 805       real     wbase(klon),            zmfuub(klon)
 806       real     upbl(klon)
 807       real     dx(klon)
 808       real     pmfude_rate(klon,klev), pmfdde_rate(klon,klev)
 809       real     zmfuus(klon,klev),      zmfdus(klon,klev)
 810       real     zuv2(klon,klev),ztenu(klon,klev),ztenv(klon,klev)
 811       real     zmfuvb(klon),zsum12(klon),zsum22(klon)
 812       integer  ilab(klon,klev),        idtop(klon),&
 813      &         ictop0(klon),           ilwmin(klon)
 814       integer  kdpl(klon)
 815       integer  kcbot(klon),            kctop(klon),&
 816      &         ktype(klon),            lndj(klon)
 817       logical  ldcum(klon)
 818       logical  loddraf(klon),          llo1,   llo2(klon)
 819       real     scale_fac(klon), scale_fac2(klon)
 820
 821 !  local varaiables
 822       real     zcons,zcons2,zqumqe,zdqmin,zdh,zmfmax
 823       real     zalfaw,zalv,zqalv,zc5ldcp,zc4les,zhsat,zgam,zzz,zhhat
 824       real     zpbmpt,zro,zdz,zdp,zeps,zfac,wspeed
 825       integer  jl,jk,ik
 826       integer  ikb,ikt,icum,itopm2
 827       real     ztmst,ztau,zerate,zderate,zmfa
 828       real     zmfs(klon),pmean(klev),zlon
 829       real     zduten,zdvten,ztdis,pgf_u,pgf_v
 830 !-------------------------------------------
 831 !     1.    specify constants and parameters
 832 !-------------------------------------------
 833       zcons=1./(g*ztmst)
 834       zcons2=3./(g*ztmst)
 835
 836 !--------------------------------------------------------------
 837 !*    2.    initialize values at vertical grid points in 'cuini'
 838 !--------------------------------------------------------------
 839       call cuinin &
 840      &    (klon,     klev,     klevp1,   klevm1,   pten, &
 841      &     pqen,     pqsen,    puen,     pven,     pverv,&
 842      &     pgeo,     paph,     zgeoh,    ztenh,    zqenh,&
 843      &     zqsenh,   ilwmin,   ptu,      pqu,      ztd,  &
 844      &     zqd,      zuu,      zvu,      zud,      zvd,  &
 845      &     pmfu,     pmfd,     zmfus,    zmfds,    zmfuq,&
 846      &     zmfdq,    zdmfup,   zdmfdp,   zdpmel,   plu,  &
 847      &     plude,    ilab)
 848
 849 !----------------------------------
 850 !*    3.0   cloud base calculations
 851 !----------------------------------
 852 !*         (a) determine cloud base values in 'cutypen',
 853 !              and the cumulus type 1 or 2
 854 !          -------------------------------------------
 855        call cutypen &
 856      &     ( klon,     klev,     klevp1,   klevm1,     pqen,&
 857      &      ztenh,    zqenh,     zqsenh,    zgeoh,     paph,&
 858      &      phhfl,    pqhfl,       pgeo,    pqsen,      pap,&
 859      &       pten,     lndj,        ptu,      pqu,     ilab,&
 860      &      ldcum,    kcbot,     ictop0,    ktype,    wbase,    plu,   kdpl)
 861
 862 !*         (b) assign the first guess mass flux at cloud base
 863 !              ------------------------------------------
 864        do jl=1,klon
 865          zdhpbl(jl)=0.0
 866          upbl(jl) = 0.0
 867          idtop(jl)=0
 868        end do
 869
 870        do jk=2,klev
 871        do jl=1,klon
 872          if(jk.ge.kcbot(jl) .and. ldcum(jl)) then
 873             zdhpbl(jl)=zdhpbl(jl)+(alv*pqte(jl,jk)+cpd*ptte(jl,jk))&
 874      &                 *(paph(jl,jk+1)-paph(jl,jk))
 875             if(lndj(jl) .eq. 0) then
 876               wspeed = sqrt(puen(jl,jk)**2 + pven(jl,jk)**2)
 877               upbl(jl) = upbl(jl) + wspeed*(paph(jl,jk+1)-paph(jl,jk))
 878             end if
 879          end if
 880        end do
 881        end do
 882
 883       do jl=1,klon
 884         if(ldcum(jl)) then
 885            ikb=kcbot(jl)
 886            zmfmax = (paph(jl,ikb)-paph(jl,ikb-1))*zcons2
 887            if(ktype(jl) == 1) then
 888              zmfub(jl)= 0.1*zmfmax
 889            else if ( ktype(jl) == 2 ) then
 890              zqumqe = pqu(jl,ikb) + plu(jl,ikb) - zqenh(jl,ikb)
 891              zdqmin = max(0.01*zqenh(jl,ikb),1.e-10)
 892              zdh = cpd*(ptu(jl,ikb)-ztenh(jl,ikb)) + alv*zqumqe
 893              zdh = g*max(zdh,1.e5*zdqmin)
 894              if ( zdhpbl(jl) > 0. ) then
 895                zmfub(jl) = zdhpbl(jl)/zdh
 896                zmfub(jl) = min(zmfub(jl),zmfmax)
 897              else
 898                zmfub(jl) = 0.1*zmfmax
 899                ldcum(jl) = .false.
 900              end if
 901             end if
 902         else
 903            zmfub(jl) = 0.
 904         end if
 905       end do
 906 !------------------------------------------------------
 907 !*    4.0   determine cloud ascent for entraining plume
 908 !------------------------------------------------------
 909 !*    (a) do ascent in 'cuasc'in absence of downdrafts
 910 !----------------------------------------------------------
 911       call cuascn &
 912      &    (klon,     klev,     klevp1,   klevm1,   ztenh,&
 913      &     zqenh,    puen,     pven,     pten,     pqen,&
 914      &     pqsen,    pgeo,     zgeoh,    pap,      paph,&
 915      &     pqte,     pverv,    ilwmin,   ldcum,    zhcbase,&
 916      &     ktype,    ilab,     ptu,      pqu,      plu,&
 917      &     zuu,      zvu,      pmfu,     zmfub,&
 918      &     zmfus,    zmfuq,    zmful,    plude,    zdmfup,&
 919      &     kcbot,    kctop,    ictop0,   icum,     ztmst,&
 920      &     zqsenh,   zlglac,   lndj,     wup,      wbase,   kdpl,  pmfude_rate)
 921
 922 !*     (b) check cloud depth and change entrainment rate accordingly
 923 !          calculate precipitation rate (for downdraft calculation)
 924 !------------------------------------------------------------------
 925       do jl=1,klon
 926         if ( ldcum(jl) ) then
 927           ikb = kcbot(jl)
 928           itopm2 = kctop(jl)
 929           zpbmpt = paph(jl,ikb) - paph(jl,itopm2)
 930           if ( ktype(jl) == 1 .and. zpbmpt <  zdnoprc ) ktype(jl) = 2
 931           if ( ktype(jl) == 2 .and. zpbmpt >= zdnoprc ) ktype(jl) = 1
 932           ictop0(jl) = kctop(jl)
 933         end if
 934         zrfl(jl)=zdmfup(jl,1)
 935       end do
 936
 937       do jk=2,klev
 938         do jl=1,klon
 939           zrfl(jl)=zrfl(jl)+zdmfup(jl,jk)
 940         end do
 941       end do
 942
 943       do jk = 1,klev
 944       do jl = 1,klon
 945         pmfd(jl,jk) = 0.
 946         zmfds(jl,jk) = 0.
 947         zmfdq(jl,jk) = 0.
 948         zdmfdp(jl,jk) = 0.
 949         zdpmel(jl,jk) = 0.
 950       end do
 951       end do
 952
 953 !-----------------------------------------
 954 !*    6.0   cumulus downdraft calculations
 955 !-----------------------------------------
 956       if(lmfdd) then
 957 !*      (a) determine lfs in 'cudlfsn'
 958 !--------------------------------------
 959         call cudlfsn    &
 960      &    (klon,     klev,&
 961      &     kcbot,    kctop,    lndj,   ldcum,   &
 962      &     ztenh,    zqenh,    puen,     pven,   &
 963      &     pten,     pqsen,    pgeo,              &
 964      &     zgeoh,    paph,     ptu,      pqu,      plu,   &
 965      &     zuu,      zvu,      zmfub,    zrfl,             &
 966      &     ztd,      zqd,      zud,      zvd,               &
 967      &     pmfd,     zmfds,    zmfdq,    zdmfdp,             &
 968      &     idtop,    loddraf)
 969 !*     (b)  determine downdraft t,q and fluxes in 'cuddrafn'
 970 !------------------------------------------------------------
 971         call cuddrafn                                         &
 972      &   ( klon,     klev,     loddraf,                      &
 973      &     ztenh,    zqenh,    puen,     pven,                 &
 974      &     pgeo,     zgeoh,    paph,     zrfl,                  &
 975      &     ztd,      zqd,      zud,      zvd,      pmfu,         &
 976      &     pmfd,     zmfds,    zmfdq,    zdmfdp,   pmfdde_rate      )
 977 !-----------------------------------------------------------
 978       end if
 979 !
 980 !-----------------------------------------------------------------------
 981 !* 6.0          closure and clean work
 982 ! ------
 983 !-- 6.1 recalculate cloud base massflux from a cape closure
 984 !       for deep convection (ktype=1)
 985 !
 986       do jl=1,klon
 987       if(ldcum(jl) .and. ktype(jl) .eq. 1) then
 988         ikb = kcbot(jl)
 989         ikt = kctop(jl)
 990         zheat(jl)=0.0
 991         zcape(jl)=0.0
 992         zcape1(jl)=0.0
 993         zcape2(jl)=0.0
 994         zmfub1(jl)=zmfub(jl)
 995
 996         ztauc(jl)  = (zgeoh(jl,ikt)-zgeoh(jl,ikb)) / &
 997                    ((2.+ min(15.0,wup(jl)))*g)
 998         if(lndj(jl) .eq. 0) then
 999           upbl(jl) = 2.+ upbl(jl)/(paph(jl,klev+1)-paph(jl,ikb))
1000           ztaubl(jl) = (zgeoh(jl,ikb)-zgeoh(jl,klev+1))/(g*upbl(jl))
1001           ztaubl(jl) = min(300., ztaubl(jl))
1002         else
1003           ztaubl(jl) = ztauc(jl)
1004         end if
1005       end if
1006       end do
1007 !
1008       do jk = 1 , klev
1009       do jl = 1 , klon
1010         llo1 = ldcum(jl) .and. ktype(jl) .eq. 1
1011         if ( llo1 .and. jk <= kcbot(jl) .and. jk > kctop(jl) ) then
1012           ikb = kcbot(jl)
1013           zdz = pgeo(jl,jk-1)-pgeo(jl,jk)
1014           zdp = pap(jl,jk)-pap(jl,jk-1)
1015           zheat(jl) = zheat(jl) + ((pten(jl,jk-1)-pten(jl,jk)+zdz*rcpd) / &
1016                       ztenh(jl,jk)+vtmpc1*(pqen(jl,jk-1)-pqen(jl,jk))) * &
1017                       (g*(pmfu(jl,jk)+pmfd(jl,jk)))
1018           zcape1(jl) = zcape1(jl) + ((ptu(jl,jk)-ztenh(jl,jk))/ztenh(jl,jk) + &
1019                       vtmpc1*(pqu(jl,jk)-zqenh(jl,jk))-plu(jl,jk))*zdp
1020         end if
1021
1022         if ( llo1 .and. jk >= kcbot(jl) ) then
1023         if((paph(jl,klev+1)-paph(jl,kdpl(jl)))<50.e2) then
1024           zdp = paph(jl,jk+1)-paph(jl,jk)
1025           zcape2(jl) = zcape2(jl) + ztaubl(jl)* &
1026                      ((1.+vtmpc1*pqen(jl,jk))*ptte(jl,jk)+vtmpc1*pten(jl,jk)*pqte(jl,jk))*zdp
1027         end if
1028         end if
1029       end do
1030       end do
1031
1032       do jl=1,klon
1033        if(ldcum(jl).and.ktype(jl).eq.1) then
1034            ikb = kcbot(jl)
1035            ikt = kctop(jl)
1036            ztauc(jl) = max(ztmst,ztauc(jl))
1037            ztauc(jl) = max(360.,ztauc(jl))
1038            ztauc(jl) = min(10800.,ztauc(jl))
1039            ztau = ztauc(jl) * scale_fac(jl)
1040            if(nonequil) then
1041              zcape2(jl)= max(0.,zcape2(jl))
1042              zcape(jl) = max(0.,min(zcape1(jl)-zcape2(jl),5000.))
1043            else
1044              zcape(jl) = max(0.,min(zcape1(jl),5000.))
1045            end if
1046            zheat(jl) = max(1.e-4,zheat(jl))
1047            zmfub1(jl) = (zcape(jl)*zmfub(jl))/(zheat(jl)*ztau)
1048            zmfub1(jl) = max(zmfub1(jl),0.001)
1049            zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
1050            zmfub1(jl)=min(zmfub1(jl),zmfmax)
1051        end if
1052       end do
1053 !
1054 !*  6.2   recalculate convective fluxes due to effect of
1055 !         downdrafts on boundary layer moist static energy budget (ktype=2)
1056 !--------------------------------------------------------
1057        do jl=1,klon
1058          if(ldcum(jl) .and. ktype(jl) .eq. 2) then
1059            ikb=kcbot(jl)
1060            if(pmfd(jl,ikb).lt.0.0 .and. loddraf(jl)) then
1061               zeps=-pmfd(jl,ikb)/max(zmfub(jl),cmfcmin)
1062            else
1063               zeps=0.
1064            endif
1065            zqumqe=pqu(jl,ikb)+plu(jl,ikb)-  &
1066      &            zeps*zqd(jl,ikb)-(1.-zeps)*zqenh(jl,ikb)
1067            zdqmin=max(0.01*zqenh(jl,ikb),cmfcmin)
1068            zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2
1069 !  using moist static engergy closure instead of moisture closure
1070            zdh=cpd*(ptu(jl,ikb)-zeps*ztd(jl,ikb)- &
1071      &       (1.-zeps)*ztenh(jl,ikb))+alv*zqumqe
1072            zdh=g*max(zdh,1.e5*zdqmin)
1073            if(zdhpbl(jl).gt.0.)then
1074              zmfub1(jl)=zdhpbl(jl)/zdh
1075            else
1076              zmfub1(jl) = zmfub(jl)
1077            end if
1078            zmfub1(jl) = zmfub1(jl)/scale_fac2(jl)
1079            zmfub1(jl) = min(zmfub1(jl),zmfmax)
1080          end if
1081
1082 !*  6.3   mid-level convection - nothing special
1083 !---------------------------------------------------------
1084          if(ldcum(jl) .and. ktype(jl) .eq. 3 ) then
1085             zmfub1(jl) = zmfub(jl)
1086          end if
1087
1088        end do
1089
1090 !*  6.4   scaling the downdraft mass flux
1091 !---------------------------------------------------------
1092        do jk=1,klev
1093        do jl=1,klon
1094         if( ldcum(jl) ) then
1095            zfac=zmfub1(jl)/max(zmfub(jl),cmfcmin)
1096            pmfd(jl,jk)=pmfd(jl,jk)*zfac
1097            zmfds(jl,jk)=zmfds(jl,jk)*zfac
1098            zmfdq(jl,jk)=zmfdq(jl,jk)*zfac
1099            zdmfdp(jl,jk)=zdmfdp(jl,jk)*zfac
1100            pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zfac
1101         end if
1102        end do
1103        end do
1104
1105 !*  6.5   scaling the updraft mass flux
1106 ! --------------------------------------------------------
1107     do jl = 1,klon
1108       if ( ldcum(jl) ) zmfs(jl) = zmfub1(jl)/max(cmfcmin,zmfub(jl))
1109     end do
1110     do jk = 2 , klev
1111       do jl = 1,klon
1112         if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1113           ikb = kcbot(jl)
1114           if ( jk>ikb ) then
1115             zdz = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
1116             pmfu(jl,jk) = pmfu(jl,ikb)*zdz
1117           end if
1118           zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
1119           if ( pmfu(jl,jk)*zmfs(jl) > zmfmax ) then
1120             zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
1121           end if
1122         end if
1123       end do
1124     end do
1125     do jk = 2 , klev
1126       do jl = 1,klon
1127         if ( ldcum(jl) .and. jk <= kcbot(jl) .and. jk >= kctop(jl)-1 ) then
1128           pmfu(jl,jk) = pmfu(jl,jk)*zmfs(jl)
1129           zmfus(jl,jk) = zmfus(jl,jk)*zmfs(jl)
1130           zmfuq(jl,jk) = zmfuq(jl,jk)*zmfs(jl)
1131           zmful(jl,jk) = zmful(jl,jk)*zmfs(jl)
1132           zdmfup(jl,jk) = zdmfup(jl,jk)*zmfs(jl)
1133           plude(jl,jk) = plude(jl,jk)*zmfs(jl)
1134           pmfude_rate(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl)
1135         end if
1136       end do
1137     end do
1138
1139 !*    6.6  if ktype = 2, kcbot=kctop is not allowed
1140 ! ---------------------------------------------------
1141     do jl = 1,klon
1142       if ( ktype(jl) == 2 .and. &
1143            kcbot(jl) == kctop(jl) .and. kcbot(jl) >= klev-1 ) then
1144         ldcum(jl) = .false.
1145         ktype(jl) = 0
1146       end if
1147     end do
1148
1149     if ( .not. lmfscv .or. .not. lmfpen ) then
1150       do jl = 1,klon
1151         llo2(jl) = .false.
1152         if ( (.not. lmfscv .and. ktype(jl) == 2) .or. &
1153              (.not. lmfpen .and. ktype(jl) == 1) ) then
1154           llo2(jl) = .true.
1155           ldcum(jl) = .false.
1156         end if
1157       end do
1158     end if
1159
1160 !*   6.7  set downdraft mass fluxes to zero above cloud top
1161 !----------------------------------------------------
1162     do jl = 1,klon
1163       if ( loddraf(jl) .and. idtop(jl) <= kctop(jl) ) then
1164         idtop(jl) = kctop(jl) + 1
1165       end if
1166     end do
1167     do jk = 2 , klev
1168       do jl = 1,klon
1169         if ( loddraf(jl) ) then
1170           if ( jk < idtop(jl) ) then
1171             pmfd(jl,jk) = 0.
1172             zmfds(jl,jk) = 0.
1173             zmfdq(jl,jk) = 0.
1174             pmfdde_rate(jl,jk) = 0.
1175             zdmfdp(jl,jk) = 0.
1176           else if ( jk == idtop(jl) ) then
1177             pmfdde_rate(jl,jk) = 0.
1178           end if
1179         end if
1180       end do
1181     end do
1182 !----------------------------------------------------------
1183 !*    7.0      determine final convective fluxes in 'cuflx'
1184 !----------------------------------------------------------
1185        call cuflxn                                        &
1186      &  (  klon,     klev,     ztmst                       &
1187      &  ,  pten,     pqen,     pqsen,    ztenh,   zqenh     &
1188      &  ,  paph,     pap,      zgeoh,    lndj,    ldcum      &
1189      &  ,  kcbot,    kctop,    idtop,    itopm2               &
1190      &  ,  ktype,    loddraf                                   &
1191      &  ,  pmfu,     pmfd,     zmfus,    zmfds                  &
1192      &  ,  zmfuq,    zmfdq,    zmful,    plude                   &
1193      &  ,  zdmfup,   zdmfdp,   zdpmel,   zlglac                   &
1194      &  ,  prain,    pmfdde_rate, pmflxr, pmflxs )
1195
1196 ! some adjustments needed
1197     do jl=1,klon
1198       zmfs(jl) = 1.
1199       zmfuub(jl)=0.
1200     end do
1201     do jk = 2 , klev
1202       do jl = 1,klon
1203         if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
1204           zmfmax = pmfu(jl,jk)*0.98
1205           if ( pmfd(jl,jk)+zmfmax+1.e-15 < 0. ) then
1206             zmfs(jl) = min(zmfs(jl),-zmfmax/pmfd(jl,jk))
1207           end if
1208         end if
1209       end do
1210     end do
1211
1212     do jk = 2 , klev
1213       do jl = 1 , klon
1214         if ( zmfs(jl) < 1. .and. jk >= idtop(jl)-1 ) then
1215           pmfd(jl,jk) = pmfd(jl,jk)*zmfs(jl)
1216           zmfds(jl,jk) = zmfds(jl,jk)*zmfs(jl)
1217           zmfdq(jl,jk) = zmfdq(jl,jk)*zmfs(jl)
1218           pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl)
1219           zmfuub(jl) = zmfuub(jl) - (1.-zmfs(jl))*zdmfdp(jl,jk)
1220           pmflxr(jl,jk+1) = pmflxr(jl,jk+1) + zmfuub(jl)
1221           zdmfdp(jl,jk) = zdmfdp(jl,jk)*zmfs(jl)
1222         end if
1223       end do
1224     end do
1225
1226     do jk = 2 , klev - 1
1227       do jl = 1, klon
1228         if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then
1229           zerate = -pmfd(jl,jk) + pmfd(jl,jk-1) + pmfdde_rate(jl,jk)
1230           if ( zerate < 0. ) then
1231             pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk) - zerate
1232           end if
1233         end if
1234         if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1235           zerate = pmfu(jl,jk) - pmfu(jl,jk+1) + pmfude_rate(jl,jk)
1236           if ( zerate < 0. ) then
1237             pmfude_rate(jl,jk) = pmfude_rate(jl,jk) - zerate
1238           end if
1239           zdmfup(jl,jk) = pmflxr(jl,jk+1) + pmflxs(jl,jk+1) - &
1240                           pmflxr(jl,jk) - pmflxs(jl,jk)
1241           zdmfdp(jl,jk) = 0.
1242         end if
1243       end do
1244     end do
1245
1246 ! avoid negative humidities at ddraught top
1247     do jl = 1,klon
1248       if ( loddraf(jl) ) then
1249         jk = idtop(jl)
1250         ik = min(jk+1,klev)
1251         if ( zmfdq(jl,jk) < 0.3*zmfdq(jl,ik) ) then
1252             zmfdq(jl,jk) = 0.3*zmfdq(jl,ik)
1253         end if
1254       end if
1255     end do
1256
1257 ! avoid negative humidities near cloud top because gradient of precip flux
1258 ! and detrainment / liquid water flux are too large
1259     do jk = 2 , klev
1260       do jl = 1, klon
1261         if ( ldcum(jl) .and. jk >= kctop(jl)-1 .and. jk < kcbot(jl) ) then
1262           zdz = ztmst*g/(paph(jl,jk+1)-paph(jl,jk))
1263           zmfa = zmfuq(jl,jk+1) + zmfdq(jl,jk+1) - &
1264                  zmfuq(jl,jk) - zmfdq(jl,jk) + &
1265                  zmful(jl,jk+1) - zmful(jl,jk) + zdmfup(jl,jk)
1266           zmfa = (zmfa-plude(jl,jk))*zdz
1267           if ( pqen(jl,jk)+zmfa < 0. ) then
1268             plude(jl,jk) = plude(jl,jk) + 2.*(pqen(jl,jk)+zmfa)/zdz
1269           end if
1270           if ( plude(jl,jk) < 0. ) plude(jl,jk) = 0.
1271         end if
1272         if ( .not. ldcum(jl) ) pmfude_rate(jl,jk) = 0.
1273         if ( abs(pmfd(jl,jk-1)) < 1.0e-20 ) pmfdde_rate(jl,jk) = 0.
1274       end do
1275     end do
1276
1277     do jl=1,klon
1278       prsfc(jl) = pmflxr(jl,klev+1)
1279       pssfc(jl) = pmflxs(jl,klev+1)
1280     end do
1281
1282 !----------------------------------------------------------------
1283 !*    8.0      update tendencies for t and q in subroutine cudtdq
1284 !----------------------------------------------------------------
1285       call cudtdqn(klon,klev,itopm2,kctop,idtop,ldcum,loddraf, &
1286                  ztmst,paph,zgeoh,pgeo,pten,ztenh,pqen,zqenh,pqsen,     &
1287                  zlglac,plude,pmfu,pmfd,zmfus,zmfds,zmfuq,zmfdq,zmful,   &
1288                  zdmfup,zdmfdp,zdpmel,ptte,pqte,pcte)
1289 !----------------------------------------------------------------
1290 !*    9.0      update tendencies for u and u in subroutine cududv
1291 !----------------------------------------------------------------
1292       if(lmfdudv) then
1293       do jk = klev-1 , 2 , -1
1294         ik = jk + 1
1295         do jl = 1,klon
1296           if ( ldcum(jl) ) then
1297             if ( jk == kcbot(jl) .and. ktype(jl) < 3 ) then
1298               ikb = kdpl(jl)
1299               zuu(jl,jk) = puen(jl,ikb-1)
1300               zvu(jl,jk) = pven(jl,ikb-1)
1301             else if ( jk == kcbot(jl) .and. ktype(jl) == 3 ) then
1302               zuu(jl,jk) = puen(jl,jk-1)
1303               zvu(jl,jk) = pven(jl,jk-1)
1304             end if
1305             if ( jk < kcbot(jl) .and. jk >= kctop(jl) ) then
1306             if(momtrans .eq. 1)then
1307               zfac = 0.
1308               if ( ktype(jl) == 1 .or. ktype(jl) == 3 ) zfac = 2.
1309               if ( ktype(jl) == 1 .and. jk <= kctop(jl)+2 ) zfac = 3.
1310               zerate = pmfu(jl,jk) - pmfu(jl,ik) + &
1311                 (1.+zfac)*pmfude_rate(jl,jk)
1312               zderate = (1.+zfac)*pmfude_rate(jl,jk)
1313               zmfa = 1./max(cmfcmin,pmfu(jl,jk))
1314               zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
1315                 zerate*puen(jl,jk)-zderate*zuu(jl,ik))*zmfa
1316               zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
1317                 zerate*pven(jl,jk)-zderate*zvu(jl,ik))*zmfa
1318             else
1319               pgf_u = -pgcoef*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+&
1320                                    pmfu(jl,jk)*(puen(jl,jk)-puen(jl,jk-1)))
1321               pgf_v = -pgcoef*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+&
1322                                    pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1)))
1323               zerate = pmfu(jl,jk) - pmfu(jl,ik) + pmfude_rate(jl,jk)
1324               zderate = pmfude_rate(jl,jk)
1325               zmfa = 1./max(cmfcmin,pmfu(jl,jk))
1326               zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + &
1327                 zerate*puen(jl,jk)-zderate*zuu(jl,ik)+pgf_u)*zmfa
1328               zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + &
1329                 zerate*pven(jl,jk)-zderate*zvu(jl,ik)+pgf_v)*zmfa
1330             end if
1331             end if
1332           end if
1333         end do
1334       end do
1335
1336       if(lmfdd) then
1337       do jk = 3 , klev
1338         ik = jk - 1
1339         do jl = 1,klon
1340           if ( ldcum(jl) ) then
1341             if ( jk == idtop(jl) ) then
1342               zud(jl,jk) = 0.5*(zuu(jl,jk)+puen(jl,ik))
1343               zvd(jl,jk) = 0.5*(zvu(jl,jk)+pven(jl,ik))
1344             else if ( jk > idtop(jl) ) then
1345               zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pmfdde_rate(jl,jk)
1346               zmfa = 1./min(-cmfcmin,pmfd(jl,jk))
1347               zud(jl,jk) = (zud(jl,ik)*pmfd(jl,ik) - &
1348                 zerate*puen(jl,ik)+pmfdde_rate(jl,jk)*zud(jl,ik))*zmfa
1349               zvd(jl,jk) = (zvd(jl,ik)*pmfd(jl,ik) - &
1350                 zerate*pven(jl,ik)+pmfdde_rate(jl,jk)*zvd(jl,ik))*zmfa
1351             end if
1352           end if
1353         end do
1354       end do
1355       end if
1356 !   --------------------------------------------------
1357 !   rescale massfluxes for stability in Momentum
1358 !------------------------------------------------------------------------
1359       zmfs(:) = 1.
1360       do jk = 2 , klev
1361         do jl = 1, klon
1362           if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1363             zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons
1364             if ( pmfu(jl,jk) > zmfmax .and. jk >= kctop(jl) ) then
1365               zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk))
1366             end if
1367           end if
1368         end do
1369       end do
1370       do jk = 1 , klev
1371         do jl = 1, klon
1372           zmfuus(jl,jk) = pmfu(jl,jk)
1373           zmfdus(jl,jk) = pmfd(jl,jk)
1374           if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1375             zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl)
1376             zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl)
1377           end if
1378         end do
1379       end do
1380 !*  9.1          update u and v in subroutine cududvn
1381 !-------------------------------------------------------------------
1382      do jk = 1 , klev
1383         do jl = 1, klon
1384           ztenu(jl,jk) = pvom(jl,jk)
1385           ztenv(jl,jk) = pvol(jl,jk)
1386         end do
1387       end do
1388
1389       call cududvn(klon,klev,itopm2,ktype,kcbot,kctop, &
1390                   ldcum,ztmst,paph,puen,pven,zmfuus,zmfdus,zuu,  &
1391                   zud,zvu,zvd,pvom,pvol)
1392
1393 !  calculate KE dissipation
1394       do jl = 1, klon
1395         zsum12(jl) = 0.
1396         zsum22(jl) = 0.
1397       end do
1398         do jk = 1 , klev
1399           do jl = 1, klon
1400             zuv2(jl,jk) = 0.
1401             if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then
1402               zdz = (paph(jl,jk+1)-paph(jl,jk))
1403               zduten = pvom(jl,jk) - ztenu(jl,jk)
1404               zdvten = pvol(jl,jk) - ztenv(jl,jk)
1405               zuv2(jl,jk) = sqrt(zduten**2+zdvten**2)
1406               zsum22(jl) = zsum22(jl) + zuv2(jl,jk)*zdz
1407               zsum12(jl) = zsum12(jl) - &
1408                 (puen(jl,jk)*zduten+pven(jl,jk)*zdvten)*zdz
1409             end if
1410           end do
1411         end do
1412         do jk = 1 , klev
1413           do jl = 1, klon
1414             if ( ldcum(jl) .and. jk>=kctop(jl)-1 ) then
1415               ztdis = rcpd*zsum12(jl)*zuv2(jl,jk)/max(1.e-15,zsum22(jl))
1416               ptte(jl,jk) = ptte(jl,jk) + ztdis
1417             end if
1418           end do
1419         end do
1420
1421       end if
1422
1423 !----------------------------------------------------------------------
1424 !*   10.           IN CASE THAT EITHER DEEP OR SHALLOW IS SWITCHED OFF
1425 ! NEED TO SET SOME VARIABLES A POSTERIORI TO ZERO
1426 ! ---------------------------------------------------
1427     if ( .not. lmfscv .or. .not. lmfpen ) then
1428       do jk = 2 , klev
1429         do jl = 1, klon
1430           if ( llo2(jl) .and. jk >= kctop(jl)-1 ) then
1431             ptu(jl,jk) = pten(jl,jk)
1432             pqu(jl,jk) = pqen(jl,jk)
1433             plu(jl,jk) = 0.
1434             pmfude_rate(jl,jk) = 0.
1435             pmfdde_rate(jl,jk) = 0.
1436           end if
1437         end do
1438       end do
1439       do jl = 1, klon
1440         if ( llo2(jl) ) then
1441           kctop(jl) = klev - 1
1442           kcbot(jl) = klev - 1
1443         end if
1444       end do
1445     end if
1446
1447       return
1448       end subroutine cumastrn
1449
1450 !**********************************************
1451 !    level 3  subroutine cuinin
1452 !**********************************************
1453 !
1454       subroutine cuinin &
1455      &    (klon,     klev,     klevp1,   klevm1,   pten,&
1456      &     pqen,     pqsen,    puen,     pven,     pverv,&
1457      &     pgeo,     paph,     pgeoh,    ptenh,    pqenh,&
1458      &     pqsenh,   klwmin,   ptu,      pqu,      ptd,&
1459      &     pqd,      puu,      pvu,      pud,      pvd,&
1460      &     pmfu,     pmfd,     pmfus,    pmfds,    pmfuq,&
1461      &     pmfdq,    pdmfup,   pdmfdp,   pdpmel,   plu,&
1462      &     plude,    klab)
1463       implicit none
1464 !      m.tiedtke         e.c.m.w.f.     12/89
1465 !***purpose
1466 !   -------
1467 !          this routine interpolates large-scale fields of t,q etc.
1468 !          to half levels (i.e. grid for massflux scheme),
1469 !          and initializes values for updrafts and downdrafts
1470 !***interface
1471 !   ---------
1472 !          this routine is called from *cumastr*.
1473 !***method.
1474 !  --------
1475 !          for extrapolation to half levels see tiedtke(1989)
1476 !***externals
1477 !   ---------
1478 !          *cuadjtq* to specify qs at half levels
1479 ! ----------------------------------------------------------------
1480       integer  klon,klev,klevp1,klevm1
1481       real     pten(klon,klev),        pqen(klon,klev),&
1482      &         puen(klon,klev),        pven(klon,klev),&
1483      &         pqsen(klon,klev),       pverv(klon,klev),&
1484      &         pgeo(klon,klev),        pgeoh(klon,klevp1),&
1485      &         paph(klon,klevp1),      ptenh(klon,klev),&
1486      &         pqenh(klon,klev),       pqsenh(klon,klev)
1487       real     ptu(klon,klev),         pqu(klon,klev),&
1488      &         ptd(klon,klev),         pqd(klon,klev),&
1489      &         puu(klon,klev),         pud(klon,klev),&
1490      &         pvu(klon,klev),         pvd(klon,klev),&
1491      &         pmfu(klon,klev),        pmfd(klon,klev),&
1492      &         pmfus(klon,klev),       pmfds(klon,klev),&
1493      &         pmfuq(klon,klev),       pmfdq(klon,klev),&
1494      &         pdmfup(klon,klev),      pdmfdp(klon,klev),&
1495      &         plu(klon,klev),         plude(klon,klev)
1496       real     zwmax(klon),            zph(klon), &
1497      &         pdpmel(klon,klev)
1498       integer  klab(klon,klev),        klwmin(klon)
1499       logical  loflag(klon)
1500 !  local variables
1501       integer  jl,jk
1502       integer  icall,ik
1503       real     zzs
1504 !------------------------------------------------------------
1505 !*    1.       specify large scale parameters at half levels
1506 !*             adjust temperature fields if staticly unstable
1507 !*             find level of maximum vertical velocity
1508 ! -----------------------------------------------------------
1509       do jk=2,klev
1510       do jl=1,klon
1511         ptenh(jl,jk)=(max(cpd*pten(jl,jk-1)+pgeo(jl,jk-1), &
1512      &             cpd*pten(jl,jk)+pgeo(jl,jk))-pgeoh(jl,jk))*rcpd
1513         pqenh(jl,jk) = pqen(jl,jk-1)
1514         pqsenh(jl,jk)= pqsen(jl,jk-1)
1515         zph(jl)=paph(jl,jk)
1516         loflag(jl)=.true.
1517       end do
1518
1519       if ( jk >= klev-1 .or. jk < 2 ) cycle
1520       ik=jk
1521       icall=0
1522       call cuadjtqn(klon,klev,ik,zph,ptenh,pqsenh,loflag,icall)
1523       do jl=1,klon
1524         pqenh(jl,jk)=min(pqen(jl,jk-1),pqsen(jl,jk-1)) &
1525      &            +(pqsenh(jl,jk)-pqsen(jl,jk-1))
1526         pqenh(jl,jk)=max(pqenh(jl,jk),0.)
1527       end do
1528       end do
1529
1530       do jl=1,klon
1531         ptenh(jl,klev)=(cpd*pten(jl,klev)+pgeo(jl,klev)- &
1532      &                pgeoh(jl,klev))*rcpd
1533         pqenh(jl,klev)=pqen(jl,klev)
1534         ptenh(jl,1)=pten(jl,1)
1535         pqenh(jl,1)=pqen(jl,1)
1536         klwmin(jl)=klev
1537         zwmax(jl)=0.
1538       end do
1539
1540       do jk=klevm1,2,-1
1541       do jl=1,klon
1542         zzs=max(cpd*ptenh(jl,jk)+pgeoh(jl,jk), &
1543      &        cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))
1544         ptenh(jl,jk)=(zzs-pgeoh(jl,jk))*rcpd
1545       end do
1546       end do
1547
1548       do jk=klev,3,-1
1549       do jl=1,klon
1550         if(pverv(jl,jk).lt.zwmax(jl)) then
1551            zwmax(jl)=pverv(jl,jk)
1552            klwmin(jl)=jk
1553         end if
1554       end do
1555       end do
1556 !-----------------------------------------------------------
1557 !*    2.0      initialize values for updrafts and downdrafts
1558 !-----------------------------------------------------------
1559       do jk=1,klev
1560       ik=jk-1
1561       if(jk.eq.1) ik=1
1562       do jl=1,klon
1563       ptu(jl,jk)=ptenh(jl,jk)
1564       ptd(jl,jk)=ptenh(jl,jk)
1565       pqu(jl,jk)=pqenh(jl,jk)
1566       pqd(jl,jk)=pqenh(jl,jk)
1567       plu(jl,jk)=0.
1568       puu(jl,jk)=puen(jl,ik)
1569       pud(jl,jk)=puen(jl,ik)
1570       pvu(jl,jk)=pven(jl,ik)
1571       pvd(jl,jk)=pven(jl,ik)
1572       klab(jl,jk)=0
1573       end do
1574       end do
1575       return
1576       end subroutine cuinin
1577
1578 !---------------------------------------------------------
1579 !  level 3 souroutines
1580 !--------------------------------------------------------
1581       subroutine cutypen &
1582      &   (  klon,    klev,     klevp1,   klevm1,     pqen,&
1583      &     ptenh,   pqenh,     pqsenh,    pgeoh,     paph,&
1584      &       hfx,     qfx,       pgeo,    pqsen,      pap,&
1585      &      pten,    lndj,       cutu,     cuqu,    culab,&
1586      &     ldcum,   cubot,      cutop,    ktype,    wbase,    culu,   kdpl  )
1587 !      zhang & wang      iprc           2011-2013
1588 !***purpose.
1589 !   --------
1590 !          to produce first guess updraught for cu-parameterizations
1591 !          calculates condensation level, and sets updraught base variables and
1592 !          first guess cloud type
1593 !***interface
1594 !   ---------
1595 !          this routine is called from *cumastr*.
1596 !          input are environm. values of t,q,p,phi at half levels.
1597 !          it returns cloud types as follows;
1598 !                 ktype=1 for deep cumulus
1599 !                 ktype=2 for shallow cumulus
1600 !***method.
1601 !  --------
1602 !          based on a simplified updraught equation
1603 !            partial(hup)/partial(z)=eta(h - hup)
1604 !            eta is the entrainment rate for test parcel
1605 !            h stands for dry static energy or the total water specific humidity
1606 !            references: christian jakob, 2003: a new subcloud model for
1607 !            mass-flux convection schemes
1608 !                        influence on triggering, updraft properties, and model
1609 !                        climate, mon.wea.rev.
1610 !                        131, 2765-2778
1611 !            and
1612 !                        ifs documentation - cy36r1,cy38r1
1613 !***input variables:
1614 !       ptenh [ztenh] - environment temperature on half levels. (cuini)
1615 !       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
1616 !       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
1617 !       paph - pressure of half levels. (mssflx)
1618 !       rho  - density of the lowest model level
1619 !       qfx  - net upward moisture flux at the surface (kg/m^2/s)
1620 !       hfx  - net upward heat flux at the surface (w/m^2)
1621 !***variables output by cutype:
1622 !       ktype - convection type - 1: penetrative  (cumastr)
1623 !                                 2: stratocumulus (cumastr)
1624 !                                 3: mid-level  (cuasc)
1625 !       information for updraft parcel (ptu,pqu,plu,kcbot,klab,kdpl...)
1626 ! ----------------------------------------------------------------
1627 !-------------------------------------------------------------------
1628       implicit none
1629 !-------------------------------------------------------------------
1630       integer  klon, klev, klevp1, klevm1
1631       real     ptenh(klon,klev),       pqenh(klon,klev),&
1632      &         pqsen(klon,klev),       pqsenh(klon,klev),&
1633      &         pgeoh(klon,klevp1),     paph(klon,klevp1),&
1634      &         pap(klon,klev),         pqen(klon,klev)
1635       real     pten(klon,klev)
1636       real     ptu(klon,klev),pqu(klon,klev),plu(klon,klev)
1637       real     pgeo(klon,klev)
1638       integer  klab(klon,klev)
1639       integer  kctop(klon),kcbot(klon)
1640
1641       real     qfx(klon),hfx(klon)
1642       real     zph(klon)
1643       integer  lndj(klon)
1644       logical  loflag(klon), deepflag(klon), resetflag(klon)
1645
1646 ! output variables
1647       real     cutu(klon,klev), cuqu(klon,klev), culu(klon,klev)
1648       integer  culab(klon,klev)
1649       real     wbase(klon)
1650       integer  ktype(klon),cubot(klon),cutop(klon),kdpl(klon)
1651       logical  ldcum(klon)
1652
1653 ! local variables
1654       real     zqold(klon)
1655       real     rho, part1, part2, root, conw, deltt, deltq
1656       real     eta(klon),dz(klon),coef(klon)
1657       real     dhen(klon,klev), dh(klon,klev)
1658       real     plude(klon,klev)
1659       real     kup(klon,klev)
1660       real     vptu(klon,klev),vten(klon,klev)
1661       real     zbuo(klon,klev),abuoy(klon,klev)
1662
1663       real     zz,zdken,zdq
1664       real     fscale,crirh1,pp
1665       real     atop1,atop2,abot
1666       real     tmix,zmix,qmix,pmix
1667       real     zlglac,dp
1668       integer  nk,is,ikb,ikt
1669
1670       real     zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp
1671       real     zpdifftop, zpdiffbot
1672       integer  zcbase(klon), itoppacel(klon)
1673       integer  jl,jk,ik,icall,levels
1674       logical  needreset, lldcum(klon)
1675 !--------------------------------------------------------------
1676       do jl=1,klon
1677         kcbot(jl)=klev
1678         kctop(jl)=klev
1679         kdpl(jl) =klev
1680         ktype(jl)=0
1681         wbase(jl)=0.
1682         ldcum(jl)=.false.
1683       end do
1684
1685 !-----------------------------------------------------------
1686 ! let's do test,and check the shallow convection first
1687 ! the first level is klev
1688 ! define deltat and deltaq
1689 !-----------------------------------------------------------
1690       do jk=1,klev
1691       do jl=1,klon
1692           plu(jl,jk)=culu(jl,jk)  ! parcel liquid water
1693           ptu(jl,jk)=cutu(jl,jk)  ! parcel temperature
1694           pqu(jl,jk)=cuqu(jl,jk)  ! parcel specific humidity
1695           klab(jl,jk)=culab(jl,jk)
1696            dh(jl,jk)=0.0  ! parcel dry static energy
1697          dhen(jl,jk)=0.0  ! environment dry static energy
1698           kup(jl,jk)=0.0  ! updraught kinetic energy for parcel
1699          vptu(jl,jk)=0.0  ! parcel virtual temperature considering water-loading
1700          vten(jl,jk)=0.0  ! environment virtual temperature
1701          zbuo(jl,jk)=0.0  ! parcel buoyancy
1702          abuoy(jl,jk)=0.0
1703       end do
1704       end do
1705
1706       do jl=1,klon
1707          zqold(jl) = 0.
1708          lldcum(jl) = .false.
1709          loflag(jl) = .true.
1710       end do
1711
1712 ! check the levels from lowest level to second top level
1713       do jk=klevm1,2,-1
1714
1715 ! define the variables at the first level
1716       if(jk .eq. klevm1) then
1717       do jl=1,klon
1718         rho=pap(jl,klev)/ &
1719      &         (rd*(pten(jl,klev)*(1.+vtmpc1*pqen(jl,klev))))
1720         part1 = 1.5*0.4*pgeo(jl,klev)/ &
1721      &              (rho*pten(jl,klev))
1722         part2 = -hfx(jl)*rcpd-vtmpc1*pten(jl,klev)*qfx(jl)
1723         root  = 0.001-part1*part2
1724         if(part2 .lt. 0.) then
1725            conw  = 1.2*(root)**t13
1726            deltt = max(1.5*hfx(jl)/(rho*cpd*conw),0.)
1727            deltq = max(1.5*qfx(jl)/(rho*conw),0.)
1728            kup(jl,klev) = 0.5*(conw**2)
1729            pqu(jl,klev)= pqenh(jl,klev) + deltq
1730           dhen(jl,klev)= pgeoh(jl,klev) + ptenh(jl,klev)*cpd
1731            dh(jl,klev) = dhen(jl,klev)  + deltt*cpd
1732           ptu(jl,klev) = (dh(jl,klev)-pgeoh(jl,klev))*rcpd
1733           vptu(jl,klev)=ptu(jl,klev)*(1.+vtmpc1*pqu(jl,klev))
1734           vten(jl,klev)=ptenh(jl,klev)*(1.+vtmpc1*pqenh(jl,klev))
1735           zbuo(jl,klev)=(vptu(jl,klev)-vten(jl,klev))/vten(jl,klev)
1736           klab(jl,klev) = 1
1737         else
1738           loflag(jl) = .false.
1739         end if
1740       end do
1741       end if
1742
1743       is=0
1744       do jl=1,klon
1745          if(loflag(jl))then
1746             is=is+1
1747          endif
1748       enddo
1749       if(is.eq.0) exit
1750
1751 ! the next levels, we use the variables at the first level as initial values
1752       do jl=1,klon
1753       if(loflag(jl)) then
1754         eta(jl) = 0.8/(pgeo(jl,jk)*zrg)+2.e-4
1755         dz(jl)  = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1756         coef(jl)= 0.5*eta(jl)*dz(jl)
1757         dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
1758         dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
1759      &              +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
1760         pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
1761      &              +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
1762         ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
1763         zqold(jl) = pqu(jl,jk)
1764         zph(jl)=paph(jl,jk)
1765       end if
1766       end do
1767 ! check if the parcel is saturated
1768       ik=jk
1769       icall=1
1770       call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1771       do jl=1,klon
1772         if( loflag(jl) ) then
1773           zdq = max((zqold(jl) - pqu(jl,jk)),0.)
1774           plu(jl,jk) = plu(jl,jk+1) + zdq
1775           zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
1776                       (1.-foealfa(ptu(jl,jk+1))))
1777           plu(jl,jk) = min(plu(jl,jk),5.e-3)
1778           dh(jl,jk) =  pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
1779 ! compute the updraft speed
1780           vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
1781                         ralfdcp*zlglac
1782           vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
1783           zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
1784           abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
1785           atop1 = 1.0 - 2.*coef(jl)
1786           atop2 = 2.0*dz(jl)*abuoy(jl,jk)
1787           abot =  1.0 + 2.*coef(jl)
1788           kup(jl,jk)  = (atop1*kup(jl,jk+1) + atop2) / abot
1789
1790 ! let's find the exact cloud base
1791          if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
1792               ik = jk + 1
1793               zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
1794               zqsu = min(0.5,zqsu)
1795               zcor = 1./(1.-vtmpc1*zqsu)
1796               zqsu = zqsu*zcor
1797               zdq = min(0.,pqu(jl,ik)-zqsu)
1798               zalfaw = foealfa(ptu(jl,ik))
1799               zfacw = c5les/((ptu(jl,ik)-c4les)**2)
1800               zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
1801               zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
1802               zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
1803               zcor = 1./(1.-vtmpc1*zesdp)
1804               zdqsdt = zfac*zcor*zqsu
1805               zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
1806               zdp = zdq/(zdqsdt*zdtdp)
1807               zcbase(jl) = paph(jl,ik) + zdp
1808 ! chose nearest half level as cloud base (jk or jk+1)
1809               zpdifftop = zcbase(jl) - paph(jl,jk)
1810               zpdiffbot = paph(jl,jk+1) - zcbase(jl)
1811               if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
1812                 ikb = min(klev-1,jk+1)
1813                 klab(jl,ikb) = 2
1814                 klab(jl,jk) = 2
1815                 kcbot(jl) = ikb
1816                 plu(jl,jk+1) = 1.0e-8
1817               else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
1818                 klab(jl,jk) = 2
1819                 kcbot(jl) = jk
1820               end if
1821           end if
1822
1823           if(kup(jl,jk) .lt. 0.)then
1824             loflag(jl) = .false.
1825             if(plu(jl,jk+1) .gt. 0.) then
1826               kctop(jl) = jk
1827               lldcum(jl) = .true.
1828             else
1829               lldcum(jl) = .false.
1830             end if
1831           else
1832             if(plu(jl,jk) .gt. 0.)then
1833               klab(jl,jk)=2
1834             else
1835               klab(jl,jk)=1
1836             end if
1837           end if
1838         end if
1839       end do
1840
1841       end do ! end all the levels
1842
1843       do jl=1,klon
1844         ikb = kcbot(jl)
1845         ikt = kctop(jl)
1846         if(paph(jl,ikb) - paph(jl,ikt) > zdnoprc) lldcum(jl) = .false.
1847         if(lldcum(jl)) then
1848             ktype(jl) = 2
1849             ldcum(jl) = .true.
1850             wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.))
1851             cubot(jl) = ikb
1852             cutop(jl) = ikt
1853             kdpl(jl)  = klev
1854         else
1855             cutop(jl) = -1
1856             cubot(jl) = -1
1857             kdpl(jl)  = klev - 1
1858             ldcum(jl) = .false.
1859             wbase(jl) = 0.
1860         end if
1861       end do
1862
1863       do jk=klev,1,-1
1864        do jl=1,klon
1865            ikt = kctop(jl)
1866            if(jk .ge. ikt)then
1867              culab(jl,jk) = klab(jl,jk)
1868              cutu(jl,jk)  = ptu(jl,jk)
1869              cuqu(jl,jk)  = pqu(jl,jk)
1870              culu(jl,jk)  = plu(jl,jk)
1871            end if
1872        end do
1873       end do
1874
1875 !-----------------------------------------------------------
1876 ! next, let's check the deep convection
1877 ! the first level is klevm1-1
1878 ! define deltat and deltaq
1879 !----------------------------------------------------------
1880 ! we check the parcel starting level by level
1881 ! assume the mix-layer is 60hPa
1882       deltt = 0.2
1883       deltq = 1.0e-4
1884       do jl=1,klon
1885         deepflag(jl) = .false.
1886       end do
1887
1888       do jk=klev,1,-1
1889        do jl=1,klon
1890          if((paph(jl,klev+1)-paph(jl,jk)) .lt. 350.e2) itoppacel(jl) = jk
1891        end do
1892       end do
1893
1894       do levels=klevm1-1,klev/2+1,-1 ! loop starts
1895         do jk=1,klev
1896           do jl=1,klon
1897              plu(jl,jk)=0.0  ! parcel liquid water
1898              ptu(jl,jk)=0.0  ! parcel temperature
1899              pqu(jl,jk)=0.0  ! parcel specific humidity
1900              dh(jl,jk)=0.0   ! parcel dry static energy
1901              dhen(jl,jk)=0.0  ! environment dry static energy
1902              kup(jl,jk)=0.0   ! updraught kinetic energy for parcel
1903              vptu(jl,jk)=0.0  ! parcel virtual temperature consideringwater-loading
1904              vten(jl,jk)=0.0  ! environment virtual temperature
1905              abuoy(jl,jk)=0.0
1906              zbuo(jl,jk)=0.0
1907              klab(jl,jk)=0
1908           end do
1909         end do
1910
1911         do jl=1,klon
1912            kcbot(jl)    =  levels
1913            kctop(jl)    =  levels
1914            zqold(jl)    = 0.
1915            lldcum(jl)   = .false.
1916            resetflag(jl)= .false.
1917            loflag(jl)   = (.not. deepflag(jl)) .and. (levels.ge.itoppacel(jl))
1918         end do
1919
1920 ! start the inner loop to search the deep convection points
1921       do jk=levels,2,-1
1922         is=0
1923         do jl=1,klon
1924          if(loflag(jl))then
1925             is=is+1
1926          endif
1927         enddo
1928         if(is.eq.0) exit
1929
1930 ! define the variables at the departure level
1931         if(jk .eq. levels) then
1932           do jl=1,klon
1933           if(loflag(jl)) then
1934             if((paph(jl,klev+1)-paph(jl,jk)) < 60.e2) then
1935               tmix=0.
1936               qmix=0.
1937               zmix=0.
1938               pmix=0.
1939               do nk=jk+2,jk,-1
1940               if(pmix < 50.e2) then
1941                 dp = paph(jl,nk) - paph(jl,nk-1)
1942                 tmix=tmix+dp*ptenh(jl,nk)
1943                 qmix=qmix+dp*pqenh(jl,nk)
1944                 zmix=zmix+dp*pgeoh(jl,nk)
1945                 pmix=pmix+dp
1946               end if
1947               end do
1948               tmix=tmix/pmix
1949               qmix=qmix/pmix
1950               zmix=zmix/pmix
1951             else
1952               tmix=ptenh(jl,jk+1)
1953               qmix=pqenh(jl,jk+1)
1954               zmix=pgeoh(jl,jk+1)
1955             end if
1956
1957             pqu(jl,jk+1) = qmix + deltq
1958             dhen(jl,jk+1)= zmix + tmix*cpd
1959             dh(jl,jk+1)  = dhen(jl,jk+1) + deltt*cpd
1960             ptu(jl,jk+1) = (dh(jl,jk+1)-pgeoh(jl,jk+1))*rcpd
1961             kup(jl,jk+1) = 0.5
1962             klab(jl,jk+1)= 1
1963             vptu(jl,jk+1)=ptu(jl,jk+1)*(1.+vtmpc1*pqu(jl,jk+1))
1964             vten(jl,jk+1)=ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))
1965             zbuo(jl,jk+1)=(vptu(jl,jk+1)-vten(jl,jk+1))/vten(jl,jk+1)
1966           end if
1967           end do
1968         end if
1969
1970 ! the next levels, we use the variables at the first level as initial values
1971         do jl=1,klon
1972            if(loflag(jl)) then
1973 ! define the fscale
1974              fscale = min(1.,(pqsen(jl,jk)/pqsen(jl,levels))**3)
1975              eta(jl) = 1.75e-3*fscale
1976              dz(jl)  = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
1977              coef(jl)= 0.5*eta(jl)*dz(jl)
1978              dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk)
1979              dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))&
1980      &              +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl))
1981              pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))&
1982      &              +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl))
1983              ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd
1984              zqold(jl) = pqu(jl,jk)
1985              zph(jl)=paph(jl,jk)
1986            end if
1987         end do
1988 ! check if the parcel is saturated
1989         ik=jk
1990         icall=1
1991         call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
1992
1993       do jl=1,klon
1994         if( loflag(jl) ) then
1995           zdq = max((zqold(jl) - pqu(jl,jk)),0.)
1996           plu(jl,jk) = plu(jl,jk+1) + zdq
1997           zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - &
1998                       (1.-foealfa(ptu(jl,jk+1))))
1999           plu(jl,jk) = 0.5*plu(jl,jk)
2000           dh(jl,jk) =  pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac)
2001 ! compute the updraft speed
2002           vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+&
2003                         ralfdcp*zlglac
2004           vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2005           zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk)
2006           abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g
2007           atop1 = 1.0 - 2.*coef(jl)
2008           atop2 = 2.0*dz(jl)*abuoy(jl,jk)
2009           abot =  1.0 + 2.*coef(jl)
2010           kup(jl,jk)  = (atop1*kup(jl,jk+1) + atop2) / abot
2011 ! let's find the exact cloud base
2012           if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then
2013               ik = jk + 1
2014               zqsu = foeewm(ptu(jl,ik))/paph(jl,ik)
2015               zqsu = min(0.5,zqsu)
2016               zcor = 1./(1.-vtmpc1*zqsu)
2017               zqsu = zqsu*zcor
2018               zdq = min(0.,pqu(jl,ik)-zqsu)
2019               zalfaw = foealfa(ptu(jl,ik))
2020               zfacw = c5les/((ptu(jl,ik)-c4les)**2)
2021               zfaci = c5ies/((ptu(jl,ik)-c4ies)**2)
2022               zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci
2023               zesdp = foeewm(ptu(jl,ik))/paph(jl,ik)
2024               zcor = 1./(1.-vtmpc1*zesdp)
2025               zdqsdt = zfac*zcor*zqsu
2026               zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik))
2027               zdp = zdq/(zdqsdt*zdtdp)
2028               zcbase(jl) = paph(jl,ik) + zdp
2029 ! chose nearest half level as cloud base (jk or jk+1)
2030               zpdifftop = zcbase(jl) - paph(jl,jk)
2031               zpdiffbot = paph(jl,jk+1) - zcbase(jl)
2032               if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then
2033                 ikb = min(klev-1,jk+1)
2034                 klab(jl,ikb) = 2
2035                 klab(jl,jk) = 2
2036                 kcbot(jl) = ikb
2037                 plu(jl,jk+1) = 1.0e-8
2038               else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then
2039                 klab(jl,jk) = 2
2040                 kcbot(jl) = jk
2041               end if
2042           end if
2043
2044           if(kup(jl,jk) .lt. 0.)then
2045             loflag(jl) = .false.
2046             if(plu(jl,jk+1) .gt. 0.) then
2047               kctop(jl) = jk
2048               lldcum(jl) = .true.
2049             else
2050               lldcum(jl) = .false.
2051             end if
2052           else
2053             if(plu(jl,jk) .gt. 0.)then
2054               klab(jl,jk)=2
2055             else
2056               klab(jl,jk)=1
2057             end if
2058           end if
2059         end if
2060       end do
2061
2062       end do ! end all the levels
2063
2064       needreset = .false.
2065       do jl=1,klon
2066         ikb = kcbot(jl)
2067         ikt = kctop(jl)
2068         if(paph(jl,ikb) - paph(jl,ikt) < zdnoprc) lldcum(jl) = .false.
2069         if(lldcum(jl)) then
2070          ktype(jl)    = 1
2071          ldcum(jl)    = .true.
2072          deepflag(jl) = .true.
2073          wbase(jl)    = sqrt(max(2.*kup(jl,ikb),0.))
2074          cubot(jl)    = ikb
2075          cutop(jl)    = ikt
2076          kdpl(jl)     = levels+1
2077          needreset    = .true.
2078          resetflag(jl)= .true.
2079         end if
2080       end do
2081
2082       if(needreset) then
2083       do jk=klev,1,-1
2084         do jl=1,klon
2085           if(resetflag(jl)) then
2086             ikt = kctop(jl)
2087             ikb = kdpl(jl)
2088             if(jk .le. ikb .and. jk .ge. ikt )then
2089              culab(jl,jk) = klab(jl,jk)
2090              cutu(jl,jk)  = ptu(jl,jk)
2091              cuqu(jl,jk)  = pqu(jl,jk)
2092              culu(jl,jk)  = plu(jl,jk)
2093             else
2094              culab(jl,jk) = 1
2095              cutu(jl,jk)  = ptenh(jl,jk)
2096              cuqu(jl,jk)  = pqenh(jl,jk)
2097              culu(jl,jk)  = 0.
2098             end if
2099             if ( jk .lt. ikt ) culab(jl,jk) = 0
2100           end if
2101         end do
2102       end do
2103       end if
2104
2105       end do ! end all cycles
2106
2107       return
2108       end subroutine cutypen
2109
2110 !-----------------------------------------------------------------
2111 !    level 3 subroutines 'cuascn'
2112 !-----------------------------------------------------------------
2113       subroutine cuascn &
2114      &    (klon,     klev,     klevp1,   klevm1,   ptenh,&
2115      &     pqenh,    puen,     pven,     pten,     pqen,&
2116      &     pqsen,    pgeo,     pgeoh,    pap,      paph,&
2117      &     pqte,     pverv,    klwmin,   ldcum,    phcbase,&
2118      &     ktype,    klab,     ptu,      pqu,      plu,&
2119      &     puu,      pvu,      pmfu,     pmfub,    &
2120      &     pmfus,    pmfuq,    pmful,    plude,    pdmfup,&
2121      &     kcbot,    kctop,    kctop0,   kcum,     ztmst,&
2122      &     pqsenh,   plglac,   lndj,     wup,      wbase,   kdpl, pmfude_rate)
2123       implicit none
2124 !     this routine does the calculations for cloud ascents
2125 !     for cumulus parameterization
2126 !     m.tiedtke         e.c.m.w.f.     7/86 modif.  12/89
2127 !     y.wang            iprc           11/01 modif.
2128 !     c.zhang           iprc           05/12 modif.
2129 !***purpose.
2130 !   --------
2131 !          to produce cloud ascents for cu-parametrization
2132 !          (vertical profiles of t,q,l,u and v and corresponding
2133 !           fluxes as well as precipitation rates)
2134 !***interface
2135 !   ---------
2136 !          this routine is called from *cumastr*.
2137 !***method.
2138 !  --------
2139 !          lift surface air dry-adiabatically to cloud base
2140 !          and then calculate moist ascent for
2141 !          entraining/detraining plume.
2142 !          entrainment and detrainment rates differ for
2143 !          shallow and deep cumulus convection.
2144 !          in case there is no penetrative or shallow convection
2145 !          check for possibility of mid level convection
2146 !          (cloud base values calculated in *cubasmc*)
2147 !***externals
2148 !   ---------
2149 !          *cuadjtqn* adjust t and q due to condensation in ascent
2150 !          *cuentrn*  calculate entrainment/detrainment rates
2151 !          *cubasmcn* calculate cloud base values for midlevel convection
2152 !***reference
2153 !   ---------
2154 !          (tiedtke,1989)
2155 !***input variables:
2156 !       ptenh [ztenh] - environ temperature on half levels. (cuini)
2157 !       pqenh [zqenh] - env. specific humidity on half levels. (cuini)
2158 !       puen - environment wind u-component. (mssflx)
2159 !       pven - environment wind v-component. (mssflx)
2160 !       pten - environment temperature. (mssflx)
2161 !       pqen - environment specific humidity. (mssflx)
2162 !       pqsen - environment saturation specific humidity. (mssflx)
2163 !       pgeo - geopotential. (mssflx)
2164 !       pgeoh [zgeoh] - geopotential on half levels, (mssflx)
2165 !       pap - pressure in pa.  (mssflx)
2166 !       paph - pressure of half levels. (mssflx)
2167 !       pqte - moisture convergence (delta q/delta t). (mssflx)
2168 !       pverv - large scale vertical velocity (omega). (mssflx)
2169 !       klwmin [ilwmin] - level of minimum omega. (cuini)
2170 !       klab [ilab] - level label - 1: sub-cloud layer.
2171 !                                   2: condensation level (cloud base)
2172 !       pmfub [zmfub] - updraft mass flux at cloud base. (cumastr)
2173 !***variables modified by cuasc:
2174 !       ldcum - logical denoting profiles. (cubase)
2175 !       ktype - convection type - 1: penetrative  (cumastr)
2176 !                                 2: stratocumulus (cumastr)
2177 !                                 3: mid-level  (cuasc)
2178 !       ptu - cloud temperature.
2179 !       pqu - cloud specific humidity.
2180 !       plu - cloud liquid water (moisture condensed out)
2181 !       puu [zuu] - cloud momentum u-component.
2182 !       pvu [zvu] - cloud momentum v-component.
2183 !       pmfu - updraft mass flux.
2184 !       pmfus [zmfus] - updraft flux of dry static energy. (cubasmc)
2185 !       pmfuq [zmfuq] - updraft flux of specific humidity.
2186 !       pmful [zmful] - updraft flux of cloud liquid water.
2187 !       plude - liquid water returned to environment by detrainment.
2188 !       pdmfup [zmfup] -
2189 !       kcbot - cloud base level. (cubase)
2190 !       kctop - cloud top level
2191 !       kctop0 [ictop0] - estimate of cloud top. (cumastr)
2192 !       kcum [icum] - flag to control the call
2193
2194       integer  klev,klon,klevp1,klevm1
2195       real     ptenh(klon,klev),       pqenh(klon,klev), &
2196      &         puen(klon,klev),        pven(klon,klev),&
2197      &         pten(klon,klev),        pqen(klon,klev),&
2198      &         pgeo(klon,klev),        pgeoh(klon,klevp1),&
2199      &         pap(klon,klev),         paph(klon,klevp1),&
2200      &         pqsen(klon,klev),       pqte(klon,klev),&
2201      &         pverv(klon,klev),       pqsenh(klon,klev)
2202       real     ptu(klon,klev),         pqu(klon,klev),&
2203      &         puu(klon,klev),         pvu(klon,klev),&
2204      &         pmfu(klon,klev),        zph(klon),&
2205      &         pmfub(klon),            &
2206      &         pmfus(klon,klev),       pmfuq(klon,klev),&
2207      &         plu(klon,klev),         plude(klon,klev),&
2208      &         pmful(klon,klev),       pdmfup(klon,klev)
2209       real     zdmfen(klon),           zdmfde(klon),&
2210      &         zmfuu(klon),            zmfuv(klon),&
2211      &         zpbase(klon),           zqold(klon)
2212       real     phcbase(klon),          zluold(klon)
2213       real     zprecip(klon),          zlrain(klon,klev)
2214       real     zbuo(klon,klev),        kup(klon,klev)
2215       real     wup(klon)
2216       real     wbase(klon),            zodetr(klon,klev)
2217       real     plglac(klon,klev)
2218
2219       real     eta(klon),dz(klon)
2220
2221       integer  klwmin(klon),           ktype(klon),&
2222      &         klab(klon,klev),        kcbot(klon),&
2223      &         kctop(klon),            kctop0(klon)
2224       integer  lndj(klon)
2225       logical  ldcum(klon),            loflag(klon)
2226       logical  llo2,llo3,              llo1(klon)
2227
2228       integer  kdpl(klon)
2229       real     zoentr(klon),           zdpmean(klon)
2230       real     pdmfen(klon,klev),      pmfude_rate(klon,klev)
2231 ! local variables
2232       integer  jl,jk
2233       integer  ikb,icum,itopm2,ik,icall,is,kcum,jlm,jll
2234       integer  jlx(klon)
2235       real     ztmst,zcons2,zfacbuo,zprcdgw,z_cwdrag,z_cldmax,z_cwifrac,z_cprc2
2236       real     zmftest,zmfmax,zqeen,zseen,zscde,zqude
2237       real     zmfusk,zmfuqk,zmfulk
2238       real     zbc,zbe,zkedke,zmfun,zwu,zprcon,zdt,zcbf,zzco
2239       real     zlcrit,zdfi,zc,zd,zint,zlnew,zvw,zvi,zalfaw,zrold
2240       real     zrnew,zz,zdmfeu,zdmfdu,dp
2241       real     zfac,zbuoc,zdkbuo,zdken,zvv,zarg,zchange,zxe,zxs,zdshrd
2242       real     atop1,atop2,abot
2243 !--------------------------------
2244 !*    1.       specify parameters
2245 !--------------------------------
2246       zcons2=3./(g*ztmst)
2247       zfacbuo = 0.5/(1.+0.5)
2248       zprcdgw = cprcon*zrg
2249       z_cldmax = 5.e-3
2250       z_cwifrac = 0.5
2251       z_cprc2 = 0.5
2252       z_cwdrag = (3.0/8.0)*0.506/0.2
2253 !---------------------------------
2254 !     2.        set default values
2255 !---------------------------------
2256       llo3 = .false.
2257       do jl=1,klon
2258         zluold(jl)=0.
2259         wup(jl)=0.
2260         zdpmean(jl)=0.
2261         zoentr(jl)=0.
2262         if(.not.ldcum(jl)) then
2263           ktype(jl)=0
2264           kcbot(jl) = -1
2265           pmfub(jl) = 0.
2266           pqu(jl,klev) = 0.
2267         end if
2268       end do
2269
2270  ! initialize variout quantities
2271       do jk=1,klev
2272       do jl=1,klon
2273           if(jk.ne.kcbot(jl)) plu(jl,jk)=0.
2274           pmfu(jl,jk)=0.
2275           pmfus(jl,jk)=0.
2276           pmfuq(jl,jk)=0.
2277           pmful(jl,jk)=0.
2278           plude(jl,jk)=0.
2279           plglac(jl,jk)=0.
2280           pdmfup(jl,jk)=0.
2281           zlrain(jl,jk)=0.
2282           zbuo(jl,jk)=0.
2283           kup(jl,jk)=0.
2284           pdmfen(jl,jk) = 0.
2285           pmfude_rate(jl,jk) = 0.
2286           if(.not.ldcum(jl).or.ktype(jl).eq.3) klab(jl,jk)=0
2287           if(.not.ldcum(jl).and.paph(jl,jk).lt.4.e4) kctop0(jl)=jk
2288       end do
2289       end do
2290
2291       do jl = 1,klon
2292         if ( ktype(jl) == 3 ) ldcum(jl) = .false.
2293       end do
2294 !------------------------------------------------
2295 !     3.0      initialize values at cloud base level
2296 !------------------------------------------------
2297       do jl=1,klon
2298         kctop(jl)=kcbot(jl)
2299         if(ldcum(jl)) then
2300           ikb = kcbot(jl)
2301           kup(jl,ikb) = 0.5*wbase(jl)**2
2302           pmfu(jl,ikb) = pmfub(jl)
2303           pmfus(jl,ikb) = pmfub(jl)*(cpd*ptu(jl,ikb)+pgeoh(jl,ikb))
2304           pmfuq(jl,ikb) = pmfub(jl)*pqu(jl,ikb)
2305           pmful(jl,ikb) = pmfub(jl)*plu(jl,ikb)
2306         end if
2307       end do
2308 !
2309 !-----------------------------------------------------------------
2310 !     4.       do ascent: subcloud layer (klab=1) ,clouds (klab=2)
2311 !              by doing first dry-adiabatic ascent and then
2312 !              by adjusting t,q and l accordingly in *cuadjtqn*,
2313 !              then check for buoyancy and set flags accordingly
2314 !-----------------------------------------------------------------
2315 !
2316       do jk=klevm1,3,-1
2317 !             specify cloud base values for midlevel convection
2318 !             in *cubasmc* in case there is not already convection
2319 ! ---------------------------------------------------------------------
2320       ik=jk
2321       call cubasmcn&
2322      &    (klon,     klev,     klevm1,   ik,      pten,&
2323      &     pqen,     pqsen,    puen,     pven,    pverv,&
2324      &     pgeo,     pgeoh,    ldcum,   ktype,   klab,  zlrain,&
2325      &     pmfu,     pmfub,    kcbot,   ptu,&
2326      &     pqu,      plu,      puu,     pvu,      pmfus,&
2327      &     pmfuq,    pmful,    pdmfup  )
2328       is = 0
2329       jlm = 0
2330       do jl = 1,klon
2331         loflag(jl) = .false.
2332         zprecip(jl) = 0.
2333         llo1(jl) = .false.
2334         is = is + klab(jl,jk+1)
2335         if ( klab(jl,jk+1) == 0 ) klab(jl,jk) = 0
2336         if ( (ldcum(jl) .and. klab(jl,jk+1) == 2) .or.   &
2337              (ktype(jl) == 3 .and. klab(jl,jk+1) == 1) ) then
2338           loflag(jl) = .true.
2339           jlm = jlm + 1
2340           jlx(jlm) = jl
2341         end if
2342         zph(jl) = paph(jl,jk)
2343         if ( ktype(jl) == 3 .and. jk == kcbot(jl) ) then
2344           zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
2345           if ( pmfub(jl) > zmfmax ) then
2346             zfac = zmfmax/pmfub(jl)
2347             pmfu(jl,jk+1) = pmfu(jl,jk+1)*zfac
2348             pmfus(jl,jk+1) = pmfus(jl,jk+1)*zfac
2349             pmfuq(jl,jk+1) = pmfuq(jl,jk+1)*zfac
2350             pmfub(jl) = zmfmax
2351           end if
2352           pmfub(jl)=min(pmfub(jl),zmfmax)
2353         end if
2354       end do
2355
2356       if(is.gt.0) llo3 = .true.
2357 !
2358 !*     specify entrainment rates in *cuentr*
2359 ! -------------------------------------
2360       ik=jk
2361       call cuentrn(klon,klev,ik,kcbot,ldcum,llo3, &
2362                   pgeoh,pmfu,zdmfen,zdmfde)
2363 !
2364 !      do adiabatic ascent for entraining/detraining plume
2365       if(llo3) then
2366 ! -------------------------------------------------------
2367 !
2368         do jl = 1,klon
2369           zqold(jl) = 0.
2370         end do
2371         do jll = 1 , jlm
2372           jl = jlx(jll)
2373           zdmfde(jl) = min(zdmfde(jl),0.75*pmfu(jl,jk+1))
2374           if ( jk == kcbot(jl) ) then
2375             zoentr(jl) = -1.75e-3*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - &
2376                          1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg
2377             zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk+1)
2378           end if
2379           if ( jk < kcbot(jl) ) then
2380             zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2
2381             zxs = max(pmfu(jl,jk+1)-zmfmax,0.)
2382             wup(jl) = wup(jl) + kup(jl,jk+1)*(pap(jl,jk+1)-pap(jl,jk))
2383             zdpmean(jl) = zdpmean(jl) + pap(jl,jk+1) - pap(jl,jk)
2384             zdmfen(jl) = zoentr(jl)
2385             if ( ktype(jl) >= 2 ) then
2386               zdmfen(jl) = 2.0*zdmfen(jl)
2387               zdmfde(jl) = zdmfen(jl)
2388             end if
2389             zdmfde(jl) = zdmfde(jl) * &
2390                          (1.6-min(1.,pqen(jl,jk)/pqsen(jl,jk)))
2391             zmftest = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2392             zchange = max(zmftest-zmfmax,0.)
2393             zxe = max(zchange-zxs,0.)
2394             zdmfen(jl) = zdmfen(jl) - zxe
2395             zchange = zchange - zxe
2396             zdmfde(jl) = zdmfde(jl) + zchange
2397           end if
2398           pdmfen(jl,jk) = zdmfen(jl) - zdmfde(jl)
2399           pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2400           zqeen = pqenh(jl,jk+1)*zdmfen(jl)
2401           zseen = (cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*zdmfen(jl)
2402           zscde = (cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1))*zdmfde(jl)
2403           zqude = pqu(jl,jk+1)*zdmfde(jl)
2404           plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2405           zmfusk = pmfus(jl,jk+1) + zseen - zscde
2406           zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude
2407           zmfulk = pmful(jl,jk+1) - plude(jl,jk)
2408           plu(jl,jk) = zmfulk*(1./max(cmfcmin,pmfu(jl,jk)))
2409           pqu(jl,jk) = zmfuqk*(1./max(cmfcmin,pmfu(jl,jk)))
2410           ptu(jl,jk) = (zmfusk * &
2411             (1./max(cmfcmin,pmfu(jl,jk)))-pgeoh(jl,jk))*rcpd
2412           ptu(jl,jk) = max(100.,ptu(jl,jk))
2413           ptu(jl,jk) = min(400.,ptu(jl,jk))
2414           zqold(jl) = pqu(jl,jk)
2415           zlrain(jl,jk) = zlrain(jl,jk+1)*(pmfu(jl,jk+1)-zdmfde(jl)) * &
2416                           (1./max(cmfcmin,pmfu(jl,jk)))
2417           zluold(jl) = plu(jl,jk)
2418         end do
2419 ! reset to environmental values if below departure level
2420         do jl = 1,klon
2421           if ( jk > kdpl(jl) ) then
2422             ptu(jl,jk) = ptenh(jl,jk)
2423             pqu(jl,jk) = pqenh(jl,jk)
2424             plu(jl,jk) = 0.
2425             zluold(jl) = plu(jl,jk)
2426           end if
2427         end do
2428 !*             do corrections for moist ascent
2429 !*             by adjusting t,q and l in *cuadjtq*
2430 !------------------------------------------------
2431       ik=jk
2432       icall=1
2433 !
2434       if ( jlm > 0 ) then
2435         call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall)
2436       end if
2437 ! compute the upfraft speed in cloud layer
2438         do jll = 1 , jlm
2439           jl = jlx(jll)
2440           if ( pqu(jl,jk) /= zqold(jl) ) then
2441             plglac(jl,jk) = plu(jl,jk) * &
2442                            ((1.-foealfa(ptu(jl,jk)))- &
2443                             (1.-foealfa(ptu(jl,jk+1))))
2444             ptu(jl,jk) = ptu(jl,jk) + ralfdcp*plglac(jl,jk)
2445           end if
2446         end do
2447         do jll = 1 , jlm
2448           jl = jlx(jll)
2449           if ( pqu(jl,jk) /= zqold(jl) ) then
2450             klab(jl,jk) = 2
2451             plu(jl,jk) = plu(jl,jk) + zqold(jl) - pqu(jl,jk)
2452             zbc = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk+1) - &
2453               zlrain(jl,jk+1))
2454             zbe = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2455             zbuo(jl,jk) = zbc - zbe
2456 ! set flags for the case of midlevel convection
2457             if ( ktype(jl) == 3 .and. klab(jl,jk+1) == 1 ) then
2458               if ( zbuo(jl,jk) > -0.5 ) then
2459                 ldcum(jl) = .true.
2460                 kctop(jl) = jk
2461                 kup(jl,jk) = 0.5
2462               else
2463                 klab(jl,jk) = 0
2464                 pmfu(jl,jk) = 0.
2465                 plude(jl,jk) = 0.
2466                 plu(jl,jk) = 0.
2467               end if
2468             end if
2469             if ( klab(jl,jk+1) == 2 ) then
2470               if ( zbuo(jl,jk) < 0. ) then
2471                 ptenh(jl,jk) = 0.5*(pten(jl,jk)+pten(jl,jk-1))
2472                 pqenh(jl,jk) = 0.5*(pqen(jl,jk)+pqen(jl,jk-1))
2473                 zbuo(jl,jk) = zbc - ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk))
2474               end if
2475               zbuoc = (zbuo(jl,jk) / &
2476                 (ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)))+zbuo(jl,jk+1) / &
2477                 (ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))))*0.5
2478               zdkbuo = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zfacbuo*zbuoc
2479 ! mixing and "pressure" gradient term in upper troposphere
2480               if ( zdmfen(jl) > 0. ) then
2481                 zdken = min(1.,(1.+z_cwdrag)*zdmfen(jl) / &
2482                         max(cmfcmin,pmfu(jl,jk+1)))
2483               else
2484                 zdken = min(1.,(1.+z_cwdrag)*zdmfde(jl) / &
2485                         max(cmfcmin,pmfu(jl,jk+1)))
2486               end if
2487               kup(jl,jk) = (kup(jl,jk+1)*(1.-zdken)+zdkbuo) / &
2488                            (1.+zdken)
2489               if ( zbuo(jl,jk) < 0. ) then
2490                 zkedke = kup(jl,jk)/max(1.e-10,kup(jl,jk+1))
2491                 zkedke = max(0.,min(1.,zkedke))
2492                 zmfun = sqrt(zkedke)*pmfu(jl,jk+1)
2493                 zdmfde(jl) = max(zdmfde(jl),pmfu(jl,jk+1)-zmfun)
2494                 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2495                 pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl)
2496               end if
2497               if ( zbuo(jl,jk) > -0.2  ) then
2498                 ikb = kcbot(jl)
2499                 zoentr(jl) = 1.75e-3*(0.3-(min(1.,pqen(jl,jk-1) /    &
2500                   pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * &
2501                   zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3
2502                 zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk)
2503               else
2504                 zoentr(jl) = 0.
2505               end if
2506 ! erase values if below departure level
2507               if ( jk > kdpl(jl) ) then
2508                 pmfu(jl,jk) = pmfu(jl,jk+1)
2509                 kup(jl,jk) = 0.5
2510               end if
2511               if ( kup(jl,jk) > 0. .and. pmfu(jl,jk) > 0. ) then
2512                 kctop(jl) = jk
2513                 llo1(jl) = .true.
2514               else
2515                 klab(jl,jk) = 0
2516                 pmfu(jl,jk) = 0.
2517                 kup(jl,jk) = 0.
2518                 zdmfde(jl) = pmfu(jl,jk+1)
2519                 plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2520               end if
2521 ! save detrainment rates for updraught
2522               if ( pmfu(jl,jk+1) > 0. ) pmfude_rate(jl,jk) = zdmfde(jl)
2523             end if
2524           else if ( ktype(jl) == 2 .and. pqu(jl,jk) == zqold(jl) ) then
2525             klab(jl,jk) = 0
2526             pmfu(jl,jk) = 0.
2527             kup(jl,jk) = 0.
2528             zdmfde(jl) = pmfu(jl,jk+1)
2529             plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl)
2530             pmfude_rate(jl,jk) = zdmfde(jl)
2531           end if
2532         end do
2533
2534         do jl = 1,klon
2535           if ( llo1(jl) ) then
2536 ! conversions only proceeds if plu is greater than a threshold liquid water
2537 ! content of 0.3 g/kg over water and 0.5 g/kg over land to prevent precipitation
2538 ! generation from small water contents.
2539             if ( lndj(jl).eq.1 ) then
2540               zdshrd = 5.e-4
2541             else
2542               zdshrd = 3.e-4
2543             end if
2544             ikb=kcbot(jl)
2545             if ( plu(jl,jk) > zdshrd )then
2546               zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk+1))))
2547               zprcon = zprcdgw/(0.75*zwu)
2548 ! PARAMETERS FOR BERGERON-FINDEISEN PROCESS (T < -5C)
2549               zdt = min(rtber-rtice,max(rtber-ptu(jl,jk),0.))
2550               zcbf = 1. + z_cprc2*sqrt(zdt)
2551               zzco = zprcon*zcbf
2552               zlcrit = zdshrd/zcbf
2553               zdfi = pgeoh(jl,jk) - pgeoh(jl,jk+1)
2554               zc = (plu(jl,jk)-zluold(jl))
2555               zarg = (plu(jl,jk)/zlcrit)**2
2556               if ( zarg < 25.0 ) then
2557                 zd = zzco*(1.-exp(-zarg))*zdfi
2558               else
2559                 zd = zzco*zdfi
2560               end if
2561               zint = exp(-zd)
2562               zlnew = zluold(jl)*zint + zc/zd*(1.-zint)
2563               zlnew = max(0.,min(plu(jl,jk),zlnew))
2564               zlnew = min(z_cldmax,zlnew)
2565               zprecip(jl) = max(0.,zluold(jl)+zc-zlnew)
2566               pdmfup(jl,jk) = zprecip(jl)*pmfu(jl,jk)
2567               zlrain(jl,jk) = zlrain(jl,jk) + zprecip(jl)
2568               plu(jl,jk) = zlnew
2569             end if
2570           end if
2571         end do
2572         do jl = 1, klon
2573           if ( llo1(jl) ) then
2574             if ( zlrain(jl,jk) > 0. ) then
2575               zvw = 21.18*zlrain(jl,jk)**0.2
2576               zvi = z_cwifrac*zvw
2577               zalfaw = foealfa(ptu(jl,jk))
2578               zvv = zalfaw*zvw + (1.-zalfaw)*zvi
2579               zrold = zlrain(jl,jk) - zprecip(jl)
2580               zc = zprecip(jl)
2581               zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk))))
2582               zd = zvv/zwu
2583               zint = exp(-zd)
2584               zrnew = zrold*zint + zc/zd*(1.-zint)
2585               zrnew = max(0.,min(zlrain(jl,jk),zrnew))
2586               zlrain(jl,jk) = zrnew
2587             end if
2588           end if
2589         end do
2590         do jll = 1 , jlm
2591           jl = jlx(jll)
2592           pmful(jl,jk) = plu(jl,jk)*pmfu(jl,jk)
2593           pmfus(jl,jk) = (cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk)
2594           pmfuq(jl,jk) = pqu(jl,jk)*pmfu(jl,jk)
2595         end do
2596       end if
2597     end do
2598 !----------------------------------------------------------------------
2599 ! 5.       final calculations
2600 ! ------------------
2601       do jl = 1,klon
2602        if ( kctop(jl) == -1 ) ldcum(jl) = .false.
2603         kcbot(jl) = max(kcbot(jl),kctop(jl))
2604         if ( ldcum(jl) ) then
2605           wup(jl) = max(1.e-2,wup(jl)/max(1.,zdpmean(jl)))
2606           wup(jl) = sqrt(2.*wup(jl))
2607         end if
2608       end do
2609
2610       return
2611       end subroutine cuascn
2612 !---------------------------------------------------------
2613 !  level 3 souroutines
2614 !--------------------------------------------------------
2615       subroutine cudlfsn   &
2616      &    (klon,     klev,  &
2617      &     kcbot,    kctop,    lndj,   ldcum,   &
2618      &     ptenh,    pqenh,    puen,     pven,     &
2619      &     pten,     pqsen,    pgeo,                &
2620      &     pgeoh,    paph,     ptu,      pqu,      plu,&
2621      &     puu,      pvu,      pmfub,    prfl,          &
2622      &     ptd,      pqd,      pud,      pvd,            &
2623      &     pmfd,     pmfds,    pmfdq,    pdmfdp,          &
2624      &     kdtop,    lddraf)
2625
2626 !          this routine calculates level of free sinking for
2627 !          cumulus downdrafts and specifies t,q,u and v values
2628
2629 !          m.tiedtke         e.c.m.w.f.    12/86 modif. 12/89
2630
2631 !          purpose.
2632 !          --------
2633 !          to produce lfs-values for cumulus downdrafts
2634 !          for massflux cumulus parameterization
2635
2636 !          interface
2637 !          ---------
2638 !          this routine is called from *cumastr*.
2639 !          input are environmental values of t,q,u,v,p,phi
2640 !          and updraft values t,q,u and v and also
2641 !          cloud base massflux and cu-precipitation rate.
2642 !          it returns t,q,u and v values and massflux at lfs.
2643
2644 !          method.
2645
2646 !          check for negative buoyancy of air of equal parts of
2647 !          moist environmental air and cloud air.
2648
2649 !     parameter     description                                   units
2650 !     ---------     -----------                                   -----
2651 !     input parameters (integer):
2652
2653 !    *klon*         number of grid points per packet
2654 !    *klev*         number of levels
2655 !    *kcbot*        cloud base level
2656 !    *kctop*        cloud top level
2657
2658 !    input parameters (logical):
2659
2660 !    *lndj*       land sea mask (1 for land)
2661 !    *ldcum*        flag: .true. for convective points
2662
2663 !    input parameters (real):
2664
2665 !    *ptenh*        env. temperature (t+1) on half levels          k
2666 !    *pqenh*        env. spec. humidity (t+1) on half levels     kg/kg
2667 !    *puen*         provisional environment u-velocity (t+1)      m/s
2668 !    *pven*         provisional environment v-velocity (t+1)      m/s
2669 !    *pten*         provisional environment temperature (t+1)       k
2670 !    *pqsen*        environment spec. saturation humidity (t+1)   kg/kg
2671 !    *pgeo*         geopotential                                  m2/s2
2672 !    *pgeoh*        geopotential on half levels                  m2/s2
2673 !    *paph*         provisional pressure on half levels           pa
2674 !    *ptu*          temperature in updrafts                        k
2675 !    *pqu*          spec. humidity in updrafts                   kg/kg
2676 !    *plu*          liquid water content in updrafts             kg/kg
2677 !    *puu*          u-velocity in updrafts                        m/s
2678 !    *pvu*          v-velocity in updrafts                        m/s
2679 !    *pmfub*        massflux in updrafts at cloud base           kg/(m2*s)
2680
2681 !    updated parameters (real):
2682
2683 !    *prfl*         precipitation rate                           kg/(m2*s)
2684
2685 !    output parameters (real):
2686
2687 !    *ptd*          temperature in downdrafts                      k
2688 !    *pqd*          spec. humidity in downdrafts                 kg/kg
2689 !    *pud*          u-velocity in downdrafts                      m/s
2690 !    *pvd*          v-velocity in downdrafts                      m/s
2691 !    *pmfd*         massflux in downdrafts                       kg/(m2*s)
2692 !    *pmfds*        flux of dry static energy in downdrafts       j/(m2*s)
2693 !    *pmfdq*        flux of spec. humidity in downdrafts         kg/(m2*s)
2694 !    *pdmfdp*       flux difference of precip. in downdrafts     kg/(m2*s)
2695
2696 !    output parameters (integer):
2697
2698 !    *kdtop*        top level of downdrafts
2699
2700 !    output parameters (logical):
2701
2702 !    *lddraf*       .true. if downdrafts exist
2703
2704 !          externals
2705 !          ---------
2706 !          *cuadjtq* for calculating wet bulb t and q at lfs
2707 !----------------------------------------------------------------------
2708       implicit none
2709
2710       integer  klev,klon
2711       real     ptenh(klon,klev),       pqenh(klon,klev), &
2712      &         puen(klon,klev),        pven(klon,klev),  &
2713      &         pten(klon,klev),        pqsen(klon,klev),  &
2714      &         pgeo(klon,klev),                       &
2715      &         pgeoh(klon,klev+1),     paph(klon,klev+1),&
2716      &         ptu(klon,klev),         pqu(klon,klev),   &
2717      &         puu(klon,klev),         pvu(klon,klev),   &
2718      &         plu(klon,klev),                          &
2719      &         pmfub(klon),            prfl(klon)
2720
2721       real     ptd(klon,klev),         pqd(klon,klev),   &
2722      &         pud(klon,klev),         pvd(klon,klev),    &
2723      &         pmfd(klon,klev),        pmfds(klon,klev),  &
2724      &         pmfdq(klon,klev),       pdmfdp(klon,klev)
2725       integer  kcbot(klon),            kctop(klon),       &
2726      &         kdtop(klon),            ikhsmin(klon)
2727       logical  ldcum(klon),                              &
2728      &         lddraf(klon)
2729       integer  lndj(klon)
2730
2731       real     ztenwb(klon,klev),      zqenwb(klon,klev), &
2732      &         zcond(klon),            zph(klon),         &
2733      &         zhsmin(klon)
2734       logical  llo2(klon)
2735 ! local variables
2736       integer  jl,jk
2737       integer  is,ik,icall,ike
2738       real     zhsk,zttest,zqtest,zbuo,zmftop
2739
2740 !----------------------------------------------------------------------
2741
2742 !     1.           set default values for downdrafts
2743 !                  ---------------------------------
2744       do jl=1,klon
2745         lddraf(jl)=.false.
2746         kdtop(jl)=klev+1
2747         ikhsmin(jl)=klev+1
2748         zhsmin(jl)=1.e8
2749       enddo
2750 !----------------------------------------------------------------------
2751
2752 !     2.           determine level of free sinking:
2753 !                  downdrafts shall start at model level of minimum
2754 !                  of saturation moist static energy or below
2755 !                  respectively
2756
2757 !                  for every point and proceed as follows:
2758
2759 !                    (1) determine level of minimum of hs
2760 !                    (2) determine wet bulb environmental t and q
2761 !                    (3) do mixing with cumulus cloud air
2762 !                    (4) check for negative buoyancy
2763 !                    (5) if buoyancy>0 repeat (2) to (4) for next
2764 !                        level below
2765
2766 !                  the assumption is that air of downdrafts is mixture
2767 !                  of 50% cloud air + 50% environmental air at wet bulb
2768 !                  temperature (i.e. which became saturated due to
2769 !                  evaporation of rain and cloud water)
2770 !                  ----------------------------------------------------
2771       do jk=3,klev-2
2772          do jl=1,klon
2773            zhsk=cpd*pten(jl,jk)+pgeo(jl,jk) +  &
2774      &         foelhm(pten(jl,jk))*pqsen(jl,jk)
2775            if(zhsk .lt. zhsmin(jl)) then
2776               zhsmin(jl) = zhsk
2777               ikhsmin(jl)= jk
2778            end if
2779          end do
2780       end do
2781
2782
2783       ike=klev-3
2784       do jk=3,ike
2785
2786 !     2.1          calculate wet-bulb temperature and moisture
2787 !                  for environmental air in *cuadjtq*
2788 !                  -------------------------------------------
2789         is=0
2790         do jl=1,klon
2791           ztenwb(jl,jk)=ptenh(jl,jk)
2792           zqenwb(jl,jk)=pqenh(jl,jk)
2793           zph(jl)=paph(jl,jk)
2794           llo2(jl)=ldcum(jl).and.prfl(jl).gt.0..and..not.lddraf(jl).and.   &
2795      &     (jk.lt.kcbot(jl).and.jk.gt.kctop(jl)).and. jk.ge.ikhsmin(jl)
2796           if(llo2(jl))then
2797             is=is+1
2798           endif
2799         enddo
2800         if(is.eq.0) cycle
2801
2802         ik=jk
2803         icall=2
2804         call cuadjtqn                                           &
2805      &   ( klon, klev, ik, zph, ztenwb, zqenwb, llo2, icall)
2806
2807 !     2.2          do mixing of cumulus and environmental air
2808 !                  and check for negative buoyancy.
2809 !                  then set values for downdraft at lfs.
2810 !                  ----------------------------------------
2811         do jl=1,klon
2812           if(llo2(jl)) then
2813             zttest=0.5*(ptu(jl,jk)+ztenwb(jl,jk))
2814             zqtest=0.5*(pqu(jl,jk)+zqenwb(jl,jk))
2815             zbuo=zttest*(1.+vtmpc1  *zqtest)-                    &
2816      &       ptenh(jl,jk)*(1.+vtmpc1  *pqenh(jl,jk))
2817             zcond(jl)=pqenh(jl,jk)-zqenwb(jl,jk)
2818             zmftop=-cmfdeps*pmfub(jl)
2819             if(zbuo.lt.0..and.prfl(jl).gt.10.*zmftop*zcond(jl)) then
2820               kdtop(jl)=jk
2821               lddraf(jl)=.true.
2822               ptd(jl,jk)=zttest
2823               pqd(jl,jk)=zqtest
2824               pmfd(jl,jk)=zmftop
2825               pmfds(jl,jk)=pmfd(jl,jk)*(cpd*ptd(jl,jk)+pgeoh(jl,jk))
2826               pmfdq(jl,jk)=pmfd(jl,jk)*pqd(jl,jk)
2827               pdmfdp(jl,jk-1)=-0.5*pmfd(jl,jk)*zcond(jl)
2828               prfl(jl)=prfl(jl)+pdmfdp(jl,jk-1)
2829             endif
2830           endif
2831         enddo
2832
2833       enddo
2834
2835       return
2836       end subroutine cudlfsn
2837
2838 !---------------------------------------------------------
2839 !  level 3 souroutines
2840 !--------------------------------------------------------
2841 !**********************************************
2842 !       subroutine cuddrafn
2843 !**********************************************
2844        subroutine cuddrafn                                 &
2845      &   ( klon,     klev,    lddraf                      &
2846      &   , ptenh,    pqenh,    puen,     pven             &
2847      &   , pgeo,     pgeoh,    paph,     prfl             &
2848      &   , ptd,      pqd,      pud,      pvd,      pmfu   &
2849      &   , pmfd,     pmfds,    pmfdq,    pdmfdp,   pmfdde_rate     )
2850
2851 !          this routine calculates cumulus downdraft descent
2852
2853 !          m.tiedtke         e.c.m.w.f.    12/86 modif. 12/89
2854
2855 !          purpose.
2856 !          --------
2857 !          to produce the vertical profiles for cumulus downdrafts
2858 !          (i.e. t,q,u and v and fluxes)
2859
2860 !          interface
2861 !          ---------
2862
2863 !          this routine is called from *cumastr*.
2864 !          input is t,q,p,phi,u,v at half levels.
2865 !          it returns fluxes of s,q and evaporation rate
2866 !          and u,v at levels where downdraft occurs
2867
2868 !          method.
2869 !          --------
2870 !          calculate moist descent for entraining/detraining plume by
2871 !          a) moving air dry-adiabatically to next level below and
2872 !          b) correcting for evaporation to obtain saturated state.
2873
2874 !     parameter     description                                   units
2875 !     ---------     -----------                                   -----
2876 !     input parameters (integer):
2877
2878 !    *klon*         number of grid points per packet
2879 !    *klev*         number of levels
2880
2881 !    input parameters (logical):
2882
2883 !    *lddraf*       .true. if downdrafts exist
2884
2885 !    input parameters (real):
2886
2887 !    *ptenh*        env. temperature (t+1) on half levels          k
2888 !    *pqenh*        env. spec. humidity (t+1) on half levels     kg/kg
2889 !    *puen*         provisional environment u-velocity (t+1)      m/s
2890 !    *pven*         provisional environment v-velocity (t+1)      m/s
2891 !    *pgeo*         geopotential                                  m2/s2
2892 !    *pgeoh*        geopotential on half levels                  m2/s2
2893 !    *paph*         provisional pressure on half levels           pa
2894 !    *pmfu*         massflux updrafts                           kg/(m2*s)
2895
2896 !    updated parameters (real):
2897
2898 !    *prfl*         precipitation rate                           kg/(m2*s)
2899
2900 !    output parameters (real):
2901
2902 !    *ptd*          temperature in downdrafts                      k
2903 !    *pqd*          spec. humidity in downdrafts                 kg/kg
2904 !    *pud*          u-velocity in downdrafts                      m/s
2905 !    *pvd*          v-velocity in downdrafts                      m/s
2906 !    *pmfd*         massflux in downdrafts                       kg/(m2*s)
2907 !    *pmfds*        flux of dry static energy in downdrafts       j/(m2*s)
2908 !    *pmfdq*        flux of spec. humidity in downdrafts         kg/(m2*s)
2909 !    *pdmfdp*       flux difference of precip. in downdrafts     kg/(m2*s)
2910
2911 !          externals
2912 !          ---------
2913 !          *cuadjtq* for adjusting t and q due to evaporation in
2914 !          saturated descent
2915 !----------------------------------------------------------------------
2916       implicit none
2917
2918       integer  klev,klon
2919       real     ptenh(klon,klev),       pqenh(klon,klev),   &
2920      &         puen(klon,klev),        pven(klon,klev),    &
2921      &         pgeoh(klon,klev+1),     paph(klon,klev+1),  &
2922      &         pgeo(klon,klev),        pmfu(klon,klev)
2923
2924       real     ptd(klon,klev),         pqd(klon,klev),     &
2925      &         pud(klon,klev),         pvd(klon,klev),     &
2926      &         pmfd(klon,klev),        pmfds(klon,klev),   &
2927      &         pmfdq(klon,klev),       pdmfdp(klon,klev),  &
2928      &         prfl(klon)
2929       real     pmfdde_rate(klon,klev)
2930       logical  lddraf(klon)
2931
2932       real     zdmfen(klon),           zdmfde(klon),       &
2933      &         zcond(klon),            zoentr(klon),       &
2934      &         zbuoy(klon)
2935       real     zph(klon)
2936       logical  llo2(klon)
2937       logical  llo1
2938 ! local variables
2939       integer  jl,jk
2940       integer  is,ik,icall,ike, itopde(klon)
2941       real     zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp
2942       real     zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk
2943
2944 !----------------------------------------------------------------------
2945 !     1.           calculate moist descent for cumulus downdraft by
2946 !                     (a) calculating entrainment/detrainment rates,
2947 !                         including organized entrainment dependent on
2948 !                         negative buoyancy and assuming
2949 !                         linear decrease of massflux in pbl
2950 !                     (b) doing moist descent - evaporative cooling
2951 !                         and moistening is calculated in *cuadjtq*
2952 !                     (c) checking for negative buoyancy and
2953 !                         specifying final t,q,u,v and downward fluxes
2954 !                    -------------------------------------------------
2955       do jl=1,klon
2956         zoentr(jl)=0.
2957         zbuoy(jl)=0.
2958         zdmfen(jl)=0.
2959         zdmfde(jl)=0.
2960       enddo
2961
2962       do jk=klev,1,-1
2963        do jl=1,klon
2964          pmfdde_rate(jl,jk) = 0.
2965          if((paph(jl,klev+1)-paph(jl,jk)).lt. 60.e2) itopde(jl) = jk
2966        end do
2967       end do
2968
2969       do jk=3,klev
2970         is=0
2971         do jl=1,klon
2972           zph(jl)=paph(jl,jk)
2973           llo2(jl)=lddraf(jl).and.pmfd(jl,jk-1).lt.0.
2974           if(llo2(jl)) then
2975             is=is+1
2976           endif
2977         enddo
2978         if(is.eq.0) cycle
2979
2980         do jl=1,klon
2981           if(llo2(jl)) then
2982             zentr = entrdd*pmfd(jl,jk-1)*(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
2983             zdmfen(jl)=zentr
2984             zdmfde(jl)=zentr
2985           endif
2986         enddo
2987
2988         do jl=1,klon
2989           if(llo2(jl)) then
2990           if(jk.gt.itopde(jl)) then
2991             zdmfen(jl)=0.
2992             zdmfde(jl)=pmfd(jl,itopde(jl))*                    &
2993      &       (paph(jl,jk)-paph(jl,jk-1))/                  &
2994      &       (paph(jl,klev+1)-paph(jl,itopde(jl)))
2995           endif
2996           endif
2997         enddo
2998
2999         do jl=1,klon
3000           if(llo2(jl)) then
3001           if(jk.le.itopde(jl)) then
3002             zdz=-(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg
3003             zzentr=zoentr(jl)*zdz*pmfd(jl,jk-1)
3004             zdmfen(jl)=zdmfen(jl)+zzentr
3005             zdmfen(jl)=max(zdmfen(jl),0.3*pmfd(jl,jk-1))
3006             zdmfen(jl)=max(zdmfen(jl),-0.75*pmfu(jl,jk)-   &
3007    &         (pmfd(jl,jk-1)-zdmfde(jl)))
3008             zdmfen(jl)=min(zdmfen(jl),0.)
3009           endif
3010           endif
3011         enddo
3012
3013         do jl=1,klon
3014           if(llo2(jl)) then
3015             pmfd(jl,jk)=pmfd(jl,jk-1)+zdmfen(jl)-zdmfde(jl)
3016             zseen=(cpd*ptenh(jl,jk-1)+pgeoh(jl,jk-1))*zdmfen(jl)
3017             zqeen=pqenh(jl,jk-1)*zdmfen(jl)
3018             zsdde=(cpd*ptd(jl,jk-1)+pgeoh(jl,jk-1))*zdmfde(jl)
3019             zqdde=pqd(jl,jk-1)*zdmfde(jl)
3020             zmfdsk=pmfds(jl,jk-1)+zseen-zsdde
3021             zmfdqk=pmfdq(jl,jk-1)+zqeen-zqdde
3022             pqd(jl,jk)=zmfdqk*(1./min(-cmfcmin,pmfd(jl,jk)))
3023             ptd(jl,jk)=(zmfdsk*(1./min(-cmfcmin,pmfd(jl,jk)))-&
3024      &                  pgeoh(jl,jk))*rcpd
3025             ptd(jl,jk)=min(400.,ptd(jl,jk))
3026             ptd(jl,jk)=max(100.,ptd(jl,jk))
3027             zcond(jl)=pqd(jl,jk)
3028           endif
3029         enddo
3030
3031         ik=jk
3032         icall=2
3033         call cuadjtqn(klon, klev, ik, zph, ptd, pqd, llo2, icall )
3034
3035         do jl=1,klon
3036           if(llo2(jl)) then
3037             zcond(jl)=zcond(jl)-pqd(jl,jk)
3038             zbuo=ptd(jl,jk)*(1.+vtmpc1  *pqd(jl,jk))-          &
3039      &      ptenh(jl,jk)*(1.+vtmpc1  *pqenh(jl,jk))
3040             if(prfl(jl).gt.0..and.pmfu(jl,jk).gt.0.) then
3041               zrain=prfl(jl)/pmfu(jl,jk)
3042               zbuo=zbuo-ptd(jl,jk)*zrain
3043             endif
3044             if(zbuo.ge.0 .or. prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then
3045               pmfd(jl,jk)=0.
3046               zbuo=0.
3047             endif
3048             pmfds(jl,jk)=(cpd*ptd(jl,jk)+pgeoh(jl,jk))*pmfd(jl,jk)
3049             pmfdq(jl,jk)=pqd(jl,jk)*pmfd(jl,jk)
3050             zdmfdp=-pmfd(jl,jk)*zcond(jl)
3051             pdmfdp(jl,jk-1)=zdmfdp
3052             prfl(jl)=prfl(jl)+zdmfdp
3053
3054 ! compute organized entrainment for use at next level
3055             zbuoyz=zbuo/ptenh(jl,jk)
3056             zbuoyz=min(zbuoyz,0.0)
3057             zdz=-(pgeo(jl,jk-1)-pgeo(jl,jk))
3058             zbuoy(jl)=zbuoy(jl)+zbuoyz*zdz
3059             zoentr(jl)=g*zbuoyz*0.5/(1.+zbuoy(jl))
3060             pmfdde_rate(jl,jk) = -zdmfde(jl)
3061           endif
3062         enddo
3063
3064       enddo
3065
3066       return
3067       end subroutine cuddrafn
3068 !---------------------------------------------------------
3069 !  level 3 souroutines
3070 !--------------------------------------------------------
3071        subroutine cuflxn         &
3072      &  (  klon,     klev,     ztmst &
3073      &  ,  pten,     pqen,     pqsen,    ptenh,    pqenh   &
3074      &  ,  paph,     pap,      pgeoh,    lndj,   ldcum   &
3075      &  ,  kcbot,    kctop,    kdtop,    ktopm2            &
3076      &  ,  ktype,    lddraf                                &
3077      &  ,  pmfu,     pmfd,     pmfus,    pmfds             &
3078      &  ,  pmfuq,    pmfdq,    pmful,    plude             &
3079      &  ,  pdmfup,   pdmfdp,   pdpmel,   plglac            &
3080      &  ,  prain,    pmfdde_rate, pmflxr, pmflxs )
3081
3082 !          m.tiedtke         e.c.m.w.f.     7/86 modif. 12/89
3083
3084 !          purpose
3085 !          -------
3086
3087 !          this routine does the final calculation of convective
3088 !          fluxes in the cloud layer and in the subcloud layer
3089
3090 !          interface
3091 !          ---------
3092 !          this routine is called from *cumastr*.
3093
3094
3095 !     parameter     description                                   units
3096 !     ---------     -----------                                   -----
3097 !     input parameters (integer):
3098
3099 !    *klon*         number of grid points per packet
3100 !    *klev*         number of levels
3101 !    *kcbot*        cloud base level
3102 !    *kctop*        cloud top level
3103 !    *kdtop*        top level of downdrafts
3104
3105 !    input parameters (logical):
3106
3107 !    *lndj*       land sea mask (1 for land)
3108 !    *ldcum*        flag: .true. for convective points
3109
3110 !    input parameters (real):
3111
3112 !    *ptsphy*       time step for the physics                       s
3113 !    *pten*         provisional environment temperature (t+1)       k
3114 !    *pqen*         provisional environment spec. humidity (t+1)  kg/kg
3115 !    *pqsen*        environment spec. saturation humidity (t+1)   kg/kg
3116 !    *ptenh*        env. temperature (t+1) on half levels           k
3117 !    *pqenh*        env. spec. humidity (t+1) on half levels      kg/kg
3118 !    *paph*         provisional pressure on half levels            pa
3119 !    *pap*          provisional pressure on full levels            pa
3120 !    *pgeoh*        geopotential on half levels                   m2/s2
3121
3122 !    updated parameters (integer):
3123
3124 !    *ktype*        set to zero if ldcum=.false.
3125
3126 !    updated parameters (logical):
3127
3128 !    *lddraf*       set to .false. if ldcum=.false. or kdtop<kctop
3129
3130 !    updated parameters (real):
3131
3132 !    *pmfu*         massflux in updrafts                          kg/(m2*s)
3133 !    *pmfd*         massflux in downdrafts                        kg/(m2*s)
3134 !    *pmfus*        flux of dry static energy in updrafts          j/(m2*s)
3135 !    *pmfds*        flux of dry static energy in downdrafts        j/(m2*s)
3136 !    *pmfuq*        flux of spec. humidity in updrafts            kg/(m2*s)
3137 !    *pmfdq*        flux of spec. humidity in downdrafts          kg/(m2*s)
3138 !    *pmful*        flux of liquid water in updrafts              kg/(m2*s)
3139 !    *plude*        detrained liquid water                        kg/(m3*s)
3140 !    *pdmfup*       flux difference of precip. in updrafts        kg/(m2*s)
3141 !    *pdmfdp*       flux difference of precip. in downdrafts      kg/(m2*s)
3142
3143 !    output parameters (real):
3144
3145 !    *pdpmel*       change in precip.-fluxes due to melting       kg/(m2*s)
3146 !    *plglac*       flux of frozen cloud water in updrafts        kg/(m2*s)
3147 !    *pmflxr*       convective rain flux                          kg/(m2*s)
3148 !    *pmflxs*       convective snow flux                          kg/(m2*s)
3149 !    *prain*        total precip. produced in conv. updrafts      kg/(m2*s)
3150 !                   (no evaporation in downdrafts)
3151
3152 !          externals
3153 !          ---------
3154 !          none
3155 !----------------------------------------------------------------------
3156       implicit none
3157
3158       integer  klev,klon,ktopm2
3159       real     pten(klon,klev),        ptenh(klon,klev),          &
3160      &         pqen(klon,klev),        pqsen(klon,klev),          &
3161      &         pqenh(klon,klev),       pap(klon,klev),            &
3162      &         paph(klon,klev+1),      pgeoh(klon,klev+1),        &
3163      &         plglac(klon,klev)
3164
3165       real     pmfu(klon,klev),        pmfd(klon,klev),           &
3166      &         pmfus(klon,klev),       pmfds(klon,klev),          &
3167      &         pmfuq(klon,klev),       pmfdq(klon,klev),          &
3168      &         pdmfup(klon,klev),      pdmfdp(klon,klev),         &
3169      &         pdpmel(klon,klev),      prain(klon),               &
3170      &         pmful(klon,klev),       plude(klon,klev),          &
3171      &         pmflxr(klon,klev+1),    pmflxs(klon,klev+1)
3172       real     pmfdde_rate(klon,klev)
3173       integer  kcbot(klon),            kctop(klon),               &
3174      &         kdtop(klon),            ktype(klon)
3175       logical  lddraf(klon),                                &
3176      &         ldcum(klon)
3177       integer  lndj(klon)
3178 ! local variables
3179       integer  jl,jk
3180       integer  is,ik,icall,ike,ikb
3181       real     ztmst,ztaumel,zcons1a,zcons1,zcons2,zcucov,zcpecons
3182       real     zalfaw,zrfl,zdrfl1,zrnew,zrmin,zrfln,zdrfl,zdenom
3183       real     zpdr,zpds,zzp,zfac,zsnmlt
3184       real     rhevap(klon)
3185       integer  idbas(klon)
3186       logical  llddraf
3187 !--------------------------------------------------------------------
3188 !*             specify constants
3189
3190       ztaumel=18000.
3191       zcons1a=cpd/(alf*g*ztaumel)
3192       zcons2=3./(g*ztmst)
3193       zcucov=0.05
3194       zcpecons=5.44e-4/g
3195
3196 !*    1.0          determine final convective fluxes
3197 !                  ---------------------------------
3198       do jl=1,klon
3199         prain(jl)=0.
3200         if(.not.ldcum(jl).or.kdtop(jl).lt.kctop(jl)) lddraf(jl)=.false.
3201         if(.not.ldcum(jl)) ktype(jl)=0
3202         idbas(jl) = klev
3203         if(lndj(jl) .eq. 1) then
3204           rhevap(jl) = 0.7
3205         else
3206           rhevap(jl) = 0.9
3207         end if
3208       enddo
3209
3210       ktopm2= 2
3211       do jk=ktopm2,klev
3212         ikb = min(jk+1,klev)
3213         do jl=1,klon
3214           pmflxr(jl,jk) = 0.
3215           pmflxs(jl,jk) = 0.
3216           pdpmel(jl,jk) = 0.
3217           if(ldcum(jl).and.jk.ge.kctop(jl)) then
3218             pmfus(jl,jk)=pmfus(jl,jk)-pmfu(jl,jk)*           &
3219      &       (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
3220             pmfuq(jl,jk)=pmfuq(jl,jk)-pmfu(jl,jk)*pqenh(jl,jk)
3221             plglac(jl,jk)=pmfu(jl,jk)*plglac(jl,jk)
3222             llddraf = lddraf(jl) .and. jk >= kdtop(jl)
3223             if ( llddraf .and.jk.ge.kdtop(jl)) then
3224               pmfds(jl,jk) = pmfds(jl,jk)-pmfd(jl,jk) * &
3225                            (cpd*ptenh(jl,jk)+pgeoh(jl,jk))
3226               pmfdq(jl,jk) = pmfdq(jl,jk)-pmfd(jl,jk)*pqenh(jl,jk)
3227             else
3228               pmfd(jl,jk) = 0.
3229               pmfds(jl,jk) = 0.
3230               pmfdq(jl,jk) = 0.
3231               pdmfdp(jl,jk-1) = 0.
3232             end if
3233             if ( llddraf .and. pmfd(jl,jk) < 0. .and. &
3234                abs(pmfd(jl,ikb)) < 1.e-20 ) then
3235                idbas(jl) = jk
3236             end if
3237           else
3238             pmfu(jl,jk)=0.
3239             pmfd(jl,jk)=0.
3240             pmfus(jl,jk)=0.
3241             pmfds(jl,jk)=0.
3242             pmfuq(jl,jk)=0.
3243             pmfdq(jl,jk)=0.
3244             pmful(jl,jk)=0.
3245             plglac(jl,jk)=0.
3246             pdmfup(jl,jk-1)=0.
3247             pdmfdp(jl,jk-1)=0.
3248             plude(jl,jk-1)=0.
3249           endif
3250         enddo
3251       enddo
3252
3253       do jl=1,klon
3254         pmflxr(jl,klev+1) = 0.
3255         pmflxs(jl,klev+1) = 0.
3256       end do
3257       do jl=1,klon
3258         if(ldcum(jl)) then
3259           ikb=kcbot(jl)
3260           ik=ikb+1
3261           zzp=((paph(jl,klev+1)-paph(jl,ik))/                  &
3262      &     (paph(jl,klev+1)-paph(jl,ikb)))
3263           if(ktype(jl).eq.3) then
3264             zzp=zzp**2
3265           endif
3266           pmfu(jl,ik)=pmfu(jl,ikb)*zzp
3267           pmfus(jl,ik)=(pmfus(jl,ikb)-                         &
3268      &         foelhm(ptenh(jl,ikb))*pmful(jl,ikb))*zzp
3269           pmfuq(jl,ik)=(pmfuq(jl,ikb)+pmful(jl,ikb))*zzp
3270           pmful(jl,ik)=0.
3271         endif
3272       enddo
3273
3274       do jk=ktopm2,klev
3275         do jl=1,klon
3276           if(ldcum(jl).and.jk.gt.kcbot(jl)+1) then
3277             ikb=kcbot(jl)+1
3278             zzp=((paph(jl,klev+1)-paph(jl,jk))/                &
3279      &       (paph(jl,klev+1)-paph(jl,ikb)))
3280             if(ktype(jl).eq.3) then
3281               zzp=zzp**2
3282             endif
3283             pmfu(jl,jk)=pmfu(jl,ikb)*zzp
3284             pmfus(jl,jk)=pmfus(jl,ikb)*zzp
3285             pmfuq(jl,jk)=pmfuq(jl,ikb)*zzp
3286             pmful(jl,jk)=0.
3287           endif
3288           ik = idbas(jl)
3289           llddraf = lddraf(jl) .and. jk > ik .and. ik < klev
3290           if ( llddraf .and. ik == kcbot(jl)+1 ) then
3291            zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ik)))
3292            if ( ktype(jl) == 3 ) zzp = zzp*zzp
3293             pmfd(jl,jk) = pmfd(jl,ik)*zzp
3294             pmfds(jl,jk) = pmfds(jl,ik)*zzp
3295             pmfdq(jl,jk) = pmfdq(jl,ik)*zzp
3296             pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
3297            else if ( llddraf .and. ik /= kcbot(jl)+1 .and. jk == ik+1 ) then
3298             pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk))
3299            end if
3300         enddo
3301       enddo
3302 !*    2.            calculate rain/snow fall rates
3303 !*                  calculate melting of snow
3304 !*                  calculate evaporation of precip
3305 !                   -------------------------------
3306
3307         do jk=ktopm2,klev
3308         do jl=1,klon
3309           if(ldcum(jl) .and. jk >=kctop(jl)-1 ) then
3310             prain(jl)=prain(jl)+pdmfup(jl,jk)
3311             if(pmflxs(jl,jk).gt.0..and.pten(jl,jk).gt.tmelt) then
3312               zcons1=zcons1a*(1.+0.5*(pten(jl,jk)-tmelt))
3313               zfac=zcons1*(paph(jl,jk+1)-paph(jl,jk))
3314               zsnmlt=min(pmflxs(jl,jk),zfac*(pten(jl,jk)-tmelt))
3315               pdpmel(jl,jk)=zsnmlt
3316               pqsen(jl,jk)=foeewm(pten(jl,jk)-zsnmlt/zfac)/pap(jl,jk)
3317             endif
3318             zalfaw=foealfa(pten(jl,jk))
3319           !
3320           ! No liquid precipitation above melting level
3321           !
3322             if ( pten(jl,jk) < tmelt .and. zalfaw > 0. ) then
3323               plglac(jl,jk) = plglac(jl,jk)+zalfaw*(pdmfup(jl,jk)+pdmfdp(jl,jk))
3324               zalfaw = 0.
3325             end if
3326             pmflxr(jl,jk+1)=pmflxr(jl,jk)+zalfaw*               &
3327      &       (pdmfup(jl,jk)+pdmfdp(jl,jk))+pdpmel(jl,jk)
3328             pmflxs(jl,jk+1)=pmflxs(jl,jk)+(1.-zalfaw)*          &
3329      &       (pdmfup(jl,jk)+pdmfdp(jl,jk))-pdpmel(jl,jk)
3330             if(pmflxr(jl,jk+1)+pmflxs(jl,jk+1).lt.0.0) then
3331               pdmfdp(jl,jk)=-(pmflxr(jl,jk)+pmflxs(jl,jk)+pdmfup(jl,jk))
3332               pmflxr(jl,jk+1)=0.0
3333               pmflxs(jl,jk+1)=0.0
3334               pdpmel(jl,jk)  =0.0
3335             else if ( pmflxr(jl,jk+1) < 0. ) then
3336               pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
3337               pmflxr(jl,jk+1) = 0.
3338             else if ( pmflxs(jl,jk+1) < 0. ) then
3339               pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
3340               pmflxs(jl,jk+1) = 0.
3341             end if
3342           endif
3343         enddo
3344       enddo
3345       do jk=ktopm2,klev
3346         do jl=1,klon
3347           if(ldcum(jl).and.jk.ge.kcbot(jl)) then
3348             zrfl=pmflxr(jl,jk)+pmflxs(jl,jk)
3349             if(zrfl.gt.1.e-20) then
3350               zdrfl1=zcpecons*max(0.,pqsen(jl,jk)-pqen(jl,jk))*zcucov* &
3351      &         (sqrt(paph(jl,jk)/paph(jl,klev+1))/5.09e-3*             &
3352      &         zrfl/zcucov)**0.5777*                                   &
3353      &         (paph(jl,jk+1)-paph(jl,jk))
3354               zrnew=zrfl-zdrfl1
3355               zrmin=zrfl-zcucov*max(0.,rhevap(jl)*pqsen(jl,jk) &
3356      &              -pqen(jl,jk)) *zcons2*(paph(jl,jk+1)-paph(jl,jk))
3357               zrnew=max(zrnew,zrmin)
3358               zrfln=max(zrnew,0.)
3359               zdrfl=min(0.,zrfln-zrfl)
3360               zdenom=1./max(1.e-20,pmflxr(jl,jk)+pmflxs(jl,jk))
3361               zalfaw=foealfa(pten(jl,jk))
3362               if ( pten(jl,jk) < tmelt ) zalfaw = 0.
3363               zpdr=zalfaw*pdmfdp(jl,jk)
3364               zpds=(1.-zalfaw)*pdmfdp(jl,jk)
3365               pmflxr(jl,jk+1)=pmflxr(jl,jk)+zpdr         &
3366      &         +pdpmel(jl,jk)+zdrfl*pmflxr(jl,jk)*zdenom
3367               pmflxs(jl,jk+1)=pmflxs(jl,jk)+zpds         &
3368      &         -pdpmel(jl,jk)+zdrfl*pmflxs(jl,jk)*zdenom
3369               pdmfup(jl,jk)=pdmfup(jl,jk)+zdrfl
3370               if ( pmflxr(jl,jk+1)+pmflxs(jl,jk+1) < 0. ) then
3371                 pdmfup(jl,jk) = pdmfup(jl,jk)-(pmflxr(jl,jk+1)+pmflxs(jl,jk+1))
3372                 pmflxr(jl,jk+1) = 0.
3373                 pmflxs(jl,jk+1) = 0.
3374                 pdpmel(jl,jk)   = 0.
3375               else if ( pmflxr(jl,jk+1) < 0. ) then
3376                 pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1)
3377                 pmflxr(jl,jk+1) = 0.
3378               else if ( pmflxs(jl,jk+1) < 0. ) then
3379                 pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1)
3380                 pmflxs(jl,jk+1) = 0.
3381               end if
3382             else
3383               pmflxr(jl,jk+1)=0.0
3384               pmflxs(jl,jk+1)=0.0
3385               pdmfdp(jl,jk)=0.0
3386               pdpmel(jl,jk)=0.0
3387             endif
3388           endif
3389         enddo
3390       enddo
3391
3392       return
3393       end subroutine cuflxn
3394 !---------------------------------------------------------
3395 !  level 3 souroutines
3396 !--------------------------------------------------------
3397      subroutine cudtdqn(klon,klev,ktopm2,kctop,kdtop,ldcum, &
3398                      lddraf,ztmst,paph,pgeoh,pgeo,pten,ptenh,pqen,  &
3399                      pqenh,pqsen,plglac,plude,pmfu,pmfd,pmfus,pmfds, &
3400                      pmfuq,pmfdq,pmful,pdmfup,pdmfdp,pdpmel,ptent,ptenq,pcte)
3401     implicit none
3402     integer  klon,klev,ktopm2
3403     integer  kctop(klon),  kdtop(klon)
3404     logical  ldcum(klon),  lddraf(klon)
3405     real     ztmst
3406     real     paph(klon,klev+1), pgeoh(klon,klev+1)
3407     real     pgeo(klon,klev),   pten(klon,klev), &
3408              pqen(klon,klev),   ptenh(klon,klev),&
3409              pqenh(klon,klev),  pqsen(klon,klev),&
3410              plglac(klon,klev), plude(klon,klev)
3411     real     pmfu(klon,klev),   pmfd(klon,klev),&
3412              pmfus(klon,klev),  pmfds(klon,klev),&
3413              pmfuq(klon,klev),  pmfdq(klon,klev),&
3414              pmful(klon,klev),  pdmfup(klon,klev),&
3415              pdpmel(klon,klev), pdmfdp(klon,klev)
3416     real     ptent(klon,klev),  ptenq(klon,klev)
3417     real     pcte(klon,klev)
3418
3419 ! local variables
3420     integer  jk , ik , jl
3421     real     zalv , zzp
3422     real     zdtdt(klon,klev) , zdqdt(klon,klev) , zdp(klon,klev)
3423     !*    1.0          SETUP AND INITIALIZATIONS
3424     ! -------------------------
3425     do jk = 1 , klev
3426       do jl = 1, klon
3427         if ( ldcum(jl) ) then
3428           zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
3429         end if
3430       end do
3431     end do
3432     !-----------------------------------------------------------------------
3433     !*    2.0          COMPUTE TENDENCIES
3434     ! ------------------
3435     do jk = ktopm2 , klev
3436       if ( jk < klev ) then
3437         do jl = 1,klon
3438           if ( ldcum(jl) ) then
3439             zalv = foelhm(pten(jl,jk))
3440             zdtdt(jl,jk) = zdp(jl,jk)*rcpd * &
3441               (pmfus(jl,jk+1)-pmfus(jl,jk)+pmfds(jl,jk+1) - &
3442                pmfds(jl,jk)+alf*plglac(jl,jk)-alf*pdpmel(jl,jk) - &
3443                zalv*(pmful(jl,jk+1)-pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk)))
3444             zdqdt(jl,jk) = zdp(jl,jk)*(pmfuq(jl,jk+1) - &
3445               pmfuq(jl,jk)+pmfdq(jl,jk+1)-pmfdq(jl,jk)+pmful(jl,jk+1) - &
3446               pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk))
3447           end if
3448         end do
3449       else
3450         do jl = 1,klon
3451           if ( ldcum(jl) ) then
3452             zalv = foelhm(pten(jl,jk))
3453             zdtdt(jl,jk) = -zdp(jl,jk)*rcpd * &
3454               (pmfus(jl,jk)+pmfds(jl,jk)+alf*pdpmel(jl,jk) - &
3455                zalv*(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)+plude(jl,jk)))
3456             zdqdt(jl,jk) = -zdp(jl,jk)*(pmfuq(jl,jk) + plude(jl,jk) + &
3457               pmfdq(jl,jk)+(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)))
3458           end if
3459         end do
3460       end if
3461     end do
3462   !---------------------------------------------------------------
3463   !*  3.0          UPDATE TENDENCIES
3464   !   -----------------
3465     do jk = ktopm2 , klev
3466      do jl = 1, klon
3467        if ( ldcum(jl) ) then
3468          ptent(jl,jk) = ptent(jl,jk) + zdtdt(jl,jk)
3469          ptenq(jl,jk) = ptenq(jl,jk) + zdqdt(jl,jk)
3470          pcte(jl,jk)  = zdp(jl,jk)*plude(jl,jk)
3471        end if
3472      end do
3473     end do
3474
3475     return
3476   end subroutine cudtdqn
3477 !---------------------------------------------------------
3478 !  level 3 souroutines
3479 !--------------------------------------------------------
3480     subroutine cududvn(klon,klev,ktopm2,ktype,kcbot,kctop,ldcum,  &
3481                     ztmst,paph,puen,pven,pmfu,pmfd,puu,pud,pvu,pvd,ptenu, &
3482                     ptenv)
3483     implicit none
3484     integer  klon,klev,ktopm2
3485     integer  ktype(klon), kcbot(klon), kctop(klon)
3486     logical  ldcum(klon)
3487     real     ztmst
3488     real     paph(klon,klev+1)
3489     real     puen(klon,klev),  pven(klon,klev),&
3490              pmfu(klon,klev),  pmfd(klon,klev),&
3491              puu(klon,klev),   pud(klon,klev),&
3492              pvu(klon,klev),   pvd(klon,klev)
3493     real     ptenu(klon,klev), ptenv(klon,klev)
3494
3495 !local variables
3496     real     zuen(klon,klev) , zven(klon,klev) , zmfuu(klon,klev), &
3497              zmfdu(klon,klev), zmfuv(klon,klev), zmfdv(klon,klev)
3498
3499     integer  ik , ikb , jk , jl
3500     real     zzp, zdtdt
3501
3502     real     zdudt(klon,klev), zdvdt(klon,klev), zdp(klon,klev)
3503 !
3504     do jk = 1 , klev
3505       do jl = 1, klon
3506         if ( ldcum(jl) ) then
3507           zuen(jl,jk) = puen(jl,jk)
3508           zven(jl,jk) = pven(jl,jk)
3509           zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk))
3510         end if
3511       end do
3512     end do
3513 !----------------------------------------------------------------------
3514 !*    1.0          CALCULATE FLUXES AND UPDATE U AND V TENDENCIES
3515 ! ----------------------------------------------
3516     do jk = ktopm2 , klev
3517       ik = jk - 1
3518       do jl = 1,klon
3519         if ( ldcum(jl) ) then
3520           zmfuu(jl,jk) = pmfu(jl,jk)*(puu(jl,jk)-zuen(jl,ik))
3521           zmfuv(jl,jk) = pmfu(jl,jk)*(pvu(jl,jk)-zven(jl,ik))
3522           zmfdu(jl,jk) = pmfd(jl,jk)*(pud(jl,jk)-zuen(jl,ik))
3523           zmfdv(jl,jk) = pmfd(jl,jk)*(pvd(jl,jk)-zven(jl,ik))
3524         end if
3525       end do
3526     end do
3527     ! linear fluxes below cloud
3528       do jk = ktopm2 , klev
3529         do jl = 1, klon
3530           if ( ldcum(jl) .and. jk > kcbot(jl) ) then
3531             ikb = kcbot(jl)
3532             zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb)))
3533             if ( ktype(jl) == 3 ) zzp = zzp*zzp
3534             zmfuu(jl,jk) = zmfuu(jl,ikb)*zzp
3535             zmfuv(jl,jk) = zmfuv(jl,ikb)*zzp
3536             zmfdu(jl,jk) = zmfdu(jl,ikb)*zzp
3537             zmfdv(jl,jk) = zmfdv(jl,ikb)*zzp
3538           end if
3539         end do
3540       end do
3541 !----------------------------------------------------------------------
3542 !*    2.0          COMPUTE TENDENCIES
3543 ! ------------------
3544     do jk = ktopm2 , klev
3545       if ( jk < klev ) then
3546         ik = jk + 1
3547         do jl = 1,klon
3548           if ( ldcum(jl) ) then
3549             zdudt(jl,jk) = zdp(jl,jk) * &
3550                           (zmfuu(jl,ik)-zmfuu(jl,jk)+zmfdu(jl,ik)-zmfdu(jl,jk))
3551             zdvdt(jl,jk) = zdp(jl,jk) * &
3552                           (zmfuv(jl,ik)-zmfuv(jl,jk)+zmfdv(jl,ik)-zmfdv(jl,jk))
3553           end if
3554         end do
3555       else
3556         do jl = 1,klon
3557           if ( ldcum(jl) ) then
3558             zdudt(jl,jk) = -zdp(jl,jk)*(zmfuu(jl,jk)+zmfdu(jl,jk))
3559             zdvdt(jl,jk) = -zdp(jl,jk)*(zmfuv(jl,jk)+zmfdv(jl,jk))
3560           end if
3561         end do
3562       end if
3563     end do
3564 !---------------------------------------------------------------------
3565 !*  3.0        UPDATE TENDENCIES
3566 !   -----------------
3567     do jk = ktopm2 , klev
3568       do jl = 1, klon
3569         if ( ldcum(jl) ) then
3570           ptenu(jl,jk) = ptenu(jl,jk) + zdudt(jl,jk)
3571           ptenv(jl,jk) = ptenv(jl,jk) + zdvdt(jl,jk)
3572         end if
3573       end do
3574     end do
3575 !----------------------------------------------------------------------
3576     return
3577     end subroutine cududvn
3578 !---------------------------------------------------------
3579 !  level 3 souroutines
3580 !--------------------------------------------------------
3581       subroutine cuadjtqn                 &
3582      &    (klon, klev, kk, psp, pt, pq, ldflag,  kcall)
3583 !          m.tiedtke         e.c.m.w.f.     12/89
3584 !          purpose.
3585 !          --------
3586 !          to produce t,q and l values for cloud ascent
3587
3588 !          interface
3589 !          ---------
3590 !          this routine is called from subroutines:
3591 !              *cond*     (t and q at condensation level)
3592 !              *cubase*   (t and q at condensation level)
3593 !              *cuasc*    (t and q at cloud levels)
3594 !              *cuini*    (environmental t and qs values at half levels)
3595 !          input are unadjusted t and q values,
3596 !          it returns adjusted values of t and q
3597
3598 !     parameter     description                                   units
3599 !     ---------     -----------                                   -----
3600 !     input parameters (integer):
3601
3602 !    *klon*         number of grid points per packet
3603 !    *klev*         number of levels
3604 !    *kk*           level
3605 !    *kcall*        defines calculation as
3606 !                      kcall=0  env. t and qs in*cuini*
3607 !                      kcall=1  condensation in updrafts  (e.g. cubase, cuasc)
3608 !                      kcall=2  evaporation in downdrafts (e.g. cudlfs,cuddraf)
3609 !     input parameters (real):
3610
3611 !    *psp*          pressure                                        pa
3612
3613 !     updated parameters (real):
3614
3615 !    *pt*           temperature                                     k
3616 !    *pq*           specific humidity                             kg/kg
3617 !          externals
3618 !          ---------
3619 !          for condensation calculations.
3620 !          the tables are initialised in *suphec*.
3621
3622 !----------------------------------------------------------------------
3623
3624       implicit none
3625
3626       integer  klev,klon
3627       real     pt(klon,klev),          pq(klon,klev),  &
3628      &         psp(klon)
3629       logical  ldflag(klon)
3630 ! local variables
3631       integer  jl,jk
3632       integer  isum,kcall,kk
3633       real     zqmax,zqsat,zcor,zqp,zcond,zcond1,zl,zi,zf
3634 !----------------------------------------------------------------------
3635 !     1.           define constants
3636 !                  ----------------
3637       zqmax=0.5
3638
3639 !     2.           calculate condensation and adjust t and q accordingly
3640 !                  -----------------------------------------------------
3641
3642       if ( kcall == 1 ) then
3643       do jl = 1,klon
3644         if ( ldflag(jl) ) then
3645           zqp = 1./psp(jl)
3646           zl = 1./(pt(jl,kk)-c4les)
3647           zi = 1./(pt(jl,kk)-c4ies)
3648           zqsat = c2es*(foealfa(pt(jl,kk))*exp(c3les*(pt(jl,kk)-tmelt)*zl) + &
3649                 (1.-foealfa(pt(jl,kk)))*exp(c3ies*(pt(jl,kk)-tmelt)*zi))
3650           zqsat = zqsat*zqp
3651           zqsat = min(0.5,zqsat)
3652           zcor = 1. - vtmpc1*zqsat
3653           zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
3654                (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
3655           zcond = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
3656           if ( zcond > 0. ) then
3657             pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
3658             pq(jl,kk) = pq(jl,kk) - zcond
3659             zl = 1./(pt(jl,kk)-c4les)
3660             zi = 1./(pt(jl,kk)-c4ies)
3661             zqsat = c2es*(foealfa(pt(jl,kk)) * &
3662               exp(c3les*(pt(jl,kk)-tmelt)*zl)+(1.-foealfa(pt(jl,kk))) * &
3663               exp(c3ies*(pt(jl,kk)-tmelt)*zi))
3664             zqsat = zqsat*zqp
3665             zqsat = min(0.5,zqsat)
3666             zcor = 1. - vtmpc1*zqsat
3667             zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + &
3668                  (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2
3669             zcond1 = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf)
3670             if ( abs(zcond) < 1.e-20 ) zcond1 = 0.
3671             pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3672             pq(jl,kk) = pq(jl,kk) - zcond1
3673           end if
3674         end if
3675       end do
3676       elseif ( kcall == 2 ) then
3677       do jl = 1,klon
3678         if ( ldflag(jl) ) then
3679           zqp = 1./psp(jl)
3680           zqsat = foeewm(pt(jl,kk))*zqp
3681           zqsat = min(0.5,zqsat)
3682           zcor = 1./(1.-vtmpc1*zqsat)
3683           zqsat = zqsat*zcor
3684           zcond = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3685           zcond = min(zcond,0.)
3686           pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond
3687           pq(jl,kk) = pq(jl,kk) - zcond
3688           zqsat = foeewm(pt(jl,kk))*zqp
3689           zqsat = min(0.5,zqsat)
3690           zcor = 1./(1.-vtmpc1*zqsat)
3691           zqsat = zqsat*zcor
3692           zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3693           if ( abs(zcond) < 1.e-20 ) zcond1 = min(zcond1,0.)
3694           pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3695           pq(jl,kk) = pq(jl,kk) - zcond1
3696         end if
3697       end do
3698       else if ( kcall == 0 ) then
3699       do jl = 1,klon
3700         zqp = 1./psp(jl)
3701         zqsat = foeewm(pt(jl,kk))*zqp
3702         zqsat = min(0.5,zqsat)
3703         zcor = 1./(1.-vtmpc1*zqsat)
3704         zqsat = zqsat*zcor
3705         zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3706         pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3707         pq(jl,kk) = pq(jl,kk) - zcond1
3708         zqsat = foeewm(pt(jl,kk))*zqp
3709         zqsat = min(0.5,zqsat)
3710         zcor = 1./(1.-vtmpc1*zqsat)
3711         zqsat = zqsat*zcor
3712         zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk)))
3713         pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1
3714         pq(jl,kk) = pq(jl,kk) - zcond1
3715       end do
3716       end if
3717
3718       return
3719       end subroutine cuadjtqn
3720 !---------------------------------------------------------
3721 !  level 4 souroutines
3722 !--------------------------------------------------------
3723       subroutine cubasmcn &
3724      &    (klon,     klev,     klevm1,  kk,     pten,&
3725      &     pqen,     pqsen,    puen,    pven,   pverv,&
3726      &     pgeo,     pgeoh,    ldcum,   ktype,  klab,  plrain,&
3727      &     pmfu,     pmfub,    kcbot,   ptu,&
3728      &     pqu,      plu,      puu,     pvu,    pmfus,&
3729      &     pmfuq,    pmful,    pdmfup  )
3730       implicit none
3731 !      m.tiedtke         e.c.m.w.f.     12/89
3732 !      c.zhang           iprc           05/2012
3733 !***purpose.
3734 !   --------
3735 !          this routine calculates cloud base values
3736 !          for midlevel convection
3737 !***interface
3738 !   ---------
3739 !          this routine is called from *cuasc*.
3740 !          input are environmental values t,q etc
3741 !          it returns cloudbase values for midlevel convection
3742 !***method.
3743 !   -------
3744 !          s. tiedtke (1989)
3745 !***externals
3746 !   ---------
3747 !          none
3748 ! ----------------------------------------------------------------
3749       real     pten(klon,klev),        pqen(klon,klev),&
3750      &         puen(klon,klev),        pven(klon,klev),&
3751      &         pqsen(klon,klev),       pverv(klon,klev),&
3752      &         pgeo(klon,klev),        pgeoh(klon,klev+1)
3753       real     ptu(klon,klev),         pqu(klon,klev),&
3754      &         puu(klon,klev),         pvu(klon,klev),&
3755      &         plu(klon,klev),         pmfu(klon,klev),&
3756      &         pmfub(klon),   &
3757      &         pmfus(klon,klev),       pmfuq(klon,klev),&
3758      &         pmful(klon,klev),       pdmfup(klon,klev),&
3759      &         plrain(klon,klev)
3760       integer  ktype(klon),            kcbot(klon),&
3761      &         klab(klon,klev)
3762       logical  ldcum(klon)
3763 ! local variabels
3764       integer  jl,kk,klev,klon,klevp1,klevm1
3765       real     zzzmb
3766 !--------------------------------------------------------
3767 !*    1.      calculate entrainment and detrainment rates
3768 ! -------------------------------------------------------
3769          do jl=1,klon
3770           if(.not.ldcum(jl) .and. klab(jl,kk+1).eq.0) then
3771             if(lmfmid .and. pqen(jl,kk) .gt. 0.80*pqsen(jl,kk).and. &
3772               pgeo(jl,kk)*zrg .gt. 5.0e2  .and.  &
3773      &        pgeo(jl,kk)*zrg .lt. 1.0e4 )  then
3774             ptu(jl,kk+1)=(cpd*pten(jl,kk)+pgeo(jl,kk)-pgeoh(jl,kk+1))&
3775      &                          *rcpd
3776             pqu(jl,kk+1)=pqen(jl,kk)
3777             plu(jl,kk+1)=0.
3778             zzzmb=max(cmfcmin,-pverv(jl,kk)*zrg)
3779             zzzmb=min(zzzmb,cmfcmax)
3780             pmfub(jl)=zzzmb
3781             pmfu(jl,kk+1)=pmfub(jl)
3782             pmfus(jl,kk+1)=pmfub(jl)*(cpd*ptu(jl,kk+1)+pgeoh(jl,kk+1))
3783             pmfuq(jl,kk+1)=pmfub(jl)*pqu(jl,kk+1)
3784             pmful(jl,kk+1)=0.
3785             pdmfup(jl,kk+1)=0.
3786             kcbot(jl)=kk
3787             klab(jl,kk+1)=1
3788             plrain(jl,kk+1)=0.0
3789             ktype(jl)=3
3790           end if
3791          end if
3792         end do
3793       return
3794       end subroutine cubasmcn
3795 !---------------------------------------------------------
3796 !  level 4 souroutines
3797 !---------------------------------------------------------
3798       subroutine cuentrn(klon,klev,kk,kcbot,ldcum,ldwork, &
3799                     pgeoh,pmfu,pdmfen,pdmfde)
3800        implicit none
3801        integer  klon,klev,kk
3802        integer  kcbot(klon)
3803        logical  ldcum(klon)
3804        logical  ldwork
3805        real  pgeoh(klon,klev+1)
3806        real  pmfu(klon,klev)
3807        real  pdmfen(klon)
3808        real  pdmfde(klon)
3809        logical  llo1
3810        integer  jl
3811        real  zdz , zmf
3812        real  zentr(klon)
3813     !
3814     !* 1. CALCULATE ENTRAINMENT AND DETRAINMENT RATES
3815     ! -------------------------------------------
3816     if ( ldwork ) then
3817       do jl = 1,klon
3818         pdmfen(jl) = 0.
3819         pdmfde(jl) = 0.
3820         zentr(jl) = 0.
3821       end do
3822       !
3823       !*  1.1 SPECIFY ENTRAINMENT RATES
3824       !   -------------------------
3825       do jl = 1, klon
3826         if ( ldcum(jl) ) then
3827           zdz = (pgeoh(jl,kk)-pgeoh(jl,kk+1))*zrg
3828           zmf = pmfu(jl,kk+1)*zdz
3829           llo1 = kk < kcbot(jl)
3830           if ( llo1 ) then
3831             pdmfen(jl) = zentr(jl)*zmf
3832             pdmfde(jl) = 0.75e-4*zmf
3833           end if
3834         end if
3835       end do
3836     end if
3837     end subroutine cuentrn
3838 !--------------------------------------------------------
3839 ! external functions
3840 !------------------------------------------------------
3841       real function foealfa(tt)
3842 !     foealfa is calculated to distinguish the three cases:
3843 !
3844 !                       foealfa=1            water phase
3845 !                       foealfa=0            ice phase
3846 !                       0 < foealfa < 1      mixed phase
3847 !
3848 !               input : tt = temperature
3849 !
3850         implicit none
3851         real tt
3852          foealfa = min(1.,((max(rtice,min(rtwat,tt))-rtice) &
3853      &  /(rtwat-rtice))**2)
3854
3855         return
3856       end function foealfa
3857
3858       real function foelhm(tt)
3859         implicit none
3860         real tt
3861         foelhm = foealfa(tt)*alv + (1.-foealfa(tt))*als
3862       return
3863       end function foelhm
3864
3865       real function foeewm(tt)
3866         implicit none
3867         real tt
3868         foeewm  = c2es * &
3869      &     (foealfa(tt)*exp(c3les*(tt-tmelt)/(tt-c4les))+ &
3870      &     (1.-foealfa(tt))*exp(c3ies*(tt-tmelt)/(tt-c4ies)))
3871       return
3872       end function foeewm
3873
3874       real function foedem(tt)
3875         implicit none
3876         real tt
3877         foedem  = foealfa(tt)*r5alvcp*(1./(tt-c4les)**2)+ &
3878      &              (1.-foealfa(tt))*r5alscp*(1./(tt-c4ies)**2)
3879       return
3880       end function foedem
3881
3882       real function foeldcpm(tt)
3883         implicit none
3884         real tt
3885         foeldcpm = foealfa(tt)*ralvdcp+ &
3886      &        (1.-foealfa(tt))*ralsdcp
3887       return
3888       end function  foeldcpm
3889
3890 end module module_cu_ntiedtke
3891