/wrfv2_fire/phys/module_mp_HWRF.F
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
- !WRF:MODEL_MP:PHYSICS
- !
- MODULE module_mp_HWRF
- !-----------------------------------------------------------------------
- !-- The following changes were made on 24 July 2006.
- ! (1) All known version 2.1 dependencies were removed from the
- ! operational WRF NMM model code (search for "!HWRF")
- ! (2) Incorporated code changes from the GFDL model (search for "!GFDL")
- !-----------------------------------------------------------------------
- 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_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=.45e-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),PRIVATE,SAVE :: TBPVS,TBPVS0
- REAL, PRIVATE,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 &
- & ,XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, XMR3=1.e6*DMR3
- INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3
- REAL:: WC
- !
- ! ======================================================================
- !--- Important tunable parameters that are exported to other modules
- !GFDL * RHgrd - generic reference to the threshold relative humidity for
- !GFDL the onset of condensation
- !GFDL (new) * RHgrd_in - "RHgrd" for the inner domain
- !GFDL (new) * RHgrd_out - "RHgrd" for the outer domain
- !HWRF 6/11/2010 mod - use lower RHgrd_out for p < 850 hPa
- ! * 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 0.45 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
- ! * PRINT_diag - for extended model diagnostics (code currently commented out)
- ! * PRINT_err - for run-time prints when water budgets are not conserved (for debugging)
- ! ======================================================================
- REAL, PUBLIC,PARAMETER :: &
- ! & RHgrd=1. &
- & RHgrd_in=1. & !GFDL
- & ,RHgrd_out=0.975 & !GFDL
- & ,P_RHgrd_out=850.E2 & !HWRF 6/11/2010
- & ,T_ICE=-40. & !GFDL
- !GFDL & ,T_ICE=-30. &
- & ,T_ICEK=T0C+T_ICE &
- & ,T_ICE_init=-5. &
- & ,NLImax=5.E3 &
- & ,NLImin=1.E3 &
- & ,N0r0=8.E6 &
- & ,N0rmin=1.E4 &
- & ,NCW=60.E6 & !GFDL
- !HWRF & ,NCW=100.E6 &
- & ,FLARGE1=1. &
- & ,FLARGE2=.2
- ! LOGICAL, PARAMETER :: PRINT_diag=.FALSE. !GFDL
- LOGICAL, PARAMETER :: PRINT_err=.TRUE. !GFDL** => eventually set this to .false.
- !--- Other public variables passed to other routines:
- REAL,PUBLIC,SAVE :: QAUT0
- REAL, PUBLIC,DIMENSION(MDImin:MDImax) :: MASSI
- !
- !
- CONTAINS
- !-----------------------------------------------------------------------
- !-----------------------------------------------------------------------
- SUBROUTINE ETAMP_NEW_HWRF (itimestep,DT,DX,DY,GID,RAINNC,RAINNCV, & !GID
- & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qt, & !gopal's doing
- & LOWLYR,SR, &
- & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- & QC,QR,QI, &
- & ids,ide, jds,jde, kds,kde, &
- & ims,ime, jms,jme, kms,kme, &
- & its,ite, jts,jte, kts,kte )
- !HWRF SUBROUTINE ETAMP_NEW (itimestep,DT,DX,DY, &
- !HWRF & dz8w,rho_phy,p_phy,pi_phy,th_phy,qv,qc, &
- !HWRF & LOWLYR,SR, &
- !HWRF & F_ICE_PHY,F_RAIN_PHY,F_RIMEF_PHY, &
- !HWRF & mp_restart_state,tbpvs_state,tbpvs0_state, &
- !HWRF & RAINNC,RAINNCV, &
- !HWRF & ids,ide, jds,jde, kds,kde, &
- !HWRF & ims,ime, jms,jme, kms,kme, &
- !HWRF & 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,GID ! GID gopal's doing
- 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,qc,qr,qi
- 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
- !
- !HWRF REAL,DIMENSION(*),INTENT(INOUT) :: MP_RESTART_STATE
- !
- !HWRF 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
- !
- !-----------------------------------------------------------------------
- !**********************************************************************
- !-----------------------------------------------------------------------
- !
- !HWRF MY_GROWTH(MY_T1:MY_T2)=MP_RESTART_STATE(MY_T1:MY_T2)
- !HWRF!
- !HWRF C1XPVS0=MP_RESTART_STATE(MY_T2+1)
- !HWRF C2XPVS0=MP_RESTART_STATE(MY_T2+2)
- !HWRF C1XPVS =MP_RESTART_STATE(MY_T2+3)
- !HWRF C2XPVS =MP_RESTART_STATE(MY_T2+4)
- !HWRF CIACW =MP_RESTART_STATE(MY_T2+5)
- !HWRF CIACR =MP_RESTART_STATE(MY_T2+6)
- !HWRF CRACW =MP_RESTART_STATE(MY_T2+7)
- !HWRF CRAUT =MP_RESTART_STATE(MY_T2+8)
- !HWRF!
- !HWRF TBPVS(1:NX) =TBPVS_STATE(1:NX)
- !HWRF 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
- !-- 6/11/2010: Update QT, F_ice, F_rain arrays
- DO j = jts,jte
- DO k = kts,kte
- DO i = its,ite
- QT(I,K,J)=QC(I,K,J)+QR(I,K,J)+QI(I,K,J)
- IF (QI(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., QI(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)/(QR(I,K,J)+QC(I,K,J))
- ENDIF
- ENDDO
- ENDDO
- ENDDO
- !-----------------------------------------------------------------------
- CALL EGCP01DRV(GID,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)
- QI(I,K,J)=0.
- QR(I,K,J)=0.
- QC(I,K,J)=0.
- IF(F_ICE_PHY(I,K,J)>=1.)THEN
- QI(I,K,J)=WC
- ELSEIF(F_ICE_PHY(I,K,J)<=0.)THEN
- QC(I,K,J)=WC
- ELSE
- QI(I,K,J)=F_ICE_PHY(I,K,J)*WC
- QC(I,K,J)=WC-QI(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
- !
- !HWRF MP_RESTART_STATE(MY_T1:MY_T2)=MY_GROWTH(MY_T1:MY_T2)
- !HWRF MP_RESTART_STATE(MY_T2+1)=C1XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+2)=C2XPVS0
- !HWRF MP_RESTART_STATE(MY_T2+3)=C1XPVS
- !HWRF MP_RESTART_STATE(MY_T2+4)=C2XPVS
- !HWRF MP_RESTART_STATE(MY_T2+5)=CIACW
- !HWRF MP_RESTART_STATE(MY_T2+6)=CIACR
- !HWRF MP_RESTART_STATE(MY_T2+7)=CRACW
- !HWRF MP_RESTART_STATE(MY_T2+8)=CRAUT
- !HWRF!
- !HWRF TBPVS_STATE(1:NX) =TBPVS(1:NX)
- !HWRF TBPVS0_STATE(1:NX)=TBPVS0(1:NX)
- !-----------------------------------------------------------------------
- END SUBROUTINE ETAMP_NEW_HWRF
- !-----------------------------------------------------------------------
- SUBROUTINE EGCP01DRV(GID, & !GID gopal's doing
- & 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
- ! 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
- INTEGER,INTENT(IN ) :: GID ! grid%id gopal's doing
- 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
- !-----------------------------------------------------------------------
- !
- !HWRF - Below are directives in the operational code that have been removed,
- ! where "TEMP_DEX" has been replaced with "I,J,L" and "TEMP_DIMS" has
- ! been replaced with "its:ite,jts:jte,kts:kte"
- !HWRF#define CACHE_FRIENDLY_MP_ETANEW
- !HWRF#ifdef CACHE_FRIENDLY_MP_ETANEW
- !HWRF# define TEMP_DIMS kts:kte,its:ite,jts:jte
- !HWRF# define TEMP_DEX L,I,J
- !HWRF#else
- !HWRF# define TEMP_DIMS its:ite,jts:jte,kts:kte
- !HWRF# define TEMP_DEX I,J,L
- !HWRF#endif
- !HWRF!
- INTEGER :: LSFC,I,J,I_index,J_index,L,K,KFLIP
- !HWRF REAL,DIMENSION(TEMP_DIMS) :: CWM,T,Q,TRAIN,TLATGS,P
- REAL,DIMENSION(its:ite,jts:jte,kts:kte) :: &
- & 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, RHC_col, DPCOL !GFDL
- 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(I,J,L)=CWM_PHY(I,KFLIP,J)
- T(I,J,L)=T_PHY(I,KFLIP,J)
- Q(I,J,L)=Q_PHY(I,KFLIP,J)
- P(I,J,L)=P_PHY(I,KFLIP,J)
- TLATGS(I,J,L)=TLATGS_PHY(I,KFLIP,J)
- TRAIN(I,J,L)=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)
- !
- IF (CWM(I,J,1) .LE. EPSQ) CWM(I,J,1)=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(I,J,L)
- !
- !--- Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
- !
- THICK_col(L)=DPCOL(L)*RGRAV
- T_col(L)=T(I,J,L)
- TC=T_col(L)-T0C
- QV_col(L)=max(EPSQ, Q(I,J,L))
- IF (CWM(I,J,L) .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(I,J,L)
- 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
- RimeF_col(L)=F_RimeF(L,I,J)
- QI_col(L)=QI
- QR_col(L)=QR
- QW_col(L)=QW
- !GFDL => New. Added RHC_col to allow for height- and grid-dependent values for
- !GFDL the relative humidity threshold for condensation ("RHgrd")
- !6/11/2010 mod - Use lower RHgrd_out threshold for < 850 hPa
- !------------------------------------------------------------
- IF(GID .EQ. 1 .AND. P_col(L)<P_RHgrd_out) THEN ! gopal's doing based on GFDL
- RHC_col(L)=RHgrd_out
- ELSE
- RHC_col(L)=RHgrd_in
- ENDIF
- !------------------------------------------------------------
- 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, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT ) !GFDL
- !
- !#######################################################################
- !
- !
- !--- Update storage arrays
- !
- DO L=1,LSFC
- TRAIN(I,J,L)=(T_col(L)-T(I,J,L))/DTPH
- TLATGS(I,J,L)=T_col(L)-T(I,J,L)
- T(I,J,L)=T_col(L)
- Q(I,J,L)=QV_col(L)
- CWM(I,J,L)=WC_col(L)
- !
- !--- REAL*4 array storage
- !
- F_RimeF(L,I,J)=MAX(1., RimeF_col(L))
- 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.
- ELSE
- F_ice(L,I,J)=MAX( 0., MIN(1., QI_col(L)/WC_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(I,J,L)
- T_PHY(I,KFLIP,J)=T(I,J,L)
- Q_PHY(I,KFLIP,J)=Q(I,J,L)
- TLATGS_PHY(I,KFLIP,J)=TLATGS(I,J,L)
- TRAIN_PHY(I,KFLIP,J)=TRAIN(I,J,L)
- 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, RHC_col, KTS,KTE,NSTATS,QMAX,QTOT) !GFDL
- !
- !###############################################################################
- !###############################################################################
- !
- !-------------------------------------------------------------------------------
- !----- 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)
- ! RHC_col - vertical column of threshold relative humidity for onset of condensation (ratio) !GFDL
- !
- !
- ! 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, RHC_col !GFDL
- !
- !-------------------------------------------------------------------------
- !-------------- 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 450 microns
- ! (0.45 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.
- !
- !-----------------------------------------------------------------------
- !--- Local variables
- !-----------------------------------------------------------------------
- !
- REAL EMAIRI, N0r, NLICE, NSmICE, RHgrd !GFDL
- 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)
- RHgrd=RHC_col(L) ! Threshold relative humidity for the onset of condensation
- !
- !-----------------------------------------------------------------------
- !--- 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=1000.*FPVS0(TK) ! Saturation vapor pressure w/r/t water
- QSW=EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water
- QSI = QSW ! Initialize variable
- WS=QSW ! General saturation mixing ratio (water/ice)
- IF (TC .LT. 0.) THEN
- ESI=1000.*FPVS(TK) ! 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
- !GFDL XLV2=XLV*XLV*QSW*TK2/RV ! Lv**2*Qsw/(Rv*TK**2)
- XLV2=XLV*XLV*QSWgrd*TK2/RV ! Lv**2*QSWgrd/(Rv*TK**2) !GFDL
- 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)
- !GFDL => turned on in GFDL code, but not here => FLARGE=1.0
- 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
- Vri…
Large files files are truncated, but you can click here to view the full file