#undef DEBUG #define WRF_PORT #define MODAL_AERO ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov module cldwat !----------------------------------------------------------------------- ! ! Purpose: Prognostic cloud water data and methods. ! ! Public interfaces: ! ! inimc -- Initialize constants ! pcond -- Calculate prognostic condensate ! ! cldwat_fice -- calculate fraction of condensate in ice phase (radiation partitioning) ! ! Author: P. Rasch, with Modifications by Minghua Zhang ! January 2010, modified by J. Kay to add precip fluxes for COSP simulator ! Ported to WRF by William.Gustafson@pnl.gov, Nov. 2009 ! updated to CESM1_0_1, Nov. 2010 ! !----------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 #ifndef WRF_PORT use spmd_utils, only: masterproc use ppgrid, only: pcols, pver, pverp #endif use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, & cp, epsqs, ttrice #ifndef WRF_PORT use abortutils, only: endrun use cam_logfile, only: iulog #else use module_cam_support, only: masterproc, & pcols, pver, pverp, & endrun, & iulog #endif implicit none !----------------------------------------------------------------------- ! PUBLIC: Make default data and interfaces private !----------------------------------------------------------------------- private save #ifndef WRF_PORT public inimc, pcond ! Public interfaces #endif public cldwat_fice ! Public interfaces public cldwat_readnl integer, public:: ktop ! Level above 10 hPa #ifndef WRF_PORT real(r8),public :: icritc ! threshold for autoconversion of cold ice real(r8),public :: icritw ! threshold for autoconversion of warm ice !!$ real(r8),public,parameter:: conke = 1.e-6 ! tunable constant for evaporation of precip !!$ real(r8),public,parameter:: conke = 2.e-6 ! tunable constant for evaporation of precip real(r8),public :: conke ! tunable constant for evaporation of precip real(r8),public :: r3lcrit ! critical radius where liq conversion begins #else ! Currently, the WRF_PORT bipasses the namelist initialization of these ! tunable parameters. We are hard-coding them to the default values for ! the fv 0.23x0.31 grid. One might choose to implement an option to set ! these via the WRF Registry in the future. real(r8),public :: icritw = 2.0e-4 real(r8),public :: icritc = 45.0e-6 real(r8),public :: conke = 5.0e-6 #endif !----------------------------------------------------------------------- ! PRIVATE: Everything else is private to this module !----------------------------------------------------------------------- real(r8), private:: tmax_fice! max temperature for cloud ice formation real(r8), private:: tmin_fice! min temperature for cloud ice formation real(r8), private:: tmax_fsnow ! max temperature for transition to convective snow real(r8), private:: tmin_fsnow ! min temperature for transition to convective snow real(r8), private:: rhonot ! air density at surface real(r8), private:: t0 ! Freezing temperature real(r8), private:: cldmin ! assumed minimum cloud amount real(r8), private:: small ! small number compared to unity real(r8), private:: c ! constant for graupel like snow cm**(1-d)/s real(r8), private:: d ! constant for graupel like snow real(r8), private:: esi ! collection efficient for ice by snow real(r8), private:: esw ! collection efficient for water by snow real(r8), private:: nos ! particles snow / cm**4 real(r8), private:: pi ! Mathematical constant real(r8), private:: gravit ! Gravitational acceleration at surface real(r8), private:: rh2o real(r8), private:: prhonos real(r8), private:: thrpd ! numerical three added to d real(r8), private:: gam3pd ! gamma function on (3+d) real(r8), private:: gam4pd ! gamma function on (4+d) real(r8), private:: rhoi ! ice density real(r8), private:: rhos ! snow density real(r8), private:: rhow ! water density real(r8), private:: mcon01 ! constants used in cloud microphysics real(r8), private:: mcon02 ! constants used in cloud microphysics real(r8), private:: mcon03 ! constants used in cloud microphysics real(r8), private:: mcon04 ! constants used in cloud microphysics real(r8), private:: mcon05 ! constants used in cloud microphysics real(r8), private:: mcon06 ! constants used in cloud microphysics real(r8), private:: mcon07 ! constants used in cloud microphysics real(r8), private:: mcon08 ! constants used in cloud microphysics integer, private :: k1mb ! index of the eta level near 1 mb ! Parameters used in findmcnew real(r8) :: capnsi ! sea ice cloud particles / cm3 real(r8) :: capnc ! cold and oceanic cloud particles / cm3 real(r8) :: capnw ! warm continental cloud particles / cm3 real(r8) :: kconst ! const for terminal velocity (stokes regime) real(r8) :: effc ! collection efficiency real(r8) :: alpha ! ratio of 3rd moment radius to 2nd real(r8) :: capc ! constant for autoconversion real(r8) :: convfw ! constant used for fall velocity calculation real(r8) :: cracw ! constant used for rain accreting water real(r8) :: critpr ! critical precip rate collection efficiency changes real(r8) :: ciautb ! coefficient of autoconversion of ice (1/s) #ifdef DEBUG integer, private,parameter :: nlook = 1 ! Number of points to examine integer, private :: ilook(nlook) ! Longitude index to examine integer, private :: latlook(nlook) ! Latitude index to examine integer, private :: lchnklook(nlook) ! Chunk index to examine integer, private :: icollook(nlook) ! Column index to examine #endif ! Private data real(r8), parameter :: unset_r8 = huge(1.0_r8) contains !=============================================================================== subroutine cldwat_readnl(nlfile) #ifndef WRF_PORT use namelist_utils, only: find_group_name use units, only: getunit, freeunit use mpishorthand #endif character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input #ifndef WRF_PORT ! Namelist variables real(r8) :: cldwat_icritw = unset_r8 ! icritw = threshold for autoconversion of warm ice real(r8) :: cldwat_icritc = unset_r8 ! icritc = threshold for autoconversion of cold ice real(r8) :: cldwat_conke = unset_r8 ! conke = tunable constant for evaporation of precip real(r8) :: cldwat_r3lcrit = unset_r8 ! r3lcrit = critical radius where liq conversion begins ! Local variables integer :: unitn, ierr character(len=*), parameter :: subname = 'cldwat_readnl' namelist /cldwat_nl/ cldwat_icritw, cldwat_icritc, cldwat_conke, cldwat_r3lcrit !----------------------------------------------------------------------------- if (masterproc) then unitn = getunit() open( unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'cldwat_nl', status=ierr) if (ierr == 0) then read(unitn, cldwat_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') end if end if close(unitn) call freeunit(unitn) ! set local variables icritw = cldwat_icritw icritc = cldwat_icritc conke = cldwat_conke r3lcrit = cldwat_r3lcrit end if #ifdef SPMD ! Broadcast namelist variables call mpibcast(icritw, 1, mpir8, 0, mpicom) call mpibcast(icritc, 1, mpir8, 0, mpicom) call mpibcast(conke, 1, mpir8, 0, mpicom) call mpibcast(r3lcrit, 1, mpir8, 0, mpicom) #endif #endif end subroutine cldwat_readnl !================================================================================================ subroutine cldwat_fice(ncol, t, fice, fsnow) ! ! Compute the fraction of the total cloud water which is in ice phase. ! The fraction depends on temperature only. ! This is the form that was used for radiation, the code came from cldefr originally ! ! Author: B. A. Boville Sept 10, 2002 ! modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection ) !----------------------------------------------------------------------- use physconst, only: tmelt implicit none ! Arguments integer, intent(in) :: ncol ! number of active columns real(r8), intent(in) :: t(pcols,pver) ! temperature real(r8), intent(out) :: fice(pcols,pver) ! Fractional ice content within cloud real(r8), intent(out) :: fsnow(pcols,pver) ! Fractional snow content for convection ! Local variables integer :: i,k ! loop indexes !----------------------------------------------------------------------- !BSINGH - This is a temporary fix for uninitialized tmax_* and tmin_* variables !This can be fixed by calling the inimc subroutine in physics init calls tmax_fice = tmelt - 10._r8 tmin_fice = tmax_fice - 30._r8 tmax_fsnow = tmelt tmin_fsnow = tmelt - 5._r8 ! Define fractional amount of cloud that is ice do k=1,pver do i=1,ncol ! If warmer than tmax then water phase if (t(i,k) > tmax_fice) then fice(i,k) = 0.0_r8 ! If colder than tmin then ice phase else if (t(i,k) < tmin_fice) then fice(i,k) = 1.0_r8 ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax else fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice) end if ! snow fraction partitioning ! If warmer than tmax then water phase if (t(i,k) > tmax_fsnow) then fsnow(i,k) = 0.0_r8 ! If colder than tmin then ice phase else if (t(i,k) < tmin_fsnow) then fsnow(i,k) = 1.0_r8 ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax else fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow) end if end do end do return end subroutine cldwat_fice #ifndef WRF_PORT subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox) !----------------------------------------------------------------------- ! ! Purpose: ! initialize constants for the prognostic condensate ! ! Author: P. Rasch, April 1997 ! !----------------------------------------------------------------------- use pmgrid, only: plev, plevp use dycore, only: dycore_is, get_resolution use hycoef, only: hypm use phys_control, only: phys_getopts integer k real(r8), intent(in) :: tmeltx real(r8), intent(in) :: rhonotx real(r8), intent(in) :: gravitx real(r8), intent(in) :: rh2ox #ifdef UNICOSMP real(r8) signgam ! variable required by cray gamma function external gamma #endif character(len=16) :: microp_scheme ! Get microphysics option call phys_getopts( microp_scheme_out = microp_scheme ) ! Set following for all physics packages tmax_fice = tmeltx - 10._r8 !! tmax_fice = tmeltx !! tmin_fice = tmax_fice - 20. tmin_fice = tmax_fice - 30._r8 tmax_fsnow = tmeltx tmin_fsnow = tmeltx - 5._r8 ! Set remaining for RK microphysics if( microp_scheme .eq. 'RK' ) then rhonot = rhonotx ! air density at surface (gm/cm3) gravit = gravitx rh2o = rh2ox rhos = .1_r8 ! assumed snow density (gm/cm3) rhow = 1._r8 ! water density rhoi = 1._r8 ! ice density esi = 1.0_r8 ! collection efficient for ice by snow esw = 0.1_r8 ! collection efficient for water by snow t0 = tmeltx ! approximate freezing temp cldmin = 0.02_r8 ! assumed minimum cloud amount small = 1.e-22_r8 ! a small number compared to unity c = 152.93_r8 ! constant for graupel like snow cm**(1-d)/s d = 0.25_r8 ! constant for graupel like snow nos = 3.e-2_r8 ! particles snow / cm**4 pi = 4._r8*atan(1.0_r8) prhonos = pi*rhos*nos thrpd = 3._r8 + d if (d==0.25_r8) then gam3pd = 2.549256966718531_r8 ! only right for d = 0.25 gam4pd = 8.285085141835282_r8 else #ifdef UNICOSMP call gamma(3._r8+d, signgam, gam3pd) gam3pd = sign(exp(gam3pd),signgam) call gamma(4._r8+d, signgam, gam4pd) gam4pd = sign(exp(gam4pd),signgam) write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd #ifdef WRF_PORT call wrf_message(iulog) #endif #else write(iulog,*) ' can only use d ne 0.25 on a cray ' #ifdef WRF_PORT call wrf_message(iulog) #endif stop #endif endif mcon01 = pi*nos*c*gam3pd/4._r8 mcon02 = 1._r8/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._r8))) mcon03 = -(0.5_r8+d/4._r8) mcon04 = 4._r8/(4._r8+d) mcon05 = (3+d)/(4+d) mcon06 = (3+d)/4._r8 mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06 mcon08 = -0.5_r8/(4._r8+d) ! find the level about 1mb, we wont do the microphysics above this level k1mb = 1 do k=1,pver-1 if (hypm(k) < 1.e2_r8 .and. hypm(k+1) >= 1.e2_r8) then if (1.e2_r8-hypm(k) < hypm(k+1)-1.e2_r8) then k1mb = k else k1mb = k + 1 end if goto 20 end if end do if (masterproc) then write(iulog,*)'inimc: model levels bracketing 1 mb not found' #ifdef WRF_PORT call wrf_message(iulog) #endif end if ! call endrun k1mb = 1 20 if( masterproc ) write(iulog,*)'inimc: model level nearest 1 mb is',k1mb,'which is',hypm(k1mb),'pascals' #ifdef WRF_PORT call wrf_message(iulog) #endif if( masterproc ) write(iulog,*) 'cloud water initialization by inimc complete ' #ifdef WRF_PORT call wrf_message(iulog) #endif ! Initialize parameters used by findmcnew capnw = 400._r8 ! warm continental cloud particles / cm3 capnc = 150._r8 ! cold and oceanic cloud particles / cm3 ! capnsi = 40._r8 ! sea ice cloud particles density / cm3 capnsi = 75._r8 ! sea ice cloud particles density / cm3 kconst = 1.18e6_r8 ! const for terminal velocity ! effc = 1._r8 ! autoconv collection efficiency following boucher 96 ! effc = .55*0.05_r8 ! autoconv collection efficiency following baker 93 effc = 0.55_r8 ! autoconv collection efficiency following tripoli and cotton ! effc = 0._r8 ! turn off warm-cloud autoconv alpha = 1.1_r8**4 capc = pi**(-.333_r8)*kconst*effc *(0.75_r8)**(1.333_r8)*alpha ! constant for autoconversion ! critical precip rate at which we assume the collector drops can change the ! drop size enough to enhance the auto-conversion process (mm/day) critpr = 0.5_r8 convfw = 1.94_r8*2.13_r8*sqrt(rhow*1000._r8*9.81_r8*2.7e-4_r8) ! liquid microphysics ! cracw = 6_r8 ! beheng cracw = .884_r8*sqrt(9.81_r8/(rhow*1000._r8*2.7e-4_r8)) ! tripoli and cotton ! ice microphysics ciautb = 5.e-4_r8 if ( masterproc ) then write(iulog,*)'tuning parameters cldwat: icritw',icritw,'icritc',icritc,'conke',conke,'r3lcrit',r3lcrit #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'tuning parameters cldwat: capnw',capnw,'capnc',capnc,'capnsi',capnsi,'kconst',kconst #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'tuning parameters cldwat: effc',effc,'alpha',alpha,'capc',capc #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*)'tuning parameters cldwat: critpr',critpr,'convfw',convfw,'cracw',cracw,'ciautb',ciautb #ifdef WRF_PORT call wrf_message(iulog) #endif endif endif return end subroutine inimc subroutine pcond (lchnk ,ncol , & tn ,ttend ,qn ,qtend ,omega , & cwat ,p ,pdel ,cldn ,fice , fsnow, & cme ,prodprec,prodsnow,evapprec,evapsnow,evapheat, prfzheat, & meltheat,precip ,snowab ,deltat ,fwaut , & fsaut ,fracw ,fsacw ,fsaci ,lctend , & rhdfda ,rhu00 ,seaicef, zi ,ice2pr, liq2pr, & liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio) !----------------------------------------------------------------------- ! ! Purpose: ! The public interface to the cloud water parameterization ! returns tendencies to water vapor, temperature and cloud water variables ! ! For basic method ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3 ! model climate using diagnosed and ! predicted condensate parameterizations, 1998, J. Clim., 11, ! pp1587---1614. ! ! For important modifications to improve the method of determining ! condensation/evaporation see Zhang et al (2001, in preparation) ! ! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson ! B. A. Boville (latent heat of fusion) !----------------------------------------------------------------------- use wv_saturation, only: vqsatd use cam_control_mod, only: nlvdry ! !--------------------------------------------------------------------- ! ! Input Arguments ! integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: fice(pcols,pver) ! fraction of cwat that is ice real(r8), intent(in) :: fsnow(pcols,pver) ! fraction of rain that freezes to snow real(r8), intent(in) :: cldn(pcols,pver) ! new value of cloud fraction (fraction) real(r8), intent(in) :: cwat(pcols,pver) ! cloud water (kg/kg) real(r8), intent(in) :: omega(pcols,pver) ! vert pressure vel (Pa/s) real(r8), intent(in) :: p(pcols,pver) ! pressure (K) real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness (Pa) real(r8), intent(in) :: qn(pcols,pver) ! new water vapor (kg/kg) real(r8), intent(in) :: qtend(pcols,pver) ! mixing ratio tend (kg/kg/s) real(r8), intent(in) :: tn(pcols,pver) ! new temperature (K) real(r8), intent(in) :: ttend(pcols,pver) ! temp tendencies (K/s) real(r8), intent(in) :: deltat ! time step to advance solution over real(r8), intent(in) :: lctend(pcols,pver) ! cloud liquid water tendencies ====wlin real(r8), intent(in) :: rhdfda(pcols,pver) ! dG(a)/da, rh=G(a), when rh>u00 ====wlin real(r8), intent(in) :: rhu00 (pcols,pver) ! Rhlim for cloud ====wlin real(r8), intent(in) :: seaicef(pcols) ! sea ice fraction (fraction) real(r8), intent(in) :: zi(pcols,pverp) ! layer interfaces (m) real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m) ! ! Output Arguments ! real(r8), intent(out) :: cme (pcols,pver) ! rate of cond-evap of condensate (1/s) real(r8), intent(out) :: prodprec(pcols,pver) ! rate of conversion of condensate to precip (1/s) real(r8), intent(out) :: evapprec(pcols,pver) ! rate of evaporation of falling precip (1/s) real(r8), intent(out) :: evapsnow(pcols,pver) ! rate of evaporation of falling snow (1/s) real(r8), intent(out) :: evapheat(pcols,pver) ! heating rate due to evaporation of precip (W/kg) real(r8), intent(out) :: prfzheat(pcols,pver) ! heating rate due to freezing of precip (W/kg) real(r8), intent(out) :: meltheat(pcols,pver) ! heating rate due to snow melt (W/kg) real(r8), intent(out) :: precip(pcols) ! rate of precipitation (kg / (m**2 * s)) real(r8), intent(out) :: snowab(pcols) ! rate of snow (kg / (m**2 * s)) real(r8), intent(out) :: ice2pr(pcols,pver) ! rate of conversion of ice to precip real(r8), intent(out) :: liq2pr(pcols,pver) ! rate of conversion of liquid to precip real(r8), intent(out) :: liq2snow(pcols,pver) ! rate of conversion of liquid to snow real(r8), intent(out) :: rkflxprc(pcols,pverp) ! grid-box mean RK flux_large_scale_cloud_rain+snow at interfaces (kg m^-2 s^-1) real(r8), intent(out) :: rkflxsnw(pcols,pverp) ! grid-box mean RK flux_large_scale_cloud_snow at interfaces (kg m^-2 s^-1) ! intent(out)s here for pcond to pass to stratiform.F90 to be addflded/outflded real(r8), intent(out) :: pracwo(pcols,pver) ! accretion of cloud water by rain (1/s) real(r8), intent(out) :: psacwo(pcols,pver) ! accretion of cloud water by snow (1/s) real(r8), intent(out) :: psacio(pcols,pver) ! accretion of cloud ice by snow (1/s) real(r8) nice2pr ! rate of conversion of ice to snow real(r8) nliq2pr ! rate of conversion of liquid to precip real(r8) nliq2snow ! rate of conversion of liquid to snow real(r8) prodsnow(pcols,pver) ! rate of production of snow ! ! Local workspace ! real(r8) :: precab(pcols) ! rate of precipitation (kg / (m**2 * s)) integer i ! work variable integer iter ! #iterations for precipitation calculation integer k ! work variable integer l ! work variable real(r8) cldm(pcols) ! mean cloud fraction over the time step real(r8) cldmax(pcols) ! max cloud fraction above real(r8) coef(pcols) ! conversion time scale for condensate to rain real(r8) cwm(pcols) ! cwat mixing ratio at midpoint of time step real(r8) cwn(pcols) ! cwat mixing ratio at end real(r8) denom ! work variable real(r8) dqsdt ! change in sat spec. hum. wrt temperature real(r8) es(pcols) ! sat. vapor pressure real(r8) fracw(pcols,pver) ! relative importance of collection of liquid by rain real(r8) fsaci(pcols,pver) ! relative importance of collection of ice by snow real(r8) fsacw(pcols,pver) ! relative importance of collection of liquid by snow real(r8) fsaut(pcols,pver) ! relative importance of ice auto conversion real(r8) fwaut(pcols,pver) ! relative importance of warm cloud autoconversion real(r8) gamma(pcols) ! d qs / dT real(r8) icwc(pcols) ! in-cloud water content (kg/kg) real(r8) mincld ! a small cloud fraction to avoid / zero real(r8) omeps ! 1 minus epsilon real(r8),parameter ::omsm=0.99999_r8 ! a number just less than unity (for rounding) real(r8) prprov(pcols) ! provisional value of precip at btm of layer real(r8) prtmp ! work variable real(r8) q(pcols,pver) ! mixing ratio before time step ignoring condensate real(r8) qs(pcols) ! spec. hum. of water vapor real(r8) qsn, esn ! work variable real(r8) qsp(pcols,pver) ! sat pt mixing ratio real(r8) qtl(pcols) ! tendency which would saturate the grid box in deltat real(r8) qtmp, ttmp ! work variable real(r8) relhum1(pcols) ! relative humidity real(r8) relhum(pcols) ! relative humidity !!$ real(r8) tc ! crit temp of transition to ice real(r8) t(pcols,pver) ! temp before time step ignoring condensate real(r8) tsp(pcols,pver) ! sat pt temperature real(r8) pol ! work variable real(r8) cdt ! work variable real(r8) wtthick ! work variable ! Extra local work space for cloud scheme modification real(r8) cpohl !Cp/Hlatv real(r8) hlocp !Hlatv/Cp real(r8) dto2 !0.5*deltat (delta=2.0*dt) real(r8) calpha(pcols) !alpha of new C - E scheme formulation real(r8) cbeta (pcols) !beta of new C - E scheme formulation real(r8) cbetah(pcols) !beta_hat at saturation portion real(r8) cgamma(pcols) !gamma of new C - E scheme formulation real(r8) cgamah(pcols) !gamma_hat at saturation portion real(r8) rcgama(pcols) !gamma/gamma_hat real(r8) csigma(pcols) !sigma of new C - E scheme formulation real(r8) cmec1 (pcols) !c1 of new C - E scheme formulation real(r8) cmec2 (pcols) !c2 of new C - E scheme formulation real(r8) cmec3 (pcols) !c3 of new C - E scheme formulation real(r8) cmec4 (pcols) !c4 of new C - E scheme formulation real(r8) cmeres(pcols) !residual cond of over-sat after cme and evapprec real(r8) ctmp !a scalar representation of cmeres real(r8) clrh2o ! Ratio of latvap to water vapor gas const real(r8) ice(pcols,pver) ! ice mixing ratio real(r8) liq(pcols,pver) ! liquid mixing ratio real(r8) rcwn(pcols,2,pver), rliq(pcols,2,pver), rice(pcols,2,pver) real(r8) cwnsave(pcols,2,pver), cmesave(pcols,2,pver) real(r8) prodprecsave(pcols,2,pver) logical error_found ! !------------------------------------------------------------ ! clrh2o = hlatv/rh2o ! Ratio of latvap to water vapor gas const omeps = 1.0_r8 - epsqs #ifdef PERGRO mincld = 1.e-4_r8 iter = 1 ! number of times to iterate the precipitation calculation #else mincld = 1.e-4_r8 iter = 2 #endif ! omsm = 0.99999 cpohl = cp/hlatv hlocp = hlatv/cp dto2=0.5_r8*deltat ! ! Constant for computing rate of evaporation of precipitation: ! !!$ conke = 1.e-5 !!$ conke = 1.e-6 ! ! initialize a few single level fields ! do i = 1,ncol precip(i) = 0.0_r8 precab(i) = 0.0_r8 snowab(i) = 0.0_r8 cldmax(i) = 0.0_r8 end do ! ! initialize multi-level fields ! do k = 1,pver do i = 1,ncol q(i,k) = qn(i,k) t(i,k) = tn(i,k) ! q(i,k)=qn(i,k)-qtend(i,k)*deltat ! t(i,k)=tn(i,k)-ttend(i,k)*deltat end do end do cme (:ncol,:) = 0._r8 evapprec(:ncol,:) = 0._r8 prodprec(:ncol,:) = 0._r8 evapsnow(:ncol,:) = 0._r8 prodsnow(:ncol,:) = 0._r8 evapheat(:ncol,:) = 0._r8 meltheat(:ncol,:) = 0._r8 prfzheat(:ncol,:) = 0._r8 ice2pr(:ncol,:) = 0._r8 liq2pr(:ncol,:) = 0._r8 liq2snow(:ncol,:) = 0._r8 fwaut(:ncol,:) = 0._r8 fsaut(:ncol,:) = 0._r8 fracw(:ncol,:) = 0._r8 fsacw(:ncol,:) = 0._r8 fsaci(:ncol,:) = 0._r8 rkflxprc(:ncol,:) = 0._r8 rkflxsnw(:ncol,:) = 0._r8 pracwo(:ncol,:) = 0._r8 psacwo(:ncol,:) = 0._r8 psacio(:ncol,:) = 0._r8 ! ! find the wet bulb temp and saturation value ! for the provisional t and q without condensation ! call findsp (lchnk, ncol, qn, tn, p, tsp, qsp) do 800 k = k1mb,pver call vqsatd (t(1,k), p(1,k), es, qs, gamma, ncol) do i = 1,ncol relhum(i) = q(i,k)/qs(i) ! cldm(i) = max(cldn(i,k),mincld) ! ! the max cloud fraction above this level ! cldmax(i) = max(cldmax(i), cldm(i)) ! define the coefficients for C - E calculation calpha(i) = 1.0_r8/qs(i) cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl cbetah(i) = 1.0_r8/qs(i)*gamma(i)*cpohl cgamma(i) = calpha(i)+hlatv*cbeta(i)/cp cgamah(i) = calpha(i)+hlatv*cbetah(i)/cp rcgama(i) = cgamma(i)/cgamah(i) if(cldm(i) > mincld) then icwc(i) = max(0._r8,cwat(i,k)/cldm(i)) else icwc(i) = 0.0_r8 endif !PJR the above logic give zero icwc with nonzero cwat, dont like it! !PJR generates problems with csigma !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds ! icwc(i) = max(1.e-8_r8,cwat(i,k)/cldm(i)) ! ! initial guess of evaporation, will be updated within iteration ! evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) & *(1._r8 - min(relhum(i),1._r8)) ! ! zero cmeres before iteration for each level ! cmeres(i)=0.0_r8 end do do i = 1,ncol ! ! fractions of ice at this level ! !!$ tc = t(i,k) - t0 !!$ fice(i,k) = max(0._r8,min(-tc*0.05,1.0_r8)) ! ! calculate the cooling due to a phase change of the rainwater ! from above ! if (t(i,k) >= t0) then meltheat(i,k) = -hlatf * snowab(i) * gravit/pdel(i,k) snowab(i) = 0._r8 else meltheat(i,k) = 0._r8 endif end do ! ! calculate cme and formation of precip. ! ! The cloud microphysics is highly nonlinear and coupled with cme ! Both rain processes and cme are calculated iteratively. ! do 100 l = 1,iter do i = 1,ncol ! ! calculation of cme has 4 scenarios ! ================================== ! if(relhum(i) > rhu00(i,k)) then ! 1. whole grid saturation ! ======================== if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i) ! 2. fractional saturation ! ======================== else if (rhdfda(i,k) .eq. 0. .and. icwc(i).eq.0.) then write (iulog,*) ' cldwat.F90: empty rh cloud ', i, k, lchnk #ifdef WRF_PORT call wrf_message(iulog) #endif write (iulog,*) ' relhum, iter ', relhum(i), l, rhu00(i,k), cldm(i), cldn(i,k) #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun () endif csigma(i) = 1.0_r8/(rhdfda(i,k)+cgamma(i)*icwc(i)) cmec1(i) = (1.0_r8-cldm(i))*csigma(i)*rhdfda(i,k) cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_r8-rcgama(i)*cldm(i))* & csigma(i)*calpha(i)*icwc(i) cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) + & (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i) cmec4(i) = csigma(i)*cgamma(i)*icwc(i) ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k) & -cmec3(i)*ttend(i,k) + cmec4(i)*evapprec(i,k) endif ! 3. when rh < rhu00, evaporate existing cloud water ! ================================================== else if(cwat(i,k) > 0.0_r8)then ! liquid water should be evaporated but not to exceed ! saturation point. if qn > qsp, not to evaporate cwat cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat ! 4. no condensation nor evaporation ! ================================== else cme(i,k)=0.0_r8 endif end do !end loop for cme update ! Because of the finite time step, ! place a bound here not to exceed wet bulb point ! and not to evaporate more than available water ! do i = 1, ncol qtmp = qn(i,k) - cme(i,k)*deltat ! possibilities to have qtmp > qsp ! ! 1. if qn > qs(tn), it condenses; ! if after applying cme, qtmp > qsp, more condensation is applied. ! ! 2. if qn < qs, evaporation should not exceed qsp, if(qtmp > qsp(i,k)) then cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat endif ! ! if net evaporation, it should not exceed available cwat ! if(cme(i,k) < -cwat(i,k)/deltat) & cme(i,k) = -cwat(i,k)/deltat ! ! addition of residual condensation from previous step of iteration ! cme(i,k) = cme(i,k) + cmeres(i) end do ! limit cme for roundoff errors do i = 1, ncol cme(i,k) = cme(i,k)*omsm end do do i = 1,ncol ! ! as a safe limit, condensation should not reduce grid mean rh below rhu00 ! if(cme(i,k) > 0.0_r8 .and. relhum(i) > rhu00(i,k) ) & cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu00(i,k))/deltat) ! ! initial guess for cwm (mean cloud water over time step) if 1st iteration ! if(l < 2) then cwm(i) = max(cwat(i,k)+cme(i,k)*dto2, 0._r8) endif enddo ! provisional precipitation falling through model layer do i = 1,ncol !!$ prprov(i) = precab(i) + prodprec(i,k)*pdel(i,k)/gravit ! rain produced in this layer not too effective in collection process wtthick = max(0._r8,min(0.5_r8,((zi(i,k)-zi(i,k+1))/1000._r8)**2)) prprov(i) = precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit end do ! calculate conversion of condensate to precipitation by cloud microphysics call findmcnew (lchnk ,ncol , & k ,prprov ,snowab, t ,p , & cwm ,cldm ,cldmax ,fice(1,k),coef , & fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), & seaicef, snowh, pracwo(1,k), psacwo(1,k), psacio(1,k)) ! ! calculate the precip rate ! error_found = .false. do i = 1,ncol if (cldm(i) > 0) then ! ! first predict the cloud water ! cdt = coef(i)*deltat if (cdt > 0.01_r8) then pol = cme(i,k)/coef(i) ! production over loss cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol) else cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt)) endif ! ! now back out the tendency of net rain production ! prodprec(i,k) = max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat) else prodprec(i,k) = 0.0_r8 cwn(i) = 0._r8 endif ! provisional calculation of conversion terms ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k)) liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k)) !old liq2snow(i,k) = prodprec(i,k)*fsacw(i,k) ! revision suggested by Jim McCaa ! it controls the amount of snow hitting the sfc ! by forcing a lot of conversion of cloud liquid to snow phase ! it might be better done later by an explicit representation of ! rain accreting ice (and freezing), or by an explicit freezing of raindrops liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k)) ! bounds nice2pr = min(ice2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*fice(i,k)/deltat) nliq2pr = min(liq2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*(1._r8-fice(i,k))/deltat) ! write(iulog,*) ' prodprec ', i, k, prodprec(i,k) ! write(iulog,*) ' nliq2pr, nice2pr ', nliq2pr, nice2pr if (liq2pr(i,k).ne.0._r8) then nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k) ! correction else nliq2snow = liq2snow(i,k) endif ! avoid roundoff problems generating negatives nliq2snow = nliq2snow*omsm nliq2pr = nliq2pr*omsm nice2pr = nice2pr*omsm ! final estimates of conversion to precip and snow prodprec(i,k) = (nliq2pr + nice2pr) prodsnow(i,k) = (nice2pr + nliq2snow) rcwn(i,l,k) = cwat(i,k) + (cme(i,k)- prodprec(i,k))*deltat rliq(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)*(1._r8-fice(i,k)) - nliq2pr * deltat rice(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)* fice(i,k) - nice2pr *deltat ! Save for sanity check later... ! Putting sanity checks inside loops 100 and 800 screws up the ! IBM compiler for reasons as yet unknown. TBH cwnsave(i,l,k) = cwn(i) cmesave(i,l,k) = cme(i,k) prodprecsave(i,l,k) = prodprec(i,k) ! End of save for sanity check later... ! final version of condensate to precip terms liq2pr(i,k) = nliq2pr liq2snow(i,k) = nliq2snow ice2pr(i,k) = nice2pr cwn(i) = rcwn(i,l,k) ! ! update any remaining provisional values ! cwm(i) = (cwn(i) + cwat(i,k))*0.5_r8 ! ! update in cloud water ! if(cldm(i) > mincld) then icwc(i) = cwm(i)/cldm(i) else icwc(i) = 0.0_r8 endif !PJR the above logic give zero icwc with nonzero cwat, dont like it! !PJR generates problems with csigma !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds ! icwc(i) = max(1.e-8_r8,cwm(i)/cldm(i)) end do ! end of do i = 1,ncol ! ! calculate provisional value of cloud water for ! evaporation of precipitate (evapprec) calculation ! do i = 1,ncol qtmp = qn(i,k) - cme(i,k)*deltat ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k) & + (hlatv + hlatf*fice(i,k)) * cme(i,k) ) esn = estblf(ttmp) qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8) qtl(i) = max((qsn - qtmp)/deltat,0._r8) relhum1(i) = qtmp/qsn end do ! do i = 1,ncol #ifdef PERGRO evapprec(i,k) = conke*(1._r8 - max(cldm(i),mincld))* & sqrt(precab(i))*(1._r8 - min(relhum1(i),1._r8)) #else evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) & *(1._r8 - min(relhum1(i),1._r8)) #endif ! ! limit the evaporation to the amount which is entering the box ! or saturates the box ! prtmp = precab(i)*gravit/pdel(i,k) evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm #ifdef PERGRO ! zeroing needed for pert growth evapprec(i,k) = 0._r8 #endif ! ! Partition evaporation of precipitate between rain and snow using ! the fraction of snow falling into the box. Determine the heating ! due to evaporation. Note that evaporation is positive (loss of precip, ! gain of vapor) and that heating is negative. if (evapprec(i,k) > 0._r8) then evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i) evapheat(i,k) = -hlatv * evapprec(i,k) - hlatf * evapsnow(i,k) else evapsnow(i,k) = 0._r8 evapheat(i,k) = 0._r8 end if ! Account for the latent heat of fusion for liquid drops collected by falling snow prfzheat(i,k) = hlatf * liq2snow(i,k) end do ! now remove the residual of any over-saturation. Normally, ! the oversaturated water vapor should have been removed by ! cme formulation plus constraints by wet bulb tsp/qsp ! as computed above. However, because of non-linearity, ! addition of (cme-evapprec) to update t and q may still cause ! a very small amount of over saturation. It is called a ! residual of over-saturation because theoretically, cme ! should have taken care of all of large scale condensation. ! do i = 1,ncol qtmp = qn(i,k)-(cme(i,k)-evapprec(i,k))*deltat ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k) & + (hlatv + hlatf*fice(i,k)) * cme(i,k) ) esn = estblf(ttmp) qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8) ! !Upper stratosphere and mesosphere, qsn calculated !above may be negative. Here just to skip it instead !of resetting it to 1 as in aqsat ! if(qtmp > qsn .and. qsn > 0) then !calculate dqsdt, a more precise calculation !which taking into account different range of T !can be found in aqsatd.F. Here follows !cond.F to calculate it. ! denom = (p(i,k)-omeps*esn)*ttmp*ttmp dqsdt = clrh2o*qsn*p(i,k)/denom ! !now extra condensation to bring air to just saturation ! ctmp = (qtmp-qsn)/(1._r8+hlocp*dqsdt)/deltat cme(i,k) = cme(i,k)+ctmp ! ! save residual on cmeres to addtion to cme on entering next iteration ! cme exit here contain the residual but overrided if back to iteration ! cmeres(i) = ctmp else cmeres(i) = 0.0_r8 endif end do 100 continue ! end of do l = 1,iter ! ! precipitation ! do i = 1,ncol precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k)) precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k)) if(precab(i).lt.0._r8) precab(i)=0._r8 ! snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k)) snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k)) ! If temperature above freezing, all precip is rain flux. if temperature below freezing, all precip is snow flux. rkflxprc(i,k+1) = precab(i) !! making this consistent with other precip fluxes. prc = rain + snow !!rkflxprc(i,k+1) = precab(i) - snowab(i) rkflxsnw(i,k+1) = snowab(i) !!$ if ((precab(i)) < 1.e-10) then !!$ precab(i) = 0. !!$ snowab(i) = 0. !!$ endif end do 800 continue ! level loop (k=1,pver) ! begin sanity checks error_found = .false. do k = k1mb,pver do l = 1,iter do i = 1,ncol if (abs(rcwn(i,l,k)).lt.1.e-300_r8) rcwn(i,l,k) = 0._r8 if (abs(rliq(i,l,k)).lt.1.e-300_r8) rliq(i,l,k) = 0._r8 if (abs(rice(i,l,k)).lt.1.e-300_r8) rice(i,l,k) = 0._r8 if (rcwn(i,l,k).lt.0._r8) error_found = .true. if (rliq(i,l,k).lt.0._r8) error_found = .true. if (rice(i,l,k).lt.0._r8) error_found = .true. enddo enddo enddo if (error_found) then do k = k1mb,pver do l = 1,iter do i = 1,ncol if (rcwn(i,l,k).lt.0._r8) then write(iulog,*) ' prob with neg rcwn1 ', rcwn(i,l,k), & cwnsave(i,l,k) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' cwat, cme*deltat, prodprec*deltat ', & cwat(i,k), cmesave(i,l,k)*deltat, & prodprecsave(i,l,k)*deltat, & (cmesave(i,l,k)-prodprecsave(i,l,k))*deltat #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun('PCOND') endif if (rliq(i,l,k).lt.0._r8) then write(iulog,*) ' prob with neg rliq1 ', rliq(i,l,k) #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun('PCOND') endif if (rice(i,l,k).lt.0._r8) then write(iulog,*) ' prob with neg rice ', rice(i,l,k) #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun('PCOND') endif enddo enddo enddo end if ! end sanity checks return end subroutine pcond !############################################################################## subroutine findmcnew (lchnk ,ncol , & k ,precab ,snowab, t ,p , & cwm ,cldm ,cldmax ,fice ,coef , & fwaut ,fsaut ,fracw ,fsacw ,fsaci , & seaicef ,snowh, pracwo, psacwo, psacio ) !----------------------------------------------------------------------- ! ! Purpose: ! calculate the conversion of condensate to precipitate ! ! Method: ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3 ! model climate using diagnosed and ! predicted condensate parameterizations, 1998, J. Clim., 11, ! pp1587---1614. ! ! Author: P. Rasch ! !----------------------------------------------------------------------- use phys_grid, only: get_rlat_all_p use comsrf, only: landm ! ! input args ! integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns integer, intent(in) :: k ! level index real(r8), intent(in) :: precab(pcols) ! rate of precipitation from above (kg / (m**2 * s)) real(r8), intent(in) :: t(pcols,pver) ! temperature (K) real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa) real(r8), intent(in) :: cldm(pcols) ! cloud fraction real(r8), intent(in) :: cldmax(pcols) ! max cloud fraction above this level real(r8), intent(in) :: cwm(pcols) ! condensate mixing ratio (kg/kg) real(r8), intent(in) :: fice(pcols) ! fraction of cwat that is ice real(r8), intent(in) :: seaicef(pcols) ! sea ice fraction real(r8), intent(in) :: snowab(pcols) ! rate of snow from above (kg / (m**2 * s)) real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m) ! output arguments real(r8), intent(out) :: coef(pcols) ! conversion rate (1/s) real(r8), intent(out) :: fwaut(pcols) ! relative importance of liquid autoconversion (a diagnostic) real(r8), intent(out) :: fsaut(pcols) ! relative importance of ice autoconversion (a diagnostic) real(r8), intent(out) :: fracw(pcols) ! relative importance of rain accreting liquid (a diagnostic) real(r8), intent(out) :: fsacw(pcols) ! relative importance of snow accreting liquid (a diagnostic) real(r8), intent(out) :: fsaci(pcols) ! relative importance of snow accreting ice (a diagnostic) real(r8), intent(out) :: pracwo(pcols) ! accretion of cloud water by rain (1/s) real(r8), intent(out) :: psacwo(pcols) ! accretion of cloud water by snow (1/s) real(r8), intent(out) :: psacio(pcols) ! accretion of cloud ice by snow (1/s) ! work variables integer i integer ii integer ind(pcols) integer ncols real(r8), parameter :: degrad = 57.296_r8 ! divide by this to convert degrees to radians real(r8) capn ! local cloud particles / cm3 real(r8) capnoice ! local cloud particles when not over sea ice / cm3 real(r8) ciaut ! coefficient of autoconversion of ice (1/s) real(r8) cldloc(pcols) ! non-zero amount of cloud real(r8) cldpr(pcols) ! assumed cloudy volume occupied by rain and cloud real(r8) con1 ! work constant real(r8) con2 ! work constant real(r8) csacx ! constant used for snow accreting liquid or ice !!$ real(r8) dtice ! interval for transition from liquid to ice real(r8) icemr(pcols) ! in-cloud ice mixing ratio real(r8) icrit ! threshold for autoconversion of ice real(r8) liqmr(pcols) ! in-cloud liquid water mixing ratio real(r8) pracw ! rate of rain accreting water real(r8) prlloc(pcols) ! local rain flux in mm/day real(r8) prscgs(pcols) ! local snow amount in cgs units real(r8) psaci ! rate of collection of ice by snow (lin et al 1983) real(r8) psacw ! rate of collection of liquid by snow (lin et al 1983) real(r8) psaut ! rate of autoconversion of ice condensate real(r8) ptot ! total rate of conversion real(r8) pwaut ! rate of autoconversion of liquid condensate real(r8) r3l ! volume radius real(r8) rainmr(pcols) ! in-cloud rain mixing ratio real(r8) rat1 ! work constant real(r8) rat2 ! work constant !!$ real(r8) rdtice ! recipricol of dtice real(r8) rho(pcols) ! density (mks units) real(r8) rhocgs ! density (cgs units) real(r8) rlat(pcols) ! latitude (radians) real(r8) snowfr ! fraction of precipate existing as snow real(r8) totmr(pcols) ! in-cloud total condensate mixing ratio real(r8) vfallw ! fall speed of precipitate as liquid real(r8) wp ! weight factor used in calculating pressure dep of autoconversion real(r8) wsi ! weight factor for sea ice real(r8) wt ! fraction of ice real(r8) wland ! fraction of land ! real(r8) csaci ! real(r8) csacw ! real(r8) cwaut ! real(r8) efact ! real(r8) lamdas ! real(r8) lcrit ! real(r8) rcwm ! real(r8) r3lc2 ! real(r8) snowmr(pcols) ! real(r8) vfalls real(8) ftot ! inline statement functions real(r8) heavy, heavym, a1, a2, heavyp, heavymp heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2)) ! heavyside function heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2)) ! modified heavyside function ! ! New heavyside functions to perhaps address error growth problems ! heavyp(a1,a2) = a1/(a2+a1+1.e-36_r8) heavymp(a1,a2) = (a1+0.01_r8*a2)/(a2+a1+1.e-36_r8) ! ! find all the points where we need to do the microphysics ! and set the output variables to zero ! ncols = 0 do i = 1,ncol coef(i) = 0._r8 fwaut(i) = 0._r8 fsaut(i) = 0._r8 fracw(i) = 0._r8 fsacw(i) = 0._r8 fsaci(i) = 0._r8 liqmr(i) = 0._r8 rainmr(i) = 0._r8 if (cwm(i) > 1.e-20_r8) then ncols = ncols + 1 ind(ncols) = i endif end do !cdir nodep !DIR$ CONCURRENT do ii = 1,ncols i = ind(ii) ! ! the local cloudiness at this level ! cldloc(i) = max(cldmin,cldm(i)) ! ! a weighted mean between max cloudiness above, and this layer ! cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_r8) ! ! decompose the suspended condensate into ! an incloud liquid and ice phase component ! totmr(i) = cwm(i)/cldloc(i) icemr(i) = totmr(i)*fice(i) liqmr(i) = totmr(i)*(1._r8-fice(i)) ! ! density ! rho(i) = p(i,k)/(287._r8*t(i,k)) rhocgs = rho(i)*1.e-3_r8 ! density in cgs units ! ! decompose the precipitate into a liquid and ice phase ! if (t(i,k) > t0) then vfallw = convfw/sqrt(rho(i)) rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i)) snowfr = 0 ! snowmr(i) else snowfr = 1 rainmr(i) = 0._r8 endif ! rainmr(i) = (precab(i)-snowab(i))/(rho(i)*vfallw*cldpr(i)) ! ! local snow amount in cgs units ! prscgs(i) = precab(i)/cldpr(i)*0.1_r8*snowfr ! prscgs(i) = snowab(i)/cldpr(i)*0.1 ! ! local rain amount in mm/day ! prlloc(i) = precab(i)*86400._r8/cldpr(i) end do con1 = 1._r8/(1.333_r8*pi)**0.333_r8 * 0.01_r8 ! meters ! ! calculate the conversion terms ! call get_rlat_all_p(lchnk, ncol, rlat) !cdir nodep !DIR$ CONCURRENT do ii = 1,ncols i = ind(ii) rhocgs = rho(i)*1.e-3_r8 ! density in cgs units ! ! exponential temperature factor ! ! efact = exp(0.025*(t(i,k)-t0)) ! ! some temperature dependent constants ! !!$ wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice)) wt = fice(i) icrit = icritc*wt + icritw*(1-wt) ! ! jrm Reworked droplet number concentration algorithm ! Start with pressure-dependent value appropriate for continental air ! Note: reltab has a temperature dependence here capn = capnw + (capnc-capnw) * min(1._r8,max(0._r8,1.0_r8-(p(i,k)-0.8_r8*p(i,pver))/(0.2_r8*p(i,pver)))) ! Modify for snow depth over land capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,snowh(i)*10._r8)) ! Ramp between polluted value over land to clean value over ocean. capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,1.0_r8-landm(i,lchnk))) ! Ramp between the resultant value and a sea ice value in the presence of ice. capn = capn + (capnsi-capn) * min(1.0_r8,max(0.0_r8,seaicef(i))) ! end jrm ! #ifdef DEBUG2 if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then if (i == ilook(1)) then write(iulog,*) ' findmcnew: lat, k, seaicef, landm, wp, capnoice, capn ', & lat(i), k, seaicef(i), landm(i,lat(i)), wp, capnoice, capn #ifdef WRF_PORT call wrf_message(iulog) #endif endif endif #endif ! ! useful terms in following calculations ! rat1 = rhocgs/rhow rat2 = liqmr(i)/capn con2 = (rat1*rat2)**0.333_r8 ! ! volume radius ! ! r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters r3l = con1*con2 ! ! critical threshold for autoconversion if modified for mixed phase ! clouds to mimic a bergeron findeisen process ! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i))) ! ! autoconversion of liquid ! ! cwaut = 2.e-4 ! cwaut = 1.e-3 ! lcrit = 2.e-4 ! lcrit = 5.e-4 ! pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut ! ! pwaut is following tripoli and cotton (and many others) ! we reduce the autoconversion below critpr, because these are regions where ! the drop size distribution is likely to imply much smaller collector drops than ! those relevant for a cloud distribution corresponding to the value of effc = 0.55 ! suggested by cotton (see austin 1995 JAS, baker 1993) ! easy to follow form ! pwaut = capc*liqmr(i)**2*rhocgs/rhow ! $ *(liqmr(i)*rhocgs/(rhow*capn))**(.333) ! $ *heavy(r3l,r3lcrit) ! $ *max(0.10_r8,min(1._r8,prlloc(i)/critpr)) ! somewhat faster form #define HEAVYNEW #ifdef HEAVYNEW !#ifdef PERGRO pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * & max(0.10_r8,min(1._r8,prlloc(i)/critpr)) #else pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* & max(0.10_r8,min(1._r8,prlloc(i)/critpr)) #endif ! ! autoconversion of ice ! ! ciaut = ciautb*efact ciaut = ciautb ! psaut = capc*totmr(i)**2*rhocgs/rhoi ! $ *(totmr(i)*rhocgs/(rhoi*capn))**(.333) ! ! autoconversion of ice condensate ! #ifdef PERGRO psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut #else psaut = max(0._r8,icemr(i)-icrit)*ciaut #endif ! ! collection of liquid by rain ! ! pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994) pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton) pracwo(i)=pracw !! pracw = 0. ! ! the following lines calculate the slope parameter and snow mixing ratio ! from the precip rate using the equations found in lin et al 83 ! in the most natural form, but it is expensive, so after some tedious ! algebraic manipulation you can use the cheaper form found below ! vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs) ! $ *0.01 ! convert from cm/s to m/s ! snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i)) ! snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04 ! lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25 ! csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd) ! ! coefficient for collection by snow independent of phase ! csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05 ! ! collection of liquid by snow (lin et al 1983) ! psacw = csacx*liqmr(i)*esw #ifdef PERGRO ! this is necessary for pergro psacw = 0._r8 #endif psacwo(i)=psacw ! ! collection of ice by snow (lin et al 1983) ! psaci = csacx*icemr(i)*esi ! psacio(i)=psaci ! total conversion of condensate to precipitate ! ptot = pwaut + psaut + pracw + psacw + psaci ! ! the recipricol of cloud water amnt (or zero if no cloud water) ! ! rcwm = totmr(i)/(max(totmr(i),small)**2) ! ! turn the tendency back into a loss rate (1/seconds) ! if (totmr(i) > 0._r8) then coef(i) = ptot/totmr(i) else coef(i) = 0._r8 endif if (ptot.gt.0._r8) then fwaut(i) = pwaut/ptot fsaut(i) = psaut/ptot fracw(i) = pracw/ptot fsacw(i) = psacw/ptot fsaci(i) = psaci/ptot else fwaut(i) = 0._r8 fsaut(i) = 0._r8 fracw(i) = 0._r8 fsacw(i) = 0._r8 fsaci(i) = 0._r8 endif ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i) ! if (abs(ftot-1._r8).gt.1.e-14_r8.and.ftot.ne.0._r8) then ! write(iulog,*) ' something is wrong in findmcnew ', ftot, & ! fwaut(i),fsaut(i),fracw(i),fsacw(i),fsaci(i) ! write(iulog,*) ' unscaled ', ptot, & ! pwaut,psaut,pracw,psacw,psaci ! write(iulog,*) ' totmr, liqmr, icemr ', totmr(i), liqmr(i), icemr(i) ! call endrun() ! endif end do #ifdef DEBUG i = icollook(nlook) if (lchnk == lchnklook(nlook) ) then write(iulog,*) write(iulog,*) '------', k, i, lchnk write(iulog,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4_r8 write(iulog,*) ' frac: waut,saut,racw,sacw,saci ', & fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i) endif #endif return end subroutine findmcnew !############################################################################## subroutine findsp (lchnk, ncol, q, t, p, tsp, qsp) !----------------------------------------------------------------------- ! ! Purpose: ! find the wet bulb temperature for a given t and q ! in a longitude height section ! wet bulb temp is the temperature and spec humidity that is ! just saturated and has the same enthalpy ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q ! ! Method: ! a Newton method is used ! first guess uses an algorithm provided by John Petch from the UKMO ! we exclude points where the physical situation is unrealistic ! e.g. where the temperature is outside the range of validity for the ! saturation vapor pressure, or where the water vapor pressure ! exceeds the ambient pressure, or the saturation specific humidity is ! unrealistic ! ! Author: P. Rasch ! !----------------------------------------------------------------------- ! ! input arguments ! integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: q(pcols,pver) ! water vapor (kg/kg) real(r8), intent(in) :: t(pcols,pver) ! temperature (K) real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa) ! ! output arguments ! real(r8), intent(out) :: tsp(pcols,pver) ! saturation temp (K) real(r8), intent(out) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg) ! ! local variables ! integer i ! work variable integer k ! work variable logical lflg ! work variable integer iter ! work variable integer l ! work variable logical :: error_found real(r8) omeps ! 1 minus epsilon real(r8) trinv ! work variable real(r8) es ! sat. vapor pressure real(r8) desdt ! change in sat vap pressure wrt temperature ! real(r8) desdp ! change in sat vap pressure wrt pressure real(r8) dqsdt ! change in sat spec. hum. wrt temperature real(r8) dgdt ! work variable real(r8) g ! work variable real(r8) weight(pcols) ! work variable real(r8) hlatsb ! (sublimation) real(r8) hlatvp ! (vaporization) real(r8) hltalt(pcols,pver) ! lat. heat. of vap. real(r8) tterm ! work var. real(r8) qs ! spec. hum. of water vapor real(r8) tc ! crit temp of transition to ice ! work variables real(r8) t1, q1, dt, dq real(r8) dtm, dqm real(r8) qvd, a1, tmp real(r8) rair real(r8) r1b, c1, c2, c3 real(r8) denom real(r8) dttol real(r8) dqtol integer doit(pcols) real(r8) enin(pcols), enout(pcols) real(r8) tlim(pcols) omeps = 1.0_r8 - epsqs trinv = 1.0_r8/ttrice a1 = 7.5_r8*log(10._r8) rair = 287.04_r8 c3 = rair*a1/cp dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration ! tmin = 173.16 ! the coldest temperature we can deal with ! ! max number of times to iterate the calculation iter = 8 ! do k = k1mb,pver ! ! first guess on the wet bulb temperature ! do i = 1,ncol #ifdef DEBUG if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then write(iulog,*) ' ' #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k) #ifdef WRF_PORT call wrf_message(iulog) #endif endif #endif ! limit the temperature range to that relevant to the sat vap pres tables #if ( ! defined WACCM_PHYS ) tlim(i) = min(max(t(i,k),173._r8),373._r8) #else tlim(i) = min(max(t(i,k),128._r8),373._r8) #endif es = estblf(tlim(i)) denom = p(i,k) - omeps*es qs = epsqs*es/denom doit(i) = 0 enout(i) = 1._r8 ! make sure a meaningful calculation is possible if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then ! ! Saturation specific humidity ! qs = min(epsqs*es/denom,1._r8) ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! tc = tlim(i) - t0 lflg = (tc >= -ttrice .and. tc < 0.0_r8) weight(i) = min(-tc*trinv,1.0_r8) hlatsb = hlatv + weight(i)*hlatf hlatvp = hlatv - 2369.0_r8*tc if (tlim(i) < t0) then hltalt(i,k) = hlatsb else hltalt(i,k) = hlatvp end if enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k) ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch) tmp = q(i,k) - qs c1 = hltalt(i,k)*c3 c2 = (tlim(i) + 36._r8)**2 r1b = c2/(c2 + c1*qs) qvd = r1b*tmp tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd) #ifdef DEBUG if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then write(iulog,*) ' relative humidity ', q(i,k)/qs #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' first guess ', tsp(i,k) #ifdef WRF_PORT call wrf_message(iulog) #endif endif #endif es = estblf(tsp(i,k)) qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8) else doit(i) = 1 tsp(i,k) = tlim(i) qsp(i,k) = q(i,k) enin(i) = 1._r8 endif end do ! end do i ! ! now iterate on first guess ! do l = 1, iter dtm = 0 dqm = 0 do i = 1,ncol if (doit(i) == 0) then es = estblf(tsp(i,k)) ! ! Saturation specific humidity ! qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8) ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! tc = tsp(i,k) - t0 lflg = (tc >= -ttrice .and. tc < 0.0_r8) weight(i) = min(-tc*trinv,1.0_r8) hlatsb = hlatv + weight(i)*hlatf hlatvp = hlatv - 2369.0_r8*tc if (tsp(i,k) < t0) then hltalt(i,k) = hlatsb else hltalt(i,k) = hlatvp end if if (lflg) then tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5)))) else tterm = 0.0_r8 end if desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt ! g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k) g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)) dgdt = -(cp + hltalt(i,k)*dqsdt) t1 = tsp(i,k) - g/dgdt dt = abs(t1 - tsp(i,k))/t1 tsp(i,k) = max(t1,tmin) es = estblf(tsp(i,k)) q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8) dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8) qsp(i,k) = q1 #ifdef DEBUG if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then write(iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g #ifdef WRF_PORT call wrf_message(iulog) #endif endif #endif dtm = max(dtm,dt) dqm = max(dqm,dq) ! if converged at this point, exclude it from more iterations if (dt < dttol .and. dq < dqtol) then doit(i) = 2 endif enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k) ! bail out if we are too near the end of temp range #if ( ! defined WACCM_PHYS ) if (tsp(i,k) < 174.16_r8) then #else if (tsp(i,k) < 130.16_r8) then #endif doit(i) = 4 endif else endif end do ! do i = 1,ncol if (dtm < dttol .and. dqm < dqtol) then go to 10 endif end do ! do l = 1,iter 10 continue error_found = .false. if (dtm > dttol .or. dqm > dqtol) then do i = 1,ncol if (doit(i) == 0) error_found = .true. end do if (error_found) then do i = 1,ncol if (doit(i) == 0) then write(iulog,*) ' findsp not converging at point i, k ', i, k #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i) #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun ('FINDSP') endif end do endif endif do i = 1,ncol if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then error_found = .true. endif end do if (error_found) then do i = 1,ncol if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then write(iulog,*) ' the enthalpy is not conserved for point ', & i, k, enin(i), enout(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i) #ifdef WRF_PORT call wrf_message(iulog) #endif write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i) #ifdef WRF_PORT call wrf_message(iulog) #endif call endrun ('FINDSP') endif end do endif end do ! level loop (k=1,pver) return end subroutine findsp #endif end module cldwat