| blob | b44b8ec42e1c315224d68f737d6c9a1fa2e8d88b |
1 #undef DEBUG
2 #define WRF_PORT
3 #define MODAL_AERO
4 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
5 module cldwat
6 !-----------------------------------------------------------------------
7 !
8 ! Purpose: Prognostic cloud water data and methods.
9 !
10 ! Public interfaces:
11 !
12 ! inimc -- Initialize constants
13 ! pcond -- Calculate prognostic condensate
14 !
15 ! cldwat_fice -- calculate fraction of condensate in ice phase (radiation partitioning)
16 !
17 ! Author: P. Rasch, with Modifications by Minghua Zhang
18 ! January 2010, modified by J. Kay to add precip fluxes for COSP simulator
19 ! Ported to WRF by William.Gustafson@pnl.gov, Nov. 2009
20 ! updated to CESM1_0_1, Nov. 2010
21 !
22 !-----------------------------------------------------------------------
23 use shr_kind_mod, only: r8 => shr_kind_r8
24 #ifndef WRF_PORT
25 use spmd_utils, only: masterproc
26 use ppgrid, only: pcols, pver, pverp
27 #endif
28 use wv_saturation, only: estblf, hlatv, tmin, hlatf, rgasv, pcf, &
29 cp, epsqs, ttrice
30 #ifndef WRF_PORT
31 use abortutils, only: endrun
32 use cam_logfile, only: iulog
33 #else
34 use module_cam_support, only: masterproc, &
35 pcols, pver, pverp, &
36 endrun, &
37 iulog
38 #endif
40 implicit none
42 !-----------------------------------------------------------------------
43 ! PUBLIC: Make default data and interfaces private
44 !-----------------------------------------------------------------------
45 private
46 save
47 #ifndef WRF_PORT
48 public inimc, pcond ! Public interfaces
49 #endif
50 public cldwat_fice ! Public interfaces
51 public cldwat_readnl
52 integer, public:: ktop ! Level above 10 hPa
53 #ifndef WRF_PORT
54 real(r8),public :: icritc ! threshold for autoconversion of cold ice
55 real(r8),public :: icritw ! threshold for autoconversion of warm ice
56 !!$ real(r8),public,parameter:: conke = 1.e-6 ! tunable constant for evaporation of precip
57 !!$ real(r8),public,parameter:: conke = 2.e-6 ! tunable constant for evaporation of precip
58 real(r8),public :: conke ! tunable constant for evaporation of precip
59 real(r8),public :: r3lcrit ! critical radius where liq conversion begins
60 #else
61 ! Currently, the WRF_PORT bipasses the namelist initialization of these
62 ! tunable parameters. We are hard-coding them to the default values for
63 ! the fv 0.23x0.31 grid. One might choose to implement an option to set
64 ! these via the WRF Registry in the future.
65 real(r8),public :: icritw = 2.0e-4
66 real(r8),public :: icritc = 45.0e-6
67 real(r8),public :: conke = 5.0e-6
69 #endif
71 !-----------------------------------------------------------------------
72 ! PRIVATE: Everything else is private to this module
73 !-----------------------------------------------------------------------
74 real(r8), private:: tmax_fice! max temperature for cloud ice formation
75 real(r8), private:: tmin_fice! min temperature for cloud ice formation
76 real(r8), private:: tmax_fsnow ! max temperature for transition to convective snow
77 real(r8), private:: tmin_fsnow ! min temperature for transition to convective snow
78 real(r8), private:: rhonot ! air density at surface
79 real(r8), private:: t0 ! Freezing temperature
80 real(r8), private:: cldmin ! assumed minimum cloud amount
81 real(r8), private:: small ! small number compared to unity
82 real(r8), private:: c ! constant for graupel like snow cm**(1-d)/s
83 real(r8), private:: d ! constant for graupel like snow
84 real(r8), private:: esi ! collection efficient for ice by snow
85 real(r8), private:: esw ! collection efficient for water by snow
86 real(r8), private:: nos ! particles snow / cm**4
87 real(r8), private:: pi ! Mathematical constant
88 real(r8), private:: gravit ! Gravitational acceleration at surface
89 real(r8), private:: rh2o
90 real(r8), private:: prhonos
91 real(r8), private:: thrpd ! numerical three added to d
92 real(r8), private:: gam3pd ! gamma function on (3+d)
93 real(r8), private:: gam4pd ! gamma function on (4+d)
94 real(r8), private:: rhoi ! ice density
95 real(r8), private:: rhos ! snow density
96 real(r8), private:: rhow ! water density
97 real(r8), private:: mcon01 ! constants used in cloud microphysics
98 real(r8), private:: mcon02 ! constants used in cloud microphysics
99 real(r8), private:: mcon03 ! constants used in cloud microphysics
100 real(r8), private:: mcon04 ! constants used in cloud microphysics
101 real(r8), private:: mcon05 ! constants used in cloud microphysics
102 real(r8), private:: mcon06 ! constants used in cloud microphysics
103 real(r8), private:: mcon07 ! constants used in cloud microphysics
104 real(r8), private:: mcon08 ! constants used in cloud microphysics
106 integer, private :: k1mb ! index of the eta level near 1 mb
108 ! Parameters used in findmcnew
109 real(r8) :: capnsi ! sea ice cloud particles / cm3
110 real(r8) :: capnc ! cold and oceanic cloud particles / cm3
111 real(r8) :: capnw ! warm continental cloud particles / cm3
112 real(r8) :: kconst ! const for terminal velocity (stokes regime)
113 real(r8) :: effc ! collection efficiency
114 real(r8) :: alpha ! ratio of 3rd moment radius to 2nd
115 real(r8) :: capc ! constant for autoconversion
116 real(r8) :: convfw ! constant used for fall velocity calculation
117 real(r8) :: cracw ! constant used for rain accreting water
118 real(r8) :: critpr ! critical precip rate collection efficiency changes
119 real(r8) :: ciautb ! coefficient of autoconversion of ice (1/s)
121 #ifdef DEBUG
122 integer, private,parameter :: nlook = 1 ! Number of points to examine
123 integer, private :: ilook(nlook) ! Longitude index to examine
124 integer, private :: latlook(nlook) ! Latitude index to examine
125 integer, private :: lchnklook(nlook) ! Chunk index to examine
126 integer, private :: icollook(nlook) ! Column index to examine
127 #endif
128 ! Private data
129 real(r8), parameter :: unset_r8 = huge(1.0_r8)
131 contains
132 !===============================================================================
133 subroutine cldwat_readnl(nlfile)
134 #ifndef WRF_PORT
135 use namelist_utils, only: find_group_name
136 use units, only: getunit, freeunit
137 use mpishorthand
138 #endif
140 character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
141 #ifndef WRF_PORT
142 ! Namelist variables
143 real(r8) :: cldwat_icritw = unset_r8 ! icritw = threshold for autoconversion of warm ice
144 real(r8) :: cldwat_icritc = unset_r8 ! icritc = threshold for autoconversion of cold ice
145 real(r8) :: cldwat_conke = unset_r8 ! conke = tunable constant for evaporation of precip
146 real(r8) :: cldwat_r3lcrit = unset_r8 ! r3lcrit = critical radius where liq conversion begins
148 ! Local variables
149 integer :: unitn, ierr
150 character(len=*), parameter :: subname = 'cldwat_readnl'
152 namelist /cldwat_nl/ cldwat_icritw, cldwat_icritc, cldwat_conke, cldwat_r3lcrit
154 !-----------------------------------------------------------------------------
156 if (masterproc) then
157 unitn = getunit()
158 open( unitn, file=trim(nlfile), status='old' )
159 call find_group_name(unitn, 'cldwat_nl', status=ierr)
160 if (ierr == 0) then
161 read(unitn, cldwat_nl, iostat=ierr)
162 if (ierr /= 0) then
163 call endrun(subname // ':: ERROR reading namelist')
164 end if
165 end if
166 close(unitn)
167 call freeunit(unitn)
169 ! set local variables
170 icritw = cldwat_icritw
171 icritc = cldwat_icritc
172 conke = cldwat_conke
173 r3lcrit = cldwat_r3lcrit
175 end if
179 #ifdef SPMD
180 ! Broadcast namelist variables
181 call mpibcast(icritw, 1, mpir8, 0, mpicom)
182 call mpibcast(icritc, 1, mpir8, 0, mpicom)
183 call mpibcast(conke, 1, mpir8, 0, mpicom)
184 call mpibcast(r3lcrit, 1, mpir8, 0, mpicom)
185 #endif
186 #endif
187 end subroutine cldwat_readnl
189 !================================================================================================
190 subroutine cldwat_fice(ncol, t, fice, fsnow)
191 !
192 ! Compute the fraction of the total cloud water which is in ice phase.
193 ! The fraction depends on temperature only.
194 ! This is the form that was used for radiation, the code came from cldefr originally
195 !
196 ! Author: B. A. Boville Sept 10, 2002
197 ! modified: PJR 3/13/03 (added fsnow to ascribe snow production for convection )
198 !-----------------------------------------------------------------------
199 use physconst, only: tmelt
200 implicit none
202 ! Arguments
203 integer, intent(in) :: ncol ! number of active columns
204 real(r8), intent(in) :: t(pcols,pver) ! temperature
206 real(r8), intent(out) :: fice(pcols,pver) ! Fractional ice content within cloud
207 real(r8), intent(out) :: fsnow(pcols,pver) ! Fractional snow content for convection
209 ! Local variables
210 integer :: i,k ! loop indexes
212 !-----------------------------------------------------------------------
214 !BSINGH - This is a temporary fix for uninitialized tmax_* and tmin_* variables
215 !This can be fixed by calling the inimc subroutine in physics init calls
216 tmax_fice = tmelt - 10._r8
217 tmin_fice = tmax_fice - 30._r8
218 tmax_fsnow = tmelt
219 tmin_fsnow = tmelt - 5._r8
222 ! Define fractional amount of cloud that is ice
223 do k=1,pver
224 do i=1,ncol
226 ! If warmer than tmax then water phase
227 if (t(i,k) > tmax_fice) then
228 fice(i,k) = 0.0_r8
230 ! If colder than tmin then ice phase
231 else if (t(i,k) < tmin_fice) then
232 fice(i,k) = 1.0_r8
234 ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
235 else
236 fice(i,k) =(tmax_fice - t(i,k)) / (tmax_fice - tmin_fice)
237 end if
239 ! snow fraction partitioning
241 ! If warmer than tmax then water phase
242 if (t(i,k) > tmax_fsnow) then
243 fsnow(i,k) = 0.0_r8
245 ! If colder than tmin then ice phase
246 else if (t(i,k) < tmin_fsnow) then
247 fsnow(i,k) = 1.0_r8
249 ! Otherwise mixed phase, with ice fraction decreasing linearly from tmin to tmax
250 else
251 fsnow(i,k) =(tmax_fsnow - t(i,k)) / (tmax_fsnow - tmin_fsnow)
252 end if
254 end do
255 end do
257 return
258 end subroutine cldwat_fice
259 #ifndef WRF_PORT
260 subroutine inimc( tmeltx, rhonotx, gravitx, rh2ox)
261 !-----------------------------------------------------------------------
262 !
263 ! Purpose:
264 ! initialize constants for the prognostic condensate
265 !
266 ! Author: P. Rasch, April 1997
267 !
268 !-----------------------------------------------------------------------
269 use pmgrid, only: plev, plevp
270 use dycore, only: dycore_is, get_resolution
271 use hycoef, only: hypm
272 use phys_control, only: phys_getopts
274 integer k
275 real(r8), intent(in) :: tmeltx
276 real(r8), intent(in) :: rhonotx
277 real(r8), intent(in) :: gravitx
278 real(r8), intent(in) :: rh2ox
280 #ifdef UNICOSMP
281 real(r8) signgam ! variable required by cray gamma function
282 external gamma
283 #endif
284 character(len=16) :: microp_scheme
286 ! Get microphysics option
288 call phys_getopts( microp_scheme_out = microp_scheme )
290 ! Set following for all physics packages
292 tmax_fice = tmeltx - 10._r8
293 !! tmax_fice = tmeltx
294 !! tmin_fice = tmax_fice - 20.
295 tmin_fice = tmax_fice - 30._r8
296 tmax_fsnow = tmeltx
297 tmin_fsnow = tmeltx - 5._r8
299 ! Set remaining for RK microphysics
301 if( microp_scheme .eq. 'RK' ) then
302 rhonot = rhonotx ! air density at surface (gm/cm3)
303 gravit = gravitx
304 rh2o = rh2ox
305 rhos = .1_r8 ! assumed snow density (gm/cm3)
306 rhow = 1._r8 ! water density
307 rhoi = 1._r8 ! ice density
308 esi = 1.0_r8 ! collection efficient for ice by snow
309 esw = 0.1_r8 ! collection efficient for water by snow
310 t0 = tmeltx ! approximate freezing temp
311 cldmin = 0.02_r8 ! assumed minimum cloud amount
312 small = 1.e-22_r8 ! a small number compared to unity
313 c = 152.93_r8 ! constant for graupel like snow cm**(1-d)/s
314 d = 0.25_r8 ! constant for graupel like snow
315 nos = 3.e-2_r8 ! particles snow / cm**4
316 pi = 4._r8*atan(1.0_r8)
317 prhonos = pi*rhos*nos
318 thrpd = 3._r8 + d
319 if (d==0.25_r8) then
320 gam3pd = 2.549256966718531_r8 ! only right for d = 0.25
321 gam4pd = 8.285085141835282_r8
322 else
323 #ifdef UNICOSMP
324 call gamma(3._r8+d, signgam, gam3pd)
325 gam3pd = sign(exp(gam3pd),signgam)
326 call gamma(4._r8+d, signgam, gam4pd)
327 gam4pd = sign(exp(gam4pd),signgam)
328 write(iulog,*) ' d, gamma(3+d), gamma(4+d) =', gam3pd, gam4pd
329 #ifdef WRF_PORT
330 call wrf_message(iulog)
331 #endif
332 #else
333 write(iulog,*) ' can only use d ne 0.25 on a cray '
334 #ifdef WRF_PORT
335 call wrf_message(iulog)
336 #endif
337 stop
338 #endif
339 endif
340 mcon01 = pi*nos*c*gam3pd/4._r8
341 mcon02 = 1._r8/(c*gam4pd*sqrt(rhonot)/(6*prhonos**(d/4._r8)))
342 mcon03 = -(0.5_r8+d/4._r8)
343 mcon04 = 4._r8/(4._r8+d)
344 mcon05 = (3+d)/(4+d)
345 mcon06 = (3+d)/4._r8
346 mcon07 = mcon01*sqrt(rhonot)*mcon02**mcon05/prhonos**mcon06
347 mcon08 = -0.5_r8/(4._r8+d)
349 ! find the level about 1mb, we wont do the microphysics above this level
350 k1mb = 1
351 do k=1,pver-1
352 if (hypm(k) < 1.e2_r8 .and. hypm(k+1) >= 1.e2_r8) then
353 if (1.e2_r8-hypm(k) < hypm(k+1)-1.e2_r8) then
354 k1mb = k
355 else
356 k1mb = k + 1
357 end if
358 goto 20
359 end if
360 end do
361 if (masterproc) then
362 write(iulog,*)'inimc: model levels bracketing 1 mb not found'
363 #ifdef WRF_PORT
364 call wrf_message(iulog)
365 #endif
366 end if
367 ! call endrun
368 k1mb = 1
369 20 if( masterproc ) write(iulog,*)'inimc: model level nearest 1 mb is',k1mb,'which is',hypm(k1mb),'pascals'
370 #ifdef WRF_PORT
371 call wrf_message(iulog)
372 #endif
374 if( masterproc ) write(iulog,*) 'cloud water initialization by inimc complete '
375 #ifdef WRF_PORT
376 call wrf_message(iulog)
377 #endif
379 ! Initialize parameters used by findmcnew
380 capnw = 400._r8 ! warm continental cloud particles / cm3
381 capnc = 150._r8 ! cold and oceanic cloud particles / cm3
382 ! capnsi = 40._r8 ! sea ice cloud particles density / cm3
383 capnsi = 75._r8 ! sea ice cloud particles density / cm3
385 kconst = 1.18e6_r8 ! const for terminal velocity
387 ! effc = 1._r8 ! autoconv collection efficiency following boucher 96
388 ! effc = .55*0.05_r8 ! autoconv collection efficiency following baker 93
389 effc = 0.55_r8 ! autoconv collection efficiency following tripoli and cotton
390 ! effc = 0._r8 ! turn off warm-cloud autoconv
391 alpha = 1.1_r8**4
392 capc = pi**(-.333_r8)*kconst*effc *(0.75_r8)**(1.333_r8)*alpha ! constant for autoconversion
394 ! critical precip rate at which we assume the collector drops can change the
395 ! drop size enough to enhance the auto-conversion process (mm/day)
396 critpr = 0.5_r8
397 convfw = 1.94_r8*2.13_r8*sqrt(rhow*1000._r8*9.81_r8*2.7e-4_r8)
399 ! liquid microphysics
400 ! cracw = 6_r8 ! beheng
401 cracw = .884_r8*sqrt(9.81_r8/(rhow*1000._r8*2.7e-4_r8)) ! tripoli and cotton
403 ! ice microphysics
404 ciautb = 5.e-4_r8
406 if ( masterproc ) then
407 write(iulog,*)'tuning parameters cldwat: icritw',icritw,'icritc',icritc,'conke',conke,'r3lcrit',r3lcrit
408 #ifdef WRF_PORT
409 call wrf_message(iulog)
410 #endif
411 write(iulog,*)'tuning parameters cldwat: capnw',capnw,'capnc',capnc,'capnsi',capnsi,'kconst',kconst
412 #ifdef WRF_PORT
413 call wrf_message(iulog)
414 #endif
415 write(iulog,*)'tuning parameters cldwat: effc',effc,'alpha',alpha,'capc',capc
416 #ifdef WRF_PORT
417 call wrf_message(iulog)
418 #endif
419 write(iulog,*)'tuning parameters cldwat: critpr',critpr,'convfw',convfw,'cracw',cracw,'ciautb',ciautb
420 #ifdef WRF_PORT
421 call wrf_message(iulog)
422 #endif
423 endif
424 endif
426 return
427 end subroutine inimc
429 subroutine pcond (lchnk ,ncol , &
430 tn ,ttend ,qn ,qtend ,omega , &
431 cwat ,p ,pdel ,cldn ,fice , fsnow, &
432 cme ,prodprec,prodsnow,evapprec,evapsnow,evapheat, prfzheat, &
433 meltheat,precip ,snowab ,deltat ,fwaut , &
434 fsaut ,fracw ,fsacw ,fsaci ,lctend , &
435 rhdfda ,rhu00 ,seaicef, zi ,ice2pr, liq2pr, &
436 liq2snow, snowh, rkflxprc, rkflxsnw, pracwo, psacwo, psacio)
438 !-----------------------------------------------------------------------
439 !
440 ! Purpose:
441 ! The public interface to the cloud water parameterization
442 ! returns tendencies to water vapor, temperature and cloud water variables
443 !
444 ! For basic method
445 ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
446 ! model climate using diagnosed and
447 ! predicted condensate parameterizations, 1998, J. Clim., 11,
448 ! pp1587---1614.
449 !
450 ! For important modifications to improve the method of determining
451 ! condensation/evaporation see Zhang et al (2001, in preparation)
452 !
453 ! Authors: M. Zhang, W. Lin, P. Rasch and J.E. Kristjansson
454 ! B. A. Boville (latent heat of fusion)
455 !-----------------------------------------------------------------------
456 use wv_saturation, only: vqsatd
457 use cam_control_mod, only: nlvdry
458 !
459 !---------------------------------------------------------------------
460 !
461 ! Input Arguments
462 !
463 integer, intent(in) :: lchnk ! chunk identifier
464 integer, intent(in) :: ncol ! number of atmospheric columns
466 real(r8), intent(in) :: fice(pcols,pver) ! fraction of cwat that is ice
467 real(r8), intent(in) :: fsnow(pcols,pver) ! fraction of rain that freezes to snow
468 real(r8), intent(in) :: cldn(pcols,pver) ! new value of cloud fraction (fraction)
469 real(r8), intent(in) :: cwat(pcols,pver) ! cloud water (kg/kg)
470 real(r8), intent(in) :: omega(pcols,pver) ! vert pressure vel (Pa/s)
471 real(r8), intent(in) :: p(pcols,pver) ! pressure (K)
472 real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness (Pa)
473 real(r8), intent(in) :: qn(pcols,pver) ! new water vapor (kg/kg)
474 real(r8), intent(in) :: qtend(pcols,pver) ! mixing ratio tend (kg/kg/s)
475 real(r8), intent(in) :: tn(pcols,pver) ! new temperature (K)
476 real(r8), intent(in) :: ttend(pcols,pver) ! temp tendencies (K/s)
477 real(r8), intent(in) :: deltat ! time step to advance solution over
478 real(r8), intent(in) :: lctend(pcols,pver) ! cloud liquid water tendencies ====wlin
479 real(r8), intent(in) :: rhdfda(pcols,pver) ! dG(a)/da, rh=G(a), when rh>u00 ====wlin
480 real(r8), intent(in) :: rhu00 (pcols,pver) ! Rhlim for cloud ====wlin
481 real(r8), intent(in) :: seaicef(pcols) ! sea ice fraction (fraction)
482 real(r8), intent(in) :: zi(pcols,pverp) ! layer interfaces (m)
483 real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m)
484 !
485 ! Output Arguments
486 !
487 real(r8), intent(out) :: cme (pcols,pver) ! rate of cond-evap of condensate (1/s)
488 real(r8), intent(out) :: prodprec(pcols,pver) ! rate of conversion of condensate to precip (1/s)
489 real(r8), intent(out) :: evapprec(pcols,pver) ! rate of evaporation of falling precip (1/s)
490 real(r8), intent(out) :: evapsnow(pcols,pver) ! rate of evaporation of falling snow (1/s)
491 real(r8), intent(out) :: evapheat(pcols,pver) ! heating rate due to evaporation of precip (W/kg)
492 real(r8), intent(out) :: prfzheat(pcols,pver) ! heating rate due to freezing of precip (W/kg)
493 real(r8), intent(out) :: meltheat(pcols,pver) ! heating rate due to snow melt (W/kg)
494 real(r8), intent(out) :: precip(pcols) ! rate of precipitation (kg / (m**2 * s))
495 real(r8), intent(out) :: snowab(pcols) ! rate of snow (kg / (m**2 * s))
496 real(r8), intent(out) :: ice2pr(pcols,pver) ! rate of conversion of ice to precip
497 real(r8), intent(out) :: liq2pr(pcols,pver) ! rate of conversion of liquid to precip
498 real(r8), intent(out) :: liq2snow(pcols,pver) ! rate of conversion of liquid to snow
499 real(r8), intent(out) :: rkflxprc(pcols,pverp) ! grid-box mean RK flux_large_scale_cloud_rain+snow at interfaces (kg m^-2 s^-1)
500 real(r8), intent(out) :: rkflxsnw(pcols,pverp) ! grid-box mean RK flux_large_scale_cloud_snow at interfaces (kg m^-2 s^-1)
501 ! intent(out)s here for pcond to pass to stratiform.F90 to be addflded/outflded
502 real(r8), intent(out) :: pracwo(pcols,pver) ! accretion of cloud water by rain (1/s)
503 real(r8), intent(out) :: psacwo(pcols,pver) ! accretion of cloud water by snow (1/s)
504 real(r8), intent(out) :: psacio(pcols,pver) ! accretion of cloud ice by snow (1/s)
506 real(r8) nice2pr ! rate of conversion of ice to snow
507 real(r8) nliq2pr ! rate of conversion of liquid to precip
508 real(r8) nliq2snow ! rate of conversion of liquid to snow
509 real(r8) prodsnow(pcols,pver) ! rate of production of snow
511 !
512 ! Local workspace
513 !
514 real(r8) :: precab(pcols) ! rate of precipitation (kg / (m**2 * s))
515 integer i ! work variable
516 integer iter ! #iterations for precipitation calculation
517 integer k ! work variable
518 integer l ! work variable
520 real(r8) cldm(pcols) ! mean cloud fraction over the time step
521 real(r8) cldmax(pcols) ! max cloud fraction above
522 real(r8) coef(pcols) ! conversion time scale for condensate to rain
523 real(r8) cwm(pcols) ! cwat mixing ratio at midpoint of time step
524 real(r8) cwn(pcols) ! cwat mixing ratio at end
525 real(r8) denom ! work variable
526 real(r8) dqsdt ! change in sat spec. hum. wrt temperature
527 real(r8) es(pcols) ! sat. vapor pressure
528 real(r8) fracw(pcols,pver) ! relative importance of collection of liquid by rain
529 real(r8) fsaci(pcols,pver) ! relative importance of collection of ice by snow
530 real(r8) fsacw(pcols,pver) ! relative importance of collection of liquid by snow
531 real(r8) fsaut(pcols,pver) ! relative importance of ice auto conversion
532 real(r8) fwaut(pcols,pver) ! relative importance of warm cloud autoconversion
533 real(r8) gamma(pcols) ! d qs / dT
534 real(r8) icwc(pcols) ! in-cloud water content (kg/kg)
535 real(r8) mincld ! a small cloud fraction to avoid / zero
536 real(r8) omeps ! 1 minus epsilon
537 real(r8),parameter ::omsm=0.99999_r8 ! a number just less than unity (for rounding)
538 real(r8) prprov(pcols) ! provisional value of precip at btm of layer
539 real(r8) prtmp ! work variable
540 real(r8) q(pcols,pver) ! mixing ratio before time step ignoring condensate
541 real(r8) qs(pcols) ! spec. hum. of water vapor
542 real(r8) qsn, esn ! work variable
543 real(r8) qsp(pcols,pver) ! sat pt mixing ratio
544 real(r8) qtl(pcols) ! tendency which would saturate the grid box in deltat
545 real(r8) qtmp, ttmp ! work variable
546 real(r8) relhum1(pcols) ! relative humidity
547 real(r8) relhum(pcols) ! relative humidity
548 !!$ real(r8) tc ! crit temp of transition to ice
549 real(r8) t(pcols,pver) ! temp before time step ignoring condensate
550 real(r8) tsp(pcols,pver) ! sat pt temperature
551 real(r8) pol ! work variable
552 real(r8) cdt ! work variable
553 real(r8) wtthick ! work variable
555 ! Extra local work space for cloud scheme modification
557 real(r8) cpohl !Cp/Hlatv
558 real(r8) hlocp !Hlatv/Cp
559 real(r8) dto2 !0.5*deltat (delta=2.0*dt)
560 real(r8) calpha(pcols) !alpha of new C - E scheme formulation
561 real(r8) cbeta (pcols) !beta of new C - E scheme formulation
562 real(r8) cbetah(pcols) !beta_hat at saturation portion
563 real(r8) cgamma(pcols) !gamma of new C - E scheme formulation
564 real(r8) cgamah(pcols) !gamma_hat at saturation portion
565 real(r8) rcgama(pcols) !gamma/gamma_hat
566 real(r8) csigma(pcols) !sigma of new C - E scheme formulation
567 real(r8) cmec1 (pcols) !c1 of new C - E scheme formulation
568 real(r8) cmec2 (pcols) !c2 of new C - E scheme formulation
569 real(r8) cmec3 (pcols) !c3 of new C - E scheme formulation
570 real(r8) cmec4 (pcols) !c4 of new C - E scheme formulation
571 real(r8) cmeres(pcols) !residual cond of over-sat after cme and evapprec
572 real(r8) ctmp !a scalar representation of cmeres
573 real(r8) clrh2o ! Ratio of latvap to water vapor gas const
574 real(r8) ice(pcols,pver) ! ice mixing ratio
575 real(r8) liq(pcols,pver) ! liquid mixing ratio
576 real(r8) rcwn(pcols,2,pver), rliq(pcols,2,pver), rice(pcols,2,pver)
577 real(r8) cwnsave(pcols,2,pver), cmesave(pcols,2,pver)
578 real(r8) prodprecsave(pcols,2,pver)
579 logical error_found
580 !
581 !------------------------------------------------------------
582 !
583 clrh2o = hlatv/rh2o ! Ratio of latvap to water vapor gas const
584 omeps = 1.0_r8 - epsqs
585 #ifdef PERGRO
586 mincld = 1.e-4_r8
587 iter = 1 ! number of times to iterate the precipitation calculation
588 #else
589 mincld = 1.e-4_r8
590 iter = 2
591 #endif
592 ! omsm = 0.99999
593 cpohl = cp/hlatv
594 hlocp = hlatv/cp
595 dto2=0.5_r8*deltat
596 !
597 ! Constant for computing rate of evaporation of precipitation:
598 !
599 !!$ conke = 1.e-5
600 !!$ conke = 1.e-6
601 !
602 ! initialize a few single level fields
603 !
604 do i = 1,ncol
605 precip(i) = 0.0_r8
606 precab(i) = 0.0_r8
607 snowab(i) = 0.0_r8
608 cldmax(i) = 0.0_r8
609 end do
610 !
611 ! initialize multi-level fields
612 !
613 do k = 1,pver
614 do i = 1,ncol
615 q(i,k) = qn(i,k)
616 t(i,k) = tn(i,k)
617 ! q(i,k)=qn(i,k)-qtend(i,k)*deltat
618 ! t(i,k)=tn(i,k)-ttend(i,k)*deltat
619 end do
620 end do
621 cme (:ncol,:) = 0._r8
622 evapprec(:ncol,:) = 0._r8
623 prodprec(:ncol,:) = 0._r8
624 evapsnow(:ncol,:) = 0._r8
625 prodsnow(:ncol,:) = 0._r8
626 evapheat(:ncol,:) = 0._r8
627 meltheat(:ncol,:) = 0._r8
628 prfzheat(:ncol,:) = 0._r8
629 ice2pr(:ncol,:) = 0._r8
630 liq2pr(:ncol,:) = 0._r8
631 liq2snow(:ncol,:) = 0._r8
632 fwaut(:ncol,:) = 0._r8
633 fsaut(:ncol,:) = 0._r8
634 fracw(:ncol,:) = 0._r8
635 fsacw(:ncol,:) = 0._r8
636 fsaci(:ncol,:) = 0._r8
637 rkflxprc(:ncol,:) = 0._r8
638 rkflxsnw(:ncol,:) = 0._r8
640 pracwo(:ncol,:) = 0._r8
641 psacwo(:ncol,:) = 0._r8
642 psacio(:ncol,:) = 0._r8
643 !
644 ! find the wet bulb temp and saturation value
645 ! for the provisional t and q without condensation
646 !
647 call findsp (lchnk, ncol, qn, tn, p, tsp, qsp)
648 do 800 k = k1mb,pver
649 call vqsatd (t(1,k), p(1,k), es, qs, gamma, ncol)
650 do i = 1,ncol
651 relhum(i) = q(i,k)/qs(i)
652 !
653 cldm(i) = max(cldn(i,k),mincld)
654 !
655 ! the max cloud fraction above this level
656 !
657 cldmax(i) = max(cldmax(i), cldm(i))
659 ! define the coefficients for C - E calculation
661 calpha(i) = 1.0_r8/qs(i)
662 cbeta (i) = q(i,k)/qs(i)**2*gamma(i)*cpohl
663 cbetah(i) = 1.0_r8/qs(i)*gamma(i)*cpohl
664 cgamma(i) = calpha(i)+hlatv*cbeta(i)/cp
665 cgamah(i) = calpha(i)+hlatv*cbetah(i)/cp
666 rcgama(i) = cgamma(i)/cgamah(i)
668 if(cldm(i) > mincld) then
669 icwc(i) = max(0._r8,cwat(i,k)/cldm(i))
670 else
671 icwc(i) = 0.0_r8
672 endif
673 !PJR the above logic give zero icwc with nonzero cwat, dont like it!
674 !PJR generates problems with csigma
675 !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
676 ! icwc(i) = max(1.e-8_r8,cwat(i,k)/cldm(i))
678 !
679 ! initial guess of evaporation, will be updated within iteration
680 !
681 evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
682 *(1._r8 - min(relhum(i),1._r8))
684 !
685 ! zero cmeres before iteration for each level
686 !
687 cmeres(i)=0.0_r8
689 end do
690 do i = 1,ncol
691 !
692 ! fractions of ice at this level
693 !
694 !!$ tc = t(i,k) - t0
695 !!$ fice(i,k) = max(0._r8,min(-tc*0.05,1.0_r8))
696 !
697 ! calculate the cooling due to a phase change of the rainwater
698 ! from above
699 !
700 if (t(i,k) >= t0) then
701 meltheat(i,k) = -hlatf * snowab(i) * gravit/pdel(i,k)
702 snowab(i) = 0._r8
703 else
704 meltheat(i,k) = 0._r8
705 endif
706 end do
708 !
709 ! calculate cme and formation of precip.
710 !
711 ! The cloud microphysics is highly nonlinear and coupled with cme
712 ! Both rain processes and cme are calculated iteratively.
713 !
714 do 100 l = 1,iter
716 do i = 1,ncol
718 !
719 ! calculation of cme has 4 scenarios
720 ! ==================================
721 !
722 if(relhum(i) > rhu00(i,k)) then
724 ! 1. whole grid saturation
725 ! ========================
726 if(relhum(i) >= 0.999_r8 .or. cldm(i) >= 0.999_r8 ) then
727 cme(i,k)=(calpha(i)*qtend(i,k)-cbetah(i)*ttend(i,k))/cgamah(i)
729 ! 2. fractional saturation
730 ! ========================
731 else
732 if (rhdfda(i,k) .eq. 0. .and. icwc(i).eq.0.) then
733 write (iulog,*) ' cldwat.F90: empty rh cloud ', i, k, lchnk
734 #ifdef WRF_PORT
735 call wrf_message(iulog)
736 #endif
737 write (iulog,*) ' relhum, iter ', relhum(i), l, rhu00(i,k), cldm(i), cldn(i,k)
738 #ifdef WRF_PORT
739 call wrf_message(iulog)
740 #endif
741 call endrun ()
742 endif
743 csigma(i) = 1.0_r8/(rhdfda(i,k)+cgamma(i)*icwc(i))
744 cmec1(i) = (1.0_r8-cldm(i))*csigma(i)*rhdfda(i,k)
745 cmec2(i) = cldm(i)*calpha(i)/cgamah(i)+(1.0_r8-rcgama(i)*cldm(i))* &
746 csigma(i)*calpha(i)*icwc(i)
747 cmec3(i) = cldm(i)*cbetah(i)/cgamah(i) + &
748 (cbeta(i)-rcgama(i)*cldm(i)*cbetah(i))*csigma(i)*icwc(i)
749 cmec4(i) = csigma(i)*cgamma(i)*icwc(i)
751 ! Q=C-E=-C1*Al + C2*Aq - C3* At + C4*Er
753 cme(i,k) = -cmec1(i)*lctend(i,k) + cmec2(i)*qtend(i,k) &
754 -cmec3(i)*ttend(i,k) + cmec4(i)*evapprec(i,k)
755 endif
757 ! 3. when rh < rhu00, evaporate existing cloud water
758 ! ==================================================
759 else if(cwat(i,k) > 0.0_r8)then
760 ! liquid water should be evaporated but not to exceed
761 ! saturation point. if qn > qsp, not to evaporate cwat
762 cme(i,k)=-min(max(0._r8,qsp(i,k)-qn(i,k)),cwat(i,k))/deltat
764 ! 4. no condensation nor evaporation
765 ! ==================================
766 else
767 cme(i,k)=0.0_r8
768 endif
771 end do !end loop for cme update
773 ! Because of the finite time step,
774 ! place a bound here not to exceed wet bulb point
775 ! and not to evaporate more than available water
776 !
777 do i = 1, ncol
778 qtmp = qn(i,k) - cme(i,k)*deltat
780 ! possibilities to have qtmp > qsp
781 !
782 ! 1. if qn > qs(tn), it condenses;
783 ! if after applying cme, qtmp > qsp, more condensation is applied.
784 !
785 ! 2. if qn < qs, evaporation should not exceed qsp,
787 if(qtmp > qsp(i,k)) then
788 cme(i,k) = cme(i,k) + (qtmp-qsp(i,k))/deltat
789 endif
791 !
792 ! if net evaporation, it should not exceed available cwat
793 !
794 if(cme(i,k) < -cwat(i,k)/deltat) &
795 cme(i,k) = -cwat(i,k)/deltat
796 !
797 ! addition of residual condensation from previous step of iteration
798 !
799 cme(i,k) = cme(i,k) + cmeres(i)
801 end do
803 ! limit cme for roundoff errors
804 do i = 1, ncol
805 cme(i,k) = cme(i,k)*omsm
806 end do
808 do i = 1,ncol
809 !
810 ! as a safe limit, condensation should not reduce grid mean rh below rhu00
811 !
812 if(cme(i,k) > 0.0_r8 .and. relhum(i) > rhu00(i,k) ) &
813 cme(i,k) = min(cme(i,k), (qn(i,k)-qs(i)*rhu00(i,k))/deltat)
814 !
815 ! initial guess for cwm (mean cloud water over time step) if 1st iteration
816 !
817 if(l < 2) then
818 cwm(i) = max(cwat(i,k)+cme(i,k)*dto2, 0._r8)
819 endif
821 enddo
823 ! provisional precipitation falling through model layer
824 do i = 1,ncol
825 !!$ prprov(i) = precab(i) + prodprec(i,k)*pdel(i,k)/gravit
826 ! rain produced in this layer not too effective in collection process
827 wtthick = max(0._r8,min(0.5_r8,((zi(i,k)-zi(i,k+1))/1000._r8)**2))
828 prprov(i) = precab(i) + wtthick*prodprec(i,k)*pdel(i,k)/gravit
829 end do
831 ! calculate conversion of condensate to precipitation by cloud microphysics
832 call findmcnew (lchnk ,ncol , &
833 k ,prprov ,snowab, t ,p , &
834 cwm ,cldm ,cldmax ,fice(1,k),coef , &
835 fwaut(1,k),fsaut(1,k),fracw(1,k),fsacw(1,k),fsaci(1,k), &
836 seaicef, snowh, pracwo(1,k), psacwo(1,k), psacio(1,k))
838 !
839 ! calculate the precip rate
840 !
841 error_found = .false.
842 do i = 1,ncol
843 if (cldm(i) > 0) then
844 !
845 ! first predict the cloud water
846 !
847 cdt = coef(i)*deltat
848 if (cdt > 0.01_r8) then
849 pol = cme(i,k)/coef(i) ! production over loss
850 cwn(i) = max(0._r8,(cwat(i,k)-pol)*exp(-cdt)+ pol)
851 else
852 cwn(i) = max(0._r8,(cwat(i,k) + cme(i,k)*deltat)/(1+cdt))
853 endif
854 !
855 ! now back out the tendency of net rain production
856 !
857 prodprec(i,k) = max(0._r8,cme(i,k)-(cwn(i)-cwat(i,k))/deltat)
858 else
859 prodprec(i,k) = 0.0_r8
860 cwn(i) = 0._r8
861 endif
863 ! provisional calculation of conversion terms
864 ice2pr(i,k) = prodprec(i,k)*(fsaut(i,k)+fsaci(i,k))
865 liq2pr(i,k) = prodprec(i,k)*(fwaut(i,k)+fsacw(i,k)+fracw(i,k))
866 !old liq2snow(i,k) = prodprec(i,k)*fsacw(i,k)
868 ! revision suggested by Jim McCaa
869 ! it controls the amount of snow hitting the sfc
870 ! by forcing a lot of conversion of cloud liquid to snow phase
871 ! it might be better done later by an explicit representation of
872 ! rain accreting ice (and freezing), or by an explicit freezing of raindrops
873 liq2snow(i,k) = max(prodprec(i,k)*fsacw(i,k), fsnow(i,k)*liq2pr(i,k))
875 ! bounds
876 nice2pr = min(ice2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*fice(i,k)/deltat)
877 nliq2pr = min(liq2pr(i,k),(cwat(i,k)+cme(i,k)*deltat)*(1._r8-fice(i,k))/deltat)
878 ! write(iulog,*) ' prodprec ', i, k, prodprec(i,k)
879 ! write(iulog,*) ' nliq2pr, nice2pr ', nliq2pr, nice2pr
880 if (liq2pr(i,k).ne.0._r8) then
881 nliq2snow = liq2snow(i,k)*nliq2pr/liq2pr(i,k) ! correction
882 else
883 nliq2snow = liq2snow(i,k)
884 endif
886 ! avoid roundoff problems generating negatives
887 nliq2snow = nliq2snow*omsm
888 nliq2pr = nliq2pr*omsm
889 nice2pr = nice2pr*omsm
891 ! final estimates of conversion to precip and snow
892 prodprec(i,k) = (nliq2pr + nice2pr)
893 prodsnow(i,k) = (nice2pr + nliq2snow)
895 rcwn(i,l,k) = cwat(i,k) + (cme(i,k)- prodprec(i,k))*deltat
896 rliq(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)*(1._r8-fice(i,k)) - nliq2pr * deltat
897 rice(i,l,k) = (cwat(i,k) + cme(i,k)*deltat)* fice(i,k) - nice2pr *deltat
899 ! Save for sanity check later...
900 ! Putting sanity checks inside loops 100 and 800 screws up the
901 ! IBM compiler for reasons as yet unknown. TBH
902 cwnsave(i,l,k) = cwn(i)
903 cmesave(i,l,k) = cme(i,k)
904 prodprecsave(i,l,k) = prodprec(i,k)
905 ! End of save for sanity check later...
907 ! final version of condensate to precip terms
908 liq2pr(i,k) = nliq2pr
909 liq2snow(i,k) = nliq2snow
910 ice2pr(i,k) = nice2pr
912 cwn(i) = rcwn(i,l,k)
913 !
914 ! update any remaining provisional values
915 !
916 cwm(i) = (cwn(i) + cwat(i,k))*0.5_r8
917 !
918 ! update in cloud water
919 !
920 if(cldm(i) > mincld) then
921 icwc(i) = cwm(i)/cldm(i)
922 else
923 icwc(i) = 0.0_r8
924 endif
925 !PJR the above logic give zero icwc with nonzero cwat, dont like it!
926 !PJR generates problems with csigma
927 !PJR set the icwc to a very small number, so we can start from zero cloud cover and make some clouds
928 ! icwc(i) = max(1.e-8_r8,cwm(i)/cldm(i))
930 end do ! end of do i = 1,ncol
932 !
933 ! calculate provisional value of cloud water for
934 ! evaporation of precipitate (evapprec) calculation
935 !
936 do i = 1,ncol
937 qtmp = qn(i,k) - cme(i,k)*deltat
938 ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k) &
939 + (hlatv + hlatf*fice(i,k)) * cme(i,k) )
940 esn = estblf(ttmp)
941 qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
942 qtl(i) = max((qsn - qtmp)/deltat,0._r8)
943 relhum1(i) = qtmp/qsn
944 end do
945 !
946 do i = 1,ncol
947 #ifdef PERGRO
948 evapprec(i,k) = conke*(1._r8 - max(cldm(i),mincld))* &
949 sqrt(precab(i))*(1._r8 - min(relhum1(i),1._r8))
950 #else
951 evapprec(i,k) = conke*(1._r8 - cldm(i))*sqrt(precab(i)) &
952 *(1._r8 - min(relhum1(i),1._r8))
953 #endif
954 !
955 ! limit the evaporation to the amount which is entering the box
956 ! or saturates the box
957 !
958 prtmp = precab(i)*gravit/pdel(i,k)
959 evapprec(i,k) = min(evapprec(i,k), prtmp, qtl(i))*omsm
960 #ifdef PERGRO
961 ! zeroing needed for pert growth
962 evapprec(i,k) = 0._r8
963 #endif
964 !
965 ! Partition evaporation of precipitate between rain and snow using
966 ! the fraction of snow falling into the box. Determine the heating
967 ! due to evaporation. Note that evaporation is positive (loss of precip,
968 ! gain of vapor) and that heating is negative.
969 if (evapprec(i,k) > 0._r8) then
970 evapsnow(i,k) = evapprec(i,k) * snowab(i) / precab(i)
971 evapheat(i,k) = -hlatv * evapprec(i,k) - hlatf * evapsnow(i,k)
972 else
973 evapsnow(i,k) = 0._r8
974 evapheat(i,k) = 0._r8
975 end if
976 ! Account for the latent heat of fusion for liquid drops collected by falling snow
977 prfzheat(i,k) = hlatf * liq2snow(i,k)
978 end do
980 ! now remove the residual of any over-saturation. Normally,
981 ! the oversaturated water vapor should have been removed by
982 ! cme formulation plus constraints by wet bulb tsp/qsp
983 ! as computed above. However, because of non-linearity,
984 ! addition of (cme-evapprec) to update t and q may still cause
985 ! a very small amount of over saturation. It is called a
986 ! residual of over-saturation because theoretically, cme
987 ! should have taken care of all of large scale condensation.
988 !
990 do i = 1,ncol
991 qtmp = qn(i,k)-(cme(i,k)-evapprec(i,k))*deltat
992 ttmp = tn(i,k) + deltat/cp * ( meltheat(i,k) + evapheat(i,k) + prfzheat(i,k) &
993 + (hlatv + hlatf*fice(i,k)) * cme(i,k) )
994 esn = estblf(ttmp)
995 qsn = min(epsqs*esn/(p(i,k) - omeps*esn),1._r8)
996 !
997 !Upper stratosphere and mesosphere, qsn calculated
998 !above may be negative. Here just to skip it instead
999 !of resetting it to 1 as in aqsat
1000 !
1001 if(qtmp > qsn .and. qsn > 0) then
1002 !calculate dqsdt, a more precise calculation
1003 !which taking into account different range of T
1004 !can be found in aqsatd.F. Here follows
1005 !cond.F to calculate it.
1006 !
1007 denom = (p(i,k)-omeps*esn)*ttmp*ttmp
1008 dqsdt = clrh2o*qsn*p(i,k)/denom
1009 !
1010 !now extra condensation to bring air to just saturation
1011 !
1012 ctmp = (qtmp-qsn)/(1._r8+hlocp*dqsdt)/deltat
1013 cme(i,k) = cme(i,k)+ctmp
1014 !
1015 ! save residual on cmeres to addtion to cme on entering next iteration
1016 ! cme exit here contain the residual but overrided if back to iteration
1017 !
1018 cmeres(i) = ctmp
1019 else
1020 cmeres(i) = 0.0_r8
1021 endif
1022 end do
1024 100 continue ! end of do l = 1,iter
1026 !
1027 ! precipitation
1028 !
1029 do i = 1,ncol
1030 precip(i) = precip(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
1031 precab(i) = precab(i) + pdel(i,k)/gravit * (prodprec(i,k) - evapprec(i,k))
1032 if(precab(i).lt.0._r8) precab(i)=0._r8
1033 ! snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodprec(i,k)*fice(i,k) - evapsnow(i,k))
1034 snowab(i) = snowab(i) + pdel(i,k)/gravit * (prodsnow(i,k) - evapsnow(i,k))
1036 ! If temperature above freezing, all precip is rain flux. if temperature below freezing, all precip is snow flux.
1037 rkflxprc(i,k+1) = precab(i) !! making this consistent with other precip fluxes. prc = rain + snow
1038 !!rkflxprc(i,k+1) = precab(i) - snowab(i)
1039 rkflxsnw(i,k+1) = snowab(i)
1041 !!$ if ((precab(i)) < 1.e-10) then
1042 !!$ precab(i) = 0.
1043 !!$ snowab(i) = 0.
1044 !!$ endif
1045 end do
1046 800 continue ! level loop (k=1,pver)
1048 ! begin sanity checks
1049 error_found = .false.
1050 do k = k1mb,pver
1051 do l = 1,iter
1052 do i = 1,ncol
1053 if (abs(rcwn(i,l,k)).lt.1.e-300_r8) rcwn(i,l,k) = 0._r8
1054 if (abs(rliq(i,l,k)).lt.1.e-300_r8) rliq(i,l,k) = 0._r8
1055 if (abs(rice(i,l,k)).lt.1.e-300_r8) rice(i,l,k) = 0._r8
1056 if (rcwn(i,l,k).lt.0._r8) error_found = .true.
1057 if (rliq(i,l,k).lt.0._r8) error_found = .true.
1058 if (rice(i,l,k).lt.0._r8) error_found = .true.
1059 enddo
1060 enddo
1061 enddo
1062 if (error_found) then
1063 do k = k1mb,pver
1064 do l = 1,iter
1065 do i = 1,ncol
1066 if (rcwn(i,l,k).lt.0._r8) then
1067 write(iulog,*) ' prob with neg rcwn1 ', rcwn(i,l,k), &
1068 cwnsave(i,l,k)
1069 #ifdef WRF_PORT
1070 call wrf_message(iulog)
1071 #endif
1072 write(iulog,*) ' cwat, cme*deltat, prodprec*deltat ', &
1073 cwat(i,k), cmesave(i,l,k)*deltat, &
1074 prodprecsave(i,l,k)*deltat, &
1075 (cmesave(i,l,k)-prodprecsave(i,l,k))*deltat
1076 #ifdef WRF_PORT
1077 call wrf_message(iulog)
1078 #endif
1079 call endrun('PCOND')
1080 endif
1081 if (rliq(i,l,k).lt.0._r8) then
1082 write(iulog,*) ' prob with neg rliq1 ', rliq(i,l,k)
1083 #ifdef WRF_PORT
1084 call wrf_message(iulog)
1085 #endif
1086 call endrun('PCOND')
1087 endif
1088 if (rice(i,l,k).lt.0._r8) then
1089 write(iulog,*) ' prob with neg rice ', rice(i,l,k)
1090 #ifdef WRF_PORT
1091 call wrf_message(iulog)
1092 #endif
1093 call endrun('PCOND')
1094 endif
1095 enddo
1096 enddo
1097 enddo
1098 end if
1099 ! end sanity checks
1101 return
1102 end subroutine pcond
1104 !##############################################################################
1106 subroutine findmcnew (lchnk ,ncol , &
1107 k ,precab ,snowab, t ,p , &
1108 cwm ,cldm ,cldmax ,fice ,coef , &
1109 fwaut ,fsaut ,fracw ,fsacw ,fsaci , &
1110 seaicef ,snowh, pracwo, psacwo, psacio )
1112 !-----------------------------------------------------------------------
1113 !
1114 ! Purpose:
1115 ! calculate the conversion of condensate to precipitate
1116 !
1117 ! Method:
1118 ! See: Rasch, P. J, and J. E. Kristjansson, A Comparison of the CCM3
1119 ! model climate using diagnosed and
1120 ! predicted condensate parameterizations, 1998, J. Clim., 11,
1121 ! pp1587---1614.
1122 !
1123 ! Author: P. Rasch
1124 !
1125 !-----------------------------------------------------------------------
1126 use phys_grid, only: get_rlat_all_p
1127 use comsrf, only: landm
1128 !
1129 ! input args
1130 !
1131 integer, intent(in) :: lchnk ! chunk identifier
1132 integer, intent(in) :: ncol ! number of atmospheric columns
1133 integer, intent(in) :: k ! level index
1135 real(r8), intent(in) :: precab(pcols) ! rate of precipitation from above (kg / (m**2 * s))
1136 real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
1137 real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa)
1138 real(r8), intent(in) :: cldm(pcols) ! cloud fraction
1139 real(r8), intent(in) :: cldmax(pcols) ! max cloud fraction above this level
1140 real(r8), intent(in) :: cwm(pcols) ! condensate mixing ratio (kg/kg)
1141 real(r8), intent(in) :: fice(pcols) ! fraction of cwat that is ice
1142 real(r8), intent(in) :: seaicef(pcols) ! sea ice fraction
1143 real(r8), intent(in) :: snowab(pcols) ! rate of snow from above (kg / (m**2 * s))
1144 real(r8), intent(in) :: snowh(pcols) ! Snow depth over land, water equivalent (m)
1146 ! output arguments
1147 real(r8), intent(out) :: coef(pcols) ! conversion rate (1/s)
1148 real(r8), intent(out) :: fwaut(pcols) ! relative importance of liquid autoconversion (a diagnostic)
1149 real(r8), intent(out) :: fsaut(pcols) ! relative importance of ice autoconversion (a diagnostic)
1150 real(r8), intent(out) :: fracw(pcols) ! relative importance of rain accreting liquid (a diagnostic)
1151 real(r8), intent(out) :: fsacw(pcols) ! relative importance of snow accreting liquid (a diagnostic)
1152 real(r8), intent(out) :: fsaci(pcols) ! relative importance of snow accreting ice (a diagnostic)
1153 real(r8), intent(out) :: pracwo(pcols) ! accretion of cloud water by rain (1/s)
1154 real(r8), intent(out) :: psacwo(pcols) ! accretion of cloud water by snow (1/s)
1155 real(r8), intent(out) :: psacio(pcols) ! accretion of cloud ice by snow (1/s)
1158 ! work variables
1160 integer i
1161 integer ii
1162 integer ind(pcols)
1163 integer ncols
1165 real(r8), parameter :: degrad = 57.296_r8 ! divide by this to convert degrees to radians
1166 real(r8) capn ! local cloud particles / cm3
1167 real(r8) capnoice ! local cloud particles when not over sea ice / cm3
1168 real(r8) ciaut ! coefficient of autoconversion of ice (1/s)
1169 real(r8) cldloc(pcols) ! non-zero amount of cloud
1170 real(r8) cldpr(pcols) ! assumed cloudy volume occupied by rain and cloud
1171 real(r8) con1 ! work constant
1172 real(r8) con2 ! work constant
1173 real(r8) csacx ! constant used for snow accreting liquid or ice
1174 !!$ real(r8) dtice ! interval for transition from liquid to ice
1175 real(r8) icemr(pcols) ! in-cloud ice mixing ratio
1176 real(r8) icrit ! threshold for autoconversion of ice
1177 real(r8) liqmr(pcols) ! in-cloud liquid water mixing ratio
1178 real(r8) pracw ! rate of rain accreting water
1179 real(r8) prlloc(pcols) ! local rain flux in mm/day
1180 real(r8) prscgs(pcols) ! local snow amount in cgs units
1181 real(r8) psaci ! rate of collection of ice by snow (lin et al 1983)
1182 real(r8) psacw ! rate of collection of liquid by snow (lin et al 1983)
1183 real(r8) psaut ! rate of autoconversion of ice condensate
1184 real(r8) ptot ! total rate of conversion
1185 real(r8) pwaut ! rate of autoconversion of liquid condensate
1186 real(r8) r3l ! volume radius
1187 real(r8) rainmr(pcols) ! in-cloud rain mixing ratio
1188 real(r8) rat1 ! work constant
1189 real(r8) rat2 ! work constant
1190 !!$ real(r8) rdtice ! recipricol of dtice
1191 real(r8) rho(pcols) ! density (mks units)
1192 real(r8) rhocgs ! density (cgs units)
1193 real(r8) rlat(pcols) ! latitude (radians)
1194 real(r8) snowfr ! fraction of precipate existing as snow
1195 real(r8) totmr(pcols) ! in-cloud total condensate mixing ratio
1196 real(r8) vfallw ! fall speed of precipitate as liquid
1197 real(r8) wp ! weight factor used in calculating pressure dep of autoconversion
1198 real(r8) wsi ! weight factor for sea ice
1199 real(r8) wt ! fraction of ice
1200 real(r8) wland ! fraction of land
1202 ! real(r8) csaci
1203 ! real(r8) csacw
1204 ! real(r8) cwaut
1205 ! real(r8) efact
1206 ! real(r8) lamdas
1207 ! real(r8) lcrit
1208 ! real(r8) rcwm
1209 ! real(r8) r3lc2
1210 ! real(r8) snowmr(pcols)
1211 ! real(r8) vfalls
1213 real(8) ftot
1215 ! inline statement functions
1216 real(r8) heavy, heavym, a1, a2, heavyp, heavymp
1217 heavy(a1,a2) = max(0._r8,sign(1._r8,a1-a2)) ! heavyside function
1218 heavym(a1,a2) = max(0.01_r8,sign(1._r8,a1-a2)) ! modified heavyside function
1219 !
1220 ! New heavyside functions to perhaps address error growth problems
1221 !
1222 heavyp(a1,a2) = a1/(a2+a1+1.e-36_r8)
1223 heavymp(a1,a2) = (a1+0.01_r8*a2)/(a2+a1+1.e-36_r8)
1225 !
1226 ! find all the points where we need to do the microphysics
1227 ! and set the output variables to zero
1228 !
1229 ncols = 0
1230 do i = 1,ncol
1231 coef(i) = 0._r8
1232 fwaut(i) = 0._r8
1233 fsaut(i) = 0._r8
1234 fracw(i) = 0._r8
1235 fsacw(i) = 0._r8
1236 fsaci(i) = 0._r8
1237 liqmr(i) = 0._r8
1238 rainmr(i) = 0._r8
1239 if (cwm(i) > 1.e-20_r8) then
1240 ncols = ncols + 1
1241 ind(ncols) = i
1242 endif
1243 end do
1245 !cdir nodep
1246 !DIR$ CONCURRENT
1247 do ii = 1,ncols
1248 i = ind(ii)
1249 !
1250 ! the local cloudiness at this level
1251 !
1252 cldloc(i) = max(cldmin,cldm(i))
1253 !
1254 ! a weighted mean between max cloudiness above, and this layer
1255 !
1256 cldpr(i) = max(cldmin,(cldmax(i)+cldm(i))*0.5_r8)
1257 !
1258 ! decompose the suspended condensate into
1259 ! an incloud liquid and ice phase component
1260 !
1261 totmr(i) = cwm(i)/cldloc(i)
1262 icemr(i) = totmr(i)*fice(i)
1263 liqmr(i) = totmr(i)*(1._r8-fice(i))
1264 !
1265 ! density
1266 !
1267 rho(i) = p(i,k)/(287._r8*t(i,k))
1268 rhocgs = rho(i)*1.e-3_r8 ! density in cgs units
1269 !
1270 ! decompose the precipitate into a liquid and ice phase
1271 !
1272 if (t(i,k) > t0) then
1273 vfallw = convfw/sqrt(rho(i))
1274 rainmr(i) = precab(i)/(rho(i)*vfallw*cldpr(i))
1275 snowfr = 0
1276 ! snowmr(i)
1277 else
1278 snowfr = 1
1279 rainmr(i) = 0._r8
1280 endif
1281 ! rainmr(i) = (precab(i)-snowab(i))/(rho(i)*vfallw*cldpr(i))
1282 !
1283 ! local snow amount in cgs units
1284 !
1285 prscgs(i) = precab(i)/cldpr(i)*0.1_r8*snowfr
1286 ! prscgs(i) = snowab(i)/cldpr(i)*0.1
1287 !
1288 ! local rain amount in mm/day
1289 !
1290 prlloc(i) = precab(i)*86400._r8/cldpr(i)
1291 end do
1293 con1 = 1._r8/(1.333_r8*pi)**0.333_r8 * 0.01_r8 ! meters
1294 !
1295 ! calculate the conversion terms
1296 !
1297 call get_rlat_all_p(lchnk, ncol, rlat)
1299 !cdir nodep
1300 !DIR$ CONCURRENT
1301 do ii = 1,ncols
1302 i = ind(ii)
1303 rhocgs = rho(i)*1.e-3_r8 ! density in cgs units
1304 !
1305 ! exponential temperature factor
1306 !
1307 ! efact = exp(0.025*(t(i,k)-t0))
1308 !
1309 ! some temperature dependent constants
1310 !
1311 !!$ wt = min(1._r8,max(0._r8,(t0-t(i,k))*rdtice))
1312 wt = fice(i)
1313 icrit = icritc*wt + icritw*(1-wt)
1314 !
1315 ! jrm Reworked droplet number concentration algorithm
1316 ! Start with pressure-dependent value appropriate for continental air
1317 ! Note: reltab has a temperature dependence here
1318 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))))
1319 ! Modify for snow depth over land
1320 capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,snowh(i)*10._r8))
1321 ! Ramp between polluted value over land to clean value over ocean.
1322 capn = capn + (capnc-capn) * min(1.0_r8,max(0.0_r8,1.0_r8-landm(i,lchnk)))
1323 ! Ramp between the resultant value and a sea ice value in the presence of ice.
1324 capn = capn + (capnsi-capn) * min(1.0_r8,max(0.0_r8,seaicef(i)))
1325 ! end jrm
1326 !
1327 #ifdef DEBUG2
1328 if ( (lat(i) == latlook(1)) .or. (lat(i) == latlook(2)) ) then
1329 if (i == ilook(1)) then
1330 write(iulog,*) ' findmcnew: lat, k, seaicef, landm, wp, capnoice, capn ', &
1331 lat(i), k, seaicef(i), landm(i,lat(i)), wp, capnoice, capn
1332 #ifdef WRF_PORT
1333 call wrf_message(iulog)
1334 #endif
1335 endif
1336 endif
1337 #endif
1339 !
1340 ! useful terms in following calculations
1341 !
1342 rat1 = rhocgs/rhow
1343 rat2 = liqmr(i)/capn
1344 con2 = (rat1*rat2)**0.333_r8
1345 !
1346 ! volume radius
1347 !
1348 ! r3l = (rhocgs*liqmr(i)/(1.333*pi*capn*rhow))**0.333 * 0.01 ! meters
1349 r3l = con1*con2
1350 !
1351 ! critical threshold for autoconversion if modified for mixed phase
1352 ! clouds to mimic a bergeron findeisen process
1353 ! r3lc2 = r3lcrit*(1.-0.5*fice(i)*(1-fice(i)))
1354 !
1355 ! autoconversion of liquid
1356 !
1357 ! cwaut = 2.e-4
1358 ! cwaut = 1.e-3
1359 ! lcrit = 2.e-4
1360 ! lcrit = 5.e-4
1361 ! pwaut = max(0._r8,liqmr(i)-lcrit)*cwaut
1362 !
1363 ! pwaut is following tripoli and cotton (and many others)
1364 ! we reduce the autoconversion below critpr, because these are regions where
1365 ! the drop size distribution is likely to imply much smaller collector drops than
1366 ! those relevant for a cloud distribution corresponding to the value of effc = 0.55
1367 ! suggested by cotton (see austin 1995 JAS, baker 1993)
1369 ! easy to follow form
1370 ! pwaut = capc*liqmr(i)**2*rhocgs/rhow
1371 ! $ *(liqmr(i)*rhocgs/(rhow*capn))**(.333)
1372 ! $ *heavy(r3l,r3lcrit)
1373 ! $ *max(0.10_r8,min(1._r8,prlloc(i)/critpr))
1374 ! somewhat faster form
1375 #define HEAVYNEW
1376 #ifdef HEAVYNEW
1377 !#ifdef PERGRO
1378 pwaut = capc*liqmr(i)**2*rat1*con2*heavymp(r3l,r3lcrit) * &
1379 max(0.10_r8,min(1._r8,prlloc(i)/critpr))
1380 #else
1381 pwaut = capc*liqmr(i)**2*rat1*con2*heavym(r3l,r3lcrit)* &
1382 max(0.10_r8,min(1._r8,prlloc(i)/critpr))
1383 #endif
1384 !
1385 ! autoconversion of ice
1386 !
1387 ! ciaut = ciautb*efact
1388 ciaut = ciautb
1389 ! psaut = capc*totmr(i)**2*rhocgs/rhoi
1390 ! $ *(totmr(i)*rhocgs/(rhoi*capn))**(.333)
1391 !
1392 ! autoconversion of ice condensate
1393 !
1394 #ifdef PERGRO
1395 psaut = heavyp(icemr(i),icrit)*icemr(i)*ciaut
1396 #else
1397 psaut = max(0._r8,icemr(i)-icrit)*ciaut
1398 #endif
1399 !
1400 ! collection of liquid by rain
1401 !
1402 ! pracw = cracw*rho(i)*liqmr(i)*rainmr(i) !(beheng 1994)
1403 pracw = cracw*rho(i)*sqrt(rho(i))*liqmr(i)*rainmr(i) !(tripoli and cotton)
1405 pracwo(i)=pracw
1407 !! pracw = 0.
1408 !
1409 ! the following lines calculate the slope parameter and snow mixing ratio
1410 ! from the precip rate using the equations found in lin et al 83
1411 ! in the most natural form, but it is expensive, so after some tedious
1412 ! algebraic manipulation you can use the cheaper form found below
1413 ! vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs)
1414 ! $ *0.01 ! convert from cm/s to m/s
1415 ! snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i))
1416 ! snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04
1417 ! lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25
1418 ! csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd)
1419 !
1420 ! coefficient for collection by snow independent of phase
1421 !
1422 csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05
1424 !
1425 ! collection of liquid by snow (lin et al 1983)
1426 !
1427 psacw = csacx*liqmr(i)*esw
1428 #ifdef PERGRO
1429 ! this is necessary for pergro
1430 psacw = 0._r8
1431 #endif
1433 psacwo(i)=psacw
1435 !
1436 ! collection of ice by snow (lin et al 1983)
1437 !
1438 psaci = csacx*icemr(i)*esi
1439 !
1440 psacio(i)=psaci
1442 ! total conversion of condensate to precipitate
1443 !
1444 ptot = pwaut + psaut + pracw + psacw + psaci
1445 !
1446 ! the recipricol of cloud water amnt (or zero if no cloud water)
1447 !
1448 ! rcwm = totmr(i)/(max(totmr(i),small)**2)
1449 !
1450 ! turn the tendency back into a loss rate (1/seconds)
1451 !
1452 if (totmr(i) > 0._r8) then
1453 coef(i) = ptot/totmr(i)
1454 else
1455 coef(i) = 0._r8
1456 endif
1458 if (ptot.gt.0._r8) then
1459 fwaut(i) = pwaut/ptot
1460 fsaut(i) = psaut/ptot
1461 fracw(i) = pracw/ptot
1462 fsacw(i) = psacw/ptot
1463 fsaci(i) = psaci/ptot
1464 else
1465 fwaut(i) = 0._r8
1466 fsaut(i) = 0._r8
1467 fracw(i) = 0._r8
1468 fsacw(i) = 0._r8
1469 fsaci(i) = 0._r8
1470 endif
1472 ftot = fwaut(i)+fsaut(i)+fracw(i)+fsacw(i)+fsaci(i)
1473 ! if (abs(ftot-1._r8).gt.1.e-14_r8.and.ftot.ne.0._r8) then
1474 ! write(iulog,*) ' something is wrong in findmcnew ', ftot, &
1475 ! fwaut(i),fsaut(i),fracw(i),fsacw(i),fsaci(i)
1476 ! write(iulog,*) ' unscaled ', ptot, &
1477 ! pwaut,psaut,pracw,psacw,psaci
1478 ! write(iulog,*) ' totmr, liqmr, icemr ', totmr(i), liqmr(i), icemr(i)
1479 ! call endrun()
1480 ! endif
1481 end do
1482 #ifdef DEBUG
1483 i = icollook(nlook)
1484 if (lchnk == lchnklook(nlook) ) then
1485 write(iulog,*)
1486 write(iulog,*) '------', k, i, lchnk
1487 write(iulog,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4_r8
1488 write(iulog,*) ' frac: waut,saut,racw,sacw,saci ', &
1489 fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i)
1490 endif
1491 #endif
1493 return
1494 end subroutine findmcnew
1496 !##############################################################################
1498 subroutine findsp (lchnk, ncol, q, t, p, tsp, qsp)
1499 !-----------------------------------------------------------------------
1500 !
1501 ! Purpose:
1502 ! find the wet bulb temperature for a given t and q
1503 ! in a longitude height section
1504 ! wet bulb temp is the temperature and spec humidity that is
1505 ! just saturated and has the same enthalpy
1506 ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q
1507 ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q
1508 !
1509 ! Method:
1510 ! a Newton method is used
1511 ! first guess uses an algorithm provided by John Petch from the UKMO
1512 ! we exclude points where the physical situation is unrealistic
1513 ! e.g. where the temperature is outside the range of validity for the
1514 ! saturation vapor pressure, or where the water vapor pressure
1515 ! exceeds the ambient pressure, or the saturation specific humidity is
1516 ! unrealistic
1517 !
1518 ! Author: P. Rasch
1519 !
1520 !-----------------------------------------------------------------------
1521 !
1522 ! input arguments
1523 !
1524 integer, intent(in) :: lchnk ! chunk identifier
1525 integer, intent(in) :: ncol ! number of atmospheric columns
1527 real(r8), intent(in) :: q(pcols,pver) ! water vapor (kg/kg)
1528 real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
1529 real(r8), intent(in) :: p(pcols,pver) ! pressure (Pa)
1530 !
1531 ! output arguments
1532 !
1533 real(r8), intent(out) :: tsp(pcols,pver) ! saturation temp (K)
1534 real(r8), intent(out) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg)
1535 !
1536 ! local variables
1537 !
1538 integer i ! work variable
1539 integer k ! work variable
1540 logical lflg ! work variable
1541 integer iter ! work variable
1542 integer l ! work variable
1543 logical :: error_found
1545 real(r8) omeps ! 1 minus epsilon
1546 real(r8) trinv ! work variable
1547 real(r8) es ! sat. vapor pressure
1548 real(r8) desdt ! change in sat vap pressure wrt temperature
1549 ! real(r8) desdp ! change in sat vap pressure wrt pressure
1550 real(r8) dqsdt ! change in sat spec. hum. wrt temperature
1551 real(r8) dgdt ! work variable
1552 real(r8) g ! work variable
1553 real(r8) weight(pcols) ! work variable
1554 real(r8) hlatsb ! (sublimation)
1555 real(r8) hlatvp ! (vaporization)
1556 real(r8) hltalt(pcols,pver) ! lat. heat. of vap.
1557 real(r8) tterm ! work var.
1558 real(r8) qs ! spec. hum. of water vapor
1559 real(r8) tc ! crit temp of transition to ice
1561 ! work variables
1562 real(r8) t1, q1, dt, dq
1563 real(r8) dtm, dqm
1564 real(r8) qvd, a1, tmp
1565 real(r8) rair
1566 real(r8) r1b, c1, c2, c3
1567 real(r8) denom
1568 real(r8) dttol
1569 real(r8) dqtol
1570 integer doit(pcols)
1571 real(r8) enin(pcols), enout(pcols)
1572 real(r8) tlim(pcols)
1574 omeps = 1.0_r8 - epsqs
1575 trinv = 1.0_r8/ttrice
1576 a1 = 7.5_r8*log(10._r8)
1577 rair = 287.04_r8
1578 c3 = rair*a1/cp
1579 dtm = 0._r8 ! needed for iter=0 blowup with f90 -ei
1580 dqm = 0._r8 ! needed for iter=0 blowup with f90 -ei
1581 dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
1582 dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
1583 ! tmin = 173.16 ! the coldest temperature we can deal with
1584 !
1585 ! max number of times to iterate the calculation
1586 iter = 8
1587 !
1588 do k = k1mb,pver
1590 !
1591 ! first guess on the wet bulb temperature
1592 !
1593 do i = 1,ncol
1595 #ifdef DEBUG
1596 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
1597 write(iulog,*) ' '
1598 #ifdef WRF_PORT
1599 call wrf_message(iulog)
1600 #endif
1601 write(iulog,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
1602 #ifdef WRF_PORT
1603 call wrf_message(iulog)
1604 #endif
1605 endif
1606 #endif
1607 ! limit the temperature range to that relevant to the sat vap pres tables
1608 #if ( ! defined WACCM_PHYS )
1609 tlim(i) = min(max(t(i,k),173._r8),373._r8)
1610 #else
1611 tlim(i) = min(max(t(i,k),128._r8),373._r8)
1612 #endif
1613 es = estblf(tlim(i))
1614 denom = p(i,k) - omeps*es
1615 qs = epsqs*es/denom
1616 doit(i) = 0
1617 enout(i) = 1._r8
1618 ! make sure a meaningful calculation is possible
1619 if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
1620 !
1621 ! Saturation specific humidity
1622 !
1623 qs = min(epsqs*es/denom,1._r8)
1624 !
1625 ! "generalized" analytic expression for t derivative of es
1626 ! accurate to within 1 percent for 173.16 < t < 373.16
1627 !
1628 ! Weighting of hlat accounts for transition from water to ice
1629 ! polynomial expression approximates difference between es over
1630 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1631 ! -40): required for accurate estimate of es derivative in transition
1632 ! range from ice to water also accounting for change of hlatv with t
1633 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1634 !
1635 tc = tlim(i) - t0
1636 lflg = (tc >= -ttrice .and. tc < 0.0_r8)
1637 weight(i) = min(-tc*trinv,1.0_r8)
1638 hlatsb = hlatv + weight(i)*hlatf
1639 hlatvp = hlatv - 2369.0_r8*tc
1640 if (tlim(i) < t0) then
1641 hltalt(i,k) = hlatsb
1642 else
1643 hltalt(i,k) = hlatvp
1644 end if
1645 enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)
1647 ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
1648 tmp = q(i,k) - qs
1649 c1 = hltalt(i,k)*c3
1650 c2 = (tlim(i) + 36._r8)**2
1651 r1b = c2/(c2 + c1*qs)
1652 qvd = r1b*tmp
1653 tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
1654 #ifdef DEBUG
1655 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
1656 write(iulog,*) ' relative humidity ', q(i,k)/qs
1657 #ifdef WRF_PORT
1658 call wrf_message(iulog)
1659 #endif
1660 write(iulog,*) ' first guess ', tsp(i,k)
1661 #ifdef WRF_PORT
1662 call wrf_message(iulog)
1663 #endif
1664 endif
1665 #endif
1666 es = estblf(tsp(i,k))
1667 qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
1668 else
1669 doit(i) = 1
1670 tsp(i,k) = tlim(i)
1671 qsp(i,k) = q(i,k)
1672 enin(i) = 1._r8
1673 endif
1674 end do ! end do i
1675 !
1676 ! now iterate on first guess
1677 !
1678 do l = 1, iter
1679 dtm = 0
1680 dqm = 0
1681 do i = 1,ncol
1682 if (doit(i) == 0) then
1683 es = estblf(tsp(i,k))
1684 !
1685 ! Saturation specific humidity
1686 !
1687 qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
1688 !
1689 ! "generalized" analytic expression for t derivative of es
1690 ! accurate to within 1 percent for 173.16 < t < 373.16
1691 !
1692 ! Weighting of hlat accounts for transition from water to ice
1693 ! polynomial expression approximates difference between es over
1694 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1695 ! -40): required for accurate estimate of es derivative in transition
1696 ! range from ice to water also accounting for change of hlatv with t
1697 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1698 !
1699 tc = tsp(i,k) - t0
1700 lflg = (tc >= -ttrice .and. tc < 0.0_r8)
1701 weight(i) = min(-tc*trinv,1.0_r8)
1702 hlatsb = hlatv + weight(i)*hlatf
1703 hlatvp = hlatv - 2369.0_r8*tc
1704 if (tsp(i,k) < t0) then
1705 hltalt(i,k) = hlatsb
1706 else
1707 hltalt(i,k) = hlatvp
1708 end if
1709 if (lflg) then
1710 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
1711 else
1712 tterm = 0.0_r8
1713 end if
1714 desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
1715 dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
1716 ! g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
1717 g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
1718 dgdt = -(cp + hltalt(i,k)*dqsdt)
1719 t1 = tsp(i,k) - g/dgdt
1720 dt = abs(t1 - tsp(i,k))/t1
1721 tsp(i,k) = max(t1,tmin)
1722 es = estblf(tsp(i,k))
1723 q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
1724 dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
1725 qsp(i,k) = q1
1726 #ifdef DEBUG
1727 if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
1728 write(iulog,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
1729 #ifdef WRF_PORT
1730 call wrf_message(iulog)
1731 #endif
1732 endif
1733 #endif
1734 dtm = max(dtm,dt)
1735 dqm = max(dqm,dq)
1736 ! if converged at this point, exclude it from more iterations
1737 if (dt < dttol .and. dq < dqtol) then
1738 doit(i) = 2
1739 endif
1740 enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
1741 ! bail out if we are too near the end of temp range
1742 #if ( ! defined WACCM_PHYS )
1743 if (tsp(i,k) < 174.16_r8) then
1744 #else
1745 if (tsp(i,k) < 130.16_r8) then
1746 #endif
1747 doit(i) = 4
1748 endif
1749 else
1750 endif
1751 end do ! do i = 1,ncol
1753 if (dtm < dttol .and. dqm < dqtol) then
1754 go to 10
1755 endif
1757 end do ! do l = 1,iter
1758 10 continue
1760 error_found = .false.
1761 if (dtm > dttol .or. dqm > dqtol) then
1762 do i = 1,ncol
1763 if (doit(i) == 0) error_found = .true.
1764 end do
1765 if (error_found) then
1766 do i = 1,ncol
1767 if (doit(i) == 0) then
1768 write(iulog,*) ' findsp not converging at point i, k ', i, k
1769 #ifdef WRF_PORT
1770 call wrf_message(iulog)
1771 #endif
1772 write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
1773 #ifdef WRF_PORT
1774 call wrf_message(iulog)
1775 #endif
1776 write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
1777 #ifdef WRF_PORT
1778 call wrf_message(iulog)
1779 #endif
1780 call endrun ('FINDSP')
1781 endif
1782 end do
1783 endif
1784 endif
1785 do i = 1,ncol
1786 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
1787 error_found = .true.
1788 endif
1789 end do
1790 if (error_found) then
1791 do i = 1,ncol
1792 if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
1793 write(iulog,*) ' the enthalpy is not conserved for point ', &
1794 i, k, enin(i), enout(i)
1795 #ifdef WRF_PORT
1796 call wrf_message(iulog)
1797 #endif
1798 write(iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
1799 #ifdef WRF_PORT
1800 call wrf_message(iulog)
1801 #endif
1802 write(iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
1803 #ifdef WRF_PORT
1804 call wrf_message(iulog)
1805 #endif
1806 call endrun ('FINDSP')
1807 endif
1808 end do
1809 endif
1811 end do ! level loop (k=1,pver)
1813 return
1814 end subroutine findsp
1815 #endif
1816 end module cldwat