!--------------------------------------------------------------------------------------------------
!  DOXYGEN
!> @brief GRID-SCALE MICROPHYSICAL PROCESSES - CONDENSATION AND PRECIPITATION
!> @details MERGES ORIGINAL GSCOND AND 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.
!!
!! 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 AND 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.
!!
!! NOTE:  CODE IS CURRENTLY SET UP W/O THREADING
!> @author ORIGINATOR - FERRIER 
!> @date 01-08-?? \n
!> @author LUCCI 
!> @date 18-03-20 \n
!> @version V1.1.0
!> @details MODERNIZATION OF THE CODE, INCLUDING:
!!                      * F77 TO F90/F95
!!                      * INDENTATION & UNIFORMIZATION CODE
!!                      * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!!                      * DOCUMENTATION WITH DOXYGEN
!!                      * OPENMP FUNCTIONALITY 
!<
!> @param[in] DTPH          - PHYSICS TIME STEP (S)
!> @param[in] I_INDEX       - I INDEX
!> @param[in] J_INDEX       - J INDEX
!> @param[in] LSFC          - ETA LEVEL OF LEVEL ABOVE SURFACE, GROUND
!> @param[in] P_COL         - VERTICAL COLUMN OF MODEL PRESSURE (PA)
!> @param[in] THICK_COL     - VERTICAL COLUMN OF MODEL MASS THICKNESS (DENSITY*HEIGHT INCREMENT)
!> @param[out] ARAIN        - ACCUMULATED RAINFALL AT THE SURFACE (KG)
!> @param[out] ASNOW        - ACCUMULATED SNOWFALL AT THE SURFACE (KG)
!> @param[inout] QI_COL     - VERTICAL COLUMN OF MODEL ICE MIXING RATIO (KG/KG)
!> @param[inout] QR_COL     - VERTICAL COLUMN OF MODEL RAIN RATIO (KG/KG)
!> @param[inout] QV_COL     - VERTICAL COLUMN OF MODEL WATER VAPOR SPECIFIC HUMIDITY (KG/KG)
!> @param[inout] QW_COL     - VERTICAL COLUMN OF MODEL CLOUD WATER MIXING RATIO (KG/KG)
!> @param[inout] RIMEF_COL  - VERTICAL COLUMN OF RIME FACTOR FOR ICE IN MODEL (RATIO, DEFINED BELOW)
!> @param[inout] T_COL      - VERTICAL COLUMN OF MODEL TEMPERATURE (DEG K)
!> @param[inout] WC_COL     - VERTICAL COLUMN OF MODEL MIXING RATIO OF TOTAL CONDENSATE (KG/KG)
!> @details <b>Use Module:</b>
!! @arg @c CMICRO_CONS
!! @arg @c CMICRO_STATS
!! @arg @c CMY600
!! @arg @c F77KINDS
!! @arg @c IACCR_TABLES
!! @arg @c IMASS_TABLES
!! @arg @c IRATE_TABLES
!! @arg @c IRIME_TABLES
!! @arg @c IVENT_TABLES
!! @arg @c MPPSTAFF
!! @arg @c PARMETA
!! @arg @c RACCR_TABLES
!! @arg @c RMASS_TABLES
!! @arg @c RVELR_TABLES
!! @arg @c RVENT_TABLES
!> @details <b>Driver:</b> 
!! @arg @c GSMDRIVE
!> @details <b>Calls:</b>
!--------------------------------------------------------------------------------------------------
    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)
!--------------------------------------------------------------------------------------------------
! SUBROUTINE GSMCOLUMN
!
! SUBPROGRAM: GSMCOLUMN - GRID-SCALE MICROPHYSICAL PROCESSES - CONDENSATION AND PRECIPITATION
! PROGRAMMER: FERRIER
! ORG: W/NP22
! DATE: 01-08-??
!
! ABSTRACT: 
! MERGES ORIGINAL GSCOND AND 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.
!
! 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 AND 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.
!
! NOTE:  CODE IS CURRENTLY SET UP W/O THREADING
! 
! PROGRAM HISTORY LOG:
! 01-08-??  FERRIER    - ORIGINATOR
! 18-01-15  LUCCI      - MODERNIZATION OF THE CODE, INCLUDING:
!                        * F77 TO F90/F95
!                        * INDENTATION & UNIFORMIZATION CODE
!                        * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!                        * DOCUMENTATION WITH DOXYGEN
!                        * OPENMP FUNCTIONALITY
!
! 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)
! THICK_COL  - VERTICAL COLUMN OF MODEL MASS THICKNESS (DENSITY*HEIGHT INCREMENT)
!
! OUTPUT ARGUMENT LIST:
! ARAIN      - ACCUMULATED RAINFALL AT THE SURFACE (KG)
! ASNOW      - ACCUMULATED SNOWFALL AT THE SURFACE (KG)
!
! INPUT/OUTPUT ARGUMENT LIST:
! 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)
! WC_COL     - VERTICAL COLUMN OF MODEL MIXING RATIO OF TOTAL CONDENSATE (KG/KG)
!
! USE MODULES: CMICRO_CONS
!              CMICRO_STATS
!              CMY600
!              F77KINDS
!              GLB_TABLE
!              IACCR_TABLES 
!              IMASS_TABLES
!              IRATE_TABLES
!              IRIME_TABLES
!              IVENT_TABLES
!              MAPPINGS
!              MPPCOM
!              PARMETA
!              RACCR_TABLES
!              RRATE_TABLES
!              RVELR_TABLES
!              RVENT_TABLES
!              TEMPCOM
!              TOPO
!
! DRIVER     : GSMDRIVE
!
! CALLS      : -----
!--------------------------------------------------------------------------------------------------
    USE CMICRO_CONS
    USE CMICRO_STATS
    USE CMY600
    USE F77KINDS
    USE IACCR_TABLES
    USE IMASS_TABLES
    USE IRATE_TABLES
    USE IRIME_TABLES
    USE IVENT_TABLES
    USE MPPSTAFF
    USE PARMETA
    USE RACCR_TABLES
    USE RMASS_TABLES
    USE RVELR_TABLES
    USE RVENT_TABLES
!
    IMPLICIT NONE
!-------------------------------------- 
! ARRAYS AND CONSTANTS IN ARGUMENT LIST
!--------------------------------------
!    REAL   (KIND=R4)                                                      , INTENT(OUT)         ::&
    REAL   (KIND=R4)                                                                            ::&
    & ARAIN   , ASNOW      
!
!    REAL   (KIND=R4)    , DIMENSION(LM)                                   , INTENT(IN)          ::&
    REAL   (KIND=R4)    , DIMENSION(LM)                                                         ::& 
    & P_COL   , THICK_COL 
!
!    REAL   (KIND=R4)    , DIMENSION(LM)                                   , INTENT(INOUT)       ::&
    REAL   (KIND=R4)    , DIMENSION(LM)                                                         ::& 
    & QI_COL  , QR_COL  , QV_COL  , QW_COL  , RIMEF_COL, T_COL  , WC_COL
!
!    INTEGER(KIND=I4)                                                      , INTENT(IN)          ::&
    INTEGER(KIND=I4)                                                                            ::&
    & I_INDEX , J_INDEX , LSFC
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
!    REAL   (KIND=R4)                                                      , INTENT(IN)          ::&
    REAL   (KIND=R4)                                                                            ::&
    & DTPH 
!--------------------------------------------------------------------------------------------------
! KEY PARAMETERS, LOCAL VARIABLES, AND 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   (KIND=R4)    , PARAMETER :: CP         = 1004.6
    REAL   (KIND=R4)    , PARAMETER :: EPSQ       =    1.E-20      
    REAL   (KIND=R4)    , PARAMETER :: EPSW       =    1.E-12    
    REAL   (KIND=R4)    , PARAMETER :: RD         =  287.04    
    REAL   (KIND=R4)    , PARAMETER :: RHOL       = 1000.   
    REAL   (KIND=R4)    , PARAMETER :: RV         =  461.5 
    REAL   (KIND=R4)    , PARAMETER :: T0C        =  273.15        
    REAL   (KIND=R4)    , PARAMETER :: XLS        =    2.834E6    
    REAL   (KIND=R4)    , PARAMETER :: EPS        = RD         / RV     
    REAL   (KIND=R4)    , PARAMETER :: NLIMAX     =   20.E3   
    REAL   (KIND=R4)    , PARAMETER :: NLIMIN     =  100.    
    REAL   (KIND=R4)    , PARAMETER :: T_ICE      =  -10.        
    REAL   (KIND=R4)    , PARAMETER :: T_ICE_INIT =   -5. 
    REAL   (KIND=R4)    , PARAMETER :: TOLER      =    5.E-7   
    REAL   (KIND=R4)    , PARAMETER :: CLIMIT     =   10.      * EPSQ      
    REAL   (KIND=R4)    , PARAMETER :: CLIMIT1    = -CLIMIT 
    REAL   (KIND=R4)    , PARAMETER :: EPS1       = RV         / RD - 1.     
    REAL   (KIND=R4)    , PARAMETER :: RCP        =    1.      / CP      
    REAL   (KIND=R4)    , PARAMETER :: RCPRV      = RCP        / RV  
    REAL   (KIND=R4)    , PARAMETER :: RRHOL      =    1.      / RHOL  
    REAL   (KIND=R4)    , PARAMETER :: XLS1       = XLS        * RCP   
    REAL   (KIND=R4)    , PARAMETER :: XLS2       = XLS        * XLS * RCPRV 
    REAL   (KIND=R4)    , PARAMETER :: XLS3       = XLS        * XLS / RV
    REAL   (KIND=R4)    , PARAMETER :: C1         =    1.      / 3.      
    REAL   (KIND=R4)    , PARAMETER :: C2         =    1.      / 6.       
    REAL   (KIND=R4)    , PARAMETER :: C3         =    3.31    / 6.     
    REAL   (KIND=R4)    , PARAMETER :: DMR1       =     .1E-3        
    REAL   (KIND=R4)    , PARAMETER :: DMR2       =     .2E-3     
    REAL   (KIND=R4)    , PARAMETER :: DMR3       =     .32E-3   
    REAL   (KIND=R4)    , PARAMETER :: N0R0       =    8.E6      
    REAL   (KIND=R4)    , PARAMETER :: N0RMIN     =    1.E4    
    REAL   (KIND=R4)    , PARAMETER :: N0S0       =    4.E6         
    REAL   (KIND=R4)    , PARAMETER :: RHO0       =    1.194    
    REAL   (KIND=R4)    , PARAMETER :: XMR1       =    1.E6    * DMR1 
    REAL   (KIND=R4)    , PARAMETER :: XMR2       =    1.E6    * DMR2 
    REAL   (KIND=R4)    , PARAMETER :: XMR3       =    1.E6    * DMR3 
    REAL   (KIND=R4)    , PARAMETER :: XRATIO     =     .025
!
    INTEGER(KIND=I4)    , PARAMETER :: MDR1       = XMR1
    INTEGER(KIND=I4)    , PARAMETER :: MDR2       = XMR2
    INTEGER(KIND=I4)    , PARAMETER :: 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   (KIND=R4)    , PARAMETER :: BLEND = 1.
!----------------------------------------------------  
! THIS VARIABLE IS FOR DEBUGGING PURPOSES (IF .TRUE.)
!---------------------------------------------------- 
    LOGICAL(KIND=L4)    , PARAMETER :: PRINT_DIAG = .TRUE. 
!----------------  
! LOCAL VARIABLES
!----------------
    REAL   (KIND=R4)                                                                            ::&
    & EMAIRI  , N0R     , NLICE   , NSMICE
!
    LOGICAL(KIND=L4)                                                                            ::&
    & CLEAR   , ICE_LOGICAL       , DBG_LOGICAL           , DIAG_PRINT    , RAIN_LOGICAL
!------------------------
! IMPLICIT NONE VARIABLES
!------------------------
    REAL   (KIND=R4)                                                                            ::&
    & TK       , TC      , PP        , QV       , WV       , ESW     , FPVS0   , QSW       ,      &
    & WS       , ESI     , FPVS      , QSI      , QSWGRD   , QSIGRD  , WSGRD   , RHO       ,      &
    & DTRHO    , BLDTRH  , THICK     , ARAINNEW , ASNOWNEW , QI      , WC      , RRHO      ,      &
    & QINEW    , QR      , QRNEW     , QW       , QWNEW    , PCOND   , PIDEP   , PIACW     ,      &
    & PIACWI   , PIACWR  , PIACR     , PICND    , PIEVP    , PIMLT   , PRAUT   , PRACW     ,      &
    & PREVP    , XLV     , XLF       , XLV1     , XLF1     , TK2     , XLV2    , DENOMW    ,      &
    & TFACTOR  , DYNVIS  , THERM_COND, DIFFUS   , GAMMAS   , GAMMAR  , QLICE   , QTICE     ,      &
    & WVQW     , FLARGE  , FSMALL    , XSIMASS  , TOT_ICE  , PILOSS  , RIMEF1  , TOT_ICENEW,      &
    & VRIMEF   , VEL_INC , VSNOW     , XLIMASS  , FLIMASS  , DUM     , DUM1    , DUM2      ,      &
    & XRF      , XLI     , DLI       , QW0      , FWS      , DEPOSIT , DWVI    , DENOMI    ,      &
    & PIDEP_MAX, SFACTOR , ABI       , VENTIL   , VENTIS   , DIDEP   , CONDENSE, DENOMWI   ,      &
    & DENOMF   , TCC     , AIEVP     , DIEVP    , QSW0     , DWV0    , TOT_RAIN, VRAIN1    ,      &
    & QTRAIN   , PRLOSS  , RQR       , RR       , DWVR     , RFACTOR ,                            &
    & ABW      , VENTR   , CREVP     , FWR      , FIR      , DELW    , DELT    , TNEW      ,      &
    & DELV     , WVNEW   , DELI      , RIMEF    , VRAIN2   , WCNEW   , QT      , QTNEW     ,      &
    & BUDGET   , QSWNEW  , QSINEW    , WSNEW    , DELR     , TOT_RAINNEW
!
    INTEGER(KIND=I4)                                                                            ::&
    & L        , LBEF    , IPASS     , IXS      , IXRF     , INDEX_MY, IDR     , ITDX      ,      &
    & INDEXR1  , INDEXR  , INDEXS  
!---------------- 
! 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
!-------------------------------------------------------------------------------------------------
!CHOU        IF (QV_COL(L) <= EPSQ .AND. WC_COL(L) <= EPSQ) GOTO 10
        IF ((QV_COL(L) <= EPSQ .AND. WC_COL(L) <= EPSQ) .OR. (P_COL(L) <= 2500.)) GOTO 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
!-------------------------------------------------------------------------------   
! JBF: Tentando achar o erro. 
        ESW = 1000. * FPVS0(TK)          ! SATURATION VAPOR PRESSURE W/R/T WATER
!        if (mype ==0) write(*,*) "GSMCOLUMN", "FPVS0(TK)=", FPVS0(TK) 
!        if (mype ==0) write(*,*) "GSMCOLUMN", "PP=", P_COL(L), "L=", L
!        if (mype ==0) write(*,*) "GSMCOLUMN", "TK=", T_COL(L), "L=", L

        QSW = EPS   * ESW / (PP-ESW)     ! SATURATION MIXING RATIO W/R/T WATER
        WS  = QSW                        ! GENERAL SATURATION MIXING RATIO (WATER/ICE)
        QSI = QSW                        ! SATURATION MIXING RATIO W/R/T 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)
        END IF
!----------------------------------------------     
! 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.
        END IF
!-----------------------------------------------------------------------------------------    
! 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.
        END IF
!----------------------------------------------------------------- 
! LOOP TO THE END IF IN CLEAR, SUBSATURATED AIR FREE OF CONDENSATE 
!----------------------------------------------------------------- 
        IF (CLEAR) GOTO 10
!-------------------------------------------------------------------------------------------------- 
! INITIALIZE RHO, THICK AND 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
!--------------------------------------------- 
!    
!------------------------------------------ 
! 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.
        END IF
!-----------------------------  
! 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 AND CLOUD WATER
!        
            IF (TC >= 0. .OR. WVQW < QSIGRD) THEN
!-----------------------------------------------------------------------              
! ELIMINATE SMALL ICE PARTICLE CONTRIBUTIONS FOR MELTING AND 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
            
            END IF  ! 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))
                        END IF
!
                    END IF             ! 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)
!
                        END IF ! END IF (XLI <= MASSI(MDIMIN) )
                    END IF ! END IF ( (NLICE >= NLIMIN .AND. NLICE <= NLIMAX) 
                END DO     ! END DO IPASS=0,1
            END IF         ! END IF (QI <= CLIMIT .AND. ASNOW <= CLIMIT)
        END IF             ! 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
            END IF                            ! END IF (QLICE > CLIMIT)
        END IF                                ! 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. ) GOTO 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<T_ICE (-10C) IF SUPERSATURATED.
! SOURCES OF ICE DUE TO NUCLEATION AND CONVECTIVE DETRAINMENT ARE EITHER POORLY UNDERSTOOD, POORLY
! RESOLVED AT TYPICAL NWP RESOLUTIONS, OR ARE NOT REPRESENTED 
! (E.G., NO DETRAINED CONDENSATE IN BMJ CU SCHEME).
!--------------------------------------------------------------------------------------------------    
            PCOND = -QW
            DUM1  =  TK + XLV1 * PCOND                 ! UPDATED (DUMMY) TEMPERATURE (DEG K)
            DUM2  =  WV + QW                           ! UPDATED (DUMMY) WATER VAPOR MIXING RATIO
            DUM   =  1000. * FPVS(DUM1)                ! UPDATED (DUMMY) SATURATION VAPOR PRESSURE W/R/T ICE
            DUM   =  RHgrd * EPS * DUM / (PP - DUM)    ! UPDATED (DUMMY) SATURATION MIXING RATIO W/R/T ICE
!
            IF (DUM2 > 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)
                END IF
!            
            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 AND 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)
            END IF       ! END IF (QTICE > 0.)
!        
        END IF           ! 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
                END IF ! END IF (DUM2 < DUM1)
            END IF     ! END IF (PIACWI == 0. .AND. PIDEP == 0.)
        END IF         ! 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
        END IF
 !   
        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)
            END IF
!
            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)
            END IF            ! 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
            END IF           ! END IF (DUM > QTICE)
        END IF               ! 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
                END IF              ! END IF (RR <= RR_DRMIN) ETC.
                VRAIN1 = GAMMAR * VRAIN(INDEXR)
            END IF                  ! 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))
            END IF
!        
            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
                    END IF
!
                    PREVP = MAX(DUM, PRLOSS)
!
                ELSE IF (QW > CLIMIT) THEN
                    FWR   = CRACW * GAMMAR * N0R * ACCRR(INDEXR)
                    PRACW = MIN(1.,FWR) * QW
                END IF       ! 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)
                    END IF   ! END IF (QLICE > CLIMIT)
!
                    DUM = PREVP - PIACR
!
                    If (DUM < PRLOSS) THEN
                        DUM1  = PRLOSS / DUM
                        PREVP = DUM1   * PREVP
                        PIACR = DUM1   * PIACR
                    END IF ! END IF (DUM < PRLOSS)
                END IF     ! END IF (TC < 0. .AND. TCC < 0.)
            END IF         ! END IF (TC < T_ICE)
        END IF             ! 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
        END IF
!
        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.
        END IF
!-------------------------------------------------    
! 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 AND 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 AND 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, MDRMIN
!              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.
            END IF
!
            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))
                    END IF
                END IF   ! 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.
                END IF
!
                ASNOWNEW = BLDTRH * FLIMASS * VSNOW * QINEW
!
                IF (ASNOW > 0. .AND. ASNOWNEW /= 0.) THEN
                    DUM = ASNOWNEW / ASNOW
                    IF (DUM < TOLER) ASNOWNEW = 0.
                END IF
            END IF     ! END IF (TOT_ICENEW <= CLIMIT)
        END IF         ! 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.
        END IF
!
        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
            END IF  ! 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.
            END IF
!
            ARAINNEW = BLDTRH * VRAIN2 * QRNEW
            IF (ARAIN > 0. .AND. ARAINNEW /= 0.) THEN
                DUM = ARAINNEW / ARAIN
                IF (DUM < TOLER) ARAINNEW = 0.
            END IF
        END IF
!    
        WCNEW = QWNEW + QRNEW + QINEW
!--------------------------------- 
! BEGIN DEBUGGING AND VERIFICATION  
!--------------------------------- 
!
!------------------------------------------------------------------------------------ 
! QT, QTNEW - TOTAL WATER (VAPOR AND CONDENSATE) BEFORE AND 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. 
!
        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
            END IF
!
            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
!        
        END IF
!------------------------------------------- 
! WATER BUDGET STATISTICS AND 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
!        
        END IF
!-------------- 
! 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
!
!
!--------------------------------------------------------------------------------------------------
!  DOXYGEN
!> @brief Inserir Brief de CONDENSE
!> @details Inserir Details de CONDENSE
!> @author ORIGINATOR - ????? 
!> @date ??-??-?? \n
!> @author LUCCI 
!> @date 18-03-20 \n
!> @version V1.1.0
!> @details MODERNIZATION OF THE CODE, INCLUDING:
!!                      * F77 TO F90/F95
!!                      * INDENTATION & UNIFORMIZATION CODE
!!                      * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!!                      * DOCUMENTATION WITH DOXYGEN
!!                      * OPENMP FUNCTIONALITY 
!<
!> @param[in] PP    - Significado de PP
!> @param[in] QW  - Significado de QW
!> @param[in] RHGRD  - Significado de RHGRD
!> @param[in] TK  - Significado de TK
!> @param[in] WV  - Significado de WV
!> @details <b>Use Module:</b>
!! @arg @c F77KINDS
!> @details <b>Driver:</b> 
!! @arg @c GSMCOLUMN
!> @details <b>Calls:</b>
!--------------------------------------------------------------------------------------------------
!---------------------------------------------------- 
! PRODUCES ACCURATE CALCULATION OF CLOUD CONDENSATION 
!---------------------------------------------------- 
    REAL FUNCTION CONDENSE(PP, QW, RHGRD, TK, WV)
!
    USE F77KINDS
!--------------------------------------------------------------------------------------------------
! THE ASAI (1965) ALGORITHM TAKES INTO CONSIDERATION THE RELEASE OF LATENT HEAT IN INCREASING THE 
! TEMPERATURE AND 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
    REAL   (KIND=HIGH_PRES), PARAMETER :: RHLIMIT  =  .001
    REAL   (KIND=HIGH_PRES), PARAMETER :: RHLIMIT1 = -RHLIMIT
!
    REAL   (KIND=R4)    , PARAMETER :: CP    = 1004.6
    REAL   (KIND=R4)    , PARAMETER :: RD    =  287.04
    REAL   (KIND=R4)    , PARAMETER :: RV    =  461.5
    REAL   (KIND=R4)    , PARAMETER :: EPS   =  RD / RV
    REAL   (KIND=R4)    , PARAMETER :: RCP   =  1. / CP
    REAL   (KIND=R4)    , PARAMETER :: 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 ICE
    WS       = RHGRD * EPS * ESW / (PP - ESW)      ! SATURATION MIXING RATIO
    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
    END DO
!
    RETURN
!
    END FUNCTION CONDENSE
!
!
!--------------------------------------------------------------------------------------------------
!  DOXYGEN
!> @brief Inserir Brief de DEPOSIT
!> @details Inserir Details de DEPOSIT
!> @author ORIGINATOR - ????? 
!> @date ??-??-?? \n
!> @author LUCCI 
!> @date 18-03-20 \n
!> @version V1.1.0
!> @details MODERNIZATION OF THE CODE, INCLUDING:
!!                      * F77 TO F90/F95
!!                      * INDENTATION & UNIFORMIZATION CODE
!!                      * REPLACEMENT OF COMMONS BLOCK FOR MODULES
!!                      * DOCUMENTATION WITH DOXYGEN
!!                      * OPENMP FUNCTIONALITY 
!<
!> @param[in] PP    - Significado de PP
!> @param[in] RHGRD  - Significado de RHGRD
!> @param[in] TDUM  - Significado de TDUM
!> @param[in] WVDUM  - Significado de WVDUM
!> @details <b>Use Module:</b>
!! @arg @c F77KINDS
!> @details <b>Driver:</b> 
!! @arg @c GSMCOLUMN
!> @details <b>Calls:</b>
!--------------------------------------------------------------------------------------------------
!------------------------------------
! CALCULATE ICE DEPOSITION AT T<T_ICE 
!------------------------------------
    REAL FUNCTION DEPOSIT(PP, RHGRD, TDUM, WVDUM)
!
    USE F77KINDS
!--------------------------------------------------------------------------------------------------
! ALSO USES THE ASAI (1965) ALGORITHM BUT USES A DIFFERENT TARGET VAPOR PRESSURE FOR THE ADJUSTMENT
!--------------------------------------------------------------------------------------------------
    INTEGER                , PARAMETER :: HIGH_PRES=SELECTED_REAL_KIND(15)
!
    REAL   (KIND=HIGH_PRES), PARAMETER :: RHLIMIT  =      .001
    REAL   (KIND=HIGH_PRES), PARAMETER :: RHLIMIT1 = -RHLIMIT
!
    REAL   (KIND=R4)    , PARAMETER :: CP       =  1004.6 
    REAL   (KIND=R4)    , PARAMETER :: RD       =   287.04
    REAL   (KIND=R4)    , PARAMETER :: RV       =   461.5
    REAL   (KIND=R4)    , PARAMETER :: XLS      =     2.834E6 
    REAL   (KIND=R4)    , PARAMETER :: EPS      = RD / RV
    REAL   (KIND=R4)    , PARAMETER :: RCP      =     1. / CP
    REAL   (KIND=R4)    , PARAMETER :: RCPRV    = RCP / RV
    REAL   (KIND=R4)    , PARAMETER :: XLS1     = XLS * RCP
    REAL   (KIND=R4)    , PARAMETER :: XLS2     = XLS * XLS * RCPRV
!
    REAL   (KIND=HIGH_PRES)                                                                     ::&
    & DEP     , SSAT
!
    ESI     = 1000. * FPVS(TDUM)             ! SATURATION VAPOR PRESS W/R/T ICE
    WS      = RHGRD * EPS * ESI / (PP - ESI) ! SATURATION MIXING RATIO
    DWV     = WVDUM - WS                     ! DEFICIT GRID-SCALE WATER VAPOR MIXING RATIO
    SSAT    = DWV   / WS                     ! SUPERSATURATION RATIO
    DEPOSIT = 0.
!
    DO WHILE (SSAT > 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
    END DO
!
    RETURN
!
    END FUNCTION DEPOSIT
