| blob | 0987df60de31fa2190a3659841ea2589db02200a |
1 ! Note: Comments starting with "!!" are directly taken from CAM's driver routine for the micro-physics (stratiform.F90)
2 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
4 !----------------------------------------------------------------------------------------------------------
5 !Possible future Modifications:
6 !1. *** IMPORTANT-WARNING *** conv_water_in_rad is hardwired to be equal to 1. It should be a namelist variable
7 !2. QME3D computation is WRONG as it doesn't include CMELIQ contribution as CMELIQ is
8 ! computed in macrophysics, which is not ported into WRF yet
9 !3. Computations such as cldfsnow,icswp etc. are not correct as CONCLD (convective
10 ! cloud fraction) variable is not avaible to microphysics
11 !4. The code ouputs variables which can be used in RRTMG radiation but RRTMG radiation doesn't
12 ! use these variables as of now
13 !5. More robust cloud fraction calculation can be included in the future versions
14 !----------------------------------------------------------------------------------------------------------
16 !----------------------------------------------------------------------------------------------------------
17 !Important Notes:
18 !1. New cloud fraction computations are added by Po-Lun Ma (PMA)
19 !2. Macrophysics processes are added by PMA into this driver
20 !----------------------------------------------------------------------------------------------------------
22 #define MODAL_AERO
23 module module_mp_cammgmp_driver
25 !!-------------------------------------------------------------------------------------------------------
26 !! Purpose:
27 !!
28 !! Provides the CAM interface to the prognostic cloud microphysics
29 !!
30 !! Author: Andrew Gettelman, Cheryl Craig, October 2010
31 !! Origin: modified from stratiform.F90
32 !! (Boville 2002, Coleman 2004, Park 2009, Kay 2010)
33 ! Ported to WRF by Balwinder.Singh@pnnl.gov
34 !!-------------------------------------------------------------------------------------------------------
36 use shr_kind_mod, only: r8=>shr_kind_r8
37 use physconst, only: gravit
38 use module_cam_support, only: pcnst =>pcnst_mp, pcols, pver, pverp
39 #if ( WRF_CHEM == 1 )
40 use module_cam_support, only: cam_mam_aerosols
41 #endif
42 use constituents, only: cnst_get_ind, cnst_name, qmin
44 #if ( WRF_CHEM == 1 )
45 use module_state_description, only: num_chem, param_first_scalar, &
46 CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM7_AQ, CBMZ_CAM_MAM3_NOAQ, &
47 CBMZ_CAM_MAM3_AQ
48 #endif
51 implicit none
52 private
53 save
55 public :: CAMMGMP_INIT
56 public :: CAMMGMP
57 #if ( WRF_CHEM == 1 )
58 public :: physics_ptend_init ! Mimics CAM's physics ptend init. Also used in wet scavaging code in WRF_CHEM
59 public :: physics_update ! Mimics CAM's physics update. Also used in wet scavaging code in WRF_CHEM
60 #endif
62 !! ------------------------- !
63 !! Private Module Parameters !
64 !! ------------------------- !
66 !! Choose either 'intermediate' ('inter') or complete ('compl') cloud microphysics
67 !! inter : Microphysics assumes 'liquid stratus frac = ice stratus frac = max( liquid stratus frac, ice stratus frac )'.
68 !! compl : Microphysics explicitly treats 'liquid stratus frac .ne. ice stratus frac'
69 !! for CAM5, only 'inter' is functional
71 character(len=5), private, parameter :: micro_treatment = 'inter'
73 logical, private :: sub_column = .false. ! True = configure microphysics for sub-columns
74 !! False = use in regular mode w/o sub-columns
76 !! -------------------------------- !
77 !! End of Private Module Parameters !
78 !! -------------------------------- !
80 integer :: &
81 ixcldliq, &! cloud liquid amount index
82 ixcldice, &! cloud ice amount index
83 ixnumliq, &! cloud liquid number index
84 ixnumice ! cloud ice water index
85 integer :: chem_opt, ptend_top_level, ptend_bot_level
87 real(r8), parameter :: qthresh_mass = 1.0e-11 ! Typical val: 1e-6 kg/kg ; Negligible val: 1e-11 kg/kg
88 real(r8), parameter :: qthresh_numl = 5.0e2 ! Typical val: 5E7 #/kg ; Negligible val: 5e2 #/kg
89 real(r8), parameter :: qthresh_numi = 1.0e-1 ! Typical val: 1E4 #/kg ; Negligible val: 1e-1 #/kg
91 !Following constants (dp1 and dp2) are obtained from cloud_fraction.F90 of CAM for highest resolution(0.23x0.31)
92 real(r8), parameter :: dp1 = 0.10_r8
93 real(r8), parameter :: dp2 = 500.0_r8
95 contains
97 subroutine CAMMGMP(itimestep, dt, p8w, p_hyd &
98 , t_phy, pi_phy, z_at_w, qfx &
99 , tke_pbl, turbtype3d, smaw3d &
100 , dlf3d, dlf2_3d, rliq2d, z_sea_level &
101 , kvh3d, ht, alt, accum_mode &
102 , aitken_mode, coarse_mode, icwmrsh3d &
103 , icwmrdp3d, shfrc3d, cmfmc3d, cmfmc2_3d &
104 , config_flags, f_ice_phy, f_rain_phy &
105 #if ( WRF_CHEM == 1 )
106 , dgnum4d, dgnumwet4d &
107 #endif
108 , ids, ide, jds, jde, kds, kde &
109 , ims, ime, jms, jme, kms, kme &
110 , its, ite, jts, jte, kts, kte &
111 !Output variables from CAMMGMP
112 , th, cldfra_old_mp, cldfra_mp &
113 ,cldfra_mp_all, lradius, iradius, cldfrai &
114 , cldfral, cldfra_conv,wsedl3d, rainnc &
115 , rainncv, snownc, snowncv,sr &
116 , qv_curr, qc_curr, qi_curr,qs_curr &
117 , qr_curr, nc3d, ni3d,ns3d,nr3d,qndrop &
118 , rh_old_mp,lcd_old_mp & !PMA- added for macrophysics
119 #if ( WRF_CHEM == 1 )
120 , chem &
121 , qme3d,prain3d,nevapr3d &
122 , rate1ord_cw2pr_st3d &
123 #endif
124 , xland,snowh )
126 !!-------------------------------------------------------- !
127 !! !
128 !! Purpose: !
129 !! !
130 !! Interface to sedimentation, detrain, cloud fraction and !
131 !! cloud macro - microphysics subroutines !
132 !! !
133 !! Author: D.B. Coleman !
134 !! Date: Apr 2004 !
135 !! !
136 ! Polun Ma added simple macrophysics scheme
137 !!-------------------------------------------------------- !
138 use shr_kind_mod, only: r8 => shr_kind_r8
139 use cldwat2m_micro, only: mmicro_pcond
140 use microp_aero, only: microp_aero_ts
141 use physconst, only: cpair, tmelt ! tmelt is added for computing saturation specific humidity for ice and liquid
143 #ifdef MODAL_AERO
144 use modal_aero_data, only: modeptr_accum => modeptr_accum_mp, &
145 numptr_amode => numptr_amode_mp , lspectype_amode => lspectype_amode_mp, &
146 specdens_amode => specdens_amode_mp, voltonumb_amode => voltonumb_amode_mp, &
147 lmassptr_amode => lmassptr_amode_mp, modeptr_aitken => modeptr_aitken_mp, &
148 modeptr_coarse => modeptr_coarse_mp, ntot_amode => ntot_amode_mp, &
149 nspec_amode => nspec_amode_mp
150 #endif
152 #if ( WRF_CHEM == 1 )
153 use module_data_cam_mam_asect, only: lptr_chem_to_q, lptr_chem_to_qqcw, factconv_chem_to_q
154 #endif
156 use wv_saturation, only: epsqs, polysvp ! Added for computing saturation specific humidity for ice and liquid
157 use conv_water, only: conv_water_4rad
158 use cldwat, only: cldwat_fice
159 use module_state_description, only: CAMZMSCHEME, CAMUWSHCUSCHEME,F_QV, F_QC, F_QI, F_QS
160 use module_configure, only: grid_config_rec_type
162 implicit none
163 !
164 ! Input arguments
165 !
166 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
167 integer, intent(in) :: itimestep !Timestep count
169 integer, intent(in) :: ids,ide, jds,jde, kds,kde
170 integer, intent(in) :: ims,ime, jms,jme, kms,kme
171 integer, intent(in) :: its,ite, jts,jte, kts,kte
173 real, intent(in) :: dt !Timestep
175 real, intent(in) :: accum_mode, aitken_mode, coarse_mode
177 real, dimension( ims:ime, jms:jme ), intent(in) :: ht !Terrain height (m)
178 real, dimension( ims:ime, jms:jme ), intent(in) :: qfx !Upward moisture flux at the surface(kg m-2 s-1)
179 real, dimension( ims:ime, jms:jme ), intent(in) :: rliq2d !Vertically-integrated reserved cloud condensate(m s-1)
180 real, dimension( ims:ime, jms:jme ), intent(in) :: xland
181 real, dimension( ims:ime, jms:jme ), intent(in) :: snowh
183 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: dlf3d !Detraining cloud water tendency from convection
184 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: dlf2_3d !dq/dt due to export of cloud water into environment by shallow convection(kg/kg/s)
185 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p8w !Hydrostatic Pressure at level interface (Pa)
186 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: p_hyd !Hydrostatic pressure(Pa)
187 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_sea_level !Height above sea level at mid-level (m)
188 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: t_phy !Temperature (K)
189 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: z_at_w !Height above sea level at layer interfaces (m)
190 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: pi_phy !Dimensionless pressure [unitless]
191 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: tke_pbl !Turbulence kinetic energy from Mellor-Yamada-Janjic (MYJ) (m^2/s^2)
192 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: kvh3d !Eddy diffusivity for heat(m^2/s)
193 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: turbtype3d !Turbulent interface types [ no unit ]
194 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: smaw3d !Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
195 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: alt !inverse density(m3/kg)
196 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: F_ICE_PHY !Fraction of ice
197 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: F_RAIN_PHY !Fraction of rain
198 !For conv_water_4rad
199 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: shfrc3d !Shallow cloud fraction
200 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cmfmc3d !Deep + Shallow Convective mass flux [ kg /s/m^2 ]
201 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: cmfmc2_3d !Shallow convective mass flux [ kg/s/m^2 ]
202 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: icwmrsh3d !Shallow cumulus in-cloud water mixing ratio (kg/m2)
203 real, dimension( ims:ime, kms:kme, jms:jme ), intent(in) :: icwmrdp3d !Deep Convection in-cloud water mixing ratio (kg/m2)
205 #if ( WRF_CHEM == 1 )
206 real, dimension( ims:ime, kms:kme, jms:jme, ntot_amode ), intent(in) :: dgnum4d, dgnumwet4d ! 4-dimensional Number mode diameters
207 #endif
208 !2d in-outs
209 real, dimension( ims:ime , jms:jme ), intent(inout) :: rainnc !grid scale precipitation (mm)
210 real, dimension( ims:ime , jms:jme ), intent(inout) :: rainncv !one time step grid scale precipitation (mm/step)
211 real, dimension( ims:ime , jms:jme ), intent(inout) :: snownc !grid scale snow and ice (mm)
212 real, dimension( ims:ime , jms:jme ), intent(inout) :: snowncv !one time step grid scale snow and ice (mm/step)
214 !3d in-outs
215 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfra_old_mp !Old Cloud fraction
216 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rh_old_mp !Old RH
217 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: lcd_old_mp !Old liquid cloud fraction
218 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qv_curr !Water vapor mixing ratio - Moisture (kg/kg)
219 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qc_curr !Cloud water mixing ratio - Cloud liq (kg/kg)
220 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qi_curr !Ice mixing ratio -Cloud ice (kg/kg)
221 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qs_curr !Snow mixing ratio -Cloud ice (kg/kg)
222 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qr_curr !Rain mixing ratio -Cloud ice (kg/kg)
223 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: ni3d !Cloud ice number concentration (#/kg)
224 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: nc3d !Cloud water number concentration (#/kg)
225 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: ns3d !Snow number concentration (#/kg)
226 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: nr3d !rain number concentration (#/kg)
227 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qndrop !droplet number mixing ratio (#/kg)
228 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: th !Potential temperature (K)
229 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: wsedl3d !Sedimentation velocity of stratiform liquid cloud droplet (m/s)
230 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfra_mp !New Cloud fraction
231 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfra_conv !New Cloud fraction
232 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfra_mp_all !New Cloud fraction
233 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: lradius
234 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: iradius
235 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfrai
236 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: cldfral
237 #if ( WRF_CHEM == 1 )
238 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: qme3d !Net condensation rate (kg/kg/s)
239 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: prain3d !Rate of conversion of condensate to precipitation (kg/kg/s)
240 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: nevapr3d !Evaporation rate of rain + snow (kg/kg/s)
241 real, dimension( ims:ime, kms:kme, jms:jme ), intent(inout) :: rate1ord_cw2pr_st3d !1st order rate for direct conversion of strat. cloud water to precip (1/s)
242 #endif
244 #if ( WRF_CHEM == 1 )
245 !4d in-outs
246 real, dimension( ims:ime, kms:kme, jms:jme, num_chem ), intent(inout) :: chem !Chem array
247 #endif
248 !2d outs
249 real, dimension( ims:ime , jms:jme ), intent( out) :: sr !one time step mass ratio of snow to total precip
252 !
253 ! Local variables
254 !
256 !Other variables
257 character*250 :: wrfmessage
258 logical :: is_first_step !Flag to know if this a first time step
259 logical :: ptend_loc_ls, ptend_all_ls !Flag for updating tendencies
260 logical :: ptend_loc_lq(pcnst),ptend_all_lq(pcnst) !Flag for updating tendencies
261 integer :: i,k,m,n,iw,jw,kw,imode,itsm1,itile_len,ktep1,kflip,ips,kcam
262 #if ( WRF_CHEM == 1 )
263 integer :: p1st,l2 !For iterating loop of NUM_CHEM for CHEM array
264 #endif
265 integer :: lchnk !Chunk identifier
266 integer :: ncol !Number of atmospheric columns
267 integer :: conv_water_in_rad
268 real (r8) :: multFrc
269 real(r8) :: dtime !Timestep in real(8) format
270 real(r8) :: dp
272 real(r8) :: phis(pcols) !Geopotential at terrain height (m2/s2)
273 real(r8) :: state_t(pcols,kte)
274 real(r8) :: state_q(pcols,kte,pcnst)
275 real(r8) :: state_s(pcols,kte)
276 real(r8) :: state_pmid(pcols,kte)
277 real(r8) :: state_pdel(pcols,kte)
278 real(r8) :: state_rpdel(pcols,kte)
279 real(r8) :: state_zm(pcols,kte)
280 real(r8) :: state_omega(pcols,kte)
282 real(r8) :: state_pint(pcols,kte+1)
284 real(r8) :: state1_t(pcols,kte)
285 real(r8) :: state1_q(pcols,kte,pcnst)
286 real(r8) :: state1_s(pcols,kte)
287 real(r8) :: state1_pmid(pcols,kte)
288 real(r8) :: state1_pdel(pcols,kte)
289 real(r8) :: state1_rpdel(pcols,kte)
290 real(r8) :: state1_zm(pcols,kte)
291 real(r8) :: state1_omega(pcols,kte)
293 real(r8) :: state1_pint(pcols,kte+1)
295 character*24 :: ptend_loc_name
296 real(r8) :: ptend_loc_s(pcols,kte)
297 real(r8) :: ptend_loc_q(pcols,kte,pcnst)
299 character*24 :: ptend_all_name
300 real(r8) :: ptend_all_s(pcols,kte)
301 real(r8) :: ptend_all_q(pcols,kte,pcnst)
303 real(r8) :: tmpnaer,tmpmaer
304 real(r8) :: cmeliq(pcols,kte) ! Rate of cond-evap of liq within the cloud
305 real(r8) :: nrout(pcols,kte) ! rain number concentration (1/m3)
306 real(r8) :: nsout(pcols,kte) ! snow number concentration (1/m3)
307 real(r8) :: sh_icwmr(pcols,kte) ! Shallow conv. cloud water
308 real(r8) :: dp_icwmr(pcols,kte) ! Deep conv. cloud water
309 real(r8) :: fice(pcols,kte) ! Ice partitioning ratio
310 real(r8) :: fsnow(pcols,kte) ! Fractional snow content for convection
311 real(r8) :: sh_frac(pcols,kte) ! Shallow convective cloud fraction
312 real(r8) :: dp_frac(pcols,kte) ! Deep convective cloud fraction
313 real(r8) :: cmfmc(pcols,kte+1)
314 real(r8) :: cmfmc2(pcols,kte+1)
316 real(r8), dimension(pcols,kte) :: esl ! liquid sat vapor pressure (pa) ! for phil suggested modifications
317 real(r8), dimension(pcols,kte) :: esi ! ice sat vapor pressure (pa)
318 real(r8), dimension(pcols,kte) :: qvs ! liquid saturation vapor mixing ratio ! for phil suggested modifications
319 real(r8), dimension(pcols,kte) :: qvi ! ice saturation vapor mixing ratio
322 !Following variables are added by PMA for Macrophysics
323 real(r8), dimension(pcols,kte) :: ttt ! Temperature (PMA)
324 real(r8), dimension(pcols,kte) :: rh_old ! RH at previous timestep (PMA)
325 real(r8), dimension(pcols,kte) :: qi_mac ! Temperature (PMA)
328 !End variables addition by PMA
331 !Variables with CAM data structure
332 real(r8) :: rliq(pcols) !Vertical integral of liquid not yet in q(ixcldliq) or vertically-integrated reserved cloud condensate (m/s)
334 !Following 6 variables are intent-out in CAM's driver for micro-physics
335 real(r8) :: prec_str(pcols) ! [Total] Sfc flux of precip from stratiform [ m/s ]
336 real(r8) :: snow_str(pcols) ! [Total] Sfc flux of snow from stratiform [ m/s ]
337 real(r8) :: prec_sed(pcols) ! Surface flux of total cloud water from sedimentation
338 real(r8) :: snow_sed(pcols) ! Surface flux of cloud ice from sedimentation
339 real(r8) :: prec_pcw(pcols) ! Sfc flux of precip from microphysics [ m/s ]
340 real(r8) :: snow_pcw(pcols) ! Sfc flux of snow from microphysics [ m/s ]
342 real(r8) :: cflx(pcols,pcnst) ! Constituent flux from surface
343 real(r8) :: dlf(pcols,kte) ! Detrained water from convection schemes
344 real(r8) :: dlf2(pcols,kte) ! Detrained water from shallow convection scheme
345 real(r8) :: rate1cld(pcols,kte) ! array to hold rate1ord_cw2pr_st from microphysics
347 ! Physics buffer fields
349 real(r8), dimension(pcols,kte) :: cld ! Total cloud fraction
350 real(r8), dimension(pcols,kte) :: ast ! Relative humidity cloud fraction
351 real(r8), dimension(pcols,kte) :: aist ! Physical ice stratus fraction
352 real(r8), dimension(pcols,kte) :: alst ! Physical liquid stratus fraction
353 real(r8), dimension(pcols,kte) :: concld ! Convective cloud fraction
354 real(r8), dimension(pcols,kte) :: qme
355 real(r8), dimension(pcols,kte) :: prain ! Total precipitation (rain + snow)
356 real(r8), dimension(pcols,kte) :: nevapr ! Evaporation of total precipitation (rain + snow)
357 real(r8), dimension(pcols,kte) :: rel ! Liquid effective drop radius (microns)
358 real(r8), dimension(pcols,kte) :: rei ! Ice effective drop size (microns)
359 real(r8), dimension(pcols,kte) :: rel2 ! Liquid effective drop radius (microns)
360 real(r8), dimension(pcols,kte) :: rei2 ! Ice effective drop size (microns)
361 real(r8), dimension(pcols,kte) :: cldo ! Old cloud fraction
362 real(r8), dimension(pcols,kte) :: wsedl ! Sedimentation velocity of liquid stratus cloud droplet [ m/s ]
363 real(r8), dimension(pcols,kte) :: rel_fn ! Ice effective drop size at fixed number (indirect effect) (microns)
364 real(r8), dimension(pcols,kte+1) :: kkvh ! Vertical eddy diffusivity
365 #ifdef MODAL_AERO
366 real(r8), target, dimension(pcols,kte,pcnst) :: qqcw ! Cloud-borne aerosol
367 real(r8), dimension(pcols,kte,ntot_amode) :: dgnumwet ! Number mode diameter
368 real(r8), dimension(pcols,kte,ntot_amode) :: dgnum ! Number mode diameter
369 real(r8), dimension(pcols,kte) :: rate1ord_cw2pr_st ! 1st order rate for direct conversion of
370 ! strat. cloud water to precip (1/s) ! rce 2010/05/01
371 #endif
374 ! physics buffer fields for COSP simulator
375 real(r8), dimension(pcols,kte+1) :: mgflxprc ! MG grid-box mean flux_large_scale_cloud_rain+snow at interfaces (kg/m2/s)
376 real(r8), dimension(pcols,kte+1) :: mgflxsnw ! MG grid-box mean flux_large_scale_cloud_snow at interfaces (kg/m2/s)
377 real(r8), dimension(pcols,kte) :: mgmrprc ! MG grid-box mean mixingratio_large_scale_cloud_rain+snow at interfaces (kg/kg)
378 real(r8), dimension(pcols,kte) :: mgmrsnw ! MG grid-box mean mixingratio_large_scale_cloud_snow at interfaces (kg/kg)
379 real(r8), dimension(pcols,kte) :: mgreffrain ! MG diagnostic rain effective radius (um)
380 real(r8), dimension(pcols,kte) :: mgreffsnow ! MG diagnostic snow effective radius (um)
381 real(r8), dimension(pcols,kte) :: cvreffliq ! convective cloud liquid effective radius (um)
382 real(r8), dimension(pcols,kte) :: cvreffice ! convective cloud ice effective radius (um)
384 ! physics buffer fields for radiation
386 real(r8), dimension(pcols,kte) :: dei ! Ice effective diameter (meters) (AG: microns?)
387 real(r8), dimension(pcols,kte) :: mu ! Size distribution shape parameter for radiation
388 real(r8), dimension(pcols,kte) :: lambdac ! Size distribution slope parameter for radiation
389 real(r8), dimension(pcols,kte) :: iciwp ! In-cloud ice water path for radiation
390 real(r8), dimension(pcols,kte) :: iclwp ! In-cloud liquid water path for radiation
392 ! For rrtm optics. specificed distribution.
394 real(r8) :: mucon ! Convective size distribution shape parameter
395 real(r8) :: dcon ! Convective size distribution effective radius (meters)
396 real(r8) :: lamcon ! Convective size distribution slope parameter (meters-1)
397 real(r8) :: deicon ! Convective ice effective diameter (meters)
399 ! Physics buffer fields
401 real(r8), dimension(pcols,kte) :: deiconv ! Ice effective diameter (microns)
402 real(r8), dimension(pcols,kte) :: muconv ! Size distribution shape parameter for radiation
403 real(r8), dimension(pcols,kte) :: lambdaconv ! Size distribution slope parameter for radiation
404 real(r8), dimension(pcols,kte) :: iciwpst ! Stratiform in-cloud ice water path for radiation
405 real(r8), dimension(pcols,kte) :: iclwpst ! Stratiform in-cloud liquid water path for radiation
406 real(r8), dimension(pcols,kte) :: iciwpconv ! Convective in-cloud ice water path for radiation
407 real(r8), dimension(pcols,kte) :: iclwpconv ! Convective in-cloud liquid water path for radiation
409 real(r8), dimension(pcols,kte+1) :: tke ! TKE from the moist PBL scheme
410 real(r8), dimension(pcols,kte+1) :: turbtype ! Turbulence type from the moist PBL scheme
411 real(r8), dimension(pcols,kte+1) :: smaw ! Instability function of momentum from the moist PBL scheme
413 ! Local variables for in-cloud water quantities adjusted for convective water
415 real(r8) allcld_ice (pcols,kte) ! All-cloud cloud ice
416 real(r8) allcld_liq (pcols,kte) ! All-cloud liquid
418 ! Snow
420 real(r8), dimension(pcols,kte) :: cldfsnow ! Cloud fraction for liquid+snow
421 real(r8), dimension(pcols,kte) :: icswp ! In-cloud snow water path
422 real(r8), dimension(pcols,kte) :: des ! Snow effective diameter (m)
423 real(r8) qsout(pcols,kte) ! Snow mixing ratio
425 real(r8) :: qrout(pcols,kte) ! Rain mixing ratio
426 real(r8) :: rflx(pcols,kte+1) ! grid-box average rain flux (kg m^-2 s^-1)
427 real(r8) :: sflx(pcols,kte+1) ! grid-box average snow flux (kg m^-2 s^-1)
428 real(r8) :: reff_rain(pcols,kte) ! rain effective radius (um)
429 real(r8) :: reff_snow(pcols,kte) ! snow effective radius (um)
431 real(r8) icecldf(pcols,kte) ! Ice cloud fraction
432 real(r8) liqcldf(pcols,kte) ! Liquid cloud fraction (combined into cloud)
433 real(r8) icecldf_out(pcols,kte) ! Ice cloud fraction
434 real(r8) liqcldf_out(pcols,kte) ! Liquid cloud fraction (combined into cloud)
436 ! Local variables for microphysics
438 real(r8) rdtime ! 1./dtime
439 real(r8) qtend(pcols,kte) ! Moisture tendencies
440 real(r8) ttend(pcols,kte) ! Temperature tendencies
441 real(r8) ltend(pcols,kte) ! Cloud liquid water tendencies
442 real(r8) evapsnow(pcols,kte) ! Local evaporation of snow
443 real(r8) prodsnow(pcols,kte) ! Local production of snow
444 real(r8) icimr(pcols,kte) ! In cloud ice mixing ratio
445 real(r8) icwmr(pcols,kte) ! In cloud water mixing ratio
446 real(r8) icimrst(pcols,kte) ! In stratus ice mixing ratio
447 real(r8) icwmrst(pcols,kte) ! In stratus water mixing ratio
448 real(r8) icimrst_out(pcols,kte) ! In stratus ice mixing ratio
449 real(r8) icwmrst_out(pcols,kte) ! In stratus water mixing ratio
450 real(r8) cmeice(pcols,kte) ! Rate of cond-evap of ice within the cloud
451 real(r8) temp(pcols)
452 real(r8) res(pcols,kte)
454 ! MG micro diagnostics
456 real(r8) qcsevap(pcols,kte) ! Evaporation of falling cloud water
457 real(r8) qisevap(pcols,kte) ! Sublimation of falling cloud ice
458 real(r8) qvres(pcols,kte) ! Residual condensation term to remove excess saturation
459 real(r8) cmeiout(pcols,kte) ! Deposition/sublimation rate of cloud ice
460 real(r8) vtrmc(pcols,kte) ! Mass-weighted cloud water fallspeed
461 real(r8) vtrmi(pcols,kte) ! Mass-weighted cloud ice fallspeed
462 real(r8) qcsedten(pcols,kte) ! Cloud water mixing ratio tendency from sedimentation
463 real(r8) qisedten(pcols,kte) ! Cloud ice mixing ratio tendency from sedimentation
465 real(r8) prao(pcols,kte)
466 real(r8) prco(pcols,kte)
467 real(r8) mnuccco(pcols,kte)
468 real(r8) mnuccto(pcols,kte)
469 real(r8) mnuccdo(pcols,kte)
470 real(r8) mnuccdohet(pcols,kte)
471 real(r8) msacwio(pcols,kte)
472 real(r8) psacwso(pcols,kte)
473 real(r8) bergso(pcols,kte)
474 real(r8) bergo(pcols,kte)
475 real(r8) melto(pcols,kte)
476 real(r8) homoo(pcols,kte)
477 real(r8) qcreso(pcols,kte)
478 real(r8) prcio(pcols,kte)
479 real(r8) praio(pcols,kte)
480 real(r8) qireso(pcols,kte)
481 real(r8) ftem(pcols,kte)
482 real(r8) mnuccro(pcols,kte)
483 real(r8) pracso (pcols,kte)
484 real(r8) meltsdt(pcols,kte)
485 real(r8) frzrdt (pcols,kte)
486 real(r8) dpdlfliq(pcols,kte)
487 real(r8) dpdlfice(pcols,kte)
488 real(r8) shdlfliq(pcols,kte)
489 real(r8) shdlfice(pcols,kte)
490 real(r8) dpdlft (pcols,kte)
491 real(r8) shdlft (pcols,kte)
493 #ifdef MODAL_AERO
494 integer l, lnum, lnumcw, lmass, lmasscw
495 #endif
497 ! Variables for MG microphysics
499 real(r8) dum1,dum2
500 real(r8) qc(pcols,kte)
501 real(r8) qi(pcols,kte)
502 real(r8) nc(pcols,kte)
503 real(r8) ni(pcols,kte)
504 real(r8) icinc(pcols,kte) ! In cloud ice number conc
505 real(r8) cdnumc(pcols) ! Vertically-integrated droplet concentration
506 real(r8) icwnc(pcols,kte) ! In cloud water number conc
507 real(r8) iwc(pcols,kte) ! Grid box average ice water content
508 real(r8) lwc(pcols,kte) ! Grid box average liquid water content
509 real(r8) effliq(pcols,kte) ! In cloud liq eff rad
510 real(r8) effice(pcols,kte) ! In cloud ice eff rad
511 real(r8) effliq_fn(pcols,kte) ! In cloud liq eff rad at fixed number concentration
512 real(r8) wsub(pcols,kte) ! Sub-grid vertical velocity (m/s)
513 real(r8) wsubi(pcols,kte) ! Sub-grid vertical velocity for ice (m/s)
515 ! Output from mmicro_pcond
517 real(r8) tlat(pcols,kte)
518 real(r8) qvlat(pcols,kte)
519 real(r8) qcten(pcols,kte)
520 real(r8) qiten(pcols,kte)
521 real(r8) ncten(pcols,kte)
522 real(r8) niten(pcols,kte)
523 real(r8) effc(pcols,kte)
524 real(r8) effc_fn(pcols,kte) ! Liquid effective radius at fixed number (for indirect calc)
525 real(r8) effi(pcols,kte)
526 real(r8) prect(pcols)
527 real(r8) preci(pcols)
529 ! Output from microp_aero_ts for aerosol actication
530 real(r8) naai(pcols,kte) !ice nucleation number
531 real(r8) naai_hom(pcols,kte) !ice nucleation number (homogeneous)
532 real(r8) npccn(pcols,kte) !liquid activation number tendency
533 real(r8) rndst(pcols,kte,4)
534 real(r8) nacon(pcols,kte,4)
536 ! Averaging arrays for effective radius and number....
538 real(r8) efiout(pcols,kte)
539 real(r8) efcout(pcols,kte)
540 real(r8) ncout(pcols,kte)
541 real(r8) niout(pcols,kte)
543 real(r8) freqi(pcols,kte)
544 real(r8) freql(pcols,kte)
546 ! Average cloud top radius & number
548 real(r8) ctrel(pcols)
549 real(r8) ctrei(pcols)
550 real(r8) ctnl(pcols)
551 real(r8) ctni(pcols)
552 real(r8) fcti(pcols)
553 real(r8) fctl(pcols)
555 ! Gather mass mixing ratio for all aerosols affecting the climate
557 integer :: naer_all
558 real(r8) :: aermmr1(pcols,kte)
559 real(r8), allocatable :: aer_mmr(:,:,:) ! Aerosol mass mixing ratio
561 real(r8) zeros(pcols,kte)
563 real(r8) alst_mic(pcols,kte)
564 real(r8) aist_mic(pcols,kte)
566 real(r8), pointer :: fldcw(:,:)
567 real(r8) :: xland_pt
568 real(r8) :: snowh_pt
569 ! ======================================================================
572 !cmeliq is an output from Macro physics. Its assigned zero here as Macrophysics is NOT implemented yet
573 !*** IMPORTANT-WARNING ***: The following assignment will produce an incorrect value of QME3D
574 cmeliq(:,:)=0.0_r8
576 !*** IMPORTANT-WARNING *** :Should be in the namelist file. Hardwired currently
577 conv_water_in_rad = 1
579 #if ( WRF_CHEM == 1 )
580 p1st = param_first_scalar ! Obtain CHEM array's first element's index
581 #endif
583 ! Specify diameter for each mode
584 #if ( WRF_CHEM != 1 )
585 dgnumwet(:,:,:) = 0.1e-6_r8 !*** IMPORTANT-WARNING *** :Constant value for the whole domain (prescribed value). The value match precribe_aerosol_mixactivate, but units change
586 dgnumwet(:,:,modeptr_aitken) = 0.02e-6
587 dgnumwet(:,:,modeptr_coarse) = 1.0e-6
588 #else
589 if(chem_opt == 0 ) then
590 dgnumwet(:,:,:) = 0.1e-6_r8 !*** IMPORTANT-WARNING *** :Constant value for the whole domain (prescribed value). The value match precribe_aerosol_mixactivate, but units change
591 dgnumwet(:,:,modeptr_aitken) = 0.02e-6
592 dgnumwet(:,:,modeptr_coarse) = 1.0e-6
593 endif
594 #endif
596 !dgnum wet and dry are assumed to be same for now
597 dgnum(:,:,:) = dgnumwet(:,:,:)
599 !Time step is stored in the (r8) format in dtime
600 dtime = dt
603 !default values for radius
604 lradius(:,:,:) = 10._r8
605 iradius(:,:,:) = 70._r8
607 !Flag for first time step
608 is_first_step = .false.
609 if(itimestep == 1) then
610 is_first_step = .true.
611 cldfra_old_mp(:,:,:) = 0.0_r8
612 rate1ord_cw2pr_st(:,:) = 0.0_r8
613 cldfrai(:,:,:) = 0._r8
614 cldfral(:,:,:) = 0._r8
615 cldfra_mp(:,:,:) = 0._r8
616 cldfra_mp_all(:,:,:) = 0._r8
617 cldfra_conv(:,:,:) = 0._r8
619 if(config_flags%shcu_physics .NE. CAMUWSHCUSCHEME) call wrf_message('WARNING: sh_icwmr,cmfmc,cmfmc2 and sh_frac are set to zero in CAMMGMP')
620 if(config_flags%cu_physics .NE. CAMZMSCHEME) call wrf_message('WARNING: dp_icwmr is set to zero in CAMMGMP')
621 endif
623 ncol = pcols
624 !This subroutine requires that ncol == 1
625 if(ncol .NE. 1) then
626 call wrf_error_fatal('Number of CAM Columns (NCOL) in CAMMGMP scheme must be 1')
627 endif
629 !*** IMPORTANT-WARNING ***:Initializing variables which are **NEVER** used but required for call consistency for CAM specific subroutine calls
630 concld(:,:) =0.0_r8 !Required for computing some diagnostic variables. These diagnostic variables are
631 !commented out for now as we didn't map 'concld' with the WRF sided variables
633 state_omega(:,:) = 0.0_r8 !Required for calling mmicro_pcond and microp_aero_ts. Used in microp_aero_ts for
634 !calling dropmixnuc but **NEVER** used for any output from dropmixnuc
636 !Divide domain in chuncks and map WRF variables into CAM
637 !Loop counters are named iw,jw,kw to represent that they index WRF sided variables
638 !The CAM's side variables are kept almost same by just replacing '%' with an '_' [e.g state%t in CAM is state_t in WRF]
640 itsm1 = its - 1
641 itile_len = ite - itsm1
642 do jw = jts , jte
643 do iw = its , ite
645 xland_pt = xland(iw,jw)
646 snowh_pt = snowh(iw,jw)
648 lchnk = (jw - jts) * itile_len + (iw - itsm1) !1-D index location from a 2-D tile
649 phis(1) = ht(iw,jw) * gravit
650 ktep1 = kte + 1
652 !Flip vertically quantities computed at the mid points
653 do kw = kts, kte
654 kflip = ktep1 - kw
655 state_pmid(1,kflip) = p_hyd(iw,kw,jw) !Pressure at the mid-points (Pa) [state%pmid in CAM]
656 dp = p8w(iw,kw,jw) - p8w(iw,kw+1,jw) !Change in pressure (Pa)
657 state_pdel(1,kflip) = dp
658 state_rpdel(1,kflip) = 1.0_r8/dp !Reciprocal of pressure difference (1/Pa) [state%rpdel in CAM]
659 state_zm(1,kflip) = z_sea_level(iw,kw,jw) - ht(iw,jw) !Height above the ground at the midpoints (m) [state%zm in CAM]
660 state_t(1,kflip) = t_phy(iw,kw,jw) !Temprature at the mid points (K) [state%t in CAM]
661 state_s(1,kflip) = cpair *th(iw,kw,jw)*pi_phy(iw,kw,jw) + gravit*state_zm(1,kflip) + phis(1) ! Dry static energy (m2/s2) [state%s in CAM]-Formula tested in vertical_diffusion.F90 of CAM
663 !Following three formulas are obtained from ported CAM's ZM cumulus scheme
664 !Values of 0 cause a crash in entropy
665 multFrc = 1._r8/(1._r8 + qv_curr(iw,kw,jw))
666 state_q(1,kflip,1) = qv_curr(iw,kw,jw)*multFrc !Specific humidity [state%q(:,:,1) in CAM]
667 state_q(1,kflip,ixcldliq) = qc_curr(iw,kw,jw)*multFrc !Convert to moist mix ratio-cloud liquid [state%q(:,:,2) in CAM]
668 state_q(1,kflip,ixcldice) = qi_curr(iw,kw,jw)*multFrc !cloud ice [state%q(:,:,3) in CAM]
669 state_q(1,kflip,ixnumliq) = nc3d(iw,kw,jw)*multFrc
670 state_q(1,kflip,ixnumice) = ni3d(iw,kw,jw)*multFrc
672 qi_mac(1,kflip) = qi_curr(iw,kw,jw)*multFrc
674 !Check for negative constituents, value below 'qthresh_...' are reset to qmin
675 !For mass mixing ratio
676 if(abs(state_q(1,kflip,1)) < qthresh_mass .and. state_q(1,kflip,1) < qmin(1)) state_q(1,kflip,1) = qmin(1)
677 if(abs(state_q(1,kflip,ixcldliq)) < qthresh_mass .and. state_q(1,kflip,ixcldliq) < qmin(ixcldliq)) state_q(1,kflip,ixcldliq) = qmin(ixcldliq)
678 if(abs(state_q(1,kflip,ixcldice)) < qthresh_mass .and. state_q(1,kflip,ixcldice) < qmin(ixcldice)) state_q(1,kflip,ixcldice) = qmin(ixcldice)
680 !For number mixing ratio
681 if(abs(state_q(1,kflip,ixnumliq)) < qthresh_numl .and. state_q(1,kflip,ixnumliq) < qmin(ixnumliq)) state_q(1,kflip,ixnumliq) = qmin(ixnumliq)
682 if(abs(state_q(1,kflip,ixnumice)) < qthresh_numi .and. state_q(1,kflip,ixnumice) < qmin(ixnumice)) state_q(1,kflip,ixnumice) = qmin(ixnumice)
685 !Set negative values to qmin and generate a warning message
686 if(state_q(1,kflip,1) < 0.0_r8)then
687 !write(wrfmessage,*)'Water vapor mass mixing ratio was:', state_q(1,kflip,1),'; Reset to',qmin(1), 'in CAMMGMP driver'
688 !call wrf_message(wrfmessage)
689 state_q(1,kflip,1) = qmin(1)
690 endif
692 if(state_q(1,kflip,ixcldliq) < 0.0_r8) then
693 !write(wrfmessage,*)'Liquid cloud mass mixing ratio was:', state_q(1,kflip,ixcldliq),'; Reset to',qmin(ixcldliq), 'in CAMMGMP driver'
694 !call wrf_message(wrfmessage)
695 state_q(1,kflip,ixcldliq) = qmin(ixcldliq)
696 endif
698 if(state_q(1,kflip,ixcldice) < 0.0_r8) then
699 !write(wrfmessage,*)'Ice cloud mass mixing ratio was:', state_q(1,kflip,ixcldice),'; Reset to',qmin(ixcldice), 'in CAMMGMP driver'
700 !call wrf_message(wrfmessage)
701 state_q(1,kflip,ixcldice) = qmin(ixcldice)
702 endif
704 if(state_q(1,kflip,ixnumliq) < 0.0_r8) then
705 !write(wrfmessage,*)'Liquid droplet number was:', state_q(1,kflip,ixnumliq),'; Reset to',qmin(ixnumliq), 'in CAMMGMP driver'
706 !call wrf_message(wrfmessage)
707 state_q(1,kflip,ixnumliq) = qmin(ixnumliq)
708 endif
710 if(state_q(1,kflip,ixnumice) < 0.0_r8) then
711 !write(wrfmessage,*)'Ice crystal number was:', state_q(1,kflip,ixnumice),'; Reset to',qmin(ixnumice), 'in CAMMGMP driver'
712 !call wrf_message(wrfmessage)
713 state_q(1,kflip,ixnumice) = qmin(ixnumice)
714 endif
715 ! For prescribed aerosols
716 qqcw(1,kflip,:) = 1.e-38_r8 !used in microp_aero for dropmicnuc, presently set to a constant value
717 #if ( WRF_CHEM == 1 )
718 if(chem_opt.NE.0 .and. config_flags%CAM_MP_MAM_cpled ) then
719 !Following Do-Loop is obtained from chem/module_cam_mam_aerchem_driver.F
720 do l = p1st, num_chem
721 l2 = lptr_chem_to_q(l)
722 if ((l2 >= 1) .and. (l2 <= pcnst)) then
723 state_q(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)
724 if(state_q(1,kflip,l2) < 0.0_r8) then
725 !write(wrfmessage,*)'Chem species with index', l,' was:', state_q(1,kflip,l2),'; Reset to 0.0 in CAMMGMP driver'
726 !call wrf_message(wrfmessage)
727 state_q(1,kflip,l2) = 0.0_r8
728 endif
729 end if
730 l2 = lptr_chem_to_qqcw(l)
731 if ((l2 >= 1) .and. (l2 <= pcnst)) then
732 qqcw(1,kflip,l2) = chem(iw,kw,jw,l)*factconv_chem_to_q(l)
733 if(qqcw(1,kflip,l2)< 0.0_r8) then
734 !write(wrfmessage,*)'Chem species with index', l,' was:', qqcw(1,kflip,l2),'; Reset to 0.0 in CAMMGMP driver'
735 !call wrf_message(wrfmessage)
736 qqcw(1,kflip,l2) = 0.0_r8
737 endif
738 end if
739 end do ! l
740 !Populate dgnums appropriately
741 do imode = 1 , ntot_amode
742 if(is_first_step) then
743 dgnum(1,kflip,imode) = 0.0_r8 !Set it to zero for the first time step
744 dgnumwet(1,kflip,imode) = 0.0_r8
745 else
746 dgnum(1,kflip,imode) = dgnum4D(iw,kw,jw,imode) !Obtained from 4D arrays
747 dgnumwet(1,kflip,imode) = dgnumwet4D(iw,kw,jw,imode)
748 endif
749 enddo
750 endif
751 #endif
753 cldo(1,kflip) = cldfra_mp_all(iw,kw,jw)
755 !New code by PMA
756 rh_old(1,kflip) = rh_old_mp(iw,kw,jw)
757 dlf(1,kflip) = 0.0_r8
758 dlf2(1,kflip) = 0.0_r8
760 sh_icwmr(1,kflip) = 0.0_r8
761 sh_frac(1,kflip) = 0.0_r8
762 dp_frac(1,kflip) = 0._r8
763 if(config_flags%shcu_physics==CAMUWSHCUSCHEME) then
764 sh_icwmr(1,kflip) = icwmrsh3d(iw,kw,jw) !Shallow cumulus in-cloud water mixing ratio (kg/m2)
765 sh_frac(1,kflip) = shfrc3d(iw,kw,jw) !Shallow convective cloud fraction
766 dlf2(1,kflip) = dlf2_3d(iw,kw,jw) !PMA QC detrained from shallow
767 endif
768 dp_icwmr(1,kflip) = 0.0_r8
769 if(config_flags%cu_physics==CAMZMSCHEME)then
770 dp_icwmr(1,kflip) = icwmrdp3d(iw,kw,jw) !Deep conv. cloud water
771 dlf(1,kflip) = dlf3d(iw,kw,jw) !PMA Qc detrained from deep+shallow
772 endif
774 enddo
776 do kw = kts, kte+1
777 kflip = kte - kw + 2
779 state_pint(1,kflip) = p8w(iw,kw,jw) !Pressure at interfaces [state%pint in CAM]
780 tke(1,kflip) = tke_pbl(iw,kw,jw) !Turbulent kinetic energy
781 kkvh(1,kflip) = kvh3d(iw,kw,jw) !Vertical diffusivity (m2/s)
782 turbtype(1,kflip) = turbtype3d(iw,kw,jw) !Turbulent interface types [ no unit ]
783 smaw(1,kflip) = smaw3d(iw,kw,jw) !Normalized Galperin instability function for momentum ( 0<= <=4.964 and 1 at neutral ) [no units]
785 cmfmc(1,kflip) = 0.0_r8
786 cmfmc2(1,kflip) = 0.0_r8
787 if(config_flags%shcu_physics==CAMUWSHCUSCHEME .or. config_flags%cu_physics == CAMZMSCHEME ) then
788 cmfmc(1,kflip) = cmfmc3d(iw,kw,jw) !Deep + Shallow Convective mass flux [ kg /s/m^2 ]
789 endif
791 if(config_flags%shcu_physics==CAMUWSHCUSCHEME ) then
792 cmfmc2(1,kflip) = cmfmc2_3d(iw,kw,jw) !Shallow convective mass flux [ kg/s/m^2 ]
793 endif
795 end do
796 do kcam = 1, kte
797 !Formulation for dp_frac is obtained from cloud_fraction.F90 of CAM
798 dp_frac(1,kcam) = max(0.0_r8,min(dp1*log(1.0_r8+dp2*(cmfmc(1,kcam+1)-cmfmc2(1,kcam+1))),0.60_r8)) ! Deep convective cloud fraction
799 end do
802 !set cflx zero for all the constituents, update the water vapor surface flux after that
803 cflx(:pcols,:) = 0.0_r8 !Surface constituent flux (kg/m2/s)
804 cflx(:pcols,1) = qfx(iw,jw) !Surface constituent flux (kg/m2/s)
805 rliq(:pcols) = rliq2d(iw,jw) !Vertically-integrated reserved cloud condensate (m/s)
808 !Following assignments mimics CAM's 'physics_state_copy' call
809 state1_pmid(:,:) = state_pmid(:,:)
810 state1_pdel(:,:) = state_pdel(:,:)
811 state1_rpdel(:,:) = state_rpdel(:,:)
812 state1_zm(:,:) = state_zm(:,:)
813 state1_t(:,:) = state_t(:,:)
814 state1_s(:,:) = state_s(:,:)
815 state1_pint(:,:) = state_pint(:,:)
816 state1_q(:,:,:) = state_q(:,:,:)
817 state1_omega(:,:) = state_omega(:,:) !*** IMPORTANT-WARNING ***:State_omega is *NOT* used for any meaningful computation at this time
819 call physics_ptend_init(ptend_loc_name,ptend_loc_q,ptend_loc_s,ptend_loc_lq,ptend_loc_ls,pcnst) !! Initialize local ptend type
820 call physics_ptend_init(ptend_all_name,ptend_all_q,ptend_all_s,ptend_all_lq,ptend_all_ls,pcnst) !! Initialize output ptend type
822 #ifdef MODAL_AERO
823 if( is_first_step) then
824 do i=1,pcnst
825 fldcw => qqcw(:,:,i)
826 if(associated(fldcw)) then
827 fldcw(1:pcols,1:pver) = 1.e-38_r8
828 end if
829 end do
830 endif
831 #endif
833 !! If first timestep, initialize heatflux....in pbuf at all time levels.
835 if( is_first_step) then
836 kkvh(:,:) = 0._r8
837 tke(:,:) = 0._r8
838 ! *** IMPORTANT-WARNING ***:turbtype and smaw are used to compute 'smm' in microp_aero_ts subroutine but 'smm' is NEVER
839 ! again used for any output from microp_aero_ts. Therefore, these variables are NOT used in any
840 ! meaningful computations
841 turbtype(:,:) = 0._r8
842 smaw(:,:) = 0._r8
843 dlf(:,:) = 0._r8
844 dlf2(:,:) = 0._r8
845 endif
847 !! Assign default size distribution parameters for no-stratiform clouds (convection only)
848 !! Also put into physics buffer for possible separate use by radiation
850 dcon = 25.e-6_r8
851 mucon = 5.3_r8
852 deicon = 50._r8
854 muconv(:,:) = mucon
855 lambdaconv(:,:) = (mucon + 1._r8)/dcon
856 deiconv(:,:) = deicon
858 !BSINGH - CALL Macrophysics
859 call Macrophysics_simple(is_first_step, pcnst, pcols, kte, ixcldliq, ixnumliq, &
860 lchnk, dtime, state1_t, state_pmid, state_q, dp_icwmr, sh_icwmr, dp_frac, &
861 sh_frac, ast,aist,alst,qi_mac,xland_pt,snowh_pt, &
862 !intent-inout
863 state1_s, state1_q, rh_old, cldo, esl, qvs, ptend_loc_name, ptend_loc_ls, &
864 ptend_loc_lq, ptend_loc_s,ptend_loc_q, ptend_all_name, ptend_all_ls, &
865 ptend_all_lq,ptend_all_s,ptend_all_q )
867 do kw = kts , kte
868 kflip = kte-kw+1
870 cldfral(iw,kw,jw) = alst(1,kflip)
871 cldfrai(iw,kw,jw) = aist(1,kflip)
872 CLDFRA_MP(iw,kw,jw) = ast(1,kflip)
873 cld(1,kflip) = ast(1,kflip)
874 CLDFRA_MP_all(iw,kw,jw) = min(ast(1,kflip)+min(0.8_r8,dp_frac(1,kflip)+sh_frac(1,kflip)),1._r8)
875 CLDFRA_CONV(iw,kw,jw) = max(min(0.8_r8,dp_frac(1,kflip)+sh_frac(1,kflip)),0._r8)
876 concld(1,kflip) = CLDFRA_CONV(iw,kw,jw)
877 cldo(1,kflip) = cldfra_old_mp(iw,kw,jw)
878 end do
880 !! ------------------------------------- !
881 !! From here, process computation begins !
882 !! ------------------------------------- !
885 !! ----------------------- !
886 !! Microp_Driver Microphysics !
887 !! ----------------------- !
889 ptend_loc_name = 'microp'
890 ptend_loc_ls = .true.
891 ptend_loc_lq(1) = .true.
892 ptend_loc_lq(ixcldliq) = .true.
893 ptend_loc_lq(ixcldice) = .true.
894 ptend_loc_lq(ixnumliq) = .true.
895 ptend_loc_lq(ixnumice) = .true.
897 #ifndef MODAL_AERO
898 call rad_cnst_get_info( 0, naero = naer_all )
899 allocate( aer_mmr( pcols, pver, naer_all ) )
900 do m = 1, naer_all
901 call rad_cnst_get_aer_mmr( 0, m, state1, pbuf, aermmr1 )
902 aer_mmr(:ncol,:,m) = aermmr1(:ncol,:)
903 enddo
904 #endif
906 #ifdef MODAL_AERO
907 #if ( WRF_CHEM == 1 )
908 do m = 1, ntot_amode
909 lnum = numptr_amode(m)
910 if( lnum > 0 ) then
911 ptend_loc_lq(lnum)= .true.
912 endif
913 do l = 1, nspec_amode(m)
914 lmass = lmassptr_amode(l,m)
915 ptend_loc_lq(lmass)= .true.
916 enddo
917 enddo
918 !Do not update dummy aerosols in state if chem_opt is 0 or CAM_MP_MAM_cpled == .false.
919 if(chem_opt == 0 .or. .NOT.config_flags%CAM_MP_MAM_cpled) then
920 do m = 6 , pcnst !First m = 1 to 5 are water species
921 ptend_loc_lq(m)= .false.
922 enddo
923 endif
924 #endif
925 #endif
927 zeros(:ncol,:pver) = 0._r8
928 qc(:ncol,:pver) = state1_q(:ncol,:pver,ixcldliq)
929 qi(:ncol,:pver) = state1_q(:ncol,:pver,ixcldice)
930 nc(:ncol,:pver) = state1_q(:ncol,:pver,ixnumliq)
931 ni(:ncol,:pver) = state1_q(:ncol,:pver,ixnumice)
933 if( micro_treatment .eq. 'inter' ) then
934 alst_mic(:ncol,:pver) = ast(:ncol,:pver)
935 aist_mic(:ncol,:pver) = ast(:ncol,:pver)
936 elseif( micro_treatment .eq. 'compl' ) then
937 alst_mic(:ncol,:pver) = alst(:ncol,:pver)
938 aist_mic(:ncol,:pver) = aist(:ncol,:pver)
939 endif
941 #if ( WRF_CHEM != 1 )
942 !Adding 7 new fields to state1_q(# mixing ratio and mass mixing ratio
943 !for 3 modes of prescribed aerosols[s04,dust(+sea salt),aitken]) array.
944 !state%q holds 5 constituents already therefore the prescribed aerosol
945 !are defined from 6th to 12th constituents
947 n = modeptr_accum !Adding 7th constituent (# mixing ratio)
948 !tmpnaer = 1000.0e6 !This value is taken from precribe_aerosol_mixactivate! 1.e9 #/kg ~= 1.e3 #/cm3
949 tmpnaer = accum_mode !172.7e6 -For ISDAC case ONLY
950 l = numptr_amode(n)
951 state1_q(:,:,l) = tmpnaer
953 m = 1 !Adding 6th constituent (mass mixing ratio)
954 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
955 l = lmassptr_amode(m,n)
956 state1_q(:,:,l) = tmpmaer
958 n = modeptr_aitken !Adding 9th constituent (# mixing ratio)
959 !tmpnaer = 300.0e6
960 tmpnaer = aitken_mode !0.0e6 -For ISDAC case ONLY
961 l = numptr_amode(n)
962 state1_q(:,:,l) = tmpnaer
964 m = 1 !Adding 8th constituent (mass mixing ratio)
965 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
966 l = lmassptr_amode(m,n)
967 state1_q(:,:,l) = tmpmaer
969 n = modeptr_coarse !Adding 12th constituent (# mixing ratio)
970 !tmpnaer = 0.2e6
971 tmpnaer = coarse_mode !1.1e6 -For ISDAC case ONLY
972 l = numptr_amode(n)
973 state1_q(:,:,l) = tmpnaer
975 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
976 m = 1 !Adding 10th constituent (mass mixing ratio)
977 l = lmassptr_amode(m,n)
978 state1_q(:,:,l) = tmpmaer*0.5
980 m = 2 !Adding 11th constituent (mass mixing ratio)
981 l = lmassptr_amode(m,n)
982 state1_q(:,:,l) = tmpmaer*0.5
983 #else
984 if(chem_opt==0 .or. .NOT. config_flags%CAM_MP_MAM_cpled) then
985 do k = 1, pver
986 n = modeptr_accum !Adding 7th constituent (# mixing ratio)
987 tmpnaer = accum_mode !172.7e6 -For ISDAC case ONLY
988 l = numptr_amode(n)
989 qqcw(1,k,l) = tmpnaer*alst(1,k)*0.8
990 state1_q(1,k,l) = tmpnaer-qqcw(1,k,l)
992 m = 1 !Adding 6th constituent (mass mixing ratio)
993 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
994 l = lmassptr_amode(m,n)
995 qqcw(1,k,l) = tmpmaer*alst(1,k)*0.8
996 state1_q(1,k,l) = tmpmaer-qqcw(1,k,l)
998 n = modeptr_aitken !Adding 9th constituent (# mixing ratio)
999 tmpnaer = aitken_mode !0.0e6 -For ISDAC case ONLY
1000 l = numptr_amode(n)
1001 qqcw(1,k,l) = tmpnaer*alst(1,k)*0.8
1002 state1_q(1,k,l) = tmpnaer-qqcw(1,k,l)
1004 m = 1 !Adding 8th constituent (mass mixing ratio)
1005 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
1006 l = lmassptr_amode(m,n)
1007 qqcw(1,k,l) = tmpmaer*alst(1,k)*0.8
1008 state1_q(1,k,l) = tmpmaer-qqcw(1,k,l)
1009 n = modeptr_coarse !Adding 12th constituent (# mixing ratio)
1010 tmpnaer = coarse_mode !1.1e6 -For ISDAC case ONLY
1011 l = numptr_amode(n)
1012 qqcw(1,k,l) = tmpnaer*alst(1,k)*0.8
1013 state1_q(1,k,l) = tmpnaer-qqcw(1,k,l)
1015 tmpmaer = specdens_amode(lspectype_amode(m,n)) * tmpnaer / voltonumb_amode(n)
1016 m = 1 !Adding 10th constituent (mass mixing ratio)
1017 l = lmassptr_amode(m,n)
1018 qqcw(1,k,l) = tmpmaer*alst(1,k)*0.8*0.5
1019 state1_q(1,k,l) = tmpmaer*0.5-qqcw(1,k,l)
1021 m = 2 !Adding 11th constituent (mass mixing ratio)
1022 l = lmassptr_amode(m,n)
1023 qqcw(1,k,l) = tmpmaer*alst(1,k)*0.8*0.5
1024 state1_q(1,k,l) = tmpmaer*0.5-qqcw(1,k,l)
1025 enddo
1026 endif
1028 #endif
1029 !! calculate aerosol activation (naai for ice, npccn for liquid) and
1030 !! dust size (rndst) and number (nacon) for contact nucleation
1032 ! Output from the following call:wsub, wsubi, naai, naai_hom, npccn, rndst,nacon
1033 ! In-out from the following call: ptend_loc_q, qqcw
1034 !*NOTE*: Added QQCW as an output
1035 call microp_aero_ts ( lchnk, ncol, dtime, state1_t, zeros, &
1036 state1_q(1,1,1), qc, qi, &
1037 nc, ni, state1_pmid, state1_pdel, ast, &
1038 alst_mic, aist_mic, &
1039 cldo, state1_pint, state1_rpdel, state1_zm, state1_omega, &
1040 #ifdef MODAL_AERO
1041 state1_q, cflx, ptend_loc_q, dgnumwet, dgnum, &
1042 #else
1043 aer_mmr, &
1044 #endif
1045 kkvh, tke, turbtype, smaw, wsub, wsubi, &
1046 naai,naai_hom, npccn, rndst,nacon,qqcw)
1048 !! Call MG Microphysics
1051 !Output from the following call:rate1cld,tlat,qvlat,qcten,qiten,
1052 ! ncten,niten,effc,effc_fn,effi,prect,preci,nevapr,evapsnow,
1053 ! prain,prodsnow,cmeice,dei,mu,lambdac,qsout,des,rflx,sflx,
1054 ! qrout,reff_rain,reff_snow,qcsevap,qisevap,qvres,cmeiout,
1055 ! vtrmc,vtrmi,qcsedten,qisedten,prao,prco,mnuccco,mnuccto,
1056 ! msacwio,psacwso,bergso,bergo,melto,homoo,qcreso,prcio,
1057 ! praio,qireso,mnuccro,pracso,meltsdt,frzrdt,mnuccdo
1059 !in-out from the following call: qc,qi,nc,ni,cldo
1060 !*NOTE*: Added nsout and nrout as an ouput
1061 call mmicro_pcond( sub_column, lchnk, ncol, dtime, state1_t, &
1062 state1_q(1,1,1), qc, qi, &
1063 nc, ni, state1_pmid, state1_pdel, ast, &
1064 alst_mic, aist_mic, &
1065 cldo, &
1066 rate1cld, &
1067 naai, npccn, rndst,nacon, &
1068 tlat, qvlat, &
1069 qcten, qiten, ncten, niten, effc, &
1070 effc_fn, effi, prect, preci, &
1071 nevapr, evapsnow, &
1072 prain, prodsnow, cmeice, dei, mu, &
1073 lambdac, qsout, des, &
1074 rflx,sflx, qrout,reff_rain,reff_snow, &
1075 qcsevap, qisevap, qvres, cmeiout, &
1076 vtrmc, vtrmi, qcsedten, qisedten, &
1077 prao, prco, mnuccco, mnuccto, msacwio, psacwso, &
1078 bergso, bergo, melto, homoo, qcreso, prcio, praio, qireso, &
1079 mnuccro, pracso, meltsdt, frzrdt , mnuccdo, nsout, nrout )
1081 do k=1,pver
1082 do i=1,ncol
1083 if (naai(i,k) .gt. 0._r8) then
1084 mnuccdohet(i,k) = mnuccdo(i,k) - (naai_hom(i,k)/naai(i,k))*mnuccdo(i,k)
1085 else
1086 mnuccdohet(i,k) = 0._r8
1087 end if
1088 end do
1089 end do
1091 mgflxprc(:ncol,:pverp) = rflx(:ncol,:pverp) + sflx(:ncol,:pverp)
1092 mgflxsnw(:ncol,:pverp) = sflx(:ncol,:pverp)
1094 mgmrprc(:ncol,:pver) = qrout(:ncol,:pver) + qsout(:ncol,:pver)
1095 mgmrsnw(:ncol,:pver) = qsout(:ncol,:pver)
1098 mgreffrain(:ncol,:pver) = reff_rain(:ncol,:pver)
1099 mgreffsnow(:ncol,:pver) = reff_snow(:ncol,:pver)
1101 !! calculate effective radius of convective liquid and ice using dcon and deicon (not used by code, not useful for COSP)
1102 cvreffliq(:ncol,:pver) = 9.0_r8 !! hard-coded as average of hard-coded values used for deep/shallow convective detrainment (near line 1502)
1103 cvreffice(:ncol,:pver) = 37.0_r8 !! hard-coded as average of hard-coded values used for deep/shallow convective detrainment (near line 1505)
1105 !! Reassign rate1 if modal aerosols
1106 #ifdef MODAL_AERO
1107 rate1ord_cw2pr_st(1:ncol,1:pver)=rate1cld(1:ncol,1:pver)
1108 #endif
1109 !! Sedimentation velocity for liquid stratus cloud droplet
1111 wsedl(:ncol,:pver) = vtrmc(:ncol,:pver)
1113 !! Nominal values for no microp_driver (convective only) cloud.
1114 !! Convert snow mixing ratio to microns
1116 do k = 1, pver
1117 do i = 1, ncol
1118 des(i,k) = des(i,k) * 1.e6_r8
1119 if( ast(i,k) .lt. 1.e-4_r8 ) then
1120 mu(i,k) = mucon
1121 lambdac(i,k) = (mucon + 1._r8)/dcon
1122 dei(i,k) = deicon
1123 endif
1124 end do
1125 end do
1127 !! Net microp_driver condensation rate
1128 !*** IMPORTANT-WARNING ***: Following computation of QME will yeild WRONG answer as CMELIQ is set to zero
1129 qme(:ncol,:pver) = cmeliq(:ncol,:pver) + cmeiout(:ncol,:pver)
1132 #ifndef MODAL_AERO
1133 deallocate(aer_mmr)
1134 #endif
1136 do k = 1, pver
1137 do i = 1, ncol
1138 ptend_loc_s(i,k) = tlat(i,k)
1139 ptend_loc_q(i,k,1) = qvlat(i,k)
1140 ptend_loc_q(i,k,ixcldliq) = qcten(i,k)
1141 ptend_loc_q(i,k,ixcldice) = qiten(i,k)
1142 ptend_loc_q(i,k,ixnumliq) = ncten(i,k)
1143 ptend_loc_q(i,k,ixnumice) = niten(i,k)
1144 enddo
1145 enddo
1147 !! For precip, accumulate only total precip in prec_pwc and snow_pwc variables.
1148 !! Other precip output varirables are set to 0
1149 prec_pcw(:ncol) = prect(:ncol)
1150 snow_pcw(:ncol) = preci(:ncol)
1151 prec_sed(:ncol) = 0._r8
1152 snow_sed(:ncol) = 0._r8
1153 prec_str(:ncol) = prec_pcw(:ncol) + prec_sed(:ncol) - rliq(:ncol)
1154 snow_str(:ncol) = snow_pcw(:ncol) + snow_sed(:ncol)
1156 !! ------------------------------- !
1157 !! Update microphysical tendencies !
1158 !! ------------------------------- !
1160 call physics_ptend_sum(ptend_loc_ls,ptend_loc_lq,ptend_loc_s,ptend_loc_q,ptend_all_ls,ptend_all_lq,ptend_all_s,ptend_all_q,pcnst)
1161 ptend_all_name = 'cldwat'
1162 call physics_update( lchnk,dtime,state1_q,ptend_loc_q,state1_s,ptend_loc_s,ptend_loc_name,ptend_loc_lq,ptend_loc_ls,pcnst)
1163 call physics_ptend_init( ptend_loc_name,ptend_loc_q,ptend_loc_s, ptend_loc_lq,ptend_loc_ls,pcnst)
1165 if( micro_treatment .eq. 'inter' ) then
1166 icecldf(:ncol,:pver) = ast(:ncol,:pver)
1167 liqcldf(:ncol,:pver) = ast(:ncol,:pver)
1168 elseif( micro_treatment .eq. 'compl' ) then
1169 icecldf(:ncol,:pver) = aist(:ncol,:pver)
1170 liqcldf(:ncol,:pver) = alst(:ncol,:pver)
1171 endif
1173 !! Effective droplet radius
1174 rel(:ncol,:pver) = effc(:ncol,:pver)
1175 rel_fn(:ncol,:pver) = effc_fn(:ncol,:pver)
1176 rei(:ncol,:pver) = effi(:ncol,:pver)
1177 rel2(:ncol,:pver) = rel(:ncol,:pver) * 0.9071_r8 !! Convect to effective volume radius assuming pgam = 8
1178 rei2(:ncol,:pver) = rei(:ncol,:pver) * 0.6057_r8 !! Convect to effective volume radius at pgam = 0 for ice
1179 !! ----------------------------------------------------------- !
1180 !! Adjust in-cloud water values to take account of convective !
1181 !! in-cloud water. It is used to calculate the values of !
1182 !! icwlp and iciwp to pass to the radiation. !
1183 !! ----------------------------------------------------------- !
1184 allcld_ice(:ncol,:pver) = 0._r8 !! Grid-avg all cloud liquid
1185 allcld_liq(:ncol,:pver) = 0._r8 !! Grid-avg all cloud ice
1186 if( conv_water_in_rad /= 0 ) then
1187 !Following formulation is obtained from macrop_driver.F90 of CAM
1188 call cldwat_fice( ncol, state1_t, fice, fsnow )
1190 !Balwinder.Singh@pnnl.gov: Changed the argument list to avoid pbuf
1191 call conv_water_4rad( lchnk, ncol, ast, sh_icwmr, dp_icwmr, &
1192 fice, sh_frac, dp_frac, conv_water_in_rad, rei, state1_pdel, &
1193 state1_q(:,:,ixcldliq), state1_q(:,:,ixcldice), &
1194 allcld_liq, allcld_ice )
1195 else
1196 allcld_liq(:ncol,:) = state1_q(:ncol,:,ixcldliq) ! Grid-ave all cloud liquid
1197 allcld_ice(:ncol,:) = state1_q(:ncol,:,ixcldice) ! " ice
1198 end if
1199 !! ------------------------------------------------------------ !
1200 !! Compute in cloud ice and liquid mixing ratios !
1201 !! Note that 'iclwp, iciwp' are used for radiation computation. !
1202 !! ------------------------------------------------------------ !
1203 do k = 1, pver
1204 do i = 1, ncol
1205 !! Limits for in-cloud mixing ratios consistent with MG microphysics
1206 !! in-cloud mixing ratio 0.0001 to 0.005 kg/kg
1207 icimr(i,k) = min( allcld_ice(i,k) / max(0.0001_r8,cld(i,k)),0.005_r8 )
1208 icwmr(i,k) = min( allcld_liq(i,k) / max(0.0001_r8,cld(i,k)),0.005_r8 )
1209 icimrst(i,k) = min( state1_q(i,k,ixcldice) / max(0.0001_r8,icecldf(i,k)),0.005_r8 )
1210 icwmrst(i,k) = min( state1_q(i,k,ixcldliq) / max(0.0001_r8,liqcldf(i,k)),0.005_r8 )
1211 icinc(i,k) = state1_q(i,k,ixnumice) / max(0.0001_r8,icecldf(i,k)) * state1_pmid(i,k) / (287.15_r8*state1_t(i,k))
1212 icwnc(i,k) = state1_q(i,k,ixnumliq) / max(0.0001_r8,liqcldf(i,k)) * state1_pmid(i,k) / (287.15_r8*state1_t(i,k))
1213 iwc(i,k) = allcld_ice(i,k) * state1_pmid(i,k) / (287.15_r8*state1_t(i,k))
1214 lwc(i,k) = allcld_liq(i,k) * state1_pmid(i,k) / (287.15_r8*state1_t(i,k))
1215 effliq(i,k) = effc(i,k)
1216 effliq_fn(i,k) = effc_fn(i,k)
1217 effice(i,k) = effi(i,k)
1218 !! Calculate total cloud water paths in each layer
1219 iciwp(i,k) = icimr(i,k) * state1_pdel(i,k) / gravit
1220 iclwp(i,k) = icwmr(i,k) * state1_pdel(i,k) / gravit
1221 !! Calculate microp_driver cloud water paths in each layer
1222 !! Note: uses stratiform cloud fraction!
1223 iciwpst(i,k) = min(state1_q(i,k,ixcldice)/max(0.0001_r8,ast(i,k)),0.005_r8) * state1_pdel(i,k) / gravit
1224 iclwpst(i,k) = min(state1_q(i,k,ixcldliq)/max(0.0001_r8,ast(i,k)),0.005_r8) * state1_pdel(i,k) / gravit
1225 !! Calculate convective in-cloud LWP.
1226 !Commented the following two diagnostics as 'concld' is not mapped correctly with variables on WRF side
1227 !iclwpconv(i,k) = max(allcld_liq(i,k) - state_q(i,k,ixcldliq),0._r8)/max(0.0001_r8,concld(i,k))
1228 !iciwpconv(i,k) = max(allcld_ice(i,k) - state_q(i,k,ixcldice),0._r8)/max(0.0001_r8,concld(i,k))
1229 !! ------------------------------ !
1230 !! Adjust cloud fraction for snow !
1231 !! ------------------------------ !
1232 cldfsnow(i,k) = cld(i,k)
1233 !! If cloud and only ice ( no convective cloud or ice ), then set to 0.
1234 if( ( cldfsnow(i,k) .gt. 1.e-4_r8 ) .and. &
1235 ( concld(i,k) .lt. 1.e-4_r8 ) .and. &
1236 ( state1_q(i,k,ixcldliq) .lt. 1.e-10_r8 ) ) then
1237 !Commented the following diagnostics as 'concld' is not mapped correctly with variables on WRF side
1238 !cldfsnow(i,k) = 0._r8
1239 endif
1240 !! If no cloud and snow, then set to 0.25
1241 if( ( cldfsnow(i,k) .lt. 1.e-4_r8 ) .and. ( qsout(i,k) .gt. 1.e-6_r8 ) ) then
1242 !Commented the following diagnostics as 'concld' is not mapped correctly with variables on WRF side
1243 !cldfsnow(i,k) = 0.25_r8
1244 endif
1245 !! Calculate in-cloud snow water path
1246 !Commented the following diagnostics as 'concld' is not mapped correctly with variables on WRF side
1247 !icswp(i,k) = qsout(i,k) / max( 0.0001_r8, cldfsnow(i,k) ) * state1_pdel(i,k) / gravit
1248 enddo
1249 enddo
1251 !! --------------------- !
1252 !! History Output Fields !
1253 !! --------------------- !
1255 !! Column droplet concentration
1257 do i = 1, ncol
1258 cdnumc(i) = 0._r8
1259 do k = 1, pver
1260 cdnumc(i) = cdnumc(i) + state1_q(i,k,ixnumliq)*state1_pdel(i,k)/gravit
1261 end do
1262 end do
1264 !! Averaging for new output fields
1266 efcout(:,:) = 10._r8
1267 efiout(:,:) = 25._r8
1268 ncout(:,:) = 0._r8
1269 niout(:,:) = 0._r8
1270 freql(:,:) = 0._r8
1271 freqi(:,:) = 0._r8
1272 liqcldf_out(:,:) = 0._r8
1273 icecldf_out(:,:) = 0._r8
1274 icwmrst_out(:,:) = 0._r8
1275 icimrst_out(:,:) = 0._r8
1276 do k = 1, pver
1277 do i = 1, ncol
1278 if( liqcldf(i,k) .gt. 0.01_r8 .and. icwmrst(i,k) .gt. 5.e-5_r8 ) then
1279 efcout(i,k) = effc(i,k)
1280 ncout(i,k) = icwnc(i,k)
1281 freql(i,k) = 1._r8
1282 liqcldf_out(i,k) = liqcldf(i,k)
1283 icwmrst_out(i,k) = icwmrst(i,k)
1284 endif
1285 if( icecldf(i,k) .gt. 0.01_r8 .and. icimrst(i,k) .gt. 1.e-6_r8 ) then
1286 efiout(i,k) = effi(i,k)
1287 niout(i,k) = icinc(i,k)
1288 freqi(i,k) = 1._r8
1289 icecldf_out(i,k) = icecldf(i,k)
1290 icimrst_out(i,k) = icimrst(i,k)
1291 endif
1292 end do
1293 end do
1295 !! Cloud top effective radius and number.
1297 fcti(:) = 0._r8
1298 fctl(:) = 0._r8
1299 ctrel(:) = 0._r8
1300 ctrei(:) = 0._r8
1301 ctnl(:) = 0._r8
1302 ctni(:) = 0._r8
1303 do i = 1, ncol
1304 do k = 1, pver
1305 if( liqcldf(i,k) .gt. 0.01_r8 .and. icwmrst(i,k) .gt. 1.e-7_r8 ) then
1306 ctrel(i) = effc(i,k)
1307 ctnl(i) = icwnc(i,k)
1308 fctl(i) = 1._r8
1309 exit
1310 endif
1311 if( icecldf(i,k) .gt. 0.01_r8 .and. icimrst(i,k) .gt. 1.e-7_r8 ) then
1312 ctrei(i) = effi(i,k)
1313 ctni(i) = icinc(i,k)
1314 fcti(i) = 1._r8
1315 exit
1316 endif
1317 enddo
1318 enddo
1320 !Following update finally updates the state variables
1321 call physics_update( lchnk,dtime,state_q,ptend_all_q,state_s,ptend_all_s,ptend_all_name,ptend_all_lq,ptend_all_ls,pcnst)
1323 !Post processing of the output from CAMMGMP
1324 do kw=kts,kte
1326 kflip = kte-kw+1
1327 !Following equation are derived following UWPBL and CAMZM schemes
1328 qv_curr(iw,kw,jw) = state_q(1,kflip,1) / (1.0_r8 - state_q(1,kflip,1))
1329 multFrc = 1._r8 + qv_curr(iw,kw,jw)
1331 qc_curr(iw,kw,jw) = state_q(1,kflip,2) * multFrc
1332 qi_curr(iw,kw,jw) = state_q(1,kflip,3) * multFrc
1333 qs_curr(iw,kw,jw) = qsout(1,kflip) * multFrc
1334 qr_curr(iw,kw,jw) = qrout(1,kflip) * multFrc
1336 nc3d(iw,kw,jw) = state_q(1,kflip,4) * multFrc
1337 qndrop(iw,kw,jw) = nc3d(iw,kw,jw) !QNDROP is used in RRTMG radiation scheme
1338 ni3d(iw,kw,jw) = state_q(1,kflip,5) * multFrc
1339 !nr3d ans ns3d are in the units of #/kg but nrout and nsout are in the units of #/m3.
1340 !Therefore these variables are multiplied by ALT (inverse density(m3/kg))
1341 ns3d(iw,kw,jw) = nsout(1,kflip) * alt(iw,kw,jw) * multFrc
1342 nr3d(iw,kw,jw) = nrout(1,kflip) * alt(iw,kw,jw) * multFrc
1344 th(iw,kw,jw) = (state_s(1,kflip)-gravit*z_sea_level(iw,kw,jw))/(cpair*pi_phy(iw,kw,jw))
1345 wsedl3d(iw,kw,jw) = wsedl(1,kflip)
1346 cldfra_old_mp(iw,kw,jw) = ast(1,kflip)
1348 !Update RH - PMA
1349 ttt(1,kflip) = (state_s(1,kflip)-gravit*z_sea_level(iw,kw,jw))/cpair
1350 esl(1,kflip) = polysvp(ttt(1,kflip),0)
1351 qvs(1,kflip) = max(1.e-30_r8,epsqs*esl(1,kflip)/(state_pmid(1,kflip)-(1._r8-epsqs)*esl(1,kflip)))
1352 rh_old_mp(iw,kw,jw) = max(0._r8,state_q(1,kflip,1) / qvs(1,kflip))
1353 lcd_old_mp(iw,kw,jw) = alst(1,kflip)
1354 lradius(iw,kw,jw) = efcout(1,kflip)
1355 iradius(iw,kw,jw) = efiout(1,kflip)
1357 #if ( WRF_CHEM == 1 )
1358 if(chem_opt .NE. 0 .and. config_flags%CAM_MP_MAM_cpled ) then
1359 do l = p1st, num_chem
1360 l2 = lptr_chem_to_q(l)
1361 if ((l2 >= 1) .and. (l2 <= pcnst)) then
1362 chem(iw,kw,jw,l) = state_q(1,kflip,l2)/factconv_chem_to_q(l)
1363 end if
1364 l2 = lptr_chem_to_qqcw(l)
1365 if ((l2 >= 1) .and. (l2 <= pcnst)) then
1366 chem(iw,kw,jw,l) = qqcw(1,kflip,l2)/factconv_chem_to_q(l)
1367 end if
1368 end do ! l
1369 endif
1371 qme3d(iw,kw,jw) = qme(1,kflip)
1372 prain3d(iw,kw,jw) = prain(1,kflip)
1373 nevapr3d(iw,kw,jw) = nevapr(1,kflip)
1374 rate1ord_cw2pr_st3d(iw,kw,jw) = rate1ord_cw2pr_st(1,kflip)
1375 #endif
1376 end do
1378 ! Precipitation(for each time step) is stored in RAINNCV variable of WRF. Following precipitation
1379 ! calculation is taken from 'PRECL'variable of CAM. In 'cam_diagnostic.F90' file of CAM, PRECL is
1380 ! computed as PRECL=PREC_SED+PREC_PCW.
1381 ! PREC_PCW and PREC_SED are output from CAMMGMP (or stratiform.F90 in CAM). Here, we store PREC_SED+PREC_PCW
1382 ! (i.e. PRECL) in RAINNCV variable of WRF. Units of PRECL are m/s. We convert it to mm (millimeter)
1383 ! by multiplying it by 1.e3 and multiplying it by 'DT'
1385 rainncv(iw,jw) = (prec_sed(1) + prec_pcw(1))*dtime*1.0e3_r8 ! in mm
1386 rainnc (iw,jw) = rainnc(iw,jw)+ rainncv(iw,jw)
1388 !Similar to RAINNCV, SNOWNCV is computed following CAM's SNOWL variable in 'cam_diagnostic.F90'
1389 snowncv(iw,jw) = (snow_sed(1) + snow_pcw(1))*dtime*1.0e3_r8 ! in mm
1390 snownc (iw,jw) = snownc(iw,jw) + snowncv(iw,jw)
1392 sr(iw,jw) = snowncv(iw,jw)/(rainncv(iw,jw) + 1.E-12_r8)
1395 enddo !iw loop
1396 enddo !jw loop
1397 end subroutine CAMMGMP
1402 !---------------------------------------------------------------------------------------------
1403 subroutine physics_update(state_lchnk,dt,state_q,ptend_q,state_s,ptend_s,ptend_name,ptend_lq,ptend_ls,pcnst_in)
1404 ! This subroutine partially mimics CAM's subroutine by similar name in physics_type.F90 module
1405 ! Please note that it is implimented just to serve the purpose of porting CAMMGMP scheme
1406 !
1407 !Author: Balwinder.Singh@pnl.gov
1408 !---------------------------------------------------------------------------------------------
1410 implicit none
1411 integer , intent(in) :: state_lchnk, pcnst_in
1412 real(r8), intent(in) :: dt
1414 character*24, intent(inout) :: ptend_name
1416 logical , intent(inout) :: ptend_ls
1417 logical , intent(inout) :: ptend_lq(pcnst_in)
1418 real(r8), intent(inout) :: ptend_s(pcols,pver)
1419 real(r8), intent(inout) :: ptend_q(pcols,pver,pcnst_in)
1421 real(r8), intent(inout) :: state_q(pcols,pver,pcnst_in)
1422 real(r8), intent(inout) :: state_s(pcols,pver)
1424 character*40 :: name ! param and tracer name for qneg3
1425 character(len=16) :: microp_scheme
1426 integer :: m,i,k,ncol
1427 integer :: ixcldice, ixcldliq, ixnumice, ixnumliq
1430 microp_scheme = 'MG'
1431 ptend_top_level = 1
1432 ptend_bot_level = pver
1433 ncol = pcols
1435 ! Update dry static energy
1436 if(ptend_ls) then
1437 do k = ptend_top_level, ptend_bot_level
1438 do i = 1, ncol
1439 state_s(i,k) = state_s(i,k) + ptend_s(i,k) * dt
1440 end do
1441 end do
1442 end if
1444 ! Update constituents, all schemes use time split q: no tendency kept
1445 call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
1446 call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
1447 ! Check for number concentration of cloud liquid and cloud ice (if not present
1448 ! the indices will be set to -1)
1449 call cnst_get_ind('NUMICE', ixnumice, abort=.false.)
1450 call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.)
1452 do m = 1, pcnst_in
1453 if(ptend_lq(m)) then
1454 do k = ptend_top_level, ptend_bot_level
1455 do i = 1,ncol
1456 state_q(i,k,m) = state_q(i,k,m) + ptend_q(i,k,m) * dt
1457 end do
1458 end do
1460 ! now test for mixing ratios which are too small
1461 ! don't call qneg3 for number concentration variables
1462 if (m .ne. ixnumice .and. m .ne. ixnumliq) then
1463 name = trim(ptend_name) // '/' // trim(cnst_name(m))
1464 call qneg3(trim(name), state_lchnk, ncol, pcols, pver, m, m, qmin(m), state_q(1,1,m))
1465 else
1466 do k = ptend_top_level, ptend_bot_level
1467 do i = 1,ncol
1468 ! checks for number concentration
1469 state_q(i,k,m) = max(1.e-12_r8,state_q(i,k,m))
1470 state_q(i,k,m) = min(1.e10_r8,state_q(i,k,m))
1471 end do
1472 end do
1473 end if
1475 end if
1476 end do
1478 ! special tests for cloud liquid
1479 if (ixcldliq > 1) then
1480 if(ptend_lq(ixcldliq)) then
1481 if (ptend_name == 'stratiform' .or. ptend_name == 'cldwat' ) then
1483 else if (ptend_name == 'convect_deep') then
1484 where (state_q(:ncol,:pver,ixcldliq) < 1.e-36_r8)
1485 state_q(:ncol,:pver,ixcldliq) = 0._r8
1486 end where
1487 ! also zero out number concentration
1488 if ( microp_scheme .eq. 'MG' ) then
1489 where (state_q(:ncol,:pver,ixcldliq) < 1.e-36_r8)
1490 state_q(:ncol,:pver,ixnumliq) = 0._r8
1491 end where
1492 end if
1493 end if
1494 end if
1495 end if
1497 ! special tests for cloud ice
1498 if (ixcldice > 1) then
1499 if(ptend_lq(ixcldice)) then
1500 if (ptend_name == 'stratiform' .or. ptend_name == 'cldwat' ) then
1502 else if (ptend_name == 'convect_deep') then
1503 where (state_q(:ncol,:pver,ixcldice) < 1.e-36_r8)
1504 state_q(:ncol,:pver,ixcldice) = 0._r8
1505 end where
1506 ! also zero out number concentration
1507 if ( microp_scheme .eq. 'MG' ) then
1508 where (state_q(:ncol,:pver,ixcldice) < 1.e-36_r8)
1509 state_q(:ncol,:pver,ixnumice) = 0._r8
1510 end where
1511 end if
1512 end if
1513 end if
1514 end if
1516 ! Derive new temperature and geopotential fields if heating or water tendency not 0.
1517 !if (ptend_ls .or. ptend_lq(1)) then
1518 ! call geopotential_dse( &
1519 ! state_lnpint, state_lnpmid, state_pint , state_pmid , state_pdel , state_rpdel , &
1520 ! state_s , state_q(1,1,1),state_phis , rair , gravit , cpair , &
1521 ! zvir , state_t , state_zi , state_zm , ncol )
1522 !end if
1524 ! Reset all parameterization tendency flags to false
1525 call physics_ptend_reset(ptend_name,ptend_q,ptend_s,ptend_lq,ptend_ls,pcnst_in)
1527 end subroutine physics_update
1530 !---------------------------------------------------------------------------------------------
1531 subroutine physics_ptend_init(ptend_name,ptend_q,ptend_s,ptend_lq,ptend_ls,pcnst_in)
1532 ! This subroutine partially mimics CAM's subroutine by similar name in physics_type.F90 module
1533 ! Please note that it is implimented just to serve the purpose of porting CAMMGMP scheme
1534 !Author: Balwinder.Singh@pnl.gov
1535 !---------------------------------------------------------------------------------------------
1536 implicit none
1537 integer, intent(in) :: pcnst_in
1539 character*24, intent(inout) :: ptend_name
1541 logical , intent(inout):: ptend_ls
1542 logical , intent(inout):: ptend_lq(pcnst_in)
1543 real(r8), intent(inout):: ptend_s(pcols,pver)
1544 real(r8), intent(inout):: ptend_q(pcols,pver,pcnst_in)
1546 ptend_name = "none"
1547 ptend_lq(:) = .TRUE.
1548 ptend_ls = .TRUE.
1549 call physics_ptend_reset(ptend_name,ptend_q,ptend_s,ptend_lq,ptend_ls,pcnst_in)
1550 return
1551 end subroutine physics_ptend_init
1554 !---------------------------------------------------------------------------------------------
1555 subroutine physics_ptend_reset(ptend_name,ptend_q,ptend_s,ptend_lq,ptend_ls,pcnst_in)
1556 ! This subroutine partially mimics CAM's subroutine by similar name in physics_type.F90 module
1557 ! Please note that it is implimented just to serve the purpose of porting CAMMGMP scheme
1558 !Author: Balwinder.Singh@pnl.gov
1559 !---------------------------------------------------------------------------------------------
1560 implicit none
1561 integer, intent(in) :: pcnst_in
1563 character*24, intent(inout) :: ptend_name
1565 logical , intent(inout):: ptend_ls
1566 logical , intent(inout):: ptend_lq(pcnst_in)
1567 real(r8), intent(inout):: ptend_s(pcols,pver)
1568 real(r8), intent(inout):: ptend_q(pcols,pver,pcnst_in)
1569 integer :: m
1571 if(ptend_ls) then
1572 ptend_s = 0._r8
1573 endif
1574 do m = 1, pcnst_in
1575 if(ptend_lq(m)) then
1576 ptend_q(:,:,m) = 0._r8
1577 endif
1578 end do
1580 ptend_name = "none"
1581 ptend_lq(:) = .FALSE.
1582 ptend_ls = .FALSE.
1584 ptend_top_level = 1
1585 ptend_bot_level = pver
1587 return
1588 end subroutine physics_ptend_reset
1591 !---------------------------------------------------------------------------------------------
1592 subroutine physics_ptend_sum(ptend_ls,ptend_lq,ptend_s,ptend_q,ptend_sum_ls,ptend_sum_lq,ptend_sum_s,ptend_sum_q, pcnst_in)
1593 ! This subroutine partially mimics CAM's subroutine by similar name in physics_type.F90 module
1594 ! Please note that it is implimented just to serve the purpose of porting CAMMGMP scheme
1595 !Author: Balwinder.Singh@pnl.gov
1596 !---------------------------------------------------------------------------------------------
1598 !-----------------------------------------------------------------------
1599 ! Add ptend fields to ptend_sum for ptend logical flags = .true.
1600 ! Where ptend logical flags = .false, don't change ptend_sum
1601 !-----------------------------------------------------------------------
1603 !------------------------------Arguments--------------------------------
1604 integer, intent(in) :: pcnst_in
1605 logical, intent(in) :: ptend_ls ! New parameterization flag
1606 logical, intent(in) :: ptend_lq(pcnst_in) ! New parameterization flag
1607 real(r8), intent(in) :: ptend_s(pcols,pver) ! New parameterization tendency
1608 real(r8), intent(in) :: ptend_q(pcols,pver,pcnst_in) ! New parameterization tendency
1610 logical, intent(inout) :: ptend_sum_ls ! Flag for ptend_sum
1611 logical, intent(inout) :: ptend_sum_lq(pcnst_in) ! Flag for ptend_sum
1612 real(r8), intent(inout) :: ptend_sum_s(pcols,pver) ! Sum of incoming ptend_sum and ptend
1613 real(r8), intent(inout) :: ptend_sum_q(pcols,pver,pcnst_in) ! Sum of incoming ptend_sum and ptend
1616 !---------------------------Local storage-------------------------------
1617 integer :: k,m ! constituent indices
1618 integer :: ncol,i
1620 !-----------------------------------------------------------------------
1621 ptend_top_level = 1
1622 ptend_bot_level = pver
1624 ncol = pcols
1626 if(ptend_ls) then
1627 ptend_sum_ls = .true.
1628 do i = 1, ncol
1629 do k = ptend_top_level, ptend_bot_level
1630 ptend_sum_s(i,k) = ptend_sum_s(i,k) + ptend_s(i,k)
1631 end do
1632 end do
1633 end if
1635 ! Update constituents
1636 do m = 1, pcnst_in
1637 if(ptend_lq(m)) then
1638 ptend_sum_lq(m) = .true.
1639 do i = 1,ncol
1640 do k = ptend_top_level, ptend_bot_level
1641 ptend_sum_q(i,k,m) = ptend_sum_q(i,k,m) + ptend_q(i,k,m)
1642 end do
1643 end do
1644 end if
1645 end do
1647 end subroutine physics_ptend_sum
1650 !---------------------------------------------------------------------------------------------
1651 subroutine cal_cldfra_mp_1d(cldl,cldi,qqv, qqc, qqi,tt,pp,xland_pt,snowh_pt)
1652 !Author: Po-Lun Ma
1653 !Ported by: Balwinder Singh
1654 !---------------------------------------------------------------------------------------------
1656 use wv_saturation, only: epsqs, polysvp
1657 use shr_kind_mod, only: r8=>shr_kind_r8
1658 use physconst, only: tmelt
1660 real(r8), intent(in) :: qqv
1661 real(r8), intent(in) :: qqc
1662 real(r8), intent(in) :: qqi
1663 real(r8), intent(in) :: tt
1664 real(r8), intent(in) :: pp
1665 real(r8), intent(in) :: xland_pt
1666 real(r8), intent(in) :: snowh_pt
1667 real(r8), intent(inout) :: cldl,cldi
1668 REAL(r8):: rhl,dV,esl,omeps,cldrh,qvl
1669 real(r8):: cldcr,esi,qvi,rhi,rhdif,icimr,qv,qc,qi
1671 qv=max(qqv,0._r8)
1672 qc=max(qqc,0._r8)
1673 qi=max(qqi,0._r8)
1674 cldi= 0._r8
1675 cldl= 0._r8
1678 cldrh=1._r8
1680 if (pp >= 70000._r8) then
1681 if (xland_pt < 1.1_r8 .and. snowh_pt < 0.001_r8) then !xland_pt=1: land, snowh_pt<1e-6m=1e-3kgm-2
1682 cldcr = 0.7875_r8
1683 else
1684 cldcr = 0.8875_r8
1685 endif
1686 elseif (pp <= 40000._r8) then
1687 cldcr = 0.8_r8
1688 else
1689 cldcr = (0.8875_r8 - 0.8_r8)/30000._r8*(pp-40000._r8)+0.8_r8
1690 end if
1692 cldcr=min(max(0._r8,cldcr),1._r8)
1694 dV=cldrh-cldcr
1695 omeps = 1._r8 - 0.622_r8
1696 rhdif = 0._r8
1697 icimr = 0._r8
1699 esl = polysvp(tt,0)
1700 qvl = max(1.e-12_r8,epsqs*esl/(pp-(1._r8-epsqs)*esl))
1701 rhl = (qv+qc)/qvl
1702 rhl = min(1.1_r8,max(0._r8,rhl))
1703 if( rhl .ge. 1._r8 ) then
1704 cldl = 1._r8
1705 elseif( rhl .gt. (cldrh-dV/6._r8) .and. rhl .lt. 1._r8 ) then
1706 cldl = 1._r8 - (-3._r8/sqrt(2._r8)*(rhl-cldrh)/dV)**(2._r8/3._r8)
1707 elseif( rhl .gt. (cldrh-dV) .and. rhl .le. (cldrh-dV/6._r8) ) then
1708 cldl = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* &
1709 (1._r8+(rhl-cldrh)/dV))-2._r8*3.14159_r8)))**2._r8
1710 elseif( rhl .le. (cldrh-dV) ) then
1711 cldl = 0._r8
1712 endif
1714 esi = polysvp(tt,1)
1715 if (tt.gt.tmelt) esi=esl
1716 qvi = max(1.e-12_r8,epsqs*esi/(pp-(1._r8-epsqs)*esi))
1717 rhi = (qv+qi)/qvi
1718 rhi = min(1.1_r8,max(0._r8,rhi))
1720 if (tt < 273.15_r8) then
1721 rhdif = (rhi-0.8_r8) / 0.3_r8 ! rhmaxi=1.1 - rhmini=0.8
1722 cldi = min(1._r8, max(rhdif,0._r8)**2._r8)
1723 if (qi.lt.1.e-12_r8) then
1724 cldi=0._r8
1725 else
1726 cldi=max(1.e-4_r8,cldi)
1727 endif
1729 if (qi.ge.1.e-12_r8) then
1730 icimr=qi/(max(cldi,1.e-20_r8))
1731 if (icimr.lt.1.e-7_r8) then
1732 cldi = max(0._r8,min(1._r8,qi/1.e-7_r8))
1733 endif
1734 if (icimr.gt.5.e-3_r8) then
1735 cldi = max(0._r8,min(1._r8,qi/5.e-3_r8))
1736 endif
1737 endif
1738 endif
1740 cldl=min(max(0._r8,cldl),1._r8)
1741 cldi=min(max(0._r8,cldi),1._r8)
1743 end subroutine cal_cldfra_mp_1d
1746 !---------------------------------------------------------------------------------------------
1747 subroutine Macrophysics_simple(is_first_step, pcnst, pcols, kte, ixcldliq, ixnumliq, &
1748 lchnk, dtime, state1_t, state_pmid, state_q, dp_icwmr, sh_icwmr, dp_frac, &
1749 sh_frac, ast,aist,alst,qi_mac,xland_pt,snowh_pt, &
1750 !intent-inout
1751 state1_s, state1_q, rh_old, cldo, esl, qvs, ptend_loc_name, ptend_loc_ls, &
1752 ptend_loc_lq, ptend_loc_s,ptend_loc_q, ptend_all_name, ptend_all_ls, &
1753 ptend_all_lq,ptend_all_s,ptend_all_q )
1755 ! This subroutine implements macrophysical processes
1756 !
1757 ! Author: Po-Lun Ma
1758 ! Ported by Balwinder Singh
1759 !-------------------------------------------------------------------------------------------
1761 !PMA
1762 !-------------------------------------------------------------------------------------------
1763 !This simple scheme does the followings:
1764 !(1) QC detrained from convections are converted to (a) cloud water and (b) cloud ice in the
1765 ! cloudy part; and (c) vapor in the cloud-free part then latent heating from the phase
1766 ! change
1767 !(2) RH change, then apportion the change of Qv to cloudy part and cloud-free part, then
1768 ! latent heating from phase change
1769 !!(3) If a grid is super-saturated, condense to Qc and update latent heating
1770 !(4) Brutal positive moisture fixer
1771 !
1772 !The idea is to take the essential formulations from Macrophysics in CAM5, but simplify the
1773 !computation. When the full CAM5 macrophysics is implemented, this scheme should be removed
1774 !Cloud fraction calculation is in this module
1775 !-------------------------------------------------------------------------------------------
1777 use wv_saturation, only: epsqs, polysvp
1778 use physconst, only: latvap,cpair
1780 !Subroutine arguments
1782 !Intent -ins
1783 logical, intent(in) :: is_first_step
1784 integer, intent(in) :: pcnst, pcols, kte, ixcldliq, ixnumliq, lchnk
1785 real(r8), intent(in) :: dtime,xland_pt,snowh_pt
1787 real(r8), dimension(pcols,kte), intent(in) :: state_pmid
1788 real(r8), dimension(pcols,kte), intent(in) :: dp_icwmr
1789 real(r8), dimension(pcols,kte), intent(in) :: sh_icwmr
1790 real(r8), dimension(pcols,kte), intent(in) :: dp_frac
1791 real(r8), dimension(pcols,kte), intent(in) :: sh_frac
1792 real(r8), dimension(pcols,kte), intent(in) :: qi_mac
1794 real(r8), dimension(pcols,kte,pcnst), intent(in) :: state_q
1796 !Intents-inouts
1797 character*24, intent(inout) :: ptend_loc_name, ptend_all_name
1799 logical, intent(inout) :: ptend_loc_ls, ptend_all_ls !Flag for updating tendencies
1800 logical, dimension(pcnst),intent(inout) :: ptend_loc_lq !Flag for updating tendencies
1801 logical, dimension(pcnst),intent(inout) :: ptend_all_lq !Flag for updating tendencies
1803 real(r8), dimension(pcols,kte),intent(inout) :: ptend_loc_s
1804 real(r8), dimension(pcols,kte),intent(inout) :: ptend_all_s
1805 real(r8), dimension(pcols,kte),intent(inout) :: rh_old ! RH at previous timestep (PMA)
1806 real(r8), dimension(pcols,kte),intent(inout) :: ast
1807 real(r8), dimension(pcols,kte),intent(inout) :: alst
1808 real(r8), dimension(pcols,kte),intent(inout) :: aist
1809 real(r8), dimension(pcols,kte),intent(inout) :: cldo ! Old cloud fraction
1810 real(r8), dimension(pcols,kte),intent(inout) :: esl ! liquid sat vapor pressure (pa) ! for phil suggested modifications
1811 real(r8), dimension(pcols,kte),intent(inout) :: qvs ! liquid saturation vapor mixing ratio ! for phil suggested modifications
1812 real(r8), dimension(pcols,kte),intent(inout) :: state1_s
1813 real(r8), dimension(pcols,kte),intent(inout) :: state1_t
1815 real(r8), dimension(pcols,kte,pcnst), intent(inout) :: ptend_loc_q
1816 real(r8), dimension(pcols,kte,pcnst), intent(inout) :: ptend_all_q
1817 real(r8), dimension(pcols,kte,pcnst), intent(inout) :: state1_q
1820 !Local workspace
1821 integer :: k
1822 real(r8), dimension(pcols,kte) :: cv_frac ! Total convective cloud fraction (PMA)
1823 real(r8), dimension(pcols,kte) :: st_frac ! Total stratiform cloud fraction (PMA)
1824 real(r8), dimension(pcols,kte) :: ttt ! Temperature (PMA)
1825 real(r8), dimension(pcols,kte) :: rh_curr ! RH at current timestep (PMA)
1826 real(r8), dimension(pcols,kte) :: del_rh ! RH difference between 2 timesteps (PMA)
1827 real(r8), dimension(pcols,kte) :: qc_nc ! non-convective LWC (PMA)
1828 real(r8), dimension(pcols,kte) :: del_qq ! Qv difference between 2 timesteps (PMA)
1829 real(r8), dimension(pcols,kte) :: del_cf ! Qv difference between 2 timesteps (PMA)
1830 real(r8), dimension(pcols,kte) :: rh_temp !
1832 real(r8), dimension(pcols) :: qsll !
1833 real(r8), dimension(pcols) :: tsp !
1834 real(r8) :: cldcr
1835 real(r8) :: factcc,fact_adj
1838 character*250 :: wrfmessage
1841 !
1842 !restore flags to true
1843 !
1845 ptend_loc_name = 'Macro1'
1846 ptend_loc_ls = .true.
1847 ptend_loc_lq(1) = .true.
1848 ptend_loc_lq(ixcldliq) = .true.
1849 ptend_loc_lq(ixnumliq) = .true.
1851 factcc = 0.1_r8
1852 fact_adj = 0.35_r8
1853 cldcr = 0.8875_r8
1854 !
1855 !Simple computation for condensation/evaporation
1856 !
1857 do k=1,pver
1858 if (state_pmid(1,k) >= 70000._r8) then
1859 if (xland_pt < 1.1_r8 .and. snowh_pt < 0.001_r8) then !xland_pt=1: land, snowh_pt<1e-6m=1e-3kgm-2
1860 cldcr = 0.7875_r8
1861 else
1862 cldcr = 0.8875_r8
1863 endif
1864 elseif (state_pmid(1,k)<= 40000._r8) then
1865 cldcr = 0.8_r8
1866 else
1867 cldcr = (0.8875_r8 - 0.8_r8)/30000._r8*(state_pmid(1,k)-40000._r8)+0.8_r8
1868 end if
1870 cldcr=min(max(0._r8,cldcr),1._r8)
1872 call cal_cldfra_mp_1d(alst(1,k),aist(1,k),state1_q(1,k,1), state1_q(1,k,ixcldliq), qi_mac(1,k),state1_t(1,k),state_pmid(1,k),xland_pt,snowh_pt)
1873 qsll = 0._r8
1874 tsp = 0._r8
1875 del_qq(1,k) = 0._r8
1876 rh_temp(1,k) = 0._r8
1877 cv_frac(1,k) = max(0._r8,min(0.8_r8,dp_frac(1,k)+sh_frac(1,k)))
1878 alst(1,k) = (1._r8-cv_frac(1,k))*alst(1,k)
1879 aist(1,k) = (1._r8-cv_frac(1,k))*aist(1,k)
1880 ast(1,k) = max(alst(1,k),aist(1,k))
1881 st_frac(1,k) = min(1._r8,ast(1,k)+cv_frac(1,k))
1882 esl(1,k) = polysvp(state1_t(1,k),0)
1883 qvs(1,k) = max(1.e-30_r8,epsqs*esl(1,k)/(state_pmid(1,k)-(1._r8-epsqs)*esl(1,k)))
1884 rh_curr(1,k) = max(0._r8,state1_q(1,k,1) / qvs(1,k))
1885 qc_nc(1,k) = state1_q(1,k,ixcldliq)
1887 if( is_first_step ) then
1888 rh_old(1,k) = rh_curr(1,k)
1889 cldo(1,k) = min(1._r8,ast(1,k)+cv_frac(1,k))
1890 endif
1892 !
1893 !if super-saturate, condense delta RH
1894 !
1896 if (rh_curr(1,k) > 1._r8) then
1897 call findsp_water(lchnk,pcols,state1_q(1,k,1),state1_t(1,k),state_pmid(1,k),tsp,qsll)
1898 del_qq(1,k) = min(3.e-3_r8,max(0._r8,(state1_q(1,k,1)-qsll(1))*0.999_r8))
1899 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k) /dtime
1900 ptend_loc_s(1,k) = del_qq(1,k) * latvap/dtime !condensational heating
1901 ptend_loc_q(1,k,1) = -del_qq(1,k)/dtime
1903 !
1904 !if not super-saturate, compare RH diff between two timesteps,
1905 !then condense/evaporate the cloudy part
1906 !
1907 else if (rh_curr(1,k) > cldcr .and. rh_curr(1,k) <= 1._r8 ) then
1908 del_cf(1,k) = st_frac(1,k)-cldo(1,k)
1909 if( del_cf(1,k) < 0._r8 ) then
1910 del_qq(1,k) = qc_nc(1,k)*del_cf(1,k)/cldo(1,k)*factcc
1911 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k) /dtime
1912 ptend_loc_s(1,k) = del_qq(1,k) * latvap/dtime !evaporative cooling
1913 ptend_loc_q(1,k,1) = -del_qq(1,k) /dtime
1914 ptend_loc_q(1,k,ixnumliq) = ptend_loc_q(1,k,ixcldliq)*state1_q(1,k,ixnumliq)/max(1.e-30_r8,state1_q(1,k,ixcldliq))
1915 ptend_loc_q(1,k,ixnumliq) = min(0._r8,max(ptend_loc_q(1,k,ixnumliq),-state1_q(1,k,ixnumliq)/dtime))
1916 else if( del_cf(1,k) > 0._r8 .and. cldo(1,k) > 1.e-8_r8) then
1917 del_qq(1,k) = min(state1_q(1,k,1)*0.1_r8,min(3.e-3_r8*st_frac(1,k),qc_nc(1,k)*del_cf(1,k)/cldo(1,k)))*factcc
1918 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k) /dtime
1919 ptend_loc_s(1,k) = del_qq(1,k) * latvap/dtime !condensational heating
1920 ptend_loc_q(1,k,1) = -del_qq(1,k) /dtime
1921 ptend_loc_q(1,k,ixnumliq) = 0._r8
1922 else if( del_cf(1,k) > 0._r8 .and. cldo(1,k) < 1.e-8_r8) then
1923 del_qq(1,k) = min(state1_q(1,k,1)*0.1_r8,2.e-5_r8*st_frac(1,k))*factcc
1924 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k) /dtime
1925 ptend_loc_s(1,k) = del_qq(1,k) * latvap/dtime !condensational heating
1926 ptend_loc_q(1,k,1) = -del_qq(1,k) /dtime
1927 ptend_loc_q(1,k,ixnumliq) = 0._r8
1928 else
1929 del_qq(1,k) = 0._r8
1930 ptend_loc_q(1,k,ixcldliq) = 0._r8
1931 ptend_loc_s(1,k) = 0._r8
1932 ptend_loc_q(1,k,1) = 0._r8
1933 ptend_loc_q(1,k,ixnumliq) = 0._r8
1934 end if
1935 else if (rh_curr(1,k) < cldcr .and. cldo(1,k) <1.e-2_r8 .and. qc_nc(1,k) >1.e-12_r8) then !RH<cldcr, evaporate qc
1936 call findsp_water(lchnk,pcols,state1_q(1,k,1),state1_t(1,k),state_pmid(1,k),tsp,qsll)
1937 del_qq(1,k)= max(min(qc_nc(1,k), (qsll(1)-state1_q(1,k,1))),0._r8)*0.99_r8
1938 ptend_loc_q(1,k,ixcldliq) = -del_qq(1,k) /dtime
1939 ptend_loc_s(1,k) = -del_qq(1,k) * latvap/dtime !evaporative cooling
1940 ptend_loc_q(1,k,1) = del_qq(1,k)/dtime
1941 ptend_loc_q(1,k,ixnumliq) = ptend_loc_q(1,k,ixcldliq)*state1_q(1,k,ixnumliq)/max(1.e-30_r8,state1_q(1,k,ixcldliq))
1942 ptend_loc_q(1,k,ixnumliq) = min(0._r8,max(ptend_loc_q(1,k,ixnumliq),-state1_q(1,k,ixnumliq)/dtime))
1943 else if (rh_curr(1,k) < cldcr .and. cldo(1,k) >=1.e-2_r8 .and. qc_nc(1,k) >1.e-12_r8) then !RH<cldcr, evaporate qc
1944 call findsp_water(lchnk,pcols,state1_q(1,k,1),state1_t(1,k),state_pmid(1,k),tsp,qsll)
1945 del_qq(1,k)= max(min(qc_nc(1,k), qsll(1)-state1_q(1,k,1)),0._r8)*factcc
1946 ptend_loc_q(1,k,ixcldliq) = -del_qq(1,k) /dtime
1947 ptend_loc_s(1,k) = -del_qq(1,k) * latvap/dtime !evaporative cooling
1948 ptend_loc_q(1,k,1) = del_qq(1,k)/dtime
1949 ptend_loc_q(1,k,ixnumliq) = ptend_loc_q(1,k,ixcldliq)*state1_q(1,k,ixnumliq)/max(1.e-30_r8,state1_q(1,k,ixcldliq))
1950 ! ptend_loc_q(1,k,ixnumliq) = ptend_loc_q(1,k,ixcldliq)*state1_q(1,k,ixnumliq)/max(1.e-30_r8,qc_nc(1,k))
1951 ptend_loc_q(1,k,ixnumliq) = min(0._r8,max(ptend_loc_q(1,k,ixnumliq),-state1_q(1,k,ixnumliq)/dtime))
1952 else
1953 ptend_loc_q(1,k,ixcldliq) = 0._r8
1954 ptend_loc_s(1,k) = 0._r8
1955 ptend_loc_q(1,k,1) = 0._r8
1956 ptend_loc_q(1,k,ixnumliq) = 0._r8
1957 end if
1958 state1_t(1,k)=state1_t(1,k)+ptend_loc_s(1,k)*dtime/cpair
1959 end do !k loop
1961 !
1962 !Update state1 variables and tendencies
1963 !
1965 call physics_ptend_sum(ptend_loc_ls,ptend_loc_lq,ptend_loc_s,ptend_loc_q,ptend_all_ls,ptend_all_lq,ptend_all_s,ptend_all_q,pcnst)
1966 ptend_all_name = 'cldwat'
1967 call physics_update( lchnk,dtime,state1_q,ptend_loc_q,state1_s,ptend_loc_s,ptend_loc_name,ptend_loc_lq,ptend_loc_ls,pcnst)
1968 call physics_ptend_init( ptend_loc_name,ptend_loc_q,ptend_loc_s, ptend_loc_lq,ptend_loc_ls,pcnst)
1970 !
1971 !restore flags to true
1972 !
1973 ptend_loc_name = 'Macro2'
1974 ptend_loc_ls = .true.
1975 ptend_loc_lq(1) = .true.
1976 ptend_loc_lq(ixcldliq) = .true.
1977 ptend_loc_lq(ixnumliq) = .true.
1979 !
1980 !Some limiters
1981 !
1983 do k = 1 , pver
1984 del_qq(1,k)=0._r8
1985 call cal_cldfra_mp_1d(alst(1,k),aist(1,k),state1_q(1,k,1), state1_q(1,k,ixcldliq), qi_mac(1,k),state1_t(1,k),state_pmid(1,k),xland_pt,snowh_pt)
1986 alst(1,k) = (1._r8-cv_frac(1,k))*alst(1,k)
1987 aist(1,k) = (1._r8-cv_frac(1,k))*aist(1,k)
1988 ast(1,k) = max(alst(1,k),aist(1,k))
1989 st_frac(1,k) = alst(1,k)
1991 qc_nc(1,k) = max(0._r8, state1_q(1,k,ixcldliq) - &
1992 max(0._r8,(dp_frac(1,k)*dp_icwmr(1,k) +sh_frac(1,k)*sh_icwmr(1,k))))
1993 !
1994 !Case 1: unrealistic cloud: evaporate qc under no cldfra
1995 !
1996 ! THIS SHOULD NO LONGER HAPPEN
1997 !
1998 ! if (st_frac(1,k) < 1.e-8_r8 ) then
1999 ! if ( qc_nc(1,k) > 0._r8 ) then
2000 ! ptend_loc_q(1,k,ixcldliq) = -(qc_nc(1,k))/dtime*0.99_r8
2001 ! ptend_loc_q(1,k,1)= qc_nc(1,k)/dtime*0.99_r8
2002 ! ptend_loc_s(1,k) =-(qc_nc(1,k)) * latvap /dtime*0.99_r8 !evaporative cooling
2003 ! ptend_loc_q(1,k,ixnumliq) = -state1_q(1,k,ixnumliq)/dtime*0.99_r8
2004 ! ptend_loc_q(1,k,ixnumliq) = min(0._r8,-0.99_r8*state1_q(1,k,ixnumliq)/dtime)
2005 ! end if
2007 if (st_frac(1,k) > 1.e-5_r8 ) then
2008 ! else !if stratiform cldfra not zero
2009 !
2010 !Case 2: dense cloud: evaporate some stratiform qc until qc_max
2011 !
2012 if (qc_nc(1,k) > 3.e-3_r8 * st_frac(1,k)) then
2013 del_qq(1,k) = qc_nc(1,k)-3.e-3_r8* st_frac(1,k)*fact_adj
2014 ptend_loc_q(1,k,ixcldliq) = -del_qq(1,k)/dtime
2015 ptend_loc_q(1,k,1) = del_qq(1,k)/dtime
2016 ptend_loc_s(1,k) = -del_qq(1,k) * latvap /dtime !evaporative cooling
2017 ptend_loc_q(1,k,ixnumliq) = ptend_loc_q(1,k,ixcldliq)*state1_q(1,k,ixnumliq)/max(1.e-20_r8,state1_q(1,k,ixcldliq))
2018 ptend_loc_q(1,k,ixnumliq) = min(0._r8,max(ptend_loc_q(1,k,ixnumliq),-state1_q(1,k,ixnumliq)/dtime))
2020 !
2021 !Case 3: empty cloud: condense stratiform qc until qc_min when cldfra is not zero
2022 !
2023 else if (qc_nc(1,k) < 2.e-5_r8 * st_frac(1,k)) then
2024 del_qq(1,k) = 2.e-5_r8 *st_frac(1,k)-qc_nc(1,k)
2025 del_qq(1,k) = max(0._r8,min(del_qq(1,k), state1_q(1,k,1)*0.1_r8))*fact_adj
2026 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k)/dtime
2027 ptend_loc_q(1,k,1) = -del_qq(1,k)/dtime
2028 ptend_loc_s(1,k) = del_qq(1,k)* latvap / dtime !condensational heating
2029 end if
2030 end if
2031 state1_t(1,k)=state1_t(1,k)+ptend_loc_s(1,k)*dtime/cpair
2032 end do
2033 !
2034 !Update state1 variables and tendencies
2035 !
2036 call physics_ptend_sum(ptend_loc_ls,ptend_loc_lq,ptend_loc_s,ptend_loc_q,ptend_all_ls,ptend_all_lq,ptend_all_s,ptend_all_q,pcnst)
2037 ptend_all_name = 'cldwat'
2038 call physics_update( lchnk,dtime,state1_q,ptend_loc_q,state1_s,ptend_loc_s,ptend_loc_name,ptend_loc_lq,ptend_loc_ls,pcnst)
2039 call physics_ptend_init( ptend_loc_name,ptend_loc_q,ptend_loc_s, ptend_loc_lq,ptend_loc_ls,pcnst)
2041 go to 2000
2042 !
2043 !restore flags to true
2044 !
2046 ptend_loc_name = 'Macro3'
2047 ptend_loc_ls = .true.
2048 ptend_loc_lq(1) = .true.
2049 ptend_loc_lq(ixcldliq) = .true.
2051 !
2052 !Last checker to remove supersaturation
2053 !
2054 do k=1,pver
2055 qsll = 0._r8
2056 tsp = 0._r8
2057 del_qq(1,k) = 0._r8
2058 esl(1,k) = polysvp(state1_t(1,k),0)
2059 qvs(1,k) = max(1.e-30_r8,epsqs*esl(1,k)/(state_pmid(1,k)-(1._r8-epsqs)*esl(1,k)))
2060 rh_curr(1,k) = state1_q(1,k,1) / qvs(1,k)
2062 if (rh_curr(1,k) > 1._r8) then
2063 call findsp_water(lchnk,pcols,state1_q(1,k,1),state1_t(1,k),state_pmid(1,k),tsp,qsll)
2064 del_qq(1,k)=min(2.e-2_r8,max(0._r8,(state1_q(1,k,1)-qsll(1))*0.999_r8))
2065 ptend_loc_q(1,k,ixcldliq) = del_qq(1,k) /dtime
2066 ptend_loc_s(1,k) = del_qq(1,k) * latvap/dtime !condensational heating
2067 ptend_loc_q(1,k,1) = -del_qq(1,k)/dtime
2068 end if
2069 state1_t(1,k)=state1_t(1,k)+ptend_loc_s(1,k)*dtime/cpair
2070 end do
2072 !Update state1 variables and tendencies
2073 !
2074 call physics_ptend_sum(ptend_loc_ls,ptend_loc_lq,ptend_loc_s,ptend_loc_q,ptend_all_ls,ptend_all_lq,ptend_all_s,ptend_all_q,pcnst)
2075 ptend_all_name = 'cldwat'
2076 call physics_update( lchnk,dtime,state1_q,ptend_loc_q,state1_s,ptend_loc_s,ptend_loc_name,ptend_loc_lq,ptend_loc_ls,pcnst)
2077 call physics_ptend_init( ptend_loc_name,ptend_loc_q,ptend_loc_s, ptend_loc_lq,ptend_loc_ls,pcnst)
2079 !
2080 !positive state1 variables
2081 !
2082 2000 continue
2084 do k = 1 , pver
2085 state1_q(1,k,1) = max(qmin(1),state1_q(1,k,1))
2086 state1_q(1,k,ixcldliq) = max(qmin(ixcldliq),state1_q(1,k,ixcldliq))
2087 state1_q(1,k,ixcldice) = max(qmin(ixcldice),state1_q(1,k,ixcldice))
2088 state1_q(1,k,ixnumliq) = max(qmin(ixnumliq),state1_q(1,k,ixnumliq))
2089 state1_q(1,k,ixnumice) = max(qmin(ixnumice),state1_q(1,k,ixnumice))
2090 call cal_cldfra_mp_1d(alst(1,k),aist(1,k),state1_q(1,k,1), state1_q(1,k,ixcldliq), qi_mac(1,k),state1_t(1,k),state_pmid(1,k),xland_pt,snowh_pt)
2091 ! ast(1,k)=max(alst(1,k),aist(1,k))
2092 aist(1,k) = (1._r8-cv_frac(1,k))*aist(1,k)
2093 alst(1,k) = (1._r8-cv_frac(1,k))*alst(1,k)
2094 ast(1,k)=max(alst(1,k),aist(1,k))
2095 end do
2096 !
2097 !PMA<<<
2098 !
2099 end subroutine Macrophysics_simple
2101 subroutine CAMMGMP_INIT(ixcldliq_in,ixcldice_in &
2102 ,ixnumliq_in,ixnumice_in,chem_opt_in &
2103 ,ids, ide, jds, jde, kds, kde &
2104 ,ims, ime, jms, jme, kms, kme &
2105 ,its, ite, jts, jte, kts, kte )
2107 !!-------------------------------------------- !
2108 !! !
2109 !! Initialize the cloud water parameterization !
2110 !! !
2111 !!-------------------------------------------- !
2112 use microp_aero, only: ini_microp_aero
2113 use cldwat2m_micro, only: ini_micro
2114 #ifdef MODAL_AERO
2115 use ndrop, only: activate_init
2116 #endif
2118 !!-----------------------------------------------------------------------
2119 implicit none
2120 integer, intent(in) :: ixcldliq_in, ixcldice_in, ixnumliq_in, ixnumice_in
2121 integer, intent(in) :: chem_opt_in
2122 integer, intent(in) :: ids, ide, jds, jde, kds, kde &
2123 ,ims, ime, jms, jme, kms, kme &
2124 ,its, ite, jts, jte, kts, kte
2126 !Local variables
2127 character(len=1000) :: msg
2128 integer :: jtf,ktf,itf
2130 chem_opt = chem_opt_in
2132 jtf = min(jte,jde-1)
2133 ktf = min(kte,kde-1)
2134 itf = min(ite,ide-1)
2136 !Map CAM veritcal level variables
2137 pver = ktf - kts + 1
2138 pverp = pver + 1
2140 !!-----------------------------------------------------------------------
2141 ixcldliq = ixcldliq_in
2142 ixcldice = ixcldice_in
2143 ixnumliq = ixnumliq_in
2144 ixnumice = ixnumice_in
2145 !! Initialization routine for cloud macrophysics and microphysics
2147 call ini_micro
2148 call ini_microp_aero
2150 #ifdef MODAL_AERO
2151 #if ( WRF_CHEM != 1 )
2152 !When WRF_CHEM is 1, activate_init is called from MODULE_CAM_MAM_INIT after initializing aerosols
2153 call activate_init
2154 #else
2155 !BSINGH:02/01/2013: Sanity check for Non-MAM simulations
2156 if(.NOT.cam_mam_aerosols .AND. chem_opt .NE. 0) then
2157 write(msg,*)'CAMMGMP DRIVER - cammgmp is valid for only MAM aerosols ', &
2158 '(chem_opts:',CBMZ_CAM_MAM3_NOAQ,CBMZ_CAM_MAM3_AQ,CBMZ_CAM_MAM7_NOAQ,CBMZ_CAM_MAM7_AQ ,')'
2159 call wrf_error_fatal( msg )
2160 endif
2161 if(chem_opt == 0)then
2162 !NOTE*: Not tested yet for other CHEM packages except MAM packages
2163 call activate_init
2164 endif
2165 #endif
2166 #endif
2167 return
2168 end subroutine CAMMGMP_INIT
2169 subroutine findsp_water (lchnk, ncol, q, t, p, tsp, qsp)
2170 !-----------------------------------------------------------------------
2171 !
2172 ! Purpose:
2173 ! find the wet bulb temperature for a given t and q
2174 ! in a longitude height section
2175 ! wet bulb temp is the temperature and spec humidity that is
2176 ! just saturated and has the same enthalpy
2177 ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
2178 ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
2179 !
2180 ! Method:
2181 ! a Newton method is used
2182 ! first guess uses an algorithm provided by John Petch from the UKMO
2183 ! we exclude points where the physical situation is unrealistic
2184 ! e.g. where the temperature is outside the range of validity for the
2185 ! saturation vapor pressure, or where the water vapor pressure
2186 ! exceeds the ambient pressure, or the saturation specific humidity is
2187 ! unrealistic
2188 !
2189 ! Author: P. Rasch
2190 !
2191 !-----------------------------------------------------------------------
2193 use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, &
2194 cp, epsqs, ttrice,polysvp
2196 !
2197 ! input arguments
2198 !
2199 integer, intent(in) :: lchnk ! chunk identifier
2200 integer, intent(in) :: ncol ! number of atmospheric columns
2202 real(r8), intent(in) :: q(pcols) ! water vapor (kg/kg)
2203 real(r8), intent(in) :: t(pcols) ! temperature (K)
2204 real(r8), intent(in) :: p(pcols) ! pressure (Pa)
2205 !
2206 ! output arguments
2207 !
2208 real(r8), intent(out) :: tsp(pcols) ! saturation temp (K)
2209 real(r8), intent(out) :: qsp(pcols) ! saturation mixing ratio (kg/kg)
2210 !
2211 ! local variables
2212 !
2213 character*250 :: iulog
2214 integer i ! work variable
2215 ! integer k ! work variable
2216 logical lflg ! work variable
2217 integer iter ! work variable
2218 integer l ! work variable
2219 logical :: error_found
2221 real(r8) omeps ! 1 minus epsilon
2222 real(r8) trinv ! work variable
2223 real(r8) es ! sat. vapor pressure
2224 real(r8) desdt ! change in sat vap pressure wrt temperature
2225 ! real(r8) desdp ! change in sat vap pressure wrt pressure
2226 real(r8) dqsdt ! change in sat spec. hum. wrt temperature
2227 real(r8) dgdt ! work variable
2228 real(r8) g ! work variable
2229 real(r8) weight(pcols) ! work variable
2230 real(r8) hlatsb ! (sublimation)
2231 real(r8) hlatvp ! (vaporization)
2232 real(r8) hltalt(pcols) ! lat. heat. of vap.
2233 real(r8) tterm ! work var.
2234 real(r8) qs ! spec. hum. of water vapor
2235 real(r8) tc ! crit temp of transition to ice
2236 real(r8) tt0 ! Freezing temperature. Needed for findsp
2238 ! work variables
2239 real(r8) t1, q1, dt, dq
2240 real(r8) dtm, dqm
2241 real(r8) qvd, a1, tmp
2242 real(r8) rair
2243 real(r8) r1b, c1, c2, c3
2244 real(r8) denom
2245 real(r8) dttol
2246 real(r8) dqtol
2247 integer doit(pcols)
2248 real(r8) enin(pcols), enout(pcols)
2249 real(r8) tlim(pcols)
2251 omeps = 1.0_r8 - epsqs
2252 ! trinv = 1.0_r8/ttrice ! put just before it is used, in case ttrice = 0.0 !++xl
2253 a1 = 7.5_r8*log(10._r8)
2254 rair = 287.04_r8
2255 c3 = rair*a1/cp
2256 dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei
2257 dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei
2258 dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
2259 dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
2260 ! tmin = 173.16 ! the coldest temperature we can deal with
2261 tt0 = 273.15_r8
2262 !
2263 ! max number of times to iterate the calculation
2264 iter = 8
2265 !
2267 ! Sungsu comment out k loop
2268 ! do k = 1,pver
2270 !
2271 ! first guess on the wet bulb temperature
2272 !
2273 do i = 1,ncol
2275 !#ifdef DEBUG
2276 ! if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
2277 ! write (iulog,*) ' '
2278 ! write (iulog,*) ' level, t, q, p', t(i), q(i), p(i)
2279 ! endif
2280 !#endif
2281 ! limit the temperature range to that relevant to the sat vap pres tables
2282 !#if ( ! defined WACCM_MOZART )
2283 ! tlim(i) = min(max(t(i),173._r8),373._r8)
2284 !#else
2285 tlim(i) = min(max(t(i),128._r8),373._r8)
2286 !#endif
2287 ! es = estblf(tlim(i))
2288 es = polysvp(tlim(i),0) !++xl
2289 denom = p(i) - omeps*es
2290 qs = epsqs*es/denom
2291 doit(i) = 0
2292 enout(i) = 1._r8
2293 ! make sure a meaningful calculation is possible
2294 if (p(i) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
2295 !
2296 ! Saturation specific humidity
2297 !
2298 qs = min(epsqs*es/denom,1._r8)
2299 !
2300 ! "generalized" analytic expression for t derivative of es
2301 ! accurate to within 1 percent for 173.16 < t < 373.16
2303 !++xl
2304 ! Weighting of hlat accounts for transition from water to ice
2305 ! polynomial expression approximates difference between es over
2306 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
2307 ! -40): required for accurate estimate of es derivative in transition
2308 ! range from ice to water also accounting for change of hlatv with t
2309 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
2310 !
2311 ! tc = tlim(i) - tt0
2312 ! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
2313 ! weight(i) = min(-tc*trinv,1.0_r8)
2314 ! hlatsb = hlatv + weight(i)*hlatf
2315 ! hlatvp = hlatv - 2369.0_r8*tc
2316 ! if (tlim(i) < tt0) then
2317 ! hltalt(i) = hlatsb
2318 ! else
2319 ! hltalt(i) = hlatvp
2320 ! end if
2321 !
2322 ! No icephs or water to ice transition
2323 !
2324 hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
2325 hlatsb = hlatv
2326 if (tlim(i) < tt0) then
2327 hltalt(i) = hlatsb
2328 else
2329 hltalt(i) = hlatvp
2330 end if
2331 !--xl
2332 enin(i) = cp*tlim(i) + hltalt(i)*q(i)
2334 ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
2335 tmp = q(i) - qs
2336 c1 = hltalt(i)*c3
2337 c2 = (tlim(i) + 36._r8)**2
2338 r1b = c2/(c2 + c1*qs)
2339 qvd = r1b*tmp
2340 tsp(i) = tlim(i) + ((hltalt(i)/cp)*qvd)
2341 !#ifdef DEBUG
2342 ! if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
2343 ! write (iulog,*) ' relative humidity ', q(i)/qs
2344 ! write (iulog,*) ' first guess ', tsp(i)
2345 ! endif
2346 !#endif
2347 ! es = estblf(tsp(i))
2348 es = polysvp(tsp(i),0) !++xl
2349 qsp(i) = min(epsqs*es/(p(i) - omeps*es),1._r8)
2350 else
2351 doit(i) = 1
2352 tsp(i) = tlim(i)
2353 qsp(i) = q(i)
2354 enin(i) = 1._r8
2355 endif
2356 end do ! end do i
2357 !
2358 ! now iterate on first guess
2359 !
2360 do l = 1, iter
2361 dtm = 0
2362 dqm = 0
2363 do i = 1,ncol
2364 if (doit(i) == 0) then
2365 ! es = estblf(tsp(i))
2366 es = polysvp(tsp(i),0) !++xl
2367 !
2368 ! Saturation specific humidity
2369 !
2370 qs = min(epsqs*es/(p(i) - omeps*es),1._r8)
2371 !
2372 ! "generalized" analytic expression for t derivative of es
2373 ! accurate to within 1 percent for 173.16 < t < 373.16
2374 !
2375 !++xl
2376 ! Weighting of hlat accounts for transition from water to ice
2377 ! polynomial expression approximates difference between es over
2378 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
2379 ! -40): required for accurate estimate of es derivative in transition
2380 ! range from ice to water also accounting for change of hlatv with t
2381 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
2382 !
2383 ! tc = tsp(i) - tt0
2384 ! lflg = (tc >= -ttrice .and. tc < 0.0_r8)
2385 ! weight(i) = min(-tc*trinv,1.0_r8)
2386 ! hlatsb = hlatv + weight(i)*hlatf
2387 ! hlatvp = hlatv - 2369.0_r8*tc
2388 ! if (tsp(i) < tt0) then
2389 ! hltalt(i) = hlatsb
2390 ! else
2391 ! hltalt(i) = hlatvp
2392 ! end if
2393 ! if (lflg) then
2394 ! tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
2395 ! else
2396 ! tterm = 0.0_r8
2397 ! end if
2398 ! desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i)) + tterm*trinv
2399 !
2400 ! No icephs or water to ice transition
2401 !
2402 hlatvp = hlatv - 2369.0*(tsp(i)-tt0)
2403 hlatsb = hlatv
2404 if (tsp(i) < tt0) then
2405 hltalt(i) = hlatsb
2406 else
2407 hltalt(i) = hlatvp
2408 end if
2409 desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i))
2410 !--xl
2411 dqsdt = (epsqs + omeps*qs)/(p(i) - omeps*es)*desdt
2412 ! g = cp*(tlim(i)-tsp(i)) + hltalt(i)*q(i)- hltalt(i)*qsp(i)
2413 g = enin(i) - (cp*tsp(i) + hltalt(i)*qsp(i))
2414 dgdt = -(cp + hltalt(i)*dqsdt)
2415 t1 = tsp(i) - g/dgdt
2416 dt = abs(t1 - tsp(i))/t1
2417 tsp(i) = max(t1,tmin)
2418 ! es = estblf(tsp(i))
2419 es = polysvp(tsp(i),0) !++xl
2420 q1 = min(epsqs*es/(p(i) - omeps*es),1._r8)
2421 dq = abs(q1 - qsp(i))/max(q1,1.e-12_r8)
2422 qsp(i) = q1
2423 !#ifdef DEBUG
2424 ! if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
2425 ! write (iulog,*) ' rel chg lev, iter, t, q ', l, dt, dq, g
2426 ! endif
2427 !#endif
2428 dtm = max(dtm,dt)
2429 dqm = max(dqm,dq)
2430 ! if converged at this point, exclude it from more iterations
2431 if (dt < dttol .and. dq < dqtol) then
2432 doit(i) = 2
2433 endif
2434 enout(i) = cp*tsp(i) + hltalt(i)*qsp(i)
2435 ! bail out if we are too near the end of temp range
2436 !#if ( ! defined WACCM_MOZART )
2437 if (tsp(i) < 174.16_r8) then
2438 !#else
2439 ! if (tsp(i) < 130.16_r8) then
2440 !#endif
2441 doit(i) = 4
2442 endif
2443 else
2444 endif
2445 end do ! do i = 1,ncol
2447 if (dtm < dttol .and. dqm < dqtol) then
2448 go to 10
2449 endif
2451 end do ! do l = 1,iter
2452 10 continue
2454 error_found = .false.
2455 if (dtm > dttol .or. dqm > dqtol) then
2456 do i = 1,ncol
2457 if (doit(i) == 0) error_found = .true.
2458 end do
2459 if (error_found) then
2460 do i = 1,ncol
2461 if (doit(i) == 0) then
2462 write (iulog,*) ' findsp not converging at point i ', i
2463 write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
2464 write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
2465 call wrf_error_fatal ('FINDSP')
2466 endif
2467 end do
2468 endif
2469 endif
2470 do i = 1,ncol
2471 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
2472 error_found = .true.
2473 endif
2474 end do
2475 if (error_found) then
2476 do i = 1,ncol
2477 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
2478 write (iulog,*) ' the enthalpy is not conserved for point ', &
2479 i, enin(i), enout(i)
2480 write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
2481 write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
2482 call wrf_error_fatal ('FINDSP')
2483 endif
2484 end do
2485 endif
2487 ! end do ! level loop (k=1,pver)
2489 return
2490 end subroutine findsp_water
2491 end module module_mp_cammgmp_driver