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