!> \file cs_conv.F90 !! This file contains the Chikira-Sugiyama Convection scheme. module cs_conv !--------------------------------------------------------------------------------- ! Purpose: ! !>--------------------------------------------------------------------------------- ! Purpose: ! !> Interface for Chikira-Sugiyama convection scheme !! !! Author: Minoru Chikira !--------------------------------------------------------------------------------- ! use machine , only : kind_phys use physcons, only : cp => con_cp, grav => con_g, & & rair => con_rd, rvap => con_rv, & & cliq => con_cliq, cvap => con_cvap, & & epsv => con_eps, epsvm1 => con_epsm1, & & epsvt => con_fvirt, & & el => con_hvap, emelt => con_hfus, t0c => con_t0c use funcphys, only : fpvs ! this is saturation vapor pressure in funcphys.f implicit none private ! Make default type private to the module real(kind_phys), parameter :: zero=0.0d0, one=1.0d0, half=0.5d0 real(kind_phys), parameter :: cpoel=cp/el, cpoesub=cp/(el+emelt), esubocp=1.0/cpoesub, & elocp=el/cp, oneocp=one/cp, gocp=grav/cp, gravi=one/grav,& emeltocp=emelt/cp, cpoemelt=cp/emelt, epsln=1.e-10_kind_phys real(kind_phys), parameter :: fact1=(cvap-cliq)/rvap, fact2=el/rvap-fact1*t0c !< to calculate d(qs)/dT logical, parameter :: adjustp=.true. ! logical, parameter :: adjustp=.false. ! Tuning parameters set from namelist ! ! real(kind_phys), parameter, public :: CLMD = 0.60, & !< entrainment efficiency (now thru argument) real(kind_phys), parameter, public :: & PA=0.15, & !< factor for buoyancy to affect updraft velocity CPRES = 0.55, & !< pressure factor for momentum transport ALP0 = 5.0e7, & !< alpha parameter in prognostic closure ! ALP0 = 8.0e7, & !< alpha parameter in prognostic closure ! CLMP = (one-CLMD)*(PA+PA), & ! CLMDPA = CLMD*PA, & spblmin=0.05, & !< minimum cloudbase height in p/ps spblmax=0.30, & !< maximum cloudbase height in p/ps ! spblcrit=0.03, & !< minimum cloudbase height in p/ps ! spblcrit=0.035,& !< minimum cloudbase height in p/ps ! spblcrit=0.025,& !< minimum cloudbase height in p/ps cincrit= -150.0 ! cincrit= -120.0 ! cincrit= -100.0 !DD precz0 and preczh control partitioning of water between detrainment !DD and precipitation. Decrease for more precip real(kind_phys), public :: precz0, preczh, clmd, clmp, clmdpa ! ! Private data ! real(kind_phys), parameter :: unset_kind_phys = -999._kind_phys ! missing value ! integer :: iulog !< unit to write debugging and diagnostic output !DD Note - see if I can find corresponding variable in a GFS module ! ! Shared variables ! integer, parameter :: ITI = 2, ITL = 3 !< index of ice and liquid water integer, save, dimension(50) :: IMFXR !< 0: mass fixer is not applied !! tracers which may become negative !! values e.g. subgrid-PDFs !! 1: mass fixer is applied, total mass !! may change through cumulus scheme !! e.g. moisture, liquid cloud, ice !! cloud, aerosols !! 2: mass fixer is applied, total mass !! never change through cumulus scheme !! e.g. CO2 ! PUBLIC: interfaces ! public cs_conv_run ! CS scheme main driver contains !>\defgroup cs_scheme Chikira-Sugiyama Cumulus Scheme Module !> \brief The subroutine contains the main driver for Chikira-Sugiyama convective scheme. !! !! \author Minoru Chikira !! !! History: !! - Jun 26 2014 : D. Dazlich - Modified for GFS !! - Apr 10 2015 : S. Moorthi - check for allocatable arrays and fix argument for cbmfx !! - Oct 2015 : D. Dazlich - Add computation of updraft area fraction (sigma) for !! diagnostic purposes. !! - Aug 2016 : D. Dazlich - Create flux form of tendencies and multiply by !! Arakawa-Wu functions of sigma !! - Sep 2016 : S. Moorthi - found two bugs - cleanup and some optimization !! - Oct 2016 : S. Moorthi - added sigma affects on tracers and CUMFLX and CUMDET !! made many cosmetic changes !! - Nov 2016 : S. Moorthi - further optimization and cleanup and several bug fixes !! - April 2017 : S. moorthi - many changes including removing elam and making gcym !! a function of cloud type. This makes it possible for !! AW affect propagate to other routines such as CUMUPR !! - Apr 12, 2017 : S. Moorthi - Added flx_form logical and relevant code to compute AW !! without flux form when false. !! - May 17, 2017 : S. Moorthi - Added routine CUMSBW for just momentum change !! in advective form !! - Sep 08, 2017 : D. Dazlich - tracers in flux form for AW !! - Nov 2017 : S. Moorthi - fix some bugs and fix fluxform for tracers !! - Nov 22 2017 : S. Moorthi - add kcnv array to identify points where deep convection !! operates - 0 - no convection 1 - with convection !! - Jan 30 2018 : S. Moorthi - fixed sigmad dimension error in CUMDWN and an error when adjustp=.true. !! - May 2018 : S. Moorthi - modified cumup to compute total workfunction (positive plus negative) !! and negative part only and to let a particular ensemble exist only if !! the ratio of negative to total is less than some prescribed percent. !! Also, added an extra iteration in this k loop. Reduced some memory. !! - June 2018 : S. Moorthi - the output mass fluxes ud_mf, dd_mf and dt_mf are over time step delta !! !! \b Arakawa-Wu \b implemtation: !! for background, consult An Introduction to the !! General Circulation of the Atmosphere, Randall, chapter six. !! Traditional parameterizations compute tendencies like those in eq 103, 105 and 106. !! Because Arakawa-Wu applies different functions to different components to the !! terms within these equations, it requires the terms used in alternate eqns 91 - 93. !! The code required to compute these terms is added within, and the appropriate !! functions of updraft area fraction (sigma) are applied. Thus, AW requires three !! steps: !! -# computation of the updraft area fraction !! -# alternative representation of the tendency terms !! -# application of functions of sigma to the alternative tendency terms !! here, and in gbphys to the large-scale microphysics tendencies. !! !! The bulk of AW is implemented within subroutine CS_CUMLUS(), and the routines it calls. !! !! !! JLS NOTE: The convective mass fluxes (dt_mf, dd_mf and ud_mf) passed in and out of cs_conv have not been multiplied by !! the timestep (kg/m2/sec) as they are in all other convective schemes. EMC is aware of this problem, !! and in the future will be fixing this discrepancy. In the meantime, CCPP will use the same mass flux standard_name !! and long_name as the other convective schemes, where the units are in kg/m2. (Aug 2018) !! !! \section arg_table_cs_conv_run Argument Table !! \htmlinclude cs_conv_run.html !! !! \section general_cs_conv CS Convection Scheme General Algorithm !> @{ subroutine cs_conv_run( IJSDIM , KMAX , ntracp1 , NN, & NTR , nctp , & !DD dimensions otspt , lat , kdt , & t , q , rain1 , clw , & zm , zi , pap , paph , & delta , delti , ud_mf , dd_mf , dt_mf, & u , v , fscav , fswtr, & cbmfx , mype , wcbmaxm , precz0in, preczhin, & clmdin , sigma , do_aw , do_awdd , flx_form, & lprnt , ipr, kcnv, & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & ! for coupling to MG microphysics CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE, & mp_phys,errmsg,errflg) implicit none ! ! input arguments ! INTEGER, INTENT(IN) :: IJSDIM, KMAX, ntracp1, nn, NTR, mype, nctp, mp_phys, kdt, lat !! DD, for GFS, pass in logical, intent(in) :: otspt(:,:) ! otspt(:,1) - on/off switch for tracer transport by updraft and ! downdraft. should not include subgrid PDF and turbulence ! otspt(:,2) - on/off switch for tracer transport by subsidence ! should include subgrid PDF and turbulence real(kind_phys), intent(inout) :: t(:,:) ! temperature at mid-layer (K) real(kind_phys), intent(inout) :: q(:,:) ! water vapor array including moisture (kg/kg) real(kind_phys), intent(inout) :: clw(:,:,:) ! tracer array including cloud condensate (kg/kg) real(kind_phys), intent(in) :: pap(:,:) ! pressure at mid-layer (Pa) real(kind_phys), intent(in) :: paph(:,:) ! pressure at boundaries (Pa) real(kind_phys), intent(in) :: zm(:,:) ! geopotential at mid-layer (m) real(kind_phys), intent(in) :: zi(:,:) ! geopotential at boundaries (m) real(kind_phys), intent(in) :: fscav(:), fswtr(:), wcbmaxm(:) real(kind_phys), intent(in) :: precz0in, preczhin, clmdin ! added for cs_convr real(kind_phys), intent(inout) :: u(:,:) ! zonal wind at mid-layer (m/s) real(kind_phys), intent(inout) :: v(:,:) ! meridional wind at mid-layer (m/s) real(kind_phys), intent(in) :: DELTA ! physics time step real(kind_phys), intent(in) :: DELTI ! dynamics time step (model time increment in seconds) logical, intent(in) :: do_aw, do_awdd, flx_form ! ! modified arguments ! real(kind_phys), intent(inout) :: CBMFX(:,:) ! cloud base mass flux (kg/m2/s) ! ! output arguments ! ! updraft, downdraft, and detrainment mass flux (kg/m2/s) real(kind_phys), intent(inout), dimension(:,:) :: ud_mf, dd_mf, dt_mf real(kind_phys), intent(out) :: rain1(:) ! lwe thickness of deep convective precipitation amount (m) ! GJF* These variables are conditionally allocated depending on whether the ! Morrison-Gettelman microphysics is used, so they must be declared ! using assumed shape. real(kind_phys), intent(out), dimension(:,:) :: qlcn, qicn, w_upi,cnv_mfd, & cnv_dqldt, clcn, cnv_fice, & cnv_ndrop, cnv_nice, cf_upi ! *GJF logical, intent(in) :: lprnt integer, intent(in) :: ipr integer, intent(inout) :: kcnv(:) ! zero if no deep convection and 1 otherwise character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg !DDsigma - output added for AW sigma diagnostics ! interface sigma and vertical velocity by cloud type (1=sfc) ! real(kind_phys), intent(out), dimension(:,:,:) :: sigmai, vverti real(kind_phys), intent(out), dimension(:,:) :: sigma ! sigma sigma totaled over cloud type - on interfaces (1=sfc) ! sigma terms in eq 91 and 92 ! real(kind_phys), dimension(IJSDIM,KMAX) :: sfluxterm, qvfluxterm, condterm !DDsigma ! ! output arguments of CS_CUMLUS ! real(kind_phys), dimension(IJSDIM,KMAX,nctp) :: vverti real(kind_phys) GTT(IJSDIM,KMAX) !< temperature tendency [K/s] real(kind_phys) GTQ(IJSDIM,KMAX,NTR) !< tracer tendency [kg/kg/s] real(kind_phys) GTU(IJSDIM,KMAX) !< zonal velocity tendency [m/s2] real(kind_phys) GTV(IJSDIM,KMAX) !< meridional velocity tendency [m/s2] real(kind_phys) GTPRP(IJSDIM,KMAX) !< precipitation (including snowfall) flux at interfaces [kg/m2/s] real(kind_phys) GSNWP(IJSDIM,KMAX) !< snowfall flux at interfaces [kg/m2/s] integer KT(IJSDIM,nctp) !< cloud top index for each cloud type real(kind_phys) :: cape(IJSDIM) !< convective available potential energy (J/kg) real(kind_phys) :: prec(IJSDIM) !< precipitation at surface (including snowfall) (kg/m2/s) real(kind_phys) :: snow(IJSDIM) !< snowfall at surface (kg/m2/s) ! ! input arguments of CS_CUMLUS ! real(kind_phys) GDT(IJSDIM,KMAX) !< temperature [K] real(kind_phys) GDQ(IJSDIM,KMAX,NTR) !< tracers including moisture [kg/kg] !DDsigmadiag real(kind_phys) GDU(IJSDIM,KMAX) !< zonal wind [m/s] real(kind_phys) GDV(IJSDIM,KMAX) !< meridional wind [m/s] real(kind_phys) GDTM(IJSDIM,KMAX+1) !< temperature at boundaries of layers [K] real(kind_phys) GDP(IJSDIM,KMAX) !< pressure [Pa] real(kind_phys) GDPM(IJSDIM,KMAX+1) !< pressure at boundaries of layers [Pa] real(kind_phys) GDZ(IJSDIM,KMAX) !< altitude [m] real(kind_phys) GDZM(IJSDIM,KMAX+1) !< altitude at boundaries of layers [m] real(kind_phys) delp(IJSDIM,KMAX) !< pressure difference between layers [Pa] real(kind_phys) delpi(IJSDIM,KMAX) !< grav/delp ! ! local variables ! !DD real(kind_phys) :: zs(IJSDIM) !< surface height [m] integer KTMAX(IJSDIM) !< max of KT real(kind_phys) :: ftintm, wrk, wrk1, tem integer i, k, n, ISTS, IENS, kp1 !DD borrowed from RAS to go form total condensate to ice/water separately ! parameter (tf=130.16, tcr=160.16, tcrf=1.0/(tcr-tf),tcl=2.0) ! parameter (tf=230.16, tcr=260.16, tcrf=1.0/(tcr-tf)) real(kind_phys), parameter :: tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf), tcl=2.0 logical, save :: first=.true. ! Initialize CCPP error handling variables errmsg = '' errflg = 0 ! lprnt = kdt == 1 .and. mype == 38 ! ipr = 43 precz0 = precz0in preczh = preczhin clmd = clmdin CLMP = (one-CLMD)*(PA+PA) CLMDPA = CLMD*PA ! if (first) then do i=1,ntr IMFXR(i) = 0 enddo ! IMFXR(1) = 1 ! IMFXR(ITL) = 1 ! IMFXR(ITI) = 1 first = .false. endif ! ISTS = 1 IENS = IJSDIM do k=1,KMAX+1 do i=1,IJSDIM GDZM(i,k) = zi(i,k) * gravi GDPM(i,k) = paph(i,k) enddo enddo do k=1,KMAX do i=1,IJSDIM GDT(i,k) = t(i,k) GDU(i,k) = u(i,k) GDV(i,k) = v(i,k) GDZ(i,k) = zm(i,k) * gravi GDP(i,k) = pap(i,k) GDQ(i,k,1) = q(i,k) delp(i,k) = paph(i,k) - paph(i,k+1) delpi(i,k) = grav / delp(i,k) enddo enddo !DD following adapted from ras !> -# Following the Relaxed Arakawa Schubert Scheme (RAS; !! Moorthi and Suarez 1992 \cite moorthi_and_suarez_1992 ), !! separate total condensate between ice and water. !! The ratio of cloud ice to cloud water is determined by a linear function !! of temperature: !!\f[ !! F_i(T)= (T_2-T)/(T_2-T_1) !!\f] !! where T is temperature, and\f$T_1\f$ and \f$T_2\f$ are set as tcf=263.16 !! and tf= 233.16 if (clw(1,1,2) <= -999.0) then ! input ice/water are together do k=1,kmax do i=1,IJSDIM tem = clw(i,k,1) * MAX(ZERO, MIN(ONE, (TCR-t(i,k))*TCRF)) clw(i,k,2) = clw(i,k,1) - tem clw(i,k,1) = tem enddo enddo endif !DD end ras adaptation do k=1,kmax do i=1,ijsdim tem = min(clw(i,k,1), 0.0) wrk = min(clw(i,k,2), 0.0) clw(i,k,1) = clw(i,k,1) - tem clw(i,k,2) = clw(i,k,2) - wrk gdq(i,k,1) = gdq(i,k,1) + tem + wrk enddo enddo ! if (lprnt) write(0,*)'in cs clw1b=',clw(ipr,:,1),' kdt=',kdt ! if (lprnt) write(0,*)'in cs clw2b=',clw(ipr,:,2),' kdt=',kdt do n=2,NTR do k=1,KMAX do i=1,IJSDIM GDQ(i,k,n) = clw(i,k,n-1) enddo enddo enddo ! if (lprnt) write(0,*)' incs tke=',gdq(ipr,1:25,ntr) ! !*************************************************************************************** ! !> -# Calculate temperature at interfaces ! DO K=2,KMAX DO I=ISTS,IENS wrk = one / GDP(I,K) wrk1 = one / LOG(GDP(I,K-1)*wrk) FTINTM = wrk1 * LOG(GDPM(I,K)*wrk) GDTM(I,K) = FTINTM*GDT(I,K-1) + (one-FTINTM)*GDT(I,K) ENDDO ENDDO DO I=ISTS,IENS GDTM(I,KMAX+1) = GDT(I,KMAX) GDTM(I,1) = GDT(I,1) ! Is this a good approximation ? - Moorthi ENDDO !> -# Initialize the sigma diagnostics do n=1,nctp do k=1,kmax do i=ists,iens vverti(i,k,n) = zero enddo enddo enddo do k=1,kmax do i=ists,iens sigma(i,k) = zero enddo enddo ! !> -# Call cs_cumlus() for the main CS cumulus parameterization call CS_CUMLUS (IJSDIM, IJSDIM, KMAX , NTR , & !DD dimensions otspt(1:ntr,1), otspt(1:ntr,2), & lprnt , ipr , & GTT , GTQ , GTU , GTV , & ! output dt_mf , & ! output GTPRP , GSNWP , ud_mf , & ! output dd_mf , cape , KT , & ! output CBMFX , & ! modified GDT , GDQ , GDU , GDV , & ! input GDTM , & ! input GDP , GDPM , GDZ , GDZM , & ! input delp , delpi , & DELTA , DELTI , ISTS , IENS, mype,& ! input fscav, fswtr, wcbmaxm, nctp, & sigma, vverti, & ! input/output !DDsigma do_aw, do_awdd, flx_form) ! ! !DD detrainment has to be added in for GFS ! ! if (lprnt) write(0,*)' aft cs_cum gtqi=',gtq(ipr,:,2) ! if (lprnt) write(0,*)' aft cs_cum gtql=',gtq(ipr,:,3) do n=2,NTR do k=1,KMAX do i=1,IJSDIM clw(i,k,n-1) = max(zero, GDQ(i,k,n) + GTQ(i,k,n) * delta) enddo enddo enddo ! if (lprnt) write(0,*)' aftcs_cum tkein=',gdq(ipr,1:25,ntr),' delta=',delta ! if (lprnt) write(0,*)' aftcs_cum tke=',clw(ipr,1:25,ntr-1) ! if (lprnt) write(0,*)'in cs clw1a=',clw(ipr,:,1),' kdt=',kdt ! if (lprnt) write(0,*)'in cs clw2a=',clw(ipr,:,2),' kdt=',kdt ! do k=1,KMAX do i=1,IJSDIM q(i,k) = max(zero, GDQ(i,k,1) + GTQ(i,k,1) * delta) t(i,k) = GDT(i,k) + GTT(i,k) * delta u(i,k) = GDU(i,k) + GTU(i,k) * delta v(i,k) = GDV(i,k) + GTV(i,k) * delta enddo enddo if (mp_phys == 10) then ! for 2M microphysics, always output these variables if (do_aw) then do k=1,KMAX kp1 = min(k+1,kmax) do i=1,IJSDIM qicn(i,k) = max(0.0, clw(i,k,1)-gdq(i,k,2)) qlcn(i,k) = max(0.0, clw(i,k,2)-gdq(i,k,3)) wrk = qicn(i,k) + qlcn(i,k) if (wrk > 1.0e-12) then cnv_fice(i,k) = qicn(i,k) / wrk else cnv_fice(i,k) = 0.0 endif ! ! CNV_MFD(i,k) = dt_mf(i,k) * (1.0/delta) CNV_MFD(i,k) = dt_mf(i,k) CNV_DQLDT(i,k) = wrk / delta ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 cf_upi(i,k) = max(0.0, min(1.0, 0.5*(sigma(i,k)+sigma(i,kp1)))) CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft !! clcn(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)/delta),0.25)) w_upi(i,k) = 0.0 !! w_upi(i,k) = ud_mf(i,k)*(t(i,k)+epsvt*gdq(i,k,1)) * rair & !! / (delta*max(cf_upi(i,k),1.e-12)*gdp(i,k)) enddo enddo do k=1,KMAX do i=1,ijsdim do n=1,nctp w_upi(i,k) = w_upi(i,k) + vverti(i,k,n) enddo if (sigma(i,k) > 1.0e-10) then w_upi(i,k) = w_upi(i,k) / sigma(i,k) else w_upi(i,k) = 0.0 endif enddo enddo else do k=1,KMAX do i=1,IJSDIM qicn(i,k) = max(0.0, clw(i,k,1)-gdq(i,k,2)) qlcn(i,k) = max(0.0, clw(i,k,2)-gdq(i,k,3)) cnv_fice(i,k) = qicn(i,k) / max(1.0e-10,qicn(i,k)+qlcn(i,k)) ! ! CNV_MFD(i,k) = dt_mf(i,k) * (1/delta) CNV_MFD(i,k) = dt_mf(i,k) CNV_DQLDT(i,k) = (qicn(i,k)+qlcn(i,k)) / delta ! CNV_PRC3(i,k) = 0.0 CNV_NDROP(i,k) = 0.0 CNV_NICE(i,k) = 0.0 cf_upi(i,k) = max(0.0,min(0.01*log(1.0+500*ud_mf(i,k)),0.25)) ! & 500*ud_mf(i,k)),0.60)) CLCN(i,k) = cf_upi(i,k) !downdraft is below updraft w_upi(i,k) = ud_mf(i,k)*(t(i,k)+epsvt*gdq(i,k,1)) * rair & / (max(cf_upi(i,k),1.e-12)*gdp(i,k)) enddo enddo endif endif !**************************************************************************** KTMAX = 1 do n=1,nctp do i=1,IJSDIM KTMAX(i) = max(KTMAX(i), KT(i,n)) enddo enddo ! do i=1,IJSDIM prec(i) = GTPRP(i,1) snow(i) = GSNWP(i,1) if (prec(i)+snow(i) > 0.0) then kcnv(i) = 1 else kcnv(i) = 0 endif enddo !> -# Multiply mass fluxes by the time step do k=1,kmax do i=1,ijsdim ud_mf(i,k) = ud_mf(i,k) * delta dd_mf(i,k) = dd_mf(i,k) * delta dt_mf(i,k) = dt_mf(i,k) * delta enddo enddo ! rain1(:) = prec(:) * (delta*0.001) ! Convert prec flux kg/m2/sec to rain1 in m do i= 1, IJSDIM rain1(i) = prec(i) * (delta*0.001) enddo ! if (lprnt) then ! write(0,*)' aft cs_cum prec=',prec(ipr),'GTPRP=',GTPRP(ipr,1) ! endif ! if (do_aw) then ! call moist_bud(ijsdim,ijsdim,im,kmax,mype,kdt,grav,delta,delp,prec & ! , gdq(1,1,1), gdq(1,1,2), gdq(1,1,3) & ! , q,clw(1,1,1),clw(1,1,2),'cs_conv_aw') ! endif end subroutine cs_conv_run !> @} !************************************************************************ !>\ingroup cs_scheme !! Main subroutine for the cumulus parameterization with !! state-dependent entrainment rate developed by Minoru Chikira. !! !! - This routine works as the prognostic Arakawa-Schubert scheme !! if OPT_ASMODE is specified. !! - Specify OPT_NS02 to use entrainment rate of Neggers et al. (2002) !! - Specify OPT_CUMBGT to check water and energy budget. !! - Specify OPT_CUMCHK to check range of output values. !! !! History(yy/mm/dd): !! - 08/09/19(chikira) MIROC4.1 !! - 08/10/30(hiro) CMT modified !! - 08/11/11(chikira) Neggers et al. (2002) !! - 08/12/3 (chikira) downdraft detrainment modified !! - 08/12/3 (chikira) COSP output !! - 09/02/24(chikira) fix convective inhibition !! - 09/04/16(hiro) CMIP5 output (cbasep,ctopp) !! - 09/09/03(yokohata) COSP !! - 10/11/19(toshi) small bug fix !! - 14/02/07(chikira) CUMDWN bug fix, CMT modified !!\section gen_cs_cumlus CSAW cs_cumlus General Algorithm !> @{ SUBROUTINE CS_CUMLUS (im , IJSDIM, KMAX , NTR , & !DD dimensions otspt1, otspt2, lprnt , ipr , & GTT , GTQ , GTU , GTV , & ! output CMDET , & ! output GTPRP , GSNWP , GMFX0 , & ! output GMFX1 , CAPE , KT , & ! output CBMFX , & ! modified GDT , GDQ , GDU , GDV , & ! input GDTM , & ! input GDP , GDPM , GDZ , GDZM , & ! input delp , delpi , & DELTA , DELTI , ISTS , IENS, mype,& ! input fscav, fswtr, wcbmaxm, nctp, & ! sigma, vverti, & ! input/output !DDsigma do_aw, do_awdd, flx_form ) ! IMPLICIT NONE Integer, parameter :: ntrq=4 ! starting index for tracers INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, mype, nctp, ipr !! DD, for GFS, pass in logical, intent(in) :: do_aw, do_awdd, flx_form ! switch to apply Arakawa-Wu to the tendencies logical, intent(in) :: otspt1(ntr), otspt2(ntr), lprnt ! ! [OUTPUT] REAL(kind_phys), INTENT(OUT) :: GTT (IJSDIM, KMAX ) ! heating rate REAL(kind_phys), INTENT(OUT) :: GTQ (IJSDIM, KMAX, NTR) ! change in q REAL(kind_phys), INTENT(OUT) :: GTU (IJSDIM, KMAX ) ! tendency of u REAL(kind_phys), INTENT(OUT) :: GTV (IJSDIM, KMAX ) ! tendency of v REAL(kind_phys), INTENT(OUT) :: CMDET (IJSDIM, KMAX ) ! detrainment mass flux ! assuming there is no flux at the top of the atmospherea - Moorthi REAL(kind_phys), INTENT(OUT) :: GTPRP (IJSDIM, KMAX ) ! rain+snow flux REAL(kind_phys), INTENT(OUT) :: GSNWP (IJSDIM, KMAX ) ! snowfall flux REAL(kind_phys), INTENT(OUT) :: GMFX0 (IJSDIM, KMAX ) ! updraft mass flux REAL(kind_phys), INTENT(OUT) :: GMFX1 (IJSDIM, KMAX ) ! downdraft mass flux REAL(kind_phys), INTENT(OUT) :: CAPE (IJSDIM ) INTEGER , INTENT(OUT) :: KT (IJSDIM, NCTP ) ! cloud top ! ! [MODIFIED] REAL(kind_phys), INTENT(INOUT) :: CBMFX (IM, NCTP) ! cloud base mass flux !DDsigma - output added for AW sigma diagnostics ! sigma and vert. velocity as a function of cloud type (1==sfc) real(kind_phys), intent(out), dimension(IM,KMAX) :: sigma !sigma totaled over cloud type - on interfaces (1=sfc) real(kind_phys), intent(out), dimension(IM,KMAX,nctp) :: vverti ! for computing AW flux form of tendencies ! The tendencies are summed over all cloud types ! real(kind_phys), intent(out), dimension(IM,KMAX) :: & !DDsigmadiag real(kind_phys), allocatable, dimension(:,:) :: sfluxterm, qvfluxterm,& ! tendencies of DSE and water vapor due to eddy mass flux qlfluxterm, qifluxterm,& ! tendencies of cloud water and cloud ice due to eddy mass flux ! The fluxes are for an individual cloud type and reused. ! condtermt, condtermq are eddy flux of temperature and water vapor condtermt, condtermq, frzterm, & prectermq, prectermfrz real(kind_phys), allocatable, dimension(:,:,:) :: trfluxterm ! tendencies of tracers due to eddy mass flux ! ! [INPUT] REAL(kind_phys), INTENT(IN) :: GDT (IJSDIM, KMAX ) ! temperature T REAL(kind_phys), INTENT(IN) :: GDQ (IJSDIM, KMAX, NTR) ! humidity, tracer !DDsigmadiag REAL(kind_phys), INTENT(IN) :: GDU (IJSDIM, KMAX ) ! westerly u REAL(kind_phys), INTENT(IN) :: GDV (IJSDIM, KMAX ) ! southern wind v REAL(kind_phys), INTENT(IN) :: GDTM (IJSDIM, KMAX+1 ) ! temperature T REAL(kind_phys), INTENT(IN) :: GDP (IJSDIM, KMAX ) ! pressure P REAL(kind_phys), INTENT(IN) :: GDPM (IJSDIM, KMAX+1 ) ! pressure (half lev) REAL(kind_phys), INTENT(IN) :: GDZ (IJSDIM, KMAX ) ! altitude REAL(kind_phys), INTENT(IN) :: GDZM (IJSDIM, KMAX+1 ) ! altitude REAL(kind_phys), INTENT(IN) :: DELTA ! delta(t) (dynamics) REAL(kind_phys), INTENT(IN) :: DELTI ! delta(t) (internal variable) INTEGER, INTENT(IN) :: ISTS, IENS ! array range real(kind_phys), intent(in) :: fscav(ntr), fswtr(ntr), wcbmaxm(ijsdim) ! ! [INTERNAL WORK] REAL(kind_phys), allocatable :: GPRCC (:, :) ! rainfall ! REAL(kind_phys) GPRCC (IJSDIM, NTR) ! rainfall ! REAL(kind_phys) GSNWC (IJSDIM) ! snowfall ! REAL(kind_phys) CUMCLW(IJSDIM, KMAX) ! cloud water in cumulus ! REAL(kind_phys) CUMFRC(IJSDIM) ! cumulus cloud fraction ! ! REAL(kind_phys) GTCFRC(IJSDIM, KMAX) ! change in cloud fraction ! REAL(kind_phys) FLIQC (IJSDIM, KMAX) ! liquid ratio in cumulus ! ! REAL(kind_phys) GDCFRC(IJSDIM, KMAX) ! cloud fraction ! REAL(kind_phys) GDW (IJSDIM, KMAX) ! total water REAL(kind_phys) DELP (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) REAL(kind_phys) GDQS (IJSDIM, KMAX) ! saturate moisture REAL(kind_phys) FDQS (IJSDIM, KMAX) REAL(kind_phys) GAM (IJSDIM, KMAX) REAL(kind_phys) GDS (IJSDIM, KMAX) ! dry static energy REAL(kind_phys) GDH (IJSDIM, KMAX) ! moist static energy REAL(kind_phys) GDHS (IJSDIM, KMAX) ! saturate MSE ! REAL(kind_phys) GCYM (IJSDIM, KMAX, NCTP)! norm. mass flux (half lev) REAL(kind_phys) GCHB (IJSDIM) ! cloud base MSE-Li*Qi REAL(kind_phys) GCWB (IJSDIM) ! cloud base total water REAL(kind_phys) GCUB (IJSDIM) ! cloud base U REAL(kind_phys) GCVB (IJSDIM) ! cloud base V REAL(kind_phys) GCIB (IJSDIM) ! cloud base ice REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr) ! cloud base tracer REAL(kind_phys) GCYT (IJSDIM, NCTP) ! norm. mass flux @top REAL(kind_phys) GCHT (IJSDIM, NCTP) ! cloud top MSE REAL(kind_phys) GCQT (IJSDIM, NCTP) ! cloud top q REAL(kind_phys) GCwT (IJSDIM) ! cloud top total water REAL(kind_phys) GCUT (IJSDIM, NCTP) ! cloud top U REAL(kind_phys) GCVT (IJSDIM, NCTP) ! cloud top V REAL(kind_phys) GCLT (IJSDIM, NCTP) ! cloud top cloud water REAL(kind_phys) GCIT (IJSDIM, NCTP) ! cloud top cloud ice REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP) ! cloud top tracer REAL(kind_phys) GTPRT (IJSDIM, NCTP) ! precipitation/M REAL(kind_phys) GCLZ (IJSDIM, KMAX) ! cloud liquid for each CTP REAL(kind_phys) GCIZ (IJSDIM, KMAX) ! cloud ice for each CTP ! REAL(kind_phys) ACWF (IJSDIM, NCTP) ! cloud work function REAL(kind_phys) ACWF (IJSDIM ) ! cloud work function REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) ! precipitation REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) ! snowfall REAL(kind_phys) GTPRC0(IJSDIM) ! precip. before evap. REAL(kind_phys) GMFLX (IJSDIM, KMAX) ! mass flux (updraft+downdraft) REAL(kind_phys) QLIQ (IJSDIM, KMAX) ! total cloud liquid REAL(kind_phys) QICE (IJSDIM, KMAX) ! total cloud ice REAL(kind_phys) GPRCI (IJSDIM, KMAX) ! rainfall generation REAL(kind_phys) GSNWI (IJSDIM, KMAX) ! snowfall generation REAL(kind_phys) GPRCP (IJSDIM, KMAX) ! rainfall flux ! REAL(kind_phys) GTEVP (IJSDIM, KMAX) ! evaporation+sublimation REAL(kind_phys) GMDD (IJSDIM, KMAX) ! downdraft mass flux ! REAL(kind_phys) CUMHGT(IJSDIM, NCTP) ! cloud top height ! REAL(kind_phys) CTOPP (IJSDIM) ! cloud top pressure REAL(kind_phys) GDZTR (IJSDIM) ! tropopause height ! REAL(kind_phys) FLIQOU(IJSDIM, KMAX) ! liquid ratio in cumulus INTEGER KB (IJSDIM) INTEGER KSTRT (IJSDIM) ! tropopause level REAL(kind_phys) GAMX REAL(kind_phys) CIN (IJSDIM) INTEGER JBUOY (IJSDIM) REAL(kind_phys) DELZ, BUOY, DELWC, DELER REAL(kind_phys) WCBX (IJSDIM) ! REAL(kind_phys) ERMR (NCTP) ! entrainment rate (ASMODE) ! SAVE ERMR INTEGER KTMX (NCTP) ! max of cloud top INTEGER KTMXT ! max of cloud top ! REAL(kind_phys) TIMED REAL(kind_phys) GDCLDX, GDMU2X, GDMU3X ! ! REAL(kind_phys) HBGT (IJSDIM) ! imbalance in column heat ! REAL(kind_phys) WBGT (IJSDIM) ! imbalance in column water !DDsigma begin local work variables - all on model interfaces (sfc=1) REAL(kind_phys) lamdai ! lamda for cloud type ctp REAL(kind_phys) gdqm, gdlm, gdim ! water vapor REAL(kind_phys) gdtrm(ntrq:ntr) ! tracer ! the following are new arguments to cumup to get them out for AW REAL(kind_phys) wcv (IJSDIM, KMAX) ! in-cloud vertical velocity REAL(kind_phys) GCTM (IJSDIM, KMAX) ! cloud T (half lev) !DDsigmadiag make output REAL(kind_phys) GCQM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCwM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCiM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GClM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GChM (IJSDIM, KMAX) ! cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCtrM (IJSDIM, KMAX, ntrq:ntr) ! cloud tracer (half lev) !DDsigmadiag make output ! eddy flux profiles for dse, water vapor, cloud water, cloud ice REAL(kind_phys), dimension(Kmax+1) :: sfluxtem, qvfluxtem, qlfluxtem, qifluxtem REAL(kind_phys), dimension(Kmax+1,ntrq:ntr) :: trfluxtem ! tracer ! tendency profiles - condensation heating, condensation moistening, heating due to ! freezing, total precip production, frozen precip production REAL(kind_phys), dimension(ijsdim,Kmax) :: dtcondtem, dqcondtem, dtfrztem, dqprectem,& ! Moorthi dfrzprectem, lamdaprod !< product of (1+lamda) through cloud type ctp REAL(kind_phys), dimension(ijsdim,Kmax) :: dtevap, dqevap, dtmelt, dtsubl ! factor to modify precip rate to force conservation of water. With bug fixes it's ! not doing anything now. REAL(kind_phys), dimension(ijsdim) :: moistening_aw real(kind_phys), dimension(ijsdim,kmax) :: gctbl, gcqbl,gcwbl, gcqlbl, gcqibl, & !DDsigmadiag updraft profiles below cloud Base sigmad ! downdraft area fraction real(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: gctrbl !DDsigmadiag tracer updraft profiles below cloud Base ! rhs_q, rhs_h are residuals of condensed water, MSE budgets to compute condensation, ! and heating due to freezing real(kind_phys) :: rhs_q, rhs_h, fsigma, sigmai, delpinv ! real(kind_phys) :: rhs_q, rhs_h, sftem, qftem, qlftem, qiftem, & ! fsigma ! factor to reduce mass flux terms (1-sigma**2) for AW !DDsigma end local work variables ! ! profiles of heating due to precip evaporation, melting and sublimation, and the ! evap, melting and sublimation rates. REAL(kind_phys), allocatable, dimension(:,:) :: dtdwn, & ! t tendency downdraft detrainment dqvdwn, & ! qv tendency downdraft detrainment dqldwn, & ! ql tendency downdraft detrainment dqidwn ! qi tendency downdraft detrainment REAL(kind_phys), allocatable, dimension(:,:,:) :: dtrdwn ! tracer tendency downdraft detrainment !DDsigma end local work variables ! ! [INTERNAL PARM] REAL(kind_phys), parameter :: WCBMIN = zero ! min. of updraft velocity at cloud base !M REAL(kind_phys) :: WCBMAX = 1.4_kind_phys ! max. of updraft velocity at cloud base !M wcbas commented by Moorthi since it is not used !M REAL(kind_phys) :: WCBAS = 2._kind_phys ! updraft velocity**2 at cloud base (ASMODE) !M REAL(kind_phys) :: ERAMIN = 1.e-5_kind_phys ! min. of entrainment rate ! used only in OPT_ASMODE !M REAL(kind_phys) :: ERAMAX = 2.e-3_kind_phys ! max. of entrainment rate ! used only in OPT_ASMODE LOGICAL :: OINICB = .false. ! set 0.d0 to CBMFX when .true. ! REAL(kind_phys) :: VARMIN = 1.e-13_kind_phys ! minimum of PDF variance ! REAL(kind_phys) :: VARMAX = 5.e-7_kind_phys ! maximum of PDF variance ! REAL(kind_phys) :: SKWMAX = 0.566_kind_phys ! maximum of PDF skewness REAL(kind_phys) :: PSTRMX = 400.e2_kind_phys ! max P of tropopause REAL(kind_phys) :: PSTRMN = 50.e2_kind_phys ! min P of tropopause REAL(kind_phys) :: GCRSTR = 1.e-4_kind_phys ! crit. dT/dz tropopause real(kind=kind_phys) :: tem, esat, mflx_e, cbmfl, tem1, tem2, tem3 INTEGER :: KBMX, I, K, CTP, ierr, n, kp1, km1, kk, kbi, l, l1 ! LOGICAL, SAVE :: OFIRST = .TRUE. ! called first time? ! IF (OFIRST) THEN OFIRST = .FALSE. IF (OINICB) THEN CBMFX = zero ENDIF ENDIF ! do n=1,ntr do k=1,kmax do i=1,ijsdim gtq(i,k,n) = zero enddo enddo enddo do k=1,kmax do i=1,ijsdim gtt(i,k) = zero gtu(i,k) = zero gtv(i,k) = zero gmflx(i,k) = zero gmfx0(i,k) = zero gprci(i,k) = zero gsnwi(i,k) = zero qliq(i,k) = zero qice(i,k) = zero ! gtcfrc(i,k) = zero ! cumclw(i,k) = zero ! fliqc(i,k) = zero sigma(i,k) = zero enddo enddo if (flx_form) then allocate(sfluxterm(ijsdim,kmax), qvfluxterm(ijsdim,kmax), qlfluxterm(ijsdim,kmax), & qifluxterm(ijsdim,kmax), condtermt(ijsdim,kmax), condtermq(ijsdim,kmax), & frzterm(ijsdim,kmax), prectermq(ijsdim,kmax), prectermfrz(ijsdim,kmax), & dtdwn(ijsdim,kmax), dqvdwn(ijsdim,kmax), dqldwn(ijsdim,kmax), & dqidwn(ijsdim,kmax), trfluxterm(ijsdim,kmax,ntrq:ntr), & dtrdwn(ijsdim,kmax,ntrq:ntr)) do k=1,kmax do i=1,ijsdim sfluxterm(i,k) = zero qvfluxterm(i,k) = zero qlfluxterm(i,k) = zero qifluxterm(i,k) = zero condtermt(i,k) = zero condtermq(i,k) = zero frzterm(i,k) = zero prectermq(i,k) = zero prectermfrz(i,k) = zero dtdwn(i,k) = zero dqvdwn(i,k) = zero dqldwn(i,k) = zero dqidwn(i,k) = zero cmdet(i,k) = zero enddo enddo do n = ntrq,ntr do k=1,kmax do i=1,ijsdim trfluxterm(i,k,n) = zero dtrdwn(i,k,n) = zero enddo enddo enddo endif do i=1,ijsdim ! gprcc(i,:) = zero gtprc0(i) = zero ! hbgt(i) = zero ! wbgt(i) = zero gdztr(i) = zero kstrt(i) = kmax enddo do k=1,kmax do i=1,ijsdim GDW(i,k) = GDQ(i,k,1) + GDQ(i,k,ITL) + GDQ(i,k,iti) enddo enddo !> -# Compute layer saturate moisture \f$Q_i\f$(GDQS) and !! saturate moist static energy (GDHS; see Appendix B in !! Chikira and Sugiyama (2010) \cite Chikira_2010) DO K=1,KMAX DO I=ISTS,IENS esat = min(gdp(i,k), fpvs(gdt(i,k))) GDQS(I,K) = min(EPSV*esat/max(gdp(i,k)+epsvm1*esat, 1.0e-10), 0.1) tem = one / GDT(I,K) FDQS(I,K) = GDQS(I,K) * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT GAM (I,K) = ELOCP*FDQS(I,K) GDS (I,K) = CP*GDT(I,K) + GRAV*GDZ(I,K) ! layer dry static energy GDH (I,K) = GDS(I,K) + EL*GDQ(I,K,1) ! layer moist static energy GDHS(I,K) = GDS(I,K) + EL*GDQS(I,K) ! layer sat. moist static energy ENDDO ENDDO ! ! < tropopause > ! !> -# Compute tropopause height (GDZTR) DO K=1,KMAX kp1 = k + 1 DO I=ISTS,IENS GAMX = (GDTM(I,KP1)-GDTM(I,K)) / (GDZM(I,KP1)-GDZM(I,K)) IF ((GDP(I,K) < PSTRMX .AND. GAMX > GCRSTR) .OR. GDP(I,K) < PSTRMN) THEN KSTRT(I) = MIN(K, KSTRT(I)) ENDIF ENDDO ENDDO DO I=ISTS,IENS K = KSTRT(I) GDZTR(I) = GDZM(I,K) ENDDO ! !DDsigma - arguments added to get subcloud profiles in updraft ! so AW eddy flux tendencies can be computed !> -# Call cumbas() to compute cloud base properties CALL CUMBAS(IJSDIM, KMAX , & !DD dimensions KB , GCYM(1,1,1) , KBMX , & ! output ntr , ntrq , & GCHB , GCWB , GCUB , GCVB , & ! output GCIB , gctrb, & ! output GDH , GDW , GDHS , GDQS , & ! input GDQ(1,1,iti) , GDU , GDV , GDZM , & ! input GDPM , FDQS , GAM , & ! input lprnt, ipr, & ISTS , IENS , & !) ! input gctbl, gcqbl,gdq,gcwbl, gcqlbl, gcqibl, gctrbl) ! sub cloud tendencies ! !> -# Compute CAPE and CIN ! DO I=ISTS,IENS CAPE(i) = zero CIN(i) = zero JBUOY(i) = 0 enddo DO K=2,KMAX DO I=ISTS,IENS if (kb(i) > 0) then IF (K >= KB(I)) THEN BUOY = (GDH(I,1)-GDHS(I,K)) / ((one+ELOCP*FDQS(I,K)) * CP*GDT(I,K)) ELSE BUOY = (GDS(I,1)-GDS(I,K)) / (CP*GDT(I,K)) END IF IF (BUOY > zero .AND. JBUOY(I) >= -1) THEN CAPE(I) = CAPE(I) + BUOY * GRAV * (GDZM(I,K+1) - GDZM(I,K)) JBUOY(I) = 2 ELSEIF (BUOY < zero .AND. JBUOY(I) /= 2) THEN CIN(I) = CIN(I) - BUOY * GRAV * (GDZM(I,K+1) - GDZM(I,K)) JBUOY(I) = -1 ENDIF endif ENDDO ENDDO DO I=ISTS,IENS IF (JBUOY(I) /= 2) CIN(I) = -999.D0 if (cin(i) < cincrit) kb(i) = -1 ENDDO !DDsigma some initialization before summing over cloud type !> -# Initialize variables before summing over cloud types do k=1,kmax ! Moorthi do i=1,ijsdim lamdaprod(i,k) = one enddo enddo do ctp=2,nctp do k=1,kmax do i=1,ijsdim gcym(i,k,ctp) = gcym(i,k,1) enddo enddo enddo DO CTP=1,NCTP ! loop over cloud types tem = ctp / DBLE(NCTP) do i=1,ijsdim DELWC = tem * (WCBMAXm(i) - WCBMIN) WCBX(I) = DELWC * DELWC enddo do k=1,kmax ! Moorthi do i=1,ijsdim dqcondtem(i,k) = zero dqprectem(i,k) = zero dfrzprectem(i,k) = zero dtfrztem(i,k) = zero dtcondtem(i,k) = zero enddo enddo ! getting more incloud profiles of variables to compute eddy flux tendencies ! and condensation rates ! DH* GNU crashes - check all arguments to CUMUP for their dimensions ! before and after CUMUP (i.e. here), and inside the routine, in ! particular: gctm, gcqm, gcwm, gchm, gcwt, gclm, gcim,gctrm ! also, inside, check that no reads/writes out of bounds occur *DH !> -# Call cumup() to compute in-cloud properties CALL CUMUP(IJSDIM, KMAX, NTR, ntrq, & !DD dimensions ACWF , & ! output GCLZ , GCIZ , GPRCIZ , GSNWIZ, & ! output GCYT(1,CTP) , GCHT(1,CTP) , GCQT (1,CTP), & ! output GCLT(1,CTP) , GCIT(1,CTP) , GTPRT(1,CTP), & ! output GCUT(1,CTP) , GCVT(1,CTP) , gctrt(1,ntrq:ntr,ctp), & ! output KT (1,CTP) , KTMX(CTP) , & ! output GCYM(1,1,CTP) , & ! modified wcv , & ! !DD-sigma new output GCHB , GCWB , GCUB , GCVB , & ! input !DDsigmadiag GCIB , gctrb , & ! input GDU , GDV , GDH , GDW , & ! input GDHS , GDQS , GDT , GDTM , & ! input GDQ , GDQ(1,1,iti) , GDZ , GDZM , & ! input GDPM , FDQS , GAM , GDZTR , & ! input CPRES , WCBX , & ! input KB , CTP , ISTS , IENS , & ! input gctm , gcqm, gcwm, gchm, gcwt, gclm, gcim, gctrm, & ! additional incloud profiles and cloud top total water lprnt , ipr ) ! !> -# Call cumbmx() to compute cloud base mass flux CALL CUMBMX(IJSDIM, KMAX, & !DD dimensions CBMFX(1,CTP), & ! modified ACWF , GCYT(1,CTP), GDZM , & ! input GDW , GDQS , DELP , & ! input KT (1,CTP), KTMX(CTP) , KB , & ! input DELTI , ISTS , IENS ) !DDsigma - begin sigma computation ! At this point cbmfx is updated and we have everything we need to compute sigma do i=ISTS,IENS if (flx_form) then !> -# Initialize eddy fluxes for cloud types do k=1,kmax+1 sfluxtem(k) = zero qvfluxtem(k) = zero qlfluxtem(k) = zero qifluxtem(k) = zero enddo do n=ntrq,ntr ! tracers do k=1,kmax+1 trfluxtem(k,n) = zero enddo enddo endif cbmfl = cbmfx(i,ctp) kk = kt(i,ctp) ! cloud top index if(cbmfl > zero) then ! this should avoid zero wcv in the denominator kbi = kb(i) ! cloud base index do k=kbi,kk ! loop from cloud base to cloud top km1 = k - 1 rhs_h = zero rhs_q = zero !> -# Interpolate environment variables to layer interface GDQM = half * (GDQ(I,K,1) + GDQ(I,KM1,1)) ! as computed in cumup ! GDwM = half * (GDw(I,K) + GDw(I,KM1 )) GDlM = half * (GDQ(I,K,3) + GDQ(I,KM1,3)) GDiM = half * (GDQ(I,K,2) + GDQ(I,KM1,2)) do n = ntrq,NTR GDtrM(n) = half * (GDQ(I,K,n) + GDQ(I,KM1,n)) ! as computed in cumup enddo mflx_e = gcym(i,k,ctp) * cbmfl ! mass flux at level k for cloud ctp if (do_aw) then !> -# Compute lamda for a cloud type and then updraft area fraction !! (sigmai) following Equations 23 and 12 of !! Arakawa and Wu (2013) \cite arakawa_and_wu_2013 , respectively lamdai = mflx_e * rair * gdtm(i,k)*(one+epsvt*gdqm) & / (gdpm(i,k)*wcv(i,k)) lamdaprod(i,k) = lamdaprod(i,k) * (one+lamdai) ! vverti(i,k,ctp) = wcv(i,k) ! sigmai(i,k,ctp) = lamdai / lamdaprod(i,k) ! sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai(i,k,ctp))) sigmai = lamdai / lamdaprod(i,k) sigma(i,k) = max(zero, min(one, sigma(i,k) + sigmai)) vverti(i,k,ctp) = sigmai * wcv(i,k) else sigma(i,k) = 0.0 endif if (flx_form) then ! fsigma = 1.0 ! no aw effect, comment following lines to undo AW fsigma = one - sigma(i,k) !> -# Compute tendencies based on mass flux and condensation ! fsigma is the AW reduction of flux tendencies if(k == kbi) then do l=2,kbi ! compute eddy fluxes below cloud base tem = - fsigma * gcym(i,l,ctp) * cbmfl ! first get environment variables at layer interface l1 = l - 1 GDQM = half * (GDQ(I,l,1) + GDQ(I,l1,1)) GDlM = half * (GDQ(I,l,3) + GDQ(I,l1,3)) GDiM = half * (GDQ(I,l,2) + GDQ(I,l1,2)) !! GDwM = half * (GDw(I,l) + GDw(I,l1)) do n = ntrq,NTR GDtrM(n) = half * (GDQ(I,l,n) + GDQ(I,l1,n)) ! as computed in cumup enddo ! flux = mass flux * (updraft variable minus environment variable) !centered differences sfluxtem(l) = tem * (gdtm(i,l)-gctbl(i,l)) qvfluxtem(l) = tem * (gdqm-gcqbl(i,l)) qlfluxtem(l) = tem * (gdlm-gcqlbl(i,l)) qifluxtem(l) = tem * (gdim-gcqibl(i,l)) do n = ntrq,NTR trfluxtem(l,n) = tem * (gdtrm(n)-gctrbl(i,l,n)) enddo ! if(lprnt .and. i == ipr) write(0,*)' l=',l,' kbi=',kbi,' tem =', tem,' trfluxtem=',trfluxtem(l,ntr),& ! ' gdtrm=',gdtrm(ntr),' gctrbl=',gctrbl(i,l,ntr),' gq=',GDQ(I,l,ntr),GDQ(I,l1,ntr),' l1=',l1,' ctp=',ctp,& ! ' fsigma=',fsigma,' gcym=',gcym(i,l,ctp),' cbmfl=',cbmfl,' sigma=',sigma(i,k) ! The following commented out by Moorthi on April 13, 2018 because tke below ! cloud base becomes too large otherwise when shoc is used !upstream - This better matches what the original CS tendencies do ! sfluxtem(l) = tem * (gdt(i,l)+gocp*(gdz(i,l)-gdzm(i,l))-gctbl(i,l)) ! qvfluxtem(l) = tem * (gdq(i,l,1)-gcqbl(i,l)) ! qlfluxtem(l) = tem * (gdq(i,l,3)-gcqlbl(i,l)) ! qifluxtem(l) = tem * (gdq(i,l,2)-gcqibl(i,l)) ! do n = ntrq,NTR ! trfluxtem(l,n) = tem * (gdq(i,l,n)-gctrbl(i,l,n)) ! enddo enddo else ! flux = mass flux * (updraft variable minus environment variable) tem = - fsigma * mflx_e !centered sfluxtem(k) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k)) qvfluxtem(k) = tem * (gdqm-gcqm(i,k)) qlfluxtem(k) = tem * (gdlm-gclm(i,k)) qifluxtem(k) = tem * (gdim-gcim(i,k)) do n = ntrq,NTR trfluxtem(k,n) = tem * (gdtrm(n)-gctrm(i,k,n)) enddo !upstream - This better matches what the original CS tendencies do ! if(k < kk) then ! sfluxtem(k) = tem * (gdt(i,k)+gocp*gdz(i,k)-gctm(i,k)) ! qvfluxtem(k) = tem * (gdq(i,k,1)-gcqm(i,k)) ! qlfluxtem(k) = tem * (gdq(i,k,3)-gclm(i,k)) ! qifluxtem(k) = tem * (gdq(i,k,2)-gcim(i,k)) ! do n = ntrq,NTR ! trfluxtem(k,n) = tem * (gdq(i,k,n)-gctrm(i,k,n)) ! enddo ! if(lprnt .and. i == ipr) write(0,*)' k=',k,' kbi=',kbi,' tem =', tem,' kk=',kk,& ! ' gctrm=',gctrm(i,k,ntr),' gdq=',gdq(I,k,ntr),' gctrm=',gctrm(I,k,ntr),' ctp=',ctp,& ! ' fsigma=',fsigma,' mflx_e=',mflx_e,' trfluxtemk=',trfluxtem(k,ntr),' sigma=',sigma(i,k) ! else ! centered at top of cloud ! sfluxtem(k) = tem * (gdtm(i,k)+gocp*gdzm(i,k)-gctm(i,k)) ! qvfluxtem(k) = tem * (gdqm-gcqm(i,k)) ! qlfluxtem(k) = tem * (gdlm-gclm(i,k)) ! qifluxtem(k) = tem * (gdim-gcim(i,k)) ! do n = ntrq,NTR ! trfluxtem(k,n) = tem * (gdtrm(n)-gctrm(i,k,n)) ! enddo ! endif ! if(lprnt .and. i == ipr) write(0,*)' k=',k,' kbi=',kbi,' tem =', tem,' kk=',kk,& ! ' gctrm=',gctrm(i,k,ntr),' gdtrm=',gdtrm(ntr),' gctrm=',gctrm(I,k,ntr),' ctp=',ctp,& ! ' fsigma=',fsigma,' mflx_e=',mflx_e,' trfluxtemk=',trfluxtem(k,ntr),' sigma=',sigma(i,k) ! the condensation terms - these come from the MSE and condensed water budgets for ! an entraining updraft ! if(k > kb(i)) then ! comment for test ! if(k <= kk) then ! Moorthi ! if(k < kt(i,ctp)) then ! rhs_h = cbmfl*(gcym(i,k)*gchm(i,k) - (gcym(i,km1)*gchm(i,km1) & ! + GDH(I,Km1 )*(gcym(i,k)-gcym(i,km1))) ) ! rhs_q = cbmfl*(gcym(i,k)*(gcwm(i,k)-gcqm(i,k)) & ! - (gcym(i,km1)*(gcwm(i,km1)-gcqm(i,km1)) & ! + (GDw( I,Km1 )-gdq(i,km1,1))*(gcym(i,k)-gcym(i,km1))) ) ! tem = cbmfl * (one - sigma(i,k)) tem = cbmfl * (one - 0.5*(sigma(i,k)+sigma(i,km1))) tem1 = gcym(i,k,ctp) * (one - sigma(i,k)) tem2 = gcym(i,km1,ctp) * (one - sigma(i,km1)) rhs_h = cbmfl * (tem1*gchm(i,k) - (tem2*gchm(i,km1) & + GDH(I,Km1)*(tem1-tem2)) ) rhs_q = cbmfl * (tem1*(gcwm(i,k)-gcqm(i,k)) & - (tem2*(gcwm(i,km1)-gcqm(i,km1)) & + (GDw(I,Km1)-gdq(i,km1,1))*(tem1-tem2)) ) ! ELSE ! rhs_h = cbmfl*(gcht(i,ctp) - (gcym(i,k-1)*gchm(i,k-1) + GDH( I,K-1 )*(gcyt(i,ctp)-gcym(i,k-1))) ) ! rhs_q = cbmfl*((gcwt(i)-gcqt(i,ctp)) - (gcym(i,k-1)*(gcwm(i,k-1)-gcqm(i,k-1)) + (GDw( I,K-1 )-gdq(i,k-1,1))*(gcyt(i,ctp)-gcym(i,k-1))) ) ! endif !> -# Compute condensation, total precipitation production, frozen precipitation production, !! heating due to freezing, and total temperature tendency due to in-cloud microphysics dqcondtem(i,km1) = -rhs_q ! condensation ! dqprectem(i,km1) = cbmfl * (GPRCIZ(i,k) + GSNWIZ(i,k)) dqprectem(i,km1) = tem * (GPRCIZ(i,k) + GSNWIZ(i,k)) ! total precip production ! dfrzprectem(i,km1) = cbmfl * GSNWIZ(i,k) dfrzprectem(i,km1) = tem * GSNWIZ(i,k) ! production of frozen precip dtfrztem(i,km1) = rhs_h*oneocp ! heating due to freezing dtcondtem(i,km1) = - elocp * dqcondtem(i,km1) + dtfrztem(i,km1) endif ! if(k > kbi) then endif ! if (flx_form) enddo ! end of k=kbi,kk loop endif ! end of if(cbmfl > zero) ! get tendencies by difference of fluxes, sum over cloud type if (flx_form) then do k = 1,kk delpinv = delpi(i,k) !> -# Sum single cloud microphysical tendencies over all cloud types condtermt(i,k) = condtermt(i,k) + dtcondtem(i,k) * delpinv condtermq(i,k) = condtermq(i,k) + dqcondtem(i,k) * delpinv prectermq(i,k) = prectermq(i,k) + dqprectem(i,k) * delpinv prectermfrz(i,k) = prectermfrz(i,k) + dfrzprectem(i,k) * delpinv frzterm(i,k) = frzterm(i,k) + dtfrztem(i,k) * delpinv !> -# Compute flux tendencies and vertical flux divergence sfluxterm(i,k) = sfluxterm(i,k) - (sfluxtem(k+1) - sfluxtem(k)) * delpinv qvfluxterm(i,k) = qvfluxterm(i,k) - (qvfluxtem(k+1) - qvfluxtem(k)) * delpinv qlfluxterm(i,k) = qlfluxterm(i,k) - (qlfluxtem(k+1) - qlfluxtem(k)) * delpinv qifluxterm(i,k) = qifluxterm(i,k) - (qifluxtem(k+1) - qifluxtem(k)) * delpinv do n = ntrq,ntr trfluxterm(i,k,n) = trfluxterm(i,k,n) - (trfluxtem(k+1,n) - trfluxtem(k,n)) * delpinv enddo ! if (lprnt .and. i == ipr) write(0,*)' k=',k,' trfluxtem=',trfluxtem(k+1,ntr),trfluxtem(k,ntr),& ! ' ctp=',ctp,' trfluxterm=',trfluxterm(i,k,ntr) enddo endif ! if (flx_form) enddo ! end of i loop ! do i=ists,iens if (cbmfx(i,ctp) > zero) then tem = one - sigma(i,kt(i,ctp)) gcyt(i,ctp) = tem * gcyt(i,ctp) gtprt(i,ctp) = tem * gtprt(i,ctp) gclt(i,ctp) = tem * gclt(i,ctp) gcht(i,ctp) = tem * gcht(i,ctp) gcqt(i,ctp) = tem * gcqt(i,ctp) gcit(i,ctp) = tem * gcit(i,ctp) if (.not. flx_form) then do n = ntrq,ntr gctrt(i,n,ctp) = tem * gctrt(i,n,ctp) enddo end if gcut(i,ctp) = tem * gcut(i,ctp) gcvt(i,ctp) = tem * gcvt(i,ctp) do k=1,kmax kk = kb(i) if (k < kk) then tem = one - sigma(i,kk) tem1 = tem else tem = one - sigma(i,k) tem1 = one - 0.5*(sigma(i,k)+sigma(i,k-1)) endif gcym(i,k,ctp) = tem * gcym(i,k,ctp) gprciz(i,k) = tem1 * gprciz(i,k) gsnwiz(i,k) = tem1 * gsnwiz(i,k) gclz(i,k) = tem1 * gclz(i,k) gciz(i,k) = tem1 * gciz(i,k) enddo endif enddo ! !> -# Call cumflx() to compute cloud mass flux and precipitation CALL CUMFLX(IM , IJSDIM, KMAX , & !DD dimensions GMFX0 , GPRCI , GSNWI , CMDET, & ! output QLIQ , QICE , GTPRC0, & ! output CBMFX(1,CTP) , GCYM(1,1,ctp), GPRCIZ , GSNWIZ , & ! input GTPRT(1,CTP) , GCLZ , GCIZ , GCYT(1,ctp),& ! input KB , KT(1,CTP) , KTMX(CTP) , & ! input ISTS , IENS ) ! input ENDDO ! end of cloud type ctp loop !> -# Compute net updraft mass flux for all clouds do k=1,kmax do i=ists,iens GMFLX(I,k) = GMFX0(I,k) ! contains net updraft mass flux for all clouds enddo enddo KTMXT = 3 DO CTP=1,NCTP IF (KTMX(CTP) > KTMXT) KTMXT = KTMX(CTP) ENDDO ! DO K=1,KTMXT ! DO I=ISTS,IENS ! CUMCLW(I,K) = QLIQ(I,K) + QICE(I,K) ! IF (CUMCLW(I,K) > zero) THEN ! FLIQC(I,K) = QLIQ(I,K) / CUMCLW(I,K) ! ENDIF ! ENDDO ! ENDDO ! ! Cumulus Cloudiness ! CALL CUMCLD(IJSDIM, KMAX , & !DD dimensions ! CUMCLW, QLIQ , QICE , FLIQC , & ! modified ! CUMFRC, & ! output ! GMFLX , KTMXT , ISTS , IENS ) ! input ! ! - Call cumdet() to compute cloud detrainment heating. if (.not. flx_form) then CALL CUMDET(im , IJSDIM, KMAX , NTR , ntrq , & !DD dimensions GTT , GTQ , GTU , GTV , & ! modified GDH , GDQ , GDU , GDV , & ! input ! GTT , GTQ , GTCFRC, GTU , GTV , & ! modified ! GDH , GDQ , GDCFRC, GDU , GDV , & ! input CBMFX , GCYT , DELPI , GCHT , GCQT , & ! input GCLT , GCIT , GCUT , GCVT , GDQ(1,1,iti),& ! input gctrt , & KT , ISTS , IENS, nctp ) ! input endif !for now area fraction of the downdraft is zero, it will be computed ! within cumdwn and applied there ! Get AW downdraft eddy flux and microphysical tendencies out of downdraft code. do k=1,kmax do i=ists,iens sigmad(i,k) = zero enddo enddo !> -# Call cumdwn() to compute cumulus downdraft and assocated melt, freeze !! and evaporation CALL CUMDWN(IM , IJSDIM, KMAX , NTR , ntrq , & ! DD dimensions GTT , GTQ , GTU , GTV , & ! modified GMFLX , & ! modified updraft+downdraft flux GPRCP , GSNWP , GTEVP , GMDD , & ! output GPRCI , GSNWI , & ! input GDH , GDW , GDQ , GDQ(1,1,iti) , & ! input GDQS , GDS , GDHS , GDT , & ! input GDU , GDV , GDZ , & ! input GDZM , FDQS , DELP , DELPI , & ! input sigmad, do_aw , do_awdd, flx_form, & ! DDsigma input dtmelt, dtevap, dtsubl, & ! DDsigma input dtdwn , dqvdwn, dqldwn, dqidwn, & ! DDsigma input dtrdwn, & KB , KTMXT , ISTS , IENS ) ! input ! here we substitute the AW tendencies into tendencies to be passed out ! if (do_aw) then ! do k=1,kmax ! do i=ists,iens ! sigma(i,k) = sigma(i,k) + sigmad(i,k) ! enddo ! enddo !> -# Call cumsbw() to compute cloud subsidence heating if (.not. flx_form) then ! Cloud Subsidence Heating ! -----------------------= CALL CUMSBH(IM , IJSDIM, KMAX , NTR , ntrq , & !DD dimensions GTT , GTQ , & ! modified GTU , GTV , & ! modified GDH , GDQ , GDQ(1,1,iti) , & ! input GDU , GDV , & ! input DELPI , GMFLX , GMFX0 , & ! input KTMXT , CPRES , kb, ISTS , IENS ) ! input else CALL CUMSBW(IM , IJSDIM, KMAX , & !DD dimensions GTU , GTV , & ! modified GDU , GDV , & ! input DELPI , GMFLX , GMFX0 , & ! input KTMXT , CPRES , kb, ISTS , IENS ) ! input endif ! ! for now the following routines appear to be of no consequence to AW - DD ! if (.not. flx_form) then ! Tracer Updraft properties ! ------------- allocate (gprcc(ijsdim,ntr)) do n=1,ntr do i=1,ijsdim gprcc(i,n) = zero enddo enddo CALL CUMUPR(im , IJSDIM, KMAX , NTR , & !DD dimensions GTQ , GPRCC , & ! modified GDQ , CBMFX , & ! input GCYM , GCYT , GCQT , GCLT , GCIT , & ! input GTPRT , GTEVP , GTPRC0, & ! input KB , KBMX , KT , KTMX , KTMXT , & ! input DELPI , OTSPT1, ISTS , IENS, & ! input fscav , fswtr, nctp) ! ! Tracer Change due to Downdraft ! --------------- CALL CUMDNR(im ,IJSDIM , KMAX , NTR , & !DD dimensions GTQ , & ! modified GDQ , GMDD , DELPI , & ! input KTMXT , OTSPT1, ISTS , IENS ) ! input !! !! Tracer change due to Subsidence !! --------------- !! This will be done by cumsbh, now DD 20170907 ! CALL CUMSBR(im , IJSDIM, KMAX , NTR , & !DD dimensions ! GTQ , & ! modified ! GDQ , DELPI , & ! input ! GMFLX , KTMXT , OTSPT2, & ! input ! ISTS , IENS ) ! input endif ! if this tracer not advected zero it out DO n = ntrq,NTR if (.not. OTSPT2(n)) then DO K=1,KMAX DO I=ISTS,IENS gtq(i,k,n) = 0.0 ENDDO ENDDO endif ENDDO ! if(do_aw .and. flx_form) then ! compute AW tendencies !> -# Compute AW tendencies of T, ql and qi if(flx_form) then ! compute AW tendencies ! AW lump all heating together, compute qv term do k=1,kmax do i=ists,iens dqevap(i,k) = - dtevap(i,k)*cpoel - dtsubl(i,k)*cpoesub dtevap(i,k) = dtevap(i,k) + dtsubl(i,k) dtsubl(i,k) = zero enddo enddo do i=1,ijsdim moistening_aw(i) = zero enddo tem2 = one / delta DO K=1,KMAX DO I=ISTS,IENS tem = frzterm(i,k)*cpoEMELT - prectermfrz(i,k) gtt(i,k) = dtdwn(i,k) + sfluxterm(i,k) + condtermt(i,k) & + dtmelt(i,k) + dtevap(i,k) gtq(i,k,1) = dqvdwn(i,k) + qvfluxterm(i,k) + condtermq(i,k) & + dqevap(i,k) gtq(i,k,itl) = dqldwn(i,k) + qlfluxterm(i,k) - condtermq(i,k) & - prectermq(i,k) - tem gtq(i,k,iti) = dqidwn(i,k) + qifluxterm(i,k) + tem ! detrainment terms get zeroed ! gtldet(i,k) = zero ! gtidet(i,k) = zero tem1 = - gdq(i,k,itl)*tem2 if (gtq(i,k,itl) < tem1) then tem3 = gtq(i,k,itl) - tem1 gtq(i,k,1) = gtq(i,k,1) + tem3 gtq(i,k,itl) = tem1 gtt(i,k) = gtt(i,k) - elocp*tem3 endif tem1 = - gdq(i,k,iti)*tem2 if (gtq(i,k,iti) < tem1) then tem3 = gtq(i,k,iti) - tem1 gtq(i,k,1) = gtq(i,k,1) + tem3 gtq(i,k,iti) = tem1 gtt(i,k) = gtt(i,k) - esubocp*tem3 endif tem1 = - gdq(i,k,1)*tem2 if (gtq(i,k,1) < tem1) then gtt(i,k) = gtt(i,k) + elocp*(gtq(i,k,1)-tem1) gtq(i,k,1) = tem1 endif ! column-integrated total water tendency - to be used to impose water conservation moistening_aw(i) = moistening_aw(i) & + (gtq(i,k,1)+gtq(i,k,itl)+gtq(i,k,iti)) * delp(i,k) * gravi ENDDO ENDDO ! replace tracer tendency only if to be advected. DO n = ntrq,NTR if (OTSPT2(n)) then DO K=1,KMAX DO I=ISTS,IENS gtq(i,k,n) = dtrdwn(i,k,n) + trfluxterm(i,k,n) ENDDO ENDDO endif ENDDO ! if (lprnt) write(0,*)' endcs_cum gtq=',gtq(ipr,1:25,ntr) ! if (lprnt) write(0,*)' endcs_cum trfluxterm=',trfluxterm(ipr,1:25,ntr) endif ! if (flx_form) !!!! this section may need adjustment for cloud ice and water with flux_form ! ! do k=1,kmax ! do i=ISTS,IENS ! GTQ(I,k,ITI) = GTQI(I,k) ! enddo ! enddo ! !> -# Call cumfxr() for tracer mass fixer without detrainment CALL CUMFXR(IM , IJSDIM, KMAX , NTR , & !DD dimensions GTQ , & ! modified GDQ , DELP , DELTA , KTMXT , IMFXR, & ! input ISTS , IENS ) ! input ! ! do k=1,kmax ! do i=ISTS,IENS ! GTQL(I,k) = GTQ(I,k,ITL) + GTLDET(I,k) + GTIDET(I,k) ! enddo ! enddo ! ! Tracer mass fixer with detrainment ! CALL CUMFXR1(IM , IJSDIM, KMAX , & !DD dimensions ! GTQL , & ! modified ! GDQ(1,1,ITL), DELP, DELTA, KTMXT, IMFXR(ITL), & ! input ! ISTS , IENS ) ! input ! !!!!! end fixer section ! DO K=1,KMAX ! DO I=ISTS, IENS ! GTLDET(I,k) = GTQL(I,k) - GTQ(I,k,ITL) - GTIDET(I,k) ! tendencies of subgrid PDF (turned off) ! GDCLDX = GDCFRC( I,K ) + GTCFRC( I,K )*DELTA ! GDCLDX = MIN( MAX( GDCLDX, 0.D0 ), one ) ! GTCFRC( I,K ) = ( GDCLDX - GDCFRC( I,K ) )/DELTA ! ! GDMU2X = GDQ( I,K,IMU2 ) + GTQ( I,K,IMU2 )*DELTA ! GDMU2X = MIN( MAX( GDMU2X,VARMIN ),VARMAX ) ! GDMU3X = GDQ( I,K,IMU3 ) + GTQ( I,K,IMU3 )*DELTA ! GDMU3X = MIN( MAX( GDMU3X,-SKWMAX ),SKWMAX ) ! GTQ( I,K,IMU2 ) = ( GDMU2X - GDQ( I,K,IMU2 ))/DELTA ! GTQ( I,K,IMU3 ) = ( GDMU3X - GDQ( I,K,IMU3 ))/DELTA ! ! tem = DELP(I,K)*GRAVI ! HBGT(I) = HBGT(I) + (CP*GTT(I,K) + EL*GTQ(I,K,1) & ! - EMELT*GTQ(I,K,ITI)) * tem ! - EMELT*(GTQ(I,K,ITI)+GTIDET(I,K))) * tem ! WBGT(I) = WBGT(I) + (GTQ(I,K,1) + GTQ(I,K,ITL) + GTQ(I,K,ITI)) * tem ! + GTLDET(I,K) + GTIDET(I,K)) * tem ! ENDDO ! ENDDO ! ! DO I=ISTS,IENS ! HBGT(I) = HBGT(I) - EMELT*GSNWC(I) ! WBGT(I) = WBGT(I) + GPRCC(I,1) + GSNWC(I) ! CTOPP(I) = 1.D6 ! ENDDO ! ! The following commented out because they are unused ! DO CTP=1,NCTP ! DO I=ISTS, IENS ! kk = kt(i,ctp) ! IF (KK > KB(I) ) THEN ! CUMHGT(I,CTP) = GDZ(I,KK) ! CTOPP(I) = MIN(CTOPP(I), GDP(I,KK)) ! ELSE ! CUMHGT (I,CTP) = -999.D0 ! ENDIF ! ENDDO ! ENDDO ! DO I=ISTS,IENS ! IF(CTOPP(I) >= 1.D6) THEN ! CTOPP(I) = -999.D0 ! ENDIF ! ENDDO ! !> -# Ensures conservation of water. !In fact, no adjustment of the precip ! is occuring now which is a good sign! DD if(flx_form .and. adjustp) then DO I = ISTS, IENS if(gprcp(i,1)+gsnwp(i,1) > 1.e-12_kind_phys) then moistening_aw(i) = - moistening_aw(i) / (gprcp(i,1)+gsnwp(i,1)) else moistening_aw(i) = 1.0 endif ENDDO do k=1,kmax DO I = ISTS, IENS gprcp(i,k) = max(0.0, gprcp(i,k) * moistening_aw(i)) gsnwp(i,k) = max(0.0, gsnwp(i,k) * moistening_aw(i)) ENDDO enddo endif ! ! do i=ISTS,IENS ! GPRCC(I,1) = GPRCP(I,1) ! GSNWC(I ) = GSNWP(I,1) ! enddo do k=1,kmax do i=ISTS,IENS GTPRP(I,k) = GPRCP(I,k) + GSNWP(I,k) enddo enddo ! !DD provide GFS with a separate downdraft mass flux DO K=1,KMAX DO I=ISTS,IENS GMFX1(I,K) = GMFX0(I,K) - GMFLX(I,K) ENDDO ENDDO ! if (flx_form) then deallocate(sfluxterm, qvfluxterm, qlfluxterm, qifluxterm,& condtermt, condtermq, frzterm, prectermq, & prectermfrz, dtdwn, dqvdwn, dqldwn, & dqidwn, trfluxterm, dtrdwn) endif if (allocated(gprcc)) deallocate(gprcc) ! END SUBROUTINE CS_CUMLUS !> @} !*********************************************************************** !>\ingroup cs_scheme !! This subroutine calculates cloud base properties. SUBROUTINE CUMBAS & ! cloud base ( IJSDIM, KMAX , & !DD dimensions KB , GCYM , KBMX , & ! output ntr , ntrq , & GCHB , GCWB , GCUB , GCVB , & ! output GCIB , gctrb, & ! output GDH , GDW , GDHS , GDQS , & ! input GDQI , GDU , GDV , GDZM , & ! input GDPM , FDQS , GAM , & ! input lprnt, ipr, & ISTS , IENS , gctbl, gcqbl ,gdq, & gcwbl, gcqlbl, gcqibl, gctrbl ) ! input !DDsigmadiag add updraft profiles below cloud base ! ! IMPLICIT NONE ! integer, parameter :: crtrh=0.80 integer, parameter :: crtrh=0.70 INTEGER, INTENT(IN) :: IJSDIM, KMAX , ntr, ntrq ! DD, for GFS, pass in integer ipr logical lprnt ! ! [OUTPUT] INTEGER KB (IJSDIM) ! cloud base REAL(kind_phys) GCYM (IJSDIM, KMAX) ! norm. mass flux (half lev) INTEGER KBMX REAL(kind_phys) GCHB (IJSDIM) ! cloud base MSE REAL(kind_phys) GCWB (IJSDIM) ! cloud base total water REAL(kind_phys) GCUB (IJSDIM) ! cloud base U REAL(kind_phys) GCVB (IJSDIM) ! cloud base V REAL(kind_phys) GCIB (IJSDIM) ! cloud base ice REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr) ! cloud base tracer !DDsigma added to arglist for AW, subcloud updraft profiles: temperature, water vapor ! total water, cloud water, and cloud ice respectively REAL(kind_phys), dimension(ijsdim,kmax) :: gctbl, gcqbl, gcwbl, gcqlbl, gcqibl !>DDsigmadiag REAL(kind_phys), dimension(ijsdim,kmax,ntrq:ntr) :: gctrbl !DDsigmadiag ! ! [INPUT] REAL(kind_phys) GDH (IJSDIM, KMAX) ! moist static energy REAL(kind_phys) GDW (IJSDIM, KMAX) ! total water REAL(kind_phys) GDq (IJSDIM, KMAX, ntr) ! water vapor and tracer REAL(kind_phys) GDHS (IJSDIM, KMAX) ! saturate MSE REAL(kind_phys) GDQS (IJSDIM, KMAX) ! saturate humidity REAL(kind_phys) GDQI (IJSDIM, KMAX) ! cloud ice REAL(kind_phys) GDU (IJSDIM, KMAX) ! u-velocity REAL(kind_phys) GDV (IJSDIM, KMAX) ! v-velocity REAL(kind_phys) GDZM (IJSDIM, KMAX+1) ! Altitude (half lev) REAL(kind_phys) GDPM (IJSDIM, KMAX+1) ! pressure (half lev) REAL(kind_phys) FDQS (IJSDIM, KMAX) REAL(kind_phys) GAM (IJSDIM, KMAX) INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) CBASE (IJSDIM) ! one over cloud base height ! REAL(kind_phys) CBASEP(IJSDIM) ! cloud base pressure REAL(kind_phys) DELZ, GAMX, wrk ! REAL(kind_phys) DELZ, QSL, GAMX, wrk ! REAL(kind_phys), dimension(ijsdim,kmax) :: gchbl !DDsigmadiag real(kind_phys), dimension(ijsdim) :: gcqb, tx1, spbl, qsl INTEGER I, K, kp1, n ! ! [INTERNAL PARM] INTEGER :: KMAXM1 INTEGER :: KLCLB !! LCL base level INTEGER :: KCB !! fix cloud bottom INTEGER :: KBMAX !! cloud base max INTEGER :: KBOFS !! cloud base offset KMAXM1 = KMAX-1 KLCLB = 1 ! LCL base level KCB = 0 ! fix cloud bottom ! KBMAX = KMAXM1 ! cloud base max KBMAX = KMAX/2 ! cloud base max KBOFS = 0 ! cloud base offset ! do k=1,kmax do i=ists,iens GCYM(I,k) = zero enddo enddo ! IF (KCB > 0) THEN DO I=ISTS,IENS KB(I) = KCB ENDDO ELSE DO I=ISTS,IENS ! KB(I) = KBMAX KB(I) = -1 tx1(i) = one / gdpm(i,1) ENDDO DO K=KLCLB+1,KBMAX-1 DO I=ISTS,IENS GAMX = FDQS(I,K) / (one+GAM(I,K)) * oneocp QSL(i) = GDQS(I,K) + GAMX * (GDH(I,KLCLB)-GDHS(I,K)) spbl(i) = one - gdpm(i,k) * tx1(i) IF (GDW(I,KLCLB) >= QSL(i) .and. kb(i) < 0 & .and. spbl(i) >= spblmin) THEN ! .and. spbl(i) >= spblcrit .and. spbl(i) < spblcrit*10.0) THEN KB(I) = K + KBOFS ENDIF ENDDO ENDDO DO K=KLCLB+1,KBMAX-1 DO I=ISTS,IENS spbl(i) = one - gdpm(i,k) * tx1(i) IF (kb(i) > k .and. spbl(i) > spblmax) THEN KB(I) = K ENDIF ENDDO ENDDO ! DO K=KBMAX-1,KLCLB+1,-1 ! DO I=ISTS,IENS ! GAMX = FDQS(I,K) / (one+GAM(I,K)) * oneocp ! QSL(i) = GDQS(I,K) + GAMX * (GDH(I,KLCLB)-GDHS(I,K)) ! spbl(i) = one - gdpm(i,k) * tx1(i) ! IF (GDW(I,KLCLB) >= QSL(i) .and. spbl(i) >= spblcrit & ! .and. spbl(i) < spblcrit*6.0) THEN ! .and. spbl(i) < spblcrit*8.0) THEN ! KB(I) = K + KBOFS ! ENDIF ! ENDDO ! if(lprnt) write(0,*)' k=',k,' gdh1=',gdh(ipr,klclb),' gdhs=',gdhs(ipr,k),' kb=',kb(ipr) & ! ,' GDQS=',GDQS(ipr,k),' GDW=',GDW(ipr,KLCLB),' gdpm=',gdpm(ipr,k),' spbl=',spbl(ipr),' qsl=',qsl(ipr) ! ENDDO ENDIF ! do i=ists,iens tx1(i) = zero qsl(i) = zero enddo do k=1,kbmax do i=ists,iens if (k < kb(i)) then tx1(i) = tx1(i) + gdw(i,k) * (GDZM(i,k+1)-GDZM(i,k)) qsl(i) = qsl(i) + gdqs(i,k) * (GDZM(i,k+1)-GDZM(i,k)) endif enddo enddo do i=ists,iens if (qsl(i) > zero) tx1(i) = tx1(i) / qsl(i) if (tx1(i) < crtrh) kb(i) = -1 enddo ! KBMX = 1 DO I=ISTS,IENS KBMX = MAX(KBMX, KB(I)) if (kb(i) > 0) then CBASE (I) = one / (GDZM(I,KB(I)) - GDZM(I,1)) ! CBASEP(I) = GDPM(I,KB(I)) endif ENDDO ! DO K=2,KBMX DO I=ISTS,IENS IF (K <= KB(I)) THEN ! GCYM(I,K) = sqrt((GDZM(I,K)-GDZM(I,1))*CBASE(i)) GCYM(I,K) = (GDZM(I,K)-GDZM(I,1))*CBASE(i) ENDIF ENDDO ENDDO ! DO I=ISTS,IENS GCHB(I) = zero GCWB(I) = zero GCUB(I) = zero GCVB(I) = zero GCIB(I) = zero GCQB(I) = zero ENDDO do n = ntrq,ntr DO I=ISTS,IENS GCtrB(I,n) = zero enddo enddo do k=1,kmax do i=ists,iens ! GChbl(i,k) = zero gcqbl(i,k) = zero gcqlbl(i,k) = zero gcqibl(i,k) = zero gctbl(i,k) = zero gcwbl(i,k) = zero enddo enddo ! do n=ntrq,ntr ! do k=1,kmax ! do i=ists,iens ! gtrqbl(i,k,n) = zero ! enddo ! enddo ! enddo ! DO K=1,KBMX kp1 = min(k+1, kmax) DO I=ISTS,IENS IF (K < KB(I)) THEN DELZ = GCYM(I,Kp1) - GCYM(I,K) GCHB(I) = GCHB(I) + DELZ * GDH (I,K) GCWB(I) = GCWB(I) + DELZ * GDW (I,K) GCUB(I) = GCUB(I) + DELZ * GDU (I,K) GCVB(I) = GCVB(I) + DELZ * GDV (I,K) GCIB(I) = GCIB(I) + DELZ * GDQI(I,K) GCqB(I) = GCqB(I) + DELZ * GDQ (I,K,1) ! do n = ntrq,ntr ! GCtrB(I,n) = GCtrB(I,n) + DELZ * GDQ (I,K,n) ! enddo ! get subcloud profiles to pass out and do AW eddy flux tendencies ! removing the normalized mass flux weighting wrk = one / gcym(i,kp1) ! gchbl(i,kp1) = gchb(i) * wrk gcqbl(i,kp1) = gcqb(i) * wrk gcqibl(i,kp1) = gcib(i) * wrk gcwbl(i,kp1) = gcwb(i) * wrk gcqlbl(i,kp1) = gcwbl(i,kp1) - (gcqibl(i,kp1)+gcqbl(i,kp1)) ! gctbl(i,kp1) = (gchbl(i,kp1) - grav*gdzm(i,kp1) - el*gcqbl(i,kp1)) * oneocp gctbl(i,kp1) = (gchb(i)*wrk - grav*gdzm(i,kp1) - el*gcqbl(i,kp1)) * oneocp ! tracers do n=ntrq,ntr GCtrB(I,n) = GCtrB(I,n) + DELZ * GDQ (I,K,n) GCtrBl(I,kp1,n) = gctrb(i,n) * wrk enddo ENDIF ENDDO ENDDO ! END SUBROUTINE CUMBAS !*********************************************************************** !>\ingroup cs_scheme !! This subroutine calculates in-cloud properties. SUBROUTINE CUMUP & !! in-cloud properties ( IJSDIM, KMAX , NTR , ntrq , & !DD dimensions ACWF , & ! output GCLZ , GCIZ , GPRCIZ, GSNWIZ, & ! output GCYT , GCHT , GCQT , & ! output GCLT , GCIT , GTPRT , & ! output GCUT , GCVT , gctrt , & ! output KT , KTMX , & ! output GCYM , & ! modified wcv , & ! !DDsigma new output GCHB , GCWB , GCUB , GCVB , & ! input !DDsigmadiag GCIB , gctrb , & ! input GDU , GDV , GDH , GDW , & ! input GDHS , GDQS , GDT , GDTM , & ! input GDQ , GDQI , GDZ , GDZM , & ! input GDPM , FDQS , GAM , GDZTR , & ! input CPRES , WCB , & ! input ! CPRES , WCB , ERMR , & ! input KB , CTP , ISTS , IENS, & ! input gctm , gcqm , gcwm , gchm, gcwt,& gclm, gcim , gctrm , lprnt, ipr ) ! !DD AW the above line of arguments were previously local, and often scalars. ! Dimensions were added to them to save profiles for each grid point. ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, KMAX, NTR, ipr , ntrq ! DD, for GFS, pass in logical :: lprnt ! ! [OUTPUT] REAL(kind_phys) ACWF (IJSDIM) !< cloud work function REAL(kind_phys) GCLZ (IJSDIM, KMAX) !< cloud liquid water*eta REAL(kind_phys) GCIZ (IJSDIM, KMAX) !< cloud ice*eta REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) !< rain generation*eta REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) !< snow generation*eta REAL(kind_phys) GCYT (IJSDIM) !< norm. mass flux @top REAL(kind_phys) GCHT (IJSDIM) !< cloud top MSE*eta REAL(kind_phys) GCQT (IJSDIM) !< cloud top moisture*eta REAL(kind_phys) GCLT (IJSDIM) !< cloud top liquid water*eta REAL(kind_phys) GCIT (IJSDIM) !< cloud top ice*eta REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr) !< cloud top tracer*eta REAL(kind_phys) GTPRT (IJSDIM) !< cloud top (rain+snow)*eta REAL(kind_phys) GCUT (IJSDIM) !< cloud top u*eta REAL(kind_phys) GCVT (IJSDIM) !< cloud top v*eta REAL(kind_phys) GCwT (IJSDIM) !< cloud top v*eta INTEGER KT (IJSDIM) !< cloud top INTEGER KTMX !< max of cloud top REAL(kind_phys) WCV (IJSDIM, KMAX) !< updraft velocity (half lev) !DD sigma make output ! ! [MODIFIED] REAL(kind_phys) GCYM (IJSDIM, KMAX) !< norm. mass flux ! ! [INPUT] REAL(kind_phys) GCHB (IJSDIM) !< cloud base Moist Static Energy REAL(kind_phys) GCWB (IJSDIM) !< cloud base total water REAL(kind_phys) GCUB (IJSDIM) !< cloud base U REAL(kind_phys) GCVB (IJSDIM) !< cloud base V REAL(kind_phys) GCIB (IJSDIM) !< cloud base ice REAL(kind_phys) GCtrB (IJSDIM,ntrq:ntr) !< cloud base tracers REAL(kind_phys) GDU (IJSDIM, KMAX) !< U REAL(kind_phys) GDV (IJSDIM, KMAX) !< V REAL(kind_phys) GDH (IJSDIM, KMAX) !< moist static energy REAL(kind_phys) GDW (IJSDIM, KMAX) !< total water REAL(kind_phys) GDHS (IJSDIM, KMAX) !< saturation MSE REAL(kind_phys) GDQS (IJSDIM, KMAX) !< saturation q REAL(kind_phys) GDT (IJSDIM, KMAX) !< T REAL(kind_phys) GDTM (IJSDIM, KMAX+1) !< T (half lev) REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR) !< q !!DDsigmadiag REAL(kind_phys) GDQI (IJSDIM, KMAX) !< cloud ice REAL(kind_phys) GDZ (IJSDIM, KMAX) !< z REAL(kind_phys) GDZM (IJSDIM, KMAX+1) !< z (half lev) REAL(kind_phys) GDPM (IJSDIM, KMAX+1) !< p (half lev) REAL(kind_phys) FDQS (IJSDIM, KMAX) REAL(kind_phys) GAM (IJSDIM, KMAX) REAL(kind_phys) GDZTR (IJSDIM) !< tropopause height REAL(kind_phys) CPRES !< pres. fac. for cum. fric. REAL(kind_phys) WCB(ijsdim) !< cloud base updraft velocity**2 ! REAL(kind_phys) ERMR !< entrainment rate (ASMODE) INTEGER KB (IJSDIM) INTEGER CTP, ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) ACWFK (IJSDIM,KMAX) !< cloud work function REAL(kind_phys) ACWFN (IJSDIM,KMAX) !< negative part of cloud work function REAL(kind_phys) myGCHt !< cloud top h *eta (half lev) REAL(kind_phys) GCHMZ (IJSDIM, KMAX) !< cloud h *eta (half lev) REAL(kind_phys) GCWMZ (IJSDIM, KMAX) !< cloud Qt*eta (half lev) REAL(kind_phys) GCUMZ (IJSDIM, KMAX) !< cloud U *eta (half lev) REAL(kind_phys) GCVMZ (IJSDIM, KMAX) !< cloud V *eta (half lev) REAL(kind_phys) GCqMZ (IJSDIM ) !< cloud qv*eta (half lev) REAL(kind_phys) GCIMZ (IJSDIM, KMAX) !< cloud Qi*eta (half lev) REAL(kind_phys) GCtrMZ(IJSDIM, KMAX,ntrq:ntr)!< cloud tracer*eta (half lev) REAL(kind_phys) GTPRMZ(IJSDIM, KMAX) !< rain+snow *eta (half lev) ! REAL(kind_phys) BUOY (IJSDIM, KMAX) !< buoyancy REAL(kind_phys) BUOYM (IJSDIM, KMAX) !< buoyancy (half lev) REAL(kind_phys) WCM (IJSDIM ) !< updraft velocity**2 (half lev) ! REAL(kind_phys) WCM (IJSDIM, KMAX) !< updraft velocity**2 (half lev) !DD sigma make output REAL(kind_phys) WCV ( IJSDIM, KMAX+1 ) !! updraft velocity (half lev) REAL(kind_phys) GCY (IJSDIM, KMAX) !< norm. mass flux ! REAL(kind_phys) ELAR (IJSDIM, KMAX) !< entrainment rate REAL(kind_phys) ELAR !< entrainment rate at mid layer ! REAL(kind_phys) GCHM (IJSDIM, KMAX) !< cloud MSE (half lev) REAL(kind_phys) GCWM (IJSDIM, KMAX) !< cloud Qt (half lev) !DDsigmadiag REAL(kind_phys) GCTM (IJSDIM, KMAX) !< cloud T (half lev) !DDsigmadiag make output REAL(kind_phys) GCQM (IJSDIM, KMAX) !< cloud q (half lev) !DDsigmadiag make output REAL(kind_phys) GCLM (IJSDIM, KMAX) !< cloud liquid ( half lev) REAL(kind_phys) GCIM (IJSDIM, KMAX) !< cloud ice (half lev) REAL(kind_phys) GCUM (IJSDIM, KMAX) !< cloud U (half lev) REAL(kind_phys) GCVM (IJSDIM, KMAX) !< cloud V (half lev) REAL(kind_phys) GCtrM (IJSDIM, KMAX,ntrq:ntr) !< cloud tracer (half lev) ! REAL(kind_phys), dimension(IJSDIM) :: WCM_, ELARM1, GDZMKB REAL(kind_phys) GDQSM, GDHSM, GDQM, GDSM, GDCM, FDQSM, GCCM, & DELZ, ELADZ, DCTM , CPGMI, DELC, FICE, ELARM2,GCCMZ, & PRECR, GTPRIZ, DELZL, GCCT, DCT, WCVX, PRCZH, wrk INTEGER K, I, kk, km1, kp1, n ! CHARACTER CTNUM*2 ! !DD#ifdef OPT_CUMBGT !DD REAL(kind_phys) HBGT (IJSDIM) ! heat budget !DD REAL(kind_phys) WBGT (IJSDIM) ! water budget !DD REAL(kind_phys) PBGT (IJSDIM) ! precipitation budget !DD REAL(kind_phys) MBGT (IJSDIM) ! mass budget !DD REAL(kind_phys) GTPRX (IJSDIM) ! (rain+snow)*eta at top !DD REAL(kind_phys) GSNWT (IJSDIM) ! cloud top snow*eta !DD REAL(kind_phys) HBMX, WBMX, PBMX, MBMX !DD SAVE HBMX, WBMX, PBMX, MBMX !DD#endif ! ! [INTERNAL PARAM] REAL(kind_phys), parameter :: ZTREF = 1._kind_phys, ztrefi = one/ztref, & ELAMIN = zero, ELAMAX = 4.e-3 ! min and max entrainment rate REAL(kind_phys) :: PB = 1.0_kind_phys !m REAL(kind_phys) :: TAUZ = 5.0e3_kind_phys REAL(kind_phys) :: TAUZ = 1.0e4_kind_phys !m REAL(kind_phys) :: ELMD = 2.4e-3 ! for Neggers and Siebesma (2002) !m REAL(kind_phys) :: ELAMAX = 5.e-3 ! max. of entrainment rate ! REAL(kind_phys) :: WCCRT = zero !m REAL(kind_phys) :: WCCRT = 0.01 REAL(kind_phys) :: WCCRT = 1.0e-6_kind_phys, wvcrt=1.0e-3_kind_phys REAL(kind_phys) :: TSICE = 268.15_kind_phys ! compatible with macrop_driver REAL(kind_phys) :: TWICE = 238.15_kind_phys ! compatible with macrop_driver ! REAL(kind_phys) :: wfn_neg = 0.1 REAL(kind_phys) :: wfn_neg = 0.15 ! REAL(kind_phys) :: wfn_neg = 0.25 ! REAL(kind_phys) :: wfn_neg = 0.30 ! REAL(kind_phys) :: wfn_neg = 0.35 REAL(kind_phys) :: esat, tem ! REAL(kind_phys) :: esat, tem, rhs_h, rhs_q ! REAL(kind_phys) Z ! altitude REAL(kind_phys) ZH ! scale height REAL(kind_phys) T ! temperature ! ! ! Note: iteration is not made to diagnose cloud ice for simplicity ! do i=ists,iens ACWF (I) = zero GCYT (I) = zero GCHT (I) = zero GCQT (I) = zero GCLT (I) = zero GCIT (I) = zero GTPRT(I) = zero GCUT (I) = zero GCVT (I) = zero GCwT (I) = zero enddo do k=1,kmax do i=ists,iens ACWFK (I,k) = unset_kind_phys ACWFN (I,k) = unset_kind_phys GCLZ (I,k) = zero GCIZ (I,k) = zero GPRCIZ(I,k) = zero GSNWIZ(I,k) = zero ! GCHMZ (I,k) = zero GCWMZ (I,k) = zero GCIMZ (I,k) = zero GCUMZ (I,k) = zero GCVMZ (I,k) = zero GTPRMZ(I,k) = zero ! BUOY (I,k) = unset_kind_phys BUOYM (I,k) = unset_kind_phys WCV (I,k) = unset_kind_phys GCY (I,k) = unset_kind_phys ! GCHM (I,k) = unset_kind_phys GCWM (I,k) = unset_kind_phys GCTM (I,k) = unset_kind_phys GCQM (I,k) = unset_kind_phys GCLM (I,k) = unset_kind_phys GCIM (I,k) = unset_kind_phys GCUM (I,k) = unset_kind_phys GCVM (I,k) = unset_kind_phys enddo enddo do i=ists,iens GCqMZ(I) = zero WCM(I) = unset_kind_phys WCM_(I) = zero enddo ! tracers do n=ntrq,ntr do i=ists,iens GCtrT(I,n) = zero enddo do k=1,kmax do i=ists,iens GCTRM(I,k,n) = unset_kind_phys enddo enddo enddo ! DO I=ISTS,IENS ! if (kb(i) > 0) then ! GDZMKB(I) = GDZM(I,KB(I)) ! cloud base height ! endif ! ENDDO ! ! < cloud base properties > ! DO I=ISTS,IENS K = KB(I) if (k > 0) then GDZMKB(I) = GDZM(I,K) ! cloud base height GCHM(I,K) = GCHB(I) GCWM(I,K) = GCWB(I) WCM_(I) = WCB(i) GCUM(I,K) = GCUB(I) GCVM(I,K) = GCVB(I) do n = ntrq,ntr GCtrM(I,K,n) = GCtrB(I,n) enddo ! esat = min(gdpm(i,k), fpvs(gdtm(i,k))) GDQSM = min(EPSV*esat/max(gdpm(i,k)+epsvm1*esat, 1.0e-10), 0.1) gdsm = CP*GDTM(I,K) + GRAV*GDZMKB(I) ! dse GDHSM = gdsm + EL*GDQSM ! saturated mse ! FDQSM = FDQSAT(GDTM(I,K), GDQSM) tem = one / GDTM(I,K) FDQSM = GDQSM * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT ! tem = one / (CP+EL*FDQSM) DCTM = (GCHB(I) - GDHSM) * tem GCQM(I,K) = min(GDQSM + FDQSM*DCTM, GCWM(I,K)) GCCM = MAX(GCWM(I,K)-GCQM(I,K), zero) ! GCTM(I,K) = GDT(I,K) + DCTM ! old ! GCTM(I,K) = (GCHB(I) - gdsm - el*gcqm(i,k)) * oneocp + dctm ! new GCTM(I,K) = (GCHB(I) - grav*gdzm(i,k) - el*gcqm(i,k)) * oneocp + dctm ! new ! GCIM(I,K) = FRICE(GCTM(I,K)) * GCCM ! cloud base ice GCLM(I,K) = MAX(GCCM-GCIM(I,K), zero) ! cloud base liquid GCHM(I,K) = GCHM(I,K) + EMELT * (GCIM(I,K)-GCIB(I)) DCTM = (GCHM(I,K) - GDHSM) * tem ! GCTM(I,K) = dctm + GDT(I,K) + gocp*gdzm(i,k) !DD old AW convert to DSE GCTM(I,K) = dctm + (GCHB(I) - el*gcqm(i,k)) * oneocp ! new, make dse ! GDQM = half * (GDQ(I,K,1) + GDQ(I,K-1,1)) GDCM = half * (GDQ(I,K,ITL) + GDQI(I,K) & + GDQ(I,K-1,ITL) + GDQI(I,K-1)) ! BUOYM(I,K) = (DCTM*tem + EPSVT*(GCQM(I,K)-GDQM) - GCCM + GDCM )*GRAV ! ACWFK(I,K) = zero ACWFN(I,K) = zero ! !DD#ifdef OPT_ASMODE !DD ELARM1(I) = ERMR !DD#elif defined OPT_NS02 !DD ELARM1(I) = ELMD / SQRT(WCM(I,K)) !DD#else ! ELARM1(I) = CLMD*PA*BUOYM(I,K)/WCM(I,K) ! ELARM1(I) = min(max(CLMD*PA*BUOYM(I,K)/WCM_(I), ELAMIN), ELAMAX) !DD#endif ! ELARM1(I) = MIN(MAX(ELARM1(I), ELAMIN), ELAMAX) ! GCHMZ(I,K) = GCHM(I,K) GCWMZ(I,K) = GCWM(I,K) GCUMZ(I,K) = GCUM(I,K) GCVMZ(I,K) = GCVM(I,K) GCqMZ(I) = GCqM(I,K) GCIMZ(I,K) = GCIM(I,K) do n = ntrq,ntr GCtrMZ(I,K,n) = GCtrM(I,K,n) enddo endif ENDDO ! ! < in-cloud properties > ! DO K=3,KMAX km1 = k - 1 DO I=ISTS,IENS IF (kb(i) > 0 .and. K > KB(I) .AND. WCM_(I) > WCCRT) THEN WCV(I,KM1) = SQRT(MAX(WCM_(I), zero)) DELZ = GDZM(I,K) - GDZM(I,KM1) ELARM1(I) = min(max(CLMDPA*BUOYM(I,KM1)/WCM_(I), ELAMIN), ELAMAX) GCYM(I,K) = GCYM(I,KM1) * EXP(ELARM1(I)*DELZ) ELADZ = GCYM(I,K) - GCYM(I,KM1) ! GCHMZ(I,K) = GCHMZ(I,KM1) + GDH(I,KM1)*ELADZ GCWMZ(I,K) = GCWMZ(I,KM1) + GDW(I,KM1)*ELADZ ! esat = min(gdpm(i,k), fpvs(gdtm(i,k))) GDQSM = min(EPSV*esat/max(gdpm(i,k)+epsvm1*esat, 1.0e-10), 0.1) GDHSM = CP*GDTM(I,K ) + GRAV*GDZM(I,K) + EL*GDQSM ! FDQSM = FDQSAT(GDTM(I,K), GDQSM) tem = one / GDTM(I,K) FDQSM = GDQSM * tem * (fact1 + fact2*tem) ! calculate d(qs)/dT CPGMI = one / (CP + EL*FDQSM) PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) PRECR = FPREC(GDZM(I,K)-GDZMKB(I), PRCZH ) ! wrk = one / GCYM(I,K) DCTM = (GCHMZ(I,K)*wrk - GDHSM) * CPGMI GCQMZ(i) = min((GDQSM+FDQSM*DCTM)*GCYM(I,K), GCWMZ(I,K)) GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) GTPRMZ(I,K) = MAX(GTPRMZ(I,K), GTPRMZ(I,KM1)) GCCMZ = GCWMZ(I,K) - GCQMZ(i) - GTPRMZ(I,K ) DELC = MIN(GCCMZ, zero) GCCMZ = GCCMZ - DELC GCQMZ(i) = GCQMZ(i) + DELC ! FICE = FRICE(GDTM(I,K)+DCTM ) GCIMZ(I,K) = FICE * GCCMZ GSNWIZ(I,KM1) = FICE * (GTPRMZ(I,K)-GTPRMZ(I,KM1)) GCHMZ(I,K) = GCHMZ(I,K) + EMELT * (GCIMZ(I,K) + GSNWIZ(I,KM1) & - GCIMZ(I,KM1) - GDQI(I,KM1)*ELADZ) DCTM = (GCHMZ(I,K)*wrk - GDHSM) * CPGMI ! GDQM = half * (GDQ(I,K,1) + GDQ(I,KM1,1)) GDCM = half * (GDQ(I,K,ITL) + GDQI(I,K) & + GDQ(I,KM1,ITL) + GDQI(I,KM1)) GCQM(I,K) = wrk * GCQMZ(i) GCCM = wrk * GCCMZ ! BUOYM(I,K) = (DCTM*tem + EPSVT*(GCQM(I,K)-GDQM)-GCCM+GDCM) * GRAV BUOY(I,KM1) = half * (BUOYM(I,K)+BUOYM(I,KM1)) ! !DD#ifdef OPT_ASMODE !DD WCM(I,K ) & !DD = (WCM_(I) + 2.D0*PA*DELZ*BUOY(I,KM1) ) & !DD / (one + 2.D0*PB*DELZ*ERMR) !DD#elif OPT_NS02 !DD WCM(I,K ) = WCM_(I ) & !DD + 2.D0*DELZ*(PA*BUOYM(I,KM1)-ELMD*WCV(I,KM1)) !DD WCM(I,K ) = MAX(WCM(I,K ), zero ) !DD WCVX = SQRT(half*(WCM(I,K )+WCM_(I))) !DD WCM(I,K) = WCM_(I) + 2.D0*DELZ*(PA*BUOY(I,KM1)-ELMD*WCVX) !DD#else IF (BUOY(I,KM1) > zero) THEN WCM(I) = (WCM_(I) + CLMP*DELZ*BUOY(I,KM1)) / (one + DELZ/TAUZ) ELSE WCM(I) = (WCM_(I) + PA*(DELZ+DELZ)*BUOY(I,KM1) ) & / (one + DELZ/TAUZ + (DELZ+DELZ)*ELAMIN ) ENDIF !DD#endif ! !DD#ifdef OPT_ASMODE !DD ELARM2 = ERMR !DD#elif OPT_NS02 !DD ELARM2 = ELMD/SQRT(MAX(WCM(I,K), EPSln)) !DD#else ! ELARM2 = CLMD*PA*BUOYM(I,K) / MAX(WCM(I), EPSln) !DD#endif if (WCM(I) > zero) then ELARM2 = min(max(CLMDPA*BUOYM(I,K)/WCM(I),ELAMIN), ELAMAX) else ELARM2 = zero endif ELAR = half * (ELARM1(I) + ELARM2) GCYM(I,K) = GCYM(I,KM1) * EXP(ELAR*DELZ) ELADZ = GCYM(I,K) - GCYM(I,KM1) ! GCHMZ(I,K) = GCHMZ(I,KM1) + GDH(I,KM1)*ELADZ GCWMZ(I,K) = GCWMZ(I,KM1) + GDW(I,KM1)*ELADZ GCUMZ(I,K) = GCUMZ(I,KM1) + GDU(I,KM1)*ELADZ GCVMZ(I,K) = GCVMZ(I,KM1) + GDV(I,KM1)*ELADZ do n = ntrq,ntr GCtrMZ(I,K,n) = GCtrMZ(I,KM1,n) + GDq(I,KM1,n)*ELADZ enddo ! wrk = one / GCYM(I,K) DCTM = (GCHMZ(I,K)*wrk - GDHSM) * CPGMI GCQMZ(i) = min((GDQSM+FDQSM*DCTM)*GCYM(I,K), GCWMZ(I,K)) GTPRMZ(I,K) = PRECR * (GCWMZ(I,K)-GCQMZ(i)) GTPRMZ(I,K) = MAX(GTPRMZ(I,K), GTPRMZ(I,KM1)) GCCMZ = GCWMZ(I,K) - GCQMZ(i) - GTPRMZ(I,K) DELC = MIN(GCCMZ, zero) GCCMZ = GCCMZ - DELC GCQMZ(i) = GCQMZ(i) + DELC GCCM = wrk * GCCMZ GCQM(I,K) = wrk * GCQMZ(i) ! FICE = FRICE(GDTM(I,K)+DCTM ) GCIMZ(I,K) = FICE*GCCMZ GCIM(I,K) = GCIMZ(I,K)*wrk GCLM(I,K) = MAX(GCCM-GCIM(I,K), zero) GTPRIZ = GTPRMZ(I,K) - GTPRMZ(I,KM1) GSNWIZ(I,KM1) = FICE*GTPRIZ GPRCIZ(I,KM1) = (one-FICE )*GTPRIZ GCHMZ(I,K) = GCHMZ(I,K) + EMELT*(GCIMZ(I,K) + GSNWIZ(I,KM1) & - GCIMZ(I,KM1) - GDQI(I,KM1)*ELADZ ) GCHM(I,K) = GCHMZ(I,K)*wrk DCTM = (GCHM(I,K)-GDHSM) * CPGMI ! GCTM(I,K) = dctm + GDTM(I,K) + gocp*gdzm(i,k) ! old, make dse GCTM(I,K) = dctm + (GCHM(I,K) - el*gcqm(i,k)) * oneocp ! new, make dse ! GCWM(I,K) = GCWMZ(I,K) * wrk GCUM(I,K) = GCUMZ(I,K) * wrk GCVM(I,K) = GCVMZ(I,K) * wrk do n = ntrq,ntr GCtrM(I,K,n) = GCtrMZ(I,K,n) * wrk enddo DELZL = GDZ(I,KM1)-GDZM(I,KM1) GCY (I,KM1) = GCYM(I,KM1) * EXP(ELAR*DELZL) GCLZ(I,KM1) = half * (GCLM(I,K) + GCLM(I,KM1)) * GCY(I,KM1) GCIZ(I,KM1) = half * (GCIM(I,K) + GCIM(I,KM1)) * GCY(I,KM1) ! BUOYM(I,K) = (DCTM*tem + EPSVT*(GCQM(I,K)-GDQM)-GCCM+GDCM) * GRAV BUOY(I,KM1) = half * (BUOYM(I,K)+BUOYM(I,KM1)) ! IF (BUOY(I,KM1) > zero) THEN WCM(I) = (WCM_(I) + CLMP*DELZ*BUOY(I,KM1)) / (one + DELZ/TAUZ) ELSE WCM(I) = (WCM_(I) + PA*(DELZ+DELZ)*BUOY(I,KM1) ) & / (one + DELZ/TAUZ + (DELZ+DELZ)*ELAMIN ) ENDIF ! ! IF (BUOY(I,KM1) > zero) THEN ! ACWF(I) = ACWF(I) + BUOY(I,KM1)*GCY(I,KM1)*DELZ ! ENDIF ! ACWF(I) = ACWF(I) + BUOY(I,KM1)*GCY(I,KM1)*DELZ !!! wrk = BUOY(I,KM1)*GCY(I,KM1)*DELZ !!! ACWFK(I,K) = ACWFK(I,KM1) + wrk !!! ACWFN(I,K) = ACWFN(I,KM1) - min(wrk,0.0) ! ACWFN(I,K) = ACWFN(I,KM1) + abs(min(wrk,0.0)) ! wrk = BUOY(I,KM1)*GCY(I,KM1)*DELZ ACWFK(I,K) = ACWFK(I,KM1) + wrk ACWFN(I,K) = ACWFN(I,KM1) - min(wrk,0.0) WCM_(I) = WCM(I) ! if (lprnt .and. i == ipr) write(0,*) ' in cumup k=',k,' km1=',km1,' WCM_=',WCM_(I),' gcy=',gcy(i,km1),' buoym=',buoym(i,km1) ENDIF ! IF (K > KB(I) .AND. WCM_(I) > WCCRT) THEN ENDDO ENDDO ! ! < find cloud top > ! DO I=ISTS,IENS KT(I) = -1 enddo DO K=KMAX,2,-1 DO I=ISTS,IENS if (kb(i) > 0 .and. k > kb(i) .and. ACWFK(I,K) > 1.0e-10) then wrk = ACWFN(I,K) / ACWFK(I,K) IF (KT(I) == -1 .and. wrk < wfn_neg .AND. WCV(I,K) > WVCRT) THEN KT(I) = K ACWF(I) = ACWFK(I,K) ENDIF endif ENDDO ENDDO ! if (lprnt .and. kt(ipr) > 0) write(0,*) ' in cumup kt=',kt(ipr),' gcy=',gcy(ipr,kt(ipr)) ! KTMX = 2 DO I=ISTS,IENS kt(i) = min(kt(i), kmax-1) KTMX = max(ktmx, KT(I)) ENDDO ! DO I=ISTS,IENS kk = max(1, kt(i)+1) do k=kk,kmax GCYM (I,K) = zero GCLZ (I,K) = zero GCIZ (I,K) = zero GPRCIZ(I,K) = zero GSNWIZ(I,K) = zero enddo ENDDO ! if (lprnt .and. kt(ipr) > 0) write(0,*) ' in cumup2 kt=',kt(ipr),' gcy=',gcy(ipr,kt(ipr)) ! ! < cloud top properties > ! DO I=ISTS,IENS IF (kb(i) > 0 .and. KT(I) > kb(i)) THEN K = KT(I) kp1 = k + 1 GCYT(I) = GCY(I,K) ELADZ = GCYT(I) - GCYM(I,K) ! GCHT(I) = GCHMZ(I,K) + GDH(I,K)*ELADZ GCWT(i) = GCWMZ(I,K) + GDW(I,K)*ELADZ GCUT(I) = GCUMZ(I,K) + GDU(I,K)*ELADZ GCVT(I) = GCVMZ(I,K) + GDV(I,K)*ELADZ do n = ntrq,NTR GCtrT(I,n) = GCtrMZ(I,K,n) + GDq(I,K,n)*ELADZ enddo ! wrk = one / gcyt(i) DCT = (GCHT(I)*wrk - GDHS(I,K)) / (CP*(one + GAM(I,K))) GCQT(I) = min((GDQS(I,K) + FDQS(I,K)*DCT) * GCYT(I), GCWT(i)) PRCZH = PRECZH * MIN(GDZTR(I)*ZTREFI, one) GTPRT(I) = FPREC(GDZ(I,K)-GDZMKB(I), PRCZH) * (GCWT(i)-GCQT(I)) GTPRT(I) = MAX(GTPRT(I), GTPRMZ(I,K)) GCCT = GCWT(i) - GCQT(I) - GTPRT(I) DELC = MIN(GCCT, zero) GCCT = GCCT - DELC GCQT(I) = GCQT(I) + DELC ! FICE = FRICE(GDT(I,K)+DCT) GCIT(I) = FICE*GCCT GCLT(I) = (one-FICE) * GCCT GTPRIZ = GTPRT(I) - GTPRMZ(I,K) GPRCIZ(I,K) = (one-FICE) * GTPRIZ GSNWIZ(I,K) = FICE * GTPRIZ GCHT(I) = GCHT(I) & + EMELT * (GCIT(I) + GSNWIZ(I,K) - GCIMZ(I,K) - GDQI(I,K)*ELADZ) ! GCUT(I) = GCUT(I)*(one-CPRES) + GCY(I,K)*GDU(I,K)*CPRES GCVT(I) = GCVT(I)*(one-CPRES) + GCY(I,K)*GDV(I,K)*CPRES do n = ntrq,NTR ! GCtrT(I,n) = GCtrT(I,n)*(one-CPRES) + GCY(I,K)*GDq(I,K,n)*CPRES GCtrT(I,n) = GCtrT(I,n) + GCY(I,K)*GDq(I,K,n) enddo GCLZ(I,K) = GCLT(I) GCIZ(I,K) = GCIT(I) !DD AW get the cloud top values denormalized and put into profile mygcht = gcht(I) - el*(gcwt(i) - gcqt(i)) gctm(i,kp1) = wrk * (mygcht - el*gcqt(i)) * oneocp !Moorthi gcqm(i,kp1) = gcqt(i) gcqm(i,kp1) = gcqt(i)*wrk ! check this - oct17 2016 gcim(i,kp1) = gcit(i)*wrk gclm(i,kp1) = gclt(i)*wrk do n = ntrq,NTR gctrm(i,kp1,n) = gctrt(i,n)*wrk enddo ! ENDIF ENDDO ! !DD#ifdef OPT_CUMBGT /* budget check */ !DD HBGT ( ISTS:IENS ) = 0.D0 !DD WBGT ( ISTS:IENS ) = 0.D0 !DD PBGT ( ISTS:IENS ) = 0.D0 !DD MBGT ( ISTS:IENS ) = 0.D0 !DD GTPRX( ISTS:IENS ) = 0.D0 !DD GSNWT( ISTS:IENS ) = 0.D0 !DD! !DD IF ( CTP .EQ. 1 ) THEN !DD HBMX = 0.D0 !DD WBMX = 0.D0 !DD PBMX = 0.D0 !DD MBMX = 0.D0 !DD END IF !DD! !DD DO K = 2, KMAX !DD DO I = ISTS, IENS !DD IF ( K .GE. KB( I ) .AND. K .LT. KT( I ) ) THEN !DD ELADZ = GCYM( I,K+1 ) - GCYM( I,K ) !DD DELZ = GDZM( I,K+1 ) - GDZM( I,K ) !DD HBGT( I ) = HBGT( I ) + ( GDH( I,K )-EMELT*GDQI( I,K ) )*ELADZ !DD WBGT( I ) = WBGT( I ) + GDW( I,K )*ELADZ !DD MBGT( I ) = MBGT( I ) + ELAM( I,K )*DELZ !DD GTPRX( I ) = GTPRX( I ) + GPRCIZ( I,K ) + GSNWIZ( I,K ) !DD GSNWT( I ) = GSNWT( I ) + GSNWIZ( I,K ) !DD END IF !DD END DO !DD END DO !DD! !DD DO I = ISTS, IENS !DD IF ( KT( I ) .GT. KB( I ) ) THEN !DD ELADZ = GCYT( I ) - GCYM( I,KT(I) ) !DD DELZ = GDZ( I,KT(I) )-GDZM( I,KT(I) ) !DD GTPRX( I ) = GTPRX( I ) + GPRCIZ( I,KT(I) ) + GSNWIZ( I,KT(I) ) !DD GSNWT( I ) = GSNWT( I ) + GSNWIZ( I,KT(I) ) !DD HBGT( I ) = HBGT( I ) + GCHB( I ) - EMELT*GCIB( I ) & !DD + ( GDH( I,KT(I) )-EMELT*GDQI( I,KT(I) ) ) *ELADZ & !DD - ( GCHT(I)-EMELT*( GCIT(I)+GSNWT(I) ) ) !DD WBGT( I ) = WBGT( I ) & !DD + GCWB( I ) + GDW( I,KT(I) )*ELADZ & !DD - GCQT( I ) - GCLT( I ) - GCIT( I ) & !DD - GTPRT( I ) !DD MBGT( I ) = MBGT( I ) + one + ELAM( I,KT(I) )*DELZ & !DD - GCYT( I ) !DD PBGT( I ) = GTPRT( I ) - GTPRX( I ) !DD! !DD IF ( ABS( HBGT(I) ) .GT. ABS( HBMX ) ) HBMX = HBGT(I) !DD IF ( ABS( WBGT(I) ) .GT. ABS( WBMX ) ) WBMX = WBGT(I) !DD IF ( ABS( PBGT(I) ) .GT. ABS( PBMX ) ) PBMX = PBGT(I) !DD IF ( ABS( MBGT(I) ) .GT. ABS( MBMX ) ) MBMX = MBGT(I) !DD END IF !DD END DO !DD! !DD IF ( CTP .EQ. NCTP ) THEN !DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): energy imbalance =', HBMX !DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): water imbalance =', WBMX !DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): precipitation imbalance =', PBMX !DD WRITE( iulog,* ) '### CUMUP(rank=',irank,'): mass imbalance =', MBMX !DD END IF !DD#endif ! ! WRITE( CTNUM, '(I2.2)' ) CTP ! contains pure function FPREC(Z,ZH) implicit none real(kind_phys), intent(in) :: Z real(kind_phys), intent(in) :: ZH real(kind_phys) :: FPREC FPREC = MIN(MAX(one-EXP(-(Z-PRECZ0)/ZH), zero), one) end function FPREC pure function FRICE(T) implicit none real(kind_phys), intent(in) :: T real(kind_phys) :: FRICE FRICE = MIN(MAX((TSICE-T)/(TSICE-TWICE), zero), one) end function FRICE END SUBROUTINE CUMUP !*********************************************************************** !>\ingroup cs_scheme !! This subroutine computes cloud base mass flux. SUBROUTINE CUMBMX & !! cloud base mass flux ( IJSDIM, KMAX , & !DD dimensions CBMFX , & ! modified ACWF , GCYT , GDZM , & ! input GDW , GDQS , DELP , & ! input KT , KTMX , KB , & ! input DELT , ISTS , IENS ) ! input ! ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, KMAX ! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) CBMFX (IJSDIM) !< cloud base mass flux ! ! [INPUT] REAL(kind_phys) ACWF (IJSDIM) !< cloud work function REAL(kind_phys) GCYT (IJSDIM) !< norm mass flux @top REAL(kind_phys) GDZM (IJSDIM, KMAX+1) !< height REAL(kind_phys) GDW (IJSDIM, KMAX) !< total water REAL(kind_phys) GDQS (IJSDIM, KMAX) !< saturate humidity REAL(kind_phys) DELP (IJSDIM, KMAX) !< delt pressure INTEGER KT (IJSDIM) !< cloud top INTEGER KTMX !< max. of cloud top INTEGER KB (IJSDIM) !< cloud base REAL(kind_phys) DELT !< time step INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys), dimension(ijsdim) :: QX, QSX, RHM INTEGER I, K REAL(kind_phys) ALP, FMAX1, wrk ! ! [INTERNAL PARAM] REAL(kind_phys) :: FMAX = 1.5e-2_kind_phys ! maximum flux ! REAL(kind_phys) :: RHMCRT = zero ! critical val. of cloud mean RH ! REAL(kind_phys) :: RHMCRT = 0.25_kind_phys ! critical val. of cloud mean RH REAL(kind_phys) :: RHMCRT = 0.50_kind_phys ! critical val. of cloud mean RH REAL(kind_phys) :: ALP1 = zero REAL(kind_phys) :: TAUD = 1.e3_kind_phys ! REAL(kind_phys) :: TAUD = 6.e2_kind_phys REAL(kind_phys) :: ZFMAX = 3.5e3_kind_phys REAL(kind_phys) :: ZDFMAX = 5.e2_kind_phys ! REAL(kind_phys) :: FMAXP = 2._kind_phys ! do i=ists,iens qx(i) = zero qsx(i) = zero enddo ! DO K=1,KTMX DO I=ISTS,IENS IF (kb(i) > 0 .and. K >= KB(I) .AND. K <= KT(I)) THEN QX (I) = QX (I) + GDW (I,K) * DELP(I,K) QSX(I) = QSX(I) + GDQS(I,K) * DELP(I,K) ENDIF ENDDO ENDDO DO I=ISTS,IENS RHM(I) = min(one, max(zero, QX(I)/MAX(QSX(I),EPSln))) ENDDO ! wrk = one + delt/(taud+taud) DO I=ISTS,IENS k = kb(i) IF (k > 0 .and. KT(I) > K .AND. RHM(I) >= RHMCRT) THEN cbmfx(i) = max(cbmfx(i), zero) ALP = ALP0 + ALP1 * (GDZM(I,KT(I))-GDZM(I,K)) FMAX1 = (one - TANH((GDZM(I,1)-ZFMAX)/ZDFMAX)) * half ! FMAX1 = FMAX * FMAX1**FMAXP FMAX1 = FMAX * FMAX1 * FMAX1 ! CBMFX(I) = CBMFX(I) + MAX(ACWF(I), zero)/(ALP+ALP)*DELT ! CBMFX(I) = CBMFX(I) / (one + DELT/(TAUD+TAUD)) CBMFX(I) = (CBMFX(I) + ACWF(I)*(delt/(ALP+ALP))) * wrk CBMFX(I) = MIN(max(CBMFX(I), zero), FMAX1/GCYT(I)) ELSE CBMFX(I) = zero ENDIF ENDDO ! END SUBROUTINE CUMBMX !*********************************************************************** !>\ingroup cs_scheme !! This subroutine computes cloud mass flux & precip. SUBROUTINE CUMFLX & !! cloud mass flux ( IM , IJSDIM, KMAX , & !DD dimensions GMFLX , GPRCI , GSNWI , CMDET, & ! output QLIQ , QICE , GTPRC0, & ! output CBMFX , GCYM , GPRCIZ, GSNWIZ, & ! input GTPRT , GCLZ , GCIZ , GCYT, & ! input KB , KT , KTMX , & ! input ISTS , IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, KMAX, IM !! DD, for GFS, pass in ! ! [OUTPUT] REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux REAL(kind_phys) CMDET (IJSDIM, KMAX) !< detrainment mass flux REAL(kind_phys) GPRCI (IJSDIM, KMAX) !< rainfall generation REAL(kind_phys) GSNWI (IJSDIM, KMAX) !< snowfall generation REAL(kind_phys) QLIQ (IJSDIM, KMAX) !< cloud liquid REAL(kind_phys) QICE (IJSDIM, KMAX) !< cloud ice REAL(kind_phys) GTPRC0(IJSDIM) !< precip. before evap. ! ! [INPUT] REAL(kind_phys) CBMFX (IJSDIM) !< cloud base mass flux REAL(kind_phys) GCYM (IJSDIM, KMAX) !< normalized mass flux REAL(kind_phys) GCYT (IJSDIM) !< detraining mass flux REAL(kind_phys) GPRCIZ(IJSDIM, KMAX) !< precipitation/M REAL(kind_phys) GSNWIZ(IJSDIM, KMAX) !< snowfall/M REAL(kind_phys) GTPRT (IJSDIM) !< rain+snow @top REAL(kind_phys) GCLZ (IJSDIM, KMAX) !< cloud liquid/M REAL(kind_phys) GCIZ (IJSDIM, KMAX) !< cloud ice/M INTEGER KB (IJSDIM) !< cloud base INTEGER KT (IJSDIM) !< cloud top INTEGER KTMX !< max of cloud top INTEGER ISTS, IENS, I, K ! DO K=1,KTMX DO I=ISTS,IENS if (kb(i) > 0) then GMFLX(I,K) = GMFLX(I,K) + CBMFX(I) * GCYM(I,K) GPRCI(I,K) = GPRCI(I,K) + CBMFX(I) * GPRCIZ(I,K) GSNWI(I,K) = GSNWI(I,K) + CBMFX(I) * GSNWIZ(I,K) QLIQ(I,K) = QLIQ (I,K) + CBMFX(I) * GCLZ(I,K) QICE(I,K) = QICE (I,K) + CBMFX(I) * GCIZ(I,K) endif ENDDO ENDDO ! DO I= ISTS,IENS k = kt(i) if (kb(i) > 0 .and. k > kb(i)) then GTPRC0(I) = GTPRC0(I) + CBMFX(I) * GTPRT(I) CMDET(I,K) = CMDET(I,K) + CBMFX(I) * GCYT(I) endif ENDDO ! END SUBROUTINE CUMFLX !*********************************************************************** !>\ingroup cs_scheme !! This subroutine calculates cloud detrainment heating. SUBROUTINE CUMDET & !! detrainment ( im , IJSDIM, KMAX , NTR , ntrq , & !DD dimensions GTT , GTQ , GTU , GTV , & ! modified GDH , GDQ , GDU , GDV , & ! input ! GTT , GTQ , GTCFRC, GTU , GTV , & ! modified ! GDH , GDQ , GDCFRC, GDU , GDV , & ! input CBMFX , GCYT , DELPI , GCHT , GCQT , & ! input GCLT , GCIT , GCUT , GCVT , GDQI , & ! input gctrt, & KT , ISTS , IENS , nctp ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, nctp, ntrq !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTT (IJSDIM, KMAX) !< temperature tendency REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) !< moisture tendency ! REAL(kind_phys) GTCFRC(IJSDIM, KMAX) !< cloud fraction tendency REAL(kind_phys) GTU (IJSDIM, KMAX) !< u tendency REAL(kind_phys) GTV (IJSDIM, KMAX) !< v tendency ! ! [INPUT] REAL(kind_phys) GDH (IJSDIM, KMAX) !< moist static energy REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR) !< humidity qv ! REAL(kind_phys) GDCFRC(IJSDIM, KMAX) !< cloud fraction REAL(kind_phys) GDU (IJSDIM, KMAX) REAL(kind_phys) GDV (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) REAL(kind_phys) CBMFX (IM, NCTP) !< cloud base mass flux REAL(kind_phys) GCYT (IJSDIM, NCTP) !< detraining mass flux REAL(kind_phys) GCHT (IJSDIM, NCTP) !< detraining MSE REAL(kind_phys) GCQT (IJSDIM, NCTP) !< detraining qv REAL(kind_phys) GCLT (IJSDIM, NCTP) !< detraining ql REAL(kind_phys) GCIT (IJSDIM, NCTP) !< detraining qi REAL(kind_phys) GCtrT (IJSDIM, ntrq:ntr, NCTP)!< detraining tracer REAL(kind_phys) GCUT (IJSDIM, NCTP) !< detraining u REAL(kind_phys) GCVT (IJSDIM, NCTP) !< detraining v REAL(kind_phys) GDQI (IJSDIM, KMAX) !< cloud ice INTEGER KT (IJSDIM, NCTP) !< cloud top INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) GTHCI, GTQVCI, GTXCI integer I, K, CTP, kk,n ! DO CTP=1,NCTP DO I=ISTS,IENS K = KT(I,CTP) IF (K > 0) THEN GTXCI = DELPI(I,K)*CBMFX(I,CTP) GTHCI = GTXCI * (GCHT(I,CTP) - GCYT(I,CTP)*GDH(I,K)) GTQVCI = GTXCI * (GCQT(I,CTP) - GCYT(I,CTP)*GDQ(I,K,1)) ! GTQ(I,K,1) = GTQ(I,K,1) + GTQVCI GTT(I,K) = GTT(I,K) + (GTHCI - EL*GTQVCI) * oneocp ! ql tendency by detrainment is treated by stratiform scheme GTQ(I,K,ITL) = GTQ(I,K,ITL) + GTXCI * (GCLT(I,CTP) - GCYT(I,CTP)*GDQ(I,K,ITL)) ! qi tendency by detrainment is treated by stratiform scheme GTQ(I,K,ITI) = GTQ(I,K,ITI) + GTXCI * (GCIT(I,CTP) - GCYT(I,CTP)*GDQI(I,K)) do n = ntrq,NTR GTQ(I,K,n) = GTQ(I,K,n) + GTXCI * (GCtrT(I,n,CTP) - GCYT(I,CTP)*GDQ(I,K,n)) enddo ! GTCFRC(I,K) = GTCFRC(I,K) + GTXCI * (GCYT(I,CTP) - GCYT(I,CTP)*GDCFRC(I,K)) GTU(I,K) = GTU(I,K) + GTXCI * (GCUT(I,CTP) - GCYT(I,CTP)*GDU(I,K)) GTV(I,K) = GTV(I,K) + GTXCI * (GCVT(I,CTP) - GCYT(I,CTP)*GDV(I,K)) ENDIF ENDDO ENDDO ! END SUBROUTINE CUMDET !*********************************************************************** !>\ingroup cs_scheme SUBROUTINE CUMSBH & !! adiabat. descent ( IM , IJSDIM, KMAX , NTR, ntrq, & !DD dimensions GTT , GTQ , & ! modified GTU , GTV , & ! modified GDH , GDQ , GDQI , & ! input GDU , GDV , & ! input DELPI , GMFLX , GMFX0 , & ! input KTMX , CPRES , KB, ISTS , IENS ) ! input ! ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, IM, KMAX, NTR, ntrq !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTT (IJSDIM, KMAX) !< Temperature tendency REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) !< Moisture etc tendency REAL(kind_phys) GTU (IJSDIM, KMAX) !< u tendency REAL(kind_phys) GTV (IJSDIM, KMAX) !< v tendency ! ! [INPUT] REAL(kind_phys) GDH (IJSDIM, KMAX) REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR) !< humidity etc REAL(kind_phys) GDQI (IJSDIM, KMAX) REAL(kind_phys) GDU (IJSDIM, KMAX) REAL(kind_phys) GDV (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux (updraft+downdraft) REAL(kind_phys) GMFX0 (IJSDIM, KMAX) !< mass flux (updraft only) INTEGER KB(IJSDIM) !< cloud base index - negative means no convection INTEGER KTMX REAL(kind_phys) CPRES !< pressure factor for cumulus friction INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) SBH0, SBQ0, SBL0, SBI0, SBC0, SBS0, & SBH1, SBQ1, SBL1, SBI1, SBC1, SBS1, FX1, & SBU0, SBV0, SBU1, SBV1, GTHCI, GTQVCI, & GTQLCI, GTQICI, GTM2CI, GTM3CI, wrk, wrk1 REAL(kind_phys) FX(ISTS:IENS) REAL(kind_phys), dimension(IJSDIM, KMAX) :: GTLSBH, GTISBH integer :: I, K, KM, KP, n ! ! FX = zero do k=1,kmax do i=ists,iens GTLSBH(i,k) = zero GTISBH(i,k) = zero enddo enddo ! DO K=KTMX,1,-1 KM = MAX(K-1, 1) KP = MIN(K+1, KMAX) DO I=ISTS,IENS if (kb(i) > 0) then SBH0 = GMFLX(I,KP) * (GDH(I,KP)-GDH(I,K)) SBQ0 = GMFLX(I,KP) * (GDQ(I,KP,1)-GDQ(I,K,1)) SBL0 = GMFLX(I,KP) * (GDQ(I,KP,ITL )-GDQ(I,K,ITL)) SBI0 = GMFLX(I,KP) * (GDQI(I,KP)-GDQI(I,K)) SBU0 = GMFLX(I,KP) * (GDU(I,KP)-GDU(I,K)) & - GMFX0(I,KP) * (GDU(I,KP)-GDU(I,K))*CPRES SBV0 = GMFLX(I,KP) * (GDV(I,KP)-GDV(I,K)) & - GMFX0(I,KP) * (GDV(I,KP)-GDV(I,K))*CPRES ! SBH1 = GMFLX(I,K) * (GDH(I,K)-GDH(I,KM)) SBQ1 = GMFLX(I,K) * (GDQ(I,K,1)-GDQ(I,KM,1)) SBL1 = GMFLX(I,K) * (GDQ(I,K,ITL)-GDQ(I,KM,ITL)) SBI1 = GMFLX(I,K) * (GDQI(I,K)-GDQI(I,KM)) SBU1 = GMFLX(I,K) * (GDU(I,K)-GDU(I,KM)) & - GMFX0(I,K) * (GDU(I,K)-GDU(I,KM))*CPRES SBV1 = GMFLX(I,K) * (GDV(I,K)-GDV(I,KM)) & - GMFX0(I,K) * (GDV(I,K)-GDV(I,KM))*CPRES ! IF (GMFLX(I,K) > GMFLX(I,KP)) THEN FX1 = half ELSE FX1 = zero ENDIF ! wrk = DELPI(I,K) wrk1 = one - FX(I) GTHCI = wrk * (wrk1*SBH0 + FX1 *SBH1) GTQVCI = wrk * (wrk1*SBQ0 + FX1 *SBQ1) GTQLCI = wrk * (wrk1*SBL0 + FX1 *SBL1) GTQICI = wrk * (wrk1*SBI0 + FX1 *SBI1) ! GTT (I,K ) = GTT(I,K) + (GTHCI-EL*GTQVCI)*oneocp GTQ (I,K,1 ) = GTQ(I,K,1) + GTQVCI GTQ (I,K,ITL) = GTQ(I,K,ITL) + GTQLCI GTQ (I,K,ITI) = GTQ(I,K,ITI) + GTQICI GTU (I,K) = GTU(I,K) + wrk * (wrk1*SBU0 + FX1*SBU1) GTV (I,K) = GTV(I,K) + wrk * (wrk1*SBV0 + FX1*SBV1) DO n = ntrq, ntr GTQ (I,K,n) = GTQ(I,K,n) + wrk & * ( wrk1 * (GMFLX(I,KP) * (GDQ(I,KP,n)-GDQ(I,K ,n))) & + FX1 * (GMFLX(I,K ) * (GDQ(I,K ,n)-GDQ(I,KM,n))) ) ENDDO GTLSBH(I,K) = GTQLCI GTISBH(I,K) = GTQICI ! ! SBC0 = GMFLX(I,K+1) * (GDQ(I,KP,IMU2)-GDQ(I,K,IMU2)) ! SBS0 = GMFLX(I,K+1) * (GDQ(I,KP,IMU3)-GDQ(I,K,IMU3)) ! SBC1 = GMFLX(I,K) * (GDQ(I,K,IMU2)-GDQ(I,KM,IMU2)) ! SBS1 = GMFLX(I,K) * (GDQ(I,K,IMU3)-GDQ(I,KM,IMU3)) ! GTM2CI = DELPI(I,K) * (( one-FX(I))*SBC0 + FX1 *SBC1) ! GTM3CI = DELPI(I,K) * ((one-FX(I))*SBS0 + FX1 *SBS1) ! GTQ(I,K,IMU2) = GTQ(I,K,IMU2) + GTM2CI ! GTQ(I,K,IMU3) = GTQ(I,K,IMU3) + GTM3CI ! FX(I) = FX1 endif enddo enddo ! END SUBROUTINE CUMSBH !*********************************************************************** ! !*********************************************************************** !>\ingroup cs_scheme !! This subroutine calculate cloud subsidence heating. SUBROUTINE CUMSBW & !! adiabat. descent ( IM , IJSDIM, KMAX , & !DD dimensions GTU , GTV , & ! modified GDU , GDV , & ! input DELPI , GMFLX , GMFX0 , & ! input KTMX , CPRES , KB, ISTS , IENS ) ! input ! ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, IM, KMAX!! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTU (IJSDIM, KMAX) !< u tendency REAL(kind_phys) GTV (IJSDIM, KMAX) !< v tendency ! ! [INPUT] REAL(kind_phys) GDU (IJSDIM, KMAX) REAL(kind_phys) GDV (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux (updraft+downdraft) REAL(kind_phys) GMFX0 (IJSDIM, KMAX) !< mass flux (updraft only) INTEGER KB(IJSDIM) !< cloud base index - negative means no convection INTEGER KTMX, ISTS, IENS REAL(kind_phys) CPRES !< pressure factor for cumulus friction ! ! [INTERNAL WORK] REAL(kind_phys) FX1, SBU0, SBV0, SBU1, SBV1, wrk, wrk1 REAL(kind_phys) FX(ISTS:IENS) integer :: I, K, KM, KP ! ! FX = zero ! DO K=KTMX,1,-1 KM = MAX(K-1, 1) KP = MIN(K+1, KMAX) DO I=ISTS,IENS if (kb(i) > 0) then SBU0 = GMFLX(I,KP) * (GDU(I,KP)-GDU(I,K)) & - GMFX0(I,KP) * (GDU(I,KP)-GDU(I,K))*CPRES SBV0 = GMFLX(I,KP) * (GDV(I,KP)-GDV(I,K)) & - GMFX0(I,KP) * (GDV(I,KP)-GDV(I,K))*CPRES ! SBU1 = GMFLX(I,K) * (GDU(I,K)-GDU(I,KM)) & - GMFX0(I,K) * (GDU(I,K)-GDU(I,KM))*CPRES SBV1 = GMFLX(I,K) * (GDV(I,K)-GDV(I,KM)) & - GMFX0(I,K) * (GDV(I,K)-GDV(I,KM))*CPRES ! IF (GMFLX(I,K) > GMFLX(I,KP)) THEN FX1 = half ELSE FX1 = zero ENDIF ! wrk = DELPI(I,K) wrk1 = one - FX(I) ! GTU(I,K) = GTU(I,K) + wrk * (wrk1*SBU0 + FX1*SBU1) GTV(I,K) = GTV(I,K) + wrk * (wrk1*SBV0 + FX1*SBV1) ! FX(I) = FX1 endif enddo enddo ! END SUBROUTINE CUMSBW !*********************************************************************** !>\ingroup cs_scheme !! This subroution calculates freeze, melt and evaporation in cumulus downdraft. SUBROUTINE CUMDWN & ! Freeze & Melt & Evaporation ( IM , IJSDIM, KMAX , NTR , ntrq, & !DD dimensions GTT , GTQ , GTU , GTV , & ! modified GMFLX , & ! modified GPRCP , GSNWP , GTEVP , GMDD , & ! output GPRCI , GSNWI , & ! input GDH , GDW , GDQ , GDQI , & ! input GDQS , GDS , GDHS , GDT , & ! input GDU , GDV , GDZ , & ! input GDZM , FDQS , DELP , & ! input DELPI , & sigmad, do_aw , do_awdd, flx_form, & !DDsigma input gtmelt, gtevap, gtsubl, & !DDsigma input dtdwn , dqvdwn, dqldwn, dqidwn, & !DDsigma input dtrdwn, & KB , KTMX , ISTS , IENS ) ! input ! ! DD AW : modify to get eddy fluxes and microphysical tendencies for AW ! IMPLICIT NONE INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR, ntrq ! DD, for GFS, pass in logical, intent(in) :: do_aw, do_awdd, flx_form ! ! [MODIFY] REAL(kind_phys) GTT (IJSDIM, KMAX) !< Temperature tendency REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) !< Moisture etc tendency REAL(kind_phys) GTU (IJSDIM, KMAX) !< u tendency REAL(kind_phys) GTV (IJSDIM, KMAX) !< v tendency REAL(kind_phys) GMFLX (IJSDIM, KMAX) !< mass flux ! ! [OUTPUT] REAL(kind_phys) GPRCP (IJSDIM, KMAX) !< rainfall flux REAL(kind_phys) GSNWP (IJSDIM, KMAX) !< snowfall flux REAL(kind_phys) GTEVP (IJSDIM, KMAX) !< evaporation+sublimation REAL(kind_phys) GMDD (IJSDIM, KMAX) !< downdraft mass flux !AW microphysical tendencies REAL(kind_phys) gtmelt (IJSDIM, KMAX) !< t tendency ice-liq REAL(kind_phys) gtevap (IJSDIM, KMAX) !< t tendency liq-vapor REAL(kind_phys) gtsubl (IJSDIM, KMAX) !< t tendency ice-vapor !AW eddy flux tendencies REAL(kind_phys) dtdwn (IJSDIM, KMAX) !< t tendency downdraft detrainment REAL(kind_phys) dqvdwn (IJSDIM, KMAX) !< qv tendency downdraft detrainment REAL(kind_phys) dqldwn (IJSDIM, KMAX) !< ql tendency downdraft detrainment REAL(kind_phys) dqidwn (IJSDIM, KMAX) !< qi tendency downdraft detrainment REAL(kind_phys) dtrdwn (IJSDIM, KMAX, ntrq:ntr) !< tracer tendency downdraft detrainment ! AW downdraft area fraction (assumed zero for now) REAL(kind_phys) sigmad (IJSDIM,KMAX) !< DDsigma cloud downdraft area fraction ! [INPUT] REAL(kind_phys) GPRCI (IJSDIM, KMAX) !< rainfall generation REAL(kind_phys) GSNWI (IJSDIM, KMAX) !< snowfall generation REAL(kind_phys) GDH (IJSDIM, KMAX) !< moist static energy REAL(kind_phys) GDW (IJSDIM, KMAX) !< total water REAL(kind_phys) GDQ (IJSDIM, KMAX, NTR) !< humidity etc REAL(kind_phys) GDQI (IJSDIM, KMAX) !< cloud ice REAL(kind_phys) GDQS (IJSDIM, KMAX) !< saturate humidity REAL(kind_phys) GDS (IJSDIM, KMAX) !< dry static energy REAL(kind_phys) GDHS (IJSDIM, KMAX) !< saturate moist static energy REAL(kind_phys) GDT (IJSDIM, KMAX) !< air temperature T REAL(kind_phys) GDU (IJSDIM, KMAX) !< u-velocity REAL(kind_phys) GDV (IJSDIM, KMAX) !< v-velocity REAL(kind_phys) GDZ (IJSDIM, KMAX) !< altitude REAL(kind_phys) GDZM (IJSDIM, KMAX+1) !< altitude (half lev) REAL(kind_phys) FDQS (IJSDIM, KMAX) REAL(kind_phys) DELP (IJSDIM, KMAX) REAL(kind_phys) DELPI (IJSDIM, KMAX) INTEGER KB (IJSDIM) INTEGER KTMX, ISTS, IENS ! ! [INTERNAL WORK] ! Note: Some variables have 3-dimensions for the purpose of budget check. REAL(kind_phys) EVAPD (IJSDIM, KMAX) !< evap. in downdraft REAL(kind_phys) SUBLD (IJSDIM, KMAX) !< subl. in downdraft REAL(kind_phys) EVAPE (IJSDIM, KMAX) !< evap. in environment REAL(kind_phys) SUBLE (IJSDIM, KMAX) !< subl. in environment REAL(kind_phys) EVAPX (IJSDIM, KMAX) !< evap. env. to DD REAL(kind_phys) SUBLX (IJSDIM, KMAX) !< subl. env. to DD REAL(kind_phys) GMDDE (IJSDIM, KMAX) !< downdraft entrainment REAL(kind_phys) SNMLT (IJSDIM, KMAX) !< melt - freeze REAL(kind_phys) GCHDD (IJSDIM, KMAX) !< MSE detrainment REAL(kind_phys) GCWDD (IJSDIM, KMAX) !< water detrainment REAL(kind_phys) GTTEV (IJSDIM, KMAX) !< T tendency by evaporation REAL(kind_phys) GTQEV (IJSDIM, KMAX) !< q tendency by evaporation REAL(kind_phys) GCHD (ISTS:IENS) !< downdraft MSE REAL(kind_phys) GCWD (ISTS:IENS) !< downdraft q ! profiles of downdraft variables for AW flux tendencies REAL(kind_phys) GCdseD(ISTS:IENS, KMAX) !< downdraft dse REAL(kind_phys) GCqvD (ISTS:IENS, KMAX) !< downdraft qv ! REAL(kind_phys) GCqlD (ISTS:IENS, KMAX) !< downdraft ql ! REAL(kind_phys) GCqiD (ISTS:IENS, KMAX) !< downdraft qi REAL(kind_phys) GCtrD (ISTS:IENS, ntrq:ntr) !< downdraft tracer REAL(kind_phys) GCUD (ISTS:IENS) !< downdraft u REAL(kind_phys) GCVD (ISTS:IENS) !< downdraft v REAL(kind_phys) FSNOW (ISTS:IENS) REAL(kind_phys) GMDDD (ISTS:IENS) REAL(kind_phys) GDTW, GCHX, GCTX, GCQSX, GTPRP, EVSU, GTEVE, LVIC, & DQW, DTW, GDQW, DZ, GCSD, FDET, GDHI, GMDDX, & GMDDMX, GCHDX, GCWDX, GCUDD, GCVDD, GTHCI, GTQVCI, & wrk, wrk1, wrk2, wrk3, wrk4, tx1, & WMX, HMX, DDWMX, DDHMX, dp_above, dp_below, fsigma, & fmelt, fevp, wrkn, gctrdd(ntrq:ntr) !M REAL(kind_phys) GTHCI, GTQVCI, GTQLCI, GTQICI, GTUCI, GTVCI !DD#ifdef OPT_CUMBGT ! Water, energy, downdraft water and downdraft energy budgets ! REAL(kind_phys), dimension(ISTS:IENS) :: WBGT, HBGT, DDWBGT, DDHBGT integer ij, i, k, kp1, n !DD#endif ! ! [INTERNAL PARM] REAL(kind_phys), parameter :: TWSNOW = 273.15_kind_phys !< wet-bulb temp. rain/snow REAL(kind_phys), parameter :: FTMLT = 4._kind_phys !< temp. factor for melt REAL(kind_phys), parameter :: GMFLXC = 5.e-2_kind_phys !< critical mass flux REAL(kind_phys), parameter :: VTERMS = 2._kind_phys !< terminal velocity of snowflake ! REAL(kind_phys), parameter :: MELTAU = 10._kind_phys !< melting timescale REAL(kind_phys), parameter :: MELTAU = 20._kind_phys !< melting timescale ! Moorthi june 30, 2017 ! ! REAL(kind_phys), parameter :: EVAPR = 0.4_kind_phys !< evaporation factor ! Moorthi June 28, 2017 REAL(kind_phys), parameter :: EVAPR = 0.3_kind_phys !< evaporation factor ! REAL(kind_phys), parameter :: EVAPR = 0._kind_phys !< evaporation factor REAL(kind_phys), parameter :: REVPDD = 1._kind_phys !< max rate of DD to evapolation REAL(kind_phys), parameter :: RDDR = 5.e-4_kind_phys !< DD rate (T0 R0 W0)^-1 ! REAL(kind_phys), parameter :: RDDR = 0._kind_phys !< DD rate (T0 R0 W0)^-1 REAL(kind_phys), parameter :: RDDMX = 0.5_kind_phys !< norm. flux of downdraft REAL(kind_phys), parameter :: VTERM = 5._kind_phys !< term. vel. of precip. ! REAL(kind_phys), parameter :: VTERM = 4._kind_phys !< term. vel. of precip. ! Moorthi June 28, 2017 REAL(kind_phys), parameter :: EVATAU = 2._kind_phys !< evaporation/sublimation timescale REAL(kind_phys), parameter :: ZDMIN = 5.e2_kind_phys !< min altitude of downdraft detrainment real(kind_phys), parameter :: evapovtrm=EVAPR/VTERM !NOTE ! downdraft area ffraction still needs to be computed for AW, assumed zero for now, ! as passed in. ! ! Note: It is assumed that condensate evaporates in downdraft air. ! do k=1,kmax do i=ists,iens GPRCP (I,k) = zero GSNWP (I,k) = zero GMDD (I,k) = zero GTEVP (I,k) = zero EVAPD (I,k) = zero SUBLD (I,k) = zero EVAPE (I,k) = zero SUBLE (I,k) = zero EVAPX (I,k) = zero SUBLX (I,k) = zero GMDDE (I,k) = zero SNMLT (I,k) = zero GCHDD (I,k) = zero GCWDD (I,k) = zero GTTEV (I,k) = zero GTQEV (I,k) = zero GCdseD(I,k) = zero GCqvD (I,k) = zero ! GCqlD (I,k) = zero ! GCqiD (I,k) = zero gtevap(I,k) = zero gtmelt(I,k) = zero gtsubl(I,k) = zero enddo enddo ! testing on oct 17 2016 if (flx_form) then if (.not. do_awdd) then do k=1,kmax do i=ists,iens if (kb(i) > 0) then dtdwn (i,k) = gtt(i,k) dqvdwn(i,k) = gtq(i,k,1) dqldwn(i,k) = gtq(i,k,itl) dqidwn(i,k) = gtq(i,k,iti) endif enddo enddo do n=ntrq,ntr do k=1,kmax do i=ists,iens if (kb(i) > 0) then dtrdwn(i,k,n) = gtq(i,k,n) endif enddo enddo enddo else do k=1,kmax do i=ists,iens dtdwn (I,k) = zero dqvdwn(I,k) = zero dqldwn(I,k) = zero dqidwn(I,k) = zero enddo enddo do n=ntrq,ntr do k=1,kmax do i=ists,iens dtrdwn(i,k,n) = zero enddo enddo enddo endif endif ! do i=ists,iens GCHD(I) = zero GCWD(I) = zero GCUD(I) = zero GCVD(I) = zero enddo do n=ntrq,ntr do i=ists,iens GCtrD (I,n) = zero enddo enddo ! DO K=KTMX,1,-1 ! loop A kp1 = min(k+1,kmax) ! ! < precipitation melt & freeze > ! DO I=ISTS,IENS if (kb(i) > 0) then GTPRP = GPRCP(I,KP1) + GSNWP(I,KP1) IF (GTPRP > zero) THEN FSNOW(I) = GSNWP(I,KP1) / GTPRP ELSE FSNOW(I) = zero ENDIF LVIC = ELocp + EMELTocp*FSNOW(I) GDTW = GDT(I,K) - LVIC*(GDQS(I,K) - GDQ(I,K,1)) & / (one + LVIC*FDQS(I,K)) IF (GDTW < TWSNOW) THEN GSNWP(I,K) = GSNWP(I,KP1) + GPRCI(I,K) + GSNWI(I,K) GTTEV(I,K) = EMELToCP * GPRCI(I,K) * DELPI(I,K) SNMLT(I,K) = -GPRCI(I,K) ELSE DZ = GDZM(I,KP1) - GDZM(I,K) FMELT = (one + FTMLT*(GDTW - TWSNOW)) & * (one - TANH(GMFLX(I,KP1)/GMFLXC)) & * (one - TANH(VTERMS*MELTAU/DZ)) SNMLT(I,K) = GSNWP(I,KP1)*max(min(FMELT, one), zero) GSNWP(I,K) = GSNWP(I,KP1)+GSNWI(I,K) - SNMLT(I,K) GPRCP(I,K) = GPRCP(I,KP1)+GPRCI(I,K) + SNMLT(I,K) GTTEV(I,K) = -EMELToCP * SNMLT(I,K) * DELPI(I,K) ENDIF !DD heating rate due to precip melting for AW gtmelt(i,k) = gtmelt(i,k) + GTTEV(I,K) endif ENDDO ! ! < downdraft > ! DO I=ISTS,IENS ! loop B if (kb(i) > 0) then wrk = delpi(i,k) wrk1 = oneocp * wrk DZ = GDZM(I,KP1) - GDZM(I,K) FEVP = (one - TANH(EVATAU*VTERM/DZ)) IF (GMDD(I,KP1) > zero) THEN GCHX = GCHD(I) / GMDD(I,KP1) GCTX = GDT(I,K) + (GCHX-GDHS(I,K)) / (CP+EL*FDQS(I,K)) GCQSX = GDQS(I,K) + FDQS(I,K) * (GCTX - GDT(I,K)) GCQSX = GCQSX*GMDD(I,KP1) EVSU = MAX(GCQSX-GCWD(I), zero) * FEVP GTPRP = GPRCP(I,K) + GSNWP(I,K) IF (GTPRP > zero) THEN FSNOW(I) = GSNWP(I,K) / GTPRP ELSE FSNOW(I) = zero ENDIF EVAPD(I,K) = min(EVSU*(one-FSNOW(I)), GPRCP(I,K)) SUBLD(I,K) = min(EVSU*FSNOW(I), GSNWP(I,K)) GPRCP(I,K) = GPRCP(I,K) - EVAPD(I,K) GSNWP(I,K) = GSNWP(I,K) - SUBLD(I,K) ! temperature tendencies due to evaporation and sublimation of precip ! This is within downdraft gtevap(i,k) = gtevap(i,k) - elocp * evapd(i,k) * wrk gtsubl(i,k) = gtsubl(i,k) - esubocp * subld(i,k) * wrk GCWD(I) = GCWD(I) + EVAPD(I,K) + SUBLD(I,K) GCHD(I) = GCHD(I) - EMELT*SUBLD(I,K) ENDIF GMDD(I,K) = GMDD(I,KP1) ! LVIC = ELocp + EMELTocp*FSNOW(I) DQW = (GDQS(I,K) - GDW(I,K)) / (one + LVIC*FDQS(I,K)) DQW = MAX(DQW, zero) DTW = LVIC*DQW GDQW = GDW(I,K) + DQW*FEVP ! EVSU = min(one, EVAPOVTRM*DQW*DZ*FEVP) EVAPE(I,K) = EVSU*GPRCP(I,K) SUBLE(I,K) = EVSU*GSNWP(I,K) GTEVP(I,K) = EVAPD(I,K) + SUBLD(I,K) + EVAPE(I,K) + SUBLE(I,K) ! GTPRP = GPRCP(I,K) + GSNWP(I,K) GPRCP(I,K) = GPRCP(I,K) - EVAPE(I,K) GSNWP(I,K) = GSNWP(I,K) - SUBLE(I,K) ! additional temperature tendencies due to evaporation and sublimation of precip ! This is outside of downdraft gtevap(i,k) = gtevap(i,k) - el*evape(i,k) * wrk1 gtsubl(i,k) = gtsubl(i,k) - (el+emelt)*suble(i,k) * wrk1 ! GMDDD(I) = zero IF (GDZ(I,K)-GDZM(I,1) > ZDMIN) THEN GTEVE = EVAPE(I,K) + SUBLE(I,K) GMDDMX = REVPDD*GTEVE/MAX(DQW, 1.D-10) GMDDE(I,K) = RDDR * (DTW*GTPRP*DELP(I,K)) GMDDE(I,K) = MAX(MIN(GMDDE(I,K), GMDDMX), zero) GMDDX = GMDD(I,KP1) + GMDDE(I,K) EVSU = GMDDE(I,K)*DQW*FEVP IF (GTEVE > zero) THEN FSNOW(I) = SUBLE(I,K) / GTEVE ELSE FSNOW(I) = zero END IF EVAPX(I,K) = (one-FSNOW(I)) * EVSU SUBLX(I,K) = FSNOW(I) * EVSU ! IF (GMDDX > zero) THEN GDHI = GDH(I,K) - EMELT*GDQI(I,K) GCHDX = GCHD(I) + GDHI*GMDDE(I,K) - EMELT*SUBLX(I,K) GCWDX = GCWD(I) + GDQW*GMDDE(I,K) GCSD = (GCHDX - EL*GCWDX) / GMDDX IF (GCSD < GDS(I,K)) THEN GCHD(I) = GCHDX GCWD(I) = GCWDX GCUD(I) = GCUD(I) + GDU(I,K)*GMDDE(I,K) GCVD(I) = GCVD(I) + GDV(I,K)*GMDDE(I,K) do n = ntrq,ntr GCtrD(I,n) = GCtrD(I,n) + GDq(I,K,n)*GMDDE(I,K) enddo GMDD(I,K) = GMDDX EVAPE(I,K) = EVAPE(I,K) - EVAPX(I,K) SUBLE(I,K) = SUBLE(I,K) - SUBLX(I,K) EVAPD(I,K) = EVAPD(I,K) + EVAPX(I,K) SUBLD(I,K) = SUBLD(I,K) + SUBLX(I,K) GMDDD(I) = zero ELSE GMDDE(I,K) = zero GMDDD(I) = GMDD(I,KP1) ENDIF ENDIF ELSE GMDDD(I) = DZ / (GDZM(I,KP1)-GDZM(I,1)) * GMDD(I,KP1) ENDIF ! GMDDD(I) = MAX(GMDDD(I), GMDD(I,K)-RDDMX*GMFLX(I,K)) ! IF (GMDDD(I) > zero) THEN FDET = GMDDD(I)/GMDD(I,K) GCHDD(I,K) = FDET*GCHD(I) GCWDD(I,K) = FDET*GCWD(I) GCUDD = FDET*GCUD(I) GCVDD = FDET*GCVD(I) do n = ntrq,ntr GCtrDD(n) = FDET*GCtrD(I,n) enddo ! GTHCI = wrk * (GCHDD(I,K) - GMDDD(I)*GDH(I,K)) GTQVCI = wrk * (GCWDD(I,K) - GMDDD(I)*GDQ(I,K,1)) ! GTT (I,K) = GTT(I,K) + (GTHCI - EL*GTQVCI)*oneoCP GTQ (I,K,1) = GTQ(I,K,1) + GTQVCI GTQ (I,K,ITL) = GTQ(I,K,ITL) - wrk * GMDDD(I)*GDQ(I,K,ITL) GTQ (I,K,ITI) = GTQ(I,K,ITI) - wrk * GMDDD(I)*GDQI(I,K) do n = ntrq,ntr GTQ (I,K,n) = GTQ(I,K,n) + wrk * (GCtrDD(n) - GMDDD(I)*GDQ(I,K,n)) GCtrD(I,n) = GCtrD(I,n) - GCtrDD(n) enddo GTU (I,K) = GTU(I,K) + wrk * (GCUDD - GMDDD(I)*GDU(I,K)) GTV (I,K) = GTV(I,K) + wrk * (GCVDD - GMDDD(I)*GDV(I,K)) ! GCHD(I) = GCHD(I) - GCHDD(I,K) GCWD(I) = GCWD(I) - GCWDD(I,K) GCUD(I) = GCUD(I) - GCUDD GCVD(I) = GCVD(I) - GCVDD GMDD(I,K) = GMDD(I,K) - GMDDD(I) ENDIF GCdseD(I,K) = GCHD(I) - el*GCWD(I) GCqvD (I,K) = GCWD(I) endif ENDDO ! loop B ! ENDDO ! loop A ! DO K=1,KTMX kp1 = min(k+1,kmax) DO I=ISTS,IENS if (kb(i) > 0) then wrk = DELPI(I,k) tx1 = DELPI(I,kp1) GTTEV(I,K) = GTTEV(I,K) - wrk & * (ELocp*EVAPE(I,K)+(ELocp+EMELTocp)*SUBLE(I,K)) GTT(I,K) = GTT(I,K) + GTTEV(I,K) ! GTQEV(I,K) = GTQEV(I,K) + (EVAPE(I,K)+SUBLE(I,K)) * wrk GTQ(I,K,1) = GTQ(I,K,1) + GTQEV(I,K) ! GMFLX(I,K) = GMFLX(I,K) - GMDD(I,K) ! AW tendencies due to vertical divergence of eddy fluxes if (k > 1 .and. flx_form) then fsigma = one - sigmad(i,kp1) dp_below = wrk * (one - sigmad(i,k)) dp_above = tx1 * (one - sigmad(i,kp1)) wrk1 = gmdd(i,kp1) * (gdt(i,k)+gocp*gdz(i,k)) - gcdsed(i,kp1)*oneocp wrk2 = gmdd(i,kp1) * gdq(i,k,1) - gcqvd(i,kp1) wrk3 = gmdd(i,kp1) * gdq(i,k,itl) wrk4 = gmdd(i,kp1) * gdqi(i,k) dtdwn(i,k) = dtdwn(i,k) + dp_below * wrk1 dqvdwn(i,k) = dqvdwn(i,k) + dp_below * wrk2 dqldwn(i,k) = dqldwn(i,k) + dp_below * wrk3 ! gcqld=0 - gcqld(i,k)) dqidwn(i,k) = dqidwn(i,k) + dp_below * wrk4 ! gcqid=0 - gcqid(i,k)) dtdwn(i,kp1) = dtdwn(i,kp1) - dp_above * wrk1 dqvdwn(i,kp1) = dqvdwn(i,kp1) - dp_above * wrk2 dqldwn(i,kp1) = dqldwn(i,kp1) - dp_above * wrk3 ! gcqld=0 - gcqld(i,k)) dqidwn(i,kp1) = dqidwn(i,kp1) - dp_above * wrk4 ! gcqid=0 - gcqid(i,k)) do n = ntrq,ntr wrkn = gmdd(i,kp1) * gdq(i,k,n) dtrdwn(i,k,n) = dtrdwn(i,k,n) + dp_below * wrkn dtrdwn(i,kp1,n) = dtrdwn(i,kp1,n) - dp_above * wrkn enddo endif endif ENDDO ! end of i loop ENDDO ! end of k loop ! if (.not. do_awdd .and. flx_form) then do k=1,kmax do i=ists,iens if (kb(i) > 0) then dtdwn(i,k) = gtt(i,k) - dtdwn(i,k) dqvdwn(i,k) = gtq(i,k,1) - dqvdwn(i,k) dqldwn(i,k) = gtq(i,k,itl) - dqldwn(i,k) dqidwn(i,k) = gtq(i,k,iti) - dqidwn(i,k) endif enddo enddo do n=ntrq,ntr do k=1,kmax do i=ists,iens if (kb(i) > 0) then dtrdwn(i,k,n) = gtq(i,k,n) - dtrdwn(i,k,n) endif enddo enddo enddo endif ! END SUBROUTINE CUMDWN !*********************************************************************** !>\ingroup cs_scheme !> This subroutine computes cumulus cloudiness. SUBROUTINE CUMCLD & !! cloudiness ( IJSDIM, KMAX , & !DD dimensions CUMCLW, QLIQ , QICE , FLIQC , & ! modified CUMFRC, & ! output GMFLX , KTMX , ISTS, IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, KMAX ! DD, for GFS, pass in ! ! [OUTPUT] REAL(kind_phys) CUMFRC(IJSDIM) !< cumulus cloud fraction ! ! [MODIFY] REAL(kind_phys) CUMCLW(IJSDIM, KMAX) !< cloud water in cumulus REAL(kind_phys) QLIQ (IJSDIM, KMAX) !< cloud liquid REAL(kind_phys) QICE (IJSDIM, KMAX) !< cloud ice REAL(kind_phys) FLIQC (IJSDIM, KMAX) !< liquid ratio in cumulus ! ! [INPUT] REAL(kind_phys) GMFLX (IJSDIM, KMAX) ! cumulus mass flux INTEGER KTMX INTEGER ISTS, IENS ! ! [WORK] INTEGER I, K REAL(kind_phys) CUMF, QC, wrk ! ! [INTERNAL PARAM] REAL(kind_phys), parameter :: CMFMIN = 2.e-3_kind_phys, &!< Mc->cloudiness CMFMAX = 3.e-1_kind_phys, &!< Mc->cloudiness CLMIN = 1.e-3_kind_phys, &!< cloudiness Min. CLMAX = 0.1_kind_phys, &!< cloudiness Max. FACLW = 0.1_kind_phys, &!< Mc->CLW FACLF = (CLMAX-CLMIN)/LOG(CMFMAX/CMFMIN) ! CUMFRC(ISTS:IENS) = zero DO K=1,KTMX DO I=ISTS,IENS CUMFRC(I) = MAX(CUMFRC(I), GMFLX(I,K)) ENDDO ENDDO DO I=ISTS,IENS IF (CUMFRC(I) > zero) THEN CUMF = LOG(MAX(CUMFRC(I), CMFMIN)/CMFMIN) CUMFRC(I) = MIN(FACLF*CUMF+CLMIN, CLMAX) ENDIF ENDDO ! DO K=1,KTMX DO I=ISTS,IENS IF (GMFLX(I,K) > zero) THEN wrk = FACLW / GMFLX(I,K) * CUMFRC(I) QLIQ (I,K) = wrk * QLIQ(I,K) QICE (I,K) = wrk * QICE(I,K) CUMCLW(I,K) = wrk * CUMCLW(I,K) QC = QLIQ(I,K) + QICE(I,K) IF (QC > zero) THEN FLIQC(I,K) = QLIQ(I,K) / QC ENDIF ENDIF ENDDO ENDDO ! END SUBROUTINE CUMCLD !*********************************************************************** !>\ingroup cs_scheme !! This subroutine calculates SUBROUTINE CUMUPR & !! Tracer Updraft ( im , IJSDIM, KMAX , NTR , & !DD dimensions GTR , GPRCC , & ! modified GDR , CBMFX , & ! input GCYM , GCYT , GCQT , GCLT , GCIT , & ! input GTPRT , GTEVP , GTPRC0, & ! input KB , KBMX , KT , KTMX , KTMXT , & ! input DELPI , OTSPT , ISTS , IENS, & ! input fscav, fswtr, nctp) ! IMPLICIT NONE INTEGER, INTENT(IN) :: im, IJSDIM, KMAX, NTR, nctp !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) REAL(kind_phys) GPRCC (IJSDIM, NTR) ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) REAL(kind_phys) CBMFX (IM, NCTP) REAL(kind_phys) GCYM (IJSDIM, KMAX, nctp) REAL(kind_phys) GCYT (IJSDIM, NCTP) REAL(kind_phys) GCQT (IJSDIM, NCTP) REAL(kind_phys) GCLT (IJSDIM, NCTP) REAL(kind_phys) GCIT (IJSDIM, NCTP) REAL(kind_phys) GTPRT (IJSDIM, NCTP) REAL(kind_phys) GTEVP (IJSDIM, KMAX) REAL(kind_phys) GTPRC0(IJSDIM) !! precip. before evap. real(kind_phys) fscav(ntr), fswtr(ntr) INTEGER KB (IJSDIM ) INTEGER KBMX INTEGER KT (IJSDIM, NCTP) INTEGER KTMX (NCTP) INTEGER KTMXT REAL(kind_phys) DELPI (IJSDIM, KMAX) LOGICAL OTSPT (NTR) !! transport with this routine? INTEGER ISTS, IENS ! ! [INTERNAL WORK] INTEGER I, K, LT, TP, CTP REAL(kind_phys) :: GCRTD, SCAV, GCWT, GPRCR, evpf, cbmfxl REAL(kind_phys), dimension(ists:iens) :: GCRB, GCRT, DR, gtprc0i ! REAL(kind_phys), dimension(ists:iens,kmax) :: DGCB, DZ, RDZM, EVPF ! REAL(kind_phys), dimension(ists:iens,nctp) :: DZT, RGCWT, MASK1, MASK2 REAL(kind_phys), dimension(ists:iens,nctp) :: RGCWT, MASK1 ! do i=ists,iens if (gtprc0(i) > zero) then gtprc0i(i) = one / gtprc0(i) else gtprc0i(i) = zero endif enddo DO CTP=1,NCTP DO I=ISTS,IENS K = KT(I,CTP) ! GCWT = GCQT(I,CTP) + GCLT(I,CTP) + GCIT(I,CTP) RGCWT(I,CTP) = zero IF (GCWT > zero) THEN RGCWT(I,CTP) = one / GCWT ENDIF ! MASK1(I,CTP) = zero IF (kb(i) > 0 .and. K > KB(I)) THEN MASK1(I,CTP) = one ENDIF ! MASK2(I,CTP) = zero ! IF (CBMFX(I,CTP) > zero) then ! MASK2(I,CTP) = one ! ENDIF ENDDO ENDDO ! DO LT=1,NTR ! outermost tracer LT loop ! IF (OTSPT(LT)) THEN DO CTP=1,NCTP DO I=ISTS,IENS GCRB(i) = zero DR(i) = zero enddo DO K=1,KBMX DO I=ISTS,IENS IF (kb(i) > 0 .and. K < KB(I)) THEN GCRB(I) = GCRB(I) + (GCYM(I,K+1,ctp)-GCYM(I,K,ctp))* GDR(I,K,LT) ENDIF ENDDO ENDDO ! DO K=2,KTMX(CTP) DO I=ISTS,IENS IF (kb(i) > 0 .and. K >= KB(I) .AND. K < KT(I,CTP)) THEN DR(I) = DR(I) + (GCYM(I,K+1,ctp)-GCYM(I,K,ctp)) * GDR(I,K,LT) ENDIF ENDDO ENDDO ! DO I=ISTS,IENS K = KT(I,CTP) if (kb(i) > 0 .and. k > kb(i)) then DR(I) = DR(I) + (GCYT(I,CTP) - GCYM(I,K,ctp)) * GDR (I,K,LT) & * MASK1(I,CTP) GCRT(I) = (GCRB(I) + DR(I)) * MASK1(I,CTP) ! SCAV = FSCAV(LT)*GTPRT(I,CTP) + FSWTR(LT)*GTPRT(I,CTP)*RGCWT(I,CTP) SCAV = MIN(SCAV, one) GCRTD = GCRT(I) * (one - SCAV) cbmfxl = max(zero, CBMFX(I,CTP)) GPRCR = SCAV * GCRT(I) * CBMFXl GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * CBMFXl & * (GCRTD - GCYT(I,CTP) * GDR(I,K,LT)) GPRCC(I,LT) = GPRCC(I,LT) + GPRCR ! GPRCR = SCAV * GCRT(I) * CBMFX(I,CTP) ! GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * CBMFX(I,CTP) & ! * (GCRTD - GCYT(I,CTP) * GDR(I,K,LT)) * MASK2(I,CTP) ! GPRCC(I,LT) = GPRCC(I,LT) + GPRCR * MASK2(I,CTP) endif ENDDO ENDDO ! DO K=KTMXT,1,-1 DO I=ISTS,IENS evpf = GTEVP(i,k) * gtprc0i(i) GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * GPRCC(I,LT) * EVPF GPRCC(I,LT) = GPRCC(I,LT) * (one - EVPF) ! GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) * GPRCC(I,LT) * EVPF(I,K) ! GPRCC(I,LT) = GPRCC(I,LT) * (one - EVPF(I,K)) ENDDO ENDDO ! ENDIF ! ENDDO ! outermost tracer LT loop ! END SUBROUTINE CUMUPR !*********************************************************************** !>\ingroup cs_scheme SUBROUTINE CUMDNR & !! Tracer Downdraft ( IM , IJSDIM, KMAX , NTR , & !DD dimensions GTR , & ! modified GDR , GMDD , DELPI , & ! input KTMX , OTSPT , ISTS , IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! Temperature tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) REAL(kind_phys) GMDD (IJSDIM, KMAX) ! downdraft mass flux REAL(kind_phys) DELPI (IJSDIM, KMAX ) LOGICAL OTSPT (NTR) INTEGER KTMX, ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) GCRD (ISTS:IENS) ! downdraft q REAL(kind_phys) GMDDE, GMDDD, GCRDD INTEGER I, K, LT, kp1 ! ! DO LT=1,NTR IF (OTSPT(LT)) THEN GCRD = zero DO K=KTMX,1,-1 kp1 = min(k+1,kmax) DO I=ISTS,IENS GMDDE = GMDD(I,K) - GMDD(I,KP1) IF (GMDDE >= zero) THEN GCRD(I) = GCRD(I) + GDR(I,K,LT)*GMDDE ELSEIF (GMDD(I,KP1) > zero) THEN GMDDD = - GMDDE GCRDD = GMDDD/GMDD(I,KP1) * GCRD(I) GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) & * (GCRDD - GMDDD*GDR(I,K,LT)) GCRD(I) = GCRD(I) - GCRDD ENDIF ENDDO ENDDO ENDIF ENDDO ! END SUBROUTINE CUMDNR !*********************************************************************** !>\ingroup cs_scheme SUBROUTINE CUMSBR & !! Tracer Subsidence ( IM , IJSDIM, KMAX , NTR , & !DD dimensions GTR , & ! modified GDR , DELPI , & ! input GMFLX , KTMX , OTSPT , & ! input ISTS, IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) !! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) !! tracer REAL(kind_phys) DELPI (IJSDIM, KMAX) REAL(kind_phys) GMFLX (IJSDIM, KMAX) !! mass flux INTEGER KTMX LOGICAL OTSPT (NTR) !! tracer transport on/off INTEGER ISTS, IENS ! ! [INTERNAL WORK] INTEGER I, K, KM, KP, LT REAL(kind_phys) SBR0, SBR1, FX1 REAL(kind_phys) FX(ISTS:IENS) ! DO LT=1,NTR IF (OTSPT(LT)) THEN DO I=ISTS,IENS FX(I) = zero enddo DO K=KTMX,1,-1 KM = MAX(K-1, 1) KP = MIN(K+1, KMAX) DO I=ISTS,IENS SBR0 = GMFLX(I,KP) * (GDR(I,KP,LT) - GDR(I,K,LT)) SBR1 = GMFLX(I,K) * (GDR(I,K,LT) - GDR(I,KM,LT)) IF (GMFLX(I,K) > GMFLX(I,KP)) THEN FX1 = half ELSE FX1 = zero END IF GTR(I,K,LT) = GTR(I,K,LT) + DELPI(I,K) & * ((one-FX(I))*SBR0 + FX1*SBR1) FX(I) = FX1 ENDDO ENDDO ENDIF ENDDO ! END SUBROUTINE CUMSBR !********************************************************************* !>\ingroup cs_scheme !! This subroutine calculates tracer mass fixer without detrainment. SUBROUTINE CUMFXR & ! Tracer mass fixer ( IM , IJSDIM, KMAX , NTR , & !DD dimensions GTR , & ! modified GDR , DELP , DELTA , KTMX , IMFXR , & ! input ISTS , IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX, NTR !! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX, NTR) ! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX, NTR) ! tracer REAL(kind_phys) DELP (IJSDIM, KMAX) REAL(kind_phys) DELTA ! time step INTEGER KTMX INTEGER IMFXR (NTR) ! 0: mass fixer is not applied ! tracers which may become negative values ! e.g. subgrid-PDFs ! 1: mass fixer is applied, total mass may change through cumulus scheme ! e.g. moisture, liquid cloud, ice cloud, aerosols ! 2: mass fixer is applied, total mass never change through cumulus scheme ! e.g. CO2 !DD add new CASE ! 3: just fill holes, no attempt to conserve INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) GDR1 REAL(kind_phys) GDR2 (ISTS:IENS, KMAX) REAL(kind_phys), dimension(ISTS:IENS) :: TOT0, TOT1, TRAT REAL(kind_phys) FWAT INTEGER I, K, LT ! ! Attention: tracers are forced to be positive unless IMFXR=0. ! DO LT=1,NTR SELECT CASE (IMFXR(LT)) CASE (0) CYCLE CASE (1) FWAT = one CASE (2) FWAT = zero CASE (3) CASE DEFAULT EXIT END SELECT ! DO I=ISTS,IENS TOT0(I) = zero TOT1(I) = zero TRAT(I) = one enddo ! DO K=KTMX,1,-1 DO I=ISTS,IENS IF (GTR(I,K,LT) /= zero) THEN GDR1 = GDR(I,K,LT) + DELTA*GTR(I,K,LT) GDR2(I,K) = MAX(GDR1, zero) GDR1 = GDR1 * FWAT + GDR(I,K,LT)*(one - FWAT) TOT0(I) = TOT0(I) + GDR1 *(DELP(I,K)*GRAVI) TOT1(I) = TOT1(I) + GDR2(I,K)*(DELP(I,K)*GRAVI) ENDIF ENDDO ENDDO ! if(imfxr(LT) .ne. 3) then DO I=ISTS,IENS IF (TOT1(I) > zero ) THEN TRAT(I) = MAX(TOT0(I), zero) / TOT1(I) ENDIF ENDDO endif ! DO K=KTMX,1,-1 DO I=ISTS,IENS IF (GTR(I,K,LT) /= zero ) THEN GDR2(I,K ) = GDR2(I,K)*TRAT(I) GTR (I,K,LT) = (GDR2(I,K)-GDR(I,K,LT)) / DELTA ENDIF ENDDO ENDDO ! ENDDO ! LT-loop ! END SUBROUTINE CUMFXR !********************************************************************* !>\ingroup cs_scheme SUBROUTINE CUMFXR1 & ! Tracer mass fixer ( IM , IJSDIM, KMAX , & !DD dimensions GTR , & ! modified GDR , DELP , DELTA , KTMX , IMFXR , & ! input ISTS , IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IM, IJSDIM, KMAX ! DD, for GFS, pass in ! ! [MODIFY] REAL(kind_phys) GTR (IJSDIM, KMAX) ! tracer tendency ! ! [INPUT] REAL(kind_phys) GDR (IJSDIM, KMAX) ! tracer REAL(kind_phys) DELP (IJSDIM, KMAX) REAL(kind_phys) DELTA ! time step INTEGER KTMX INTEGER IMFXR ! 0: mass fixer is not applied ! tracers which may become negative values ! e.g. subgrid-PDFs ! 1: mass fixer is applied, total mass may change through cumulus scheme ! e.g. moisture, liquid cloud, ice cloud, aerosols ! 2: mass fixer is applied, total mass never change through cumulus scheme ! e.g. CO2 INTEGER ISTS, IENS ! ! [INTERNAL WORK] REAL(kind_phys) GDR1 REAL(kind_phys) GDR2 (ISTS:IENS, KMAX) REAL(kind_phys), dimension(ISTS:IENS) :: TOT0, TOT1, TRAT REAL(kind_phys) FWAT INTEGER I, K ! ! Attention: tracers are forced to be positive unless IMFXR=0. ! SELECT CASE (IMFXR) CASE (0) RETURN CASE (1) FWAT = one CASE (2) FWAT = zero CASE DEFAULT RETURN END SELECT ! DO I=ISTS,IENS TOT0(I) = zero TOT1(I) = zero enddo ! DO K=KTMX,1,-1 DO I=ISTS,IENS IF (GTR(I,K) /= zero) THEN GDR1 = GDR(I,K) + DELTA*GTR(I,K) GDR2(I,K) = MAX(GDR1, zero) GDR1 = GDR1*FWAT + GDR(I,K)*(one - FWAT) TOT0(I) = TOT0(I) + GDR1 *(DELP(I,K)*GRAVI) TOT1(I) = TOT1(I) + GDR2(I,K)*(DELP(I,K)*GRAVI) ENDIF ENDDO ENDDO ! DO I=ISTS,IENS IF (TOT1(I) > zero) THEN TRAT(I) = MAX(TOT0(I), zero) / TOT1(I) ELSE TRAT(I) = one ENDIF ENDDO ! DO K=KTMX,1,-1 DO I=ISTS,IENS IF (GTR(I,K) /= zero) THEN GDR2(I,K) = GDR2(I,K)*TRAT(I) GTR (I,K) = (GDR2(I,K)-GDR(I,K)) / DELTA ENDIF ENDDO ENDDO ! END SUBROUTINE CUMFXR1 !********************************************************************* !>\ingroup cs_scheme SUBROUTINE CUMCHK & ! check range of output values ( IJSDIM, KMAX , NTR , & !DD dimensions GTT , GTQ , GTU , GTV , & ! input GTCFRC, GPRCC , GSNWC , CUMCLW, & ! input CUMFRC, FLIQC , GTPRP , & ! input ISTS , IENS ) ! input ! IMPLICIT NONE INTEGER, INTENT(IN) :: IJSDIM, KMAX, NTR ! DD, for GFS, pass in ! ! [INPUT] REAL(kind_phys) GTT (IJSDIM, KMAX) ! heating rate REAL(kind_phys) GTQ (IJSDIM, KMAX, NTR) ! change in q REAL(kind_phys) GTU (IJSDIM, KMAX) ! tendency of u REAL(kind_phys) GTV (IJSDIM, KMAX) ! tendency of v REAL(kind_phys) GPRCC (IJSDIM, NTR ) ! rainfall REAL(kind_phys) GSNWC (IJSDIM) ! snowfall REAL(kind_phys) CUMCLW(IJSDIM, KMAX) ! cloud water in cumulus REAL(kind_phys) CUMFRC(IJSDIM) ! cumulus cloud fraction REAL(kind_phys) GTCFRC(IJSDIM, KMAX) ! change in cloud fraction REAL(kind_phys) FLIQC (IJSDIM, KMAX) ! liquid ratio in cumulus REAL(kind_phys) GTPRP (IJSDIM, KMAX) ! rain+snow flux ! INTEGER ISTS, IENS ! ! [INTERNAL WORK] INTEGER I, K ! ! [INTERNAL PARM] REAL(kind_phys) :: GTTMAX = 1.e-2_kind_phys REAL(kind_phys) :: GTQVMAX = 1.e-4_kind_phys REAL(kind_phys) :: GTQLMAX = 1.e-5_kind_phys REAL(kind_phys) :: GTUMAX = 1.e-2_kind_phys REAL(kind_phys) :: GTVMAX = 1.e-2_kind_phys REAL(kind_phys) :: GTCFMAX = 1.e-3_kind_phys REAL(kind_phys) :: PRCCMAX = 1.e-2_kind_phys REAL(kind_phys) :: SNWCMAX = 1.e-2_kind_phys REAL(kind_phys) :: CLWMAX = 1.e-3_kind_phys REAL(kind_phys) :: TPRPMAX = 1.e-2_kind_phys REAL(kind_phys) :: GTQIMAX = 1.e-5_kind_phys !REAL(kind_phys) :: GTM2MAX = 1._kind_phys !REAL(kind_phys) :: GTM3MAX = 1._kind_phys ! DO K=1,KMAX DO I=ISTS, IENS IF (ABS(GTT(I,K)) > GTTMAX) THEN WRITE(iulog,*) '### CUMCHK: GTT(',I,',',K,')=',GTT(I,K) ENDIF IF (ABS(GTQ(I,K,1) ) > GTQVMAX) THEN WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',1 )=', GTQ(I,K,1) ENDIF IF (ABS(GTQ(I,K,ITL)) > GTQLMAX) THEN WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',ITL )=', GTQ(I,K,ITL) ENDIF IF (ABS(GTU(I,K)) > GTUMAX) THEN WRITE(iulog,*) '### CUMCHK: GTU(',I,',',K,')=',GTU(I,K) END IF IF (ABS(GTV(I,K)) > GTVMAX) THEN WRITE(iulog,*) '### CUMCHK: GTV(',I,',',K,')=',GTV(I,K) ENDIF IF (ABS(GTCFRC(I,K)) > GTCFMAX) THEN WRITE(iulog,*) '### CUMCHK: GTCFRC(',I,',',K,')=', GTCFRC(I,K) ENDIF IF (CUMCLW(I,K) > CLWMAX .OR. CUMCLW(I,K) < zero) THEN WRITE(iulog,*) '### CUMCHK: CUMCLW(',I,',',K,')=', CUMCLW(I,K) ENDIF IF (FLIQC(I,K) > one .OR. FLIQC(I,K) < zero) THEN WRITE(iulog,*) '### CUMCHK: FLIQC(',I,',',K,')=', FLIQC(I,K) ENDIF IF (GTPRP(I,K) > TPRPMAX .OR. GTPRP(I,K) < zero) THEN WRITE(iulog,*) '### CUMCHK: GTPRP(',I,',',K,')=', GTPRP(I,K) ENDIF IF (ABS(GTQ(I,K,ITI)) > GTQIMAX) THEN WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',ITI )=', GTQ(I,K,ITI) ENDIF ! IF (ABS(GTQ(I,K,IMU2) ) > GTM2MAX) THEN ! WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',IMU2 )=', GTQ(I,K,IMU2) ! ENDIF ! IF (ABS(GTQ(I,K,IMU3)) > GTM3MAX) THEN ! WRITE(iulog,*) '### CUMCHK: GTQ(',I,',',K,',IMU3 )=', GTQ(I,K,IMU3) ! ENDIF ENDDO ENDDO ! DO I=ISTS,IENS IF (GPRCC(I,1) > PRCCMAX .OR. GPRCC(I,1) < zero) THEN WRITE(iulog,*) '### CUMCHK: GPRCC(',I,')=',GPRCC(I,1) END IF IF (GSNWC(I) > SNWCMAX .OR. GSNWC(I) < zero) THEN WRITE(iulog,*) '### CUMCHK: GSNWC(',I,')=',GSNWC(I) END IF IF (CUMFRC(I) > one .OR. CUMFRC(I) < zero) THEN WRITE(iulog,*) '### CUMCHK: CUMFRC(',I,')=',CUMFRC(I) ENDIF ENDDO ! END SUBROUTINE CUMCHK !*********************************************************************** end module cs_conv