! =======================================================================
! cloud micro - physics package for gfdl global cloud resolving model
! the algorithms are originally derived from lin et al 1983. most of the
! key elements have been simplified / improved. this code at this stage
! bears little to no similarity to the original lin mp in zetac.
! therefore, it is best to be called gfdl micro - physics (gfdl mp) .
! developer: shian-jiann lin, linjiong zhou
! =======================================================================
module gfdl_cloud_microphys_mod
! use mpp_mod, only: stdlog, mpp_pe, mpp_root_pe, mpp_clock_id, &
! mpp_clock_begin, mpp_clock_end, clock_routine, &
! input_nml_file
! use diag_manager_mod, only: register_diag_field, send_data
! use time_manager_mod, only: time_type, get_time
! use constants_mod, only: grav, rdgas, rvgas, cp_air, hlv, hlf, pi => pi_8
! use fms_mod, only: write_version_number, open_namelist_file, &
! check_nml_error, file_exist, close_file
implicit none
private
public gfdl_cloud_microphys_driver, gfdl_cloud_microphys_init, gfdl_cloud_microphys_end
public wqs1, wqs2, qs_blend, wqsat_moist, wqsat2_moist
public qsmith_init, qsmith, es2_table1d, es3_table1d, esw_table1d
public setup_con, wet_bulb
public cloud_diagnosis
real :: missing_value = - 1.e10
logical :: module_is_initialized = .false.
logical :: qsmith_tables_initialized = .false.
character (len = 17) :: mod_name = 'gfdl_cloud_microphys'
real, parameter :: grav = 9.80665 ! gfs: acceleration due to gravity
real, parameter :: rdgas = 287.05 ! gfs: gas constant for dry air
real, parameter :: rvgas = 461.50 ! gfs: gas constant for water vapor
real, parameter :: cp_air = 1004.6 ! gfs: heat capacity of dry air at constant pressure
real, parameter :: hlv = 2.5e6 ! gfs: latent heat of evaporation
real, parameter :: hlf = 3.3358e5 ! gfs: latent heat of fusion
real, parameter :: pi = 3.1415926535897931 ! gfs: ratio of circle circumference to diameter
! real, parameter :: rdgas = 287.04 ! gfdl: gas constant for dry air
! real, parameter :: cp_air = rdgas * 7. / 2. ! 1004.675, heat capacity of dry air at constant pressure
real, parameter :: cp_vap = 4.0 * rvgas ! 1846.0, heat capacity of water vapore at constnat pressure
! real, parameter :: cv_air = 717.56 ! satoh value
real, parameter :: cv_air = cp_air - rdgas ! 717.55, heat capacity of dry air at constant volume
! real, parameter :: cv_vap = 1410.0 ! emanuel value
real, parameter :: cv_vap = 3.0 * rvgas ! 1384.5, heat capacity of water vapor at constant volume
! the following two are from emanuel's book "atmospheric convection"
! real, parameter :: c_ice = 2106.0 ! heat capacity of ice at 0 deg c: c = c_ice + 7.3 * (t - tice)
! real, parameter :: c_liq = 4190.0 ! heat capacity of water at 0 deg c
real, parameter :: c_ice = 1972.0 ! gfdl: heat capacity of ice at - 15 deg c
real, parameter :: c_liq = 4185.5 ! gfdl: heat capacity of water at 15 deg c
! real, parameter :: c_liq = 4218.0 ! ifs: heat capacity of liquid at 0 deg c
real, parameter :: eps = rdgas / rvgas ! 0.6219934995
real, parameter :: zvir = rvgas / rdgas - 1. ! 0.6077338443
real, parameter :: t_ice = 273.16 ! freezing temperature
real, parameter :: table_ice = 273.16 ! freezing point for qs table
! real, parameter :: e00 = 610.71 ! gfdl: saturation vapor pressure at 0 deg c
real, parameter :: e00 = 611.21 ! ifs: saturation vapor pressure at 0 deg c
real, parameter :: dc_vap = cp_vap - c_liq ! - 2339.5, isobaric heating / cooling
real, parameter :: dc_ice = c_liq - c_ice ! 2213.5, isobaric heating / colling
real, parameter :: hlv0 = hlv ! gfs: evaporation latent heat coefficient at 0 deg c
! real, parameter :: hlv0 = 2.501e6 ! emanuel appendix - 2
real, parameter :: hlf0 = hlf ! gfs: fussion latent heat coefficient at 0 deg c
! real, parameter :: hlf0 = 3.337e5 ! emanuel
real, parameter :: lv0 = hlv0 - dc_vap * t_ice! 3.13905782e6, evaporation latent heat coefficient at 0 deg k
real, parameter :: li00 = hlf0 - dc_ice * t_ice! - 2.7105966e5, fussion latend heat coefficient at 0 deg k
real, parameter :: d2ice = dc_vap + dc_ice ! - 126, isobaric heating / cooling
real, parameter :: li2 = lv0 + li00 ! 2.86799816e6, sublimation latent heat coefficient at 0 deg k
real, parameter :: qrmin = 1.e-8 ! min value for ???
real, parameter :: qvmin = 1.e-20 ! min value for water vapor (treated as zero)
real, parameter :: qcmin = 1.e-12 ! min value for cloud condensates
real, parameter :: vr_min = 1.e-3 ! min fall speed for rain
real, parameter :: vf_min = 1.e-5 ! min fall speed for cloud ice, snow, graupel
real, parameter :: dz_min = 1.e-2 ! use for correcting flipped height
real, parameter :: sfcrho = 1.2 ! surface air density
real, parameter :: rhor = 1.e3 ! density of rain water, lin83
real :: cracs, csacr, cgacr, cgacs, csacw, craci, csaci, cgacw, cgaci, cracw ! constants for accretions
real :: acco (3, 4) ! constants for accretions
real :: cssub (5), cgsub (5), crevp (5), cgfr (2), csmlt (5), cgmlt (5)
real :: es0, ces0
real :: pie, rgrav, fac_rc
real :: c_air, c_vap
real :: lati, latv, lats, lat2, lcp, icp, tcp ! used in bigg mechanism and wet bulk
real :: d0_vap ! the same as dc_vap, except that cp_vap can be cp_vap or cv_vap
real :: lv00 ! the same as lv0, except that cp_vap can be cp_vap or cv_vap
! cloud microphysics switchers
integer :: icloud_f = 0 ! cloud scheme
integer :: irain_f = 0 ! cloud water to rain auto conversion scheme
logical :: de_ice = .false. ! to prevent excessive build - up of cloud ice from external sources
logical :: sedi_transport = .true. ! transport of momentum in sedimentation
logical :: do_sedi_w = .false. ! transport of vertical motion in sedimentation
logical :: do_sedi_heat = .true. ! transport of heat in sedimentation
logical :: prog_ccn = .false. ! do prognostic ccn (yi ming's method)
logical :: do_qa = .true. ! do inline cloud fraction
logical :: rad_snow = .true. ! consider snow in cloud fraciton calculation
logical :: rad_graupel = .true. ! consider graupel in cloud fraction calculation
logical :: rad_rain = .true. ! consider rain in cloud fraction calculation
logical :: fix_negative = .false. ! fix negative water species
logical :: do_setup = .true. ! setup constants and parameters
logical :: p_nonhydro = .false. ! perform hydrosatic adjustment on air density
real, allocatable :: table (:), table2 (:), table3 (:), tablew (:)
real, allocatable :: des (:), des2 (:), des3 (:), desw (:)
logical :: tables_are_initialized = .false.
! logical :: master
! integer :: id_rh, id_vtr, id_vts, id_vtg, id_vti, id_rain, id_snow, id_graupel, &
! id_ice, id_prec, id_cond, id_var, id_droplets
! integer :: gfdl_mp_clock ! clock for timing of driver routine
real, parameter :: dt_fr = 8. ! homogeneous freezing of all cloud water at t_wfr - dt_fr
! minimum temperature water can exist (moore & molinero nov. 2011, nature)
! dt_fr can be considered as the error bar
real :: p_min = 100. ! minimum pressure (pascal) for mp to operate
! slj, the following parameters are for cloud - resolving resolution: 1 - 5 km
! qi0_crt = 0.8e-4
! qs0_crt = 0.6e-3
! c_psaci = 0.1
! c_pgacs = 0.1
! -----------------------------------------------------------------------
! namelist parameters
! -----------------------------------------------------------------------
real :: cld_min = 0.05 ! minimum cloud fraction
real :: tice = 273.16 ! set tice = 165. to trun off ice - phase phys (kessler emulator)
real :: t_min = 178. ! min temp to freeze - dry all water vapor
real :: t_sub = 184. ! min temp for sublimation of cloud ice
real :: mp_time = 150. ! maximum micro - physics time step (sec)
! relative humidity increment
real :: rh_inc = 0.25 ! rh increment for complete evaporation of cloud water and cloud ice
real :: rh_inr = 0.25 ! rh increment for minimum evaporation of rain
real :: rh_ins = 0.25 ! rh increment for sublimation of snow
! conversion time scale
real :: tau_r2g = 900. ! rain freezing during fast_sat
real :: tau_smlt = 900. ! snow melting
real :: tau_g2r = 600. ! graupel melting to rain
real :: tau_imlt = 600. ! cloud ice melting
real :: tau_i2s = 1000. ! cloud ice to snow auto - conversion
real :: tau_l2r = 900. ! cloud water to rain auto - conversion
real :: tau_v2l = 150. ! water vapor to cloud water (condensation)
real :: tau_l2v = 300. ! cloud water to water vapor (evaporation)
real :: tau_g2v = 900. ! grapuel sublimation
real :: tau_v2g = 21600. ! grapuel deposition -- make it a slow process
! horizontal subgrid variability
real :: dw_land = 0.20 ! base value for subgrid deviation / variability over land
real :: dw_ocean = 0.10 ! base value for ocean
! prescribed ccn
real :: ccn_o = 90. ! ccn over ocean (cm^ - 3)
real :: ccn_l = 270. ! ccn over land (cm^ - 3)
real :: rthresh = 10.0e-6 ! critical cloud drop radius (micro m)
! -----------------------------------------------------------------------
! wrf / wsm6 scheme: qi_gen = 4.92e-11 * (1.e3 * exp (0.1 * tmp)) ** 1.33
! optimized: qi_gen = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp)))
! qi_gen ~ 4.808e-7 at 0 c; 1.818e-6 at - 10 c, 9.82679e-5 at - 40c
! the following value is constructed such that qc_crt = 0 at zero c and @ - 10c matches
! wrf / wsm6 ice initiation scheme; qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den
! -----------------------------------------------------------------------
real :: sat_adj0 = 0.90 ! adjustment factor (0: no, 1: full) during fast_sat_adj
real :: qc_crt = 5.0e-8 ! mini condensate mixing ratio to allow partial cloudiness
real :: qi_lim = 1. ! cloud ice limiter to prevent large ice build up
real :: ql_mlt = 2.0e-3 ! max value of cloud water allowed from melted cloud ice
real :: qs_mlt = 1.0e-6 ! max cloud water due to snow melt
real :: ql_gen = 1.0e-3 ! max cloud water generation during remapping step if fast_sat_adj = .t.
real :: qi_gen = 1.82e-6 ! max cloud ice generation during remapping step
! cloud condensate upper bounds: "safety valves" for ql & qi
real :: ql0_max = 2.0e-3 ! max cloud water value (auto converted to rain)
real :: qi0_max = 1.0e-4 ! max cloud ice value (by other sources)
real :: qi0_crt = 1.0e-4 ! cloud ice to snow autoconversion threshold (was 1.e-4)
! qi0_crt is highly dependent on horizontal resolution
real :: qr0_crt = 1.0e-4 ! rain to snow or graupel / hail threshold
! lfo used * mixing ratio * = 1.e-4 (hail in lfo)
real :: qs0_crt = 1.0e-3 ! snow to graupel density threshold (0.6e-3 in purdue lin scheme)
real :: c_paut = 0.55 ! autoconversion cloud water to rain (use 0.5 to reduce autoconversion)
real :: c_psaci = 0.02 ! accretion: cloud ice to snow (was 0.1 in zetac)
real :: c_piacr = 5.0 ! accretion: rain to ice:
real :: c_cracw = 0.9 ! rain accretion efficiency
real :: c_pgacs = 2.0e-3 ! snow to graupel "accretion" eff. (was 0.1 in zetac)
! decreasing clin to reduce csacw (so as to reduce cloud water --- > snow)
real :: alin = 842.0 ! "a" in lin1983
real :: clin = 4.8 ! "c" in lin 1983, 4.8 -- > 6. (to ehance ql -- > qs)
! fall velocity tuning constants:
logical :: const_vi = .false. ! if .t. the constants are specified by v * _fac
logical :: const_vs = .false. ! if .t. the constants are specified by v * _fac
logical :: const_vg = .false. ! if .t. the constants are specified by v * _fac
logical :: const_vr = .false. ! if .t. the constants are specified by v * _fac
! good values:
real :: vi_fac = 1. ! if const_vi: 1 / 3
real :: vs_fac = 1. ! if const_vs: 1.
real :: vg_fac = 1. ! if const_vg: 2.
real :: vr_fac = 1. ! if const_vr: 4.
! upper bounds of fall speed (with variable speed option)
real :: vi_max = 0.5 ! max fall speed for ice
real :: vs_max = 5.0 ! max fall speed for snow
real :: vg_max = 8.0 ! max fall speed for graupel
real :: vr_max = 12. ! max fall speed for rain
! cloud microphysics switchers
logical :: fast_sat_adj = .false. ! has fast saturation adjustments
logical :: z_slope_liq = .true. ! use linear mono slope for autocconversions
logical :: z_slope_ice = .false. ! use linear mono slope for autocconversions
logical :: use_ccn = .false. ! must be true when prog_ccn is false
logical :: use_ppm = .false. ! use ppm fall scheme
logical :: mono_prof = .true. ! perform terminal fall with mono ppm scheme
logical :: mp_print = .false. ! cloud microphysics debugging printout
! real :: global_area = - 1.
real :: log_10, tice0, t_wfr
! -----------------------------------------------------------------------
! namelist
! -----------------------------------------------------------------------
namelist / gfdl_cloud_microphysics_nml / &
mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, mp_print
public &
mp_time, t_min, t_sub, tau_r2g, tau_smlt, tau_g2r, dw_land, dw_ocean, &
vi_fac, vr_fac, vs_fac, vg_fac, ql_mlt, do_qa, fix_negative, vi_max, &
vs_max, vg_max, vr_max, qs_mlt, qs0_crt, qi_gen, ql0_max, qi0_max, &
qi0_crt, qr0_crt, fast_sat_adj, rh_inc, rh_ins, rh_inr, const_vi, &
const_vs, const_vg, const_vr, use_ccn, rthresh, ccn_l, ccn_o, qc_crt, &
tau_g2v, tau_v2g, sat_adj0, c_piacr, tau_imlt, tau_v2l, tau_l2v, &
tau_i2s, tau_l2r, qi_lim, ql_gen, c_paut, c_psaci, c_pgacs, &
z_slope_liq, z_slope_ice, prog_ccn, c_cracw, alin, clin, tice, &
rad_snow, rad_graupel, rad_rain, cld_min, use_ppm, mono_prof, &
do_sedi_heat, sedi_transport, do_sedi_w, de_ice, icloud_f, irain_f, mp_print
contains
! -----------------------------------------------------------------------
! the driver of the gfdl cloud microphysics
! -----------------------------------------------------------------------
!subroutine gfdl_cloud_microphys_driver (qv, ql, qr, qi, qs, qg, qa, qn, &
! qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, &
! pt_dt, pt, w, uin, vin, udt, vdt, dz, delp, area, dt_in, &
! land, rain, snow, ice, graupel, &
! hydrostatic, phys_hydrostatic, &
! iis, iie, jjs, jje, kks, kke, ktop, kbot, time)
subroutine gfdl_cloud_microphys_driver (qv, ql, qr, qi, qs, qg, qa, qn, &
qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, pt_dt, pt, w, &
uin, vin, udt, vdt, dz, delp, area, dt_in, land, rain, snow, ice, &
graupel, hydrostatic, phys_hydrostatic, iis, iie, jjs, jje, kks, &
kke, ktop, kbot, seconds)
implicit none
logical, intent (in) :: hydrostatic, phys_hydrostatic
integer, intent (in) :: iis, iie, jjs, jje ! physics window
integer, intent (in) :: kks, kke ! vertical dimension
integer, intent (in) :: ktop, kbot ! vertical compute domain
integer, intent (in) :: seconds
real, intent (in) :: dt_in ! physics time step
real, intent (in), dimension (:, :) :: area ! cell area
real, intent (in), dimension (:, :) :: land ! land fraction
real, intent (in), dimension (:, :, :) :: delp, dz, uin, vin
real, intent (in), dimension (:, :, :) :: pt, qv, ql, qr, qg, qa, qn
real, intent (inout), dimension (:, :, :) :: qi, qs
real, intent (inout), dimension (:, :, :) :: pt_dt, qa_dt, udt, vdt, w
real, intent (inout), dimension (:, :, :) :: qv_dt, ql_dt, qr_dt
real, intent (inout), dimension (:, :, :) :: qi_dt, qs_dt, qg_dt
real, intent (out), dimension (:, :) :: rain, snow, ice, graupel
! logical :: used
real :: mpdt, rdt, dts, convt, tot_prec
integer :: i, j, k
integer :: is, ie, js, je ! physics window
integer :: ks, ke ! vertical dimension
integer :: days, ntimes
real, dimension (iie - iis + 1, jje - jjs + 1) :: prec_mp, prec1, cond, w_var, rh0
real, dimension (iie - iis + 1, jje - jjs + 1, kke - kks + 1) :: vt_r, vt_s, vt_g, vt_i, qn2
real, dimension (size (pt, 1), size (pt, 3)) :: m2_rain, m2_sol
real :: allmax
is = 1
js = 1
ks = 1
ie = iie - iis + 1
je = jje - jjs + 1
ke = kke - kks + 1
! call mpp_clock_begin (gfdl_mp_clock)
! -----------------------------------------------------------------------
! define heat capacity of dry air and water vapor based on hydrostatical property
! -----------------------------------------------------------------------
if (phys_hydrostatic .or. hydrostatic) then
c_air = cp_air
c_vap = cp_vap
p_nonhydro = .false.
else
c_air = cv_air
c_vap = cv_vap
p_nonhydro = .true.
endif
d0_vap = c_vap - c_liq
lv00 = hlv0 - d0_vap * t_ice
if (hydrostatic) do_sedi_w = .false.
! -----------------------------------------------------------------------
! define latent heat coefficient used in wet bulb and bigg mechanism
! -----------------------------------------------------------------------
latv = hlv
lati = hlf
lats = latv + lati
lat2 = lats * lats
lcp = latv / cp_air
icp = lati / cp_air
tcp = (latv + lati) / cp_air
! tendency zero out for am moist processes should be done outside the driver
! -----------------------------------------------------------------------
! define cloud microphysics sub time step
! -----------------------------------------------------------------------
mpdt = min (dt_in, mp_time)
rdt = 1. / dt_in
ntimes = nint (dt_in / mpdt)
! small time step:
dts = dt_in / real (ntimes)
! call get_time (time, seconds, days)
! -----------------------------------------------------------------------
! initialize precipitation
! -----------------------------------------------------------------------
do j = js, je
do i = is, ie
graupel (i, j) = 0.
rain (i, j) = 0.
snow (i, j) = 0.
ice (i, j) = 0.
cond (i, j) = 0.
enddo
enddo
! -----------------------------------------------------------------------
! major cloud microphysics
! -----------------------------------------------------------------------
do j = js, je
call mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, qg, &
qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
rain (:, j), snow (:, j), graupel (:, j), ice (:, j), m2_rain, &
m2_sol, cond (:, j), area (:, j), land (:, j), udt, vdt, pt_dt, &
qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, w_var, vt_r, &
vt_s, vt_g, vt_i, qn2)
enddo
! -----------------------------------------------------------------------
! no clouds allowed above ktop
! -----------------------------------------------------------------------
if (ks < ktop) then
do k = ks, ktop
if (do_qa) then
do j = js, je
do i = is, ie
qa_dt (i, j, k) = 0.
enddo
enddo
else
do j = js, je
do i = is, ie
! qa_dt (i, j, k) = - qa (i, j, k) * rdt
qa_dt (i, j, k) = 0. ! gfs
enddo
enddo
endif
enddo
endif
! -----------------------------------------------------------------------
! diagnostic output
! -----------------------------------------------------------------------
! if (id_vtr > 0) then
! used = send_data (id_vtr, vt_r, time, is_in = iis, js_in = jjs)
! endif
! if (id_vts > 0) then
! used = send_data (id_vts, vt_s, time, is_in = iis, js_in = jjs)
! endif
! if (id_vtg > 0) then
! used = send_data (id_vtg, vt_g, time, is_in = iis, js_in = jjs)
! endif
! if (id_vti > 0) then
! used = send_data (id_vti, vt_i, time, is_in = iis, js_in = jjs)
! endif
! if (id_droplets > 0) then
! used = send_data (id_droplets, qn2, time, is_in = iis, js_in = jjs)
! endif
! if (id_var > 0) then
! used = send_data (id_var, w_var, time, is_in = iis, js_in = jjs)
! endif
! convert to mm / day
convt = 86400. * rdt * rgrav
do j = js, je
do i = is, ie
rain (i, j) = rain (i, j) * convt
snow (i, j) = snow (i, j) * convt
ice (i, j) = ice (i, j) * convt
graupel (i, j) = graupel (i, j) * convt
prec_mp (i, j) = rain (i, j) + snow (i, j) + ice (i, j) + graupel (i, j)
enddo
enddo
! if (id_cond > 0) then
! do j = js, je
! do i = is, ie
! cond (i, j) = cond (i, j) * rgrav
! enddo
! enddo
! used = send_data (id_cond, cond, time, is_in = iis, js_in = jjs)
! endif
! if (id_snow > 0) then
! used = send_data (id_snow, snow, time, iis, jjs)
! used = send_data (id_snow, snow, time, is_in = iis, js_in = jjs)
! if (mp_print .and. seconds == 0) then
! tot_prec = g_sum (snow, is, ie, js, je, area, 1)
! if (master) write (*, *) 'mean snow = ', tot_prec
! endif
! endif
!
! if (id_graupel > 0) then
! used = send_data (id_graupel, graupel, time, iis, jjs)
! used = send_data (id_graupel, graupel, time, is_in = iis, js_in = jjs)
! if (mp_print .and. seconds == 0) then
! tot_prec = g_sum (graupel, is, ie, js, je, area, 1)
! if (master) write (*, *) 'mean graupel = ', tot_prec
! endif
! endif
!
! if (id_ice > 0) then
! used = send_data (id_ice, ice, time, iis, jjs)
! used = send_data (id_ice, ice, time, is_in = iis, js_in = jjs)
! if (mp_print .and. seconds == 0) then
! tot_prec = g_sum (ice, is, ie, js, je, area, 1)
! if (master) write (*, *) 'mean ice_mp = ', tot_prec
! endif
! endif
!
! if (id_rain > 0) then
! used = send_data (id_rain, rain, time, iis, jjs)
! used = send_data (id_rain, rain, time, is_in = iis, js_in = jjs)
! if (mp_print .and. seconds == 0) then
! tot_prec = g_sum (rain, is, ie, js, je, area, 1)
! if (master) write (*, *) 'mean rain = ', tot_prec
! endif
! endif
!
! if (id_rh > 0) then !not used?
! used = send_data (id_rh, rh0, time, iis, jjs)
! used = send_data (id_rh, rh0, time, is_in = iis, js_in = jjs)
! endif
!
!
! if (id_prec > 0) then
! used = send_data (id_prec, prec_mp, time, iis, jjs)
! used = send_data (id_prec, prec_mp, time, is_in = iis, js_in = jjs)
! endif
! if (mp_print) then
! prec1 (:, :) = prec1 (:, :) + prec_mp (:, :)
! if (seconds == 0) then
! prec1 (:, :) = prec1 (:, :) * dt_in / 86400.
! tot_prec = g_sum (prec1, is, ie, js, je, area, 1)
! if (master) write (*, *) 'daily prec_mp = ', tot_prec
! prec1 (:, :) = 0.
! endif
! endif
! call mpp_clock_end (gfdl_mp_clock)
end subroutine gfdl_cloud_microphys_driver
! -----------------------------------------------------------------------
! gfdl cloud microphysics, major program
! lin et al., 1983, jam, 1065 - 1092, and
! rutledge and hobbs, 1984, jas, 2949 - 2972
! terminal fall is handled lagrangianly by conservative fv algorithm
! pt: temperature (k)
! 6 water species:
! 1) qv: water vapor (kg / kg)
! 2) ql: cloud water (kg / kg)
! 3) qr: rain (kg / kg)
! 4) qi: cloud ice (kg / kg)
! 5) qs: snow (kg / kg)
! 6) qg: graupel (kg / kg)
! -----------------------------------------------------------------------
subroutine mpdrv (hydrostatic, uin, vin, w, delp, pt, qv, ql, qr, qi, qs, &
qg, qa, qn, dz, is, ie, js, je, ks, ke, ktop, kbot, j, dt_in, ntimes, &
rain, snow, graupel, ice, m2_rain, m2_sol, cond, area1, land, &
u_dt, v_dt, pt_dt, qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt, qa_dt, &
w_var, vt_r, vt_s, vt_g, vt_i, qn2)
implicit none
logical, intent (in) :: hydrostatic
integer, intent (in) :: j, is, ie, js, je, ks, ke
integer, intent (in) :: ntimes, ktop, kbot
real, intent (in) :: dt_in
real, intent (in), dimension (is:) :: area1, land
real, intent (in), dimension (is:, js:, ks:) :: uin, vin, delp, pt, dz
real, intent (in), dimension (is:, js:, ks:) :: qv, ql, qr, qg, qa, qn
real, intent (inout), dimension (is:, js:, ks:) :: qi, qs
real, intent (inout), dimension (is:, js:, ks:) :: u_dt, v_dt, w, pt_dt, qa_dt
real, intent (inout), dimension (is:, js:, ks:) :: qv_dt, ql_dt, qr_dt, qi_dt, qs_dt, qg_dt
real, intent (inout), dimension (is:) :: rain, snow, ice, graupel, cond
real, intent (out), dimension (is:, js:) :: w_var
real, intent (out), dimension (is:, js:, ks:) :: vt_r, vt_s, vt_g, vt_i, qn2
real, intent (out), dimension (is:, ks:) :: m2_rain, m2_sol
real, dimension (ktop:kbot) :: qvz, qlz, qrz, qiz, qsz, qgz, qaz
real, dimension (ktop:kbot) :: vtiz, vtsz, vtgz, vtrz
real, dimension (ktop:kbot) :: dp0, dp1, dz0, dz1
real, dimension (ktop:kbot) :: qv0, ql0, qr0, qi0, qs0, qg0, qa0
real, dimension (ktop:kbot) :: t0, den, den0, tz, p1, denfac
real, dimension (ktop:kbot) :: ccn, c_praut, m1_rain, m1_sol, m1
real, dimension (ktop:kbot) :: u0, v0, u1, v1, w1
real :: cpaut, rh_adj, rh_rain
real :: r1, s1, i1, g1, rdt, ccn0
real :: dt_rain, dts
real :: s_leng, t_land, t_ocean, h_var
real :: cvm, tmp, omq
real :: dqi, qio, qin
integer :: i, k, n
dts = dt_in / real (ntimes)
dt_rain = dts * 0.5
rdt = 1. / dt_in
! -----------------------------------------------------------------------
! use local variables
! -----------------------------------------------------------------------
do i = is, ie
do k = ktop, kbot
qiz (k) = qi (i, j, k)
qsz (k) = qs (i, j, k)
enddo
! -----------------------------------------------------------------------
! this is to prevent excessive build - up of cloud ice from external sources
! -----------------------------------------------------------------------
if (de_ice) then
do k = ktop, kbot
qio = qiz (k) - dt_in * qi_dt (i, j, k) ! original qi before phys
qin = max (qio, qi0_max) ! adjusted value
if (qiz (k) > qin) then
qsz (k) = qsz (k) + qiz (k) - qin
qiz (k) = qin
dqi = (qin - qio) * rdt ! modified qi tendency
qs_dt (i, j, k) = qs_dt (i, j, k) + qi_dt (i, j, k) - dqi
qi_dt (i, j, k) = dqi
qi (i, j, k) = qiz (k)
qs (i, j, k) = qsz (k)
endif
enddo
endif
do k = ktop, kbot
t0 (k) = pt (i, j, k)
tz (k) = t0 (k)
dp1 (k) = delp (i, j, k)
dp0 (k) = dp1 (k) ! moist air mass * grav
! -----------------------------------------------------------------------
! convert moist mixing ratios to dry mixing ratios
! -----------------------------------------------------------------------
qvz (k) = qv (i, j, k)
qlz (k) = ql (i, j, k)
qrz (k) = qr (i, j, k)
qgz (k) = qg (i, j, k)
! dp1: dry air_mass
! dp1 (k) = dp1 (k) * (1. - (qvz (k) + qlz (k) + qrz (k) + qiz (k) + qsz (k) + qgz (k)))
dp1 (k) = dp1 (k) * (1. - qvz (k)) ! gfs
omq = dp0 (k) / dp1 (k)
qvz (k) = qvz (k) * omq
qlz (k) = qlz (k) * omq
qrz (k) = qrz (k) * omq
qiz (k) = qiz (k) * omq
qsz (k) = qsz (k) * omq
qgz (k) = qgz (k) * omq
qa0 (k) = qa (i, j, k)
qaz (k) = 0.
dz0 (k) = dz (i, j, k)
den0 (k) = - dp1 (k) / (grav * dz0 (k)) ! density of dry air
p1 (k) = den0 (k) * rdgas * t0 (k) ! dry air pressure
! -----------------------------------------------------------------------
! save a copy of old value for computing tendencies
! -----------------------------------------------------------------------
qv0 (k) = qvz (k)
ql0 (k) = qlz (k)
qr0 (k) = qrz (k)
qi0 (k) = qiz (k)
qs0 (k) = qsz (k)
qg0 (k) = qgz (k)
! -----------------------------------------------------------------------
! for sedi_momentum
! -----------------------------------------------------------------------
m1 (k) = 0.
u0 (k) = uin (i, j, k)
v0 (k) = vin (i, j, k)
u1 (k) = u0 (k)
v1 (k) = v0 (k)
enddo
if (do_sedi_w) then
do k = ktop, kbot
w1 (k) = w (i, j, k)
enddo
endif
! -----------------------------------------------------------------------
! calculate cloud condensation nuclei (ccn)
! the following is based on klein eq. 15
! -----------------------------------------------------------------------
cpaut = c_paut * 0.104 * grav / 1.717e-5
if (prog_ccn) then
do k = ktop, kbot
! convert # / cc to # / m^3
ccn (k) = qn (i, j, k) * 1.e6
c_praut (k) = cpaut * (ccn (k) * rhor) ** (- 1. / 3.)
enddo
use_ccn = .false.
else
ccn0 = (ccn_l * land (i) + ccn_o * (1. - land (i))) * 1.e6
if (use_ccn) then
! -----------------------------------------------------------------------
! ccn is formulted as ccn = ccn_surface * (den / den_surface)
! -----------------------------------------------------------------------
ccn0 = ccn0 * rdgas * tz (kbot) / p1 (kbot)
endif
tmp = cpaut * (ccn0 * rhor) ** (- 1. / 3.)
do k = ktop, kbot
c_praut (k) = tmp
ccn (k) = ccn0
enddo
endif
! -----------------------------------------------------------------------
! calculate horizontal subgrid variability
! total water subgrid deviation in horizontal direction
! default area dependent form: use dx ~ 100 km as the base
! -----------------------------------------------------------------------
s_leng = sqrt (sqrt (area1 (i) / 1.e10))
t_land = dw_land * s_leng
t_ocean = dw_ocean * s_leng
h_var = t_land * land (i) + t_ocean * (1. - land (i))
h_var = min (0.20, max (0.01, h_var))
! if (id_var > 0) w_var (i, j) = h_var
! -----------------------------------------------------------------------
! relative humidity increment
! -----------------------------------------------------------------------
rh_adj = 1. - h_var - rh_inc
rh_rain = max (0.35, rh_adj - rh_inr) ! rh_inr = 0.25
! -----------------------------------------------------------------------
! fix all negative water species
! -----------------------------------------------------------------------
if (fix_negative) &
call neg_adj (ktop, kbot, tz, dp1, qvz, qlz, qrz, qiz, qsz, qgz)
m2_rain (i, :) = 0.
m2_sol (i, :) = 0.
do n = 1, ntimes
! -----------------------------------------------------------------------
! define air density based on hydrostatical property
! -----------------------------------------------------------------------
if (p_nonhydro) then
do k = ktop, kbot
dz1 (k) = dz0 (k)
den (k) = den0 (k) ! dry air density remains the same
denfac (k) = sqrt (sfcrho / den (k))
enddo
else
do k = ktop, kbot
dz1 (k) = dz0 (k) * tz (k) / t0 (k) ! hydrostatic balance
den (k) = den0 (k) * dz0 (k) / dz1 (k)
denfac (k) = sqrt (sfcrho / den (k))
enddo
endif
! -----------------------------------------------------------------------
! time - split warm rain processes: 1st pass
! -----------------------------------------------------------------------
call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
rain (i) = rain (i) + r1
do k = ktop, kbot
m2_rain (i, k) = m2_rain (i, k) + m1_rain (k)
m1 (k) = m1 (k) + m1_rain (k)
enddo
! -----------------------------------------------------------------------
! sedimentation of cloud ice, snow, and graupel
! -----------------------------------------------------------------------
call fall_speed (ktop, kbot, den, qsz, qiz, qgz, qlz, tz, vtsz, vtiz, vtgz)
call terminal_fall (dts, ktop, kbot, tz, qvz, qlz, qrz, qgz, qsz, qiz, &
dz1, dp1, den, vtgz, vtsz, vtiz, r1, g1, s1, i1, m1_sol, w1)
rain (i) = rain (i) + r1 ! from melted snow & ice that reached the ground
snow (i) = snow (i) + s1
graupel (i) = graupel (i) + g1
ice (i) = ice (i) + i1
! -----------------------------------------------------------------------
! heat transportation during sedimentation
! -----------------------------------------------------------------------
if (do_sedi_heat) &
call sedi_heat (ktop, kbot, dp1, m1_sol, dz1, tz, qvz, qlz, qrz, qiz, &
qsz, qgz, c_ice)
! -----------------------------------------------------------------------
! time - split warm rain processes: 2nd pass
! -----------------------------------------------------------------------
call warm_rain (dt_rain, ktop, kbot, dp1, dz1, tz, qvz, qlz, qrz, qiz, qsz, &
qgz, den, denfac, ccn, c_praut, rh_rain, vtrz, r1, m1_rain, w1, h_var)
rain (i) = rain (i) + r1
do k = ktop, kbot
m2_rain (i, k) = m2_rain (i, k) + m1_rain (k)
m2_sol (i, k) = m2_sol (i, k) + m1_sol (k)
m1 (k) = m1 (k) + m1_rain (k) + m1_sol (k)
enddo
! -----------------------------------------------------------------------
! ice - phase microphysics
! -----------------------------------------------------------------------
call icloud (ktop, kbot, tz, p1, qvz, qlz, qrz, qiz, qsz, qgz, dp1, den, &
denfac, vtsz, vtgz, vtrz, qaz, rh_adj, rh_rain, dts, h_var)
enddo
! convert units from Pa*kg/kg to kg/m^2/s
m2_rain (i, :) = m2_rain (i, :) * rdt * rgrav
m2_sol (i, :) = m2_sol (i, :) * rdt * rgrav
! -----------------------------------------------------------------------
! momentum transportation during sedimentation
! note: dp1 is dry mass; dp0 is the old moist (total) mass
! -----------------------------------------------------------------------
if (sedi_transport) then
do k = ktop + 1, kbot
u1 (k) = (dp0 (k) * u1 (k) + m1 (k - 1) * u1 (k - 1)) / (dp0 (k) + m1 (k - 1))
v1 (k) = (dp0 (k) * v1 (k) + m1 (k - 1) * v1 (k - 1)) / (dp0 (k) + m1 (k - 1))
u_dt (i, j, k) = u_dt (i, j, k) + (u1 (k) - u0 (k)) * rdt
v_dt (i, j, k) = v_dt (i, j, k) + (v1 (k) - v0 (k)) * rdt
enddo
endif
if (do_sedi_w) then
do k = ktop, kbot
w (i, j, k) = w1 (k)
enddo
endif
! -----------------------------------------------------------------------
! update moist air mass (actually hydrostatic pressure)
! convert to dry mixing ratios
! -----------------------------------------------------------------------
do k = ktop, kbot
omq = dp1 (k) / dp0 (k)
qv_dt (i, j, k) = qv_dt (i, j, k) + rdt * (qvz (k) - qv0 (k)) * omq
ql_dt (i, j, k) = ql_dt (i, j, k) + rdt * (qlz (k) - ql0 (k)) * omq
qr_dt (i, j, k) = qr_dt (i, j, k) + rdt * (qrz (k) - qr0 (k)) * omq
qi_dt (i, j, k) = qi_dt (i, j, k) + rdt * (qiz (k) - qi0 (k)) * omq
qs_dt (i, j, k) = qs_dt (i, j, k) + rdt * (qsz (k) - qs0 (k)) * omq
qg_dt (i, j, k) = qg_dt (i, j, k) + rdt * (qgz (k) - qg0 (k)) * omq
cvm = c_air + qvz (k) * c_vap + (qrz (k) + qlz (k)) * c_liq + (qiz (k) + qsz (k) + qgz (k)) * c_ice
pt_dt (i, j, k) = pt_dt (i, j, k) + rdt * (tz (k) - t0 (k)) * cvm / cp_air
enddo
! -----------------------------------------------------------------------
! update cloud fraction tendency
! -----------------------------------------------------------------------
do k = ktop, kbot
if (do_qa) then
qa_dt (i, j, k) = 0.
else
qa_dt (i, j, k) = qa_dt (i, j, k) + rdt * (qaz (k) / real (ntimes) - qa0 (k))
endif
enddo
! -----------------------------------------------------------------------
! fms diagnostics:
! -----------------------------------------------------------------------
! if (id_cond > 0) then
! do k = ktop, kbot ! total condensate
! cond (i) = cond (i) + dp1 (k) * (qlz (k) + qrz (k) + qsz (k) + qiz (k) + qgz (k))
! enddo
! endif
!
! if (id_vtr > 0) then
! do k = ktop, kbot
! vt_r (i, j, k) = vtrz (k)
! enddo
! endif
!
! if (id_vts > 0) then
! do k = ktop, kbot
! vt_s (i, j, k) = vtsz (k)
! enddo
! endif
!
! if (id_vtg > 0) then
! do k = ktop, kbot
! vt_g (i, j, k) = vtgz (k)
! enddo
! endif
!
! if (id_vts > 0) then
! do k = ktop, kbot
! vt_i (i, j, k) = vtiz (k)
! enddo
! endif
!
! if (id_droplets > 0) then
! do k = ktop, kbot
! qn2 (i, j, k) = ccn (k)
! enddo
! endif
enddo
end subroutine mpdrv
! -----------------------------------------------------------------------
! sedimentation of heat
! -----------------------------------------------------------------------
subroutine sedi_heat (ktop, kbot, dm, m1, dz, tz, qv, ql, qr, qi, qs, qg, cw)
implicit none
! input q fields are dry mixing ratios, and dm is dry air mass
integer, intent (in) :: ktop, kbot
real, intent (in), dimension (ktop:kbot) :: dm, m1, dz, qv, ql, qr, qi, qs, qg
real, intent (inout), dimension (ktop:kbot) :: tz
real, intent (in) :: cw ! heat capacity
real, dimension (ktop:kbot) :: dgz, cvn
real :: tmp
integer :: k
do k = ktop, kbot
dgz (k) = - 0.5 * grav * dz (k) ! > 0
cvn (k) = dm (k) * (cv_air + qv (k) * cv_vap + (qr (k) + ql (k)) * &
c_liq + (qi (k) + qs (k) + qg (k)) * c_ice)
enddo
! -----------------------------------------------------------------------
! sjl, july 2014
! assumption: the ke in the falling condensates is negligible compared to the potential energy
! that was unaccounted for. local thermal equilibrium is assumed, and the loss in pe is transformed
! into internal energy (to heat the whole grid box)
! backward time - implicit upwind transport scheme:
! dm here is dry air mass
! -----------------------------------------------------------------------
k = ktop
tmp = cvn (k) + m1 (k) * cw
tz (k) = (tmp * tz (k) + m1 (k) * dgz (k)) / tmp
! -----------------------------------------------------------------------
! implicit algorithm: can't be vectorized
! needs an inner i - loop for vectorization
! -----------------------------------------------------------------------
do k = ktop + 1, kbot
tz (k) = ((cvn (k) + cw * (m1 (k) - m1 (k - 1))) * tz (k) + m1 (k - 1) * &
cw * tz (k - 1) + dgz (k) * (m1 (k - 1) + m1 (k))) / (cvn (k) + cw * m1 (k))
enddo
end subroutine sedi_heat
! -----------------------------------------------------------------------
! warm rain cloud microphysics
! -----------------------------------------------------------------------
subroutine warm_rain (dt, ktop, kbot, dp, dz, tz, qv, ql, qr, qi, qs, qg, &
den, denfac, ccn, c_praut, rh_rain, vtr, r1, m1_rain, w1, h_var)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: dt ! time step (s)
real, intent (in) :: rh_rain, h_var
real, intent (in), dimension (ktop:kbot) :: dp, dz, den
real, intent (in), dimension (ktop:kbot) :: denfac, ccn, c_praut
real, intent (inout), dimension (ktop:kbot) :: tz, vtr
real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qi, qs, qg
real, intent (inout), dimension (ktop:kbot) :: m1_rain, w1
real, intent (out) :: r1
real, parameter :: so3 = 7. / 3.
real, dimension (ktop:kbot) :: dl, dm
real, dimension (ktop:kbot + 1) :: ze, zt
real :: sink, dq, qc0, qc
real :: qden
real :: zs = 0.
real :: dt5
integer :: k
! fall velocity constants:
real, parameter :: vconr = 2503.23638966667
real, parameter :: normr = 25132741228.7183
real, parameter :: thr = 1.e-8
logical :: no_fall
dt5 = 0.5 * dt
! -----------------------------------------------------------------------
! terminal speed of rain
! -----------------------------------------------------------------------
m1_rain (:) = 0.
call check_column (ktop, kbot, qr, no_fall)
if (no_fall) then
vtr (:) = vf_min
r1 = 0.
else
! -----------------------------------------------------------------------
! fall speed of rain
! -----------------------------------------------------------------------
if (const_vr) then
vtr (:) = vr_fac ! ifs_2016: 4.0
else
do k = ktop, kbot
qden = qr (k) * den (k)
if (qr (k) < thr) then
vtr (k) = vr_min
else
vtr (k) = vr_fac * vconr * sqrt (min (10., sfcrho / den (k))) * &
exp (0.2 * log (qden / normr))
vtr (k) = min (vr_max, max (vr_min, vtr (k)))
endif
enddo
endif
ze (kbot + 1) = zs
do k = kbot, ktop, - 1
ze (k) = ze (k + 1) - dz (k) ! dz < 0
enddo
! -----------------------------------------------------------------------
! evaporation and accretion of rain for the first 1 / 2 time step
! -----------------------------------------------------------------------
! if (.not. fast_sat_adj) &
call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
if (do_sedi_w) then
do k = ktop, kbot
dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
enddo
endif
! -----------------------------------------------------------------------
! mass flux induced by falling rain
! -----------------------------------------------------------------------
if (use_ppm) then
zt (ktop) = ze (ktop)
do k = ktop + 1, kbot
zt (k) = ze (k) - dt5 * (vtr (k - 1) + vtr (k))
enddo
zt (kbot + 1) = zs - dt * vtr (kbot)
do k = ktop, kbot
if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
enddo
call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qr, r1, m1_rain, mono_prof)
else
call implicit_fall (dt, ktop, kbot, ze, vtr, dp, qr, r1, m1_rain)
endif
! -----------------------------------------------------------------------
! vertical velocity transportation during sedimentation
! -----------------------------------------------------------------------
if (do_sedi_w) then
w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_rain (ktop) * vtr (ktop)) / (dm (ktop) - m1_rain (ktop))
do k = ktop + 1, kbot
w1 (k) = (dm (k) * w1 (k) - m1_rain (k - 1) * vtr (k - 1) + m1_rain (k) * vtr (k)) &
/ (dm (k) + m1_rain (k - 1) - m1_rain (k))
enddo
endif
! -----------------------------------------------------------------------
! heat transportation during sedimentation
! -----------------------------------------------------------------------
if (do_sedi_heat) &
call sedi_heat (ktop, kbot, dp, m1_rain, dz, tz, qv, ql, qr, qi, qs, qg, c_liq)
! -----------------------------------------------------------------------
! evaporation and accretion of rain for the remaing 1 / 2 time step
! -----------------------------------------------------------------------
call revap_racc (ktop, kbot, dt5, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
endif
! -----------------------------------------------------------------------
! auto - conversion
! assuming linear subgrid vertical distribution of cloud water
! following lin et al. 1994, mwr
! -----------------------------------------------------------------------
if (irain_f /= 0) then
! -----------------------------------------------------------------------
! no subgrid varaibility
! -----------------------------------------------------------------------
do k = ktop, kbot
qc0 = fac_rc * ccn (k)
if (tz (k) > t_wfr) then
if (use_ccn) then
! -----------------------------------------------------------------------
! ccn is formulted as ccn = ccn_surface * (den / den_surface)
! -----------------------------------------------------------------------
qc = qc0
else
qc = qc0 / den (k)
endif
dq = ql (k) - qc
if (dq > 0.) then
sink = min (dq, dt * c_praut (k) * den (k) * exp (so3 * log (ql (k))))
ql (k) = ql (k) - sink
qr (k) = qr (k) + sink
endif
endif
enddo
else
! -----------------------------------------------------------------------
! with subgrid varaibility
! -----------------------------------------------------------------------
call linear_prof (kbot - ktop + 1, ql (ktop), dl (ktop), z_slope_liq, h_var)
do k = ktop, kbot
qc0 = fac_rc * ccn (k)
if (tz (k) > t_wfr + dt_fr) then
dl (k) = min (max (1.e-6, dl (k)), 0.5 * ql (k))
! --------------------------------------------------------------------
! as in klein's gfdl am2 stratiform scheme (with subgrid variations)
! --------------------------------------------------------------------
if (use_ccn) then
! --------------------------------------------------------------------
! ccn is formulted as ccn = ccn_surface * (den / den_surface)
! --------------------------------------------------------------------
qc = qc0
else
qc = qc0 / den (k)
endif
dq = 0.5 * (ql (k) + dl (k) - qc)
! --------------------------------------------------------------------
! dq = dl if qc == q_minus = ql - dl
! dq = 0 if qc == q_plus = ql + dl
! --------------------------------------------------------------------
if (dq > 0.) then ! q_plus > qc
! --------------------------------------------------------------------
! revised continuous form: linearly decays (with subgrid dl) to zero at qc == ql + dl
! --------------------------------------------------------------------
sink = min (1., dq / dl (k)) * dt * c_praut (k) * den (k) * exp (so3 * log (ql (k)))
ql (k) = ql (k) - sink
qr (k) = qr (k) + sink
endif
endif
enddo
endif
end subroutine warm_rain
! -----------------------------------------------------------------------
! evaporation of rain
! -----------------------------------------------------------------------
subroutine revap_racc (ktop, kbot, dt, tz, qv, ql, qr, qi, qs, qg, den, denfac, rh_rain, h_var)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: dt ! time step (s)
real, intent (in) :: rh_rain, h_var
real, intent (in), dimension (ktop:kbot) :: den, denfac
real, intent (inout), dimension (ktop:kbot) :: tz, qv, qr, ql, qi, qs, qg
real, dimension (ktop:kbot) :: lhl, cvm, q_liq, q_sol, lcpk
real :: dqv, qsat, dqsdt, evap, t2, qden, q_plus, q_minus, sink
real :: qpz, dq, dqh, tin
integer :: k
do k = ktop, kbot
if (tz (k) > t_wfr .and. qr (k) > qrmin) then
! -----------------------------------------------------------------------
! define heat capacity and latent heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
q_liq (k) = ql (k) + qr (k)
q_sol (k) = qi (k) + qs (k) + qg (k)
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
lcpk (k) = lhl (k) / cvm (k)
tin = tz (k) - lcpk (k) * ql (k) ! presence of clouds suppresses the rain evap
qpz = qv (k) + ql (k)
qsat = wqs2 (tin, den (k), dqsdt)
dqh = max (ql (k), h_var * max (qpz, qcmin))
dqh = min (dqh, 0.2 * qpz) ! new limiter
dqv = qsat - qv (k) ! use this to prevent super - sat the gird box
q_minus = qpz - dqh
q_plus = qpz + dqh
! -----------------------------------------------------------------------
! qsat must be > q_minus to activate evaporation
! qsat must be < q_plus to activate accretion
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! rain evaporation
! -----------------------------------------------------------------------
if (dqv > qvmin .and. qsat > q_minus) then
if (qsat > q_plus) then
dq = qsat - qpz
else
! -----------------------------------------------------------------------
! q_minus < qsat < q_plus
! dq == dqh if qsat == q_minus
! -----------------------------------------------------------------------
dq = 0.25 * (q_minus - qsat) ** 2 / dqh
endif
qden = qr (k) * den (k)
t2 = tin * tin
evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * &
exp (0.725 * log (qden))) / (crevp (4) * t2 + crevp (5) * qsat * den (k))
evap = min (qr (k), dt * evap, dqv / (1. + lcpk (k) * dqsdt))
! -----------------------------------------------------------------------
! alternative minimum evap in dry environmental air
! sink = min (qr (k), dim (rh_rain * qsat, qv (k)) / (1. + lcpk (k) * dqsdt))
! evap = max (evap, sink)
! -----------------------------------------------------------------------
qr (k) = qr (k) - evap
qv (k) = qv (k) + evap
q_liq (k) = q_liq (k) - evap
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) - evap * lhl (k) / cvm (k)
endif
! -----------------------------------------------------------------------
! accretion: pracc
! -----------------------------------------------------------------------
! if (qr (k) > qrmin .and. ql (k) > 1.e-7 .and. qsat < q_plus) then
if (qr (k) > qrmin .and. ql (k) > 1.e-6 .and. qsat < q_minus) then
sink = dt * denfac (k) * cracw * exp (0.95 * log (qr (k) * den (k)))
sink = sink / (1. + sink) * ql (k)
ql (k) = ql (k) - sink
qr (k) = qr (k) + sink
endif
endif ! warm - rain
enddo
end subroutine revap_racc
! -----------------------------------------------------------------------
! definition of vertical subgrid variability
! used for cloud ice and cloud water autoconversion
! qi -- > ql & ql -- > qr
! edges: qe == qbar + / - dm
! -----------------------------------------------------------------------
subroutine linear_prof (km, q, dm, z_var, h_var)
implicit none
integer, intent (in) :: km
real, intent (in) :: q (km), h_var
real, intent (out) :: dm (km)
logical, intent (in) :: z_var
real :: dq (km)
integer :: k
if (z_var) then
do k = 2, km
dq (k) = 0.5 * (q (k) - q (k - 1))
enddo
dm (1) = 0.
! -----------------------------------------------------------------------
! use twice the strength of the positive definiteness limiter (lin et al 1994)
! -----------------------------------------------------------------------
do k = 2, km - 1
dm (k) = 0.5 * min (abs (dq (k) + dq (k + 1)), 0.5 * q (k))
if (dq (k) * dq (k + 1) <= 0.) then
if (dq (k) > 0.) then ! local max
dm (k) = min (dm (k), dq (k), - dq (k + 1))
else
dm (k) = 0.
endif
endif
enddo
dm (km) = 0.
! -----------------------------------------------------------------------
! impose a presumed background horizontal variability that is proportional to the value itself
! -----------------------------------------------------------------------
do k = 1, km
dm (k) = max (dm (k), qvmin, h_var * q (k))
enddo
else
do k = 1, km
dm (k) = max (qvmin, h_var * q (k))
enddo
endif
end subroutine linear_prof
! =======================================================================
! ice cloud microphysics processes
! bulk cloud micro - physics; processes splitting
! with some un - split sub - grouping
! time implicit (when possible) accretion and autoconversion
! author: shian - jiann lin, gfdl
! =======================================================================
subroutine icloud (ktop, kbot, tzk, p1, qvk, qlk, qrk, qik, qsk, qgk, dp1, &
den, denfac, vts, vtg, vtr, qak, rh_adj, rh_rain, dts, h_var)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in), dimension (ktop:kbot) :: p1, dp1, den, denfac, vts, vtg, vtr
real, intent (inout), dimension (ktop:kbot) :: tzk, qvk, qlk, qrk, qik, qsk, qgk, qak
real, intent (in) :: rh_adj, rh_rain, dts, h_var
real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, di, lhl, lhi
real, dimension (ktop:kbot) :: cvm, q_liq, q_sol
real :: rdts, fac_g2v, fac_v2g, fac_i2s, fac_imlt
real :: tz, qv, ql, qr, qi, qs, qg, melt
real :: pracs, psacw, pgacw, psacr, pgacr, pgaci, praci, psaci
real :: pgmlt, psmlt, pgfr, pgaut, psaut, pgsub
real :: tc, tsq, dqs0, qden, qim, qsm
real :: dt5, factor, sink, qi_crt
real :: tmp, qsw, qsi, dqsdt, dq
real :: dtmp, qc, q_plus, q_minus
integer :: k
dt5 = 0.5 * dts
rdts = 1. / dts
! -----------------------------------------------------------------------
! define conversion scalar / factor
! -----------------------------------------------------------------------
fac_i2s = 1. - exp (- dts / tau_i2s)
fac_g2v = 1. - exp (- dts / tau_g2v)
fac_v2g = 1. - exp (- dts / tau_v2g)
fac_imlt = 1. - exp (- dt5 / tau_imlt)
! -----------------------------------------------------------------------
! define heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
do k = ktop, kbot
lhi (k) = li00 + dc_ice * tzk (k)
q_liq (k) = qlk (k) + qrk (k)
q_sol (k) = qik (k) + qsk (k) + qgk (k)
cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
icpk (k) = lhi (k) / cvm (k)
enddo
! -----------------------------------------------------------------------
! sources of cloud ice: pihom, cold rain, and the sat_adj
! (initiation plus deposition)
! sources of snow: cold rain, auto conversion + accretion (from cloud ice)
! sat_adj (deposition; requires pre - existing snow) ; initial snow comes from auto conversion
! -----------------------------------------------------------------------
do k = ktop, kbot
if (tzk (k) > tice .and. qik (k) > qcmin) then
! -----------------------------------------------------------------------
! pimlt: instant melting of cloud ice
! -----------------------------------------------------------------------
melt = min (qik (k), fac_imlt * (tzk (k) - tice) / icpk (k))
tmp = min (melt, dim (ql_mlt, qlk (k))) ! max ql amount
qlk (k) = qlk (k) + tmp
qrk (k) = qrk (k) + melt - tmp
qik (k) = qik (k) - melt
q_liq (k) = q_liq (k) + melt
q_sol (k) = q_sol (k) - melt
cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tzk (k) = tzk (k) - melt * lhi (k) / cvm (k)
elseif (tzk (k) < t_wfr .and. qlk (k) > qcmin) then
! -----------------------------------------------------------------------
! pihom: homogeneous freezing of cloud water into cloud ice
! this is the 1st occurance of liquid water freezing in the split mp process
! -----------------------------------------------------------------------
dtmp = t_wfr - tzk (k)
factor = min (1., dtmp / dt_fr)
sink = min (qlk (k) * factor, dtmp / icpk (k))
qi_crt = qi_gen * min (qi_lim, 0.1 * (tice - tzk (k))) / den (k)
tmp = min (sink, dim (qi_crt, qik (k)))
qlk (k) = qlk (k) - sink
qsk (k) = qsk (k) + sink - tmp
qik (k) = qik (k) + tmp
q_liq (k) = q_liq (k) - sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qvk (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tzk (k) = tzk (k) + sink * lhi (k) / cvm (k)
endif
enddo
! -----------------------------------------------------------------------
! vertical subgrid variability
! -----------------------------------------------------------------------
call linear_prof (kbot - ktop + 1, qik (ktop), di (ktop), z_slope_ice, h_var)
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
do k = ktop, kbot
lhl (k) = lv00 + d0_vap * tzk (k)
lhi (k) = li00 + dc_ice * tzk (k)
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
enddo
do k = ktop, kbot
! -----------------------------------------------------------------------
! do nothing above p_min
! -----------------------------------------------------------------------
if (p1 (k) < p_min) cycle
tz = tzk (k)
qv = qvk (k)
ql = qlk (k)
qi = qik (k)
qr = qrk (k)
qs = qsk (k)
qg = qgk (k)
pgacr = 0.
pgacw = 0.
tc = tz - tice
if (tc .ge. 0.) then
! -----------------------------------------------------------------------
! melting of snow
! -----------------------------------------------------------------------
dqs0 = ces0 / p1 (k) - qv
if (qs > qcmin) then
! -----------------------------------------------------------------------
! psacw: accretion of cloud water by snow
! only rate is used (for snow melt) since tc > 0.
! -----------------------------------------------------------------------
if (ql > qrmin) then
factor = denfac (k) * csacw * exp (0.8125 * log (qs * den (k)))
psacw = factor / (1. + dts * factor) * ql ! rate
else
psacw = 0.
endif
! -----------------------------------------------------------------------
! psacr: accretion of rain by melted snow
! pracs: accretion of snow by rain
! -----------------------------------------------------------------------
if (qr > qrmin) then
psacr = min (acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), &
den (k)), qr * rdts)
pracs = acr3d (vtr (k), vts (k), qs, qr, cracs, acco (1, 1), den (k))
else
psacr = 0.
pracs = 0.
endif
! -----------------------------------------------------------------------
! total snow sink:
! psmlt: snow melt (due to rain accretion)
! -----------------------------------------------------------------------
psmlt = max (0., smlt (tc, dqs0, qs * den (k), psacw, psacr, csmlt, &
den (k), denfac (k)))
sink = min (qs, dts * (psmlt + pracs), tc / icpk (k))
qs = qs - sink
! sjl, 20170321:
tmp = min (sink, dim (qs_mlt, ql)) ! max ql due to snow melt
ql = ql + tmp
qr = qr + sink - tmp
! qr = qr + sink
! sjl, 20170321:
q_liq (k) = q_liq (k) + sink
q_sol (k) = q_sol (k) - sink
cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz = tz - sink * lhi (k) / cvm (k)
tc = tz - tice
endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhi (k) = li00 + dc_ice * tz
icpk (k) = lhi (k) / cvm (k)
! -----------------------------------------------------------------------
! melting of graupel
! -----------------------------------------------------------------------
if (qg > qcmin .and. tc > 0.) then
! -----------------------------------------------------------------------
! pgacr: accretion of rain by graupel
! -----------------------------------------------------------------------
if (qr > qrmin) &
pgacr = min (acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), &
den (k)), rdts * qr)
! -----------------------------------------------------------------------
! pgacw: accretion of cloud water by graupel
! -----------------------------------------------------------------------
qden = qg * den (k)
if (ql > qrmin) then
factor = cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden)))
pgacw = factor / (1. + dts * factor) * ql ! rate
endif
! -----------------------------------------------------------------------
! pgmlt: graupel melt
! -----------------------------------------------------------------------
pgmlt = dts * gmlt (tc, dqs0, qden, pgacw, pgacr, cgmlt, den (k))
pgmlt = min (max (0., pgmlt), qg, tc / icpk (k))
qg = qg - pgmlt
qr = qr + pgmlt
q_liq (k) = q_liq (k) + pgmlt
q_sol (k) = q_sol (k) - pgmlt
cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz = tz - pgmlt * lhi (k) / cvm (k)
endif
else
! -----------------------------------------------------------------------
! cloud ice proc:
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! psaci: accretion of cloud ice by snow
! -----------------------------------------------------------------------
if (qi > 3.e-7) then ! cloud ice sink terms
if (qs > 1.e-7) then
! -----------------------------------------------------------------------
! sjl added (following lin eq. 23) the temperature dependency
! to reduce accretion, use esi = exp (0.05 * tc) as in hong et al 2004
! -----------------------------------------------------------------------
factor = dts * denfac (k) * csaci * exp (0.05 * tc + 0.8125 * log (qs * den (k)))
psaci = factor / (1. + factor) * qi
else
psaci = 0.
endif
! -----------------------------------------------------------------------
! pasut: autoconversion: cloud ice -- > snow
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! similar to lfo 1983: eq. 21 solved implicitly
! threshold from wsm6 scheme, hong et al 2004, eq (13) : qi0_crt ~0.8e-4
! -----------------------------------------------------------------------
qim = qi0_crt / den (k)
! -----------------------------------------------------------------------
! assuming linear subgrid vertical distribution of cloud ice
! the mismatch computation following lin et al. 1994, mwr
! -----------------------------------------------------------------------
if (const_vi) then
tmp = fac_i2s
else
tmp = fac_i2s * exp (0.025 * tc)
endif
di (k) = max (di (k), qrmin)
q_plus = qi + di (k)
if (q_plus > (qim + qrmin)) then
if (qim > (qi - di (k))) then
dq = (0.25 * (q_plus - qim) ** 2) / di (k)
else
dq = qi - qim
endif
psaut = tmp * dq
else
psaut = 0.
endif
! -----------------------------------------------------------------------
! sink is no greater than 75% of qi
! -----------------------------------------------------------------------
sink = min (0.75 * qi, psaci + psaut)
qi = qi - sink
qs = qs + sink
! -----------------------------------------------------------------------
! pgaci: accretion of cloud ice by graupel
! -----------------------------------------------------------------------
if (qg > 1.e-6) then
! -----------------------------------------------------------------------
! factor = dts * cgaci / sqrt (den (k)) * exp (0.05 * tc + 0.875 * log (qg * den (k)))
! simplified form: remove temp dependency & set the exponent "0.875" -- > 1
! -----------------------------------------------------------------------
factor = dts * cgaci * sqrt (den (k)) * qg
pgaci = factor / (1. + factor) * qi
qi = qi - pgaci
qg = qg + pgaci
endif
endif
! -----------------------------------------------------------------------
! cold - rain proc:
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! rain to ice, snow, graupel processes:
! -----------------------------------------------------------------------
tc = tz - tice
if (qr > 1.e-7 .and. tc < 0.) then
! -----------------------------------------------------------------------
! * sink * terms to qr: psacr + pgfr
! source terms to qs: psacr
! source terms to qg: pgfr
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! psacr accretion of rain by snow
! -----------------------------------------------------------------------
if (qs > 1.e-7) then ! if snow exists
psacr = dts * acr3d (vts (k), vtr (k), qr, qs, csacr, acco (1, 2), den (k))
else
psacr = 0.
endif
! -----------------------------------------------------------------------
! pgfr: rain freezing -- > graupel
! -----------------------------------------------------------------------
pgfr = dts * cgfr (1) / den (k) * (exp (- cgfr (2) * tc) - 1.) * &
exp (1.75 * log (qr * den (k)))
! -----------------------------------------------------------------------
! total sink to qr
! -----------------------------------------------------------------------
sink = psacr + pgfr
factor = min (sink, qr, - tc / icpk (k)) / max (sink, qrmin)
psacr = factor * psacr
pgfr = factor * pgfr
sink = psacr + pgfr
qr = qr - sink
qs = qs + psacr
qg = qg + pgfr
q_liq (k) = q_liq (k) - sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz = tz + sink * lhi (k) / cvm (k)
endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhi (k) = li00 + dc_ice * tz
icpk (k) = lhi (k) / cvm (k)
! -----------------------------------------------------------------------
! graupel production terms:
! -----------------------------------------------------------------------
if (qs > 1.e-7) then
! -----------------------------------------------------------------------
! accretion: snow -- > graupel
! -----------------------------------------------------------------------
if (qg > qrmin) then
sink = dts * acr3d (vtg (k), vts (k), qs, qg, cgacs, acco (1, 4), den (k))
else
sink = 0.
endif
! -----------------------------------------------------------------------
! autoconversion snow -- > graupel
! -----------------------------------------------------------------------
qsm = qs0_crt / den (k)
if (qs > qsm) then
factor = dts * 1.e-3 * exp (0.09 * (tz - tice))
sink = sink + factor / (1. + factor) * (qs - qsm)
endif
sink = min (qs, sink)
qs = qs - sink
qg = qg + sink
endif ! snow existed
if (qg > 1.e-7 .and. tz < tice0) then
! -----------------------------------------------------------------------
! pgacw: accretion of cloud water by graupel
! -----------------------------------------------------------------------
if (ql > 1.e-6) then
qden = qg * den (k)
factor = dts * cgacw * qden / sqrt (den (k) * sqrt (sqrt (qden)))
pgacw = factor / (1. + factor) * ql
else
pgacw = 0.
endif
! -----------------------------------------------------------------------
! pgacr: accretion of rain by graupel
! -----------------------------------------------------------------------
if (qr > 1.e-6) then
pgacr = min (dts * acr3d (vtg (k), vtr (k), qr, qg, cgacr, acco (1, 3), &
den (k)), qr)
else
pgacr = 0.
endif
sink = pgacr + pgacw
factor = min (sink, dim (tice, tz) / icpk (k)) / max (sink, qrmin)
pgacr = factor * pgacr
pgacw = factor * pgacw
sink = pgacr + pgacw
qg = qg + sink
qr = qr - pgacr
ql = ql - pgacw
q_liq (k) = q_liq (k) - sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz = tz + sink * lhi (k) / cvm (k)
endif
endif
tzk (k) = tz
qvk (k) = qv
qlk (k) = ql
qik (k) = qi
qrk (k) = qr
qsk (k) = qs
qgk (k) = qg
enddo
! -----------------------------------------------------------------------
! subgrid cloud microphysics
! -----------------------------------------------------------------------
call subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tzk, qvk, &
qlk, qrk, qik, qsk, qgk, qak, h_var, rh_rain)
end subroutine icloud
! =======================================================================
! temperature sentive high vertical resolution processes
! =======================================================================
subroutine subgrid_z_proc (ktop, kbot, p1, den, denfac, dts, rh_adj, tz, qv, &
ql, qr, qi, qs, qg, qa, h_var, rh_rain)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in), dimension (ktop:kbot) :: p1, den, denfac
real, intent (in) :: dts, rh_adj, h_var, rh_rain
real, intent (inout), dimension (ktop:kbot) :: tz, qv, ql, qr, qi, qs, qg, qa
real, dimension (ktop:kbot) :: lcpk, icpk, tcpk, tcp3, lhl, lhi
real, dimension (ktop:kbot) :: cvm, q_liq, q_sol, q_cond
real :: fac_v2l, fac_l2v
real :: pidep, qi_crt
! -----------------------------------------------------------------------
! qstar over water may be accurate only down to - 80 deg c with ~10% uncertainty
! must not be too large to allow psc
! -----------------------------------------------------------------------
real :: rh, rqi, tin, qsw, qsi, qpz, qstar
real :: dqsdt, dwsdt, dq, dq0, factor, tmp
real :: q_plus, q_minus, dt_evap, dt_pisub
real :: evap, sink, tc, pisub, q_adj, dtmp
real :: pssub, pgsub, tsq, qden, fac_g2v, fac_v2g
integer :: k
if (fast_sat_adj) then
dt_evap = 0.5 * dts
else
dt_evap = dts
endif
! -----------------------------------------------------------------------
! define conversion scalar / factor
! -----------------------------------------------------------------------
fac_v2l = 1. - exp (- dt_evap / tau_v2l)
fac_l2v = 1. - exp (- dt_evap / tau_l2v)
fac_g2v = 1. - exp (- dts / tau_g2v)
fac_v2g = 1. - exp (- dts / tau_v2g)
! -----------------------------------------------------------------------
! define heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
do k = ktop, kbot
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
q_liq (k) = ql (k) + qr (k)
q_sol (k) = qi (k) + qs (k) + qg (k)
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr))
enddo
do k = ktop, kbot
if (p1 (k) < p_min) cycle
! -----------------------------------------------------------------------
! instant deposit all water vapor to cloud ice when temperature is super low
! -----------------------------------------------------------------------
if (tz (k) < t_min) then
sink = dim (qv (k), 1.e-7)
qv (k) = qv (k) - sink
qi (k) = qi (k) + sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k)
if (.not. do_qa) qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover
cycle
endif
! -----------------------------------------------------------------------
! update heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
tcp3 (k) = lcpk (k) + icpk (k) * min (1., dim (tice, tz (k)) / (tice - t_wfr))
! -----------------------------------------------------------------------
! instant evaporation / sublimation of all clouds if rh < rh_adj -- > cloud free
! -----------------------------------------------------------------------
qpz = qv (k) + ql (k) + qi (k)
tin = tz (k) - (lhl (k) * (ql (k) + qi (k)) + lhi (k) * qi (k)) / (c_air + &
qpz * c_vap + qr (k) * c_liq + (qs (k) + qg (k)) * c_ice)
if (tin > t_sub + 6.) then
rh = qpz / iqs1 (tin, den (k))
if (rh < rh_adj) then ! qpz / rh_adj < qs
tz (k) = tin
qv (k) = qpz
ql (k) = 0.
qi (k) = 0.
cycle ! cloud free
endif
endif
! -----------------------------------------------------------------------
! cloud water < -- > vapor adjustment:
! -----------------------------------------------------------------------
qsw = wqs2 (tz (k), den (k), dwsdt)
dq0 = qsw - qv (k)
if (dq0 > 0.) then
! SJL 20170703 added ql factor to prevent the situation of high ql and low RH
! factor = min (1., fac_l2v * sqrt (max (0., ql (k)) / 1.e-5) * 10. * dq0 / qsw)
! factor = fac_l2v
! factor = 1
factor = min (1., fac_l2v * (10. * dq0 / qsw)) ! the rh dependent factor = 1 at 90%
evap = min (ql (k), factor * dq0 / (1. + tcp3 (k) * dwsdt))
else ! condensate all excess vapor into cloud water
! -----------------------------------------------------------------------
! evap = fac_v2l * dq0 / (1. + tcp3 (k) * dwsdt)
! sjl, 20161108
! -----------------------------------------------------------------------
evap = dq0 / (1. + tcp3 (k) * dwsdt)
endif
qv (k) = qv (k) + evap
ql (k) = ql (k) - evap
q_liq (k) = q_liq (k) - evap
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) - evap * lhl (k) / cvm (k)
! -----------------------------------------------------------------------
! update heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
lhi (k) = li00 + dc_ice * tz (k)
icpk (k) = lhi (k) / cvm (k)
! -----------------------------------------------------------------------
! enforce complete freezing below - 48 c
! -----------------------------------------------------------------------
dtmp = t_wfr - tz (k) ! [ - 40, - 48]
if (dtmp > 0. .and. ql (k) > qcmin) then
sink = min (ql (k), ql (k) * dtmp * 0.125, dtmp / icpk (k))
ql (k) = ql (k) - sink
qi (k) = qi (k) + sink
q_liq (k) = q_liq (k) - sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) + sink * lhi (k) / cvm (k)
endif
! -----------------------------------------------------------------------
! update heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
lhi (k) = li00 + dc_ice * tz (k)
icpk (k) = lhi (k) / cvm (k)
! -----------------------------------------------------------------------
! bigg mechanism
! -----------------------------------------------------------------------
if (fast_sat_adj) then
dt_pisub = 0.5 * dts
else
dt_pisub = dts
tc = tice - tz (k)
if (ql (k) > qrmin .and. tc > 0.) then
sink = 3.3333e-10 * dts * (exp (0.66 * tc) - 1.) * den (k) * ql (k) * ql (k)
sink = min (ql (k), tc / icpk (k), sink)
ql (k) = ql (k) - sink
qi (k) = qi (k) + sink
q_liq (k) = q_liq (k) - sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) + sink * lhi (k) / cvm (k)
endif ! significant ql existed
endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
! -----------------------------------------------------------------------
! sublimation / deposition of ice
! -----------------------------------------------------------------------
if (tz (k) < tice) then
qsi = iqs2 (tz (k), den (k), dqsdt)
dq = qv (k) - qsi
sink = dq / (1. + tcpk (k) * dqsdt)
if (qi (k) > qrmin) then
! eq 9, hong et al. 2004, mwr
! for a and b, see dudhia 1989: page 3103 eq (b7) and (b8)
pidep = dt_pisub * dq * 349138.78 * exp (0.875 * log (qi (k) * den (k))) &
/ (qsi * den (k) * lat2 / (0.0243 * rvgas * tz (k) ** 2) + 4.42478e4)
else
pidep = 0.
endif
if (dq > 0.) then ! vapor - > ice
tmp = tice - tz (k)
! 20160912: the following should produce more ice at higher altitude
! qi_crt = 4.92e-11 * exp (1.33 * log (1.e3 * exp (0.1 * tmp))) / den (k)
qi_crt = qi_gen * min (qi_lim, 0.1 * tmp) / den (k)
sink = min (sink, max (qi_crt - qi (k), pidep), tmp / tcpk (k))
else ! ice -- > vapor
pidep = pidep * min (1., dim (tz (k), t_sub) * 0.2)
sink = max (pidep, sink, - qi (k))
endif
qv (k) = qv (k) - sink
qi (k) = qi (k) + sink
q_sol (k) = q_sol (k) + sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) + sink * (lhl (k) + lhi (k)) / cvm (k)
endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
! -----------------------------------------------------------------------
! sublimation / deposition of snow
! this process happens for all temp rage
! -----------------------------------------------------------------------
if (qs (k) > qrmin) then
qsi = iqs2 (tz (k), den (k), dqsdt)
qden = qs (k) * den (k)
tmp = exp (0.65625 * log (qden))
tsq = tz (k) * tz (k)
dq = (qsi - qv (k)) / (1. + tcpk (k) * dqsdt)
pssub = cssub (1) * tsq * (cssub (2) * sqrt (qden) + cssub (3) * tmp * &
sqrt (denfac (k))) / (cssub (4) * tsq + cssub (5) * qsi * den (k))
pssub = (qsi - qv (k)) * dts * pssub
if (pssub > 0.) then ! qs -- > qv, sublimation
pssub = min (pssub * min (1., dim (tz (k), t_sub) * 0.2), qs (k))
else
if (tz (k) > tice) then
pssub = 0. ! no deposition
else
pssub = max (pssub, dq, (tz (k) - tice) / tcpk (k))
endif
endif
qs (k) = qs (k) - pssub
qv (k) = qv (k) + pssub
q_sol (k) = q_sol (k) - pssub
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) - pssub * (lhl (k) + lhi (k)) / cvm (k)
endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
tcpk (k) = lcpk (k) + icpk (k)
! -----------------------------------------------------------------------
! simplified 2 - way grapuel sublimation - deposition mechanism
! -----------------------------------------------------------------------
if (qg (k) > qrmin) then
qsi = iqs2 (tz (k), den (k), dqsdt)
dq = (qv (k) - qsi) / (1. + tcpk (k) * dqsdt)
pgsub = (qv (k) / qsi - 1.) * qg (k)
if (pgsub > 0.) then ! deposition
if (tz (k) > tice) then
pgsub = 0. ! no deposition
else
pgsub = min (fac_v2g * pgsub, 0.2 * dq, ql (k) + qr (k), &
(tice - tz (k)) / tcpk (k))
endif
else ! submilation
pgsub = max (fac_g2v * pgsub, dq) * min (1., dim (tz (k), t_sub) * 0.1)
endif
qg (k) = qg (k) + pgsub
qv (k) = qv (k) - pgsub
q_sol (k) = q_sol (k) + pgsub
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) + pgsub * (lhl (k) + lhi (k)) / cvm (k)
endif
#ifdef USE_MIN_EVAP
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
lcpk (k) = lhl (k) / cvm (k)
! -----------------------------------------------------------------------
! * minimum evap of rain in dry environmental air
! -----------------------------------------------------------------------
if (qr (k) > qcmin) then
qsw = wqs2 (tz (k), den (k), dqsdt)
sink = min (qr (k), dim (rh_rain * qsw, qv (k)) / (1. + lcpk (k) * dqsdt))
qv (k) = qv (k) + sink
qr (k) = qr (k) - sink
q_liq (k) = q_liq (k) - sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) - sink * lhl (k) / cvm (k)
endif
#endif
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
lhl (k) = lv00 + d0_vap * tz (k)
cvm (k) = c_air + (qv (k) + q_liq (k) + q_sol (k)) * c_vap
lcpk (k) = lhl (k) / cvm (k)
! -----------------------------------------------------------------------
! compute cloud fraction
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! combine water species
! -----------------------------------------------------------------------
if (do_qa) cycle
if (rad_snow) then
q_sol (k) = qi (k) + qs (k)
else
q_sol (k) = qi (k)
endif
if (rad_rain) then
q_liq (k) = ql (k) + qr (k)
else
q_liq (k) = ql (k)
endif
q_cond (k) = q_liq (k) + q_sol (k)
qpz = qv (k) + q_cond (k) ! qpz is conserved
! -----------------------------------------------------------------------
! use the "liquid - frozen water temperature" (tin) to compute saturated specific humidity
! -----------------------------------------------------------------------
tin = tz (k) - (lcpk (k) * q_cond (k) + icpk (k) * q_sol (k)) ! minimum temperature
! tin = tz (k) - ((lv00 + d0_vap * tz (k)) * q_cond (k) + &
! (li00 + dc_ice * tz (k)) * q_sol (k)) / (c_air + qpz * c_vap)
! -----------------------------------------------------------------------
! determine saturated specific humidity
! -----------------------------------------------------------------------
if (tin <= t_wfr) then
! ice phase:
qstar = iqs1 (tin, den (k))
elseif (tin >= tice) then
! liquid phase:
qstar = wqs1 (tin, den (k))
else
! mixed phase:
qsi = iqs1 (tin, den (k))
qsw = wqs1 (tin, den (k))
if (q_cond (k) > 3.e-6) then
rqi = q_sol (k) / q_cond (k)
else
! -----------------------------------------------------------------------
! mostly liquid water q_cond (k) at initial cloud development stage
! -----------------------------------------------------------------------
rqi = (tice - tin) / (tice - t_wfr)
endif
qstar = rqi * qsi + (1. - rqi) * qsw
endif
! -----------------------------------------------------------------------
! assuming subgrid linear distribution in horizontal; this is effectively a smoother for the
! binary cloud scheme
! -----------------------------------------------------------------------
if (qpz > qrmin) then
! partial cloudiness by pdf:
dq = max (qcmin, h_var * qpz)
q_plus = qpz + dq ! cloud free if qstar > q_plus
q_minus = qpz - dq
if (qstar < q_minus) then
qa (k) = qa (k) + 1. ! air fully saturated; 100 % cloud cover
elseif (qstar < q_plus .and. q_cond (k) > qc_crt) then
qa (k) = qa (k) + (q_plus - qstar) / (dq + dq) ! partial cloud cover
! qa (k) = sqrt (qa (k) + (q_plus - qstar) / (dq + dq))
endif
endif
enddo
end subroutine subgrid_z_proc
! =======================================================================
! rain evaporation
! =======================================================================
subroutine revap_rac1 (hydrostatic, is, ie, dt, tz, qv, ql, qr, qi, qs, qg, den, hvar)
implicit none
logical, intent (in) :: hydrostatic
integer, intent (in) :: is, ie
real, intent (in) :: dt ! time step (s)
real, intent (in), dimension (is:ie) :: den, hvar, qi, qs, qg
real, intent (inout), dimension (is:ie) :: tz, qv, qr, ql
real, dimension (is:ie) :: lcp2, denfac, q_liq, q_sol, cvm, lhl
real :: dqv, qsat, dqsdt, evap, qden, q_plus, q_minus, sink
real :: tin, t2, qpz, dq, dqh
integer :: i
! -----------------------------------------------------------------------
! define latend heat coefficient
! -----------------------------------------------------------------------
do i = is, ie
lhl (i) = lv00 + d0_vap * tz (i)
q_liq (i) = ql (i) + qr (i)
q_sol (i) = qi (i) + qs (i) + qg (i)
cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
lcp2 (i) = lhl (i) / cvm (i)
! denfac (i) = sqrt (sfcrho / den (i))
enddo
do i = is, ie
if (qr (i) > qrmin .and. tz (i) > t_wfr) then
qpz = qv (i) + ql (i)
tin = tz (i) - lcp2 (i) * ql (i) ! presence of clouds suppresses the rain evap
qsat = wqs2 (tin, den (i), dqsdt)
dqh = max (ql (i), hvar (i) * max (qpz, qcmin))
dqv = qsat - qv (i)
q_minus = qpz - dqh
q_plus = qpz + dqh
! -----------------------------------------------------------------------
! qsat must be > q_minus to activate evaporation
! qsat must be < q_plus to activate accretion
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! rain evaporation
! -----------------------------------------------------------------------
if (dqv > qvmin .and. qsat > q_minus) then
if (qsat > q_plus) then
dq = qsat - qpz
else
! q_minus < qsat < q_plus
! dq == dqh if qsat == q_minus
dq = 0.25 * (q_minus - qsat) ** 2 / dqh
endif
qden = qr (i) * den (i)
t2 = tin * tin
evap = crevp (1) * t2 * dq * (crevp (2) * sqrt (qden) + crevp (3) * exp (0.725 * log (qden))) &
/ (crevp (4) * t2 + crevp (5) * qsat * den (i))
evap = min (qr (i), dt * evap, dqv / (1. + lcp2 (i) * dqsdt))
qr (i) = qr (i) - evap
qv (i) = qv (i) + evap
q_liq (i) = q_liq (i) - evap
cvm (i) = c_air + qv (i) * c_vap + q_liq (i) * c_liq + q_sol (i) * c_ice
tz (i) = tz (i) - evap * lhl (i) / cvm (i)
endif
! -----------------------------------------------------------------------
! accretion: pracc
! -----------------------------------------------------------------------
if (qr (i) > qrmin .and. ql (i) > 1.e-8 .and. qsat < q_plus) then
denfac (i) = sqrt (sfcrho / den (i))
sink = dt * denfac (i) * cracw * exp (0.95 * log (qr (i) * den (i)))
sink = sink / (1. + sink) * ql (i)
ql (i) = ql (i) - sink
qr (i) = qr (i) + sink
endif
endif
enddo
end subroutine revap_rac1
! =======================================================================
! compute terminal fall speed
! consider cloud ice, snow, and graupel's melting during fall
! =======================================================================
subroutine terminal_fall (dtm, ktop, kbot, tz, qv, ql, qr, qg, qs, qi, dz, dp, &
den, vtg, vts, vti, r1, g1, s1, i1, m1_sol, w1)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: dtm ! time step (s)
real, intent (in), dimension (ktop:kbot) :: vtg, vts, vti, den, dp, dz
real, intent (inout), dimension (ktop:kbot) :: qv, ql, qr, qg, qs, qi, tz, m1_sol, w1
real, intent (out) :: r1, g1, s1, i1
real, dimension (ktop:kbot + 1) :: ze, zt
real :: qsat, dqsdt, dt5, evap, dtime
real :: factor, frac
real :: tmp, precip, tc, sink
real, dimension (ktop:kbot) :: lcpk, icpk, cvm, q_liq, q_sol, lhl, lhi
real, dimension (ktop:kbot) :: m1, dm
real :: zs = 0.
real :: fac_imlt
integer :: k, k0, m
logical :: no_fall
dt5 = 0.5 * dtm
fac_imlt = 1. - exp (- dt5 / tau_imlt)
! -----------------------------------------------------------------------
! define heat capacity and latend heat coefficient
! -----------------------------------------------------------------------
do k = ktop, kbot
m1_sol (k) = 0.
lhl (k) = lv00 + d0_vap * tz (k)
lhi (k) = li00 + dc_ice * tz (k)
q_liq (k) = ql (k) + qr (k)
q_sol (k) = qi (k) + qs (k) + qg (k)
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
lcpk (k) = lhl (k) / cvm (k)
icpk (k) = lhi (k) / cvm (k)
enddo
! -----------------------------------------------------------------------
! find significant melting level
! -----------------------------------------------------------------------
k0 = kbot
do k = ktop, kbot - 1
if (tz (k) > tice) then
k0 = k
exit
endif
enddo
! -----------------------------------------------------------------------
! melting of cloud_ice (before fall) :
! -----------------------------------------------------------------------
do k = k0, kbot
tc = tz (k) - tice
if (qi (k) > qcmin .and. tc > 0.) then
sink = min (qi (k), fac_imlt * tc / icpk (k))
tmp = min (sink, dim (ql_mlt, ql (k)))
ql (k) = ql (k) + tmp
qr (k) = qr (k) + sink - tmp
qi (k) = qi (k) - sink
q_liq (k) = q_liq (k) + sink
q_sol (k) = q_sol (k) - sink
cvm (k) = c_air + qv (k) * c_vap + q_liq (k) * c_liq + q_sol (k) * c_ice
tz (k) = tz (k) - sink * lhi (k) / cvm (k)
tc = tz (k) - tice
endif
enddo
! -----------------------------------------------------------------------
! turn off melting when cloud microphysics time step is small
! -----------------------------------------------------------------------
if (dtm < 60.) k0 = kbot
! sjl, turn off melting of falling cloud ice, snow and graupel
k0 = kbot
! sjl, turn off melting of falling cloud ice, snow and graupel
ze (kbot + 1) = zs
do k = kbot, ktop, - 1
ze (k) = ze (k + 1) - dz (k) ! dz < 0
enddo
zt (ktop) = ze (ktop)
! -----------------------------------------------------------------------
! update capacity heat and latend heat coefficient
! -----------------------------------------------------------------------
do k = k0, kbot
lhi (k) = li00 + dc_ice * tz (k)
icpk (k) = lhi (k) / cvm (k)
enddo
! -----------------------------------------------------------------------
! melting of falling cloud ice into rain
! -----------------------------------------------------------------------
call check_column (ktop, kbot, qi, no_fall)
if (vi_fac < 1.e-5 .or. no_fall) then
i1 = 0.
else
do k = ktop + 1, kbot
zt (k) = ze (k) - dt5 * (vti (k - 1) + vti (k))
enddo
zt (kbot + 1) = zs - dtm * vti (kbot)
do k = ktop, kbot
if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
enddo
if (k0 < kbot) then
do k = kbot - 1, k0, - 1
if (qi (k) > qrmin) then
do m = k + 1, kbot
if (zt (k + 1) >= ze (m)) exit
if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
dtime = min (1.0, (ze (m) - ze (m + 1)) / (max (vr_min, vti (k)) * tau_imlt))
sink = min (qi (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
tmp = min (sink, dim (ql_mlt, ql (m)))
ql (m) = ql (m) + tmp
qr (m) = qr (m) - tmp + sink
tz (m) = tz (m) - sink * icpk (m)
qi (k) = qi (k) - sink * dp (m) / dp (k)
endif
enddo
endif
enddo
endif
if (do_sedi_w) then
do k = ktop, kbot
dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
enddo
endif
if (use_ppm) then
call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qi, i1, m1_sol, mono_prof)
else
call implicit_fall (dtm, ktop, kbot, ze, vti, dp, qi, i1, m1_sol)
endif
if (do_sedi_w) then
w1 (ktop) = (dm (ktop) * w1 (ktop) + m1_sol (ktop) * vti (ktop)) / (dm (ktop) - m1_sol (ktop))
do k = ktop + 1, kbot
w1 (k) = (dm (k) * w1 (k) - m1_sol (k - 1) * vti (k - 1) + m1_sol (k) * vti (k)) &
/ (dm (k) + m1_sol (k - 1) - m1_sol (k))
enddo
endif
endif
! -----------------------------------------------------------------------
! melting of falling snow into rain
! -----------------------------------------------------------------------
r1 = 0.
call check_column (ktop, kbot, qs, no_fall)
if (no_fall) then
s1 = 0.
else
do k = ktop + 1, kbot
zt (k) = ze (k) - dt5 * (vts (k - 1) + vts (k))
enddo
zt (kbot + 1) = zs - dtm * vts (kbot)
do k = ktop, kbot
if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
enddo
if (k0 < kbot) then
do k = kbot - 1, k0, - 1
if (qs (k) > qrmin) then
do m = k + 1, kbot
if (zt (k + 1) >= ze (m)) exit
dtime = min (dtm, (ze (m) - ze (m + 1)) / (vr_min + vts (k)))
if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
dtime = min (1.0, dtime / tau_smlt)
sink = min (qs (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
tz (m) = tz (m) - sink * icpk (m)
qs (k) = qs (k) - sink * dp (m) / dp (k)
if (zt (k) < zs) then
r1 = r1 + sink * dp (m) ! precip as rain
else
! qr source here will fall next time step (therefore, can evap)
qr (m) = qr (m) + sink
endif
endif
if (qs (k) < qrmin) exit
enddo
endif
enddo
endif
if (do_sedi_w) then
do k = ktop, kbot
dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
enddo
endif
if (use_ppm) then
call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qs, s1, m1, mono_prof)
else
call implicit_fall (dtm, ktop, kbot, ze, vts, dp, qs, s1, m1)
endif
do k = ktop, kbot
m1_sol (k) = m1_sol (k) + m1 (k)
enddo
if (do_sedi_w) then
w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vts (ktop)) / (dm (ktop) - m1 (ktop))
do k = ktop + 1, kbot
w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vts (k - 1) + m1 (k) * vts (k)) &
/ (dm (k) + m1 (k - 1) - m1 (k))
enddo
endif
endif
! ----------------------------------------------
! melting of falling graupel into rain
! ----------------------------------------------
call check_column (ktop, kbot, qg, no_fall)
if (no_fall) then
g1 = 0.
else
do k = ktop + 1, kbot
zt (k) = ze (k) - dt5 * (vtg (k - 1) + vtg (k))
enddo
zt (kbot + 1) = zs - dtm * vtg (kbot)
do k = ktop, kbot
if (zt (k + 1) >= zt (k)) zt (k + 1) = zt (k) - dz_min
enddo
if (k0 < kbot) then
do k = kbot - 1, k0, - 1
if (qg (k) > qrmin) then
do m = k + 1, kbot
if (zt (k + 1) >= ze (m)) exit
dtime = min (dtm, (ze (m) - ze (m + 1)) / vtg (k))
if (zt (k) < ze (m + 1) .and. tz (m) > tice) then
dtime = min (1., dtime / tau_g2r)
sink = min (qg (k) * dp (k) / dp (m), dtime * (tz (m) - tice) / icpk (m))
tz (m) = tz (m) - sink * icpk (m)
qg (k) = qg (k) - sink * dp (m) / dp (k)
if (zt (k) < zs) then
r1 = r1 + sink * dp (m)
else
qr (m) = qr (m) + sink
endif
endif
if (qg (k) < qrmin) exit
enddo
endif
enddo
endif
if (do_sedi_w) then
do k = ktop, kbot
dm (k) = dp (k) * (1. + qv (k) + ql (k) + qr (k) + qi (k) + qs (k) + qg (k))
enddo
endif
if (use_ppm) then
call lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, qg, g1, m1, mono_prof)
else
call implicit_fall (dtm, ktop, kbot, ze, vtg, dp, qg, g1, m1)
endif
do k = ktop, kbot
m1_sol (k) = m1_sol (k) + m1 (k)
enddo
if (do_sedi_w) then
w1 (ktop) = (dm (ktop) * w1 (ktop) + m1 (ktop) * vtg (ktop)) / (dm (ktop) - m1 (ktop))
do k = ktop + 1, kbot
w1 (k) = (dm (k) * w1 (k) - m1 (k - 1) * vtg (k - 1) + m1 (k) * vtg (k)) &
/ (dm (k) + m1 (k - 1) - m1 (k))
enddo
endif
endif
end subroutine terminal_fall
! =======================================================================
! check if water species large enough to fall
! =======================================================================
subroutine check_column (ktop, kbot, q, no_fall)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: q (ktop:kbot)
logical, intent (out) :: no_fall
integer :: k
no_fall = .true.
do k = ktop, kbot
if (q (k) > qrmin) then
no_fall = .false.
exit
endif
enddo
end subroutine check_column
! =======================================================================
! time - implicit monotonic scheme
! developed by sj lin, 2016
! =======================================================================
subroutine implicit_fall (dt, ktop, kbot, ze, vt, dp, q, precip, m1)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: dt
real, intent (in), dimension (ktop:kbot + 1) :: ze
real, intent (in), dimension (ktop:kbot) :: vt, dp
real, intent (inout), dimension (ktop:kbot) :: q
real, intent (out), dimension (ktop:kbot) :: m1
real, intent (out) :: precip
real, dimension (ktop:kbot) :: dz, qm, dd
integer :: k
do k = ktop, kbot
dz (k) = ze (k) - ze (k + 1)
dd (k) = dt * vt (k)
q (k) = q (k) * dp (k)
enddo
! -----------------------------------------------------------------------
! sedimentation: non - vectorizable loop
! -----------------------------------------------------------------------
qm (ktop) = q (ktop) / (dz (ktop) + dd (ktop))
do k = ktop + 1, kbot
qm (k) = (q (k) + dd (k - 1) * qm (k - 1)) / (dz (k) + dd (k))
enddo
! -----------------------------------------------------------------------
! qm is density at this stage
! -----------------------------------------------------------------------
do k = ktop, kbot
qm (k) = qm (k) * dz (k)
enddo
! -----------------------------------------------------------------------
! output mass fluxes: non - vectorizable loop
! -----------------------------------------------------------------------
m1 (ktop) = q (ktop) - qm (ktop)
do k = ktop + 1, kbot
m1 (k) = m1 (k - 1) + q (k) - qm (k)
enddo
precip = m1 (kbot)
! -----------------------------------------------------------------------
! update:
! -----------------------------------------------------------------------
do k = ktop, kbot
q (k) = qm (k) / dp (k)
enddo
end subroutine implicit_fall
! =======================================================================
! lagrangian scheme
! developed by sj lin, ????
! =======================================================================
subroutine lagrangian_fall_ppm (ktop, kbot, zs, ze, zt, dp, q, precip, m1, mono)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in) :: zs
logical, intent (in) :: mono
real, intent (in), dimension (ktop:kbot + 1) :: ze, zt
real, intent (in), dimension (ktop:kbot) :: dp
! m1: flux
real, intent (inout), dimension (ktop:kbot) :: q, m1
real, intent (out) :: precip
real, dimension (ktop:kbot) :: qm, dz
real :: a4 (4, ktop:kbot)
real :: pl, pr, delz, esl
integer :: k, k0, n, m
real, parameter :: r3 = 1. / 3., r23 = 2. / 3.
! -----------------------------------------------------------------------
! density:
! -----------------------------------------------------------------------
do k = ktop, kbot
dz (k) = zt (k) - zt (k + 1) ! note: dz is positive
q (k) = q (k) * dp (k)
a4 (1, k) = q (k) / dz (k)
qm (k) = 0.
enddo
! -----------------------------------------------------------------------
! construct vertical profile with zt as coordinate
! -----------------------------------------------------------------------
call cs_profile (a4 (1, ktop), dz (ktop), kbot - ktop + 1, mono)
k0 = ktop
do k = ktop, kbot
do n = k0, kbot
if (ze (k) <= zt (n) .and. ze (k) >= zt (n + 1)) then
pl = (zt (n) - ze (k)) / dz (n)
if (zt (n + 1) <= ze (k + 1)) then
! entire new grid is within the original grid
pr = (zt (n) - ze (k + 1)) / dz (n)
qm (k) = a4 (2, n) + 0.5 * (a4 (4, n) + a4 (3, n) - a4 (2, n)) * (pr + pl) - &
a4 (4, n) * r3 * (pr * (pr + pl) + pl ** 2)
qm (k) = qm (k) * (ze (k) - ze (k + 1))
k0 = n
goto 555
else
qm (k) = (ze (k) - zt (n + 1)) * (a4 (2, n) + 0.5 * (a4 (4, n) + &
a4 (3, n) - a4 (2, n)) * (1. + pl) - a4 (4, n) * (r3 * (1. + pl * (1. + pl))))
if (n < kbot) then
do m = n + 1, kbot
! locate the bottom edge: ze (k + 1)
if (ze (k + 1) < zt (m + 1)) then
qm (k) = qm (k) + q (m)
else
delz = zt (m) - ze (k + 1)
esl = delz / dz (m)
qm (k) = qm (k) + delz * (a4 (2, m) + 0.5 * esl * &
(a4 (3, m) - a4 (2, m) + a4 (4, m) * (1. - r23 * esl)))
k0 = m
goto 555
endif
enddo
endif
goto 555
endif
endif
enddo
555 continue
enddo
m1 (ktop) = q (ktop) - qm (ktop)
do k = ktop + 1, kbot
m1 (k) = m1 (k - 1) + q (k) - qm (k)
enddo
precip = m1 (kbot)
! convert back to * dry * mixing ratio:
! dp must be dry air_mass (because moist air mass will be changed due to terminal fall) .
do k = ktop, kbot
q (k) = qm (k) / dp (k)
enddo
end subroutine lagrangian_fall_ppm
subroutine cs_profile (a4, del, km, do_mono)
implicit none
integer, intent (in) :: km ! vertical dimension
real, intent (in) :: del (km)
logical, intent (in) :: do_mono
real, intent (inout) :: a4 (4, km)
real, parameter :: qp_min = 1.e-6
real :: gam (km)
real :: q (km + 1)
real :: d4, bet, a_bot, grat, pmp, lac
real :: pmp_1, lac_1, pmp_2, lac_2
real :: da1, da2, a6da
integer :: k
logical extm (km)
grat = del (2) / del (1) ! grid ratio
bet = grat * (grat + 0.5)
q (1) = (2. * grat * (grat + 1.) * a4 (1, 1) + a4 (1, 2)) / bet
gam (1) = (1. + grat * (grat + 1.5)) / bet
do k = 2, km
d4 = del (k - 1) / del (k)
bet = 2. + 2. * d4 - gam (k - 1)
q (k) = (3. * (a4 (1, k - 1) + d4 * a4 (1, k)) - q (k - 1)) / bet
gam (k) = d4 / bet
enddo
a_bot = 1. + d4 * (d4 + 1.5)
q (km + 1) = (2. * d4 * (d4 + 1.) * a4 (1, km) + a4 (1, km - 1) - a_bot * q (km)) &
/ (d4 * (d4 + 0.5) - a_bot * gam (km))
do k = km, 1, - 1
q (k) = q (k) - gam (k) * q (k + 1)
enddo
! -----------------------------------------------------------------------
! apply constraints
! -----------------------------------------------------------------------
do k = 2, km
gam (k) = a4 (1, k) - a4 (1, k - 1)
enddo
! -----------------------------------------------------------------------
! apply large - scale constraints to all fields if not local max / min
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! top:
! -----------------------------------------------------------------------
q (1) = max (q (1), 0.)
q (2) = min (q (2), max (a4 (1, 1), a4 (1, 2)))
q (2) = max (q (2), min (a4 (1, 1), a4 (1, 2)), 0.)
! -----------------------------------------------------------------------
! interior:
! -----------------------------------------------------------------------
do k = 3, km - 1
if (gam (k - 1) * gam (k + 1) > 0.) then
q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k)))
q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k)))
else
if (gam (k - 1) > 0.) then
! there exists a local max
q (k) = max (q (k), min (a4 (1, k - 1), a4 (1, k)))
else
! there exists a local min
q (k) = min (q (k), max (a4 (1, k - 1), a4 (1, k)))
q (k) = max (q (k), 0.0)
endif
endif
enddo
! -----------------------------------------------------------------------
! bottom :
! -----------------------------------------------------------------------
q (km) = min (q (km), max (a4 (1, km - 1), a4 (1, km)))
q (km) = max (q (km), min (a4 (1, km - 1), a4 (1, km)), 0.)
! q (km + 1) = max (q (km + 1), 0.)
! -----------------------------------------------------------------------
! f (s) = al + s * [ (ar - al) + a6 * (1 - s) ] (0 <= s <= 1)
! -----------------------------------------------------------------------
do k = 1, km - 1
a4 (2, k) = q (k)
a4 (3, k) = q (k + 1)
enddo
do k = 2, km - 1
if (gam (k) * gam (k + 1) > 0.0) then
extm (k) = .false.
else
extm (k) = .true.
endif
enddo
if (do_mono) then
do k = 3, km - 2
if (extm (k)) then
! positive definite constraint only if true local extrema
if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then
a4 (2, k) = a4 (1, k)
a4 (3, k) = a4 (1, k)
endif
else
a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k))
if (abs (a4 (4, k)) > abs (a4 (2, k) - a4 (3, k))) then
! check within the smooth region if subgrid profile is non - monotonic
pmp_1 = a4 (1, k) - 2.0 * gam (k + 1)
lac_1 = pmp_1 + 1.5 * gam (k + 2)
a4 (2, k) = min (max (a4 (2, k), min (a4 (1, k), pmp_1, lac_1)), &
max (a4 (1, k), pmp_1, lac_1))
pmp_2 = a4 (1, k) + 2.0 * gam (k)
lac_2 = pmp_2 - 1.5 * gam (k - 1)
a4 (3, k) = min (max (a4 (3, k), min (a4 (1, k), pmp_2, lac_2)), &
max (a4 (1, k), pmp_2, lac_2))
endif
endif
enddo
else
do k = 3, km - 2
if (extm (k)) then
if (a4 (1, k) < qp_min .or. extm (k - 1) .or. extm (k + 1)) then
a4 (2, k) = a4 (1, k)
a4 (3, k) = a4 (1, k)
endif
endif
enddo
endif
do k = 1, km - 1
a4 (4, k) = 6. * a4 (1, k) - 3. * (a4 (2, k) + a4 (3, k))
enddo
k = km - 1
if (extm (k)) then
a4 (2, k) = a4 (1, k)
a4 (3, k) = a4 (1, k)
a4 (4, k) = 0.
else
da1 = a4 (3, k) - a4 (2, k)
da2 = da1 ** 2
a6da = a4 (4, k) * da1
if (a6da < - da2) then
a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k))
a4 (3, k) = a4 (2, k) - a4 (4, k)
elseif (a6da > da2) then
a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k))
a4 (2, k) = a4 (3, k) - a4 (4, k)
endif
endif
call cs_limiters (km - 1, a4)
! -----------------------------------------------------------------------
! bottom layer:
! -----------------------------------------------------------------------
a4 (2, km) = a4 (1, km)
a4 (3, km) = a4 (1, km)
a4 (4, km) = 0.
end subroutine cs_profile
subroutine cs_limiters (km, a4)
implicit none
integer, intent (in) :: km
real, intent (inout) :: a4 (4, km) ! ppm array
real, parameter :: r12 = 1. / 12.
integer :: k
! -----------------------------------------------------------------------
! positive definite constraint
! -----------------------------------------------------------------------
do k = 1, km
if (abs (a4 (3, k) - a4 (2, k)) < - a4 (4, k)) then
if ((a4 (1, k) + 0.25 * (a4 (3, k) - a4 (2, k)) ** 2 / a4 (4, k) + a4 (4, k) * r12) < 0.) then
if (a4 (1, k) < a4 (3, k) .and. a4 (1, k) < a4 (2, k)) then
a4 (3, k) = a4 (1, k)
a4 (2, k) = a4 (1, k)
a4 (4, k) = 0.
elseif (a4 (3, k) > a4 (2, k)) then
a4 (4, k) = 3. * (a4 (2, k) - a4 (1, k))
a4 (3, k) = a4 (2, k) - a4 (4, k)
else
a4 (4, k) = 3. * (a4 (3, k) - a4 (1, k))
a4 (2, k) = a4 (3, k) - a4 (4, k)
endif
endif
endif
enddo
end subroutine cs_limiters
! =======================================================================
! calculation of vertical fall speed
! =======================================================================
subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in), dimension (ktop:kbot) :: den, qs, qi, qg, ql, tk
real, intent (out), dimension (ktop:kbot) :: vts, vti, vtg
! fall velocity constants:
real, parameter :: thi = 1.0e-8 ! cloud ice threshold for terminal fall
real, parameter :: thg = 1.0e-8
real, parameter :: ths = 1.0e-8
real, parameter :: aa = - 4.14122e-5
real, parameter :: bb = - 0.00538922
real, parameter :: cc = - 0.0516344
real, parameter :: dd = 0.00216078
real, parameter :: ee = 1.9714
! marshall - palmer constants
real, parameter :: vcons = 6.6280504
real, parameter :: vcong = 87.2382675
real, parameter :: norms = 942477796.076938
real, parameter :: normg = 5026548245.74367
real, dimension (ktop:kbot) :: qden, tc, rhof
real :: vi0
integer :: k
! -----------------------------------------------------------------------
! marshall - palmer formula
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! try the local air density -- for global model; the true value could be
! much smaller than sfcrho over high mountains
! -----------------------------------------------------------------------
do k = ktop, kbot
rhof (k) = sqrt (min (10., sfcrho / den (k)))
enddo
! -----------------------------------------------------------------------
! ice:
! -----------------------------------------------------------------------
if (const_vi) then
vti (:) = vi_fac
else
! -----------------------------------------------------------------------
! use deng and mace (2008, grl), which gives smaller fall speed than hd90 formula
! -----------------------------------------------------------------------
vi0 = 0.01 * vi_fac
do k = ktop, kbot
if (qi (k) < thi) then ! this is needed as the fall - speed maybe problematic for small qi
vti (k) = vf_min
else
tc (k) = tk (k) - tice
vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee
vti (k) = vi0 * exp (log_10 * vti (k))
vti (k) = min (vi_max, max (vf_min, vti (k)))
endif
enddo
endif
! -----------------------------------------------------------------------
! snow:
! -----------------------------------------------------------------------
if (const_vs) then
vts (:) = vs_fac ! 1. ifs_2016
else
do k = ktop, kbot
if (qs (k) < ths) then
vts (k) = vf_min
else
vts (k) = vs_fac * vcons * rhof (k) * exp (0.0625 * log (qs (k) * den (k) / norms))
vts (k) = min (vs_max, max (vf_min, vts (k)))
endif
enddo
endif
! -----------------------------------------------------------------------
! graupel:
! -----------------------------------------------------------------------
if (const_vg) then
vtg (:) = vg_fac ! 2.
else
do k = ktop, kbot
if (qg (k) < thg) then
vtg (k) = vf_min
else
vtg (k) = vg_fac * vcong * rhof (k) * sqrt (sqrt (sqrt (qg (k) * den (k) / normg)))
vtg (k) = min (vg_max, max (vf_min, vtg (k)))
endif
enddo
endif
end subroutine fall_speed
! =======================================================================
! setup gfdl cloud microphysics parameters
! =======================================================================
subroutine setupm
implicit none
real :: gcon, cd, scm3, pisq, act (8)
real :: vdifu, tcond
real :: visk
real :: ch2o, hltf
real :: hlts, hltc, ri50
real, parameter :: gam263 = 1.456943, gam275 = 1.608355, gam290 = 1.827363, &
gam325 = 2.54925, gam350 = 3.323363, gam380 = 4.694155, &
gam425 = 8.285063, gam450 = 11.631769, gam480 = 17.837789, &
gam625 = 184.860962, gam680 = 496.604067
! intercept parameters
real, parameter :: rnzr = 8.0e6 ! lin83
real, parameter :: rnzs = 3.0e6 ! lin83
real, parameter :: rnzg = 4.0e6 ! rh84
! density parameters
real, parameter :: rhos = 0.1e3 ! lin83 (snow density; 1 / 10 of water)
real, parameter :: rhog = 0.4e3 ! rh84 (graupel density)
real, parameter :: acc (3) = (/ 5.0, 2.0, 0.5 /)
real den_rc
integer :: i, k
pie = 4. * atan (1.0)
! s. klein's formular (eq 16) from am2
fac_rc = (4. / 3.) * pie * rhor * rthresh ** 3
if (prog_ccn) then
! if (master) write (*, *) 'prog_ccn option is .t.'
else
den_rc = fac_rc * ccn_o * 1.e6
! if (master) write (*, *) 'mp: for ccn_o = ', ccn_o, 'ql_rc = ', den_rc
den_rc = fac_rc * ccn_l * 1.e6
! if (master) write (*, *) 'mp: for ccn_l = ', ccn_l, 'ql_rc = ', den_rc
endif
vdifu = 2.11e-5
tcond = 2.36e-2
visk = 1.259e-5
hlts = 2.8336e6
hltc = 2.5e6
hltf = 3.336e5
ch2o = 4.1855e3
ri50 = 1.e-4
pisq = pie * pie
scm3 = (visk / vdifu) ** (1. / 3.)
cracs = pisq * rnzr * rnzs * rhos
csacr = pisq * rnzr * rnzs * rhor
cgacr = pisq * rnzr * rnzg * rhor
cgacs = pisq * rnzg * rnzs * rhos
cgacs = cgacs * c_pgacs
! act: 1 - 2:racs (s - r) ; 3 - 4:sacr (r - s) ;
! 5 - 6:gacr (r - g) ; 7 - 8:gacs (s - g)
act (1) = pie * rnzs * rhos
act (2) = pie * rnzr * rhor
act (6) = pie * rnzg * rhog
act (3) = act (2)
act (4) = act (1)
act (5) = act (2)
act (7) = act (1)
act (8) = act (6)
do i = 1, 3
do k = 1, 4
acco (i, k) = acc (i) / (act (2 * k - 1) ** ((7 - i) * 0.25) * act (2 * k) ** (i * 0.25))
enddo
enddo
gcon = 40.74 * sqrt (sfcrho) ! 44.628
csacw = pie * rnzs * clin * gam325 / (4. * act (1) ** 0.8125)
! decreasing csacw to reduce cloud water --- > snow
craci = pie * rnzr * alin * gam380 / (4. * act (2) ** 0.95)
csaci = csacw * c_psaci
cgacw = pie * rnzg * gam350 * gcon / (4. * act (6) ** 0.875)
! cgaci = cgacw * 0.1
! sjl, may 28, 2012
cgaci = cgacw * 0.05
! sjl, may 28, 2012
cracw = craci ! cracw = 3.27206196043822
cracw = c_cracw * cracw
! subl and revp: five constants for three separate processes
cssub (1) = 2. * pie * vdifu * tcond * rvgas * rnzs
cgsub (1) = 2. * pie * vdifu * tcond * rvgas * rnzg
crevp (1) = 2. * pie * vdifu * tcond * rvgas * rnzr
cssub (2) = 0.78 / sqrt (act (1))
cgsub (2) = 0.78 / sqrt (act (6))
crevp (2) = 0.78 / sqrt (act (2))
cssub (3) = 0.31 * scm3 * gam263 * sqrt (clin / visk) / act (1) ** 0.65625
cgsub (3) = 0.31 * scm3 * gam275 * sqrt (gcon / visk) / act (6) ** 0.6875
crevp (3) = 0.31 * scm3 * gam290 * sqrt (alin / visk) / act (2) ** 0.725
cssub (4) = tcond * rvgas
cssub (5) = hlts ** 2 * vdifu
cgsub (4) = cssub (4)
crevp (4) = cssub (4)
cgsub (5) = cssub (5)
crevp (5) = hltc ** 2 * vdifu
cgfr (1) = 20.e2 * pisq * rnzr * rhor / act (2) ** 1.75
cgfr (2) = 0.66
! smlt: five constants (lin et al. 1983)
csmlt (1) = 2. * pie * tcond * rnzs / hltf
csmlt (2) = 2. * pie * vdifu * rnzs * hltc / hltf
csmlt (3) = cssub (2)
csmlt (4) = cssub (3)
csmlt (5) = ch2o / hltf
! gmlt: five constants
cgmlt (1) = 2. * pie * tcond * rnzg / hltf
cgmlt (2) = 2. * pie * vdifu * rnzg * hltc / hltf
cgmlt (3) = cgsub (2)
cgmlt (4) = cgsub (3)
cgmlt (5) = ch2o / hltf
es0 = 6.107799961e2 ! ~6.1 mb
ces0 = eps * es0
end subroutine setupm
! =======================================================================
! initialization of gfdl cloud microphysics
! =======================================================================
!subroutine gfdl_cloud_microphys_init (id, jd, kd, axes, time)
subroutine gfdl_cloud_microphys_init (me, master, nlunit, logunit, fn_nml)
implicit none
integer, intent (in) :: me
integer, intent (in) :: master
integer, intent (in) :: nlunit
integer, intent (in) :: logunit
character (len = 64), intent (in) :: fn_nml
integer :: ios
logical :: exists
! integer, intent (in) :: id, jd, kd
! integer, intent (in) :: axes (4)
! type (time_type), intent (in) :: time
! integer :: unit, io, ierr, k, logunit
! logical :: flag
! real :: tmp, q1, q2
! master = (mpp_pe () .eq.mpp_root_pe ())
!#ifdef internal_file_nml
! read (input_nml_file, nml = gfdl_cloud_microphys_nml, iostat = io)
! ierr = check_nml_error (io, 'gfdl_cloud_microphys_nml')
!#else
! if (file_exist ('input.nml')) then
! unit = open_namelist_file ()
! io = 1
! do while (io .ne. 0)
! read (unit, nml = gfdl_cloud_microphys_nml, iostat = io, end = 10)
! ierr = check_nml_error (io, 'gfdl_cloud_microphys_nml')
! enddo
!10 call close_file (unit)
! endif
!#endif
! call write_version_number ('gfdl_cloud_microphys_mod', version)
! logunit = stdlog ()
inquire (file = trim (fn_nml), exist = exists)
if (.not. exists) then
write (6, *) 'gfdl - mp :: namelist file: ', trim (fn_nml), ' does not exist'
stop
else
open (unit = nlunit, file = fn_nml, readonly, status = 'old', iostat = ios)
endif
rewind (nlunit)
read (nlunit, nml = gfdl_cloud_microphysics_nml)
close (nlunit)
! write version number and namelist to log file
if (me == master) then
write (logunit, *) " ================================================================== "
write (logunit, *) "gfdl_cloud_microphys_mod"
write (logunit, nml = gfdl_cloud_microphysics_nml)
endif
if (do_setup) then
call setup_con
call setupm
do_setup = .false.
endif
log_10 = log (10.)
tice0 = tice - 0.01
t_wfr = tice - 40.0 ! supercooled water can exist down to - 48 c, which is the "absolute"
! if (master) write (logunit, nml = gfdl_cloud_microphys_nml)
!
! id_vtr = register_diag_field (mod_name, 'vt_r', axes (1:3), time, &
! 'rain fall speed', 'm / s', missing_value = missing_value)
! id_vts = register_diag_field (mod_name, 'vt_s', axes (1:3), time, &
! 'snow fall speed', 'm / s', missing_value = missing_value)
! id_vtg = register_diag_field (mod_name, 'vt_g', axes (1:3), time, &
! 'graupel fall speed', 'm / s', missing_value = missing_value)
! id_vti = register_diag_field (mod_name, 'vt_i', axes (1:3), time, &
! 'ice fall speed', 'm / s', missing_value = missing_value)
! id_droplets = register_diag_field (mod_name, 'droplets', axes (1:3), time, &
! 'droplet number concentration', '# / m3', missing_value = missing_value)
! id_rh = register_diag_field (mod_name, 'rh_lin', axes (1:2), time, &
! 'relative humidity', 'n / a', missing_value = missing_value)
! id_rain = register_diag_field (mod_name, 'rain_lin', axes (1:2), time, &
! 'rain_lin', 'mm / day', missing_value = missing_value)
! id_snow = register_diag_field (mod_name, 'snow_lin', axes (1:2), time, &
! 'snow_lin', 'mm / day', missing_value = missing_value)
! id_graupel = register_diag_field (mod_name, 'graupel_lin', axes (1:2), time, &
! 'graupel_lin', 'mm / day', missing_value = missing_value)
! id_ice = register_diag_field (mod_name, 'ice_lin', axes (1:2), time, &
! 'ice_lin', 'mm / day', missing_value = missing_value)
! id_prec = register_diag_field (mod_name, 'prec_lin', axes (1:2), time, &
! 'prec_lin', 'mm / day', missing_value = missing_value)
! if (master) write (*, *) 'prec_lin diagnostics initialized.', id_prec
! id_cond = register_diag_field (mod_name, 'cond_lin', axes (1:2), time, &
! 'total condensate', 'kg / m ** 2', missing_value = missing_value)
! id_var = register_diag_field (mod_name, 'var_lin', axes (1:2), time, &
! 'subgrid variance', 'n / a', missing_value = missing_value)
! call qsmith_init
! testing the water vapor tables
! if (mp_debug .and. master) then
! write (*, *) 'testing water vapor tables in gfdl_cloud_microphys'
! tmp = tice - 90.
! do k = 1, 25
! q1 = wqsat_moist (tmp, 0., 1.e5)
! q2 = qs1d_m (tmp, 0., 1.e5)
! write (*, *) nint (tmp - tice), q1, q2, 'dq = ', q1 - q2
! tmp = tmp + 5.
! enddo
! endif
! if (master) write (*, *) 'gfdl_cloud_micrphys diagnostics initialized.'
! gfdl_mp_clock = mpp_clock_id ('gfdl_cloud_microphys', grain = clock_routine)
module_is_initialized = .true.
end subroutine gfdl_cloud_microphys_init
! =======================================================================
! end of gfdl cloud microphysics
! =======================================================================
subroutine gfdl_cloud_microphys_end
implicit none
deallocate (table)
deallocate (table2)
deallocate (table3)
deallocate (tablew)
deallocate (des)
deallocate (des2)
deallocate (des3)
deallocate (desw)
tables_are_initialized = .false.
end subroutine gfdl_cloud_microphys_end
! =======================================================================
! qsmith table initialization
! =======================================================================
subroutine setup_con
implicit none
! master = (mpp_pe () .eq.mpp_root_pe ())
rgrav = 1. / grav
if (.not. qsmith_tables_initialized) call qsmith_init
qsmith_tables_initialized = .true.
end subroutine setup_con
! =======================================================================
! accretion function (lin et al. 1983)
! =======================================================================
real function acr3d (v1, v2, q1, q2, c, cac, rho)
implicit none
real, intent (in) :: v1, v2, c, rho
real, intent (in) :: q1, q2 ! mixing ratio!!!
real, intent (in) :: cac (3)
real :: t1, s1, s2
! integer :: k
!
! real :: a
!
! a = 0.0
! do k = 1, 3
! a = a + cac (k) * ((q1 * rho) ** ((7 - k) * 0.25) * (q2 * rho) ** (k * 0.25))
! enddo
! acr3d = c * abs (v1 - v2) * a / rho
! optimized
t1 = sqrt (q1 * rho)
s1 = sqrt (q2 * rho)
s2 = sqrt (s1) ! s1 = s2 ** 2
acr3d = c * abs (v1 - v2) * q1 * s2 * (cac (1) * t1 + cac (2) * sqrt (t1) * s2 + cac (3) * s1)
end function acr3d
! =======================================================================
! melting of snow function (lin et al. 1983)
! note: psacw and psacr must be calc before smlt is called
! =======================================================================
real function smlt (tc, dqs, qsrho, psacw, psacr, c, rho, rhofac)
implicit none
real, intent (in) :: tc, dqs, qsrho, psacw, psacr, c (5), rho, rhofac
smlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qsrho) + &
c (4) * qsrho ** 0.65625 * sqrt (rhofac)) + c (5) * tc * (psacw + psacr)
end function smlt
! =======================================================================
! melting of graupel function (lin et al. 1983)
! note: pgacw and pgacr must be calc before gmlt is called
! =======================================================================
real function gmlt (tc, dqs, qgrho, pgacw, pgacr, c, rho)
implicit none
real, intent (in) :: tc, dqs, qgrho, pgacw, pgacr, c (5), rho
gmlt = (c (1) * tc / rho - c (2) * dqs) * (c (3) * sqrt (qgrho) + &
c (4) * qgrho ** 0.6875 / rho ** 0.25) + c (5) * tc * (pgacw + pgacr)
end function gmlt
! =======================================================================
! initialization
! prepare saturation water vapor pressure tables
! =======================================================================
subroutine qsmith_init
implicit none
integer, parameter :: length = 2621
integer :: i
if (.not. tables_are_initialized) then
! master = (mpp_pe () .eq. mpp_root_pe ())
! if (master) print *, ' gfdl mp: initializing qs tables'
! debug code
! print *, mpp_pe (), allocated (table), allocated (table2), &
! allocated (table3), allocated (tablew), allocated (des), &
! allocated (des2), allocated (des3), allocated (desw)
! end debug code
! generate es table (dt = 0.1 deg. c)
allocate (table (length))
allocate (table2 (length))
allocate (table3 (length))
allocate (tablew (length))
allocate (des (length))
allocate (des2 (length))
allocate (des3 (length))
allocate (desw (length))
call qs_table (length)
call qs_table2 (length)
call qs_table3 (length)
call qs_tablew (length)
do i = 1, length - 1
des (i) = max (0., table (i + 1) - table (i))
des2 (i) = max (0., table2 (i + 1) - table2 (i))
des3 (i) = max (0., table3 (i + 1) - table3 (i))
desw (i) = max (0., tablew (i + 1) - tablew (i))
enddo
des (length) = des (length - 1)
des2 (length) = des2 (length - 1)
des3 (length) = des3 (length - 1)
desw (length) = desw (length - 1)
tables_are_initialized = .true.
endif
end subroutine qsmith_init
! =======================================================================
! compute the saturated specific humidity for table ii
! =======================================================================
real function wqs1 (ta, den)
implicit none
! pure water phase; universal dry / moist formular using air density
! input "den" can be either dry or moist air density
real, intent (in) :: ta, den
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = tablew (it) + (ap1 - it) * desw (it)
wqs1 = es / (rvgas * ta * den)
end function wqs1
! =======================================================================
! compute the gradient of saturated specific humidity for table ii
! =======================================================================
real function wqs2 (ta, den, dqdt)
implicit none
! pure water phase; universal dry / moist formular using air density
! input "den" can be either dry or moist air density
real, intent (in) :: ta, den
real, intent (out) :: dqdt
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
if (.not. tables_are_initialized) call qsmith_init
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = tablew (it) + (ap1 - it) * desw (it)
wqs2 = es / (rvgas * ta * den)
it = ap1 - 0.5
! finite diff, del_t = 0.1:
dqdt = 10. * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) / (rvgas * ta * den)
end function wqs2
! =======================================================================
! compute wet buld temperature
! =======================================================================
real function wet_bulb (q, t, den)
implicit none
real, intent (in) :: t, q, den
real :: qs, tp, dqdt
wet_bulb = t
qs = wqs2 (wet_bulb, den, dqdt)
tp = 0.5 * (qs - q) / (1. + lcp * dqdt) * lcp
wet_bulb = wet_bulb - tp
! tp is negative if super - saturated
if (tp > 0.01) then
qs = wqs2 (wet_bulb, den, dqdt)
tp = (qs - q) / (1. + lcp * dqdt) * lcp
wet_bulb = wet_bulb - tp
endif
end function wet_bulb
! =======================================================================
! compute the saturated specific humidity for table iii
! =======================================================================
real function iqs1 (ta, den)
implicit none
! water - ice phase; universal dry / moist formular using air density
! input "den" can be either dry or moist air density
real, intent (in) :: ta, den
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = table2 (it) + (ap1 - it) * des2 (it)
iqs1 = es / (rvgas * ta * den)
end function iqs1
! =======================================================================
! compute the gradient of saturated specific humidity for table iii
! =======================================================================
real function iqs2 (ta, den, dqdt)
implicit none
! water - ice phase; universal dry / moist formular using air density
! input "den" can be either dry or moist air density
real, intent (in) :: ta, den
real, intent (out) :: dqdt
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = table2 (it) + (ap1 - it) * des2 (it)
iqs2 = es / (rvgas * ta * den)
it = ap1 - 0.5
dqdt = 10. * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) / (rvgas * ta * den)
end function iqs2
! =======================================================================
! compute the gradient of saturated specific humidity for table iii
! =======================================================================
real function qs1d_moist (ta, qv, pa, dqdt)
implicit none
real, intent (in) :: ta, pa, qv
real, intent (out) :: dqdt
real :: es, ap1, tmin, eps10
integer :: it
tmin = table_ice - 160.
eps10 = 10. * eps
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = table2 (it) + (ap1 - it) * des2 (it)
qs1d_moist = eps * es * (1. + zvir * qv) / pa
it = ap1 - 0.5
dqdt = eps10 * (des2 (it) + (ap1 - it) * (des2 (it + 1) - des2 (it))) * (1. + zvir * qv) / pa
end function qs1d_moist
! =======================================================================
! compute the gradient of saturated specific humidity for table ii
! =======================================================================
real function wqsat2_moist (ta, qv, pa, dqdt)
implicit none
real, intent (in) :: ta, pa, qv
real, intent (out) :: dqdt
real :: es, ap1, tmin, eps10
integer :: it
tmin = table_ice - 160.
eps10 = 10. * eps
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = tablew (it) + (ap1 - it) * desw (it)
wqsat2_moist = eps * es * (1. + zvir * qv) / pa
it = ap1 - 0.5
dqdt = eps10 * (desw (it) + (ap1 - it) * (desw (it + 1) - desw (it))) * (1. + zvir * qv) / pa
end function wqsat2_moist
! =======================================================================
! compute the saturated specific humidity for table ii
! =======================================================================
real function wqsat_moist (ta, qv, pa)
implicit none
real, intent (in) :: ta, pa, qv
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = tablew (it) + (ap1 - it) * desw (it)
wqsat_moist = eps * es * (1. + zvir * qv) / pa
end function wqsat_moist
! =======================================================================
! compute the saturated specific humidity for table iii
! =======================================================================
real function qs1d_m (ta, qv, pa)
implicit none
real, intent (in) :: ta, pa, qv
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = table2 (it) + (ap1 - it) * des2 (it)
qs1d_m = eps * es * (1. + zvir * qv) / pa
end function qs1d_m
! =======================================================================
! computes the difference in saturation vapor * density * between water and ice
! =======================================================================
real function d_sat (ta, den)
implicit none
real, intent (in) :: ta, den
real :: es_w, es_i, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es_w = tablew (it) + (ap1 - it) * desw (it)
es_i = table2 (it) + (ap1 - it) * des2 (it)
d_sat = dim (es_w, es_i) / (rvgas * ta * den) ! take positive difference
end function d_sat
! =======================================================================
! compute the saturated water vapor pressure for table ii
! =======================================================================
real function esw_table (ta)
implicit none
real, intent (in) :: ta
real :: ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
esw_table = tablew (it) + (ap1 - it) * desw (it)
end function esw_table
! =======================================================================
! compute the saturated water vapor pressure for table iii
! =======================================================================
real function es2_table (ta)
implicit none
real, intent (in) :: ta
real :: ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (ta, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es2_table = table2 (it) + (ap1 - it) * des2 (it)
end function es2_table
! =======================================================================
! compute the saturated water vapor pressure for table ii
! =======================================================================
subroutine esw_table1d (ta, es, n)
implicit none
integer, intent (in) :: n
real, intent (in) :: ta (n)
real, intent (out) :: es (n)
real :: ap1, tmin
integer :: i, it
tmin = table_ice - 160.
do i = 1, n
ap1 = 10. * dim (ta (i), tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es (i) = tablew (it) + (ap1 - it) * desw (it)
enddo
end subroutine esw_table1d
! =======================================================================
! compute the saturated water vapor pressure for table iii
! =======================================================================
subroutine es2_table1d (ta, es, n)
implicit none
integer, intent (in) :: n
real, intent (in) :: ta (n)
real, intent (out) :: es (n)
real :: ap1, tmin
integer :: i, it
tmin = table_ice - 160.
do i = 1, n
ap1 = 10. * dim (ta (i), tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es (i) = table2 (it) + (ap1 - it) * des2 (it)
enddo
end subroutine es2_table1d
! =======================================================================
! compute the saturated water vapor pressure for table iv
! =======================================================================
subroutine es3_table1d (ta, es, n)
implicit none
integer, intent (in) :: n
real, intent (in) :: ta (n)
real, intent (out) :: es (n)
real :: ap1, tmin
integer :: i, it
tmin = table_ice - 160.
do i = 1, n
ap1 = 10. * dim (ta (i), tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es (i) = table3 (it) + (ap1 - it) * des3 (it)
enddo
end subroutine es3_table1d
! =======================================================================
! saturation water vapor pressure table ii
! 1 - phase table
! =======================================================================
subroutine qs_tablew (n)
implicit none
integer, intent (in) :: n
real :: delt = 0.1
real :: tmin, tem, fac0, fac1, fac2
integer :: i
tmin = table_ice - 160.
! -----------------------------------------------------------------------
! compute es over water
! -----------------------------------------------------------------------
do i = 1, n
tem = tmin + delt * real (i - 1)
fac0 = (tem - t_ice) / (tem * t_ice)
fac1 = fac0 * lv0
fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas
tablew (i) = e00 * exp (fac2)
enddo
end subroutine qs_tablew
! =======================================================================
! saturation water vapor pressure table iii
! 2 - phase table
! =======================================================================
subroutine qs_table2 (n)
implicit none
integer, intent (in) :: n
real :: delt = 0.1
real :: tmin, tem0, tem1, fac0, fac1, fac2
integer :: i, i0, i1
tmin = table_ice - 160.
do i = 1, n
tem0 = tmin + delt * real (i - 1)
fac0 = (tem0 - t_ice) / (tem0 * t_ice)
if (i <= 1600) then
! -----------------------------------------------------------------------
! compute es over ice between - 160 deg c and 0 deg c.
! -----------------------------------------------------------------------
fac1 = fac0 * li2
fac2 = (d2ice * log (tem0 / t_ice) + fac1) / rvgas
else
! -----------------------------------------------------------------------
! compute es over water between 0 deg c and 102 deg c.
! -----------------------------------------------------------------------
fac1 = fac0 * lv0
fac2 = (dc_vap * log (tem0 / t_ice) + fac1) / rvgas
endif
table2 (i) = e00 * exp (fac2)
enddo
! -----------------------------------------------------------------------
! smoother around 0 deg c
! -----------------------------------------------------------------------
i0 = 1600
i1 = 1601
tem0 = 0.25 * (table2 (i0 - 1) + 2. * table (i0) + table2 (i0 + 1))
tem1 = 0.25 * (table2 (i1 - 1) + 2. * table (i1) + table2 (i1 + 1))
table2 (i0) = tem0
table2 (i1) = tem1
end subroutine qs_table2
! =======================================================================
! saturation water vapor pressure table iv
! 2 - phase table with " - 2 c" as the transition point
! =======================================================================
subroutine qs_table3 (n)
implicit none
integer, intent (in) :: n
real :: delt = 0.1
real :: esbasw, tbasw, esbasi, tmin, tem, aa, b, c, d, e
real :: tem0, tem1
integer :: i, i0, i1
esbasw = 1013246.0
tbasw = table_ice + 100.
esbasi = 6107.1
tmin = table_ice - 160.
do i = 1, n
tem = tmin + delt * real (i - 1)
! if (i <= 1600) then
if (i <= 1580) then ! change to - 2 c
! -----------------------------------------------------------------------
! compute es over ice between - 160 deg c and 0 deg c.
! see smithsonian meteorological tables page 350.
! -----------------------------------------------------------------------
aa = - 9.09718 * (table_ice / tem - 1.)
b = - 3.56654 * alog10 (table_ice / tem)
c = 0.876793 * (1. - tem / table_ice)
e = alog10 (esbasi)
table3 (i) = 0.1 * 10 ** (aa + b + c + e)
else
! -----------------------------------------------------------------------
! compute es over water between - 2 deg c and 102 deg c.
! see smithsonian meteorological tables page 350.
! -----------------------------------------------------------------------
aa = - 7.90298 * (tbasw / tem - 1.)
b = 5.02808 * alog10 (tbasw / tem)
c = - 1.3816e-7 * (10 ** ((1. - tem / tbasw) * 11.344) - 1.)
d = 8.1328e-3 * (10 ** ((tbasw / tem - 1.) * (- 3.49149)) - 1.)
e = alog10 (esbasw)
table3 (i) = 0.1 * 10 ** (aa + b + c + d + e)
endif
enddo
! -----------------------------------------------------------------------
! smoother around - 2 deg c
! -----------------------------------------------------------------------
i0 = 1580
i1 = 1581
tem0 = 0.25 * (table3 (i0 - 1) + 2. * table (i0) + table3 (i0 + 1))
tem1 = 0.25 * (table3 (i1 - 1) + 2. * table (i1) + table3 (i1 + 1))
table3 (i0) = tem0
table3 (i1) = tem1
end subroutine qs_table3
! =======================================================================
! compute the saturated specific humidity for table
! note: this routine is based on "moist" mixing ratio
! =======================================================================
real function qs_blend (t, p, q)
implicit none
real, intent (in) :: t, p, q
real :: es, ap1, tmin
integer :: it
tmin = table_ice - 160.
ap1 = 10. * dim (t, tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es = table (it) + (ap1 - it) * des (it)
qs_blend = eps * es * (1. + zvir * q) / p
end function qs_blend
! =======================================================================
! saturation water vapor pressure table i
! 3 - phase table
! =======================================================================
subroutine qs_table (n)
implicit none
integer, intent (in) :: n
real :: delt = 0.1
real :: tmin, tem, esh20
real :: wice, wh2o, fac0, fac1, fac2
real :: esupc (200)
integer :: i
tmin = table_ice - 160.
! -----------------------------------------------------------------------
! compute es over ice between - 160 deg c and 0 deg c.
! -----------------------------------------------------------------------
do i = 1, 1600
tem = tmin + delt * real (i - 1)
fac0 = (tem - t_ice) / (tem * t_ice)
fac1 = fac0 * li2
fac2 = (d2ice * log (tem / t_ice) + fac1) / rvgas
table (i) = e00 * exp (fac2)
enddo
! -----------------------------------------------------------------------
! compute es over water between - 20 deg c and 102 deg c.
! -----------------------------------------------------------------------
do i = 1, 1221
tem = 253.16 + delt * real (i - 1)
fac0 = (tem - t_ice) / (tem * t_ice)
fac1 = fac0 * lv0
fac2 = (dc_vap * log (tem / t_ice) + fac1) / rvgas
esh20 = e00 * exp (fac2)
if (i <= 200) then
esupc (i) = esh20
else
table (i + 1400) = esh20
endif
enddo
! -----------------------------------------------------------------------
! derive blended es over ice and supercooled water between - 20 deg c and 0 deg c
! -----------------------------------------------------------------------
do i = 1, 200
tem = 253.16 + delt * real (i - 1)
wice = 0.05 * (table_ice - tem)
wh2o = 0.05 * (tem - 253.16)
table (i + 1400) = wice * table (i + 1400) + wh2o * esupc (i)
enddo
end subroutine qs_table
! =======================================================================
! compute the saturated specific humidity and the gradient of saturated specific humidity
! input t in deg k, p in pa; p = rho rdry tv, moist pressure
! =======================================================================
subroutine qsmith (im, km, ks, t, p, q, qs, dqdt)
implicit none
integer, intent (in) :: im, km, ks
real, intent (in), dimension (im, km) :: t, p, q
real, intent (out), dimension (im, km) :: qs
real, intent (out), dimension (im, km), optional :: dqdt
real :: eps10, ap1, tmin
real, dimension (im, km) :: es
integer :: i, k, it
tmin = table_ice - 160.
eps10 = 10. * eps
if (.not. tables_are_initialized) then
call qsmith_init
endif
do k = ks, km
do i = 1, im
ap1 = 10. * dim (t (i, k), tmin) + 1.
ap1 = min (2621., ap1)
it = ap1
es (i, k) = table (it) + (ap1 - it) * des (it)
qs (i, k) = eps * es (i, k) * (1. + zvir * q (i, k)) / p (i, k)
enddo
enddo
if (present (dqdt)) then
do k = ks, km
do i = 1, im
ap1 = 10. * dim (t (i, k), tmin) + 1.
ap1 = min (2621., ap1) - 0.5
it = ap1
dqdt (i, k) = eps10 * (des (it) + (ap1 - it) * (des (it + 1) - des (it))) * (1. + zvir * q (i, k)) / p (i, k)
enddo
enddo
endif
end subroutine qsmith
! =======================================================================
! fix negative water species
! this is designed for 6 - class micro - physics schemes
! =======================================================================
subroutine neg_adj (ktop, kbot, pt, dp, qv, ql, qr, qi, qs, qg)
implicit none
integer, intent (in) :: ktop, kbot
real, intent (in), dimension (ktop:kbot) :: dp
real, intent (inout), dimension (ktop:kbot) :: pt, qv, ql, qr, qi, qs, qg
real, dimension (ktop:kbot) :: lcpk, icpk
real :: dq, cvm
integer :: k
! -----------------------------------------------------------------------
! define heat capacity and latent heat coefficient
! -----------------------------------------------------------------------
do k = ktop, kbot
cvm = c_air + qv (k) * c_vap + (qr (k) + ql (k)) * c_liq + (qi (k) + qs (k) + qg (k)) * c_ice
lcpk (k) = (lv00 + d0_vap * pt (k)) / cvm
icpk (k) = (li00 + dc_ice * pt (k)) / cvm
enddo
do k = ktop, kbot
! -----------------------------------------------------------------------
! ice phase:
! -----------------------------------------------------------------------
! if cloud ice < 0, borrow from snow
if (qi (k) < 0.) then
qs (k) = qs (k) + qi (k)
qi (k) = 0.
endif
! if snow < 0, borrow from graupel
if (qs (k) < 0.) then
qg (k) = qg (k) + qs (k)
qs (k) = 0.
endif
! if graupel < 0, borrow from rain
if (qg (k) < 0.) then
qr (k) = qr (k) + qg (k)
pt (k) = pt (k) - qg (k) * icpk (k) ! heating
qg (k) = 0.
endif
! -----------------------------------------------------------------------
! liquid phase:
! -----------------------------------------------------------------------
! if rain < 0, borrow from cloud water
if (qr (k) < 0.) then
ql (k) = ql (k) + qr (k)
qr (k) = 0.
endif
! if cloud water < 0, borrow from water vapor
if (ql (k) < 0.) then
qv (k) = qv (k) + ql (k)
pt (k) = pt (k) - ql (k) * lcpk (k) ! heating
ql (k) = 0.
endif
enddo
! -----------------------------------------------------------------------
! fix water vapor; borrow from below
! -----------------------------------------------------------------------
do k = ktop, kbot - 1
if (qv (k) < 0.) then
qv (k + 1) = qv (k + 1) + qv (k) * dp (k) / dp (k + 1)
qv (k) = 0.
endif
enddo
! -----------------------------------------------------------------------
! bottom layer; borrow from above
! -----------------------------------------------------------------------
if (qv (kbot) < 0. .and. qv (kbot - 1) > 0.) then
dq = min (- qv (kbot) * dp (kbot), qv (kbot - 1) * dp (kbot - 1))
qv (kbot - 1) = qv (kbot - 1) - dq / dp (kbot - 1)
qv (kbot) = qv (kbot) + dq / dp (kbot)
endif
end subroutine neg_adj
! =======================================================================
! compute global sum
! quick local sum algorithm
! =======================================================================
!real function g_sum (p, ifirst, ilast, jfirst, jlast, area, mode)
!
! use mpp_mod, only: mpp_sum
!
! implicit none
!
! integer, intent (in) :: ifirst, ilast, jfirst, jlast
! integer, intent (in) :: mode ! if == 1 divided by area
!
! real, intent (in), dimension (ifirst:ilast, jfirst:jlast) :: p, area
!
! integer :: i, j
!
! real :: gsum
!
! if (global_area < 0.) then
! global_area = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! global_area = global_area + area (i, j)
! enddo
! enddo
! call mpp_sum (global_area)
! endif
!
! gsum = 0.
! do j = jfirst, jlast
! do i = ifirst, ilast
! gsum = gsum + p (i, j) * area (i, j)
! enddo
! enddo
! call mpp_sum (gsum)
!
! if (mode == 1) then
! g_sum = gsum / global_area
! else
! g_sum = gsum
! endif
!
!end function g_sum
! =======================================================================
! interpolate to a prescribed height
! =======================================================================
subroutine interpolate_z (is, ie, js, je, km, zl, hgt, a3, a2)
implicit none
integer, intent (in) :: is, ie, js, je, km
real, intent (in), dimension (is:ie, js:je, km) :: a3
real, intent (in), dimension (is:ie, js:je, km + 1) :: hgt ! hgt (k) > hgt (k + 1)
real, intent (in) :: zl
real, intent (out), dimension (is:ie, js:je) :: a2
real, dimension (km) :: zm ! middle layer height
integer :: i, j, k
!$omp parallel do default (none) shared (is, ie, js, je, km, hgt, zl, a2, a3) private (zm)
do j = js, je
do i = is, ie
do k = 1, km
zm (k) = 0.5 * (hgt (i, j, k) + hgt (i, j, k + 1))
enddo
if (zl >= zm (1)) then
a2 (i, j) = a3 (i, j, 1)
elseif (zl <= zm (km)) then
a2 (i, j) = a3 (i, j, km)
else
do k = 1, km - 1
if (zl <= zm (k) .and. zl >= zm (k + 1)) then
a2 (i, j) = a3 (i, j, k) + (a3 (i, j, k + 1) - a3 (i, j, k)) * (zm (k) - zl) / (zm (k) - zm (k + 1))
exit
endif
enddo
endif
enddo
enddo
end subroutine interpolate_z
! =======================================================================
! radius of cloud species diagnosis
! =======================================================================
subroutine cloud_diagnosis (is, ie, js, je, den, qw, qi, qr, qs, qg, t, &
qcw, qci, qcr, qcs, qcg, rew, rei, rer, res, reg)
implicit none
integer, intent (in) :: is, ie, js, je
real, intent (in), dimension (is:ie, js:je) :: den, t
real, intent (in), dimension (is:ie, js:je) :: qw, qi, qr, qs, qg ! units: kg / kg
real, intent (out), dimension (is:ie, js:je) :: qcw, qci, qcr, qcs, qcg ! units: kg / m^3
real, intent (out), dimension (is:ie, js:je) :: rew, rei, rer, res, reg ! units: micron
integer :: i, j
real :: lambdar, lambdas, lambdag
real :: rhow = 1.0e3, rhor = 1.0e3, rhos = 1.0e2, rhog = 4.0e2
real :: n0r = 8.0e6, n0s = 3.0e6, n0g = 4.0e6
real :: alphar = 0.8, alphas = 0.25, alphag = 0.5
real :: gammar = 17.837789, gammas = 8.2850630, gammag = 11.631769
real :: qmin = 1.0e-5, ccn = 1.0e8, beta = 1.22
! real :: rewmin = 1.0, rewmax = 25.0
! real :: reimin = 10.0, reimax = 300.0
! real :: rermin = 25.0, rermax = 225.0
! real :: resmin = 300, resmax = 1000.0
! real :: regmin = 1000.0, regmax = 1.0e5
real :: rewmin = 5.0, rewmax = 10.0
real :: reimin = 10.0, reimax = 150.0
real :: rermin = 0.0, rermax = 10000.0
real :: resmin = 0.0, resmax = 10000.0
real :: regmin = 0.0, regmax = 10000.0
do j = js, je
do i = is, ie
! -----------------------------------------------------------------------
! cloud water (martin et al., 1994)
! -----------------------------------------------------------------------
if (qw (i, j) .gt. qmin) then
qcw (i, j) = den (i, j) * qw (i, j)
rew (i, j) = exp (1.0 / 3.0 * log ((3 * qcw (i, j)) / (4 * pi * rhow * ccn))) * 1.0e6
rew (i, j) = max (rewmin, min (rewmax, rew (i, j)))
else
qcw (i, j) = 0.0
rew (i, j) = rewmin
endif
! -----------------------------------------------------------------------
! cloud ice (heymsfield and mcfarquhar, 1996)
! -----------------------------------------------------------------------
if (qi (i, j) .gt. qmin) then
qci (i, j) = den (i, j) * qi (i, j)
if (t (i, j) - tice .lt. - 50) then
rei (i, j) = beta / 9.917 * exp ((1 - 0.891) * log (1.0e3 * qci (i, j))) * 1.0e3
elseif (t (i, j) - tice .lt. - 40) then
rei (i, j) = beta / 9.337 * exp ((1 - 0.920) * log (1.0e3 * qci (i, j))) * 1.0e3
elseif (t (i, j) - tice .lt. - 30) then
rei (i, j) = beta / 9.208 * exp ((1 - 0.945) * log (1.0e3 * qci (i, j))) * 1.0e3
else
rei (i, j) = beta / 9.387 * exp ((1 - 0.969) * log (1.0e3 * qci (i, j))) * 1.0e3
endif
rei (i, j) = max (reimin, min (reimax, rei (i, j)))
else
qci (i, j) = 0.0
rei (i, j) = reimin
endif
! -----------------------------------------------------------------------
! rain (lin et al., 1983)
! -----------------------------------------------------------------------
if (qr (i, j) .gt. qmin) then
qcr (i, j) = den (i, j) * qr (i, j)
lambdar = exp (0.25 * log (pi * rhor * n0r / qcr (i, j)))
rer (i, j) = 0.5 * exp (log (gammar / 6) / alphar) / lambdar * 1.0e6
rer (i, j) = max (rermin, min (rermax, rer (i, j)))
else
qcr (i, j) = 0.0
rer (i, j) = rermin
endif
! -----------------------------------------------------------------------
! snow (lin et al., 1983)
! -----------------------------------------------------------------------
if (qs (i, j) .gt. qmin) then
qcs (i, j) = den (i, j) * qs (i, j)
lambdas = exp (0.25 * log (pi * rhos * n0s / qcs (i, j)))
res (i, j) = 0.5 * exp (log (gammas / 6) / alphas) / lambdas * 1.0e6
res (i, j) = max (resmin, min (resmax, res (i, j)))
else
qcs (i, j) = 0.0
res (i, j) = resmin
endif
! -----------------------------------------------------------------------
! graupel (lin et al., 1983)
! -----------------------------------------------------------------------
if (qg (i, j) .gt. qmin) then
qcg (i, j) = den (i, j) * qg (i, j)
lambdag = exp (0.25 * log (pi * rhog * n0g / qcg (i, j)))
reg (i, j) = 0.5 * exp (log (gammag / 6) / alphag) / lambdag * 1.0e6
reg (i, j) = max (regmin, min (regmax, reg (i, j)))
else
qcg (i, j) = 0.0
reg (i, j) = regmin
endif
enddo
enddo
end subroutine cloud_diagnosis
end module gfdl_cloud_microphys_mod