module cldwat2m_micro 2,13
!---------------------------------------------------------------------------------
! Purpose:
! CAM Interface for microphysics
!
! Author: Andrew Gettelman, Hugh Morrison.
! Contributions from: Xiaohong Liu and Steve Ghan
! December 2005-May 2010
! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
! for questions contact Hugh Morrison, Andrew Gettelman
! e-mail: morrison@ucar.edu, andrew@ucar.edu
!--Pull out aerosols & activation from code? Probably to complex for now.
!---------------------------------------------------------------------------------
use shr_kind_mod
, only: r8=>shr_kind_r8
use spmd_utils, only: masterproc
use ppgrid, only: pcols, pver, pverp
use physconst, only: gravit, rair, tmelt, cpair, rh2o, r_universal, mwh2o, rhoh2o
use physconst, only: latvap, latice
use abortutils, only: endrun
use error_function, only: erf,erfc
use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice, &
vqsatd2, vqsatd2_single,polysvp
use cam_history, only: addfld, add_default, phys_decomp, outfld
use cam_logfile, only: iulog
use rad_constituents, only: rad_cnst_get_clim_info, rad_cnst_get_clim_aer_props
use phys_control, only: phys_getopts
use cldwat2m_macro, only: rhmini, rhmaxi
implicit none
private
save
logical, public :: liu_in = .true. ! True = Liu et al 2007 Ice nucleation
! False = cooper fixed ice nucleation (MG2008)
public :: ini_micro, mmicro_pcond,gamma,derf
! Private module data
logical :: wsubTKE ! If .true. (.false.), use UW's TKE (kkvh) for computing wsub
!constants remaped
real(r8), private:: g !gravity
real(r8), private:: mw !molecular weight of water
real(r8), private:: r !Dry air Gas constant
real(r8), private:: rv !water vapor gas contstant
real(r8), private:: rr !universal gas constant
real(r8), private:: cpp !specific heat of dry air
real(r8), private:: rhow !density of liquid water
real(r8), private:: xxlv ! latent heat of vaporization
real(r8), private:: xlf !latent heat of freezing
real(r8), private:: xxls !latent heat of sublimation
real(r8), private:: rhosn ! bulk density snow
real(r8), private:: rhoi ! bulk density ice
real(r8), private:: ac,bc,as,bs,ai,bi,ar,br !fall speed parameters
real(r8), private:: ci,di !ice mass-diameter relation parameters
real(r8), private:: cs,ds !snow mass-diameter relation parameters
real(r8), private:: cr,dr !drop mass-diameter relation parameters
real(r8), private:: f1s,f2s !ventilation param for snow
real(r8), private:: Eii !collection efficiency aggregation of ice
real(r8), private:: Ecc !collection efficiency
real(r8), private:: Ecr !collection efficiency cloud droplets/rain
real(r8), private:: f1r,f2r !ventilation param for rain
real(r8), private:: DCS !autoconversion size threshold
real(r8), private:: qsmall !min mixing ratio
real(r8), private:: bimm,aimm !immersion freezing
real(r8), private:: rhosu !typical 850mn air density
real(r8), private:: mi0 ! new crystal mass
real(r8), private:: rin ! radius of contact nuclei
real(r8), private:: qcvar ! 1/relative variance of sub-grid qc
real(r8), private:: pi ! pi
real(r8), private:: rn_dst1, rn_dst2, rn_dst3, rn_dst4 !dust number mean radius for contact freezing
! hm, add additional constants to help speed up code
real(r8), private:: cons1
real(r8), private:: cons2
real(r8), private:: cons3
real(r8), private:: cons4
real(r8), private:: cons5
real(r8), private:: cons6
real(r8), private:: cons7
real(r8), private:: cons8
real(r8), private:: cons9
real(r8), private:: cons10
real(r8), private:: cons11
real(r8), private:: cons12
real(r8), private:: cons13
real(r8), private:: cons14
real(r8), private:: cons15
real(r8), private:: cons16
real(r8), private:: cons17
real(r8), private:: cons18
real(r8), private:: cons19
real(r8), private:: cons20
real(r8), private:: cons21
real(r8), private:: cons22
real(r8), private:: cons23
real(r8), private:: cons24
real(r8), private:: cons25
real(r8), private:: cons26
real(r8), private:: cons27
real(r8), private:: cons28
real(r8), private:: lammini
real(r8), private:: lammaxi
real(r8), private:: lamminr
real(r8), private:: lammaxr
real(r8), private:: lammins
real(r8), private:: lammaxs
! parameters for snow/rain fraction for convective clouds
real(r8), private, parameter :: tmax_fsnow = tmelt ! max temperature for transition to convective snow
real(r8), private, parameter :: tmin_fsnow = tmelt-5._r8 ! min temperature for transition to convective snow
!needed for findsp
real(r8), private:: tt0 ! Freezing temperature
! activate parameters
integer, private:: psat
parameter (psat=6) ! number of supersaturations to calc ccn concentration
real(r8), private:: aten
real(r8), allocatable, private:: alogsig(:) ! natl log of geometric standard dev of aerosol
real(r8), allocatable, private:: exp45logsig(:)
real(r8), allocatable, private:: argfactor(:)
real(r8), allocatable, private:: amcube(:) ! cube of dry mode radius (m)
real(r8), allocatable, private:: smcrit(:) ! critical supersatuation for activation
real(r8), allocatable, private:: lnsm(:) ! ln(smcrit)
real(r8), private:: amcubesulfate(pcols) ! cube of dry mode radius (m) of sulfate
real(r8), private:: smcritsulfate(pcols) ! critical supersatuation for activation of sulfate
real(r8), allocatable, private:: amcubefactor(:) ! factors for calculating mode radius
real(r8), allocatable, private:: smcritfactor(:) ! factors for calculating critical supersaturation
real(r8), private:: super(psat)
real(r8), private:: alogten,alog2,alog3,alogaten
real(r8), private, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
(/0.02,0.05,0.1,0.2,0.5,1.0/)
real(r8), allocatable, private:: ccnfact(:,:)
real(r8), allocatable, private:: f1(:),f2(:) ! abdul-razzak functions of width
real(r8), private:: third, sixth,zero
real(r8), private:: sq2, sqpi
integer :: naer_all ! number of aerosols affecting climate
integer :: idxsul = -1 ! index in aerosol list for sulfate
integer :: idxdst1 = -1 ! index in aerosol list for dust1
integer :: idxdst2 = -1 ! index in aerosol list for dust2
integer :: idxdst3 = -1 ! index in aerosol list for dust3
integer :: idxdst4 = -1 ! index in aerosol list for dust4
integer :: idxbcphi = -1 ! index in aerosol list for Soot (BCPHIL)
! aerosol properties
character(len=20), allocatable :: aername(:)
real(r8), allocatable :: dryrad_aer(:)
real(r8), allocatable :: density_aer(:)
real(r8), allocatable :: hygro_aer(:)
real(r8), allocatable :: dispersion_aer(:)
real(r8), allocatable :: num_to_mass_aer(:)
real(r8), private:: csmin,csmax,minrefl,mindbz
contains
!===============================================================================
subroutine ini_micro 1,65
!-----------------------------------------------------------------------
!
! Purpose:
! initialize constants for the morrison microphysics
! called from stratiform.F90
!
! Author: Andrew Gettelman Dec 2005
!
!-----------------------------------------------------------------------
use cloud_fraction, only: cldfrc_getparams
integer k
integer l,m, iaer
real(r8) surften ! surface tension of water w/respect to air (N/m)
real(r8) arg
real(r8) derf
character(len=16) :: eddy_scheme = ' '
logical :: history_microphysics ! output variables for microphysics diagnostics package
! Query the PBL eddy scheme
call phys_getopts(eddy_scheme_out = eddy_scheme, &
history_microphysics_out = history_microphysics )
wsubTKE = .false.
if (trim(eddy_scheme) .eq. 'diag_TKE') wsubTKE = .true.
! Access the physical properties of the aerosols that are affecting the climate
! by using routines from the rad_constituents module.
call rad_cnst_get_clim_info(naero=naer_all)
allocate( &
aername(naer_all), &
dryrad_aer(naer_all), &
density_aer(naer_all), &
hygro_aer(naer_all), &
dispersion_aer(naer_all), &
num_to_mass_aer(naer_all) )
do iaer = 1, naer_all
call rad_cnst_get_clim_aer_props(iaer, &
aername = aername(iaer), &
dryrad_aer = dryrad_aer(iaer), &
density_aer = density_aer(iaer), &
hygro_aer = hygro_aer(iaer), &
dispersion_aer = dispersion_aer(iaer), &
num_to_mass_aer = num_to_mass_aer(iaer) )
! Look for sulfate aerosol in this list (Bulk aerosol only)
if (trim(aername(iaer)) == 'SULFATE') idxsul = iaer
if (trim(aername(iaer)) == 'DUST1') idxdst1 = iaer
if (trim(aername(iaer)) == 'DUST2') idxdst2 = iaer
if (trim(aername(iaer)) == 'DUST3') idxdst3 = iaer
if (trim(aername(iaer)) == 'DUST4') idxdst4 = iaer
if (trim(aername(iaer)) == 'BCPHIL') idxbcphi = iaer
!addfield calls for outputting aerosol number concentration
call addfld(aername(iaer)//'_m3', 'm-3', pver, 'A', 'aerosol number concentration', phys_decomp)
end do
if (masterproc) then
write(iulog,*) 'ini_micro: iaer, name, dryrad, density, hygro, dispersion, num_to_mass'
do iaer = 1, naer_all
write(iulog,*) iaer, aername(iaer), dryrad_aer(iaer), density_aer(iaer), hygro_aer(iaer), &
dispersion_aer(iaer), num_to_mass_aer(iaer)
end do
if (idxsul < 1) then
write(iulog,*) 'ini_micro: SULFATE aerosol properties NOT FOUND'
else
write(iulog,*) 'ini_micro: SULFATE aerosol properties FOUND at index ',idxsul
end if
end if
call addfld ('CCN1 ','#/cm3 ',pver, 'A','CCN concentration at S=0.02%',phys_decomp)
call addfld ('CCN2 ','#/cm3 ',pver, 'A','CCN concentration at S=0.05%',phys_decomp)
call addfld ('CCN3 ','#/cm3 ',pver, 'A','CCN concentration at S=0.1%',phys_decomp)
call addfld ('CCN4 ','#/cm3 ',pver, 'A','CCN concentration at S=0.2%',phys_decomp)
call addfld ('CCN5 ','#/cm3 ',pver, 'A','CCN concentration at S=0.5%',phys_decomp)
call addfld ('CCN6 ','#/cm3 ',pver, 'A','CCN concentration at S=1.0%',phys_decomp)
! diagnostic precip
call addfld ('QRAIN ','kg/kg ',pver, 'A','Diagnostic grid-mean rain mixing ratio' ,phys_decomp)
call addfld ('QSNOW ','kg/kg ',pver, 'A','Diagnostic grid-mean snow mixing ratio' ,phys_decomp)
call addfld ('NRAIN ','m-3 ',pver, 'A','Diagnostic grid-mean rain number conc' ,phys_decomp)
call addfld ('NSNOW ','m-3 ',pver, 'A','Diagnostic grid-mean snow number conc' ,phys_decomp)
! size of precip
call addfld ('RERCLD ','m ',pver, 'A','Diagnostic effective radius of Liquid Cloud and Rain' ,phys_decomp)
call addfld ('DSNOW ','m ',pver, 'A','Diagnostic grid-mean snow diameter' ,phys_decomp)
! precip flux variables for ciccs (instrument simulator package)
call addfld ('MGFLXPRC ','kg/m2/s ',pver+1, 'A','Diagnostic grid-mean rain flux at layer interface', phys_decomp)
call addfld ('MGFLXSNW ','kg/m2/s ',pver+1, 'A','Diagnostic grid-mean snow flux at layer interface', phys_decomp)
! diagnostic radar reflectivity, cloud-averaged
call addfld ('REFL ','DBz ',pver, 'A','94 GHz radar reflectivity' ,phys_decomp)
call addfld ('AREFL ','DBz ',pver, 'A','Average 94 GHz radar reflectivity' ,phys_decomp)
call addfld ('FREFL ','fraction ',pver, 'A','Fractional occurance of radar reflectivity' ,phys_decomp)
call addfld ('CSRFL ','DBz ',pver, 'A','94 GHz radar reflectivity (CloudSat thresholds)' ,phys_decomp)
call addfld ('ACSRFL ','DBz ',pver, 'A','Average 94 GHz radar reflectivity (CloudSat thresholds)' ,phys_decomp)
call addfld ('FCSRFL ','fraction ',pver, 'A','Fractional occurance of radar reflectivity (CloudSat thresholds)' &
,phys_decomp)
call addfld ('AREFLZ ','mm^6/m^3 ',pver, 'A','Average 94 GHz radar reflectivity' ,phys_decomp)
call addfld ('NCAL ','#/m3 ',pver, 'A','Number Concentation Activated for Liquid',phys_decomp)
call addfld ('NCAI ','#/m3 ',pver, 'A','Number Concentation Activated for Ice',phys_decomp)
call addfld ('NIHF ','#/m3 ',pver, 'A','Activated Ice Number Concentation due to homogenous freezing',phys_decomp)
call addfld ('NIDEP ','#/m3 ',pver, 'A','Activated Ice Number Concentation due to deposition nucleation',phys_decomp)
call addfld ('NIIMM ','#/m3 ',pver, 'A','Activated Ice Number Concentation due to immersion freezing',phys_decomp)
call addfld ('NIMEY ','#/m3 ',pver, 'A','Activated Ice Number Concentation due to meyers deposition',phys_decomp)
!--ag
! Average rain and snow mixing ratio (Q), number (N) and diameter (D), with frequency
call addfld ('AQRAIN ','kg/kg ',pver, 'A','Average rain mixing ratio' ,phys_decomp)
call addfld ('AQSNOW ','kg/kg ',pver, 'A','Average snow mixing ratio' ,phys_decomp)
call addfld ('ANRAIN ','m-3 ',pver, 'A','Average rain number conc' ,phys_decomp)
call addfld ('ANSNOW ','m-3 ',pver, 'A','Average snow number conc' ,phys_decomp)
call addfld ('ADRAIN ','Micron ',pver, 'A','Average rain effective Diameter' ,phys_decomp)
call addfld ('ADSNOW ','Micron ',pver, 'A','Average snow effective Diameter' ,phys_decomp)
call addfld ('FREQR ','fraction ',pver, 'A','Fractional occurance of rain' ,phys_decomp)
call addfld ('FREQS ','fraction ',pver, 'A','Fractional occurance of snow' ,phys_decomp)
if ( history_microphysics) then
call add_default ('AQSNOW ', 1, ' ')
call add_default ('CCN3 ', 1, ' ')
call add_default ('FREQR ', 1, ' ')
call add_default ('FREQS ', 1, ' ')
call add_default ('AQRAIN ', 1, ' ')
call add_default ('AQSNOW ', 1, ' ')
call add_default ('ANRAIN ', 1, ' ')
call add_default ('ANSNOW ', 1, ' ')
end if
!declarations for morrison codes (transforms variable names)
g= gravit !gravity
! mw = mwh2o / 1000._r8 !molecular weight of water
r= rair !Dry air Gas constant: note units(phys_constants are in J/K/kmol)
rv= rh2o !water vapor gas contstant
rr = r_universal !universal gas constant
cpp = cpair !specific heat of dry air
rhow = rhoh2o !density of liquid water
! latent heats
xxlv = latvap ! latent heat vaporization
xlf = latice ! latent heat freezing
xxls = xxlv + xlf ! latent heat of sublimation
! parameters below from Reisner et al. (1998)
! density parameters (kg/m3)
rhosn = 100._r8 ! bulk density snow
rhoi = 500._r8 ! bulk density ice
rhow = 1000._r8 ! bulk density liquid
! fall speed parameters, V = aD^b
! V is in m/s
! droplets
ac = 3.e7_r8
bc = 2._r8
! snow
as = 11.72_r8
bs = 0.41_r8
! cloud ice
ai = 700._r8
bi = 1._r8
! rain
ar = 841.99667_r8
br = 0.8_r8
! particle mass-diameter relationship
! currently we assume spherical particles for cloud ice/snow
! m = cD^d
pi= 3.1415927_r8
! cloud ice mass-diameter relationship
ci = rhoi*pi/6._r8
di = 3._r8
! snow mass-diameter relationship
cs = rhosn*pi/6._r8
ds = 3._r8
! drop mass-diameter relationship
cr = rhow*pi/6._r8
dr = 3._r8
! ventilation parameters for snow
! hall and prupacher
f1s = 0.86_r8
f2s = 0.28_r8
! collection efficiency, aggregation of cloud ice and snow
Eii = 0.1_r8
! collection efficiency, accretion of cloud water by rain
Ecr = 1.0_r8
! ventilation constants for rain
f1r = 0.78_r8
f2r = 0.32_r8
! autoconversion size threshold for cloud ice to snow (m)
Dcs = 325.e-6_r8
! smallest mixing ratio considered in microphysics
qsmall = 1.e-18_r8
! immersion freezing parameters, bigg 1953
bimm = 100._r8
aimm = 0.66_r8
! contact freezing due to dust
! dust number mean radius (m), Zender et al JGR 2003 assuming number mode radius of 0.6 micron, sigma=2
rn_dst1=0.258e-6_r8
rn_dst2=0.717e-6_r8
rn_dst3=1.576e-6_r8
rn_dst4=3.026e-6_r8
! typical air density at 850 mb
rhosu = 85000._r8/(rair * tmelt)
! mass of new crystal due to aerosol freezing and growth (kg)
mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)
! radius of contact nuclei aerosol (m)
rin = 0.1e-6_r8
! 1 / relative variance of sub-grid cloud water distribution
! see morrison and gettelman, 2007, J. Climate for details
qcvar = 2._r8
! freezing temperature
tt0=273.15_r8
! set parameters for droplet activation, following abdul-razzak and ghan 2000, JGR
! mathematical constants
zero=0._r8
third=1./3._r8
sixth=1./6._r8
sq2=sqrt(2._r8)
pi=4._r8*atan(1.0_r8)
sqpi=sqrt(pi)
surften=0.076_r8
aten=2.*mwh2o*surften/(r_universal*tt0*rhoh2o)
alogaten=log(aten)
alog2=log(2._r8)
alog3=log(3._r8)
super(:)=0.01*supersat(:)
allocate(alogsig(naer_all), &
exp45logsig(naer_all), &
argfactor(naer_all), &
amcube(naer_all), &
smcrit(naer_all), &
lnsm(naer_all), &
amcubefactor(naer_all), &
smcritfactor(naer_all), &
ccnfact(psat,naer_all), &
f1(naer_all), &
f2(naer_all) )
do m=1,naer_all
! use only if width of size distribution is prescribed
alogsig(m)=log(dispersion_aer(m))
exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
argfactor(m)=2./(3.*sqrt(2.)*alogsig(m))
f1(m)=0.5*exp(2.5*alogsig(m)*alogsig(m))
f2(m)=1.+0.25*alogsig(m)
amcubefactor(m)=3._r8/(4._r8*pi*exp45logsig(m)*density_aer(m))
smcritfactor(m)=2._r8*aten*sqrt(aten/(27._r8*max(1.e-10_r8,hygro_aer(m))))
amcube(m)=amcubefactor(m)/(num_to_mass_aer(m))
if(hygro_aer(m).gt.1.e-10) then
smcrit(m)=smcritfactor(m)/sqrt(amcube(m))
else
smcrit(m)=100.
endif
lnsm(m)=log(smcrit(m))
do l=1,psat
arg=argfactor(m)*log(smcrit(m)/super(l))
if(arg<2)then
if(arg<-2)then
ccnfact(l,m)=1.e-6
else
ccnfact(l,m)=1.e-6_r8*0.5_r8*erfc(arg)
endif
else
ccnfact(l,m) = 0.
endif
enddo
end do
!Range of cloudsat reflectivities (dBz)
csmin= -30._r8
csmax= 26._r8
mindbz = -99._r8
! minrefl = 10._r8**(mindbz/10._r8) hard code below
minrefl = 1.26e-10_r8
! hm, add constants to help speed up code
cons1=gamma(1._r8+di)
cons2=gamma(qcvar+2.47_r8)
cons3=gamma(qcvar)
cons4=gamma(1._r8+br)
cons5=gamma(4._r8+br)
cons6=gamma(1._r8+ds)
cons7=gamma(1._r8+bs)
cons8=gamma(4._r8+bs)
cons9=gamma(qcvar+2._r8)
cons10=gamma(qcvar+1._r8)
cons11=gamma(3._r8+bs)
cons12=gamma(qcvar+1.15_r8)
cons13=gamma(5._r8/2._r8+br/2._r8)
cons14=gamma(5._r8/2._r8+bs/2._r8)
cons15=gamma(qcvar+bc/3._r8)
cons16=gamma(1._r8+bi)
cons17=gamma(4._r8+bi)
cons18=qcvar**2.47_r8
cons19=qcvar**2
cons20=qcvar**1.15_r8
cons21=qcvar**(bc/3._r8)
cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
cons23=dcs**3
cons24=dcs**2
cons25=dcs**bs
! cons26=pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)
cons27=xxlv**2
cons28=xxls**2
lammaxi = 1._r8/10.e-6_r8
lammini = 1._r8/(2._r8*dcs)
lammaxr = 1._r8/20.e-6_r8
lamminr = 1._r8/500.e-6_r8
lammaxs = 1._r8/10.e-6_r8
lammins = 1._r8/2000.e-6_r8
return
end subroutine ini_micro
!===============================================================================
!microphysics routine for each timestep goes here...
subroutine mmicro_pcond ( & 1,77
lchnk, ncol, deltatin, tn, ttend, &
qn, qtend, cwtend, qc, qi, &
nc, ni, p, pdel, cldn, &
liqcldf, icecldf, &
cldo, pint, rpdel, zm, omega, &
#ifdef MODAL_AERO
qaer, cflx, qaertend, qqcw,dgnumwet,dgnum, &
rate1ord_cw2pr_st, & ! rce 2010/05/01
#else
aer_mmr, &
#endif
rhdfda, rhu00, fice, tlat, qvlat, &
qctend, qitend, nctend, nitend, effc, &
effc_fn, effi, prect, preci, kkvh, &
tke, turbtype, smaw, wsub, wsubi, nevapr, evapsnow, &
prain, prodsnow, cmeout, deffi, pgamrad, &
lamcrad,qsout,dsout, &
qcsevap,qisevap,qvres,cmeiout, &
vtrmc,vtrmi,qcsedten,qisedten, &
prao,prco,mnuccco,mnuccto,msacwio,psacwso,&
bergso,bergo,melto,homoo,qcreso,prcio,praio,qireso,&
mnuccro,pracso,meltsdt,frzrdt)
!Author: Hugh Morrison, Andrew Gettelman, NCAR
! e-mail: morrison@ucar.edu, andrew@ucar.edu
use wv_saturation, only: vqsatd, vqsatd_water
use constituents, only: pcnst
#ifdef MODAL_AERO
use ndrop, only: dropmixnuc
use modal_aero_data, only: numptr_amode, modeptr_accum, modeptr_coarse, modeptr_aitken, &
lptr_dust_a_amode,lptr_nacl_a_amode,maxd_amode
#endif
! input arguments
integer, intent(in) :: lchnk
integer, intent(in) :: ncol
real(r8), intent(in) :: deltatin ! time step (s)
real(r8), intent(in) :: tn(pcols,pver) ! input temperature (K)
real(r8), intent(in) :: ttend(pcols,pver) ! non-microphysical temperature tendency (K/s)
real(r8), intent(in) :: qn(pcols,pver) ! input h20 vapor mixing ratio (kg/kg)
real(r8), intent(in) :: qtend(pcols,pver) ! non-microphysical qv tendency (1/s)
real(r8), intent(in) :: cwtend(pcols,pver) ! non-microphysical tendency of cloud water (1/s)
! note: all input cloud variables are grid-averaged
real(r8), intent(inout) :: qc(pcols,pver) ! cloud water mixing ratio (kg/kg)
real(r8), intent(inout) :: qi(pcols,pver) ! cloud ice mixing ratio (kg/kg)
real(r8), intent(inout) :: nc(pcols,pver) ! cloud water number conc (1/kg)
real(r8), intent(inout) :: ni(pcols,pver) ! cloud ice number conc (1/kg)
real(r8), intent(in) :: p(pcols,pver) ! air pressure (pa)
real(r8), intent(in) :: pdel(pcols,pver) ! pressure difference across level (pa)
real(r8), intent(in) :: cldn(pcols,pver) ! cloud fraction
real(r8), intent(in) :: icecldf(pcols,pver) ! ice cloud fraction
real(r8), intent(in) :: liqcldf(pcols,pver) ! liquid cloud fraction
real(r8), intent(inout) :: cldo(pcols,pver) ! old cloud fraction
real(r8), intent(in) :: pint(pcols,pverp) ! air pressure layer interfaces (pa)
real(r8), intent(in) :: rpdel(pcols,pver) ! inverse pressure difference across level (pa)
real(r8), intent(in) :: zm(pcols,pver) ! geopotential height of model levels (m)
real(r8), intent(in) :: omega(pcols,pver) ! vertical velocity (Pa/s)
#ifdef MODAL_AERO
real(r8), intent(in) :: qaer(pcols,pver,pcnst) ! aerosol number and mass mixing ratios
real(r8), intent(in) :: cflx(pcols,pcnst) ! constituent flux from surface
real(r8), intent(inout) :: qaertend(pcols,pver,pcnst) ! qaer tendency (1/s)
real(r8), intent(inout) :: qqcw(pcols,pver,pcnst) ! cloud-borne aerosol
real(r8), intent(in) :: dgnumwet(pcols,pver,maxd_amode) ! aerosol mode diameter
real(r8), intent(in) :: dgnum(pcols,pver,maxd_amode) ! aerosol mode dry diameter
real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for direct cw to precip conversion ! rce 2010/05/01
#else
real(r8), intent(in) :: aer_mmr(:,:,:) ! aerosol mass mixing ratio
#endif
real(r8), intent(in) :: rhdfda(pcols,pver) ! dA/dRH
real(r8), intent(in) :: rhu00(pcols,pver) ! threshold rh for cloud
real(r8), intent(in) :: fice(pcols,pver) ! fraction of ice/liquid
real(r8), intent(in) :: kkvh(pcols,pver+1) ! vertical eddy diff coef (m2 s-1)
real(r8), intent(in) :: tke(pcols,pver+1) ! TKE from the UW PBL scheme (m2 s-2)
real(r8), intent(in) :: turbtype(pcols,pver+1) ! Turbulence type at each interface
real(r8), intent(in) :: smaw(pcols,pver+1) ! Normalized instability function of momentum. 0 <= <= 4.964 and 1 at neutral condition.
! output arguments
real(r8), intent(out) :: tlat(pcols,pver) ! latent heating rate (W/kg)
real(r8), intent(out) :: qvlat(pcols,pver) ! microphysical tendency qv (1/s)
real(r8), intent(out) :: qctend(pcols,pver) ! microphysical tendency qc (1/s)
real(r8), intent(out) :: qitend(pcols,pver) ! microphysical tendency qi (1/s)
real(r8), intent(out) :: nctend(pcols,pver) ! microphysical tendency nc (1/(kg*s))
real(r8), intent(out) :: nitend(pcols,pver) ! microphysical tendency ni (1/(kg*s))
real(r8), intent(out) :: effc(pcols,pver) ! droplet effective radius (micron)
real(r8), intent(out) :: effc_fn(pcols,pver) ! droplet effective radius, assuming nc = 1.e8 kg-1
real(r8), intent(out) :: effi(pcols,pver) ! cloud ice effective radius (micron)
real(r8), intent(out) :: prect(pcols) ! surface precip rate (m/s)
real(r8), intent(out) :: preci(pcols) ! cloud ice/snow precip rate (m/s)
real(r8), intent(out) :: wsub(pcols,pver) ! diagnosed sub-grid vertical velocity st. dev. (m/s)
real(r8), intent(out) :: wsubi(pcols,pver) ! diagnosed sub-grid vertical velocity ice (m/s)
real(r8), intent(out) :: nevapr(pcols,pver) ! evaporation rate of rain + snow
real(r8), intent(out) :: evapsnow(pcols,pver)! sublimation rate of snow
real(r8), intent(out) :: prain(pcols,pver) ! production of rain + snow
real(r8), intent(out) :: prodsnow(pcols,pver)! production of snow
real(r8), intent(out) :: cmeout(pcols,pver) ! evap/sub of cloud
real(r8), intent(out) :: deffi(pcols,pver) ! ice effective diameter for optics (radiation)
real(r8), intent(out) :: pgamrad(pcols,pver) ! ice gamma parameter for optics (radiation)
real(r8), intent(out) :: lamcrad(pcols,pver) ! slope of droplet distribution for optics (radiation)
real(r8), intent(out) :: qcsevap(pcols,pver) ! cloud water evaporation due to sedimentation
real(r8), intent(out) :: qisevap(pcols,pver) ! cloud ice sublimation due to sublimation
real(r8), intent(out) :: qvres(pcols,pver) ! residual condensation term to ensure RH < 100%
real(r8), intent(out) :: cmeiout(pcols,pver) ! grid-mean cloud ice sub/dep
real(r8), intent(out) :: vtrmc(pcols,pver) ! mass-weighted cloud water fallspeed
real(r8), intent(out) :: vtrmi(pcols,pver) ! mass-weighted cloud ice fallspeed
real(r8), intent(out) :: qcsedten(pcols,pver) ! qc sedimentation tendency
real(r8), intent(out) :: qisedten(pcols,pver) ! qi sedimentation tendency
! microphysical process rates for output (mixing ratio tendencies)
real(r8), intent(out) :: prao(pcols,pver) ! accretion of cloud by rain
real(r8), intent(out) :: prco(pcols,pver) ! autoconversion of cloud to rain
real(r8), intent(out) :: mnuccco(pcols,pver) ! mixing rat tend due to immersion freezing
real(r8), intent(out) :: mnuccto(pcols,pver) ! mixing ratio tend due to contact freezing
real(r8), intent(out) :: msacwio(pcols,pver) ! mixing ratio tend due to H-M splintering
real(r8), intent(out) :: psacwso(pcols,pver) ! collection of cloud water by snow
real(r8), intent(out) :: bergso(pcols,pver) ! bergeron process on snow
real(r8), intent(out) :: bergo(pcols,pver) ! bergeron process on cloud ice
real(r8), intent(out) :: melto(pcols,pver) ! melting of cloud ice
real(r8), intent(out) :: homoo(pcols,pver) ! homogeneos freezign cloud water
real(r8), intent(out) :: qcreso(pcols,pver) ! residual cloud condensation due to removal of excess supersat
real(r8), intent(out) :: prcio(pcols,pver) ! autoconversion of cloud ice to snow
real(r8), intent(out) :: praio(pcols,pver) ! accretion of cloud ice by snow
real(r8), intent(out) :: qireso(pcols,pver) ! residual ice deposition due to removal of excess supersat
real(r8), intent(out) :: mnuccro(pcols,pver) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
real(r8), intent(out) :: pracso (pcols,pver) ! mixing ratio tendency due to accretion of rain by snow (1/s)
real(r8), intent(out) :: meltsdt(pcols,pver) ! latent heating rate due to melting of snow (W/kg)
real(r8), intent(out) :: frzrdt (pcols,pver) ! latent heating rate due to homogeneous freezing of rain (W/kg)
real(r8) :: tkem(pcols,pver) ! Layer-mean TKE
real(r8) :: smm(pcols,pver) ! Layer-mean instability function
! local workspace
! all units mks unless otherwise stated
! temporary variables for sub-stepping
real(r8) :: t1(pcols,pver)
real(r8) :: q1(pcols,pver)
real(r8) :: qc1(pcols,pver)
real(r8) :: qi1(pcols,pver)
real(r8) :: nc1(pcols,pver)
real(r8) :: ni1(pcols,pver)
real(r8) :: tlat1(pcols,pver)
real(r8) :: qvlat1(pcols,pver)
real(r8) :: qctend1(pcols,pver)
real(r8) :: qitend1(pcols,pver)
real(r8) :: nctend1(pcols,pver)
real(r8) :: nitend1(pcols,pver)
real(r8) :: prect1(pcols)
real(r8) :: preci1(pcols)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(r8) :: deltat ! sub-time step (s)
real(r8) :: omsm ! number near unity for round-off issues
real(r8) :: dto2 ! dt/2 (s)
real(r8) :: mincld ! minimum allowed cloud fraction
real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg)
real(r8) :: t(pcols,pver) ! temperature (K)
real(r8) :: rho(pcols,pver) ! air density (kg m-3)
real(r8) :: dv(pcols,pver) ! diffusivity of water vapor in air
real(r8) :: mu(pcols,pver) ! viscocity of air
real(r8) :: sc(pcols,pver) ! schmidt number
real(r8) :: kap(pcols,pver) ! thermal conductivity of air
real(r8) :: rhof(pcols,pver) ! air density correction factor for fallspeed
real(r8) :: cldmax(pcols,pver) ! precip fraction assuming maximum overlap
real(r8) :: cldm(pcols,pver) ! cloud fraction
real(r8) :: icldm(pcols,pver) ! ice cloud fraction
real(r8) :: lcldm(pcols,pver) ! liq cloud fraction
real(r8) :: icwc(pcols) ! in cloud water content (liquid+ice)
real(r8) :: calpha(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cbeta(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cbetah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cgamma(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cgamah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: rcgama(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cmec1(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cmec2(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cmec3(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: cmec4(pcols) ! parameter for cond/evap (Zhang et al. 2003)
real(r8) :: qtmp ! dummy qv
real(r8) :: dum ! temporary dummy variable
real(r8) :: cme(pcols,pver) ! total (liquid+ice) cond/evap rate of cloud
real(r8) :: cmei(pcols,pver) ! dep/sublimation rate of cloud ice
real(r8) :: cwml(pcols,pver) ! cloud water mixing ratio
real(r8) :: cwmi(pcols,pver) ! cloud ice mixing ratio
real(r8) :: nnuccd(pver) ! ice nucleation rate from deposition/cond.-freezing
real(r8) :: mnuccd(pver) ! mass tendency from ice nucleation
real(r8) :: qcld ! total cloud water
real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
real(r8) :: nctend_mixnuc(pcols,pver)
real(r8) :: arg ! argument of erfc
#ifdef MODAL_AERO
real(r8) :: qcsinksum_rate1ord(pver) ! sum over iterations of cw to precip sink ! rce 2010/05/01
real(r8) :: qcsum_rate1ord(pver) ! sum over iterations of cloud water ! rce 2010/05/01
#endif
real(r8) :: alpha
real(r8) :: dum1,dum2 !general dummy variables
real(r8) :: nihf2,niimm2,nidep2,nimey2 !dummy variables for activated ice nuclei by diffferent processes
real(r8) :: npccn(pver) ! droplet activation rate
real(r8) :: qcic(pcols,pver) ! in-cloud cloud liquid mixing ratio
real(r8) :: qiic(pcols,pver) ! in-cloud cloud ice mixing ratio
real(r8) :: qniic(pcols,pver) ! in-precip snow mixing ratio
real(r8) :: qric(pcols,pver) ! in-precip rain mixing ratio
real(r8) :: ncic(pcols,pver) ! in-cloud droplet number conc
real(r8) :: niic(pcols,pver) ! in-cloud cloud ice number conc
real(r8) :: nsic(pcols,pver) ! in-precip snow number conc
real(r8) :: nric(pcols,pver) ! in-precip rain number conc
real(r8) :: lami(pver) ! slope of cloud ice size distr
real(r8) :: n0i(pver) ! intercept of cloud ice size distr
real(r8) :: lamc(pver) ! slope of cloud liquid size distr
real(r8) :: n0c(pver) ! intercept of cloud liquid size distr
real(r8) :: lams(pver) ! slope of snow size distr
real(r8) :: n0s(pver) ! intercept of snow size distr
real(r8) :: lamr(pver) ! slope of rain size distr
real(r8) :: n0r(pver) ! intercept of rain size distr
real(r8) :: cdist1(pver) ! size distr parameter to calculate droplet freezing
! combined size of precip & cloud drops
real(r8) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud
real(r8) :: arcld(pcols,pver) ! averaging control flag
real(r8) :: Actmp !area cross section of drops
real(r8) :: Artmp !area cross section of rain
real(r8) :: pgam(pver) ! spectral width parameter of droplet size distr
real(r8) :: lammax ! maximum allowed slope of size distr
real(r8) :: lammin ! minimum allowed slope of size distr
real(r8) :: nacnt ! number conc of contact ice nuclei
real(r8) :: mnuccc(pver) ! mixing ratio tendency due to freezing of cloud water
real(r8) :: nnuccc(pver) ! number conc tendency due to freezing of cloud water
real(r8) :: mnucct(pver) ! mixing ratio tendency due to contact freezing of cloud water
real(r8) :: nnucct(pver) ! number conc tendency due to contact freezing of cloud water
real(r8) :: msacwi(pver) ! mixing ratio tendency due to HM ice multiplication
real(r8) :: nsacwi(pver) ! number conc tendency due to HM ice multiplication
real(r8) :: prc(pver) ! qc tendency due to autoconversion of cloud droplets
real(r8) :: nprc(pver) ! number conc tendency due to autoconversion of cloud droplets
real(r8) :: nprc1(pver) ! qr tendency due to autoconversion of cloud droplets
real(r8) :: nsagg(pver) ! ns tendency due to self-aggregation of snow
real(r8) :: dc0 ! mean size droplet size distr
real(r8) :: ds0 ! mean size snow size distr (area weighted)
real(r8) :: eci ! collection efficiency for riming of snow by droplets
real(r8) :: psacws(pver) ! mixing rat tendency due to collection of droplets by snow
real(r8) :: npsacws(pver) ! number conc tendency due to collection of droplets by snow
real(r8) :: uni ! number-weighted cloud ice fallspeed
real(r8) :: umi ! mass-weighted cloud ice fallspeed
real(r8) :: uns(pver) ! number-weighted snow fallspeed
real(r8) :: ums(pver) ! mass-weighted snow fallspeed
real(r8) :: unr(pver) ! number-weighted rain fallspeed
real(r8) :: umr(pver) ! mass-weighted rain fallspeed
real(r8) :: unc ! number-weighted cloud droplet fallspeed
real(r8) :: umc ! mass-weighted cloud droplet fallspeed
real(r8) :: pracs(pver) ! mixing rat tendency due to collection of rain by snow
real(r8) :: npracs(pver) ! number conc tendency due to collection of rain by snow
real(r8) :: mnuccr(pver) ! mixing rat tendency due to freezing of rain
real(r8) :: nnuccr(pver) ! number conc tendency due to freezing of rain
real(r8) :: pra(pver) ! mixing rat tendnency due to accretion of droplets by rain
real(r8) :: npra(pver) ! nc tendnency due to accretion of droplets by rain
real(r8) :: nragg(pver) ! nr tendency due to self-collection of rain
real(r8) :: prci(pver) ! mixing rat tendency due to autoconversion of cloud ice to snow
real(r8) :: nprci(pver) ! number conc tendency due to autoconversion of cloud ice to snow
real(r8) :: prai(pver) ! mixing rat tendency due to accretion of cloud ice by snow
real(r8) :: nprai(pver) ! number conc tendency due to accretion of cloud ice by snow
real(r8) :: qvs ! liquid saturation vapor mixing ratio
real(r8) :: qvi ! ice saturation vapor mixing ratio
real(r8) :: dqsdt ! change of sat vapor mixing ratio with temperature
real(r8) :: dqsidt ! change of ice sat vapor mixing ratio with temperature
real(r8) :: ab ! correction factor for rain evap to account for latent heat
real(r8) :: qclr ! water vapor mixing ratio in clear air
real(r8) :: abi ! correction factor for snow sublimation to account for latent heat
real(r8) :: epss ! 1/ sat relaxation timescale for snow
real(r8) :: epsr ! 1/ sat relaxation timescale for rain
real(r8) :: pre(pver) ! rain mixing rat tendency due to evaporation
real(r8) :: prds(pver) ! snow mixing rat tendency due to sublimation
real(r8) :: qce ! dummy qc for conservation check
real(r8) :: qie ! dummy qi for conservation check
real(r8) :: nce ! dummy nc for conservation check
real(r8) :: nie ! dummy ni for conservation check
real(r8) :: ratio ! parameter for conservation check
real(r8) :: dumc(pcols,pver) ! dummy in-cloud qc
real(r8) :: dumnc(pcols,pver) ! dummy in-cloud nc
real(r8) :: dumi(pcols,pver) ! dummy in-cloud qi
real(r8) :: dumni(pcols,pver) ! dummy in-cloud ni
real(r8) :: dums(pcols,pver) ! dummy in-cloud snow mixing rat
real(r8) :: dumns(pcols,pver) ! dummy in-cloud snow number conc
real(r8) :: dumr(pcols,pver) ! dummy in-cloud rain mixing rat
real(r8) :: dumnr(pcols,pver) ! dummy in-cloud rain number conc
! below are parameters for cloud water and cloud ice sedimentation calculations
real(r8) :: fr(pver)
real(r8) :: fnr(pver)
real(r8) :: fc(pver)
real(r8) :: fnc(pver)
real(r8) :: fi(pver)
real(r8) :: fni(pver)
real(r8) :: fs(pver)
real(r8) :: fns(pver)
real(r8) :: faloutr(pver)
real(r8) :: faloutnr(pver)
real(r8) :: faloutc(pver)
real(r8) :: faloutnc(pver)
real(r8) :: falouti(pver)
real(r8) :: faloutni(pver)
real(r8) :: falouts(pver)
real(r8) :: faloutns(pver)
real(r8) :: faltndr
real(r8) :: faltndnr
real(r8) :: faltndc
real(r8) :: faltndnc
real(r8) :: faltndi
real(r8) :: faltndni
real(r8) :: faltnds
real(r8) :: faltndns
real(r8) :: faltndqie
real(r8) :: faltndqce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(r8) :: relhum(pcols,pver) ! relative humidity
real(r8) :: csigma(pcols) ! parameter for cond/evap of cloud water/ice
real(r8) :: rgvm ! max fallspeed for all species
real(r8) :: arn(pcols,pver) ! air density corrected rain fallspeed parameter
real(r8) :: asn(pcols,pver) ! air density corrected snow fallspeed parameter
real(r8) :: acn(pcols,pver) ! air density corrected cloud droplet fallspeed parameter
real(r8) :: ain(pcols,pver) ! air density corrected cloud ice fallspeed parameter
real(r8) :: nsubi(pver) ! evaporation of cloud ice number
real(r8) :: nsubc(pver) ! evaporation of droplet number
real(r8) :: nsubs(pver) ! evaporation of snow number
real(r8) :: nsubr(pver) ! evaporation of rain number
real(r8) :: mtime ! factor to account for droplet activation timescale
real(r8) :: dz(pcols,pver) ! height difference across model vertical level
!fice variable
real(r8) :: nfice(pcols,pver)
!add variables for rain and snow flux at layer interfaces
real(r8) :: rflx(pcols,pver+1)
real(r8) :: sflx(pcols,pver+1)
!! also add precip flux variables for sub-stepping
real(r8) :: rflx1(pcols,pver+1)
real(r8) :: sflx1(pcols,pver+1)
! returns from function/subroutine calls
real(r8) :: tsp(pcols,pver) ! saturation temp (K)
real(r8) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg)
real(r8) :: qsphy(pcols,pver) ! saturation mixing ratio (kg/kg): hybrid rh
real(r8) :: qs(pcols) ! liquid-ice weighted sat mixing rat (kg/kg)
real(r8) :: es(pcols) ! liquid-ice weighted sat vapor press (pa)
real(r8) :: esl(pcols,pver) ! liquid sat vapor pressure (pa)
real(r8) :: esi(pcols,pver) ! ice sat vapor pressure (pa)
real(r8) :: gammas(pcols) ! parameter for cond/evap of cloud water
! sum of source/sink terms for diagnostic precip
real(r8) :: qnitend(pcols,pver) ! snow mixing ratio source/sink term
real(r8) :: nstend(pcols,pver) ! snow number concentration source/sink term
real(r8) :: qrtend(pcols,pver) ! rain mixing ratio source/sink term
real(r8) :: nrtend(pcols,pver) ! rain number concentration source/sink term
real(r8) :: qrtot ! vertically-integrated rain mixing rat source/sink term
real(r8) :: nrtot ! vertically-integrated rain number conc source/sink term
real(r8) :: qstot ! vertically-integrated snow mixing rat source/sink term
real(r8) :: nstot ! vertically-integrated snow number conc source/sink term
! new terms for Bergeron process
real(r8) :: dumnnuc ! provisional ice nucleation rate (for calculating bergeron)
real(r8) :: ninew ! provisional cloud ice number conc (for calculating bergeron)
real(r8) :: qinew ! provisional cloud ice mixing ratio (for calculating bergeron)
real(r8) :: qvl ! liquid sat mixing ratio
real(r8) :: epsi ! 1/ sat relaxation timecale for cloud ice
real(r8) :: prd ! provisional deposition rate of cloud ice at water sat
real(r8) :: berg(pcols,pver) ! mixing rat tendency due to bergeron process for cloud ice
real(r8) :: bergs(pver) ! mixing rat tendency due to bergeron process for snow
!bergeron terms
real(r8) :: bergtsf !bergeron timescale to remove all liquid
real(r8) :: rhin !modified RH for vapor deposition
! new aerosol variables
real(r8), allocatable :: naermod(:) ! aerosol number concentration (/m3)
real(r8), allocatable :: naer2(:,:,:) ! new aerosol number concentration (/m3)
real(r8), allocatable :: maerosol(:,:) ! aerosol mass conc (kg/m3)
real(r8) naer(pcols)
real(r8) ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat
character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
! diagnostic rain/snow for output to history
! values are in-precip (local) !!!!
real(r8) :: qrout(pcols,pver) ! rain mixing ratio (kg/kg)
real(r8) :: nrout(pcols,pver) ! rain number concentration (1/m3)
real(r8) :: nsout(pcols,pver) ! snow number concentration (1/m3)
real(r8), intent(out) :: dsout(pcols,pver) ! snow diameter (m)
real(r8), intent(out) :: qsout(pcols,pver) ! snow mixing ratio (kg/kg)
!averageed rain/snow for history
real(r8) :: qrout2(pcols,pver)
real(r8) :: qsout2(pcols,pver)
real(r8) :: nrout2(pcols,pver)
real(r8) :: nsout2(pcols,pver)
real(r8) :: freqs(pcols,pver)
real(r8) :: freqr(pcols,pver)
real(r8) :: dumfice
real(r8) :: drout2(pcols,pver) ! mean rain particle diameter (m)
real(r8) :: dsout2(pcols,pver) ! mean snow particle diameter (m)
!ice nucleation, droplet activation
real(r8) :: dum2i(pcols,pver) ! number conc of ice nuclei available (1/kg)
real(r8) :: dum2l(pcols,pver) ! number conc of CCN (1/kg)
real(r8) :: ncmax
real(r8) :: nimax
!output fields for number conc
real(r8) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3)
real(r8) :: ncal(pcols,pver) ! output number conc of CCN (1/m3)
!output for ice nucleation
real(r8) :: nimey(pcols,pver) !output number conc of ice nuclei due to meyers deposition (1/m3)
real(r8) :: nihf(pcols,pver) !output number conc of ice nuclei due to heterogenous freezing (1/m3)
real(r8) :: nidep(pcols,pver) !output number conc of ice nuclei due to deoposion nucleation (hetero nuc) (1/m3)
real(r8) :: niimm(pcols,pver) !output number conc of ice nuclei due to immersion freezing (hetero nuc) (1/m3)
! loop array variables
integer i,k,nstep,n, l
integer ii,kk, m
integer, allocatable :: ntype(:)
integer ftrue ! integer used to determine cloud depth
! loop variables for sub-step solution
integer iter,it,ltrue(pcols)
! used in contact freezing via dust particles
real(r8) tcnt, viscosity, mfp
real(r8) slip1, slip2, slip3, slip4
real(r8) dfaer1, dfaer2, dfaer3, dfaer4
real(r8) nacon1,nacon2,nacon3,nacon4
! used in ice effective radius
real(r8) bbi, cci, ak, iciwc, rvi
! used in Bergeron processe and water vapor deposition
real(r8) Tk, deles, Aprpr, Bprpr, Cice, qi0, Crate, qidep
! mean cloud fraction over the time step
real(r8) cldmw(pcols,pver)
! used in secondary ice production
real(r8) ni_secp
! variabels to check for RH after rain evap
real(r8) :: esn
real(r8) :: qsn
real(r8) :: ttmp
real(r8) :: refl(pcols,pver) ! analytic radar reflectivity
real(r8) :: rainrt(pcols,pver) ! rain rate for reflectivity calculation
real(r8) :: rainrt1(pcols,pver)
real(r8) :: csrfl(pcols,pver) !cloudsat reflectivity
real(r8) :: arefl(pcols,pver) !average reflectivity will zero points outside valid range
real(r8) :: acsrfl(pcols,pver) !cloudsat average
real(r8) :: frefl(pcols,pver)
real(r8) :: fcsrfl(pcols,pver)
real(r8) :: areflz(pcols,pver) !average reflectivity in z.
real(r8) :: tmp
real(r8) dmc,ssmc,dstrn ! variables for modal scheme.
real(r8), parameter :: cdnl = 0.e6_r8 ! cloud droplet number limiter
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! initialize output fields for number conc qand ice nucleation
ncai(1:ncol,1:pver)=0._r8
ncal(1:ncol,1:pver)=0._r8
nimey(1:ncol,1:pver)=0._r8
nihf(1:ncol,1:pver)=0._r8
nidep(1:ncol,1:pver)=0._r8
niimm(1:ncol,1:pver)=0._r8
!Initialize rain size
rercld(1:ncol,1:pver)=0._r8
arcld(1:ncol,1:pver)=0._r8
! TKE initialization
tkem(1:ncol,1:pver)=0._r8
smm(1:ncol,1:pver)=0._r8
!initialize radiation output variables
pgamrad(1:ncol,1:pver)=0._r8 ! liquid gamma parameter for optics (radiation)
lamcrad(1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
deffi (1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
!initialize radiation output variables
!initialize water vapor tendency term output
qcsevap(1:ncol,1:pver)=0._r8
qisevap(1:ncol,1:pver)=0._r8
qvres (1:ncol,1:pver)=0._r8
cmeiout (1:ncol,1:pver)=0._r8
vtrmc (1:ncol,1:pver)=0._r8
vtrmi (1:ncol,1:pver)=0._r8
qcsedten (1:ncol,1:pver)=0._r8
qisedten (1:ncol,1:pver)=0._r8
prao(1:ncol,1:pver)=0._r8
prco(1:ncol,1:pver)=0._r8
mnuccco(1:ncol,1:pver)=0._r8
mnuccto(1:ncol,1:pver)=0._r8
msacwio(1:ncol,1:pver)=0._r8
psacwso(1:ncol,1:pver)=0._r8
bergso(1:ncol,1:pver)=0._r8
bergo(1:ncol,1:pver)=0._r8
melto(1:ncol,1:pver)=0._r8
homoo(1:ncol,1:pver)=0._r8
qcreso(1:ncol,1:pver)=0._r8
prcio(1:ncol,1:pver)=0._r8
praio(1:ncol,1:pver)=0._r8
qireso(1:ncol,1:pver)=0._r8
mnuccro(1:ncol,1:pver)=0._r8
pracso (1:ncol,1:pver)=0._r8
meltsdt(1:ncol,1:pver)=0._r8
frzrdt (1:ncol,1:pver)=0._r8
allocate(naermod(naer_all), &
naer2(pcols,pver,naer_all), &
maerosol(1,naer_all), &
ntype(naer_all))
! assign variable deltat for sub-stepping...
deltat=deltatin
! parameters for scheme
omsm=0.99999_r8
dto2=0.5_r8*deltat
mincld=0.0001_r8
! initialize multi-level fields
q(1:ncol,1:pver)=qn(1:ncol,1:pver)
t(1:ncol,1:pver)=tn(1:ncol,1:pver)
! More refined computation of sub-grid vertical velocity
! Set to be zero at the surface by initialization.
if( wsubTKE ) then
do i = 1, ncol
do k = 1, pver
tkem(i,k) = 0.5_r8 * ( tke(i,k) + tke(i,k+1) )
smm(i,k) = 0.5_r8 * ( smaw(i,k) + smaw(i,k+1) )
if( turbtype(i,k) .eq. 3._r8 ) then ! Bottom entrainment interface
tkem(i,k) = 0.5_r8 * tke(i,k+1)
smm(i,k) = 0.5_r8 * smaw(i,k+1)
elseif( turbtype(i,k+1) .eq. 4._r8 ) then ! Top entrainment interface
tkem(i,k) = 0.5_r8 * tke(i,k)
smm(i,k) = 0.5_r8 * smaw(i,k)
endif
smm(i,k) = 0.259_r8*smm(i,k)
smm(i,k) = max(smm(i,k), 0.4743_r8)
end do
end do
endif
do i=1,ncol
ftrue=0
do k=1,pver
! get sub-grid vertical velocity from diff coef.
! following morrison et al. 2005, JAS
! assume mixing length of 30 m
dum=(kkvh(i,k)+kkvh(i,k+1))/2._r8/30._r8
! use maximum sub-grid vertical vel of 10 m/s
dum=min(dum,10._r8)
! set wsub to value at current vertical level
wsub(i,k)=dum
! new computation
if( wsubTKE ) then
wsub(i,k) = sqrt(0.5_r8*(tke(i,k)+tke(i,k+1))*(2._r8/3._r8))
wsub(i,k) = min(wsub(i,k),10._r8)
endif
wsubi(i,k) = max(0.001_r8,wsub(i,k))
wsubi(i,k) = min(wsubi(i,k),0.2_r8)
wsub(i,k) = max(0.20_r8,wsub(i,k))
end do
end do
! initialize time-varying parameters
do k=1,pver
do i=1,ncol
rho(i,k)=p(i,k)/(r*t(i,k))
dv(i,k) = 8.794E-5_r8*t(i,k)**1.81_r8/p(i,k)
mu(i,k) = 1.496E-6_r8*t(i,k)**1.5_r8/ &
(t(i,k)+120._r8)
sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
kap(i,k) = 1.414e3_r8*1.496e-6_r8*t(i,k)** &
1.5_r8/(t(i,k)+120._r8)
! air density adjustment for fallspeed parameters
! includes air density correction factor to the
! power of 0.54 following Heymsfield and Bansemer 2007
rhof(i,k)=(rhosu/rho(i,k))**0.54
arn(i,k)=ar*rhof(i,k)
asn(i,k)=as*rhof(i,k)
acn(i,k)=ac*rhof(i,k)
ain(i,k)=ai*rhof(i,k)
! get dz from dp and hydrostatic approx
! keep dz positive (define as layer k-1 - layer k)
dz(i,k)= pdel(i,k)/(rho(i,k)*g)
end do
end do
! initialization
t1(1:ncol,1:pver) = t(1:ncol,1:pver)
q1(1:ncol,1:pver) = q(1:ncol,1:pver)
qc1(1:ncol,1:pver) = qc(1:ncol,1:pver)
qi1(1:ncol,1:pver) = qi(1:ncol,1:pver)
nc1(1:ncol,1:pver) = nc(1:ncol,1:pver)
ni1(1:ncol,1:pver) = ni(1:ncol,1:pver)
! initialize tendencies to zero
tlat1(1:ncol,1:pver)=0._r8
qvlat1(1:ncol,1:pver)=0._r8
qctend1(1:ncol,1:pver)=0._r8
qitend1(1:ncol,1:pver)=0._r8
nctend1(1:ncol,1:pver)=0._r8
nitend1(1:ncol,1:pver)=0._r8
! initialize precip output
qrout(1:ncol,1:pver)=0._r8
qsout(1:ncol,1:pver)=0._r8
nrout(1:ncol,1:pver)=0._r8
nsout(1:ncol,1:pver)=0._r8
dsout(1:ncol,1:pver)=0._r8
! initialize variables for trop_mozart
nevapr(1:ncol,1:pver) = 0._r8
evapsnow(1:ncol,1:pver) = 0._r8
prain(1:ncol,1:pver) = 0._r8
prodsnow(1:ncol,1:pver) = 0._r8
cmeout(1:ncol,1:pver) = 0._r8
! for refl calc
rainrt1(1:ncol,1:pver) = 0._r8
! initialize precip fraction and output tendencies
cldmax(1:ncol,1:pver)=mincld
!initialize aerosol number
naer2(1:ncol,1:pver,:)=0._r8
dum2l(1:ncol,1:pver)=0._r8
dum2i(1:ncol,1:pver)=0._r8
! initialize avg precip rate
prect1(1:ncol)=0._r8
preci1(1:ncol)=0._r8
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!Get humidity and saturation vapor pressures
! hm: tsp, qsp no longer used below
! call findsp1 (lchnk,ncol,qn,tn,p,tsp,qsphy)
! call findsp1_water (lchnk,ncol,qn,tn,p,tsp,qsp)
do k=1,pver
! find wet bulk temperature and saturation value for provisional t and q without
! condensation
call vqsatd_water(t(1,k),p(1,k),es,qs,gammas,ncol) ! use rhw
do i=1,ncol
esl(i,k)=polysvp(t(i,k),0)
esi(i,k)=polysvp(t(i,k),1)
! hm fix, make sure when above freezing that esi=esl, not active yet
if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
relhum(i,k)=q(i,k)/qs(i)
! get cloud fraction, check for minimum
cldm(i,k)=max(cldn(i,k),mincld)
cldmw(i,k)=max(cldn(i,k),mincld)
icldm(i,k)=max(icecldf(i,k),mincld)
lcldm(i,k)=max(liqcldf(i,k),mincld)
! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)
nfice(i,k)=0._r8
dumfice=qc(i,k)+qi(i,k)
if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then
nfice(i,k)=qi(i,k)/dumfice
endif
if (t(i,k).lt.tmelt - 5._r8) then
if (liu_in) then
#ifndef MODAL_AERO
do m=1,naer_all
ntype(m)=1
maerosol(1,m)=aer_mmr(i,k,m)*rho(i,k)
! For bulk aerosols only:
! set number nucleated for sulfate based on Lohmann et al 2000 (JGR) eq 2
! Na=340.*(massSO4)^0.58 where Na=cm-3 and massSO4=ug/m3
! convert units to Na [m-3] and SO4 [kgm-3]
! Na(m-3)= 1.e6 cm3 m-3 Na(cm-3)=340. *(massSO4[kg/m3]*1.e9ug/kg)^0.58
! or Na(m-3)= 1.e6* 340.*(1.e9ug/kg)^0.58 * (massSO4[kg/m3])^0.58
if(m .eq. idxsul) then
naer2(i,k,m)= 5.64259e13_r8 * maerosol(1,m)**0.58
else
naer2(i,k,m)=maerosol(1,m)*num_to_mass_aer(m)
endif
enddo
#endif
! get ice nucleation rate
call nucleati(wsubi(i,k),t(i,k),relhum(i,k),icldm(i,k),qc(i,k),nfice(i,k),rho(i,k), &
#ifdef MODAL_AERO
qaer(i,k,:)*rho(i,k),dgnum(i,k,:),1,naer_all,dum2,nihf2,niimm2,nidep2,nimey2)
#else
naer2(i,k,:),1,naer_all,dum2,nihf2,niimm2,nidep2,nimey2)
#endif
else
! cooper curve (factor of 1000 is to convert from L-1 to m-3)
dum2=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
dum2=min(dum,208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1
endif
dumnnuc=(dum2-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
dumnnuc=max(dumnnuc,0._r8)
! get provisional ni and qi after nucleation in order to calculate
! Bergeron process below
ninew=ni(i,k)+dumnnuc*deltat
qinew=qi(i,k)+dumnnuc*deltat*mi0
!T>268
else
ninew=ni(i,k)
qinew=qi(i,k)
end if
! Initialize CME components
cme(i,k) = 0._r8
cmei(i,k)=0._r8
!-------------------------------------------------------------------
!Bergeron process
! make sure to initialize bergeron process to zero
berg(i,k)=0._r8
prd = 0._r8
!condensation loop.
! get in-cloud qi and ni after nucleation
if (icldm(i,k) .gt. 0._r8) then
qiic(i,k)=qinew/icldm(i,k)
niic(i,k)=ninew/icldm(i,k)
else
qiic(i,k)=0._r8
niic(i,k)=0._r8
endif
!if T < 0 C then bergeron.
if (t(i,k).lt.273.15) then
!if ice exists
if (qi(i,k).gt.qsmall) then
bergtsf = 0._r8 ! bergeron time scale (fraction of timestep)
qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
dqsidt = xxls*qvi/(rv*t(i,k)**2)
abi = 1._r8+dqsidt*xxls/cpp
! get ice size distribution parameters
if (qiic(i,k).ge.qsmall) then
lami(k) = (cons1*ci* &
niic(i,k)/qiic(i,k))**(1._r8/di)
n0i(k) = niic(i,k)*lami(k)
! check for slope
! adjust vars
if (lami(k).lt.lammini) then
lami(k) = lammini
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
else if (lami(k).gt.lammaxi) then
lami(k) = lammaxi
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
end if
epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k) &
/(lami(k)*lami(k))
!if liquid exists
if (qc(i,k).gt. qsmall) then
!begin bergeron process
! do bergeron (vapor deposition with RHw=1)
! code to find berg (a rate) goes here
! calculate Bergeron process
prd = epsi*(qvl-qvi)/abi
else
prd = 0._r8
end if
! multiply by cloud fraction
prd = prd*min(icldm(i,k),lcldm(i,k))
! transfer of existing cloud liquid to ice
berg(i,k)=max(0._r8,prd)
end if !end liquid exists bergeron
if (berg(i,k).gt.0._r8) then
bergtsf=max(0._r8,(qc(i,k)/berg(i,k))/deltat)
if(bergtsf.lt.1._r8) berg(i,k) = max(0._r8,qc(i,k)/deltat)
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (bergtsf.lt.1._r8.or.icldm(i,k).gt.lcldm(i,k)) then
if (qiic(i,k).ge.qsmall) then
! first case is for case when liquid water is present, but is completely depleted in time step, i.e., bergrsf > 0 but < 1
if (qc(i,k).ge.qsmall) then
rhin = (1.0_r8 + relhum(i,k)) / 2._r8
if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
prd = epsi*(rhin*qvl-qvi)/abi
! multiply by cloud fraction assuming liquid/ice maximum overlap
prd = prd*min(icldm(i,k),lcldm(i,k))
! add to cmei
cmei(i,k) = cmei(i,k) + (prd * (1._r8- bergtsf))
end if ! rhin
end if ! qc > qsmall
! second case is for pure ice cloud, either no liquid, or icldm > lcldm
if (qc(i,k).lt.qsmall.or.icldm(i,k).gt.lcldm(i,k)) then
! note: for case of no liquid, need to set liquid cloud fraction to zero
! store liquid cloud fraction in 'dum'
if (qc(i,k).lt.qsmall) then
dum=0._r8
else
dum=lcldm(i,k)
end if
! set RH to grid-mean value for pure ice cloud
rhin = relhum(i,k)
if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
prd = epsi*(rhin*qvl-qvi)/abi
! multiply by relevant cloud fraction for pure ice cloud
! assuming maximum overlap of liquid/ice
prd = prd*max((icldm(i,k)-dum),0._r8)
cmei(i,k) = cmei(i,k) + prd
end if ! rhin
end if ! qc or icldm > lcldm
end if ! qiic
end if ! bergtsf or icldm > lcldm
! if deposition, it should not reduce grid mean rhi below 1.0
if(cmei(i,k) > 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) > 1._r8 ) &
cmei(i,k)=min(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat)
end if !end ice exists loop
!this ends temperature < 0. loop
!-------------------------------------------------------------------
end if !
!..............................................................
! evaporation should not exceed available water
if ((-berg(i,k)).lt.-qc(i,k)/deltat) &
berg(i,k) = max(qc(i,k)/deltat,0._r8)
!sublimation process...
if ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall ) then
qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
dqsidt = xxls*qvi/(rv*t(i,k)**2)
abi = 1._r8+dqsidt*xxls/cpp
! get ice size distribution parameters
lami(k) = (cons1*ci* &
niic(i,k)/qiic(i,k))**(1._r8/di)
n0i(k) = niic(i,k)*lami(k)
! check for slope
! adjust vars
if (lami(k).lt.lammini) then
lami(k) = lammini
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
else if (lami(k).gt.lammaxi) then
lami(k) = lammaxi
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
end if
epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k) &
/(lami(k)*lami(k))
! modify for ice fraction below
prd = epsi*(relhum(i,k)*qvl-qvi)/abi * icldm(i,k)
cmei(i,k)=min(prd,0._r8)
endif
! sublimation should not exceed available ice
if (cmei(i,k).lt.-qi(i,k)/deltat) &
cmei(i,k)=-qi(i,k)/deltat
! sublimation should not increase grid mean rhi above 1.0
if(cmei(i,k) < 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) < 1._r8 ) &
cmei(i,k)=min(0._r8,max(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat))
! limit cmei due for roundoff error
cmei(i,k)=cmei(i,k)*omsm
! conditional for ice nucleation
if (t(i,k).lt.(tmelt - 5._r8)) then
if ( liu_in ) then
! using Liu et al. (2007) ice nucleation with hooks into simulated aerosol
! ice nucleation rate (dum2) has already been calculated above and state unchanged since then.
dum2i(i,k)=dum2
nihf(i,k)=nihf2
niimm(i,k)=niimm2
nidep(i,k)=nidep2
nimey(i,k)=nimey2
else
! cooper curve (factor of 1000 is to convert from L-1 to m-3)
dum2i(i,k)=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
dum2i(i,k)=min(dum2i(i,k),208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1
endif
else
dum2i(i,k)=0._r8
end if
#ifndef MODAL_AERO
! more general treatment with hooks into simulated aerosol
!
do m=1,naer_all
ntype(m)=1
maerosol(1,m)=aer_mmr(i,k,m)*rho(i,k)
! set number concentration for sulfate based on Lohmann et al 2000 (JGR) eq 2
! Na=340.*(massSO4)^0.58 where Na=cm-3 and massSO4=ug/m3
! convert units to Na [m-3] and SO4 [kgm-3]
! Na(m-3)= 1.e6 cm3 m-3 Na(cm-3)=340. *(massSO4[kg/m3]*1.e9ug/kg)^0.58
! or Na(m-3)= 1.e6* 340.*(1.e9ug/kg)^0.58 * (massSO4[kg/m3])^0.58
!
if(m .eq. idxSUL) then
naer2(i,k,m)= 5.64259e13_r8 * maerosol(1,m)**0.58
else
naer2(i,k,m)=maerosol(1,m)*num_to_mass_aer(m)
endif
enddo
if (qc(i,k).ge.qsmall) then
! get droplet activation rate
call activate(wsub(i,k),t(i,k),rho(i,k), &
naer2(i,k,:), naer_all,naer_all, maerosol, &
dispersion_aer,hygro_aer, density_aer, dum2)
dum2l(i,k) = dum2
else
dum2l(i,k) = 0._r8
end if
#endif
end do ! i loop
end do ! k loop
#ifdef MODAL_AERO
! partition cloud fraction into liquid water part
do k=1,pver
do i=1,ncol
qcld=qc(i,k)+qi(i,k)
if(qcld.gt.qsmall)then
lcldn(i,k)=cldn(i,k)*qc(i,k)/qcld
lcldo(i,k)=cldo(i,k)*qc(i,k)/qcld
else
lcldn(i,k)=0._r8
lcldo(i,k)=0._r8
endif
enddo
enddo
call dropmixnuc(lchnk, ncol, nc, nctend_mixnuc, t, omega, &
p, pint, pdel, rpdel, zm, kkvh, wsub, lcldn, lcldo, &
qaer, cflx, qaertend, qqcw, &
deltat)
cldo(:ncol,:)=cldn(:ncol,:)
#endif
!! initialize sub-step precip flux variables
do i=1,ncol
!! flux is zero at top interface, so these should stay as 0.
rflx1(i,1)=0._r8
sflx1(i,1)=0._r8
do k=1,pver
! initialize normal and sub-step precip flux variables
rflx1(i,k+1)=0._r8
sflx1(i,k+1)=0._r8
end do ! i loop
end do ! k loop
!! initialize final precip flux variables.
do i=1,ncol
!! flux is zero at top interface, so these should stay as 0.
rflx(i,1)=0._r8
sflx(i,1)=0._r8
do k=1,pver
! initialize normal and sub-step precip flux variables
rflx(i,k+1)=0._r8
sflx(i,k+1)=0._r8
end do ! i loop
end do ! k loop
do i=1,ncol
ltrue(i)=0
do k=1,pver
! skip microphysical calculations if no cloud water
if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall.or.cmei(i,k).ge.qsmall) ltrue(i)=1
end do
end do
! assign number of sub-steps to iter
! use 2 sub-steps, following tests described in MG2008
iter = 2
! get sub-step time step
deltat=deltat/real(iter)
! since activation/nucleation processes are fast, need to take into account
! factor mtime = mixing timescale in cloud / model time step
! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity
! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk
#ifdef MODAL_AERO
mtime=1._r8
rate1ord_cw2pr_st(:,:)=0._r8 ! rce 2010/05/01
#else
mtime=deltat/1200._r8
#endif
!!!! skip calculations if no cloud water
do i=1,ncol
if (ltrue(i).eq.0) then
tlat(i,1:pver)=0._r8
qvlat(i,1:pver)=0._r8
qctend(i,1:pver)=0._r8
qitend(i,1:pver)=0._r8
qnitend(i,1:pver)=0._r8
qrtend(i,1:pver)=0._r8
nctend(i,1:pver)=0._r8
nitend(i,1:pver)=0._r8
nrtend(i,1:pver)=0._r8
nstend(i,1:pver)=0._r8
prect(i)=0._r8
preci(i)=0._r8
qniic(i,1:pver)=0._r8
qric(i,1:pver)=0._r8
nsic(i,1:pver)=0._r8
nric(i,1:pver)=0._r8
rainrt(i,1:pver)=0._r8
goto 300
end if
#ifdef MODAL_AERO
qcsinksum_rate1ord(1:pver)=0._r8 ! rce 2010/05/01
qcsum_rate1ord(1:pver)=0._r8 ! rce 2010/05/01
#endif
!!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!.....................................................................................................
do it=1,iter
! initialize sub-step microphysical tendencies
tlat(i,1:pver)=0._r8
qvlat(i,1:pver)=0._r8
qctend(i,1:pver)=0._r8
qitend(i,1:pver)=0._r8
qnitend(i,1:pver)=0._r8
qrtend(i,1:pver)=0._r8
nctend(i,1:pver)=0._r8
nitend(i,1:pver)=0._r8
nrtend(i,1:pver)=0._r8
nstend(i,1:pver)=0._r8
! initialize diagnostic precipitation to zero
qniic(i,1:pver)=0._r8
qric(i,1:pver)=0._r8
nsic(i,1:pver)=0._r8
nric(i,1:pver)=0._r8
rainrt(i,1:pver)=0._r8
! begin new i,k loop, calculate new cldmax after adjustment to cldm above
! initialize vertically-integrated rain and snow tendencies
qrtot = 0._r8
nrtot = 0._r8
qstot = 0._r8
nstot = 0._r8
! initialize precip at surface
prect(i)=0._r8
preci(i)=0._r8
do k=1,pver
! set cwml and cwmi to current qc and qi
cwml(i,k)=qc(i,k)
cwmi(i,k)=qi(i,k)
! initialize precip fallspeeds to zero
ums(k)=0._r8
uns(k)=0._r8
umr(k)=0._r8
unr(k)=0._r8
! calculate precip fraction based on maximum overlap assumption
if (k.eq.1) then
cldmax(i,k)=cldm(i,k)
else
! if rain or snow mix ratio is smaller than
! threshold, then set cldmax to cloud fraction at current level
if (qric(i,k-1).ge.qsmall.or.qniic(i,k-1).ge.qsmall) then
cldmax(i,k)=max(cldmax(i,k-1),cldm(i,k))
else
cldmax(i,k)=cldm(i,k)
end if
end if
! decrease in number concentration due to sublimation/evap
! divide by cloud fraction to get in-cloud decrease
! don't reduce Nc due to bergeron process
if (cmei(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. cldm(i,k) > mincld) then
nsubi(k)=cmei(i,k)/qi(i,k)*ni(i,k)/cldm(i,k)
else
nsubi(k)=0._r8
end if
nsubc(k)=0._r8
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
if (dum2i(i,k).gt.0._r8.and.t(i,k).lt.(tmelt - 5._r8).and. &
relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then
!if NCAI > 0. then set numice = ncai (as before)
!note: this is gridbox averaged
nnuccd(k)=(dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
nnuccd(k)=max(nnuccd(k),0._r8)
nimax = dum2i(i,k)*icldm(i,k)
!Calc mass of new particles using new crystal mass...
!also this will be multiplied by mtime as nnuccd is...
mnuccd(k) = nnuccd(k) * mi0
! add mnuccd to cmei....
cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime
! limit cmei
qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
dqsidt = xxls*qvi/(rv*t(i,k)**2)
abi = 1._r8+dqsidt*xxls/cpp
cmei(i,k)=min(cmei(i,k),(q(i,k)-qvi)/abi/deltat)
! limit for roundoff error
cmei(i,k)=cmei(i,k)*omsm
else
nnuccd(k)=0._r8
nimax = 0._r8
mnuccd(k) = 0._r8
end if
!c............................................................................
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
! for microphysical process calculations
! units are kg/kg for mixing ratio, 1/kg for number conc
! limit in-cloud values to 0.005 kg/kg
qcic(i,k)=min(cwml(i,k)/lcldm(i,k),5.e-3_r8)
qiic(i,k)=min(cwmi(i,k)/icldm(i,k),5.e-3_r8)
ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)
if (qc(i,k) - berg(i,k)*deltat.lt.qsmall) then
qcic(i,k)=0._r8
ncic(i,k)=0._r8
if (qc(i,k)-berg(i,k)*deltat.lt.0._r8) then
berg(i,k)=qc(i,k)/deltat*omsm
end if
end if
if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.qsmall) then
qiic(i,k)=0._r8
niic(i,k)=0._r8
if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.0._r8) then
cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm
end if
end if
! add to cme output
cmeout(i,k) = cmeout(i,k)+cmei(i,k)
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! droplet activation
! calculate potential for droplet activation if cloud water is present
! formulation from Abdul-Razzak and Ghan (2000) and Abdul-Razzak et al. (1998), AR98
! assume aerosols already activated are equal to number of existing droplets for simplicity
! multiply by cloud fraction to obtain grid-average tendency
#ifdef MODAL_AERO
if (qcic(i,k).ge.qsmall) then
npccn(k) = nctend_mixnuc(i,k)
npccn(k) = max(0._r8,npccn(k))
dum2l(i,k)=(nc(i,k)+npccn(k)*deltat)/cldm(i,k)
dum2l(i,k)=max(dum2l(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
ncmax = dum2l(i,k)*cldm(i,k)
else
npccn(k)=0._r8
dum2l(i,k)=0._r8
ncmax = 0._r8
end if
#else
if (qcic(i,k).ge.qsmall) then
npccn(k) = (dum2l(i,k)-nc(i,k)/lcldm(i,k))/deltat*lcldm(i,k)
! make sure number activated > 0
npccn(k) = max(0._r8,npccn(k))
ncmax = dum2l(i,k)*lcldm(i,k)
else
npccn(k)=0._r8
ncmax = 0._r8
end if
#endif
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! get size distribution parameters based on in-cloud cloud water/ice
! these calculations also ensure consistency between number and mixing ratio
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!......................................................................
! cloud ice
if (qiic(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
niic(i,k)=min(niic(i,k),qiic(i,k)*1.e20_r8)
lami(k) = (cons1*ci* &
niic(i,k)/qiic(i,k))**(1._r8/di)
n0i(k) = niic(i,k)*lami(k)
! check for slope
! adjust vars
if (lami(k).lt.lammini) then
lami(k) = lammini
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
niic(i,k) = n0i(k)/lami(k)
else if (lami(k).gt.lammaxi) then
lami(k) = lammaxi
n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
niic(i,k) = n0i(k)/lami(k)
end if
else
lami(k) = 0._r8
n0i(k) = 0._r8
end if
if (qcic(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
ncic(i,k)=min(ncic(i,k),qcic(i,k)*1.e20_r8)
ncic(i,k)=max(ncic(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm
! get pgam from fit to observations of martin et al. 1994
! pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8/rho(i,k))+0.2714_r8
! hm fix, not active yet
pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
pgam(k)=1._r8/(pgam(k)**2)-1._r8
pgam(k)=max(pgam(k),2._r8)
pgam(k)=min(pgam(k),15._r8)
! calculate lamc
lamc(k) = (pi/6._r8*rhow*ncic(i,k)*gamma(pgam(k)+4._r8)/ &
(qcic(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
! lammin, 50 micron diameter max mean size
lammin = (pgam(k)+1._r8)/50.e-6_r8
lammax = (pgam(k)+1._r8)/2.e-6_r8
if (lamc(k).lt.lammin) then
lamc(k) = lammin
ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
gamma(pgam(k)+1._r8)/ &
(pi*rhow*gamma(pgam(k)+4._r8))
else if (lamc(k).gt.lammax) then
lamc(k) = lammax
ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
gamma(pgam(k)+1._r8)/ &
(pi*rhow*gamma(pgam(k)+4._r8))
end if
! parameter to calculate droplet freezing
cdist1(k) = ncic(i,k)/gamma(pgam(k)+1._r8)
else
lamc(k) = 0._r8
cdist1(k) = 0._r8
end if
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! begin micropysical process calculations
!.................................................................
! autoconversion of cloud liquid water to rain
! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
! minimum qc of 1 x 10^-8 prevents floating point error
if (qcic(i,k).ge.1.e-8_r8) then
! nprc is increase in rain number conc due to autoconversion
! nprc1 is decrease in cloud droplet conc due to autoconversion
! assume exponential sub-grid distribution of qc, resulting in additional
! factor related to qcvar below
prc(k) = cons2/(cons3*cons18)*1350._r8*qcic(i,k)**2.47_r8* &
(ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
! nprc(k) = prc(k)/(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
nprc(k) = prc(k)/cons22
nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
else
prc(k)=0._r8
nprc(k)=0._r8
nprc1(k)=0._r8
end if
! add autoconversion to precip from above to get provisional rain mixing ratio
! and number concentration (qric and nric)
! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)
dum=0.45_r8
dum1=0.45_r8
if (k.eq.1) then
qric(i,k)=prc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
nric(i,k)=nprc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
else
if (qric(i,k-1).ge.qsmall) then
dum=umr(k-1)
dum1=unr(k-1)
end if
! no autoconversion of rain number if rain/snow falling from above
! this assumes that new drizzle drops formed by autoconversion are rapidly collected
! by the existing rain/snow particles from above
if (qric(i,k-1).ge.1.e-9_r8.or.qniic(i,k-1).ge.1.e-9_r8) then
nprc(k)=0._r8
end if
qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*((pra(k-1)+prc(k))*lcldm(i,k)+(pre(k-1)-pracs(k-1)-mnuccr(k-1))*cldmax(i,k))))&
/(dum*rho(i,k)*cldmax(i,k))
nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*(nprc(k)*lcldm(i,k)+(nsubr(k-1)-npracs(k-1)-nnuccr(k-1)+nragg(k-1))*cldmax(i,k))))&
/(dum1*rho(i,k)*cldmax(i,k))
end if
!.......................................................................
! Autoconversion of cloud ice to snow
! similar to Ferrier (1994)
if (t(i,k).le.273.15_r8.and.qiic(i,k).ge.qsmall) then
! note: assumes autoconversion timescale of 180 sec
nprci(k) = n0i(k)/(lami(k)*180._r8)*exp(-lami(k)*dcs)
prci(k) = pi*rhoi*n0i(k)/(6._r8*180._r8)* &
(cons23/lami(k)+3._r8*cons24/lami(k)**2+ &
6._r8*dcs/lami(k)**3+6._r8/lami(k)**4)*exp(-lami(k)*dcs)
else
prci(k)=0._r8
nprci(k)=0._r8
end if
! add autoconversion to flux from level above to get provisional snow mixing ratio
! and number concentration (qniic and nsic)
dum=(asn(i,k)*cons25)
dum1=(asn(i,k)*cons25)
if (k.eq.1) then
qniic(i,k)=prci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
nsic(i,k)=nprci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
else
if (qniic(i,k-1).ge.qsmall) then
dum=ums(k-1)
dum1=uns(k-1)
end if
qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*((prci(k)+prai(k-1)+psacws(k-1)+bergs(k-1))*icldm(i,k)+(prds(k-1)+ &
pracs(k-1)+mnuccr(k-1))*cldmax(i,k))))&
/(dum*rho(i,k)*cldmax(i,k))
nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+nnuccr(k-1))*cldmax(i,k))))&
/(dum1*rho(i,k)*cldmax(i,k))
end if
! if precip mix ratio is zero so should number concentration
if (qniic(i,k).lt.qsmall) then
qniic(i,k)=0._r8
nsic(i,k)=0._r8
end if
if (qric(i,k).lt.qsmall) then
qric(i,k)=0._r8
nric(i,k)=0._r8
end if
! make sure number concentration is a positive number to avoid
! taking root of negative later
nric(i,k)=max(nric(i,k),0._r8)
nsic(i,k)=max(nsic(i,k),0._r8)
!.......................................................................
! get size distribution parameters for precip
!......................................................................
! rain
if (qric(i,k).ge.qsmall) then
lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
n0r(k) = nric(i,k)*lamr(k)
! check for slope
! adjust vars
if (lamr(k).lt.lamminr) then
lamr(k) = lamminr
n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
nric(i,k) = n0r(k)/lamr(k)
else if (lamr(k).gt.lammaxr) then
lamr(k) = lammaxr
n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
nric(i,k) = n0r(k)/lamr(k)
end if
! provisional rain number and mass weighted mean fallspeed (m/s)
unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
else
lamr(k) = 0._r8
n0r(k) = 0._r8
umr(k) = 0._r8
unr(k) = 0._r8
end if
!......................................................................
! snow
if (qniic(i,k).ge.qsmall) then
lams(k) = (cons6*cs*nsic(i,k)/ &
qniic(i,k))**(1._r8/ds)
n0s(k) = nsic(i,k)*lams(k)
! check for slope
! adjust vars
if (lams(k).lt.lammins) then
lams(k) = lammins
n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
nsic(i,k) = n0s(k)/lams(k)
else if (lams(k).gt.lammaxs) then
lams(k) = lammaxs
n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
nsic(i,k) = n0s(k)/lams(k)
end if
! provisional snow number and mass weighted mean fallspeed (m/s)
ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
else
lams(k) = 0._r8
n0s(k) = 0._r8
ums(k) = 0._r8
uns(k) = 0._r8
end if
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! heterogeneous freezing of cloud water
if (qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then
! immersion freezing (Bigg, 1953)
mnuccc(k) = cons9/(cons3*cons19)* &
pi*pi/36._r8*rhow* &
cdist1(k)*gamma(7._r8+pgam(k))* &
bimm*exp(aimm*(273.15_r8-t(i,k)))/ &
lamc(k)**3/lamc(k)**3
nnuccc(k) = cons10/(cons3*qcvar)* &
pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
*bimm* &
exp(aimm*(273.15_r8-t(i,k)))/lamc(k)**3
! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
tcnt=(270.16_r8-t(i,k))**1.3_r8
viscosity=1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8 ! Viscosity (kg/m/s)
mfp=2.0_r8*viscosity/(p(i,k) & ! Mean free path (m)
*sqrt(8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i,k))))
slip1=1.0_r8+(mfp/rn_dst1)*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rn_dst1/mfp))))! Slip correction factor
slip2=1.0_r8+(mfp/rn_dst2)*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rn_dst2/mfp))))
slip3=1.0_r8+(mfp/rn_dst3)*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rn_dst3/mfp))))
slip4=1.0_r8+(mfp/rn_dst4)*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rn_dst4/mfp))))
dfaer1=1.381e-23_r8*t(i,k)*slip1/(6._r8*pi*viscosity*rn_dst1) ! aerosol diffusivity (m2/s)
dfaer2=1.381e-23_r8*t(i,k)*slip2/(6._r8*pi*viscosity*rn_dst2)
dfaer3=1.381e-23_r8*t(i,k)*slip3/(6._r8*pi*viscosity*rn_dst3)
dfaer4=1.381e-23_r8*t(i,k)*slip4/(6._r8*pi*viscosity*rn_dst4)
nacon1=0.0_r8
nacon2=0.0_r8
nacon3=0.0_r8
nacon4=0.0_r8
!use size '3' for dust coarse mode...
!scale by dust fraction in coarse mode
#ifdef MODAL_AERO
dmc= qaer(i,k,lptr_dust_a_amode(modeptr_coarse))
ssmc=qaer(i,k,lptr_nacl_a_amode(modeptr_coarse))
if (dmc.gt.0.0_r8) then
nacon3=dmc/(ssmc+dmc) * qaer(i,k,numptr_amode(modeptr_coarse))*rho(i,k)*tcnt
else
nacon3=0._r8
endif
!also redefine parameters based on size...
dstrn=0.5_r8*dgnumwet(i,k,modeptr_coarse)
if (dstrn.le.0._r8) then
dstrn=rn_dst3
endif
slip3=1.0_r8+(mfp/dstrn)*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*dstrn/mfp))))
dfaer3=1.381e-23_r8*t(i,k)*slip3/(6._r8*pi*viscosity*dstrn)
#else
if (idxdst1.gt.0) then
nacon1=naer2(i,k,idxdst1)*tcnt *0.0_r8
endif
if (idxdst2.gt.0) then
nacon2=naer2(i,k,idxdst2)*tcnt ! 1/m3
endif
if (idxdst3.gt.0) then
nacon3=naer2(i,k,idxdst3)*tcnt
endif
if (idxdst4.gt.0) then
nacon4=naer2(i,k,idxdst4)*tcnt
endif
#endif
mnucct(k) = cons10/(cons3*qcvar)* &
(dfaer1*nacon1+dfaer2*nacon2+dfaer3*nacon3+dfaer4*nacon4)*pi*pi/3._r8*rhow* &
cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
nnucct(k) = (dfaer1*nacon1+dfaer2*nacon2+dfaer3*nacon3+dfaer4*nacon4)*2._r8*pi* &
cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
! make sure number of droplets frozen does not exceed available ice nuclei concentration
! this prevents 'runaway' droplet freezing
if (nnuccc(k)*lcldm(i,k).gt.nnuccd(k)) then
dum=(nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
! scale mixing ratio of droplet freezing with limit
mnuccc(k)=mnuccc(k)*dum
nnuccc(k)=nnuccd(k)/lcldm(i,k)
end if
else
mnuccc(k)=0._r8
nnuccc(k)=0._r8
mnucct(k)=0._r8
nnucct(k)=0._r8
end if
!.......................................................................
! snow self-aggregation from passarelli, 1978, used by reisner, 1998
! this is hard-wired for bs = 0.4 for now
! ignore self-collection of cloud ice
if (qniic(i,k).ge.qsmall .and. t(i,k).le.273.15_r8) then
nsagg(k) = -1108._r8*asn(i,k)*Eii* &
pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)*rho(i,k)** &
((2._r8+bs)/3._r8)*qniic(i,k)**((2._r8+bs)/3._r8)* &
(nsic(i,k)*rho(i,k))**((4._r8-bs)/3._r8)/ &
(4._r8*720._r8*rho(i,k))
else
nsagg(k)=0._r8
end if
!.......................................................................
! accretion of cloud droplets onto snow/graupel
! here use continuous collection equation with
! simple gravitational collection kernel
! ignore collisions between droplets/cloud ice
! since minimum size ice particle for accretion is 50 - 150 micron
! ignore collision of snow with droplets above freezing
if (qniic(i,k).ge.qsmall .and. t(i,k).le.tmelt .and. &
qcic(i,k).ge.qsmall) then
! put in size dependent collection efficiency
! mean diameter of snow is area-weighted, since
! accretion is function of crystal geometric area
! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
dc0 = (pgam(k)+1._r8)/lamc(k)
ds0 = 1._r8/lams(k)
dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0)
eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
eci = max(eci,0._r8)
eci = min(eci,1._r8)
! no impact of sub-grid distribution of qc since psacws
! is linear in qc
psacws(k) = pi/4._r8*asn(i,k)*qcic(i,k)*rho(i,k)* &
n0s(k)*Eci*cons11/ &
lams(k)**(bs+3._r8)
npsacws(k) = pi/4._r8*asn(i,k)*ncic(i,k)*rho(i,k)* &
n0s(k)*Eci*cons11/ &
lams(k)**(bs+3._r8)
else
psacws(k)=0._r8
npsacws(k)=0._r8
end if
! add secondary ice production due to accretion of droplets by snow
! (Hallet-Mossop process) (from Cotton et al., 1986)
if((t(i,k).lt.270.16_r8) .and. (t(i,k).ge.268.16_r8)) then
ni_secp = 3.5e8_r8*(270.16_r8-t(i,k))/2.0_r8*psacws(k)
nsacwi(k) = ni_secp
msacwi(k) = min(ni_secp*mi0,psacws(k))
else if((t(i,k).lt.268.16_r8) .and. (t(i,k).ge.265.16_r8)) then
ni_secp = 3.5e8_r8*(t(i,k)-265.16_r8)/3.0_r8*psacws(k)
nsacwi(k) = ni_secp
msacwi(k) = min(ni_secp*mi0,psacws(k))
else
ni_secp = 0.0_r8
nsacwi(k) = 0.0_r8
msacwi(k) = 0.0_r8
endif
psacws(k) = max(0.0_r8,psacws(k)-ni_secp*mi0)
!.......................................................................
! accretion of rain water by snow
! formula from ikawa and saito, 1991, used by reisner et al., 1998
if (qric(i,k).ge.1.e-8_r8 .and. qniic(i,k).ge.1.e-8_r8 .and. &
t(i,k).le.273.15_r8) then
pracs(k) = pi*pi*ecr*(((1.2_r8*umr(k)-0.95_r8*ums(k))**2+ &
0.08_r8*ums(k)*umr(k))**0.5_r8*rhow*rho(i,k)* &
n0r(k)*n0s(k)* &
(5._r8/(lamr(k)**6*lams(k))+ &
2._r8/(lamr(k)**5*lams(k)**2)+ &
0.5_r8/(lamr(k)**4*lams(k)**3)))
npracs(k) = pi/2._r8*rho(i,k)*ecr*(1.7_r8*(unr(k)-uns(k))**2+ &
0.3_r8*unr(k)*uns(k))**0.5_r8*n0r(k)*n0s(k)* &
(1._r8/(lamr(k)**3*lams(k))+ &
1._r8/(lamr(k)**2*lams(k)**2)+ &
1._r8/(lamr(k)*lams(k)**3))
else
pracs(k)=0._r8
npracs(k)=0._r8
end if
!.......................................................................
! heterogeneous freezing of rain drops
! follows from Bigg (1953)
if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then
mnuccr(k) = 20._r8*pi*pi*rhow*nric(i,k)*bimm* &
exp(aimm*(273.15_r8-t(i,k)))/lamr(k)**3 &
/lamr(k)**3
nnuccr(k) = pi*nric(i,k)*bimm* &
exp(aimm*(273.15_r8-t(i,k)))/lamr(k)**3
else
mnuccr(k)=0._r8
nnuccr(k)=0._r8
end if
!.......................................................................
! accretion of cloud liquid water by rain
! formula from Khrouditnov and Kogan (2000)
! gravitational collection kernel, droplet fall speed neglected
if (qric(i,k).ge.qsmall .and. qcic(i,k).ge.qsmall) then
! include sub-grid distribution of cloud water
pra(k) = cons12/(cons3*cons20)* &
67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
else
pra(k)=0._r8
npra(k)=0._r8
end if
!.......................................................................
! Self-collection of rain drops
! from Beheng(1994)
if (qric(i,k).ge.qsmall) then
nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
else
nragg(k)=0._r8
end if
!.......................................................................
! Accretion of cloud ice by snow
! For this calculation, it is assumed that the Vs >> Vi
! and Ds >> Di for continuous collection
if (qniic(i,k).ge.qsmall.and.qiic(i,k).ge.qsmall &
.and.t(i,k).le.273.15_r8) then
prai(k) = pi/4._r8*asn(i,k)*qiic(i,k)*rho(i,k)* &
n0s(k)*Eii*cons11/ &
lams(k)**(bs+3._r8)
nprai(k) = pi/4._r8*asn(i,k)*niic(i,k)* &
rho(i,k)*n0s(k)*Eii*cons11/ &
lams(k)**(bs+3._r8)
else
prai(k)=0._r8
nprai(k)=0._r8
end if
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calculate evaporation/sublimation of rain and snow
! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
! in-cloud condensation/deposition of rain and snow is neglected
! except for transfer of cloud water to snow through bergeron process
! initialize evap/sub tendncies
pre(k)=0._r8
prds(k)=0._r8
! evaporation of rain
! only calculate if there is some precip fraction > cloud fraction
if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8.or.cldmax(i,k).gt.lcldm(i,k)) then
! set temporary cloud fraction to zero if cloud water + ice is very small
! this will ensure that evaporation/sublimation of precip occurs over
! entire grid cell, since min cloud fraction is specified otherwise
if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8) then
dum=0._r8
else
dum=lcldm(i,k)
end if
! saturation vapor pressure
esn=polysvp(t(i,k),0)
qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
! recalculate saturation vapor pressure for liquid and ice
esl(i,k)=esn
esi(i,k)=polysvp(t(i,k),1)
! hm fix, make sure when above freezing that esi=esl, not active yet
if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
! calculate q for out-of-cloud region
qclr=(q(i,k)-dum*qsn)/(1._r8-dum)
if (qric(i,k).ge.qsmall) then
qvs=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
dqsdt = xxlv*qvs/(rv*t(i,k)**2)
ab = 1._r8+dqsdt*xxlv/cpp
epsr = 2._r8*pi*n0r(k)*rho(i,k)*Dv(i,k)* &
(f1r/(lamr(k)*lamr(k))+ &
f2r*(arn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
sc(i,k)**(1._r8/3._r8)*cons13/ &
(lamr(k)**(5._r8/2._r8+br/2._r8)))
pre(k) = epsr*(qclr-qvs)/ab
! only evaporate in out-of-cloud region
! and distribute across cldmax
pre(k)=min(pre(k)*(cldmax(i,k)-dum),0._r8)
pre(k)=pre(k)/cldmax(i,k)
end if
! sublimation of snow
if (qniic(i,k).ge.qsmall) then
qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
dqsidt = xxls*qvi/(rv*t(i,k)**2)
abi = 1._r8+dqsidt*xxls/cpp
epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
(f1s/(lams(k)*lams(k))+ &
f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
sc(i,k)**(1._r8/3._r8)*cons14/ &
(lams(k)**(5._r8/2._r8+bs/2._r8)))
prds(k) = epss*(qclr-qvi)/abi
! only sublimate in out-of-cloud region and distribute over cldmax
prds(k)=min(prds(k)*(cldmax(i,k)-dum),0._r8)
prds(k)=prds(k)/cldmax(i,k)
end if
! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
! get updated RH at end of time step based on cloud water/ice condensation/evap
qtmp=q(i,k)-(cmei(i,k)+(pre(k)+prds(k))*cldmax(i,k))*deltat
ttmp=t(i,k)+((pre(k)*cldmax(i,k))*xxlv+ &
(cmei(i,k)+prds(k)*cldmax(i,k))*xxls)*deltat/cpp
!limit range of temperatures!
ttmp=max(180._r8,min(ttmp,323._r8))
esn=polysvp(ttmp,0) ! use rhw to allow ice supersaturation
qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
! modify precip evaporation rate if q > qsat
if (qtmp.gt.qsn) then
if (pre(k)+prds(k).lt.-1.e-20) then
dum1=pre(k)/(pre(k)+prds(k))
! recalculate q and t after cloud water cond but without precip evap
qtmp=q(i,k)-(cmei(i,k))*deltat
ttmp=t(i,k)+(cmei(i,k)*xxls)*deltat/cpp
esn=polysvp(ttmp,0) ! use rhw to allow ice supersaturation
qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2))
dum=min(dum,0._r8)
! modify rates if needed, divide by cldmax to get local (in-precip) value
pre(k)=dum*dum1/deltat/cldmax(i,k)
! do separately using RHI for prds....
esn=polysvp(ttmp,1) ! use rhi to allow ice supersaturation
qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2))
dum=min(dum,0._r8)
! modify rates if needed, divide by cldmax to get local (in-precip) value
prds(k)=dum*(1._r8-dum1)/deltat/cldmax(i,k)
end if
end if
end if
! bergeron process - evaporation of droplets and deposition onto snow
if (qniic(i,k).ge.qsmall.and.qcic(i,k).ge.qsmall.and.t(i,k).lt.tmelt) then
qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
qvs=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
dqsidt = xxls*qvi/(rv*t(i,k)**2)
abi = 1._r8+dqsidt*xxls/cpp
epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
(f1s/(lams(k)*lams(k))+ &
f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
sc(i,k)**(1._r8/3._r8)*cons14/ &
(lams(k)**(5._r8/2._r8+bs/2._r8)))
bergs(k)=epss*(qvs-qvi)/abi
else
bergs(k)=0._r8
end if
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! conservation to ensure no negative values of cloud water/precipitation
! in case microphysical process rates are large
! make sure and use end-of-time step values for cloud water, ice, due
! condensation/deposition
! note: for check on conservation, processes are multiplied by omsm
! to prevent problems due to round off error
! include mixing timescale (mtime)
qce=(qc(i,k) - berg(i,k)*deltat)
nce=(nc(i,k)+npccn(k)*deltat*mtime)
qie=(qi(i,k)+(cmei(i,k)+berg(i,k))*deltat)
nie=(ni(i,k)+nnuccd(k)*deltat*mtime)
! conservation of qc
dum = (prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+ &
psacws(k)+bergs(k))*lcldm(i,k)*deltat
if (dum.gt.qce) then
ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm
prc(k) = prc(k)*ratio
pra(k) = pra(k)*ratio
mnuccc(k) = mnuccc(k)*ratio
mnucct(k) = mnucct(k)*ratio
msacwi(k) = msacwi(k)*ratio
psacws(k) = psacws(k)*ratio
bergs(k) = bergs(k)*ratio
end if
! conservation of nc
dum = (nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+ &
npsacws(k)-nsubc(k))*lcldm(i,k)*deltat
if (dum.gt.nce) then
ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+&
npsacws(k)-nsubc(k))*lcldm(i,k))*omsm
nprc1(k) = nprc1(k)*ratio
npra(k) = npra(k)*ratio
nnuccc(k) = nnuccc(k)*ratio
nnucct(k) = nnucct(k)*ratio !++xl
npsacws(k) = npsacws(k)*ratio
nsubc(k)=nsubc(k)*ratio
end if
! conservation of qi
dum = ((-mnuccc(k)-mnucct(k)-msacwi(k))*lcldm(i,k)+(prci(k)+ &
prai(k))*icldm(i,k))*deltat
if (dum.gt.qie) then
ratio = (qie/deltat+(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm
prci(k) = prci(k)*ratio
prai(k) = prai(k)*ratio
end if
! conservation of ni
dum = ((-nnucct(k)-nsacwi(k))*lcldm(i,k)+(nprci(k)+ &
nprai(k)-nsubi(k))*icldm(i,k))*deltat
if (dum.gt.nie) then
ratio = (nie/deltat+(nnucct(k)+nsacwi(k))*lcldm(i,k))/ &
((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
nprci(k) = nprci(k)*ratio
nprai(k) = nprai(k)*ratio
nsubi(k) = nsubi(k)*ratio
end if
! for preciptiation conservation, use logic that vertical integral
! of tendency from current level to top of model (i.e., qrtot) cannot be negative
! conservation of rain mixing rat
if (((prc(k)+pra(k))*lcldm(i,k)+(-mnuccr(k)+pre(k)-pracs(k))*&
cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot.lt.0._r8) then
if (-pre(k)+pracs(k)+mnuccr(k).ge.qsmall) then
ratio = (qrtot/(dz(i,k)*rho(i,k))+(prc(k)+pra(k))*lcldm(i,k))/&
((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm
pre(k) = pre(k)*ratio
pracs(k) = pracs(k)*ratio
mnuccr(k) = mnuccr(k)*ratio
end if
end if
! conservation of nr
! for now neglect evaporation of nr
nsubr(k)=0._r8
if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k)&
+nragg(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot.lt.0._r8) then
if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k).ge.qsmall) then
ratio = (nrtot/(dz(i,k)*rho(i,k))+nprc(k)*lcldm(i,k))/&
((-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k))*cldmax(i,k))*omsm
nsubr(k) = nsubr(k)*ratio
npracs(k) = npracs(k)*ratio
nnuccr(k) = nnuccr(k)*ratio
nragg(k) = nragg(k)*ratio
end if
end if
! conservation of snow mix ratio
if (((bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+(pracs(k)+&
mnuccr(k)+prds(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot.lt.0._r8) then
if (-prds(k).ge.qsmall) then
ratio = (qstot/(dz(i,k)*rho(i,k))+(bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+&
(pracs(k)+mnuccr(k))*cldmax(i,k))/(-prds(k)*cldmax(i,k))*omsm
prds(k) = prds(k)*ratio
end if
end if
! conservation of ns
! calculate loss of number due to sublimation
! for now neglect sublimation of ns
nsubs(k)=0._r8
if ((nprci(k)*icldm(i,k)+(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))*&
dz(i,k)*rho(i,k)+nstot.lt.0._r8) then
if (-nsubs(k)-nsagg(k).ge.qsmall) then
ratio = (nstot/(dz(i,k)*rho(i,k))+nprci(k)*icldm(i,k)+&
nnuccr(k)*cldmax(i,k))/((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm
nsubs(k) = nsubs(k)*ratio
nsagg(k) = nsagg(k)*ratio
end if
end if
! get tendencies due to microphysical conversion processes
! note: tendencies are multiplied by appropaiate cloud/precip
! fraction to get grid-scale values
! note: cmei is already grid-average values
qvlat(i,k) = qvlat(i,k)- &
(pre(k)+prds(k))*cldmax(i,k)-cmei(i,k)
tlat(i,k) = tlat(i,k)+((pre(k)*cldmax(i,k)) &
*xxlv+(prds(k)*cldmax(i,k)+cmei(i,k))*xxls+ &
((bergs(k)+psacws(k)+mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(mnuccr(k)+ &
pracs(k))*cldmax(i,k)+berg(i,k))*xlf)
qctend(i,k) = qctend(i,k)+ &
(-pra(k)-prc(k)-mnuccc(k)-mnucct(k)-msacwi(k)- &
psacws(k)-bergs(k))*lcldm(i,k)-berg(i,k)
qitend(i,k) = qitend(i,k)+ &
(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(-prci(k)- &
prai(k))*icldm(i,k)+cmei(i,k)+berg(i,k)
qrtend(i,k) = qrtend(i,k)+ &
(pra(k)+prc(k))*lcldm(i,k)+(pre(k)-pracs(k)- &
mnuccr(k))*cldmax(i,k)
qnitend(i,k) = qnitend(i,k)+ &
(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(prds(k)+ &
pracs(k)+mnuccr(k))*cldmax(i,k)
! add output for cmei (accumulate)
cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)
! assign variables for trop_mozart, these are grid-average
! evaporation/sublimation is stored here as positive term
evapsnow(i,k) = evapsnow(i,k)-prds(k)*cldmax(i,k)
nevapr(i,k) = nevapr(i,k)-pre(k)*cldmax(i,k)
! change to make sure prain is positive: do not remove snow from
! prain used for wet deposition
prain(i,k) = prain(i,k)+(pra(k)+prc(k))*lcldm(i,k)+(-pracs(k)- &
mnuccr(k))*cldmax(i,k)
prodsnow(i,k) = prodsnow(i,k)+(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(&
pracs(k)+mnuccr(k))*cldmax(i,k)
#ifdef MODAL_AERO
! rce 2010/05/01
! following are used to calculate 1st order conversion rate of cloud water
! to rain and snow (1/s), for later use in aerosol wet removal routine
! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
! used to calculate pra, prc, ... in this routine
! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
! (no cloud ice or bergeron terms)
! qcsum_rate1ord = sum over iterations{ qc used in calculation of the transfer terms }
qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k) + (pra(k)+prc(k)+psacws(k))*lcldm(i,k) ! rce 2010/05/01
qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k) ! rce 2010/05/01
#endif
! microphysics output, note this is grid-averaged
prao(i,k)=prao(i,k)+pra(k)*lcldm(i,k)
prco(i,k)=prco(i,k)+prc(k)*lcldm(i,k)
mnuccco(i,k)=mnuccco(i,k)+mnuccc(k)*lcldm(i,k)
mnuccto(i,k)=mnuccto(i,k)+mnucct(k)*lcldm(i,k)
msacwio(i,k)=msacwio(i,k)+msacwi(k)*lcldm(i,k)
psacwso(i,k)=psacwso(i,k)+psacws(k)*lcldm(i,k)
bergso(i,k)=bergso(i,k)+bergs(k)*lcldm(i,k)
bergo(i,k)=bergo(i,k)+berg(i,k)
prcio(i,k)=prcio(i,k)+prci(k)*icldm(i,k)
praio(i,k)=praio(i,k)+prai(k)*icldm(i,k)
mnuccro(i,k)=mnuccro(i,k)+mnuccr(k)*cldmax(i,k)
pracso (i,k)=pracso (i,k)+pracs (k)*cldmax(i,k)
! multiply activation/nucleation by mtime to account for fast timescale
nctend(i,k) = nctend(i,k)+ npccn(k)*mtime+&
(-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k) &
-npra(k)-nprc1(k))*lcldm(i,k)
nitend(i,k) = nitend(i,k)+ nnuccd(k)*mtime+ &
(nnucct(k)+nsacwi(k))*lcldm(i,k)+(nsubi(k)-nprci(k)- &
nprai(k))*icldm(i,k)
nstend(i,k) = nstend(i,k)+(nsubs(k)+ &
nsagg(k)+nnuccr(k))*cldmax(i,k)+nprci(k)*icldm(i,k)
nrtend(i,k) = nrtend(i,k)+ &
nprc(k)*lcldm(i,k)+(nsubr(k)-npracs(k)-nnuccr(k) &
+nragg(k))*cldmax(i,k)
! make sure that nc and ni at advanced time step do not exceed
! maximum (existing N + source terms*dt), which is possible due to
! fast nucleation timescale
if (nctend(i,k).gt.0._r8.and.nc(i,k)+nctend(i,k)*deltat.gt.ncmax) then
nctend(i,k)=max(0._r8,(ncmax-nc(i,k))/deltat)
end if
if (nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax) then
nitend(i,k)=max(0._r8,(nimax-ni(i,k))/deltat)
end if
! get final values for precipitation q and N, based on
! flux of precip from above, source/sink term, and terminal fallspeed
! see eq. 15-16 in MG2008
! rain
if (qric(i,k).ge.qsmall) then
if (k.eq.1) then
qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(k)
nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(k)
else
qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(k)*rho(i,k)*cldmax(i,k))
nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(k)*rho(i,k)*cldmax(i,k))
end if
else
qric(i,k)=0._r8
nric(i,k)=0._r8
end if
! snow
if (qniic(i,k).ge.qsmall) then
if (k.eq.1) then
qniic(i,k)=qnitend(i,k)*dz(i,k)/cldmax(i,k)/ums(k)
nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(k)
else
qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*qnitend(i,k)))/(ums(k)*rho(i,k)*cldmax(i,k))
nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
(rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(k)*rho(i,k)*cldmax(i,k))
end if
else
qniic(i,k)=0._r8
nsic(i,k)=0._r8
end if
! calculate precipitation flux at surface
! divide by density of water to get units of m/s
prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+&
qnitend(i,k)*dz(i,k)*rho(i,k))/rhow
preci(i) = preci(i)+qnitend(i,k)*dz(i,k)*rho(i,k)/rhow
! convert rain rate from m/s to mm/hr
rainrt(i,k)=qric(i,k)*rho(i,k)*umr(k)/rhow*3600._r8*1000._r8
! vertically-integrated precip source/sink terms (note: grid-averaged)
qrtot = max(qrtot+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
qstot = max(qstot+qnitend(i,k)*dz(i,k)*rho(i,k),0._r8)
nrtot = max(nrtot+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
nstot = max(nstot+nstend(i,k)*dz(i,k)*rho(i,k),0._r8)
! calculate melting and freezing of precip
! melt snow at +2 C
if (t(i,k)+tlat(i,k)/cp*deltat > 275.15_r8) then
if (qstot > 0._r8) then
! make sure melting snow doesn't reduce temperature below threshold
dum = -xlf/cp*qstot/(dz(i,k)*rho(i,k))
if (t(i,k)+tlat(i,k)/cp*deltat+dum.lt.275.15_r8) then
dum = (t(i,k)+tlat(i,k)/cp*deltat-275.15_r8)*cp/xlf
dum = dum/(xlf/cp*qstot/(dz(i,k)*rho(i,k)))
dum = max(0._r8,dum)
dum = min(1._r8,dum)
else
dum = 1._r8
end if
qric(i,k)=qric(i,k)+dum*qniic(i,k)
nric(i,k)=nric(i,k)+dum*nsic(i,k)
qniic(i,k)=(1._r8-dum)*qniic(i,k)
nsic(i,k)=(1._r8-dum)*nsic(i,k)
! heating tendency
tmp=-xlf*dum*qstot/(dz(i,k)*rho(i,k))
meltsdt(i,k)=meltsdt(i,k) + tmp
tlat(i,k)=tlat(i,k)+tmp
qrtot=qrtot+dum*qstot
nrtot=nrtot+dum*nstot
qstot=(1._r8-dum)*qstot
nstot=(1._r8-dum)*nstot
preci(i)=(1._r8-dum)*preci(i)
end if
end if
! freeze all rain at -5C for Arctic
if (t(i,k)+tlat(i,k)/cp*deltat < (tmelt - 5._r8)) then
if (qrtot > 0._r8) then
! make sure freezing rain doesn't increase temperature above threshold
dum = xlf/cp*qrtot/(dz(i,k)*rho(i,k))
if (t(i,k)+tlat(i,k)/cp*deltat+dum.gt.(tmelt - 5._r8)) then
dum = -(t(i,k)+tlat(i,k)/cp*deltat-(tmelt-5._r8))*cp/xlf
dum = dum/(xlf/cp*qrtot/(dz(i,k)*rho(i,k)))
dum = max(0._r8,dum)
dum = min(1._r8,dum)
else
dum = 1._r8
end if
qniic(i,k)=qniic(i,k)+dum*qric(i,k)
nsic(i,k)=nsic(i,k)+dum*nric(i,k)
qric(i,k)=(1._r8-dum)*qric(i,k)
nric(i,k)=(1._r8-dum)*nric(i,k)
! heating tendency
tmp = xlf*dum*qrtot/(dz(i,k)*rho(i,k))
frzrdt(i,k)=frzrdt(i,k) + tmp
tlat(i,k)=tlat(i,k)+tmp
qstot=qstot+dum*qrtot
qrtot=(1._r8-dum)*qrtot
nstot=nstot+dum*nrtot
nrtot=(1._r8-dum)*nrtot
preci(i)=preci(i)+dum*(prect(i)-preci(i))
end if
end if
! if rain/snow mix ratio is zero so should number concentration
if (qniic(i,k).lt.qsmall) then
qniic(i,k)=0._r8
nsic(i,k)=0._r8
end if
if (qric(i,k).lt.qsmall) then
qric(i,k)=0._r8
nric(i,k)=0._r8
end if
! make sure number concentration is a positive number to avoid
! taking root of negative
nric(i,k)=max(nric(i,k),0._r8)
nsic(i,k)=max(nsic(i,k),0._r8)
!.......................................................................
! get size distribution parameters for fallspeed calculations
!......................................................................
! rain
if (qric(i,k).ge.qsmall) then
lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
n0r(k) = nric(i,k)*lamr(k)
! check for slope
! change lammax and lammin for rain and snow
! adjust vars
if (lamr(k).lt.lamminr) then
lamr(k) = lamminr
n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
nric(i,k) = n0r(k)/lamr(k)
else if (lamr(k).gt.lammaxr) then
lamr(k) = lammaxr
n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
nric(i,k) = n0r(k)/lamr(k)
end if
! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
else
lamr(k) = 0._r8
n0r(k) = 0._r8
umr(k)=0._r8
unr(k)=0._r8
end if
!calculate mean size of combined rain and snow
if (lamr(k).gt.0._r8) then
Artmp = n0r(k) * pi / (2 * lamr(k)**3._r8)
else
Artmp = 0._r8
endif
if (lamc(k).gt.0._r8) then
Actmp = cdist1(k) * pi * gamma(pgam(k)+3._r8)/(4._r8 * lamc(k)**2._r8)
else
Actmp = 0._r8
endif
if (Actmp.gt.0_r8.or.Artmp.gt.0) then
rercld(i,k)=rercld(i,k) + 3._r8 *(qric(i,k) + qcic(i,k)) / (4._r8 * rhow * (Actmp + Artmp))
arcld(i,k)=arcld(i,k)+1._r8
endif
!......................................................................
! snow
if (qniic(i,k).ge.qsmall) then
lams(k) = (cons6*cs*nsic(i,k)/ &
qniic(i,k))**(1._r8/ds)
n0s(k) = nsic(i,k)*lams(k)
! check for slope
! adjust vars
if (lams(k).lt.lammins) then
lams(k) = lammins
n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
nsic(i,k) = n0s(k)/lams(k)
else if (lams(k).gt.lammaxs) then
lams(k) = lammaxs
n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
nsic(i,k) = n0s(k)/lams(k)
end if
! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
else
lams(k) = 0._r8
n0s(k) = 0._r8
ums(k) = 0._r8
uns(k) = 0._r8
end if
!c........................................................................
! sum over sub-step for average process rates
! convert rain/snow q and N for output to history, note,
! output is for gridbox average
qrout(i,k)=qrout(i,k)+qric(i,k)*cldmax(i,k)
qsout(i,k)=qsout(i,k)+qniic(i,k)*cldmax(i,k)
nrout(i,k)=nrout(i,k)+nric(i,k)*rho(i,k)*cldmax(i,k)
nsout(i,k)=nsout(i,k)+nsic(i,k)*rho(i,k)*cldmax(i,k)
tlat1(i,k)=tlat1(i,k)+tlat(i,k)
qvlat1(i,k)=qvlat1(i,k)+qvlat(i,k)
qctend1(i,k)=qctend1(i,k)+qctend(i,k)
qitend1(i,k)=qitend1(i,k)+qitend(i,k)
nctend1(i,k)=nctend1(i,k)+nctend(i,k)
nitend1(i,k)=nitend1(i,k)+nitend(i,k)
t(i,k)=t(i,k)+tlat(i,k)*deltat/cpp
q(i,k)=q(i,k)+qvlat(i,k)*deltat
qc(i,k)=qc(i,k)+qctend(i,k)*deltat
qi(i,k)=qi(i,k)+qitend(i,k)*deltat
nc(i,k)=nc(i,k)+nctend(i,k)*deltat
ni(i,k)=ni(i,k)+nitend(i,k)*deltat
rainrt1(i,k)=rainrt1(i,k)+rainrt(i,k)
!divide rain radius over substeps for average
if (arcld(i,k) .gt. 0._r8) then
rercld(i,k)=rercld(i,k)/arcld(i,k)
end if
!calculate precip fluxes and adding them to summing sub-stepping variables
!! flux is zero at top interface
rflx(i,1)=0.0_r8
sflx(i,1)=0.0_r8
!! calculating the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s)
rflx(i,k+1)=qrout(i,k)*rho(i,k)*umr(k)
sflx(i,k+1)=qsout(i,k)*rho(i,k)*ums(k)
!! add to summing sub-stepping variable
rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1)
sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1)
!c........................................................................
end do ! k loop
prect1(i)=prect1(i)+prect(i)
preci1(i)=preci1(i)+preci(i)
end do ! it loop, sub-step
#ifdef MODAL_AERO ! rce 2010/05/01
do k = 1, pver
rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)/max(qcsum_rate1ord(k),1.0e-30_r8)
end do
#endif
300 continue ! continue if no cloud water
end do ! i loop
! convert dt from sub-step back to full time step
deltat=deltat*real(iter)
!c.............................................................................
do i=1,ncol
! skip all calculations if no cloud water
if (ltrue(i).eq.0) then
do k=1,pver
! assign default values for effective radius
effc(i,k)=10._r8
effi(i,k)=25._r8
effc_fn(i,k)=10._r8
lamcrad(i,k)=0._r8
pgamrad(i,k)=0._r8
deffi(i,k)=0._r8
end do
goto 500
end if
! initialize nstep for sedimentation sub-steps
nstep = 1
! divide precip rate by number of sub-steps to get average over time step
prect(i)=prect1(i)/real(iter)
preci(i)=preci1(i)/real(iter)
do k=1,pver
! assign variables back to start-of-timestep values before updating after sub-steps
t(i,k)=t1(i,k)
q(i,k)=q1(i,k)
qc(i,k)=qc1(i,k)
qi(i,k)=qi1(i,k)
nc(i,k)=nc1(i,k)
ni(i,k)=ni1(i,k)
! divide microphysical tendencies by number of sub-steps to get average over time step
tlat(i,k)=tlat1(i,k)/real(iter)
qvlat(i,k)=qvlat1(i,k)/real(iter)
qctend(i,k)=qctend1(i,k)/real(iter)
qitend(i,k)=qitend1(i,k)/real(iter)
nctend(i,k)=nctend1(i,k)/real(iter)
nitend(i,k)=nitend1(i,k)/real(iter)
rainrt(i,k)=rainrt1(i,k)/real(iter)
! divide by number of sub-steps to find final values
rflx(i,k+1)=rflx1(i,k+1)/real(iter)
sflx(i,k+1)=sflx1(i,k+1)/real(iter)
! divide output precip q and N by number of sub-steps to get average over time step
qrout(i,k)=qrout(i,k)/real(iter)
qsout(i,k)=qsout(i,k)/real(iter)
nrout(i,k)=nrout(i,k)/real(iter)
nsout(i,k)=nsout(i,k)/real(iter)
! divide trop_mozart variables by number of sub-steps to get average over time step
nevapr(i,k) = nevapr(i,k)/real(iter)
evapsnow(i,k) = evapsnow(i,k)/real(iter)
prain(i,k) = prain(i,k)/real(iter)
prodsnow(i,k) = prodsnow(i,k)/real(iter)
cmeout(i,k) = cmeout(i,k)/real(iter)
cmeiout(i,k) = cmeiout(i,k)/real(iter)
meltsdt(i,k) = meltsdt(i,k)/real(iter)
frzrdt (i,k) = frzrdt (i,k)/real(iter)
! microphysics output
prao(i,k)=prao(i,k)/real(iter)
prco(i,k)=prco(i,k)/real(iter)
mnuccco(i,k)=mnuccco(i,k)/real(iter)
mnuccto(i,k)=mnuccto(i,k)/real(iter)
msacwio(i,k)=msacwio(i,k)/real(iter)
psacwso(i,k)=psacwso(i,k)/real(iter)
bergso(i,k)=bergso(i,k)/real(iter)
bergo(i,k)=bergo(i,k)/real(iter)
prcio(i,k)=prcio(i,k)/real(iter)
praio(i,k)=praio(i,k)/real(iter)
mnuccro(i,k)=mnuccro(i,k)/real(iter)
pracso (i,k)=pracso (i,k)/real(iter)
! modify to include snow. in prain & evap (diagnostic here: for wet dep)
nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
prain(i,k) = prain(i,k) + prodsnow(i,k)
!--ag
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calculate sedimentation for cloud water and ice
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! update in-cloud cloud mixing ratio and number concentration
! with microphysical tendencies to calculate sedimentation, assign to dummy vars
! note: these are in-cloud values***, hence we divide by cloud fraction
dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)
! obtain new slope parameter to avoid possible singularity
if (dumi(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
lami(k) = (cons1*ci* &
dumni(i,k)/dumi(i,k))**(1._r8/di)
lami(k)=max(lami(k),lammini)
lami(k)=min(lami(k),lammaxi)
else
lami(k)=0._r8
end if
if (dumc(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
! add lower limit to in-cloud number concentration
dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
! pgam(k)=0.0005714_r8*(dumnc(i,k)/1.e6_r8/rho(i,k))+0.2714_r8
! hm fix, not active yet
pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
pgam(k)=1._r8/(pgam(k)**2)-1._r8
pgam(k)=max(pgam(k),2._r8)
pgam(k)=min(pgam(k),15._r8)
lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
(dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
lammin = (pgam(k)+1._r8)/50.e-6_r8
lammax = (pgam(k)+1._r8)/2.e-6_r8
lamc(k)=max(lamc(k),lammin)
lamc(k)=min(lamc(k),lammax)
else
lamc(k)=0._r8
end if
! calculate number and mass weighted fall velocity for droplets
! include effects of sub-grid distribution of cloud water
if (dumc(i,k).ge.qsmall) then
! unc = cons15/(cons3*cons21)* &
! acn(i,k)*gamma(1._r8+bc+pgam(k))/ &
! (lamc(k)**bc*gamma(pgam(k)+1._r8))
! umc = cons15/(cons3*cons21)* &
! acn(i,k)*gamma(4._r8+bc+pgam(k))/ &
! (lamc(k)**bc*gamma(pgam(k)+4._r8))
! hm fix, remove sub-grid qc from fallspeed calculation
unc = &
acn(i,k)*gamma(1._r8+bc+pgam(k))/ &
(lamc(k)**bc*gamma(pgam(k)+1._r8))
umc = &
acn(i,k)*gamma(4._r8+bc+pgam(k))/ &
(lamc(k)**bc*gamma(pgam(k)+4._r8))
! fallspeed for output
vtrmc(i,k)=umc
else
umc = 0._r8
unc = 0._r8
end if
! calculate number and mass weighted fall velocity for cloud ice
if (dumi(i,k).ge.qsmall) then
uni = ain(i,k)*cons16/lami(k)**bi
umi = ain(i,k)*cons17/(6._r8*lami(k)**bi)
uni=min(uni,1.2_r8*rhof(i,k))
umi=min(umi,1.2_r8*rhof(i,k))
! fallspeed
vtrmi(i,k)=umi
else
umi = 0._r8
uni = 0._r8
end if
fi(k) = g*rho(i,k)*umi
fni(k) = g*rho(i,k)*uni
fc(k) = g*rho(i,k)*umc
fnc(k) = g*rho(i,k)*unc
! calculate number of split time steps to ensure courant stability criteria
! for sedimentation calculations
rgvm = max(fi(k),fc(k),fni(k),fnc(k))
nstep = max(int(rgvm*deltat/pdel(i,k)+1._r8),nstep)
! redefine dummy variables - sedimentation is calculated over grid-scale
! quantities to ensure conservation
dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)
if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
end do !!! vertical loop
do n = 1,nstep !! loop over sub-time step to ensure stability
do k = 1,pver
falouti(k) = fi(k)*dumi(i,k)
faloutni(k) = fni(k)*dumni(i,k)
faloutc(k) = fc(k)*dumc(i,k)
faloutnc(k) = fnc(k)*dumnc(i,k)
end do
! top of model
k = 1
faltndi = falouti(k)/pdel(i,k)
faltndni = faloutni(k)/pdel(i,k)
faltndc = faloutc(k)/pdel(i,k)
faltndnc = faloutnc(k)/pdel(i,k)
! add fallout terms to microphysical tendencies
qitend(i,k) = qitend(i,k)-faltndi/nstep
nitend(i,k) = nitend(i,k)-faltndni/nstep
qctend(i,k) = qctend(i,k)-faltndc/nstep
nctend(i,k) = nctend(i,k)-faltndnc/nstep
! sedimentation tendencies for output
qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
qisedten(i,k)=qisedten(i,k)-faltndi/nstep
dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
do k = 2,pver
! for cloud liquid and ice, if cloud fraction increases with height
! then add flux from above to both vapor and cloud water of current level
! this means that flux entering clear portion of cell from above evaporates
! instantly
dum=lcldm(i,k)/lcldm(i,k-1)
dum=min(dum,1._r8)
dum1=icldm(i,k)/icldm(i,k-1)
dum1=min(dum1,1._r8)
faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k)
faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k)
faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k)
faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k)
faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k)
faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k)
! add fallout terms to eulerian tendencies
qitend(i,k) = qitend(i,k)-faltndi/nstep
nitend(i,k) = nitend(i,k)-faltndni/nstep
qctend(i,k) = qctend(i,k)-faltndc/nstep
nctend(i,k) = nctend(i,k)-faltndnc/nstep
! sedimentation tendencies for output
qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
qisedten(i,k)=qisedten(i,k)-faltndi/nstep
! add terms to to evap/sub of cloud water
qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep
! for output
qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep
qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep
! for output
qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep
tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep
tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep
dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
Fni(K)=MAX(Fni(K)/pdel(i,K),Fni(K-1)/pdel(i,K-1))*pdel(i,K)
FI(K)=MAX(FI(K)/pdel(i,K),FI(K-1)/pdel(i,K-1))*pdel(i,K)
fnc(k)=max(fnc(k)/pdel(i,k),fnc(k-1)/pdel(i,k-1))*pdel(i,k)
Fc(K)=MAX(Fc(K)/pdel(i,K),Fc(K-1)/pdel(i,K-1))*pdel(i,K)
end do !! k loop
! units below are m/s
! cloud water/ice sedimentation flux at surface
! is added to precip flux at surface to get total precip (cloud + precip water)
! rate
prect(i) = prect(i)+(faloutc(pver)+falouti(pver)) &
/g/nstep/1000._r8
preci(i) = preci(i)+(falouti(pver)) &
/g/nstep/1000._r8
end do !! nstep loop
! end sedimentation
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! get new update for variables that includes sedimentation tendency
! note : here dum variables are grid-average, NOT in-cloud
do k=1,pver
dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)
if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
! calculate instantaneous processes (melting, homogeneous freezing)
if (t(i,k)+tlat(i,k)/cp*deltat > tmelt) then
if (dumi(i,k) > 0._r8) then
! limit so that melting does not push temperature below freezing
dum = -dumi(i,k)*xlf/cp
if (t(i,k)+tlat(i,k)/cp*deltat+dum.lt.tmelt) then
dum = (t(i,k)+tlat(i,k)/cp*deltat-tmelt)*cp/xlf
dum = dum/dumi(i,k)*xlf/cp
dum = max(0._r8,dum)
dum = min(1._r8,dum)
else
dum = 1._r8
end if
qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat
! for output
melto(i,k)=dum*dumi(i,k)/deltat
! assume melting ice produces droplet
! mean volume radius of 8 micron
nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ &
(4._r8*pi*5.12e-16_r8*rhow)
qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat
nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat
tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat
end if
end if
! homogeneously freeze droplets at -40 C
if (t(i,k)+tlat(i,k)/cp*deltat < 233.15_r8) then
if (dumc(i,k) > 0._r8) then
! limit so that freezing does not push temperature above threshold
dum = dumc(i,k)*xlf/cp
if (t(i,k)+tlat(i,k)/cp*deltat+dum.gt.233.15_r8) then
dum = -(t(i,k)+tlat(i,k)/cp*deltat-233.15_r8)*cp/xlf
dum = dum/dumc(i,k)*xlf/cp
dum = max(0._r8,dum)
dum = min(1._r8,dum)
else
dum = 1._r8
end if
qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat
! for output
homoo(i,k)=dum*dumc(i,k)/deltat
! assume 25 micron mean volume radius of homogeneously frozen droplets
! consistent with size of detrained ice in stratiform.F90
nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* &
500._r8)/deltat
qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat
nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat
tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat
end if
end if
! remove any excess over-saturation, which is possible due to non-linearity when adding
! together all microphysical processes
! follow code similar to old CAM scheme
qtmp=q(i,k)+qvlat(i,k)*deltat
ttmp=t(i,k)+tlat(i,k)/cpp*deltat
esn = polysvp(ttmp,0) ! use rhw to allow ice supersaturation
qsn = min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
if (qtmp > qsn .and. qsn > 0) then
! expression below is approximate since there may be ice deposition
dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat
! add to output cme
cmeout(i,k) = cmeout(i,k)+dum
! now add to tendencies, partition between liquid and ice based on temperature
if (ttmp > 268.15_r8) then
dum1=0.0_r8
! now add to tendencies, partition between liquid and ice based on te
else if (ttmp < 238.15_r8) then
dum1=1.0_r8
else
dum1=(268.15_r8-ttmp)/30._r8
end if
dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
*qsn/(cpp*rv*ttmp**2))/deltat
qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
! for output
qcreso(i,k)=dum*(1._r8-dum1)
qitend(i,k)=qitend(i,k)+dum*dum1
qireso(i,k)=dum*dum1
qvlat(i,k)=qvlat(i,k)-dum
! for output
qvres(i,k)=-dum
tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
end if
!...............................................................................
! calculate effective radius for pass to radiation code
! if no cloud water, default value is 10 micron for droplets,
! 25 micron for cloud ice
! update cloud variables after instantaneous processes to get effective radius
! variables are in-cloud to calculate size dist parameters
dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)
! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
dumc(i,k)=min(dumc(i,k),5.e-3_r8)
dumi(i,k)=min(dumi(i,k),5.e-3_r8)
!...................
! cloud ice effective radius
if (dumi(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
lami(k) = (cons1*ci* &
dumni(i,k)/dumi(i,k))**(1._r8/di)
if (lami(k).lt.lammini) then
lami(k) = lammini
n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
niic(i,k) = n0i(k)/lami(k)
! adjust number conc if needed to keep mean size in reasonable range
nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
else if (lami(k).gt.lammaxi) then
lami(k) = lammaxi
n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
niic(i,k) = n0i(k)/lami(k)
! adjust number conc if needed to keep mean size in reasonable range
nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
end if
effi(i,k) = 1.5_r8/lami(k)*1.e6_r8
else
effi(i,k) = 25._r8
end if
!...................
! cloud droplet effective radius
if (dumc(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! set tendency to ensure minimum droplet concentration
! after update by microphysics, except when lambda exceeds bounds on mean drop
! size or if there is no cloud water
if (dumnc(i,k).lt.cdnl/rho(i,k)) then
nctend(i,k)=(cdnl/rho(i,k)*cldm(i,k)-nc(i,k))/deltat
end if
dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! pgam(k)=0.0005714_r8*(dumnc(i,k)/1.e6_r8/rho(i,k))+0.2714_r8
! hm fix, not active yet
pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
pgam(k)=1._r8/(pgam(k)**2)-1._r8
pgam(k)=max(pgam(k),2._r8)
pgam(k)=min(pgam(k),15._r8)
lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
(dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
lammin = (pgam(k)+1._r8)/50.e-6_r8
lammax = (pgam(k)+1._r8)/2.e-6_r8
if (lamc(k).lt.lammin) then
lamc(k) = lammin
ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
gamma(pgam(k)+1._r8)/ &
(pi*rhow*gamma(pgam(k)+4._r8))
! adjust number conc if needed to keep mean size in reasonable range
nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
else if (lamc(k).gt.lammax) then
lamc(k) = lammax
ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
gamma(pgam(k)+1._r8)/ &
(pi*rhow*gamma(pgam(k)+4._r8))
! adjust number conc if needed to keep mean size in reasonable range
nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
end if
effc(i,k) = &
gamma(pgam(k)+4._r8)/ &
gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
!assign output fields for shape here
lamcrad(i,k)=lamc(k)
pgamrad(i,k)=pgam(k)
else
effc(i,k) = 10._r8
lamcrad(i,k)=0._r8
pgamrad(i,k)=0._r8
end if
! ice effective diameter for david mitchell's optics
deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8
!!! recalculate effective radius for constant number, in order to separate
! first and second indirect effects
! assume constant number of 10^8 kg-1
dumnc(i,k)=1.e8
if (dumc(i,k).ge.qsmall) then
! pgam(k)=0.0005714_r8*(dumnc(i,k)/1.e6_r8/rho(i,k))+0.2714_r8
! hm fix, not active yet
pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
pgam(k)=1._r8/(pgam(k)**2)-1._r8
pgam(k)=max(pgam(k),2._r8)
pgam(k)=min(pgam(k),15._r8)
lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
(dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
lammin = (pgam(k)+1._r8)/50.e-6_r8
lammax = (pgam(k)+1._r8)/2.e-6_r8
if (lamc(k).lt.lammin) then
lamc(k) = lammin
else if (lamc(k).gt.lammax) then
lamc(k) = lammax
end if
! effc_fn(i,k) = gamma(qcvar+1._r8/3._r8)/(gamma(qcvar)*qcvar**(1._r8/3._r8))* &
! gamma(pgam(k)+4._r8)/ &
! gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
! hm, modify, don't include sub-grid qc
effc_fn(i,k) = &
gamma(pgam(k)+4._r8)/ &
gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
else
effc_fn(i,k) = 10._r8
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
end do ! vertical k loop
500 continue
do k=1,pver
! if updated q (after microphysics) is zero, then ensure updated n is also zero
if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat
if (qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat
end do
end do ! i loop
! ccn concentration as diagnostic
#ifndef MODAL_AERO
do k=1,pver
ccn(:ncol,k,:) = 0.
do m=1,naer_all
if(m.eq.idxsul)then
! Lohmann treatment for sulfate has variable size distribution
do i=1,ncol
if (naer2(i,k,m).gt.0._r8) then
amcubesulfate(i)=amcubefactor(m)*aer_mmr(i,k,m)*rho(i,k)/(naer2(i,k,m))
smcritsulfate(i)=smcritfactor(m)/sqrt(amcubesulfate(i))
else
smcritsulfate(i)=1._r8
endif
end do
end if
do l=1,psat
if(m.eq.idxsul)then
do i=1,ncol
arg=argfactor(m)*log(smcritsulfate(i)/super(l))
if(arg<2)then
if(arg<-2)then
ccnfact(l,m)=1.0e-6_r8
else
ccnfact(l,m)=0.5e-6_r8*erfc(arg)
endif
else
ccnfact(l,m)=0.0_r8
endif
ccn(i,k,l)=ccn(i,k,l)+naer2(i,k,m)*ccnfact(l,m)
end do
else
ccn(:ncol,k,l)=ccn(:ncol,k,l)+naer2(:ncol,k,m)*ccnfact(l,m)
endif
enddo
enddo
enddo
do l=1,psat
call outfld(ccn_name(l), ccn(1,1,l) , pcols, lchnk )
enddo
#endif
do l=1,naer_all
call outfld(aername(l)//'_m3', naer2(1,1,l), pcols, lchnk)
enddo
! hm add rain/snow mixing ratio and number concentration as diagnostic
call outfld('QRAIN',qrout, pcols, lchnk)
call outfld('QSNOW',qsout, pcols, lchnk)
call outfld('NRAIN',nrout, pcols, lchnk)
call outfld('NSNOW',nsout, pcols, lchnk)
! add snow ouptut
do i = 1,ncol
do k=1,pver
if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
dsout(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
endif
end do
end do
call outfld('DSNOW',dsout, pcols, lchnk)
! add precip fluxes as output fields
call outfld('MGFLXPRC',rflx, pcols, lchnk)
call outfld('MGFLXSNW',sflx, pcols, lchnk)
! analytic radar reflectivity
! formulas from Matthew Shupe, NOAA/CERES
! *****note: radar reflectivity is local (in-precip average)
! units of mm^6/m^3
do i = 1,ncol
do k=1,pver
if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall) then
dum=((qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
/(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/cldmax(i,k)
else
dum=0._r8
end if
if (qi(i,k)+qitend(i,k)*deltat.ge.qsmall) then
dum1=((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/cldmax(i,k)
else
dum1=0._r8
end if
if (qsout(i,k).ge.qsmall) then
dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
end if
refl(i,k)=dum+dum1
! add rain rate, but for 37 GHz formulation instead of 94 GHz
! formula approximated from data of Matrasov (2007)
! rainrt is the rain rate in mm/hr
! reflectivity (dum) is in DBz
if (rainrt(i,k).ge.0.001) then
dum=log10(rainrt(i,k)**6._r8)+16._r8
! convert from DBz to mm^6/m^3
dum = 10._r8**(dum/10._r8)
else
! don't include rain rate in R calculation for values less than 0.001 mm/hr
dum=0.
end if
! add to refl
refl(i,k)=refl(i,k)+dum
!output reflectivity in Z.
areflz(i,k)=refl(i,k)
! convert back to DBz
if (refl(i,k).gt.minrefl) then
refl(i,k)=10._r8*log10(refl(i,k))
else
refl(i,k)=-9999._r8
end if
!set averaging flag
if (refl(i,k).gt.mindbz) then
arefl(i,k)=refl(i,k)
frefl(i,k)=1.0_r8
else
arefl(i,k)=0._r8
areflz(i,k)=0._r8
frefl(i,k)=0._r8
end if
! bound cloudsat reflectivity
csrfl(i,k)=min(csmax,refl(i,k))
!set averaging flag
if (csrfl(i,k).gt.csmin) then
acsrfl(i,k)=refl(i,k)
fcsrfl(i,k)=1.0_r8
else
acsrfl(i,k)=0._r8
fcsrfl(i,k)=0._r8
end if
end do
end do
call outfld('REFL',refl, pcols, lchnk)
call outfld('AREFL',arefl, pcols, lchnk)
call outfld('AREFLZ',areflz, pcols, lchnk)
call outfld('FREFL',frefl, pcols, lchnk)
call outfld('CSRFL',csrfl, pcols, lchnk)
call outfld('ACSRFL',acsrfl, pcols, lchnk)
call outfld('FCSRFL',fcsrfl, pcols, lchnk)
call outfld('RERCLD',rercld, pcols, lchnk)
! averaging for snow and rain number and diameter
qrout2(:,:)=0._r8
qsout2(:,:)=0._r8
nrout2(:,:)=0._r8
nsout2(:,:)=0._r8
drout2(:,:)=0._r8
dsout2(:,:)=0._r8
freqs(:,:)=0._r8
freqr(:,:)=0._r8
do i = 1,ncol
do k=1,pver
if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
qrout2(i,k)=qrout(i,k)
nrout2(i,k)=nrout(i,k)
drout2(i,k)=(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)
freqr(i,k)=1._r8
endif
if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
qsout2(i,k)=qsout(i,k)
nsout2(i,k)=nsout(i,k)
dsout2(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
freqs(i,k)=1._r8
endif
end do
end do
! output activated liquid and ice (convert from #/kg -> #/m3)
do i = 1,ncol
do k=1,pver
ncai(i,k)=dum2i(i,k)*rho(i,k)
ncal(i,k)=dum2l(i,k)*rho(i,k)
nihf(i,k)=nihf(i,k)*rho(i,k)
niimm(i,k)=niimm(i,k)*rho(i,k)
nidep(i,k)=nidep(i,k)*rho(i,k)
nimey(i,k)=nimey(i,k)*rho(i,k)
end do
end do
call outfld('NCAL',ncal, pcols,lchnk)
call outfld('NCAI',ncai, pcols,lchnk)
call outfld('NIHF',nihf, pcols,lchnk)
call outfld('NIIMM',niimm, pcols,lchnk)
call outfld('NIDEP',nidep, pcols,lchnk)
call outfld('NIMEY',nimey, pcols,lchnk)
!add averaged output fields.
call outfld('AQRAIN',qrout2, pcols,lchnk)
call outfld('AQSNOW',qsout2, pcols,lchnk)
call outfld('ANRAIN',nrout2, pcols,lchnk)
call outfld('ANSNOW',nsout2, pcols,lchnk)
call outfld('ADRAIN',drout2, pcols,lchnk)
call outfld('ADSNOW',dsout2, pcols,lchnk)
call outfld('FREQR',freqr, pcols,lchnk)
call outfld('FREQS',freqs, pcols,lchnk)
!redefine fice here....
nfice(:,:)=0._r8
do k=1,pver
do i=1,ncol
dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
dumfice=qsout(i,k) + qrout(i,k) + dumc(i,k) + dumi(i,k)
if (dumfice.gt.0._r8) then
nfice(i,k)=(qsout(i,k) + dumi(i,k))/dumfice
endif
enddo
enddo
call outfld('FICE',nfice, pcols, lchnk)
deallocate( &
naermod, &
naer2, &
maerosol, &
ntype )
return
end subroutine mmicro_pcond
!##############################################################################
subroutine findsp1 (lchnk, ncol, q, t, p, tsp, qsp),6
!-----------------------------------------------------------------------
!
! 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 = 1,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,*) ' '
write (iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
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),128._r8),373._r8)
#else
tlim(i) = min(max(t(i,k),173._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) - tt0
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) < tt0) 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
write (iulog,*) ' first guess ', tsp(i,k)
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) - tt0
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) < tt0) 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
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) < 130.16_r8) then
#else
if (tsp(i,k) < 174.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
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
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)
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
call endrun ('FINDSP')
endif
end do
endif
end do ! level loop (k=1,pver)
return
end subroutine findsp1
subroutine findsp1_water (lchnk, ncol, q, t, p, tsp, qsp),6
!-----------------------------------------------------------------------
!
! 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
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
!
! max number of times to iterate the calculation
iter = 8
!
do k = 1,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,*) ' '
write (iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
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),128._r8),373._r8)
#else
tlim(i) = min(max(t(i,k),173._r8),373._r8)
#endif
es = polysvp(tlim(i),0)
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
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
hlatsb = hlatv
if (tlim(i) < tt0) then
hltalt(i,k) = hlatsb
else
hltalt(i,k) = hlatvp
end if
!--xl
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
write (iulog,*) ' first guess ', tsp(i,k)
endif
#endif
es = polysvp(tsp(i,k),0)
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 = polysvp(tsp(i,k),0)
!
! 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
!
!
! No icephs or water to ice transition
!
hlatvp = hlatv - 2369.0*(tsp(i,k)-tt0)
hlatsb = hlatv
if (tsp(i,k) < tt0) then
hltalt(i,k) = hlatsb
else
hltalt(i,k) = hlatvp
end if
desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k))
dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
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 = polysvp(tsp(i,k),0)
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
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) < 130.16_r8) then
#else
if (tsp(i,k) < 174.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
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
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)
write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
call endrun ('FINDSP')
endif
end do
endif
end do ! level loop (k=1,pver)
return
end subroutine findsp1_water
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! error function in single precision
!
! Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
! You may use, copy, modify this code for any purpose and
! without fee. You may distribute this ORIGINAL package.
function derf(x) 1,1
implicit real (a - h, o - z)
real(r8) a,b,x
dimension a(0 : 64), b(0 : 64)
integer i,k
data (a(i), i = 0, 12) / &
0.00000000005958930743d0, -0.00000000113739022964d0, &
0.00000001466005199839d0, -0.00000016350354461960d0, &
0.00000164610044809620d0, -0.00001492559551950604d0, &
0.00012055331122299265d0, -0.00085483269811296660d0, &
0.00522397762482322257d0, -0.02686617064507733420d0, &
0.11283791670954881569d0, -0.37612638903183748117d0, &
1.12837916709551257377d0 /
data (a(i), i = 13, 25) / &
0.00000000002372510631d0, -0.00000000045493253732d0, &
0.00000000590362766598d0, -0.00000006642090827576d0, &
0.00000067595634268133d0, -0.00000621188515924000d0, &
0.00005103883009709690d0, -0.00037015410692956173d0, &
0.00233307631218880978d0, -0.01254988477182192210d0, &
0.05657061146827041994d0, -0.21379664776456006580d0, &
0.84270079294971486929d0 /
data (a(i), i = 26, 38) / &
0.00000000000949905026d0, -0.00000000018310229805d0, &
0.00000000239463074000d0, -0.00000002721444369609d0, &
0.00000028045522331686d0, -0.00000261830022482897d0, &
0.00002195455056768781d0, -0.00016358986921372656d0, &
0.00107052153564110318d0, -0.00608284718113590151d0, &
0.02986978465246258244d0, -0.13055593046562267625d0, &
0.67493323603965504676d0 /
data (a(i), i = 39, 51) / &
0.00000000000382722073d0, -0.00000000007421598602d0, &
0.00000000097930574080d0, -0.00000001126008898854d0, &
0.00000011775134830784d0, -0.00000111992758382650d0, &
0.00000962023443095201d0, -0.00007404402135070773d0, &
0.00050689993654144881d0, -0.00307553051439272889d0, &
0.01668977892553165586d0, -0.08548534594781312114d0, &
0.56909076642393639985d0 /
data (a(i), i = 52, 64) / &
0.00000000000155296588d0, -0.00000000003032205868d0, &
0.00000000040424830707d0, -0.00000000471135111493d0, &
0.00000005011915876293d0, -0.00000048722516178974d0, &
0.00000430683284629395d0, -0.00003445026145385764d0, &
0.00024879276133931664d0, -0.00162940941748079288d0, &
0.00988786373932350462d0, -0.05962426839442303805d0, &
0.49766113250947636708d0 /
data (b(i), i = 0, 12) / &
-0.00000000029734388465d0, 0.00000000269776334046d0, &
-0.00000000640788827665d0, -0.00000001667820132100d0, &
-0.00000021854388148686d0, 0.00000266246030457984d0, &
0.00001612722157047886d0, -0.00025616361025506629d0, &
0.00015380842432375365d0, 0.00815533022524927908d0, &
-0.01402283663896319337d0, -0.19746892495383021487d0,&
0.71511720328842845913d0 /
data (b(i), i = 13, 25) / &
-0.00000000001951073787d0, -0.00000000032302692214d0, &
0.00000000522461866919d0, 0.00000000342940918551d0, &
-0.00000035772874310272d0, 0.00000019999935792654d0, &
0.00002687044575042908d0, -0.00011843240273775776d0, &
-0.00080991728956032271d0, 0.00661062970502241174d0, &
0.00909530922354827295d0, -0.20160072778491013140d0, &
0.51169696718727644908d0 /
data (b(i), i = 26, 38) / &
0.00000000003147682272d0, -0.00000000048465972408d0, &
0.00000000063675740242d0, 0.00000003377623323271d0, &
-0.00000015451139637086d0, -0.00000203340624738438d0,&
0.00001947204525295057d0, 0.00002854147231653228d0, &
-0.00101565063152200272d0, 0.00271187003520095655d0, &
0.02328095035422810727d0, -0.16725021123116877197d0, &
0.32490054966649436974d0 /
data (b(i), i = 39, 51) / &
0.00000000002319363370d0, -0.00000000006303206648d0, &
-0.00000000264888267434d0, 0.00000002050708040581d0, &
0.00000011371857327578d0, -0.00000211211337219663d0, &
0.00000368797328322935d0, 0.00009823686253424796d0, &
-0.00065860243990455368d0, -0.00075285814895230877d0,&
0.02585434424202960464d0, -0.11637092784486193258d0, &
0.18267336775296612024d0 /
data (b(i), i = 52, 64) / &
-0.00000000000367789363d0, 0.00000000020876046746d0, &
-0.00000000193319027226d0, -0.00000000435953392472d0, &
0.00000018006992266137d0, -0.00000078441223763969d0, &
-0.00000675407647949153d0, 0.00008428418334440096d0, &
-0.00017604388937031815d0, -0.00239729611435071610d0, &
0.02064129023876022970d0, -0.06905562880005864105d0, &
0.09084526782065478489d0 /
w = abs(x)
if (w .lt. 2.2d0) then
t = w * w
k = int(t)
t = t - k
k = k * 13
y = ((((((((((((a(k) * t + a(k + 1)) * t + &
a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + &
a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + &
a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + &
a(k + 11)) * t + a(k + 12)) * w
else if (w .lt. 6.9d0) then
k = int(w)
t = w - k
k = 13 * (k - 2)
y = (((((((((((b(k) * t + b(k + 1)) * t + &
b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + &
b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + &
b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + &
b(k + 11)) * t + b(k + 12)
y = y * y
y = y * y
y = y * y
y = 1 - y * y
else
y = 1
end if
if (x .lt. 0) y = -y
derf = y
end function derf
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
FUNCTION GAMMA(X) 90
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!D DOUBLE PRECISION FUNCTION DGAMMA(X)
!----------------------------------------------------------------------
!
! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
! MACHINE-DEPENDENT CONSTANTS.
!
!
!*******************************************************************
!*******************************************************************
!
! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
!
! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
! GAMMA(XBIG) = BETA**MAXEXP
! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
! APPROXIMATELY BETA**MAXEXP
! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
! 1.0+EPS .GT. 1.0
! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
! 1/XMININ IS MACHINE REPRESENTABLE
!
! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
!
! BETA MAXEXP XBIG
!
! CRAY-1 (S.P.) 2 8191 966.961
! CYBER 180/855
! UNDER NOS (S.P.) 2 1070 177.803
! IEEE (IBM/XT,
! SUN, ETC.) (S.P.) 2 128 35.040
! IEEE (IBM/XT,
! SUN, ETC.) (D.P.) 2 1024 171.624
! IBM 3033 (D.P.) 16 63 57.574
! VAX D-FORMAT (D.P.) 2 127 34.844
! VAX G-FORMAT (D.P.) 2 1023 171.489
!
! XINF EPS XMININ
!
! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
! CYBER 180/855
! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
! IEEE (IBM/XT,
! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
! IEEE (IBM/XT,
! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
!
!*******************************************************************
!*******************************************************************
!
! ERROR RETURNS
!
! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
! TO BE FREE OF UNDERFLOW AND OVERFLOW.
!
!
! INTRINSIC FUNCTIONS REQUIRED ARE:
!
! INT, DBLE, EXP, LOG, REAL, SIN
!
!
! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
! (ED.), SPRINGER VERLAG, BERLIN, 1976.
!
! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
! SONS, NEW YORK, 1968.
!
! LATEST MODIFICATION: OCTOBER 12, 1989
!
! AUTHORS: W. J. CODY AND L. STOLTZ
! APPLIED MATHEMATICS DIVISION
! ARGONNE NATIONAL LABORATORY
! ARGONNE, IL 60439
!
!----------------------------------------------------------------------
INTEGER I,N
LOGICAL PARITY
real(r8) gamma
REAL(r8) &
!D DOUBLE PRECISION
C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
DIMENSION C(7),P(8),Q(8)
!----------------------------------------------------------------------
! MATHEMATICAL CONSTANTS
!----------------------------------------------------------------------
DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0_r8,0.5E0_r8,12.0E0_r8,2.0E0_r8,0.0E0_r8/, &
SQRTPI/0.9189385332046727417803297E0_r8/, &
PI/3.1415926535897932384626434E0_r8/
!D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
!D 1 SQRTPI/0.9189385332046727417803297D0/,
!D 2 PI/3.1415926535897932384626434D0/
!----------------------------------------------------------------------
! MACHINE DEPENDENT PARAMETERS
!----------------------------------------------------------------------
DATA XBIG,XMININ,EPS/35.040E0_r8,1.18E-38_r8,1.19E-7_r8/, &
XINF/3.4E38_r8/
!D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
!D 1 XINF/1.79D308/
!----------------------------------------------------------------------
! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
! APPROXIMATION OVER (1,2).
!----------------------------------------------------------------------
DATA P/-1.71618513886549492533811E+0_r8,2.47656508055759199108314E+1_r8,&
-3.79804256470945635097577E+2_r8,6.29331155312818442661052E+2_r8,&
8.66966202790413211295064E+2_r8,-3.14512729688483675254357E+4_r8,&
-3.61444134186911729807069E+4_r8,6.64561438202405440627855E+4_r8/
DATA Q/-3.08402300119738975254353E+1_r8,3.15350626979604161529144E+2_r8,&
-1.01515636749021914166146E+3_r8,-3.10777167157231109440444E+3_r8,&
2.25381184209801510330112E+4_r8,4.75584627752788110767815E+3_r8,&
-1.34659959864969306392456E+5_r8,-1.15132259675553483497211E+5_r8/
!D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
!D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
!D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
!D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
!D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
!D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
!D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
!D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
!----------------------------------------------------------------------
! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
!----------------------------------------------------------------------
DATA C/-1.910444077728E-03_r8,8.4171387781295E-04_r8, &
-5.952379913043012E-04_r8,7.93650793500350248E-04_r8,&
-2.777777777777681622553E-03_r8,8.333333333333333331554247E-02_r8,&
5.7083835261E-03_r8/
!D DATA C/-1.910444077728D-03,8.4171387781295D-04,
!D 1 -5.952379913043012D-04,7.93650793500350248D-04,
!D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
!D 3 5.7083835261D-03/
!----------------------------------------------------------------------
! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
!----------------------------------------------------------------------
CONV(I) = REAL(I,r8)
!D CONV(I) = DBLE(I)
PARITY=.FALSE.
FACT=ONE
N=0
Y=X
IF(Y.LE.ZERO)THEN
!----------------------------------------------------------------------
! ARGUMENT IS NEGATIVE
!----------------------------------------------------------------------
Y=-X
Y1=AINT(Y)
RES=Y-Y1
IF(RES.NE.ZERO)THEN
IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
FACT=-PI/SIN(PI*RES)
Y=Y+ONE
ELSE
RES=XINF
GOTO 900
ENDIF
ENDIF
!----------------------------------------------------------------------
! ARGUMENT IS POSITIVE
!----------------------------------------------------------------------
IF(Y.LT.EPS)THEN
!----------------------------------------------------------------------
! ARGUMENT .LT. EPS
!----------------------------------------------------------------------
IF(Y.GE.XMININ)THEN
RES=ONE/Y
ELSE
RES=XINF
GOTO 900
ENDIF
ELSEIF(Y.LT.TWELVE)THEN
Y1=Y
IF(Y.LT.ONE)THEN
!----------------------------------------------------------------------
! 0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
Z=Y
Y=Y+ONE
ELSE
!----------------------------------------------------------------------
! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
!----------------------------------------------------------------------
N=INT(Y)-1
Y=Y-CONV(N)
Z=Y-ONE
ENDIF
!----------------------------------------------------------------------
! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
!----------------------------------------------------------------------
XNUM=ZERO
XDEN=ONE
DO 260 I=1,8
XNUM=(XNUM+P(I))*Z
XDEN=XDEN*Z+Q(I)
260 CONTINUE
RES=XNUM/XDEN+ONE
IF(Y1.LT.Y)THEN
!----------------------------------------------------------------------
! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
RES=RES/Y1
ELSEIF(Y1.GT.Y)THEN
!----------------------------------------------------------------------
! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
!----------------------------------------------------------------------
DO 290 I=1,N
RES=RES*Y
Y=Y+ONE
290 CONTINUE
ENDIF
ELSE
!----------------------------------------------------------------------
! EVALUATE FOR ARGUMENT .GE. 12.0,
!----------------------------------------------------------------------
IF(Y.LE.XBIG)THEN
YSQ=Y*Y
SUM=C(7)
DO 350 I=1,6
SUM=SUM/YSQ+C(I)
350 CONTINUE
SUM=SUM/Y-Y+SQRTPI
SUM=SUM+(Y-HALF)*LOG(Y)
RES=EXP(SUM)
ELSE
RES=XINF
GOTO 900
ENDIF
ENDIF
!----------------------------------------------------------------------
! FINAL ADJUSTMENTS AND RETURN
!----------------------------------------------------------------------
IF(PARITY)RES=-RES
IF(FACT.NE.ONE)RES=FACT/RES
900 GAMMA=RES
!D900 DGAMMA = RES
RETURN
! ---------- LAST LINE OF GAMMA ----------
END function gamma
subroutine activate(wbar, tair, rhoair, & 1,5
na, pmode, nmode, ma, sigman, hygro, rhodry, nact)
! calculates number fraction of aerosols activated as CCN
! assumes an internal mixture within each of up to pmode multiple aerosol modes
! a gaussiam spectrum of updrafts can be treated.
! mks units
! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
use physconst, only: rair, epsilo, cpair, rh2o, latvap, gravit, &
rhoh2o, mwh2o, r_universal
use wv_saturation, only: estblf, epsqs
implicit none
! input
integer pmode,ptype ! dimension of modes, types in modes
real(r8) wbar ! grid cell mean vertical velocity (m/s)
real(r8) tair ! air temperature (K)
real(r8) rhoair ! air density (kg/m3)
real(r8) na(pmode) ! aerosol number concentration (/m3)
integer nmode ! number of aerosol modes
real(r8) ma(pmode) ! aerosol mass concentration (kg/m3)
real(r8) rhodry(pmode) ! density of aerosol material
real(r8) sigman(pmode) ! geometric standard deviation of aerosol size distribution
real(r8) hygro(pmode) ! hygroscopicity of aerosol mode
! output
real(r8) nact ! number fraction of aerosols activated
! local
! external erf,erfc
! real(r8) erf,erfc
#if (defined AIX)
#define ERF erf
#define ERFC erfc
#else
#define ERF derf
#define ERFC derfc
real(r8) derf,derfc
#endif
integer, parameter:: nx=200
integer :: maxmodes
real(r8) surften ! surface tension of water w/respect to air (N/m)
data surften/0.076/
save surften
real(r8) p0 ! reference pressure (Pa)
data p0/1013.25e2/
save p0
real(r8), allocatable :: volc(:) ! total aerosol volume concentration (m3/m3)
real(r8) tmass ! total aerosol mass concentration (g/cm3)
real(r8) rm ! number mode radius of aerosol at max supersat (cm)
real(r8) pres ! pressure (Pa)
real(r8) path ! mean free path (m)
real(r8) diff ! diffusivity (m2/s)
real(r8) conduct ! thermal conductivity (Joule/m/sec/deg)
real(r8) diff0,conduct0
real(r8) qs ! water vapor saturation mixing ratio
real(r8) dqsdt ! change in qs with temperature
real(r8) dqsdp ! change in qs with pressure
real(r8) gloc ! thermodynamic function (m2/s)
real(r8) zeta
real(r8), allocatable :: eta(:)
real(r8), allocatable :: smc(:)
real(r8) lnsmax ! ln(smax)
real(r8) alpha
real(r8) gammaloc
real(r8) beta
real(r8) sqrtg
real(r8) alogam
real(r8) rlo,rhi,xint1,xint2,xint3,xint4
real(r8) w,wnuc,wb
real(r8) alw,sqrtalw
real(r8) smax
real(r8) x,arg
real(r8) xmincoeff,xcut,volcut,surfcut
real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
real(r8) :: etafactor1,etafactor2max
real(r8),allocatable :: etafactor2(:)
real(r8) es
integer m,n
real(r8),allocatable :: amcubeloc(:)
real(r8),allocatable :: lnsmloc(:)
maxmodes = naer_all
allocate( &
volc(maxmodes), &
eta(maxmodes), &
smc(maxmodes), &
etafactor2(maxmodes), &
amcubeloc(maxmodes), &
lnsmloc(maxmodes) )
if(maxmodes<pmode)then
write(iulog,*)'maxmodes,pmode in activate =',maxmodes,pmode
call endrun('activate')
endif
nact=0._r8
if(nmode.eq.1.and.na(1).lt.1.e-20)return
if(wbar.le.0.)return
pres=rair*rhoair*tair
diff0=0.211e-4*(p0/pres)*(tair/tt0)**1.94
conduct0=(5.69+0.017*(tair-tt0))*4.186e2*1.e-5 ! convert to J/m/s/deg
es = estblf(tair)
qs = epsilo*es/(pres-(1.0_r8 - epsqs)*es)
dqsdt=latvap/(rh2o*tair*tair)*qs
alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1./(rair*tair))
gammaloc=(1+latvap/cpair*dqsdt)/(rhoair*qs)
! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
! should depend on mean radius of mode to account for gas kinetic effects
gloc=1./(rhoh2o/(diff0*rhoair*qs) &
+latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1.))
sqrtg=sqrt(gloc)
beta=4.*pi*rhoh2o*gloc*gammaloc
etafactor2max=1.e10/(alpha*wbar)**1.5 ! this should make eta big if na is very small.
do m=1,nmode
! internal mixture of aerosols
volc(m)=ma(m)/(rhodry(m)) ! only if variable size dist
if(volc(m).gt.1.e-39.and.na(m).gt.1.e-39)then
etafactor2(m)=1./(na(m)*beta*sqrtg) !fixed or variable size dist
! number mode radius (m)
amcubeloc(m)=(3.*volc(m)/(4.*pi*exp45logsig(m)*na(m))) ! only if variable size dist
smc(m)=smcrit(m) ! only for prescribed size dist
if(hygro(m).gt.1.e-10)then ! loop only if variable size dist
smc(m)=2.*aten*sqrt(aten/(27.*hygro(m)*amcubeloc(m)))
else
smc(m)=100.
endif
else
smc(m)=1.
etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
endif
lnsmloc(m)=log(smc(m)) ! only if variable size dist
enddo
! single updraft
wnuc=wbar
w=wbar
alw=alpha*wnuc
sqrtalw=sqrt(alw)
zeta=2.*sqrtalw*aten/(3.*sqrtg)
etafactor1=2.*alw*sqrtalw
do m=1,nmode
eta(m)=etafactor1*etafactor2(m)
enddo
call maxsat(zeta,eta,nmode,smc,smax)
lnsmax=log(smax)
xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
nact=0._r8
do m=1,nmode
x=2*(lnsmloc(m)-lnsmax)/(3*sq2*alogsig(m))
nact=nact+0.5*(1.-erf(x))*na(m)
enddo
nact=nact/rhoair ! convert from #/m3 to #/kg
deallocate( &
volc, &
eta, &
smc, &
etafactor2, &
amcubeloc, &
lnsmloc )
return
end subroutine activate
subroutine maxsat(zeta,eta,nmode,smc,smax) 3
! calculates maximum supersaturation for multiple
! competing aerosol modes.
! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
implicit none
! integer pmode
! parameter (pmode=naer_all)
integer nmode ! number of modes
real(r8) :: smc(:) ! critical supersaturation for number mode radius
real(r8) zeta
real(r8) :: eta(:)
real(r8) smax ! maximum supersaturation
integer m ! mode index
real(r8) sum, g1, g2
do m=1,nmode
if(zeta.gt.1.e5*eta(m).or.smc(m)*smc(m).gt.1.e5*eta(m))then
! weak forcing. essentially none activated
smax=1.e-20
else
! significant activation of this mode. calc activation all modes.
go to 1
endif
enddo
return
1 continue
sum=0
do m=1,nmode
if(eta(m).gt.1.e-20)then
g1=sqrt(zeta/eta(m))
g1=g1*g1*g1
g2=smc(m)/sqrt(eta(m)+3*zeta)
g2=sqrt(g2)
g2=g2*g2*g2
sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
else
sum=1.e20
endif
enddo
smax=1./sqrt(sum)
return
end subroutine maxsat
subroutine nucleati(wbar, tair, relhum, cldn, qc, nfice, rhoair, & 1,11
#ifdef MODAL_AERO
qaerpt, dgnum, ptype, naer_all, nuci, onihf, oniimm, onidep, onimey)
#else
na, ptype, naer_all, nuci, onihf, oniimm, onidep, onimey)
#endif
!---------------------------------------------------------------
! Purpose:
! The parameterization of ice nucleation.
!
! Method: The current method is based on Liu & Penner (2005)
! It related the ice nucleation with the aerosol number, temperature and the
! updraft velocity. It includes homogeneous freezing of sulfate, immersion
! freezing of soot, and Meyers et al. (1992) deposition nucleation
!
! Authors: Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010
!----------------------------------------------------------------
#ifdef MODAL_AERO
use modal_aero_data, only: numptr_amode, modeptr_accum, modeptr_coarse, modeptr_aitken, &
maxd_amode, sigmag_amode, &
lptr_dust_a_amode,lptr_nacl_a_amode
use constituents, only: pcnst
#endif
!-----------------------------------------------------
! Input Arguments
!
integer ptype, naer_all
real(r8), intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
real(r8), intent(in) :: tair ! temperature (K)
real(r8), intent(in) :: relhum ! relative humidity with respective to liquid
real(r8), intent(in) :: cldn ! new value of cloud fraction (fraction)
real(r8), intent(in) :: qc ! liquid water mixing ratio (kg/kg)
real(r8), intent(in) :: nfice ! ice mass fraction
real(r8), intent(in) :: rhoair ! air density (kg/m3)
#ifdef MODAL_AERO
real(r8), intent(in) :: qaerpt(pcnst) ! aerosol number and mass mixing ratios
real(r8), intent(in) :: dgnum(maxd_amode) ! aerosol mode dry diameter (m)
#else
real(r8), intent(in) :: na(naer_all) ! aerosol number concentration (/m3)
#endif
!
! Output Arguments
!
real(r8), intent(out) :: nuci ! ice number nucleated (#/kg)
real(r8), intent(out) :: onihf ! nucleated number from homogeneous freezing of so4
real(r8), intent(out) :: oniimm ! nucleated number from immersion freezing
real(r8), intent(out) :: onidep ! nucleated number from deposition nucleation
real(r8), intent(out) :: onimey ! nucleated number from deposition nucleation (meyers: mixed phase)
!
! Local workspace
!
real(r8) so4_num ! so4 aerosol number (#/cm^3)
real(r8) soot_num ! soot (hydrophilic) aerosol number (#/cm^3)
real(r8) dst1_num,dst2_num,dst3_num,dst4_num ! dust aerosol number (#/cm^3)
real(r8) dst_num ! total dust aerosol number (#/cm^3)
real(r8) nihf ! nucleated number from homogeneous freezing of so4
real(r8) niimm ! nucleated number from immersion freezing
real(r8) nidep ! nucleated number from deposition nucleation
real(r8) nimey ! nucleated number from deposition nucleation (meyers)
real(r8) n1,ni ! nucleated number
real(r8) tc,A,B,C,regm,RHw ! work variable
real(r8) esl,esi,deles ! work variable
real(r8) dst_scale
real(r8) subgrid
real(r8) dmc,ssmc ! variables for modal scheme.
so4_num=0.0_r8
soot_num=0.0_r8
dst_num=0.0_r8
dst1_num = 0.0_r8
dst2_num = 0.0_r8
dst3_num = 0.0_r8
dst4_num = 0.0_r8
!For modal aerosols, assume for the upper troposphere:
! soot = accumulation mode
! sulfate = aiken mode
! dust = coarse mode
! since modal has internal mixtures.
#ifdef MODAL_AERO
soot_num = qaerpt(numptr_amode(modeptr_accum))*1.0e-6_r8 !#/cm^3
dmc= qaerpt(lptr_dust_a_amode(modeptr_coarse))
ssmc=qaerpt(lptr_nacl_a_amode(modeptr_coarse))
if (dmc.gt.0._r8) then
dst_num=dmc/(ssmc+dmc) * qaerpt(numptr_amode(modeptr_coarse))*1.0e-6_r8 !#/cm^3
else
dst_num=0.0_r8
endif
if (dgnum(modeptr_aitken) .gt. 0._r8 ) then
so4_num = qaerpt(numptr_amode(modeptr_aitken))*1.0e-6_r8 & !#/cm^3, only allow so4 with D>0.1 um in ice nucleation
* (0.5_r8 - 0.5_r8*erf(log(0.1e-6_r8/dgnum(modeptr_aitken))/ &
(2._r8**0.5_r8*log(sigmag_amode(modeptr_aitken)))))
else
so4_num = 0.0_r8
endif
so4_num = max(0.0_r8,so4_num)
#else
if(idxsul .gt. 0) then
so4_num=na(idxsul)*1.0e-6_r8 ! #/cm^3
end if
!continue above philosophy here....
if(idxbcphi .gt. 0) then
soot_num=na(idxbcphi)*1.0e-6_r8 !#/cm^3
end if
if(idxdst1 .gt. 0) then
dst1_num=na(idxdst1) *1.0e-6_r8 !#/cm^3
end if
if(idxdst2 .gt. 0) then
dst2_num=na(idxdst2)*1.0e-6_r8 !#/cm^3
end if
if(idxdst3 .gt. 0) then
dst3_num=na(idxdst3)*1.0e-6_r8 !#/cm^3
end if
if(idxdst4 .gt. 0) then
dst4_num=na(idxdst4)*1.0e-6_r8 !#/cm^3
end if
dst_num =dst1_num+dst2_num+dst3_num+dst4_num
#endif
! no soot nucleation for now.
soot_num=0.0_r8
ni=0._r8
tc=tair-273.15_r8
! initialize
niimm=0._r8
nidep=0._r8
nihf=0._r8
if(so4_num.ge.1.0e-10_r8 .and. (soot_num+dst_num).ge.1.0e-10_r8 .and. cldn.gt.0._r8) then
!-----------------------------
! RHw parameterization for heterogeneous immersion nucleation
A = 0.0073_r8
B = 1.477_r8
C = 131.74_r8
RHw=(A*tc*tc+B*tc+C)*0.01_r8 ! RHi ~ 120-130%
subgrid = 1.2_r8
if((tc.le.-35.0_r8) .and. ((relhum*polysvp(tair,0)/polysvp(tair,1)*subgrid).ge.1.2_r8)) then ! use higher RHi threshold
A = -1.4938_r8 * log(soot_num+dst_num) + 12.884_r8
B = -10.41_r8 * log(soot_num+dst_num) - 67.69_r8
regm = A * log(wbar) + B
if(tc.gt.regm) then ! heterogeneous nucleation only
if(tc.lt.-40._r8 .and. wbar.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
niimm=0._r8
nidep=0._r8
n1=nihf
else
call hetero(tc,wbar,soot_num+dst_num,niimm,nidep)
nihf=0._r8
n1=niimm+nidep
endif
elseif (tc.lt.regm-5._r8) then ! homogeneous nucleation only
call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
niimm=0._r8
nidep=0._r8
n1=nihf
else ! transition between homogeneous and heterogeneous: interpolate in-between
if(tc.lt.-40._r8 .and. wbar.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
call hf(tc,wbar,relhum,subgrid,so4_num,nihf)
niimm=0._r8
nidep=0._r8
n1=nihf
else
call hf(regm-5._r8,wbar,relhum,subgrid,so4_num,nihf)
call hetero(regm,wbar,soot_num+dst_num,niimm,nidep)
if(nihf.le.(niimm+nidep)) then
n1=nihf
else
n1=(niimm+nidep)*((niimm+nidep)/nihf)**((tc-regm)/5._r8)
endif
endif
endif
ni=n1
endif
endif
1100 continue
! deposition/condensation nucleation in mixed clouds (-40<T<0C) (Meyers, 1992)
if(tc.lt.0._r8 .and. tc.gt.-37._r8 .and. qc.gt.1.e-12_r8) then
esl = polysvp(tair,0) ! over water in mixed clouds
esi = polysvp(tair,1) ! over ice
deles = (esl - esi)
nimey=1.e-3_r8*exp(12.96_r8*deles/esi - 0.639_r8)
else
nimey=0._r8
endif
nuci=ni+nimey
if(nuci.gt.9999._r8.or.nuci.lt.0._r8) then
write(iulog, *) 'incorrect ice nucleation number'
write(iulog, *) ni, tair, relhum, wbar, nihf, niimm, nidep,deles,esi,dst2_num,dst3_num,dst4_num
nuci=0._r8
endif
nuci=nuci*1.e+6_r8/rhoair ! change unit from #/cm3 to #/kg
onimey=nimey*1.e+6_r8/rhoair
onidep=nidep*1.e+6_r8/rhoair
oniimm=niimm*1.e+6_r8/rhoair
onihf=nihf*1.e+6_r8/rhoair
return
end subroutine nucleati
subroutine hetero(T,ww,Ns,Nis,Nid) 2
real(r8), intent(in) :: T, ww, Ns
real(r8), intent(out) :: Nis, Nid
real(r8) A11,A12,A21,A22,B11,B12,B21,B22
real(r8) A,B,C
!---------------------------------------------------------------------
! parameters
A11 = 0.0263_r8
A12 = -0.0185_r8
A21 = 2.758_r8
A22 = 1.3221_r8
B11 = -0.008_r8
B12 = -0.0468_r8
B21 = -0.2667_r8
B22 = -1.4588_r8
! ice from immersion nucleation (cm-3)
B = (A11+B11*log(Ns)) * log(ww) + (A12+B12*log(Ns))
C = A21+B21*log(Ns)
Nis = exp(A22) * Ns**B22 * exp(B*T) * ww**C
Nis = min(Nis,Ns)
Nid = 0.0_r8 ! don't include deposition nucleation for cirrus clouds when T<-37C
return
end subroutine hetero
subroutine hf(T,ww,RH,subgrid,Na,Ni) 4
real(r8), intent(in) :: T, ww, RH, subgrid, Na
real(r8), intent(out) :: Ni
real(r8) A1_fast,A21_fast,A22_fast,B1_fast,B21_fast,B22_fast
real(r8) A2_fast,B2_fast
real(r8) C1_fast,C2_fast,k1_fast,k2_fast
real(r8) A1_slow,A2_slow,B1_slow,B2_slow,B3_slow
real(r8) C1_slow,C2_slow,k1_slow,k2_slow
real(r8) regm
real(r8) A,B,C
real(r8) RHw
!---------------------------------------------------------------------
! parameters
A1_fast =0.0231_r8
A21_fast =-1.6387_r8 !(T>-64 deg)
A22_fast =-6.045_r8 !(T<=-64 deg)
B1_fast =-0.008_r8
B21_fast =-0.042_r8 !(T>-64 deg)
B22_fast =-0.112_r8 !(T<=-64 deg)
C1_fast =0.0739_r8
C2_fast =1.2372_r8
A1_slow =-0.3949_r8
A2_slow =1.282_r8
B1_slow =-0.0156_r8
B2_slow =0.0111_r8
B3_slow =0.0217_r8
C1_slow =0.120_r8
C2_slow =2.312_r8
Ni = 0.0_r8
!----------------------------
!RHw xiaohong's parameter
A = 6.0e-4_r8*log(ww)+6.6e-3_r8
B = 6.0e-2_r8*log(ww)+1.052_r8
C = 1.68_r8 *log(ww)+129.35_r8
RHw=(A*T*T+B*T+C)*0.01_r8
if((T.le.-37.0_r8) .and. ((RH*subgrid).ge.RHw)) then
regm = 6.07_r8*log(ww)-55.0_r8
if(T.ge.regm) then ! fast-growth regime
if(T.gt.-64.0_r8) then
A2_fast=A21_fast
B2_fast=B21_fast
else
A2_fast=A22_fast
B2_fast=B22_fast
endif
k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww))
k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww)
Ni = k1_fast*Na**(k2_fast)
Ni = min(Ni,Na)
else ! slow-growth regime
k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww))
k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww)
Ni = k1_slow*Na**(k2_slow)
Ni = min(Ni,Na)
endif
end if
return
end subroutine hf
!--xl
end module cldwat2m_micro