LCOV - code coverage report
Current view: | top level - physics/cam - cldwat.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 8 | 338 | 2.4 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 5 | 20.0 %

              Line data    Source code
           1             : #undef DEBUG
           2             : module cldwat
           3             : !-----------------------------------------------------------------------
           4             : !
           5             : ! Purpose: Prognostic cloud water data and methods.
           6             : !
           7             : ! Public interfaces:
           8             : !
           9             : ! inimc -- Initialize constants
          10             : ! pcond -- Calculate prognostic condensate
          11             : !
          12             : ! Author: P. Rasch, with Modifications by Minghua Zhang
          13             : ! January 2010, modified by J. Kay to add precip fluxes for COSP simulator
          14             : !
          15             : !-----------------------------------------------------------------------
          16             :    use shr_kind_mod,   only: r8 => shr_kind_r8
          17             :    use spmd_utils,     only: masterproc
          18             :    use ppgrid,         only: pcols, pver, pverp
          19             :    use physconst,      only: latvap, latice, cpair
          20             :    use cam_abortutils, only: endrun
          21             :    use cam_logfile,    only: iulog
          22             :    use ref_pres,       only: top_lev => trop_cloud_top_lev
          23             :
          24             :    implicit none
          25             :
          26             : !-----------------------------------------------------------------------
          27             : ! PUBLIC: Make default data and interfaces private
          28             : !-----------------------------------------------------------------------
          29             :    private
          30             :    save
          31             :    public inimc, pcond          ! Public interfaces
          32             :    public :: cldwat_init
          33             :    integer, public::  ktop      ! Level above 10 hPa
          34             :
          35             :    real(r8),public ::  icritc               ! threshold for autoconversion of cold ice
          36             :    real(r8),public ::  icritw               ! threshold for autoconversion of warm ice
          37             : !!$   real(r8),public,parameter::  conke  = 1.e-6    ! tunable constant for evaporation of precip
          38             : !!$   real(r8),public,parameter::  conke  =  2.e-6    ! tunable constant for evaporation of precip
          39             :    real(r8),public ::  conke                ! tunable constant for evaporation of precip
          40             :    real(r8),public ::  r3lcrit              ! critical radius where liq conversion begins
          41             :
          42             : !-----------------------------------------------------------------------
          43             : ! PRIVATE: Everything else is private to this module
          44             : !-----------------------------------------------------------------------
          45             :    real(r8), private:: rhonot   ! air density at surface
          46             :    real(r8), private:: t0       ! Freezing temperature
          47             :    real(r8), private:: cldmin   ! assumed minimum cloud amount
          48             :    real(r8), private:: small    ! small number compared to unity
          49             :    real(r8), private:: c        ! constant for graupel like snow cm**(1-d)/s
          50             :    real(r8), private:: d        ! constant for graupel like snow
          51             :    real(r8), private:: esi      ! collection efficient for ice by snow
          52             :    real(r8), private:: esw      ! collection efficient for water by snow
          53             :    real(r8), private:: nos      ! particles snow / cm**4
          54             :    real(r8), private:: pi       ! Mathematical constant
          55             :    real(r8), private:: gravit   ! Gravitational acceleration at surface
          56             :    real(r8), private:: rh2o
          57             :    real(r8), private:: prhonos
          58             :    real(r8), private:: thrpd    ! numerical three added to d
          59             :    real(r8), private:: gam3pd   ! gamma function on (3+d)
          60             :    real(r8), private:: gam4pd   ! gamma function on (4+d)
          61             :    real(r8), private:: rhoi     ! ice density
          62             :    real(r8), private:: rhos     ! snow density
          63             :    real(r8), private:: rhow     ! water density
          64             :    real(r8), private:: mcon01   ! constants used in cloud microphysics
          65             :    real(r8), private:: mcon02   ! constants used in cloud microphysics
          66             :    real(r8), private:: mcon03   ! constants used in cloud microphysics
          67             :    real(r8), private:: mcon04   ! constants used in cloud microphysics
          68             :    real(r8), private:: mcon05   ! constants used in cloud microphysics
          69             :    real(r8), private:: mcon06   ! constants used in cloud microphysics
          70             :    real(r8), private:: mcon07   ! constants used in cloud microphysics
          71             :    real(r8), private:: mcon08   ! constants used in cloud microphysics
          72             :
          73             :
          74             : ! Parameters used in findmcnew
          75             :    real(r8) :: capnsi               ! sea ice cloud particles / cm3
          76             :    real(r8) :: capnc                ! cold and oceanic cloud particles / cm3
          77             :    real(r8) :: capnw                ! warm continental cloud particles / cm3
          78             :    real(r8) :: kconst               ! const for terminal velocity (stokes regime)
          79             :    real(r8) :: effc                 ! collection efficiency
          80             :    real(r8) :: alpha                ! ratio of 3rd moment radius to 2nd
          81             :    real(r8) :: capc                 ! constant for autoconversion
          82             :    real(r8) :: convfw               ! constant used for fall velocity calculation
          83             :    real(r8) :: cracw                ! constant used for rain accreting water
          84             :    real(r8) :: critpr               ! critical precip rate collection efficiency changes
          85             :    real(r8) :: ciautb               ! coefficient of autoconversion of ice (1/s)
          86             :    real(r8) :: psrhmin              ! condensation threshold in polar stratosphere
          87             :    logical  :: do_psrhmin
          88             :
          89             : #ifdef DEBUG
          90             :    integer, private,parameter ::  nlook = 1  ! Number of points to examine
          91             :    integer, private ::  ilook(nlook)         ! Longitude index to examine
          92             :    integer, private ::  latlook(nlook)       ! Latitude index to examine
          93             :    integer, private ::  lchnklook(nlook)     ! Chunk index to examine
          94             :    integer, private ::  icollook(nlook)      ! Column index to examine
          95             : #endif
          96             :
          97             : contains
          98             : !===============================================================================
          99        1536 : subroutine cldwat_init(icritw_in, icritc_in, conke_in, r3lcrit_in, psrhmin_in, do_psrhmin_in )
         100             :
         101             :     real(r8), intent(in) :: icritw_in    !   icritw  = threshold for autoconversion of warm ice
         102             :     real(r8), intent(in) :: icritc_in    !   icritc  = threshold for autoconversion of cold ice
         103             :     real(r8), intent(in) :: conke_in     !   conke   = tunable constant for evaporation of precip
         104             :     real(r8), intent(in) :: r3lcrit_in   !   r3lcrit = critical radius where liq conversion begins
         105             :     real(r8), intent(in) :: psrhmin_in   ! condensation threadhold in polar stratosphere
         106             :     logical,  intent(in) :: do_psrhmin_in
         107             :
         108        1536 :     icritw  = icritw_in
         109        1536 :     icritc  = icritc_in
         110        1536 :     conke   = conke_in
         111        1536 :     r3lcrit = r3lcrit_in
         112        1536 :     psrhmin = psrhmin_in
         113        1536 :     do_psrhmin = do_psrhmin_in
         114             :
         115        1536 :   end subroutine cldwat_init
         116             :
         117           0 : subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox)
         118             : !-----------------------------------------------------------------------
         119             : !
         120             : ! Purpose:
         121             : ! initialize constants for the prognostic condensate
         122             : !
         123             : ! Author: P. Rasch, April 1997
         124             : !
         125             : !-----------------------------------------------------------------------
         126             :    use pmgrid,       only: plev, plevp
         127             :    use ref_pres,     only: pref_mid
         128             :
         129             :    integer k
         130             :    real(r8), intent(in) :: tmeltx
         131             :    real(r8), intent(in) :: rhonotx
         132             :    real(r8), intent(in) :: gravitx
         133             :    real(r8), intent(in) :: rh2ox
         134             :
         135             : #ifdef UNICOSMP
         136             :    real(r8) signgam              ! variable required by cray gamma function
         137             :    external gamma
         138             : #endif
         139             :
         140           0 :    rhonot = rhonotx             ! air density at surface (gm/cm3)
         141           0 :    gravit = gravitx
         142           0 :    rh2o   = rh2ox
         143           0 :    rhos = .1_r8                 ! assumed snow density (gm/cm3)
         144           0 :    rhow = 1._r8                 ! water density
         145           0 :    rhoi = 1._r8                 ! ice density
         146           0 :    esi = 1.0_r8                 ! collection efficient for ice by snow
         147           0 :    esw = 0.1_r8                 ! collection efficient for water by snow
         148           0 :    t0 = tmeltx                  ! approximate freezing temp
         149           0 :    cldmin = 0.02_r8             ! assumed minimum cloud amount
         150           0 :    small = 1.e-22_r8            ! a small number compared to unity
         151           0 :    c = 152.93_r8                ! constant for graupel like snow cm**(1-d)/s
         152           0 :    d = 0.25_r8                  ! constant for graupel like snow
         153           0 :    nos = 3.e-2_r8               ! particles snow / cm**4
         154           0 :    pi = 4._r8*atan(1.0_r8)
         155           0 :    prhonos = pi*rhos*nos
         156           0 :    thrpd = 3._r8 + d
         157             :    if (d==0.25_r8) then
         158           0 :       gam3pd = 2.549256966718531_r8 ! only right for d = 0.25
         159           0 :       gam4pd = 8.285085141835282_r8
         160             :    else
         161             : #ifdef UNICOSMP
         162             :       call gamma(3._r8+d, signgam, gam3pd)
         163             :       gam3pd = sign(exp(gam3pd),signgam)
         164             :       call gamma(4._r8+d, signgam, gam4pd)
         165             :       gam4pd = sign(exp(gam4pd),signgam)
         166             :       write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd
         167             : #else
         168             :       call endrun(' can only use d ne 0.25 on a cray ')
         169             : #endif
         170             :    endif
         171           0 :    mcon01 = pi*nos*c*gam3pd/4._r8
         172           0 :    mcon02 = 1._r8/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._r8)))
         173           0 :    mcon03 = -(0.5_r8+d/4._r8)
         174           0 :    mcon04 = 4._r8/(4._r8+d)
         175           0 :    mcon05 = (3+d)/(4+d)
         176           0 :    mcon06 = (3+d)/4._r8
         177           0 :    mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06
         178           0 :    mcon08 = -0.5_r8/(4._r8+d)
         179             :
         180           0 :    if( masterproc ) write(iulog,*) 'cloud water initialization by inimc complete '
         181             :
         182             : ! Initialize parameters used by findmcnew
         183           0 :    capnw = 400._r8              ! warm continental cloud particles / cm3
         184           0 :    capnc = 150._r8              ! cold and oceanic cloud particles / cm3
         185             : !  capnsi = 40._r8              ! sea ice cloud particles density  / cm3
         186           0 :    capnsi = 75._r8              ! sea ice cloud particles density  / cm3
         187             :
         188           0 :    kconst = 1.18e6_r8           ! const for terminal velocity
         189             :
         190             : !  effc = 1._r8                 ! autoconv collection efficiency following boucher 96
         191             : !  effc = .55*0.05_r8           ! autoconv collection efficiency following baker 93
         192           0 :    effc = 0.55_r8               ! autoconv collection efficiency following tripoli and cotton
         193             : !  effc = 0._r8    ! turn off warm-cloud autoconv
         194           0 :    alpha = 1.1_r8**4
         195           0 :    capc = pi**(-.333_r8)*kconst*effc *(0.75_r8)**(1.333_r8)*alpha  ! constant for autoconversion
         196             :
         197             : ! critical precip rate at which we assume the collector drops can change the
         198             : ! drop size enough to enhance the auto-conversion process (mm/day)
         199           0 :    critpr = 0.5_r8
         200           0 :    convfw = 1.94_r8*2.13_r8*sqrt(rhow*1000._r8*9.81_r8*2.7e-4_r8)
         201             :
         202             : ! liquid microphysics
         203             : !  cracw = 6_r8                 ! beheng
         204           0 :    cracw = .884_r8*sqrt(9.81_r8/(rhow*1000._r8*2.7e-4_r8)) ! tripoli and cotton
         205             :
         206             : ! ice microphysics
         207           0 :    ciautb = 5.e-4_r8
         208             :
         209           0 :    if ( masterproc ) then
         210           0 :       write(iulog,*)'tuning parameters cldwat: icritw',icritw,'icritc',icritc,'conke',conke,'r3lcrit',r3lcrit
         211           0 :       write(iulog,*)'tuning parameters cldwat: capnw',capnw,'capnc',capnc,'capnsi',capnsi,'kconst',kconst
         212           0 :       write(iulog,*)'tuning parameters cldwat: effc',effc,'alpha',alpha,'capc',capc
         213           0 :       write(iulog,*)'tuning parameters cldwat: critpr',critpr,'convfw',convfw,'cracw',cracw,'ciautb',ciautb
         214             :    endif
         215             :
         216           0 : end subroutine inimc
         217             :
         218           0 : subroutine pcond (lchnk   ,ncol    ,troplev ,dlat    , &
         219             :                   tn      ,ttend   ,qn      ,qtend   ,omega   , &
         220             :                   cwat    ,p       ,pdel    ,cldn    ,fice    , fsnow, &
         221             :                   cme     ,prodprec,prodsnow,evapprec,evapsnow,evapheat, prfzheat, &
         222             :                   meltheat,precip  ,snowab  ,deltat  ,fwaut   , &
         223             :                   fsaut   ,fracw   ,fsacw   ,fsaci   ,lctend  , &
         224             :                   rhdfda  ,rhu00   ,landm   ,seaicef ,zi      ,ice2pr, liq2pr, &
         225             :                   liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio)
         226             : !-----------------------------------------------------------------------
         227             : !
         228             : ! Purpose:
         229             : ! The public interface to the cloud water parameterization
         230             : ! returns tendencies to water vapor, temperature and cloud water variables
         231             : !
         232             : ! For basic method
         233             : !  See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
         234             : !  model climate using diagnosed and
         235             : !  predicted condensate parameterizations, 1998, J. Clim., 11,
         236             : !  pp1587---1614.
         237             : !
         238             : ! For important modifications to improve the method of determining
         239             : ! condensation/evaporation see Zhang et al (2001, in preparation)
         240             : !
         241             : ! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson
         242             : !          B. A. Boville (latent heat of fusion)
         243             : !-----------------------------------------------------------------------
         244           0 :    use wv_saturation, only: qsat, estblf, svp_to_qsat, findsp_vc
         245             :    use physconst, only: epsilo
         246             : !
         247             : !---------------------------------------------------------------------
         248             : !
         249             : ! Input Arguments
         250             : !
         251             :    integer,  intent(in) :: lchnk                ! chunk identifier
         252             :    integer,  intent(in) :: ncol                 ! number of atmospheric columns
         253             :    integer,  intent(in) :: troplev(pcols)       ! tropopause level
         254             :    real(r8), intent(in) :: dlat(pcols)          ! latitudes in degrees
         255             :    real(r8), intent(in) :: fice(pcols,pver)     ! fraction of cwat that is ice
         256             :    real(r8), intent(in) :: fsnow(pcols,pver)    ! fraction of rain that freezes to snow
         257             :    real(r8), intent(in) :: cldn(pcols,pver)     ! new value of cloud fraction    (fraction)
         258             :    real(r8), intent(in) :: cwat(pcols,pver)     ! cloud water (kg/kg)
         259             :    real(r8), intent(in) :: omega(pcols,pver)    ! vert pressure vel (Pa/s)
         260             :    real(r8), intent(in) :: p(pcols,pver)        ! pressure          (K)
         261             :    real(r8), intent(in) :: pdel(pcols,pver)     ! pressure thickness (Pa)
         262             :    real(r8), intent(in) :: qn(pcols,pver)       ! new water vapor    (kg/kg)
         263             :    real(r8), intent(in) :: qtend(pcols,pver)    ! mixing ratio tend  (kg/kg/s)
         264             :    real(r8), intent(in) :: tn(pcols,pver)       ! new temperature    (K)
         265             :    real(r8), intent(in) :: ttend(pcols,pver)    ! temp tendencies    (K/s)
         266             :    real(r8), intent(in) :: deltat               ! time step to advance solution over
         267             :    real(r8), intent(in) :: lctend(pcols,pver)   ! cloud liquid water tendencies   ====wlin
         268             :    real(r8), intent(in) :: rhdfda(pcols,pver)   ! dG(a)/da, rh=G(a), when rh>u00  ====wlin
         269             :    real(r8), intent(in) :: rhu00 (pcols,pver)   ! Rhlim for cloud                 ====wlin
         270             :    real(r8), intent(in) :: landm(pcols)         ! Land fraction ramped over water (fraction)
         271             :    real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction  (fraction)
         272             :    real(r8), intent(in) :: zi(pcols,pverp)      ! layer interfaces (m)
         273             :    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
         274             : !
         275             : ! Output Arguments
         276             : !
         277             :    real(r8), intent(out) :: cme     (pcols,pver) ! rate of cond-evap of condensate (1/s)
         278             :    real(r8), intent(out) :: prodprec(pcols,pver) ! rate of conversion of condensate to precip (1/s)
         279             :    real(r8), intent(out) :: evapprec(pcols,pver) ! rate of evaporation of falling precip (1/s)
         280             :    real(r8), intent(out) :: evapsnow(pcols,pver) ! rate of evaporation of falling snow (1/s)
         281             :    real(r8), intent(out) :: evapheat(pcols,pver) ! heating rate due to evaporation of precip (W/kg)
         282             :    real(r8), intent(out) :: prfzheat(pcols,pver) ! heating rate due to freezing of precip (W/kg)
         283             :    real(r8), intent(out) :: meltheat(pcols,pver) ! heating rate due to snow melt (W/kg)
         284             :    real(r8), intent(out) :: precip(pcols)        ! rate of precipitation (kg / (m**2 * s))
         285             :    real(r8), intent(out) :: snowab(pcols)        ! rate of snow (kg / (m**2 * s))
         286             :    real(r8), intent(out) :: ice2pr(pcols,pver)   ! rate of conversion of ice to precip
         287             :    real(r8), intent(out) :: liq2pr(pcols,pver)   ! rate of conversion of liquid to precip
         288             :    real(r8), intent(out) :: liq2snow(pcols,pver) ! rate of conversion of liquid to snow
         289             :    real(r8), intent(out) :: rkflxprc(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_rain+snow at interfaces (kg m^-2 s^-1)
         290             :    real(r8), intent(out) :: rkflxsnw(pcols,pverp)   ! grid-box mean RK flux_large_scale_cloud_snow at interfaces (kg m^-2 s^-1)
         291             : ! intent(out)s here for pcond to pass to stratiform.F90 to be addflded/outflded
         292             :    real(r8), intent(out) :: pracwo(pcols,pver)      ! accretion of cloud water by rain (1/s)
         293             :    real(r8), intent(out) :: psacwo(pcols,pver)      ! accretion of cloud water by snow (1/s)
         294             :    real(r8), intent(out) :: psacio(pcols,pver)      ! accretion of cloud ice by snow (1/s)
         295             :
         296             :    real(r8) nice2pr     ! rate of conversion of ice to snow
         297             :    real(r8) nliq2pr     ! rate of conversion of liquid to precip
         298             :    real(r8) nliq2snow   ! rate of conversion of liquid to snow
         299             :    real(r8) prodsnow(pcols,pver) ! rate of production of snow
         300             :
         301             : !
         302             : ! Local workspace
         303             : !
         304             :    real(r8) :: precab(pcols)        ! rate of precipitation (kg / (m**2 * s))
         305             :    integer i                 ! work variable
         306             :    integer iter              ! #iterations for precipitation calculation
         307             :    integer k                 ! work variable
         308             :    integer l                 ! work variable
         309             :
         310             :    real(r8) cldm(pcols)          ! mean cloud fraction over the time step
         311             :    real(r8) cldmax(pcols)        ! max cloud fraction above
         312             :    real(r8) coef(pcols)          ! conversion time scale for condensate to rain
         313             :    real(r8) cwm(pcols)           ! cwat mixing ratio at midpoint of time step
         314             :    real(r8) cwn(pcols)           ! cwat mixing ratio at end
         315             :    real(r8) denom                ! work variable
         316             :    real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
         317             :    real(r8) es(pcols)            ! sat. vapor pressure
         318             :    real(r8) fracw(pcols,pver)    ! relative importance of collection of liquid by rain
         319             :    real(r8) fsaci(pcols,pver)    ! relative importance of collection of ice by snow
         320             :    real(r8) fsacw(pcols,pver)    ! relative importance of collection of liquid by snow
         321             :    real(r8) fsaut(pcols,pver)    ! relative importance of ice auto conversion
         322             :    real(r8) fwaut(pcols,pver)    ! relative importance of warm cloud autoconversion
         323             :    real(r8) gamma(pcols)         ! d qs / dT
         324             :    real(r8) icwc(pcols)          ! in-cloud water content (kg/kg)
         325             :    real(r8) mincld               ! a small cloud fraction to avoid / zero
         326             :    real(r8),parameter ::omsm=0.99999_r8                 ! a number just less than unity (for rounding)
         327             :    real(r8) prprov(pcols)        ! provisional value of precip at btm of layer
         328             :    real(r8) prtmp                ! work variable
         329             :    real(r8) q(pcols,pver)        ! mixing ratio before time step ignoring condensate
         330             :    real(r8) qs(pcols)            ! spec. hum. of water vapor
         331             :    real(r8) qsn, esn             ! work variable
         332             :    real(r8) qsp(pcols,pver)      ! sat pt mixing ratio
         333             :    real(r8) qtl(pcols)           ! tendency which would saturate the grid box in deltat
         334             :    real(r8) qtmp, ttmp           ! work variable
         335             :    real(r8) relhum1(pcols)        ! relative humidity
         336             :    real(r8) relhum(pcols)        ! relative humidity
         337             : !!$   real(r8) tc                   ! crit temp of transition to ice
         338             :    real(r8) t(pcols,pver)        ! temp before time step ignoring condensate
         339             :    real(r8) tsp(pcols,pver)      ! sat pt temperature
         340             :    real(r8) pol                  ! work variable
         341             :    real(r8) cdt                  ! work variable
         342             :    real(r8) wtthick              ! work variable
         343             :
         344             : ! Extra local work space for cloud scheme modification
         345             :
         346             :    real(r8) cpohl                !cpair/Latvap
         347             :    real(r8) hlocp                !Latvap/cpair
         348             :    real(r8) dto2                 !0.5*deltat (delta=2.0*dt)
         349             :    real(r8) calpha(pcols)        !alpha of new C - E scheme formulation
         350             :    real(r8) cbeta (pcols)        !beta  of new C - E scheme formulation
         351             :    real(r8) cbetah(pcols)        !beta_hat at saturation portion
         352             :    real(r8) cgamma(pcols)        !gamma of new C - E scheme formulation
         353             :    real(r8) cgamah(pcols)        !gamma_hat at saturation portion
         354             :    real(r8) rcgama(pcols)        !gamma/gamma_hat
         355             :    real(r8) csigma(pcols)        !sigma of new C - E scheme formulation
         356             :    real(r8) cmec1 (pcols)        !c1    of new C - E scheme formulation
         357             :    real(r8) cmec2 (pcols)        !c2    of new C - E scheme formulation
         358             :    real(r8) cmec3 (pcols)        !c3    of new C - E scheme formulation
         359             :    real(r8) cmec4 (pcols)        !c4    of new C - E scheme formulation
         360             :    real(r8) cmeres(pcols)        !residual cond of over-sat after cme and evapprec
         361             :    real(r8) ctmp                 !a scalar representation of cmeres
         362             :    real(r8) clrh2o               ! Ratio of latvap to water vapor gas const
         363             :    real(r8) ice(pcols,pver)    ! ice mixing ratio
         364             :    real(r8) liq(pcols,pver)    ! liquid mixing ratio
         365             :    real(r8) rcwn(pcols,2,pver), rliq(pcols,2,pver), rice(pcols,2,pver)
         366             :    real(r8) cwnsave(pcols,2,pver), cmesave(pcols,2,pver)
         367             :    real(r8) prodprecsave(pcols,2,pver)
         368             :    logical error_found
         369             :
         370             :    real(r8) :: rhu_adj(pcols,pver)   ! adjusted rhlim for dehydration
         371             : !
         372             : !------------------------------------------------------------
         373             : !
         374           0 :    clrh2o = latvap/rh2o   ! Ratio of latvap to water vapor gas const
         375             : #ifdef PERGRO
         376             :    mincld = 1.e-4_r8
         377             :    iter = 1   ! number of times to iterate the precipitation calculation
         378             : #else
         379           0 :    mincld = 1.e-4_r8
         380           0 :    iter = 2
         381             : #endif
         382             : !   omsm = 0.99999
         383           0 :    cpohl = cpair/latvap
         384           0 :    hlocp = latvap/cpair
         385           0 :    dto2=0.5_r8*deltat
         386             : !
         387             : ! Constant for computing rate of evaporation of precipitation:
         388             : !
         389             : !!$   conke = 1.e-5
         390             : !!$   conke = 1.e-6
         391             : !
         392             : ! initialize a few single level fields
         393             : !
         394           0 :    do i = 1,ncol
         395           0 :       precip(i) = 0.0_r8
         396           0 :       precab(i) = 0.0_r8
         397           0 :       snowab(i) = 0.0_r8
         398           0 :       cldmax(i) = 0.0_r8
         399             :    end do
         400             : !
         401             : ! initialize multi-level fields
         402             : !
         403           0 :    do k = 1,pver
         404           0 :       do i = 1,ncol
         405           0 :          q(i,k) = qn(i,k)
         406           0 :          t(i,k) = tn(i,k)
         407             : !         q(i,k)=qn(i,k)-qtend(i,k)*deltat
         408             : !         t(i,k)=tn(i,k)-ttend(i,k)*deltat
         409             :     end do
         410             :    end do
         411           0 :    cme     (:ncol,:) = 0._r8
         412           0 :    evapprec(:ncol,:) = 0._r8
         413           0 :    prodprec(:ncol,:) = 0._r8
         414           0 :    evapsnow(:ncol,:) = 0._r8
         415           0 :    prodsnow(:ncol,:) = 0._r8
         416           0 :    evapheat(:ncol,:) = 0._r8
         417           0 :    meltheat(:ncol,:) = 0._r8
         418           0 :    prfzheat(:ncol,:) = 0._r8
         419           0 :    ice2pr(:ncol,:)   = 0._r8
         420           0 :    liq2pr(:ncol,:)   = 0._r8
         421           0 :    liq2snow(:ncol,:) = 0._r8
         422           0 :    fwaut(:ncol,:) = 0._r8
         423           0 :    fsaut(:ncol,:) = 0._r8
         424           0 :    fracw(:ncol,:) = 0._r8
         425           0 :    fsacw(:ncol,:) = 0._r8
         426           0 :    fsaci(:ncol,:) = 0._r8
         427           0 :    rkflxprc(:ncol,:) = 0._r8
         428           0 :    rkflxsnw(:ncol,:) = 0._r8
         429             :
         430           0 :    pracwo(:ncol,:) = 0._r8
         431           0 :    psacwo(:ncol,:) = 0._r8
         432           0 :    psacio(:ncol,:) = 0._r8
         433             :
         434             : !
         435             : ! find the wet bulb temp and saturation value
         436             : ! for the provisional t and q without condensation
         437             : !
         438           0 :    do 800 k = top_lev,pver
         439             :
         440             :       ! "True" means that ice will be taken into account.
         441             :       call findsp_vc(qn(:ncol,k), tn(:ncol,k), p(:ncol,k), .true., &
         442           0 :            tsp(:ncol,k), qsp(:ncol,k))
         443             :
         444           0 :       call qsat(t(1:ncol,k), p(1:ncol,k), es(1:ncol), qs(1:ncol), ncol, gam=gamma(1:ncol))
         445             :
         446           0 :       do i = 1,ncol
         447             : !
         448           0 :          relhum(i) = q(i,k)/qs(i)
         449             : !
         450           0 :          cldm(i) = max(cldn(i,k),mincld)
         451             : !
         452             : ! the max cloud fraction above this level
         453             : !
         454           0 :          cldmax(i) = max(cldmax(i), cldm(i))
         455             :
         456             : ! define the coefficients for C - E calculation
         457             :
         458           0 :          calpha(i) = 1.0_r8/qs(i)
         459           0 :          cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl
         460           0 :          cbetah(i) = 1.0_r8/qs(i)*gamma(i)*cpohl
         461           0 :          cgamma(i) = calpha(i)+latvap*cbeta(i)/cpair
         462           0 :          cgamah(i) = calpha(i)+latvap*cbetah(i)/cpair
         463           0 :          rcgama(i) = cgamma(i)/cgamah(i)
         464             :
         465           0 :          if(cldm(i) > mincld) then
         466           0 :             icwc(i) = max(0._r8,cwat(i,k)/cldm(i))
         467             :          else
         468           0 :             icwc(i) = 0.0_r8
         469             :          endif
         470             : !PJR the above logic give zero icwc with nonzero cwat, dont like it!
         471             : !PJR generates problems with csigma
         472             : !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
         473             : !         icwc(i) = max(1.e-8_r8,cwat(i,k)/cldm(i))
         474             :
         475             : !
         476             : ! initial guess of evaporation, will be updated within iteration
         477             : !
         478             :          evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
         479           0 :                         *(1._r8 - min(relhum(i),1._r8))
         480             :
         481             : !
         482             : ! zero cmeres before iteration for each level
         483             : !
         484           0 :          cmeres(i)=0.0_r8
         485             :
         486             :       end do
         487           0 :       do i = 1,ncol
         488             : !
         489             : ! fractions of ice at this level
         490             : !
         491             : !!$         tc = t(i,k) - t0
         492             : !!$         fice(i,k) = max(0._r8,min(-tc*0.05,1.0_r8))
         493             : !
         494             : ! calculate the cooling due to a phase change of the rainwater
         495             : ! from above
         496             : !
         497           0 :          if (t(i,k) >= t0) then
         498           0 :             meltheat(i,k) =  -latice * snowab(i) * gravit/pdel(i,k)
         499           0 :             snowab(i) = 0._r8
         500             :          else
         501           0 :             meltheat(i,k) = 0._r8
         502             :          endif
         503             :       end do
         504             :
         505             : !
         506             : ! calculate cme and formation of precip.
         507             : !
         508             : ! The cloud microphysics is highly nonlinear and coupled with cme
         509             : ! Both rain processes and cme are calculated iteratively.
         510             : !
         511           0 :       do 100 l = 1,iter
         512             :
         513           0 :         do i = 1,ncol
         514             :
         515             : !
         516             : ! calculation of cme has 4 scenarios
         517             : ! ==================================
         518             : !
         519           0 :            call relhum_min_adj( ncol, troplev, dlat, rhu00, rhu_adj )
         520             :
         521           0 :            if(relhum(i) > rhu_adj(i,k)) then
         522             :
         523             :            ! 1. whole grid saturation
         524             :            ! ========================
         525           0 :               if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then
         526           0 :                  cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)
         527             :
         528             :            ! 2. fractional saturation
         529             :            ! ========================
         530             :               else
         531           0 :                  if (rhdfda(i,k) .eq. 0._r8 .and. icwc(i).eq.0._r8) then
         532           0 :                     write (iulog,*) ' cldwat.F90:  empty rh cloud ', i, k, lchnk
         533           0 :                     write (iulog,*) ' relhum, iter ', relhum(i), l, rhu_adj(i,k), cldm(i), cldn(i,k)
         534           0 :                     call endrun ()
         535             :                  endif
         536           0 :                   csigma(i) = 1.0_r8/(rhdfda(i,k)+cgamma(i)*icwc(i))
         537           0 :                   cmec1(i) = (1.0_r8-cldm(i))*csigma(i)*rhdfda(i,k)
         538             :                   cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_r8-rcgama(i)*cldm(i))*   &
         539           0 :                              csigma(i)*calpha(i)*icwc(i)
         540             :                   cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) +  &
         541           0 :                            (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)
         542           0 :                   cmec4(i) = csigma(i)*cgamma(i)*icwc(i)
         543             :
         544             :                   ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er
         545             :
         546             :                   cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k)  &
         547           0 :                              -cmec3(i)*ttend(i,k) + cmec4(i)*evapprec(i,k)
         548             :                endif
         549             :
         550             :            ! 3. when rh < rhu00, evaporate existing cloud water
         551             :            ! ==================================================
         552           0 :            else if(cwat(i,k) > 0.0_r8)then
         553             :               ! liquid water should be evaporated but not to exceed
         554             :               ! saturation point. if qn > qsp, not to evaporate cwat
         555           0 :               cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat
         556             :
         557             :            ! 4. no condensation nor evaporation
         558             :            ! ==================================
         559             :            else
         560           0 :               cme(i,k)=0.0_r8
         561             :            endif
         562             :
         563             :
         564             :         end do    !end loop for cme update
         565             :
         566             : ! Because of the finite time step,
         567             : ! place a bound here not to exceed wet bulb point
         568             : ! and not to evaporate more than available water
         569             : !
         570           0 :          do i = 1, ncol
         571           0 :             qtmp = qn(i,k) - cme(i,k)*deltat
         572             :
         573             : ! possibilities to have qtmp > qsp
         574             : !
         575             : !   1. if qn > qs(tn), it condenses;
         576             : !      if after applying cme,  qtmp > qsp,  more condensation is applied.
         577             : !
         578             : !   2. if qn < qs, evaporation should not exceed qsp,
         579             :
         580           0 :             if(qtmp > qsp(i,k)) then
         581           0 :               cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat
         582             :             endif
         583             :
         584             : !
         585             : ! if net evaporation, it should not exceed available cwat
         586             : !
         587           0 :             if(cme(i,k) < -cwat(i,k)/deltat)  &
         588           0 :                cme(i,k) = -cwat(i,k)/deltat
         589             : !
         590             : ! addition of residual condensation from previous step of iteration
         591             : !
         592           0 :             cme(i,k) = cme(i,k) + cmeres(i)
         593             :
         594             :          end do
         595             :
         596             :          !      limit cme for roundoff errors
         597           0 :          do i = 1, ncol
         598           0 :             cme(i,k) = cme(i,k)*omsm
         599             :          end do
         600             :
         601           0 :          do i = 1,ncol
         602             : !
         603             : ! as a safe limit, condensation should not reduce grid mean rh below rhu00
         604             : !
         605           0 :            if(cme(i,k) > 0.0_r8 .and. relhum(i) > rhu_adj(i,k) )  &
         606           0 :               cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu_adj(i,k))/deltat)
         607             : !
         608             : ! initial guess for cwm (mean cloud water over time step) if 1st iteration
         609             : !
         610           0 :            if(l < 2) then
         611           0 :              cwm(i) = max(cwat(i,k)+cme(i,k)*dto2,  0._r8)
         612             :            endif
         613             :
         614             :          enddo
         615             :
         616             : ! provisional precipitation falling through model layer
         617           0 :          do i = 1,ncol
         618             : !!$            prprov(i) =  precab(i) + prodprec(i,k)*pdel(i,k)/gravit
         619             : ! rain produced in this layer not too effective in collection process
         620           0 :             wtthick = max(0._r8,min(0.5_r8,((zi(i,k)-zi(i,k+1))/1000._r8)**2))
         621           0 :             prprov(i) =  precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit
         622             :          end do
         623             :
         624             : ! calculate conversion of condensate to precipitation by cloud microphysics
         625             :          call findmcnew (lchnk   ,ncol    , &
         626             :                          k       ,prprov  ,snowab,  t       ,p        , &
         627           0 :                          cwm     ,cldm    ,cldmax  ,fice(1,k),coef    , &
         628             :                          fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &
         629           0 :                          landm, seaicef, snowh, pracwo(1,k), psacwo(1,k), psacio(1,k))
         630             :
         631             : !
         632             : ! calculate the precip rate
         633             : !
         634           0 :          error_found = .false.
         635           0 :          do i = 1,ncol
         636           0 :             if (cldm(i) > 0) then
         637             : !
         638             : ! first predict the cloud water
         639             : !
         640           0 :                cdt = coef(i)*deltat
         641           0 :                if (cdt > 0.01_r8) then
         642           0 :                   pol = cme(i,k)/coef(i) ! production over loss
         643           0 :                   cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)
         644             :                else
         645           0 :                   cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))
         646             :                endif
         647             : !
         648             : ! now back out the tendency of net rain production
         649             : !
         650           0 :                prodprec(i,k) =  max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)
         651             :             else
         652           0 :                prodprec(i,k) = 0.0_r8
         653           0 :                cwn(i) = 0._r8
         654             :             endif
         655             :
         656             :             ! provisional calculation of conversion terms
         657           0 :             ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k))
         658           0 :             liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k))
         659             : !old        liq2snow(i,k) = prodprec(i,k)*fsacw(i,k)
         660             :
         661             : !           revision suggested by Jim McCaa
         662             : !           it controls the amount of snow hitting the sfc
         663             : !           by forcing a lot of conversion of cloud liquid to snow phase
         664             : !           it might be better done later by an explicit representation of
         665             : !           rain accreting ice (and freezing), or by an explicit freezing of raindrops
         666           0 :             liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k))
         667             :
         668             :             ! bounds
         669           0 :             nice2pr = min(ice2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*fice(i,k)/deltat)
         670           0 :             nliq2pr = min(liq2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*(1._r8-fice(i,k))/deltat)
         671             : !            write(iulog,*) ' prodprec ', i, k, prodprec(i,k)
         672             : !            write(iulog,*) ' nliq2pr, nice2pr ', nliq2pr, nice2pr
         673           0 :             if (liq2pr(i,k).ne.0._r8) then
         674           0 :                nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k)   ! correction
         675             :             else
         676             :                nliq2snow = liq2snow(i,k)
         677             :             endif
         678             :
         679             : !           avoid roundoff problems generating negatives
         680           0 :             nliq2snow = nliq2snow*omsm
         681           0 :             nliq2pr = nliq2pr*omsm
         682           0 :             nice2pr = nice2pr*omsm
         683             :
         684             : !           final estimates of conversion to precip and snow
         685           0 :             prodprec(i,k) = (nliq2pr + nice2pr)
         686           0 :             prodsnow(i,k) = (nice2pr + nliq2snow)
         687             :
         688           0 :             rcwn(i,l,k) =  cwat(i,k) + (cme(i,k)-   prodprec(i,k))*deltat
         689           0 :             rliq(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)*(1._r8-fice(i,k)) - nliq2pr * deltat
         690           0 :             rice(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)* fice(i,k)      -    nice2pr                     *deltat
         691             :
         692             : !           Save for sanity check later...
         693             : !           Putting sanity checks inside loops 100 and 800 screws up the
         694             : !           IBM compiler for reasons as yet unknown.  TBH
         695           0 :             cwnsave(i,l,k)      = cwn(i)
         696           0 :             cmesave(i,l,k)      = cme(i,k)
         697           0 :             prodprecsave(i,l,k) = prodprec(i,k)
         698             : !           End of save for sanity check later...
         699             :
         700             : !           final version of condensate to precip terms
         701           0 :             liq2pr(i,k) = nliq2pr
         702           0 :             liq2snow(i,k) = nliq2snow
         703           0 :             ice2pr(i,k) = nice2pr
         704             :
         705           0 :             cwn(i) = rcwn(i,l,k)
         706             : !
         707             : ! update any remaining  provisional values
         708             : !
         709           0 :             cwm(i) = (cwn(i) + cwat(i,k))*0.5_r8
         710             : !
         711             : ! update in cloud water
         712             : !
         713           0 :             if(cldm(i) > mincld) then
         714           0 :                icwc(i) = cwm(i)/cldm(i)
         715             :             else
         716           0 :                icwc(i) = 0.0_r8
         717             :             endif
         718             : !PJR the above logic give zero icwc with nonzero cwat, dont like it!
         719             : !PJR generates problems with csigma
         720             : !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
         721             : !         icwc(i) = max(1.e-8_r8,cwm(i)/cldm(i))
         722             :
         723             :          end do              ! end of do i = 1,ncol
         724             :
         725             : !
         726             : ! calculate provisional value of cloud water for
         727             : ! evaporation of precipitate (evapprec) calculation
         728             : !
         729           0 :       do i = 1,ncol
         730           0 :          qtmp = qn(i,k) - cme(i,k)*deltat
         731             :          ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k)       &
         732           0 :               + (latvap + latice*fice(i,k)) * cme(i,k) )
         733           0 :          esn = estblf(ttmp)
         734           0 :          qsn = svp_to_qsat(esn, p(i,k))
         735           0 :          qtl(i) = max((qsn - qtmp)/deltat,0._r8)
         736           0 :          relhum1(i) = qtmp/qsn
         737             :       end do
         738             : !
         739           0 :       do i = 1,ncol
         740             : #ifdef PERGRO
         741             :          evapprec(i,k) = conke*(1._r8 - max(cldm(i),mincld))* &
         742             :                          sqrt(precab(i))*(1._r8 - min(relhum1(i),1._r8))
         743             : #else
         744           0 :          evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
         745           0 :                          *(1._r8 - min(relhum1(i),1._r8))
         746             : #endif
         747             : !
         748             : ! limit the evaporation to the amount which is entering the box
         749             : ! or saturates the box
         750             : !
         751           0 :          prtmp = precab(i)*gravit/pdel(i,k)
         752           0 :          evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm
         753             : #ifdef PERGRO
         754             : !           zeroing needed for pert growth
         755             :          evapprec(i,k) = 0._r8
         756             : #endif
         757             : !
         758             : ! Partition evaporation of precipitate between rain and snow using
         759             : ! the fraction of snow falling into the box. Determine the heating
         760             : ! due to evaporation. Note that evaporation is positive (loss of precip,
         761             : ! gain of vapor) and that heating is negative.
         762           0 :          if (evapprec(i,k) > 0._r8) then
         763           0 :             evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i)
         764           0 :             evapheat(i,k) = -latvap * evapprec(i,k) - latice * evapsnow(i,k)
         765             :          else
         766           0 :             evapsnow(i,k) = 0._r8
         767           0 :             evapheat(i,k) = 0._r8
         768             :          end if
         769             : ! Account for the latent heat of fusion for liquid drops collected by falling snow
         770           0 :          prfzheat(i,k) = latice * liq2snow(i,k)
         771             :       end do
         772             :
         773             : ! now remove the residual of any over-saturation. Normally,
         774             : ! the oversaturated water vapor should have been removed by
         775             : ! cme formulation plus constraints by wet bulb tsp/qsp
         776             : ! as computed above. However, because of non-linearity,
         777             : ! addition of (cme-evapprec) to update t and q may still cause
         778             : ! a very small amount of over saturation. It is called a
         779             : ! residual of over-saturation because theoretically, cme
         780             : ! should have taken care of all of large scale condensation.
         781             : !
         782             :
         783           0 :        do i = 1,ncol
         784           0 :           qtmp = qn(i,k)-(cme(i,k)-evapprec(i,k))*deltat
         785             :           ttmp = tn(i,k) + deltat/cpair * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k)      &
         786           0 :               + (latvap + latice*fice(i,k)) * cme(i,k) )
         787             :
         788           0 :           call qsat(ttmp, p(i,k), esn, qsn, dqsdt=dqsdt)
         789             :
         790           0 :           if( qtmp > qsn ) then
         791             :              !
         792             :              !now extra condensation to bring air to just saturation
         793             :              !
         794           0 :              ctmp = (qtmp-qsn)/(1._r8+hlocp*dqsdt)/deltat
         795           0 :              cme(i,k) = cme(i,k)+ctmp
         796             : !
         797             : ! save residual on cmeres to addtion to cme on entering next iteration
         798             : ! cme exit here contain the residual but overrided if back to iteration
         799             : !
         800           0 :              cmeres(i) = ctmp
         801             :           else
         802           0 :              cmeres(i) = 0.0_r8
         803             :           endif
         804             :        end do
         805             :
         806           0 :  100 continue              ! end of do l = 1,iter
         807             :
         808             : !
         809             : ! precipitation
         810             : !
         811           0 :       do i = 1,ncol
         812           0 :          precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         813           0 :          precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
         814           0 :          if(precab(i).lt.0._r8) precab(i)=0._r8
         815             : !         snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k))
         816           0 :          snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k))
         817             :
         818             :          ! If temperature above freezing, all precip is rain flux.  if temperature below freezing, all precip is snow flux.
         819           0 :          rkflxprc(i,k+1) = precab(i)   !! making this consistent with other precip fluxes.  prc = rain + snow
         820             :          !!rkflxprc(i,k+1) = precab(i) - snowab(i)
         821           0 :          rkflxsnw(i,k+1) = snowab(i)
         822             :
         823             : !!$         if ((precab(i)) < 1.e-10) then
         824             : !!$            precab(i) = 0.
         825             : !!$            snowab(i) = 0.
         826             : !!$         endif
         827             :       end do
         828           0 :  800 continue                ! level loop (k=1,pver)
         829             :
         830             : ! begin sanity checks
         831           0 :    error_found = .false.
         832           0 :    do k = top_lev,pver
         833           0 :       do l = 1,iter
         834           0 :          do i = 1,ncol
         835           0 :             if (abs(rcwn(i,l,k)).lt.1.e-300_r8) rcwn(i,l,k) = 0._r8
         836           0 :             if (abs(rliq(i,l,k)).lt.1.e-300_r8) rliq(i,l,k) = 0._r8
         837           0 :             if (abs(rice(i,l,k)).lt.1.e-300_r8) rice(i,l,k) = 0._r8
         838           0 :             if (rcwn(i,l,k).lt.0._r8) error_found = .true.
         839           0 :             if (rliq(i,l,k).lt.0._r8) error_found = .true.
         840           0 :             if (rice(i,l,k).lt.0._r8) error_found = .true.
         841             :          enddo
         842             :       enddo
         843             :    enddo
         844           0 :    if (error_found) then
         845           0 :       do k = top_lev,pver
         846           0 :          do l = 1,iter
         847           0 :             do i = 1,ncol
         848           0 :                if (rcwn(i,l,k).lt.0._r8) then
         849           0 :                   write(iulog,*) ' prob with neg rcwn1 ', rcwn(i,l,k),  &
         850           0 :                      cwnsave(i,l,k)
         851           0 :                   write(iulog,*) ' cwat, cme*deltat, prodprec*deltat ', &
         852           0 :                      cwat(i,k), cmesave(i,l,k)*deltat,               &
         853           0 :                      prodprecsave(i,l,k)*deltat,                     &
         854           0 :                      (cmesave(i,l,k)-prodprecsave(i,l,k))*deltat
         855           0 :                   call endrun('PCOND')
         856             :                endif
         857           0 :                if (rliq(i,l,k).lt.0._r8) then
         858           0 :                   write(iulog,*) ' prob with neg rliq1 ', rliq(i,l,k)
         859           0 :                   call endrun('PCOND')
         860             :                endif
         861           0 :                if (rice(i,l,k).lt.0._r8) then
         862           0 :                   write(iulog,*) ' prob with neg rice ', rice(i,l,k)
         863           0 :                   call endrun('PCOND')
         864             :                endif
         865             :             enddo
         866             :          enddo
         867             :       enddo
         868             :    end if
         869             : ! end sanity checks
         870             :
         871           0 :    return
         872             : end subroutine pcond
         873             :
         874             : !##############################################################################
         875             :
         876           0 : subroutine findmcnew (lchnk   ,ncol    , &
         877             :                       k       ,precab  ,snowab,  t       ,p       , &
         878             :                       cwm     ,cldm    ,cldmax  ,fice    ,coef    , &
         879             :                       fwaut   ,fsaut   ,fracw   ,fsacw   ,fsaci   , &
         880             :                       landm   ,seaicef ,snowh   ,pracwo  ,psacwo  ,psacio )
         881             :
         882             : !-----------------------------------------------------------------------
         883             : !
         884             : ! Purpose:
         885             : ! calculate the conversion of condensate to precipitate
         886             : !
         887             : ! Method:
         888             : ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
         889             : !  model climate using diagnosed and
         890             : !  predicted condensate parameterizations, 1998, J. Clim., 11,
         891             : !  pp1587---1614.
         892             : !
         893             : ! Author: P. Rasch
         894             : !
         895             : !-----------------------------------------------------------------------
         896             :    use phys_grid, only: get_rlat_all_p
         897             : !
         898             : ! input args
         899             : !
         900             :    integer, intent(in) :: lchnk                 ! chunk identifier
         901             :    integer, intent(in) :: ncol                  ! number of atmospheric columns
         902             :    integer, intent(in) :: k                     ! level index
         903             :
         904             :    real(r8), intent(in) :: precab(pcols)        ! rate of precipitation from above (kg / (m**2 * s))
         905             :    real(r8), intent(in) :: t(pcols,pver)        ! temperature       (K)
         906             :    real(r8), intent(in) :: p(pcols,pver)        ! pressure          (Pa)
         907             :    real(r8), intent(in) :: cldm(pcols)          ! cloud fraction
         908             :    real(r8), intent(in) :: cldmax(pcols)        ! max cloud fraction above this level
         909             :    real(r8), intent(in) :: cwm(pcols)           ! condensate mixing ratio (kg/kg)
         910             :    real(r8), intent(in) :: fice(pcols)          ! fraction of cwat that is ice
         911             :    real(r8), intent(in) :: landm(pcols)         ! Land fraction ramped over water
         912             :    real(r8), intent(in) :: seaicef(pcols)       ! sea ice fraction
         913             :    real(r8), intent(in) :: snowab(pcols)        ! rate of snow from above (kg / (m**2 * s))
         914             :    real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
         915             :
         916             : ! output arguments
         917             :    real(r8), intent(out) :: coef(pcols)          ! conversion rate (1/s)
         918             :    real(r8), intent(out) :: fwaut(pcols)         ! relative importance of liquid autoconversion (a diagnostic)
         919             :    real(r8), intent(out) :: fsaut(pcols)         ! relative importance of ice autoconversion    (a diagnostic)
         920             :    real(r8), intent(out) :: fracw(pcols)         ! relative importance of rain accreting liquid (a diagnostic)
         921             :    real(r8), intent(out) :: fsacw(pcols)         ! relative importance of snow accreting liquid (a diagnostic)
         922             :    real(r8), intent(out) :: fsaci(pcols)         ! relative importance of snow accreting ice    (a diagnostic)
         923             :    real(r8), intent(out) :: pracwo(pcols)        ! accretion of cloud water by rain (1/s)
         924             :    real(r8), intent(out) :: psacwo(pcols)        ! accretion of cloud water by snow (1/s)
         925             :    real(r8), intent(out) :: psacio(pcols)        ! accretion of cloud ice by snow (1/s)
         926             :
         927             :
         928             : ! work variables
         929             :
         930             :    integer i
         931             :    integer ii
         932             :    integer ind(pcols)
         933             :    integer ncols
         934             :
         935             :    real(r8), parameter :: degrad = 57.296_r8 ! divide by this to convert degrees to radians
         936             :    real(r8) capn                 ! local cloud particles / cm3
         937             :    real(r8) capnoice             ! local cloud particles when not over sea ice / cm3
         938             :    real(r8) ciaut                ! coefficient of autoconversion of ice (1/s)
         939             :    real(r8) cldloc(pcols)        ! non-zero amount of cloud
         940             :    real(r8) cldpr(pcols)         ! assumed cloudy volume occupied by rain and cloud
         941             :    real(r8) con1                 ! work constant
         942             :    real(r8) con2                 ! work constant
         943             :    real(r8) csacx                ! constant used for snow accreting liquid or ice
         944             : !!$   real(r8) dtice                ! interval for transition from liquid to ice
         945             :    real(r8) icemr(pcols)         ! in-cloud ice mixing ratio
         946             :    real(r8) icrit                ! threshold for autoconversion of ice
         947             :    real(r8) liqmr(pcols)         ! in-cloud liquid water mixing ratio
         948             :    real(r8) pracw                ! rate of rain accreting water
         949             :    real(r8) prlloc(pcols)        ! local rain flux in mm/day
         950             :    real(r8) prscgs(pcols)        ! local snow amount in cgs units
         951             :    real(r8) psaci                ! rate of collection of ice by snow (lin et al 1983)
         952             :    real(r8) psacw                ! rate of collection of liquid by snow (lin et al 1983)
         953             :    real(r8) psaut                ! rate of autoconversion of ice condensate
         954             :    real(r8) ptot                 ! total rate of conversion
         955             :    real(r8) pwaut                ! rate of autoconversion of liquid condensate
         956             :    real(r8) r3l                  ! volume radius
         957             :    real(r8) rainmr(pcols)        ! in-cloud rain mixing ratio
         958             :    real(r8) rat1                 ! work constant
         959             :    real(r8) rat2                 ! work constant
         960             : !!$   real(r8) rdtice               ! recipricol of dtice
         961             :    real(r8) rho(pcols)           ! density (mks units)
         962             :    real(r8) rhocgs               ! density (cgs units)
         963             :    real(r8) rlat(pcols)          ! latitude (radians)
         964             :    real(r8) snowfr               ! fraction of precipate existing as snow
         965             :    real(r8) totmr(pcols)         ! in-cloud total condensate mixing ratio
         966             :    real(r8) vfallw               ! fall speed of precipitate as liquid
         967             :    real(r8) wp                   ! weight factor used in calculating pressure dep of autoconversion
         968             :    real(r8) wsi                  ! weight factor for sea ice
         969             :    real(r8) wt                   ! fraction of ice
         970             :    real(r8) wland                ! fraction of land
         971             :
         972             : !      real(r8) csaci
         973             : !      real(r8) csacw
         974             : !      real(r8) cwaut
         975             : !      real(r8) efact
         976             : !      real(r8) lamdas
         977             : !      real(r8) lcrit
         978             : !      real(r8) rcwm
         979             : !      real(r8) r3lc2
         980             : !      real(r8) snowmr(pcols)
         981             : !      real(r8) vfalls
         982             :
         983             :    real(r8) ftot
         984             :
         985             : !     inline statement functions
         986             :    real(r8) heavy, heavym, a1, a2, heavyp, heavymp
         987             :    heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2))  ! heavyside function
         988             :    heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2))  ! modified heavyside function
         989             : !
         990             : ! New heavyside functions to perhaps address error growth problems
         991             : !
         992             :    heavyp(a1,a2) = a1/(a2+a1+1.e-36_r8)
         993             :    heavymp(a1,a2) = (a1+0.01_r8*a2)/(a2+a1+1.e-36_r8)
         994             :
         995             : !
         996             : ! find all the points where we need to do the microphysics
         997             : ! and set the output variables to zero
         998             : !
         999           0 :    ncols = 0
        1000           0 :    do i = 1,ncol
        1001           0 :       coef(i) = 0._r8
        1002           0 :       fwaut(i) = 0._r8
        1003           0 :       fsaut(i) = 0._r8
        1004           0 :       fracw(i) = 0._r8
        1005           0 :       fsacw(i) = 0._r8
        1006           0 :       fsaci(i) = 0._r8
        1007           0 :       liqmr(i) = 0._r8
        1008           0 :       rainmr(i) = 0._r8
        1009           0 :       if (cwm(i) > 1.e-20_r8) then
        1010           0 :          ncols = ncols + 1
        1011           0 :          ind(ncols) = i
        1012             :       endif
        1013             :    end do
        1014             :
        1015           0 :    do ii = 1,ncols
        1016           0 :       i = ind(ii)
        1017             : !
        1018             : ! the local cloudiness at this level
        1019             : !
        1020           0 :       cldloc(i) = max(cldmin,cldm(i))
        1021             : !
        1022             : ! a weighted mean between max cloudiness above, and this layer
        1023             : !
        1024           0 :       cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_r8)
        1025             : !
        1026             : ! decompose the suspended condensate into
        1027             : ! an incloud liquid and ice phase component
        1028             : !
        1029           0 :       totmr(i) = cwm(i)/cldloc(i)
        1030           0 :       icemr(i) = totmr(i)*fice(i)
        1031           0 :       liqmr(i) = totmr(i)*(1._r8-fice(i))
        1032             : !
        1033             : ! density
        1034             : !
        1035           0 :       rho(i) = p(i,k)/(287._r8*t(i,k))
        1036           0 :       rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
        1037             : !
        1038             : ! decompose the precipitate into a liquid and ice phase
        1039             : !
        1040           0 :       if (t(i,k) > t0) then
        1041           0 :          vfallw = convfw/sqrt(rho(i))
        1042           0 :          rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i))
        1043           0 :          snowfr = 0
        1044             : !        snowmr(i)
        1045             :       else
        1046           0 :          snowfr = 1
        1047           0 :          rainmr(i) = 0._r8
        1048             :       endif
        1049             : !     rainmr(i) = (precab(i)-snowab(i))/(rho(i)*vfallw*cldpr(i))
        1050             : !
        1051             : ! local snow amount in cgs units
        1052             : !
        1053           0 :       prscgs(i) = precab(i)/cldpr(i)*0.1_r8*snowfr
        1054             : !     prscgs(i) = snowab(i)/cldpr(i)*0.1
        1055             : !
        1056             : ! local rain amount in mm/day
        1057             : !
        1058           0 :       prlloc(i) = precab(i)*86400._r8/cldpr(i)
        1059             :    end do
        1060             :
        1061           0 :    con1 = 1._r8/(1.333_r8*pi)**0.333_r8 * 0.01_r8 ! meters
        1062             : !
        1063             : ! calculate the conversion terms
        1064             : !
        1065           0 :    call get_rlat_all_p(lchnk, ncol, rlat)
        1066             :
        1067           0 :    do ii = 1,ncols
        1068           0 :       i = ind(ii)
        1069           0 :       rhocgs = rho(i)*1.e-3_r8     ! density in cgs units
        1070             : !
        1071             : ! exponential temperature factor
        1072             : !
        1073             : !        efact = exp(0.025*(t(i,k)-t0))
        1074             : !
        1075             : ! some temperature dependent constants
        1076             : !
        1077             : !!$      wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice))
        1078           0 :       wt = fice(i)
        1079           0 :       icrit = icritc*wt + icritw*(1-wt)
        1080             : !
        1081             : ! jrm Reworked droplet number concentration algorithm
        1082             :       ! Start with pressure-dependent value appropriate for continental air
        1083             :       ! Note: reltab has a temperature dependence here
        1084           0 :       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))))
        1085             :       ! Modify for snow depth over land
        1086           0 :       capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,snowh(i)*10._r8))
        1087             :       ! Ramp between polluted value over land to clean value over ocean.
        1088           0 :       capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,1.0_r8-landm(i)))
        1089             :       ! Ramp between the resultant value and a sea ice value in the presence of ice.
        1090           0 :       capn = capn + (capnsi-capn) * min(1.0_r8,max(0.0_r8,seaicef(i)))
        1091             : ! end jrm
        1092             : !
        1093             : #ifdef DEBUG2
        1094             :       if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then
        1095             :          if (i == ilook(1)) then
        1096             :             write(iulog,*) ' findmcnew: lat, k, seaicef, landm, wp, capnoice, capn ', &
        1097             :                  lat(i), k, seaicef(i), landm(i), wp, capnoice, capn
        1098             :          endif
        1099             :       endif
        1100             : #endif
        1101             :
        1102             : !
        1103             : ! useful terms in following calculations
        1104             : !
        1105           0 :       rat1 = rhocgs/rhow
        1106           0 :       rat2 = liqmr(i)/capn
        1107           0 :       con2 = (rat1*rat2)**0.333_r8
        1108             : !
        1109             : ! volume radius
        1110             : !
        1111             : !        r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters
        1112           0 :       r3l = con1*con2
        1113             : !
        1114             : ! critical threshold for autoconversion if modified for mixed phase
        1115             : ! clouds to mimic a bergeron findeisen process
        1116             : ! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i)))
        1117             : !
        1118             : ! autoconversion of liquid
        1119             : !
        1120             : !        cwaut = 2.e-4
        1121             : !        cwaut = 1.e-3
        1122             : !        lcrit = 2.e-4
        1123             : !        lcrit = 5.e-4
        1124             : !        pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut
        1125             : !
        1126             : ! pwaut is following tripoli and cotton (and many others)
        1127             : ! we reduce the autoconversion below critpr, because these are regions where
        1128             : ! the drop size distribution is likely to imply much smaller collector drops than
        1129             : ! those relevant for a cloud distribution corresponding to the value of effc = 0.55
        1130             : ! suggested by cotton (see austin 1995 JAS, baker 1993)
        1131             :
        1132             : ! easy to follow form
        1133             : !        pwaut = capc*liqmr(i)**2*rhocgs/rhow
        1134             : !    $           *(liqmr(i)*rhocgs/(rhow*capn))**(.333)
        1135             : !    $           *heavy(r3l,r3lcrit)
        1136             : !    $           *max(0.10_r8,min(1._r8,prlloc(i)/critpr))
        1137             : ! somewhat faster form
        1138             : #define HEAVYNEW
        1139             : #ifdef HEAVYNEW
        1140             : !#ifdef PERGRO
        1141             :       pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * &
        1142           0 :               max(0.10_r8,min(1._r8,prlloc(i)/critpr))
        1143             : #else
        1144             :       pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* &
        1145             :               max(0.10_r8,min(1._r8,prlloc(i)/critpr))
        1146             : #endif
        1147             : !
        1148             : ! autoconversion of ice
        1149             : !
        1150             : !        ciaut = ciautb*efact
        1151           0 :       ciaut = ciautb
        1152             : !        psaut = capc*totmr(i)**2*rhocgs/rhoi
        1153             : !     $           *(totmr(i)*rhocgs/(rhoi*capn))**(.333)
        1154             : !
        1155             : ! autoconversion of ice condensate
        1156             : !
        1157             : #ifdef PERGRO
        1158             :       psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut
        1159             : #else
        1160           0 :       psaut = max(0._r8,icemr(i)-icrit)*ciaut
        1161             : #endif
        1162             : !
        1163             : ! collection of liquid by rain
        1164             : !
        1165             : !        pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994)
        1166           0 :       pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton)
        1167             :
        1168           0 :       pracwo(i)=pracw
        1169             :
        1170             : !!      pracw = 0.
        1171             : !
        1172             : ! the following lines calculate the slope parameter and snow mixing ratio
        1173             : ! from the precip rate using the equations found in lin et al 83
        1174             : ! in the most natural form, but it is expensive, so after some tedious
        1175             : ! algebraic manipulation you can use the cheaper form found below
        1176             : !            vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs)
        1177             : !     $               *0.01   ! convert from cm/s to m/s
        1178             : !            snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i))
        1179             : !            snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04
        1180             : !            lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25
        1181             : !            csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd)
        1182             : !
        1183             : ! coefficient for collection by snow independent of phase
        1184             : !
        1185           0 :       csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05
        1186             :
        1187             : !
        1188             : ! collection of liquid by snow (lin et al 1983)
        1189             : !
        1190           0 :       psacw = csacx*liqmr(i)*esw
        1191             : #ifdef PERGRO
        1192             : ! this is necessary for pergro
        1193             :       psacw = 0._r8
        1194             : #endif
        1195             :
        1196           0 :       psacwo(i)=psacw
        1197             :
        1198             : !
        1199             : ! collection of ice by snow (lin et al 1983)
        1200             : !
        1201           0 :       psaci = csacx*icemr(i)*esi
        1202             : !
        1203           0 :       psacio(i)=psaci
        1204             :
        1205             : ! total conversion of condensate to precipitate
        1206             : !
        1207           0 :       ptot = pwaut + psaut + pracw + psacw + psaci
        1208             : !
        1209             : ! the recipricol of cloud water amnt (or zero if no cloud water)
        1210             : !
        1211             : !         rcwm =  totmr(i)/(max(totmr(i),small)**2)
        1212             : !
        1213             : ! turn the tendency back into a loss rate (1/seconds)
        1214             : !
        1215           0 :       if (totmr(i) > 0._r8) then
        1216           0 :          coef(i) = ptot/totmr(i)
        1217             :       else
        1218           0 :          coef(i) = 0._r8
        1219             :       endif
        1220             :
        1221           0 :       if (ptot.gt.0._r8) then
        1222           0 :          fwaut(i) = pwaut/ptot
        1223           0 :          fsaut(i) = psaut/ptot
        1224           0 :          fracw(i) = pracw/ptot
        1225           0 :          fsacw(i) = psacw/ptot
        1226           0 :          fsaci(i) = psaci/ptot
        1227             :       else
        1228           0 :          fwaut(i) = 0._r8
        1229           0 :          fsaut(i) = 0._r8
        1230           0 :          fracw(i) = 0._r8
        1231           0 :          fsacw(i) = 0._r8
        1232           0 :          fsaci(i) = 0._r8
        1233             :       endif
        1234             :
        1235           0 :       ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i)
        1236             : !      if (abs(ftot-1._r8).gt.1.e-14_r8.and.ftot.ne.0._r8) then
        1237             : !         write(iulog,*) ' something is wrong in findmcnew ', ftot, &
        1238             : !              fwaut(i),fsaut(i),fracw(i),fsacw(i),fsaci(i)
        1239             : !         write(iulog,*) ' unscaled ', ptot, &
        1240             : !              pwaut,psaut,pracw,psacw,psaci
        1241             : !         write(iulog,*) ' totmr, liqmr, icemr ', totmr(i), liqmr(i), icemr(i)
        1242             : !         call endrun()
        1243             : !      endif
        1244             :    end do
        1245             : #ifdef DEBUG
        1246             :    i = icollook(nlook)
        1247             :    if (lchnk == lchnklook(nlook) ) then
        1248             :       write(iulog,*)
        1249             :       write(iulog,*) '------', k, i, lchnk
        1250             :       write(iulog,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4_r8
        1251             :       write(iulog,*) ' frac: waut,saut,racw,sacw,saci ', &
        1252             :            fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i)
        1253             :    endif
        1254             : #endif
        1255             :
        1256           0 :    return
        1257           0 : end subroutine findmcnew
        1258             :
        1259             : !-----------------------------------------------------------------------------
        1260             : ! Sets rhu to a different value poleward of +/- 50 deg latitude and
        1261             : ! levels above the tropopause if cldwat_polstrat_rhmin is specified
        1262             : ! ** This is used only for special waccm/cam-chem cases with cam4 physics **
        1263             : !-----------------------------------------------------------------------------
        1264           0 : subroutine relhum_min_adj( ncol, troplev, dlat, rhu, rhu_adj )
        1265             :
        1266             :   integer,  intent(in)  :: ncol
        1267             :   integer,  intent(in)  :: troplev(:)
        1268             :   real(r8), intent(in)  :: dlat(:)        ! latitudes in degrees
        1269             :   real(r8), intent(in)  :: rhu(:,:)
        1270             :   real(r8), intent(out) :: rhu_adj(:,:)
        1271             :
        1272             :   integer :: i,k
        1273             :
        1274           0 :   rhu_adj(:,:) = rhu(:,:)
        1275           0 :   if ( .not.do_psrhmin ) return
        1276             :
        1277           0 :   do k = 1,pver
        1278           0 :      do i = 1,ncol
        1279           0 :         if ((k .lt. troplev(i)) .and. &
        1280           0 :              ( abs( dlat(i) ) .gt. 50._r8 ) ) then
        1281           0 :            rhu_adj(i,k) = psrhmin
        1282             :         endif
        1283             :      enddo
        1284             :   enddo
        1285             :
        1286           0 : end subroutine relhum_min_adj
        1287             :
        1288             : end module cldwat

Generated by: LCOV version 1.14