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