[Dart-dev] [3442] DART/trunk/models/PBL_1d:
Updates from Josh Hacker to bring the PBL source up to date.
nancy at ucar.edu
nancy at ucar.edu
Wed Jul 2 07:48:03 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080702/8d861000/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/models/PBL_1d/src/module_mp_etanew.F90
===================================================================
--- DART/trunk/models/PBL_1d/src/module_mp_etanew.F90 2008-07-01 23:43:29 UTC (rev 3441)
+++ DART/trunk/models/PBL_1d/src/module_mp_etanew.F90 2008-07-02 13:48:02 UTC (rev 3442)
@@ -153,2441 +153,8 @@
!-----------------------------------------------------------------------
!**********************************************************************
!-----------------------------------------------------------------------
-!
- 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
-
-!-----------------------------------------------------------------------
-
- 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
+ STOP 'ETANEW NOT IMPLEMENTED'
!
- 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
- 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
- !
- 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(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 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.
-!
-!--- 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=1000.*FPVS0(TK) ! 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=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
- 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
- !
- !--- Reduce excessive accumulation of ice at upper levels
- ! associated with strong grid-resolved ascent
- !
- !--- Force NLICE to be between NLImin and NLImax
- !
- DUM=MAX(NLImin, MIN(NLImax, NLICE) )
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list