module_mp_etanew.F
References to this file elsewhere.
1 !WRF:MODEL_MP:PHYSICS
2 !
3 MODULE module_mp_etanew
4 !
5 !-----------------------------------------------------------------------
6 REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, &
7 & CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
8 & RFmax, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
9 & RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax
10 !
11 INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35
12 REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH
13 !
14 REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, &
15 & DelDMI=1.e-6,XMImin=1.e6*DMImin
16 INTEGER, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536, &
17 & MDImin=XMImin, MDImax=XMImax
18 REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: &
19 & ACCRI,SDENS,VSNOWI,VENTI1,VENTI2
20 !
21 REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=.45e-3, &
22 & DelDMR=1.e-6,XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax
23 INTEGER, PRIVATE,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax
24 REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: &
25 & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
26 !
27 INTEGER, PRIVATE,PARAMETER :: Nrime=40
28 REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF
29 !
30 INTEGER,PARAMETER :: NX=7501
31 REAL, PARAMETER :: XMIN=180.0,XMAX=330.0
32 REAL, DIMENSION(NX),PRIVATE,SAVE :: TBPVS,TBPVS0
33 REAL, PRIVATE,SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS
34 !
35 REAL, PRIVATE,PARAMETER :: &
36 !--- Physical constants follow:
37 & CP=1004.6, EPSQ=1.E-12, GRAV=9.806, RHOL=1000., RD=287.04 &
38 & ,RV=461.5, T0C=273.15, XLS=2.834E6 &
39 !--- Derived physical constants follow:
40 & ,EPS=RD/RV, EPS1=RV/RD-1., EPSQ1=1.001*EPSQ &
41 & ,RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV, RRHOL=1./RHOL &
42 & ,XLS1=XLS*RCP, XLS2=XLS*XLS*RCPRV, XLS3=XLS*XLS/RV &
43 !--- Constants specific to the parameterization follow:
44 !--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
45 & ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT &
46 & ,C1=1./3. &
47 & ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3 &
48 & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3
49 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
50 !
51 ! ======================================================================
52 !--- Important tunable parameters that are exported to other modules
53 ! * RHgrd - threshold relative humidity for onset of condensation
54 ! * T_ICE - temperature (C) threshold at which all remaining liquid water
55 ! is glaciated to ice
56 ! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
57 ! * NLImax - maximum number concentrations (m**-3) of large ice (snow/graupel/sleet)
58 ! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
59 ! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 0.45 mm
60 ! * N0rmin - minimum intercept (m**-4) for rain drops
61 ! * NCW - number concentrations of cloud droplets (m**-3)
62 ! * FLARGE1, FLARGE2 - number fraction of large ice to total (large+snow) ice
63 ! at T>0C and in presence of sublimation (FLARGE1), otherwise in
64 ! presence of ice saturated/supersaturated conditions
65 ! ======================================================================
66 REAL, PUBLIC,PARAMETER :: &
67 & RHgrd=1. &
68 & ,T_ICE=-30. &
69 & ,T_ICEK=T0C+T_ICE &
70 & ,T_ICE_init=-5. &
71 & ,NLImax=5.E3 &
72 & ,NLImin=1.E3 &
73 & ,N0r0=8.E6 &
74 & ,N0rmin=1.E4 &
75 & ,NCW=100.E6 &
76 & ,FLARGE1=1. &
77 & ,FLARGE2=.2
78 !--- Other public variables passed to other routines:
79 REAL,PUBLIC,SAVE :: QAUT0
80 REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
81 !
82 !
83 CONTAINS
84
85 !-----------------------------------------------------------------------
86 !-----------------------------------------------------------------------
87 SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
88 & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, &
89 & LOWLYR,SR, &
90 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
91 & QC,QR,QS, &
92 & mp_restart_state,tbpvs_state,tbpvs0_state, &
93 & RAINNC,RAINNCV, &
94 & ids,ide, jds,jde, kds,kde, &
95 & ims,ime, jms,jme, kms,kme, &
96 & its,ite, jts,jte, kts,kte )
97 !-----------------------------------------------------------------------
98 IMPLICIT NONE
99 !-----------------------------------------------------------------------
100 INTEGER, PARAMETER :: ITLO=-60, ITHI=40
101
102 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
103 & ,IMS,IME,JMS,JME,KMS,KME &
104 & ,ITS,ITE,JTS,JTE,KTS,KTE &
105 & ,ITIMESTEP
106
107 REAL, INTENT(IN) :: DT,DX,DY
108 REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
109 & dz8w,p_phy,pi_phy,rho_phy
110 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
111 & th_phy,qv,qt
112 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
113 & qc,qr,qs
114 REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
115 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
116 REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
117 & RAINNC,RAINNCV
118 REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
119 !
120 REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
121 !
122 REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
123 !
124 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
125
126 !-----------------------------------------------------------------------
127 ! LOCAL VARS
128 !-----------------------------------------------------------------------
129
130 ! NSTATS,QMAX,QTOT are diagnostic vars
131
132 INTEGER,DIMENSION(ITLO:ITHI,4) :: NSTATS
133 REAL, DIMENSION(ITLO:ITHI,5) :: QMAX
134 REAL, DIMENSION(ITLO:ITHI,22):: QTOT
135
136 ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
137 ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
138 ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
139
140 ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
141 ! the microphysics scheme. Instead, they will be used by Eta precip
142 ! assimilation.
143
144 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
145 & TLATGS_PHY,TRAIN_PHY
146 REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
147 REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
148
149 INTEGER :: I,J,K,KFLIP
150 REAL :: WC
151 !
152 !-----------------------------------------------------------------------
153 !**********************************************************************
154 !-----------------------------------------------------------------------
155 !
156 MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
157 !
158 C1XPVS0=MP_RESTART_STATE(MY_T2+1)
159 C2XPVS0=MP_RESTART_STATE(MY_T2+2)
160 C1XPVS =MP_RESTART_STATE(MY_T2+3)
161 C2XPVS =MP_RESTART_STATE(MY_T2+4)
162 CIACW =MP_RESTART_STATE(MY_T2+5)
163 CIACR =MP_RESTART_STATE(MY_T2+6)
164 CRACW =MP_RESTART_STATE(MY_T2+7)
165 CRAUT =MP_RESTART_STATE(MY_T2+8)
166 !
167 TBPVS(1:NX) =TBPVS_STATE(1:NX)
168 TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
169 !
170 DO j = jts,jte
171 DO k = kts,kte
172 DO i = its,ite
173 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
174 qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
175 ENDDO
176 ENDDO
177 ENDDO
178
179 ! initial diagnostic variables and data assimilation vars
180 ! (will need to delete this part in the future)
181
182 DO k = 1,4
183 DO i = ITLO,ITHI
184 NSTATS(i,k)=0.
185 ENDDO
186 ENDDO
187
188 DO k = 1,5
189 DO i = ITLO,ITHI
190 QMAX(i,k)=0.
191 ENDDO
192 ENDDO
193
194 DO k = 1,22
195 DO i = ITLO,ITHI
196 QTOT(i,k)=0.
197 ENDDO
198 ENDDO
199
200 ! initial data assimilation vars (will need to delete this part in the future)
201
202 DO j = jts,jte
203 DO k = kts,kte
204 DO i = its,ite
205 TLATGS_PHY (i,k,j)=0.
206 TRAIN_PHY (i,k,j)=0.
207 ENDDO
208 ENDDO
209 ENDDO
210
211 DO j = jts,jte
212 DO i = its,ite
213 ACPREC(i,j)=0.
214 APREC (i,j)=0.
215 PREC (i,j)=0.
216 SR (i,j)=0.
217 ENDDO
218 ENDDO
219
220 !-----------------------------------------------------------------------
221
222 CALL EGCP01DRV(DT,LOWLYR, &
223 & APREC,PREC,ACPREC,SR,NSTATS,QMAX,QTOT, &
224 & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
225 & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
226 & ids,ide, jds,jde, kds,kde, &
227 & ims,ime, jms,jme, kms,kme, &
228 & its,ite, jts,jte, kts,kte )
229 !-----------------------------------------------------------------------
230
231 DO j = jts,jte
232 DO k = kts,kte
233 DO i = its,ite
234 th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
235 qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
236 WC=qt(I,K,J)
237 QS(I,K,J)=0.
238 QR(I,K,J)=0.
239 QC(I,K,J)=0.
240 IF(F_ICE_PHY(I,K,J)>=1.)THEN
241 QS(I,K,J)=WC
242 ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
243 QC(I,K,J)=WC
244 ELSE
245 QS(I,K,J)=F_ICE_PHY(I,K,J)*WC
246 QC(I,K,J)=WC-QS(I,K,J)
247 ENDIF
248 !
249 IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
250 IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
251 QR(I,K,J)=QC(I,K,J)
252 QC(I,K,J)=0.
253 ELSE
254 QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
255 QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
256 ENDIF
257 ENDIF
258 ENDDO
259 ENDDO
260 ENDDO
261 !
262 ! update rain (from m to mm)
263
264 DO j=jts,jte
265 DO i=its,ite
266 RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
267 RAINNCV(i,j)=APREC(i,j)*1000.
268 ENDDO
269 ENDDO
270 !
271 MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
272 MP_RESTART_STATE(MY_T2+1)=C1XPVS0
273 MP_RESTART_STATE(MY_T2+2)=C2XPVS0
274 MP_RESTART_STATE(MY_T2+3)=C1XPVS
275 MP_RESTART_STATE(MY_T2+4)=C2XPVS
276 MP_RESTART_STATE(MY_T2+5)=CIACW
277 MP_RESTART_STATE(MY_T2+6)=CIACR
278 MP_RESTART_STATE(MY_T2+7)=CRACW
279 MP_RESTART_STATE(MY_T2+8)=CRAUT
280 !
281 TBPVS_STATE(1:NX) =TBPVS(1:NX)
282 TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
283
284 !-----------------------------------------------------------------------
285
286 END SUBROUTINE ETAMP_NEW
287
288 !-----------------------------------------------------------------------
289
290 SUBROUTINE EGCP01DRV( &
291 & DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
292 & NSTATS,QMAX,QTOT, &
293 & dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
294 & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
295 & ids,ide, jds,jde, kds,kde, &
296 & ims,ime, jms,jme, kms,kme, &
297 & its,ite, jts,jte, kts,kte)
298 !-----------------------------------------------------------------------
299 ! DTPH Physics time step (s)
300 ! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
301 ! Q_PHY Mixing ratio of water vapor. kg/kg
302 ! F_RAIN_PHY Fraction of rain.
303 ! F_ICE_PHY Fraction of ice.
304 ! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
305 !
306 !TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
307 !micrphysics sechme. Instead, they will be used by Eta precip assimilation.
308 !
309 !NSTATS,QMAX,QTOT are used for diagnosis purposes.
310 !
311 !-----------------------------------------------------------------------
312 !--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
313 ! and/or ZHAO's scheme in Eta and are not required by this microphysics
314 ! scheme itself.
315 !--- NSTATS,QMAX,QTOT are used for diagnosis purposes only. They will be
316 ! printed out when PRINT_diag is true.
317 !
318 !-----------------------------------------------------------------------
319 IMPLICIT NONE
320 !-----------------------------------------------------------------------
321 !
322 INTEGER, PARAMETER :: ITLO=-60, ITHI=40
323 LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
324 ! VARIABLES PASSED IN/OUT
325 INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
326 & ,ims,ime, jms,jme, kms,kme &
327 & ,its,ite, jts,jte, kts,kte
328 REAL,INTENT(IN) :: DTPH
329 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
330 INTEGER,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
331 REAL,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
332 REAL,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
333 REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
334 & APREC,PREC,ACPREC,SR
335 REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
336 REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
337 & dz8w,P_PHY,RHO_PHY
338 REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
339 & CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
340 & ,Q_PHY,TRAIN_PHY
341 !
342 !-----------------------------------------------------------------------
343 !LOCAL VARIABLES
344 !-----------------------------------------------------------------------
345 !
346 #define CACHE_FRIENDLY_MP_ETANEW
347 #ifdef CACHE_FRIENDLY_MP_ETANEW
348 # define TEMP_DIMS kts:kte,its:ite,jts:jte
349 # define TEMP_DEX L,I,J
350 #else
351 # define TEMP_DIMS its:ite,jts:jte,kts:kte
352 # define TEMP_DEX I,J,L
353 #endif
354 !
355 INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
356 REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
357 REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
358 INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
359 REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
360 REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
361 RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL
362 REAL,DIMENSION(2) :: PRECtot,PRECmax
363 !-----------------------------------------------------------------------
364 !
365 DO J=JTS,JTE
366 DO I=ITS,ITE
367 LMH(I,J) = KTE-LOWLYR(I,J)+1
368 ENDDO
369 ENDDO
370
371 DO 98 J=JTS,JTE
372 DO 98 I=ITS,ITE
373 DO L=KTS,KTE
374 KFLIP=KTE+1-L
375 CWM(TEMP_DEX)=CWM_PHY(I,KFLIP,J)
376 T(TEMP_DEX)=T_PHY(I,KFLIP,J)
377 Q(TEMP_DEX)=Q_PHY(I,KFLIP,J)
378 P(TEMP_DEX)=P_PHY(I,KFLIP,J)
379 TLATGS(TEMP_DEX)=TLATGS_PHY(I,KFLIP,J)
380 TRAIN(TEMP_DEX)=TRAIN_PHY(I,KFLIP,J)
381 F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
382 F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
383 F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
384 ENDDO
385 98 CONTINUE
386
387 DO 100 J=JTS,JTE
388 DO 100 I=ITS,ITE
389 LSFC=LMH(I,J) ! "L" of surface
390 !
391 DO K=KTS,KTE
392 KFLIP=KTE+1-K
393 DPCOL(K)=RHO_PHY(I,KFLIP,J)*GRAV*dz8w(I,KFLIP,J)
394 ENDDO
395 !
396 !
397 !--- Initialize column data (1D arrays)
398 !
399 L=1
400 IF (CWM(TEMP_DEX) .LE. EPSQ) CWM(TEMP_DEX)=EPSQ
401 F_ice(1,I,J)=1.
402 F_rain(1,I,J)=0.
403 F_RimeF(1,I,J)=1.
404 DO L=1,LSFC
405 !
406 !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
407 !
408 P_col(L)=P(TEMP_DEX)
409 !
410 !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
411 !
412 THICK_col(L)=DPCOL(L)*RGRAV
413 T_col(L)=T(TEMP_DEX)
414 TC=T_col(L)-T0C
415 QV_col(L)=max(EPSQ, Q(TEMP_DEX))
416 IF (CWM(TEMP_DEX) .LE. EPSQ1) THEN
417 WC_col(L)=0.
418 IF (TC .LT. T_ICE) THEN
419 F_ice(L,I,J)=1.
420 ELSE
421 F_ice(L,I,J)=0.
422 ENDIF
423 F_rain(L,I,J)=0.
424 F_RimeF(L,I,J)=1.
425 ELSE
426 WC_col(L)=CWM(TEMP_DEX)
427 ENDIF
428 !
429 !--- Determine composition of condensate in terms of
430 ! cloud water, ice, & rain
431 !
432 WC=WC_col(L)
433 QI=0.
434 QR=0.
435 QW=0.
436 Fice=F_ice(L,I,J)
437 Frain=F_rain(L,I,J)
438 IF (Fice .GE. 1.) THEN
439 QI=WC
440 ELSE IF (Fice .LE. 0.) THEN
441 QW=WC
442 ELSE
443 QI=Fice*WC
444 QW=WC-QI
445 ENDIF
446 IF (QW.GT.0. .AND. Frain.GT.0.) THEN
447 IF (Frain .GE. 1.) THEN
448 QR=QW
449 QW=0.
450 ELSE
451 QR=Frain*QW
452 QW=QW-QR
453 ENDIF
454 ENDIF
455 RimeF_col(L)=F_RimeF(L,I,J) ! (real)
456 QI_col(L)=QI
457 QR_col(L)=QR
458 QW_col(L)=QW
459 ENDDO
460 !
461 !#######################################################################
462 !
463 !--- Perform the microphysical calculations in this column
464 !
465 I_index=I
466 J_index=J
467 CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
468 & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
469 & THICK_col, WC_col,KTS,KTE,NSTATS,QMAX,QTOT )
470
471
472 !
473 !#######################################################################
474 !
475 !
476 !--- Update storage arrays
477 !
478 DO L=1,LSFC
479 TRAIN(TEMP_DEX)=(T_col(L)-T(TEMP_DEX))/DTPH
480 TLATGS(TEMP_DEX)=T_col(L)-T(TEMP_DEX)
481 T(TEMP_DEX)=T_col(L)
482 Q(TEMP_DEX)=QV_col(L)
483 CWM(TEMP_DEX)=WC_col(L)
484 !
485 !--- REAL*4 array storage
486 !
487 F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
488 IF (QI_col(L) .LE. EPSQ) THEN
489 F_ice(L,I,J)=0.
490 IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
491 ELSE
492 F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
493 ENDIF
494 IF (QR_col(L) .LE. EPSQ) THEN
495 DUM=0
496 ELSE
497 DUM=QR_col(L)/(QR_col(L)+QW_col(L))
498 ENDIF
499 F_rain(L,I,J)=DUM
500 !
501 ENDDO
502 !
503 !--- Update accumulated precipitation statistics
504 !
505 !--- Surface precipitation statistics; SR is fraction of surface
506 ! precipitation (if >0) associated with snow
507 !
508 APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
509 PREC(I,J)=PREC(I,J)+APREC(I,J)
510 ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
511 IF(APREC(I,J) .LT. 1.E-8) THEN
512 SR(I,J)=0.
513 ELSE
514 SR(I,J)=RRHOL*ASNOW/APREC(I,J)
515 ENDIF
516 !
517 !--- Debug statistics
518 !
519 IF (PRINT_diag) THEN
520 PRECtot(1)=PRECtot(1)+ARAIN
521 PRECtot(2)=PRECtot(2)+ASNOW
522 PRECmax(1)=MAX(PRECmax(1), ARAIN)
523 PRECmax(2)=MAX(PRECmax(2), ASNOW)
524 ENDIF
525
526
527 !#######################################################################
528 !#######################################################################
529 !
530 100 CONTINUE ! End "I" & "J" loops
531 DO 101 J=JTS,JTE
532 DO 101 I=ITS,ITE
533 DO L=KTS,KTE
534 KFLIP=KTE+1-L
535 CWM_PHY(I,KFLIP,J)=CWM(TEMP_DEX)
536 T_PHY(I,KFLIP,J)=T(TEMP_DEX)
537 Q_PHY(I,KFLIP,J)=Q(TEMP_DEX)
538 TLATGS_PHY(I,KFLIP,J)=TLATGS(TEMP_DEX)
539 TRAIN_PHY(I,KFLIP,J)=TRAIN(TEMP_DEX)
540 F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
541 F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
542 F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
543 ENDDO
544 101 CONTINUE
545 END SUBROUTINE EGCP01DRV
546 !
547 !
548 !###############################################################################
549 ! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
550 ! (1) Represents sedimentation by preserving a portion of the precipitation
551 ! through top-down integration from cloud-top. Modified procedure to
552 ! Zhao and Carr (1997).
553 ! (2) Microphysical equations are modified to be less sensitive to time
554 ! steps by use of Clausius-Clapeyron equation to account for changes in
555 ! saturation mixing ratios in response to latent heating/cooling.
556 ! (3) Prevent spurious temperature oscillations across 0C due to
557 ! microphysics.
558 ! (4) Uses lookup tables for: calculating two different ventilation
559 ! coefficients in condensation and deposition processes; accretion of
560 ! cloud water by precipitation; precipitation mass; precipitation rate
561 ! (and mass-weighted precipitation fall speeds).
562 ! (5) Assumes temperature-dependent variation in mean diameter of large ice
563 ! (Houze et al., 1979; Ryan et al., 1996).
564 ! -> 8/22/01: This relationship has been extended to colder temperatures
565 ! to parameterize smaller large-ice particles down to mean sizes of MDImin,
566 ! which is 50 microns reached at -55.9C.
567 ! (6) Attempts to differentiate growth of large and small ice, mainly for
568 ! improved transition from thin cirrus to thick, precipitating ice
569 ! anvils.
570 ! -> 8/22/01: This feature has been diminished by effectively adjusting to
571 ! ice saturation during depositional growth at temperatures colder than
572 ! -10C. Ice sublimation is calculated more explicitly. The logic is
573 ! that sources of are either poorly understood (e.g., nucleation for NWP)
574 ! or are not represented in the Eta model (e.g., detrainment of ice from
575 ! convection). Otherwise the model is too wet compared to the radiosonde
576 ! observations based on 1 Feb - 18 March 2001 retrospective runs.
577 ! (7) Top-down integration also attempts to treat mixed-phase processes,
578 ! allowing a mixture of ice and water. Based on numerous observational
579 ! studies, ice growth is based on nucleation at cloud top &
580 ! subsequent growth by vapor deposition and riming as the ice particles
581 ! fall through the cloud. Effective nucleation rates are a function
582 ! of ice supersaturation following Meyers et al. (JAM, 1992).
583 ! -> 8/22/01: The simulated relative humidities were far too moist compared
584 ! to the rawinsonde observations. This feature has been substantially
585 ! diminished, limited to a much narrower temperature range of 0 to -10C.
586 ! (8) Depositional growth of newly nucleated ice is calculated for large time
587 ! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
588 ! using their ice crystal masses calculated after 600 s of growth in water
589 ! saturated conditions. The growth rates are normalized by time step
590 ! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
591 ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
592 ! (9) Ice precipitation rates can increase due to increase in response to
593 ! cloud water riming due to (a) increased density & mass of the rimed
594 ! ice, and (b) increased fall speeds of rimed ice.
595 ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
596 !###############################################################################
597 !###############################################################################
598 !
599 SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, &
600 & LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
601 & THICK_col, WC_col ,KTS,KTE,NSTATS,QMAX,QTOT)
602 !
603 !###############################################################################
604 !###############################################################################
605 !
606 !-------------------------------------------------------------------------------
607 !----- NOTE: Code is currently set up w/o threading!
608 !-------------------------------------------------------------------------------
609 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
610 ! . . .
611 ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
612 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
613 ! PRGRMMR: Jin (Modification for WRF structure)
614 !-------------------------------------------------------------------------------
615 ! ABSTRACT:
616 ! * Merges original GSCOND & PRECPD subroutines.
617 ! * Code has been substantially streamlined and restructured.
618 ! * Exchange between water vapor & small cloud condensate is calculated using
619 ! the original Asai (1965, J. Japan) algorithm. See also references to
620 ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
621 ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
622 ! parameterization.
623 !-------------------------------------------------------------------------------
624 !
625 ! USAGE:
626 ! * CALL EGCP01COLUMN FROM SUBROUTINE EGCP01DRV
627 !
628 ! INPUT ARGUMENT LIST:
629 ! DTPH - physics time step (s)
630 ! I_index - I index
631 ! J_index - J index
632 ! LSFC - Eta level of level above surface, ground
633 ! P_col - vertical column of model pressure (Pa)
634 ! QI_col - vertical column of model ice mixing ratio (kg/kg)
635 ! QR_col - vertical column of model rain ratio (kg/kg)
636 ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
637 ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
638 ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
639 ! T_col - vertical column of model temperature (deg K)
640 ! THICK_col - vertical column of model mass thickness (density*height increment)
641 ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
642 !
643 !
644 ! OUTPUT ARGUMENT LIST:
645 ! ARAIN - accumulated rainfall at the surface (kg)
646 ! ASNOW - accumulated snowfall at the surface (kg)
647 ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
648 ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
649 ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
650 ! QI_col - vertical column of model ice mixing ratio (kg/kg)
651 ! QR_col - vertical column of model rain ratio (kg/kg)
652 ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
653 ! T_col - vertical column of model temperature (deg K)
654 !
655 ! OUTPUT FILES:
656 ! NONE
657 !
658 ! Subprograms & Functions called:
659 ! * Real Function CONDENSE - cloud water condensation
660 ! * Real Function DEPOSIT - ice deposition (not sublimation)
661 !
662 ! UNIQUE: NONE
663 !
664 ! LIBRARY: NONE
665 !
666 ! COMMON BLOCKS:
667 ! CMICRO_CONS - key constants initialized in GSMCONST
668 ! CMICRO_STATS - accumulated and maximum statistics
669 ! CMY_GROWTH - lookup table for growth of ice crystals in
670 ! water saturated conditions (Miller & Young, 1979)
671 ! IVENT_TABLES - lookup tables for ventilation effects of ice
672 ! IACCR_TABLES - lookup tables for accretion rates of ice
673 ! IMASS_TABLES - lookup tables for mass content of ice
674 ! IRATE_TABLES - lookup tables for precipitation rates of ice
675 ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
676 ! RVENT_TABLES - lookup tables for ventilation effects of rain
677 ! RACCR_TABLES - lookup tables for accretion rates of rain
678 ! RMASS_TABLES - lookup tables for mass content of rain
679 ! RVELR_TABLES - lookup tables for fall speeds of rain
680 ! RRATE_TABLES - lookup tables for precipitation rates of rain
681 !
682 ! ATTRIBUTES:
683 ! LANGUAGE: FORTRAN 90
684 ! MACHINE : IBM SP
685 !
686 !
687 !-------------------------------------------------------------------------
688 !--------------- Arrays & constants in argument list ---------------------
689 !-------------------------------------------------------------------------
690 !
691 IMPLICIT NONE
692 !
693 INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
694 REAL,INTENT(INOUT) :: ARAIN, ASNOW
695 REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col, QI_col,QR_col &
696 & ,QV_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col
697 !
698 !-------------------------------------------------------------------------
699 !-------------- Common blocks for microphysical statistics ---------------
700 !-------------------------------------------------------------------------
701 !
702 !-------------------------------------------------------------------------
703 !--------- Common blocks for constants initialized in GSMCONST ----------
704 !
705 INTEGER, PARAMETER :: ITLO=-60, ITHI=40
706 INTEGER,INTENT(INOUT) :: NSTATS(ITLO:ITHI,4)
707 REAL,INTENT(INOUT) :: QMAX(ITLO:ITHI,5),QTOT(ITLO:ITHI,22)
708 !
709 !-------------------------------------------------------------------------
710 !--------------- Common blocks for various lookup tables -----------------
711 !
712 !--- Discretized growth rates of small ice crystals after their nucleation
713 ! at 1 C intervals from -1 C to -35 C, based on calculations by Miller
714 ! and Young (1979, JAS) after 600 s of growth. Resultant growth rates
715 ! are multiplied by physics time step in GSMCONST.
716 !
717 !-------------------------------------------------------------------------
718 !
719 !--- Mean ice-particle diameters varying from 50 microns to 1000 microns
720 ! (1 mm), assuming an exponential size distribution.
721 !
722 !---- Meaning of the following arrays:
723 ! - mdiam - mean diameter (m)
724 ! - VENTI1 - integrated quantity associated w/ ventilation effects
725 ! (capacitance only) for calculating vapor deposition onto ice
726 ! - VENTI2 - integrated quantity associated w/ ventilation effects
727 ! (with fall speed) for calculating vapor deposition onto ice
728 ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
729 ! - MASSI - integrated quantity associated w/ ice mass
730 ! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
731 ! precipitation rates
732 !
733 !
734 !-------------------------------------------------------------------------
735 !
736 !--- VEL_RF - velocity increase of rimed particles as functions of crude
737 ! particle size categories (at 0.1 mm intervals of mean ice particle
738 ! sizes) and rime factor (different values of Rime Factor of 1.1**N,
739 ! where N=0 to Nrime).
740 !
741 !-------------------------------------------------------------------------
742 !
743 !--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
744 ! (0.45 mm), assuming an exponential size distribution.
745 !
746 !-------------------------------------------------------------------------
747 !------- Key parameters, local variables, & important comments ---------
748 !-----------------------------------------------------------------------
749 !
750 !--- TOLER => Tolerance or precision for accumulated precipitation
751 !
752 REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, Xratio=.025
753 !
754 !--- If BLEND=1:
755 ! precipitation (large) ice amounts are estimated at each level as a
756 ! blend of ice falling from the grid point above and the precip ice
757 ! present at the start of the time step (see TOT_ICE below).
758 !--- If BLEND=0:
759 ! precipitation (large) ice amounts are estimated to be the precip
760 ! ice present at the start of the time step.
761 !
762 !--- Extended to include sedimentation of rain on 2/5/01
763 !
764 REAL, PARAMETER :: BLEND=1.
765 !
766 !--- This variable is for debugging purposes (if .true.)
767 !
768 LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
769 !
770 !-----------------------------------------------------------------------
771 !--- Local variables
772 !-----------------------------------------------------------------------
773 !
774 REAL EMAIRI, N0r, NLICE, NSmICE
775 LOGICAL CLEAR, ICE_logical, DBG_logical, RAIN_logical
776 INTEGER :: IDR,INDEX_MY,INDEXR,INDEXR1,INDEXS,IPASS,ITDX,IXRF, &
777 & IXS,LBEF,L
778 !
779 REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
780 & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
781 & DENOMI,DENOMW,DENOMWI,DIDEP, &
782 & DIEVP,DIFFUS,DLI,DTPH,DTRHO,DUM,DUM1, &
783 & DUM2,DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLARGE,FLIMASS, &
784 & FSMALL,FWR,FWS,GAMMAR,GAMMAS, &
785 & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
786 & PIEVP,PILOSS,PIMLT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
787 & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
788 & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,QV,QW,QW0,QWnew, &
789 & RFACTOR,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
790 & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
791 & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
792 & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
793 & WC,WCnew,WSgrd,WS,WSnew,WV,WVnew,WVQW, &
794 & XLF,XLF1,XLI,XLV,XLV1,XLV2,XLIMASS,XRF,XSIMASS
795 !
796 !#######################################################################
797 !########################## Begin Execution ############################
798 !#######################################################################
799 !
800 !
801 ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
802 ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
803 !
804 !-----------------------------------------------------------------------
805 !------------ Loop from top (L=1) to surface (L=LSFC) ------------------
806 !-----------------------------------------------------------------------
807 !
808
809 DO 10 L=1,LSFC
810
811 !--- Skip this level and go to the next lower level if no condensate
812 ! and very low specific humidities
813 !
814 IF (QV_col(L).LE.EPSQ .AND. WC_col(L).LE.EPSQ) GO TO 10
815 !
816 !-----------------------------------------------------------------------
817 !------------ Proceed with cloud microphysics calculations -------------
818 !-----------------------------------------------------------------------
819 !
820 TK=T_col(L) ! Temperature (deg K)
821 TC=TK-T0C ! Temperature (deg C)
822 PP=P_col(L) ! Pressure (Pa)
823 QV=QV_col(L) ! Specific humidity of water vapor (kg/kg)
824 WV=QV/(1.-QV) ! Water vapor mixing ratio (kg/kg)
825 WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
826 !
827 !-----------------------------------------------------------------------
828 !--- Moisture variables below are mixing ratios & not specifc humidities
829 !-----------------------------------------------------------------------
830 !
831 CLEAR=.TRUE.
832 !
833 !--- This check is to determine grid-scale saturation when no condensate is present
834 !
835 ESW=1000.*FPVS0(TK) ! Saturation vapor pressure w/r/t water
836 QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
837 WS=QSW ! General saturation mixing ratio (water/ice)
838 IF (TC .LT. 0.) THEN
839 ESI=1000.*FPVS(TK) ! Saturation vapor pressure w/r/t ice
840 QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
841 WS=QSI ! General saturation mixing ratio (water/ice)
842 ENDIF
843 !
844 !--- Effective grid-scale Saturation mixing ratios
845 !
846 QSWgrd=RHgrd*QSW
847 QSIgrd=RHgrd*QSI
848 WSgrd=RHgrd*WS
849 !
850 !--- Check if air is subsaturated and w/o condensate
851 !
852 IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
853 !
854 !--- Check if any rain is falling into layer from above
855 !
856 IF (ARAIN .GT. CLIMIT) THEN
857 CLEAR=.FALSE.
858 ELSE
859 ARAIN=0.
860 ENDIF
861 !
862 !--- Check if any ice is falling into layer from above
863 !
864 !--- NOTE that "SNOW" in variable names is synonomous with
865 ! large, precipitation ice particles
866 !
867 IF (ASNOW .GT. CLIMIT) THEN
868 CLEAR=.FALSE.
869 ELSE
870 ASNOW=0.
871 ENDIF
872 !
873 !-----------------------------------------------------------------------
874 !-- Loop to the end if in clear, subsaturated air free of condensate ---
875 !-----------------------------------------------------------------------
876 !
877 IF (CLEAR) GO TO 10
878 !
879 !-----------------------------------------------------------------------
880 !--------- Initialize RHO, THICK & microphysical processes -------------
881 !-----------------------------------------------------------------------
882 !
883 !
884 !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
885 ! (see pp. 63-65 in Fleagle & Businger, 1963)
886 !
887 RHO=PP/(RD*TK*(1.+EPS1*QV)) ! Air density (kg/m**3)
888 RRHO=1./RHO ! Reciprocal of air density
889 DTRHO=DTPH*RHO ! Time step * air density
890 BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density
891 THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
892 !
893 ARAINnew=0. ! Updated accumulated rainfall
894 ASNOWnew=0. ! Updated accumulated snowfall
895 QI=QI_col(L) ! Ice mixing ratio
896 QInew=0. ! Updated ice mixing ratio
897 QR=QR_col(L) ! Rain mixing ratio
898 QRnew=0. ! Updated rain ratio
899 QW=QW_col(L) ! Cloud water mixing ratio
900 QWnew=0. ! Updated cloud water ratio
901 !
902 PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
903 PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
904 PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
905 PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
906 PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
907 PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
908 PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
909 PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
910 PIMLT=0. ! Melting ice (kg/kg; >0)
911 PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
912 PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
913 PREVP=0. ! Rain evaporation (kg/kg; <0)
914 !
915 !--- Double check input hydrometeor mixing ratios
916 !
917 ! DUM=WC-(QI+QW+QR)
918 ! DUM1=ABS(DUM)
919 ! DUM2=TOLER*MIN(WC, QI+QW+QR)
920 ! IF (DUM1 .GT. DUM2) THEN
921 ! WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
922 ! & ' L=',L
923 ! WRITE(6,"(4(a12,g11.4,1x))")
924 ! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
925 ! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
926 ! ENDIF
927 !
928 !***********************************************************************
929 !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
930 !***********************************************************************
931 !
932 !--- Calculate a few variables, which are used more than once below
933 !
934 !--- Latent heat of vaporization as a function of temperature from
935 ! Bolton (1980, JAS)
936 !
937 XLV=3.148E6-2370*TK ! Latent heat of vaporization (Lv)
938 XLF=XLS-XLV ! Latent heat of fusion (Lf)
939 XLV1=XLV*RCP ! Lv/Cp
940 XLF1=XLF*RCP ! Lf/Cp
941 TK2=1./(TK*TK) ! 1./TK**2
942 XLV2=XLV*XLV*QSW*TK2/RV ! Lv**2*Qsw/(Rv*TK**2)
943 DENOMW=1.+XLV2*RCP ! Denominator term, Clausius-Clapeyron correction
944 !
945 !--- Basic thermodynamic quantities
946 ! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
947 ! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
948 ! * DIFFUS - diffusivity of water vapor [ m**2/s ]
949 !
950 TFACTOR=TK**1.5/(TK+120.)
951 DYNVIS=1.496E-6*TFACTOR
952 THERM_COND=2.116E-3*TFACTOR
953 DIFFUS=8.794E-5*TK**1.81/PP
954 !
955 !--- Air resistance term for the fall speed of ice following the
956 ! basic research by Heymsfield, Kajikawa, others
957 !
958 GAMMAS=(1.E5/PP)**C1
959 !
960 !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
961 !
962 GAMMAR=(RHO0/RHO)**.4
963 !
964 !----------------------------------------------------------------------
965 !------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
966 !----------------------------------------------------------------------
967 !
968 !--- Determine if conditions supporting ice are present
969 !
970 IF (TC.LT.0. .OR. QI.GT. EPSQ .OR. ASNOW.GT.CLIMIT) THEN
971 ICE_logical=.TRUE.
972 ELSE
973 ICE_logical=.FALSE.
974 QLICE=0.
975 QTICE=0.
976 ENDIF
977 !
978 !--- Determine if rain is present
979 !
980 RAIN_logical=.FALSE.
981 IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
982 !
983 IF (ICE_logical) THEN
984 !
985 !--- IMPORTANT: Estimate time-averaged properties.
986 !
987 !---
988 ! * FLARGE - ratio of number of large ice to total (large & small) ice
989 ! * FSMALL - ratio of number of small ice crystals to large ice particles
990 ! -> Small ice particles are assumed to have a mean diameter of 50 microns.
991 ! * XSIMASS - used for calculating small ice mixing ratio
992 !---
993 ! * TOT_ICE - total mass (small & large) ice before microphysics,
994 ! which is the sum of the total mass of large ice in the
995 ! current layer and the input flux of ice from above
996 ! * PILOSS - greatest loss (<0) of total (small & large) ice by
997 ! sublimation, removing all of the ice falling from above
998 ! and the ice within the layer
999 ! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
1000 ! ice mass to the unrimed ice mass (>=1)
1001 ! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
1002 ! * VSNOW - Fall speed of rimed snow w/ air resistance correction
1003 ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
1004 ! * XLIMASS - used for calculating large ice mixing ratio
1005 ! * FLIMASS - mass fraction of large ice
1006 ! * QTICE - time-averaged mixing ratio of total ice
1007 ! * QLICE - time-averaged mixing ratio of large ice
1008 ! * NLICE - time-averaged number concentration of large ice
1009 ! * NSmICE - number concentration of small ice crystals at current level
1010 !---
1011 !--- Assumed number fraction of large ice particles to total (large & small)
1012 ! ice particles, which is based on a general impression of the literature.
1013 !
1014 WVQW=WV+QW ! Water vapor & cloud water
1015 !
1016
1017
1018 IF (TC.GE.0. .OR. WVQW.LT.QSIgrd) THEN
1019 !
1020 !--- Eliminate small ice particle contributions for melting & sublimation
1021 !
1022 FLARGE=FLARGE1
1023 ELSE
1024 !
1025 !--- Enhanced number of small ice particles during depositional growth
1026 ! (effective only when 0C > T >= T_ice [-10C] )
1027 !
1028 FLARGE=FLARGE2
1029 !
1030 !--- Larger number of small ice particles due to rime splintering
1031 !
1032 IF (TC.GE.-8. .AND. TC.LE.-3.) FLARGE=.5*FLARGE
1033 !
1034 ENDIF ! End IF (TC.GE.0. .OR. WVQW.LT.QSIgrd)
1035 FSMALL=(1.-FLARGE)/FLARGE
1036 XSIMASS=RRHO*MASSI(MDImin)*FSMALL
1037 IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT) THEN
1038 INDEXS=MDImin
1039 TOT_ICE=0.
1040 PILOSS=0.
1041 RimeF1=1.
1042 VrimeF=1.
1043 VEL_INC=GAMMAS
1044 VSNOW=0.
1045 EMAIRI=THICK
1046 XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
1047 FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
1048 QLICE=0.
1049 QTICE=0.
1050 NLICE=0.
1051 NSmICE=0.
1052 ELSE
1053 !
1054 !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
1055 ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
1056 ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
1057 !
1058 DUM=XMImax*EXP(.0536*TC)
1059 INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
1060 TOT_ICE=THICK*QI+BLEND*ASNOW
1061 PILOSS=-TOT_ICE/THICK
1062 LBEF=MAX(1,L-1)
1063 DUM1=RimeF_col(LBEF)
1064 DUM2=RimeF_col(L)
1065 RimeF1=(DUM2*THICK*QI+DUM1*BLEND*ASNOW)/TOT_ICE
1066 RimeF1=MIN(RimeF1, RFmax)
1067 DO IPASS=0,1
1068 IF (RimeF1 .LE. 1.) THEN
1069 RimeF1=1.
1070 VrimeF=1.
1071 ELSE
1072 IXS=MAX(2, MIN(INDEXS/100, 9))
1073 XRF=10.492*ALOG(RimeF1)
1074 IXRF=MAX(0, MIN(INT(XRF), Nrime))
1075 IF (IXRF .GE. Nrime) THEN
1076 VrimeF=VEL_RF(IXS,Nrime)
1077 ELSE
1078 VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
1079 & (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
1080 ENDIF
1081 ENDIF ! End IF (RimeF1 .LE. 1.)
1082 VEL_INC=GAMMAS*VrimeF
1083 VSNOW=VEL_INC*VSNOWI(INDEXS)
1084 EMAIRI=THICK+BLDTRH*VSNOW
1085 XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
1086 FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
1087 QTICE=TOT_ICE/EMAIRI
1088 QLICE=FLIMASS*QTICE
1089 NLICE=QLICE/XLIMASS
1090 NSmICE=Fsmall*NLICE
1091 !
1092 IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) &
1093 & .OR. IPASS.EQ.1) THEN
1094 EXIT
1095 ELSE
1096 !
1097 !--- Reduce excessive accumulation of ice at upper levels
1098 ! associated with strong grid-resolved ascent
1099 !
1100 !--- Force NLICE to be between NLImin and NLImax
1101 !
1102 DUM=MAX(NLImin, MIN(NLImax, NLICE) )
1103 XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
1104 IF (XLI .LE. MASSI(MDImin) ) THEN
1105 INDEXS=MDImin
1106 ELSE IF (XLI .LE. MASSI(450) ) THEN
1107 DLI=9.5885E5*XLI**.42066 ! DLI in microns
1108 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1109 ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
1110 DLI=3.9751E6*XLI**.49870 ! DLI in microns
1111 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1112 ELSE
1113 INDEXS=MDImax
1114 !
1115 !--- 8/22/01: Increase density of large ice if maximum limits
1116 ! are reached for number concentration (NLImax) and mean size
1117 ! (MDImax). Done to increase fall out of ice.
1118 !
1119 IF (DUM .GE. NLImax) &
1120 & RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
1121 ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
1122 ! WRITE(6,"(4(a12,g11.4,1x))")
1123 ! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
1124 ! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
1125 ! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
1126 ENDIF ! End IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) ...
1127 ENDDO ! End DO IPASS=0,1
1128 ENDIF ! End IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT)
1129 ENDIF ! End IF (ICE_logical)
1130 !
1131 !----------------------------------------------------------------------
1132 !--------------- Calculate individual processes -----------------------
1133 !----------------------------------------------------------------------
1134 !
1135 !--- Cloud water autoconversion to rain and collection by rain
1136 !
1137 IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
1138 !
1139 !--- QW0 could be modified based on land/sea properties,
1140 ! presence of convection, etc. This is why QAUT0 and CRAUT
1141 ! are passed into the subroutine as externally determined
1142 ! parameters. Can be changed in the future if desired.
1143 !
1144 QW0=QAUT0*RRHO
1145 PRAUT=MAX(0., QW-QW0)*CRAUT
1146 IF (QLICE .GT. EPSQ) THEN
1147 !
1148 !--- Collection of cloud water by large ice particles ("snow")
1149 ! PIACWI=PIACW for riming, PIACWI=0 for shedding
1150 !
1151 FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1)
1152 PIACW=FWS*QW
1153 IF (TC .LT. 0.) PIACWI=PIACW ! Large ice riming
1154 ENDIF ! End IF (QLICE .GT. EPSQ)
1155 ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
1156 !
1157 !----------------------------------------------------------------------
1158 !--- Loop around some of the ice-phase processes if no ice should be present
1159 !----------------------------------------------------------------------
1160 !
1161 IF (ICE_logical .EQV. .FALSE.) GO TO 20
1162 !
1163 !--- Now the pretzel logic of calculating ice deposition
1164 !
1165 IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ)) THEN
1166 !
1167 !--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
1168 ! Sources of ice due to nucleation and convective detrainment are
1169 ! either poorly understood, poorly resolved at typical NWP
1170 ! resolutions, or are not represented (e.g., no detrained
1171 ! condensate in BMJ Cu scheme).
1172 !
1173 PCOND=-QW
1174 DUM1=TK+XLV1*PCOND ! Updated (dummy) temperature (deg K)
1175 DUM2=WV+QW ! Updated (dummy) water vapor mixing ratio
1176 DUM=1000.*FPVS(DUM1) ! Updated (dummy) saturation vapor pressure w/r/t ice
1177 DUM=RHgrd*EPS*DUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
1178 IF (DUM2 .GT. DUM) PIDEP=DEPOSIT (PP, DUM1, DUM2)
1179 DWVi=0. ! Used only for debugging
1180 !
1181 ELSE IF (TC .LT. 0.) THEN
1182 !
1183 !--- These quantities are handy for ice deposition/sublimation
1184 ! PIDEP_max - max deposition or minimum sublimation to ice saturation
1185 !
1186 DENOMI=1.+XLS2*QSI*TK2
1187 DWVi=MIN(WVQW,QSW)-QSI
1188 PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
1189 IF (QTICE .GT. 0.) THEN
1190 !
1191 !--- Calculate ice deposition/sublimation
1192 ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1193 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1194 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1195 ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
1196 ! VENTI2 - m**2/s**.5 ; DIDEP - unitless
1197 !
1198 SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1199 ABI=1./(RHO*XLS3*QSI*TK2/THERM_COND+1./DIFFUS)
1200 !
1201 !--- VENTIL - Number concentration * ventilation factors for large ice
1202 !--- VENTIS - Number concentration * ventilation factors for small ice
1203 !
1204 !--- Variation in the number concentration of ice with time is not
1205 ! accounted for in these calculations (could be in the future).
1206 !
1207 VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE
1208 VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE
1209 DIDEP=ABI*(VENTIL+VENTIS)*DTPH
1210 !
1211 !--- Account for change in water vapor supply w/ time
1212 !
1213 IF (DIDEP .GE. Xratio)then
1214 DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
1215 endif
1216 IF (DWVi .GT. 0.) THEN
1217 PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
1218 ELSE IF (DWVi .LT. 0.) THEN
1219 PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
1220 ENDIF
1221 !
1222 ELSE IF (WVQW.GT.QSI .AND. TC.LE.T_ICE_init) THEN
1223 !
1224 !--- Ice nucleation in near water-saturated conditions. Ice crystal
1225 ! growth during time step calculated using Miller & Young (1979, JAS).
1226 !--- These deposition rates could drive conditions below water saturation,
1227 ! which is the basis of these calculations. Intended to approximate
1228 ! more complex & computationally intensive calculations.
1229 !
1230 INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
1231 !
1232 !--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
1233 !
1234 !--- DUM2 is the number of ice crystals nucleated at water-saturated
1235 ! conditions based on Meyers et al. (JAM, 1992).
1236 !
1237 !--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
1238 ! if DUM2 values are increased in future experiments
1239 !
1240 DUM1=QSW/QSI-1.
1241 DUM2=1.E3*EXP(12.96*DUM1-.639)
1242 PIDEP=MIN(PIDEP_max, DUM2*MY_GROWTH(INDEX_MY)*RRHO)
1243 !
1244 ENDIF ! End IF (QTICE .GT. 0.)
1245 !
1246 ENDIF ! End IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ))
1247 !
1248 !------------------------------------------------------------------------
1249 !
1250 20 CONTINUE ! Jump here if conditions for ice are not present
1251
1252
1253 !
1254 !------------------------------------------------------------------------
1255 !
1256 !--- Cloud water condensation
1257 !
1258 IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd)) THEN
1259 IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.) THEN
1260 PCOND=CONDENSE (PP, QW, TK, WV)
1261 ELSE
1262 !
1263 !--- Modify cloud condensation in response to ice processes
1264 !
1265 DUM=XLV*QSWgrd*RCPRV*TK2
1266 DENOMWI=1.+XLS*DUM
1267 DENOMF=XLF*DUM
1268 DUM=MAX(0., PIDEP)
1269 PCOND=(WV-QSWgrd-DENOMWI*DUM-DENOMF*PIACWI)/DENOMW
1270 DUM1=-QW
1271 DUM2=PCOND-PIACW
1272 IF (DUM2 .LT. DUM1) THEN
1273 !
1274 !--- Limit cloud water sinks
1275 !
1276 DUM=DUM1/DUM2
1277 PCOND=DUM*PCOND
1278 PIACW=DUM*PIACW
1279 PIACWI=DUM*PIACWI
1280 ENDIF ! End IF (DUM2 .LT. DUM1)
1281 ENDIF ! End IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.)
1282 ENDIF ! End IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd))
1283 !
1284 !--- Limit freezing of accreted rime to prevent temperature oscillations,
1285 ! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
1286 !
1287 TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI
1288 IF (TCC .GT. 0.) THEN
1289 PIACWI=0.
1290 TCC=TC+XLV1*PCOND+XLS1*PIDEP
1291 ENDIF
1292 IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
1293 !
1294 !--- Calculate melting and evaporation/condensation
1295 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1296 ! VENTIL - m**-2 ; VENTI1 - m ;
1297 ! VENTI2 - m**2/s**.5 ; CIEVP - /s
1298 !
1299 SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1300 VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))
1301 AIEVP=VENTIL*DIFFUS*DTPH
1302 IF (AIEVP .LT. Xratio) THEN
1303 DIEVP=AIEVP
1304 ELSE
1305 DIEVP=1.-EXP(-AIEVP)
1306 ENDIF
1307 QSW0=EPS*ESW0/(PP-ESW0)
1308 DWV0=MIN(WV,QSW)-QSW0
1309 DUM=QW+PCOND
1310 IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
1311 !
1312 !--- Evaporation from melting snow (sink of snow) or shedding
1313 ! of water condensed onto melting snow (source of rain)
1314 !
1315 DUM=DWV0*DIEVP
1316 PIEVP=MAX( MIN(0., DUM), PILOSS)
1317 PICND=MAX(0., DUM)
1318 ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
1319 PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF
1320 !
1321 !--- Limit melting to prevent temperature oscillations across 0C
1322 !
1323 DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
1324 PIMLT=MIN(PIMLT, DUM1)
1325 !
1326 !--- Limit loss of snow by melting (>0) and evaporation
1327 !
1328 DUM=PIEVP-PIMLT
1329 IF (DUM .LT. PILOSS) THEN
1330 DUM1=PILOSS/DUM
1331 PIMLT=PIMLT*DUM1
1332 PIEVP=PIEVP*DUM1
1333 ENDIF ! End IF (DUM .GT. QTICE)
1334 ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
1335 !
1336 !--- IMPORTANT: Estimate time-averaged properties.
1337 !
1338 ! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
1339 ! the total mass of rain in the current layer and the input
1340 ! flux of rain from above
1341 ! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction)
1342 ! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
1343 ! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
1344 ! above and the rain within the layer
1345 ! * RQR - rain content (kg/m**3)
1346 ! * INDEXR - mean size of rain drops to the nearest 1 micron in size
1347 ! * N0r - intercept of rain size distribution (typically 10**6 m**-4)
1348 !
1349 TOT_RAIN=0.
1350 VRAIN1=0.
1351 QTRAIN=0.
1352 PRLOSS=0.
1353 RQR=0.
1354 N0r=0.
1355 INDEXR1=0
1356 INDEXR=MDRmin
1357 IF (RAIN_logical) THEN
1358 IF (ARAIN .LE. 0.) THEN
1359 INDEXR=MDRmin
1360 VRAIN1=0.
1361 ELSE
1362 !
1363 !--- INDEXR (related to mean diameter) & N0r could be modified
1364 ! by land/sea properties, presence of convection, etc.
1365 !
1366 !--- Rain rate normalized to a density of 1.194 kg/m**3
1367 !
1368 RR=ARAIN/(DTPH*GAMMAR)
1369 !
1370 IF (RR .LE. RR_DRmin) THEN
1371 !
1372 !--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
1373 ! instead vary N0r with rain rate
1374 !
1375 INDEXR=MDRmin
1376 ELSE IF (RR .LE. RR_DR1) THEN
1377 !
1378 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1379 ! for mean diameters (Dr) between 0.05 and 0.10 mm:
1380 ! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
1381 ! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
1382 ! RR in kg/(m**2*s)
1383 ! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
1384 !
1385 INDEXR=INT( 1.123E3*RR**.1947 + .5 )
1386 INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
1387 ELSE IF (RR .LE. RR_DR2) THEN
1388 !
1389 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1390 ! for mean diameters (Dr) between 0.10 and 0.20 mm:
1391 ! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
1392 ! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
1393 ! RR in kg/(m**2*s)
1394 ! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
1395 !
1396 INDEXR=INT( 1.225E3*RR**.2017 + .5 )
1397 INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
1398 ELSE IF (RR .LE. RR_DR3) THEN
1399 !
1400 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1401 ! for mean diameters (Dr) between 0.20 and 0.32 mm:
1402 ! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
1403 ! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
1404 ! RR in kg/(m**2*s)
1405 ! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
1406 !
1407 INDEXR=INT( 1.3006E3*RR**.2083 + .5 )
1408 INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
1409 ELSE IF (RR .LE. RR_DRmax) THEN
1410 !
1411 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1412 ! for mean diameters (Dr) between 0.32 and 0.45 mm:
1413 ! V(Dr)=944.8*Dr**.6636, V in m/s and Dr in m
1414 ! RR = PI*1000.*N0r0*944.8*Dr**(4+.6636) = 2.3745e13*Dr**4.6636,
1415 ! RR in kg/(m**2*s)
1416 ! Dr (m) = 1.355e-3*RR**.2144 -> Dr (microns) = 1.355e3*RR**.2144
1417 !
1418 INDEXR=INT( 1.355E3*RR**.2144 + .5 )
1419 INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
1420 ELSE
1421 !
1422 !--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
1423 ! instead vary N0r with rain rate
1424 !
1425 INDEXR=MDRmax
1426 ENDIF ! End IF (RR .LE. RR_DRmin) etc.
1427 VRAIN1=GAMMAR*VRAIN(INDEXR)
1428 ENDIF ! End IF (ARAIN .LE. 0.)
1429 INDEXR1=INDEXR ! For debugging only
1430 TOT_RAIN=THICK*QR+BLEND*ARAIN
1431 QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
1432 PRLOSS=-TOT_RAIN/THICK
1433 RQR=RHO*QTRAIN
1434 !
1435 !--- RQR - time-averaged rain content (kg/m**3)
1436 !
1437 IF (RQR .LE. RQR_DRmin) THEN
1438 N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
1439 INDEXR=MDRmin
1440 ELSE IF (RQR .GE. RQR_DRmax) THEN
1441 N0r=CN0r_DMRmax*RQR
1442 INDEXR=MDRmax
1443 ELSE
1444 N0r=N0r0
1445 INDEXR=MAX( XMRmin, MIN(CN0r0*RQR**.25, XMRmax) )
1446 ENDIF
1447 !
1448 IF (TC .LT. T_ICE) THEN
1449 PIACR=-PRLOSS
1450 ELSE
1451 DWVr=WV-PCOND-QSW
1452 DUM=QW+PCOND
1453 IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
1454 !
1455 !--- Rain evaporation
1456 !
1457 ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1458 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1459 !
1460 ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ;
1461 ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ;
1462 ! CREVP - unitless
1463 !
1464 RFACTOR=GAMMAR**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1465 ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
1466 !
1467 !--- Note that VENTR1, VENTR2 lookup tables do not include the
1468 ! 1/Davg multiplier as in the ice tables
1469 !
1470 VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR))
1471 CREVP=ABW*VENTR*DTPH
1472 IF (CREVP .LT. Xratio) THEN
1473 DUM=DWVr*CREVP
1474 ELSE
1475 DUM=DWVr*(1.-EXP(-CREVP*DENOMW))/DENOMW
1476 ENDIF
1477 PREVP=MAX(DUM, PRLOSS)
1478 ELSE IF (QW .GT. EPSQ) THEN
1479 FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR)
1480 PRACW=MIN(1.,FWR)*QW
1481 ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
1482 !
1483 IF (TC.LT.0. .AND. TCC.LT.0.) THEN
1484 !
1485 !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
1486 ! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
1487 !
1488 DUM=.001*FLOAT(INDEXR)
1489 DUM=(EXP(ABFR*TC)-1.)*DUM*DUM*DUM*DUM*DUM*DUM*DUM
1490 PIACR=MIN(CBFR*N0r*RRHO*DUM, QTRAIN)
1491 IF (QLICE .GT. EPSQ) THEN
1492 !
1493 !--- Freezing of rain by collisions w/ large ice
1494 !
1495 DUM=GAMMAR*VRAIN(INDEXR)
1496 DUM1=DUM-VSNOW
1497 !
1498 !--- DUM2 - Difference in spectral fall speeds of rain and
1499 ! large ice, parameterized following eq. (48) on p. 112 of
1500 ! Murakami (J. Meteor. Soc. Japan, 1990)
1501 !
1502 DUM2=(DUM1*DUM1+.04*DUM*VSNOW)**.5
1503 DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS &
1504 & +.5E-12*INDEXS*INDEXS
1505 FIR=MIN(1., CIACR*NLICE*DUM1*DUM2)
1506 !
1507 !--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
1508 !
1509 PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
1510 ENDIF ! End IF (QLICE .GT. EPSQ)
1511 DUM=PREVP-PIACR
1512 If (DUM .LT. PRLOSS) THEN
1513 DUM1=PRLOSS/DUM
1514 PREVP=DUM1*PREVP
1515 PIACR=DUM1*PIACR
1516 ENDIF ! End If (DUM .LT. PRLOSS)
1517 ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
1518 ENDIF ! End IF (TC .LT. T_ICE)
1519 ENDIF ! End IF (RAIN_logical)
1520 !
1521 !----------------------------------------------------------------------
1522 !---------------------- Main Budget Equations -------------------------
1523 !----------------------------------------------------------------------
1524 !
1525 !
1526 !-----------------------------------------------------------------------
1527 !--- Update fields, determine characteristics for next lower layer ----
1528 !-----------------------------------------------------------------------
1529 !
1530 !--- Carefully limit sinks of cloud water
1531 !
1532 DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
1533 IF (DUM1 .GT. QW) THEN
1534 DUM=QW/DUM1
1535 PIACW=DUM*PIACW
1536 PIACWI=DUM*PIACWI
1537 PRAUT=DUM*PRAUT
1538 PRACW=DUM*PRACW
1539 IF (PCOND .LT. 0.) PCOND=DUM*PCOND
1540 ENDIF
1541 PIACWR=PIACW-PIACWI ! TC >= 0C
1542 !
1543 !--- QWnew - updated cloud water mixing ratio
1544 !
1545 DELW=PCOND-PIACW-PRAUT-PRACW
1546 QWnew=QW+DELW
1547 IF (QWnew .LE. EPSQ) QWnew=0.
1548 IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
1549 DUM=QWnew/QW
1550 IF (DUM .LT. TOLER) QWnew=0.
1551 ENDIF
1552 !
1553 !--- Update temperature and water vapor mixing ratios
1554 !
1555 DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
1556 & +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT)
1557 Tnew=TK+DELT
1558 !
1559 DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
1560 WVnew=WV+DELV
1561 !
1562 !--- Update ice mixing ratios
1563 !
1564 !---
1565 ! * TOT_ICEnew - total mass (small & large) ice after microphysics,
1566 ! which is the sum of the total mass of large ice in the
1567 ! current layer and the flux of ice out of the grid box below
1568 ! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
1569 ! rimed) ice mass to the unrimed ice mass (>=1)
1570 ! * QInew - updated mixing ratio of total (large & small) ice in layer
1571 ! -> TOT_ICEnew=QInew*THICK+BLDTRH*QLICEnew*VSNOW
1572 ! -> But QLICEnew=QInew*FLIMASS, so
1573 ! -> TOT_ICEnew=QInew*(THICK+BLDTRH*FLIMASS*VSNOW)
1574 ! * ASNOWnew - updated accumulation of snow at bottom of grid cell
1575 !---
1576 !
1577 DELI=0.
1578 RimeF=1.
1579 IF (ICE_logical) THEN
1580 DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
1581 TOT_ICEnew=TOT_ICE+THICK*DELI
1582 IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
1583 DUM=TOT_ICEnew/TOT_ICE
1584 IF (DUM .LT. TOLER) TOT_ICEnew=0.
1585 ENDIF
1586 IF (TOT_ICEnew .LE. CLIMIT) THEN
1587 TOT_ICEnew=0.
1588 RimeF=1.
1589 QInew=0.
1590 ASNOWnew=0.
1591 ELSE
1592 !
1593 !--- Update rime factor if appropriate
1594 !
1595 DUM=PIACWI+PIACR
1596 IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ) THEN
1597 RimeF=RimeF1
1598 ELSE
1599 !
1600 !--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
1601 ! DUM1 - Total ice mass, rimed & unrimed
1602 ! DUM2 - Estimated mass of *unrimed* ice
1603 !
1604 DUM1=TOT_ICE+THICK*(PIDEP+DUM)
1605 DUM2=TOT_ICE/RimeF1+THICK*PIDEP
1606 IF (DUM2 .LE. 0.) THEN
1607 RimeF=RFmax
1608 ELSE
1609 RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
1610 ENDIF
1611 ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
1612 QInew=TOT_ICEnew/(THICK+BLDTRH*FLIMASS*VSNOW)
1613 IF (QInew .LE. EPSQ) QInew=0.
1614 IF (QI.GT.0. .AND. QInew.NE.0.) THEN
1615 DUM=QInew/QI
1616 IF (DUM .LT. TOLER) QInew=0.
1617 ENDIF
1618 ASNOWnew=BLDTRH*FLIMASS*VSNOW*QInew
1619 IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
1620 DUM=ASNOWnew/ASNOW
1621 IF (DUM .LT. TOLER) ASNOWnew=0.
1622 ENDIF
1623 ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
1624 ENDIF ! End IF (ICE_logical)
1625
1626
1627 !
1628 !--- Update rain mixing ratios
1629 !
1630 !---
1631 ! * TOT_RAINnew - total mass of rain after microphysics
1632 ! current layer and the input flux of ice from above
1633 ! * VRAIN2 - time-averaged fall speed of rain in grid and below
1634 ! (with air resistance correction)
1635 ! * QRnew - updated rain mixing ratio in layer
1636 ! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
1637 ! * ARAINnew - updated accumulation of rain at bottom of grid cell
1638 !---
1639 !
1640 DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
1641 TOT_RAINnew=TOT_RAIN+THICK*DELR
1642 IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
1643 DUM=TOT_RAINnew/TOT_RAIN
1644 IF (DUM .LT. TOLER) TOT_RAINnew=0.
1645 ENDIF
1646 IF (TOT_RAINnew .LE. CLIMIT) THEN
1647 TOT_RAINnew=0.
1648 VRAIN2=0.
1649 QRnew=0.
1650 ARAINnew=0.
1651 ELSE
1652 !
1653 !--- 1st guess time-averaged rain rate at bottom of grid box
1654 !
1655 RR=TOT_RAINnew/(DTPH*GAMMAR)
1656 !
1657 !--- Use same algorithm as above for calculating mean drop diameter
1658 ! (IDR, in microns), which is used to estimate the time-averaged
1659 ! fall speed of rain drops at the bottom of the grid layer. This
1660 ! isn't perfect, but the alternative is solving a transcendental
1661 ! equation that is numerically inefficient and nasty to program
1662 ! (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
1663 !
1664 IF (RR .LE. RR_DRmin) THEN
1665 IDR=MDRmin
1666 ELSE IF (RR .LE. RR_DR1) THEN
1667 IDR=INT( 1.123E3*RR**.1947 + .5 )
1668 IDR=MAX( MDRmin, MIN(IDR, MDR1) )
1669 ELSE IF (RR .LE. RR_DR2) THEN
1670 IDR=INT( 1.225E3*RR**.2017 + .5 )
1671 IDR=MAX( MDR1, MIN(IDR, MDR2) )
1672 ELSE IF (RR .LE. RR_DR3) THEN
1673 IDR=INT( 1.3006E3*RR**.2083 + .5 )
1674 IDR=MAX( MDR2, MIN(IDR, MDR3) )
1675 ELSE IF (RR .LE. RR_DRmax) THEN
1676 IDR=INT( 1.355E3*RR**.2144 + .5 )
1677 IDR=MAX( MDR3, MIN(IDR, MDRmax) )
1678 ELSE
1679 IDR=MDRmax
1680 ENDIF ! End IF (RR .LE. RR_DRmin)
1681 VRAIN2=GAMMAR*VRAIN(IDR)
1682 QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
1683 IF (QRnew .LE. EPSQ) QRnew=0.
1684 IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
1685 DUM=QRnew/QR
1686 IF (DUM .LT. TOLER) QRnew=0.
1687 ENDIF
1688 ARAINnew=BLDTRH*VRAIN2*QRnew
1689 IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
1690 DUM=ARAINnew/ARAIN
1691 IF (DUM .LT. TOLER) ARAINnew=0.
1692 ENDIF
1693 ENDIF
1694 !
1695 WCnew=QWnew+QRnew+QInew
1696 !
1697 !----------------------------------------------------------------------
1698 !-------------- Begin debugging & verification ------------------------
1699 !----------------------------------------------------------------------
1700 !
1701 !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
1702 !
1703
1704
1705 QT=THICK*(WV+WC)+ARAIN+ASNOW
1706 QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
1707 BUDGET=QT-QTnew
1708 !
1709 !--- Additional check on budget preservation, accounting for truncation effects
1710 !
1711 DBG_logical=.FALSE.
1712 ! DUM=ABS(BUDGET)
1713 ! IF (DUM .GT. TOLER) THEN
1714 ! DUM=DUM/MIN(QT, QTnew)
1715 ! IF (DUM .GT. TOLER) DBG_logical=.TRUE.
1716 ! ENDIF
1717 !!
1718 ! DUM=(RHgrd+.001)*QSInew
1719 ! IF ( (QWnew.GT.EPSQ) .OR. QRnew.GT.EPSQ .OR. WVnew.GT.DUM)
1720 ! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE.
1721 !
1722 ! IF (TC.GT.5. .AND. QInew.GT.EPSQ) DBG_logical=.TRUE.
1723 !
1724 IF ((WVnew.LT.EPSQ .OR. DBG_logical) .AND. PRINT_diag) THEN
1725 !
1726 WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index,&
1727 & ' L=',L,' LSFC=',LSFC
1728 !
1729 ESW=1000.*FPVS0(Tnew)
1730 QSWnew=EPS*ESW/(PP-ESW)
1731 IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
1732 ESI=1000.*FPVS(Tnew)
1733 QSInew=EPS*ESI/(PP-ESI)
1734 ELSE
1735 QSI=QSW
1736 QSInew=QSWnew
1737 ENDIF
1738 WSnew=QSInew
1739 WRITE(6,"(4(a12,g11.4,1x))") &
1740 & '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
1741 & '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
1742 & 'RHgrd=',RHgrd, &
1743 & '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
1744 & 'RHInew=',WVnew/QSInew, &
1745 & '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
1746 & '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
1747 & '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
1748 & '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
1749 & '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
1750 & 'ASNOWnew=',ASNOWnew, &
1751 & '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
1752 & 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
1753 & '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
1754 !
1755 WRITE(6,"(4(a12,g11.4,1x))") &
1756 & '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
1757 & '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
1758 & '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
1759 & '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
1760 & PIMLT, &
1761 & '{} PIACR=',PIACR
1762 !
1763 IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
1764 & '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
1765 & 'VSNOW=',VSNOW, &
1766 & '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
1767 & 'FLIMASS=',FLIMASS, &
1768 & '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
1769 & 'QTICE=',QTICE, &
1770 & '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
1771 & 'EMAIRI=',EMAIRI, &
1772 & '{} RimeF=',RimeF
1773 !
1774 IF (TOT_RAIN.GT.0. .OR. TOT_RAINnew.GT.0.) &
1775 & WRITE(6,"(4(a12,g11.4,1x))") &
1776 & '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
1777 & 'GAMMAR=',GAMMAR,'N0r=',N0r, &
1778 & '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
1779 & '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
1780 & 'VOLR2=',THICK+BLDTRH*VRAIN2
1781 !
1782 IF (PRAUT .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
1783 !
1784 IF (PRACW .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
1785 !
1786 IF (PIACR .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
1787 !
1788 DUM=PIMLT+PICND-PREVP-PIEVP
1789 IF (DUM.GT.0. .or. DWVi.NE.0.) &
1790 & WRITE(6,"(4(a12,g11.4,1x))") &
1791 & '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
1792 & 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
1793 !
1794 IF (PREVP .LT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1795 & '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
1796 & '{} DWVr=',DWVr,'DENOMW=',DENOMW
1797 !
1798 IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
1799 & WRITE(6,"(4(a12,g11.4,1x))") &
1800 & '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
1801 & 'SFACTOR=',SFACTOR, &
1802 & '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1803 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1804 & '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
1805 !
1806 IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
1807 & WRITE(6,"(4(a12,g11.4,1x))") &
1808 & '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
1809 & 'DUM2=',PCOND-PIACW
1810 !
1811 IF (FWS .GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1812 & '{} FWS=',FWS
1813 !
1814 DUM=PIMLT+PICND-PIEVP
1815 IF (DUM.GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1816 & '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1817 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1818 & '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
1819 !
1820 ENDIF
1821
1822
1823 !
1824 !-----------------------------------------------------------------------
1825 !--------------- Water budget statistics & maximum values --------------
1826 !-----------------------------------------------------------------------
1827 !
1828 IF (PRINT_diag) THEN
1829 ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
1830 IF (QInew .GT. EPSQ) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
1831 IF (QInew.GT.EPSQ .AND. QRnew+QWnew.GT.EPSQ) &
1832 & NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
1833 IF (QWnew .GT. EPSQ) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
1834 IF (QRnew .GT. EPSQ) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
1835 !
1836 QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
1837 QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
1838 QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
1839 QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
1840 QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
1841 QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
1842 QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
1843 QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
1844 !
1845 QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
1846 QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
1847 QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
1848 QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
1849 QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
1850 QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
1851 QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
1852 QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
1853 QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
1854 QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
1855 QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
1856 QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
1857 !
1858 QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
1859 QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
1860 QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
1861 QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
1862 QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
1863 QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
1864 IF (QInew .GT. 0.) &
1865 & QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
1866 !
1867 ENDIF
1868 !
1869 !----------------------------------------------------------------------
1870 !------------------------- Update arrays ------------------------------
1871 !----------------------------------------------------------------------
1872 !
1873
1874
1875 T_col(L)=Tnew ! Updated temperature
1876 !
1877 QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity
1878 WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio
1879 QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio
1880 QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio
1881 QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio
1882 RimeF_col(L)=RimeF ! Updated rime factor
1883 ASNOW=ASNOWnew ! Updated accumulated snow
1884 ARAIN=ARAINnew ! Updated accumulated rain
1885 !
1886 !#######################################################################
1887 !
1888 10 CONTINUE ! ##### End "L" loop through model levels #####
1889
1890
1891 !
1892 !#######################################################################
1893 !
1894 !-----------------------------------------------------------------------
1895 !--------------------------- Return to GSMDRIVE -----------------------
1896 !-----------------------------------------------------------------------
1897 !
1898 CONTAINS
1899 !#######################################################################
1900 !--------- Produces accurate calculation of cloud condensation ---------
1901 !#######################################################################
1902 !
1903 REAL FUNCTION CONDENSE (PP, QW, TK, WV)
1904 !
1905 !---------------------------------------------------------------------------------
1906 !------ The Asai (1965) algorithm takes into consideration the release of ------
1907 !------ latent heat in increasing the temperature & in increasing the ------
1908 !------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
1909 !---------------------------------------------------------------------------------
1910 !
1911 IMPLICIT NONE
1912 !
1913 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1914 REAL (KIND=HIGH_PRES), PARAMETER :: &
1915 & RHLIMIT=.001, RHLIMIT1=-RHLIMIT
1916 REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
1917 !
1918 REAL,INTENT(IN) :: QW,PP,WV,TK
1919 REAL WVdum,Tdum,XLV2,DWV,WS,ESW,XLV1,XLV
1920 integer nsteps
1921 !
1922 !-----------------------------------------------------------------------
1923 !
1924 !--- LV (T) is from Bolton (JAS, 1980)
1925 !
1926 XLV=3.148E6-2370.*TK
1927 XLV1=XLV*RCP
1928 XLV2=XLV*XLV*RCPRV
1929 Tdum=TK
1930 WVdum=WV
1931 WCdum=QW
1932 ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water
1933 WS=RHgrd*EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
1934 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1935 SSAT=DWV/WS ! Supersaturation ratio
1936 CONDENSE=0.
1937 nsteps = 0
1938 DO WHILE ((SSAT.LT.RHLIMIT1 .AND. WCdum.GT.EPSQ) &
1939 & .OR. SSAT.GT.RHLIMIT)
1940 nsteps = nsteps + 1
1941 COND=DWV/(1.+XLV2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
1942 COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
1943 Tdum=Tdum+XLV1*COND ! Updated temperature
1944 WVdum=WVdum-COND ! Updated water vapor mixing ratio
1945 WCdum=WCdum+COND ! Updated cloud water mixing ratio
1946 CONDENSE=CONDENSE+COND ! Total cloud water condensation
1947 ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water
1948 WS=RHgrd*EPS*ESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
1949 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1950 SSAT=DWV/WS ! Grid-scale supersaturation ratio
1951 ENDDO
1952 !
1953 END FUNCTION CONDENSE
1954 !
1955 !#######################################################################
1956 !---------------- Calculate ice deposition at T<T_ICE ------------------
1957 !#######################################################################
1958 !
1959 REAL FUNCTION DEPOSIT (PP, Tdum, WVdum)
1960 !
1961 !--- Also uses the Asai (1965) algorithm, but uses a different target
1962 ! vapor pressure for the adjustment
1963 !
1964 IMPLICIT NONE
1965 !
1966 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1967 REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
1968 & RHLIMIT1=-RHLIMIT
1969 REAL (KIND=HIGH_PRES) :: DEP, SSAT
1970 !
1971 real,INTENT(IN) :: PP
1972 real,INTENT(INOUT) :: WVdum,Tdum
1973 real ESI,WS,DWV
1974 !
1975 !-----------------------------------------------------------------------
1976 !
1977 ESI=1000.*FPVS(Tdum) ! Saturation vapor press w/r/t ice
1978 WS=RHgrd*EPS*ESI/(PP-ESI) ! Saturation mixing ratio
1979 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1980 SSAT=DWV/WS ! Supersaturation ratio
1981 DEPOSIT=0.
1982 DO WHILE (SSAT.GT.RHLIMIT .OR. SSAT.LT.RHLIMIT1)
1983 !
1984 !--- Note that XLVS2=LS*LV/(CP*RV)=LV*WS/(RV*T*T)*(LS/CP*DEP1),
1985 ! where WS is the saturation mixing ratio following Clausius-
1986 ! Clapeyron (see Asai,1965; Young,1993,p.405)
1987 !
1988 DEP=DWV/(1.+XLS2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
1989 Tdum=Tdum+XLS1*DEP ! Updated temperature
1990 WVdum=WVdum-DEP ! Updated ice mixing ratio
1991 DEPOSIT=DEPOSIT+DEP ! Total ice deposition
1992 ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice
1993 WS=RHgrd*EPS*ESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
1994 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1995 SSAT=DWV/WS ! Grid-scale supersaturation ratio
1996 ENDDO
1997 !
1998 END FUNCTION DEPOSIT
1999 !
2000 END SUBROUTINE EGCP01COLUMN
2001
2002 !#######################################################################
2003 !------- Initialize constants & lookup tables for microphysics ---------
2004 !#######################################################################
2005 !
2006
2007 ! SH 0211/2002
2008
2009 !-----------------------------------------------------------------------
2010 SUBROUTINE etanewinit (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
2011 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
2012 & MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
2013 & ALLOWED_TO_READ, &
2014 & IDS,IDE,JDS,JDE,KDS,KDE, &
2015 & IMS,IME,JMS,JME,KMS,KME, &
2016 & ITS,ITE,JTS,JTE,KTS,KTE )
2017 !-----------------------------------------------------------------------
2018 !-------------------------------------------------------------------------------
2019 !--- SUBPROGRAM DOCUMENTATION BLOCK
2020 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
2021 ! Jin ORG: W/NP22 DATE: January 2002
2022 ! (Modification for WRF structure)
2023 !
2024 !-------------------------------------------------------------------------------
2025 ! ABSTRACT:
2026 ! * Reads various microphysical lookup tables used in COLUMN_MICRO
2027 ! * Lookup tables were created "offline" and are read in during execution
2028 ! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
2029 !-------------------------------------------------------------------------------
2030 !
2031 ! USAGE: CALL etanewinit FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
2032 !
2033 ! INPUT ARGUMENT LIST:
2034 ! DTPH - physics time step (s)
2035 !
2036 ! OUTPUT ARGUMENT LIST:
2037 ! NONE
2038 !
2039 ! OUTPUT FILES:
2040 ! NONE
2041 !
2042 ! SUBROUTINES:
2043 ! MY_GROWTH_RATES - lookup table for growth of nucleated ice
2044 ! GPVS - lookup tables for saturation vapor pressure (water, ice)
2045 !
2046 ! UNIQUE: NONE
2047 !
2048 ! LIBRARY: NONE
2049 !
2050 ! COMMON BLOCKS:
2051 ! CMICRO_CONS - constants used in GSMCOLUMN
2052 ! CMY600 - lookup table for growth of ice crystals in
2053 ! water saturated conditions (Miller & Young, 1979)
2054 ! IVENT_TABLES - lookup tables for ventilation effects of ice
2055 ! IACCR_TABLES - lookup tables for accretion rates of ice
2056 ! IMASS_TABLES - lookup tables for mass content of ice
2057 ! IRATE_TABLES - lookup tables for precipitation rates of ice
2058 ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
2059 ! MAPOT - Need lat/lon grid resolution
2060 ! RVENT_TABLES - lookup tables for ventilation effects of rain
2061 ! RACCR_TABLES - lookup tables for accretion rates of rain
2062 ! RMASS_TABLES - lookup tables for mass content of rain
2063 ! RVELR_TABLES - lookup tables for fall speeds of rain
2064 ! RRATE_TABLES - lookup tables for precipitation rates of rain
2065 !
2066 ! ATTRIBUTES:
2067 ! LANGUAGE: FORTRAN 90
2068 ! MACHINE : IBM SP
2069 !
2070 !-----------------------------------------------------------------------
2071 !
2072 !
2073 !-----------------------------------------------------------------------
2074 IMPLICIT NONE
2075 !-----------------------------------------------------------------------
2076 !-------------------------------------------------------------------------
2077 !-------------- Parameters & arrays for lookup tables --------------------
2078 !-------------------------------------------------------------------------
2079 !
2080 !--- Common block of constants used in column microphysics
2081 !
2082 !WRF
2083 ! real DLMD,DPHD
2084 !WRF
2085 !
2086 !-----------------------------------------------------------------------
2087 !--- Parameters & data statement for local calculations
2088 !-----------------------------------------------------------------------
2089 !
2090 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
2091 !
2092 ! VARIABLES PASSED IN
2093 integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2094 & ,IMS,IME,JMS,JME,KMS,KME &
2095 & ,ITS,ITE,JTS,JTE,KTS,KTE
2096 !WRF
2097 INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
2098 !
2099 real, INTENT(IN) :: DELX,DELY
2100 real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
2101 real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
2102 real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT) :: &
2103 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
2104 INTEGER, PARAMETER :: ITLO=-60, ITHI=40
2105 ! integer,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
2106 ! real,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
2107 ! real,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
2108 ! real,INTENT(INOUT) :: PRECtot(2),PRECmax(2)
2109 real,INTENT(IN) :: DT,GSMDT
2110 LOGICAL,INTENT(IN) :: allowed_to_read,restart
2111 !
2112 !-----------------------------------------------------------------------
2113 ! LOCAL VARIABLES
2114 !-----------------------------------------------------------------------
2115 REAL :: BBFR,DTPH,PI,DX,Thour_print
2116 INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
2117 INTEGER :: etampnew_unit1
2118 LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
2119 LOGICAL :: opened
2120 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
2121 CHARACTER*80 errmess
2122 !
2123 !-----------------------------------------------------------------------
2124 !
2125 JTF=MIN0(JTE,JDE-1)
2126 KTF=MIN0(KTE,KDE-1)
2127 ITF=MIN0(ITE,IDE-1)
2128 !
2129 DO J=JTS,JTF
2130 DO I=ITS,ITF
2131 LOWLYR(I,J)=1
2132 ENDDO
2133 ENDDO
2134 !
2135 IF(.NOT.RESTART)THEN
2136 DO J = jts,jte
2137 DO K = kts,kte
2138 DO I= its,ite
2139 F_ICE_PHY(i,k,j)=0.
2140 F_RAIN_PHY(i,k,j)=0.
2141 F_RIMEF_PHY(i,k,j)=1.
2142 ENDDO
2143 ENDDO
2144 ENDDO
2145 ENDIF
2146 !
2147 !-----------------------------------------------------------------------
2148 IF(ALLOWED_TO_READ)THEN
2149 !-----------------------------------------------------------------------
2150 !
2151 DX=((DELX)**2+(DELY)**2)**.5/1000. ! Model resolution at equator (km)
2152 DX=MIN(100., MAX(5., DX) )
2153 !
2154 !-- Relative humidity threshold for the onset of grid-scale condensation
2155 !!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
2156 !! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
2157 ! RHgrd=0.90+.08*((100.-DX)/95.)**.5
2158 !
2159 DTPH=MAX(GSMDT*60.,DT)
2160 DTPH=NINT(DTPH/DT)*DT
2161 !
2162 !--- Create lookup tables for saturation vapor pressure w/r/t water & ice
2163 !
2164 CALL GPVS
2165 !
2166 !--- Read in various lookup tables
2167 !
2168 IF ( wrf_dm_on_monitor() ) THEN
2169 DO i = 31,99
2170 INQUIRE ( i , OPENED = opened )
2171 IF ( .NOT. opened ) THEN
2172 etampnew_unit1 = i
2173 GOTO 2061
2174 ENDIF
2175 ENDDO
2176 etampnew_unit1 = -1
2177 2061 CONTINUE
2178 ENDIF
2179 !
2180 CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
2181 !
2182 IF ( etampnew_unit1 < 0 ) THEN
2183 CALL wrf_error_fatal ( 'module_mp_etanew: etanewinit: Can not find '// &
2184 'unused fortran unit to read in lookup table.' )
2185 ENDIF
2186 !
2187 IF ( wrf_dm_on_monitor() ) THEN
2188 !!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
2189 OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
2190 & FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
2191 !
2192 READ(etampnew_unit1) VENTR1
2193 READ(etampnew_unit1) VENTR2
2194 READ(etampnew_unit1) ACCRR
2195 READ(etampnew_unit1) MASSR
2196 READ(etampnew_unit1) VRAIN
2197 READ(etampnew_unit1) RRATE
2198 READ(etampnew_unit1) VENTI1
2199 READ(etampnew_unit1) VENTI2
2200 READ(etampnew_unit1) ACCRI
2201 READ(etampnew_unit1) MASSI
2202 READ(etampnew_unit1) VSNOWI
2203 READ(etampnew_unit1) VEL_RF
2204 ! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
2205 CLOSE (etampnew_unit1)
2206 ENDIF
2207 !
2208 CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
2209 CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
2210 CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
2211 CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
2212 CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
2213 CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
2214 CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
2215 CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
2216 CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
2217 CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
2218 CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
2219 CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
2220 !
2221 !--- Calculates coefficients for growth rates of ice nucleated in water
2222 ! saturated conditions, scaled by physics time step (lookup table)
2223 !
2224 CALL MY_GROWTH_RATES (DTPH)
2225 ! CALL MY_GROWTH_RATES (DTPH,MY_GROWTH)
2226 !
2227 PI=ACOS(-1.)
2228 !
2229 !--- Constants associated with Biggs (1953) freezing of rain, as parameterized
2230 ! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
2231 !
2232 ABFR=-0.66
2233 BBFR=100.
2234 CBFR=20.*PI*PI*BBFR*RHOL*1.E-21
2235 !
2236 !--- CIACW is used in calculating riming rates
2237 ! The assumed effective collection efficiency of cloud water rimed onto
2238 ! ice is =0.5 below:
2239 !
2240 CIACW=DTPH*0.25*PI*0.5*(1.E5)**C1
2241 !
2242 !--- CIACR is used in calculating freezing of rain colliding with large ice
2243 ! The assumed collection efficiency is 1.0
2244 !
2245 CIACR=PI*DTPH
2246 !
2247 !--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
2248 ! * Four different functional relationships of mean drop diameter as
2249 ! a function of rain rate (RR), derived based on simple fits to
2250 ! mass-weighted fall speeds of rain as functions of mean diameter
2251 ! from the lookup tables.
2252 !
2253 RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
2254 RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
2255 RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
2256 RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
2257 RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of .45 mm
2258 !
2259 RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
2260 RQR_DR1=N0r0*MASSR(MDR1) ! Rain content for mean drop diameter of .10 mm
2261 RQR_DR2=N0r0*MASSR(MDR2) ! Rain content for mean drop diameter of .20 mm
2262 RQR_DR3=N0r0*MASSR(MDR3) ! Rain content for mean drop diameter of .32 mm
2263 RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of .45 mm
2264 C_N0r0=PI*RHOL*N0r0
2265 CN0r0=1.E6/C_N0r0**.25
2266 CN0r_DMRmin=1./(PI*RHOL*DMRmin**4)
2267 CN0r_DMRmax=1./(PI*RHOL*DMRmax**4)
2268 !
2269 !--- CRACW is used in calculating collection of cloud water by rain (an
2270 ! assumed collection efficiency of 1.0)
2271 !
2272 CRACW=DTPH*0.25*PI*1.0
2273 !
2274 ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
2275 RFmax=1.1**Nrime ! Maximum rime factor allowed
2276 !
2277 !------------------------------------------------------------------------
2278 !--------------- Constants passed through argument list -----------------
2279 !------------------------------------------------------------------------
2280 !
2281 !--- Important parameters for self collection (autoconversion) of
2282 ! cloud water to rain.
2283 !
2284 !--- CRAUT is proportional to the rate that cloud water is converted by
2285 ! self collection to rain (autoconversion rate)
2286 !
2287 CRAUT=1.-EXP(-1.E-3*DTPH)
2288 !
2289 !--- QAUT0 is the threshold cloud content for autoconversion to rain
2290 ! needed for droplets to reach a diameter of 20 microns (following
2291 ! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM)
2292 !--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
2293 ! of 300, 200, and 100 cm**-3, respectively
2294 !
2295 QAUT0=PI*RHOL*NCW*(20.E-6)**3/6.
2296 !
2297 !--- For calculating snow optical depths by considering bulk density of
2298 ! snow based on emails from Q. Fu (6/27-28/01), where optical
2299 ! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
2300 ! is effective radius, and DENS is the bulk density of snow.
2301 !
2302 ! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
2303 ! T = 1.5*1.E3*SWPrad/(Reff*DENS)
2304 !
2305 ! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
2306 !
2307 ! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
2308 !
2309 DO I=MDImin,MDImax
2310 SDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I)
2311 ENDDO
2312 !
2313 Thour_print=-DTPH/3600.
2314
2315 ! SH 0211/2002
2316 ! IF (PRINT_diag) THEN
2317
2318 !-------- Total and maximum quantities
2319 !
2320 ! NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts
2321 ! QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass
2322 ! QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass
2323 ! PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h)
2324 ! PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface
2325 ! ENDIF
2326
2327 !wrf
2328 IF(.NOT.RESTART)THEN
2329 MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
2330 MP_RESTART_STATE(MY_T2+1)=C1XPVS0
2331 MP_RESTART_STATE(MY_T2+2)=C2XPVS0
2332 MP_RESTART_STATE(MY_T2+3)=C1XPVS
2333 MP_RESTART_STATE(MY_T2+4)=C2XPVS
2334 MP_RESTART_STATE(MY_T2+5)=CIACW
2335 MP_RESTART_STATE(MY_T2+6)=CIACR
2336 MP_RESTART_STATE(MY_T2+7)=CRACW
2337 MP_RESTART_STATE(MY_T2+8)=CRAUT
2338 TBPVS_STATE(1:NX) =TBPVS(1:NX)
2339 TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
2340 ENDIF
2341
2342 ENDIF ! Allowed_to_read
2343
2344 RETURN
2345 !
2346 !-----------------------------------------------------------------------
2347 !
2348 9061 CONTINUE
2349 WRITE( errmess , '(A,I4)' ) &
2350 'module_mp_etanew: error opening ETAMPNEW_DATA on unit ' &
2351 &, etampnew_unit1
2352 CALL wrf_error_fatal(errmess)
2353 !
2354 !-----------------------------------------------------------------------
2355 END SUBROUTINE etanewinit
2356 !
2357 SUBROUTINE MY_GROWTH_RATES (DTPH)
2358 ! SUBROUTINE MY_GROWTH_RATES (DTPH,MY_GROWTH)
2359 !
2360 !--- Below are tabulated values for the predicted mass of ice crystals
2361 ! after 600 s of growth in water saturated conditions, based on
2362 ! calculations from Miller and Young (JAS, 1979). These values are
2363 ! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
2364 ! Young (1993). Values at temperatures colder than -27C were
2365 ! assumed to be invariant with temperature.
2366 !
2367 !--- Used to normalize Miller & Young (1979) calculations of ice growth
2368 ! over large time steps using their tabulated values at 600 s.
2369 ! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
2370 !
2371 IMPLICIT NONE
2372 !
2373 REAL,INTENT(IN) :: DTPH
2374 !
2375 REAL DT_ICE
2376 REAL,DIMENSION(35) :: MY_600
2377 !WRF
2378 !
2379 !-----------------------------------------------------------------------
2380 DATA MY_600 / &
2381 & 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
2382 & 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
2383 & 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
2384 & 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
2385 & 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
2386 & 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
2387 & 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
2388 !
2389 !-----------------------------------------------------------------------
2390 !
2391 DT_ICE=(DTPH/600.)**1.5
2392 MY_GROWTH=DT_ICE*MY_600
2393 !
2394 !-----------------------------------------------------------------------
2395 !
2396 END SUBROUTINE MY_GROWTH_RATES
2397 !
2398 !-----------------------------------------------------------------------
2399 !--------- Old GFS saturation vapor pressure lookup tables -----------
2400 !-----------------------------------------------------------------------
2401 !
2402 SUBROUTINE GPVS
2403 ! ******************************************************************
2404 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2405 ! . . .
2406 ! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
2407 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2408 !
2409 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
2410 ! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
2411 ! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
2412 ! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
2413 ! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
2414 !
2415 ! PROGRAM HISTORY LOG:
2416 ! 91-05-07 IREDELL
2417 ! 94-12-30 IREDELL EXPAND TABLE
2418 ! 96-02-19 HONG ICE EFFECT
2419 ! 01-11-29 JIN MODIFIED FOR WRF
2420 !
2421 ! USAGE: CALL GPVS
2422 !
2423 ! SUBPROGRAMS CALLED:
2424 ! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
2425 !
2426 ! COMMON BLOCKS:
2427 ! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
2428 !
2429 ! ATTRIBUTES:
2430 ! LANGUAGE: FORTRAN 90
2431 !
2432 !$$$
2433 IMPLICIT NONE
2434 real :: X,XINC,T
2435 integer :: JX
2436 !----------------------------------------------------------------------
2437 XINC=(XMAX-XMIN)/(NX-1)
2438 C1XPVS=1.-XMIN/XINC
2439 C2XPVS=1./XINC
2440 C1XPVS0=1.-XMIN/XINC
2441 C2XPVS0=1./XINC
2442 !
2443 DO JX=1,NX
2444 X=XMIN+(JX-1)*XINC
2445 T=X
2446 TBPVS(JX)=FPVSX(T)
2447 TBPVS0(JX)=FPVSX0(T)
2448 ENDDO
2449 !
2450 END SUBROUTINE GPVS
2451 !-----------------------------------------------------------------------
2452 !***********************************************************************
2453 !-----------------------------------------------------------------------
2454 REAL FUNCTION FPVS(T)
2455 !-----------------------------------------------------------------------
2456 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2457 ! . . .
2458 ! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
2459 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2460 !
2461 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
2462 ! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
2463 ! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
2464 ! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
2465 ! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
2466 ! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
2467 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2468 !
2469 ! PROGRAM HISTORY LOG:
2470 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2471 ! 94-12-30 IREDELL EXPAND TABLE
2472 ! 96-02-19 HONG ICE EFFECT
2473 ! 01-11-29 JIN MODIFIED FOR WRF
2474 !
2475 ! USAGE: PVS=FPVS(T)
2476 !
2477 ! INPUT ARGUMENT LIST:
2478 ! T - REAL TEMPERATURE IN KELVIN
2479 !
2480 ! OUTPUT ARGUMENT LIST:
2481 ! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2482 !
2483 ! ATTRIBUTES:
2484 ! LANGUAGE: FORTRAN 90
2485 !$$$
2486 IMPLICIT NONE
2487 real,INTENT(IN) :: T
2488 real XJ
2489 integer :: JX
2490 !-----------------------------------------------------------------------
2491 XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
2492 JX=MIN(XJ,NX-1.)
2493 FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
2494 !
2495 END FUNCTION FPVS
2496 !-----------------------------------------------------------------------
2497 !-----------------------------------------------------------------------
2498 REAL FUNCTION FPVS0(T)
2499 !-----------------------------------------------------------------------
2500 IMPLICIT NONE
2501 real,INTENT(IN) :: T
2502 real :: XJ1
2503 integer :: JX1
2504 !-----------------------------------------------------------------------
2505 XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
2506 JX1=MIN(XJ1,NX-1.)
2507 FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
2508 !
2509 END FUNCTION FPVS0
2510 !-----------------------------------------------------------------------
2511 !***********************************************************************
2512 !-----------------------------------------------------------------------
2513 REAL FUNCTION FPVSX(T)
2514 !-----------------------------------------------------------------------
2515 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2516 ! . . .
2517 ! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
2518 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2519 !
2520 ! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
2521 ! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
2522 ! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
2523 ! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
2524 ! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
2525 ! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
2526 ! TO GET THE FORMULA
2527 ! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2528 ! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
2529 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2530 !
2531 ! PROGRAM HISTORY LOG:
2532 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2533 ! 94-12-30 IREDELL EXACT COMPUTATION
2534 ! 96-02-19 HONG ICE EFFECT
2535 ! 01-11-29 JIN MODIFIED FOR WRF
2536 !
2537 ! USAGE: PVS=FPVSX(T)
2538 ! REFERENCE: EMANUEL(1994),116-117
2539 !
2540 ! INPUT ARGUMENT LIST:
2541 ! T - REAL TEMPERATURE IN KELVIN
2542 !
2543 ! OUTPUT ARGUMENT LIST:
2544 ! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2545 !
2546 ! ATTRIBUTES:
2547 ! LANGUAGE: FORTRAN 90
2548 !$$$
2549 IMPLICIT NONE
2550 !-----------------------------------------------------------------------
2551 real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2552 , CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
2553 , CICE=2.1060E+3,HSUB=2.8340E+6
2554 !
2555 real, parameter :: PSATK=PSAT*1.E-3
2556 real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2557 real, parameter :: DLDTI=CVAP-CICE &
2558 , XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
2559 real T,TR
2560 !-----------------------------------------------------------------------
2561 TR=TTP/T
2562 !
2563 IF(T.GE.TTP)THEN
2564 FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2565 ELSE
2566 FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
2567 ENDIF
2568 !
2569 END FUNCTION FPVSX
2570 !-----------------------------------------------------------------------
2571 !-----------------------------------------------------------------------
2572 REAL FUNCTION FPVSX0(T)
2573 !-----------------------------------------------------------------------
2574 IMPLICIT NONE
2575 real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2576 , CLIQ=4.1855E+3,CVAP=1.8460E+3 &
2577 , CICE=2.1060E+3,HSUB=2.8340E+6
2578 real,PARAMETER :: PSATK=PSAT*1.E-3
2579 real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2580 real,PARAMETER :: DLDTI=CVAP-CICE &
2581 , XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
2582 real :: T,TR
2583 !-----------------------------------------------------------------------
2584 TR=TTP/T
2585 FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2586 !
2587 END FUNCTION FPVSX0
2588 !
2589 END MODULE module_mp_etanew