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=-40. &
69 & ,T_ICEK=T0C+T_ICE &
70 & ,T_ICE_init=-15. &
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 IF (TC < 0) THEN
1097 XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
1098 IF (XLI .LE. MASSI(MDImin) ) THEN
1099 INDEXS=MDImin
1100 ELSE IF (XLI .LE. MASSI(450) ) THEN
1101 DLI=9.5885E5*XLI**.42066 ! DLI in microns
1102 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1103 ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
1104 DLI=3.9751E6*XLI**.49870 ! DLI in microns
1105 INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
1106 ELSE
1107 INDEXS=MDImax
1108 ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
1109 ENDIF ! End IF (TC < 0)
1110 !
1111 !--- Reduce excessive accumulation of ice at upper levels
1112 ! associated with strong grid-resolved ascent
1113 !
1114 !--- Force NLICE to be between NLImin and NLImax
1115 !
1116 !
1117 !--- 8/22/01: Increase density of large ice if maximum limits
1118 ! are reached for number concentration (NLImax) and mean size
1119 ! (MDImax). Done to increase fall out of ice.
1120 !
1121 DUM=MAX(NLImin, MIN(NLImax, NLICE) )
1122 IF (DUM.GE.NLImax .AND. INDEXS.GE.MDImax) &
1123 & RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
1124 ! WRITE(6,"(4(a12,g11.4,1x))")
1125 ! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
1126 ! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
1127 ! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
1128 ENDIF ! End IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) ...
1129 ENDDO ! End DO IPASS=0,1
1130 ENDIF ! End IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT)
1131 ENDIF ! End IF (ICE_logical)
1132 !
1133 !----------------------------------------------------------------------
1134 !--------------- Calculate individual processes -----------------------
1135 !----------------------------------------------------------------------
1136 !
1137 !--- Cloud water autoconversion to rain and collection by rain
1138 !
1139 IF (QW.GT.EPSQ .AND. TC.GE.T_ICE) THEN
1140 !
1141 !--- QW0 could be modified based on land/sea properties,
1142 ! presence of convection, etc. This is why QAUT0 and CRAUT
1143 ! are passed into the subroutine as externally determined
1144 ! parameters. Can be changed in the future if desired.
1145 !
1146 QW0=QAUT0*RRHO
1147 PRAUT=MAX(0., QW-QW0)*CRAUT
1148 IF (QLICE .GT. EPSQ) THEN
1149 !
1150 !--- Collection of cloud water by large ice particles ("snow")
1151 ! PIACWI=PIACW for riming, PIACWI=0 for shedding
1152 !
1153 FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1)
1154 PIACW=FWS*QW
1155 IF (TC .LT. 0.) PIACWI=PIACW ! Large ice riming
1156 ENDIF ! End IF (QLICE .GT. EPSQ)
1157 ENDIF ! End IF (QW.GT.EPSQ .AND. TC.GE.T_ICE)
1158 !
1159 !----------------------------------------------------------------------
1160 !--- Loop around some of the ice-phase processes if no ice should be present
1161 !----------------------------------------------------------------------
1162 !
1163 IF (ICE_logical .EQV. .FALSE.) GO TO 20
1164 !
1165 !--- Now the pretzel logic of calculating ice deposition
1166 !
1167 IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ)) THEN
1168 !
1169 !--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
1170 ! Sources of ice due to nucleation and convective detrainment are
1171 ! either poorly understood, poorly resolved at typical NWP
1172 ! resolutions, or are not represented (e.g., no detrained
1173 ! condensate in BMJ Cu scheme).
1174 !
1175 PCOND=-QW
1176 DUM1=TK+XLV1*PCOND ! Updated (dummy) temperature (deg K)
1177 DUM2=WV+QW ! Updated (dummy) water vapor mixing ratio
1178 DUM=1000.*FPVS(DUM1) ! Updated (dummy) saturation vapor pressure w/r/t ice
1179 DUM=RHgrd*EPS*DUM/(PP-DUM) ! Updated (dummy) saturation mixing ratio w/r/t ice
1180 IF (DUM2 .GT. DUM) PIDEP=DEPOSIT (PP, DUM1, DUM2)
1181 DWVi=0. ! Used only for debugging
1182 !
1183 ELSE IF (TC .LT. 0.) THEN
1184 !
1185 !--- These quantities are handy for ice deposition/sublimation
1186 ! PIDEP_max - max deposition or minimum sublimation to ice saturation
1187 !
1188 DENOMI=1.+XLS2*QSI*TK2
1189 DWVi=MIN(WVQW,QSW)-QSI
1190 PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
1191 IF (QTICE .GT. 0.) THEN
1192 !
1193 !--- Calculate ice deposition/sublimation
1194 ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1195 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1196 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1197 ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ;
1198 ! VENTI2 - m**2/s**.5 ; DIDEP - unitless
1199 !
1200 SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1201 ABI=1./(RHO*XLS3*QSI*TK2/THERM_COND+1./DIFFUS)
1202 !
1203 !--- VENTIL - Number concentration * ventilation factors for large ice
1204 !--- VENTIS - Number concentration * ventilation factors for small ice
1205 !
1206 !--- Variation in the number concentration of ice with time is not
1207 ! accounted for in these calculations (could be in the future).
1208 !
1209 VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE
1210 VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE
1211 DIDEP=ABI*(VENTIL+VENTIS)*DTPH
1212 !
1213 !--- Account for change in water vapor supply w/ time
1214 !
1215 IF (DIDEP .GE. Xratio)then
1216 DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
1217 endif
1218 IF (DWVi .GT. 0.) THEN
1219 PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
1220 ELSE IF (DWVi .LT. 0.) THEN
1221 PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
1222 ENDIF
1223 !
1224 ELSE IF (WVQW.GT.QSI .AND. TC.LE.T_ICE_init) THEN
1225 !
1226 !--- Ice nucleation in near water-saturated conditions. Ice crystal
1227 ! growth during time step calculated using Miller & Young (1979, JAS).
1228 !--- These deposition rates could drive conditions below water saturation,
1229 ! which is the basis of these calculations. Intended to approximate
1230 ! more complex & computationally intensive calculations.
1231 !
1232 INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
1233 !
1234 !--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
1235 !
1236 !--- DUM2 is the number of ice crystals nucleated at water-saturated
1237 ! conditions based on Meyers et al. (JAM, 1992).
1238 !
1239 !--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
1240 ! if DUM2 values are increased in future experiments
1241 !
1242 DUM1=QSW/QSI-1.
1243 DUM2=1.E3*EXP(12.96*DUM1-.639)
1244 PIDEP=MIN(PIDEP_max, DUM2*MY_GROWTH(INDEX_MY)*RRHO)
1245 !
1246 ENDIF ! End IF (QTICE .GT. 0.)
1247 !
1248 ENDIF ! End IF (TC.LT.T_ICE .AND. (WV.GT.QSIgrd .OR. QW.GT.EPSQ))
1249 !
1250 !------------------------------------------------------------------------
1251 !
1252 20 CONTINUE ! Jump here if conditions for ice are not present
1253
1254
1255 !
1256 !------------------------------------------------------------------------
1257 !
1258 !--- Cloud water condensation
1259 !
1260 IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd)) THEN
1261 IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.) THEN
1262 PCOND=CONDENSE (PP, QW, TK, WV)
1263 ELSE
1264 !
1265 !--- Modify cloud condensation in response to ice processes
1266 !
1267 DUM=XLV*QSWgrd*RCPRV*TK2
1268 DENOMWI=1.+XLS*DUM
1269 DENOMF=XLF*DUM
1270 DUM=MAX(0., PIDEP)
1271 PCOND=(WV-QSWgrd-DENOMWI*DUM-DENOMF*PIACWI)/DENOMW
1272 DUM1=-QW
1273 DUM2=PCOND-PIACW
1274 IF (DUM2 .LT. DUM1) THEN
1275 !
1276 !--- Limit cloud water sinks
1277 !
1278 DUM=DUM1/DUM2
1279 PCOND=DUM*PCOND
1280 PIACW=DUM*PIACW
1281 PIACWI=DUM*PIACWI
1282 ENDIF ! End IF (DUM2 .LT. DUM1)
1283 ENDIF ! End IF (PIACWI.EQ.0. .AND. PIDEP.EQ.0.)
1284 ENDIF ! End IF (TC.GE.T_ICE .AND. (QW.GT.EPSQ .OR. WV.GT.QSWgrd))
1285 !
1286 !--- Limit freezing of accreted rime to prevent temperature oscillations,
1287 ! a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
1288 !
1289 TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI
1290 IF (TCC .GT. 0.) THEN
1291 PIACWI=0.
1292 TCC=TC+XLV1*PCOND+XLS1*PIDEP
1293 ENDIF
1294 IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical) THEN
1295 !
1296 !--- Calculate melting and evaporation/condensation
1297 ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ;
1298 ! VENTIL - m**-2 ; VENTI1 - m ;
1299 ! VENTI2 - m**2/s**.5 ; CIEVP - /s
1300 !
1301 SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1302 VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))
1303 AIEVP=VENTIL*DIFFUS*DTPH
1304 IF (AIEVP .LT. Xratio) THEN
1305 DIEVP=AIEVP
1306 ELSE
1307 DIEVP=1.-EXP(-AIEVP)
1308 ENDIF
1309 QSW0=EPS*ESW0/(PP-ESW0)
1310 DWV0=MIN(WV,QSW)-QSW0
1311 DUM=QW+PCOND
1312 IF (WV.LT.QSW .AND. DUM.LE.EPSQ) THEN
1313 !
1314 !--- Evaporation from melting snow (sink of snow) or shedding
1315 ! of water condensed onto melting snow (source of rain)
1316 !
1317 DUM=DWV0*DIEVP
1318 PIEVP=MAX( MIN(0., DUM), PILOSS)
1319 PICND=MAX(0., DUM)
1320 ENDIF ! End IF (WV.LT.QSW .AND. DUM.LE.EPSQ)
1321 PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF
1322 !
1323 !--- Limit melting to prevent temperature oscillations across 0C
1324 !
1325 DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
1326 PIMLT=MIN(PIMLT, DUM1)
1327 !
1328 !--- Limit loss of snow by melting (>0) and evaporation
1329 !
1330 DUM=PIEVP-PIMLT
1331 IF (DUM .LT. PILOSS) THEN
1332 DUM1=PILOSS/DUM
1333 PIMLT=PIMLT*DUM1
1334 PIEVP=PIEVP*DUM1
1335 ENDIF ! End IF (DUM .GT. QTICE)
1336 ENDIF ! End IF (TC.GT.0. .AND. TCC.GT.0. .AND. ICE_logical)
1337 !
1338 !--- IMPORTANT: Estimate time-averaged properties.
1339 !
1340 ! * TOT_RAIN - total mass of rain before microphysics, which is the sum of
1341 ! the total mass of rain in the current layer and the input
1342 ! flux of rain from above
1343 ! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction)
1344 ! * QTRAIN - time-averaged mixing ratio of rain (kg/kg)
1345 ! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from
1346 ! above and the rain within the layer
1347 ! * RQR - rain content (kg/m**3)
1348 ! * INDEXR - mean size of rain drops to the nearest 1 micron in size
1349 ! * N0r - intercept of rain size distribution (typically 10**6 m**-4)
1350 !
1351 TOT_RAIN=0.
1352 VRAIN1=0.
1353 QTRAIN=0.
1354 PRLOSS=0.
1355 RQR=0.
1356 N0r=0.
1357 INDEXR=MDRmin
1358 INDEXR1=INDEXR !-- For debugging only
1359 IF (RAIN_logical) THEN
1360 IF (ARAIN .LE. 0.) THEN
1361 INDEXR=MDRmin
1362 VRAIN1=0.
1363 ELSE
1364 !
1365 !--- INDEXR (related to mean diameter) & N0r could be modified
1366 ! by land/sea properties, presence of convection, etc.
1367 !
1368 !--- Rain rate normalized to a density of 1.194 kg/m**3
1369 !
1370 RR=ARAIN/(DTPH*GAMMAR)
1371 !
1372 IF (RR .LE. RR_DRmin) THEN
1373 !
1374 !--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
1375 ! instead vary N0r with rain rate
1376 !
1377 INDEXR=MDRmin
1378 ELSE IF (RR .LE. RR_DR1) THEN
1379 !
1380 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1381 ! for mean diameters (Dr) between 0.05 and 0.10 mm:
1382 ! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
1383 ! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
1384 ! RR in kg/(m**2*s)
1385 ! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
1386 !
1387 INDEXR=INT( 1.123E3*RR**.1947 + .5 )
1388 INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
1389 ELSE IF (RR .LE. RR_DR2) THEN
1390 !
1391 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1392 ! for mean diameters (Dr) between 0.10 and 0.20 mm:
1393 ! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
1394 ! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
1395 ! RR in kg/(m**2*s)
1396 ! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
1397 !
1398 INDEXR=INT( 1.225E3*RR**.2017 + .5 )
1399 INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
1400 ELSE IF (RR .LE. RR_DR3) THEN
1401 !
1402 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1403 ! for mean diameters (Dr) between 0.20 and 0.32 mm:
1404 ! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
1405 ! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
1406 ! RR in kg/(m**2*s)
1407 ! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
1408 !
1409 INDEXR=INT( 1.3006E3*RR**.2083 + .5 )
1410 INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
1411 ELSE IF (RR .LE. RR_DRmax) THEN
1412 !
1413 !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
1414 ! for mean diameters (Dr) between 0.32 and 0.45 mm:
1415 ! V(Dr)=944.8*Dr**.6636, V in m/s and Dr in m
1416 ! RR = PI*1000.*N0r0*944.8*Dr**(4+.6636) = 2.3745e13*Dr**4.6636,
1417 ! RR in kg/(m**2*s)
1418 ! Dr (m) = 1.355e-3*RR**.2144 -> Dr (microns) = 1.355e3*RR**.2144
1419 !
1420 INDEXR=INT( 1.355E3*RR**.2144 + .5 )
1421 INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
1422 ELSE
1423 !
1424 !--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
1425 ! instead vary N0r with rain rate
1426 !
1427 INDEXR=MDRmax
1428 ENDIF ! End IF (RR .LE. RR_DRmin) etc.
1429 VRAIN1=GAMMAR*VRAIN(INDEXR)
1430 ENDIF ! End IF (ARAIN .LE. 0.)
1431 INDEXR1=INDEXR ! For debugging only
1432 TOT_RAIN=THICK*QR+BLEND*ARAIN
1433 QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
1434 PRLOSS=-TOT_RAIN/THICK
1435 RQR=RHO*QTRAIN
1436 !
1437 !--- RQR - time-averaged rain content (kg/m**3)
1438 !
1439 IF (RQR .LE. RQR_DRmin) THEN
1440 N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
1441 INDEXR=MDRmin
1442 ELSE IF (RQR .GE. RQR_DRmax) THEN
1443 N0r=CN0r_DMRmax*RQR
1444 INDEXR=MDRmax
1445 ELSE
1446 N0r=N0r0
1447 INDEXR=MAX( XMRmin, MIN(CN0r0*RQR**.25, XMRmax) )
1448 ENDIF
1449 !
1450 IF (TC .LT. T_ICE) THEN
1451 PIACR=-PRLOSS
1452 ELSE
1453 DWVr=WV-PCOND-QSW
1454 DUM=QW+PCOND
1455 IF (DWVr.LT.0. .AND. DUM.LE.EPSQ) THEN
1456 !
1457 !--- Rain evaporation
1458 !
1459 ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
1460 ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
1461 !
1462 ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ;
1463 ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ;
1464 ! CREVP - unitless
1465 !
1466 RFACTOR=GAMMAR**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
1467 ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
1468 !
1469 !--- Note that VENTR1, VENTR2 lookup tables do not include the
1470 ! 1/Davg multiplier as in the ice tables
1471 !
1472 VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR))
1473 CREVP=ABW*VENTR*DTPH
1474 IF (CREVP .LT. Xratio) THEN
1475 DUM=DWVr*CREVP
1476 ELSE
1477 DUM=DWVr*(1.-EXP(-CREVP*DENOMW))/DENOMW
1478 ENDIF
1479 PREVP=MAX(DUM, PRLOSS)
1480 ELSE IF (QW .GT. EPSQ) THEN
1481 FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR)
1482 PRACW=MIN(1.,FWR)*QW
1483 ENDIF ! End IF (DWVr.LT.0. .AND. DUM.LE.EPSQ)
1484 !
1485 IF (TC.LT.0. .AND. TCC.LT.0.) THEN
1486 !
1487 !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
1488 ! - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
1489 !
1490 DUM=.001*FLOAT(INDEXR)
1491 DUM=(EXP(ABFR*TC)-1.)*DUM*DUM*DUM*DUM*DUM*DUM*DUM
1492 PIACR=MIN(CBFR*N0r*RRHO*DUM, QTRAIN)
1493 IF (QLICE .GT. EPSQ) THEN
1494 !
1495 !--- Freezing of rain by collisions w/ large ice
1496 !
1497 DUM=GAMMAR*VRAIN(INDEXR)
1498 DUM1=DUM-VSNOW
1499 !
1500 !--- DUM2 - Difference in spectral fall speeds of rain and
1501 ! large ice, parameterized following eq. (48) on p. 112 of
1502 ! Murakami (J. Meteor. Soc. Japan, 1990)
1503 !
1504 DUM2=(DUM1*DUM1+.04*DUM*VSNOW)**.5
1505 DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS &
1506 & +.5E-12*INDEXS*INDEXS
1507 FIR=MIN(1., CIACR*NLICE*DUM1*DUM2)
1508 !
1509 !--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
1510 !
1511 PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
1512 ENDIF ! End IF (QLICE .GT. EPSQ)
1513 DUM=PREVP-PIACR
1514 If (DUM .LT. PRLOSS) THEN
1515 DUM1=PRLOSS/DUM
1516 PREVP=DUM1*PREVP
1517 PIACR=DUM1*PIACR
1518 ENDIF ! End If (DUM .LT. PRLOSS)
1519 ENDIF ! End IF (TC.LT.0. .AND. TCC.LT.0.)
1520 ENDIF ! End IF (TC .LT. T_ICE)
1521 ENDIF ! End IF (RAIN_logical)
1522 !
1523 !----------------------------------------------------------------------
1524 !---------------------- Main Budget Equations -------------------------
1525 !----------------------------------------------------------------------
1526 !
1527 !
1528 !-----------------------------------------------------------------------
1529 !--- Update fields, determine characteristics for next lower layer ----
1530 !-----------------------------------------------------------------------
1531 !
1532 !--- Carefully limit sinks of cloud water
1533 !
1534 DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
1535 IF (DUM1 .GT. QW) THEN
1536 DUM=QW/DUM1
1537 PIACW=DUM*PIACW
1538 PIACWI=DUM*PIACWI
1539 PRAUT=DUM*PRAUT
1540 PRACW=DUM*PRACW
1541 IF (PCOND .LT. 0.) PCOND=DUM*PCOND
1542 ENDIF
1543 PIACWR=PIACW-PIACWI ! TC >= 0C
1544 !
1545 !--- QWnew - updated cloud water mixing ratio
1546 !
1547 DELW=PCOND-PIACW-PRAUT-PRACW
1548 QWnew=QW+DELW
1549 IF (QWnew .LE. EPSQ) QWnew=0.
1550 IF (QW.GT.0. .AND. QWnew.NE.0.) THEN
1551 DUM=QWnew/QW
1552 IF (DUM .LT. TOLER) QWnew=0.
1553 ENDIF
1554 !
1555 !--- Update temperature and water vapor mixing ratios
1556 !
1557 DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
1558 & +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT)
1559 Tnew=TK+DELT
1560 !
1561 DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
1562 WVnew=WV+DELV
1563 !
1564 !--- Update ice mixing ratios
1565 !
1566 !---
1567 ! * TOT_ICEnew - total mass (small & large) ice after microphysics,
1568 ! which is the sum of the total mass of large ice in the
1569 ! current layer and the flux of ice out of the grid box below
1570 ! * RimeF - Rime Factor, which is the mass ratio of total (unrimed &
1571 ! rimed) ice mass to the unrimed ice mass (>=1)
1572 ! * QInew - updated mixing ratio of total (large & small) ice in layer
1573 ! -> TOT_ICEnew=QInew*THICK+BLDTRH*QLICEnew*VSNOW
1574 ! -> But QLICEnew=QInew*FLIMASS, so
1575 ! -> TOT_ICEnew=QInew*(THICK+BLDTRH*FLIMASS*VSNOW)
1576 ! * ASNOWnew - updated accumulation of snow at bottom of grid cell
1577 !---
1578 !
1579 DELI=0.
1580 RimeF=1.
1581 IF (ICE_logical) THEN
1582 DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
1583 TOT_ICEnew=TOT_ICE+THICK*DELI
1584 IF (TOT_ICE.GT.0. .AND. TOT_ICEnew.NE.0.) THEN
1585 DUM=TOT_ICEnew/TOT_ICE
1586 IF (DUM .LT. TOLER) TOT_ICEnew=0.
1587 ENDIF
1588 IF (TOT_ICEnew .LE. CLIMIT) THEN
1589 TOT_ICEnew=0.
1590 RimeF=1.
1591 QInew=0.
1592 ASNOWnew=0.
1593 ELSE
1594 !
1595 !--- Update rime factor if appropriate
1596 !
1597 DUM=PIACWI+PIACR
1598 IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ) THEN
1599 RimeF=RimeF1
1600 ELSE
1601 !
1602 !--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
1603 ! DUM1 - Total ice mass, rimed & unrimed
1604 ! DUM2 - Estimated mass of *unrimed* ice
1605 !
1606 DUM1=TOT_ICE+THICK*(PIDEP+DUM)
1607 DUM2=TOT_ICE/RimeF1+THICK*PIDEP
1608 IF (DUM2 .LE. 0.) THEN
1609 RimeF=RFmax
1610 ELSE
1611 RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
1612 ENDIF
1613 ENDIF ! End IF (DUM.LE.EPSQ .AND. PIDEP.LE.EPSQ)
1614 QInew=TOT_ICEnew/(THICK+BLDTRH*FLIMASS*VSNOW)
1615 IF (QInew .LE. EPSQ) QInew=0.
1616 IF (QI.GT.0. .AND. QInew.NE.0.) THEN
1617 DUM=QInew/QI
1618 IF (DUM .LT. TOLER) QInew=0.
1619 ENDIF
1620 ASNOWnew=BLDTRH*FLIMASS*VSNOW*QInew
1621 IF (ASNOW.GT.0. .AND. ASNOWnew.NE.0.) THEN
1622 DUM=ASNOWnew/ASNOW
1623 IF (DUM .LT. TOLER) ASNOWnew=0.
1624 ENDIF
1625 ENDIF ! End IF (TOT_ICEnew .LE. CLIMIT)
1626 ENDIF ! End IF (ICE_logical)
1627
1628
1629 !
1630 !--- Update rain mixing ratios
1631 !
1632 !---
1633 ! * TOT_RAINnew - total mass of rain after microphysics
1634 ! current layer and the input flux of ice from above
1635 ! * VRAIN2 - time-averaged fall speed of rain in grid and below
1636 ! (with air resistance correction)
1637 ! * QRnew - updated rain mixing ratio in layer
1638 ! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
1639 ! * ARAINnew - updated accumulation of rain at bottom of grid cell
1640 !---
1641 !
1642 DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
1643 TOT_RAINnew=TOT_RAIN+THICK*DELR
1644 IF (TOT_RAIN.GT.0. .AND. TOT_RAINnew.NE.0.) THEN
1645 DUM=TOT_RAINnew/TOT_RAIN
1646 IF (DUM .LT. TOLER) TOT_RAINnew=0.
1647 ENDIF
1648 IF (TOT_RAINnew .LE. CLIMIT) THEN
1649 TOT_RAINnew=0.
1650 VRAIN2=0.
1651 QRnew=0.
1652 ARAINnew=0.
1653 ELSE
1654 !
1655 !--- 1st guess time-averaged rain rate at bottom of grid box
1656 !
1657 RR=TOT_RAINnew/(DTPH*GAMMAR)
1658 !
1659 !--- Use same algorithm as above for calculating mean drop diameter
1660 ! (IDR, in microns), which is used to estimate the time-averaged
1661 ! fall speed of rain drops at the bottom of the grid layer. This
1662 ! isn't perfect, but the alternative is solving a transcendental
1663 ! equation that is numerically inefficient and nasty to program
1664 ! (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
1665 !
1666 IF (RR .LE. RR_DRmin) THEN
1667 IDR=MDRmin
1668 ELSE IF (RR .LE. RR_DR1) THEN
1669 IDR=INT( 1.123E3*RR**.1947 + .5 )
1670 IDR=MAX( MDRmin, MIN(IDR, MDR1) )
1671 ELSE IF (RR .LE. RR_DR2) THEN
1672 IDR=INT( 1.225E3*RR**.2017 + .5 )
1673 IDR=MAX( MDR1, MIN(IDR, MDR2) )
1674 ELSE IF (RR .LE. RR_DR3) THEN
1675 IDR=INT( 1.3006E3*RR**.2083 + .5 )
1676 IDR=MAX( MDR2, MIN(IDR, MDR3) )
1677 ELSE IF (RR .LE. RR_DRmax) THEN
1678 IDR=INT( 1.355E3*RR**.2144 + .5 )
1679 IDR=MAX( MDR3, MIN(IDR, MDRmax) )
1680 ELSE
1681 IDR=MDRmax
1682 ENDIF ! End IF (RR .LE. RR_DRmin)
1683 VRAIN2=GAMMAR*VRAIN(IDR)
1684 QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
1685 IF (QRnew .LE. EPSQ) QRnew=0.
1686 IF (QR.GT.0. .AND. QRnew.NE.0.) THEN
1687 DUM=QRnew/QR
1688 IF (DUM .LT. TOLER) QRnew=0.
1689 ENDIF
1690 ARAINnew=BLDTRH*VRAIN2*QRnew
1691 IF (ARAIN.GT.0. .AND. ARAINnew.NE.0.) THEN
1692 DUM=ARAINnew/ARAIN
1693 IF (DUM .LT. TOLER) ARAINnew=0.
1694 ENDIF
1695 ENDIF
1696 !
1697 WCnew=QWnew+QRnew+QInew
1698 !
1699 !----------------------------------------------------------------------
1700 !-------------- Begin debugging & verification ------------------------
1701 !----------------------------------------------------------------------
1702 !
1703 !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
1704 !
1705
1706
1707 QT=THICK*(WV+WC)+ARAIN+ASNOW
1708 QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
1709 BUDGET=QT-QTnew
1710 !
1711 !--- Additional check on budget preservation, accounting for truncation effects
1712 !
1713 DBG_logical=.FALSE.
1714 ! DUM=ABS(BUDGET)
1715 ! IF (DUM .GT. TOLER) THEN
1716 ! DUM=DUM/MIN(QT, QTnew)
1717 ! IF (DUM .GT. TOLER) DBG_logical=.TRUE.
1718 ! ENDIF
1719 !!
1720 ! DUM=(RHgrd+.001)*QSInew
1721 ! IF ( (QWnew.GT.EPSQ) .OR. QRnew.GT.EPSQ .OR. WVnew.GT.DUM)
1722 ! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE.
1723 !
1724 ! IF (TC.GT.5. .AND. QInew.GT.EPSQ) DBG_logical=.TRUE.
1725 !
1726 IF ((WVnew.LT.EPSQ .OR. DBG_logical) .AND. PRINT_diag) THEN
1727 !
1728 WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index,&
1729 & ' L=',L,' LSFC=',LSFC
1730 !
1731 ESW=1000.*FPVS0(Tnew)
1732 QSWnew=EPS*ESW/(PP-ESW)
1733 IF (TC.LT.0. .OR. Tnew .LT. 0.) THEN
1734 ESI=1000.*FPVS(Tnew)
1735 QSInew=EPS*ESI/(PP-ESI)
1736 ELSE
1737 QSI=QSW
1738 QSInew=QSWnew
1739 ENDIF
1740 WSnew=QSInew
1741 WRITE(6,"(4(a12,g11.4,1x))") &
1742 & '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
1743 & '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
1744 & 'RHgrd=',RHgrd, &
1745 & '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
1746 & 'RHInew=',WVnew/QSInew, &
1747 & '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
1748 & '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
1749 & '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
1750 & '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
1751 & '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
1752 & 'ASNOWnew=',ASNOWnew, &
1753 & '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
1754 & 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
1755 & '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
1756 !
1757 WRITE(6,"(4(a12,g11.4,1x))") &
1758 & '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
1759 & '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
1760 & '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
1761 & '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
1762 & PIMLT, &
1763 & '{} PIACR=',PIACR
1764 !
1765 IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
1766 & '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
1767 & 'VSNOW=',VSNOW, &
1768 & '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
1769 & 'FLIMASS=',FLIMASS, &
1770 & '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
1771 & 'QTICE=',QTICE, &
1772 & '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
1773 & 'EMAIRI=',EMAIRI, &
1774 & '{} RimeF=',RimeF
1775 !
1776 IF (TOT_RAIN.GT.0. .OR. TOT_RAINnew.GT.0.) &
1777 & WRITE(6,"(4(a12,g11.4,1x))") &
1778 & '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
1779 & 'GAMMAR=',GAMMAR,'N0r=',N0r, &
1780 & '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
1781 & '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
1782 & 'VOLR2=',THICK+BLDTRH*VRAIN2
1783 !
1784 IF (PRAUT .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
1785 !
1786 IF (PRACW .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
1787 !
1788 IF (PIACR .GT. 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
1789 !
1790 DUM=PIMLT+PICND-PREVP-PIEVP
1791 IF (DUM.GT.0. .or. DWVi.NE.0.) &
1792 & WRITE(6,"(4(a12,g11.4,1x))") &
1793 & '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
1794 & 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
1795 !
1796 IF (PREVP .LT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1797 & '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
1798 & '{} DWVr=',DWVr,'DENOMW=',DENOMW
1799 !
1800 IF (PIDEP.NE.0. .AND. DWVi.NE.0.) &
1801 & WRITE(6,"(4(a12,g11.4,1x))") &
1802 & '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
1803 & 'SFACTOR=',SFACTOR, &
1804 & '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1805 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1806 & '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
1807 !
1808 IF (PIDEP.GT.0. .AND. PCOND.NE.0.) &
1809 & WRITE(6,"(4(a12,g11.4,1x))") &
1810 & '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
1811 & 'DUM2=',PCOND-PIACW
1812 !
1813 IF (FWS .GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1814 & '{} FWS=',FWS
1815 !
1816 DUM=PIMLT+PICND-PIEVP
1817 IF (DUM.GT. 0.) WRITE(6,"(4(a12,g11.4,1x))") &
1818 & '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
1819 & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
1820 & '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
1821 !
1822 ENDIF
1823
1824
1825 !
1826 !-----------------------------------------------------------------------
1827 !--------------- Water budget statistics & maximum values --------------
1828 !-----------------------------------------------------------------------
1829 !
1830 IF (PRINT_diag) THEN
1831 ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
1832 IF (QInew .GT. EPSQ) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
1833 IF (QInew.GT.EPSQ .AND. QRnew+QWnew.GT.EPSQ) &
1834 & NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
1835 IF (QWnew .GT. EPSQ) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
1836 IF (QRnew .GT. EPSQ) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
1837 !
1838 QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
1839 QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
1840 QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
1841 QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
1842 QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
1843 QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
1844 QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
1845 QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
1846 !
1847 QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
1848 QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
1849 QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
1850 QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
1851 QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
1852 QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
1853 QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
1854 QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
1855 QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
1856 QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
1857 QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
1858 QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
1859 !
1860 QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
1861 QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
1862 QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
1863 QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
1864 QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
1865 QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
1866 IF (QInew .GT. 0.) &
1867 & QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
1868 !
1869 ENDIF
1870 !
1871 !----------------------------------------------------------------------
1872 !------------------------- Update arrays ------------------------------
1873 !----------------------------------------------------------------------
1874 !
1875
1876
1877 T_col(L)=Tnew ! Updated temperature
1878 !
1879 QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity
1880 WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio
1881 QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio
1882 QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio
1883 QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio
1884 RimeF_col(L)=RimeF ! Updated rime factor
1885 ASNOW=ASNOWnew ! Updated accumulated snow
1886 ARAIN=ARAINnew ! Updated accumulated rain
1887 !
1888 !#######################################################################
1889 !
1890 10 CONTINUE ! ##### End "L" loop through model levels #####
1891
1892
1893 !
1894 !#######################################################################
1895 !
1896 !-----------------------------------------------------------------------
1897 !--------------------------- Return to GSMDRIVE -----------------------
1898 !-----------------------------------------------------------------------
1899 !
1900 CONTAINS
1901 !#######################################################################
1902 !--------- Produces accurate calculation of cloud condensation ---------
1903 !#######################################################################
1904 !
1905 REAL FUNCTION CONDENSE (PP, QW, TK, WV)
1906 !
1907 !---------------------------------------------------------------------------------
1908 !------ The Asai (1965) algorithm takes into consideration the release of ------
1909 !------ latent heat in increasing the temperature & in increasing the ------
1910 !------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------
1911 !---------------------------------------------------------------------------------
1912 !
1913 IMPLICIT NONE
1914 !
1915 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1916 REAL (KIND=HIGH_PRES), PARAMETER :: &
1917 & RHLIMIT=.001, RHLIMIT1=-RHLIMIT
1918 REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum
1919 !
1920 REAL,INTENT(IN) :: QW,PP,WV,TK
1921 REAL WVdum,Tdum,XLV2,DWV,WS,ESW,XLV1,XLV
1922 integer nsteps
1923 !
1924 !-----------------------------------------------------------------------
1925 !
1926 !--- LV (T) is from Bolton (JAS, 1980)
1927 !
1928 XLV=3.148E6-2370.*TK
1929 XLV1=XLV*RCP
1930 XLV2=XLV*XLV*RCPRV
1931 Tdum=TK
1932 WVdum=WV
1933 WCdum=QW
1934 ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water
1935 WS=RHgrd*EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
1936 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1937 SSAT=DWV/WS ! Supersaturation ratio
1938 CONDENSE=0.
1939 nsteps = 0
1940 DO WHILE ((SSAT.LT.RHLIMIT1 .AND. WCdum.GT.EPSQ) &
1941 & .OR. SSAT.GT.RHLIMIT)
1942 nsteps = nsteps + 1
1943 COND=DWV/(1.+XLV2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
1944 COND=MAX(COND, -WCdum) ! Limit cloud water evaporation
1945 Tdum=Tdum+XLV1*COND ! Updated temperature
1946 WVdum=WVdum-COND ! Updated water vapor mixing ratio
1947 WCdum=WCdum+COND ! Updated cloud water mixing ratio
1948 CONDENSE=CONDENSE+COND ! Total cloud water condensation
1949 ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water
1950 WS=RHgrd*EPS*ESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water
1951 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1952 SSAT=DWV/WS ! Grid-scale supersaturation ratio
1953 ENDDO
1954 !
1955 END FUNCTION CONDENSE
1956 !
1957 !#######################################################################
1958 !---------------- Calculate ice deposition at T<T_ICE ------------------
1959 !#######################################################################
1960 !
1961 REAL FUNCTION DEPOSIT (PP, Tdum, WVdum)
1962 !
1963 !--- Also uses the Asai (1965) algorithm, but uses a different target
1964 ! vapor pressure for the adjustment
1965 !
1966 IMPLICIT NONE
1967 !
1968 INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
1969 REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, &
1970 & RHLIMIT1=-RHLIMIT
1971 REAL (KIND=HIGH_PRES) :: DEP, SSAT
1972 !
1973 real,INTENT(IN) :: PP
1974 real,INTENT(INOUT) :: WVdum,Tdum
1975 real ESI,WS,DWV
1976 !
1977 !-----------------------------------------------------------------------
1978 !
1979 ESI=1000.*FPVS(Tdum) ! Saturation vapor press w/r/t ice
1980 WS=RHgrd*EPS*ESI/(PP-ESI) ! Saturation mixing ratio
1981 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1982 SSAT=DWV/WS ! Supersaturation ratio
1983 DEPOSIT=0.
1984 DO WHILE (SSAT.GT.RHLIMIT .OR. SSAT.LT.RHLIMIT1)
1985 !
1986 !--- Note that XLVS2=LS*LV/(CP*RV)=LV*WS/(RV*T*T)*(LS/CP*DEP1),
1987 ! where WS is the saturation mixing ratio following Clausius-
1988 ! Clapeyron (see Asai,1965; Young,1993,p.405)
1989 !
1990 DEP=DWV/(1.+XLS2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan)
1991 Tdum=Tdum+XLS1*DEP ! Updated temperature
1992 WVdum=WVdum-DEP ! Updated ice mixing ratio
1993 DEPOSIT=DEPOSIT+DEP ! Total ice deposition
1994 ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice
1995 WS=RHgrd*EPS*ESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice
1996 DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio
1997 SSAT=DWV/WS ! Grid-scale supersaturation ratio
1998 ENDDO
1999 !
2000 END FUNCTION DEPOSIT
2001 !
2002 END SUBROUTINE EGCP01COLUMN
2003
2004 !#######################################################################
2005 !------- Initialize constants & lookup tables for microphysics ---------
2006 !#######################################################################
2007 !
2008
2009 ! SH 0211/2002
2010
2011 !-----------------------------------------------------------------------
2012 SUBROUTINE etanewinit (GSMDT,DT,DELX,DELY,LOWLYR,restart, &
2013 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
2014 & MP_RESTART_STATE,TBPVS_STATE,TBPVS0_STATE, &
2015 & ALLOWED_TO_READ, &
2016 & IDS,IDE,JDS,JDE,KDS,KDE, &
2017 & IMS,IME,JMS,JME,KMS,KME, &
2018 & ITS,ITE,JTS,JTE,KTS,KTE )
2019 !-----------------------------------------------------------------------
2020 !-------------------------------------------------------------------------------
2021 !--- SUBPROGRAM DOCUMENTATION BLOCK
2022 ! PRGRMMR: Ferrier ORG: W/NP22 DATE: February 2001
2023 ! Jin ORG: W/NP22 DATE: January 2002
2024 ! (Modification for WRF structure)
2025 !
2026 !-------------------------------------------------------------------------------
2027 ! ABSTRACT:
2028 ! * Reads various microphysical lookup tables used in COLUMN_MICRO
2029 ! * Lookup tables were created "offline" and are read in during execution
2030 ! * Creates lookup tables for saturation vapor pressure w/r/t water & ice
2031 !-------------------------------------------------------------------------------
2032 !
2033 ! USAGE: CALL etanewinit FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
2034 !
2035 ! INPUT ARGUMENT LIST:
2036 ! DTPH - physics time step (s)
2037 !
2038 ! OUTPUT ARGUMENT LIST:
2039 ! NONE
2040 !
2041 ! OUTPUT FILES:
2042 ! NONE
2043 !
2044 ! SUBROUTINES:
2045 ! MY_GROWTH_RATES - lookup table for growth of nucleated ice
2046 ! GPVS - lookup tables for saturation vapor pressure (water, ice)
2047 !
2048 ! UNIQUE: NONE
2049 !
2050 ! LIBRARY: NONE
2051 !
2052 ! COMMON BLOCKS:
2053 ! CMICRO_CONS - constants used in GSMCOLUMN
2054 ! CMY600 - lookup table for growth of ice crystals in
2055 ! water saturated conditions (Miller & Young, 1979)
2056 ! IVENT_TABLES - lookup tables for ventilation effects of ice
2057 ! IACCR_TABLES - lookup tables for accretion rates of ice
2058 ! IMASS_TABLES - lookup tables for mass content of ice
2059 ! IRATE_TABLES - lookup tables for precipitation rates of ice
2060 ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
2061 ! MAPOT - Need lat/lon grid resolution
2062 ! RVENT_TABLES - lookup tables for ventilation effects of rain
2063 ! RACCR_TABLES - lookup tables for accretion rates of rain
2064 ! RMASS_TABLES - lookup tables for mass content of rain
2065 ! RVELR_TABLES - lookup tables for fall speeds of rain
2066 ! RRATE_TABLES - lookup tables for precipitation rates of rain
2067 !
2068 ! ATTRIBUTES:
2069 ! LANGUAGE: FORTRAN 90
2070 ! MACHINE : IBM SP
2071 !
2072 !-----------------------------------------------------------------------
2073 !
2074 !
2075 !-----------------------------------------------------------------------
2076 IMPLICIT NONE
2077 !-----------------------------------------------------------------------
2078 !-------------------------------------------------------------------------
2079 !-------------- Parameters & arrays for lookup tables --------------------
2080 !-------------------------------------------------------------------------
2081 !
2082 !--- Common block of constants used in column microphysics
2083 !
2084 !WRF
2085 ! real DLMD,DPHD
2086 !WRF
2087 !
2088 !-----------------------------------------------------------------------
2089 !--- Parameters & data statement for local calculations
2090 !-----------------------------------------------------------------------
2091 !
2092 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
2093 !
2094 ! VARIABLES PASSED IN
2095 integer,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2096 & ,IMS,IME,JMS,JME,KMS,KME &
2097 & ,ITS,ITE,JTS,JTE,KTS,KTE
2098 !WRF
2099 INTEGER, DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: LOWLYR
2100 !
2101 real, INTENT(IN) :: DELX,DELY
2102 real,DIMENSION(*), INTENT(INOUT) :: MP_RESTART_STATE
2103 real,DIMENSION(NX), INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
2104 real,DIMENSION(ims:ime, kms:kme, jms:jme),INTENT(OUT) :: &
2105 & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
2106 INTEGER, PARAMETER :: ITLO=-60, ITHI=40
2107 ! integer,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
2108 ! real,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
2109 ! real,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
2110 ! real,INTENT(INOUT) :: PRECtot(2),PRECmax(2)
2111 real,INTENT(IN) :: DT,GSMDT
2112 LOGICAL,INTENT(IN) :: allowed_to_read,restart
2113 !
2114 !-----------------------------------------------------------------------
2115 ! LOCAL VARIABLES
2116 !-----------------------------------------------------------------------
2117 REAL :: BBFR,DTPH,PI,DX,Thour_print
2118 INTEGER :: I,IM,J,L,K,JTF,KTF,ITF
2119 INTEGER :: etampnew_unit1
2120 LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
2121 LOGICAL :: opened
2122 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
2123 CHARACTER*80 errmess
2124 !
2125 !-----------------------------------------------------------------------
2126 !
2127 JTF=MIN0(JTE,JDE-1)
2128 KTF=MIN0(KTE,KDE-1)
2129 ITF=MIN0(ITE,IDE-1)
2130 !
2131 DO J=JTS,JTF
2132 DO I=ITS,ITF
2133 LOWLYR(I,J)=1
2134 ENDDO
2135 ENDDO
2136 !
2137 IF(.NOT.RESTART)THEN
2138 DO J = jts,jte
2139 DO K = kts,kte
2140 DO I= its,ite
2141 F_ICE_PHY(i,k,j)=0.
2142 F_RAIN_PHY(i,k,j)=0.
2143 F_RIMEF_PHY(i,k,j)=1.
2144 ENDDO
2145 ENDDO
2146 ENDDO
2147 ENDIF
2148 !
2149 !-----------------------------------------------------------------------
2150 IF(ALLOWED_TO_READ)THEN
2151 !-----------------------------------------------------------------------
2152 !
2153 DX=((DELX)**2+(DELY)**2)**.5/1000. ! Model resolution at equator (km)
2154 DX=MIN(100., MAX(5., DX) )
2155 !
2156 !-- Relative humidity threshold for the onset of grid-scale condensation
2157 !!-- 9/1/01: Assume the following functional dependence for 5 - 100 km resolution:
2158 !! RHgrd=0.90 for dx=100 km, 0.98 for dx=5 km, where
2159 ! RHgrd=0.90+.08*((100.-DX)/95.)**.5
2160 !
2161 DTPH=MAX(GSMDT*60.,DT)
2162 DTPH=NINT(DTPH/DT)*DT
2163 !
2164 !--- Create lookup tables for saturation vapor pressure w/r/t water & ice
2165 !
2166 CALL GPVS
2167 !
2168 !--- Read in various lookup tables
2169 !
2170 IF ( wrf_dm_on_monitor() ) THEN
2171 DO i = 31,99
2172 INQUIRE ( i , OPENED = opened )
2173 IF ( .NOT. opened ) THEN
2174 etampnew_unit1 = i
2175 GOTO 2061
2176 ENDIF
2177 ENDDO
2178 etampnew_unit1 = -1
2179 2061 CONTINUE
2180 ENDIF
2181 !
2182 CALL wrf_dm_bcast_bytes ( etampnew_unit1 , IWORDSIZE )
2183 !
2184 IF ( etampnew_unit1 < 0 ) THEN
2185 CALL wrf_error_fatal ( 'module_mp_etanew: etanewinit: Can not find '// &
2186 'unused fortran unit to read in lookup table.' )
2187 ENDIF
2188 !
2189 IF ( wrf_dm_on_monitor() ) THEN
2190 !!was OPEN (UNIT=1,FILE="eta_micro_lookup.dat",FORM="UNFORMATTED")
2191 OPEN(UNIT=etampnew_unit1,FILE="ETAMPNEW_DATA", &
2192 & FORM="UNFORMATTED",STATUS="OLD",ERR=9061)
2193 !
2194 READ(etampnew_unit1) VENTR1
2195 READ(etampnew_unit1) VENTR2
2196 READ(etampnew_unit1) ACCRR
2197 READ(etampnew_unit1) MASSR
2198 READ(etampnew_unit1) VRAIN
2199 READ(etampnew_unit1) RRATE
2200 READ(etampnew_unit1) VENTI1
2201 READ(etampnew_unit1) VENTI2
2202 READ(etampnew_unit1) ACCRI
2203 READ(etampnew_unit1) MASSI
2204 READ(etampnew_unit1) VSNOWI
2205 READ(etampnew_unit1) VEL_RF
2206 ! read(etampnew_unit1) my_growth ! Applicable only for DTPH=180 s
2207 CLOSE (etampnew_unit1)
2208 ENDIF
2209 !
2210 CALL wrf_dm_bcast_bytes ( VENTR1 , size ( VENTR1 ) * RWORDSIZE )
2211 CALL wrf_dm_bcast_bytes ( VENTR2 , size ( VENTR2 ) * RWORDSIZE )
2212 CALL wrf_dm_bcast_bytes ( ACCRR , size ( ACCRR ) * RWORDSIZE )
2213 CALL wrf_dm_bcast_bytes ( MASSR , size ( MASSR ) * RWORDSIZE )
2214 CALL wrf_dm_bcast_bytes ( VRAIN , size ( VRAIN ) * RWORDSIZE )
2215 CALL wrf_dm_bcast_bytes ( RRATE , size ( RRATE ) * RWORDSIZE )
2216 CALL wrf_dm_bcast_bytes ( VENTI1 , size ( VENTI1 ) * RWORDSIZE )
2217 CALL wrf_dm_bcast_bytes ( VENTI2 , size ( VENTI2 ) * RWORDSIZE )
2218 CALL wrf_dm_bcast_bytes ( ACCRI , size ( ACCRI ) * RWORDSIZE )
2219 CALL wrf_dm_bcast_bytes ( MASSI , size ( MASSI ) * RWORDSIZE )
2220 CALL wrf_dm_bcast_bytes ( VSNOWI , size ( VSNOWI ) * RWORDSIZE )
2221 CALL wrf_dm_bcast_bytes ( VEL_RF , size ( VEL_RF ) * RWORDSIZE )
2222 !
2223 !--- Calculates coefficients for growth rates of ice nucleated in water
2224 ! saturated conditions, scaled by physics time step (lookup table)
2225 !
2226 CALL MY_GROWTH_RATES (DTPH)
2227 ! CALL MY_GROWTH_RATES (DTPH,MY_GROWTH)
2228 !
2229 PI=ACOS(-1.)
2230 !
2231 !--- Constants associated with Biggs (1953) freezing of rain, as parameterized
2232 ! following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
2233 !
2234 ABFR=-0.66
2235 BBFR=100.
2236 CBFR=20.*PI*PI*BBFR*RHOL*1.E-21
2237 !
2238 !--- CIACW is used in calculating riming rates
2239 ! The assumed effective collection efficiency of cloud water rimed onto
2240 ! ice is =0.5 below:
2241 !
2242 CIACW=DTPH*0.25*PI*0.5*(1.E5)**C1
2243 !
2244 !--- CIACR is used in calculating freezing of rain colliding with large ice
2245 ! The assumed collection efficiency is 1.0
2246 !
2247 CIACR=PI*DTPH
2248 !
2249 !--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
2250 ! * Four different functional relationships of mean drop diameter as
2251 ! a function of rain rate (RR), derived based on simple fits to
2252 ! mass-weighted fall speeds of rain as functions of mean diameter
2253 ! from the lookup tables.
2254 !
2255 RR_DRmin=N0r0*RRATE(MDRmin) ! RR for mean drop diameter of .05 mm
2256 RR_DR1=N0r0*RRATE(MDR1) ! RR for mean drop diameter of .10 mm
2257 RR_DR2=N0r0*RRATE(MDR2) ! RR for mean drop diameter of .20 mm
2258 RR_DR3=N0r0*RRATE(MDR3) ! RR for mean drop diameter of .32 mm
2259 RR_DRmax=N0r0*RRATE(MDRmax) ! RR for mean drop diameter of .45 mm
2260 !
2261 RQR_DRmin=N0r0*MASSR(MDRmin) ! Rain content for mean drop diameter of .05 mm
2262 RQR_DR1=N0r0*MASSR(MDR1) ! Rain content for mean drop diameter of .10 mm
2263 RQR_DR2=N0r0*MASSR(MDR2) ! Rain content for mean drop diameter of .20 mm
2264 RQR_DR3=N0r0*MASSR(MDR3) ! Rain content for mean drop diameter of .32 mm
2265 RQR_DRmax=N0r0*MASSR(MDRmax) ! Rain content for mean drop diameter of .45 mm
2266 C_N0r0=PI*RHOL*N0r0
2267 CN0r0=1.E6/C_N0r0**.25
2268 CN0r_DMRmin=1./(PI*RHOL*DMRmin**4)
2269 CN0r_DMRmax=1./(PI*RHOL*DMRmax**4)
2270 !
2271 !--- CRACW is used in calculating collection of cloud water by rain (an
2272 ! assumed collection efficiency of 1.0)
2273 !
2274 CRACW=DTPH*0.25*PI*1.0
2275 !
2276 ESW0=1000.*FPVS0(T0C) ! Saturation vapor pressure at 0C
2277 RFmax=1.1**Nrime ! Maximum rime factor allowed
2278 !
2279 !------------------------------------------------------------------------
2280 !--------------- Constants passed through argument list -----------------
2281 !------------------------------------------------------------------------
2282 !
2283 !--- Important parameters for self collection (autoconversion) of
2284 ! cloud water to rain.
2285 !
2286 !--- CRAUT is proportional to the rate that cloud water is converted by
2287 ! self collection to rain (autoconversion rate)
2288 !
2289 CRAUT=1.-EXP(-1.E-3*DTPH)
2290 !
2291 !--- QAUT0 is the threshold cloud content for autoconversion to rain
2292 ! needed for droplets to reach a diameter of 20 microns (following
2293 ! Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM)
2294 !--- QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for droplet number concentrations
2295 ! of 300, 200, and 100 cm**-3, respectively
2296 !
2297 QAUT0=PI*RHOL*NCW*(20.E-6)**3/6.
2298 !
2299 !--- For calculating snow optical depths by considering bulk density of
2300 ! snow based on emails from Q. Fu (6/27-28/01), where optical
2301 ! depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
2302 ! is effective radius, and DENS is the bulk density of snow.
2303 !
2304 ! SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
2305 ! T = 1.5*1.E3*SWPrad/(Reff*DENS)
2306 !
2307 ! See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
2308 !
2309 ! SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
2310 !
2311 DO I=MDImin,MDImax
2312 SDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I)
2313 ENDDO
2314 !
2315 Thour_print=-DTPH/3600.
2316
2317 ! SH 0211/2002
2318 ! IF (PRINT_diag) THEN
2319
2320 !-------- Total and maximum quantities
2321 !
2322 ! NSTATS=0 ! Microphysical statistics dealing w/ grid-point counts
2323 ! QMAX=0. ! Microphysical statistics dealing w/ hydrometeor mass
2324 ! QTOT=0. ! Microphysical statistics dealing w/ hydrometeor mass
2325 ! PRECmax=0. ! Maximum precip rates (rain, snow) at surface (mm/h)
2326 ! PRECtot=0. ! Total precipitation (rain, snow) accumulation at surface
2327 ! ENDIF
2328
2329 !wrf
2330 IF(.NOT.RESTART)THEN
2331 MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
2332 MP_RESTART_STATE(MY_T2+1)=C1XPVS0
2333 MP_RESTART_STATE(MY_T2+2)=C2XPVS0
2334 MP_RESTART_STATE(MY_T2+3)=C1XPVS
2335 MP_RESTART_STATE(MY_T2+4)=C2XPVS
2336 MP_RESTART_STATE(MY_T2+5)=CIACW
2337 MP_RESTART_STATE(MY_T2+6)=CIACR
2338 MP_RESTART_STATE(MY_T2+7)=CRACW
2339 MP_RESTART_STATE(MY_T2+8)=CRAUT
2340 TBPVS_STATE(1:NX) =TBPVS(1:NX)
2341 TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
2342 ENDIF
2343
2344 ENDIF ! Allowed_to_read
2345
2346 RETURN
2347 !
2348 !-----------------------------------------------------------------------
2349 !
2350 9061 CONTINUE
2351 WRITE( errmess , '(A,I4)' ) &
2352 'module_mp_etanew: error opening ETAMPNEW_DATA on unit ' &
2353 &, etampnew_unit1
2354 CALL wrf_error_fatal(errmess)
2355 !
2356 !-----------------------------------------------------------------------
2357 END SUBROUTINE etanewinit
2358 !
2359 SUBROUTINE MY_GROWTH_RATES (DTPH)
2360 ! SUBROUTINE MY_GROWTH_RATES (DTPH,MY_GROWTH)
2361 !
2362 !--- Below are tabulated values for the predicted mass of ice crystals
2363 ! after 600 s of growth in water saturated conditions, based on
2364 ! calculations from Miller and Young (JAS, 1979). These values are
2365 ! crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
2366 ! Young (1993). Values at temperatures colder than -27C were
2367 ! assumed to be invariant with temperature.
2368 !
2369 !--- Used to normalize Miller & Young (1979) calculations of ice growth
2370 ! over large time steps using their tabulated values at 600 s.
2371 ! Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
2372 !
2373 IMPLICIT NONE
2374 !
2375 REAL,INTENT(IN) :: DTPH
2376 !
2377 REAL DT_ICE
2378 REAL,DIMENSION(35) :: MY_600
2379 !WRF
2380 !
2381 !-----------------------------------------------------------------------
2382 DATA MY_600 / &
2383 & 5.5e-8, 1.4E-7, 2.8E-7, 6.E-7, 3.3E-6, &
2384 & 2.E-6, 9.E-7, 8.8E-7, 8.2E-7, 9.4e-7, &
2385 & 1.2E-6, 1.85E-6, 5.5E-6, 1.5E-5, 1.7E-5, &
2386 & 1.5E-5, 1.E-5, 3.4E-6, 1.85E-6, 1.35E-6, &
2387 & 1.05E-6, 1.E-6, 9.5E-7, 9.0E-7, 9.5E-7, &
2388 & 9.5E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7, &
2389 & 9.E-7, 9.E-7, 9.E-7, 9.E-7, 9.E-7 / ! -31 to -35 deg C
2390 !
2391 !-----------------------------------------------------------------------
2392 !
2393 DT_ICE=(DTPH/600.)**1.5
2394 MY_GROWTH=DT_ICE*MY_600
2395 !
2396 !-----------------------------------------------------------------------
2397 !
2398 END SUBROUTINE MY_GROWTH_RATES
2399 !
2400 !-----------------------------------------------------------------------
2401 !--------- Old GFS saturation vapor pressure lookup tables -----------
2402 !-----------------------------------------------------------------------
2403 !
2404 SUBROUTINE GPVS
2405 ! ******************************************************************
2406 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2407 ! . . .
2408 ! SUBPROGRAM: GPVS COMPUTE SATURATION VAPOR PRESSURE TABLE
2409 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2410 !
2411 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE TABLE AS A FUNCTION OF
2412 ! TEMPERATURE FOR THE TABLE LOOKUP FUNCTION FPVS.
2413 ! EXACT SATURATION VAPOR PRESSURES ARE CALCULATED IN SUBPROGRAM FPVSX.
2414 ! THE CURRENT IMPLEMENTATION COMPUTES A TABLE WITH A LENGTH
2415 ! OF 7501 FOR TEMPERATURES RANGING FROM 180.0 TO 330.0 KELVIN.
2416 !
2417 ! PROGRAM HISTORY LOG:
2418 ! 91-05-07 IREDELL
2419 ! 94-12-30 IREDELL EXPAND TABLE
2420 ! 96-02-19 HONG ICE EFFECT
2421 ! 01-11-29 JIN MODIFIED FOR WRF
2422 !
2423 ! USAGE: CALL GPVS
2424 !
2425 ! SUBPROGRAMS CALLED:
2426 ! (FPVSX) - INLINABLE FUNCTION TO COMPUTE SATURATION VAPOR PRESSURE
2427 !
2428 ! COMMON BLOCKS:
2429 ! COMPVS - SCALING PARAMETERS AND TABLE FOR FUNCTION FPVS.
2430 !
2431 ! ATTRIBUTES:
2432 ! LANGUAGE: FORTRAN 90
2433 !
2434 !$$$
2435 IMPLICIT NONE
2436 real :: X,XINC,T
2437 integer :: JX
2438 !----------------------------------------------------------------------
2439 XINC=(XMAX-XMIN)/(NX-1)
2440 C1XPVS=1.-XMIN/XINC
2441 C2XPVS=1./XINC
2442 C1XPVS0=1.-XMIN/XINC
2443 C2XPVS0=1./XINC
2444 !
2445 DO JX=1,NX
2446 X=XMIN+(JX-1)*XINC
2447 T=X
2448 TBPVS(JX)=FPVSX(T)
2449 TBPVS0(JX)=FPVSX0(T)
2450 ENDDO
2451 !
2452 END SUBROUTINE GPVS
2453 !-----------------------------------------------------------------------
2454 !***********************************************************************
2455 !-----------------------------------------------------------------------
2456 REAL FUNCTION FPVS(T)
2457 !-----------------------------------------------------------------------
2458 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2459 ! . . .
2460 ! SUBPROGRAM: FPVS COMPUTE SATURATION VAPOR PRESSURE
2461 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2462 !
2463 ! ABSTRACT: COMPUTE SATURATION VAPOR PRESSURE FROM THE TEMPERATURE.
2464 ! A LINEAR INTERPOLATION IS DONE BETWEEN VALUES IN A LOOKUP TABLE
2465 ! COMPUTED IN GPVS. SEE DOCUMENTATION FOR FPVSX FOR DETAILS.
2466 ! INPUT VALUES OUTSIDE TABLE RANGE ARE RESET TO TABLE EXTREMA.
2467 ! THE INTERPOLATION ACCURACY IS ALMOST 6 DECIMAL PLACES.
2468 ! ON THE CRAY, FPVS IS ABOUT 4 TIMES FASTER THAN EXACT CALCULATION.
2469 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2470 !
2471 ! PROGRAM HISTORY LOG:
2472 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2473 ! 94-12-30 IREDELL EXPAND TABLE
2474 ! 96-02-19 HONG ICE EFFECT
2475 ! 01-11-29 JIN MODIFIED FOR WRF
2476 !
2477 ! USAGE: PVS=FPVS(T)
2478 !
2479 ! INPUT ARGUMENT LIST:
2480 ! T - REAL TEMPERATURE IN KELVIN
2481 !
2482 ! OUTPUT ARGUMENT LIST:
2483 ! FPVS - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2484 !
2485 ! ATTRIBUTES:
2486 ! LANGUAGE: FORTRAN 90
2487 !$$$
2488 IMPLICIT NONE
2489 real,INTENT(IN) :: T
2490 real XJ
2491 integer :: JX
2492 !-----------------------------------------------------------------------
2493 XJ=MIN(MAX(C1XPVS+C2XPVS*T,1.),FLOAT(NX))
2494 JX=MIN(XJ,NX-1.)
2495 FPVS=TBPVS(JX)+(XJ-JX)*(TBPVS(JX+1)-TBPVS(JX))
2496 !
2497 END FUNCTION FPVS
2498 !-----------------------------------------------------------------------
2499 !-----------------------------------------------------------------------
2500 REAL FUNCTION FPVS0(T)
2501 !-----------------------------------------------------------------------
2502 IMPLICIT NONE
2503 real,INTENT(IN) :: T
2504 real :: XJ1
2505 integer :: JX1
2506 !-----------------------------------------------------------------------
2507 XJ1=MIN(MAX(C1XPVS0+C2XPVS0*T,1.),FLOAT(NX))
2508 JX1=MIN(XJ1,NX-1.)
2509 FPVS0=TBPVS0(JX1)+(XJ1-JX1)*(TBPVS0(JX1+1)-TBPVS0(JX1))
2510 !
2511 END FUNCTION FPVS0
2512 !-----------------------------------------------------------------------
2513 !***********************************************************************
2514 !-----------------------------------------------------------------------
2515 REAL FUNCTION FPVSX(T)
2516 !-----------------------------------------------------------------------
2517 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2518 ! . . .
2519 ! SUBPROGRAM: FPVSX COMPUTE SATURATION VAPOR PRESSURE
2520 ! AUTHOR: N PHILLIPS W/NP2 DATE: 30 DEC 82
2521 !
2522 ! ABSTRACT: EXACTLY COMPUTE SATURATION VAPOR PRESSURE FROM TEMPERATURE.
2523 ! THE WATER MODEL ASSUMES A PERFECT GAS, CONSTANT SPECIFIC HEATS
2524 ! FOR GAS AND LIQUID, AND NEGLECTS THE VOLUME OF THE LIQUID.
2525 ! THE MODEL DOES ACCOUNT FOR THE VARIATION OF THE LATENT HEAT
2526 ! OF CONDENSATION WITH TEMPERATURE. THE ICE OPTION IS NOT INCLUDED.
2527 ! THE CLAUSIUS-CLAPEYRON EQUATION IS INTEGRATED FROM THE TRIPLE POINT
2528 ! TO GET THE FORMULA
2529 ! PVS=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2530 ! WHERE TR IS TTP/T AND OTHER VALUES ARE PHYSICAL CONSTANTS
2531 ! THIS FUNCTION SHOULD BE EXPANDED INLINE IN THE CALLING ROUTINE.
2532 !
2533 ! PROGRAM HISTORY LOG:
2534 ! 91-05-07 IREDELL MADE INTO INLINABLE FUNCTION
2535 ! 94-12-30 IREDELL EXACT COMPUTATION
2536 ! 96-02-19 HONG ICE EFFECT
2537 ! 01-11-29 JIN MODIFIED FOR WRF
2538 !
2539 ! USAGE: PVS=FPVSX(T)
2540 ! REFERENCE: EMANUEL(1994),116-117
2541 !
2542 ! INPUT ARGUMENT LIST:
2543 ! T - REAL TEMPERATURE IN KELVIN
2544 !
2545 ! OUTPUT ARGUMENT LIST:
2546 ! FPVSX - REAL SATURATION VAPOR PRESSURE IN KILOPASCALS (CB)
2547 !
2548 ! ATTRIBUTES:
2549 ! LANGUAGE: FORTRAN 90
2550 !$$$
2551 IMPLICIT NONE
2552 !-----------------------------------------------------------------------
2553 real, parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2554 , CLIQ=4.1855E+3,CVAP= 1.8460E+3 &
2555 , CICE=2.1060E+3,HSUB=2.8340E+6
2556 !
2557 real, parameter :: PSATK=PSAT*1.E-3
2558 real, parameter :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2559 real, parameter :: DLDTI=CVAP-CICE &
2560 , XAI=-DLDTI/RV,XBI=XAI+HSUB/(RV*TTP)
2561 real T,TR
2562 !-----------------------------------------------------------------------
2563 TR=TTP/T
2564 !
2565 IF(T.GE.TTP)THEN
2566 FPVSX=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2567 ELSE
2568 FPVSX=PSATK*(TR**XAI)*EXP(XBI*(1.-TR))
2569 ENDIF
2570 !
2571 END FUNCTION FPVSX
2572 !-----------------------------------------------------------------------
2573 !-----------------------------------------------------------------------
2574 REAL FUNCTION FPVSX0(T)
2575 !-----------------------------------------------------------------------
2576 IMPLICIT NONE
2577 real,parameter :: TTP=2.7316E+2,HVAP=2.5000E+6,PSAT=6.1078E+2 &
2578 , CLIQ=4.1855E+3,CVAP=1.8460E+3 &
2579 , CICE=2.1060E+3,HSUB=2.8340E+6
2580 real,PARAMETER :: PSATK=PSAT*1.E-3
2581 real,PARAMETER :: DLDT=CVAP-CLIQ,XA=-DLDT/RV,XB=XA+HVAP/(RV*TTP)
2582 real,PARAMETER :: DLDTI=CVAP-CICE &
2583 , XAI=-DLDT/RV,XBI=XA+HSUB/(RV*TTP)
2584 real :: T,TR
2585 !-----------------------------------------------------------------------
2586 TR=TTP/T
2587 FPVSX0=PSATK*(TR**XA)*EXP(XB*(1.-TR))
2588 !
2589 END FUNCTION FPVSX0
2590 !
2591 END MODULE module_mp_etanew