CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_MP_FER_HIRES.F90
1
7! FSMALL parameters are no longer used.
8! (3) T_ICE_init=-12 deg C provides a slight delay in the initial onset of ice.
9! (4) NLImax is a function of rime factor (RF) and temperature.
10! a) For RF>10, NLImax=1.e3. Mean ice diameters can exceed the 1 mm maximum
11! size in the tables so that NLICE=NLImax=1.e3.
12! b) Otherwise, NLImax is 10 L-1 at 0C and decreasing to 5 L-1 at <=-40C.
13! NLICE>NLImax at the maximum ice diameter of 1 mm.
14! (5) Can turn off ice processes by setting T_ICE & T_ICE_init to be < -100 deg C
15! (6) Modified the homogeneous freezing of cloud water when T<T_ICE.
16! (7) Reduce the fall speeds of rimed ice by making VEL_INC a function of
17! VrimeF**1.5 and not VrimeF**2.
18! (8) RHgrd=0.98 (98% RH) for the onset of condensation, to match what's been
19! tested for many months in the NAMX. Made obsolete with change in (13).
20! (9) Rime factor is *never* adjusted when NLICE>NLImax.
21! (10) Ice deposition does not change the rime factor (RF) when RF>=10 & T>T_ICE.
22! (11) Limit GAMMAS to <=1.5 (air resistance impact on ice fall speeds)
23! (12) NSImax is maximum # conc of ice crsytals. At cold temperature NSImax is
24! calculated based on assuming 10% of total ice content is due to cloud ice.
25!
26!-- Further modifications starting on 23 July 2015
27! (13) RHgrd is passed in as an input argument so that it can vary for different
28! domains (RHgrd=0.98 for 12-km parent, 1.0 for 3-km nests)
29! (14) Use the old "PRAUT" cloud water autoconversion *threshold* (QAUT0)
30
31!-- Further modifications starting on 28 July 2015
32! (15) Added calculations for radar reflectivity and number concentrations of
33! rain (Nrain) and precipitating ice (Nsnow).
34! (16) Removed double counting of air resistance term for riming onto ice (PIACW)
35! (17) The maximum rime factor (RFmx) is now a function of MASSI(INDEXS), accounting
36! for the increase in unrimed ice particle densities as values of INDEXS
37! decrease from the maximum upper limit of 1000 microns to the lower limit of
38! 50 microns, coinciding with the assumed size of cloud ice; see lines 1128-1134.
39! (18) A new closure is used for updating the rime factor, which is described in
40! detail near lines 1643-1682. The revised code is near lines 1683-1718.
41! (19) Restructured the two-pass algorithm to be more robust, removed the HAIL
42! & LARGE_RF logical variables so that NLICE>NLImax can occur.
43! (20) Increased nsimax (see !aug27 below)
44! (21) Modified the rain sedimentation (see two !aug27 blocks below)
45! (22) NInuclei is the lower of Fletcher (1962), Cooper (1986), or NSImax.
46! (23) NLImax is no longer used or enforced. Instead, INDEXS=MDImax when RF>20,
47! else INDEXS is a function of temperature. Look for !sep10 comment.
48! (24) An override was inserted for (18), such that the rime density is not diluted
49! diluted when RF>20. Look for !sep10 comment.
50! (25) Radar reflectivity calculations were changes to reduce radar bright bands,
51! limit enhanced, mixed-phase reflectivity to RF>=20. Look for !sep10 comments.
52! (26) NLICE is not to exceed NSI_max (250 L^-1) when RF<20. Look for !sep16 comments.
53! Commented out! (28) Increase hail fall speeds using Thompson et al. (2008). Look for !sep22 comments.
54! (29) Modify NLImax, INDEXS for RF>=20. Look for !sep22 comments.
55! (30) Check on NSmICE, Vci based on whether FLIMASS<1. Look for !sep22a comments.
56! Revised in (34)! (31) Introduced RFlag logical, which if =T enforces a lower limit of drop sizes not
57! to go below INDEXRmin and N0r is adjusted. Look for !nov25 comments (corrections,
58! refinements to sep25 & nov18 versions, includes an additional fix in nov25-fix).
59! Also set INDEXRmin=500 rather than 250 microns.
60!-----------------------------------------------------------------------------
61!--- The following changes now refer to dates when those were made in 2016.
62!-----------------------------------------------------------------------------
63! (32) Convective (RF>=20, Ng~10 L^-1, RHOg~500 kg m^-3), transition (RF=10, Ng~25 L^-1,
64! RHOg~300 kg m^-3), & stratiform (RF<2) profiles are blended based on RF. !mar08
65! (33) Fixed bug in Biggs' freezing, put back in collisional drop freezing. !mar03
66! (34) Changes in (31) are revised so that INDEXRmin at and below 0C level is
67! based on a rain rate equal to the snowfall rate above the 0C level. !mar03
68! (35) Increase radar reflectivity when RF>10 and RQSnew > 2.5 g m^-3. !mar12
69! (36) !mar10 combines all elements of (32)-(35) together.
70! (37) Bug fixes for the changes in (34) and the RFLAG variable !apr18
71! (38) Revised Schumann-Ludlam limit. !apr18
72! (39) Simplified PCOND (cloud cond/evap) calculation !apr21
73! (40) Slight change in calculating RF. !apr22
74! (41) Reduce RF values for calculating mean sizes of snow, graupel, sleet/hail !apr22a
75! (42) Increase reflectivity from large, wet, high rime factor ice (graupel) by
76! assuming |Kw|**2/|Ki|**2 = 0.224 (Smith, 1984, JCAM).
77! (43) Major restructuring of code to allow N0r to vary from N0r0 !may11
78! (44) More major restructuring of code to use fixed XLS, XLV, XLF !may12
79! (45) Increased VEL_INC ~ VrimeF**2, put the enhanced graupel/hail fall speeds
80! from Thompson into the code but only in limited circumstances, restructured
81! and streamlined the INDEXS calculation, removed the upper limit for
82! for the vapor mixing ratio is at water saturation when calculating ice
83! deposition, and N0r is gradually increased for conditions supporting
84! drizzle when rain contents decrease below 0.25 g/m**3. !may17
85! (46) The may11 code changes that increase N0r0 when rain contents exceed 1 g m^-3
86! have been removed, limit the number of iterations calculating final rain
87! parameters, remove the revised N0r calculation for reflectivity. All of
88! the changes following those made in the may10 code. !may20
89! (47) Reduce the assumed # concentration of hail/sleet when RF>10 from 5 L^-1 to
90! 1 L^-1, and also reduce it for graupel when RF>5 from 10 L^-1 to 5 L^-1.
91! This is being done to try and make greater use of the Thompson graupel/hail
92! fallspeeds by having INDEXS==MDImax.
93! (48) Increased NCW from 200e6 to 300e6 for a more delayed onset of drizzle,
94! simplified drizzle algorithm to reduce/eliminate N0r bulls eyes and to allow
95! for supercooled drizzle, and set limits for 8.e6 <= N0r (m^-4) <= 1.e9 !may31
96! (49) Further restructuring of code to better define STRAT, DRZL logicals,
97! add these rain flags to mprates arrays !jun01
98! (50) Increase in reflectivity due to wet ice was commented out.
99! (51) Fixed minor bug to update INDEXR2 in the "rain_pass: do" loop. !jun13
100! (52) Final changes to Nsnow for boosting reflectivities from ice for
101! mass contents exceeding 5 g m^-3. !jun16
102! (53) Cosmetic changes only that do not affect the calculations. Removed old, unused
103! diagnostic arrays. Updated comments.
104!
105!-----------------------------------------------------------------------------
106!
108!
109!-----------------------------------------------------------------------------
110
111#ifdef MPI
112 USE mpi
113#endif
114 USE machine
115!MZ
116!MZ USE MODULE_CONSTANTS,ONLY : PI, CP, EPSQ, GRAV=>G, RHOL=>RHOWATER, &
117!MZ RD=>R_D, RV=>R_V, T0C=>TIW, EPS=>EP_2, EPS1=>EP_1, CLIQ, CICE, &
118!MZ XLV
119!MZ
120!MZ temporary values copied from module_CONSTANTS; ideally they come from host model
121!side
122 REAL, PARAMETER :: pi=3.141592653589793 ! ludolf number
123 REAL, PARAMETER :: cp=1004.6 ! spec. heat for dry air at constant pressure
124 REAL, PARAMETER :: epsq=1.e-12 ! floor value for specific humidity (kg/kg)
125 REAL, PARAMETER :: grav= 9.8060226 ! gravity
126 REAL, PARAMETER :: rhol=1000. ! density of water (kg/m3)
127 REAL, PARAMETER :: rd=287.04 ! gas constant for dry air
128 REAL, PARAMETER :: rv=461.6 ! gas constant for water vapor
129 REAL, PARAMETER :: t0c= 273.15 ! melting point
130 REAL, PARAMETER :: eps=rd/rv
131 REAL, PARAMETER :: eps1=rv/rd-1.
132 REAL, PARAMETER :: cliq = 4190. ! MZ: inconsistent value below
133 REAL, PARAMETER :: cice = 2106.
134 REAL, PARAMETER :: xlv = 2.5e6
135!-----------------------------------------------------------------------------
136 PUBLIC :: ferrier_init_hr, gpvs_hr,fpvs,fpvs0,nx
137!-----------------------------------------------------------------------------
138 REAL,PRIVATE,SAVE :: abfr, cbfr, ciacw, ciacr, c_n0r0, c_nr, crain, & !jul28
139 & CRACW, ARAUT, BRAUT, ESW0, RFmx1, ARcw, RH_NgC, RH_NgT, & !jul31 !mar08
140 & RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DR4, RR_DR5, RR_DRmax, & !may17
141 & BETA6, &
142 & RQhail, AVhail, BVhail, QAUT0 !may17
143!
144 INTEGER,PRIVATE,PARAMETER :: indexrstrmax=500 !mar03, stratiform maximum
145 REAL,PUBLIC,SAVE :: cn0r0, cn0r_dmrmin, cn0r_dmrmax, &
146 rfmax, rqr_drmax, rqr_drmin
147!
148 INTEGER, PRIVATE,PARAMETER :: my_t1=1, my_t2=35
149 REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: my_growth_nmm
150!
151 REAL, PRIVATE,PARAMETER :: dmimin=.05e-3, dmimax=1.e-3, &
152 & deldmi=1.e-6,xmimin=1.e6*dmimin
153 REAL, PUBLIC,PARAMETER :: xmimax=1.e6*dmimax, xmiexp=.0536
154 INTEGER, PUBLIC,PARAMETER :: mdimin=xmimin, mdimax=xmimax
155 REAL, ALLOCATABLE, DIMENSION(:) :: &
156 & ACCRI,VSNOWI,VENTI1,VENTI2
157 REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: sdens !-- For RRTM
158!
159 REAL, PRIVATE,PARAMETER :: dmrmin=.05e-3, dmrmax=1.0e-3, &
160 & deldmr=1.e-6, xmrmin=1.e6*dmrmin, xmrmax=1.e6*dmrmax
161 INTEGER, PUBLIC,PARAMETER :: mdrmin=xmrmin, mdrmax=xmrmax
162!
163 REAL, ALLOCATABLE, DIMENSION(:):: &
164 & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
165!
166 INTEGER, PRIVATE,PARAMETER :: nrime=40
167 REAL, ALLOCATABLE, DIMENSION(:,:) :: vel_rf
168!
169 INTEGER,PARAMETER :: nx=7501
170 REAL, PARAMETER :: xmin=180.0,xmax=330.0
171 REAL, DIMENSION(NX),PUBLIC,SAVE :: tbpvs,tbpvs0
172 REAL, PUBLIC,SAVE :: c1xpvs0,c2xpvs0,c1xpvs,c2xpvs
173!
174 REAL,DIMENSION(MY_T2+8) :: mp_restart_state
175 REAL,DIMENSION(nx) :: tbpvs_state,tbpvs0_state
176!
177 REAL, PRIVATE,PARAMETER :: cvap=1846., xlf=3.3358e+5, xls=xlv+xlf &
178 & ,epsq1=1.001*epsq, rcp=1./cp, rcprv=rcp/rv, rgrav=1./grav &
179 & ,rrhol=1./rhol, xlv1=xlv/cp, xlf1=xlf/cp, xls1=xls/cp &
180 & ,xlv2=xlv*xlv/rv, xls2=xls*xls/rv &
181 & ,xlv3=xlv*xlv*rcprv, xls3=xls*xls*rcprv &
182!--- Constants specific to the parameterization follow:
183!--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
184 & ,climit=10.*epsq, climit1=-climit &
185 & ,c1=1./3. &
186 & ,dmr1=.1e-3, dmr2=.2e-3, dmr3=.32e-3, dmr4=0.45e-3 &
187 & ,dmr5=0.67e-3 &
188 & ,xmr1=1.e6*dmr1, xmr2=1.e6*dmr2, xmr3=1.e6*dmr3 &
189 & ,xmr4=1.e6*dmr4, xmr5=1.e6*dmr5, rqrmix=0.05e-3, rqsmix=1.e-3 & !jul28 !apr27
190 & ,cdry=1.634e13, cwet=1./.224 !jul28 !apr27
191 INTEGER, PARAMETER :: mdr1=xmr1, mdr2=xmr2, mdr3=xmr3, mdr4=xmr4 &
192 & , mdr5=xmr5
193
194!-- Debug 20120111
195LOGICAL, SAVE :: warn1=.true.,warn2=.true.,warn3=.true.,warn5=.true.
196REAL, SAVE :: pwarn=75.e2, qtwarn=1.e-3
197INTEGER, PARAMETER :: max_iterations=10
198
199!
200! ======================================================================
201!--- Important tunable parameters that are exported to other modules
202! * T_ICE - temperature (C) threshold at which all remaining liquid water
203! is glaciated to ice
204! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
205!
206!-- To turn off ice processes, set T_ICE & T_ICE_init to <= -100. (i.e., -100 C)
207!
208! * NSImax - maximum number concentrations (m**-3) of small ice crystals
209! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
210! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 1.0 mm
211! * N0rmin - minimum intercept (m**-4) for rain drops
212! * NCW - number concentrations of cloud droplets (m**-3)
213! ======================================================================
214 REAL, PUBLIC,PARAMETER :: &
215 & RHgrd_in=1. &
216 &, P_RHgrd_out=850.E2 &
217 & ,T_ICE=-40. &
218 & ,T_ICEK=T0C+T_ICE &
219 & ,T_ICE_init=-12. &
220 & ,NSI_max=250.E3 &
221 & ,NLImin=1.0E3 &
222 & ,N0r0=8.E6 &
223 & ,N0rmin=1.E4 &
224!! based on Aligo's email,NCW is changed to 250E6
225 & ,ncw=250.e6
226 !HWRF & ,NCW=300.E6 !- 100.e6 (maritime), 500.e6 (continental)
227
228!--- Other public variables passed to other routines:
229 REAL, ALLOCATABLE ,DIMENSION(:) :: massi
230!
231
232 CONTAINS
233!-----------------------------------------------------------------------
234!-----------------------------------------------------------------------
235!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
236!-----------------------------------------------------------------------
237
244 SUBROUTINE fer_hires (DT,RHgrd, &
245 & prsi,p_phy,t_phy, &
246 & q,qt, &
247 & LOWLYR,SR,TRAIN_PHY, &
248 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
249 & QC,QR,QS, &
250 & RAINNC,RAINNCV, &
251 & threads, &
252 & ims,ime, lm, &
253 & d_ss, &
254 & refl_10cm,DX1 )
255!-----------------------------------------------------------------------
256 IMPLICIT NONE
257!-----------------------------------------------------------------------
258 INTEGER,INTENT(IN) :: D_SS,IMS,IME,LM,DX1
259 REAL, INTENT(IN) :: DT,RHgrd
260 INTEGER, INTENT(IN) :: THREADS
261 REAL, INTENT(IN), DIMENSION(ims:ime, lm+1):: &
262 & prsi
263 REAL, INTENT(IN), DIMENSION(ims:ime, lm):: &
264 & p_phy
265 REAL, INTENT(INOUT), DIMENSION(ims:ime, lm):: &
266 & q,qt,t_phy
267 REAL, INTENT(INOUT), DIMENSION(ims:ime, lm ):: & !Aligo Oct 23,2019: dry mixing ratio for cloud species
268 & qc,qr,qs
269 REAL, INTENT(INOUT), DIMENSION(ims:ime, lm) :: &
270 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
271 REAL, INTENT(OUT), DIMENSION(ims:ime, lm) :: &
272 & refl_10cm
273 REAL, INTENT(INOUT), DIMENSION(ims:ime) :: &
274 & RAINNC,RAINNCV
275 REAL, INTENT(OUT), DIMENSION(ims:ime):: SR
276 REAL, INTENT(OUT), DIMENSION( ims:ime, lm ) :: &
277 & TRAIN_PHY
278!
279 INTEGER, DIMENSION( ims:ime ),INTENT(INOUT) :: LOWLYR
280
281!-----------------------------------------------------------------------
282! LOCAL VARS
283!-----------------------------------------------------------------------
284
285! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
286! the microphysics scheme. Instead, they will be used by Eta precip
287! assimilation.
288
289 REAL, DIMENSION(ims:ime):: APREC,PREC,ACPREC
290
291 INTEGER :: I,K,KK
292 REAL :: wc, RDIS, BETA6
293!------------------------------------------------------------------------
294! For subroutine EGCP01COLUMN_hr
295!-----------------------------------------------------------------------
296 INTEGER :: LSFC,I_index,J_index,L
297 INTEGER,DIMENSION(ims:ime) :: LMH
298 REAL :: TC,QI,QRdum,QW,Fice,Frain,DUM,ASNOW,ARAIN
299 REAL,DIMENSION(lm) :: P_col,Q_col,T_col,WC_col, &
300 RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL,pcond1d, &
301 pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d, &
302 pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col, &
303 NR_col,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d, &
304 INDEXS1d,INDEXR1d,RFlag1d,RHC_col
305!
306!-----------------------------------------------------------------------
307!**********************************************************************
308!-----------------------------------------------------------------------
309!
310
311! MZ: HWRF start
312!----------
313!2015-03-30, recalculate some constants which may depend on phy time step
314 CALL my_growth_rates_nmm_hr (dt)
315
316!--- CIACW is used in calculating riming rates
317! The assumed effective collection efficiency of cloud water rimed onto
318! ice is =0.5 below:
319!
320 ciacw=dt*0.25*pi*0.5*(1.e5)**c1
321!
322!--- CIACR is used in calculating freezing of rain colliding with large ice
323! The assumed collection efficiency is 1.0
324!
325 ciacr=pi*dt
326!
327!--- CRACW is used in calculating collection of cloud water by rain (an
328! assumed collection efficiency of 1.0)
329!
330 cracw=dt*0.25*pi*1.0
331!
332!-- See comments in subroutine etanewhr_init starting with variable RDIS=
333!
334!-- Relative dispersion == standard deviation of droplet spectrum / mean radius
335! (see pp 1542-1543, Liu & Daum, JAS, 2004)
336 rdis=0.5 !-- relative dispersion of droplet spectrum
337 beta6=( (1.+3.*rdis*rdis)*(1.+4.*rdis*rdis)*(1.+5.*rdis*rdis)/ &
338 & ((1.+rdis*rdis)*(1.+2.*rdis*rdis) ) )
339
340 braut=dt*1.1e10*beta6/ncw
341
342!! END OF adding, 2015-03-30
343!-----------
344! MZ: HWRF end
345!
346
347 DO i = ims,ime
348 acprec(i)=0.
349 aprec(i)=0.
350 prec(i)=0.
351 sr(i)=0.
352 ENDDO
353
354 DO k = 1,lm
355 DO i = ims,ime
356 train_phy(i,k)=0.
357 ENDDO
358 ENDDO
359
360!-----------------------------------------------------------------------
361!-- Start of original driver for EGCP01COLUMN_hr
362!-----------------------------------------------------------------------
363!
364 DO i=ims,ime
365 lsfc=lm-lowlyr(i)+1 ! "L" of surface
366 DO k=1,lm
367 dpcol(k)=prsi(i,k)-prsi(i,k+1)
368 ENDDO
369!
370!--- Initialize column data (1D arrays)
371!
372 l=lm
373!-- qt = CWM, total condensate
374 IF (qt(i,l) .LE. epsq) qt(i,l)=epsq
375 f_ice_phy(i,l)=1.
376 f_rain_phy(i,l)=0.
377 f_rimef_phy(i,l)=1.
378 do l=lm,1,-1
379 p_col(l)=p_phy(i,l)
380!
381!--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
382!
383 thick_col(l)=dpcol(l)*rgrav
384 t_col(l)=t_phy(i,l)
385 tc=t_col(l)-t0c
386 q_col(l)=max(epsq, q(i,l))
387 IF (qt(i,l) .LE. epsq1) THEN
388 wc_col(l)=0.
389 IF (tc .LT. t_ice) THEN
390 f_ice_phy(i,l)=1.
391 ELSE
392 f_ice_phy(i,l)=0.
393 ENDIF
394 f_rain_phy(i,l)=0.
395 f_rimef_phy(i,l)=1.
396 ELSE
397 wc_col(l)=qt(i,l)
398
399!-- Debug 20120111
400! TC==TC will fail if NaN, preventing unnecessary error messages
401IF (wc_col(l)>qtwarn .AND. p_col(l)<pwarn .AND. tc==tc) THEN
402 WRITE(0,*) 'WARN4: >1 g/kg condensate in stratosphere; I,L,TC,P,QT=', &
403 i,l,tc,.01*p_col(l),1000.*wc_col(l)
404 qtwarn=max(wc_col(l),10.*qtwarn)
405 pwarn=min(p_col(l),0.5*pwarn)
406ENDIF
407!-- TC/=TC will pass if TC is NaN
408IF (warn5 .AND. tc/=tc) THEN
409 WRITE(0,*) 'WARN5: NaN temperature; I,L,P=',i,l,.01*p_col(l)
410 warn5=.false.
411ENDIF
412
413 ENDIF
414 IF (t_ice<=-100.) f_ice_phy(i,l)=0.
415! !
416! !--- Determine composition of condensate in terms of
417! ! cloud water, ice, & rain
418! !
419 wc=wc_col(l)
420 qi=0.
421 qrdum=0.
422 qw=0.
423 fice=f_ice_phy(i,l)
424 frain=f_rain_phy(i,l)
425!
426 IF (fice .GE. 1.) THEN
427 qi=wc
428 ELSE IF (fice .LE. 0.) THEN
429 qw=wc
430 ELSE
431 qi=fice*wc
432 qw=wc-qi
433 ENDIF
434!
435 IF (qw.GT.0. .AND. frain.GT.0.) THEN
436 IF (frain .GE. 1.) THEN
437 qrdum=qw
438 qw=0.
439 ELSE
440 qrdum=frain*qw
441 qw=qw-qrdum
442 ENDIF
443 ENDIF
444 IF (qi .LE. 0.) f_rimef_phy(i,l)=1.
445 rimef_col(l)=f_rimef_phy(i,l) ! (real)
446 qi_col(l)=qi
447 qr_col(l)=qrdum
448 qw_col(l)=qw
449!GFDL => New. Added RHC_col to allow for height- and grid-dependent values for
450!GFDL the relative humidity threshold for condensation ("RHgrd")
451!6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa
452!mz 05/06/2020 - 10km
453!------------------------------------------------------------
454 IF(dx1 .GE. 10000 .AND. p_col(l)<p_rhgrd_out) THEN ! gopal's doing based on GFDL
455 rhc_col(l)=rhgrd
456 ELSE
457 rhc_col(l)=rhgrd_in
458 ENDIF
459
460 ENDDO
461!
462!#######################################################################
463!
464!--- Perform the microphysical calculations in this column
465!
466 i_index=i
467 j_index=1
468 CALL egcp01column_hr ( arain, asnow, dt, rhc_col, &
469 & i_index, j_index, lsfc, &
470 & p_col, qi_col, qr_col, q_col, qw_col, rimef_col, t_col, &
471 & thick_col, wc_col,lm,pcond1d,pidep1d, &
472 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
473 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, dbz_col,nr_col,ns_col, &
474 & vsnow1d,vrain11d,vrain21d,vci1d,nsmice1d,indexs1d,indexr1d, & !jul28
475 & rflag1d,dx1) !jun01
476!#######################################################################
477!
478!--- Update storage arrays
479!
480 do l=lm,1,-1
481 train_phy(i,l)=(t_col(l)-t_phy(i,l))/dt
482 t_phy(i,l)=t_col(l)
483 q(i,l)=q_col(l)
484 qt(i,l)=wc_col(l)
485!---convert 1D source/sink terms to one 4D array
486!---d_ss is the total number of source/sink terms in the 4D mprates array
487!---if d_ss=1, only 1 source/sink term is used
488!
489
490!--- REAL*4 array storage
491!
492 IF (qi_col(l) .LE. epsq) THEN
493 f_ice_phy(i,l)=0.
494 IF (t_col(l) .LT. t_icek) f_ice_phy(i,l)=1.
495 f_rimef_phy(i,l)=1.
496 ELSE
497 f_ice_phy(i,l)=max( 0., min(1., qi_col(l)/wc_col(l)) )
498 f_rimef_phy(i,l)=max(1., rimef_col(l))
499 ENDIF
500 IF (qr_col(l) .LE. epsq) THEN
501 dum=0
502 ELSE
503 dum=qr_col(l)/(qr_col(l)+qw_col(l))
504 ENDIF
505 f_rain_phy(i,l)=dum
506 refl_10cm(i,l)=dbz_col(l) !jul28
507 ENDDO
508!
509!--- Update accumulated precipitation statistics
510!
511!--- Surface precipitation statistics; SR is fraction of surface
512! precipitation (if >0) associated with snow
513!
514 aprec(i)=(arain+asnow)*rrhol ! Accumulated surface precip (depth in m) !<--- Ying
515 prec(i)=prec(i)+aprec(i)
516 acprec(i)=acprec(i)+aprec(i)
517 IF(aprec(i) .LT. 1.e-8) THEN
518 sr(i)=0.
519 ELSE
520 sr(i)=rrhol*asnow/aprec(i)
521 ENDIF
522!
523!#######################################################################
524!#######################################################################
525!
526 enddo ! End "I" loop
527!
528!-----------------------------------------------------------------------
529!-- End of original driver for EGCP01COLUMN_hr
530!-----------------------------------------------------------------------
531!
532 do k = lm, 1, -1
533 DO i = ims,ime
534 wc=qt(i,k)
535 qs(i,k)=0.
536 qr(i,k)=0.
537 qc(i,k)=0.
538!
539 IF(f_ice_phy(i,k)>=1.)THEN
540 qs(i,k)=wc
541 ELSEIF(f_ice_phy(i,k)<=0.)THEN
542 qc(i,k)=wc
543 ELSE
544 qs(i,k)=f_ice_phy(i,k)*wc
545 qc(i,k)=wc-qs(i,k)
546 ENDIF
547!
548 IF(qc(i,k)>0..AND.f_rain_phy(i,k)>0.)THEN
549 IF(f_rain_phy(i,k).GE.1.)THEN
550 qr(i,k)=qc(i,k)
551 qc(i,k)=0.
552 ELSE
553 qr(i,k)=f_rain_phy(i,k)*qc(i,k)
554 qc(i,k)=qc(i,k)-qr(i,k)
555 ENDIF
556 ENDIF
557 ENDDO !- i
558 ENDDO !- k
559!
560!- Update rain (convert from m to kg/m**2, which is also equivalent to mm depth)
561!
562 DO i=ims,ime
563 rainnc(i)=aprec(i)*1000.+rainnc(i)
564 rainncv(i)=aprec(i)*1000.
565 ENDDO
566!
567!-----------------------------------------------------------------------
568!
569 END SUBROUTINE fer_hires
570!
571!-----------------------------------------------------------------------
572!
573!###############################################################################
574! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
575! (1) Represents sedimentation by preserving a portion of the precipitation
576! through top-down integration from cloud-top. Modified procedure to
577! Zhao and Carr (1997).
578! (2) Microphysical equations are modified to be less sensitive to time
579! steps by use of Clausius-Clapeyron equation to account for changes in
580! saturation mixing ratios in response to latent heating/cooling.
581! (3) Prevent spurious temperature oscillations across 0C due to
582! microphysics.
583! (4) Uses lookup tables for: calculating two different ventilation
584! coefficients in condensation and deposition processes; accretion of
585! cloud water by precipitation; precipitation mass; precipitation rate
586! (and mass-weighted precipitation fall speeds).
587! (5) Assumes temperature-dependent variation in mean diameter of large ice
588! (Houze et al., 1979; Ryan et al., 1996).
589! -> 8/22/01: This relationship has been extended to colder temperatures
590! to parameterize smaller large-ice particles down to mean sizes of MDImin,
591! which is 50 microns reached at -55.9C.
592! (6) Attempts to differentiate growth of large and small ice, mainly for
593! improved transition from thin cirrus to thick, precipitating ice
594! anvils.
595! (7) Top-down integration also attempts to treat mixed-phase processes,
596! allowing a mixture of ice and water. Based on numerous observational
597! studies, ice growth is based on nucleation at cloud top &
598! subsequent growth by vapor deposition and riming as the ice particles
599! fall through the cloud. There are two modes of ice nucleation
600! following Meyers et al. (JAM, 1992):
601! a) Deposition & condensation freezing nucleation - eq. (2.4) when
602! air is supersaturated w/r/t ice
603! b) Contact freezing nucleation - eq. (2.6) in presence of cloud water
604! (8) Depositional growth of newly nucleated ice is calculated for large time
605! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
606! using their ice crystal masses calculated after 600 s of growth in water
607! saturated conditions. The growth rates are normalized by time step
608! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
609! (9) Ice precipitation rates can increase due to increase in response to
610! cloud water riming due to (a) increased density & mass of the rimed
611! ice, and (b) increased fall speeds of rimed ice.
612!###############################################################################
613!###############################################################################
614!
662 SUBROUTINE egcp01column_hr ( ARAIN, ASNOW, DTPH, RHC_col, &
663 & I_index, J_index, LSFC, &
664 & P_col, QI_col, QR_col, Q_col, QW_col, RimeF_col, T_col, &
665 & THICK_col, WC_col ,LM,pcond1d,pidep1d, &
666 & piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d,pimlt1d, &
667 & praut1d,pracw1d,prevp1d,pisub1d,pevap1d, DBZ_col,NR_col,NS_col, &
668 & vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d,INDEXR1d, & !jul28
669 & RFlag1d,DX1) !jun01
670!
671!###############################################################################
672!###############################################################################
673!
674!-------------------------------------------------------------------------------
675!----- NOTE: Code is currently set up w/o threading!
676!-------------------------------------------------------------------------------
677!$$$ SUBPROGRAM DOCUMENTATION BLOCK
678! . . .
679! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
680! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
681! PRGRMMR: Jin (Modification for WRF structure)
682!-------------------------------------------------------------------------------
683! ABSTRACT:
684! * Merges original GSCOND & PRECPD subroutines.
685! * Code has been substantially streamlined and restructured.
686! * Exchange between water vapor & small cloud condensate is calculated using
687! the original Asai (1965, J. Japan) algorithm. See also references to
688! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
689! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
690! parameterization.
691!-------------------------------------------------------------------------------
692!
693! USAGE:
694! * CALL EGCP01COLUMN_hr FROM SUBROUTINE EGCP01DRV
695!
696! INPUT ARGUMENT LIST:
697! DTPH - physics time step (s)
698! RHgrd - threshold relative humidity (ratio) for onset of condensation
699! I_index - I index
700! J_index - J index
701! LSFC - Eta level of level above surface, ground
702! P_col - vertical column of model pressure (Pa)
703! QI_col - vertical column of model ice mixing ratio (kg/kg)
704! QR_col - vertical column of model rain ratio (kg/kg)
705! Q_col - vertical column of model water vapor specific humidity (kg/kg)
706! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
707! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
708! T_col - vertical column of model temperature (deg K)
709! THICK_col - vertical column of model mass thickness (density*height increment)
710! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
711! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL
712!
713!
714! OUTPUT ARGUMENT LIST:
715! ARAIN - accumulated rainfall at the surface (kg)
716! ASNOW - accumulated snowfall at the surface (kg)
717! Q_col - vertical column of model water vapor specific humidity (kg/kg)
718! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
719! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
720! QI_col - vertical column of model ice mixing ratio (kg/kg)
721! QR_col - vertical column of model rain ratio (kg/kg)
722! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
723! T_col - vertical column of model temperature (deg K)
724! DBZ_col - vertical column of radar reflectivity (dBZ)
725! NR_col - vertical column of rain number concentration (m^-3)
726! NS_col - vertical column of snow number concentration (m^-3)
727!
728! OUTPUT FILES:
729! NONE
730!
731! Subprograms & Functions called:
732! * Real Function CONDENSE - cloud water condensation
733! * Real Function DEPOSIT - ice deposition (not sublimation)
734! * Integer Function GET_INDEXR - estimate the mean size of raindrops (microns)
735!
736! UNIQUE: NONE
737!
738! LIBRARY: NONE
739!
740! ATTRIBUTES:
741! LANGUAGE: FORTRAN 90
742! MACHINE : IBM SP
743!
744!-------------------------------------------------------------------------
745!--------------- Arrays & constants in argument list ---------------------
746!-------------------------------------------------------------------------
747!
748 IMPLICIT NONE
749!
750 INTEGER,INTENT(IN) :: LM,I_index, J_index, LSFC,DX1
751 REAL,INTENT(IN) :: DTPH
752 REAL,INTENT(INOUT) :: ARAIN, ASNOW
753 REAL,DIMENSION(LM),INTENT(INOUT) :: P_col, QI_col,QR_col &
754 & ,Q_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col,pcond1d &
755 & ,pidep1d,piacw1d,piacwi1d,piacwr1d,piacr1d,picnd1d,pievp1d &
756 & ,pimlt1d,praut1d,pracw1d,prevp1d,pisub1d,pevap1d,DBZ_col,NR_col &
757 & ,NS_col,vsnow1d,vrain11d,vrain21d,vci1d,NSmICE1d,INDEXS1d & !jun01
758 & ,INDEXR1d,RFlag1d,RHC_col !jun01
759!
760!--------------------------------------------------------------------------------
761!--- The following arrays are integral calculations based on the mean
762! snow/graupel diameters, which vary from 50 microns to 1000 microns
763! (1 mm) at 1-micron intervals and assume exponential size distributions.
764! The values are normalized and require being multipled by the number
765! concentration of large ice (NLICE).
766!---------------------------------------
767! - VENTI1 - integrated quantity associated w/ ventilation effects
768! (capacitance only) for calculating vapor deposition onto ice
769! - VENTI2 - integrated quantity associated w/ ventilation effects
770! (with fall speed) for calculating vapor deposition onto ice
771! - ACCRI - integrated quantity associated w/ cloud water collection by ice
772! - MASSI - integrated quantity associated w/ ice mass
773! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
774! precipitation rates
775! - VEL_RF - velocity increase of rimed particles as functions of crude
776! particle size categories (at 0.1 mm intervals of mean ice particle
777! sizes) and rime factor (different values of Rime Factor of 1.1**N,
778! where N=0 to Nrime).
779!--------------------------------------------------------------------------------
780!--- The following arrays are integral calculations based on the mean
781! rain diameters, which vary from 50 microns to 1000 microns
782! (1 mm) at 1-micron intervals and assume exponential size distributions.
783! The values are normalized and require being multiplied by the rain intercept
784! (N0r).
785!---------------------------------------
786! - VENTR1 - integrated quantity associated w/ ventilation effects
787! (capacitance only) for calculating evaporation from rain
788! - VENTR2 - integrated quantity associated w/ ventilation effects
789! (with fall speed) for calculating evaporation from rain
790! - ACCRR - integrated quantity associated w/ cloud water collection by rain
791! - MASSR - integrated quantity associated w/ rain
792! - VRAIN - mass-weighted fall speed of rain, used to calculate
793! precipitation rates
794! - RRATE - precipitation rates, which should also be equal to RHO*QR*VRAIN
795!
796!-------------------------------------------------------------------------
797!------- Key parameters, local variables, & important comments ---------
798!-----------------------------------------------------------------------
799!
800!--- TOLER => Tolerance or precision for accumulated precipitation
801!
802 REAL, PARAMETER :: TOLER=5.e-7, c2=1./6., rho0=1.194, &
803 xratio=.025, zmin=0.01, dbzmin=-20.
804!
805!--- If BLEND=1:
806! precipitation (large) ice amounts are estimated at each level as a
807! blend of ice falling from the grid point above and the precip ice
808! present at the start of the time step (see TOT_ICE below).
809!--- If BLEND=0:
810! precipitation (large) ice amounts are estimated to be the precip
811! ice present at the start of the time step.
812!
813!--- Extended to include sedimentation of rain on 2/5/01
814!
815 REAL, PARAMETER :: BLEND=1.
816!
817!--- This variable is for debugging purposes (if .true.)
818!
819 LOGICAL, PARAMETER :: PRINT_diag=.false.
820!
821!-----------------------------------------------------------------------
822!--- Local variables
823!-----------------------------------------------------------------------
824!
825 REAL :: EMAIRI, N0r, NLICE, NSmICE, NInuclei, Nrain, Nsnow, Nmix
826 REAL :: RHgrd
827 LOGICAL :: CLEAR, ICE_logical, DBG_logical, RAIN_logical, &
828 STRAT, DRZL
829 INTEGER :: INDEX_MY,INDEXR,INDEXR1,INDEXR2,INDEXS,IPASS,ITDX,IXRF,&
830 & IXS,LBEF,L,INDEXRmin,INDEXS0C,IDR !mar03 !may20
831!
832!
833 REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
834 & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
835 & DENOMI,DENOMW,DENOMWI,DIDEP, &
836 & DIEVP,DIFFUS,DLI,DTRHO,DUM,DUM1,DUM2,DUM3, &
837 & DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLIMASS, &
838 & FWR,FWS,GAMMAR,GAMMAS, &
839 & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
840 & PIEVP,PILOSS,PIMLT,PINIT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
841 & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
842 & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,Q,QW,QWnew,Rcw, &
843 & RFACTOR,RFmx,RFrime,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
844 & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
845 & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
846 & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
847 & VSNOW1,WC,WCnew,WSgrd,WS,WSnew,WV,WVnew, &
848 & XLI,XLIMASS,XRF, &
849 & NSImax,QRdum,QSmICE,QLgIce,RQLICE,VCI,TIMLT, &
850 & RQSnew,RQRnew,Zrain,Zsnow,Ztot,RHOX0C,RFnew,PSDEP,DELS !mar03 !apr22
851 REAL, SAVE :: Revised_LICE=1.e-3 !-- kg/m**3
852!
853!#######################################################################
854!########################## Begin Execution ############################
855!#######################################################################
856!
857!
858 arain=0. ! Accumulated rainfall into grid box from above (kg/m**2)
859 vrain1=0. ! Rain fall speeds into grib box from above (m/s)
860 vsnow1=0. ! Ice fall speeds into grib box from above (m/s)
861 asnow=0. ! Accumulated snowfall into grid box from above (kg/m**2)
862 indexs0c=mdimin ! Mean snow/graupel diameter just above (<0C) freezing level (height)
863 rhox0c=22.5 ! Estimated ice density at 0C (kg m^-3) !mar03
864 timlt=0. ! Total ice melting in a layer (drizzle detection)
865 strat=.false. ! Stratiform rain DSD below melting level !may11
866 drzl=.false. ! Drizzle DSD below melting level !may23
867!
868!-----------------------------------------------------------------------
869!------------ Loop from top (L=1) to surface (L=LSFC) ------------------
870!-----------------------------------------------------------------------
871!
872big_loop: DO l=lm,1,-1
873 pcond1d(l)=0.
874 pidep1d(l)=0.
875 piacw1d(l)=0.
876 piacwi1d(l)=0.
877 piacwr1d(l)=0.
878 piacr1d(l)=0.
879 picnd1d(l)=0.
880 pievp1d(l)=0.
881 pimlt1d(l)=0.
882 praut1d(l)=0.
883 pracw1d(l)=0.
884 prevp1d(l)=0.
885 pisub1d(l)=0.
886 pevap1d(l)=0.
887 vsnow1d(l)=0.
888 vrain11d(l)=0.
889 vrain21d(l)=0.
890 vci1d(l)=0.
891 nsmice1d(l)=0.
892 dbz_col(l)=dbzmin
893 nr_col(l)=0.
894 ns_col(l)=0.
895 indexr1d(l)=0.
896 indexs1d(l)=0.
897 rflag1d(l)=0. !jun01
898!
899!--- Skip this level and go to the next lower level if no condensate
900! and very low specific humidities
901!
902!--- Check if any rain is falling into layer from above
903!
904 IF (arain .GT. climit) THEN
905 clear=.false.
906 vrain1=0.
907 ELSE
908 clear=.true.
909 arain=0.
910 ENDIF
911!
912!--- Check if any ice is falling into layer from above
913!
914!--- NOTE that "SNOW" in variable names is often synonomous with
915! large, precipitation ice particles
916!
917 IF (asnow .GT. climit) THEN
918 clear=.false.
919 vsnow1=0.
920 ELSE
921 asnow=0.
922 ENDIF
923!
924!-----------------------------------------------------------------------
925!------------ Proceed with cloud microphysics calculations -------------
926!-----------------------------------------------------------------------
927!
928 tk=t_col(l) ! Temperature (deg K)
929 tc=tk-t0c ! Temperature (deg C)
930 pp=p_col(l) ! Pressure (Pa)
931 q=q_col(l) ! Specific humidity of water vapor (kg/kg)
932 wv=q/(1.-q) ! Water vapor mixing ratio (kg/kg)
933 wc=wc_col(l) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
934 rhgrd=rhc_col(l) ! Threshold relative humidity for the onset of condensation
935!
936!-----------------------------------------------------------------------
937!--- Moisture variables below are mixing ratios & not specifc humidities
938!-----------------------------------------------------------------------
939!
940!--- This check is to determine grid-scale saturation when no condensate is present
941!
942 esw=min(1000.*fpvs0(tk),0.99*pp) ! Saturation vapor pressure w/r/t water
943 qsw=eps*esw/(pp-esw) ! Saturation mixing ratio w/r/t water
944 ws=qsw ! General saturation mixing ratio (water/ice)
945 qsi=qsw ! Saturation mixing ratio w/r/t ice
946 IF (tc .LT. 0.) THEN
947 esi=min(1000.*fpvs(tk),0.99*pp) ! Saturation vapor pressure w/r/t ice
948 qsi=eps*esi/(pp-esi) ! Saturation mixing ratio w/r/t water
949 ws=qsi ! General saturation mixing ratio (water/ice)
950 ENDIF
951!
952!--- Effective grid-scale Saturation mixing ratios
953!
954 qswgrd=rhgrd*qsw
955 qsigrd=rhgrd*qsi
956 wsgrd=rhgrd*ws
957!
958!--- Check if air is subsaturated and w/o condensate
959!
960 IF (wv.GT.wsgrd .OR. wc.GT.epsq) clear=.false.
961!
962!-----------------------------------------------------------------------
963!-- Loop to the end if in clear, subsaturated air free of condensate ---
964!-----------------------------------------------------------------------
965!
966 IF (clear) THEN
967 strat=.false. !- Reset stratiform rain flag
968 drzl=.false. !- Reset drizzle flag
969 indexrmin=mdrmin !- Reset INDEXRmin
970 timlt=0. !- Reset accumulated ice melting
971 cycle big_loop
972 ENDIF
973!
974!-----------------------------------------------------------------------
975!--------- Initialize RHO, THICK & microphysical processes -------------
976!-----------------------------------------------------------------------
977!
978!
979!--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
980! (see pp. 63-65 in Fleagle & Businger, 1963)
981!
982 rho=pp/(rd*tk*(1.+eps1*q)) ! Air density (kg/m**3)
983 rrho=1./rho ! Reciprocal of air density
984 dtrho=dtph*rho ! Time step * air density
985 bldtrh=blend*dtrho ! Blend parameter * time step * air density
986 thick=thick_col(l) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
987!
988 arainnew=0. ! Updated accumulated rainfall
989 asnownew=0. ! Updated accumulated snowfall
990 qi=qi_col(l) ! Ice mixing ratio
991 qinew=0. ! Updated ice mixing ratio
992 qr=qr_col(l) ! Rain mixing ratio
993 qrnew=0. ! Updated rain ratio
994 qw=qw_col(l) ! Cloud water mixing ratio
995 qwnew=0. ! Updated cloud water ratio
996!
997 pcond=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
998 pidep=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
999 pinit=0. ! Ice initiation (part of PIDEP calculation, kg/kg)
1000 piacw=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
1001 piacwi=0. ! Growth of precip ice by riming (kg/kg; >0)
1002 piacwr=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
1003 piacr=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
1004 picnd=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
1005 pievp=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
1006 pimlt=0. ! Melting ice (kg/kg; >0)
1007 praut=0. ! Cloud water autoconversion to rain (kg/kg; >0)
1008 pracw=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
1009 prevp=0. ! Rain evaporation (kg/kg; <0)
1010 nsmice=0. ! Cloud ice number concentration (m^-3)
1011 nrain=0. ! Rain number concentration (m^-3) !jul28 begin
1012 nsnow=0. ! "Snow" number concentration (m^-3)
1013 rqrnew=0. ! Final rain content (kg/m**3)
1014 rqsnew=0. ! Final "snow" content (kg/m**3)
1015 zrain=0. ! Radar reflectivity from rain (mm**6 m-3)
1016 zsnow=0. ! Radar reflectivity from snow (mm**6 m-3)
1017 ztot=0. ! Radar reflectivity from rain+snow (mm**6 m-3)
1018 indexr=mdrmin ! Mean diameter of rain (microns)
1019 indexr1=indexr ! 1st updated mean diameter of rain (microns)
1020 indexr2=indexr ! 2nd updated mean diameter of rain (microns)
1021 n0r=0. ! 1st estimate for rain intercept (m^-4)
1022 dum1=min(0.,tc)
1023 dum=xmimax*exp(xmiexp*dum1)
1024 indexs=min(mdimax, max(mdimin, int(dum) ) ) ! 1st estimate for mean diameter of snow (microns)
1025 vci=0. ! Cloud ice fall speeds (m/s)
1026 vsnow=0. ! "Snow" (snow/graupel/sleet/hail) fall speeds (m/s)
1027 vrain2=0. ! Rain fall speeds out of bottom of grid box (m/s)
1028 rimef1=1. ! Rime Factor (ratio, >=1, defined below)
1029!
1030!--- Double check input hydrometeor mixing ratios
1031!
1032! DUM=WC-(QI+QW+QR)
1033! DUM1=ABS(DUM)
1034! DUM2=TOLER*MIN(WC, QI+QW+QR)
1035! IF (DUM1 .GT. DUM2) THEN
1036! WRITE(0,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
1037! & ' L=',L
1038! WRITE(0,"(4(a12,g11.4,1x))")
1039! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
1040! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
1041! ENDIF
1042!
1043!***********************************************************************
1044!*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
1045!***********************************************************************
1046!
1047!--- Calculate a few variables, which are used more than once below
1048!
1049!--- Latent heat of vaporization as a function of temperature from
1050! Bolton (1980, JAS)
1051!
1052 tk2=1./(tk*tk) ! 1./TK**2
1053!
1054!--- Basic thermodynamic quantities
1055! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
1056! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
1057! * DIFFUS - diffusivity of water vapor [ m**2/s ]
1058!
1059 tfactor=sqrt(tk*tk*tk)/(tk+120.)
1060 dynvis=1.496e-6*tfactor
1061 therm_cond=2.116e-3*tfactor
1062 diffus=8.794e-5*tk**1.81/pp
1063!
1064!--- Air resistance term for the fall speed of ice following the
1065! basic research by Heymsfield, Kajikawa, others
1066!
1067 gammas=min(1.5, (1.e5/pp)**c1) !-- limited to 1.5x
1068!
1069!--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
1070!
1071 gammar=(rho0/rho)**.4
1072!
1073!----------------------------------------------------------------------
1074!------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
1075!----------------------------------------------------------------------
1076!
1077!--- Determine if conditions supporting ice are present
1078!
1079 IF (tc.LT.0. .OR. qi.GT. epsq .OR. asnow.GT.climit) THEN
1080 ice_logical=.true.
1081 ELSE
1082 ice_logical=.false.
1083 qlice=0.
1084 qtice=0.
1085 ENDIF
1086 IF (t_ice <= -100.) THEN
1087 ice_logical=.false.
1088 qlice=0.
1089 qtice=0.
1090 ENDIF
1091!
1092!--- Determine if rain is present
1093!
1094 rain_logical=.false.
1095 IF (arain.GT.climit .OR. qr.GT.epsq) rain_logical=.true.
1096!
1097ice_test: IF (ice_logical) THEN
1098!
1099!--- IMPORTANT: Estimate time-averaged properties.
1100!
1101!---
1102! -> Small ice particles are assumed to have a mean diameter of 50 microns.
1103! * QSmICE - estimated mixing ratio for small cloud ice
1104!---
1105! * TOT_ICE - total mass (small & large) ice before microphysics,
1106! which is the sum of the total mass of large ice in the
1107! current layer and the input flux of ice from above
1108! * PILOSS - greatest loss (<0) of total (small & large) ice by
1109! sublimation, removing all of the ice falling from above
1110! and the ice within the layer
1111! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
1112! ice mass to the unrimed ice mass (>=1)
1113! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
1114! * VSNOW - Fall speed of rimed snow w/ air resistance correction
1115! * VCI - Fall speed of 50-micron ice crystals w/ air resistance correction
1116! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
1117! * XLIMASS - used for debugging, associated with calculating large ice mixing ratio
1118! * FLIMASS - mass fraction of large ice
1119! * QTICE - time-averaged mixing ratio of total ice
1120! * QLICE - time-averaged mixing ratio of large ice
1121! * NLICE - time-averaged number concentration of large ice
1122! * NSmICE - number concentration of small ice crystals at current level
1123! * QSmICE - mixing ratio of small ice crystals at current level
1124!---
1125!--- Assumed number fraction of large ice particles to total (large & small)
1126! ice particles, which is based on a general impression of the literature.
1127!
1128 ninuclei=0.
1129 nsmice=0.
1130 qsmice=0.
1131 rcw=0.
1132 IF (tc<0.) THEN
1133!
1134!--- Max # conc of small ice crystals based on 10% of total ice content
1135! or the parameter NSI_max
1136!
1137 nsimax=max(nsi_max, 0.1*rho*qi/massi(mdimin) ) !aug27
1138!
1139!-- Specify Fletcher, Cooper, Meyers, etc. here for ice nuclei concentrations
1140! Cooper (1986): NInuclei=MIN(5.*EXP(-0.304*TC), NSImax)
1141! Fletcher (1962): NInuclei=MIN(0.01*EXP(-0.6*TC), NSImax)
1142!
1143!aug28: The formulas below mean that Fletcher is used for >-21C and Cooper at colder
1144! temperatures. In areas of high ice contents near the tops of deep convection,
1145! the number concentrations will be determined by the lower value of the "FQi"
1146! contribution to NSImax or Cooper.
1147!
1148 ninuclei=min(0.01*exp(-0.6*tc), nsimax) !aug28 - Fletcher (1962)
1149 ninuclei=min(5.*exp(-0.304*tc), ninuclei) !aug28 - Cooper (1984)
1150 IF (qi>epsq) THEN
1151 dum=rrho*massi(mdimin)
1152 nsmice=min(ninuclei, qi/dum)
1153 qsmice=nsmice*dum
1154 ENDIF ! End IF (QI>EPSQ)
1155 ENDIF ! End IF (TC<0.)
1156 init_ice: IF (qi<=epsq .AND. asnow<=climit) THEN
1157 tot_ice=0.
1158 piloss=0.
1159 rimef1=1.
1160 vrimef=1.
1161 vel_inc=gammas
1162 vsnow=0.
1163 vsnow1=0.
1164 vci=0.
1165 emairi=thick
1166 xlimass=rimef1*massi(indexs)
1167 flimass=1.
1168 qlice=0.
1169 rqlice=0.
1170 qtice=0.
1171 nlice=0.
1172 ELSE init_ice
1173 !
1174 !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
1175 ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
1176 ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
1177 !
1178!
1179!sep10 - Start of changes described in (23) at top of code.
1180!
1181 tot_ice=thick*qi+blend*asnow
1182 piloss=-tot_ice/thick
1183 qlgice=max(0., qi-qsmice) !-- 1st-guess estimate of large ice
1184 vci=gammas*vsnowi(mdimin)
1185!
1186!-- Need to save this original value before two_pass iteration
1187!
1188 lbef=max(1,l-1)
1189 rimef1=(rimef_col(l)*thick*qlgice &
1190 & +rimef_col(lbef)*blend*asnow)/tot_ice
1191!
1192!mar08 see (32); !apr22a see (41) start - Estimate mean diameter (INDEXS in microns)
1193 IF (rimef1>2.) THEN
1194 dum3=rh_ngc*(rho*qlgice)**c1 !- convective mean diameter
1195 dum2=rh_ngt*(rho*qlgice)**c1 !- transition mean diameter
1196 IF (rimef1>=10.) THEN
1197 dum=dum3
1198 ELSE IF (rimef1>=5.) THEN
1199 dum=0.2*(rimef1-5.) !- Blend at 5<=RF<10
1200 dum=dum3*dum+dum2*(1.-dum)
1201 ELSE
1202 dum1=real(indexs) !- stratiform mean diameter
1203 dum=0.33333*(rimef1-2.) !- Blend at 2<RF<5
1204 dum=dum2*dum+dum1*(1.-dum)
1205 ENDIF
1206 indexs=min(mdimax, max(mdimin, int(dum) ) )
1207 ENDIF
1208!may11 - begin
1209 emairi=thick+bldtrh*vsnow1
1210 qlice=(thick*qlgice+blend*asnow)/emairi !- Estimate of large ice
1211!may17 - start; now calculated before the DO IPASS iteration
1212 rqlice=rho*qlice
1213 qtice=qlice+qsmice
1214 flimass=qlice/qtice
1215!may17 end
1216!may11 - end
1217!mar08 !apr22a end
1218!
1219!===========================================
1220 two_pass: DO ipass=1,2
1221!===========================================
1222!
1223!-- Prevent rime factor (RimeF1) from exceeding a maximum value, RFmx, in which
1224! the ice has an equivalent density near that of pure ice
1225!
1226 dum=1.e-6*real(indexs) !- Mean diameter in m
1227 rfmx=rfmx1*dum*dum*dum/massi(indexs)
1228 rimef1=min(rimef1, rfmx)
1229!
1230 vel_rime: IF (rimef1<=1.) THEN
1231 rimef1=1.
1232 vrimef=1.
1233 ELSE vel_rime
1234 !--- Prevent rime factor (RimeF1) from exceeding a maximum value (RFmax)
1235 rimef1=min(rimef1, rfmax)
1236 ixs=max(2, min(indexs/100, 9))
1237 xrf=10.492*alog(rimef1)
1238 ixrf=max(0, min(int(xrf), nrime))
1239 IF (ixrf .GE. nrime) THEN
1240 vrimef=vel_rf(ixs,nrime)
1241 ELSE
1242 vrimef=vel_rf(ixs,ixrf)+(xrf-float(ixrf))* &
1243 & (vel_rf(ixs,ixrf+1)-vel_rf(ixs,ixrf))
1244 ENDIF
1245 vrimef=max(1., vrimef)
1246 ENDIF vel_rime
1247!
1248!may17 - start
1249 vel_inc=gammas*vrimef !- Normal rimed ice fall speeds
1250 vsnow=vel_inc*vsnowi(indexs)
1251 IF (rimef1>=5. .AND. indexs==mdimax .AND. rqlice>rqhail) THEN
1252!- Additional increase using Thompson graupel/hail fall speeds
1253 dum=gammas*avhail*rqlice**bvhail
1254 IF (dum>vsnow) THEN
1255 vsnow=dum
1256 vel_inc=vsnow/vsnowi(indexs)
1257 ENDIF
1258 ENDIF
1259 xlimass=rimef1*massi(indexs)
1260 nlice=rqlice/xlimass
1261!
1262!sep16 - End of change described in (26) at top of code.
1263!
1264!===========================================
1265 IF (ipass>=2 .OR. &
1266 (nlice>=nlimin .AND. nlice<=nsi_max)) EXIT two_pass
1267!may17 - end
1268!===========================================
1269!
1270!--- Force NLICE to be between NLImin and NSI_max when IPASS=1
1271!
1272! IF (PRINT_diag .AND. RQLICE>Revised_LICE) DUM2=NLICE !-- For debugging (see DUM2 below)
1273 nlice=max(nlimin, min(nsi_max, nlice) )
1274!sep16 - End of changes
1275!
1276 xli=rqlice/(nlice*rimef1) !- Mean mass of unrimed ice
1277new_size: IF (xli<=massi(mdimin) ) THEN
1278 indexs=mdimin
1279 ELSE IF (xli<=massi(450) ) THEN new_size
1280 dli=9.5885e5*xli**.42066 ! DLI in microns
1281 indexs=min(mdimax, max(mdimin, int(dli) ) )
1282 ELSE IF (xli<massi(mdimax) ) THEN new_size
1283 dli=3.9751e6*xli**.49870 ! DLI in microns
1284 indexs=min(mdimax, max(mdimin, int(dli) ) )
1285 ELSE new_size
1286 indexs=mdimax
1287! IF (PRINT_diag .AND. RQLICE>Revised_LICE) THEN
1288! WRITE(0,"(5(a12,g11.4,1x))") '{$ RimeF1=',RimeF1, &
1289! & ' RHO*QLICE=',RQLICE,' TC=',TC,' NLICE=',NLICE, &
1290! & ' NLICEold=',DUM2
1291! Revised_LICE=1.2*RQLICE
1292! ENDIF
1293 ENDIF new_size
1294!===========================================
1295 ENDDO two_pass
1296!===========================================
1297 ENDIF init_ice
1298 ENDIF ice_test
1299!
1300!mar03 !may11 !jun01 - start; see (34) above
1301 IF(tc<=0.) THEN
1302 strat=.false.
1303 indexrmin=mdrmin
1304 timlt=0.
1305 indexs0c=indexs
1306 rhox0c=22.5*max(1.,min(rimef1,40.)) !- Estimated ice density at 0C (kg m^-3)
1307 ELSE ! TC>0.
1308 IF(.NOT.rain_logical) THEN
1309 strat=.false. !- Reset STRAT
1310 indexrmin=mdrmin !- Reset INDEXRmin
1311 IF(.NOT.ice_logical) timlt=0.
1312 ELSE
1313!- STRAT=T for stratiform rain
1314 IF(timlt>epsq .AND. rhox0c<=225.) THEN
1315 strat=.true.
1316 ELSE
1317 strat=.false. !- Reset STRAT
1318 indexrmin=mdrmin !- Reset INDEXRmin
1319 ENDIF
1320 IF(strat .AND. indexrmin<=mdrmin) THEN
1321 indexrmin=indexs0c*(0.001*rhox0c)**c1
1322 indexrmin=max(mdrmin, min(indexrmin, indexrstrmax) )
1323 ENDIF
1324 ENDIF
1325 ENDIF
1326!
1327 IF(strat .OR. timlt>epsq) THEN
1328 drzl=.false.
1329 ELSE
1330!- DRZL=T for drizzle (no melted ice falling from above)
1331 drzl=.true. !mar30
1332 ENDIF
1333!jun01 - end
1334!
1335!----------------------------------------------------------------------
1336!--------------- Calculate individual processes -----------------------
1337!----------------------------------------------------------------------
1338!
1339!--- Cloud water autoconversion to rain (PRAUT) and collection of cloud
1340! water by precipitation ice (PIACW)
1341!
1342 IF (qw.GT.epsq .AND. tc.GE.t_ice) THEN
1343!-- The old autoconversion threshold returns
1344 dum2=rho*qw
1345 IF (dum2>qaut0) THEN
1346!-- July 2010 version follows Liu & Daum (JAS, 2004) and Liu et al. (JAS, 2006)
1347 dum2=dum2*dum2
1348 dum=braut*dum2*qw
1349 dum1=araut*dum2
1350 praut=min(qw, dum*(1.-exp(-dum1*dum1)) )
1351 ENDIF
1352 IF (qlice .GT. epsq) THEN
1353 !
1354 !--- Collection of cloud water by large ice particles ("snow")
1355 ! PIACWI=PIACW for riming, PIACWI=0 for shedding
1356 !
1357 fws=min(1., ciacw*vel_inc*nlice*accri(indexs) ) !jul28 (16)
1358 piacw=fws*qw
1359 IF (tc<0.) THEN
1360 piacwi=piacw !- Large ice riming
1361 rcw=arcw*dum2**c1 !- Cloud droplet radius in microns
1362 ENDIF
1363 ENDIF ! End IF (QLICE .GT. EPSQ)
1364 ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
1365!
1366!----------------------------------------------------------------------
1367!--- Calculate homogeneous freezing of cloud water (PIACW, PIACWI) and
1368! ice deposition (PIDEP), which also includes ice initiation (PINIT)
1369!
1370ice_only: IF (tc.LT.t_ice .AND. (wv.GT.qswgrd .OR. qw.GT.epsq)) THEN
1371 !
1372 !--- Adjust to ice saturation at T<T_ICE (-40C) if saturated w/r/t water
1373 ! or if cloud water is present (homogeneous glaciation).
1374 !
1375 piacw=qw
1376 piacwi=piacw
1377 rcw=0. ! Homogeneous freezing of cloud water adds to
1378 ! cloud ice, do not use to update RF for large ice
1379 dum1=tk+xlf1*piacw ! Updated (dummy) temperature (deg K)
1380 dum2=wv ! Updated (dummy) water vapor mixing ratio
1381 dum=min(1000.*fpvs(dum1),0.99*pp) ! Updated (dummy) saturation vapor pressure w/r/t ice
1382 dum=rhgrd*eps*dum/(pp-dum) ! Updated (dummy) saturation mixing ratio w/r/t ice
1383
1384!-- Debug 20120111
1385IF (warn1 .AND. dum1<xmin) THEN
1386 WRITE(0,*) 'WARN1: Water saturation T<180K; I,J,L,TC,P=', &
1387 i_index,j_index,l,dum1-t0c,.01*pp
1388 warn1=.false.
1389ENDIF
1390
1391!-- Adjust to ice saturation if homogeneous freezing occurs
1392 IF (dum2>dum) &
1393 & pidep=deposit(pp,dum1,dum2,rhgrd,i_index,j_index,l)
1394
1395 dwvi=0. ! Used only for debugging
1396 !
1397 ELSE IF (tc<0.) THEN ice_only
1398 !
1399 !--- These quantities are handy for ice deposition/sublimation
1400 ! PIDEP_max - max deposition or minimum sublimation to ice saturation
1401 !
1402 denomi=1.+xls3*qsi*tk2
1403! DWVi=MIN(WV,QSW)-QSIgrd !may17
1404 dwvi=wv-qsigrd !may17 !-- Speed up ice deposition
1405 pidep_max=max(piloss, dwvi/denomi)
1406 didep=0. !-- Vapor deposition/sublimation onto existing ice
1407 IF (qtice .GT. 0.) THEN
1408 !
1409 !--- Calculate ice deposition/sublimation
1410 ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1411 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1412 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1413 ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
1414 ! VENTI2 - m**2/s**.5 ; DIDEP - unitless
1415 !
1416 abi=1./(rho*xls2*qsi*tk2/therm_cond+1./diffus)
1417 !
1418 !--- VENTIL - Number concentration * ventilation factors for large ice
1419 !--- VENTIS - Number concentration * ventilation factors for small ice
1420 !
1421 !--- Variation in the number concentration of ice with time is not
1422 ! accounted for in these calculations (could be in the future).
1423 !
1424 dum=(rho/(diffus*diffus*dynvis))**c2
1425 sfactor=sqrt(gammas)*dum !-- GAMMAS (RF=1 for cloud ice)
1426 ventis=(venti1(mdimin)+sfactor*venti2(mdimin))*nsmice
1427 sfactor=sqrt(vel_inc)*dum
1428 ventil=(venti1(indexs)+sfactor*venti2(indexs))*nlice
1429 didep=abi*(ventil+ventis)*dtph
1430 ENDIF !-- IF (QTICE .GT. 0.) THEN
1431 !
1432 !--- WARNING: the Meyers et al. stuff is not used!
1433 !--- Two modes of ice nucleation from Meyers et al. (JAM, 1992):
1434 ! (1) Deposition & condensation freezing nucleation - eq. (2.4),
1435 ! requires ice supersaturation
1436 ! (2) Contact freezing nucleation - eq. (2.6), requires presence
1437 ! of cloud water
1438 !
1439 !--- Ice crystal growth during the physics time step is calculated using
1440 ! Miller & Young (1979, JAS) and is represented by MY_GROWTH(INDEX_MY),
1441 ! where INDEX_MY is tabulated air temperatures between -1 and -35 C.
1442 ! The original Miller & Young (MY) calculations only went down to -30C,
1443 ! so a fixed value is assumed at temperatures colder than -30C.
1444 !
1445 !--- Sensitivity test using:
1446 ! - Fletcher (1962) for ice initiation
1447 ! - Can occur only when there is water supersaturation and T<-12C.
1448 ! - Maximum cloud ice number concentration of 100 per liter
1449 !
1450 IF (wv>qswgrd .AND. tc<t_ice_init .AND. nsmice<ninuclei) THEN
1451 index_my=max(my_t1, min( int(.5-tc), my_t2 ) )
1452 !-- Only initiate small ice to get to NInuclei number concentrations
1453 dum=max(0., ninuclei-nsmice)
1454 pinit=max(0., dum*my_growth_nmm(index_my)*rrho)
1455 ENDIF
1456 !
1457 !--- Calculate PIDEP, but also account for limited water vapor supply
1458 !
1459 IF (dwvi>0.) THEN
1460 pidep=min(dwvi*didep+pinit, pidep_max)
1461 ELSE IF (dwvi<0.) THEN
1462 pidep=max(dwvi*didep, pidep_max)
1463 ENDIF
1464 !
1465 ENDIF ice_only
1466!
1467!------------------------------------------------------------------------
1468!--- Cloud water condensation
1469!
1470 IF (tc>=t_ice .AND. (qw>epsq .OR. wv>qswgrd)) THEN
1471!apr21 - start; see (39) above
1472 dum1=qw-piacwi
1473 dum2=tk+xls1*pidep+xlf1*piacwi
1474 dum3=wv-pidep
1475 pcond=condense(pp,dum1,dum2,dum3,rhgrd,i_index,j_index,l)
1476!apr21 - end; see (39) above
1477 ENDIF
1478!
1479!--- Limit freezing of accreted rime to prevent temperature oscillations,
1480! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
1481!
1482 tcc=tc+xlv1*pcond+xls1*pidep+xlf1*piacwi
1483 IF (tcc>0.) THEN
1484 piacwi=0.
1485 rcw=0. !apr18 see (38)
1486 pidep=0. !apr18 see (38)
1487 tcc=tc+xlv1*pcond !apr18 see (38)
1488 ENDIF
1489!
1490 qsw0=0.
1491 dwv0=0.
1492 IF (tc.GT.0. .AND. tcc.GT.0. .AND. ice_logical) THEN
1493 !
1494 !--- Calculate melting and evaporation/condensation
1495 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1496 ! VENTIL - m**-2 ; VENTI1 - m ;
1497 ! VENTI2 - m**2/s**.5 ; CIEVP - /s
1498 !
1499 sfactor=sqrt(vel_inc)*(rho/(diffus*diffus*dynvis))**c2
1500 ventil=nlice*(venti1(indexs)+sfactor*venti2(indexs))
1501 aievp=ventil*diffus*dtph
1502 IF (aievp .LT. xratio) THEN
1503 dievp=aievp
1504 ELSE
1505 dievp=1.-exp(-aievp)
1506 ENDIF
1507 dum=qw+pcond
1508 IF (wv.LT.qsw .AND. dum.LE.epsq) THEN
1509 !
1510 !--- Evaporation from melting snow (sink of snow) or shedding
1511 ! of water condensed onto melting snow (source of rain)
1512 !
1513 dum=min(esw0, 0.99*pp) !- Limit on ESW0 at low pressures
1514 qsw0=max(epsq, eps*dum/(pp-dum) ) !- Constrain QSW0 to be >=EPSQ
1515 dwv0=min(wv,qsw)-qsw0
1516 dum=dwv0*dievp
1517 pievp=max( min(0., dum), piloss)
1518 picnd=max(0., dum)
1519 ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
1520 pimlt=therm_cond*tcc*ventil*rrho*dtph/xlf
1521 !
1522 !--- Limit melting to prevent temperature oscillations across 0C
1523 !
1524 dum1=max( 0., (tcc+xlv1*pievp)/xlf1 )
1525 pimlt=min(pimlt, dum1)
1526 !
1527 !--- Limit loss of snow by melting (>0) and evaporation
1528 !
1529 dum=pievp-pimlt
1530 IF (dum .LT. piloss) THEN
1531 dum1=piloss/dum
1532 pimlt=pimlt*dum1
1533 pievp=pievp*dum1
1534 ENDIF ! End IF (DUM .GT. QTICE)
1535 ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
1536!
1537!--- IMPORTANT: Estimate time-averaged properties.
1538!
1539! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
1540! the total mass of rain in the current layer and the input
1541! flux of rain from above
1542! * VRAIN1 - fall speed of rain into grid from above
1543! * VRAIN2 - fall speed of rain out of grid box to the level below
1544! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
1545! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
1546! above and the rain within the layer
1547! * RQR - rain content (kg/m**3)
1548! * INDEXR - mean size of rain drops to the nearest 1 micron in size
1549! * N0r - intercept of rain size distribution (typically 10**6 m**-4)
1550!
1551 tot_rain=0.
1552 qtrain=0.
1553 prloss=0.
1554 rqr=0.
1555do_rain: IF (rain_logical) THEN
1556 tot_rain=thick*qr+blend*arain !aug27 begin mods
1557 prloss=-tot_rain/thick
1558!may11 - start
1559 qtrain=tot_rain/(thick+bldtrh*vrain1) !- Rain mixing ratio
1560 rqr=rho*qtrain !- Rain content
1561 IF(strat .AND. rqr>1.e-3) strat=.false. !may31
1562 IF(drzl .AND. pimlt>epsq) drzl=.false. !jun01
1563!
1564dsd1: IF (rqr<=rqr_drmin) THEN
1565!-- Extremely light rain
1566 indexr=mdrmin
1567 n0r=rqr/massr(indexr)
1568 ELSE IF (drzl) THEN dsd1
1569!-- Drizzle - gradually increase N0r for rain contents <0.5 g m^-3 !may31
1570 dum=max(1.0, 0.5e-3/rqr)
1571 n0r=min(1.e9, n0r0*dum*sqrt(dum) ) !- 8.e6 <= N0r <= 1.e9
1572 dum1=rqr/(pi*rhol*n0r)
1573 dum=1.e6*sqrt(sqrt(dum1))
1574 indexr=max(xmrmin, min(dum, xmrmax) )
1575 n0r=rqr/massr(indexr)
1576!may20 - start
1577 ELSE IF (rqr>=rqr_drmax) THEN dsd1
1578!-- Extremely heavy rain
1579 indexr=mdrmax
1580 n0r=rqr/massr(indexr)
1581 ELSE dsd1
1582!-- Light to heavy rain
1583 n0r=n0r0
1584 dum=cn0r0*sqrt(sqrt(rqr))
1585 indexr=max(xmrmin, min(dum, xmrmax) )
1586 IF (strat .AND. indexr<indexrmin) THEN
1587!-- Stratiform rain
1588 indexr=indexrmin
1589 n0r=rqr/massr(indexr)
1590 ENDIF
1591 ENDIF dsd1
1592!may20 - end
1593 vrain2=gammar*vrain(indexr) !-- VRAIN2 is used below
1594!may11 - end
1595!
1596 IF (tc .LT. t_ice) THEN
1597 piacr=-prloss
1598 ELSE
1599 dwvr=wv-pcond-qsw
1600 dum=qw+pcond
1601 IF (dwvr.LT.0. .AND. dum.LE.epsq) THEN
1602 !
1603 !--- Rain evaporation
1604 !
1605 ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1606 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1607 !
1608 ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ;
1609 ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ;
1610 ! CREVP - unitless
1611 !
1612 rfactor=sqrt(gammar)*(rho/(diffus*diffus*dynvis))**c2
1613 abw=1./(rho*xlv2*qsw*tk2/therm_cond+1./diffus)
1614 !
1615 !--- Note that VENTR1, VENTR2 lookup tables do not include the
1616 ! 1/Davg multiplier as in the ice tables
1617 !
1618 ventr=n0r*(ventr1(indexr)+rfactor*ventr2(indexr))
1619 crevp=abw*ventr*dtph
1620 prevp=max(dwvr*crevp, prloss)
1621 ELSE IF (qw .GT. epsq) THEN
1622 fwr=cracw*gammar*n0r*accrr(indexr)
1623 pracw=min(1.,fwr)*qw
1624 ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
1625 !
1626 IF (ice_logical .AND. tc<0. .AND. tcc<0.) THEN
1627 !
1628 !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
1629 ! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
1630 !
1631 dum=.001*float(indexr)
1632 dum=(exp(abfr*tc)-1.)*dum*dum*dum*dum*dum*dum*dum
1633 piacr=min(cbfr*n0r*rrho*dum, qtrain)
1634!- Start changes listed above in (33) - no collisional freezing of rain
1635 fir=0.
1636
1637 IF (qlice .GT. epsq) THEN
1638 !
1639 !--- Freezing of rain by collisions w/ large ice
1640 !
1641 dum=gammar*vrain(indexr)
1642 dum1=dum-vsnow
1643 !
1644 !--- DUM2 - Difference in spectral fall speeds of rain and
1645 ! large ice, parameterized following eq. (48) on p. 112 of
1646 ! Murakami (J. Meteor. Soc. Japan, 1990)
1647 !
1648 dum2=sqrt(dum1*dum1+.04*dum*vsnow)
1649 dum1=5.e-12*indexr*indexr+2.e-12*indexr*indexs &
1650 & +.5e-12*indexs*indexs
1651 fir=min(1., ciacr*nlice*dum1*dum2)
1652 !
1653 !--- Future? Should COLLECTION BY SMALL ICE BE INCLUDED???
1654 !
1655 piacr=min(piacr+fir*qtrain, qtrain)
1656 ENDIF ! End IF (QLICE .GT. EPSQ)
1657 dum=prevp-piacr
1658 If (dum .LT. prloss) THEN
1659 dum1=prloss/dum
1660 prevp=dum1*prevp
1661 piacr=dum1*piacr
1662 ENDIF ! End If (DUM .LT. PRLOSS)
1663 ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
1664 ENDIF ! End IF (TC .LT. T_ICE)
1665 ENDIF do_rain ! End IF (RAIN_logical)
1666!
1667 indexr1=indexr
1668!
1669!----------------------------------------------------------------------
1670!---------------------- Main Budget Equations -------------------------
1671!----------------------------------------------------------------------
1672!
1673!
1674!-----------------------------------------------------------------------
1675!--- Update fields, determine characteristics for next lower layer ----
1676!-----------------------------------------------------------------------
1677!
1678!--- Carefully limit sinks of cloud water
1679!
1680 dum1=piacw+praut+pracw-min(0.,pcond)
1681 IF (dum1 .GT. qw) THEN
1682 dum=qw/dum1
1683 piacw=dum*piacw
1684 piacwi=dum*piacwi
1685 praut=dum*praut
1686 pracw=dum*pracw
1687 IF (pcond .LT. 0.) pcond=dum*pcond
1688 ENDIF
1689 piacwr=piacw-piacwi ! TC >= 0C
1690!
1691!--- QWnew - updated cloud water mixing ratio
1692!
1693 delw=pcond-piacw-praut-pracw
1694 qwnew=qw+delw
1695 IF (qwnew .LE. epsq) qwnew=0.
1696 IF (qw.GT.0. .AND. qwnew.NE.0.) THEN
1697 dum=qwnew/qw
1698 IF (dum .LT. toler) qwnew=0.
1699 ENDIF
1700!
1701!--- Update temperature and water vapor mixing ratios
1702!
1703 delt= xlv1*(pcond+pievp+picnd+prevp) &
1704 & +xls1*pidep+xlf1*(piacwi+piacr-pimlt)
1705 tnew=tk+delt
1706!
1707 delv=-pcond-pidep-pievp-picnd-prevp
1708 wvnew=wv+delv
1709!
1710!--- Update ice mixing ratios
1711!
1712!---
1713! * TOT_ICEnew - total mass (small & large) ice after microphysics,
1714! which is the sum of the total mass of ice in the
1715! layer and the flux of ice out of the grid box below
1716! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
1717! rimed) ice mass to the unrimed ice mass (>=1)
1718! * QInew - updated mixing ratio of total (large & small) ice in layer
1719! * QLICEnew=FLIMASS*QInew, an estimate of the updated large ice mixing ratio
1720!
1721! -> TOT_ICEnew=QInew*THICK+QLICEnew*BLDTRH*VSNOW+(QInew-QLICEnew)*BLDTRH*VCI
1722! =QInew*THICK+QInew*FLIMASS*BLDTRH*VSNOW+QInew*(1.-FLIMASS)*BLDTRH*VCI
1723! =QInew*THICK+QInew*BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI)
1724! =QInew*(THICK+BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI))
1725! -> Rearranging this equation gives:
1726! QInew=TOT_ICEnew/(THICK+BLDTRH*(FLIMASS*VSNOW+(1.-FLIMASS)*VCI))
1727!
1728! * ASNOWnew - updated accumulation of snow and cloud ice at bottom of grid cell
1729! -> ASNOWnew=QInew*BLDTRH*(FLIMASS+VSNOW+(1.-FLIMASS)*VCI)
1730!---
1731!
1732 deli=0.
1733 rimef=1.
1734 IF (ice_logical) THEN
1735 deli=pidep+pievp+piacwi+piacr-pimlt
1736 tot_icenew=tot_ice+thick*deli
1737 IF (tot_ice.GT.0. .AND. tot_icenew.NE.0.) THEN
1738 dum=tot_icenew/tot_ice
1739 IF (dum .LT. toler) tot_icenew=0.
1740 ENDIF
1741 IF (tot_icenew .LE. climit) THEN
1742 tot_icenew=0.
1743 rimef=1.
1744 qinew=0.
1745 asnownew=0.
1746 ELSE
1747!
1748!--- Update rime factor if appropriate
1749!
1750!apr22 - start; see (40) above
1751 IF (pinit<=epsq) THEN
1752 psdep=max(0., pidep)
1753 ELSE
1754!-- Assume vapor deposition is onto cloud ice if PINT>0
1755 psdep=0.
1756 ENDIF
1757 dels=piacwi+piacr+psdep
1758 IF (dels<=epsq .OR. tnew>=t0c) THEN
1759 rimef=rimef1
1760 ELSE
1761!
1762!-----------------------------------------------------------------------
1763! RimeF is the new, updated rime factor; RimeF1 is the existing rime factor
1764!
1765! From near line 1115:
1766! RimeF1=(RimeF_col(L)*THICK*QLgICE+RimeF_col(LBEF)*BLEND*ASNOW)/TOT_ICE
1767! RimeF1*TOT_ICE=RimeF_col(L)*THICK*QLgICE+RimeF_col(LBEF)*BLEND*ASNOW
1768!
1769! TOT_ICEnew=TOT_ICE+THICK*DELI
1770! TOT_ICEnew=TOT_ICE+THICK*(PIDEP+PIEVP+PIACWI+PIACR-PIMLT)
1771!
1772! But the following processes do not affect the rime factor (ice density):
1773! 1) PIEVP, evaporation from melting ice
1774! 2) PIDEP<0, sublimation of ice
1775! 3) PIMLT, melting of ice because it is shed to rain
1776!
1777! So the final version is
1778! TOT_ICEnew=TOT_ICE+THICK*DELS,
1779! DELS=PSDEP+PIACWI+PIACR,
1780! PSDEP=MAX(0., PIDEP) only if PINIT<=EPSQ
1781!
1782!-----------------------------------------------------------------------
1783! Rime factor values associated with different processes.
1784!-----------------------------------------------------------------------
1785! 1) RF=1 for growth by vapor deposition.
1786! 2) RFmx, the RF associated with an ice density of 900 kg m^-3, the assumed
1787! properties for the freezing of rain to ice (PIACR).
1788! 3) RFrime, the RF associated with the density of rimed ice, the assumed
1789! properties for the riming of cloud water onto ice (PIACWI).
1790!
1791! The rime factor (RFnew) associated with the various microphysical processes is:
1792! RFnew=(1.*PSDEP+RFrime*PIACWI+RFmx*PIACR)/DELS,
1793!
1794! The updated rime factor (RimeF) becomes:
1795! RimeF*TOT_ICEnew=RimeF1*TOT_ICE+THICK*DELS*RFnew,
1796! RimeF=(RimeF1*TOT_ICE+THICK*DELS*RFnew)/TOT_ICEnew
1797!-----------------------------------------------------------------------
1798!
1799!-- Calculate density of rimed ice following Heymsfield and Pflaum (1985), where
1800! Rime density = 300.*(-Rcw*Vimpact/Tsfc)**0.44,
1801! in which Rcw is the mean diameter of the cloud droplets, Vimpact=VSNOW, and
1802! Tsfc (surface temperature of the particle, deg C) is approximated by TC.
1803!
1804! Here the calculations are extended to temperatures as warm as -0.5C, whereas
1805! the original study only looked at ice particles whose surface temperature were
1806! colder than -5C. Rime densities will vary from 170 to 900 kg m^-3
1807! (Straka, 2009 textbook; p. 308).
1808!
1809 dum=1.e-6*real(indexs) !- Mean diameter in m
1810 dum1=dum*dum*dum/massi(indexs) !- Used to estimate ice densities
1811 rfmx=rfmx1*dum1 !- Rime factor for density of pure ice (PIACR)
1812 IF (piacwi>0. .AND. rcw>0.) THEN
1813 dum=-rcw*vsnow/min(-0.5,tc)
1814 dum=min(12.14, max(0.275, dum) )
1815 dum=300.*dum**0.44 !- Initial rime density
1816 dum=min(900., max(170., dum) ) !- Final rime density, kg m^-3
1817 rfrime=pi*dum*dum1 !- Rime factor for the density of rimed ice
1818 ELSE
1819 rfrime=1. !- Homogeneous freezing of cloud water does not
1820 !- modify RF, contributes to cloud ice
1821 ENDIF
1822!
1823 rfnew=(psdep+rfrime*piacwi+rfmx*piacr)/dels
1824 dum2=tot_ice+thick*dels
1825 dum3=rimef1*tot_ice+rfnew*thick*dels
1826 IF (dum2<=climit) THEN
1827 rimef=rimef1
1828 ELSE
1829 rimef=min(rfmx, max(1., dum3/dum2) )
1830 ENDIF
1831!apr22 - end; see (40) above
1832 ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
1833!
1834 dum1=bldtrh*(flimass*vsnow+(1.-flimass)*vci)
1835 qinew=tot_icenew/(thick+dum1)
1836 IF (qinew .LE. epsq) qinew=0.
1837 IF (qi.GT.0. .AND. qinew.NE.0.) THEN
1838 dum=qinew/qi
1839 IF (dum .LT. toler) qinew=0.
1840 ENDIF
1841 asnownew=qinew*dum1
1842 IF (asnow.GT.0. .AND. asnownew.NE.0.) THEN
1843 dum=asnownew/asnow
1844 IF (dum .LT. toler) asnownew=0.
1845 ENDIF
1846 rqsnew=flimass*rho*qinew
1847 IF (rqsnew>0.) THEN !jul28 begin
1848!mar03 - Use mean diameter for snow/graupel (INDEXS) from above
1849 nsnow=rqsnew/(rimef*massi(indexs))
1850 IF (rqsnew>0.0025 .AND. rimef>5.) THEN
1851!mar12 - Blend to Nsnow=1.e3 for 2.5 g m^-3 < RQSnew <= 5 g m^-3 when RF>5
1852 IF (rqsnew<=0.005) THEN
1853 dum=min(1., 400.*rqsnew-1.) !- Blend at 2.5 < RQSnew (g m^-3) < 5
1854 nsnow=1.e3*dum+nsnow*(1.-dum) !- Final, blended Nsnow
1855 ELSE
1856!-- When RF>5, this branch will produce 56.1, 62.2, & 69.1 dBZ reflectivities
1857! when RQSnew reaches 5, 7.5, and 10 g m^-3 (respectively) due to
1858! Nsnow being reduced to 1, 0.55, and 0.2 L^-1 (respectively). A 60 dBZ is reached
1859! when RQSnew=6.6 g m^-3 & Nsnow=0.712 L^-1 (solving a quadratic eqn).
1860 dum=180.*(rqsnew-0.005) !- Steadily decrease Nsnow at > 5 g m^-3
1861 nsnow=1.e3*(1.-min(dum,0.8)) !- from 1 L^-1 down to 0.2 L^-1 at >=10 g m^-3
1862 ENDIF
1863 ENDIF
1864 zsnow=cdry*rqsnew*rqsnew/nsnow
1865 ENDIF !jul28 end
1866 ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
1867 ENDIF ! End IF (ICE_logical)
1868!
1869!--- Update rain mixing ratios
1870!
1871!---
1872! * TOT_RAINnew - total mass of rain after microphysics
1873! current layer and the input flux of ice from above
1874! * VRAIN2 - time-averaged fall speed of rain in grid and below
1875! * QRdum - first-guess estimate (dummy) rain mixing ratio in layer
1876! (uses rain fall speed from grid box above, VRAIN1)
1877! * QRnew - updated rain mixing ratio in layer
1878! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
1879! * ARAINnew - updated accumulation of rain at bottom of grid cell
1880!---
1881!
1882 timlt=timlt+pimlt*thick !may31
1883 delr=praut+pracw+piacwr-piacr+pimlt+prevp+picnd
1884 tot_rainnew=tot_rain+thick*delr
1885 IF (tot_rain.GT.0. .AND. tot_rainnew.NE.0.) THEN
1886 dum=tot_rainnew/tot_rain
1887 IF (dum .LT. toler) tot_rainnew=0.
1888 ENDIF
1889update_rn: IF (tot_rainnew .LE. climit) THEN
1890 tot_rainnew=0.
1891 vrain2=0.
1892 qrnew=0.
1893 arainnew=0.
1894 ELSE update_rn
1895!may11 - start
1896 qrnew=tot_rainnew/(thick+bldtrh*vrain2)
1897 rqrnew=rho*qrnew
1898!may20 - start - made slight changes to speed up algorithm, plus remove
1899! block that increases N0r for heavy rain.
1900 IF(strat .AND. rqrnew>1.e-3) strat=.false. !may31
1901 IF(drzl .AND. timlt>epsq) drzl=.false. !jun01
1902!
1903!-- Iterate for consistent QRnew, RQRnew, N0r, INDEXR2, VRAIN2
1904!
1905 indexr2=indexr1
1906rain_pass: DO ipass=1,3
1907dsd2: IF (rqrnew<=rqr_drmin) THEN
1908!-- Extremely light rain
1909 indexr=mdrmin
1910 n0r=rqrnew/massr(indexr)
1911 ELSE IF (drzl) THEN dsd2
1912!-- Drizzle - gradually increase N0r for rain contents <0.5 g m^-3 !may31
1913 dum=max(1.0, 0.5e-3/rqrnew)
1914 n0r=min(1.e9, n0r0*dum*sqrt(dum) ) !- 8.e6 <= N0r <= 1.e9
1915 dum1=rqrnew/(pi*rhol*n0r)
1916 dum=1.e6*sqrt(sqrt(dum1))
1917 indexr=max(xmrmin, min(dum, xmrmax) )
1918 n0r=rqrnew/massr(indexr)
1919 ELSE IF (rqrnew>=rqr_drmax) THEN
1920!-- Extremely heavy rain
1921 indexr=mdrmax
1922 n0r=rqrnew/massr(indexr)
1923 ELSE dsd2
1924!-- Light to heavy rain
1925 n0r=n0r0
1926 dum=cn0r0*sqrt(sqrt(rqrnew))
1927 indexr=max(xmrmin, min(dum, xmrmax) )
1928 IF (strat .AND. indexr<indexrmin) THEN
1929!-- Stratiform rain
1930 indexr=indexrmin
1931 n0r=rqrnew/massr(indexr)
1932 ENDIF
1933 ENDIF dsd2
1934 vrain2=gammar*vrain(indexr)
1935 qrnew=tot_rainnew/(thick+bldtrh*vrain2)
1936 rqrnew=rho*qrnew
1937 idr=abs(indexr-indexr2)
1938 indexr2=indexr !jun13
1939 IF(idr<=5) EXIT rain_pass !- Agreement within 5 microns
1940!may20 - end
1941 ENDDO rain_pass
1942!may11 - end
1943 IF (qrnew .LE. epsq) qrnew=0.
1944 IF (qr.GT.0. .AND. qrnew.NE.0.) THEN
1945 dum=qrnew/qr
1946 IF (dum .LT. toler) qrnew=0.
1947 ENDIF
1948 arainnew=bldtrh*vrain2*qrnew
1949 IF (arain.GT.0. .AND. arainnew.NE.0.) THEN
1950 dum=arainnew/arain
1951 IF (dum .LT. toler) arainnew=0.
1952 ENDIF
1953 rqrnew=rho*qrnew
1954!may11 - start
1955 IF (rqrnew>epsq) THEN
1956!may20 - remove may11 change in N0r calculation
1957 nrain=n0r*1.e-6*real(indexr2)
1958 zrain=crain*rqrnew*rqrnew/nrain !-- Rain reflectivity
1959 ENDIF
1960!may11 - end
1961 ENDIF update_rn
1962!
1963 wcnew=qwnew+qrnew+qinew
1964!
1965!----------------------------------------------------------------------
1966!-------------- Begin debugging & verification ------------------------
1967!----------------------------------------------------------------------
1968!
1969!--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
1970!
1971 qt=thick*(wv+wc)+arain+asnow
1972 qtnew=thick*(wvnew+wcnew)+arainnew+asnownew
1973 budget=qt-qtnew
1974!
1975!--- Additional check on budget preservation, accounting for truncation effects
1976!
1977 dbg_logical=.false.
1978 dum=abs(budget)
1979 IF (dum .GT. toler) THEN
1980 dum=dum/min(qt, qtnew)
1981 IF (dum .GT. toler) dbg_logical=.true.
1982 ENDIF
1983!
1984 IF (print_diag) THEN
1985 esw=min(1000.*fpvs0(tnew),0.99*pp)
1986 qswnew=(rhgrd+0.001)*eps*esw/(pp-esw)
1987 IF( (qwnew>epsq .OR. qrnew>epsq .OR. wvnew>qswnew) &
1988 & .AND. tc<t_ice) dbg_logical=.true.
1989 ENDIF
1990!
1991 IF ((wvnew.LT.epsq .OR. dbg_logical) .AND. print_diag) THEN
1992 !
1993 WRITE(0,"(/2(a,i4),2(a,i2))") '{} i=',i_index,' j=',j_index,&
1994 & ' L=',l,' LSFC=',lsfc
1995 !
1996 esw=min(1000.*fpvs0(tnew),0.99*pp)
1997 qswnew=eps*esw/(pp-esw)
1998 IF (tc.LT.0. .OR. tnew .LT. 0.) THEN
1999 esi=min(1000.*fpvs(tnew),0.99*pp)
2000 qsinew=eps*esi/(pp-esi)
2001 ELSE
2002 qsi=qsw
2003 qsinew=qswnew
2004 ENDIF
2005 wsnew=qsinew
2006 WRITE(0,"(4(a12,g11.4,1x))") &
2007 & '{} TCold=',tc,'TCnew=',tnew-t0c,'P=',.01*pp,'RHO=',rho, &
2008 & '{} THICK=',thick,'RHold=',wv/ws,'RHnew=',wvnew/wsnew, &
2009 & 'RHgrd=',rhgrd, &
2010 & '{} RHWold=',wv/qsw,'RHWnew=',wvnew/qswnew,'RHIold=',wv/qsi, &
2011 & 'RHInew=',wvnew/qsinew, &
2012 & '{} QSWold=',qsw,'QSWnew=',qswnew,'QSIold=',qsi,'QSInew=',qsinew, &
2013 & '{} WSold=',ws,'WSnew=',wsnew,'WVold=',wv,'WVnew=',wvnew, &
2014 & '{} WCold=',wc,'WCnew=',wcnew,'QWold=',qw,'QWnew=',qwnew, &
2015 & '{} QIold=',qi,'QInew=',qinew,'QRold=',qr,'QRnew=',qrnew, &
2016 & '{} ARAINold=',arain,'ARAINnew=',arainnew,'ASNOWold=',asnow, &
2017 & 'ASNOWnew=',asnownew, &
2018 & '{} TOT_RAIN=',tot_rain,'TOT_RAINnew=',tot_rainnew, &
2019 & 'TOT_ICE=',tot_ice,'TOT_ICEnew=',tot_icenew, &
2020 & '{} BUDGET=',budget,'QTold=',qt,'QTnew=',qtnew
2021 !
2022 WRITE(0,"(4(a12,g11.4,1x))") &
2023 & '{} DELT=',delt,'DELV=',delv,'DELW=',delw,'DELI=',deli, &
2024 & '{} DELR=',delr,'PCOND=',pcond,'PIDEP=',pidep,'PIEVP=',pievp, &
2025 & '{} PICND=',picnd,'PREVP=',prevp,'PRAUT=',praut,'PRACW=',pracw, &
2026 & '{} PIACW=',piacw,'PIACWI=',piacwi,'PIACWR=',piacwr,'PIMLT=', &
2027 & pimlt, &
2028 & '{} PIACR=',piacr
2029 !
2030 WRITE(0,"(4(a15,L2))") &
2031 & '{} ICE_logical=',ice_logical,'RAIN_logical=',rain_logical, &
2032 & 'STRAT=',strat,'DRZL=',drzl
2033 !
2034 IF (ice_logical) WRITE(0,"(4(a12,g11.4,1x))") &
2035 & '{} RimeF1=',rimef1,'GAMMAS=',gammas,'VrimeF=',vrimef, &
2036 & 'VSNOW=',vsnow, &
2037 & '{} QSmICE=',qsmice,'XLIMASS=',xlimass,'RQLICE=',rqlice, &
2038 & 'QTICE=',qtice, &
2039 & '{} NLICE=',nlice,'NSmICE=',nsmice,'PILOSS=',piloss, &
2040 & 'EMAIRI=',emairi, &
2041 & '{} RimeF=',rimef,'VCI=',vci,'INDEXS=',float(indexs), &
2042 & 'FLIMASS=',flimass, &
2043 & '{} Nsnow=',nsnow,'Zsnow=',zsnow,'TIMLT=',timlt,'RHOX0C=',rhox0c, &
2044 & '{} INDEXS0C=',float(indexs0c)
2045 !
2046 IF (rain_logical) WRITE(0,"(4(a12,g11.4,1x))") &
2047 & '{} INDEXR1=',float(indexr1),'INDEXR=',float(indexr), &
2048 & 'GAMMAR=',gammar,'N0r=',n0r, &
2049 & '{} VRAIN1=',vrain1,'VRAIN2=',vrain2,'QTRAIN=',qtrain,'RQR=',rqr, &
2050 & '{} PRLOSS=',prloss,'VOLR1=',thick+bldtrh*vrain1, &
2051 & 'VOLR2=',thick+bldtrh*vrain2,'INDEXR2=',indexr2, &
2052 & '{} Nrain=',nrain,'Zrain=',zrain
2053 !
2054 IF (pracw .GT. 0.) WRITE(0,"(a12,g11.4,1x)") '{} FWR=',fwr
2055 !
2056 IF (piacr .GT. 0.) WRITE(0,"(a12,g11.4,1x)") '{} FIR=',fir
2057 !
2058 dum=pimlt+picnd-prevp-pievp
2059 IF (dum.GT.0. .or. dwvi.NE.0.) &
2060 & WRITE(0,"(4(a12,g11.4,1x))") &
2061 & '{} TFACTOR=',tfactor,'DYNVIS=',dynvis, &
2062 & 'THERM_CON=',therm_cond,'DIFFUS=',diffus
2063 !
2064 IF (prevp .LT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2065 & '{} RFACTOR=',rfactor,'ABW=',abw,'VENTR=',ventr,'CREVP=',crevp, &
2066 & '{} DWVr=',dwvr,'DENOMW=',denomw
2067 !
2068 IF (pidep.NE.0. .AND. dwvi.NE.0.) &
2069 & WRITE(0,"(4(a12,g11.4,1x))") &
2070 & '{} DWVi=',dwvi,'DENOMI=',denomi,'PIDEP_max=',pidep_max, &
2071 & 'SFACTOR=',sfactor, &
2072 & '{} ABI=',abi,'VENTIL=',ventil,'VENTIL1=',venti1(indexs), &
2073 & 'VENTIL2=',sfactor*venti2(indexs), &
2074 & '{} VENTIS=',ventis,'DIDEP=',didep
2075 !
2076 IF (pidep.GT.0. .AND. pcond.NE.0.) &
2077 & WRITE(0,"(4(a12,g11.4,1x))") &
2078 & '{} DENOMW=',denomw,'DENOMWI=',denomwi,'DENOMF=',denomf, &
2079 & 'DUM2=',pcond-piacw
2080 !
2081 IF (fws .GT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2082 & '{} FWS=',fws
2083 !
2084 dum=pimlt+picnd-pievp
2085 IF (dum.GT. 0.) WRITE(0,"(4(a12,g11.4,1x))") &
2086 & '{} SFACTOR=',sfactor,'VENTIL=',ventil,'VENTIL1=',venti1(indexs), &
2087 & 'VENTIL2=',sfactor*venti2(indexs), &
2088 & '{} AIEVP=',aievp,'DIEVP=',dievp,'QSW0=',qsw0,'DWV0=',dwv0
2089 !
2090 ENDIF
2091!
2092!----------------------------------------------------------------------
2093!------------------------- Update arrays ------------------------------
2094!----------------------------------------------------------------------
2095!
2096 t_col(l)=tnew ! Update temperature
2097 q_col(l)=max(epsq, wvnew/(1.+wvnew)) ! Update specific humidity
2098 wc_col(l)=max(epsq, wcnew) ! Update total condensate mixing ratio
2099 qi_col(l)=max(epsq, qinew) ! Update ice mixing ratio
2100 qr_col(l)=max(epsq, qrnew) ! Update rain mixing ratio
2101 qw_col(l)=max(epsq, qwnew) ! Update cloud water mixing ratio
2102 rimef_col(l)=rimef ! Update rime factor
2103 nr_col(l)=nrain ! Update rain number concentration !jul28 begin
2104 ns_col(l)=nsnow ! Update snow number concentration
2105!aligo IF (RQRnew>RQRmix .AND. RQSnew>RQSmix .and. RimeF>10.) THEN
2106! Zsnow=Zsnow*Cwet !apr27
2107!aligo ENDIF
2108 ztot=max(ztot, zrain+zsnow)
2109 ztot=zrain+zsnow
2110 IF (ztot>zmin) dbz_col(l)=10.*alog10(ztot) ! Update radar reflectivity
2111 asnow=asnownew ! Update accumulated snow
2112 vsnow1=vsnow ! Update total ice fall speed
2113 arain=arainnew ! Update accumulated rain
2114 vrain1=vrain2 ! Update rain fall speed
2115!
2116!-- Assign microphysical processes and fall speeds to 1D array
2117!
2118 if(pcond.gt.0)then
2119 pcond1d(l)=pcond
2120 elseif(pcond.lt.0)then
2121 pevap1d(l)=pcond
2122 endif
2123 if(pidep.gt.0)then
2124 pidep1d(l)=pidep
2125 elseif(pidep.lt.0)then
2126 pisub1d(l)=pidep
2127 endif
2128 piacw1d(l)=piacw
2129 piacwi1d(l)=piacwi
2130 piacwr1d(l)=piacwr
2131 piacr1d(l)=piacr
2132 picnd1d(l)=picnd
2133 pievp1d(l)=pievp
2134 pimlt1d(l)=pimlt
2135 praut1d(l)=praut
2136 pracw1d(l)=pracw
2137 prevp1d(l)=prevp
2138 if (qinew>epsq) then
2139 vsnow1d(l)=vsnow
2140!sep22a - Start changes
2141 if (flimass<1.) then
2142 vci1d(l)=vci
2143 nsmice1d(l)=nsmice
2144 endif
2145!sep22a - End changes
2146 indexs1d(l)=float(indexs)
2147 endif
2148 if (qrnew>epsq) then
2149 vrain11d(l)=vrain1
2150 vrain21d(l)=vrain2
2151 indexr1d(l)=float(indexr2)
2152!jun01 - start
2153 idr=0
2154 IF(strat) idr=1
2155 IF(drzl) idr=idr+2
2156 rflag1d(l)=idr
2157!jun01 - end
2158 endif
2159!
2160!#######################################################################
2161!
2162 ENDDO big_loop
2163!
2164!#######################################################################
2165!
2166!-----------------------------------------------------------------------
2167!--------------------------- Return to GSMDRIVE -----------------------
2168!-----------------------------------------------------------------------
2169!
2170 CONTAINS
2171!#######################################################################
2172!--------- Produces accurate calculation of cloud condensation ---------
2173!#######################################################################
2174!
2176 REAL function condense (pp,qw,tk,wv,rhgrd,i,j,l)
2177!
2178!---------------------------------------------------------------------------------
2179!------ The Asai (1965) algorithm takes into consideration the release of ------
2180!------ latent heat in increasing the temperature & in increasing the ------
2181!------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
2182!---------------------------------------------------------------------------------
2183!
2184 IMPLICIT NONE
2185!
2186 REAL (kind=kind_phys), PARAMETER :: &
2187 & rhlimit=.001, rhlimit1=-rhlimit
2188 REAL (kind=kind_phys) :: cond, ssat, wcdum
2189!
2190 REAL,INTENT(IN) :: qw,pp,wv,tk,rhgrd
2191 REAL wvdum,tdum,dwv,ws,esw
2192
2193integer,INTENT(IN) :: i,j,l !-- Debug 20120111
2194integer :: niter
2195real :: dwvinp,ssatinp
2196
2197!
2198!-----------------------------------------------------------------------
2199!
2200 tdum=tk
2201 wvdum=wv
2202 wcdum=qw
2203 esw=min(1000.*fpvs0(tdum),0.99*pp) ! Saturation vapor press w/r/t water
2204 ws=rhgrd*eps*esw/(pp-esw) ! Saturation mixing ratio w/r/t water
2205 dwv=wvdum-ws ! Deficit grid-scale water vapor mixing ratio
2206 ssat=dwv/ws ! Supersaturation ratio
2207 condense=0.
2208
2209dwvinp=dwv !-- Debug 20120111
2210ssatinp=ssat
2211
2212 DO niter=1,max_iterations
2213 cond=dwv/(1.+xlv3*ws/(tdum*tdum)) ! Asai (1965, J. Japan)
2214 cond=max(cond, -wcdum) ! Limit cloud water evaporation
2215 tdum=tdum+xlv1*cond ! Updated temperature
2216 wvdum=wvdum-cond ! Updated water vapor mixing ratio
2217 wcdum=wcdum+cond ! Updated cloud water mixing ratio
2218 condense=condense+cond ! Total cloud water condensation
2219 esw=min(1000.*fpvs0(tdum),0.99*pp) ! Updated saturation vapor press w/r/t water
2220 ws=rhgrd*eps*esw/(pp-esw) ! Updated saturation mixing ratio w/r/t water
2221 dwv=wvdum-ws ! Deficit grid-scale water vapor mixing ratio
2222 ssat=dwv/ws ! Grid-scale supersaturation ratio
2223 IF (ssat>=rhlimit1 .AND. ssat<=rhlimit) EXIT !-- Exit if SSAT is near 0
2224 IF (ssat<rhlimit1 .AND. wcdum<=epsq) EXIT !-- Exit if SSAT<0 & no cloud water
2225 ENDDO
2226
2227!-- Debug 20120111
2228IF (niter>max_iterations) THEN
2229!-- Too many iterations - indicates possible numerical instability
2230 IF (warn3) THEN
2231 write(0,*) 'WARN3: Too many iterations in function CONDENSE; ', &
2232 ' I,J,L,TC,SSAT,QW,DWV=',i,j,l,tk-t0c,ssatinp,1000.*qw,dwvinp
2233 warn3=.false.
2234 ENDIF
2235 ssat=0.
2236 condense=dwvinp
2237ENDIF
2238
2239!
2240 END FUNCTION condense
2241!
2242!#######################################################################
2243!---------------- Calculate ice deposition at T<T_ICE ------------------
2244!#######################################################################
2245!
2247 REAL function deposit (pp,tdum,wvdum,rhgrd,i,j,l) !-- Debug 20120111
2248!
2249!--- Also uses the Asai (1965) algorithm, but uses a different target
2250! vapor pressure for the adjustment
2251!
2252 use machine, only: high_pres => kind_dbl_prec
2253 IMPLICIT NONE
2254!
2255 !INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
2256 REAL (kind=high_pres), PARAMETER :: rhlimit=.001, &
2257 & rhlimit1=-rhlimit
2258 REAL (kind=high_pres) :: dep, ssat
2259!
2260 real,INTENT(IN) :: pp,rhgrd
2261 real,INTENT(INOUT) :: wvdum,tdum
2262 real esi,ws,dwv
2263
2264integer,INTENT(IN) :: i,j,l !-- Debug 20120111
2265integer :: niter
2266real :: tinp,dwvinp,ssatinp
2267
2268!
2269!-----------------------------------------------------------------------
2270!
2271 esi=min(1000.*fpvs(tdum),0.99*pp) ! Saturation vapor press w/r/t ice
2272 ws=rhgrd*eps*esi/(pp-esi) ! Saturation mixing ratio
2273 dwv=wvdum-ws ! Deficit grid-scale water vapor mixing ratio
2274 ssat=dwv/ws ! Supersaturation ratio
2275 deposit=0.
2276
2277tinp=tdum !-- Debug 20120111
2278dwvinp=dwv
2279ssatinp=ssat
2280
2281 DO niter=1,max_iterations
2282 dep=dwv/(1.+xls3*ws/(tdum*tdum)) ! Asai (1965, J. Japan)
2283 tdum=tdum+xls1*dep ! Updated temperature
2284 wvdum=wvdum-dep ! Updated ice mixing ratio
2285 deposit=deposit+dep ! Total ice deposition
2286 esi=min(1000.*fpvs(tdum),0.99*pp) ! Updated saturation vapor press w/r/t ice
2287 ws=rhgrd*eps*esi/(pp-esi) ! Updated saturation mixing ratio w/r/t ice
2288 dwv=wvdum-ws ! Deficit grid-scale water vapor mixing ratio
2289 ssat=dwv/ws ! Grid-scale supersaturation ratio
2290 IF (ssat>=rhlimit1 .AND. ssat<=rhlimit) EXIT !-- Exit if SSAT is near 0
2291 ENDDO
2292
2293!-- Debug 20120111
2294IF (niter>max_iterations) THEN
2295!-- Too many iterations - indicates possible numerical instability
2296 IF (warn2) THEN
2297 write(0,*) 'WARN2: Too many iterations in function DEPOSIT; ', &
2298 ' I,J,L,TC,SSAT,DWV=',i,j,l,tinp-t0c,ssatinp,dwvinp
2299 warn2=.false.
2300 ENDIF
2301 ssat=0.
2302 deposit=dwvinp
2303ENDIF
2304
2305!
2306 END FUNCTION deposit
2307!
2308!#######################################################################
2309!--- Used to calculate the mean size of rain drops (INDEXR) in microns
2310!#######################################################################
2311!
2312!-- NOTE: this function is not used in this version.
2313!
2314 INTEGER FUNCTION get_indexr(RR)
2315 IMPLICIT NONE
2316 REAL, INTENT(IN) :: rr
2317 IF (rr .LE. rr_drmin) THEN
2318!
2319!--- Assume fixed mean diameter of rain (0.05 mm) for low rain rates,
2320! instead vary N0r with rain rate
2321!
2322 get_indexr=mdrmin
2323 ELSE IF (rr .LE. rr_dr1) THEN
2324!
2325!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2326! for mean diameters (Dr) between 0.05 and 0.10 mm:
2327! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
2328! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
2329! RR in kg/(m**2*s)
2330! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
2331!
2332 get_indexr=int( 1.123e3*rr**.1947 + .5 )
2333 get_indexr=max( mdrmin, min(get_indexr, mdr1) )
2334 ELSE IF (rr .LE. rr_dr2) THEN
2335!
2336!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2337! for mean diameters (Dr) between 0.10 and 0.20 mm:
2338! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
2339! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
2340! RR in kg/(m**2*s)
2341! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
2342!
2343 get_indexr=int( 1.225e3*rr**.2017 + .5 )
2344 get_indexr=max( mdr1, min(get_indexr, mdr2) )
2345 ELSE IF (rr .LE. rr_dr3) THEN
2346!
2347!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2348! for mean diameters (Dr) between 0.20 and 0.32 mm:
2349! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
2350! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
2351! RR in kg/(m**2*s)
2352! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
2353!
2354 get_indexr=int( 1.3006e3*rr**.2083 + .5 )
2355 get_indexr=max( mdr2, min(get_indexr, mdr3) )
2356 ELSE IF (rr .LE. rr_dr4) THEN
2357!
2358!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2359! for mean diameters (Dr) between 0.32 and 0.45 mm:
2360! V(Dr)=963.0*Dr**.666, V in m/s and Dr in m
2361! RR = PI*1000.*N0r0*963.0*Dr**(4+.666) = 2.4205e13*Dr**4.666,
2362! RR in kg/(m**2*s)
2363! Dr (m) = 1.354e-3*RR**.2143 -> Dr (microns) = 1.354e3*RR**.2143
2364!
2365 get_indexr=int( 1.354e3*rr**.2143 + .5 )
2366 get_indexr=max( mdr3, min(get_indexr, mdr4) )
2367 ELSE IF (rr .LE. rr_dr5) THEN
2368!
2369!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2370! for mean diameters (Dr) between 0.45 and 0.675 mm:
2371! V(Dr)=309.0*Dr**.5185, V in m/s and Dr in m
2372! RR = PI*1000.*N0r0*309.0*Dr**(4+.5185) = 7.766e12*Dr**4.5185,
2373! RR in kg/(m**2*s)
2374! Dr (m) = 1.404e-3*RR**.2213 -> Dr (microns) = 1.404e3*RR**.2213
2375!
2376 get_indexr=int( 1.404e3*rr**.2213 + .5 )
2377 get_indexr=max( mdr4, min(get_indexr, mdr5) )
2378 ELSE IF (rr .LE. rr_drmax) THEN
2379!
2380!--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
2381! for mean diameters (Dr) between 0.675 and 1.0 mm:
2382! V(Dr)=85.81Dr**.343, V in m/s and Dr in m
2383! RR = PI*1000.*N0r0*85.81*Dr**(4+.343) = 2.1566e12*Dr**4.343,
2384! RR in kg/(m**2*s)
2385! Dr (m) = 1.4457e-3*RR**.2303 -> Dr (microns) = 1.4457e3*RR**.2303
2386!
2387 get_indexr=int( 1.4457e3*rr**.2303 + .5 )
2388 get_indexr=max( mdr5, min(get_indexr, mdrmax) )
2389 ELSE
2390!
2391!--- Assume fixed mean diameter of rain (1.0 mm) for high rain rates,
2392! instead vary N0r with rain rate
2393!
2394 get_indexr=mdrmax
2395 ENDIF ! End IF (RR .LE. RR_DRmin) etc.
2396!
2397 END FUNCTION get_indexr
2398!
2399 END SUBROUTINE egcp01column_hr
2400!#######################################################################
2401!------- Initialize constants & lookup tables for microphysics ---------
2402!#######################################################################
2403!
2404
2405! SH 0211/2002
2406
2407!-----------------------------------------------------------------------
2409 SUBROUTINE ferrier_init_hr (GSMDT,MPI_COMM_COMP,MPIRANK,MPIROOT,THREADS, &
2410 errmsg,errflg)
2411!-----------------------------------------------------------------------
2412!-------------------------------------------------------------------------------
2413!--- SUBPROGRAM DOCUMENTATION BLOCK
2414! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
2415! Jin ORG: W/NP22 DATE: January 2002
2416! (Modification for WRF structure)
2417!
2418!-------------------------------------------------------------------------------
2419! ABSTRACT:
2420! * Reads various microphysical lookup tables used in COLUMN_MICRO
2421! * Lookup tables were created "offline" and are read in during execution
2422! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
2423!-------------------------------------------------------------------------------
2424!
2425! USAGE: CALL FERRIER_INIT_hr FROM SUBROUTINE PHYSICS_INITIALIZE
2426!
2427! INPUT ARGUMENT LIST:
2428! DTPH - physics time step (s)
2429!
2430! OUTPUT ARGUMENT LIST:
2431! NONE
2432!
2433! OUTPUT FILES:
2434! NONE
2435!
2436! SUBROUTINES:
2437! MY_GROWTH_RATES_NMM_hr - lookup table for growth of nucleated ice
2438! GPVS_hr - lookup tables for saturation vapor pressure (water, ice)
2439!
2440! UNIQUE: NONE
2441!
2442! LIBRARY: NONE
2443!
2444! ATTRIBUTES:
2445! LANGUAGE: FORTRAN 90
2446! MACHINE : IBM SP
2447!
2448!-----------------------------------------------------------------------
2449!
2450#ifdef MPI
2451 use mpi
2452#endif
2453 IMPLICIT NONE
2454!
2455!-------------------------------------------------------------------------
2456!-------------- Parameters & arrays for lookup tables --------------------
2457!-------------------------------------------------------------------------
2458!
2459!-----------------------------------------------------------------------
2460!--- Parameters & data statement for local calculations
2461!-----------------------------------------------------------------------
2462!
2463 INTEGER, PARAMETER :: mdr1=xmr1, mdr2=xmr2, mdr3=xmr3
2464!
2465! VARIABLES PASSED IN
2466 REAL, INTENT(IN) :: gsmdt
2467 INTEGER, INTENT(IN) :: mpirank
2468 INTEGER, INTENT(IN) :: mpiroot
2469 INTEGER, INTENT(IN) :: mpi_comm_comp
2470 INTEGER, INTENT(IN) :: threads
2471 CHARACTER(LEN=*), INTENT(OUT) :: errmsg
2472 INTEGER, INTENT(OUT) :: errflg
2473!
2474!-----------------------------------------------------------------------
2475! LOCAL VARIABLES
2476!-----------------------------------------------------------------------
2477 REAL :: bbfr,dtph,thour_print,rdis,beta6
2478 INTEGER :: i,j,l,k
2479 INTEGER :: etampnew_unit1
2480 LOGICAL :: opened
2481 INTEGER :: irtn,rc
2482 CHARACTER*80 errmess
2483 INTEGER :: ierr, good
2484 LOGICAL :: lexist,lopen, force_read_ferhires
2485!
2486!-----------------------------------------------------------------------
2487!
2488 dtph=gsmdt !-- Time step in s
2489!
2490!--- Create lookup tables for saturation vapor pressure w/r/t water & ice
2491
2492
2493!
2494 CALL gpvs_hr
2495!
2496!zhang:
2497 if (.NOT. ALLOCATED(ventr1)) ALLOCATE(ventr1(mdrmin:mdrmax))
2498 if (.NOT. ALLOCATED(ventr2)) ALLOCATE(ventr2(mdrmin:mdrmax))
2499 if (.NOT. ALLOCATED(accrr)) ALLOCATE(accrr(mdrmin:mdrmax))
2500 if (.NOT. ALLOCATED(massr)) ALLOCATE(massr(mdrmin:mdrmax))
2501 if (.NOT. ALLOCATED(vrain)) ALLOCATE(vrain(mdrmin:mdrmax))
2502 if (.NOT. ALLOCATED(rrate)) ALLOCATE(rrate(mdrmin:mdrmax))
2503 if (.NOT. ALLOCATED(venti1)) ALLOCATE(venti1(mdimin:mdimax))
2504 if (.NOT. ALLOCATED(venti2)) ALLOCATE(venti2(mdimin:mdimax))
2505 if (.NOT. ALLOCATED(accri)) ALLOCATE(accri(mdimin:mdimax))
2506 if (.NOT. ALLOCATED(massi)) ALLOCATE(massi(mdimin:mdimax))
2507 if (.NOT. ALLOCATED(vsnowi)) ALLOCATE(vsnowi(mdimin:mdimax))
2508 if (.NOT. ALLOCATED(vel_rf)) ALLOCATE(vel_rf(2:9,0:nrime))
2509
2510#ifdef MPI
2511 call mpi_barrier(mpi_comm_comp,ierr)
2512#endif
2513
2514 only_root_reads: if (mpirank==mpiroot) then
2515 force_read_ferhires = .true.
2516 good = 0
2517 INQUIRE(file="DETAMPNEW_DATA.expanded_rain_LE",exist=lexist)
2518
2519 IF (lexist) THEN
2520 OPEN(63,file="DETAMPNEW_DATA.expanded_rain_LE", &
2521 & form="UNFORMATTED",status="OLD",err=1234)
2522!
2523!sms$serial begin
2524 READ(63, err=1234) ventr1
2525 READ(63, err=1234) ventr2
2526 READ(63, err=1234) accrr
2527 READ(63, err=1234) massr
2528 READ(63, err=1234) vrain
2529 READ(63, err=1234) rrate
2530 READ(63, err=1234) venti1
2531 READ(63, err=1234) venti2
2532 READ(63, err=1234) accri
2533 READ(63, err=1234) massi
2534 READ(63, err=1234) vsnowi
2535 READ(63, err=1234) vel_rf
2536!sms$serial end
2537 good = 1
25381234 CONTINUE
2539 IF ( good .NE. 1 ) THEN
2540 INQUIRE(63,opened=lopen)
2541 IF (lopen) THEN
2542 IF( force_read_ferhires ) THEN
2543 errmsg = "Error reading DETAMPNEW_DATA.expanded_rain_LE. Aborting because force_read_ferhires is .true."
2544 errflg = 1
2545 return
2546 ENDIF
2547 CLOSE(63)
2548 ELSE
2549 IF( force_read_ferhires ) THEN
2550 errmsg = "Error opening DETAMPNEW_DATA.expanded_rain_LE. Aborting because force_read_ferhires is .true."
2551 errflg = 1
2552 return
2553 ENDIF
2554 ENDIF
2555 ELSE
2556 INQUIRE(63,opened=lopen)
2557 IF (lopen) THEN
2558 CLOSE(63)
2559 ENDIF
2560 ENDIF
2561 ELSE
2562 IF( force_read_ferhires ) THEN
2563 errmsg = "Non-existent DETAMPNEW_DATA.expanded_rain_LE. Aborting because force_read_ferhires is .true."
2564 errflg = 1
2565 return
2566 ENDIF
2567 ENDIF
2568 endif only_root_reads
2569!
2570#ifdef MPI
2571 CALL mpi_bcast(ventr1,SIZE(ventr1),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2572 CALL mpi_bcast(ventr2,SIZE(ventr2),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2573 CALL mpi_bcast(accrr, SIZE(accrr), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2574 CALL mpi_bcast(massr, SIZE(massr), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2575 CALL mpi_bcast(vrain, SIZE(vrain), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2576 CALL mpi_bcast(rrate, SIZE(rrate), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2577 CALL mpi_bcast(venti1,SIZE(venti1),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2578 CALL mpi_bcast(venti2,SIZE(venti2),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2579 CALL mpi_bcast(accri, SIZE(accri), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2580 CALL mpi_bcast(massi, SIZE(massi), mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2581 CALL mpi_bcast(vsnowi,SIZE(vsnowi),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2582 CALL mpi_bcast(vel_rf,SIZE(vel_rf),mpi_double_precision,mpiroot,mpi_comm_comp,irtn)
2583#endif
2584
2585!
2586!--- Calculates coefficients for growth rates of ice nucleated in water
2587! saturated conditions, scaled by physics time step (lookup table)
2588!
2589
2590 CALL my_growth_rates_nmm_hr (dtph)
2591!
2592!
2593!--- Constants associated with Biggs (1953) freezing of rain, as parameterized
2594! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
2595!
2596 abfr=-0.66
2597 bbfr=100.
2598 cbfr=20.*pi*pi*bbfr*rhol*1.e-21*dtph !mar03 - Bug fix in (33); include DTPH
2599!
2600!--- CIACW is used in calculating riming rates
2601! The assumed effective collection efficiency of cloud water rimed onto
2602! ice is =0.5 below:
2603!
2604 ciacw=0.5*dtph*0.25*pi !jul28 (16)
2605!
2606!--- CIACR is used in calculating freezing of rain colliding with large ice
2607! The assumed collection efficiency is 1.0
2608!
2609 ciacr=pi*dtph
2610!
2611!--- Based on rain lookup tables for mean diameters from 0.05 to 1.0 mm
2612! * Four different functional relationships of mean drop diameter as
2613! a function of rain rate (RR), derived based on simple fits to
2614! mass-weighted fall speeds of rain as functions of mean diameter
2615! from the lookup tables.
2616!
2617 rr_drmin=n0r0*rrate(mdrmin) ! RR for mean drop diameter of .05 mm
2618 rr_dr1=n0r0*rrate(mdr1) ! RR for mean drop diameter of .10 mm
2619 rr_dr2=n0r0*rrate(mdr2) ! RR for mean drop diameter of .20 mm
2620 rr_dr3=n0r0*rrate(mdr3) ! RR for mean drop diameter of .32 mm
2621 rr_dr4=n0r0*rrate(mdr4) ! RR for mean drop diameter of .45 mm
2622 rr_dr5=n0r0*rrate(mdr5) ! RR for mean drop diameter of .675 mm
2623 rr_drmax=n0r0*rrate(mdrmax) ! RR for mean drop diameter of 1.0 mm
2624!
2625 rqr_drmin=n0r0*massr(mdrmin) ! Rain content for mean drop diameter of .05 mm
2626 rqr_drmax=n0r0*massr(mdrmax) ! Rain content for mean drop diameter of 1.0 mm
2627 c_nr=1./(pi*rhol) !jul28
2628 crain=720.e18*c_nr*c_nr !jul28
2629 c_n0r0=pi*rhol*n0r0
2630 cn0r0=1.e6/sqrt(sqrt(c_n0r0))
2631 cn0r_dmrmin=1./(pi*rhol*dmrmin*dmrmin*dmrmin*dmrmin)
2632 cn0r_dmrmax=1./(pi*rhol*dmrmax*dmrmax*dmrmax*dmrmax)
2633!
2634!--- CRACW is used in calculating collection of cloud water by rain (an
2635! assumed collection efficiency of 1.0)
2636!
2637 cracw=dtph*0.25*pi*1.0
2638!
2639 esw0=1000.*fpvs0(t0c) ! Saturation vapor pressure at 0C
2640 rfmax=1.1**nrime ! Maximum rime factor allowed
2641 rfmx1=pi*900. ! Density near that of pure ice, 900 kg m^-3
2642!
2643!------------------------------------------------------------------------
2644!--------------- Constants passed through argument list -----------------
2645!------------------------------------------------------------------------
2646!
2647!--- Important parameters for self collection (autoconversion) of
2648! cloud water to rain.
2649!
2650!-- Relative dispersion == standard deviation of droplet spectrum / mean radius
2651! (see pp 1542-1543, Liu & Daum, JAS, 2004)
2652 rdis=0.5 !-- relative dispersion of droplet spectrum
2653 beta6=( (1.+3.*rdis*rdis)*(1.+4.*rdis*rdis)*(1.+5.*rdis*rdis)/ &
2654 & ((1.+rdis*rdis)*(1.+2.*rdis*rdis) ) )
2655!-- Kappa=1.1e10 g^-2 cm^3 s^-1 after eq. (8b) on p.1105 of Liu et al. (JAS, 2006)
2656! => More extensive units conversion than can be described here to go from
2657! eq. (13) in Liu et al. (JAS, 2006) to what's programmed below. Note that
2658! the units used throughout the paper are in cgs units!
2659!
2660 araut=1.03e19/(ncw*sqrt(ncw))
2661 braut=dtph*1.1e10*beta6/ncw
2662!
2663!--- QAUT0 is the *OLD* threshold cloud content for autoconversion to rain
2664! needed for droplets to reach a diameter of 20 microns (following
2665! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM). It's no longer
2666! used in this version, but the value is passed into radiation in case
2667! a ball park estimate is needed.
2668!--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
2669! of 300, 200, and 100 cm**-3, respectively
2670!
2671 qaut0=pi*rhol*ncw*(20.e-6)**3/6. !-- legacy
2672!
2673!--- For calculating cloud droplet radius in microns, Rcw
2674!
2675 arcw=1.e6*(0.75/(pi*ncw*rhol))**c1
2676!
2677!may20 - start
2678!
2679!- An explanation for the following settings assumed for "hail" or frozen drops (RF>10):
2680! RH_NgC=PI*500.*5.E3
2681! RHOg=500 kg m^-3, Ng=5.e3 m^-3 => "hail" content >7.85 g m^-3 for INDEXS=MDImax
2682!- or -
2683! RH_NgC=PI*500.*1.E3
2684! RHOg=500 kg m^-3, Ng=1.e3 m^-3 => "hail" content >1.57 g m^-3 for INDEXS=MDImax
2685!
2686 rh_ngc=pi*500.*1.e3 !- RHOg=500 kg m^-3, Ng=1.e3 m^-3
2687 rqhail=rh_ngc*(1.e-3)**3 !- Threshold "hail" content ! becomes 1.57 g m^-3
2688 bvhail=0.82*c1 !- Exponent for Thompson graupel
2689 avhail=1353.*(1./rh_ngc)**bvhail !- 1353 ~ constant for Thompson graupel
2690 rh_ngc=1.e6*(1./rh_ngc)**c1 !mar08 (convection, convert to microns)
2691!
2692!- An explanation for the following settings assumed for graupel (RF>5):
2693! RH_NgT=PI*300.*10.E3
2694! RHOg=300 kg m^-3, Ng=10.e3 m^-3 => "graupel" content must exceed 9.43 g m^-3 for INDEXS=MDImax
2695!- or -
2696! RH_NgT=PI*300.*5.E3
2697! RHOg=300 kg m^-3, Ng=5.e3 m^-3 => "graupel" content must exceed 4.71 g m^-3 for INDEXS=MDImax
2698!
2699 rh_ngt=pi*300.*5.e3 !- RHOg=300 kg m^-3, Ng=5.e3 m^-3
2700 rh_ngt=1.e6*(1./rh_ngt)**c1 !mar08 (transition, convert to microns)
2701!may20 - end
2702!
2703!--- For calculating snow optical depths by considering bulk density of
2704! snow based on emails from Q. Fu (6/27-28/01), where optical
2705! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
2706! is effective radius, and DENS is the bulk density of snow.
2707!
2708! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
2709! T = 1.5*1.E3*SWPrad/(Reff*DENS)
2710!
2711! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
2712!
2713! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
2714!
2715 DO i=mdimin,mdimax
2716 sdens(i)=pi*1.5e-15*float(i*i*i)/massi(i)
2717 ENDDO
2718!
2719 thour_print=-dtph/3600.
2720!
2721
2722 RETURN
2723!
2724!-----------------------------------------------------------------------
2725 END SUBROUTINE ferrier_init_hr
2726!
2728 SUBROUTINE my_growth_rates_nmm_hr (DTPH)
2729!
2730!--- Below are tabulated values for the predicted mass of ice crystals
2731! after 600 s of growth in water saturated conditions, based on
2732! calculations from Miller and Young (JAS, 1979). These values are
2733! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
2734! Young (1993). Values at temperatures colder than -27C were
2735! assumed to be invariant with temperature.
2736!
2737!--- Used to normalize Miller & Young (1979) calculations of ice growth
2738! over large time steps using their tabulated values at 600 s.
2739! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
2740!
2741 IMPLICIT NONE
2742!
2743 REAL,INTENT(IN) :: DTPH
2744!
2745 REAL DT_ICE
2746 REAL,DIMENSION(35) :: MY_600
2747!WRF
2748!
2749!-----------------------------------------------------------------------
2750!-- 20090714: These values are in g and need to be converted to kg below
2751 DATA my_600 / &
2752 & 5.5e-8, 1.4e-7, 2.8e-7, 6.e-7, 3.3e-6, &
2753 & 2.e-6, 9.e-7, 8.8e-7, 8.2e-7, 9.4e-7, &
2754 & 1.2e-6, 1.85e-6, 5.5e-6, 1.5e-5, 1.7e-5, &
2755 & 1.5e-5, 1.e-5, 3.4e-6, 1.85e-6, 1.35e-6, &
2756 & 1.05e-6, 1.e-6, 9.5e-7, 9.0e-7, 9.5e-7, &
2757 & 9.5e-7, 9.e-7, 9.e-7, 9.e-7, 9.e-7, &
2758 & 9.e-7, 9.e-7, 9.e-7, 9.e-7, 9.e-7 / ! -31 to -35 deg C
2759!
2760!-----------------------------------------------------------------------
2761!
2762 dt_ice=(dtph/600.)**1.5
2763 my_growth_nmm=dt_ice*my_600*1.e-3 !-- 20090714: Convert from g to kg
2764!
2765!-----------------------------------------------------------------------
2766!
2767 END SUBROUTINE my_growth_rates_nmm_hr
2768!
2769!-----------------------------------------------------------------------
2770!--------- Old GFS saturation vapor pressure lookup tables -----------
2771!-----------------------------------------------------------------------
2772!
2774 SUBROUTINE gpvs_hr
2775! ******************************************************************
2776!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2777! . . .
2778! SUBPROGRAM: GPVS_hr COMPUTE SATURATION VAPOR PRESSURE TABLE
2779! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2780!
2781! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
2782! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
2783! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
2784! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
2785! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
2786!
2787! PROGRAM HISTORY LOG:
2788! 91-05-07 IREDELL
2789! 94-12-30 IREDELL EXPAND TABLE
2790! 96-02-19 HONG ICE EFFECT
2791! 01-11-29 JIN MODIFIED FOR WRF
2792!
2793! USAGE: CALL GPVS_hr
2794!
2795! SUBPROGRAMS CALLED:
2796! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
2797!
2798! ATTRIBUTES:
2799! LANGUAGE: FORTRAN 90
2800!
2801!$$$
2802 IMPLICIT NONE
2803 real :: x,xinc,t
2804 integer :: jx
2805!----------------------------------------------------------------------
2806 xinc=(xmax-xmin)/(nx-1)
2807 c1xpvs=1.-xmin/xinc
2808 c2xpvs=1./xinc
2809 c1xpvs0=1.-xmin/xinc
2810 c2xpvs0=1./xinc
2811!
2812 DO jx=1,nx
2813 x=xmin+(jx-1)*xinc
2814 t=x
2815 tbpvs(jx)=fpvsx(t)
2816 tbpvs0(jx)=fpvsx0(t)
2817 ENDDO
2818!
2819 END SUBROUTINE gpvs_hr
2820!-----------------------------------------------------------------------
2821!***********************************************************************
2822!-----------------------------------------------------------------------
2823 REAL function fpvs(t)
2824!-----------------------------------------------------------------------
2825!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2826! . . .
2827! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
2828! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2829!
2830! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
2831! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
2832! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
2833! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
2834! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
2835! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
2836! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2837!
2838! PROGRAM HISTORY LOG:
2839! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2840! 94-12-30 IREDELL EXPAND TABLE
2841! 96-02-19 HONG ICE EFFECT
2842! 01-11-29 JIN MODIFIED FOR WRF
2843!
2844! USAGE: PVS=FPVS(T)
2845!
2846! INPUT ARGUMENT LIST:
2847! T - REAL TEMPERATURE IN KELVIN
2848!
2849! OUTPUT ARGUMENT LIST:
2850! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2851!
2852! ATTRIBUTES:
2853! LANGUAGE: FORTRAN 90
2854!$$$
2855 IMPLICIT NONE
2856 real,INTENT(IN) :: t
2857 real xj
2858 integer :: jx
2859!-----------------------------------------------------------------------
2860 IF (t>=xmin .AND. t<=xmax) THEN
2861 xj=min(max(c1xpvs+c2xpvs*t,1.),float(nx))
2862 jx=min(xj,nx-1.)
2863 fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
2864 ELSE IF (t>xmax) THEN
2865!-- Magnus Tetens formula for water saturation (Murray, 1967)
2866! (saturation vapor pressure in kPa)
2867 fpvs=0.61078*exp(17.2694*(t-273.16)/(t-35.86))
2868 ELSE
2869!-- Magnus Tetens formula for ice saturation(Murray, 1967)
2870! (saturation vapor pressure in kPa)
2871 fpvs=0.61078*exp(21.8746*(t-273.16)/(t-7.66))
2872 ENDIF
2873!
2874 END FUNCTION fpvs
2875!-----------------------------------------------------------------------
2876!-----------------------------------------------------------------------
2877 REAL function fpvs0(t)
2878!-----------------------------------------------------------------------
2879 IMPLICIT NONE
2880 real,INTENT(IN) :: t
2881 real :: xj1
2882 integer :: jx1
2883!-----------------------------------------------------------------------
2884 IF (t>=xmin .AND. t<=xmax) THEN
2885 xj1=min(max(c1xpvs0+c2xpvs0*t,1.),float(nx))
2886 jx1=min(xj1,nx-1.)
2887 fpvs0=tbpvs0(jx1)+(xj1-jx1)*(tbpvs0(jx1+1)-tbpvs0(jx1))
2888 ELSE
2889!-- Magnus Tetens formula for water saturation (Murray, 1967)
2890! (saturation vapor pressure in kPa)
2891 fpvs0=0.61078*exp(17.2694*(t-273.16)/(t-35.86))
2892 ENDIF
2893!
2894 END FUNCTION fpvs0
2895!-----------------------------------------------------------------------
2896!***********************************************************************
2897!-----------------------------------------------------------------------
2898 REAL function fpvsx(t)
2899!-----------------------------------------------------------------------
2900!$$$ SUBPROGRAM DOCUMENTATION BLOCK
2901! . . .
2902! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
2903! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2904!
2905! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
2906! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
2907! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
2908! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
2909! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
2910! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
2911! TO GET THE FORMULA
2912! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2913! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
2914! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2915!
2916! PROGRAM HISTORY LOG:
2917! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2918! 94-12-30 IREDELL EXACT COMPUTATION
2919! 96-02-19 HONG ICE EFFECT
2920! 01-11-29 JIN MODIFIED FOR WRF
2921!
2922! USAGE: PVS=FPVSX(T)
2923! REFERENCE: EMANUEL(1994),116-117
2924!
2925! INPUT ARGUMENT LIST:
2926! T - REAL TEMPERATURE IN KELVIN
2927!
2928! OUTPUT ARGUMENT LIST:
2929! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2930!
2931! ATTRIBUTES:
2932! LANGUAGE: FORTRAN 90
2933!$$$
2934 IMPLICIT NONE
2935!-----------------------------------------------------------------------
2936 real, parameter :: ttp=2.7316e+2,hvap=2.5000e+6,psat=6.1078e+2 &
2937 , cliq=4.1855e+3,cvap= 1.8460e+3 &
2938 , cice=2.1060e+3,hsub=2.8340e+6
2939!
2940 real, parameter :: psatk=psat*1.e-3
2941 real, parameter :: dldt=cvap-cliq,xa=-dldt/rv,xb=xa+hvap/(rv*ttp)
2942 real, parameter :: dldti=cvap-cice &
2943 , xai=-dldti/rv,xbi=xai+hsub/(rv*ttp)
2944 real t,tr
2945!-----------------------------------------------------------------------
2946 tr=ttp/t
2947!
2948 IF(t.GE.ttp)THEN
2949 fpvsx=psatk*(tr**xa)*exp(xb*(1.-tr))
2950 ELSE
2951 fpvsx=psatk*(tr**xai)*exp(xbi*(1.-tr))
2952 ENDIF
2953!
2954 END FUNCTION fpvsx
2955!-----------------------------------------------------------------------
2956!-----------------------------------------------------------------------
2957 REAL function fpvsx0(t)
2958!-----------------------------------------------------------------------
2959 IMPLICIT NONE
2960 real,parameter :: ttp=2.7316e+2,hvap=2.5000e+6,psat=6.1078e+2 &
2961 , cliq=4.1855e+3,cvap=1.8460e+3 &
2962 , cice=2.1060e+3,hsub=2.8340e+6
2963 real,PARAMETER :: psatk=psat*1.e-3
2964 real,PARAMETER :: dldt=cvap-cliq,xa=-dldt/rv,xb=xa+hvap/(rv*ttp)
2965 real,PARAMETER :: dldti=cvap-cice &
2966 , xai=-dldt/rv,xbi=xa+hsub/(rv*ttp)
2967 real :: t,tr
2968!-----------------------------------------------------------------------
2969 tr=ttp/t
2970 fpvsx0=psatk*(tr**xa)*exp(xb*(1.-tr))
2971!
2972 END FUNCTION fpvsx0
2973
2974 SUBROUTINE ferhires_finalize()
2975
2976 IMPLICIT NONE
2977
2978 if (ALLOCATED(ventr1)) DEALLOCATE(ventr1)
2979 if (ALLOCATED(ventr2)) DEALLOCATE(ventr2)
2980 if (ALLOCATED(accrr)) DEALLOCATE(accrr)
2981 if (ALLOCATED(massr)) DEALLOCATE(massr)
2982 if (ALLOCATED(vrain)) DEALLOCATE(vrain)
2983 if (ALLOCATED(rrate)) DEALLOCATE(rrate)
2984 if (ALLOCATED(venti1)) DEALLOCATE(venti1)
2985 if (ALLOCATED(venti2)) DEALLOCATE(venti2)
2986 if (ALLOCATED(accri)) DEALLOCATE(accri)
2987 if (ALLOCATED(massi)) DEALLOCATE(massi)
2988 if (ALLOCATED(vsnowi)) DEALLOCATE(vsnowi)
2989 if (ALLOCATED(vel_rf)) DEALLOCATE(vel_rf)
2990
2991 END SUBROUTINE ferhires_finalize
2992
2993!
2994 END MODULE module_mp_fer_hires
subroutine, public ferrier_init_hr(gsmdt, mpi_comm_comp, mpirank, mpiroot, threads, errmsg, errflg)
subroutine my_growth_rates_nmm_hr(dtph)
subroutine egcp01column_hr(arain, asnow, dtph, rhc_col, i_index, j_index, lsfc, p_col, qi_col, qr_col, q_col, qw_col, rimef_col, t_col, thick_col, wc_col, lm, pcond1d, pidep1d, piacw1d, piacwi1d, piacwr1d, piacr1d, picnd1d, pievp1d, pimlt1d, praut1d, pracw1d, prevp1d, pisub1d, pevap1d, dbz_col, nr_col, ns_col, vsnow1d, vrain11d, vrain21d, vci1d, nsmice1d, indexs1d, indexr1d, rflag1d, dx1)
This is the grid-scale microphysical processes of Ferrier-Aligo microphysics scheme (i....
subroutine fer_hires(dt, rhgrd, prsi, p_phy, t_phy, q, qt, lowlyr, sr, train_phy, f_ice_phy, f_rain_phy, f_rimef_phy, qc, qr, qs, rainnc, rainncv, threads, ims, ime, lm, d_ss, refl_10cm, dx1)
This is the driver scheme of Ferrier-Aligo microphysics scheme. NOTE: The only differences between FE...
real function condense(pp, qw, tk, wv, rhgrd, i, j, l)
real function deposit(pp, tdum, wvdum, rhgrd, i, j, l)