/wrfv2_fire/phys/module_mp_etanew.F
FORTRAN Legacy | 2675 lines | 1273 code | 68 blank | 1334 comment | 1 complexity | 8ae618e33a201896565cf54034d178fc MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !WRF:MODEL_MP:PHYSICS
- !
- MODULE module_mp_etanew
- !
- !-----------------------------------------------------------------------
- REAL,PRIVATE,SAVE :: ABFR, CBFR, CIACW, CIACR, C_N0r0, &
- & CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
- & RFmax, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
- & RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DR4, &
- & RR_DR5, RR_DRmax
- !
- INTEGER, PRIVATE,PARAMETER :: MY_T1=1, MY_T2=35
- REAL,PRIVATE,DIMENSION(MY_T1:MY_T2),SAVE :: MY_GROWTH
- !
- REAL, PRIVATE,PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, &
- & DelDMI=1.e-6,XMImin=1.e6*DMImin
- INTEGER, PUBLIC,PARAMETER :: XMImax=1.e6*DMImax, XMIexp=.0536, &
- & MDImin=XMImin, MDImax=XMImax
- REAL, PRIVATE,DIMENSION(MDImin:MDImax) :: &
- & ACCRI,SDENS,VSNOWI,VENTI1,VENTI2
- !
- REAL, PRIVATE,PARAMETER :: DMRmin=.05e-3, DMRmax=1.e-3, &
- & DelDMR=1.e-6,XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax
- INTEGER, PRIVATE,PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax
- REAL, PRIVATE,DIMENSION(MDRmin:MDRmax):: &
- & ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
- !
- INTEGER, PRIVATE,PARAMETER :: Nrime=40
- REAL, DIMENSION(2:9,0:Nrime),PRIVATE,SAVE :: VEL_RF
- !
- INTEGER,PARAMETER :: NX=7501
- REAL, PARAMETER :: XMIN=180.0,XMAX=330.0
- REAL, DIMENSION(NX),SAVE :: TBPVS,TBPVS0
- REAL, SAVE :: C1XPVS0,C2XPVS0,C1XPVS,C2XPVS
- !
- REAL, PRIVATE,PARAMETER :: &
- !--- Physical constants follow:
- & CP=1004.6, EPSQ=1.E-12, GRAV=9.806, RHOL=1000., RD=287.04 &
- & ,RV=461.5, T0C=273.15, XLS=2.834E6 &
- !--- Derived physical constants follow:
- & ,EPS=RD/RV, EPS1=RV/RD-1., EPSQ1=1.001*EPSQ &
- & ,RCP=1./CP, RCPRV=RCP/RV, RGRAV=1./GRAV, RRHOL=1./RHOL &
- & ,XLS1=XLS*RCP, XLS2=XLS*XLS*RCPRV, XLS3=XLS*XLS/RV &
- !--- Constants specific to the parameterization follow:
- !--- CLIMIT/CLIMIT1 are lower limits for treating accumulated precipitation
- & ,CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT &
- & ,C1=1./3. &
- & ,DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, DMR4=0.45E-3 &
- & ,DMR5=0.67E-3 &
- & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3 &
- & ,XMR4=1.e6*DMR4, XMR5=1.e6*DMR5
- !
- INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3, MDR4=XMR4 &
- & , MDR5=XMR5
- !
- ! ======================================================================
- !--- Important tunable parameters that are exported to other modules
- ! * RHgrd - threshold relative humidity for onset of condensation
- ! * T_ICE - temperature (C) threshold at which all remaining liquid water
- ! is glaciated to ice
- ! * T_ICE_init - maximum temperature (C) at which ice nucleation occurs
- ! * NLImax - maximum number concentrations (m**-3) of large ice (snow/graupel/sleet)
- ! * NLImin - minimum number concentrations (m**-3) of large ice (snow/graupel/sleet)
- ! * N0r0 - assumed intercept (m**-4) of rain drops if drop diameters are between 0.2 and 1.0 mm
- ! * N0rmin - minimum intercept (m**-4) for rain drops
- ! * NCW - number concentrations of cloud droplets (m**-3)
- ! * FLARGE1, FLARGE2 - number fraction of large ice to total (large+snow) ice
- ! at T>0C and in presence of sublimation (FLARGE1), otherwise in
- ! presence of ice saturated/supersaturated conditions
- ! ======================================================================
- REAL, PUBLIC,PARAMETER :: &
- & RHgrd=1. &
- & ,T_ICE=-40. &
- & ,T_ICEK=T0C+T_ICE &
- & ,T_ICE_init=-5. &
- & ,NLImax=5.E3 &
- & ,NLImin=1.E3 &
- & ,N0r0=8.E6 &
- & ,N0rmin=1.E4 &
- & ,NCW=250.E6 &
- & ,FLARGE1=1. &
- & ,FLARGE2=.2
- !--- Other public variables passed to other routines:
- REAL,PUBLIC,SAVE :: QAUT0
- REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
- !
- !
- CONTAINS
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
- & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, &
- & LOWLYR,SR, &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- & QC,QR,QS, &
- & mp_restart_state,tbpvs_state,tbpvs0_state, &
- & RAINNC,RAINNCV, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
- & ,IMS,IME,JMS,JME,KMS,KME &
- & ,ITS,ITE,JTS,JTE,KTS,KTE &
- & ,ITIMESTEP
- REAL, INTENT(IN) :: DT,DX,DY
- REAL, INTENT(IN), DIMENSION(ims:ime, kms:kme, jms:jme):: &
- & dz8w,p_phy,pi_phy,rho_phy
- REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme):: &
- & th_phy,qv,qt
- REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
- & qc,qr,qs
- REAL, INTENT(INOUT), DIMENSION(ims:ime, kms:kme, jms:jme ) :: &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY
- REAL, INTENT(INOUT), DIMENSION(ims:ime,jms:jme) :: &
- & RAINNC,RAINNCV
- REAL, INTENT(OUT), DIMENSION(ims:ime,jms:jme):: SR
- !
- REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
- !
- REAL,DIMENSION(nx),INTENT(INOUT) :: TBPVS_STATE,TBPVS0_STATE
- !
- INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
- !-----------------------------------------------------------------------
- ! LOCAL VARS
- !-----------------------------------------------------------------------
- ! NSTATS,QMAX,QTOT are diagnostic vars
- INTEGER,DIMENSION(ITLO:ITHI,4) :: NSTATS
- REAL, DIMENSION(ITLO:ITHI,5) :: QMAX
- REAL, DIMENSION(ITLO:ITHI,22):: QTOT
- ! SOME VARS WILL BE USED FOR DATA ASSIMILATION (DON'T NEED THEM NOW).
- ! THEY ARE TREATED AS LOCAL VARS, BUT WILL BECOME STATE VARS IN THE
- ! FUTURE. SO, WE DECLARED THEM AS MEMORY SIZES FOR THE FUTURE USE
- ! TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related
- ! the microphysics scheme. Instead, they will be used by Eta precip
- ! assimilation.
- REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
- & TLATGS_PHY,TRAIN_PHY
- REAL, DIMENSION(ims:ime,jms:jme):: APREC,PREC,ACPREC
- REAL, DIMENSION(its:ite, kts:kte, jts:jte):: t_phy
- INTEGER :: I,J,K,KFLIP
- REAL :: WC
- !
- !-----------------------------------------------------------------------
- !**********************************************************************
- !-----------------------------------------------------------------------
- !
- MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
- !
- C1XPVS0=MP_RESTART_STATE(MY_T2+1)
- C2XPVS0=MP_RESTART_STATE(MY_T2+2)
- C1XPVS =MP_RESTART_STATE(MY_T2+3)
- C2XPVS =MP_RESTART_STATE(MY_T2+4)
- CIACW =MP_RESTART_STATE(MY_T2+5)
- CIACR =MP_RESTART_STATE(MY_T2+6)
- CRACW =MP_RESTART_STATE(MY_T2+7)
- CRAUT =MP_RESTART_STATE(MY_T2+8)
- !
- TBPVS(1:NX) =TBPVS_STATE(1:NX)
- TBPVS0(1:NX)=TBPVS0_STATE(1:NX)
- !
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
- qv(i,k,j)=qv(i,k,j)/(1.+qv(i,k,j)) !Convert to specific humidity
- ENDDO
- ENDDO
- ENDDO
- ! initial diagnostic variables and data assimilation vars
- ! (will need to delete this part in the future)
- DO k = 1,4
- DO i = ITLO,ITHI
- NSTATS(i,k)=0.
- ENDDO
- ENDDO
- DO k = 1,5
- DO i = ITLO,ITHI
- QMAX(i,k)=0.
- ENDDO
- ENDDO
- DO k = 1,22
- DO i = ITLO,ITHI
- QTOT(i,k)=0.
- ENDDO
- ENDDO
- ! initial data assimilation vars (will need to delete this part in the future)
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- TLATGS_PHY (i,k,j)=0.
- TRAIN_PHY (i,k,j)=0.
- ENDDO
- ENDDO
- ENDDO
- DO j = jts,jte
- DO i = its,ite
- ACPREC(i,j)=0.
- APREC (i,j)=0.
- PREC (i,j)=0.
- SR (i,j)=0.
- ENDDO
- ENDDO
- !-- NOTE: ARW QT has been advected, while QR, QS and QC have not
- !
- !-- Update QT, F_ice, F_rain arrays for WRF NMM only
- #if (NMM_CORE==1)
- !
- !-- NOTE: The total ice array in this code is "QS" because the vast
- ! majority of the ice mass is in the form of snow, and using
- ! the "QS" array should result in better coupling with the
- ! Dudhia SW package. NMM calls microphysics after other
- ! physics, so use updated QR, QS and QC to update QT array.
- !
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QS(I,K,J)
- IF (QS(I,K,J) <= EPSQ) THEN
- F_ICE_PHY(I,K,J)=0.
- IF (T_PHY(I,K,J) < T_ICEK) F_ICE_PHY(I,K,J)=1.
- ELSE
- F_ICE_PHY(I,K,J)=MAX( 0., MIN(1., QS(I,K,J)/QT(I,K,J) ) )
- ENDIF
- IF (QR(I,K,J) <= EPSQ) THEN
- F_RAIN_PHY(I,K,J)=0.
- ELSE
- F_RAIN_PHY(I,K,J)=QR(I,K,J)/(QC(I,K,J)+QR(I,K,J))
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- #endif
- !-----------------------------------------------------------------------
- CALL EGCP01DRV(DT,LOWLYR, &
- & APREC,PREC,ACPREC,SR,NSTATS,QMAX,QTOT, &
- & dz8w,rho_phy,qt,t_phy,qv,F_ICE_PHY,P_PHY, &
- & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte )
- !-----------------------------------------------------------------------
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- th_phy(i,k,j) = t_phy(i,k,j)/pi_phy(i,k,j)
- qv(i,k,j)=qv(i,k,j)/(1.-qv(i,k,j)) !Convert to mixing ratio
- WC=qt(I,K,J)
- QS(I,K,J)=0.
- QR(I,K,J)=0.
- QC(I,K,J)=0.
- IF(F_ICE_PHY(I,K,J)>=1.)THEN
- QS(I,K,J)=WC
- ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
- QC(I,K,J)=WC
- ELSE
- QS(I,K,J)=F_ICE_PHY(I,K,J)*WC
- QC(I,K,J)=WC-QS(I,K,J)
- ENDIF
- !
- IF(QC(I,K,J)>0..AND.F_RAIN_PHY(I,K,J)>0.)THEN
- IF(F_RAIN_PHY(I,K,J).GE.1.)THEN
- QR(I,K,J)=QC(I,K,J)
- QC(I,K,J)=0.
- ELSE
- QR(I,K,J)=F_RAIN_PHY(I,K,J)*QC(I,K,J)
- QC(I,K,J)=QC(I,K,J)-QR(I,K,J)
- ENDIF
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- !
- ! update rain (from m to mm)
- DO j=jts,jte
- DO i=its,ite
- RAINNC(i,j)=APREC(i,j)*1000.+RAINNC(i,j)
- RAINNCV(i,j)=APREC(i,j)*1000.
- ENDDO
- ENDDO
- !
- MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
- MP_RESTART_STATE(MY_T2+1)=C1XPVS0
- MP_RESTART_STATE(MY_T2+2)=C2XPVS0
- MP_RESTART_STATE(MY_T2+3)=C1XPVS
- MP_RESTART_STATE(MY_T2+4)=C2XPVS
- MP_RESTART_STATE(MY_T2+5)=CIACW
- MP_RESTART_STATE(MY_T2+6)=CIACR
- MP_RESTART_STATE(MY_T2+7)=CRACW
- MP_RESTART_STATE(MY_T2+8)=CRAUT
- !
- TBPVS_STATE(1:NX) =TBPVS(1:NX)
- TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
- !-----------------------------------------------------------------------
- END SUBROUTINE ETAMP_NEW
- !-----------------------------------------------------------------------
- SUBROUTINE EGCP01DRV( &
- & DTPH,LOWLYR,APREC,PREC,ACPREC,SR, &
- & NSTATS,QMAX,QTOT, &
- & dz8w,RHO_PHY,CWM_PHY,T_PHY,Q_PHY,F_ICE_PHY,P_PHY, &
- & F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY,TRAIN_PHY, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte)
- !-----------------------------------------------------------------------
- ! DTPH Physics time step (s)
- ! CWM_PHY (qt) Mixing ratio of the total condensate. kg/kg
- ! Q_PHY Mixing ratio of water vapor. kg/kg
- ! F_RAIN_PHY Fraction of rain.
- ! F_ICE_PHY Fraction of ice.
- ! F_RIMEF_PHY Mass ratio of rimed ice (rime factor).
- !
- !TLATGS_PHY,TRAIN_PHY,APREC,PREC,ACPREC,SR are not directly related the
- !micrphysics sechme. Instead, they will be used by Eta precip assimilation.
- !
- !NSTATS,QMAX,QTOT are used for diagnosis purposes.
- !
- !-----------------------------------------------------------------------
- !--- Variables APREC,PREC,ACPREC,SR are calculated for precip assimilation
- ! and/or ZHAO's scheme in Eta and are not required by this microphysics
- ! scheme itself.
- !--- NSTATS,QMAX,QTOT are used for diagnosis purposes only. They will be
- ! printed out when PRINT_diag is true.
- !
- !-----------------------------------------------------------------------
- IMPLICIT NONE
- !-----------------------------------------------------------------------
- !
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
- ! VARIABLES PASSED IN/OUT
- INTEGER,INTENT(IN ) :: ids,ide, jds,jde, kds,kde &
- & ,ims,ime, jms,jme, kms,kme &
- & ,its,ite, jts,jte, kts,kte
- REAL,INTENT(IN) :: DTPH
- INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(INOUT) :: LOWLYR
- INTEGER,DIMENSION(ITLO:ITHI,4),INTENT(INOUT) :: NSTATS
- REAL,DIMENSION(ITLO:ITHI,5),INTENT(INOUT) :: QMAX
- REAL,DIMENSION(ITLO:ITHI,22),INTENT(INOUT) :: QTOT
- REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
- & APREC,PREC,ACPREC,SR
- REAL,DIMENSION( its:ite, kts:kte, jts:jte ),INTENT(INOUT) :: t_phy
- REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(IN) :: &
- & dz8w,P_PHY,RHO_PHY
- REAL,DIMENSION( ims:ime, kms:kme, jms:jme ),INTENT(INOUT) :: &
- & CWM_PHY, F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY,TLATGS_PHY &
- & ,Q_PHY,TRAIN_PHY
- !
- !-----------------------------------------------------------------------
- !LOCAL VARIABLES
- !-----------------------------------------------------------------------
- !
- #define CACHE_FRIENDLY_MP_ETANEW
- #ifdef CACHE_FRIENDLY_MP_ETANEW
- # define TEMP_DIMS kts:kte,its:ite,jts:jte
- # define TEMP_DEX L,I,J
- #else
- # define TEMP_DIMS its:ite,jts:jte,kts:kte
- # define TEMP_DEX I,J,L
- #endif
- !
- INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
- REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
- REAL,DIMENSION(kts:kte,its:ite,jts:jte) :: F_ice,F_rain,F_RimeF
- INTEGER,DIMENSION(its:ite,jts:jte) :: LMH
- REAL :: TC,WC,QI,QR,QW,Fice,Frain,DUM,ASNOW,ARAIN
- REAL,DIMENSION(kts:kte) :: P_col,Q_col,T_col,QV_col,WC_col, &
- RimeF_col,QI_col,QR_col,QW_col, THICK_col,DPCOL
- REAL,DIMENSION(2) :: PRECtot,PRECmax
- !-----------------------------------------------------------------------
- !
- DO J=JTS,JTE
- DO I=ITS,ITE
- LMH(I,J) = KTE-LOWLYR(I,J)+1
- ENDDO
- ENDDO
- DO 98 J=JTS,JTE
- DO 98 I=ITS,ITE
- DO L=KTS,KTE
- KFLIP=KTE+1-L
- CWM(TEMP_DEX)=CWM_PHY(I,KFLIP,J)
- T(TEMP_DEX)=T_PHY(I,KFLIP,J)
- Q(TEMP_DEX)=Q_PHY(I,KFLIP,J)
- P(TEMP_DEX)=P_PHY(I,KFLIP,J)
- TLATGS(TEMP_DEX)=TLATGS_PHY(I,KFLIP,J)
- TRAIN(TEMP_DEX)=TRAIN_PHY(I,KFLIP,J)
- F_ice(L,I,J)=F_ice_PHY(I,KFLIP,J)
- F_rain(L,I,J)=F_rain_PHY(I,KFLIP,J)
- F_RimeF(L,I,J)=F_RimeF_PHY(I,KFLIP,J)
- ENDDO
- 98 CONTINUE
-
- DO 100 J=JTS,JTE
- DO 100 I=ITS,ITE
- LSFC=LMH(I,J) ! "L" of surface
- !
- DO K=KTS,KTE
- KFLIP=KTE+1-K
- DPCOL(K)=RHO_PHY(I,KFLIP,J)*GRAV*dz8w(I,KFLIP,J)
- ENDDO
- !
- !
- !--- Initialize column data (1D arrays)
- !
- L=1
- IF (CWM(TEMP_DEX) .LE. EPSQ) CWM(TEMP_DEX)=EPSQ
- F_ice(1,I,J)=1.
- F_rain(1,I,J)=0.
- F_RimeF(1,I,J)=1.
- DO L=1,LSFC
- !
- !--- Pressure (Pa) = (Psfc-Ptop)*(ETA/ETA_sfc)+Ptop
- !
- P_col(L)=P(TEMP_DEX)
- !
- !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
- !
- THICK_col(L)=DPCOL(L)*RGRAV
- T_col(L)=T(TEMP_DEX)
- TC=T_col(L)-T0C
- QV_col(L)=max(EPSQ, Q(TEMP_DEX))
- IF (CWM(TEMP_DEX) .LE. EPSQ1) THEN
- WC_col(L)=0.
- IF (TC .LT. T_ICE) THEN
- F_ice(L,I,J)=1.
- ELSE
- F_ice(L,I,J)=0.
- ENDIF
- F_rain(L,I,J)=0.
- F_RimeF(L,I,J)=1.
- ELSE
- WC_col(L)=CWM(TEMP_DEX)
- ENDIF
- !
- !--- Determine composition of condensate in terms of
- ! cloud water, ice, & rain
- !
- WC=WC_col(L)
- QI=0.
- QR=0.
- QW=0.
- Fice=F_ice(L,I,J)
- Frain=F_rain(L,I,J)
- IF (Fice .GE. 1.) THEN
- QI=WC
- ELSE IF (Fice .LE. 0.) THEN
- QW=WC
- ELSE
- QI=Fice*WC
- QW=WC-QI
- ENDIF
- IF (QW.GT.0. .AND. Frain.GT.0.) THEN
- IF (Frain .GE. 1.) THEN
- QR=QW
- QW=0.
- ELSE
- QR=Frain*QW
- QW=QW-QR
- ENDIF
- ENDIF
- IF (QI .LE. 0.) F_RimeF(L,I,J)=1.
- RimeF_col(L)=F_RimeF(L,I,J) ! (real)
- QI_col(L)=QI
- QR_col(L)=QR
- QW_col(L)=QW
- ENDDO
- !
- !#######################################################################
- !
- !--- Perform the microphysical calculations in this column
- !
- I_index=I
- J_index=J
- CALL EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
- & P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
- & THICK_col, WC_col,KTS,KTE,NSTATS,QMAX,QTOT )
- !
- !#######################################################################
- !
- !
- !--- Update storage arrays
- !
- DO L=1,LSFC
- TRAIN(TEMP_DEX)=(T_col(L)-T(TEMP_DEX))/DTPH
- TLATGS(TEMP_DEX)=T_col(L)-T(TEMP_DEX)
- T(TEMP_DEX)=T_col(L)
- Q(TEMP_DEX)=QV_col(L)
- CWM(TEMP_DEX)=WC_col(L)
- !
- !--- REAL*4 array storage
- !
- IF (QI_col(L) .LE. EPSQ) THEN
- F_ice(L,I,J)=0.
- IF (T_col(L) .LT. T_ICEK) F_ice(L,I,J)=1.
- F_RimeF(L,I,J)=1.
- ELSE
- F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_col(L)) )
- F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
- ENDIF
- IF (QR_col(L) .LE. EPSQ) THEN
- DUM=0
- ELSE
- DUM=QR_col(L)/(QR_col(L)+QW_col(L))
- ENDIF
- F_rain(L,I,J)=DUM
- !
- ENDDO
- !
- !--- Update accumulated precipitation statistics
- !
- !--- Surface precipitation statistics; SR is fraction of surface
- ! precipitation (if >0) associated with snow
- !
- APREC(I,J)=(ARAIN+ASNOW)*RRHOL ! Accumulated surface precip (depth in m) !<--- Ying
- PREC(I,J)=PREC(I,J)+APREC(I,J)
- ACPREC(I,J)=ACPREC(I,J)+APREC(I,J)
- IF(APREC(I,J) .LT. 1.E-8) THEN
- SR(I,J)=0.
- ELSE
- SR(I,J)=RRHOL*ASNOW/APREC(I,J)
- ENDIF
- !
- !--- Debug statistics
- !
- IF (PRINT_diag) THEN
- PRECtot(1)=PRECtot(1)+ARAIN
- PRECtot(2)=PRECtot(2)+ASNOW
- PRECmax(1)=MAX(PRECmax(1), ARAIN)
- PRECmax(2)=MAX(PRECmax(2), ASNOW)
- ENDIF
- !#######################################################################
- !#######################################################################
- !
- 100 CONTINUE ! End "I" & "J" loops
- DO 101 J=JTS,JTE
- DO 101 I=ITS,ITE
- DO L=KTS,KTE
- KFLIP=KTE+1-L
- CWM_PHY(I,KFLIP,J)=CWM(TEMP_DEX)
- T_PHY(I,KFLIP,J)=T(TEMP_DEX)
- Q_PHY(I,KFLIP,J)=Q(TEMP_DEX)
- TLATGS_PHY(I,KFLIP,J)=TLATGS(TEMP_DEX)
- TRAIN_PHY(I,KFLIP,J)=TRAIN(TEMP_DEX)
- F_ice_PHY(I,KFLIP,J)=F_ice(L,I,J)
- F_rain_PHY(I,KFLIP,J)=F_rain(L,I,J)
- F_RimeF_PHY(I,KFLIP,J)=F_RimeF(L,I,J)
- ENDDO
- 101 CONTINUE
- END SUBROUTINE EGCP01DRV
- !
- !
- !###############################################################################
- ! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
- ! (1) Represents sedimentation by preserving a portion of the precipitation
- ! through top-down integration from cloud-top. Modified procedure to
- ! Zhao and Carr (1997).
- ! (2) Microphysical equations are modified to be less sensitive to time
- ! steps by use of Clausius-Clapeyron equation to account for changes in
- ! saturation mixing ratios in response to latent heating/cooling.
- ! (3) Prevent spurious temperature oscillations across 0C due to
- ! microphysics.
- ! (4) Uses lookup tables for: calculating two different ventilation
- ! coefficients in condensation and deposition processes; accretion of
- ! cloud water by precipitation; precipitation mass; precipitation rate
- ! (and mass-weighted precipitation fall speeds).
- ! (5) Assumes temperature-dependent variation in mean diameter of large ice
- ! (Houze et al., 1979; Ryan et al., 1996).
- ! -> 8/22/01: This relationship has been extended to colder temperatures
- ! to parameterize smaller large-ice particles down to mean sizes of MDImin,
- ! which is 50 microns reached at -55.9C.
- ! (6) Attempts to differentiate growth of large and small ice, mainly for
- ! improved transition from thin cirrus to thick, precipitating ice
- ! anvils.
- ! -> 8/22/01: This feature has been diminished by effectively adjusting to
- ! ice saturation during depositional growth at temperatures colder than
- ! -10C. Ice sublimation is calculated more explicitly. The logic is
- ! that sources of are either poorly understood (e.g., nucleation for NWP)
- ! or are not represented in the Eta model (e.g., detrainment of ice from
- ! convection). Otherwise the model is too wet compared to the radiosonde
- ! observations based on 1 Feb - 18 March 2001 retrospective runs.
- ! (7) Top-down integration also attempts to treat mixed-phase processes,
- ! allowing a mixture of ice and water. Based on numerous observational
- ! studies, ice growth is based on nucleation at cloud top &
- ! subsequent growth by vapor deposition and riming as the ice particles
- ! fall through the cloud. Effective nucleation rates are a function
- ! of ice supersaturation following Meyers et al. (JAM, 1992).
- ! -> 8/22/01: The simulated relative humidities were far too moist compared
- ! to the rawinsonde observations. This feature has been substantially
- ! diminished, limited to a much narrower temperature range of 0 to -10C.
- ! (8) Depositional growth of newly nucleated ice is calculated for large time
- ! steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
- ! using their ice crystal masses calculated after 600 s of growth in water
- ! saturated conditions. The growth rates are normalized by time step
- ! assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
- ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
- ! (9) Ice precipitation rates can increase due to increase in response to
- ! cloud water riming due to (a) increased density & mass of the rimed
- ! ice, and (b) increased fall speeds of rimed ice.
- ! -> 8/22/01: This feature has been effectively limited to 0 to -10C.
- !###############################################################################
- !###############################################################################
- !
- SUBROUTINE EGCP01COLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, &
- & LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
- & THICK_col, WC_col ,KTS,KTE,NSTATS,QMAX,QTOT)
- !
- !###############################################################################
- !###############################################################################
- !
- !-------------------------------------------------------------------------------
- !----- NOTE: Code is currently set up w/o threading!
- !-------------------------------------------------------------------------------
- !$$$ SUBPROGRAM DOCUMENTATION BLOCK
- ! . . .
- ! SUBPROGRAM: Grid-scale microphysical processes - condensation & precipitation
- ! PRGRMMR: Ferrier ORG: W/NP22 DATE: 08-2001
- ! PRGRMMR: Jin (Modification for WRF structure)
- !-------------------------------------------------------------------------------
- ! ABSTRACT:
- ! * Merges original GSCOND & PRECPD subroutines.
- ! * Code has been substantially streamlined and restructured.
- ! * Exchange between water vapor & small cloud condensate is calculated using
- ! the original Asai (1965, J. Japan) algorithm. See also references to
- ! Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
- ! (1989, MWR). This algorithm replaces the Sundqvist et al. (1989, MWR)
- ! parameterization.
- !-------------------------------------------------------------------------------
- !
- ! USAGE:
- ! * CALL EGCP01COLUMN FROM SUBROUTINE EGCP01DRV
- !
- ! INPUT ARGUMENT LIST:
- ! DTPH - physics time step (s)
- ! I_index - I index
- ! J_index - J index
- ! LSFC - Eta level of level above surface, ground
- ! P_col - vertical column of model pressure (Pa)
- ! QI_col - vertical column of model ice mixing ratio (kg/kg)
- ! QR_col - vertical column of model rain ratio (kg/kg)
- ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
- ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
- ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
- ! T_col - vertical column of model temperature (deg K)
- ! THICK_col - vertical column of model mass thickness (density*height increment)
- ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
- !
- !
- ! OUTPUT ARGUMENT LIST:
- ! ARAIN - accumulated rainfall at the surface (kg)
- ! ASNOW - accumulated snowfall at the surface (kg)
- ! QV_col - vertical column of model water vapor specific humidity (kg/kg)
- ! WC_col - vertical column of model mixing ratio of total condensate (kg/kg)
- ! QW_col - vertical column of model cloud water mixing ratio (kg/kg)
- ! QI_col - vertical column of model ice mixing ratio (kg/kg)
- ! QR_col - vertical column of model rain ratio (kg/kg)
- ! RimeF_col - vertical column of rime factor for ice in model (ratio, defined below)
- ! T_col - vertical column of model temperature (deg K)
- !
- ! OUTPUT FILES:
- ! NONE
- !
- ! Subprograms & Functions called:
- ! * Real Function CONDENSE - cloud water condensation
- ! * Real Function DEPOSIT - ice deposition (not sublimation)
- !
- ! UNIQUE: NONE
- !
- ! LIBRARY: NONE
- !
- ! COMMON BLOCKS:
- ! CMICRO_CONS - key constants initialized in GSMCONST
- ! CMICRO_STATS - accumulated and maximum statistics
- ! CMY_GROWTH - lookup table for growth of ice crystals in
- ! water saturated conditions (Miller & Young, 1979)
- ! IVENT_TABLES - lookup tables for ventilation effects of ice
- ! IACCR_TABLES - lookup tables for accretion rates of ice
- ! IMASS_TABLES - lookup tables for mass content of ice
- ! IRATE_TABLES - lookup tables for precipitation rates of ice
- ! IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
- ! RVENT_TABLES - lookup tables for ventilation effects of rain
- ! RACCR_TABLES - lookup tables for accretion rates of rain
- ! RMASS_TABLES - lookup tables for mass content of rain
- ! RVELR_TABLES - lookup tables for fall speeds of rain
- ! RRATE_TABLES - lookup tables for precipitation rates of rain
- !
- ! ATTRIBUTES:
- ! LANGUAGE: FORTRAN 90
- ! MACHINE : IBM SP
- !
- !
- !-------------------------------------------------------------------------
- !--------------- Arrays & constants in argument list ---------------------
- !-------------------------------------------------------------------------
- !
- IMPLICIT NONE
- !
- INTEGER,INTENT(IN) :: KTS,KTE,I_index, J_index, LSFC
- REAL,INTENT(INOUT) :: ARAIN, ASNOW
- REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: P_col, QI_col,QR_col &
- & ,QV_col ,QW_col, RimeF_col, T_col, THICK_col,WC_col
- !
- !-------------------------------------------------------------------------
- !-------------- Common blocks for microphysical statistics ---------------
- !-------------------------------------------------------------------------
- !
- !-------------------------------------------------------------------------
- !--------- Common blocks for constants initialized in GSMCONST ----------
- !
- INTEGER, PARAMETER :: ITLO=-60, ITHI=40
- INTEGER,INTENT(INOUT) :: NSTATS(ITLO:ITHI,4)
- REAL,INTENT(INOUT) :: QMAX(ITLO:ITHI,5),QTOT(ITLO:ITHI,22)
- !
- !-------------------------------------------------------------------------
- !--------------- Common blocks for various lookup tables -----------------
- !
- !--- Discretized growth rates of small ice crystals after their nucleation
- ! at 1 C intervals from -1 C to -35 C, based on calculations by Miller
- ! and Young (1979, JAS) after 600 s of growth. Resultant growth rates
- ! are multiplied by physics time step in GSMCONST.
- !
- !-------------------------------------------------------------------------
- !
- !--- Mean ice-particle diameters varying from 50 microns to 1000 microns
- ! (1 mm), assuming an exponential size distribution.
- !
- !---- Meaning of the following arrays:
- ! - mdiam - mean diameter (m)
- ! - VENTI1 - integrated quantity associated w/ ventilation effects
- ! (capacitance only) for calculating vapor deposition onto ice
- ! - VENTI2 - integrated quantity associated w/ ventilation effects
- ! (with fall speed) for calculating vapor deposition onto ice
- ! - ACCRI - integrated quantity associated w/ cloud water collection by ice
- ! - MASSI - integrated quantity associated w/ ice mass
- ! - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
- ! precipitation rates
- !
- !
- !-------------------------------------------------------------------------
- !
- !--- VEL_RF - velocity increase of rimed particles as functions of crude
- ! particle size categories (at 0.1 mm intervals of mean ice particle
- ! sizes) and rime factor (different values of Rime Factor of 1.1**N,
- ! where N=0 to Nrime).
- !
- !-------------------------------------------------------------------------
- !
- !--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 1000 microns
- ! (1 mm), assuming an exponential size distribution.
- !
- !-------------------------------------------------------------------------
- !------- Key parameters, local variables, & important comments ---------
- !-----------------------------------------------------------------------
- !
- !--- TOLER => Tolerance or precision for accumulated precipitation
- !
- REAL, PARAMETER :: TOLER=5.E-7, C2=1./6., RHO0=1.194, Xratio=.025
- !
- !--- If BLEND=1:
- ! precipitation (large) ice amounts are estimated at each level as a
- ! blend of ice falling from the grid point above and the precip ice
- ! present at the start of the time step (see TOT_ICE below).
- !--- If BLEND=0:
- ! precipitation (large) ice amounts are estimated to be the precip
- ! ice present at the start of the time step.
- !
- !--- Extended to include sedimentation of rain on 2/5/01
- !
- REAL, PARAMETER :: BLEND=1.
- !
- !--- This variable is for debugging purposes (if .true.)
- !
- LOGICAL, PARAMETER :: PRINT_diag=.FALSE.
- !
- !-----------------------------------------------------------------------
- !--- Local variables
- !-----------------------------------------------------------------------
- !
- REAL EMAIRI, N0r, NLICE, NSmICE
- LOGICAL CLEAR, ICE_logical, DBG_logical, RAIN_logical
- INTEGER :: IDR,INDEX_MY,INDEXR,INDEXR1,INDEXS,IPASS,ITDX,IXRF, &
- & IXS,LBEF,L
- !
- REAL :: ABI,ABW,AIEVP,ARAINnew,ASNOWnew,BLDTRH,BUDGET, &
- & CREVP,DELI,DELR,DELT,DELV,DELW,DENOMF, &
- & DENOMI,DENOMW,DENOMWI,DIDEP, &
- & DIEVP,DIFFUS,DLI,DTPH,DTRHO,DUM,DUM1, &
- & DUM2,DWV0,DWVI,DWVR,DYNVIS,ESI,ESW,FIR,FLARGE,FLIMASS, &
- & FSMALL,FWR,FWS,GAMMAR,GAMMAS, &
- & PCOND,PIACR,PIACW,PIACWI,PIACWR,PICND,PIDEP,PIDEP_max, &
- & PIEVP,PILOSS,PIMLT,PP,PRACW,PRAUT,PREVP,PRLOSS, &
- & QI,QInew,QLICE,QR,QRnew,QSI,QSIgrd,QSInew,QSW,QSW0, &
- & QSWgrd,QSWnew,QT,QTICE,QTnew,QTRAIN,QV,QW,QW0,QWnew, &
- & RFACTOR,RHO,RIMEF,RIMEF1,RQR,RR,RRHO,SFACTOR, &
- & TC,TCC,TFACTOR,THERM_COND,THICK,TK,TK2,TNEW, &
- & TOT_ICE,TOT_ICEnew,TOT_RAIN,TOT_RAINnew, &
- & VEL_INC,VENTR,VENTIL,VENTIS,VRAIN1,VRAIN2,VRIMEF,VSNOW, &
- & WC,WCnew,WSgrd,WS,WSnew,WV,WVnew,WVQW, &
- & XLF,XLF1,XLI,XLV,XLV1,XLV2,XLIMASS,XRF,XSIMASS
- !
- !#######################################################################
- !########################## Begin Execution ############################
- !#######################################################################
- !
- !
- ARAIN=0. ! Accumulated rainfall into grid box from above (kg/m**2)
- ASNOW=0. ! Accumulated snowfall into grid box from above (kg/m**2)
- !
- !-----------------------------------------------------------------------
- !------------ Loop from top (L=1) to surface (L=LSFC) ------------------
- !-----------------------------------------------------------------------
- !
- DO 10 L=1,LSFC
- !--- Skip this level and go to the next lower level if no condensate
- ! and very low specific humidities
- !
- IF (QV_col(L).LE.EPSQ .AND. WC_col(L).LE.EPSQ) GO TO 10
- !
- !-----------------------------------------------------------------------
- !------------ Proceed with cloud microphysics calculations -------------
- !-----------------------------------------------------------------------
- !
- TK=T_col(L) ! Temperature (deg K)
- TC=TK-T0C ! Temperature (deg C)
- PP=P_col(L) ! Pressure (Pa)
- QV=QV_col(L) ! Specific humidity of water vapor (kg/kg)
- WV=QV/(1.-QV) ! Water vapor mixing ratio (kg/kg)
- WC=WC_col(L) ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
- !
- !-----------------------------------------------------------------------
- !--- Moisture variables below are mixing ratios & not specifc humidities
- !-----------------------------------------------------------------------
- !
- CLEAR=.TRUE.
- !
- !--- This check is to determine grid-scale saturation when no condensate is present
- !
- ESW=MIN(1000.*FPVS0(TK),0.99*PP) ! Saturation vapor pressure w/r/t water
- QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
- WS=QSW ! General saturation mixing ratio (water/ice)
- IF (TC .LT. 0.) THEN
- ESI=MIN(1000.*FPVS(TK),0.99*PP) ! Saturation vapor pressure w/r/t ice
- QSI=EPS*ESI/(PP-ESI) ! Saturation mixing ratio w/r/t water
- WS=QSI ! General saturation mixing ratio (water/ice)
- ENDIF
- !
- !--- Effective grid-scale Saturation mixing ratios
- !
- QSWgrd=RHgrd*QSW
- QSIgrd=RHgrd*QSI
- WSgrd=RHgrd*WS
- !
- !--- Check if air is subsaturated and w/o condensate
- !
- IF (WV.GT.WSgrd .OR. WC.GT.EPSQ) CLEAR=.FALSE.
- !
- !--- Check if any rain is falling into layer from above
- !
- IF (ARAIN .GT. CLIMIT) THEN
- CLEAR=.FALSE.
- ELSE
- ARAIN=0.
- ENDIF
- !
- !--- Check if any ice is falling into layer from above
- !
- !--- NOTE that "SNOW" in variable names is synonomous with
- ! large, precipitation ice particles
- !
- IF (ASNOW .GT. CLIMIT) THEN
- CLEAR=.FALSE.
- ELSE
- ASNOW=0.
- ENDIF
- !
- !-----------------------------------------------------------------------
- !-- Loop to the end if in clear, subsaturated air free of condensate ---
- !-----------------------------------------------------------------------
- !
- IF (CLEAR) GO TO 10
- !
- !-----------------------------------------------------------------------
- !--------- Initialize RHO, THICK & microphysical processes -------------
- !-----------------------------------------------------------------------
- !
- !
- !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
- ! (see pp. 63-65 in Fleagle & Businger, 1963)
- !
- RHO=PP/(RD*TK*(1.+EPS1*QV)) ! Air density (kg/m**3)
- RRHO=1./RHO ! Reciprocal of air density
- DTRHO=DTPH*RHO ! Time step * air density
- BLDTRH=BLEND*DTRHO ! Blend parameter * time step * air density
- THICK=THICK_col(L) ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
- !
- ARAINnew=0. ! Updated accumulated rainfall
- ASNOWnew=0. ! Updated accumulated snowfall
- QI=QI_col(L) ! Ice mixing ratio
- QInew=0. ! Updated ice mixing ratio
- QR=QR_col(L) ! Rain mixing ratio
- QRnew=0. ! Updated rain ratio
- QW=QW_col(L) ! Cloud water mixing ratio
- QWnew=0. ! Updated cloud water ratio
- !
- PCOND=0. ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
- PIDEP=0. ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
- PIACW=0. ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
- PIACWI=0. ! Growth of precip ice by riming (kg/kg; >0)
- PIACWR=0. ! Shedding of accreted cloud water to form rain (kg/kg; >0)
- PIACR=0. ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
- PICND=0. ! Condensation (>0) onto wet, melting ice (kg/kg)
- PIEVP=0. ! Evaporation (<0) from wet, melting ice (kg/kg)
- PIMLT=0. ! Melting ice (kg/kg; >0)
- PRAUT=0. ! Cloud water autoconversion to rain (kg/kg; >0)
- PRACW=0. ! Cloud water collection (accretion) by rain (kg/kg; >0)
- PREVP=0. ! Rain evaporation (kg/kg; <0)
- !
- !--- Double check input hydrometeor mixing ratios
- !
- ! DUM=WC-(QI+QW+QR)
- ! DUM1=ABS(DUM)
- ! DUM2=TOLER*MIN(WC, QI+QW+QR)
- ! IF (DUM1 .GT. DUM2) THEN
- ! WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
- ! & ' L=',L
- ! WRITE(6,"(4(a12,g11.4,1x))")
- ! & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
- ! & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
- ! ENDIF
- !
- !***********************************************************************
- !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
- !***********************************************************************
- !
- !--- Calculate a few variables, which are used more than once below
- !
- !--- Latent heat of vaporization as a function of temperature from
- ! Bolton (1980, JAS)
- !
- XLV=3.148E6-2370*TK ! Latent heat of vaporization (Lv)
- XLF=XLS-XLV ! Latent heat of fusion (Lf)
- XLV1=XLV*RCP ! Lv/Cp
- XLF1=XLF*RCP ! Lf/Cp
- TK2=1./(TK*TK) ! 1./TK**2
- XLV2=XLV*XLV*QSW*TK2/RV ! Lv**2*Qsw/(Rv*TK**2)
- DENOMW=1.+XLV2*RCP ! Denominator term, Clausius-Clapeyron correction
- !
- !--- Basic thermodynamic quantities
- ! * DYNVIS - dynamic viscosity [ kg/(m*s) ]
- ! * THERM_COND - thermal conductivity [ J/(m*s*K) ]
- ! * DIFFUS - diffusivity of water vapor [ m**2/s ]
- !
- TFACTOR=TK**1.5/(TK+120.)
- DYNVIS=1.496E-6*TFACTOR
- THERM_COND=2.116E-3*TFACTOR
- DIFFUS=8.794E-5*TK**1.81/PP
- !
- !--- Air resistance term for the fall speed of ice following the
- ! basic research by Heymsfield, Kajikawa, others
- !
- GAMMAS=(1.E5/PP)**C1
- !
- !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
- !
- GAMMAR=(RHO0/RHO)**.4
- !
- !----------------------------------------------------------------------
- !------------- IMPORTANT MICROPHYSICS DECISION TREE -----------------
- !----------------------------------------------------------------------
- !
- !--- Determine if conditions supporting ice are present
- !
- IF (TC.LT.0. .OR. QI.GT. EPSQ .OR. ASNOW.GT.CLIMIT) THEN
- ICE_logical=.TRUE.
- ELSE
- ICE_logical=.FALSE.
- QLICE=0.
- QTICE=0.
- ENDIF
- !
- !--- Determine if rain is present
- !
- RAIN_logical=.FALSE.
- IF (ARAIN.GT.CLIMIT .OR. QR.GT.EPSQ) RAIN_logical=.TRUE.
- !
- IF (ICE_logical) THEN
- !
- !--- IMPORTANT: Estimate time-averaged properties.
- !
- !---
- ! * FLARGE - ratio of number of large ice to total (large & small) ice
- ! * FSMALL - ratio of number of small ice crystals to large ice particles
- ! -> Small ice particles are assumed to have a mean diameter of 50 microns.
- ! * XSIMASS - used for calculating small ice mixing ratio
- !---
- ! * TOT_ICE - total mass (small & large) ice before microphysics,
- ! which is the sum of the total mass of large ice in the
- ! current layer and the input flux of ice from above
- ! * PILOSS - greatest loss (<0) of total (small & large) ice by
- ! sublimation, removing all of the ice falling from above
- ! and the ice within the layer
- ! * RimeF1 - Rime Factor, which is the mass ratio of total (unrimed & rimed)
- ! ice mass to the unrimed ice mass (>=1)
- ! * VrimeF - the velocity increase due to rime factor or melting (ratio, >=1)
- ! * VSNOW - Fall speed of rimed snow w/ air resistance correction
- ! * EMAIRI - equivalent mass of air associated layer and with fall of snow into layer
- ! * XLIMASS - used for calculating large ice mixing ratio
- ! * FLIMASS - mass fraction of large ice
- ! * QTICE - time-averaged mixing ratio of total ice
- ! * QLICE - time-averaged mixing ratio of large ice
- ! * NLICE - time-averaged number concentration of large ice
- ! * NSmICE - number concentration of small ice crystals at current level
- !---
- !--- Assumed number fraction of large ice particles to total (large & small)
- ! ice particles, which is based on a general impression of the literature.
- !
- WVQW=WV+QW ! Water vapor & cloud water
- !
- IF (TC.GE.0. .OR. WVQW.LT.QSIgrd) THEN
- !
- !--- Eliminate small ice particle contributions for melting & sublimation
- !
- FLARGE=FLARGE1
- ELSE
- !
- !--- Enhanced number of small ice particles during depositional growth
- ! (effective only when 0C > T >= T_ice [-10C] )
- !
- FLARGE=FLARGE2
- !
- !--- Larger number of small ice particles due to rime splintering
- !
- IF (TC.GE.-8. .AND. TC.LE.-3.) FLARGE=.5*FLARGE
- !
- ENDIF ! End IF (TC.GE.0. .OR. WVQW.LT.QSIgrd)
- FSMALL=(1.-FLARGE)/FLARGE
- XSIMASS=RRHO*MASSI(MDImin)*FSMALL
- IF (QI.LE.EPSQ .AND. ASNOW.LE.CLIMIT) THEN
- INDEXS=MDImin
- TOT_ICE=0.
- PILOSS=0.
- RimeF1=1.
- VrimeF=1.
- VEL_INC=GAMMAS
- VSNOW=0.
- EMAIRI=THICK
- XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
- FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
- QLICE=0.
- QTICE=0.
- NLICE=0.
- NSmICE=0.
- ELSE
- !
- !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
- ! converted from Fig. 5 plot of LAMDAs. Similar set of relationships
- ! also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
- !
- DUM=XMImax*EXP(.0536*TC)
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
- TOT_ICE=THICK*QI+BLEND*ASNOW
- PILOSS=-TOT_ICE/THICK
- LBEF=MAX(1,L-1)
- DUM1=RimeF_col(LBEF)
- DUM2=RimeF_col(L)
- RimeF1=(DUM2*THICK*QI+DUM1*BLEND*ASNOW)/TOT_ICE
- RimeF1=MIN(RimeF1, RFmax)
- DO IPASS=0,1
- IF (RimeF1 .LE. 1.) THEN
- RimeF1=1.
- VrimeF=1.
- ELSE
- IXS=MAX(2, MIN(INDEXS/100, 9))
- XRF=10.492*ALOG(RimeF1)
- IXRF=MAX(0, MIN(INT(XRF), Nrime))
- IF (IXRF .GE. Nrime) THEN
- VrimeF=VEL_RF(IXS,Nrime)
- ELSE
- VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
- & (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
- ENDIF
- ENDIF ! End IF (RimeF1 .LE. 1.)
- VEL_INC=GAMMAS*VrimeF
- VSNOW=VEL_INC*VSNOWI(INDEXS)
- EMAIRI=THICK+BLDTRH*VSNOW
- XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
- FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
- QTICE=TOT_ICE/EMAIRI
- QLICE=FLIMASS*QTICE
- NLICE=QLICE/XLIMASS
- NSmICE=Fsmall*NLICE
- !
- IF ( (NLICE.GE.NLImin .AND. NLICE.LE.NLImax) &
- & .OR. IPASS.EQ.1) THEN
- EXIT
- ELSE
- IF (TC < 0) THEN
- XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
- IF (XLI .LE. MASSI(MDImin) ) THEN
- INDEXS=MDImin
- ELSE IF (XLI .LE. MASSI(450) ) THEN
- DLI=9.5885E5*XLI**.42066 ! DLI in microns
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
- ELSE IF (XLI .LE. MASSI(MDImax) ) THEN
- DLI=3.9751E6*XLI**.49870 ! DLI in microns
- INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
- ELSE
- INDEXS=MDImax
- ENDIF ! End IF (XLI .LE. MASSI(MDImin) )
- ENDIF ! End IF (TC < 0)
- !
- !--- Reduce excessive accumulation of ice at upper levels
- ! associated with strong grid-resolved ascent
- !
- !--- Force NLICE to be between NLImin and NLImax
- !
- !
- !--- 8/22/01: Increase density of large ice if maximum limits
- ! are reached for number concentration (NLImax) and mean size
- ! (MDImax). Done to increase fall out of ice.
- !
- DUM=MAX(NLImin, MIN(NLImax, NLICE) )
- IF (DUM.GE.NLImax .AND. INDEXS.GE.MDImax) &
- & RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
- ! WRITE(6,"(4(a12,g11.4,1x))")
- ! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
- ! & '{…
Large files files are truncated, but you can click here to view the full file