diff --git a/physics/CONV/nTiedtke/cu_ntiedtke.F90 b/physics/CONV/nTiedtke/cu_ntiedtke.F90 index ada38c6f5..1de9de72b 100644 --- a/physics/CONV/nTiedtke/cu_ntiedtke.F90 +++ b/physics/CONV/nTiedtke/cu_ntiedtke.F90 @@ -14,27 +14,29 @@ module cu_ntiedtke ! this also requires redefining derived constants in the ! parameter section below use physcons, only:rd=>con_rd, rv=>con_rv, g=>con_g, & - & cpd=>con_cp, alv=>con_hvap, alf=>con_hfus + & cpd=>con_cp, alv=>con_hvap, alf=>con_hfus implicit none - real(kind=kind_phys),private :: rcpd,vtmpc1,tmelt,als,t13, & - c1es,c2es,c3les,c3ies,c4les,c4ies,c5les,c5ies,zrg - - real(kind=kind_phys),private :: rovcp,r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp,rtwat,rtber,rtice - real(kind=kind_phys),private :: entrdd,cmfcmax,cmfcmin,cmfdeps,zdnoprc,cprcon - integer,private :: momtrans,p650 - + real(kind=kind_phys),private :: rcpd,vtmpc1,als, & + c2es,c5les,c5ies,zrg + + real(kind=kind_phys),private :: rovcp,r5alvcp,r5alscp,ralvdcp,ralsdcp,ralfdcp + + real(kind=kind_phys),parameter:: t13 = 0.333333333 + 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. parameter( & - t13 = 0.333333333,& rcpd=1.0/cpd, & - tmelt=273.16, & zrg=1.0/g, & - c1es=610.78, & c2es=c1es*rd/rv, & - c3les=17.2693882, & - c3ies=21.875, & - c4les=35.86, & - c4ies=7.66, & als = alv+alf, & c5les=c3les*(tmelt-c4les), & c5ies=c3ies*(tmelt-c4ies), & @@ -43,62 +45,74 @@ module cu_ntiedtke ralvdcp=alv*rcpd, & ralsdcp=als*rcpd, & ralfdcp=alf*rcpd, & - rtwat=tmelt, & - rtber=tmelt-5., & - rtice=tmelt-23., & vtmpc1=rv/rd-1.0, & rovcp = rd*rcpd ) -! -! entrdd: average entrainment & detrainment rate for downdrafts + +! momtrans: momentum transport method ( 1 = IFS40r1 method; 2 = new method ) +! ------- + integer,parameter:: momtrans = 2 + +! entrdd: average turbulent entrainment & detrainment rate for downdrafts (Eq. 6.15 IFS Cy48r1) ! ------ -! - parameter(entrdd = 2.0e-4) -! + real(kind=kind_phys),parameter:: entrdd = 2.0e-4 + ! cmfcmax: maximum massflux value allowed for updrafts etc ! ------- -! - parameter(cmfcmax = 1.0) -! + real(kind=kind_phys),parameter:: cmfcmax = 1.0 + ! cmfcmin: minimum massflux value (for safety) ! ------- -! - parameter(cmfcmin = 1.e-10) -! + real(kind=kind_phys),parameter:: cmfcmin = 1.e-10 + ! cmfdeps: fractional massflux for downdrafts at lfs ! ------- -! - parameter(cmfdeps = 0.30) + real(kind=kind_phys),parameter:: cmfdeps = 0.30 ! zdnoprc: deep cloud is thicker than this height (Unit:Pa) -! - parameter(zdnoprc = 2.0e4) +! NRL changed from 2.0e4 to 1.5e4 as a result of NEPTUNE tuning experiments, +! see https://github.nrlmry.navy.mil/NEPTUNE/ccpp-physics/pull/28 ! ------- -! + real(kind=kind_phys),parameter:: zdnoprc = 1.5e4 + ! cprcon: coefficient from cloud water to rain water -! - parameter(cprcon = 1.4e-3) ! ------- -! -! momtrans: momentum transport method -! ( 1 = IFS40r1 method; 2 = new method ) -! - parameter(momtrans = 2 ) + real(kind=kind_phys),parameter:: cprcon = 1.4e-3 + +! pgcoef: 0.7 to 1.0 is good depends on the basin ! ------- -! - logical :: isequil + real(kind=kind_phys),parameter:: pgcoef = 0.7 + +! entorg: organized updraft entrainment scaling factor (Eq. 6.7 IFS Cy48r1) +! NRL changed from 1.75e-3 to 2.10e-3 as a result of NEPTUNE tuning experiments, +! see https://github.nrlmry.navy.mil/NEPTUNE/ccpp-physics/pull/28 +! ------- + real(kind=kind_phys),parameter:: entorg = 2.1e-3 ! exp 2.4, 2.1, and 1.4, orig. 1.75 + +! detturb: turbulent detrainment scaling factor (Eq. 6.8 IFS Cy48r1) +! ------- + real(kind=kind_phys),parameter:: detturb = 0.75e-4 + ! isequil: representing equilibrium and nonequilibrium convection ! ( .false. [default]; .true. [experimental]. Ref. Bechtold et al. 2014 JAS ) -! - parameter(isequil = .false. ) +! Note for the diurnal simulation of precipitaton +! When isequil = .true., the CAPE is relaxed toward to a value from PBL +! It can improve the diurnal precipitation over land. +! ------- + logical,parameter:: isequil = .false. ! !-------------------- ! switches for deep, mid, shallow convections, downdraft, and momemtum transport ! ------------------ - logical :: lmfpen,lmfmid,lmfscv,lmfdd,lmfdudv - parameter(lmfpen=.true.,lmfmid=.true.,lmfscv=.true.,lmfdd=.true.,lmfdudv=.true.) -!-------------------- + logical,parameter:: lmfpen = .true. + logical,parameter:: lmfmid = .true. + logical,parameter:: lmfscv = .true. + logical,parameter:: lmfdd = .true. + logical,parameter:: lmfdudv = .true. + + +!================================================================================================================= !#################### end of variables definition########################## -!----------------------------------------------------------------------- +!================================================================================================================= ! contains !> \brief Brief description of the subroutine @@ -112,16 +126,16 @@ subroutine cu_ntiedtke_init(imfshalcnv, imfshalcnv_ntiedtke, imfdeepcnv, & implicit none integer, intent(in) :: imfshalcnv, imfshalcnv_ntiedtke - integer, intent(in) :: imfdeepcnv, imfdeepcnv_ntiedtke + integer, intent(in) :: imfdeepcnv, imfdeepcnv_ntiedtke integer, intent(in) :: mpirank integer, intent(in) :: mpiroot character(len=*), intent( out) :: errmsg integer, intent( out) :: errflg - + ! initialize ccpp error handling variables errmsg = '' errflg = 0 - + ! DH* temporary if (mpirank==mpiroot) then write(0,*) ' -----------------------------------------------------------------------------------------------------------------------------' @@ -144,66 +158,92 @@ subroutine cu_ntiedtke_init(imfshalcnv, imfshalcnv_ntiedtke, imfdeepcnv, & errflg = 1 return end if - + end subroutine cu_ntiedtke_init -! Tiedtke cumulus scheme from WRF with small modifications -! This scheme includes both deep and shallow convections !=================== ! !> \section arg_table_cu_ntiedtke_run Argument Table !! \htmlinclude cu_ntiedtke_run.html !! -!----------------------------------------------------------------------- -! level 1 subroutine 'tiecnvn' -!----------------------------------------------------------------- +!================================================================================================================= +! level 1 subroutine 'cu_ntiedkte_run' subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & evap,hfx,zprecc,lmask,lq,km,dt,dx,kbot,ktop,kcnv, & ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg) -!----------------------------------------------------------------- -! this is the interface between the model and the mass -! flux convection module -!----------------------------------------------------------------- +!================================================================================================================= +! 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 +! implicit none -! in&out variables - integer, intent(in) :: lq, km, ktrac - real(kind=kind_phys), intent(in ) :: dt - integer, dimension( : ), intent(in) :: lmask - real(kind=kind_phys), dimension( : ), intent(in ) :: evap, hfx, dx - real(kind=kind_phys), dimension( :, : ), intent(inout) :: pu, pv, pt, pqv - real(kind=kind_phys), dimension( :, :), intent(in ) :: tdi, qvdi, poz, prsl, pomg - real(kind=kind_phys), dimension( :, :), intent(in ), optional :: pqvf, ptf - real(kind=kind_phys), dimension( :, : ), intent(in ) :: pzz, prsi - real(kind=kind_phys), dimension( :, :, : ), intent(inout ) :: clw - - integer, dimension( : ), intent(out) :: kbot, ktop, kcnv - real(kind=kind_phys), dimension( : ), intent(out) :: zprecc - real(kind=kind_phys), dimension (:, :), intent(out), optional :: ud_mf - real(kind=kind_phys), dimension (:, :), intent(out) :: dd_mf, dt_mf, cnvw, cnvc - +!--- input arguments: + integer, intent(in) :: lq, km, ktrac + integer, intent(in), dimension(:) :: lmask + + real(kind=kind_phys), intent(in) :: dt + real(kind=kind_phys), dimension(:), intent(in) :: evap, hfx, dx + real(kind=kind_phys), dimension(:,:), intent(in) :: tdi, qvdi, poz, prsl, pomg + real(kind=kind_phys), dimension(:,:), intent(in), optional :: pqvf, ptf + real(kind=kind_phys), dimension(:,:),intent(in) :: pzz, prsi + +!--- inout arguments: + real(kind=kind_phys), dimension(:,:,:), intent(inout ) :: clw + real(kind=kind_phys), dimension(:,:), intent(inout) :: pu, pv, pt, pqv + +!--- output arguments: + real(kind=kind_phys), dimension(:), intent(out) :: zprecc + integer, dimension(:), intent(out) :: kbot, ktop, kcnv + real(kind=kind_phys), dimension (:,:), intent(out), optional :: ud_mf + real(kind=kind_phys), dimension (:,:), intent(out) :: dd_mf, dt_mf, cnvw, cnvc + ! error messages - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + ! local variables - real(kind=kind_phys) pum1(lq,km), pvm1(lq,km), ztt(lq,km), & - & ptte(lq,km), pqte(lq,km), pvom(lq,km), pvol(lq,km), & - & pverv(lq,km), pgeo(lq,km), pap(lq,km), paph(lq,km+1) - real(kind=kind_phys) pqhfl(lq), zqq(lq,km), & - & prsfc(lq), pssfc(lq), pcte(lq,km), & - & phhfl(lq), pgeoh(lq,km+1) - real(kind=kind_phys) ztp1(lq,km), zqp1(lq,km), ztu(lq,km), zqu(lq,km),& - & zlu(lq,km), zlude(lq,km), zmfu(lq,km), zmfd(lq,km), zmfude_rate(lq,km),& - & zqsat(lq,km), zrain(lq) - real(kind=kind_phys),allocatable :: pcen(:,:,:),ptenc(:,:,:) - - integer icbot(lq), ictop(lq), ktype(lq), lndj(lq) - logical locum(lq) -! - real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt - integer i,j,k,k1,n,km1,ktracer - real(kind=kind_phys) ztpp1 - real(kind=kind_phys) zew,zqs,zcor + real(kind=kind_phys),allocatable :: pcen(:,:,:),ptenc(:,:,:) + integer,dimension(lq):: lndj + logical,dimension(lq):: locum + integer:: i,j,k + integer:: k1,n,km1,ktracer + 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,zmfude_rate,pap + real(kind=kind_phys),dimension(lq,km+1):: pgeoh,paph + +!----------------------------------------------------------------------------------------------------------------- ! ! Initialize CCPP error handling variables errmsg = '' @@ -212,6 +252,19 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, km1 = km + 1 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 @@ -222,7 +275,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, pqhfl(j)=evap(j) phhfl(j)=hfx(j) pgeoh(j,km1)=pzz(j,1) - paph(j,km1)=prsi(j,1) + paph(j,km1)=prsi(j,1) if(lmask(j).eq.1) then lndj(j)=1 else @@ -248,12 +301,12 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, pap(j,k1)=prsl(j,k) paph(j,k1)=prsi(j,k+1) tt=ztp1(j,k1) - zew = foeewm(tt) - zqs = zew/pap(j,k1) - zqs = min(0.5,zqs) - zcor = 1./(1.-vtmpc1*zqs) - zqsat(j,k1)=zqs*zcor - pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst + zew = foeewm(tt) + zqs = zew/pap(j,k1) + zqs = min(0.5,zqs) + zcor = 1./(1.-vtmpc1*zqs) + zqsat(j,k1)=zqs*zcor + pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst zqq(j,k1) =pqte(j,k1) ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst ztt(j,k1) =ptte(j,k1) @@ -291,13 +344,13 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, end do end do end if - + ! print *, "pgeo=",pgeo(1,:) ! print *, "pgeoh=",pgeoh(1,:) ! print *, "pap=",pap(1,:) ! print *, "paph=",paph(1,:) ! print *, "ztp1=",ztp1(1,:) -! print *, "zqp1=",zqp1(1,:) +! print *, "zqp1=",zqp1(1,:) ! print *, "pum1=",pum1(1,:) ! print *, "pvm1=",pvm1(1,:) ! print *, "pverv=",pverv(1,:) @@ -309,14 +362,15 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, !* 2. call 'cumastrn'(master-routine for cumulus parameterization) ! call cumastrn & - & (lq, km, km1, km-1, ztp1, & - & zqp1, pum1, pvm1, pverv, zqsat,& - & pqhfl, ztmst, pap, paph, pgeo, & - & ptte, pqte, pvom, pvol, prsfc,& - & pssfc, locum, ktracer, pcen, ptenc,& - & ktype, icbot, ictop, ztu, zqu, & - & zlu, zlude, zmfu, zmfd, zrain,& - & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx) + & (lq, km, km1, km-1, ztp1, & + & zqp1, pum1, pvm1, pverv, zqsat, & + & pqhfl, ztmst, pap, paph, pgeo, & + & ptte, pqte, pvom, pvol, prsfc, & + & pssfc, locum, ktracer, pcen, ptenc, & + & ktype, icbot, ictop, ztu, zqu, & + & zlu, zlude, zmfu, zmfd, zrain, & + & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx, & + & scale_fac, scale_fac2) ! ! to include the cloud water and cloud ice detrained from convection ! @@ -336,7 +390,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, k1 = km-k+1 do j=1,lq pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst - pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst + pqv(j,k) = zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst ud_mf(j,k)= zmfu(j,k1)*ztmst dd_mf(j,k)= -zmfd(j,k1)*ztmst dt_mf(j,k)= zmfude_rate(j,k1)*ztmst @@ -351,8 +405,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, zprecc(j)=amax1(0.0,(prsfc(j)+pssfc(j))*ztmst*0.001) kbot(j) = km-icbot(j)+1 ktop(j) = km-ictop(j)+1 + ! deep convection flag if(ktype(j).eq.1 .or. ktype(j).eq.3) then - kcnv(j)=1 + kcnv(j)=1 else kcnv(j)=0 end if @@ -367,21 +422,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, end do end do endif - ! -! Currently, vertical mixing of tracers are turned off -! if(ktrac > 2) then -! do n=1,ktrac-2 -! do k=1,km -! k1=km-k+1 -! do j=1,lq -! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst -! end do -! end do -! end do -! end if - deallocate(pcen) - deallocate(ptenc) + errmsg = 'cu_ntiedtke_run OK' + errflg = 0 ! return end subroutine cu_ntiedtke_run @@ -395,14 +438,15 @@ end subroutine cu_ntiedtke_run ! subroutine cumastrn !*********************************************************** subroutine cumastrn & - & (klon, klev, klevp1, klevm1, pten, & - & pqen, puen, pven, pverv, pqsen,& - & pqhfl, ztmst, pap, paph, pgeo, & - & ptte, pqte, pvom, pvol, prsfc,& - & pssfc, ldcum, ktrac, pcen, ptenc,& - & ktype, kcbot, kctop, ptu, pqu,& - & plu, plude, pmfu, pmfd, prain,& - & pcte, phhfl, lndj, zgeoh, pmfude_rate, dx) + & (klon, klev, klevp1, klevm1, pten, & + & pqen, puen, pven, pverv, pqsen, & + & pqhfl, ztmst, pap, paph, pgeo, & + & ptte, pqte, pvom, pvol, prsfc, & + & pssfc, ldcum, ktrac, pcen, ptenc, & + & ktype, kcbot, kctop, ptu, pqu, & + & plu, plude, pmfu, pmfd, prain, & + & pcte, phhfl, lndj, zgeoh, pmfude_rate, dx, & + & scale_fac, scale_fac2) implicit none ! !***cumastrn* master routine for cumulus massflux-scheme @@ -462,92 +506,81 @@ subroutine cumastrn & ! ---------- ! paper on massflux scheme (tiedtke,1989) !----------------------------------------------------------------- - integer klev,klon,ktrac,klevp1,klevm1 - real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),& - & puen(klon,klev), pven(klon,klev),& - & ptte(klon,klev), pqte(klon,klev),& - & pvom(klon,klev), pvol(klon,klev),& - & pqsen(klon,klev), pgeo(klon,klev),& - & pap(klon,klev), paph(klon,klevp1),& - & pverv(klon,klev), pqhfl(klon),& - & phhfl(klon) - real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& - & plu(klon,klev), plude(klon,klev),& - & pmfu(klon,klev), pmfd(klon,klev),& - & prain(klon),& - & prsfc(klon), pssfc(klon) - real(kind=kind_phys) ztenh(klon,klev), zqenh(klon,klev),& - & zgeoh(klon,klevp1), zqsenh(klon,klev),& - & ztd(klon,klev), zqd(klon,klev),& - & zmfus(klon,klev), zmfds(klon,klev),& - & zmfuq(klon,klev), zmfdq(klon,klev),& - & zdmfup(klon,klev), zdmfdp(klon,klev),& - & zmful(klon,klev), zrfl(klon),& - & zuu(klon,klev), zvu(klon,klev),& - & zud(klon,klev), zvd(klon,klev),& - & zlglac(klon,klev) - real(kind=kind_phys) pmflxr(klon,klevp1), pmflxs(klon,klevp1) - real(kind=kind_phys) zhcbase(klon),& - & zmfub(klon), zmfub1(klon),& - & zdhpbl(klon) - real(kind=kind_phys) zsfl(klon), zdpmel(klon,klev),& - & pcte(klon,klev), zcape(klon),& - & zcape1(klon), zcape2(klon),& - & ztauc(klon), ztaubl(klon),& - & zheat(klon) - real(kind=kind_phys) pcen(klon,klev,ktrac), ptenc(klon,klev,ktrac) - real(kind=kind_phys) wup(klon), zdqcv(klon) - real(kind=kind_phys) wbase(klon), zmfuub(klon) - real(kind=kind_phys) upbl(klon) - real(kind=kind_phys) dx(klon) - real(kind=kind_phys) pmfude_rate(klon,klev), pmfdde_rate(klon,klev) - real(kind=kind_phys) zmfuus(klon,klev), zmfdus(klon,klev) - real(kind=kind_phys) zmfudr(klon,klev), zmfddr(klon,klev) - real(kind=kind_phys) zuv2(klon,klev),ztenu(klon,klev),ztenv(klon,klev) - real(kind=kind_phys) zmfuvb(klon),zsum12(klon),zsum22(klon) - integer ilab(klon,klev), idtop(klon),& - & ictop0(klon), ilwmin(klon) - integer kdpl(klon) - integer kcbot(klon), kctop(klon),& - & ktype(klon), lndj(klon) - logical ldcum(klon), lldcum(klon) - logical loddraf(klon), llddraf3(klon), llo1, llo2(klon) - -! local varaiables - real(kind=kind_phys) zcons,zcons2,zqumqe,zdqmin,zdh,zmfmax - real(kind=kind_phys) zalfaw,zalv,zqalv,zc5ldcp,zc4les,zhsat,zgam,zzz,zhhat - real(kind=kind_phys) zpbmpt,zro,zdz,zdp,zeps,zfac,wspeed - integer jl,jk,ik - integer ikb,ikt,icum,itopm2 - real(kind=kind_phys) ztmst,ztau,zerate,zderate,zmfa - real(kind=kind_phys) zmfs(klon),pmean(klev),zlon - real(kind=kind_phys) zduten,zdvten,ztdis,pgf_u,pgf_v + +!--- input arguments: + integer,intent(in):: klev,klon,klevp1,klevm1 + integer,intent(in):: ktrac + 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 + logical,dimension(klon):: lldcum,llddraf3 + + 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,ktrac):: pcen, ptenc + 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):: zmfudr,zmfddr + 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) - zlon = real(klon) - do jk = klev , 1 , -1 - pmean(jk) = sum(pap(:,jk))/zlon - end do - p650 = klev-2 - do jk = klev , 3 , -1 - if ( pmean(jk)/pmean(klev)*1.013250e5 > 650.e2 ) p650 = jk - end do - !-------------------------------------------------------------- !* 2. initialize values at vertical grid points in 'cuini' !-------------------------------------------------------------- call cuinin & - & (klon, klev, klevp1, klevm1, pten, & - & pqen, pqsen, puen, pven, pverv,& - & pgeo, paph, zgeoh, ztenh, zqenh,& - & zqsenh, ilwmin, ptu, pqu, ztd, & - & zqd, zuu, zvu, zud, zvd, & - & pmfu, pmfd, zmfus, zmfds, zmfuq,& - & zmfdq, zdmfup, zdmfdp, zdpmel, plu, & + & (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) !---------------------------------- @@ -557,11 +590,12 @@ subroutine cumastrn & ! 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) + & ( 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 ! ------------------------------------------ @@ -571,41 +605,74 @@ subroutine cumastrn & idtop(jl)=0 end do + !----------------------------------------------- + ! Calculate moist static energy and kinetic + ! energy within the subcloud layer for the + ! environment + !----------------------------------------------- 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)) + ! sum subcloud layer moist static energy + 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) + 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 + !----------------------------------------------- + ! Calculate first guess cloud base mass flux + !----------------------------------------------- do jl=1,klon + if(ldcum(jl)) then - ikb=kcbot(jl) + + ikb = kcbot(jl) zmfmax = (paph(jl,ikb)-paph(jl,ikb-1))*zcons2 + !----------------------------------------------- + ! Deep convection. + ! Initial updraft mass flux is 10% of its maximum + ! value, which is determined by the layer thickness + ! and time step. + !----------------------------------------------- if(ktype(jl) == 1) then + zmfub(jl)= 0.1*zmfmax + !----------------------------------------------- + ! Shallow convection. + ! Initial updraft mass flux is determined by + ! a balance of moist static energy in the + ! boundary layer. + !----------------------------------------------- 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) + !zdh = g*max(zdh,1.e5*zdqmin) + zdh = g*max(zdh,0.75*cpd) ! limiter updated to be consistent with IFS documentation + if ( zdhpbl(jl) > 0. ) then - zmfub(jl) = zdhpbl(jl)/zdh + zmfub(jl) = zdhpbl(jl) / zdh zmfub(jl) = min(zmfub(jl),zmfmax) else zmfub(jl) = 0.1*zmfmax ldcum(jl) = .false. end if - end if + + end if + else zmfub(jl) = 0. end if + end do !------------------------------------------------------ !* 4.0 determine cloud ascent for entraining plume @@ -613,15 +680,16 @@ subroutine cumastrn & !* (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 ) + & (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) @@ -660,24 +728,24 @@ subroutine cumastrn & if(lmfdd) then !* (a) determine lfs in 'cudlfsn' !-------------------------------------- - call 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, & + & 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 ) + 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 ! @@ -685,7 +753,7 @@ subroutine cumastrn & !* 6.0 closure and clean work ! ------ !-- 6.1 recalculate cloud base massflux from a cape closure -! for deep convection (ktype=1) +! for deep convection (ktype=1) ! do jl=1,klon if(ldcum(jl) .and. ktype(jl) .eq. 1) then @@ -696,19 +764,19 @@ subroutine cumastrn & 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 + 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 if end do -! + do jk = 1 , klev do jl = 1 , klon llo1 = ldcum(jl) .and. ktype(jl) .eq. 1 @@ -727,7 +795,7 @@ subroutine cumastrn & 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 + ((1.+vtmpc1*pqen(jl,jk))*ptte(jl,jk)+vtmpc1*pten(jl,jk)*pqte(jl,jk))*zdp end if end if end do @@ -737,10 +805,10 @@ subroutine cumastrn & if(ldcum(jl).and.ktype(jl).eq.1) then ikb = kcbot(jl) ikt = kctop(jl) - ztau = ztauc(jl) * (1.+1.33e-5*dx(jl)) - ztau = max(ztmst,ztau) - ztau = max(720.,ztau) - ztau = min(10800.,ztau) + 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(isequil) then zcape2(jl)= max(0.,zcape2(jl)) zcape(jl) = max(0.,min(zcape1(jl)-zcape2(jl),5000.)) @@ -754,44 +822,52 @@ subroutine cumastrn & 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) -!-------------------------------------------------------- + ! + !* 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) + 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) + ! 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) + zdh = g*max(zdh,0.75*cpd) ! limiter updated to be consistent with IFS documentation + if(zdhpbl(jl).gt.0.)then - zmfub1(jl)=zdhpbl(jl)/zdh + 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 -!--------------------------------------------------------- + !* 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 -!--------------------------------------------------------- + !* 6.4 scaling the downdraft mass flux + !--------------------------------------------------------- do jk=1,klev do jl=1,klon if( ldcum(jl) ) then @@ -805,8 +881,8 @@ subroutine cumastrn & end do end do -!* 6.5 scaling the updraft mass flux -! -------------------------------------------------------- + !* 6.5 scaling the updraft mass flux + ! -------------------------------------------------------- do jl = 1,klon if ( ldcum(jl) ) zmfs(jl) = zmfub1(jl)/max(cmfcmin,zmfub(jl)) end do @@ -839,8 +915,8 @@ subroutine cumastrn & end do end do -!* 6.6 if ktype = 2, kcbot=kctop is not allowed -! --------------------------------------------------- + !* 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 @@ -860,8 +936,8 @@ subroutine cumastrn & end do end if -!* 6.7 set downdraft mass fluxes to zero above cloud top -!---------------------------------------------------- + !* 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 @@ -882,21 +958,19 @@ subroutine cumastrn & end if end do end do - - itopm2 = 2 !---------------------------------------------------------- !* 7.0 determine final convective fluxes in 'cuflx' !---------------------------------------------------------- - call cuflxn & - & ( klon, klev, ztmst & - & , pten, pqen, pqsen, ztenh, zqenh & - & , paph, pap, zgeoh, lndj, ldcum & - & , kcbot, kctop, idtop, itopm2 & - & , ktype, loddraf & - & , pmfu, pmfd, zmfus, zmfds & - & , zmfuq, zmfdq, zmful, plude & - & , zdmfup, zdmfdp, zdpmel, zlglac & - & , prain, pmfdde_rate, pmflxr, pmflxs ) + 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 @@ -987,9 +1061,9 @@ subroutine cumastrn & !---------------------------------------------------------------- !* 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, & + 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 @@ -1021,15 +1095,10 @@ subroutine cumastrn & zvu(jl,jk) = (zvu(jl,ik)*pmfu(jl,ik) + & zerate*pven(jl,jk)-zderate*zvu(jl,ik))*zmfa else - if(ktype(jl) == 1 .or. ktype(jl) == 3) then - pgf_u = -0.7*0.5*(pmfu(jl,ik)*(puen(jl,ik)-puen(jl,jk))+& + 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 = -0.7*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+& + pgf_v = -pgcoef*0.5*(pmfu(jl,ik)*(pven(jl,ik)-pven(jl,jk))+& pmfu(jl,jk)*(pven(jl,jk)-pven(jl,jk-1))) - else - pgf_u = 0. - pgf_v = 0. - end if zerate = pmfu(jl,jk) - pmfu(jl,ik) + pmfude_rate(jl,jk) zderate = pmfude_rate(jl,jk) zmfa = 1./max(cmfcmin,pmfu(jl,jk)) @@ -1215,13 +1284,13 @@ end subroutine cumastrn !********************************************** ! 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,& + & (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 @@ -1240,30 +1309,33 @@ subroutine cuinin & ! --------- ! *cuadjtq* to specify qs at half levels ! ---------------------------------------------------------------- - integer klon,klev,klevp1,klevm1 - real(kind=kind_phys) pten(klon,klev), pqen(klon,klev),& - & puen(klon,klev), pven(klon,klev),& - & pqsen(klon,klev), pverv(klon,klev),& - & pgeo(klon,klev), pgeoh(klon,klevp1),& - & paph(klon,klevp1), ptenh(klon,klev),& - & pqenh(klon,klev), pqsenh(klon,klev) - real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& - & ptd(klon,klev), pqd(klon,klev),& - & puu(klon,klev), pud(klon,klev),& - & pvu(klon,klev), pvd(klon,klev),& - & pmfu(klon,klev), pmfd(klon,klev),& - & pmfus(klon,klev), pmfds(klon,klev),& - & pmfuq(klon,klev), pmfdq(klon,klev),& - & pdmfup(klon,klev), pdmfdp(klon,klev),& - & plu(klon,klev), plude(klon,klev) - real(kind=kind_phys) zwmax(klon), zph(klon), & - & pdpmel(klon,klev) - integer klab(klon,klev), klwmin(klon) - logical loflag(klon) -! local variables - integer jl,jk - integer icall,ik - real(kind=kind_phys) zzs + +!--- 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 @@ -1339,14 +1411,15 @@ subroutine cuinin & end subroutine cuinin !--------------------------------------------------------- -! level 3 souroutines +! 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 ) + 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. ! -------- @@ -1372,7 +1445,7 @@ subroutine cutypen & ! climate, mon.wea.rev. ! 131, 2765-2778 ! and -! ifs documentation - cy36r1,cy38r1 +! ifs documentation - cy36r1,cy38r1 !***input variables: ! ptenh [ztenh] - environment temperature on half levels. (cuini) ! pqenh [zqenh] - env. specific humidity on half levels. (cuini) @@ -1390,52 +1463,53 @@ subroutine cutypen & !------------------------------------------------------------------- implicit none !------------------------------------------------------------------- - integer klon, klev, klevp1, klevm1 - real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev),& - & pqsen(klon,klev), pqsenh(klon,klev),& - & pgeoh(klon,klevp1), paph(klon,klevp1),& - & pap(klon,klev), pqen(klon,klev) - real(kind=kind_phys) pten(klon,klev) - real(kind=kind_phys) ptu(klon,klev),pqu(klon,klev),plu(klon,klev) - real(kind=kind_phys) pgeo(klon,klev) - integer klab(klon,klev) - integer kctop(klon),kcbot(klon) - - real(kind=kind_phys) qfx(klon),hfx(klon) - real(kind=kind_phys) zph(klon) - integer lndj(klon) - logical loflag(klon), deepflag(klon), resetflag(klon) - -! output variables - real(kind=kind_phys) cutu(klon,klev), cuqu(klon,klev), culu(klon,klev) - integer culab(klon,klev) - real(kind=kind_phys) wbase(klon) - integer ktype(klon),cubot(klon),cutop(klon),kdpl(klon) - logical ldcum(klon) -! local variables - real(kind=kind_phys) zqold(klon) - real(kind=kind_phys) rho, part1, part2, root, conw, deltt, deltq - real(kind=kind_phys) eta(klon),dz(klon),coef(klon) - real(kind=kind_phys) dhen(klon,klev), dh(klon,klev) - real(kind=kind_phys) plude(klon,klev) - real(kind=kind_phys) kup(klon,klev) - real(kind=kind_phys) vptu(klon,klev),vten(klon,klev) - real(kind=kind_phys) zbuo(klon,klev),abuoy(klon,klev) - - real(kind=kind_phys) zz,zdken,zdq - real(kind=kind_phys) fscale,crirh1,pp - real(kind=kind_phys) atop1,atop2,abot - real(kind=kind_phys) tmix,zmix,qmix,pmix - real(kind=kind_phys) zlglac,dp - integer nk,is,ikb,ikt - - real(kind=kind_phys) zqsu,zcor,zdp,zesdp,zalfaw,zfacw,zfaci,zfac,zdsdp,zdqsdt,zdtdp - real(kind=kind_phys) zpdifftop, zpdiffbot - integer zcbase(klon), itoppacel(klon) - integer jl,jk,ik,icall,levels - logical needreset, lldcum(klon) -!-------------------------------------------------------------- +!--- input arguments: + integer,intent(in):: klon,klev,klevp1,klevm1 + integer,intent(in),dimension(klon):: lndj + + real(kind=kind_phys),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 @@ -1445,24 +1519,23 @@ subroutine cutypen & ldcum(jl)=.false. end do -!----------------------------------------------------------- -! let's do test,and check the shallow convection first -! the first level is klev -! define deltat and deltaq -!----------------------------------------------------------- + !----------------------------------------------------------- + ! let's do test,and check the shallow convection first + ! the first level is klev + !----------------------------------------------------------- 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 + 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 @@ -1472,76 +1545,100 @@ subroutine cutypen & loflag(jl) = .true. end do -! check the levels from lowest level to second top level + ! check the levels from lowest level to second top level do jk=klevm1,2,-1 -! define the variables at the first level + ! 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)))) + + rho = pap(jl,klev) / & + (rd*(pten(jl,klev)*(1.+vtmpc1*pqen(jl,klev)))) hfx(jl) = hfx(jl)*rho*cpd qfx(jl) = qfx(jl)*rho part1 = 1.5*0.4*pgeo(jl,klev)/ & - & (rho*pten(jl,klev)) + (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 + 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) + 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 - + + end if ! end if(jk .eq. klevm1) + 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 + ! the next levels, we use the variables at the first level as initial values do jl=1,klon if(loflag(jl)) then - eta(jl) = 0.55/(pgeo(jl,jk)*zrg)+1.e-4 + !---------------------------------------- + ! Parcel entrainment rate for shallow convection. + ! Used to determine whether or not shallow + ! convection is triggered. Final entrainment + ! rate is calculated later. + !---------------------------------------- + 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)) + + (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)) + + (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 + ! 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 + ! 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)) @@ -1552,8 +1649,8 @@ subroutine cutypen & 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 + ! 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) @@ -1570,9 +1667,10 @@ subroutine cutypen & 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) + ! 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 @@ -1583,6 +1681,7 @@ subroutine cutypen & klab(jl,jk) = 2 kcbot(jl) = jk end if + end if if(kup(jl,jk) .lt. 0.)then @@ -1593,7 +1692,7 @@ subroutine cutypen & else lldcum(jl) = .false. end if - else + else if(plu(jl,jk) .gt. 0.)then klab(jl,jk)=2 else @@ -1606,9 +1705,13 @@ subroutine cutypen & 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 shallow cumulus is found, define some properties if(lldcum(jl)) then ktype(jl) = 2 ldcum(jl) = .true. @@ -1623,6 +1726,7 @@ subroutine cutypen & ldcum(jl) = .false. wbase(jl) = 0. end if + end do do jk=klev,1,-1 @@ -1636,16 +1740,17 @@ subroutine cutypen & 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 + + !----------------------------------------------------------- + ! 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 ! give parcel a small temperature perturbation at surface (Eq. 6.21 IFS Cy48r1) + deltq = 1.0e-4 ! give parcel a small humidity perturbation at surface (Eq. 6.21 IFS Cy48r1) + do jl=1,klon deepflag(jl) = .false. end do @@ -1656,20 +1761,21 @@ subroutine cutypen & end do end do - do levels=klevm1-1,klevm1-20,-1 ! loop starts + 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 + 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 @@ -1682,79 +1788,95 @@ subroutine cutypen & 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 + ! 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 + is = is + 1 endif enddo + if(is.eq.0) exit -! define the variables at the departure level + ! 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. + + 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 + 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) + + 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 + 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) + 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 + end if ! end if(jk .eq. levels) then -! the next levels, we use the variables at the first level as initial values + ! 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) + + ! calculate parcel entrainment rate for deep convection + fscale = min(1.,(pqsen(jl,jk)/pqsen(jl,levels))**3) ! (env. qvsat / (env. qvsat at cloud base))**3 + !eta(jl) = 1.75e-3 * (0.3-(min(1.,pqen(jl,jk) /pqsen(jl,jk))-1.)) * fscale ! entrainment rate + eta(jl) = entorg * fscale ! entrainment rate + dz(jl) = (pgeoh(jl,jk)-pgeoh(jl,jk+1)) * zrg ! convert from geopotential to height + !coef(jl) = eta(jl) * dz(jl) + coef(jl) = 0.5 * eta(jl) * dz(jl) + + ! dry static energy of environment 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)) + ! dry static energy for entraining parcel + dh(jl,jk) = (coef(jl)*(dhen(jl,jk+1)+dhen(jl,jk)) + (1.-coef(jl))*dh(jl,jk+1)) / (1.+coef(jl)) + ! mixing ratio for entraining parcel + pqu(jl,jk) = ( coef(jl)*(pqenh(jl,jk+1)+pqenh(jl,jk)) + (1.-coef(jl))*pqu(jl,jk+1) ) / (1.+coef(jl)) + ! temperature for entraining parcel ptu(jl,jk) = (dh(jl,jk)-pgeoh(jl,jk))*rcpd + zqold(jl) = pqu(jl,jk) - zph(jl)=paph(jl,jk) + zph(jl) = paph(jl,jk) + end if end do -! check if the parcel is saturated + + ! 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.) @@ -1763,18 +1885,19 @@ subroutine cutypen & (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 + ! 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 + 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) + 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 + + ! 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) @@ -1791,9 +1914,10 @@ subroutine cutypen & 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) + ! choose 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 @@ -1804,6 +1928,7 @@ subroutine cutypen & klab(jl,jk) = 2 kcbot(jl) = jk end if + end if if(kup(jl,jk) .lt. 0.)then @@ -1814,7 +1939,7 @@ subroutine cutypen & else lldcum(jl) = .false. end if - else + else if(plu(jl,jk) .gt. 0.)then klab(jl,jk)=2 else @@ -1831,7 +1956,7 @@ subroutine cutypen & ikb = kcbot(jl) ikt = kctop(jl) if(paph(jl,ikb) - paph(jl,ikt) < zdnoprc) lldcum(jl) = .false. - if(lldcum(jl)) then + if(lldcum(jl)) then ktype(jl) = 1 ldcum(jl) = .true. deepflag(jl) = .true. @@ -1876,15 +2001,17 @@ 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) + & (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 @@ -1956,68 +2083,71 @@ subroutine cuascn & ! kctop0 [ictop0] - estimate of cloud top. (cumastr) ! kcum [icum] - flag to control the call - integer klev,klon,klevp1,klevm1 - real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & - & puen(klon,klev), pven(klon,klev),& - & pten(klon,klev), pqen(klon,klev),& - & pgeo(klon,klev), pgeoh(klon,klevp1),& - & pap(klon,klev), paph(klon,klevp1),& - & pqsen(klon,klev), pqte(klon,klev),& - & pverv(klon,klev), pqsenh(klon,klev) - real(kind=kind_phys) ptu(klon,klev), pqu(klon,klev),& - & puu(klon,klev), pvu(klon,klev),& - & pmfu(klon,klev), zph(klon),& - & pmfub(klon), & - & pmfus(klon,klev), pmfuq(klon,klev),& - & plu(klon,klev), plude(klon,klev),& - & pmful(klon,klev), pdmfup(klon,klev) - real(kind=kind_phys) zdmfen(klon), zdmfde(klon),& - & zmfuu(klon), zmfuv(klon),& - & zpbase(klon), zqold(klon) - real(kind=kind_phys) phcbase(klon), zluold(klon) - real(kind=kind_phys) zprecip(klon), zlrain(klon,klev) - real(kind=kind_phys) zbuo(klon,klev), kup(klon,klev) - real(kind=kind_phys) wup(klon) - real(kind=kind_phys) wbase(klon), zodetr(klon,klev) - real(kind=kind_phys) plglac(klon,klev) - - real(kind=kind_phys) eta(klon),dz(klon) - - integer klwmin(klon), ktype(klon),& - & klab(klon,klev), kcbot(klon),& - & kctop(klon), kctop0(klon) - integer lndj(klon) - logical ldcum(klon), loflag(klon) - logical llo2,llo3, llo1(klon) - - integer kdpl(klon) - real(kind=kind_phys) zoentr(klon), zdpmean(klon) - real(kind=kind_phys) pdmfen(klon,klev), pmfude_rate(klon,klev) -! local variables - integer jl,jk - integer ikb,icum,itopm2,ik,icall,is,kcum,jlm,jll - integer jlx(klon) - real(kind=kind_phys) ztmst,zcons2,zfacbuo,zprcdgw,z_cwdrag,z_cldmax,z_cwifrac,z_cprc2 - real(kind=kind_phys) zmftest,zmfmax,zqeen,zseen,zscde,zqude - real(kind=kind_phys) zmfusk,zmfuqk,zmfulk - real(kind=kind_phys) zbc,zbe,zkedke,zmfun,zwu,zprcon,zdt,zcbf,zzco - real(kind=kind_phys) zlcrit,zdfi,zc,zd,zint,zlnew,zvw,zvi,zalfaw,zrold - real(kind=kind_phys) zrnew,zz,zdmfeu,zdmfdu,dp - real(kind=kind_phys) zfac,zbuoc,zdkbuo,zdken,zvv,zarg,zchange,zxe,zxs,zdshrd - real(kind=kind_phys) atop1,atop2,abot -!-------------------------------- -!* 1. specify parameters -!-------------------------------- - zcons2=3./(g*ztmst) - zfacbuo = 0.5/(1.+0.5) + !--- 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 -!--------------------------------- + !--------------------------------- + ! 2. set default values + !--------------------------------- llo3 = .false. do jl=1,klon zluold(jl)=0. @@ -2032,7 +2162,7 @@ subroutine cuascn & end if end do - ! initialize variout quantities + ! initialize variout quantities do jk=1,klev do jl=1,klon if(jk.ne.kcbot(jl)) plu(jl,jk)=0. @@ -2056,9 +2186,9 @@ subroutine cuascn & do jl = 1,klon if ( ktype(jl) == 3 ) ldcum(jl) = .false. end do -!------------------------------------------------ -! 3.0 initialize values at cloud base level -!------------------------------------------------ + !------------------------------------------------ + ! 3.0 initialize values at cloud base level + !------------------------------------------------ do jl=1,klon kctop(jl)=kcbot(jl) if(ldcum(jl)) then @@ -2070,25 +2200,24 @@ subroutine cuascn & 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 -!----------------------------------------------------------------- -! + !----------------------------------------------------------------- + ! 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 -! --------------------------------------------------------------------- + ! ----------------------------------------------------- + ! 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,& + & (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 @@ -2119,40 +2248,67 @@ subroutine cuascn & end do if(is.gt.0) llo3 = .true. -! -!* specify entrainment rates in *cuentr* -! ------------------------------------- - ik=jk + ! + !* 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 + ! ------------------------------------------------------- + ! 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)) + !--------------------------------------- + ! Entrainment parameter at cloud base + ! Why is it negative? + !--------------------------------------- if ( jk == kcbot(jl) ) then - zoentr(jl) = -1.75e-3*(min(1.,pqen(jl,jk)/pqsen(jl,jk)) - & + + zoentr(jl) = -entorg*(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 within the cloud layer + !--------------------------------------- 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) + + ! current level's entrainment is equal to zoentr value of level below zdmfen(jl) = zoentr(jl) + !--------------------------------------- + ! Set entrainment/detrainment rates for + ! shallow or mid-level convection. This + ! overwrites the values from the call + ! to cuentr. + !--------------------------------------- if ( ktype(jl) >= 2 ) then + ! double the entrainment rate for shallow convection zdmfen(jl) = 2.0*zdmfen(jl) + ! set turbulent detrainment equal to entrainment for shallow convection zdmfde(jl) = zdmfen(jl) end if - zdmfde(jl) = zdmfde(jl) * & - (1.6-min(1.,pqen(jl,jk)/pqsen(jl,jk))) + !--------------------------------------- + ! Multiply detrainment rate by (1.6-RH) + ! (Eq. 6.8/6.9 IFS Cy48r1) + ! For deep convection, will be value + ! from call to cuentr. + ! For shallow convection, value is set above. + !--------------------------------------- + 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.) @@ -2160,28 +2316,45 @@ subroutine cuascn & 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) + !--------------------------------------- + ! Calculate rates of change of + ! mass flux due to entrainment/detrainment + ! for state variables within updrafts + !--------------------------------------- + pdmfen(jl,jk) = zdmfen(jl) - zdmfde(jl) ! net rate = entrainment minus detrainment + pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl) ! massflux = massflux at lower level + entr - detr + + zqeen = pqenh(jl,jk+1) * zdmfen(jl) ! entr rate of env moisture into parcel + zseen = (cpd*ptenh(jl,jk+1)+pgeoh(jl,jk+1)) * zdmfen(jl) ! entr rate of env dry static energy into parcel + zscde = (cpd*ptu(jl,jk+1)+pgeoh(jl,jk+1)) * zdmfde(jl) ! detr rate of dry static energy from parcel + zqude = pqu(jl,jk+1) * zdmfde(jl) ! detr rate of moisture from parcel + plude(jl,jk) = plu(jl,jk+1) * zdmfde(jl) ! detr rate of cloud liquid water from parcel + + zmfusk = pmfus(jl,jk+1) + zseen - zscde ! net flux of updraft dry static energy at index k + zmfuqk = pmfuq(jl,jk+1) + zqeen - zqude ! net flux of updraft moisture at index k + zmfulk = pmful(jl,jk+1) - plude(jl,jk) ! net flux of updraft cloud liquid water at index k + !--------------------------------------- + ! Update updraft properties + ! due to entrainment/detrainment + !--------------------------------------- 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) + ptu(jl,jk) = max(100.,ptu(jl,jk)) ! updraft can't get colder than 100 K + ptu(jl,jk) = min(400.,ptu(jl,jk)) ! updraft can't get warmer than 400 K + + zqold(jl) = pqu(jl,jk) ! store parcel humidity, used later to determine how much + ! cloud water to condense after adjusting 'pqu' for saturation + 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 + !--------------------------------------- + ! Reset parcel values to environmental + ! values if below departure level + !--------------------------------------- do jl = 1,klon if ( jk > kdpl(jl) ) then ptu(jl,jk) = ptenh(jl,jk) @@ -2190,16 +2363,19 @@ subroutine cuascn & 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 -! + !------------------------------------------------ + ! Do corrections for moist ascent + ! by adjusting t,q and l in *cuadjtq* + ! to account for condensation + !------------------------------------------------ + ik = jk + icall = 1 ! flag for condensation in updrafts + if ( jlm > 0 ) then call cuadjtqn(klon,klev,ik,zph,ptu,pqu,loflag,icall) end if -! compute the upfraft speed in cloud layer + + ! compute the updraft speed in cloud layer do jll = 1 , jlm jl = jlx(jll) if ( pqu(jl,jk) /= zqold(jl) ) then @@ -2209,16 +2385,20 @@ subroutine cuascn & ptu(jl,jk) = ptu(jl,jk) + ralfdcp*plglac(jl,jk) end if end do + do jll = 1 , jlm jl = jlx(jll) + !------------------------------------------------ + ! If condensation has occurred + !------------------------------------------------ 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) - & + klab(jl,jk) = 2 ! we are now inside the cloud + plu(jl,jk) = plu(jl,jk) + zqold(jl) - pqu(jl,jk) ! add condensed water vapor to cloud water + zbc = ptu(jl,jk)*(1.+vtmpc1*pqu(jl,jk)-plu(jl,jk+1) - & ! parcel virtual temperature 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 + zbe = ptenh(jl,jk)*(1.+vtmpc1*pqenh(jl,jk)) ! environmental virtual temperature + zbuo(jl,jk) = zbc - zbe ! parcel buoyancy (K) + ! 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. @@ -2231,6 +2411,9 @@ subroutine cuascn & plu(jl,jk) = 0. end if end if + !---------------------------------------------- + ! If layer below is within the cloud layer + !---------------------------------------------- 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)) @@ -2241,7 +2424,8 @@ subroutine cuascn & (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 + + ! 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))) @@ -2249,31 +2433,60 @@ subroutine cuascn & 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) + !---------------------------------------------- + ! Organized detrainment for negatively buoyant + ! updraft (generally at cloud top) based + ! on the decrease of updraft velocity with height + ! (Eq. 6.12 IFS Cy48r1 without RH term) + ! + ! Is stable -> no org. entrainment. This overwrites + ! PMFU for current level which has been calculated + ! above with organised entrainment (ICON comment) + !---------------------------------------------- if ( zbuo(jl,jk) < 0. ) then zkedke = kup(jl,jk)/max(1.e-10,kup(jl,jk+1)) zkedke = max(0.,min(1.,zkedke)) - zmfun = sqrt(zkedke)*pmfu(jl,jk+1) !* (1.6-min(1.,pqen(jl,jk) / & - ! pqsen(jl,jk))) + 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) + ! mass flux = mass flux at layer below plus entr minus detr pmfu(jl,jk) = pmfu(jl,jk+1) + zdmfen(jl) - zdmfde(jl) end if + !---------------------------------------------- + ! Calculate parcel entrainment rate given + ! a sufficiently buoyant updraft, otherwise + ! set to zero (Eq. 6.7 IFS Cy48r1) + !---------------------------------------------- if ( zbuo(jl,jk) > -0.2 ) then + ! when positively buoyant, have organised entrainment + ! which increases MF with height, while detrainment + ! is pretty small and constant (ICON comment) ikb = kcbot(jl) - zoentr(jl) = 1.75e-3*(0.3-(min(1.,pqen(jl,jk-1) / & + ! zoentr is overwritten, but not used until + ! the next jk level in the loop (ICON comment) + zoentr(jl) = entorg*(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 + + ! 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 + !------------------------------------------- + ! determine convection top level; + ! the last set of criteria serves to limit + ! the overshooting of updrafts + ! through the tropopause (ICON comment) + !------------------------------------------- if ( kup(jl,jk) > 0. .and. pmfu(jl,jk) > 0. ) then kctop(jl) = jk llo1(jl) = .true. @@ -2284,7 +2497,7 @@ subroutine cuascn & zdmfde(jl) = pmfu(jl,jk+1) plude(jl,jk) = plu(jl,jk+1)*zdmfde(jl) end if -! save detrainment rates for updraught + ! 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 @@ -2296,12 +2509,14 @@ subroutine cuascn & pmfude_rate(jl,jk) = zdmfde(jl) end if end do - + !---------------------------------------------------- + ! Calculate precipitation rates + !---------------------------------------------------- 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. + ! 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 @@ -2309,10 +2524,9 @@ subroutine cuascn & end if ikb=kcbot(jl) if ( plu(jl,jk) > zdshrd )then -! if ((paph(jl,ikb)-paph(jl,jk))>zdnoprc) then zwu = min(15.0,sqrt(2.*max(0.1,kup(jl,jk+1)))) zprcon = zprcdgw/(0.75*zwu) -! PARAMETERS FOR BERGERON-FINDEISEN PROCESS (T < -5C) + ! 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 @@ -2336,6 +2550,7 @@ subroutine cuascn & end if end if end do + do jl = 1, klon if ( llo1(jl) ) then if ( zlrain(jl,jk) > 0. ) then @@ -2354,17 +2569,20 @@ subroutine cuascn & 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 -! ------------------ + + !---------------------------------------------------------------------- + ! 5. final calculations + ! ------------------ do jl = 1,klon if ( kctop(jl) == -1 ) ldcum(jl) = .false. kcbot(jl) = max(kcbot(jl),kctop(jl)) @@ -2377,164 +2595,168 @@ subroutine cuascn & return end subroutine cuascn !--------------------------------------------------------- -! level 3 souroutines +! level 3 souroutines !-------------------------------------------------------- - subroutine cudlfsn & - & (klon, klev, & - & kcbot, kctop, lndj, ldcum, & - & ptenh, pqenh, puen, pven, & - & pten, pqsen, pgeo, & - & pgeoh, paph, ptu, pqu, plu,& - & puu, pvu, pmfub, prfl, & - & ptd, pqd, pud, pvd, & - & pmfd, pmfds, pmfdq, pdmfdp, & - & kdtop, lddraf) - -! this routine calculates level of free sinking for -! cumulus downdrafts and specifies t,q,u and v values - -! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89 - -! purpose. -! -------- -! to produce lfs-values for cumulus downdrafts -! for massflux cumulus parameterization - -! interface -! --------- -! this routine is called from *cumastr*. -! input are environmental values of t,q,u,v,p,phi -! and updraft values t,q,u and v and also -! cloud base massflux and cu-precipitation rate. -! it returns t,q,u and v values and massflux at lfs. - -! method. - -! check for negative buoyancy of air of equal parts of -! moist environmental air and cloud air. - -! parameter description units -! --------- ----------- ----- -! input parameters (integer): - -! *klon* number of grid points per packet -! *klev* number of levels -! *kcbot* cloud base level -! *kctop* cloud top level - -! input parameters (logical): - -! *lndj* land sea mask (1 for land) -! *ldcum* flag: .true. for convective points - -! input parameters (real(kind=kind_phys)): - -! *ptenh* env. temperature (t+1) on half levels k -! *pqenh* env. spec. humidity (t+1) on half levels kg/kg -! *puen* provisional environment u-velocity (t+1) m/s -! *pven* provisional environment v-velocity (t+1) m/s -! *pten* provisional environment temperature (t+1) k -! *pqsen* environment spec. saturation humidity (t+1) kg/kg -! *pgeo* geopotential m2/s2 -! *pgeoh* geopotential on half levels m2/s2 -! *paph* provisional pressure on half levels pa -! *ptu* temperature in updrafts k -! *pqu* spec. humidity in updrafts kg/kg -! *plu* liquid water content in updrafts kg/kg -! *puu* u-velocity in updrafts m/s -! *pvu* v-velocity in updrafts m/s -! *pmfub* massflux in updrafts at cloud base kg/(m2*s) - -! updated parameters (real(kind=kind_phys)): - -! *prfl* precipitation rate kg/(m2*s) - -! output parameters (real(kind=kind_phys)): - -! *ptd* temperature in downdrafts k -! *pqd* spec. humidity in downdrafts kg/kg -! *pud* u-velocity in downdrafts m/s -! *pvd* v-velocity in downdrafts m/s -! *pmfd* massflux in downdrafts kg/(m2*s) -! *pmfds* flux of dry static energy in downdrafts j/(m2*s) -! *pmfdq* flux of spec. humidity in downdrafts kg/(m2*s) -! *pdmfdp* flux difference of precip. in downdrafts kg/(m2*s) - -! output parameters (integer): - -! *kdtop* top level of downdrafts - -! output parameters (logical): - -! *lddraf* .true. if downdrafts exist - -! externals -! --------- -! *cuadjtq* for calculating wet bulb t and q at lfs -!---------------------------------------------------------------------- - implicit none - - integer klev,klon - real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & - & puen(klon,klev), pven(klon,klev), & - & pten(klon,klev), pqsen(klon,klev), & - & pgeo(klon,klev), & - & pgeoh(klon,klev+1), paph(klon,klev+1),& - & ptu(klon,klev), pqu(klon,klev), & - & puu(klon,klev), pvu(klon,klev), & - & plu(klon,klev), & - & pmfub(klon), prfl(klon) - - real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), & - & pud(klon,klev), pvd(klon,klev), & - & pmfd(klon,klev), pmfds(klon,klev), & - & pmfdq(klon,klev), pdmfdp(klon,klev) - integer kcbot(klon), kctop(klon), & - & kdtop(klon), ikhsmin(klon) - logical ldcum(klon), & - & lddraf(klon) - integer lndj(klon) - - real(kind=kind_phys) ztenwb(klon,klev), zqenwb(klon,klev), & - & zcond(klon), zph(klon), & - & zhsmin(klon) - logical llo2(klon) -! local variables - integer jl,jk - integer is,ik,icall,ike - real(kind=kind_phys) zhsk,zttest,zqtest,zbuo,zmftop - -!---------------------------------------------------------------------- - -! 1. set default values for downdrafts -! --------------------------------- - do jl=1,klon - lddraf(jl)=.false. - kdtop(jl)=klev+1 + 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) -! ---------------------------------------------------- + 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) + & @@ -2547,211 +2769,215 @@ subroutine cudlfsn & 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 + 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 +! level 3 souroutines !-------------------------------------------------------- !********************************************** ! subroutine cuddrafn !********************************************** - subroutine cuddrafn & - & ( klon, klev, lddraf & - & , ptenh, pqenh, puen, pven & - & , pgeo, pgeoh, paph, prfl & - & , ptd, pqd, pud, pvd, pmfu & - & , pmfd, pmfds, pmfdq, pdmfdp, pmfdde_rate ) - -! this routine calculates cumulus downdraft descent - -! m.tiedtke e.c.m.w.f. 12/86 modif. 12/89 - -! purpose. -! -------- -! to produce the vertical profiles for cumulus downdrafts -! (i.e. t,q,u and v and fluxes) - -! interface -! --------- - -! this routine is called from *cumastr*. -! input is t,q,p,phi,u,v at half levels. -! it returns fluxes of s,q and evaporation rate -! and u,v at levels where downdraft occurs - -! method. -! -------- -! calculate moist descent for entraining/detraining plume by -! a) moving air dry-adiabatically to next level below and -! b) correcting for evaporation to obtain saturated state. - -! parameter description units -! --------- ----------- ----- -! input parameters (integer): - -! *klon* number of grid points per packet -! *klev* number of levels - -! input parameters (logical): - -! *lddraf* .true. if downdrafts exist - -! input parameters (real(kind=kind_phys)): - -! *ptenh* env. temperature (t+1) on half levels k -! *pqenh* env. spec. humidity (t+1) on half levels kg/kg -! *puen* provisional environment u-velocity (t+1) m/s -! *pven* provisional environment v-velocity (t+1) m/s -! *pgeo* geopotential m2/s2 -! *pgeoh* geopotential on half levels m2/s2 -! *paph* provisional pressure on half levels pa -! *pmfu* massflux updrafts kg/(m2*s) - -! updated parameters (real(kind=kind_phys)): - + 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 +! *pgeoh* geopotential on half levels 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(kind=kind_phys)): - -! *ptd* temperature in downdrafts k -! *pqd* spec. humidity in downdrafts kg/kg -! *pud* u-velocity in downdrafts m/s -! *pvd* v-velocity in downdrafts m/s + +! 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 -!---------------------------------------------------------------------- + +! externals +! --------- +! *cuadjtq* for adjusting t and q due to evaporation in +! saturated descent +!---------------------------------------------------------------------- implicit none - - integer klev,klon - real(kind=kind_phys) ptenh(klon,klev), pqenh(klon,klev), & - & puen(klon,klev), pven(klon,klev), & - & pgeoh(klon,klev+1), paph(klon,klev+1), & - & pgeo(klon,klev), pmfu(klon,klev) - - real(kind=kind_phys) ptd(klon,klev), pqd(klon,klev), & - & pud(klon,klev), pvd(klon,klev), & - & pmfd(klon,klev), pmfds(klon,klev), & - & pmfdq(klon,klev), pdmfdp(klon,klev), & - & prfl(klon) - real(kind=kind_phys) pmfdde_rate(klon,klev) - logical lddraf(klon) - - real(kind=kind_phys) zdmfen(klon), zdmfde(klon), & - & zcond(klon), zoentr(klon), & - & zbuoy(klon) - real(kind=kind_phys) zph(klon) - logical llo2(klon) - logical llo1 -! local variables - integer jl,jk - integer is,ik,icall,ike, itopde(klon) - real(kind=kind_phys) zentr,zdz,zzentr,zseen,zqeen,zsdde,zqdde,zdmfdp - real(kind=kind_phys) zmfdsk,zmfdqk,zbuo,zrain,zbuoyz,zmfduk,zmfdvk - -!---------------------------------------------------------------------- -! 1. calculate moist descent for cumulus downdraft by -! (a) calculating entrainment/detrainment rates, -! including organized entrainment dependent on -! negative buoyancy and assuming -! linear decrease of massflux in pbl -! (b) doing moist descent - evaporative cooling -! and moistening is calculated in *cuadjtq* -! (c) checking for negative buoyancy and -! specifying final t,q,u,v and downward fluxes -! ------------------------------------------------- - do jl=1,klon - zoentr(jl)=0. - zbuoy(jl)=0. - zdmfen(jl)=0. + +!--- 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 + 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 + 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 - + zdmfen(jl)=zentr + zdmfde(jl)=zentr + endif + enddo + do jl=1,klon if(llo2(jl)) then if(jk.gt.itopde(jl)) then @@ -2777,182 +3003,195 @@ subroutine cuddrafn & endif enddo - do jl=1,klon - if(llo2(jl)) then + 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 + 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=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 + endif + if(zbuo.ge.0 .or. prfl(jl).le.(pmfd(jl,jk)*zcond(jl))) then pmfd(jl,jk)=0. - zbuo=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) + 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)) + zoentr(jl)=g*zbuoyz*0.5/(1.+zbuoy(jl)) pmfdde_rate(jl,jk) = -zdmfde(jl) - endif + endif enddo - - enddo - - return + + enddo + + return end subroutine cuddrafn !--------------------------------------------------------- -! level 3 souroutines +! level 3 souroutines !-------------------------------------------------------- - subroutine cuflxn & - & ( klon, klev, ztmst & - & , pten, pqen, pqsen, ptenh, pqenh & - & , paph, pap, pgeoh, lndj, ldcum & - & , kcbot, kctop, kdtop, ktopm2 & - & , ktype, lddraf & - & , pmfu, pmfd, pmfus, pmfds & - & , pmfuq, pmfdq, pmful, plude & - & , pdmfup, pdmfdp, pdpmel, plglac & - & , prain, pmfdde_rate, pmflxr, pmflxs ) - -! m.tiedtke e.c.m.w.f. 7/86 modif. 12/89 - -! purpose -! ------- - -! this routine does the final calculation of convective -! fluxes in the cloud layer and in the subcloud layer - -! interface -! --------- -! this routine is called from *cumastr*. - - -! parameter description units -! --------- ----------- ----- -! input parameters (integer): - -! *klon* number of grid points per packet -! *klev* number of levels -! *kcbot* cloud base level -! *kctop* cloud top level -! *kdtop* top level of downdrafts - -! input parameters (logical): - -! *lndj* land sea mask (1 for land) -! *ldcum* flag: .true. for convective points - -! input parameters (real(kind=kind_phys)): - -! *ztmst* time step for the physics s -! *pten* provisional environment temperature (t+1) k -! *pqen* provisional environment spec. humidity (t+1) kg/kg -! *pqsen* environment spec. saturation humidity (t+1) kg/kg -! *ptenh* env. temperature (t+1) on half levels k -! *pqenh* env. spec. humidity (t+1) on half levels kg/kg -! *paph* provisional pressure on half levels pa -! *pap* provisional pressure on full levels pa -! *pgeoh* geopotential on half levels m2/s2 - -! updated parameters (integer): - -! *ktype* set to zero if ldcum=.false. - -! updated parameters (logical): - -! *lddraf* set to .false. if ldcum=.false. or kdtop