!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! !! !! GNU General Public License !! !! !! !! This file is part of the Flexible Modeling System (FMS). !! !! !! !! FMS is free software; you can redistribute it and/or modify !! !! it and are expected to follow the terms of the GNU General Public !! !! License as published by the Free Software Foundation. !! !! !! !! FMS is distributed in the hope that it will be useful, !! !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! !! GNU General Public License for more details. !! !! !! !! You should have received a copy of the GNU General Public License !! !! along with FMS; if not, write to: !! !! Free Software Foundation, Inc. !! !! 59 Temple Place, Suite 330 !! !! Boston, MA 02111-1307 USA !! !! or see: !! !! http://www.gnu.org/licenses/gpl.txt !! !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module betts_miller_mod !---------------------------------------------------------------------- !use utilities_mod, only: file_exist, error_mesg, open_file, & ! check_nml_error, get_my_pe, FATAL, & ! close_file use fms_mod, only: file_exist, error_mesg, open_namelist_file, & check_nml_error, mpp_pe, mpp_root_pe, & FATAL, close_file, write_version_number, stdlog use sat_vapor_pres_mod, only: escomp, descomp use constants_mod, only: HLv,HLs,Cp_air,Grav,rdgas,rvgas, & kappa, es0 implicit none private !--------------------------------------------------------------------- ! ---- public interfaces ---- public betts_miller, betts_miller_init, betts_miller_end !----------------------------------------------------------------------- ! ---- version number ---- character(len=128) :: version = '$Id: betts_miller.F90,v 17.0.2.1 2009/07/24 21:48:03 wfc Exp $' character(len=128) :: tagname = '$Name: mom4p1_pubrel_dec2009_nnz $' !----------------------------------------------------------------------- ! ---- local/private data ---- real, parameter :: d622 = rdgas/rvgas real, parameter :: d378 = 1.-d622 logical :: module_is_initialized=.false. !----------------------------------------------------------------------- ! --- namelist ---- real :: tau_bm=7200. real :: rhbm = .8 logical :: do_simp = .true. !logical :: do_enadjusttemp = .false. logical :: do_shallower = .false. logical :: do_changeqref = .false. logical :: do_envsat = .false. logical :: do_taucape = .false. real :: capetaubm = 900. real :: tau_min = 2400. namelist /betts_miller_nml/ tau_bm, rhbm, do_simp, & ! do_enadjusttemp, & do_shallower, do_changeqref, & do_envsat, do_taucape, capetaubm, tau_min !----------------------------------------------------------------------- ! description of namelist variables ! ! tau_bm = betts-miller relaxation timescale (seconds) ! ! rhbm = relative humidity that you're relaxing towards ! ! do_simp = do the simple method where you adjust timescales to make ! precip continuous always ! ! do_enadjusttemp = do the betts-miller (1986) way of conserving energy: ! adjust the reference temperature by a fixed amount w/ ! height. ! *** if neither of the two above are true, then we instead adjust the ! humidity profile so its precipitation is increased to balance the ! precip_t change, and hence dries more than the ! ! do_shallower = do the shallow convection scheme where it chooses a smaller ! depth such that precipitation is zero ! ! do_changeqref = do the shallow convection scheme where if changes the ! profile of both q and T in order make precip zero ! ! do_envsat = reference profile is rhbm times saturated wrt environment ! (if false, it's rhbm times parcel) ! ! do_taucape = scheme where taubm is proportional to CAPE**-1/2 ! ! capetaubm = for the above scheme, the value of CAPE for which ! tau = tau_bm ! ! tau_min = minimum relaxation time allowed for the above scheme ! ! ! !----------------------------------------------------------------------- contains !####################################################################### subroutine betts_miller (dt, tin, qin, pfull, phalf, coldT, & rain, snow, tdel, qdel, q_ref, bmflag, & klzbs, cape, cin, t_ref,invtau_bm_t,invtau_bm_q, & mask, conv) !----------------------------------------------------------------------- ! ! Betts-Miller Convection Scheme ! !----------------------------------------------------------------------- ! ! input: dt time step in seconds ! tin temperature at full model levels ! qin specific humidity of water vapor at full ! model levels ! pfull pressure at full model levels ! phalf pressure at half (interface) model levels ! coldT should precipitation be snow at this point? ! optional: ! mask optional mask (0 or 1.) ! conv logical flag; if true then no betts-miller ! adjustment is performed at that grid-point or ! model level ! ! output: rain liquid precipitation (kg/m2) ! snow frozen precipitation (kg/m2) ! tdel temperature tendency at full model levels ! qdel specific humidity tendency (of water vapor) at ! full model levels ! bmflag flag for which routines you're calling ! klzbs stored klzb values ! cape convectively available potential energy ! cin convective inhibition (this and the above are before the ! adjustment) ! invtau_bm_t temperature relaxation timescale ! invtau_bm_q humidity relaxation timescale ! !----------------------------------------------------------------------- !--------------------- interface arguments ----------------------------- real , intent(in) , dimension(:,:,:) :: tin, qin, pfull, phalf real , intent(in) :: dt logical , intent(in) , dimension(:,:):: coldT real , intent(out), dimension(:,:) :: rain,snow, bmflag, klzbs, cape, & cin, invtau_bm_t, invtau_bm_q real , intent(out), dimension(:,:,:) :: tdel, qdel, q_ref, t_ref real , intent(in) , dimension(:,:,:), optional :: mask logical, intent(in) , dimension(:,:,:), optional :: conv !----------------------------------------------------------------------- !---------------------- local data ------------------------------------- logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust logical :: avgbl real,dimension(size(tin,1),size(tin,2),size(tin,3)) :: & rin, esat, qsat, desat, dqsat, pmes, pmass real,dimension(size(tin,1),size(tin,2)) :: & hlcp, precip, precip_t real,dimension(size(tin,3)) :: eref, rpc, tpc, & tpc1, rpc1 real :: & cape1, cin1, tot, deltak, deltaq, qrefint, deltaqfrac, deltaqfrac2, & ptopfrac, es, plzb, plcl, cape2, small integer i, j, k, ix, jx, kx, klzb, ktop, klzb2 !----------------------------------------------------------------------- ! computation of precipitation by betts-miller scheme !----------------------------------------------------------------------- if (.not. module_is_initialized) call error_mesg ('betts_miller', & 'betts_miller_init has not been called.', FATAL) ix=size(tin,1) jx=size(tin,2) kx=size(tin,3) avgbl = .false. small = 1.e-10 ! calculate r rin = qin/(1.0 - qin) do i=1,ix do j=1,jx cape1 = 0. cin1 = 0. tot = 0. klzb=0 ! the bmflag is written out to show what aspects of the bm scheme is called ! bmflag = 0 is no cape, no convection ! bmflag = 1 is shallow conv, the predicted precip is less than zero ! bmflag = 2 is deep convection bmflag(i,j) = 0. tpc = tin(i,j,:) rpc = rin(i,j,:) ! calculate cape, cin, level of zero buoyancy, and parcel properties ! new code (second order in delta ln p and exact LCL calculation) call capecalcnew( kx, pfull(i,j,:), phalf(i,j,:),& cp_air, rdgas, rvgas, hlv, kappa, tin(i,j,:), & rin(i,j,:), avgbl, cape1, cin1, tpc, & rpc, klzb) ! set values for storage cape(i,j) = cape1 cin(i,j) = cin1 klzbs(i,j) = klzb if(cape1.gt.0.) then ! if((tot.gt.0.).and.(cape1.gt.0.)) then bmflag(i,j) = 1. ! reference temperature is just that of the parcel all the way up t_ref(i,j,:) = tpc do k=klzb,kx ! sets reference spec hum to a certain relative hum (change to vapor pressure, ! multiply by rhbm, then back to spec humid) if(do_envsat) then ! call establ2(es,tin(i,j,k)) call escomp(tin(i,j,k),es) es = es*rhbm ! rpc(k) = d622*es/(pfull(i,j,k) - d378*es) rpc(k) = rdgas/rvgas*es/pfull(i,j,k) q_ref(i,j,k) = rpc(k)/(1 + rpc(k)) else rpc(k) = rhbm*rpc(k) ! eref(k) = rhbm*pfull(i,j,k)*rpc(k)/(d622 + d378*rpc(k)) ! rpc(k) = d622*eref(k)/(pfull(i,j,k) - d378*eref(k)) q_ref(i,j,k) = rpc(k)/(1 + rpc(k)) endif end do ! set the tendencies to zero where you don't adjust ! set the reference profiles to be the original profiles (for diagnostic ! purposes only -- you can think of this as what you're relaxing to in ! areas above the actual convection do k=1,max(klzb-1,1) qdel(i,j,k) = 0.0 tdel(i,j,k) = 0.0 q_ref(i,j,k) = qin(i,j,k) t_ref(i,j,k) = tin(i,j,k) end do ! initialize p to zero for the loop precip(i,j) = 0. precip_t(i,j) = 0. ! makes t_bm prop to (CAPE)**-.5. Gives a relaxation time of tau_bm when ! CAPE = sqrt(capetaubm) if(do_taucape) then tau_bm = sqrt(capetaubm)*tau_bm/sqrt(cape1) if(tau_bm.lt.tau_min) tau_bm = tau_min endif do k=klzb, kx ! relax to reference profiles tdel(i,j,k) = - (tin(i,j,k) - t_ref(i,j,k))/tau_bm*dt qdel(i,j,k) = - (qin(i,j,k) - q_ref(i,j,k))/tau_bm*dt ! Precipitation can be calculated already, based on the change in q on the ! way up (this isn't altered in the energy conservation scheme). precip(i,j) = precip(i,j) - qdel(i,j,k)*(phalf(i,j,k+1)- & phalf(i,j,k))/grav precip_t(i,j)= precip_t(i,j) + cp_air/(hlv+small)*tdel(i,j,k)* & (phalf(i,j,k+1)-phalf(i,j,k))/grav end do if ((precip(i,j).gt.0.).and.(precip_t(i,j).gt.0.)) then ! If precip > 0, then correct energy. bmflag(i,j) = 2. if(precip(i,j).gt.precip_t(i,j)) then ! if the q precip is greater, then lengthen the relaxation timescale on q to ! conserve energy. qdel is therefore changed. invtau_bm_q(i,j) = precip_t(i,j)/precip(i,j)/tau_bm qdel(i,j,klzb:kx) = tau_bm*invtau_bm_q(i,j)* & qdel(i,j,klzb:kx) precip(i,j) = precip_t(i,j) invtau_bm_t(i,j) = 1./tau_bm else if(do_simp) then ! simple scheme: ! if the t precip is greater, then lengthen the relaxation timescale on t to ! conserve energy. tdel is therefore changed. invtau_bm_t(i,j) = precip(i,j)/precip_t(i,j)/tau_bm tdel(i,j,klzb:kx) = tau_bm*invtau_bm_t(i,j)* & tdel(i,j,klzb:kx) invtau_bm_q(i,j) = 1./tau_bm ! else if(do_enadjusttemp) else ! not simple scheme: shift the reference profile of temperature ! deltak is the energy correction that you make to the temperature reference ! profile deltak = 0. do k=klzb, kx ! Calculate the integrated difference in energy change within each level. deltak = deltak - (tdel(i,j,k) + hlv/cp_air*& qdel(i,j,k))* & (phalf(i,j,k+1) - phalf(i,j,k)) end do ! Divide by total pressure. deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,klzb)) ! Adjust the reference profile (uniformly with height), and correspondingly ! the temperature change. t_ref(i,j,klzb:kx) = t_ref(i,j,klzb:kx)+ & deltak*tau_bm/dt tdel(i,j,klzb:kx) = tdel(i,j,klzb:kx) + deltak endif endif else if(precip_t(i,j).gt.0.) then ! If precip < 0, then do the shallow conv routine. ! First option: do_shallower = true ! This chooses the depth of convection based on choosing the height that ! it can make precip zero, i.e., subtract off heights until that precip ! becomes positive. if (do_shallower) then ! ktop is the new top of convection. set this initially to klzb. ktop = klzb ! Work your way down until precip is positive again. do while ((precip(i,j).lt.0.).and.(ktop.le.kx)) precip(i,j) = precip(i,j) - qdel(i,j,ktop)* & (phalf(i,j,ktop) - phalf(i,j,ktop+1))/grav ktop = ktop + 1 end do ! since there will be an overshoot (precip is going to be greater than zero ! once we finish this), the actual new top of convection is somewhere between ! the current ktop, and one level above this. set ktop to the level above. ktop = ktop - 1 ! Adjust the tendencies in the places above back to zero, and the reference ! profiles back to the original t,q. if (ktop.gt.klzb) then qdel(i,j,klzb:ktop-1) = 0. q_ref(i,j,klzb:ktop-1) = qin(i,j,klzb:ktop-1) tdel(i,j,klzb:ktop-1) = 0. t_ref(i,j,klzb:ktop-1) = tin(i,j,klzb:ktop-1) end if ! Then make the change only a fraction of the new top layer so the precip is ! identically zero. ! Calculate the fractional penetration of convection through that top layer. ! This is the amount necessary to make precip identically zero. if (precip(i,j).gt.0.) then ptopfrac = precip(i,j)/(qdel(i,j,ktop)* & (phalf(i,j,ktop+1) - phalf(i,j,ktop)))*grav ! Reduce qdel in the top layer by this fraction. qdel(i,j,ktop) = ptopfrac*qdel(i,j,ktop) ! Set precip to zero precip(i,j) = 0. ! Now change the reference temperature in such a way to make the net ! heating zero. !! Reduce tdel in the top layer tdel(i,j,ktop) = ptopfrac*tdel(i,j,ktop) deltak = 0. if (ktop.lt.kx) then ! Integrate temperature tendency up to 1 level below top. do k=ktop,kx deltak = deltak + tdel(i,j,k)* & (phalf(i,j,k) - phalf(i,j,k+1)) end do ! Normalize by the pressure difference. deltak = deltak/(phalf(i,j,kx+1) - phalf(i,j,ktop)) ! Subtract this value uniformly from tdel, and make the according change to ! t_ref. do k=ktop,kx tdel(i,j,k) = tdel(i,j,k) + deltak t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt end do end if else precip(i,j) = 0. qdel(i,j,kx) = 0. q_ref(i,j,kx) = qin(i,j,kx) tdel(i,j,kx) = 0. t_ref(i,j,kx) = tin(i,j,kx) invtau_bm_t(i,j) = 0. invtau_bm_q(i,j) = 0. end if else if(do_changeqref) then ! Change the reference profile of q by a certain fraction so that precip is ! zero. This involves calculating the total integrated q_ref dp (this is the ! quantity intqref), as well as the necessary change in q_ref (this is the ! quantity deltaq). Then the fractional change in q_ref at each level (the ! quantity deltaqfrac) is 1-deltaq/intqref. (have to multiply q_ref by ! 1-deltaq/intqref at every level) Then the change in qdel is ! -deltaq/intqref*q_ref*dt/tau_bm. ! Change the reference profile of T by a uniform amount so that precip is zero. deltak = 0. deltaq = 0. qrefint = 0. do k=klzb,kx ! deltaq = a positive quantity (since int qdel is positive). It's how ! much q_ref must be changed by, in an integrated sense. The requisite ! change in qdel is this without the factors of tau_bm and dt. deltaq = deltaq - qdel(i,j,k)*tau_bm/dt* & (phalf(i,j,k) - phalf(i,j,k+1)) ! deltak = the amount tdel needs to be changed deltak = deltak + tdel(i,j,k)* & (phalf(i,j,k) - phalf(i,j,k+1)) ! qrefint = integrated value of qref qrefint = qrefint - q_ref(i,j,k)* & (phalf(i,j,k) - phalf(i,j,k+1)) end do ! Normalize deltak by total pressure. deltak = deltak /(phalf(i,j,kx+1) - phalf(i,j,klzb)) ! multiplying factor for q_ref is 1 + the ratio deltaqfrac = 1. - deltaq/qrefint ! multiplying factor for qdel adds dt/tau_bm deltaqfrac2 = - deltaq/qrefint*dt/tau_bm precip(i,j) = 0.0 do k=klzb,kx qdel(i,j,k) = qdel(i,j,k) + deltaqfrac2*q_ref(i,j,k) q_ref(i,j,k) = deltaqfrac*q_ref(i,j,k) tdel(i,j,k) = tdel(i,j,k) + deltak t_ref(i,j,k) = t_ref(i,j,k) + deltak*tau_bm/dt end do else precip(i,j) = 0. tdel(i,j,:) = 0. qdel(i,j,:) = 0. invtau_bm_t(i,j) = 0. invtau_bm_q(i,j) = 0. end if ! for cases where virtual temp predicts CAPE but precip_t < 0. else tdel(i,j,:) = 0.0 qdel(i,j,:) = 0.0 precip(i,j) = 0.0 q_ref(i,j,:) = qin(i,j,:) t_ref(i,j,:) = tin(i,j,:) invtau_bm_t(i,j) = 0. invtau_bm_q(i,j) = 0. end if ! if no CAPE, set tendencies to zero. else tdel(i,j,:) = 0.0 qdel(i,j,:) = 0.0 precip(i,j) = 0.0 q_ref(i,j,:) = qin(i,j,:) t_ref(i,j,:) = tin(i,j,:) invtau_bm_t(i,j) = 0. invtau_bm_q(i,j) = 0. end if end do end do rain = precip snow = 0. end subroutine betts_miller !####################################################################### !all new cape calculation. subroutine capecalcnew(kx,p,phalf,cp_air,rdgas,rvgas,hlv,kappa,tin,rin,& avgbl,cape,cin,tp,rp,klzb) ! ! Input: ! ! kx number of levels ! p pressure (index 1 refers to TOA, index kx refers to surface) ! phalf pressure at half levels ! cp_air specific heat of dry air ! rdgas gas constant for dry air ! rvgas gas constant for water vapor (used in Clausius-Clapeyron, ! not for virtual temperature effects, which are not considered) ! hlv latent heat of vaporization ! kappa the constant kappa ! tin temperature of the environment ! rin specific humidity of the environment ! avgbl if true, the parcel is averaged in theta and r up to its LCL ! ! Output: ! cape Convective available potential energy ! cin Convective inhibition (if there's no LFC, then this is set ! to zero) ! tp Parcel temperature (set to the environmental temperature ! where no adjustment) ! rp Parcel specific humidity (set to the environmental humidity ! where no adjustment, and set to the saturation humidity at ! the parcel temperature below the LCL) ! klzb Level of zero buoyancy ! ! Algorithm: ! Start with surface parcel. ! Calculate the lifting condensation level (uses an analytic formula and a ! lookup table). ! Average under the LCL if desired, if this is done, then a new LCL must ! be calculated. ! Calculate parcel ascent up to LZB. ! Calculate CAPE and CIN. implicit none integer, intent(in) :: kx logical, intent(in) :: avgbl real, intent(in), dimension(:) :: p, phalf, tin, rin real, intent(in) :: rdgas, rvgas, hlv, kappa, cp_air integer, intent(out) :: klzb real, intent(out), dimension(:) :: tp, rp real, intent(out) :: cape, cin integer :: k, klcl, klfc, klcl2 logical :: nocape real, dimension(kx) :: theta real :: t0, r0, es, rs, theta0, pstar, value, tlcl, & a, b, dtdlnp, d2tdlnp2, thetam, rm, tlcl2, & plcl2, plcl, plzb, small pstar = 1.e5 ! so we can run dry limit (one expression involves 1/hlv) small = 1.e-10 nocape = .true. cape = 0. cin = 0. plcl = 0. plzb = 0. klfc = 0 klcl = 0 klzb = 0 tp(1:kx) = tin(1:kx) rp(1:kx) = rin(1:kx) ! start with surface parcel t0 = tin(kx) r0 = rin(kx) ! calculate the lifting condensation level by the following: ! are you saturated to begin with? call escomp(t0,es) rs = rdgas/rvgas*es/p(kx) if (r0.ge.rs) then ! if youčre already saturated, set lcl to be the surface value. plcl = p(kx) ! the first level where youčre completely saturated. klcl = kx ! saturate out to get the parcel temp and humidity at this level ! first order (in delta T) accurate expression for change in temp tp(kx) = t0 + (r0 - rs)/(cp_air/(hlv+small) + hlv*rs/rvgas/t0**2.) call escomp(tp(kx),es) rp(kx) = rdgas/rvgas*es/p(kx) else ! if not saturated to begin with, use the analytic expression to calculate the ! exact pressure and temperature where youčre saturated. theta0 = tin(kx)*(pstar/p(kx))**kappa ! the expression that we utilize is ! log(r/theta**(1/kappa)*pstar*rvgas/rdgas/es00) = log(es/T**(1/kappa)) ! (the division by es00 is necessary because the RHS values are tabulated ! for control moisture content) ! The right hand side of this is only a function of temperature, therefore ! this is put into a lookup table to solve for temperature. if (r0.gt.0.) then value = log(theta0**(-1/kappa)*r0*pstar*rvgas/rdgas/es0) call lcltabl(value,tlcl) plcl = pstar*(tlcl/theta0)**(1/kappa) ! just in case plcl is very high up if (plcl.lt.p(1)) then plcl = p(1) tlcl = theta0*(plcl/pstar)**kappa write (*,*) 'hi lcl' end if k = kx else ! if the parcel sp hum is zero or negative, set lcl to top level plcl = p(1) tlcl = theta0*(plcl/pstar)**kappa ! write (*,*) 'zero r0', r0 do k=1,kx tp(k) = theta0*(p(k)/pstar)**kappa rp(k) = 0. ! this definition of CIN contains everything below the LCL cin = cin + rdgas*(tin(k) - tp(k))*log(phalf(k+1)/phalf(k)) end do go to 11 end if ! calculate the parcel temperature (adiabatic ascent) below the LCL. ! the mixing ratio stays the same do while (p(k).gt.plcl) tp(k) = theta0*(p(k)/pstar)**kappa call escomp(tp(k),es) rp(k) = rdgas/rvgas*es/p(k) ! this definition of CIN contains everything below the LCL cin = cin + rdgas*(tin(k) - tp(k))*log(phalf(k+1)/phalf(k)) k = k-1 end do ! first level where you're saturated at the level klcl = k if (klcl.eq.1) klcl = 2 ! do a saturated ascent to get the parcel temp at the LCL. ! use your 2nd order equation up to the pressure above. ! moist adaibat derivatives: (use the lcl values for temp, humid, and ! pressure) a = kappa*tlcl + hlv/cp_air*r0 b = hlv**2.*r0/cp_air/rvgas/tlcl**2. dtdlnp = a/(1. + b) ! first order in p ! tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl) ! second order in p (RK2) ! first get temp halfway up tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl)/2. if ((tp(klcl).lt.173.16).and.nocape) go to 11 call escomp(tp(klcl),es) rp(klcl) = rdgas/rvgas*es/(p(klcl) + plcl)*2. a = kappa*tp(klcl) + hlv/cp_air*rp(klcl) b = hlv**2./cp_air/rvgas*rp(klcl)/tp(klcl)**2. dtdlnp = a/(1. + b) ! second half of RK2 tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl) ! d2tdlnp2 = (kappa + b - 1. - b/tlcl*(hlv/rvgas/tlcl - & ! 2.)*dtdlnp)/ (1. + b)*dtdlnp - hlv*r0/cp_air/ & ! (1. + b) ! second order in p ! tp(klcl) = tlcl + dtdlnp*log(p(klcl)/plcl) + .5*d2tdlnp2*(log(& ! p(klcl)/plcl))**2. if ((tp(klcl).lt.173.16).and.nocape) go to 11 call escomp(tp(klcl),es) rp(klcl) = rdgas/rvgas*es/p(klcl) ! write (*,*) 'tp, rp klcl:kx, new', tp(klcl:kx), rp(klcl:kx) ! CAPE/CIN stuff if ((tp(klcl).lt.tin(klcl)).and.nocape) then ! if youčre not yet buoyant, then add to the CIN and continue cin = cin + rdgas*(tin(klcl) - & tp(klcl))*log(phalf(klcl+1)/phalf(klcl)) else ! if youčre buoyant, then add to cape cape = cape + rdgas*(tp(klcl) - & tin(klcl))*log(phalf(klcl+1)/phalf(klcl)) ! if itčs the first time buoyant, then set the level of free convection to k if (nocape) then nocape = .false. klfc = klcl endif end if end if ! then average the properties over the boundary layer if so desired. to give ! a new "parcel". this may not be saturated at the LCL, so make sure you get ! to a level where it is before moist adiabatic ascent! !!!! take out all the below (between the exclamation points) if no avgbl !!!! if (avgbl) then theta(klcl:kx) = tin(klcl:kx)*(pstar/p(klcl:kx))**kappa thetam = 0. rm = 0. do k=klcl,kx thetam = thetam + theta(k)*(phalf(k+1) - phalf(k)) rm = rm + rin(k)*(phalf(k+1) - phalf(k)) end do thetam = thetam/(phalf(kx+1) - phalf(klcl)) rm = rm/(phalf(kx+1) - phalf(klcl)) ! check if you're saturated at the top level. if not, then get a new LCL tp(klcl) = thetam*(p(klcl)/pstar)**kappa call escomp(tp(klcl),es) rs = rdgas/rvgas*es/p(klcl) ! if you're not saturated, get a new LCL if (rm.lt.rs) then ! reset CIN to zero. cin = 0. ! again, use the analytic expression to calculate the exact pressure and ! temperature where youčre saturated. ! the expression that we utilize is ! log(r/theta**(1/kappa)*pstar*rvgas/rdgas/es00)= log(es/T**(1/kappa)) ! (the division by es00 is necessary because the RHS values are tabulated ! for control moisture content) ! The right hand side of this is only a function of temperature, therefore ! this is put into a lookup table to solve for temperature. value = log(thetam**(-1/kappa)*rm*pstar*rvgas/rdgas/es0) call lcltabl(value,tlcl2) plcl2 = pstar*(tlcl2/thetam)**(1/kappa) ! just in case plcl is very high up if (plcl2.lt.p(1)) then plcl2 = p(1) end if k = kx ! calculate the parcel temperature (adiabatic ascent) below the LCL. ! the mixing ratio stays the same do while (p(k).gt.plcl2) tp(k) = thetam*(p(k)/pstar)**kappa call escomp(tp(k),es) rp(k) = rdgas/rvgas*es/p(k) ! this definition of CIN contains everything below the LCL cin = cin + rdgas*(tin(k) - tp(k))*log(phalf(k+1)/phalf(k)) k = k-1 end do ! first level where youčre saturated at the level klcl2 = k if (klcl2.eq.1) klcl2 = 2 ! do a saturated ascent to get the parcel temp at the LCL. ! use your 2nd order equation up to the pressure above. ! moist adaibat derivatives: (use the lcl values for temp, humid, and ! pressure) a = kappa*tlcl2 + hlv/cp_air*rm b = hlv**2.*rm/cp_air/rvgas/tlcl2**2. dtdlnp = a/(1. + b) ! first order in p ! tp(klcl2) = tlcl2 + dtdlnp*log(p(klcl2)/plcl2) ! second order in p (RK2) ! first get temp halfway up tp(klcl2) = tlcl2 + dtdlnp*log(p(klcl2)/plcl2)/2. if ((tp(klcl2).lt.173.16).and.nocape) go to 11 call escomp(tp(klcl2),es) rp(klcl2) = rdgas/rvgas*es/(p(klcl2) + plcl2)*2. a = kappa*tp(klcl2) + hlv/cp_air*rp(klcl2) b = hlv**2./cp_air/rvgas*rp(klcl2)/tp(klcl2)**2. dtdlnp = a/(1. + b) ! second half of RK2 tp(klcl2) = tlcl2 + dtdlnp*log(p(klcl2)/plcl2) ! d2tdlnp2 = (kappa + b - 1. - b/tlcl2*(hlv/rvgas/tlcl2 - & ! 2.)*dtdlnp)/ (1. + b)*dtdlnp - hlv*rm/cp_air/ & ! (1. + b) ! second order in p ! tp(klcl2) = tlcl2 + dtdlnp*log(p(klcl2)/plcl2) + & ! .5*d2tdlnp2*(log(p(klcl2)/plcl2))**2. call escomp(tp(klcl2),es) rp(klcl2) = rdgas/rvgas*es/p(klcl2) ! CAPE/CIN stuff if ((tp(klcl2).lt.tin(klcl2)).and.nocape) then ! if youčre not yet buoyant, then add to the CIN and continue cin = cin + rdgas*(tin(klcl2) - & tp(klcl2))*log(phalf(klcl2+1)/phalf(klcl2)) else ! if youčre buoyant, then add to cape cape = cape + rdgas*(tp(klcl) - & tin(klcl))*log(phalf(klcl+1)/phalf(klcl)) ! if itčs the first time buoyant, then set the level of free convection to k if (nocape) then nocape = .false. klfc = klcl2 endif end if end if end if !!!! take out all of the above (within the exclamations) if no avgbl !!!! ! then, start at the LCL, and do moist adiabatic ascent by the first order ! scheme -- 2nd order as well do k=klcl-1,1,-1 a = kappa*tp(k+1) + hlv/cp_air*rp(k+1) b = hlv**2./cp_air/rvgas*rp(k+1)/tp(k+1)**2. dtdlnp = a/(1. + b) ! first order in p ! tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1)) ! second order in p (RK2) ! first get temp halfway up tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1))/2. if ((tp(k).lt.173.16).and.nocape) go to 11 call escomp(tp(k),es) rp(k) = rdgas/rvgas*es/(p(k) + p(k+1))*2. a = kappa*tp(k) + hlv/cp_air*rp(k) b = hlv**2./cp_air/rvgas*rp(k)/tp(k)**2. dtdlnp = a/(1. + b) ! second half of RK2 tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1)) ! d2tdlnp2 = (kappa + b - 1. - b/tp(k+1)*(hlv/rvgas/tp(k+1) - & ! 2.)*dtdlnp)/(1. + b)*dtdlnp - hlv/cp_air*rp(k+1)/(1. + b) ! second order in p ! tp(k) = tp(k+1) + dtdlnp*log(p(k)/p(k+1)) + .5*d2tdlnp2*(log( & ! p(k)/p(k+1)))**2. ! if you're below the lookup table value, just presume that there's no way ! you could have cape and call it quits if ((tp(k).lt.173.16).and.nocape) go to 11 call escomp(tp(k),es) rp(k) = rdgas/rvgas*es/p(k) if ((tp(k).lt.tin(k)).and.nocape) then ! if youčre not yet buoyant, then add to the CIN and continue cin = cin + rdgas*(tin(k) - tp(k))*log(phalf(k+1)/phalf(k)) elseif((tp(k).lt.tin(k)).and.(.not.nocape)) then ! if you have CAPE, and itčs your first time being negatively buoyant, ! then set the level of zero buoyancy to k+1, and stop the moist ascent klzb = k+1 go to 11 else ! if youčre buoyant, then add to cape cape = cape + rdgas*(tp(k) - tin(k))*log(phalf(k+1)/phalf(k)) ! if itčs the first time buoyant, then set the level of free convection to k if (nocape) then nocape = .false. klfc = k endif end if end do 11 if(nocape) then ! this is if you made it through without having a LZB ! set LZB to be the top level. plzb = p(1) klzb = 0 klfc = 0 cin = 0. tp(1:kx) = tin(1:kx) rp(1:kx) = rin(1:kx) end if ! write (*,*) 'plcl, klcl, tlcl, r0 new', plcl, klcl, tlcl, r0 ! write (*,*) 'tp, rp new', tp, rp ! write (*,*) 'tp, new', tp ! write (*,*) 'tin new', tin ! write (*,*) 'klcl, klfc, klzb new', klcl, klfc, klzb end subroutine capecalcnew ! lookup table for the analytic evaluation of LCL subroutine lcltabl(value,tlcl) ! ! Table of values used to compute the temperature of the lifting condensation ! level. ! ! the expression that we utilize is ! log(r/theta**(1/kappa)*pstar*rvgas/rdgas/es00) = log(es/T**(1/kappa)) ! the RHS is tabulated for control moisture content, hence the division ! by es00 on the LHS ! ! Gives the values of the temperature for the following range: ! starts with -23, is uniformly distributed up to -10.4. There are a ! total of 127 values, and the increment is .1. ! implicit none real, intent(in) :: value real, intent(out) :: tlcl integer :: ival real, dimension(127) :: lcltable real :: v1, v2 data lcltable/ 1.7364512e+02, 1.7427449e+02, 1.7490874e+02, & 1.7554791e+02, 1.7619208e+02, 1.7684130e+02, 1.7749563e+02, & 1.7815514e+02, 1.7881989e+02, 1.7948995e+02, 1.8016539e+02, & 1.8084626e+02, 1.8153265e+02, 1.8222461e+02, 1.8292223e+02, & 1.8362557e+02, 1.8433471e+02, 1.8504972e+02, 1.8577068e+02, & 1.8649767e+02, 1.8723077e+02, 1.8797006e+02, 1.8871561e+02, & 1.8946752e+02, 1.9022587e+02, 1.9099074e+02, 1.9176222e+02, & 1.9254042e+02, 1.9332540e+02, 1.9411728e+02, 1.9491614e+02, & 1.9572209e+02, 1.9653521e+02, 1.9735562e+02, 1.9818341e+02, & 1.9901870e+02, 1.9986158e+02, 2.0071216e+02, 2.0157057e+02, & 2.0243690e+02, 2.0331128e+02, 2.0419383e+02, 2.0508466e+02, & 2.0598391e+02, 2.0689168e+02, 2.0780812e+02, 2.0873335e+02, & 2.0966751e+02, 2.1061074e+02, 2.1156316e+02, 2.1252493e+02, & 2.1349619e+02, 2.1447709e+02, 2.1546778e+02, 2.1646842e+02, & 2.1747916e+02, 2.1850016e+02, 2.1953160e+02, 2.2057364e+02, & 2.2162645e+02, 2.2269022e+02, 2.2376511e+02, 2.2485133e+02, & 2.2594905e+02, 2.2705847e+02, 2.2817979e+02, 2.2931322e+02, & 2.3045895e+02, 2.3161721e+02, 2.3278821e+02, 2.3397218e+02, & 2.3516935e+02, 2.3637994e+02, 2.3760420e+02, 2.3884238e+02, & 2.4009473e+02, 2.4136150e+02, 2.4264297e+02, 2.4393941e+02, & 2.4525110e+02, 2.4657831e+02, 2.4792136e+02, 2.4928053e+02, & 2.5065615e+02, 2.5204853e+02, 2.5345799e+02, 2.5488487e+02, & 2.5632953e+02, 2.5779231e+02, 2.5927358e+02, 2.6077372e+02, & 2.6229310e+02, 2.6383214e+02, 2.6539124e+02, 2.6697081e+02, & 2.6857130e+02, 2.7019315e+02, 2.7183682e+02, 2.7350278e+02, & 2.7519152e+02, 2.7690354e+02, 2.7863937e+02, 2.8039954e+02, & 2.8218459e+02, 2.8399511e+02, 2.8583167e+02, 2.8769489e+02, & 2.8958539e+02, 2.9150383e+02, 2.9345086e+02, 2.9542719e+02, & 2.9743353e+02, 2.9947061e+02, 3.0153922e+02, 3.0364014e+02, & 3.0577420e+02, 3.0794224e+02, 3.1014515e+02, 3.1238386e+02, & 3.1465930e+02, 3.1697246e+02, 3.1932437e+02, 3.2171609e+02, & 3.2414873e+02, 3.2662343e+02, 3.2914139e+02, 3.3170385e+02 / v1 = value if (value.lt.-23.0) v1 = -23.0 if (value.gt.-10.4) v1 = -10.4 ival = floor(10.*(v1 + 23.0)) v2 = -230. + ival v1 = 10.*v1 tlcl = (v2 + 1.0 - v1)*lcltable(ival+1) + (v1 - v2)*lcltable(ival+2) end subroutine lcltabl !####################################################################### subroutine betts_miller_init () !----------------------------------------------------------------------- ! ! initialization for betts_miller ! !----------------------------------------------------------------------- integer unit,io,ierr, logunit !----------- read namelist --------------------------------------------- if (file_exist('input.nml')) then unit = open_namelist_file ( ) ierr=1; do while (ierr /= 0) read (unit, nml=betts_miller_nml, iostat=io, end=10) ierr = check_nml_error (io,'betts_miller_nml') enddo 10 call close_file (unit) endif !---------- output namelist -------------------------------------------- call write_version_number ( version, tagname ) logunit = stdlog() if ( mpp_pe() == mpp_root_pe() ) then write (logunit,nml=betts_miller_nml) endif module_is_initialized =.true. end subroutine betts_miller_init !####################################################################### subroutine betts_miller_end () module_is_initialized =.false. end subroutine betts_miller_end end module betts_miller_mod