PageRenderTime 58ms CodeModel.GetById 19ms RepoModel.GetById 1ms app.codeStats 0ms

/wrfv2_fire/phys/module_mp_HWRF.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2653 lines | 1168 code | 55 blank | 1430 comment | 0 complexity | d9d0c6dbc1e2fdaab2f34f4034d6d2c2 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

  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
  102. !-----------------------------------------------------------------------
  103. !-----------------------------------------------------------------------
  104. SUBROUTINE ETAMP_NEW_HWRF (itimestep,DT,DX,DY,GID,RAINNC,RAINNCV, & !GID
  105. & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, & !gopal's doing
  106. & LOWLYR,SR, &
  107. & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
  108. & QC,QR,QI, &
  109. & ids,ide, jds,jde, kds,kde, &
  110. & ims,ime, jms,jme, kms,kme, &
  111. & its,ite, jts,jte, kts,kte )
  112. !HWRF SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
  113. !HWRF & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qc, &
  114. !HWRF & LOWLYR,SR, &
  115. !HWRF & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
  116. !HWRF & mp_restart_state,tbpvs_state,tbpvs0_state, &
  117. !HWRF & RAINNC,RAINNCV, &
  118. !HWRF & ids,ide, jds,jde, kds,kde, &
  119. !HWRF & ims,ime, jms,jme, kms,kme, &
  120. !HWRF & its,ite, jts,jte, kts,kte )
  121. !-----------------------------------------------------------------------
  122. IMPLICIT NONE
  123. !-----------------------------------------------------------------------
  124. INTEGER, PARAMETER :: ITLO=-60, ITHI=40
  125. INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
  126. & ,IMS,IME,JMS,JME,KMS,KME &
  127. & ,ITS,ITE,JTS,JTE,KTS,KTE &
  128. & ,ITIMESTEP,GID ! GID gopal's doing
  129. REAL, INTENT(IN) :: DT,DX,DY
  130. REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
  131. & dz8w,p_phy,pi_phy,rho_phy
  132. REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
  133. & th_phy,qv,qt,qc,qr,qi
  134. REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
  135. & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
  136. REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
  137. & RAINNC,RAINNCV
  138. REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
  139. !
  140. !HWRF REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
  141. !
  142. !HWRF REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
  143. !
  144. INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
  145. !-----------------------------------------------------------------------
  146. ! LOCAL VARS
  147. !-----------------------------------------------------------------------
  148. ! NSTATS,QMAX,QTOT are diagnostic vars
  149. INTEGER,DIMENSION(ITLO:ITHI,4) :: NSTATS
  150. REAL, DIMENSION(ITLO:ITHI,5) :: QMAX
  151. REAL, DIMENSION(ITLO:ITHI,22):: QTOT
  152. ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
  153. ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
  154. ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
  155. ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
  156. ! the microphysics scheme. Instead, they will be used by Eta precip
  157. ! assimilation.
  158. REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
  159. & TLATGS_PHY,TRAIN_PHY
  160. REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
  161. REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
  162. INTEGER :: I,J,K,KFLIP
  163. !
  164. !-----------------------------------------------------------------------
  165. !**********************************************************************
  166. !-----------------------------------------------------------------------
  167. !
  168. !HWRF MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
  169. !HWRF!
  170. !HWRF C1XPVS0=MP_RESTART_STATE(MY_T2+1)
  171. !HWRF C2XPVS0=MP_RESTART_STATE(MY_T2+2)
  172. !HWRF C1XPVS =MP_RESTART_STATE(MY_T2+3)
  173. !HWRF C2XPVS =MP_RESTART_STATE(MY_T2+4)
  174. !HWRF CIACW =MP_RESTART_STATE(MY_T2+5)
  175. !HWRF CIACR =MP_RESTART_STATE(MY_T2+6)
  176. !HWRF CRACW =MP_RESTART_STATE(MY_T2+7)
  177. !HWRF CRAUT =MP_RESTART_STATE(MY_T2+8)
  178. !HWRF!
  179. !HWRF TBPVS(1:NX) =TBPVS_STATE(1:NX)
  180. !HWRF TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
  181. !
  182. DO j = jts,jte
  183. DO k = kts,kte
  184. DO i = its,ite
  185. t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
  186. qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
  187. ENDDO
  188. ENDDO
  189. ENDDO
  190. ! initial diagnostic variables and data assimilation vars
  191. ! (will need to delete this part in the future)
  192. DO k = 1,4
  193. DO i = ITLO,ITHI
  194. NSTATS(i,k)=0.
  195. ENDDO
  196. ENDDO
  197. DO k = 1,5
  198. DO i = ITLO,ITHI
  199. QMAX(i,k)=0.
  200. ENDDO
  201. ENDDO
  202. DO k = 1,22
  203. DO i = ITLO,ITHI
  204. QTOT(i,k)=0.
  205. ENDDO
  206. ENDDO
  207. ! initial data assimilation vars (will need to delete this part in the future)
  208. DO j = jts,jte
  209. DO k = kts,kte
  210. DO i = its,ite
  211. TLATGS_PHY (i,k,j)=0.
  212. TRAIN_PHY (i,k,j)=0.
  213. ENDDO
  214. ENDDO
  215. ENDDO
  216. DO j = jts,jte
  217. DO i = its,ite
  218. ACPREC(i,j)=0.
  219. APREC (i,j)=0.
  220. PREC (i,j)=0.
  221. SR (i,j)=0.
  222. ENDDO
  223. ENDDO
  224. !-- 6/11/2010: Update QT, F_ice, F_rain arrays
  225. DO j = jts,jte
  226. DO k = kts,kte
  227. DO i = its,ite
  228. QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QI(I,K,J)
  229. IF (QI(I,K,J) <= EPSQ) THEN
  230. F_ICE_PHY(I,K,J)=0.
  231. IF (T_PHY(I,K,J) < T_ICEK) F_ICE_PHY(I,K,J)=1.
  232. ELSE
  233. F_ICE_PHY(I,K,J)=MAX( 0., MIN(1., QI(I,K,J)/QT(I,K,J) ) )
  234. ENDIF
  235. IF (QR(I,K,J) <= EPSQ) THEN
  236. F_RAIN_PHY(I,K,J)=0.
  237. ELSE
  238. F_RAIN_PHY(I,K,J)=QR(I,K,J)/(QR(I,K,J)+QC(I,K,J))
  239. ENDIF
  240. ENDDO
  241. ENDDO
  242. ENDDO
  243. !-----------------------------------------------------------------------
  244. CALL EGCP01DRV(GID,DT,LOWLYR, &
  245. & APREC,PREC,ACPREC,SR,NSTATS,QMAX,QTOT, &
  246. & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
  247. & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
  248. & ids,ide, jds,jde, kds,kde, &
  249. & ims,ime, jms,jme, kms,kme, &
  250. & its,ite, jts,jte, kts,kte )
  251. !-----------------------------------------------------------------------
  252. DO j = jts,jte
  253. DO k = kts,kte
  254. DO i = its,ite
  255. th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
  256. qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
  257. WC=qt(I,K,J)
  258. QI(I,K,J)=0.
  259. QR(I,K,J)=0.
  260. QC(I,K,J)=0.
  261. IF(F_ICE_PHY(I,K,J)>=1.)THEN
  262. QI(I,K,J)=WC
  263. ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
  264. QC(I,K,J)=WC
  265. ELSE
  266. QI(I,K,J)=F_ICE_PHY(I,K,J)*WC
  267. QC(I,K,J)=WC-QI(I,K,J)
  268. ENDIF
  269. !
  270. IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
  271. IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
  272. QR(I,K,J)=QC(I,K,J)
  273. QC(I,K,J)=0.
  274. ELSE
  275. QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
  276. QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
  277. ENDIF
  278. endif
  279. ENDDO
  280. ENDDO
  281. ENDDO
  282. !
  283. ! update rain (from m to mm)
  284. DO j=jts,jte
  285. DO i=its,ite
  286. RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
  287. RAINNCV(i,j)=APREC(i,j)*1000.
  288. ENDDO
  289. ENDDO
  290. !
  291. !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
  292. !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
  293. !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
  294. !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
  295. !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
  296. !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
  297. !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
  298. !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
  299. !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
  300. !HWRF!
  301. !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
  302. !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
  303. !-----------------------------------------------------------------------
  304. END SUBROUTINE ETAMP_NEW_HWRF
  305. !-----------------------------------------------------------------------
  306. SUBROUTINE EGCP01DRV(GID, & !GID gopal's doing
  307. & DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
  308. & NSTATS,QMAX,QTOT, &
  309. & dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
  310. & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
  311. & ids,ide, jds,jde, kds,kde, &
  312. & ims,ime, jms,jme, kms,kme, &
  313. & its,ite, jts,jte, kts,kte)
  314. !-----------------------------------------------------------------------
  315. ! DTPH Physics time step (s)
  316. ! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
  317. ! Q_PHY Mixing ratio of water vapor. kg/kg
  318. ! F_RAIN_PHY Fraction of rain.
  319. ! F_ICE_PHY Fraction of ice.
  320. ! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
  321. !
  322. !TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
  323. !micrphysics sechme. Instead, they will be used by Eta precip assimilation.
  324. !
  325. !NSTATS,QMAX,QTOT are used for diagnosis purposes.
  326. !
  327. !-----------------------------------------------------------------------
  328. !--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
  329. ! and/or ZHAO's scheme in Eta and are not required by this microphysics
  330. ! scheme itself.
  331. !--- NSTATS,QMAX,QTOT are used for diagnosis purposes only. They will be
  332. ! printed out when PRINT_diag is true.
  333. !
  334. !-----------------------------------------------------------------------
  335. IMPLICIT NONE
  336. !-----------------------------------------------------------------------
  337. !
  338. INTEGER, PARAMETER :: ITLO=-60, ITHI=40
  339. ! VARIABLES PASSED IN/OUT
  340. INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
  341. & ,ims,ime, jms,jme, kms,kme &
  342. & ,its,ite, jts,jte, kts,kte
  343. INTEGER,INTENT(IN ) :: GID ! grid%id gopal's doing
  344. REAL,INTENT(IN) :: DTPH
  345. INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
  346. INTEGER,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
  347. REAL,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
  348. REAL,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
  349. REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
  350. & APREC,PREC,ACPREC,SR
  351. REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
  352. REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
  353. & dz8w,P_PHY,RHO_PHY
  354. REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
  355. & CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
  356. & ,Q_PHY,TRAIN_PHY
  357. !
  358. !-----------------------------------------------------------------------
  359. !LOCAL VARIABLES
  360. !-----------------------------------------------------------------------
  361. !
  362. !HWRF - Below are directives in the operational code that have been removed,
  363. ! where "TEMP_DEX" has been replaced with "I,J,L" and "TEMP_DIMS" has
  364. ! been replaced with "its:ite,jts:jte,kts:kte"
  365. !HWRF#define CACHE_FRIENDLY_MP_ETANEW
  366. !HWRF#ifdef CACHE_FRIENDLY_MP_ETANEW
  367. !HWRF# define TEMP_DIMS kts:kte,its:ite,jts:jte
  368. !HWRF# define TEMP_DEX L,I,J
  369. !HWRF#else
  370. !HWRF# define TEMP_DIMS its:ite,jts:jte,kts:kte
  371. !HWRF# define TEMP_DEX I,J,L
  372. !HWRF#endif
  373. !HWRF!
  374. INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
  375. !HWRF REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
  376. REAL,DIMENSION(its:ite,jts:jte,kts:kte) :: &
  377. & CWM,T,Q,TRAIN,TLATGS,P
  378. REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
  379. INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
  380. REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
  381. REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
  382. RimeF_col,QI_col,QR_col,QW_col, THICK_col, RHC_col, DPCOL !GFDL
  383. REAL,DIMENSION(2) :: PRECtot,PRECmax
  384. !-----------------------------------------------------------------------
  385. !
  386. DO J=JTS,JTE
  387. DO I=ITS,ITE
  388. LMH(I,J) = KTE-LOWLYR(I,J)+1
  389. ENDDO
  390. ENDDO
  391. DO 98 J=JTS,JTE
  392. DO 98 I=ITS,ITE
  393. DO L=KTS,KTE
  394. KFLIP=KTE+1-L
  395. CWM(I,J,L)=CWM_PHY(I,KFLIP,J)
  396. T(I,J,L)=T_PHY(I,KFLIP,J)
  397. Q(I,J,L)=Q_PHY(I,KFLIP,J)
  398. P(I,J,L)=P_PHY(I,KFLIP,J)
  399. TLATGS(I,J,L)=TLATGS_PHY(I,KFLIP,J)
  400. TRAIN(I,J,L)=TRAIN_PHY(I,KFLIP,J)
  401. F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
  402. F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
  403. F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
  404. ENDDO
  405. 98 CONTINUE
  406. DO 100 J=JTS,JTE
  407. DO 100 I=ITS,ITE
  408. LSFC=LMH(I,J) ! "L" of surface
  409. !
  410. DO K=KTS,KTE
  411. KFLIP=KTE+1-K
  412. DPCOL(K)=RHO_PHY(I,KFLIP,J)*GRAV*dz8w(I,KFLIP,J)
  413. ENDDO
  414. !
  415. !
  416. !--- Initialize column data (1D arrays)
  417. !
  418. IF (CWM(I,J,1) .LE. EPSQ) CWM(I,J,1)=EPSQ
  419. F_ice(1,I,J)=1.
  420. F_rain(1,I,J)=0.
  421. F_RimeF(1,I,J)=1.
  422. DO L=1,LSFC
  423. !
  424. !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
  425. !
  426. P_col(L)=P(I,J,L)
  427. !
  428. !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
  429. !
  430. THICK_col(L)=DPCOL(L)*RGRAV
  431. T_col(L)=T(I,J,L)
  432. TC=T_col(L)-T0C
  433. QV_col(L)=max(EPSQ, Q(I,J,L))
  434. IF (CWM(I,J,L) .LE. EPSQ1) THEN
  435. WC_col(L)=0.
  436. IF (TC .LT. T_ICE) THEN
  437. F_ice(L,I,J)=1.
  438. ELSE
  439. F_ice(L,I,J)=0.
  440. ENDIF
  441. F_rain(L,I,J)=0.
  442. F_RimeF(L,I,J)=1.
  443. ELSE
  444. WC_col(L)=CWM(I,J,L)
  445. ENDIF
  446. !
  447. !--- Determine composition of condensate in terms of
  448. ! cloud water, ice, & rain
  449. !
  450. WC=WC_col(L)
  451. QI=0.
  452. QR=0.
  453. QW=0.
  454. Fice=F_ice(L,I,J)
  455. Frain=F_rain(L,I,J)
  456. IF (Fice .GE. 1.) THEN
  457. QI=WC
  458. ELSE IF (Fice .LE. 0.) THEN
  459. QW=WC
  460. ELSE
  461. QI=Fice*WC
  462. QW=WC-QI
  463. ENDIF
  464. IF (QW.GT.0. .AND. Frain.GT.0.) THEN
  465. IF (Frain .GE. 1.) THEN
  466. QR=QW
  467. QW=0.
  468. ELSE
  469. QR=Frain*QW
  470. QW=QW-QR
  471. ENDIF
  472. ENDIF
  473. RimeF_col(L)=F_RimeF(L,I,J)
  474. QI_col(L)=QI
  475. QR_col(L)=QR
  476. QW_col(L)=QW
  477. !GFDL => New. Added RHC_col to allow for height- and grid-dependent values for
  478. !GFDL the relative humidity threshold for condensation ("RHgrd")
  479. !6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa
  480. !------------------------------------------------------------
  481. IF(GID .EQ. 1 .AND. P_col(L)<P_RHgrd_out) THEN ! gopal's doing based on GFDL
  482. RHC_col(L)=RHgrd_out
  483. ELSE
  484. RHC_col(L)=RHgrd_in
  485. ENDIF
  486. !------------------------------------------------------------
  487. ENDDO
  488. !
  489. !#######################################################################
  490. !
  491. !--- Perform the microphysical calculations in this column
  492. !
  493. I_index=I
  494. J_index=J
  495. CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
  496. & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
  497. & THICK_col, WC_col, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT ) !GFDL
  498. !
  499. !#######################################################################
  500. !
  501. !
  502. !--- Update storage arrays
  503. !
  504. DO L=1,LSFC
  505. TRAIN(I,J,L)=(T_col(L)-T(I,J,L))/DTPH
  506. TLATGS(I,J,L)=T_col(L)-T(I,J,L)
  507. T(I,J,L)=T_col(L)
  508. Q(I,J,L)=QV_col(L)
  509. CWM(I,J,L)=WC_col(L)
  510. !
  511. !--- REAL*4 array storage
  512. !
  513. F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
  514. IF (QI_col(L) .LE. EPSQ) THEN
  515. F_ice(L,I,J)=0.
  516. IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
  517. ELSE
  518. F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
  519. ENDIF
  520. IF (QR_col(L) .LE. EPSQ) THEN
  521. DUM=0
  522. ELSE
  523. DUM=QR_col(L)/(QR_col(L)+QW_col(L))
  524. ENDIF
  525. F_rain(L,I,J)=DUM
  526. !
  527. ENDDO
  528. !
  529. !--- Update accumulated precipitation statistics
  530. !
  531. !--- Surface precipitation statistics; SR is fraction of surface
  532. ! precipitation (if >0) associated with snow
  533. !
  534. APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
  535. PREC(I,J)=PREC(I,J)+APREC(I,J)
  536. ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
  537. IF(APREC(I,J) .LT. 1.E-8) THEN
  538. SR(I,J)=0.
  539. ELSE
  540. SR(I,J)=RRHOL*ASNOW/APREC(I,J)
  541. ENDIF
  542. ! !
  543. ! !--- Debug statistics
  544. ! !
  545. ! IF (PRINT_diag) THEN
  546. ! PRECtot(1)=PRECtot(1)+ARAIN
  547. ! PRECtot(2)=PRECtot(2)+ASNOW
  548. ! PRECmax(1)=MAX(PRECmax(1), ARAIN)
  549. ! PRECmax(2)=MAX(PRECmax(2), ASNOW)
  550. ! ENDIF
  551. !#######################################################################
  552. !#######################################################################
  553. !
  554. 100 CONTINUE ! End "I" & "J" loops
  555. DO 101 J=JTS,JTE
  556. DO 101 I=ITS,ITE
  557. DO L=KTS,KTE
  558. KFLIP=KTE+1-L
  559. CWM_PHY(I,KFLIP,J)=CWM(I,J,L)
  560. T_PHY(I,KFLIP,J)=T(I,J,L)
  561. Q_PHY(I,KFLIP,J)=Q(I,J,L)
  562. TLATGS_PHY(I,KFLIP,J)=TLATGS(I,J,L)
  563. TRAIN_PHY(I,KFLIP,J)=TRAIN(I,J,L)
  564. F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
  565. F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
  566. F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
  567. ENDDO
  568. 101 CONTINUE
  569. END SUBROUTINE EGCP01DRV
  570. !
  571. !
  572. !###############################################################################
  573. ! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
  574. ! (1) Represents sedimentation by preserving a portion of the precipitation
  575. ! through top-down integration from cloud-top. Modified procedure to
  576. ! Zhao and Carr (1997).
  577. ! (2) Microphysical equations are modified to be less sensitive to time
  578. ! steps by use of Clausius-Clapeyron equation to account for changes in
  579. ! saturation mixing ratios in response to latent heating/cooling.
  580. ! (3) Prevent spurious temperature oscillations across 0C due to
  581. ! microphysics.
  582. ! (4) Uses lookup tables for: calculating two different ventilation
  583. ! coefficients in condensation and deposition processes; accretion of
  584. ! cloud water by precipitation; precipitation mass; precipitation rate
  585. ! (and mass-weighted precipitation fall speeds).
  586. ! (5) Assumes temperature-dependent variation in mean diameter of large ice
  587. ! (Houze et al., 1979; Ryan et al., 1996).
  588. ! -> 8/22/01: This relationship has been extended to colder temperatures
  589. ! to parameterize smaller large-ice particles down to mean sizes of MDImin,
  590. ! which is 50 microns reached at -55.9C.
  591. ! (6) Attempts to differentiate growth of large and small ice, mainly for
  592. ! improved transition from thin cirrus to thick, precipitating ice
  593. ! anvils.
  594. ! -> 8/22/01: This feature has been diminished by effectively adjusting to
  595. ! ice saturation during depositional growth at temperatures colder than
  596. ! -10C. Ice sublimation is calculated more explicitly. The logic is
  597. ! that sources of are either poorly understood (e.g., nucleation for NWP)
  598. ! or are not represented in the Eta model (e.g., detrainment of ice from
  599. ! convection). Otherwise the model is too wet compared to the radiosonde
  600. ! observations based on 1 Feb - 18 March 2001 retrospective runs.
  601. ! (7) Top-down integration also attempts to treat mixed-phase processes,
  602. ! allowing a mixture of ice and water. Based on numerous observational
  603. ! studies, ice growth is based on nucleation at cloud top &
  604. ! subsequent growth by vapor deposition and riming as the ice particles
  605. ! fall through the cloud. Effective nucleation rates are a function
  606. ! of ice supersaturation following Meyers et al. (JAM, 1992).
  607. ! -> 8/22/01: The simulated relative humidities were far too moist compared
  608. ! to the rawinsonde observations. This feature has been substantially
  609. ! diminished, limited to a much narrower temperature range of 0 to -10C.
  610. ! (8) Depositional growth of newly nucleated ice is calculated for large time
  611. ! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
  612. ! using their ice crystal masses calculated after 600 s of growth in water
  613. ! saturated conditions. The growth rates are normalized by time step
  614. ! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
  615. ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
  616. ! (9) Ice precipitation rates can increase due to increase in response to
  617. ! cloud water riming due to (a) increased density & mass of the rimed
  618. ! ice, and (b) increased fall speeds of rimed ice.
  619. ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
  620. !###############################################################################
  621. !###############################################################################
  622. !
  623. SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, &
  624. & LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
  625. & THICK_col, WC_col, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT) !GFDL
  626. !
  627. !###############################################################################
  628. !###############################################################################
  629. !
  630. !-------------------------------------------------------------------------------
  631. !----- NOTE: Code is currently set up w/o threading!
  632. !-------------------------------------------------------------------------------
  633. !$$$ SUBPROGRAM DOCUMENTATION BLOCK
  634. ! . . .
  635. ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
  636. ! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
  637. ! PRGRMMR: Jin (Modification for WRF structure)
  638. !-------------------------------------------------------------------------------
  639. ! ABSTRACT:
  640. ! * Merges original GSCOND & PRECPD subroutines.
  641. ! * Code has been substantially streamlined and restructured.
  642. ! * Exchange between water vapor & small cloud condensate is calculated using
  643. ! the original Asai (1965, J. Japan) algorithm. See also references to
  644. ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
  645. ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
  646. ! parameterization.
  647. !-------------------------------------------------------------------------------
  648. !
  649. ! USAGE:
  650. ! * CALL EGCP01COLUMN FROM SUBROUTINE EGCP01DRV
  651. !
  652. ! INPUT ARGUMENT LIST:
  653. ! DTPH - physics time step (s)
  654. ! I_index - I index
  655. ! J_index - J index
  656. ! LSFC - Eta level of level above surface, ground
  657. ! P_col - vertical column of model pressure (Pa)
  658. ! QI_col - vertical column of model ice mixing ratio (kg/kg)
  659. ! QR_col - vertical column of model rain ratio (kg/kg)
  660. ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
  661. ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
  662. ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
  663. ! T_col - vertical column of model temperature (deg K)
  664. ! THICK_col - vertical column of model mass thickness (density*height increment)
  665. ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
  666. ! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL
  667. !
  668. !
  669. ! OUTPUT ARGUMENT LIST:
  670. ! ARAIN - accumulated rainfall at the surface (kg)
  671. ! ASNOW - accumulated snowfall at the surface (kg)
  672. ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
  673. ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
  674. ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
  675. ! QI_col - vertical column of model ice mixing ratio (kg/kg)
  676. ! QR_col - vertical column of model rain ratio (kg/kg)
  677. ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
  678. ! T_col - vertical column of model temperature (deg K)
  679. !
  680. ! OUTPUT FILES:
  681. ! NONE
  682. !
  683. ! Subprograms & Functions called:
  684. ! * Real Function CONDENSE - cloud water condensation
  685. ! * Real Function DEPOSIT - ice deposition (not sublimation)
  686. !
  687. ! UNIQUE: NONE
  688. !
  689. ! LIBRARY: NONE
  690. !
  691. ! COMMON BLOCKS:
  692. ! CMICRO_CONS - key constants initialized in GSMCONST
  693. ! CMICRO_STATS - accumulated and maximum statistics
  694. ! CMY_GROWTH - lookup table for growth of ice crystals in
  695. ! water saturated conditions (Miller & Young, 1979)
  696. ! IVENT_TABLES - lookup tables for ventilation effects of ice
  697. ! IACCR_TABLES - lookup tables for accretion rates of ice
  698. ! IMASS_TABLES - lookup tables for mass content of ice
  699. ! IRATE_TABLES - lookup tables for precipitation rates of ice
  700. ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
  701. ! RVENT_TABLES - lookup tables for ventilation effects of rain
  702. ! RACCR_TABLES - lookup tables for accretion rates of rain
  703. ! RMASS_TABLES - lookup tables for mass content of rain
  704. ! RVELR_TABLES - lookup tables for fall speeds of rain
  705. ! RRATE_TABLES - lookup tables for precipitation rates of rain
  706. !
  707. ! ATTRIBUTES:
  708. ! LANGUAGE: FORTRAN 90
  709. ! MACHINE : IBM SP
  710. !
  711. !
  712. !-------------------------------------------------------------------------
  713. !--------------- Arrays & constants in argument list ---------------------
  714. !-------------------------------------------------------------------------
  715. !
  716. IMPLICIT NONE
  717. !
  718. INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
  719. REAL,INTENT(INOUT) :: ARAIN, ASNOW
  720. REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col, QI_col,QR_col &
  721. & ,QV_col ,QW_col, RimeF_col, T_col, THICK_col, WC_col, RHC_col !GFDL
  722. !
  723. !-------------------------------------------------------------------------
  724. !-------------- Common blocks for microphysical statistics ---------------
  725. !-------------------------------------------------------------------------
  726. !
  727. !-------------------------------------------------------------------------
  728. !--------- Common blocks for constants initialized in GSMCONST ----------
  729. !
  730. INTEGER, PARAMETER :: ITLO=-60, ITHI=40
  731. INTEGER,INTENT(INOUT) :: NSTATS(ITLO:ITHI,4)
  732. REAL,INTENT(INOUT) :: QMAX(ITLO:ITHI,5),QTOT(ITLO:ITHI,22)
  733. !
  734. !-------------------------------------------------------------------------
  735. !--------------- Common blocks for various lookup tables -----------------
  736. !
  737. !--- Discretized growth rates of small ice crystals after their nucleation
  738. ! at 1 C intervals from -1 C to -35 C, based on calculations by Miller
  739. ! and Young (1979, JAS) after 600 s of growth. Resultant growth rates
  740. ! are multiplied by physics time step in GSMCONST.
  741. !
  742. !-------------------------------------------------------------------------
  743. !
  744. !--- Mean ice-particle diameters varying from 50 microns to 1000 microns
  745. ! (1 mm), assuming an exponential size distribution.
  746. !
  747. !---- Meaning of the following arrays:
  748. ! - mdiam - mean diameter (m)
  749. ! - VENTI1 - integrated quantity associated w/ ventilation effects
  750. ! (capacitance only) for calculating vapor deposition onto ice
  751. ! - VENTI2 - integrated quantity associated w/ ventilation effects
  752. ! (with fall speed) for calculating vapor deposition onto ice
  753. ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
  754. ! - MASSI - integrated quantity associated w/ ice mass
  755. ! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
  756. ! precipitation rates
  757. !
  758. !
  759. !-------------------------------------------------------------------------
  760. !
  761. !--- VEL_RF - velocity increase of rimed particles as functions of crude
  762. ! particle size categories (at 0.1 mm intervals of mean ice particle
  763. ! sizes) and rime factor (different values of Rime Factor of 1.1**N,
  764. ! where N=0 to Nrime).
  765. !
  766. !-------------------------------------------------------------------------
  767. !
  768. !--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
  769. ! (0.45 mm), assuming an exponential size distribution.
  770. !
  771. !-------------------------------------------------------------------------
  772. !------- Key parameters, local variables, & important comments ---------
  773. !-----------------------------------------------------------------------
  774. !
  775. !--- TOLER => Tolerance or precision for accumulated precipitation
  776. !
  777. REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, Xratio=.025
  778. !
  779. !--- If BLEND=1:
  780. ! precipitation (large) ice amounts are estimated at each level as a
  781. ! blend of ice falling from the grid point above and the precip ice
  782. ! present at the start of the time step (see TOT_ICE below).
  783. !--- If BLEND=0:
  784. ! precipitation (large) ice amounts are estimated to be the precip
  785. ! ice present at the start of the time step.
  786. !
  787. !--- Extended to include sedimentation of rain on 2/5/01
  788. !
  789. REAL, PARAMETER :: BLEND=1.
  790. !
  791. !-----------------------------------------------------------------------
  792. !--- Local variables
  793. !-----------------------------------------------------------------------
  794. !
  795. REAL EMAIRI, N0r, NLICE, NSmICE, RHgrd !GFDL
  796. LOGICAL CLEAR, ICE_logical, DBG_logical, RAIN_logical
  797. INTEGER :: IDR,INDEX_MY,INDEXR,INDEXR1,INDEXS,IPASS,ITDX,IXRF, &
  798. & IXS,LBEF,L
  799. !
  800. REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
  801. & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
  802. & DENOMI,DENOMW,DENOMWI,DIDEP, &
  803. & DIEVP,DIFFUS,DLI,DTPH,DTRHO,DUM,DUM1, &
  804. & DUM2,DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLARGE,FLIMASS, &
  805. & FSMALL,FWR,FWS,GAMMAR,GAMMAS, &
  806. & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
  807. & PIEVP,PILOSS,PIMLT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
  808. & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
  809. & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,QV,QW,QW0,QWnew, &
  810. & RFACTOR,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
  811. & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
  812. & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
  813. & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
  814. & WC,WCnew,WSgrd,WS,WSnew,WV,WVnew,WVQW, &
  815. & XLF,XLF1,XLI,XLV,XLV1,XLV2,XLIMASS,XRF,XSIMASS
  816. !
  817. !#######################################################################
  818. !########################## Begin Execution ############################
  819. !#######################################################################
  820. !
  821. !
  822. ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
  823. ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
  824. !
  825. !-----------------------------------------------------------------------
  826. !------------ Loop from top (L=1) to surface (L=LSFC) ------------------
  827. !-----------------------------------------------------------------------
  828. !
  829. DO 10 L=1,LSFC
  830. !--- Skip this level and go to the next lower level if no condensate
  831. ! and very low specific humidities
  832. !
  833. IF (QV_col(L).LE.EPSQ .AND. WC_col(L).LE.EPSQ) GO TO 10
  834. !
  835. !-----------------------------------------------------------------------
  836. !------------ Proceed with cloud microphysics calculations -------------
  837. !-----------------------------------------------------------------------
  838. !
  839. TK=T_col(L) ! Temperature (deg K)
  840. TC=TK-T0C ! Temperature (deg C)
  841. PP=P_col(L) ! Pressure (Pa)
  842. QV=QV_col(L) ! Specific humidity of water vapor (kg/kg)
  843. WV=QV/(1.-QV) ! Water vapor mixing ratio (kg/kg)
  844. WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
  845. RHgrd=RHC_col(L) ! Threshold relative humidity for the onset of condensation
  846. !
  847. !-----------------------------------------------------------------------
  848. !--- Moisture variables below are mixing ratios & not specifc humidities
  849. !-----------------------------------------------------------------------
  850. !
  851. CLEAR=.TRUE.
  852. !
  853. !--- This check is to determine grid-scale saturation when no condensate is present
  854. !
  855. ESW=1000.*FPVS0(TK) ! Saturation vapor pressure w/r/t water
  856. QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
  857. QSI = QSW ! Initialize variable
  858. WS=QSW ! General saturation mixing ratio (water/ice)
  859. IF (TC .LT. 0.) THEN
  860. ESI=1000.*FPVS(TK) ! Saturation vapor pressure w/r/t ice
  861. QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
  862. WS=QSI ! General saturation mixing ratio (water/ice)
  863. ENDIF
  864. !
  865. !--- Effective grid-scale Saturation mixing ratios
  866. !
  867. QSWgrd=RHgrd*QSW
  868. QSIgrd=RHgrd*QSI
  869. WSgrd=RHgrd*WS
  870. !
  871. !--- Check if air is subsaturated and w/o condensate
  872. !
  873. IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
  874. !
  875. !--- Check if any rain is falling into layer from above
  876. !
  877. IF (ARAIN .GT. CLIMIT) THEN
  878. CLEAR=.FALSE.
  879. ELSE
  880. ARAIN=0.
  881. ENDIF
  882. !
  883. !--- Check if any ice is falling into layer from above
  884. !
  885. !--- NOTE that "SNOW" in variable names is synonomous with
  886. ! large, precipitation ice particles
  887. !
  888. IF (ASNOW .GT. CLIMIT) THEN
  889. CLEAR=.FALSE.
  890. ELSE
  891. ASNOW=0.
  892. ENDIF
  893. !
  894. !-----------------------------------------------------------------------
  895. !-- Loop to the end if in clear, subsaturated air free of condensate ---
  896. !-----------------------------------------------------------------------
  897. !
  898. IF (CLEAR) GO TO 10
  899. !
  900. !-----------------------------------------------------------------------
  901. !--------- Initialize RHO, THICK & microphysical processes -------------
  902. !-----------------------------------------------------------------------
  903. !
  904. !
  905. !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
  906. ! (see pp. 63-65 in Fleagle & Businger, 1963)
  907. !
  908. RHO=PP/(RD*TK*(1.+EPS1*QV)) ! Air density (kg/m**3)
  909. RRHO=1./RHO ! Reciprocal of air density
  910. DTRHO=DTPH*RHO ! Time step * air density
  911. BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density
  912. THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
  913. !
  914. ARAINnew=0. ! Updated accumulated rainfall
  915. ASNOWnew=0. ! Updated accumulated snowfall
  916. QI=QI_col(L) ! Ice mixing ratio
  917. QInew=0. ! Updated ice mixing ratio
  918. QR=QR_col(L) ! Rain mixing ratio
  919. QRnew=0. ! Updated rain ratio
  920. QW=QW_col(L) ! Cloud water mixing ratio
  921. QWnew=0. ! Updated cloud water ratio
  922. !
  923. PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
  924. PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
  925. PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
  926. PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
  927. PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
  928. PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
  929. PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
  930. PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
  931. PIMLT=0. ! Melting ice (kg/kg; >0)
  932. PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
  933. PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
  934. PREVP=0. ! Rain evaporation (kg/kg; <0)
  935. !
  936. !--- Double check input hydrometeor mixing ratios
  937. !
  938. ! DUM=WC-(QI+QW+QR)
  939. ! DUM1=ABS(DUM)
  940. ! DUM2=TOLER*MIN(WC, QI+QW+QR)
  941. ! IF (DUM1 .GT. DUM2) THEN
  942. ! WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
  943. ! & ' L=',L
  944. ! WRITE(6,"(4(a12,g11.4,1x))")
  945. ! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
  946. ! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
  947. ! ENDIF
  948. !
  949. !***********************************************************************
  950. !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
  951. !***********************************************************************
  952. !
  953. !--- Calculate a few variables, which are used more than once below
  954. !
  955. !--- Latent heat of vaporization as a function of temperature from
  956. ! Bolton (1980, JAS)
  957. !
  958. XLV=3.148E6-2370*TK ! Latent heat of vaporization (Lv)
  959. XLF=XLS-XLV ! Latent heat of fusion (Lf)
  960. XLV1=XLV*RCP ! Lv/Cp
  961. XLF1=XLF*RCP ! Lf/Cp
  962. TK2=1./(TK*TK) ! 1./TK**2
  963. !GFDL XLV2=XLV*XLV*QSW*TK2/RV ! Lv**2*Qsw/(Rv*TK**2)
  964. XLV2=XLV*XLV*QSWgrd*TK2/RV ! Lv**2*QSWgrd/(Rv*TK**2) !GFDL
  965. DENOMW=1.+XLV2*RCP ! Denominator term, Clausius-Clapeyron correction
  966. !
  967. !--- Basic thermodynamic quantities
  968. ! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
  969. ! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
  970. ! * DIFFUS - diffusivity of water vapor [ m**2/s ]
  971. !
  972. TFACTOR=TK**1.5/(TK+120.)
  973. DYNVIS=1.496E-6*TFACTOR
  974. THERM_COND=2.116E-3*TFACTOR
  975. DIFFUS=8.794E-5*TK**1.81/PP
  976. !
  977. !--- Air resistance term for the fall speed of ice following the
  978. ! basic research by Heymsfield, Kajikawa, others
  979. !
  980. GAMMAS=(1.E5/PP)**C1
  981. !
  982. !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
  983. !
  984. GAMMAR=(RHO0/RHO)**.4
  985. !
  986. !----------------------------------------------------------------------
  987. !------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
  988. !----------------------------------------------------------------------
  989. !
  990. !--- Determine if conditions supporting ice are present
  991. !
  992. IF (TC.LT.0. .OR. QI.GT.EPSQ .OR. ASNOW.GT.CLIMIT) THEN
  993. ICE_logical=.TRUE.
  994. ELSE
  995. ICE_logical=.FALSE.
  996. QLICE=0.
  997. QTICE=0.
  998. ENDIF
  999. !
  1000. !--- Determine if rain is present
  1001. !
  1002. RAIN_logical=.FALSE.
  1003. IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
  1004. !
  1005. IF (ICE_logical) THEN
  1006. !
  1007. !--- IMPORTANT: Estimate time-averaged properties.
  1008. !
  1009. !---
  1010. ! * FLARGE - ratio of number of large ice to total (large & small) ice
  1011. ! * FSMALL - ratio of number of small ice crystals to large ice particles
  1012. ! -> Small ice particles are assumed to have a mean diameter of 50 microns.
  1013. ! * XSIMASS - used for calculating small ice mixing ratio
  1014. !---
  1015. ! * TOT_ICE - total mass (small & large) ice before microphysics,
  1016. ! which is the sum of the total mass of large ice in the
  1017. ! current layer and the input flux of ice from above
  1018. ! * PILOSS - greatest loss (<0) of total (small & large) ice by
  1019. ! sublimation, removing all of the ice falling from above
  1020. ! and the ice within the layer
  1021. ! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
  1022. ! ice mass to the unrimed ice mass (>=1)
  1023. ! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
  1024. ! * VSNOW - Fall speed of rimed snow w/ air resistance correction
  1025. ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
  1026. ! * XLIMASS - used for calculating large ice mixing ratio
  1027. ! * FLIMASS - mass fraction of large ice
  1028. ! * QTICE - time-averaged mixing ratio of total ice
  1029. ! * QLICE - time-averaged mixing ratio of large ice
  1030. ! * NLICE - time-averaged number concentration of large ice
  1031. ! * NSmICE - number concentration of small ice crystals at current level
  1032. !---
  1033. !--- Assumed number fraction of large ice particles to total (large & small)
  1034. ! ice particles, which is based on a general impression of the literature.
  1035. !
  1036. WVQW=WV+QW ! Water vapor & cloud water
  1037. !
  1038. IF (TC.GE.0. .OR. WVQW.LT.QSIgrd) THEN
  1039. !
  1040. !--- Eliminate small ice particle contributions for melting & sublimation
  1041. !
  1042. FLARGE=FLARGE1
  1043. ELSE
  1044. !
  1045. !--- Enhanced number of small ice particles during depositional growth
  1046. ! (effective only when 0C > T >= T_ice [-10C] )
  1047. !
  1048. FLARGE=FLARGE2
  1049. !
  1050. !--- Larger number of small ice particles due to rime splintering
  1051. !
  1052. IF (TC.GE.-8. .AND. TC.LE.-3.) FLARGE=.5*FLARGE
  1053. !
  1054. ENDIF ! End IF (TC.GE.0. .OR. WVQW.LT.QSIgrd)
  1055. !GFDL => turned on in GFDL code, but not here => FLARGE=1.0
  1056. FSMALL=(1.-FLARGE)/FLARGE
  1057. XSIMASS=RRHO*MASSI(MDImin)*FSMALL
  1058. IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT) THEN
  1059. INDEXS=MDImin
  1060. TOT_ICE=0.
  1061. PILOSS=0.
  1062. RimeF1=1.
  1063. VrimeF=1.
  1064. VEL_INC=GAMMAS
  1065. VSNOW=0.
  1066. EMAIRI=THICK
  1067. XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
  1068. FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
  1069. QLICE=0.
  1070. QTICE=0.
  1071. NLICE=0.
  1072. NSmICE=0.
  1073. ELSE
  1074. !
  1075. !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
  1076. ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
  1077. ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
  1078. !
  1079. DUM=XMImax*EXP(.0536*TC)
  1080. INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
  1081. TOT_ICE=THICK*QI+BLEND*ASNOW
  1082. PILOSS=-TOT_ICE/THICK
  1083. LBEF=MAX(1,L-1)
  1084. DUM1=RimeF_col(LBEF)
  1085. DUM2=RimeF_col(L)
  1086. RimeF1=(DUM2*THICK*QI+DUM1*BLEND*ASNOW)/TOT_ICE
  1087. RimeF1=MIN(RimeF1, RFmax)
  1088. DO IPASS=0,1
  1089. IF (RimeF1 .LE. 1.) THEN
  1090. RimeF1=1.
  1091. VrimeF=1.
  1092. ELSE
  1093. IXS=MAX(2, MIN(INDEXS/100, 9))
  1094. XRF=10.492*ALOG(RimeF1)
  1095. IXRF=MAX(0, MIN(INT(XRF), Nrime))
  1096. IF (IXRF .GE. Nrime) THEN
  1097. Vri

Large files files are truncated, but you can click here to view the full file