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