!############################################################################### ! ***** 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 GSMCOLUMN ( 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 ) !############################################################################### !############################################################################### !------------------------------------------------------------------------------- !----- 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 !------------------------------------------------------------------------------- ! 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 GSMCOLUMN FROM SUBROUTINE GSMDRIVE ! * SUBROUTINE GSMDRIVE CALLED FROM MAIN PROGRAM EBU ! 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 --------------------- !------------------------------------------------------------------------- INCLUDE "parmeta.f90" INCLUDE "mpp.h" REAL :: ARAIN, ASNOW, P_col(LM), QI_col(LM), QR_col(LM), QV_col(LM), & QW_col(LM), RimeF_col(LM), T_col(LM), THICK_col(LM), & WC_col(LM) INTEGER :: I_index, J_index, LSFC !------------------------------------------------------------------------- !-------------- Common blocks for microphysical statistics --------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !--------- Common blocks for constants initialized in GSMCONST ---------- COMMON /CMICRO_CONS/ ABFR, CBFR, CIACW, CIACR, C_N0r0, & CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, & QAUT0, RFmax, RHgrd, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, & RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax !--- The following variables are for microphysical statistics INTEGER, PARAMETER :: ITLO=-60, ITHI=40 COMMON /CMICRO_STATS/ NSTATS(ITLO:ITHI,4), QMAX(ITLO:ITHI,5), & QTOT(ITLO:ITHI,22) INTEGER :: NSTATS REAL :: QMAX, QTOT !------------------------------------------------------------------------- !--------------- 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. INTEGER, PARAMETER :: MY_T1=1, MY_T2=35 COMMON /CMY600/ MY_GROWTH(MY_T1:MY_T2) REAL :: MY_GROWTH !------------------------------------------------------------------------- !--- 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 REAL, PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, DelDMI=1.e-6, & XMImin=1.e6*DMImin, XMImax=1.e6*DMImax INTEGER, PARAMETER :: MDImin=XMImin, MDImax=XMImax COMMON /IACCR_TABLES/ ACCRI(MDImin:MDImax) COMMON /IMASS_TABLES/ MASSI(MDImin:MDImax) REAL :: MASSI COMMON /IRATE_TABLES/ VSNOWI(MDImin:MDImax) COMMON /IVENT_TABLES/ VENTI1(MDImin:MDImax), VENTI2(MDImin:MDImax) !------------------------------------------------------------------------- !--- 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). INTEGER, PARAMETER :: Nrime=40 COMMON /IRIME_TABLES/ VEL_RF(2:9,0:Nrime) !------------------------------------------------------------------------- !--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns ! (0.45 mm), assuming an exponential size distribution. REAL, PARAMETER :: DMRmin=.05e-3, DMRmax=.45e-3, DelDMR=1.e-6, & XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax INTEGER, PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax COMMON /RACCR_TABLES/ ACCRR(MDRmin:MDRmax) COMMON /RMASS_TABLES/ MASSR(MDRmin:MDRmax) REAL :: MASSR COMMON /RRATE_TABLES/ RRATE(MDRmin:MDRmax) COMMON /RVELR_TABLES/ VRAIN(MDRmin:MDRmax) COMMON /RVENT_TABLES/ VENTR1(MDRmin:MDRmax), VENTR2(MDRmin:MDRmax) !------------------------------------------------------------------------- !------- Key parameters, local variables, & important comments --------- !----------------------------------------------------------------------- !--- KEY Parameters: !---- Comments on 2 August 2000 ! * EPSQ=1.E-20 is the lower limit for specific humidity and cloud ! condensate. The value of EPSQ will need to be changed in the other ! subroutines in order to make it consistent throughout the Eta code. !- NLImax - maximum number concentration of large ice crystals (20,000 /m**3, 20 per liter) !- NLImin - minimum number concentration of large ice crystals (100 /m**3, 0.1 per liter) REAL, PARAMETER :: CP=1004.6, EPSQ=1.E-20, EPSW=1.E-12, RD=287.04, & RHOL=1000., RV=461.5, T0C=273.15, XLS=2.834E6, EPS=RD/RV, & NLImax=20.E3, NLImin=100., T_ICE=-10., T_ICE_init=-5., & TOLER=5.E-7, & CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT, EPS1=RV/RD-1., RCP=1./CP, & RCPRV=RCP/RV, RRHOL=1./RHOL, XLS1=XLS*RCP, XLS2=XLS*XLS*RCPRV, & XLS3=XLS*XLS/RV, & C1=1./3., C2=1./6., C3=3.31/6., & DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, N0r0=8.E6, N0rmin=1.e4, & N0s0=4.E6, RHO0=1.194, XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, & XMR3=1.e6*DMR3, Xratio=.025 INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3 !--- 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= .TRUE. !--- Local variables REAL :: EMAIRI, N0r, NLICE, NSmICE LOGICAL :: CLEAR, ICE_logical, DBG_logical, DIAG_print, RAIN_logical !####################################################################### !########################## 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) DIAG_print= .FALSE. !----------------------------------------------------------------------- !------------ 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) <= EPSQ .AND. WC_col(L) <= 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 < 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 > WSgrd .OR. WC > EPSW) CLEAR= .FALSE. !--- Check if any rain is falling into layer from above IF (ARAIN > 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 > 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 < 0. .OR. QI > CLIMIT .OR. ASNOW > CLIMIT) THEN ICE_logical= .TRUE. ELSE ICE_logical= .FALSE. QLICE=0. QTICE=0. ENDIF !--- Determine if rain is present RAIN_logical= .FALSE. IF (ARAIN > CLIMIT .OR. QR > CLIMIT) 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 >= 0. .OR. WVQW < QSIgrd) THEN !--- Eliminate small ice particle contributions for melting & sublimation FLARGE=1. ELSE !--- Enhanced number of small ice particles during depositional growth ! (effective only when 0C > T >= T_ice [-10C] ) FLARGE=.2 !--- Larger number of small ice particles due to rime splintering IF (TC >= -8. .AND. TC <= -3.) FLARGE=.5*FLARGE ENDIF ! End IF (TC >= 0. .OR. WVQW < QSIgrd) FSMALL=(1.-FLARGE)/FLARGE XSIMASS=RRHO*MASSI(MDImin)*FSMALL IF (QI <= CLIMIT .AND. ASNOW <= 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 <= 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 >= 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 <= 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 >= NLImin .AND. NLICE <= NLImax) & .OR. IPASS == 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) ) XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1 IF (XLI <= MASSI(MDImin) ) THEN INDEXS=MDImin ELSE IF (XLI <= MASSI(450) ) THEN DLI=9.5885E5*XLI**.42066 ! DLI in microns INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) ) ELSE IF (XLI <= MASSI(MDImax) ) THEN DLI=3.9751E6*XLI**.49870 ! DLI in microns INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) ) ELSE INDEXS=MDImax !--- 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. IF (DUM >= NLImax) & RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS) ENDIF ! End IF (XLI <= MASSI(MDImin) ) ! WRITE(6,"(4(a12,g11.4,1x))") ! & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM, ! & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE, ! & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1 ENDIF ! End IF ( (NLICE >= NLImin .AND. NLICE <= NLImax) ... ENDDO ! End DO IPASS=0,1 ENDIF ! End IF (QI <= CLIMIT .AND. ASNOW <= CLIMIT) ENDIF ! End IF (ICE_logical) !---------------------------------------------------------------------- !--------------- Calculate individual processes ----------------------- !---------------------------------------------------------------------- !--- Cloud water autoconversion to rain and collection by rain IF (QW > CLIMIT .AND. TC >= T_ICE) THEN !--- QW0 could be modified based on land/sea properties, ! presence of convection, etc. This is why QAUT0 and CRAUT ! are passed into the subroutine as externally determined ! parameters. Can be changed in the future if desired. QW0=QAUT0*RRHO PRAUT=MAX(0., QW-QW0)*CRAUT IF (QLICE > CLIMIT) THEN !--- Collection of cloud water by large ice particles ("snow") ! PIACWI=PIACW for riming, PIACWI=0 for shedding FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1) PIACW=FWS*QW IF (TC < 0.) PIACWI=PIACW ! Large ice riming ENDIF ! End IF (QLICE > CLIMIT) ENDIF ! End IF (QW > CLIMIT .AND. TC >= T_ICE) !---------------------------------------------------------------------- !--- Loop around some of the ice-phase processes if no ice should be present !---------------------------------------------------------------------- IF (ICE_logical .EQV. .FALSE. ) GO TO 20 !--- Now the pretzel logic of calculating ice deposition IF (TC < T_ICE .AND. (WV > QSIgrd .OR. QW > CLIMIT)) THEN !--- Adjust to ice saturation at T DUM) PIDEP=DEPOSIT (PP, RHgrd, DUM1, DUM2) DWVi=0. ! Used only for debugging ELSE IF (TC < 0.) THEN !--- These quantities are handy for ice deposition/sublimation ! PIDEP_max - max deposition or minimum sublimation to ice saturation DENOMI=1.+XLS2*QSI*TK2 DWVi=MIN(WVQW,QSW)-QSI PIDEP_max=MAX(PILOSS, DWVi/DENOMI) IF (QTICE > 0.) THEN !--- Calculate ice deposition/sublimation ! * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5], ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS) ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ; ! VENTIL, VENTIS - m**-2 ; VENTI1 - m ; ! VENTI2 - m**2/s**.5 ; DIDEP - unitless SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 ABI=1./(RHO*XLS3*QSI*TK2/THERM_COND+1./DIFFUS) !--- VENTIL - Number concentration * ventilation factors for large ice !--- VENTIS - Number concentration * ventilation factors for small ice !--- Variation in the number concentration of ice with time is not ! accounted for in these calculations (could be in the future). VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE DIDEP=ABI*(VENTIL+VENTIS)*DTPH !--- Account for change in water vapor supply w/ time IF (DIDEP >= Xratio) & DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI IF (DWVi > 0.) THEN PIDEP=MIN(DWVi*DIDEP, PIDEP_max) ELSE IF (DWVi < 0.) THEN PIDEP=MAX(DWVi*DIDEP, PIDEP_max) ENDIF ELSE IF (WVQW > QSI .AND. TC <= T_ICE_init) THEN !--- Ice nucleation in near water-saturated conditions. Ice crystal ! growth during time step calculated using Miller & Young (1979, JAS). !--- These deposition rates could drive conditions below water saturation, ! which is the basis of these calculations. Intended to approximate ! more complex & computationally intensive calculations. INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) ) !--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions !--- DUM2 is the number of ice crystals nucleated at water-saturated ! conditions based on Meyers et al. (JAM, 1992). !--- Prevent unrealistically large ice initiation (limited by PIDEP_max) ! if DUM2 values are increased in future experiments DUM1=QSW/QSI-1. DUM2=1.E3*EXP(12.96*DUM1-.039) PIDEP=MIN(PIDEP_max, DUM2*MY_GROWTH(INDEX_MY)*RRHO) ENDIF ! End IF (QTICE > 0.) ENDIF ! End IF (TC < T_ICE .AND. (WV > QSIgrd .OR. QW > CLIMIT)) !------------------------------------------------------------------------ 20 CONTINUE ! Jump here if conditions for ice are not present !------------------------------------------------------------------------ !--- Cloud water condensation IF (TC >= T_ICE .AND. (QW > CLIMIT .OR. WV > QSWgrd)) THEN IF (PIACWI == 0. .AND. PIDEP == 0.) THEN PCOND=CONDENSE (PP, QW, RHgrd, TK, WV) ELSE !--- Modify cloud condensation in response to ice processes DUM=XLV*QSWgrd*RCPRV*TK2 DENOMWI=1.+XLS*DUM DENOMF=XLF*DUM DUM=MAX(0., PIDEP) PCOND=(WV-QSWgrd-DENOMWI*DUM-DENOMF*PIACWI)/DENOMW DUM1=-QW DUM2=PCOND-PIACW IF (DUM2 < DUM1) THEN !--- Limit cloud water sinks DUM=DUM1/DUM2 PCOND=DUM*PCOND PIACW=DUM*PIACW PIACWI=DUM*PIACWI ENDIF ! End IF (DUM2 < DUM1) ENDIF ! End IF (PIACWI == 0. .AND. PIDEP == 0.) ENDIF ! End IF (TC >= T_ICE .AND. (QW > CLIMIT .OR. WV > QSWgrd)) !--- Limit freezing of accreted rime to prevent temperature oscillations, ! a crude Schumann-Ludlam limit (p. 209 of Young, 1993). TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI IF (TCC > 0.) THEN PIACWI=0. TCC=TC+XLV1*PCOND+XLS1*PIDEP ENDIF IF (TC > 0. .AND. TCC > 0. .AND. ICE_logical) THEN !--- Calculate melting and evaporation/condensation ! * Units: SFACTOR - s**.5/m ; ABI - m**2/s ; NLICE - m**-3 ; ! VENTIL - m**-2 ; VENTI1 - m ; ! VENTI2 - m**2/s**.5 ; CIEVP - /s SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS)) AIEVP=VENTIL*DIFFUS*DTPH IF (AIEVP < Xratio) THEN DIEVP=AIEVP ELSE DIEVP=1.-EXP(-AIEVP) ENDIF QSW0=EPS*ESW0/(PP-ESW0) DWV0=MIN(WV,QSW)-QSW0 DUM=QW+PCOND IF (WV < QSW .AND. DUM <= CLIMIT) THEN !--- Evaporation from melting snow (sink of snow) or shedding ! of water condensed onto melting snow (source of rain) DUM=DWV0*DIEVP PIEVP=MAX( MIN(0., DUM), PILOSS) PICND=MAX(0., DUM) ENDIF ! End IF (WV < QSW .AND. DUM <= CLIMIT) PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF !--- Limit melting to prevent temperature oscillations across 0C DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 ) PIMLT=MIN(PIMLT, DUM1) !--- Limit loss of snow by melting (>0) and evaporation DUM=PIEVP-PIMLT IF (DUM < PILOSS) THEN DUM1=PILOSS/DUM PIMLT=PIMLT*DUM1 PIEVP=PIEVP*DUM1 ENDIF ! End IF (DUM > QTICE) ENDIF ! End IF (TC > 0. .AND. TCC > 0. .AND. ICE_logical) !--- IMPORTANT: Estimate time-averaged properties. ! * TOT_RAIN - total mass of rain before microphysics, which is the sum of ! the total mass of rain in the current layer and the input ! flux of rain from above ! * VRAIN1 - fall speed of rain into grid from above (with air resistance correction) ! * QTRAIN - time-averaged mixing ratio of rain (kg/kg) ! * PRLOSS - greatest loss (<0) of rain, removing all rain falling from ! above and the rain within the layer ! * RQR - rain content (kg/m**3) ! * INDEXR - mean size of rain drops to the nearest 1 micron in size ! * N0r - intercept of rain size distribution (typically 10**6 m**-4) TOT_RAIN=0. VRAIN1=0. QTRAIN=0. PRLOSS=0. RQR=0. N0r=0. INDEXR1=INDEXR ! For debugging only INDEXR=MDRmin IF (RAIN_logical) THEN IF (ARAIN <= 0.) THEN INDEXR=MDRmin VRAIN1=0. ELSE !--- INDEXR (related to mean diameter) & N0r could be modified ! by land/sea properties, presence of convection, etc. !--- Rain rate normalized to a density of 1.194 kg/m**3 RR=ARAIN/(DTPH*GAMMAR) IF (RR <= RR_DRmin) THEN !--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates, ! instead vary N0r with rain rate INDEXR=MDRmin ELSE IF (RR <= RR_DR1) THEN !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables ! for mean diameters (Dr) between 0.05 and 0.10 mm: ! V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m ! RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136, ! RR in kg/(m**2*s) ! Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947 INDEXR=INT( 1.123E3*RR**.1947 + .5 ) INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) ) ELSE IF (RR <= RR_DR2) THEN !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables ! for mean diameters (Dr) between 0.10 and 0.20 mm: ! V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m ! RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958, ! RR in kg/(m**2*s) ! Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017 INDEXR=INT( 1.225E3*RR**.2017 + .5 ) INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) ) ELSE IF (RR <= RR_DR3) THEN !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables ! for mean diameters (Dr) between 0.20 and 0.32 mm: ! V(Dr)=2831.*Dr**.80, V in m/s and Dr in m ! RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80, ! RR in kg/(m**2*s) ! Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083 INDEXR=INT( 1.3006E3*RR**.2083 + .5 ) INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) ) ELSE IF (RR <= RR_DRmax) THEN !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables ! for mean diameters (Dr) between 0.32 and 0.45 mm: ! V(Dr)=944.8*Dr**.6636, V in m/s and Dr in m ! RR = PI*1000.*N0r0*944.8*Dr**(4+.6636) = 2.3745e13*Dr**4.6636, ! RR in kg/(m**2*s) ! Dr (m) = 1.355e-3*RR**.2144 -> Dr (microns) = 1.355e3*RR**.2144 INDEXR=INT( 1.355E3*RR**.2144 + .5 ) INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) ) ELSE !--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates, ! instead vary N0r with rain rate INDEXR=MDRmax ENDIF ! End IF (RR <= RR_DRmin) etc. VRAIN1=GAMMAR*VRAIN(INDEXR) ENDIF ! End IF (ARAIN <= 0.) INDEXR1=INDEXR ! For debugging only TOT_RAIN=THICK*QR+BLEND*ARAIN QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1) PRLOSS=-TOT_RAIN/THICK RQR=RHO*QTRAIN !--- RQR - time-averaged rain content (kg/m**3) IF (RQR <= RQR_DRmin) THEN N0r=MAX(N0rmin, CN0r_DMRmin*RQR) INDEXR=MDRmin ELSE IF (RQR >= RQR_DRmax) THEN N0r=CN0r_DMRmax*RQR INDEXR=MDRmax ELSE N0r=N0r0 INDEXR=MAX( XMRmin, MIN(CN0r0*RQR**.25, XMRmax) ) ENDIF IF (TC < T_ICE) THEN PIACR=-PRLOSS ELSE DWVr=WV-PCOND-QSW DUM=QW+PCOND IF (DWVr < 0. .AND. DUM <= CLIMIT) THEN !--- Rain evaporation ! * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5], ! where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS) ! * Units: RFACTOR - s**.5/m ; ABW - m**2/s ; VENTR - m**-2 ; ! N0r - m**-4 ; VENTR1 - m**2 ; VENTR2 - m**3/s**.5 ; ! CREVP - unitless RFACTOR=GAMMAR**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2 ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS) !--- Note that VENTR1, VENTR2 lookup tables do not include the ! 1/Davg multiplier as in the ice tables VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR)) CREVP=ABW*VENTR*DTPH IF (CREVP < Xratio) THEN DUM=DWVr*CREVP ELSE DUM=DWVr*(1.-EXP(-CREVP*DENOMW))/DENOMW ENDIF PREVP=MAX(DUM, PRLOSS) ELSE IF (QW > CLIMIT) THEN FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR) PRACW=MIN(1.,FWR)*QW ENDIF ! End IF (DWVr < 0. .AND. DUM <= CLIMIT) IF (TC < 0. .AND. TCC < 0.) THEN !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983) PIACR=CBFR*N0r*RRHO*(EXP(ABFR*TC)-1.)*(FLOAT(INDEXR))**7 IF (QLICE > CLIMIT) THEN !--- Freezing of rain by collisions w/ large ice DUM=GAMMAR*VRAIN(INDEXR) DUM1=DUM-VSNOW !--- DUM2 - Difference in spectral fall speeds of rain and ! large ice, parameterized following eq. (48) on p. 112 of ! Murakami (J. Meteor. Soc. Japan, 1990) DUM2=(DUM1*DUM1+.04*DUM*VSNOW)**.5 DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS & +.5E-12*INDEXS*INDEXS FIR=MIN(1., CIACR*NLICE*DUM1*DUM2) !--- Future? Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED??? PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN) ENDIF ! End IF (QLICE > CLIMIT) DUM=PREVP-PIACR If (DUM < PRLOSS) THEN DUM1=PRLOSS/DUM PREVP=DUM1*PREVP PIACR=DUM1*PIACR ENDIF ! End If (DUM < PRLOSS) ENDIF ! End IF (TC < 0. .AND. TCC < 0.) ENDIF ! End IF (TC < T_ICE) ENDIF ! End IF (RAIN_logical) !---------------------------------------------------------------------- !---------------------- Main Budget Equations ------------------------- !---------------------------------------------------------------------- !----------------------------------------------------------------------- !--- Update fields, determine characteristics for next lower layer ---- !----------------------------------------------------------------------- !--- Carefully limit sinks of cloud water DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND) IF (DUM1 > QW) THEN DUM=QW/DUM1 PIACW=DUM*PIACW PIACWI=DUM*PIACWI PRAUT=DUM*PRAUT PRACW=DUM*PRACW IF (PCOND < 0.) PCOND=DUM*PCOND ENDIF PIACWR=PIACW-PIACWI ! TC >= 0C !--- QWnew - updated cloud water mixing ratio DELW=PCOND-PIACW-PRAUT-PRACW QWnew=QW+DELW IF (QWnew <= CLIMIT) QWnew=0. IF (QW > 0. .AND. QWnew /= 0.) THEN DUM=QWnew/QW IF (DUM < TOLER) QWnew=0. ENDIF !--- Update temperature and water vapor mixing ratios DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) & +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT) Tnew=TK+DELT DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP WVnew=WV+DELV !--- Update ice mixing ratios !--- ! * TOT_ICEnew - total mass (small & large) ice after microphysics, ! which is the sum of the total mass of large ice in the ! current layer and the flux of ice out of the grid box below ! * RimeF - Rime Factor, which is the mass ratio of total (unrimed & ! rimed) ice mass to the unrimed ice mass (>=1) ! * QInew - updated mixing ratio of total (large & small) ice in layer ! -> TOT_ICEnew=QInew*THICK+BLDTRH*QLICEnew*VSNOW ! -> But QLICEnew=QInew*FLIMASS, so ! -> TOT_ICEnew=QInew*(THICK+BLDTRH*FLIMASS*VSNOW) ! * ASNOWnew - updated accumulation of snow at bottom of grid cell !--- DELI=0. RimeF=1. IF (ICE_logical) THEN DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT TOT_ICEnew=TOT_ICE+THICK*DELI IF (TOT_ICE > 0. .AND. TOT_ICEnew /= 0.) THEN DUM=TOT_ICEnew/TOT_ICE IF (DUM < TOLER) TOT_ICEnew=0. ENDIF IF (TOT_ICEnew <= CLIMIT) THEN TOT_ICEnew=0. RimeF=1. QInew=0. ASNOWnew=0. ELSE !--- Update rime factor if appropriate DUM=PIACWI+PIACR IF (DUM <= CLIMIT .AND. PIDEP <= CLIMIT) THEN RimeF=RimeF1 ELSE !--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass) ! DUM1 - Total ice mass, rimed & unrimed ! DUM2 - Estimated mass of *unrimed* ice DUM1=TOT_ICE+THICK*(PIDEP+DUM) DUM2=TOT_ICE/RimeF1+THICK*PIDEP IF (DUM2 <= 0.) THEN RimeF=RFmax ELSE RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) ) ENDIF ENDIF ! End IF (DUM <= CLIMIT .AND. PIDEP <= CLIMIT) QInew=TOT_ICEnew/(THICK+BLDTRH*FLIMASS*VSNOW) IF (QInew <= CLIMIT) QInew=0. IF (QI > 0. .AND. QInew /= 0.) THEN DUM=QInew/QI IF (DUM < TOLER) QInew=0. ENDIF ASNOWnew=BLDTRH*FLIMASS*VSNOW*QInew IF (ASNOW > 0. .AND. ASNOWnew /= 0.) THEN DUM=ASNOWnew/ASNOW IF (DUM < TOLER) ASNOWnew=0. ENDIF ENDIF ! End IF (TOT_ICEnew <= CLIMIT) ENDIF ! End IF (ICE_logical) !--- Update rain mixing ratios !--- ! * TOT_RAINnew - total mass of rain after microphysics ! current layer and the input flux of ice from above ! * VRAIN2 - time-averaged fall speed of rain in grid and below ! (with air resistance correction) ! * QRnew - updated rain mixing ratio in layer ! -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2) ! * ARAINnew - updated accumulation of rain at bottom of grid cell !--- DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND TOT_RAINnew=TOT_RAIN+THICK*DELR IF (TOT_RAIN > 0. .AND. TOT_RAINnew /= 0.) THEN DUM=TOT_RAINnew/TOT_RAIN IF (DUM < TOLER) TOT_RAINnew=0. ENDIF IF (TOT_RAINnew <= CLIMIT) THEN TOT_RAINnew=0. VRAIN2=0. QRnew=0. ARAINnew=0. ELSE !--- 1st guess time-averaged rain rate at bottom of grid box RR=TOT_RAINnew/(DTPH*GAMMAR) !--- Use same algorithm as above for calculating mean drop diameter ! (IDR, in microns), which is used to estimate the time-averaged ! fall speed of rain drops at the bottom of the grid layer. This ! isnt perfect, but the alternative is solving a transcendental ! equation that is numerically inefficient and nasty to program ! (coded in earlier versions of GSMCOLUMN prior to 8-22-01). IF (RR <= RR_DRmin) THEN IDR=MDRmin ELSE IF (RR <= RR_DR1) THEN IDR=INT( 1.123E3*RR**.1947 + .5 ) IDR=MAX( MDRmin, MIN(IDR, MDR1) ) ELSE IF (RR <= RR_DR2) THEN IDR=INT( 1.225E3*RR**.2017 + .5 ) IDR=MAX( MDR1, MIN(IDR, MDR2) ) ELSE IF (RR <= RR_DR3) THEN IDR=INT( 1.3006E3*RR**.2083 + .5 ) IDR=MAX( MDR2, MIN(IDR, MDR3) ) ELSE IF (RR <= RR_DRmax) THEN IDR=INT( 1.355E3*RR**.2144 + .5 ) IDR=MAX( MDR3, MIN(IDR, MDRmax) ) ELSE IDR=MDRmax ENDIF ! End IF (RR <= RR_DRmin) VRAIN2=GAMMAR*VRAIN(IDR) QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2) IF (QRnew <= CLIMIT) QRnew=0. IF (QR > 0. .AND. QRnew /= 0.) THEN DUM=QRnew/QR IF (DUM < TOLER) QRnew=0. ENDIF ARAINnew=BLDTRH*VRAIN2*QRnew IF (ARAIN > 0. .AND. ARAINnew /= 0.) THEN DUM=ARAINnew/ARAIN IF (DUM < TOLER) ARAINnew=0. ENDIF ENDIF WCnew=QWnew+QRnew+QInew !---------------------------------------------------------------------- !-------------- Begin debugging & verification ------------------------ !---------------------------------------------------------------------- !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp. QT=THICK*(WV+WC)+ARAIN+ASNOW QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew BUDGET=QT-QTnew !--- Additional check on budget preservation, accounting for truncation effects DBG_logical= .FALSE. ! DUM=ABS(BUDGET) ! IF (DUM .GT. TOLER) THEN ! DUM=DUM/MIN(QT, QTnew) ! IF (DUM .GT. TOLER) DBG_logical=.TRUE. ! ENDIF !! ! DUM=(RHgrd+.001)*QSInew ! IF ( (QWnew.GT.CLIMIT .OR. QRnew.GT.CLIMIT .OR. WVnew.GT.DUM) ! & .AND. TC.LT.T_ICE ) DBG_logical=.TRUE. ! IF (TC.GT.5. .AND. QInew.GT.CLIMIT) DBG_logical=.TRUE. IF ((WVnew < CLIMIT .OR. DBG_logical) .AND. PRINT_diag) THEN WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index, & ' L=',L,' LSFC=',LSFC ESW=1000.*FPVS0(Tnew) QSWnew=EPS*ESW/(PP-ESW) IF (TC < 0. .OR. Tnew < 0.) THEN ESI=1000.*FPVS(Tnew) QSInew=EPS*ESI/(PP-ESI) ELSE QSI=QSW QSInew=QSWnew ENDIF WSnew=QSInew WRITE(6,"(4(a12,g11.4,1x))") & '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, & '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, & 'RHgrd=',RHgrd, & '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, & 'RHInew=',WVnew/QSInew, & '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, & '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, & '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, & '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, & '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, & 'ASNOWnew=',ASNOWnew, & '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, & 'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, & '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew WRITE(6,"(4(a12,g11.4,1x))") & '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, & '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, & '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, & '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', & PIMLT, & '{} PIACR=',PIACR IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") & '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, & 'VSNOW=',VSNOW, & '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, & 'FLIMASS=',FLIMASS, & '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, & 'QTICE=',QTICE, & '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, & 'EMAIRI=',EMAIRI, & '{} RimeF=',RimeF IF (TOT_RAIN > 0. .OR. TOT_RAINnew > 0.) & WRITE(6,"(4(a12,g11.4,1x))") & '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), & 'GAMMAR=',GAMMAR,'N0r=',N0r, & '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, & '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, & 'VOLR2=',THICK+BLDTRH*VRAIN2 IF (PRAUT > 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0 IF (PRACW > 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR IF (PIACR > 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR DUM=PIMLT+PICND-PREVP-PIEVP IF (DUM > 0. .OR. DWVi /= 0.) & WRITE(6,"(4(a12,g11.4,1x))") & '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, & 'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS IF (PREVP < 0.) WRITE(6,"(4(a12,g11.4,1x))") & '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, & '{} DWVr=',DWVr,'DENOMW=',DENOMW IF (PIDEP /= 0. .AND. DWVi /= 0.) & WRITE(6,"(4(a12,g11.4,1x))") & '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, & 'SFACTOR=',SFACTOR, & '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), & '{} VENTIS=',VENTIS,'DIDEP=',DIDEP IF (PIDEP > 0. .AND. PCOND /= 0.) & WRITE(6,"(4(a12,g11.4,1x))") & '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, & 'DUM2=',PCOND-PIACW IF (FWS > 0.) WRITE(6,"(4(a12,g11.4,1x))") & '{} FWS=',FWS DUM=PIMLT+PICND-PIEVP IF (DUM > 0.) WRITE(6,"(4(a12,g11.4,1x))") & '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), & 'VENTIL2=',SFACTOR*VENTI2(INDEXS), & '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0 ENDIF !---------------------------------------------------------------------- !-------------- Water budget statistics & maximum values -------------- !---------------------------------------------------------------------- IF (PRINT_diag) THEN ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) ) IF (QInew > CLIMIT) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1 IF (QInew > CLIMIT .AND. QRnew+QWnew > CLIMIT) & NSTATS(ITdx,2)=NSTATS(ITdx,2)+1 IF (QWnew > CLIMIT) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1 IF (QRnew > CLIMIT) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1 QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew) QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew) QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew) QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew) QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew) QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN) QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW) IF (QInew > 0.) & QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF ENDIF !---------------------------------------------------------------------- !------------------------- Update arrays ------------------------------ !---------------------------------------------------------------------- T_col(L)=Tnew ! Updated temperature QV_col(L)=max(EPSQ, WVnew/(1.+WVnew)) ! Updated specific humidity WC_col(L)=max(EPSQ, WCnew) ! Updated total condensate mixing ratio QI_col(L)=max(EPSQ, QInew) ! Updated ice mixing ratio QR_col(L)=max(EPSQ, QRnew) ! Updated rain mixing ratio QW_col(L)=max(EPSQ, QWnew) ! Updated cloud water mixing ratio RimeF_col(L)=RimeF ! Updated rime factor ASNOW=ASNOWnew ! Updated accumulated snow ARAIN=ARAINnew ! Updated accumulated rain !####################################################################### CONTINUE ! ##### End "L" loop through model levels ##### 10 END DO !####################################################################### !----------------------------------------------------------------------- !--------------------------- Return to GSMDRIVE ----------------------- !----------------------------------------------------------------------- RETURN END SUBROUTINE GSMCOLUMN !####################################################################### !--------- Produces accurate calculation of cloud condensation --------- !####################################################################### REAL FUNCTION CONDENSE (PP, QW, RHgrd, TK, WV) !--------------------------------------------------------------------------------- !------ The Asai (1965) algorithm takes into consideration the release of ------ !------ latent heat in increasing the temperature & in increasing the ------ !------ saturation mixing ratio (following the Clausius-Clapeyron eqn.). ------ !--------------------------------------------------------------------------------- INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15) REAL (KIND=HIGH_PRES), PARAMETER :: CLIMIT=1.E-20, RHLIMIT=.001, RHLIMIT1=-RHLIMIT REAL, PARAMETER :: CP=1004.6, RD=287.04, RV=461.5, EPS=RD/RV, RCP=1./CP, RCPRV=RCP/RV REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum !----------------------------------------------------------------------- !--- LV (T) is from Bolton (JAS, 1980) XLV=3.148E6-2370.*TK XLV1=XLV*RCP XLV2=XLV*XLV*RCPRV Tdum=TK WVdum=WV WCdum=QW ESW=1000.*FPVS0(Tdum) ! Saturation vapor press w/r/t water WS=RHgrd*EPS*ESW/(PP-ESW) ! Saturation mixing ratio w/r/t water DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio SSAT=DWV/WS ! Supersaturation ratio CONDENSE=0. DO WHILE ((SSAT < RHLIMIT1 .AND. WCdum > CLIMIT) & .OR. SSAT > RHLIMIT) COND=DWV/(1.+XLV2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan) COND=MAX(COND, -WCdum) ! Limit cloud water evaporation Tdum=Tdum+XLV1*COND ! Updated temperature WVdum=WVdum-COND ! Updated water vapor mixing ratio WCdum=WCdum+COND ! Updated cloud water mixing ratio CONDENSE=CONDENSE+COND ! Total cloud water condensation ESW=1000.*FPVS0(Tdum) ! Updated saturation vapor press w/r/t water WS=RHgrd*EPS*ESW/(PP-ESW) ! Updated saturation mixing ratio w/r/t water DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio SSAT=DWV/WS ! Grid-scale supersaturation ratio ENDDO RETURN END FUNCTION CONDENSE !####################################################################### !---------------- Calculate ice deposition at T RHLIMIT .OR. SSAT < RHLIMIT1) !--- Note that XLVS2=LS*LV/(CP*RV)=LV*WS/(RV*T*T)*(LS/CP*DEP1), ! where WS is the saturation mixing ratio following Clausius- ! Clapeyron (see Asai,1965; Young,1993,p.405) DEP=DWV/(1.+XLS2*WS/(Tdum*Tdum)) ! Asai (1965, J. Japan) Tdum=Tdum+XLS1*DEP ! Updated temperature WVdum=WVdum-DEP ! Updated ice mixing ratio DEPOSIT=DEPOSIT+DEP ! Total ice deposition ESI=1000.*FPVS(Tdum) ! Updated saturation vapor press w/r/t ice WS=RHgrd*EPS*ESI/(PP-ESI) ! Updated saturation mixing ratio w/r/t ice DWV=WVdum-WS ! Deficit grid-scale water vapor mixing ratio SSAT=DWV/WS ! Grid-scale supersaturation ratio ENDDO RETURN END FUNCTION DEPOSIT