!================================================================================================================= module cu_ntiedtke_common use ccpp_kind_types,only: kind_phys implicit none save real(kind=kind_phys):: alf real(kind=kind_phys):: als real(kind=kind_phys):: alv real(kind=kind_phys):: cpd real(kind=kind_phys):: g real(kind=kind_phys):: rd real(kind=kind_phys):: rv real(kind=kind_phys),parameter:: t13 = 1.0/3.0 real(kind=kind_phys),parameter:: tmelt = 273.16 real(kind=kind_phys),parameter:: c1es = 610.78 real(kind=kind_phys),parameter:: c3les = 17.2693882 real(kind=kind_phys),parameter:: c3ies = 21.875 real(kind=kind_phys),parameter:: c4les = 35.86 real(kind=kind_phys),parameter:: c4ies = 7.66 real(kind=kind_phys),parameter:: rtwat = tmelt real(kind=kind_phys),parameter:: rtber = tmelt-5. real(kind=kind_phys),parameter:: rtice = tmelt-23. integer,parameter:: momtrans = 2 real(kind=kind_phys),parameter:: entrdd = 2.0e-4 real(kind=kind_phys),parameter:: cmfcmax = 1.0 real(kind=kind_phys),parameter:: cmfcmin = 1.e-10 real(kind=kind_phys),parameter:: cmfdeps = 0.30 real(kind=kind_phys),parameter:: zdnoprc = 2.0e4 real(kind=kind_phys),parameter:: cprcon = 1.4e-3 real(kind=kind_phys),parameter:: pgcoef = 0.7 real(kind=kind_phys):: rcpd real(kind=kind_phys):: c2es real(kind=kind_phys):: c5les real(kind=kind_phys):: c5ies real(kind=kind_phys):: r5alvcp real(kind=kind_phys):: r5alscp real(kind=kind_phys):: ralvdcp real(kind=kind_phys):: ralsdcp real(kind=kind_phys):: ralfdcp real(kind=kind_phys):: vtmpc1 real(kind=kind_phys):: zrg logical,parameter:: nonequil = .true. logical,parameter:: lmfpen = .true. logical,parameter:: lmfmid = .true. logical,parameter:: lmfscv = .true. logical,parameter:: lmfdd = .true. logical,parameter:: lmfdudv = .true. !================================================================================================================= end module cu_ntiedtke_common !================================================================================================================= module cu_ntiedtke use ccpp_kind_types,only: kind_phys use cu_ntiedtke_common implicit none private public:: cu_ntiedtke_run, & cu_ntiedtke_init, & cu_ntiedtke_finalize contains !================================================================================================================= !>\section arg_table_cu_ntiedtke_init !!\html\include cu_ntiedtke_init.html !! subroutine cu_ntiedtke_init(con_cp,con_rd,con_rv,con_xlv,con_xls,con_xlf,con_grav,errmsg,errflg) !================================================================================================================= !input arguments: real(kind=kind_phys),intent(in):: & con_cp, & con_rd, & con_rv, & con_xlv, & con_xls, & con_xlf, & con_grav !--- output arguments: character(len=*),intent(out):: errmsg integer,intent(out):: errflg !----------------------------------------------------------------------------------------------------------------- alf = con_xlf als = con_xls alv = con_xlv cpd = con_cp g = con_grav rd = con_rd rv = con_rv rcpd = 1.0/con_cp c2es = c1es*con_rd/con_rv c5les = c3les*(tmelt-c4les) c5ies = c3ies*(tmelt-c4ies) r5alvcp = c5les*con_xlv*rcpd r5alscp = c5ies*con_xls*rcpd ralvdcp = con_xlv*rcpd ralsdcp = con_xls*rcpd ralfdcp = con_xlf*rcpd vtmpc1 = con_rv/con_rd-1.0 zrg = 1.0/con_grav errmsg = 'cu_ntiedtke_init OK' errflg = 0 end subroutine cu_ntiedtke_init !================================================================================================================= !>\section arg_table_cu_ntiedtke_finalize !!\html\include cu_ntiedtke_finalize.html !! subroutine cu_ntiedtke_finalize(errmsg,errflg) !================================================================================================================= !--- output arguments: character(len=*),intent(out):: errmsg integer,intent(out):: errflg !----------------------------------------------------------------------------------------------------------------- errmsg = 'cu_ntiedtke_finalize OK' errflg = 0 end subroutine cu_ntiedtke_finalize !================================================================================================================= !>\section arg_table_cu_ntiedtke_run !!\html\include cu_ntiedtke_run.html !! ! level 1 subroutine 'cu_ntiedkte_run' subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqc,pqi,pqvf,ptf,poz,pzz,pomg, & & pap,paph,evap,hfx,zprecc,lndj,lq,km,km1,dt,dx,errmsg,errflg) !================================================================================================================= ! this is the interface between the model and the mass flux convection module ! m.tiedtke e.c.m.w.f. 1989 ! j.morcrette 1992 !-------------------------------------------- ! modifications ! C. zhang & Yuqing Wang 2011-2017 ! ! modified from IPRC IRAM - yuqing wang, university of hawaii (ICTP REGCM4.4). ! ! The current version is stable. There are many updates to the old Tiedtke scheme (cu_physics=6) ! update notes: ! the new Tiedtke scheme is similar to the Tiedtke scheme used in REGCM4 and ECMWF cy40r1. ! the major differences to the old Tiedtke (cu_physics=6) scheme are, ! (a) New trigger functions for deep and shallow convections (Jakob and Siebesma 2003; ! Bechtold et al. 2004, 2008, 2014). ! (b) Non-equilibrium situations are considered in the closure for deep convection ! (Bechtold et al. 2014). ! (c) New convection time scale for the deep convection closure (Bechtold et al. 2008). ! (d) New entrainment and detrainment rates for all convection types (Bechtold et al. 2008). ! (e) New formula for the conversion from cloud water/ice to rain/snow (Sundqvist 1978). ! (f) Different way to include cloud scale pressure gradients (Gregory et al. 1997; ! Wu and Yanai 1994) ! ! other reference: tiedtke (1989, mwr, 117, 1779-1800) ! IFS documentation - cy33r1, cy37r2, cy38r1, cy40r1 ! ! Note for climate simulation of Tropical Cyclones ! This version of Tiedtke scheme was tested with YSU PBL scheme, RRTMG radation ! schemes, and WSM6 microphysics schemes, at horizontal resolution around 20 km ! Set: momtrans = 2. ! pgcoef = 0.7 to 1.0 is good depends on the basin ! nonequil = .false. ! Note for the diurnal simulation of precipitaton ! When nonequil = .true., the CAPE is relaxed toward to a value from PBL ! It can improve the diurnal precipitation over land. !--- input arguments: integer,intent(in):: lq,km,km1 integer,intent(in),dimension(:):: lndj real(kind=kind_phys),intent(in):: dt real(kind=kind_phys),intent(in),dimension(:):: dx real(kind=kind_phys),intent(in),dimension(:):: evap,hfx real(kind=kind_phys),intent(in),dimension(:,:):: pqvf,ptf real(kind=kind_phys),intent(in),dimension(:,:):: poz,pomg,pap real(kind=kind_phys),intent(in),dimension(:,:):: pzz,paph !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(:):: zprecc real(kind=kind_phys),intent(inout),dimension(:,:):: pu,pv,pt,pqv,pqc,pqi !--- output arguments: character(len=*),intent(out):: errmsg integer,intent(out):: errflg !--- local variables and arrays: logical,dimension(lq):: locum integer:: i,j,k integer,dimension(lq):: icbot,ictop,ktype real(kind=kind_phys):: ztmst,fliq,fice,ztc,zalf,tt real(kind=kind_phys):: ztpp1,zew,zqs,zcor real(kind=kind_phys):: dxref real(kind=kind_phys),dimension(lq):: pqhfl,prsfc,pssfc,phhfl,zrain real(kind=kind_phys),dimension(lq):: scale_fac,scale_fac2 real(kind=kind_phys),dimension(lq,km):: pum1,pvm1,ztt,ptte,pqte,pvom,pvol,pverv,pgeo real(kind=kind_phys),dimension(lq,km):: zqq,pcte real(kind=kind_phys),dimension(lq,km):: ztp1,zqp1,ztu,zqu,zlu,zlude,zmfu,zmfd,zqsat real(kind=kind_phys),dimension(lq,km1):: pgeoh !----------------------------------------------------------------------------------------------------------------- ! ztmst=dt ! ! set scale-dependency factor when dx is < 15 km ! dxref = 15000. do j=1,lq if (dx(j).lt.dxref) then scale_fac(j) = (1.06133+log(dxref/dx(j)))**3 scale_fac2(j) = scale_fac(j)**0.5 else scale_fac(j) = 1.+1.33e-5*dx(j) scale_fac2(j) = 1. end if end do ! ! 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)=g*pzz(j,km1) end do ! ! convert model variables for mflux scheme ! do k=1,km do j=1,lq pcte(j,k)=0.0 pvom(j,k)=0.0 pvol(j,k)=0.0 ztp1(j,k)=pt(j,k) zqp1(j,k)=pqv(j,k)/(1.0+pqv(j,k)) pum1(j,k)=pu(j,k) pvm1(j,k)=pv(j,k) pverv(j,k)=pomg(j,k) pgeo(j,k)=g*poz(j,k) pgeoh(j,k)=g*pzz(j,k) tt=ztp1(j,k) zew = foeewm(tt) zqs = zew/pap(j,k) zqs = min(0.5,zqs) zcor = 1./(1.-vtmpc1*zqs) zqsat(j,k)=zqs*zcor pqte(j,k)=pqvf(j,k) zqq(j,k) =pqte(j,k) ptte(j,k)=ptf(j,k) ztt(j,k) =ptte(j,k) end do end do ! !----------------------------------------------------------------------- !* 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, & & ktype, icbot, ictop, ztu, zqu, & & zlu, zlude, zmfu, zmfd, zrain, & & pcte, phhfl, lndj, pgeoh, dx, & & scale_fac, scale_fac2) ! ! to include the cloud water and cloud ice detrained from convection ! do k=1,km do j=1,lq if(pcte(j,k).gt.0.) then fliq=foealfa(ztp1(j,k)) fice=1.0-fliq pqc(j,k)=pqc(j,k)+fliq*pcte(j,k)*ztmst pqi(j,k)=pqi(j,k)+fice*pcte(j,k)*ztmst endif end do end do ! do k=1,km do j=1,lq pt(j,k)= ztp1(j,k)+(ptte(j,k)-ztt(j,k))*ztmst zqp1(j,k)=zqp1(j,k)+(pqte(j,k)-zqq(j,k))*ztmst pqv(j,k)=zqp1(j,k)/(1.0-zqp1(j,k)) end do end do do j=1,lq zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst) end do if (lmfdudv) then do k=1,km do j=1,lq pu(j,k)=pu(j,k)+pvom(j,k)*ztmst pv(j,k)=pv(j,k)+pvol(j,k)*ztmst end do end do endif ! errmsg = 'cu_ntiedtke_run OK' errflg = 0 ! 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, & & ktype, kcbot, kctop, ptu, pqu, & & plu, plude, pmfu, pmfd, prain, & & pcte, phhfl, lndj, zgeoh, dx, & & scale_fac, scale_fac2) 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) !----------------------------------------------------------------- !--- input arguments: integer,intent(in):: klev,klon,klevp1,klevm1 integer,intent(in),dimension(klon):: lndj real(kind=kind_phys),intent(in):: ztmst real(kind=kind_phys),intent(in),dimension(klon):: dx real(kind=kind_phys),intent(in),dimension(klon):: pqhfl,phhfl real(kind=kind_phys),intent(in),dimension(klon):: scale_fac,scale_fac2 real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqen,puen,pven,pverv real(kind=kind_phys),intent(in),dimension(klon,klev):: pap,pgeo real(kind=kind_phys),intent(in),dimension(klon,klevp1):: paph,zgeoh !--- inout arguments: integer,intent(inout),dimension(klon):: ktype,kcbot,kctop logical,intent(inout),dimension(klon):: ldcum real(kind=kind_phys),intent(inout),dimension(klon):: pqsen real(kind=kind_phys),intent(inout),dimension(klon):: prsfc,pssfc,prain real(kind=kind_phys),intent(inout),dimension(klon,klev):: pcte,ptte,pqte,pvom,pvol real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptu,pqu,plu,plude,pmfu,pmfd !--- local variables and arrays: logical:: llo1 logical,dimension(klon):: loddraf,llo2 integer:: jl,jk,ik integer:: ikb,ikt,icum,itopm2 integer,dimension(klon):: kdpl,idtop,ictop0,ilwmin integer,dimension(klon,klev):: ilab 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 real(kind=kind_phys):: zduten,zdvten,ztdis,pgf_u,pgf_v real(kind=kind_phys):: zlon real(kind=kind_phys):: ztau,zerate,zderate,zmfa real(kind=kind_phys),dimension(klon):: zmfs real(kind=kind_phys),dimension(klon):: zsfl,zcape,zcape1,zcape2,ztauc,ztaubl,zheat real(kind=kind_phys),dimension(klon):: wup,zdqcv real(kind=kind_phys),dimension(klon):: wbase,zmfuub real(kind=kind_phys),dimension(klon):: upbl real(kind=kind_phys),dimension(klon):: zhcbase,zmfub,zmfub1,zdhpbl real(kind=kind_phys),dimension(klon):: zmfuvb,zsum12,zsum22 real(kind=kind_phys),dimension(klon):: zrfl real(kind=kind_phys),dimension(klev):: pmean real(kind=kind_phys),dimension(klon,klev):: pmfude_rate,pmfdde_rate real(kind=kind_phys),dimension(klon,klev):: zdpmel real(kind=kind_phys),dimension(klon,klev):: zmfuus,zmfdus,zuv2,ztenu,ztenv real(kind=kind_phys),dimension(klon,klev):: ztenh,zqenh,zqsenh,ztd,zqd real(kind=kind_phys),dimension(klon,klev):: zmfus,zmfds,zmfuq,zmfdq,zdmfup,zdmfdp,zmful real(kind=kind_phys),dimension(klon,klev):: zuu,zvu,zud,zvd,zlglac real(kind=kind_phys),dimension(klon,klevp1):: pmflxr,pmflxs !------------------------------------------- ! 1. specify constants and parameters !------------------------------------------- zcons=1./(g*ztmst) zcons2=3./(g*ztmst) !-------------------------------------------------------------- !* 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) ztauc(jl) = max(ztmst,ztauc(jl)) ztauc(jl) = max(360.,ztauc(jl)) ztauc(jl) = min(10800.,ztauc(jl)) ztau = ztauc(jl) * scale_fac(jl) if(nonequil) 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) = zmfub1(jl)/scale_fac2(jl) 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 !---------------------------------------------------------- !* 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 pgf_u = -pgcoef*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+& pmfu(jl,jk)*(puen(jl,jk)-puen(jl,jk-1))) pgf_v = -pgcoef*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+& pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1))) 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 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 ! ---------------------------------------------------------------- !--- input arguments: integer,intent(in):: klon,klev,klevp1,klevm1 real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqen,pqsen,puen,pven real(kind=kind_phys),intent(in),dimension(klon,klev):: pgeo,pverv real(kind=kind_phys),intent(in),dimension(klon,klev+1):: paph,pgeoh !--- output arguments: integer,intent(out),dimension(klon):: klwmin integer,intent(out),dimension(klon,klev):: klab real(kind=kind_phys),intent(out),dimension(klon,klev):: ptenh,pqenh,pqsenh real(kind=kind_phys),intent(out),dimension(klon,klev):: ptu,ptd,pqu,pqd,plu real(kind=kind_phys),intent(out),dimension(klon,klev):: puu,pud,pvu,pvd !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon,klev):: pmfu,pmfd,pmfus,pmfds,pmfuq,pmfdq real(kind=kind_phys),intent(inout),dimension(klon,klev):: pdmfup,pdmfdp,plude,pdpmel !--- local variables and arrays: logical,dimension(klon):: loflag integer:: jl,jk integer:: icall,ik real(kind=kind_phys):: zzs real(kind=kind_phys),dimension(klon):: zph,zwmax !------------------------------------------------------------ !* 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 subroutines !-------------------------------------------------------- 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 !------------------------------------------------------------------- !--- input arguments: integer,intent(in):: klon,klev,klevp1,klevm1 integer,intent(in),dimension(klon):: lndj real(kind=kind_phys),intent(in),dimension(klon):: qfx,hfx real(kind=kind_phys),intent(in),dimension(klon,klev):: pap,pgeo real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqen,pqsen real(kind=kind_phys),intent(in),dimension(klon,klev):: ptenh,pqenh,pqsenh real(kind=kind_phys),intent(in),dimension(klon,klevp1):: paph,pgeoh !--- output arguments: logical,intent(out),dimension(klon):: ldcum integer,intent(out),dimension(klon):: ktype integer,intent(out),dimension(klon):: cubot,cutop,kdpl integer,intent(out),dimension(klon,klev):: culab real(kind=kind_phys),intent(out),dimension(klon):: wbase real(kind=kind_phys),intent(out),dimension(klon,klev):: cutu,cuqu,culu !--- local variables and arrays: logical:: needreset logical,dimension(klon):: lldcum logical,dimension(klon):: loflag,deepflag,resetflag integer:: jl,jk,ik,icall,levels integer:: nk,is,ikb,ikt integer,dimension(klon):: kctop,kcbot integer,dimension(klon):: zcbase,itoppacel integer,dimension(klon,klev):: klab real(kind=kind_phys):: rho,part1,part2,root,conw,deltt,deltq 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 real(kind=kind_phys):: zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp real(kind=kind_phys):: zpdifftop, zpdiffbot real(kind=kind_phys),dimension(klon):: eta,dz,coef,zqold,zph real(kind=kind_phys),dimension(klon,klev):: dh,dhen,kup,vptu,vten real(kind=kind_phys),dimension(klon,klev):: ptu,pqu,plu real(kind=kind_phys),dimension(klon,klev):: zbuo,abuoy,plude !-------------------------------------------------------------- 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)))) 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.8/(pgeo(jl,jk)*zrg)+2.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,klev/2+1,-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 !--- input arguments: integer,intent(in):: klev,klon,klevp1,klevm1 integer,intent(in),dimension(klon):: lndj integer,intent(in),dimension(klon):: klwmin integer,intent(in),dimension(klon):: kdpl real(kind=kind_phys),intent(in):: ztmst real(kind=kind_phys),intent(in),dimension(klon):: wbase real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqen,pqsen,puen,pven,pqte,pverv real(kind=kind_phys),intent(in),dimension(klon,klev):: pap,pgeo real(kind=kind_phys),intent(in),dimension(klon,klevp1):: paph,pgeoh !--- inout arguments: logical,intent(inout),dimension(klon):: ldcum integer,intent(inout):: kcum integer,intent(inout),dimension(klon):: kcbot,kctop,kctop0 integer,intent(inout),dimension(klon,klev):: klab real(kind=kind_phys),intent(inout),dimension(klon):: phcbase real(kind=kind_phys),intent(inout),dimension(klon):: pmfub real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptenh,pqenh,pqsenh real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptu,pqu,plu,puu,pvu real(kind=kind_phys),intent(inout),dimension(klon,klev):: pmfu,pmfus,pmfuq,pmful,plude,pdmfup !--- output arguments: integer,intent(out),dimension(klon):: ktype real(kind=kind_phys),intent(out),dimension(klon):: wup real(kind=kind_phys),intent(out),dimension(klon,klev):: plglac,pmfude_rate !--- local variables and arrays: logical:: llo2,llo3 logical,dimension(klon):: loflag,llo1 integer:: jl,jk integer::ikb,icum,itopm2,ik,icall,is,jlm,jll integer,dimension(klon):: jlx real(kind=kind_phys):: 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 real(kind=kind_phys),dimension(klon):: eta,dz,zoentr,zdpmean real(kind=kind_phys),dimension(klon):: zph,zdmfen,zdmfde,zmfuu,zmfuv,zpbase,zqold,zluold,zprecip real(kind=kind_phys),dimension(klon,klev):: zlrain,zbuo,kup,zodetr,pdmfen !-------------------------------- !* 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) 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 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): ! *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): ! *prfl* precipitation rate kg/(m2*s) ! output parameters (real): ! *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 !--- input arguments: integer,intent(in):: klon logical,intent(in),dimension(klon):: ldcum integer,intent(in):: klev integer,intent(in),dimension(klon):: lndj integer,intent(in),dimension(klon):: kcbot,kctop real(kind=kind_phys),intent(in),dimension(klon):: pmfub real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqsen,pgeo,puen,pven real(kind=kind_phys),intent(in),dimension(klon,klev):: ptenh,pqenh real(kind=kind_phys),intent(in),dimension(klon,klev):: ptu,pqu,puu,pvu,plu real(kind=kind_phys),intent(in),dimension(klon,klev+1):: pgeoh,paph !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon):: prfl real(kind=kind_phys),intent(inout),dimension(klon,klev):: pud,pvd !--- output arguments: logical,intent(out),dimension(klon):: lddraf integer,intent(out),dimension(klon):: kdtop real(kind=kind_phys),intent(out),dimension(klon,klev):: ptd,pqd,pmfd,pmfds,pmfdq,pdmfdp !--- local variables and arrays: logical,dimension(klon):: llo2 integer:: jl,jk integer:: is,ik,icall,ike integer,dimension(klon):: ikhsmin real(kind=kind_phys):: zhsk,zttest,zqtest,zbuo,zmftop real(kind=kind_phys),dimension(klon):: zcond,zph,zhsmin real(kind=kind_phys),dimension(klon,klev):: ztenwb,zqenwb !---------------------------------------------------------------------- ! 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): ! *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 ! *paph* provisional pressure on half levels pa ! *pmfu* massflux updrafts kg/(m2*s) ! updated parameters (real): ! *prfl* precipitation rate kg/(m2*s) ! output parameters (real): ! *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 !--- input arguments: integer,intent(in)::klon logical,intent(in),dimension(klon):: lddraf integer,intent(in)::klev real(kind=kind_phys),intent(in),dimension(klon,klev):: ptenh,pqenh,puen,pven real(kind=kind_phys),intent(in),dimension(klon,klev):: pgeo,pmfu real(kind=kind_phys),intent(in),dimension(klon,klev+1):: pgeoh,paph !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon):: prfl real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptd,pqd,pud,pvd real(kind=kind_phys),intent(inout),dimension(klon,klev):: pmfd,pmfds,pmfdq,pdmfdp !--- output arguments: real(kind=kind_phys),intent(inout),dimension(klon,klev):: pmfdde_rate !--- local variables and arrays: logical:: llo1 logical,dimension(klon):: llo2 integer:: jl,jk integer:: is,ik,icall,ike integer,dimension(klon):: itopde real(kind=kind_phys):: zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp real(kind=kind_phys):: zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk real(kind=kind_phys),dimension(klon):: zdmfen,zdmfde,zcond,zoentr,zbuoy,zph !---------------------------------------------------------------------- ! 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): ! *ptsphy* 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 subroutines !-------------------------------------------------------- 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 !--- input arguments: integer,intent(in):: klon logical,intent(in),dimension(klon):: ldcum,lddraf integer,intent(in):: klev,ktopm2 integer,intent(in),dimension(klon):: kctop,kdtop real(kind=kind_phys),intent(in):: ztmst real(kind=kind_phys),intent(in),dimension(klon,klev):: pgeo real(kind=kind_phys),intent(in),dimension(klon,klev):: pten real(kind=kind_phys),intent(in),dimension(klon,klev):: pmfu,pmfus,pmfd,pmfds real(kind=kind_phys),intent(in),dimension(klon,klev):: pmfuq,pmfdq,pmful real(kind=kind_phys),intent(in),dimension(klon,klev):: plglac,plude,pdpmel real(kind=kind_phys),intent(in),dimension(klon,klev):: pdmfup,pdmfdp real(kind=kind_phys),intent(in),dimension(klon,klev):: pqen, ptenh,pqenh,pqsen real(kind=kind_phys),intent(in),dimension(klon,klev+1):: paph,pgeoh !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptent,ptenq,pcte !--- local variables and arrays: integer:: jk ,ik ,jl real(kind=kind_phys):: zalv ,zzp real(kind=kind_phys),dimension(klon,klev):: zdtdt,zdqdt,zdp !* 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 subroutines !-------------------------------------------------------- subroutine cududvn(klon,klev,ktopm2,ktype,kcbot,kctop,ldcum, & ztmst,paph,puen,pven,pmfu,pmfd,puu,pud,pvu,pvd,ptenu, & ptenv) implicit none !--- input arguments: integer,intent(in):: klon logical,intent(in),dimension(klon):: ldcum integer,intent(in):: klev,ktopm2 integer,intent(in),dimension(klon):: ktype,kcbot,kctop real(kind=kind_phys),intent(in):: ztmst real(kind=kind_phys),intent(in),dimension(klon,klev):: pmfu,pmfd,puen,pven real(kind=kind_phys),intent(in),dimension(klon,klev):: puu,pud,pvu,pvd real(kind=kind_phys),intent(in),dimension(klon,klev+1):: paph !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon,klev):: ptenu,ptenv !--- local variables and arrays: integer:: ik,ikb,jk,jl real(kind=kind_phys):: zzp,zdtdt real(kind=kind_Phys),dimension(klon,klev):: zdudt,zdvdt,zdp real(kind=kind_phys),dimension(klon,klev):: zuen,zven,zmfuu,zmfdu,zmfuv,zmfdv ! 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 subroutines !-------------------------------------------------------- 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): ! *psp* pressure pa ! updated parameters (real): ! *pt* temperature k ! *pq* specific humidity kg/kg ! externals ! --------- ! for condensation calculations. ! the tables are initialised in *suphec*. !---------------------------------------------------------------------- implicit none !--- input arguments: integer,intent(in):: klon logical,intent(in),dimension(klon):: ldflag integer,intent(in):: kcall,kk,klev real(kind=kind_phys),intent(in),dimension(klon):: psp !--- inout arguments: real(kind=kind_phys),intent(inout),dimension(klon,klev):: pt,pq !--- local variables and arrays: integer:: jl,jk integer:: isum 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 subroutines !-------------------------------------------------------- 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 ! ---------------------------------------------------------------- !--- input arguments: integer,intent(in):: klon logical,intent(in),dimension(klon):: ldcum integer,intent(in):: kk,klev,klevm1 real(kind=kind_phys),intent(in),dimension(klon,klev):: pten,pqen,pqsen,pgeo,pverv real(kind=kind_phys),intent(in),dimension(klon,klev):: puen,pven ! not used. real(kind=kind_phys),intent(in),dimension(klon,klev):: puu,pvu ! not used. real(kind=kind_phys),intent(in),dimension(klon,klev+1):: pgeoh !--- output arguments: integer,intent(out),dimension(klon):: ktype,kcbot integer,intent(out),dimension(klon,klev):: klab real(kind=kind_phys),intent(out),dimension(klon):: pmfub real(kind=kind_phys),intent(out),dimension(klon,klev):: plrain real(kind=kind_phys),intent(out),dimension(klon,klev):: ptu,pqu,plu real(kind=kind_phys),intent(out),dimension(klon,klev):: pmfu,pmfus,pmfuq,pmful real(kind=kind_phys),intent(out),dimension(klon,klev):: pdmfup !--- local variables and arrays: integer:: jl,klevp1 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 subroutines !--------------------------------------------------------- subroutine cuentrn(klon,klev,kk,kcbot,ldcum,ldwork, & pgeoh,pmfu,pdmfen,pdmfde) implicit none !--- input arguments: logical,intent(in):: ldwork integer,intent(in):: klon logical,intent(in),dimension(klon):: ldcum integer,intent(in):: klev,kk integer,intent(in),dimension(klon):: kcbot real(kind=kind_phys),intent(in),dimension(klon,klev):: pmfu real(kind=kind_phys),intent(in),dimension(klon,klev+1):: pgeoh !--- output arguments: real(kind=kind_phys),intent(out),dimension(klon):: pdmfen real(kind=kind_phys),intent(out),dimension(klon):: pdmfde !--- local variables and arrays: logical:: llo1 integer:: jl real(kind=kind_phys):: zdz ,zmf real(kind=kind_phys),dimension(klon):: zentr ! !* 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),intent(in):: 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),intent(in):: 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),intent(in):: 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),intent(in):: 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),intent(in):: tt foeldcpm = foealfa(tt)*ralvdcp+ & & (1.-foealfa(tt))*ralsdcp return end function foeldcpm !================================================================================================================= end module cu_ntiedtke !=================================================================================================================