!> \file cu_ntiedtke.F90 !! This file contains the CCPP-compliant new Tiedtke scheme which parameterize !! Shallow, deep, and mid-level convections in the model !! Please refer to Tiedtke (1989), Bechtold et al. (2004,2008, 2014), !! Zhang et al. (2011), Zhang and Wang (2017, 2018) !! !########################################################### module cu_ntiedtke !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ use machine , only : kind_phys ! DH* TODO - replace with arguments to subroutine calls, ! this also requires redefining derived constants in the ! parameter section below use physcons, only:rd=>con_rd, rv=>con_rv, g=>con_g, & & cpd=>con_cp, alv=>con_hvap, alf=>con_hfus implicit none real(kind=kind_phys),private :: rcpd,vtmpc1,tmelt,als,t13, & c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg real(kind=kind_phys),private :: rovcp,r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp,rtwat,rtber,rtice real(kind=kind_phys),private :: entrdd,cmfcmax,cmfcmin,cmfdeps,zdnoprc,cprcon integer,private :: momtrans,p650 parameter( & t13 = 0.333333333,& rcpd=1.0/cpd, & tmelt=273.16, & zrg=1.0/g, & c1es=610.78, & c2es=c1es*rd/rv, & c3les=17.2693882, & c3ies=21.875, & c4les=35.86, & c4ies=7.66, & als = alv+alf, & c5les=c3les*(tmelt-c4les), & c5ies=c3ies*(tmelt-c4ies), & r5alvcp=c5les*alv*rcpd, & r5alscp=c5ies*als*rcpd, & ralvdcp=alv*rcpd, & ralsdcp=als*rcpd, & ralfdcp=alf*rcpd, & rtwat=tmelt, & rtber=tmelt-5., & rtice=tmelt-23., & vtmpc1=rv/rd-1.0, & rovcp = rd*rcpd ) ! ! entrdd: average entrainment & detrainment rate for downdrafts ! ------ ! parameter(entrdd = 2.0e-4) ! ! cmfcmax: maximum massflux value allowed for updrafts etc ! ------- ! parameter(cmfcmax = 1.0) ! ! cmfcmin: minimum massflux value (for safety) ! ------- ! parameter(cmfcmin = 1.e-10) ! ! cmfdeps: fractional massflux for downdrafts at lfs ! ------- ! parameter(cmfdeps = 0.30) ! zdnoprc: deep cloud is thicker than this height (Unit:Pa) ! parameter(zdnoprc = 2.0e4) ! ------- ! ! cprcon: coefficient from cloud water to rain water ! parameter(cprcon = 1.4e-3) ! ------- ! ! momtrans: momentum transport method ! ( 1 = IFS40r1 method; 2 = new method ) ! parameter(momtrans = 2 ) ! ------- ! logical :: isequil ! isequil: representing equilibrium and nonequilibrium convection ! ( .false. [default]; .true. [experimental]. Ref. Bechtold et al. 2014 JAS ) ! parameter(isequil = .false. ) ! !-------------------- ! switches for deep, mid, shallow convections, downdraft, and momemtum transport ! ------------------ logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.) !-------------------- !#################### end of variables definition########################## !----------------------------------------------------------------------- ! contains !> \brief Brief description of the subroutine !! !! \section arg_table_cu_ntiedtke_init Argument Table !! \htmlinclude cu_ntiedtke_init.html !! subroutine cu_ntiedtke_init(imfshalcnv, imfshalcnv_ntiedtke, imfdeepcnv, & imfdeepcnv_ntiedtke,mpirank, mpiroot, errmsg, errflg) implicit none integer, intent(in) :: imfshalcnv, imfshalcnv_ntiedtke integer, intent(in) :: imfdeepcnv, imfdeepcnv_ntiedtke integer, intent(in) :: mpirank integer, intent(in) :: mpiroot character(len=*), intent( out) :: errmsg integer, intent( out) :: errflg ! initialize ccpp error handling variables errmsg = '' errflg = 0 ! DH* temporary if (mpirank==mpiroot) then write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------' write(0,*) ' --- WARNING --- the CCPP New Tiedtke convection scheme is currently under development, use at your own risk --- WARNING ---' write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------' end if ! *DH temporary ! Consistency checks if (imfshalcnv/=imfshalcnv_ntiedtke) then write(errmsg,'(*(a))') 'Logic error: namelist choice of', & & ' shallow convection is different from new Tiedtke scheme' errflg = 1 return end if if (imfdeepcnv/=imfdeepcnv_ntiedtke) then write(errmsg,'(*(a))') 'Logic error: namelist choice of', & & ' deep convection is different from new Tiedtke scheme' errflg = 1 return end if end subroutine cu_ntiedtke_init ! Tiedtke cumulus scheme from WRF with small modifications ! This scheme includes both deep and shallow convections !=================== ! !> \section arg_table_cu_ntiedtke_run Argument Table !! \htmlinclude cu_ntiedtke_run.html !! !----------------------------------------------------------------------- ! level 1 subroutine 'tiecnvn' !----------------------------------------------------------------- subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & evap,hfx,zprecc,lmask,lq,km,dt,dx,kbot,ktop,kcnv, & ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg) !----------------------------------------------------------------- ! this is the interface between the model and the mass ! flux convection module !----------------------------------------------------------------- implicit none ! in&out variables integer, intent(in) :: lq, km, ktrac real(kind=kind_phys), intent(in ) :: dt integer, dimension( : ), intent(in) :: lmask real(kind=kind_phys), dimension( : ), intent(in ) :: evap, hfx, dx real(kind=kind_phys), dimension( :, : ), intent(inout) :: pu, pv, pt, pqv real(kind=kind_phys), dimension( :, :), intent(in ) :: tdi, qvdi, poz, prsl, pomg, pqvf, ptf real(kind=kind_phys), dimension( :, : ), intent(in ) :: pzz, prsi real(kind=kind_phys), dimension( :, :, : ), intent(inout ) :: clw integer, dimension( : ), intent(out) :: kbot, ktop, kcnv real(kind=kind_phys), dimension( : ), intent(out) :: zprecc real(kind=kind_phys), dimension (:, :), intent(out) :: ud_mf, dd_mf, dt_mf, cnvw, cnvc ! error messages character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! local variables real(kind=kind_phys) pum1(lq,km), pvm1(lq,km), ztt(lq,km), & & ptte(lq,km), pqte(lq,km), pvom(lq,km), pvol(lq,km), & & pverv(lq,km), pgeo(lq,km), pap(lq,km), paph(lq,km+1) real(kind=kind_phys) pqhfl(lq), zqq(lq,km), & & prsfc(lq), pssfc(lq), pcte(lq,km), & & phhfl(lq), pgeoh(lq,km+1) real(kind=kind_phys) ztp1(lq,km), zqp1(lq,km), ztu(lq,km), zqu(lq,km),& & zlu(lq,km), zlude(lq,km), zmfu(lq,km), zmfd(lq,km), zmfude_rate(lq,km),& & zqsat(lq,km), zrain(lq) real(kind=kind_phys),allocatable :: pcen(:,:,:),ptenc(:,:,:) integer icbot(lq), ictop(lq), ktype(lq), lndj(lq) logical locum(lq) ! real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt integer i,j,k,k1,n,km1,ktracer real(kind=kind_phys) ztpp1 real(kind=kind_phys) zew,zqs,zcor ! ! Initialize CCPP error handling variables errmsg = '' errflg = 0 km1 = km + 1 ztmst=dt ! ! masv flux diagnostics. ! do j=1,lq zrain(j)=0.0 locum(j)=.false. prsfc(j)=0.0 pssfc(j)=0.0 pqhfl(j)=evap(j) phhfl(j)=hfx(j) pgeoh(j,km1)=pzz(j,1) paph(j,km1)=prsi(j,1) if(lmask(j).eq.1) then lndj(j)=1 else lndj(j)=0 end if end do ! ! convert model variables for mflux scheme ! do k=1,km k1=km-k+1 do j=1,lq pcte(j,k1)=0.0 pvom(j,k1)=0.0 pvol(j,k1)=0.0 ztp1(j,k1)=pt(j,k) zqp1(j,k1)=pqv(j,k) pum1(j,k1)=pu(j,k) pvm1(j,k1)=pv(j,k) pverv(j,k1)=pomg(j,k) pgeo(j,k1)=poz(j,k) pgeoh(j,k1)=pzz(j,k+1) pap(j,k1)=prsl(j,k) paph(j,k1)=prsi(j,k+1) tt=ztp1(j,k1) zew = foeewm(tt) zqs = zew/pap(j,k1) zqs = min(0.5,zqs) zcor = 1./(1.-vtmpc1*zqs) zqsat(j,k1)=zqs*zcor pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst zqq(j,k1) =pqte(j,k1) ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst ztt(j,k1) =ptte(j,k1) ud_mf(j,k1)=0. dd_mf(j,k1)=0. dt_mf(j,k1)=0. cnvw(j,k1)=0. cnvc(j,k1)=0. end do end do if(ktrac > 2) then ktracer = ktrac - 2 allocate(pcen(lq,km,ktracer)) allocate(ptenc(lq,km,ktracer)) do n=1,ktracer do k=1,km k1=km-k+1 do j=1,lq pcen(j,k1,n) = clw(j,k,n+2) ptenc(j,k1,n)= 0. end do end do end do else ktracer = 2 allocate(pcen(lq,km,ktracer)) allocate(ptenc(lq,km,ktracer)) do n=1,ktracer do k=1,km do j=1,lq pcen(j,k,n) = 0. ptenc(j,k,n)= 0. end do end do end do end if ! print *, "pgeo=",pgeo(1,:) ! print *, "pgeoh=",pgeoh(1,:) ! print *, "pap=",pap(1,:) ! print *, "paph=",paph(1,:) ! print *, "ztp1=",ztp1(1,:) ! print *, "zqp1=",zqp1(1,:) ! print *, "pum1=",pum1(1,:) ! print *, "pvm1=",pvm1(1,:) ! print *, "pverv=",pverv(1,:) ! print *, "pqte=",pqte(1,:) ! print *, "ptte=",ptte(1,:) ! print *, "hfx=", pqhfl(1),phhfl(1),dx(1) ! !----------------------------------------------------------------------- !* 2. call 'cumastrn'(master-routine for cumulus parameterization) ! call cumastrn & & (lq, km, km1, km-1, ztp1, & & zqp1, pum1, pvm1, pverv, zqsat,& & pqhfl, ztmst, pap, paph, pgeo, & & ptte, pqte, pvom, pvol, prsfc,& & pssfc, locum, ktracer, pcen, ptenc,& & ktype, icbot, ictop, ztu, zqu, & & zlu, zlude, zmfu, zmfd, zrain,& & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx) ! ! to include the cloud water and cloud ice detrained from convection ! do k=1,km k1=km-k+1 do j=1,lq if(pcte(j,k1).gt.0.) then fliq=foealfa(ztp1(j,k1)) fice=1.0-fliq clw(j,k,2)=clw(j,k,2)+fliq*pcte(j,k1)*ztmst clw(j,k,1)=clw(j,k,1)+fice*pcte(j,k1)*ztmst endif end do end do ! do k=1,km k1 = km-k+1 do j=1,lq pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst ud_mf(j,k)= zmfu(j,k1)*ztmst dd_mf(j,k)= -zmfd(j,k1)*ztmst dt_mf(j,k)= zmfude_rate(j,k1)*ztmst cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1)) cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k)) cnvc(j,k) = min(cnvc(j,k), 0.6) cnvc(j,k) = max(cnvc(j,k), 0.0) end do end do do j=1,lq zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst*0.001) kbot(j) = km-icbot(j)+1 ktop(j) = km-ictop(j)+1 if(ktype(j).eq.1 .or. ktype(j).eq.3) then kcnv(j)=1 else kcnv(j)=0 end if end do if (lmfdudv) then do k=1,km k1=km-k+1 do j=1,lq pu(j,k)=pu(j,k)+pvom(j,k1)*ztmst pv(j,k)=pv(j,k)+pvol(j,k1)*ztmst end do end do endif ! ! Currently, vertical mixing of tracers are turned off ! if(ktrac > 2) then ! do n=1,ktrac-2 ! do k=1,km ! k1=km-k+1 ! do j=1,lq ! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst ! end do ! end do ! end do ! end if deallocate(pcen) deallocate(ptenc) ! return end subroutine cu_ntiedtke_run !############################################################# ! ! level 2 subroutines ! !############################################################# !*********************************************************** ! subroutine cumastrn !*********************************************************** subroutine cumastrn & & (klon, klev, klevp1, klevm1, pten, & & pqen, puen, pven, pverv, pqsen,& & pqhfl, ztmst, pap, paph, pgeo, & & ptte, pqte, pvom, pvol, prsfc,& & pssfc, ldcum, ktrac, pcen, ptenc,& & ktype, kcbot, kctop, ptu, pqu,& & plu, plude, pmfu, pmfd, prain,& & pcte, phhfl, lndj, zgeoh, pmfude_rate, dx) implicit none ! !***cumastrn* master routine for cumulus massflux-scheme ! m.tiedtke e.c.m.w.f. 1986/1987/1989 ! modifications ! y.wang i.p.r.c 2001 ! c.zhang 2012 !***purpose ! ------- ! this routine computes the physical tendencies of the ! prognostic variables t,q,u and v due to convective processes. ! processes considered are: convective fluxes, formation of ! precipitation, evaporation of falling rain below cloud base, ! saturated cumulus downdrafts. !***method ! ------ ! parameterization is done using a massflux-scheme. ! (1) define constants and parameters ! (2) specify values (t,q,qs...) at half levels and ! initialize updraft- and downdraft-values in 'cuinin' ! (3) calculate cloud base in 'cutypen', calculate cloud types in cutypen, ! and specify cloud base massflux ! (4) do cloud ascent in 'cuascn' in absence of downdrafts ! (5) do downdraft calculations: ! (a) determine values at lfs in 'cudlfsn' ! (b) determine moist descent in 'cuddrafn' ! (c) recalculate cloud base massflux considering the ! effect of cu-downdrafts ! (6) do final adjusments to convective fluxes in 'cuflxn', ! do evaporation in subcloud layer ! (7) calculate increments of t and q in 'cudtdqn' ! (8) calculate increments of u and v in 'cududvn' !***externals. ! ---------- ! cuinin: initializes values at vertical grid used in cu-parametr. ! cutypen: cloud bypes, 1: deep cumulus 2: shallow cumulus ! cuascn: cloud ascent for entraining plume ! cudlfsn: determines values at lfs for downdrafts ! cuddrafn:does moist descent for cumulus downdrafts ! cuflxn: final adjustments to convective fluxes (also in pbl) ! cudqdtn: updates tendencies for t and q ! cududvn: updates tendencies for u and v !***switches. ! -------- ! lmfmid=.t. midlevel convection is switched on ! lmfdd=.t. cumulus downdrafts switched on ! lmfdudv=.t. cumulus friction switched on !*** ! model parameters (defined in subroutine cuparam) ! ------------------------------------------------ ! entrdd entrainment rate for cumulus downdrafts ! cmfcmax maximum massflux value allowed for ! cmfcmin minimum massflux value (for safety) ! cmfdeps fractional massflux for downdrafts at lfs ! cprcon coefficient for conversion from cloud water to rain !***reference. ! ---------- ! paper on massflux scheme (tiedtke,1989) !----------------------------------------------------------------- integer klev,klon,ktrac,klevp1,klevm1 real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),& & puen(klon,klev), pven(klon,klev),& & ptte(klon,klev), pqte(klon,klev),& & pvom(klon,klev), pvol(klon,klev),& & pqsen(klon,klev), pgeo(klon,klev),& & pap(klon,klev), paph(klon,klevp1),& & pverv(klon,klev), pqhfl(klon),& & phhfl(klon) real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& & plu(klon,klev), plude(klon,klev),& & pmfu(klon,klev), pmfd(klon,klev),& & prain(klon),& & prsfc(klon), pssfc(klon) real(kind=kind_phys) ztenh(klon,klev), zqenh(klon,klev),& & zgeoh(klon,klevp1), zqsenh(klon,klev),& & ztd(klon,klev), zqd(klon,klev),& & zmfus(klon,klev), zmfds(klon,klev),& & zmfuq(klon,klev), zmfdq(klon,klev),& & zdmfup(klon,klev), zdmfdp(klon,klev),& & zmful(klon,klev), zrfl(klon),& & zuu(klon,klev), zvu(klon,klev),& & zud(klon,klev), zvd(klon,klev),& & zlglac(klon,klev) real(kind=kind_phys) pmflxr(klon,klevp1), pmflxs(klon,klevp1) real(kind=kind_phys) zhcbase(klon),& & zmfub(klon), zmfub1(klon),& & zdhpbl(klon) real(kind=kind_phys) zsfl(klon), zdpmel(klon,klev),& & pcte(klon,klev), zcape(klon),& & zcape1(klon), zcape2(klon),& & ztauc(klon), ztaubl(klon),& & zheat(klon) real(kind=kind_phys) pcen(klon,klev,ktrac), ptenc(klon,klev,ktrac) real(kind=kind_phys) wup(klon), zdqcv(klon) real(kind=kind_phys) wbase(klon), zmfuub(klon) real(kind=kind_phys) upbl(klon) real(kind=kind_phys) dx(klon) real(kind=kind_phys) pmfude_rate(klon,klev), pmfdde_rate(klon,klev) real(kind=kind_phys) zmfuus(klon,klev), zmfdus(klon,klev) real(kind=kind_phys) zmfudr(klon,klev), zmfddr(klon,klev) real(kind=kind_phys) zuv2(klon,klev),ztenu(klon,klev),ztenv(klon,klev) real(kind=kind_phys) zmfuvb(klon),zsum12(klon),zsum22(klon) integer ilab(klon,klev), idtop(klon),& & ictop0(klon), ilwmin(klon) integer kdpl(klon) integer kcbot(klon), kctop(klon),& & ktype(klon), lndj(klon) logical ldcum(klon), lldcum(klon) logical loddraf(klon), llddraf3(klon), llo1, llo2(klon) ! local varaiables real(kind=kind_phys) zcons,zcons2,zqumqe,zdqmin,zdh,zmfmax real(kind=kind_phys) zalfaw,zalv,zqalv,zc5ldcp,zc4les,zhsat,zgam,zzz,zhhat real(kind=kind_phys) zpbmpt,zro,zdz,zdp,zeps,zfac,wspeed integer jl,jk,ik integer ikb,ikt,icum,itopm2 real(kind=kind_phys) ztmst,ztau,zerate,zderate,zmfa real(kind=kind_phys) zmfs(klon),pmean(klev),zlon real(kind=kind_phys) zduten,zdvten,ztdis,pgf_u,pgf_v !------------------------------------------- ! 1. specify constants and parameters !------------------------------------------- zcons=1./(g*ztmst) zcons2=3./(g*ztmst) zlon = real(klon) do jk = klev , 1 , -1 pmean(jk) = sum(pap(:,jk))/zlon end do p650 = klev-2 do jk = klev , 3 , -1 if ( pmean(jk)/pmean(klev)*1.013250e5 > 650.e2 ) p650 = jk end do !-------------------------------------------------------------- !* 2. initialize values at vertical grid points in 'cuini' !-------------------------------------------------------------- call cuinin & & (klon, klev, klevp1, klevm1, pten, & & pqen, pqsen, puen, pven, pverv,& & pgeo, paph, zgeoh, ztenh, zqenh,& & zqsenh, ilwmin, ptu, pqu, ztd, & & zqd, zuu, zvu, zud, zvd, & & pmfu, pmfd, zmfus, zmfds, zmfuq,& & zmfdq, zdmfup, zdmfdp, zdpmel, plu, & & plude, ilab) !---------------------------------- !* 3.0 cloud base calculations !---------------------------------- !* (a) determine cloud base values in 'cutypen', ! and the cumulus type 1 or 2 ! ------------------------------------------- call cutypen & & ( klon, klev, klevp1, klevm1, pqen,& & ztenh, zqenh, zqsenh, zgeoh, paph,& & phhfl, pqhfl, pgeo, pqsen, pap,& & pten, lndj, ptu, pqu, ilab,& & ldcum, kcbot, ictop0, ktype, wbase, plu, kdpl) !* (b) assign the first guess mass flux at cloud base ! ------------------------------------------ do jl=1,klon zdhpbl(jl)=0.0 upbl(jl) = 0.0 idtop(jl)=0 end do do jk=2,klev do jl=1,klon if(jk.ge.kcbot(jl) .and. ldcum(jl)) then zdhpbl(jl)=zdhpbl(jl)+(alv*pqte(jl,jk)+cpd*ptte(jl,jk))& & *(paph(jl,jk+1)-paph(jl,jk)) if(lndj(jl) .eq. 0) then wspeed = sqrt(puen(jl,jk)**2 + pven(jl,jk)**2) upbl(jl) = upbl(jl) + wspeed*(paph(jl,jk+1)-paph(jl,jk)) end if end if end do end do do jl=1,klon if(ldcum(jl)) then ikb=kcbot(jl) zmfmax = (paph(jl,ikb)-paph(jl,ikb-1))*zcons2 if(ktype(jl) == 1) then zmfub(jl)= 0.1*zmfmax else if ( ktype(jl) == 2 ) then zqumqe = pqu(jl,ikb) + plu(jl,ikb) - zqenh(jl,ikb) zdqmin = max(0.01*zqenh(jl,ikb),1.e-10) zdh = cpd*(ptu(jl,ikb)-ztenh(jl,ikb)) + alv*zqumqe zdh = g*max(zdh,1.e5*zdqmin) if ( zdhpbl(jl) > 0. ) then zmfub(jl) = zdhpbl(jl)/zdh zmfub(jl) = min(zmfub(jl),zmfmax) else zmfub(jl) = 0.1*zmfmax ldcum(jl) = .false. end if end if else zmfub(jl) = 0. end if end do !------------------------------------------------------ !* 4.0 determine cloud ascent for entraining plume !------------------------------------------------------ !* (a) do ascent in 'cuasc'in absence of downdrafts !---------------------------------------------------------- call cuascn & & (klon, klev, klevp1, klevm1, ztenh,& & zqenh, puen, pven, pten, pqen,& & pqsen, pgeo, zgeoh, pap, paph,& & pqte, pverv, ilwmin, ldcum, zhcbase,& & ktype, ilab, ptu, pqu, plu,& & zuu, zvu, pmfu, zmfub,& & zmfus, zmfuq, zmful, plude, zdmfup,& & kcbot, kctop, ictop0, icum, ztmst,& & zqsenh, zlglac, lndj, wup, wbase, kdpl, pmfude_rate ) !* (b) check cloud depth and change entrainment rate accordingly ! calculate precipitation rate (for downdraft calculation) !------------------------------------------------------------------ do jl=1,klon if ( ldcum(jl) ) then ikb = kcbot(jl) itopm2 = kctop(jl) zpbmpt = paph(jl,ikb) - paph(jl,itopm2) if ( ktype(jl) == 1 .and. zpbmpt < zdnoprc ) ktype(jl) = 2 if ( ktype(jl) == 2 .and. zpbmpt >= zdnoprc ) ktype(jl) = 1 ictop0(jl) = kctop(jl) end if zrfl(jl)=zdmfup(jl,1) end do do jk=2,klev do jl=1,klon zrfl(jl)=zrfl(jl)+zdmfup(jl,jk) end do end do do jk = 1,klev do jl = 1,klon pmfd(jl,jk) = 0. zmfds(jl,jk) = 0. zmfdq(jl,jk) = 0. zdmfdp(jl,jk) = 0. zdpmel(jl,jk) = 0. end do end do !----------------------------------------- !* 6.0 cumulus downdraft calculations !----------------------------------------- if(lmfdd) then !* (a) determine lfs in 'cudlfsn' !-------------------------------------- call cudlfsn & & (klon, klev,& & kcbot, kctop, lndj, ldcum, & & ztenh, zqenh, puen, pven, & & pten, pqsen, pgeo, & & zgeoh, paph, ptu, pqu, plu, & & zuu, zvu, zmfub, zrfl, & & ztd, zqd, zud, zvd, & & pmfd, zmfds, zmfdq, zdmfdp, & & idtop, loddraf) !* (b) determine downdraft t,q and fluxes in 'cuddrafn' !------------------------------------------------------------ call cuddrafn & & ( klon, klev, loddraf, & & ztenh, zqenh, puen, pven, & & pgeo, zgeoh, paph, zrfl, & & ztd, zqd, zud, zvd, pmfu, & & pmfd, zmfds, zmfdq, zdmfdp, pmfdde_rate ) !----------------------------------------------------------- end if ! !----------------------------------------------------------------------- !* 6.0 closure and clean work ! ------ !-- 6.1 recalculate cloud base massflux from a cape closure ! for deep convection (ktype=1) ! do jl=1,klon if(ldcum(jl) .and. ktype(jl) .eq. 1) then ikb = kcbot(jl) ikt = kctop(jl) zheat(jl)=0.0 zcape(jl)=0.0 zcape1(jl)=0.0 zcape2(jl)=0.0 zmfub1(jl)=zmfub(jl) ztauc(jl) = (zgeoh(jl,ikt)-zgeoh(jl,ikb)) / & ((2.+ min(15.0,wup(jl)))*g) if(lndj(jl) .eq. 0) then upbl(jl) = 2.+ upbl(jl)/(paph(jl,klev+1)-paph(jl,ikb)) ztaubl(jl) = (zgeoh(jl,ikb)-zgeoh(jl,klev+1))/(g*upbl(jl)) ztaubl(jl) = min(300., ztaubl(jl)) else ztaubl(jl) = ztauc(jl) end if end if end do ! do jk = 1 , klev do jl = 1 , klon llo1 = ldcum(jl) .and. ktype(jl) .eq. 1 if ( llo1 .and. jk <= kcbot(jl) .and. jk > kctop(jl) ) then ikb = kcbot(jl) zdz = pgeo(jl,jk-1)-pgeo(jl,jk) zdp = pap(jl,jk)-pap(jl,jk-1) zheat(jl) = zheat(jl) + ((pten(jl,jk-1)-pten(jl,jk)+zdz*rcpd) / & ztenh(jl,jk)+vtmpc1*(pqen(jl,jk-1)-pqen(jl,jk))) * & (g*(pmfu(jl,jk)+pmfd(jl,jk))) zcape1(jl) = zcape1(jl) + ((ptu(jl,jk)-ztenh(jl,jk))/ztenh(jl,jk) + & vtmpc1*(pqu(jl,jk)-zqenh(jl,jk))-plu(jl,jk))*zdp end if if ( llo1 .and. jk >= kcbot(jl) ) then if((paph(jl,klev+1)-paph(jl,kdpl(jl)))<50.e2) then zdp = paph(jl,jk+1)-paph(jl,jk) zcape2(jl) = zcape2(jl) + ztaubl(jl)* & ((1.+vtmpc1*pqen(jl,jk))*ptte(jl,jk)+vtmpc1*pten(jl,jk)*pqte(jl,jk))*zdp end if end if end do end do do jl=1,klon if(ldcum(jl).and.ktype(jl).eq.1) then ikb = kcbot(jl) ikt = kctop(jl) ztau = ztauc(jl) * (1.+1.33e-5*dx(jl)) ztau = max(ztmst,ztau) ztau = max(720.,ztau) ztau = min(10800.,ztau) if(isequil) then zcape2(jl)= max(0.,zcape2(jl)) zcape(jl) = max(0.,min(zcape1(jl)-zcape2(jl),5000.)) else zcape(jl) = max(0.,min(zcape1(jl),5000.)) end if zheat(jl) = max(1.e-4,zheat(jl)) zmfub1(jl) = (zcape(jl)*zmfub(jl))/(zheat(jl)*ztau) zmfub1(jl) = max(zmfub1(jl),0.001) zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2 zmfub1(jl)=min(zmfub1(jl),zmfmax) end if end do ! !* 6.2 recalculate convective fluxes due to effect of ! downdrafts on boundary layer moist static energy budget (ktype=2) !-------------------------------------------------------- do jl=1,klon if(ldcum(jl) .and. ktype(jl) .eq. 2) then ikb=kcbot(jl) if(pmfd(jl,ikb).lt.0.0 .and. loddraf(jl)) then zeps=-pmfd(jl,ikb)/max(zmfub(jl),cmfcmin) else zeps=0. endif zqumqe=pqu(jl,ikb)+plu(jl,ikb)- & & zeps*zqd(jl,ikb)-(1.-zeps)*zqenh(jl,ikb) zdqmin=max(0.01*zqenh(jl,ikb),cmfcmin) zmfmax=(paph(jl,ikb)-paph(jl,ikb-1))*zcons2 ! using moist static engergy closure instead of moisture closure zdh=cpd*(ptu(jl,ikb)-zeps*ztd(jl,ikb)- & & (1.-zeps)*ztenh(jl,ikb))+alv*zqumqe zdh=g*max(zdh,1.e5*zdqmin) if(zdhpbl(jl).gt.0.)then zmfub1(jl)=zdhpbl(jl)/zdh else zmfub1(jl) = zmfub(jl) end if zmfub1(jl) = min(zmfub1(jl),zmfmax) end if !* 6.3 mid-level convection - nothing special !--------------------------------------------------------- if(ldcum(jl) .and. ktype(jl) .eq. 3 ) then zmfub1(jl) = zmfub(jl) end if end do !* 6.4 scaling the downdraft mass flux !--------------------------------------------------------- do jk=1,klev do jl=1,klon if( ldcum(jl) ) then zfac=zmfub1(jl)/max(zmfub(jl),cmfcmin) pmfd(jl,jk)=pmfd(jl,jk)*zfac zmfds(jl,jk)=zmfds(jl,jk)*zfac zmfdq(jl,jk)=zmfdq(jl,jk)*zfac zdmfdp(jl,jk)=zdmfdp(jl,jk)*zfac pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zfac end if end do end do !* 6.5 scaling the updraft mass flux ! -------------------------------------------------------- do jl = 1,klon if ( ldcum(jl) ) zmfs(jl) = zmfub1(jl)/max(cmfcmin,zmfub(jl)) end do do jk = 2 , klev do jl = 1,klon if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then ikb = kcbot(jl) if ( jk>ikb ) then zdz = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb))) pmfu(jl,jk) = pmfu(jl,ikb)*zdz end if zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2 if ( pmfu(jl,jk)*zmfs(jl) > zmfmax ) then zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk)) end if end if end do end do do jk = 2 , klev do jl = 1,klon if ( ldcum(jl) .and. jk <= kcbot(jl) .and. jk >= kctop(jl)-1 ) then pmfu(jl,jk) = pmfu(jl,jk)*zmfs(jl) zmfus(jl,jk) = zmfus(jl,jk)*zmfs(jl) zmfuq(jl,jk) = zmfuq(jl,jk)*zmfs(jl) zmful(jl,jk) = zmful(jl,jk)*zmfs(jl) zdmfup(jl,jk) = zdmfup(jl,jk)*zmfs(jl) plude(jl,jk) = plude(jl,jk)*zmfs(jl) pmfude_rate(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl) end if end do end do !* 6.6 if ktype = 2, kcbot=kctop is not allowed ! --------------------------------------------------- do jl = 1,klon if ( ktype(jl) == 2 .and. & kcbot(jl) == kctop(jl) .and. kcbot(jl) >= klev-1 ) then ldcum(jl) = .false. ktype(jl) = 0 end if end do if ( .not. lmfscv .or. .not. lmfpen ) then do jl = 1,klon llo2(jl) = .false. if ( (.not. lmfscv .and. ktype(jl) == 2) .or. & (.not. lmfpen .and. ktype(jl) == 1) ) then llo2(jl) = .true. ldcum(jl) = .false. end if end do end if !* 6.7 set downdraft mass fluxes to zero above cloud top !---------------------------------------------------- do jl = 1,klon if ( loddraf(jl) .and. idtop(jl) <= kctop(jl) ) then idtop(jl) = kctop(jl) + 1 end if end do do jk = 2 , klev do jl = 1,klon if ( loddraf(jl) ) then if ( jk < idtop(jl) ) then pmfd(jl,jk) = 0. zmfds(jl,jk) = 0. zmfdq(jl,jk) = 0. pmfdde_rate(jl,jk) = 0. zdmfdp(jl,jk) = 0. else if ( jk == idtop(jl) ) then pmfdde_rate(jl,jk) = 0. end if end if end do end do itopm2 = 2 !---------------------------------------------------------- !* 7.0 determine final convective fluxes in 'cuflx' !---------------------------------------------------------- call cuflxn & & ( klon, klev, ztmst & & , pten, pqen, pqsen, ztenh, zqenh & & , paph, pap, zgeoh, lndj, ldcum & & , kcbot, kctop, idtop, itopm2 & & , ktype, loddraf & & , pmfu, pmfd, zmfus, zmfds & & , zmfuq, zmfdq, zmful, plude & & , zdmfup, zdmfdp, zdpmel, zlglac & & , prain, pmfdde_rate, pmflxr, pmflxs ) ! some adjustments needed do jl=1,klon zmfs(jl) = 1. zmfuub(jl)=0. end do do jk = 2 , klev do jl = 1,klon if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then zmfmax = pmfu(jl,jk)*0.98 if ( pmfd(jl,jk)+zmfmax+1.e-15 < 0. ) then zmfs(jl) = min(zmfs(jl),-zmfmax/pmfd(jl,jk)) end if end if end do end do do jk = 2 , klev do jl = 1 , klon if ( zmfs(jl) < 1. .and. jk >= idtop(jl)-1 ) then pmfd(jl,jk) = pmfd(jl,jk)*zmfs(jl) zmfds(jl,jk) = zmfds(jl,jk)*zmfs(jl) zmfdq(jl,jk) = zmfdq(jl,jk)*zmfs(jl) pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl) zmfuub(jl) = zmfuub(jl) - (1.-zmfs(jl))*zdmfdp(jl,jk) pmflxr(jl,jk+1) = pmflxr(jl,jk+1) + zmfuub(jl) zdmfdp(jl,jk) = zdmfdp(jl,jk)*zmfs(jl) end if end do end do do jk = 2 , klev - 1 do jl = 1, klon if ( loddraf(jl) .and. jk >= idtop(jl)-1 ) then zerate = -pmfd(jl,jk) + pmfd(jl,jk-1) + pmfdde_rate(jl,jk) if ( zerate < 0. ) then pmfdde_rate(jl,jk) = pmfdde_rate(jl,jk) - zerate end if end if if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then zerate = pmfu(jl,jk) - pmfu(jl,jk+1) + pmfude_rate(jl,jk) if ( zerate < 0. ) then pmfude_rate(jl,jk) = pmfude_rate(jl,jk) - zerate end if zdmfup(jl,jk) = pmflxr(jl,jk+1) + pmflxs(jl,jk+1) - & pmflxr(jl,jk) - pmflxs(jl,jk) zdmfdp(jl,jk) = 0. end if end do end do ! avoid negative humidities at ddraught top do jl = 1,klon if ( loddraf(jl) ) then jk = idtop(jl) ik = min(jk+1,klev) if ( zmfdq(jl,jk) < 0.3*zmfdq(jl,ik) ) then zmfdq(jl,jk) = 0.3*zmfdq(jl,ik) end if end if end do ! avoid negative humidities near cloud top because gradient of precip flux ! and detrainment / liquid water flux are too large do jk = 2 , klev do jl = 1, klon if ( ldcum(jl) .and. jk >= kctop(jl)-1 .and. jk < kcbot(jl) ) then zdz = ztmst*g/(paph(jl,jk+1)-paph(jl,jk)) zmfa = zmfuq(jl,jk+1) + zmfdq(jl,jk+1) - & zmfuq(jl,jk) - zmfdq(jl,jk) + & zmful(jl,jk+1) - zmful(jl,jk) + zdmfup(jl,jk) zmfa = (zmfa-plude(jl,jk))*zdz if ( pqen(jl,jk)+zmfa < 0. ) then plude(jl,jk) = plude(jl,jk) + 2.*(pqen(jl,jk)+zmfa)/zdz end if if ( plude(jl,jk) < 0. ) plude(jl,jk) = 0. end if if ( .not. ldcum(jl) ) pmfude_rate(jl,jk) = 0. if ( abs(pmfd(jl,jk-1)) < 1.0e-20 ) pmfdde_rate(jl,jk) = 0. end do end do do jl=1,klon prsfc(jl) = pmflxr(jl,klev+1) pssfc(jl) = pmflxs(jl,klev+1) end do !---------------------------------------------------------------- !* 8.0 update tendencies for t and q in subroutine cudtdq !---------------------------------------------------------------- call cudtdqn(klon,klev,itopm2,kctop,idtop,ldcum,loddraf, & ztmst,paph,zgeoh,pgeo,pten,ztenh,pqen,zqenh,pqsen, & zlglac,plude,pmfu,pmfd,zmfus,zmfds,zmfuq,zmfdq,zmful, & zdmfup,zdmfdp,zdpmel,ptte,pqte,pcte) !---------------------------------------------------------------- !* 9.0 update tendencies for u and u in subroutine cududv !---------------------------------------------------------------- if(lmfdudv) then do jk = klev-1 , 2 , -1 ik = jk + 1 do jl = 1,klon if ( ldcum(jl) ) then if ( jk == kcbot(jl) .and. ktype(jl) < 3 ) then ikb = kdpl(jl) zuu(jl,jk) = puen(jl,ikb-1) zvu(jl,jk) = pven(jl,ikb-1) else if ( jk == kcbot(jl) .and. ktype(jl) == 3 ) then zuu(jl,jk) = puen(jl,jk-1) zvu(jl,jk) = pven(jl,jk-1) end if if ( jk < kcbot(jl) .and. jk >= kctop(jl) ) then if(momtrans .eq. 1)then zfac = 0. if ( ktype(jl) == 1 .or. ktype(jl) == 3 ) zfac = 2. if ( ktype(jl) == 1 .and. jk <= kctop(jl)+2 ) zfac = 3. zerate = pmfu(jl,jk) - pmfu(jl,ik) + & (1.+zfac)*pmfude_rate(jl,jk) zderate = (1.+zfac)*pmfude_rate(jl,jk) zmfa = 1./max(cmfcmin,pmfu(jl,jk)) zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + & zerate*puen(jl,jk)-zderate*zuu(jl,ik))*zmfa zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + & zerate*pven(jl,jk)-zderate*zvu(jl,ik))*zmfa else if(ktype(jl) == 1 .or. ktype(jl) == 3) then pgf_u = -0.7*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+& pmfu(jl,jk)*(puen(jl,jk)-puen(jl,jk-1))) pgf_v = -0.7*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+& pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1))) else pgf_u = 0. pgf_v = 0. end if zerate = pmfu(jl,jk) - pmfu(jl,ik) + pmfude_rate(jl,jk) zderate = pmfude_rate(jl,jk) zmfa = 1./max(cmfcmin,pmfu(jl,jk)) zuu(jl,jk) = (zuu(jl,ik)*pmfu(jl,ik) + & zerate*puen(jl,jk)-zderate*zuu(jl,ik)+pgf_u)*zmfa zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + & zerate*pven(jl,jk)-zderate*zvu(jl,ik)+pgf_v)*zmfa end if end if end if end do end do if(lmfdd) then do jk = 3 , klev ik = jk - 1 do jl = 1,klon if ( ldcum(jl) ) then if ( jk == idtop(jl) ) then zud(jl,jk) = 0.5*(zuu(jl,jk)+puen(jl,ik)) zvd(jl,jk) = 0.5*(zvu(jl,jk)+pven(jl,ik)) else if ( jk > idtop(jl) ) then zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pmfdde_rate(jl,jk) zmfa = 1./min(-cmfcmin,pmfd(jl,jk)) zud(jl,jk) = (zud(jl,ik)*pmfd(jl,ik) - & zerate*puen(jl,ik)+pmfdde_rate(jl,jk)*zud(jl,ik))*zmfa zvd(jl,jk) = (zvd(jl,ik)*pmfd(jl,ik) - & zerate*pven(jl,ik)+pmfdde_rate(jl,jk)*zvd(jl,ik))*zmfa end if end if end do end do end if ! -------------------------------------------------- ! rescale massfluxes for stability in Momentum !------------------------------------------------------------------------ zmfs(:) = 1. do jk = 2 , klev do jl = 1, klon if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons if ( pmfu(jl,jk) > zmfmax .and. jk >= kctop(jl) ) then zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk)) end if end if end do end do do jk = 1 , klev do jl = 1, klon zmfuus(jl,jk) = pmfu(jl,jk) zmfdus(jl,jk) = pmfd(jl,jk) if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl) zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl) end if end do end do !* 9.1 update u and v in subroutine cududvn !------------------------------------------------------------------- do jk = 1 , klev do jl = 1, klon ztenu(jl,jk) = pvom(jl,jk) ztenv(jl,jk) = pvol(jl,jk) end do end do call cududvn(klon,klev,itopm2,ktype,kcbot,kctop, & ldcum,ztmst,paph,puen,pven,zmfuus,zmfdus,zuu, & zud,zvu,zvd,pvom,pvol) ! calculate KE dissipation do jl = 1, klon zsum12(jl) = 0. zsum22(jl) = 0. end do do jk = 1 , klev do jl = 1, klon zuv2(jl,jk) = 0. if ( ldcum(jl) .and. jk >= kctop(jl)-1 ) then zdz = (paph(jl,jk+1)-paph(jl,jk)) zduten = pvom(jl,jk) - ztenu(jl,jk) zdvten = pvol(jl,jk) - ztenv(jl,jk) zuv2(jl,jk) = sqrt(zduten**2+zdvten**2) zsum22(jl) = zsum22(jl) + zuv2(jl,jk)*zdz zsum12(jl) = zsum12(jl) - & (puen(jl,jk)*zduten+pven(jl,jk)*zdvten)*zdz end if end do end do do jk = 1 , klev do jl = 1, klon if ( ldcum(jl) .and. jk>=kctop(jl)-1 ) then ztdis = rcpd*zsum12(jl)*zuv2(jl,jk)/max(1.e-15,zsum22(jl)) ptte(jl,jk) = ptte(jl,jk) + ztdis end if end do end do end if !---------------------------------------------------------------------- !* 10. IN CASE THAT EITHER DEEP OR SHALLOW IS SWITCHED OFF ! NEED TO SET SOME VARIABLES A POSTERIORI TO ZERO ! --------------------------------------------------- if ( .not. lmfscv .or. .not. lmfpen ) then do jk = 2 , klev do jl = 1, klon if ( llo2(jl) .and. jk >= kctop(jl)-1 ) then ptu(jl,jk) = pten(jl,jk) pqu(jl,jk) = pqen(jl,jk) plu(jl,jk) = 0. pmfude_rate(jl,jk) = 0. pmfdde_rate(jl,jk) = 0. end if end do end do do jl = 1, klon if ( llo2(jl) ) then kctop(jl) = klev - 1 kcbot(jl) = klev - 1 end if end do end if !---------------------------------------------------------------------- !* 11.0 CHEMICAL TRACER TRANSPORT ! ------------------------- if ( ktrac > 0 ) then ! transport switched off for mid-level convection do jl = 1, klon if ( ldcum(jl) .and. ktype(jl) /= 3 .and. & kcbot(jl)-kctop(jl) >= 1 ) then lldcum(jl) = .true. llddraf3(jl) = loddraf(jl) else lldcum(jl) = .false. llddraf3(jl) = .false. end if end do ! check and correct mass fluxes for CFL criterium zmfs(:) = 1. do jk = 2 , klev do jl = 1, klon if ( lldcum(jl) .and. jk >= kctop(jl) ) then zmfmax = (paph(jl,jk)-paph(jl,jk-1))*0.8*zcons if ( pmfu(jl,jk) > zmfmax ) then zmfs(jl) = min(zmfs(jl),zmfmax/pmfu(jl,jk)) end if end if end do end do do jk = 1, klev do jl = 1, klon if ( lldcum(jl) .and. jk >= kctop(jl)-1 ) then zmfuus(jl,jk) = pmfu(jl,jk)*zmfs(jl) zmfudr(jl,jk) = pmfude_rate(jl,jk)*zmfs(jl) else zmfuus(jl,jk) = 0. zmfudr(jl,jk) = 0. end if if ( llddraf3(jl) .and. jk >= idtop(jl)-1 ) then zmfdus(jl,jk) = pmfd(jl,jk)*zmfs(jl) zmfddr(jl,jk) = pmfdde_rate(jl,jk)*zmfs(jl) else zmfdus(jl,jk) = 0. zmfddr(jl,jk) = 0. end if end do end do call cuctracer(klon,klev,ktrac,kctop,idtop, & lldcum,llddraf3,ztmst,paph,zmfuus,zmfdus, & zmfudr,zmfddr,pcen,ptenc) end if return end subroutine cumastrn !********************************************** ! level 3 subroutine cuinin !********************************************** ! subroutine cuinin & & (klon, klev, klevp1, klevm1, pten,& & pqen, pqsen, puen, pven, pverv,& & pgeo, paph, pgeoh, ptenh, pqenh,& & pqsenh, klwmin, ptu, pqu, ptd,& & pqd, puu, pvu, pud, pvd,& & pmfu, pmfd, pmfus, pmfds, pmfuq,& & pmfdq, pdmfup, pdmfdp, pdpmel, plu,& & plude, klab) implicit none ! m.tiedtke e.c.m.w.f. 12/89 !***purpose ! ------- ! this routine interpolates large-scale fields of t,q etc. ! to half levels (i.e. grid for massflux scheme), ! and initializes values for updrafts and downdrafts !***interface ! --------- ! this routine is called from *cumastr*. !***method. ! -------- ! for extrapolation to half levels see tiedtke(1989) !***externals ! --------- ! *cuadjtq* to specify qs at half levels ! ---------------------------------------------------------------- integer klon,klev,klevp1,klevm1 real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),& & puen(klon,klev), pven(klon,klev),& & pqsen(klon,klev), pverv(klon,klev),& & pgeo(klon,klev), pgeoh(klon,klevp1),& & paph(klon,klevp1), ptenh(klon,klev),& & pqenh(klon,klev), pqsenh(klon,klev) real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& & ptd(klon,klev), pqd(klon,klev),& & puu(klon,klev), pud(klon,klev),& & pvu(klon,klev), pvd(klon,klev),& & pmfu(klon,klev), pmfd(klon,klev),& & pmfus(klon,klev), pmfds(klon,klev),& & pmfuq(klon,klev), pmfdq(klon,klev),& & pdmfup(klon,klev), pdmfdp(klon,klev),& & plu(klon,klev), plude(klon,klev) real(kind=kind_phys) zwmax(klon), zph(klon), & & pdpmel(klon,klev) integer klab(klon,klev), klwmin(klon) logical loflag(klon) ! local variables integer jl,jk integer icall,ik real(kind=kind_phys) zzs !------------------------------------------------------------ !* 1. specify large scale parameters at half levels !* adjust temperature fields if staticly unstable !* find level of maximum vertical velocity ! ----------------------------------------------------------- do jk=2,klev do jl=1,klon ptenh(jl,jk)=(max(cpd*pten(jl,jk-1)+pgeo(jl,jk-1), & & cpd*pten(jl,jk)+pgeo(jl,jk))-pgeoh(jl,jk))*rcpd pqenh(jl,jk) = pqen(jl,jk-1) pqsenh(jl,jk)= pqsen(jl,jk-1) zph(jl)=paph(jl,jk) loflag(jl)=.true. end do if ( jk >= klev-1 .or. jk < 2 ) cycle ik=jk icall=0 call cuadjtqn(klon,klev,ik,zph,ptenh,pqsenh,loflag,icall) do jl=1,klon pqenh(jl,jk)=min(pqen(jl,jk-1),pqsen(jl,jk-1)) & & +(pqsenh(jl,jk)-pqsen(jl,jk-1)) pqenh(jl,jk)=max(pqenh(jl,jk),0.) end do end do do jl=1,klon ptenh(jl,klev)=(cpd*pten(jl,klev)+pgeo(jl,klev)- & & pgeoh(jl,klev))*rcpd pqenh(jl,klev)=pqen(jl,klev) ptenh(jl,1)=pten(jl,1) pqenh(jl,1)=pqen(jl,1) klwmin(jl)=klev zwmax(jl)=0. end do do jk=klevm1,2,-1 do jl=1,klon zzs=max(cpd*ptenh(jl,jk)+pgeoh(jl,jk), & & cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1)) ptenh(jl,jk)=(zzs-pgeoh(jl,jk))*rcpd end do end do do jk=klev,3,-1 do jl=1,klon if(pverv(jl,jk).lt.zwmax(jl)) then zwmax(jl)=pverv(jl,jk) klwmin(jl)=jk end if end do end do !----------------------------------------------------------- !* 2.0 initialize values for updrafts and downdrafts !----------------------------------------------------------- do jk=1,klev ik=jk-1 if(jk.eq.1) ik=1 do jl=1,klon ptu(jl,jk)=ptenh(jl,jk) ptd(jl,jk)=ptenh(jl,jk) pqu(jl,jk)=pqenh(jl,jk) pqd(jl,jk)=pqenh(jl,jk) plu(jl,jk)=0. puu(jl,jk)=puen(jl,ik) pud(jl,jk)=puen(jl,ik) pvu(jl,jk)=pven(jl,ik) pvd(jl,jk)=pven(jl,ik) klab(jl,jk)=0 end do end do return end subroutine cuinin !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cutypen & & ( klon, klev, klevp1, klevm1, pqen,& & ptenh, pqenh, pqsenh, pgeoh, paph,& & hfx, qfx, pgeo, pqsen, pap,& & pten, lndj, cutu, cuqu, culab,& & ldcum, cubot, cutop, ktype, wbase, culu, kdpl ) ! zhang & wang iprc 2011-2013 !***purpose. ! -------- ! to produce first guess updraught for cu-parameterizations ! calculates condensation level, and sets updraught base variables and ! first guess cloud type !***interface ! --------- ! this routine is called from *cumastr*. ! input are environm. values of t,q,p,phi at half levels. ! it returns cloud types as follows; ! ktype=1 for deep cumulus ! ktype=2 for shallow cumulus !***method. ! -------- ! based on a simplified updraught equation ! partial(hup)/partial(z)=eta(h - hup) ! eta is the entrainment rate for test parcel ! h stands for dry static energy or the total water specific humidity ! references: christian jakob, 2003: a new subcloud model for ! mass-flux convection schemes ! influence on triggering, updraft properties, and model ! climate, mon.wea.rev. ! 131, 2765-2778 ! and ! ifs documentation - cy36r1,cy38r1 !***input variables: ! ptenh [ztenh] - environment temperature on half levels. (cuini) ! pqenh [zqenh] - env. specific humidity on half levels. (cuini) ! pgeoh [zgeoh] - geopotential on half levels, (mssflx) ! paph - pressure of half levels. (mssflx) ! rho - density of the lowest model level ! qfx - net upward moisture flux at the surface (kg/m^2/s) ! hfx - net upward heat flux at the surface (w/m^2) !***variables output by cutype: ! ktype - convection type - 1: penetrative (cumastr) ! 2: stratocumulus (cumastr) ! 3: mid-level (cuasc) ! information for updraft parcel (ptu,pqu,plu,kcbot,klab,kdpl...) ! ---------------------------------------------------------------- !------------------------------------------------------------------- implicit none !------------------------------------------------------------------- integer klon, klev, klevp1, klevm1 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev),& & pqsen(klon,klev), pqsenh(klon,klev),& & pgeoh(klon,klevp1), paph(klon,klevp1),& & pap(klon,klev), pqen(klon,klev) real(kind=kind_phys) pten(klon,klev) real(kind=kind_phys) ptu(klon,klev),pqu(klon,klev),plu(klon,klev) real(kind=kind_phys) pgeo(klon,klev) integer klab(klon,klev) integer kctop(klon),kcbot(klon) real(kind=kind_phys) qfx(klon),hfx(klon) real(kind=kind_phys) zph(klon) integer lndj(klon) logical loflag(klon), deepflag(klon), resetflag(klon) ! output variables real(kind=kind_phys) cutu(klon,klev), cuqu(klon,klev), culu(klon,klev) integer culab(klon,klev) real(kind=kind_phys) wbase(klon) integer ktype(klon),cubot(klon),cutop(klon),kdpl(klon) logical ldcum(klon) ! local variables real(kind=kind_phys) zqold(klon) real(kind=kind_phys) rho, part1, part2, root, conw, deltt, deltq real(kind=kind_phys) eta(klon),dz(klon),coef(klon) real(kind=kind_phys) dhen(klon,klev), dh(klon,klev) real(kind=kind_phys) plude(klon,klev) real(kind=kind_phys) kup(klon,klev) real(kind=kind_phys) vptu(klon,klev),vten(klon,klev) real(kind=kind_phys) zbuo(klon,klev),abuoy(klon,klev) real(kind=kind_phys) zz,zdken,zdq real(kind=kind_phys) fscale,crirh1,pp real(kind=kind_phys) atop1,atop2,abot real(kind=kind_phys) tmix,zmix,qmix,pmix real(kind=kind_phys) zlglac,dp integer nk,is,ikb,ikt real(kind=kind_phys) zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp real(kind=kind_phys) zpdifftop, zpdiffbot integer zcbase(klon), itoppacel(klon) integer jl,jk,ik,icall,levels logical needreset, lldcum(klon) !-------------------------------------------------------------- do jl=1,klon kcbot(jl)=klev kctop(jl)=klev kdpl(jl) =klev ktype(jl)=0 wbase(jl)=0. ldcum(jl)=.false. end do !----------------------------------------------------------- ! let's do test,and check the shallow convection first ! the first level is klev ! define deltat and deltaq !----------------------------------------------------------- do jk=1,klev do jl=1,klon plu(jl,jk)=culu(jl,jk) ! parcel liquid water ptu(jl,jk)=cutu(jl,jk) ! parcel temperature pqu(jl,jk)=cuqu(jl,jk) ! parcel specific humidity klab(jl,jk)=culab(jl,jk) dh(jl,jk)=0.0 ! parcel dry static energy dhen(jl,jk)=0.0 ! environment dry static energy kup(jl,jk)=0.0 ! updraught kinetic energy for parcel vptu(jl,jk)=0.0 ! parcel virtual temperature considering water-loading vten(jl,jk)=0.0 ! environment virtual temperature zbuo(jl,jk)=0.0 ! parcel buoyancy abuoy(jl,jk)=0.0 end do end do do jl=1,klon zqold(jl) = 0. lldcum(jl) = .false. loflag(jl) = .true. end do ! check the levels from lowest level to second top level do jk=klevm1,2,-1 ! define the variables at the first level if(jk .eq. klevm1) then do jl=1,klon rho=pap(jl,klev)/ & & (rd*(pten(jl,klev)*(1.+vtmpc1*pqen(jl,klev)))) hfx(jl) = hfx(jl)*rho*cpd qfx(jl) = qfx(jl)*rho part1 = 1.5*0.4*pgeo(jl,klev)/ & & (rho*pten(jl,klev)) part2 = -hfx(jl)*rcpd-vtmpc1*pten(jl,klev)*qfx(jl) root = 0.001-part1*part2 if(part2 .lt. 0.) then conw = 1.2*(root)**t13 deltt = max(1.5*hfx(jl)/(rho*cpd*conw),0.) deltq = max(1.5*qfx(jl)/(rho*conw),0.) kup(jl,klev) = 0.5*(conw**2) pqu(jl,klev)= pqenh(jl,klev) + deltq dhen(jl,klev)= pgeoh(jl,klev) + ptenh(jl,klev)*cpd dh(jl,klev) = dhen(jl,klev) + deltt*cpd ptu(jl,klev) = (dh(jl,klev)-pgeoh(jl,klev))*rcpd vptu(jl,klev)=ptu(jl,klev)*(1.+vtmpc1*pqu(jl,klev)) vten(jl,klev)=ptenh(jl,klev)*(1.+vtmpc1*pqenh(jl,klev)) zbuo(jl,klev)=(vptu(jl,klev)-vten(jl,klev))/vten(jl,klev) klab(jl,klev) = 1 else loflag(jl) = .false. end if end do end if is=0 do jl=1,klon if(loflag(jl))then is=is+1 endif enddo if(is.eq.0) exit ! the next levels, we use the variables at the first level as initial values do jl=1,klon if(loflag(jl)) then eta(jl) = 0.55/(pgeo(jl,jk)*zrg)+1.e-4 dz(jl) = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg coef(jl)= 0.5*eta(jl)*dz(jl) dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk) dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))& & +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl)) pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))& & +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl)) ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd zqold(jl) = pqu(jl,jk) zph(jl)=paph(jl,jk) end if end do ! check if the parcel is saturated ik=jk icall=1 call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall) do jl=1,klon if( loflag(jl) ) then zdq = max((zqold(jl) - pqu(jl,jk)),0.) plu(jl,jk) = plu(jl,jk+1) + zdq zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - & (1.-foealfa(ptu(jl,jk+1)))) plu(jl,jk) = min(plu(jl,jk),5.e-3) dh(jl,jk) = pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac) ! compute the updraft speed vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+& ralfdcp*zlglac vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)) zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk) abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g atop1 = 1.0 - 2.*coef(jl) atop2 = 2.0*dz(jl)*abuoy(jl,jk) abot = 1.0 + 2.*coef(jl) kup(jl,jk) = (atop1*kup(jl,jk+1) + atop2) / abot ! let's find the exact cloud base if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then ik = jk + 1 zqsu = foeewm(ptu(jl,ik))/paph(jl,ik) zqsu = min(0.5,zqsu) zcor = 1./(1.-vtmpc1*zqsu) zqsu = zqsu*zcor zdq = min(0.,pqu(jl,ik)-zqsu) zalfaw = foealfa(ptu(jl,ik)) zfacw = c5les/((ptu(jl,ik)-c4les)**2) zfaci = c5ies/((ptu(jl,ik)-c4ies)**2) zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci zesdp = foeewm(ptu(jl,ik))/paph(jl,ik) zcor = 1./(1.-vtmpc1*zesdp) zdqsdt = zfac*zcor*zqsu zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik)) zdp = zdq/(zdqsdt*zdtdp) zcbase(jl) = paph(jl,ik) + zdp ! chose nearest half level as cloud base (jk or jk+1) zpdifftop = zcbase(jl) - paph(jl,jk) zpdiffbot = paph(jl,jk+1) - zcbase(jl) if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then ikb = min(klev-1,jk+1) klab(jl,ikb) = 2 klab(jl,jk) = 2 kcbot(jl) = ikb plu(jl,jk+1) = 1.0e-8 else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then klab(jl,jk) = 2 kcbot(jl) = jk end if end if if(kup(jl,jk) .lt. 0.)then loflag(jl) = .false. if(plu(jl,jk+1) .gt. 0.) then kctop(jl) = jk lldcum(jl) = .true. else lldcum(jl) = .false. end if else if(plu(jl,jk) .gt. 0.)then klab(jl,jk)=2 else klab(jl,jk)=1 end if end if end if end do end do ! end all the levels do jl=1,klon ikb = kcbot(jl) ikt = kctop(jl) if(paph(jl,ikb) - paph(jl,ikt) > zdnoprc) lldcum(jl) = .false. if(lldcum(jl)) then ktype(jl) = 2 ldcum(jl) = .true. wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.)) cubot(jl) = ikb cutop(jl) = ikt kdpl(jl) = klev else cutop(jl) = -1 cubot(jl) = -1 kdpl(jl) = klev - 1 ldcum(jl) = .false. wbase(jl) = 0. end if end do do jk=klev,1,-1 do jl=1,klon ikt = kctop(jl) if(jk .ge. ikt)then culab(jl,jk) = klab(jl,jk) cutu(jl,jk) = ptu(jl,jk) cuqu(jl,jk) = pqu(jl,jk) culu(jl,jk) = plu(jl,jk) end if end do end do !----------------------------------------------------------- ! next, let's check the deep convection ! the first level is klevm1-1 ! define deltat and deltaq !---------------------------------------------------------- ! we check the parcel starting level by level ! assume the mix-layer is 60hPa deltt = 0.2 deltq = 1.0e-4 do jl=1,klon deepflag(jl) = .false. end do do jk=klev,1,-1 do jl=1,klon if((paph(jl,klev+1)-paph(jl,jk)) .lt. 350.e2) itoppacel(jl) = jk end do end do do levels=klevm1-1,klevm1-20,-1 ! loop starts do jk=1,klev do jl=1,klon plu(jl,jk)=0.0 ! parcel liquid water ptu(jl,jk)=0.0 ! parcel temperature pqu(jl,jk)=0.0 ! parcel specific humidity dh(jl,jk)=0.0 ! parcel dry static energy dhen(jl,jk)=0.0 ! environment dry static energy kup(jl,jk)=0.0 ! updraught kinetic energy for parcel vptu(jl,jk)=0.0 ! parcel virtual temperature consideringwater-loading vten(jl,jk)=0.0 ! environment virtual temperature abuoy(jl,jk)=0.0 zbuo(jl,jk)=0.0 klab(jl,jk)=0 end do end do do jl=1,klon kcbot(jl) = levels kctop(jl) = levels zqold(jl) = 0. lldcum(jl) = .false. resetflag(jl)= .false. loflag(jl) = (.not. deepflag(jl)) .and. (levels.ge.itoppacel(jl)) end do ! start the inner loop to search the deep convection points do jk=levels,2,-1 is=0 do jl=1,klon if(loflag(jl))then is=is+1 endif enddo if(is.eq.0) exit ! define the variables at the departure level if(jk .eq. levels) then do jl=1,klon if(loflag(jl)) then if((paph(jl,klev+1)-paph(jl,jk)) < 60.e2) then tmix=0. qmix=0. zmix=0. pmix=0. do nk=jk+2,jk,-1 if(pmix < 50.e2) then dp = paph(jl,nk) - paph(jl,nk-1) tmix=tmix+dp*ptenh(jl,nk) qmix=qmix+dp*pqenh(jl,nk) zmix=zmix+dp*pgeoh(jl,nk) pmix=pmix+dp end if end do tmix=tmix/pmix qmix=qmix/pmix zmix=zmix/pmix else tmix=ptenh(jl,jk+1) qmix=pqenh(jl,jk+1) zmix=pgeoh(jl,jk+1) end if pqu(jl,jk+1) = qmix + deltq dhen(jl,jk+1)= zmix + tmix*cpd dh(jl,jk+1) = dhen(jl,jk+1) + deltt*cpd ptu(jl,jk+1) = (dh(jl,jk+1)-pgeoh(jl,jk+1))*rcpd kup(jl,jk+1) = 0.5 klab(jl,jk+1)= 1 vptu(jl,jk+1)=ptu(jl,jk+1)*(1.+vtmpc1*pqu(jl,jk+1)) vten(jl,jk+1)=ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1)) zbuo(jl,jk+1)=(vptu(jl,jk+1)-vten(jl,jk+1))/vten(jl,jk+1) end if end do end if ! the next levels, we use the variables at the first level as initial values do jl=1,klon if(loflag(jl)) then ! define the fscale fscale = min(1.,(pqsen(jl,jk)/pqsen(jl,levels))**3) eta(jl) = 1.75e-3*fscale dz(jl) = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg coef(jl)= 0.5*eta(jl)*dz(jl) dhen(jl,jk) = pgeoh(jl,jk) + cpd*ptenh(jl,jk) dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk))& & +(1.-coef(jl))*dh(jl,jk+1))/(1.+coef(jl)) pqu(jl,jk) =(coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk))& & +(1.-coef(jl))*pqu(jl,jk+1))/(1.+coef(jl)) ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd zqold(jl) = pqu(jl,jk) zph(jl)=paph(jl,jk) end if end do ! check if the parcel is saturated ik=jk icall=1 call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall) do jl=1,klon if( loflag(jl) ) then zdq = max((zqold(jl) - pqu(jl,jk)),0.) plu(jl,jk) = plu(jl,jk+1) + zdq zlglac=zdq*((1.-foealfa(ptu(jl,jk))) - & (1.-foealfa(ptu(jl,jk+1)))) plu(jl,jk) = 0.5*plu(jl,jk) dh(jl,jk) = pgeoh(jl,jk) + cpd*(ptu(jl,jk)+ralfdcp*zlglac) ! compute the updraft speed vptu(jl,jk) = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk))+& ralfdcp*zlglac vten(jl,jk) = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)) zbuo(jl,jk) = (vptu(jl,jk) - vten(jl,jk))/vten(jl,jk) abuoy(jl,jk)=(zbuo(jl,jk)+zbuo(jl,jk+1))*0.5*g atop1 = 1.0 - 2.*coef(jl) atop2 = 2.0*dz(jl)*abuoy(jl,jk) abot = 1.0 + 2.*coef(jl) kup(jl,jk) = (atop1*kup(jl,jk+1) + atop2) / abot ! let's find the exact cloud base if ( plu(jl,jk) > 0. .and. klab(jl,jk+1) == 1 ) then ik = jk + 1 zqsu = foeewm(ptu(jl,ik))/paph(jl,ik) zqsu = min(0.5,zqsu) zcor = 1./(1.-vtmpc1*zqsu) zqsu = zqsu*zcor zdq = min(0.,pqu(jl,ik)-zqsu) zalfaw = foealfa(ptu(jl,ik)) zfacw = c5les/((ptu(jl,ik)-c4les)**2) zfaci = c5ies/((ptu(jl,ik)-c4ies)**2) zfac = zalfaw*zfacw + (1.-zalfaw)*zfaci zesdp = foeewm(ptu(jl,ik))/paph(jl,ik) zcor = 1./(1.-vtmpc1*zesdp) zdqsdt = zfac*zcor*zqsu zdtdp = rd*ptu(jl,ik)/(cpd*paph(jl,ik)) zdp = zdq/(zdqsdt*zdtdp) zcbase(jl) = paph(jl,ik) + zdp ! chose nearest half level as cloud base (jk or jk+1) zpdifftop = zcbase(jl) - paph(jl,jk) zpdiffbot = paph(jl,jk+1) - zcbase(jl) if ( zpdifftop > zpdiffbot .and. kup(jl,jk+1) > 0. ) then ikb = min(klev-1,jk+1) klab(jl,ikb) = 2 klab(jl,jk) = 2 kcbot(jl) = ikb plu(jl,jk+1) = 1.0e-8 else if ( zpdifftop <= zpdiffbot .and.kup(jl,jk) > 0. ) then klab(jl,jk) = 2 kcbot(jl) = jk end if end if if(kup(jl,jk) .lt. 0.)then loflag(jl) = .false. if(plu(jl,jk+1) .gt. 0.) then kctop(jl) = jk lldcum(jl) = .true. else lldcum(jl) = .false. end if else if(plu(jl,jk) .gt. 0.)then klab(jl,jk)=2 else klab(jl,jk)=1 end if end if end if end do end do ! end all the levels needreset = .false. do jl=1,klon ikb = kcbot(jl) ikt = kctop(jl) if(paph(jl,ikb) - paph(jl,ikt) < zdnoprc) lldcum(jl) = .false. if(lldcum(jl)) then ktype(jl) = 1 ldcum(jl) = .true. deepflag(jl) = .true. wbase(jl) = sqrt(max(2.*kup(jl,ikb),0.)) cubot(jl) = ikb cutop(jl) = ikt kdpl(jl) = levels+1 needreset = .true. resetflag(jl)= .true. end if end do if(needreset) then do jk=klev,1,-1 do jl=1,klon if(resetflag(jl)) then ikt = kctop(jl) ikb = kdpl(jl) if(jk .le. ikb .and. jk .ge. ikt )then culab(jl,jk) = klab(jl,jk) cutu(jl,jk) = ptu(jl,jk) cuqu(jl,jk) = pqu(jl,jk) culu(jl,jk) = plu(jl,jk) else culab(jl,jk) = 1 cutu(jl,jk) = ptenh(jl,jk) cuqu(jl,jk) = pqenh(jl,jk) culu(jl,jk) = 0. end if if ( jk .lt. ikt ) culab(jl,jk) = 0 end if end do end do end if end do ! end all cycles return end subroutine cutypen !----------------------------------------------------------------- ! level 3 subroutines 'cuascn' !----------------------------------------------------------------- subroutine cuascn & & (klon, klev, klevp1, klevm1, ptenh,& & pqenh, puen, pven, pten, pqen,& & pqsen, pgeo, pgeoh, pap, paph,& & pqte, pverv, klwmin, ldcum, phcbase,& & ktype, klab, ptu, pqu, plu,& & puu, pvu, pmfu, pmfub, & & pmfus, pmfuq, pmful, plude, pdmfup,& & kcbot, kctop, kctop0, kcum, ztmst,& & pqsenh, plglac, lndj, wup, wbase, kdpl, pmfude_rate) implicit none ! this routine does the calculations for cloud ascents ! for cumulus parameterization ! m.tiedtke e.c.m.w.f. 7/86 modif. 12/89 ! y.wang iprc 11/01 modif. ! c.zhang iprc 05/12 modif. !***purpose. ! -------- ! to produce cloud ascents for cu-parametrization ! (vertical profiles of t,q,l,u and v and corresponding ! fluxes as well as precipitation rates) !***interface ! --------- ! this routine is called from *cumastr*. !***method. ! -------- ! lift surface air dry-adiabatically to cloud base ! and then calculate moist ascent for ! entraining/detraining plume. ! entrainment and detrainment rates differ for ! shallow and deep cumulus convection. ! in case there is no penetrative or shallow convection ! check for possibility of mid level convection ! (cloud base values calculated in *cubasmc*) !***externals ! --------- ! *cuadjtqn* adjust t and q due to condensation in ascent ! *cuentrn* calculate entrainment/detrainment rates ! *cubasmcn* calculate cloud base values for midlevel convection !***reference ! --------- ! (tiedtke,1989) !***input variables: ! ptenh [ztenh] - environ temperature on half levels. (cuini) ! pqenh [zqenh] - env. specific humidity on half levels. (cuini) ! puen - environment wind u-component. (mssflx) ! pven - environment wind v-component. (mssflx) ! pten - environment temperature. (mssflx) ! pqen - environment specific humidity. (mssflx) ! pqsen - environment saturation specific humidity. (mssflx) ! pgeo - geopotential. (mssflx) ! pgeoh [zgeoh] - geopotential on half levels, (mssflx) ! pap - pressure in pa. (mssflx) ! paph - pressure of half levels. (mssflx) ! pqte - moisture convergence (delta q/delta t). (mssflx) ! pverv - large scale vertical velocity (omega). (mssflx) ! klwmin [ilwmin] - level of minimum omega. (cuini) ! klab [ilab] - level label - 1: sub-cloud layer. ! 2: condensation level (cloud base) ! pmfub [zmfub] - updraft mass flux at cloud base. (cumastr) !***variables modified by cuasc: ! ldcum - logical denoting profiles. (cubase) ! ktype - convection type - 1: penetrative (cumastr) ! 2: stratocumulus (cumastr) ! 3: mid-level (cuasc) ! ptu - cloud temperature. ! pqu - cloud specific humidity. ! plu - cloud liquid water (moisture condensed out) ! puu [zuu] - cloud momentum u-component. ! pvu [zvu] - cloud momentum v-component. ! pmfu - updraft mass flux. ! pmfus [zmfus] - updraft flux of dry static energy. (cubasmc) ! pmfuq [zmfuq] - updraft flux of specific humidity. ! pmful [zmful] - updraft flux of cloud liquid water. ! plude - liquid water returned to environment by detrainment. ! pdmfup [zmfup] - ! kcbot - cloud base level. (cubase) ! kctop - cloud top level ! kctop0 [ictop0] - estimate of cloud top. (cumastr) ! kcum [icum] - flag to control the call integer klev,klon,klevp1,klevm1 real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & & puen(klon,klev), pven(klon,klev),& & pten(klon,klev), pqen(klon,klev),& & pgeo(klon,klev), pgeoh(klon,klevp1),& & pap(klon,klev), paph(klon,klevp1),& & pqsen(klon,klev), pqte(klon,klev),& & pverv(klon,klev), pqsenh(klon,klev) real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& & puu(klon,klev), pvu(klon,klev),& & pmfu(klon,klev), zph(klon),& & pmfub(klon), & & pmfus(klon,klev), pmfuq(klon,klev),& & plu(klon,klev), plude(klon,klev),& & pmful(klon,klev), pdmfup(klon,klev) real(kind=kind_phys) zdmfen(klon), zdmfde(klon),& & zmfuu(klon), zmfuv(klon),& & zpbase(klon), zqold(klon) real(kind=kind_phys) phcbase(klon), zluold(klon) real(kind=kind_phys) zprecip(klon), zlrain(klon,klev) real(kind=kind_phys) zbuo(klon,klev), kup(klon,klev) real(kind=kind_phys) wup(klon) real(kind=kind_phys) wbase(klon), zodetr(klon,klev) real(kind=kind_phys) plglac(klon,klev) real(kind=kind_phys) eta(klon),dz(klon) integer klwmin(klon), ktype(klon),& & klab(klon,klev), kcbot(klon),& & kctop(klon), kctop0(klon) integer lndj(klon) logical ldcum(klon), loflag(klon) logical llo2,llo3, llo1(klon) integer kdpl(klon) real(kind=kind_phys) zoentr(klon), zdpmean(klon) real(kind=kind_phys) pdmfen(klon,klev), pmfude_rate(klon,klev) ! local variables integer jl,jk integer ikb,icum,itopm2,ik,icall,is,kcum,jlm,jll integer jlx(klon) real(kind=kind_phys) ztmst,zcons2,zfacbuo,zprcdgw,z_cwdrag,z_cldmax,z_cwifrac,z_cprc2 real(kind=kind_phys) zmftest,zmfmax,zqeen,zseen,zscde,zqude real(kind=kind_phys) zmfusk,zmfuqk,zmfulk real(kind=kind_phys) zbc,zbe,zkedke,zmfun,zwu,zprcon,zdt,zcbf,zzco real(kind=kind_phys) zlcrit,zdfi,zc,zd,zint,zlnew,zvw,zvi,zalfaw,zrold real(kind=kind_phys) zrnew,zz,zdmfeu,zdmfdu,dp real(kind=kind_phys) zfac,zbuoc,zdkbuo,zdken,zvv,zarg,zchange,zxe,zxs,zdshrd real(kind=kind_phys) atop1,atop2,abot !-------------------------------- !* 1. specify parameters !-------------------------------- zcons2=3./(g*ztmst) zfacbuo = 0.5/(1.+0.5) zprcdgw = cprcon*zrg z_cldmax = 5.e-3 z_cwifrac = 0.5 z_cprc2 = 0.5 z_cwdrag = (3.0/8.0)*0.506/0.2 !--------------------------------- ! 2. set default values !--------------------------------- llo3 = .false. do jl=1,klon zluold(jl)=0. wup(jl)=0. zdpmean(jl)=0. zoentr(jl)=0. if(.not.ldcum(jl)) then ktype(jl)=0 kcbot(jl) = -1 pmfub(jl) = 0. pqu(jl,klev) = 0. end if end do ! initialize variout quantities do jk=1,klev do jl=1,klon if(jk.ne.kcbot(jl)) plu(jl,jk)=0. pmfu(jl,jk)=0. pmfus(jl,jk)=0. pmfuq(jl,jk)=0. pmful(jl,jk)=0. plude(jl,jk)=0. plglac(jl,jk)=0. pdmfup(jl,jk)=0. zlrain(jl,jk)=0. zbuo(jl,jk)=0. kup(jl,jk)=0. pdmfen(jl,jk) = 0. pmfude_rate(jl,jk) = 0. if(.not.ldcum(jl).or.ktype(jl).eq.3) klab(jl,jk)=0 if(.not.ldcum(jl).and.paph(jl,jk).lt.4.e4) kctop0(jl)=jk end do end do do jl = 1,klon if ( ktype(jl) == 3 ) ldcum(jl) = .false. end do !------------------------------------------------ ! 3.0 initialize values at cloud base level !------------------------------------------------ do jl=1,klon kctop(jl)=kcbot(jl) if(ldcum(jl)) then ikb = kcbot(jl) kup(jl,ikb) = 0.5*wbase(jl)**2 pmfu(jl,ikb) = pmfub(jl) pmfus(jl,ikb) = pmfub(jl)*(cpd*ptu(jl,ikb)+pgeoh(jl,ikb)) pmfuq(jl,ikb) = pmfub(jl)*pqu(jl,ikb) pmful(jl,ikb) = pmfub(jl)*plu(jl,ikb) end if end do ! !----------------------------------------------------------------- ! 4. do ascent: subcloud layer (klab=1) ,clouds (klab=2) ! by doing first dry-adiabatic ascent and then ! by adjusting t,q and l accordingly in *cuadjtqn*, ! then check for buoyancy and set flags accordingly !----------------------------------------------------------------- ! do jk=klevm1,3,-1 ! specify cloud base values for midlevel convection ! in *cubasmc* in case there is not already convection ! --------------------------------------------------------------------- ik=jk call cubasmcn& & (klon, klev, klevm1, ik, pten,& & pqen, pqsen, puen, pven, pverv,& & pgeo, pgeoh, ldcum, ktype, klab, zlrain,& & pmfu, pmfub, kcbot, ptu,& & pqu, plu, puu, pvu, pmfus,& & pmfuq, pmful, pdmfup) is = 0 jlm = 0 do jl = 1,klon loflag(jl) = .false. zprecip(jl) = 0. llo1(jl) = .false. is = is + klab(jl,jk+1) if ( klab(jl,jk+1) == 0 ) klab(jl,jk) = 0 if ( (ldcum(jl) .and. klab(jl,jk+1) == 2) .or. & (ktype(jl) == 3 .and. klab(jl,jk+1) == 1) ) then loflag(jl) = .true. jlm = jlm + 1 jlx(jlm) = jl end if zph(jl) = paph(jl,jk) if ( ktype(jl) == 3 .and. jk == kcbot(jl) ) then zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2 if ( pmfub(jl) > zmfmax ) then zfac = zmfmax/pmfub(jl) pmfu(jl,jk+1) = pmfu(jl,jk+1)*zfac pmfus(jl,jk+1) = pmfus(jl,jk+1)*zfac pmfuq(jl,jk+1) = pmfuq(jl,jk+1)*zfac pmfub(jl) = zmfmax end if pmfub(jl)=min(pmfub(jl),zmfmax) end if end do if(is.gt.0) llo3 = .true. ! !* specify entrainment rates in *cuentr* ! ------------------------------------- ik=jk call cuentrn(klon,klev,ik,kcbot,ldcum,llo3, & pgeoh,pmfu,zdmfen,zdmfde) ! ! do adiabatic ascent for entraining/detraining plume if(llo3) then ! ------------------------------------------------------- ! do jl = 1,klon zqold(jl) = 0. end do do jll = 1 , jlm jl = jlx(jll) zdmfde(jl) = min(zdmfde(jl),0.75*pmfu(jl,jk+1)) if ( jk == kcbot(jl) ) then zoentr(jl) = -1.75e-3*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - & 1.)*(pgeoh(jl,jk)-pgeoh(jl,jk+1))*zrg zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk+1) end if if ( jk < kcbot(jl) ) then zmfmax = (paph(jl,jk)-paph(jl,jk-1))*zcons2 zxs = max(pmfu(jl,jk+1)-zmfmax,0.) wup(jl) = wup(jl) + kup(jl,jk+1)*(pap(jl,jk+1)-pap(jl,jk)) zdpmean(jl) = zdpmean(jl) + pap(jl,jk+1) - pap(jl,jk) zdmfen(jl) = zoentr(jl) if ( ktype(jl) >= 2 ) then zdmfen(jl) = 2.0*zdmfen(jl) zdmfde(jl) = zdmfen(jl) end if zdmfde(jl) = zdmfde(jl) * & (1.6-min(1.,pqen(jl,jk)/pqsen(jl,jk))) zmftest = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl) zchange = max(zmftest-zmfmax,0.) zxe = max(zchange-zxs,0.) zdmfen(jl) = zdmfen(jl) - zxe zchange = zchange - zxe zdmfde(jl) = zdmfde(jl) + zchange end if pdmfen(jl,jk) = zdmfen(jl) - zdmfde(jl) pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl) zqeen = pqenh(jl,jk+1)*zdmfen(jl) zseen = (cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1))*zdmfen(jl) zscde = (cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1))*zdmfde(jl) zqude = pqu(jl,jk+1)*zdmfde(jl) plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl) zmfusk = pmfus(jl,jk+1) + zseen - zscde zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude zmfulk = pmful(jl,jk+1) - plude(jl,jk) plu(jl,jk) = zmfulk*(1./max(cmfcmin,pmfu(jl,jk))) pqu(jl,jk) = zmfuqk*(1./max(cmfcmin,pmfu(jl,jk))) ptu(jl,jk) = (zmfusk * & (1./max(cmfcmin,pmfu(jl,jk)))-pgeoh(jl,jk))*rcpd ptu(jl,jk) = max(100.,ptu(jl,jk)) ptu(jl,jk) = min(400.,ptu(jl,jk)) zqold(jl) = pqu(jl,jk) zlrain(jl,jk) = zlrain(jl,jk+1)*(pmfu(jl,jk+1)-zdmfde(jl)) * & (1./max(cmfcmin,pmfu(jl,jk))) zluold(jl) = plu(jl,jk) end do ! reset to environmental values if below departure level do jl = 1,klon if ( jk > kdpl(jl) ) then ptu(jl,jk) = ptenh(jl,jk) pqu(jl,jk) = pqenh(jl,jk) plu(jl,jk) = 0. zluold(jl) = plu(jl,jk) end if end do !* do corrections for moist ascent !* by adjusting t,q and l in *cuadjtq* !------------------------------------------------ ik=jk icall=1 ! if ( jlm > 0 ) then call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall) end if ! compute the upfraft speed in cloud layer do jll = 1 , jlm jl = jlx(jll) if ( pqu(jl,jk) /= zqold(jl) ) then plglac(jl,jk) = plu(jl,jk) * & ((1.-foealfa(ptu(jl,jk)))- & (1.-foealfa(ptu(jl,jk+1)))) ptu(jl,jk) = ptu(jl,jk) + ralfdcp*plglac(jl,jk) end if end do do jll = 1 , jlm jl = jlx(jll) if ( pqu(jl,jk) /= zqold(jl) ) then klab(jl,jk) = 2 plu(jl,jk) = plu(jl,jk) + zqold(jl) - pqu(jl,jk) zbc = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk+1) - & zlrain(jl,jk+1)) zbe = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)) zbuo(jl,jk) = zbc - zbe ! set flags for the case of midlevel convection if ( ktype(jl) == 3 .and. klab(jl,jk+1) == 1 ) then if ( zbuo(jl,jk) > -0.5 ) then ldcum(jl) = .true. kctop(jl) = jk kup(jl,jk) = 0.5 else klab(jl,jk) = 0 pmfu(jl,jk) = 0. plude(jl,jk) = 0. plu(jl,jk) = 0. end if end if if ( klab(jl,jk+1) == 2 ) then if ( zbuo(jl,jk) < 0. ) then ptenh(jl,jk) = 0.5*(pten(jl,jk)+pten(jl,jk-1)) pqenh(jl,jk) = 0.5*(pqen(jl,jk)+pqen(jl,jk-1)) zbuo(jl,jk) = zbc - ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)) end if zbuoc = (zbuo(jl,jk) / & (ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)))+zbuo(jl,jk+1) / & (ptenh(jl,jk+1)*(1.+vtmpc1*pqenh(jl,jk+1))))*0.5 zdkbuo = (pgeoh(jl,jk)-pgeoh(jl,jk+1))*zfacbuo*zbuoc ! mixing and "pressure" gradient term in upper troposphere if ( zdmfen(jl) > 0. ) then zdken = min(1.,(1.+z_cwdrag)*zdmfen(jl) / & max(cmfcmin,pmfu(jl,jk+1))) else zdken = min(1.,(1.+z_cwdrag)*zdmfde(jl) / & max(cmfcmin,pmfu(jl,jk+1))) end if kup(jl,jk) = (kup(jl,jk+1)*(1.-zdken)+zdkbuo) / & (1.+zdken) if ( zbuo(jl,jk) < 0. ) then zkedke = kup(jl,jk)/max(1.e-10,kup(jl,jk+1)) zkedke = max(0.,min(1.,zkedke)) zmfun = sqrt(zkedke)*pmfu(jl,jk+1) !* (1.6-min(1.,pqen(jl,jk) / & ! pqsen(jl,jk))) zdmfde(jl) = max(zdmfde(jl),pmfu(jl,jk+1)-zmfun) plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl) pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl) end if if ( zbuo(jl,jk) > -0.2 ) then ikb = kcbot(jl) zoentr(jl) = 1.75e-3*(0.3-(min(1.,pqen(jl,jk-1) / & pqsen(jl,jk-1))-1.))*(pgeoh(jl,jk-1)-pgeoh(jl,jk)) * & zrg*min(1.,pqsen(jl,jk)/pqsen(jl,ikb))**3 zoentr(jl) = min(0.4,zoentr(jl))*pmfu(jl,jk) else zoentr(jl) = 0. end if ! erase values if below departure level if ( jk > kdpl(jl) ) then pmfu(jl,jk) = pmfu(jl,jk+1) kup(jl,jk) = 0.5 end if if ( kup(jl,jk) > 0. .and. pmfu(jl,jk) > 0. ) then kctop(jl) = jk llo1(jl) = .true. else klab(jl,jk) = 0 pmfu(jl,jk) = 0. kup(jl,jk) = 0. zdmfde(jl) = pmfu(jl,jk+1) plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl) end if ! save detrainment rates for updraught if ( pmfu(jl,jk+1) > 0. ) pmfude_rate(jl,jk) = zdmfde(jl) end if else if ( ktype(jl) == 2 .and. pqu(jl,jk) == zqold(jl) ) then klab(jl,jk) = 0 pmfu(jl,jk) = 0. kup(jl,jk) = 0. zdmfde(jl) = pmfu(jl,jk+1) plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl) pmfude_rate(jl,jk) = zdmfde(jl) end if end do do jl = 1,klon if ( llo1(jl) ) then ! conversions only proceeds if plu is greater than a threshold liquid water ! content of 0.3 g/kg over water and 0.5 g/kg over land to prevent precipitation ! generation from small water contents. if ( lndj(jl).eq.1 ) then zdshrd = 5.e-4 else zdshrd = 3.e-4 end if ikb=kcbot(jl) if ( plu(jl,jk) > zdshrd )then ! if ((paph(jl,ikb)-paph(jl,jk))>zdnoprc) then zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk+1)))) zprcon = zprcdgw/(0.75*zwu) ! PARAMETERS FOR BERGERON-FINDEISEN PROCESS (T < -5C) zdt = min(rtber-rtice,max(rtber-ptu(jl,jk),0.)) zcbf = 1. + z_cprc2*sqrt(zdt) zzco = zprcon*zcbf zlcrit = zdshrd/zcbf zdfi = pgeoh(jl,jk) - pgeoh(jl,jk+1) zc = (plu(jl,jk)-zluold(jl)) zarg = (plu(jl,jk)/zlcrit)**2 if ( zarg < 25.0 ) then zd = zzco*(1.-exp(-zarg))*zdfi else zd = zzco*zdfi end if zint = exp(-zd) zlnew = zluold(jl)*zint + zc/zd*(1.-zint) zlnew = max(0.,min(plu(jl,jk),zlnew)) zlnew = min(z_cldmax,zlnew) zprecip(jl) = max(0.,zluold(jl)+zc-zlnew) pdmfup(jl,jk) = zprecip(jl)*pmfu(jl,jk) zlrain(jl,jk) = zlrain(jl,jk) + zprecip(jl) plu(jl,jk) = zlnew end if end if end do do jl = 1, klon if ( llo1(jl) ) then if ( zlrain(jl,jk) > 0. ) then zvw = 21.18*zlrain(jl,jk)**0.2 zvi = z_cwifrac*zvw zalfaw = foealfa(ptu(jl,jk)) zvv = zalfaw*zvw + (1.-zalfaw)*zvi zrold = zlrain(jl,jk) - zprecip(jl) zc = zprecip(jl) zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk)))) zd = zvv/zwu zint = exp(-zd) zrnew = zrold*zint + zc/zd*(1.-zint) zrnew = max(0.,min(zlrain(jl,jk),zrnew)) zlrain(jl,jk) = zrnew end if end if end do do jll = 1 , jlm jl = jlx(jll) pmful(jl,jk) = plu(jl,jk)*pmfu(jl,jk) pmfus(jl,jk) = (cpd*ptu(jl,jk)+pgeoh(jl,jk))*pmfu(jl,jk) pmfuq(jl,jk) = pqu(jl,jk)*pmfu(jl,jk) end do end if end do !---------------------------------------------------------------------- ! 5. final calculations ! ------------------ do jl = 1,klon if ( kctop(jl) == -1 ) ldcum(jl) = .false. kcbot(jl) = max(kcbot(jl),kctop(jl)) if ( ldcum(jl) ) then wup(jl) = max(1.e-2,wup(jl)/max(1.,zdpmean(jl))) wup(jl) = sqrt(2.*wup(jl)) end if end do return end subroutine cuascn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cudlfsn & & (klon, klev, & & kcbot, kctop, lndj, ldcum, & & ptenh, pqenh, puen, pven, & & pten, pqsen, pgeo, & & pgeoh, paph, ptu, pqu, plu,& & puu, pvu, pmfub, prfl, & & ptd, pqd, pud, pvd, & & pmfd, pmfds, pmfdq, pdmfdp, & & kdtop, lddraf) ! this routine calculates level of free sinking for ! cumulus downdrafts and specifies t,q,u and v values ! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89 ! purpose. ! -------- ! to produce lfs-values for cumulus downdrafts ! for massflux cumulus parameterization ! interface ! --------- ! this routine is called from *cumastr*. ! input are environmental values of t,q,u,v,p,phi ! and updraft values t,q,u and v and also ! cloud base massflux and cu-precipitation rate. ! it returns t,q,u and v values and massflux at lfs. ! method. ! check for negative buoyancy of air of equal parts of ! moist environmental air and cloud air. ! parameter description units ! --------- ----------- ----- ! input parameters (integer): ! *klon* number of grid points per packet ! *klev* number of levels ! *kcbot* cloud base level ! *kctop* cloud top level ! input parameters (logical): ! *lndj* land sea mask (1 for land) ! *ldcum* flag: .true. for convective points ! input parameters (real(kind=kind_phys)): ! *ptenh* env. temperature (t+1) on half levels k ! *pqenh* env. spec. humidity (t+1) on half levels kg/kg ! *puen* provisional environment u-velocity (t+1) m/s ! *pven* provisional environment v-velocity (t+1) m/s ! *pten* provisional environment temperature (t+1) k ! *pqsen* environment spec. saturation humidity (t+1) kg/kg ! *pgeo* geopotential m2/s2 ! *pgeoh* geopotential on half levels m2/s2 ! *paph* provisional pressure on half levels pa ! *ptu* temperature in updrafts k ! *pqu* spec. humidity in updrafts kg/kg ! *plu* liquid water content in updrafts kg/kg ! *puu* u-velocity in updrafts m/s ! *pvu* v-velocity in updrafts m/s ! *pmfub* massflux in updrafts at cloud base kg/(m2*s) ! updated parameters (real(kind=kind_phys)): ! *prfl* precipitation rate kg/(m2*s) ! output parameters (real(kind=kind_phys)): ! *ptd* temperature in downdrafts k ! *pqd* spec. humidity in downdrafts kg/kg ! *pud* u-velocity in downdrafts m/s ! *pvd* v-velocity in downdrafts m/s ! *pmfd* massflux in downdrafts kg/(m2*s) ! *pmfds* flux of dry static energy in downdrafts j/(m2*s) ! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s) ! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s) ! output parameters (integer): ! *kdtop* top level of downdrafts ! output parameters (logical): ! *lddraf* .true. if downdrafts exist ! externals ! --------- ! *cuadjtq* for calculating wet bulb t and q at lfs !---------------------------------------------------------------------- implicit none integer klev,klon real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & & puen(klon,klev), pven(klon,klev), & & pten(klon,klev), pqsen(klon,klev), & & pgeo(klon,klev), & & pgeoh(klon,klev+1), paph(klon,klev+1),& & ptu(klon,klev), pqu(klon,klev), & & puu(klon,klev), pvu(klon,klev), & & plu(klon,klev), & & pmfub(klon), prfl(klon) real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), & & pud(klon,klev), pvd(klon,klev), & & pmfd(klon,klev), pmfds(klon,klev), & & pmfdq(klon,klev), pdmfdp(klon,klev) integer kcbot(klon), kctop(klon), & & kdtop(klon), ikhsmin(klon) logical ldcum(klon), & & lddraf(klon) integer lndj(klon) real(kind=kind_phys) ztenwb(klon,klev), zqenwb(klon,klev), & & zcond(klon), zph(klon), & & zhsmin(klon) logical llo2(klon) ! local variables integer jl,jk integer is,ik,icall,ike real(kind=kind_phys) zhsk,zttest,zqtest,zbuo,zmftop !---------------------------------------------------------------------- ! 1. set default values for downdrafts ! --------------------------------- do jl=1,klon lddraf(jl)=.false. kdtop(jl)=klev+1 ikhsmin(jl)=klev+1 zhsmin(jl)=1.e8 enddo !---------------------------------------------------------------------- ! 2. determine level of free sinking: ! downdrafts shall start at model level of minimum ! of saturation moist static energy or below ! respectively ! for every point and proceed as follows: ! (1) determine level of minimum of hs ! (2) determine wet bulb environmental t and q ! (3) do mixing with cumulus cloud air ! (4) check for negative buoyancy ! (5) if buoyancy>0 repeat (2) to (4) for next ! level below ! the assumption is that air of downdrafts is mixture ! of 50% cloud air + 50% environmental air at wet bulb ! temperature (i.e. which became saturated due to ! evaporation of rain and cloud water) ! ---------------------------------------------------- do jk=3,klev-2 do jl=1,klon zhsk=cpd*pten(jl,jk)+pgeo(jl,jk) + & & foelhm(pten(jl,jk))*pqsen(jl,jk) if(zhsk .lt. zhsmin(jl)) then zhsmin(jl) = zhsk ikhsmin(jl)= jk end if end do end do ike=klev-3 do jk=3,ike ! 2.1 calculate wet-bulb temperature and moisture ! for environmental air in *cuadjtq* ! ------------------------------------------- is=0 do jl=1,klon ztenwb(jl,jk)=ptenh(jl,jk) zqenwb(jl,jk)=pqenh(jl,jk) zph(jl)=paph(jl,jk) llo2(jl)=ldcum(jl).and.prfl(jl).gt.0..and..not.lddraf(jl).and. & & (jk.lt.kcbot(jl).and.jk.gt.kctop(jl)).and. jk.ge.ikhsmin(jl) if(llo2(jl))then is=is+1 endif enddo if(is.eq.0) cycle ik=jk icall=2 call cuadjtqn & & ( klon, klev, ik, zph, ztenwb, zqenwb, llo2, icall) ! 2.2 do mixing of cumulus and environmental air ! and check for negative buoyancy. ! then set values for downdraft at lfs. ! ---------------------------------------- do jl=1,klon if(llo2(jl)) then zttest=0.5*(ptu(jl,jk)+ztenwb(jl,jk)) zqtest=0.5*(pqu(jl,jk)+zqenwb(jl,jk)) zbuo=zttest*(1.+vtmpc1 *zqtest)- & & ptenh(jl,jk)*(1.+vtmpc1 *pqenh(jl,jk)) zcond(jl)=pqenh(jl,jk)-zqenwb(jl,jk) zmftop=-cmfdeps*pmfub(jl) if(zbuo.lt.0..and.prfl(jl).gt.10.*zmftop*zcond(jl)) then kdtop(jl)=jk lddraf(jl)=.true. ptd(jl,jk)=zttest pqd(jl,jk)=zqtest pmfd(jl,jk)=zmftop pmfds(jl,jk)=pmfd(jl,jk)*(cpd*ptd(jl,jk)+pgeoh(jl,jk)) pmfdq(jl,jk)=pmfd(jl,jk)*pqd(jl,jk) pdmfdp(jl,jk-1)=-0.5*pmfd(jl,jk)*zcond(jl) prfl(jl)=prfl(jl)+pdmfdp(jl,jk-1) endif endif enddo enddo return end subroutine cudlfsn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- !********************************************** ! subroutine cuddrafn !********************************************** subroutine cuddrafn & & ( klon, klev, lddraf & & , ptenh, pqenh, puen, pven & & , pgeo, pgeoh, paph, prfl & & , ptd, pqd, pud, pvd, pmfu & & , pmfd, pmfds, pmfdq, pdmfdp, pmfdde_rate ) ! this routine calculates cumulus downdraft descent ! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89 ! purpose. ! -------- ! to produce the vertical profiles for cumulus downdrafts ! (i.e. t,q,u and v and fluxes) ! interface ! --------- ! this routine is called from *cumastr*. ! input is t,q,p,phi,u,v at half levels. ! it returns fluxes of s,q and evaporation rate ! and u,v at levels where downdraft occurs ! method. ! -------- ! calculate moist descent for entraining/detraining plume by ! a) moving air dry-adiabatically to next level below and ! b) correcting for evaporation to obtain saturated state. ! parameter description units ! --------- ----------- ----- ! input parameters (integer): ! *klon* number of grid points per packet ! *klev* number of levels ! input parameters (logical): ! *lddraf* .true. if downdrafts exist ! input parameters (real(kind=kind_phys)): ! *ptenh* env. temperature (t+1) on half levels k ! *pqenh* env. spec. humidity (t+1) on half levels kg/kg ! *puen* provisional environment u-velocity (t+1) m/s ! *pven* provisional environment v-velocity (t+1) m/s ! *pgeo* geopotential m2/s2 ! *pgeoh* geopotential on half levels m2/s2 ! *paph* provisional pressure on half levels pa ! *pmfu* massflux updrafts kg/(m2*s) ! updated parameters (real(kind=kind_phys)): ! *prfl* precipitation rate kg/(m2*s) ! output parameters (real(kind=kind_phys)): ! *ptd* temperature in downdrafts k ! *pqd* spec. humidity in downdrafts kg/kg ! *pud* u-velocity in downdrafts m/s ! *pvd* v-velocity in downdrafts m/s ! *pmfd* massflux in downdrafts kg/(m2*s) ! *pmfds* flux of dry static energy in downdrafts j/(m2*s) ! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s) ! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s) ! externals ! --------- ! *cuadjtq* for adjusting t and q due to evaporation in ! saturated descent !---------------------------------------------------------------------- implicit none integer klev,klon real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & & puen(klon,klev), pven(klon,klev), & & pgeoh(klon,klev+1), paph(klon,klev+1), & & pgeo(klon,klev), pmfu(klon,klev) real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), & & pud(klon,klev), pvd(klon,klev), & & pmfd(klon,klev), pmfds(klon,klev), & & pmfdq(klon,klev), pdmfdp(klon,klev), & & prfl(klon) real(kind=kind_phys) pmfdde_rate(klon,klev) logical lddraf(klon) real(kind=kind_phys) zdmfen(klon), zdmfde(klon), & & zcond(klon), zoentr(klon), & & zbuoy(klon) real(kind=kind_phys) zph(klon) logical llo2(klon) logical llo1 ! local variables integer jl,jk integer is,ik,icall,ike, itopde(klon) real(kind=kind_phys) zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp real(kind=kind_phys) zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk !---------------------------------------------------------------------- ! 1. calculate moist descent for cumulus downdraft by ! (a) calculating entrainment/detrainment rates, ! including organized entrainment dependent on ! negative buoyancy and assuming ! linear decrease of massflux in pbl ! (b) doing moist descent - evaporative cooling ! and moistening is calculated in *cuadjtq* ! (c) checking for negative buoyancy and ! specifying final t,q,u,v and downward fluxes ! ------------------------------------------------- do jl=1,klon zoentr(jl)=0. zbuoy(jl)=0. zdmfen(jl)=0. zdmfde(jl)=0. enddo do jk=klev,1,-1 do jl=1,klon pmfdde_rate(jl,jk) = 0. if((paph(jl,klev+1)-paph(jl,jk)).lt. 60.e2) itopde(jl) = jk end do end do do jk=3,klev is=0 do jl=1,klon zph(jl)=paph(jl,jk) llo2(jl)=lddraf(jl).and.pmfd(jl,jk-1).lt.0. if(llo2(jl)) then is=is+1 endif enddo if(is.eq.0) cycle do jl=1,klon if(llo2(jl)) then zentr = entrdd*pmfd(jl,jk-1)*(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg zdmfen(jl)=zentr zdmfde(jl)=zentr endif enddo do jl=1,klon if(llo2(jl)) then if(jk.gt.itopde(jl)) then zdmfen(jl)=0. zdmfde(jl)=pmfd(jl,itopde(jl))* & & (paph(jl,jk)-paph(jl,jk-1))/ & & (paph(jl,klev+1)-paph(jl,itopde(jl))) endif endif enddo do jl=1,klon if(llo2(jl)) then if(jk.le.itopde(jl)) then zdz=-(pgeoh(jl,jk-1)-pgeoh(jl,jk))*zrg zzentr=zoentr(jl)*zdz*pmfd(jl,jk-1) zdmfen(jl)=zdmfen(jl)+zzentr zdmfen(jl)=max(zdmfen(jl),0.3*pmfd(jl,jk-1)) zdmfen(jl)=max(zdmfen(jl),-0.75*pmfu(jl,jk)- & & (pmfd(jl,jk-1)-zdmfde(jl))) zdmfen(jl)=min(zdmfen(jl),0.) endif endif enddo do jl=1,klon if(llo2(jl)) then pmfd(jl,jk)=pmfd(jl,jk-1)+zdmfen(jl)-zdmfde(jl) zseen=(cpd*ptenh(jl,jk-1)+pgeoh(jl,jk-1))*zdmfen(jl) zqeen=pqenh(jl,jk-1)*zdmfen(jl) zsdde=(cpd*ptd(jl,jk-1)+pgeoh(jl,jk-1))*zdmfde(jl) zqdde=pqd(jl,jk-1)*zdmfde(jl) zmfdsk=pmfds(jl,jk-1)+zseen-zsdde zmfdqk=pmfdq(jl,jk-1)+zqeen-zqdde pqd(jl,jk)=zmfdqk*(1./min(-cmfcmin,pmfd(jl,jk))) ptd(jl,jk)=(zmfdsk*(1./min(-cmfcmin,pmfd(jl,jk)))-& & pgeoh(jl,jk))*rcpd ptd(jl,jk)=min(400.,ptd(jl,jk)) ptd(jl,jk)=max(100.,ptd(jl,jk)) zcond(jl)=pqd(jl,jk) endif enddo ik=jk icall=2 call cuadjtqn(klon, klev, ik, zph, ptd, pqd, llo2, icall ) do jl=1,klon if(llo2(jl)) then zcond(jl)=zcond(jl)-pqd(jl,jk) zbuo=ptd(jl,jk)*(1.+vtmpc1 *pqd(jl,jk))- & & ptenh(jl,jk)*(1.+vtmpc1 *pqenh(jl,jk)) if(prfl(jl).gt.0..and.pmfu(jl,jk).gt.0.) then zrain=prfl(jl)/pmfu(jl,jk) zbuo=zbuo-ptd(jl,jk)*zrain endif if(zbuo.ge.0 .or. prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then pmfd(jl,jk)=0. zbuo=0. endif pmfds(jl,jk)=(cpd*ptd(jl,jk)+pgeoh(jl,jk))*pmfd(jl,jk) pmfdq(jl,jk)=pqd(jl,jk)*pmfd(jl,jk) zdmfdp=-pmfd(jl,jk)*zcond(jl) pdmfdp(jl,jk-1)=zdmfdp prfl(jl)=prfl(jl)+zdmfdp ! compute organized entrainment for use at next level zbuoyz=zbuo/ptenh(jl,jk) zbuoyz=min(zbuoyz,0.0) zdz=-(pgeo(jl,jk-1)-pgeo(jl,jk)) zbuoy(jl)=zbuoy(jl)+zbuoyz*zdz zoentr(jl)=g*zbuoyz*0.5/(1.+zbuoy(jl)) pmfdde_rate(jl,jk) = -zdmfde(jl) endif enddo enddo return end subroutine cuddrafn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cuflxn & & ( klon, klev, ztmst & & , pten, pqen, pqsen, ptenh, pqenh & & , paph, pap, pgeoh, lndj, ldcum & & , kcbot, kctop, kdtop, ktopm2 & & , ktype, lddraf & & , pmfu, pmfd, pmfus, pmfds & & , pmfuq, pmfdq, pmful, plude & & , pdmfup, pdmfdp, pdpmel, plglac & & , prain, pmfdde_rate, pmflxr, pmflxs ) ! m.tiedtke e.c.m.w.f. 7/86 modif. 12/89 ! purpose ! ------- ! this routine does the final calculation of convective ! fluxes in the cloud layer and in the subcloud layer ! interface ! --------- ! this routine is called from *cumastr*. ! parameter description units ! --------- ----------- ----- ! input parameters (integer): ! *klon* number of grid points per packet ! *klev* number of levels ! *kcbot* cloud base level ! *kctop* cloud top level ! *kdtop* top level of downdrafts ! input parameters (logical): ! *lndj* land sea mask (1 for land) ! *ldcum* flag: .true. for convective points ! input parameters (real(kind=kind_phys)): ! *ztmst* time step for the physics s ! *pten* provisional environment temperature (t+1) k ! *pqen* provisional environment spec. humidity (t+1) kg/kg ! *pqsen* environment spec. saturation humidity (t+1) kg/kg ! *ptenh* env. temperature (t+1) on half levels k ! *pqenh* env. spec. humidity (t+1) on half levels kg/kg ! *paph* provisional pressure on half levels pa ! *pap* provisional pressure on full levels pa ! *pgeoh* geopotential on half levels m2/s2 ! updated parameters (integer): ! *ktype* set to zero if ldcum=.false. ! updated parameters (logical): ! *lddraf* set to .false. if ldcum=.false. or kdtop= kdtop(jl) if ( llddraf .and.jk.ge.kdtop(jl)) then pmfds(jl,jk) = pmfds(jl,jk)-pmfd(jl,jk) * & (cpd*ptenh(jl,jk)+pgeoh(jl,jk)) pmfdq(jl,jk) = pmfdq(jl,jk)-pmfd(jl,jk)*pqenh(jl,jk) else pmfd(jl,jk) = 0. pmfds(jl,jk) = 0. pmfdq(jl,jk) = 0. pdmfdp(jl,jk-1) = 0. end if if ( llddraf .and. pmfd(jl,jk) < 0. .and. & abs(pmfd(jl,ikb)) < 1.e-20 ) then idbas(jl) = jk end if else pmfu(jl,jk)=0. pmfd(jl,jk)=0. pmfus(jl,jk)=0. pmfds(jl,jk)=0. pmfuq(jl,jk)=0. pmfdq(jl,jk)=0. pmful(jl,jk)=0. plglac(jl,jk)=0. pdmfup(jl,jk-1)=0. pdmfdp(jl,jk-1)=0. plude(jl,jk-1)=0. endif enddo enddo do jl=1,klon pmflxr(jl,klev+1) = 0. pmflxs(jl,klev+1) = 0. end do do jl=1,klon if(ldcum(jl)) then ikb=kcbot(jl) ik=ikb+1 zzp=((paph(jl,klev+1)-paph(jl,ik))/ & & (paph(jl,klev+1)-paph(jl,ikb))) if(ktype(jl).eq.3) then zzp=zzp**2 endif pmfu(jl,ik)=pmfu(jl,ikb)*zzp pmfus(jl,ik)=(pmfus(jl,ikb)- & & foelhm(ptenh(jl,ikb))*pmful(jl,ikb))*zzp pmfuq(jl,ik)=(pmfuq(jl,ikb)+pmful(jl,ikb))*zzp pmful(jl,ik)=0. endif enddo do jk=ktopm2,klev do jl=1,klon if(ldcum(jl).and.jk.gt.kcbot(jl)+1) then ikb=kcbot(jl)+1 zzp=((paph(jl,klev+1)-paph(jl,jk))/ & & (paph(jl,klev+1)-paph(jl,ikb))) if(ktype(jl).eq.3) then zzp=zzp**2 endif pmfu(jl,jk)=pmfu(jl,ikb)*zzp pmfus(jl,jk)=pmfus(jl,ikb)*zzp pmfuq(jl,jk)=pmfuq(jl,ikb)*zzp pmful(jl,jk)=0. endif ik = idbas(jl) llddraf = lddraf(jl) .and. jk > ik .and. ik < klev if ( llddraf .and. ik == kcbot(jl)+1 ) then zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ik))) if ( ktype(jl) == 3 ) zzp = zzp*zzp pmfd(jl,jk) = pmfd(jl,ik)*zzp pmfds(jl,jk) = pmfds(jl,ik)*zzp pmfdq(jl,jk) = pmfdq(jl,ik)*zzp pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk)) else if ( llddraf .and. ik /= kcbot(jl)+1 .and. jk == ik+1 ) then pmfdde_rate(jl,jk) = -(pmfd(jl,jk-1)-pmfd(jl,jk)) end if enddo enddo !* 2. calculate rain/snow fall rates !* calculate melting of snow !* calculate evaporation of precip ! ------------------------------- do jk=ktopm2,klev do jl=1,klon if(ldcum(jl) .and. jk >=kctop(jl)-1 ) then prain(jl)=prain(jl)+pdmfup(jl,jk) if(pmflxs(jl,jk).gt.0..and.pten(jl,jk).gt.tmelt) then zcons1=zcons1a*(1.+0.5*(pten(jl,jk)-tmelt)) zfac=zcons1*(paph(jl,jk+1)-paph(jl,jk)) zsnmlt=min(pmflxs(jl,jk),zfac*(pten(jl,jk)-tmelt)) pdpmel(jl,jk)=zsnmlt pqsen(jl,jk)=foeewm(pten(jl,jk)-zsnmlt/zfac)/pap(jl,jk) endif zalfaw=foealfa(pten(jl,jk)) ! ! No liquid precipitation above melting level ! if ( pten(jl,jk) < tmelt .and. zalfaw > 0. ) then plglac(jl,jk) = plglac(jl,jk)+zalfaw*(pdmfup(jl,jk)+pdmfdp(jl,jk)) zalfaw = 0. end if pmflxr(jl,jk+1)=pmflxr(jl,jk)+zalfaw* & & (pdmfup(jl,jk)+pdmfdp(jl,jk))+pdpmel(jl,jk) pmflxs(jl,jk+1)=pmflxs(jl,jk)+(1.-zalfaw)* & & (pdmfup(jl,jk)+pdmfdp(jl,jk))-pdpmel(jl,jk) if(pmflxr(jl,jk+1)+pmflxs(jl,jk+1).lt.0.0) then pdmfdp(jl,jk)=-(pmflxr(jl,jk)+pmflxs(jl,jk)+pdmfup(jl,jk)) pmflxr(jl,jk+1)=0.0 pmflxs(jl,jk+1)=0.0 pdpmel(jl,jk) =0.0 else if ( pmflxr(jl,jk+1) < 0. ) then pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1) pmflxr(jl,jk+1) = 0. else if ( pmflxs(jl,jk+1) < 0. ) then pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1) pmflxs(jl,jk+1) = 0. end if endif enddo enddo do jk=ktopm2,klev do jl=1,klon if(ldcum(jl).and.jk.ge.kcbot(jl)) then zrfl=pmflxr(jl,jk)+pmflxs(jl,jk) if(zrfl.gt.1.e-20) then zdrfl1=zcpecons*max(0.,pqsen(jl,jk)-pqen(jl,jk))*zcucov* & & (sqrt(paph(jl,jk)/paph(jl,klev+1))/5.09e-3* & & zrfl/zcucov)**0.5777* & & (paph(jl,jk+1)-paph(jl,jk)) zrnew=zrfl-zdrfl1 zrmin=zrfl-zcucov*max(0.,rhevap(jl)*pqsen(jl,jk) & & -pqen(jl,jk)) *zcons2*(paph(jl,jk+1)-paph(jl,jk)) zrnew=max(zrnew,zrmin) zrfln=max(zrnew,0.) zdrfl=min(0.,zrfln-zrfl) zdenom=1./max(1.e-20,pmflxr(jl,jk)+pmflxs(jl,jk)) zalfaw=foealfa(pten(jl,jk)) if ( pten(jl,jk) < tmelt ) zalfaw = 0. zpdr=zalfaw*pdmfdp(jl,jk) zpds=(1.-zalfaw)*pdmfdp(jl,jk) pmflxr(jl,jk+1)=pmflxr(jl,jk)+zpdr & & +pdpmel(jl,jk)+zdrfl*pmflxr(jl,jk)*zdenom pmflxs(jl,jk+1)=pmflxs(jl,jk)+zpds & & -pdpmel(jl,jk)+zdrfl*pmflxs(jl,jk)*zdenom pdmfup(jl,jk)=pdmfup(jl,jk)+zdrfl if ( pmflxr(jl,jk+1)+pmflxs(jl,jk+1) < 0. ) then pdmfup(jl,jk) = pdmfup(jl,jk)-(pmflxr(jl,jk+1)+pmflxs(jl,jk+1)) pmflxr(jl,jk+1) = 0. pmflxs(jl,jk+1) = 0. pdpmel(jl,jk) = 0. else if ( pmflxr(jl,jk+1) < 0. ) then pmflxs(jl,jk+1) = pmflxs(jl,jk+1)+pmflxr(jl,jk+1) pmflxr(jl,jk+1) = 0. else if ( pmflxs(jl,jk+1) < 0. ) then pmflxr(jl,jk+1) = pmflxr(jl,jk+1)+pmflxs(jl,jk+1) pmflxs(jl,jk+1) = 0. end if else pmflxr(jl,jk+1)=0.0 pmflxs(jl,jk+1)=0.0 pdmfdp(jl,jk)=0.0 pdpmel(jl,jk)=0.0 endif endif enddo enddo return end subroutine cuflxn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cudtdqn(klon,klev,ktopm2,kctop,kdtop,ldcum, & lddraf,ztmst,paph,pgeoh,pgeo,pten,ptenh,pqen, & pqenh,pqsen,plglac,plude,pmfu,pmfd,pmfus,pmfds, & pmfuq,pmfdq,pmful,pdmfup,pdmfdp,pdpmel,ptent,ptenq,pcte) implicit none integer klon,klev,ktopm2 integer kctop(klon), kdtop(klon) logical ldcum(klon), lddraf(klon) real(kind=kind_phys) ztmst real(kind=kind_phys) paph(klon,klev+1), pgeoh(klon,klev+1) real(kind=kind_phys) pgeo(klon,klev), pten(klon,klev), & pqen(klon,klev), ptenh(klon,klev),& pqenh(klon,klev), pqsen(klon,klev),& plglac(klon,klev), plude(klon,klev) real(kind=kind_phys) pmfu(klon,klev), pmfd(klon,klev),& pmfus(klon,klev), pmfds(klon,klev),& pmfuq(klon,klev), pmfdq(klon,klev),& pmful(klon,klev), pdmfup(klon,klev),& pdpmel(klon,klev), pdmfdp(klon,klev) real(kind=kind_phys) ptent(klon,klev), ptenq(klon,klev) real(kind=kind_phys) pcte(klon,klev) ! local variables integer jk , ik , jl real(kind=kind_phys) zalv , zzp real(kind=kind_phys) zdtdt(klon,klev) , zdqdt(klon,klev) , zdp(klon,klev) !* 1.0 SETUP AND INITIALIZATIONS ! ------------------------- do jk = 1 , klev do jl = 1, klon if ( ldcum(jl) ) then zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk)) end if end do end do !----------------------------------------------------------------------- !* 2.0 COMPUTE TENDENCIES ! ------------------ do jk = ktopm2 , klev if ( jk < klev ) then do jl = 1,klon if ( ldcum(jl) ) then zalv = foelhm(pten(jl,jk)) zdtdt(jl,jk) = zdp(jl,jk)*rcpd * & (pmfus(jl,jk+1)-pmfus(jl,jk)+pmfds(jl,jk+1) - & pmfds(jl,jk)+alf*plglac(jl,jk)-alf*pdpmel(jl,jk) - & zalv*(pmful(jl,jk+1)-pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk))) zdqdt(jl,jk) = zdp(jl,jk)*(pmfuq(jl,jk+1) - & pmfuq(jl,jk)+pmfdq(jl,jk+1)-pmfdq(jl,jk)+pmful(jl,jk+1) - & pmful(jl,jk)-plude(jl,jk)-pdmfup(jl,jk)-pdmfdp(jl,jk)) end if end do else do jl = 1,klon if ( ldcum(jl) ) then zalv = foelhm(pten(jl,jk)) zdtdt(jl,jk) = -zdp(jl,jk)*rcpd * & (pmfus(jl,jk)+pmfds(jl,jk)+alf*pdpmel(jl,jk) - & zalv*(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk)+plude(jl,jk))) zdqdt(jl,jk) = -zdp(jl,jk)*(pmfuq(jl,jk) + plude(jl,jk) + & pmfdq(jl,jk)+(pmful(jl,jk)+pdmfup(jl,jk)+pdmfdp(jl,jk))) end if end do end if end do !--------------------------------------------------------------- !* 3.0 UPDATE TENDENCIES ! ----------------- do jk = ktopm2 , klev do jl = 1, klon if ( ldcum(jl) ) then ptent(jl,jk) = ptent(jl,jk) + zdtdt(jl,jk) ptenq(jl,jk) = ptenq(jl,jk) + zdqdt(jl,jk) pcte(jl,jk) = zdp(jl,jk)*plude(jl,jk) end if end do end do return end subroutine cudtdqn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cududvn(klon,klev,ktopm2,ktype,kcbot,kctop,ldcum, & ztmst,paph,puen,pven,pmfu,pmfd,puu,pud,pvu,pvd,ptenu, & ptenv) implicit none integer klon,klev,ktopm2 integer ktype(klon), kcbot(klon), kctop(klon) logical ldcum(klon) real(kind=kind_phys) ztmst real(kind=kind_phys) paph(klon,klev+1) real(kind=kind_phys) puen(klon,klev), pven(klon,klev),& pmfu(klon,klev), pmfd(klon,klev),& puu(klon,klev), pud(klon,klev),& pvu(klon,klev), pvd(klon,klev) real(kind=kind_phys) ptenu(klon,klev), ptenv(klon,klev) !local variables real(kind=kind_phys) zuen(klon,klev) , zven(klon,klev) , zmfuu(klon,klev), & zmfdu(klon,klev), zmfuv(klon,klev), zmfdv(klon,klev) integer ik , ikb , jk , jl real(kind=kind_phys) zzp, zdtdt real(kind=kind_phys) zdudt(klon,klev), zdvdt(klon,klev), zdp(klon,klev) ! do jk = 1 , klev do jl = 1, klon if ( ldcum(jl) ) then zuen(jl,jk) = puen(jl,jk) zven(jl,jk) = pven(jl,jk) zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk)) end if end do end do !---------------------------------------------------------------------- !* 1.0 CALCULATE FLUXES AND UPDATE U AND V TENDENCIES ! ---------------------------------------------- do jk = ktopm2 , klev ik = jk - 1 do jl = 1,klon if ( ldcum(jl) ) then zmfuu(jl,jk) = pmfu(jl,jk)*(puu(jl,jk)-zuen(jl,ik)) zmfuv(jl,jk) = pmfu(jl,jk)*(pvu(jl,jk)-zven(jl,ik)) zmfdu(jl,jk) = pmfd(jl,jk)*(pud(jl,jk)-zuen(jl,ik)) zmfdv(jl,jk) = pmfd(jl,jk)*(pvd(jl,jk)-zven(jl,ik)) end if end do end do ! linear fluxes below cloud do jk = ktopm2 , klev do jl = 1, klon if ( ldcum(jl) .and. jk > kcbot(jl) ) then ikb = kcbot(jl) zzp = ((paph(jl,klev+1)-paph(jl,jk))/(paph(jl,klev+1)-paph(jl,ikb))) if ( ktype(jl) == 3 ) zzp = zzp*zzp zmfuu(jl,jk) = zmfuu(jl,ikb)*zzp zmfuv(jl,jk) = zmfuv(jl,ikb)*zzp zmfdu(jl,jk) = zmfdu(jl,ikb)*zzp zmfdv(jl,jk) = zmfdv(jl,ikb)*zzp end if end do end do !---------------------------------------------------------------------- !* 2.0 COMPUTE TENDENCIES ! ------------------ do jk = ktopm2 , klev if ( jk < klev ) then ik = jk + 1 do jl = 1,klon if ( ldcum(jl) ) then zdudt(jl,jk) = zdp(jl,jk) * & (zmfuu(jl,ik)-zmfuu(jl,jk)+zmfdu(jl,ik)-zmfdu(jl,jk)) zdvdt(jl,jk) = zdp(jl,jk) * & (zmfuv(jl,ik)-zmfuv(jl,jk)+zmfdv(jl,ik)-zmfdv(jl,jk)) end if end do else do jl = 1,klon if ( ldcum(jl) ) then zdudt(jl,jk) = -zdp(jl,jk)*(zmfuu(jl,jk)+zmfdu(jl,jk)) zdvdt(jl,jk) = -zdp(jl,jk)*(zmfuv(jl,jk)+zmfdv(jl,jk)) end if end do end if end do !--------------------------------------------------------------------- !* 3.0 UPDATE TENDENCIES ! ----------------- do jk = ktopm2 , klev do jl = 1, klon if ( ldcum(jl) ) then ptenu(jl,jk) = ptenu(jl,jk) + zdudt(jl,jk) ptenv(jl,jk) = ptenv(jl,jk) + zdvdt(jl,jk) end if end do end do !---------------------------------------------------------------------- return end subroutine cududvn !--------------------------------------------------------- ! level 3 souroutines !-------------------------------------------------------- subroutine cuctracer(klon,klev,ktrac,kctop,kdtop, & ldcum,lddraf,ztmst,paph,pmfu,pmfd, & pudrate,pddrate,pcen,ptenc) implicit none integer klon,klev,ktrac integer kctop(klon), kdtop(klon) logical ldcum(klon), lddraf(klon) real(kind=kind_phys) ztmst real(kind=kind_phys) paph(klon,klev+1) real(kind=kind_phys) pmfu(klon,klev) real(kind=kind_phys) pmfd(klon,klev) real(kind=kind_phys) pudrate(klon,klev) real(kind=kind_phys) pddrate(klon,klev) real(kind=kind_phys) pcen(klon,klev,ktrac) real(kind=kind_phys) ptenc(klon,klev,ktrac) !---------------------------------------------------------------------- integer ik , jk , jl , jn real(kind=kind_phys) zzp , zmfa , zerate , zposi ! ALLOCATABLE ARAYS real(kind=kind_phys) , dimension(:,:,:) , allocatable :: zcen , zcu , zcd , & ztenc , zmfc real(kind=kind_phys) , dimension(:,:) , allocatable :: zdp logical , dimension(:,:) , allocatable :: llcumask , llcumbas !---------------------------------------------------------------------- allocate (zcen(klon,klev,ktrac)) ! Half-level environmental values allocate (zcu(klon,klev,ktrac)) ! Updraft values allocate (zcd(klon,klev,ktrac)) ! Downdraft values allocate (ztenc(klon,klev,ktrac)) ! Tendency allocate (zmfc(klon,klev,ktrac)) ! Fluxes allocate (zdp(klon,klev)) ! Pressure difference allocate (llcumask(klon,klev)) ! Mask for convection ! Initialize Cumulus mask + some setups do jk = 2, klev do jl = 1, klon llcumask(jl,jk) = .false. if ( ldcum(jl) ) then zdp(jl,jk) = g/(paph(jl,jk+1)-paph(jl,jk)) if ( jk >= kctop(jl)-1 ) llcumask(jl,jk) = .true. end if end do end do !---------------------------------------------------------------------- do jn = 1 , ktrac !* 1.0 DEFINE TRACERS AT HALF LEVELS ! ----------------------------- do jk = 2 , klev ik = jk - 1 do jl = 1, klon zcen(jl,jk,jn) = pcen(jl,jk,jn) zcd(jl,jk,jn) = pcen(jl,ik,jn) zcu(jl,jk,jn) = pcen(jl,ik,jn) zmfc(jl,jk,jn) = 0. ztenc(jl,jk,jn)= 0. end do end do do jl = 1, klon zcu(jl,klev,jn) = pcen(jl,klev,jn) end do !* 2.0 COMPUTE UPDRAFT VALUES ! ---------------------- do jk = klev - 1 , 3 , -1 ik = jk + 1 do jl = 1, klon if ( llcumask(jl,jk) ) then zerate = pmfu(jl,jk) - pmfu(jl,ik) + pudrate(jl,jk) zmfa = 1./max(cmfcmin,pmfu(jl,jk)) if ( jk >= kctop(jl) ) then zcu(jl,jk,jn) = (pmfu(jl,ik)*zcu(jl,ik,jn) + & zerate*pcen(jl,jk,jn)-pudrate(jl,jk)*zcu(jl,ik,jn))*zmfa end if end if end do end do !* 3.0 COMPUTE DOWNDRAFT VALUES ! ------------------------ do jk = 3 , klev ik = jk - 1 do jl = 1, klon if ( lddraf(jl) .and. jk == kdtop(jl) ) then ! Nota: in order to avoid final negative Tracer values at LFS ! the allowed value of ZCD depends on the jump in mass flux ! at the LFS zcd(jl,jk,jn) = 0.1*zcu(jl,jk,jn) + 0.9*pcen(jl,ik,jn) else if ( lddraf(jl).and.jk>kdtop(jl) ) then zerate = -pmfd(jl,jk) + pmfd(jl,ik) + pddrate(jl,jk) zmfa = 1./min(-cmfcmin,pmfd(jl,jk)) zcd(jl,jk,jn) = (pmfd(jl,ik)*zcd(jl,ik,jn) - & zerate*pcen(jl,ik,jn)+pddrate(jl,jk)*zcd(jl,ik,jn))*zmfa end if end do end do ! In order to avoid negative Tracer at KLEV adjust ZCD jk = klev ik = jk - 1 do jl = 1, klon if ( lddraf(jl) ) then zposi = -zdp(jl,jk) *(pmfu(jl,jk)*zcu(jl,jk,jn) + & pmfd(jl,jk)*zcd(jl,jk,jn)-(pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn)) if ( pcen(jl,jk,jn)+zposi*ztmst < 0. ) then zmfa = 1./min(-cmfcmin,pmfd(jl,jk)) zcd(jl,jk,jn) = ((pmfu(jl,jk)+pmfd(jl,jk))*pcen(jl,ik,jn) - & pmfu(jl,jk)*zcu(jl,jk,jn)+pcen(jl,jk,jn) / & (ztmst*zdp(jl,jk)))*zmfa end if end if end do end do !---------------------------------------------------------------------- do jn = 1 , ktrac !* 4.0 COMPUTE FLUXES ! -------------- do jk = 2 , klev ik = jk - 1 do jl = 1, klon if ( llcumask(jl,jk) ) then zmfa = pmfu(jl,jk) + pmfd(jl,jk) zmfc(jl,jk,jn) = pmfu(jl,jk)*zcu(jl,jk,jn) + & pmfd(jl,jk)*zcd(jl,jk,jn) - zmfa*zcen(jl,ik,jn) end if end do end do !* 5.0 COMPUTE TENDENCIES = RHS ! ------------------------ do jk = 2 , klev - 1 ik = jk + 1 do jl = 1, klon if ( llcumask(jl,jk) ) then ztenc(jl,jk,jn) = zdp(jl,jk)*(zmfc(jl,ik,jn)-zmfc(jl,jk,jn)) end if end do end do jk = klev do jl = 1, klon if ( ldcum(jl) ) ztenc(jl,jk,jn) = -zdp(jl,jk)*zmfc(jl,jk,jn) end do end do !* 6.0 UPDATE TENDENCIES ! ----------------- do jn = 1, ktrac do jk = 2, klev do jl = 1, klon if ( llcumask(jl,jk) ) then ptenc(jl,jk,jn) = ptenc(jl,jk,jn)+ztenc(jl,jk,jn) end if end do end do end do !--------------------------------------------------------------------------- deallocate (llcumask) deallocate (zdp) deallocate (zmfc) deallocate (ztenc) deallocate (zcd) deallocate (zcu) deallocate (zcen) end subroutine cuctracer !--------------------------------------------------------- ! level 4 souroutines !-------------------------------------------------------- subroutine cuadjtqn & & (klon, klev, kk, psp, pt, pq, ldflag, kcall) ! m.tiedtke e.c.m.w.f. 12/89 ! purpose. ! -------- ! to produce t,q and l values for cloud ascent ! interface ! --------- ! this routine is called from subroutines: ! *cond* (t and q at condensation level) ! *cubase* (t and q at condensation level) ! *cuasc* (t and q at cloud levels) ! *cuini* (environmental t and qs values at half levels) ! input are unadjusted t and q values, ! it returns adjusted values of t and q ! parameter description units ! --------- ----------- ----- ! input parameters (integer): ! *klon* number of grid points per packet ! *klev* number of levels ! *kk* level ! *kcall* defines calculation as ! kcall=0 env. t and qs in*cuini* ! kcall=1 condensation in updrafts (e.g. cubase, cuasc) ! kcall=2 evaporation in downdrafts (e.g. cudlfs,cuddraf) ! input parameters (real(kind=kind_phys)): ! *psp* pressure pa ! updated parameters (real(kind=kind_phys)): ! *pt* temperature k ! *pq* specific humidity kg/kg ! externals ! --------- ! for condensation calculations. ! the tables are initialised in *suphec*. !---------------------------------------------------------------------- implicit none integer klev,klon real(kind=kind_phys) pt(klon,klev), pq(klon,klev), & & psp(klon) logical ldflag(klon) ! local variables integer jl,jk integer isum,kcall,kk real(kind=kind_phys) zqmax,zqsat,zcor,zqp,zcond,zcond1,zl,zi,zf !---------------------------------------------------------------------- ! 1. define constants ! ---------------- zqmax=0.5 ! 2. calculate condensation and adjust t and q accordingly ! ----------------------------------------------------- if ( kcall == 1 ) then do jl = 1,klon if ( ldflag(jl) ) then zqp = 1./psp(jl) zl = 1./(pt(jl,kk)-c4les) zi = 1./(pt(jl,kk)-c4ies) zqsat = c2es*(foealfa(pt(jl,kk))*exp(c3les*(pt(jl,kk)-tmelt)*zl) + & (1.-foealfa(pt(jl,kk)))*exp(c3ies*(pt(jl,kk)-tmelt)*zi)) zqsat = zqsat*zqp zqsat = min(0.5,zqsat) zcor = 1. - vtmpc1*zqsat zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + & (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2 zcond = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf) if ( zcond > 0. ) then pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond pq(jl,kk) = pq(jl,kk) - zcond zl = 1./(pt(jl,kk)-c4les) zi = 1./(pt(jl,kk)-c4ies) zqsat = c2es*(foealfa(pt(jl,kk)) * & exp(c3les*(pt(jl,kk)-tmelt)*zl)+(1.-foealfa(pt(jl,kk))) * & exp(c3ies*(pt(jl,kk)-tmelt)*zi)) zqsat = zqsat*zqp zqsat = min(0.5,zqsat) zcor = 1. - vtmpc1*zqsat zf = foealfa(pt(jl,kk))*r5alvcp*zl**2 + & (1.-foealfa(pt(jl,kk)))*r5alscp*zi**2 zcond1 = (pq(jl,kk)*zcor**2-zqsat*zcor)/(zcor**2+zqsat*zf) if ( abs(zcond) < 1.e-20 ) zcond1 = 0. pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1 pq(jl,kk) = pq(jl,kk) - zcond1 end if end if end do elseif ( kcall == 2 ) then do jl = 1,klon if ( ldflag(jl) ) then zqp = 1./psp(jl) zqsat = foeewm(pt(jl,kk))*zqp zqsat = min(0.5,zqsat) zcor = 1./(1.-vtmpc1*zqsat) zqsat = zqsat*zcor zcond = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk))) zcond = min(zcond,0.) pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond pq(jl,kk) = pq(jl,kk) - zcond zqsat = foeewm(pt(jl,kk))*zqp zqsat = min(0.5,zqsat) zcor = 1./(1.-vtmpc1*zqsat) zqsat = zqsat*zcor zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk))) if ( abs(zcond) < 1.e-20 ) zcond1 = min(zcond1,0.) pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1 pq(jl,kk) = pq(jl,kk) - zcond1 end if end do else if ( kcall == 0 ) then do jl = 1,klon zqp = 1./psp(jl) zqsat = foeewm(pt(jl,kk))*zqp zqsat = min(0.5,zqsat) zcor = 1./(1.-vtmpc1*zqsat) zqsat = zqsat*zcor zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk))) pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1 pq(jl,kk) = pq(jl,kk) - zcond1 zqsat = foeewm(pt(jl,kk))*zqp zqsat = min(0.5,zqsat) zcor = 1./(1.-vtmpc1*zqsat) zqsat = zqsat*zcor zcond1 = (pq(jl,kk)-zqsat)/(1.+zqsat*zcor*foedem(pt(jl,kk))) pt(jl,kk) = pt(jl,kk) + foeldcpm(pt(jl,kk))*zcond1 pq(jl,kk) = pq(jl,kk) - zcond1 end do end if return end subroutine cuadjtqn !--------------------------------------------------------- ! level 4 souroutines !-------------------------------------------------------- subroutine cubasmcn & & (klon, klev, klevm1, kk, pten,& & pqen, pqsen, puen, pven, pverv,& & pgeo, pgeoh, ldcum, ktype, klab, plrain,& & pmfu, pmfub, kcbot, ptu,& & pqu, plu, puu, pvu, pmfus,& & pmfuq, pmful, pdmfup) implicit none ! m.tiedtke e.c.m.w.f. 12/89 ! c.zhang iprc 05/2012 !***purpose. ! -------- ! this routine calculates cloud base values ! for midlevel convection !***interface ! --------- ! this routine is called from *cuasc*. ! input are environmental values t,q etc ! it returns cloudbase values for midlevel convection !***method. ! ------- ! s. tiedtke (1989) !***externals ! --------- ! none ! ---------------------------------------------------------------- real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),& & puen(klon,klev), pven(klon,klev),& & pqsen(klon,klev), pverv(klon,klev),& & pgeo(klon,klev), pgeoh(klon,klev+1) real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& & puu(klon,klev), pvu(klon,klev),& & plu(klon,klev), pmfu(klon,klev),& & pmfub(klon), & & pmfus(klon,klev), pmfuq(klon,klev),& & pmful(klon,klev), pdmfup(klon,klev),& & plrain(klon,klev) integer ktype(klon), kcbot(klon),& & klab(klon,klev) logical ldcum(klon) ! local variabels integer jl,kk,klev,klon,klevp1,klevm1 real(kind=kind_phys) zzzmb !-------------------------------------------------------- !* 1. calculate entrainment and detrainment rates ! ------------------------------------------------------- do jl=1,klon if(.not.ldcum(jl) .and. klab(jl,kk+1).eq.0) then if(lmfmid .and. pqen(jl,kk) .gt. 0.80*pqsen(jl,kk).and. & pgeo(jl,kk)*zrg .gt. 5.0e2 .and. & & pgeo(jl,kk)*zrg .lt. 1.0e4 ) then ptu(jl,kk+1)=(cpd*pten(jl,kk)+pgeo(jl,kk)-pgeoh(jl,kk+1))& & *rcpd pqu(jl,kk+1)=pqen(jl,kk) plu(jl,kk+1)=0. zzzmb=max(cmfcmin,-pverv(jl,kk)*zrg) zzzmb=min(zzzmb,cmfcmax) pmfub(jl)=zzzmb pmfu(jl,kk+1)=pmfub(jl) pmfus(jl,kk+1)=pmfub(jl)*(cpd*ptu(jl,kk+1)+pgeoh(jl,kk+1)) pmfuq(jl,kk+1)=pmfub(jl)*pqu(jl,kk+1) pmful(jl,kk+1)=0. pdmfup(jl,kk+1)=0. kcbot(jl)=kk klab(jl,kk+1)=1 plrain(jl,kk+1)=0.0 ktype(jl)=3 end if end if end do return end subroutine cubasmcn ! !--------------------------------------------------------- ! level 4 souroutines !--------------------------------------------------------- subroutine cuentrn(klon,klev,kk,kcbot,ldcum,ldwork, & pgeoh,pmfu,pdmfen,pdmfde) implicit none integer klon,klev,kk integer kcbot(klon) logical ldcum(klon) logical ldwork real(kind=kind_phys) pgeoh(klon,klev+1) real(kind=kind_phys) pmfu(klon,klev) real(kind=kind_phys) pdmfen(klon) real(kind=kind_phys) pdmfde(klon) logical llo1 integer jl real(kind=kind_phys) zdz , zmf real(kind=kind_phys) zentr(klon) ! !* 1. CALCULATE ENTRAINMENT AND DETRAINMENT RATES ! ------------------------------------------- if ( ldwork ) then do jl = 1,klon pdmfen(jl) = 0. pdmfde(jl) = 0. zentr(jl) = 0. end do ! !* 1.1 SPECIFY ENTRAINMENT RATES ! ------------------------- do jl = 1, klon if ( ldcum(jl) ) then zdz = (pgeoh(jl,kk)-pgeoh(jl,kk+1))*zrg zmf = pmfu(jl,kk+1)*zdz llo1 = kk < kcbot(jl) if ( llo1 ) then pdmfen(jl) = zentr(jl)*zmf pdmfde(jl) = 0.75e-4*zmf end if end if end do end if end subroutine cuentrn ! !-------------------------------------------------------- ! external functions !------------------------------------------------------ real(kind=kind_phys) function foealfa(tt) ! foealfa is calculated to distinguish the three cases: ! ! foealfa=1 water phase ! foealfa=0 ice phase ! 0 < foealfa < 1 mixed phase ! ! input : tt = temperature ! implicit none real(kind=kind_phys) tt foealfa = min(1.,((max(rtice,min(rtwat,tt))-rtice) & & /(rtwat-rtice))**2) return end function foealfa real(kind=kind_phys) function foelhm(tt) implicit none real(kind=kind_phys) tt foelhm = foealfa(tt)*alv + (1.-foealfa(tt))*als return end function foelhm real(kind=kind_phys) function foeewm(tt) implicit none real(kind=kind_phys) tt foeewm = c2es * & & (foealfa(tt)*exp(c3les*(tt-tmelt)/(tt-c4les))+ & & (1.-foealfa(tt))*exp(c3ies*(tt-tmelt)/(tt-c4ies))) return end function foeewm real(kind=kind_phys) function foedem(tt) implicit none real(kind=kind_phys) tt foedem = foealfa(tt)*r5alvcp*(1./(tt-c4les)**2)+ & & (1.-foealfa(tt))*r5alscp*(1./(tt-c4ies)**2) return end function foedem real(kind=kind_phys) function foeldcpm(tt) implicit none real(kind=kind_phys) tt foeldcpm = foealfa(tt)*ralvdcp+ & & (1.-foealfa(tt))*ralsdcp return end function foeldcpm real(kind=kind_phys) function foeldcp(tt) implicit none real(kind=kind_phys) tt foeldcp = foedelta(tt)*ralvdcp + (1.-foedelta(tt))*ralsdcp end function foeldcp real(kind=kind_phys) function foedelta(tt) implicit none real(kind=kind_phys) tt foedelta = max(0.,sign(1.,tt-tmelt)) end function foedelta end module cu_ntiedtke