Skip to content
Snippets Groups Projects
ice_adjust.f90 24.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • !MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
    
    !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
    
    !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
    
    !MNH_LIC for details. version 1.
    
    !-----------------------------------------------------------------
    !     ######################
          MODULE MODI_ICE_ADJUST
    !     ######################
    !
    INTERFACE
    !
    
          SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, KMI, HRAD,                    &
    
                                 HTURBDIM, OSUBG_COND, OSIGMAS, PTSTEP, PSIGQSAT,  &
    
                                 PRHODJ, PEXNREF, PSIGS,PMFCONV,PPABST,PZZ,        &
                                 PCF_MF,PRC_MF, PRI_MF,                            &
    
                                 PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR ,     &
                                 PRRT, PRRS, PRIT, PRIS, PRST,PRSS, PRGT, PRGS,    &
                                 PRHT, PRHS                                        )
    !
    INTEGER,                  INTENT(IN)    :: KKA   !near ground array index  
    INTEGER,                  INTENT(IN)    :: KKU   !uppest atmosphere array index
    INTEGER,                  INTENT(IN)    :: KKL   !vert. levels type 1=MNH -1=ARO
    INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
    INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
    CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
                                                        ! turbulence scheme
    CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
    LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                        ! Condensation
    LOGICAL                                 :: OSIGMAS  ! Switch for Sigma_s: 
                                                        ! use values computed in CONDENSATION
                                                        ! or that from turbulence scheme
    REAL,                     INTENT(IN)   :: PTSTEP    ! Double Time step
                                                        ! (single if cold start)
    REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXNREF ! Reference Exner function
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PSIGS   ! Sigma_s at time t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PMFCONV ! convective mass flux
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t        
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PZZ     ! height of model layer
    !
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction 
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRI_MF! Convective Mass Flux ice mixing ratio
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRVT    ! Water vapor m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRCT    ! Cloud water m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
    REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
                                                       ! s'rc'/2Sigma_s2 at time t+1
                                                       ! multiplied by Lambda_3
    REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRRS ! Rain water m.r. at t+1
    REAL, DIMENSION(:,:,:),   INTENT(INOUT)::  PRIS ! Cloud ice  m.r. at t+1
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRSS ! Aggregate  m.r. at t+1
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRGS ! Graupel    m.r. at t+1
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRRT ! Rain water m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRIT ! Cloud ice  m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRST ! Aggregate  m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRGT ! Graupel    m.r. at t
    REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRHT ! Hail       m.r. at t
    REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRHS ! Hail       m.r. at t+1
    !
    !
    END SUBROUTINE ICE_ADJUST
    !
    END INTERFACE
    !
    END MODULE MODI_ICE_ADJUST
    !
    !     ##########################################################################
    
          SUBROUTINE ICE_ADJUST (KKA, KKU, KKL, KRR, KMI, HRAD,                    &
    
                                 HTURBDIM, OSUBG_COND, OSIGMAS, PTSTEP, PSIGQSAT,  &
    
                                 PRHODJ, PEXNREF, PSIGS,PMFCONV,PPABST,PZZ,        &
    
                                 PCF_MF,PRC_MF,PRI_MF,                             &
                                 PRVT, PRCT, PRVS, PRCS, PTHS, PSRCS, PCLDFR ,     &
                                 PRRT, PRRS, PRIT, PRIS, PRST,PRSS, PRGT, PRGS,    &
                                 PRHT, PRHS       )
    !     #########################################################################
    !
    !!****  *ICE_ADJUST* -  compute the ajustment of water vapor in mixed-phase 
    !!                      clouds
    !!
    !!    PURPOSE
    !!    -------
    !!    The purpose of this routine is to compute the fast microphysical sources
    !!    through a saturation ajustement procedure in case of mixed-phase clouds.
    !!
    !!
    !!**  METHOD
    !!    ------
    !!    Langlois, Tellus, 1973 for the cloudless version.
    !!    When cloud water is taken into account, refer to book 1 of the
    !!    documentation.
    !!
    !!     
    !!
    !!    EXTERNAL
    !!    --------
    !!      None
    !!     
    !!
    !!    IMPLICIT ARGUMENTS
    !!    ------------------
    !!      Module MODD_CST
    !!         XP00               ! Reference pressure
    !!         XMD,XMV            ! Molar mass of dry air and molar mass of vapor
    !!         XRD,XRV            ! Gaz constant for dry air, gaz constant for vapor
    !!         XCPD,XCPV          ! Cpd (dry air), Cpv (vapor)
    !!         XCL                ! Cl (liquid)
    !!         XCI                ! Ci (ice)
    !!         XTT                ! Triple point temperature
    !!         XLVTT              ! Vaporization heat constant
    !!         XLSTT              ! Sublimation  heat constant
    !!         XALPW,XBETAW,XGAMW ! Constants for saturation vapor over liquid
    !!                            !  pressure  function 
    !!         XALPI,XBETAI,XGAMI ! Constants for saturation vapor over ice
    !!                            !  pressure  function 
    !!      Module  MODD_CONF 
    !!         CCONF
    !!      Module MODD_BUDGET:
    !!         NBUMOD 
    !!         CBUTYPE
    !!         NBUPROCCTR 
    !!         LBU_RTH    
    !!         LBU_RRV  
    !!         LBU_RRC  
    !!         LBU_RRI  
    !!
    !!
    !!    REFERENCE
    !!    ---------
    !!      Book 1 and Book2 of documentation ( routine ICE_ADJUST )
    !!      Langlois, Tellus, 1973
    !!
    !!    AUTHOR
    !!    ------
    !!      J.-P. Pinty    * Laboratoire d'Aerologie*
    !!   
    !!
    !!    MODIFICATIONS
    !!    -------------
    !!      Original    06/12/96 
    !!      M. Tomasini 27/11/00 Change CND and DEP fct of the T instead of rc and ri
    !!                           Avoid the sub- and super-saturation before the ajustment
    !!                           Avoid rc>0 if T<T00 before the ajustment
    !!      P Bechtold 12/02/02  change subgrid condensation
    !!      JP Pinty   29/11/02  add ICE2 and IC4 cases
    !!      (P. Jabouille) 27/05/04 safety test for case where esw/i(T)> pabs (~Z>40km)
    !!      J.Pergaud and S.Malardel Add EDKF case
    
    !!      (E.Perraud)  06/08   add correction to avoid ice when T >0
    
    !!      S. Riette ice for EDKF
    !!      2012-02 Y. Seity,  add possibility to run with reversed vertical levels
    
    !!      J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 
    
    !!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
    
    !!
    !-------------------------------------------------------------------------------
    !
    !*       0.    DECLARATIONS
    !              ------------
    !
    USE MODD_PARAMETERS
    USE MODD_CST
    USE MODD_CONF
    USE MODD_BUDGET
    !
    USE MODI_CONDENSATION
    USE MODI_BUDGET
    USE MODE_FMWRIT
    
    !
    IMPLICIT NONE
    !
    !
    !*       0.1   Declarations of dummy arguments :
    !
    !
    INTEGER,                  INTENT(IN)    :: KKA  !near ground array index  
    INTEGER,                  INTENT(IN)    :: KKU  !uppest atmosphere array index
    INTEGER,                  INTENT(IN)    :: KKL  !vert. levels type 1=MNH -1=ARO
    INTEGER,                  INTENT(IN)    :: KRR      ! Number of moist variables
    INTEGER,                  INTENT(IN)    :: KMI      ! Model index 
    CHARACTER*4,              INTENT(IN)    :: HTURBDIM ! Dimensionality of the
                                                        ! turbulence scheme
    CHARACTER*4,              INTENT(IN)    :: HRAD     ! Radiation scheme name
    LOGICAL,                  INTENT(IN)    :: OSUBG_COND ! Switch for Subgrid 
                                                        ! Condensation
    LOGICAL                                 :: OSIGMAS  ! Switch for Sigma_s: 
                                                        ! use values computed in CONDENSATION
                                                        ! or that from turbulence scheme
    REAL,                     INTENT(IN)   :: PTSTEP    ! Double Time step
                                                        ! (single if cold start)
    REAL,                     INTENT(IN)   :: PSIGQSAT  ! coeff applied to qsat variance contribution
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRHODJ  ! Dry density * Jacobian
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PEXNREF ! Reference Exner function
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PSIGS   ! Sigma_s at time t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PMFCONV ! convective mass flux
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PPABST  ! Absolute Pressure at t        
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PZZ     ! height of model layer
    !
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PCF_MF! Convective Mass Flux Cloud fraction 
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRC_MF! Convective Mass Flux liquid mixing ratio
    REAL, DIMENSION(:,:,:),     INTENT(IN)    :: PRI_MF! Convective Mass Flux ice mixing ratio
    !
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRVT    ! Water vapor m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(IN)   ::  PRCT    ! Cloud water m.r. at t
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRCS    ! Cloud water m.r. source
    REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
    REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PSRCS   ! Second-order flux
                                                       ! s'rc'/2Sigma_s2 at time t+1
                                                       ! multiplied by Lambda_3
    REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PCLDFR  ! Cloud fraction          
    !
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRRS ! Rain water m.r. at t+1
    REAL, DIMENSION(:,:,:),  INTENT(INOUT)::  PRIS ! Cloud ice  m.r. at t+1
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRSS ! Aggregate  m.r. at t+1
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRGS ! Graupel    m.r. at t+1
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRRT ! Rain water m.r. at t
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRIT ! Cloud ice  m.r. at t
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRST ! Aggregate  m.r. at t
    REAL, DIMENSION(:,:,:),  INTENT(IN)   ::  PRGT ! Graupel    m.r. at t
    REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRHT ! Hail       m.r. at t
    REAL, DIMENSION(:,:,:), OPTIONAL, INTENT(IN)   ::  PRHS ! Hail       m.r. at t+1
    !
    !*       0.2   Declarations of local variables :
    !
    !
    REAL  :: ZEPS  ! Mv/Md
    REAL  :: ZT00,ZT0   ! Min and max temperature for the mixed phase liquid and solid water
                        ! for the coeff CND of the barycentric mixing ratio
    REAL, DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3)) &
                             :: ZEXNS,&  ! guess of the Exner functon at t+1
                                ZT,   &  ! guess of the temperature at t+1
                                ZCPH, &  ! guess of the CPh for the mixing
                                ZLV,  &  ! guess of the Lv at t+1
                                ZLS,  &  ! guess of the Ls at t+1
          ZW1,ZW2,ZW3,ZW4,ZW5,ZW6,ZW7,&  ! Work arrays for intermediate fields
                                ZCND     ! CND=(T-T00)/(T0-T00) cf sc doc and TAO etal (89)
    !
    INTEGER             :: IIU,IJU,IKU! dimensions of dummy arrays
    INTEGER             :: IIB,IJB    ! Horz index values of the first inner mass points
    INTEGER             :: IIE,IJE    ! Horz index values of the last inner mass points
    INTEGER             :: IKB        ! K index value of the first inner mass point
    INTEGER             :: IKE        ! K index value of the last inner mass point
    INTEGER             :: JITER,ITERMAX ! iterative loop for first order adjustment
    !
    LOGICAL             :: LPRETREATMENT, LNEW_ADJUST
    !
    !-------------------------------------------------------------------------------
    !
    !*       1.     PRELIMINARIES
    !               -------------
    !
    IIU = SIZE(PEXNREF,1)
    IJU = SIZE(PEXNREF,2)
    IKU = SIZE(PEXNREF,3)
    
    CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
    
    IKB=KKA+JPVEXT*KKL
    IKE=KKU-JPVEXT*KKL
    !
    ZEPS= XMV / XMD
    !
    ITERMAX=1
    !
    LPRETREATMENT=.TRUE.     ! FALSE to retreive the previous MASDEV4_1 version
    LNEW_ADJUST  =.TRUE.     ! FALSE to retreive the previous MASDEV4_1 version
    ZT0  = XTT               ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
    ZT00 = XTT-40.           ! Usefull if LPRETREATMENT=T or LNEW_ADJUST=T
    !
    !-------------------------------------------------------------------------------
    !
    !*       2.     COMPUTE QUANTITIES WITH THE GUESS OF THE FUTURE INSTANT
    !               -------------------------------------------------------
    !
    !*       2.1    estimate the pressure at t+1
    !
    
    ZEXNS(:,:,:) = ( PPABST(:,:,:) /XP00)**(XRD/XCPD)
    
    !
    !    beginning of the iterative loop
    !
    DO JITER =1,ITERMAX
    !
    !*       2.2    compute the intermediate temperature at t+1, T*
    !  
      ZT(:,:,:) = ( PTHS(:,:,:) * PTSTEP ) * ZEXNS(:,:,:)
    !
    !*       2.3    compute the latent heat of vaporization Lv(T*) at t+1
    !                   and the latent heat of sublimation  Ls(T*) at t+1
    !
      ZLV(:,:,:) = XLVTT + ( XCPV - XCL ) * ( ZT(:,:,:) -XTT )
      ZLS(:,:,:) = XLSTT + ( XCPV - XCI ) * ( ZT(:,:,:) -XTT )
    !
    !*       2.4    compute the specific heat for moist air (Cph) at t+1
    !
    IF     ( KRR == 7 ) THEN
      ZCPH(:,:,:) = XCPD + XCPV  *PTSTEP*   PRVS(:,:,:)                             &
                         + XCL   *PTSTEP* ( PRCS(:,:,:) + PRRS(:,:,:) )             &
                         + XCI   *PTSTEP* ( PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:) &
                                                                   + PRHS(:,:,:) )
    ELSE IF( KRR == 6 ) THEN
      ZCPH(:,:,:) = XCPD + XCPV  *PTSTEP*   PRVS(:,:,:)                             &
                         + XCL   *PTSTEP* ( PRCS(:,:,:) + PRRS(:,:,:) )             &
                         + XCI   *PTSTEP* ( PRIS(:,:,:) + PRSS(:,:,:) + PRGS(:,:,:) )
    ELSE IF( KRR == 5 ) THEN
      ZCPH(:,:,:) = XCPD + XCPV  *PTSTEP*   PRVS(:,:,:)                             &
                         + XCL   *PTSTEP* ( PRCS(:,:,:) + PRRS(:,:,:) )             &
                         + XCI   *PTSTEP* ( PRIS(:,:,:) + PRSS(:,:,:) )
    ELSE IF( KRR == 3 ) THEN
      ZCPH(:,:,:) = XCPD + XCPV  *PTSTEP*   PRVS(:,:,:)               &
                         + XCL   *PTSTEP* ( PRCS(:,:,:) + PRRS(:,:,:) )
    ELSE IF( KRR == 2 ) THEN
      ZCPH(:,:,:) = XCPD + XCPV  *PTSTEP*   PRVS(:,:,:) &
                         + XCL   *PTSTEP*   PRCS(:,:,:)
    END IF
    !
    !
    !*       3.     FIRST ORDER SUBGRID CONDENSATION SCHEME
    !               ---------------------------------------
    !
      IF ( OSUBG_COND ) THEN
    !
    !
    !*       3.1    compute condensate, cloud fraction
    !
    
      !   ZW3=water vapor    ZW1=rc (INOUT)  ZW2=ri (INOUT)   PSRC= s'rci'/Sigma_s^2
        ZW3= PRVS*PTSTEP;     ZW1=PRCS*PTSTEP;  ZW2=PRIS*PTSTEP
    
        CALL CONDENSATION( IIU, IJU, IKU, IIB, IIE, IJB, IJE, IKB, IKE, KKL,               &
           PPABST, PZZ, ZT, ZW3, ZW1, ZW2, PSIGS, PMFCONV, PCLDFR, PSRCS, .TRUE., OSIGMAS, &
           PSIGQSAT )
    !
    !*       3.2    compute the variation of mixing ratio
    !
                                                          !         Rc - Rc*
        ZW1(:,:,:) = (ZW1(:,:,:)/PTSTEP) - PRCS(:,:,:)       ! Pcon = ----------
                                                          !         2 Delta t
    
        ZW2(:,:,:) = (ZW2(:,:,:)/PTSTEP) - PRIS(:,:,:)       ! idem ZW1 but for Ri
    
      ELSE
    !
    !
    !*       4.     SECOND ORDER ALL OR NOTHING CONDENSATION SCHEME
    !                            FOR MIXED-PHASE CLOUD
    !               -----------------------------------------------
    !
    !
    !*       4.1    Eventually pretreatment
    !
        IF (LPRETREATMENT) THEN
    !
    !*    compute the saturation vapor pressures at t+1
    !
    
          ZW1(:,:,:) = EXP( XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG(ZT(:,:,:)) ) ! e_sw
          ZW2(:,:,:) = EXP( XALPI - XBETAI/ZT(:,:,:) - XGAMI*ALOG(ZT(:,:,:)) ) ! e_si
          ZW1(:,:,:) =  MIN(PPABST(:,:,:)/2.,ZW1(:,:,:))   ! safety limitation
          ZW2(:,:,:) =  MIN(PPABST(:,:,:)/2.,ZW2(:,:,:))   ! safety limitation
    !
    !*    compute the saturation mixing ratios at t+1
    !
          ZW3(:,:,:) =  ZW1(:,:,:) * ZEPS     &
    
                     / (  PPABST(:,:,:)  - ZW1(:,:,:) )     ! r_sw
    
          ZW4(:,:,:) =  ZW2(:,:,:) * ZEPS     &
    
                     / (  PPABST(:,:,:)  - ZW2(:,:,:) )     ! r_si
    
    !
          WHERE(PRVS(:,:,:)*PTSTEP.LT.ZW4(:,:,:).AND. PRCS(:,:,:).GT.0..AND. ZT(:,:,:).LT.XTT)
    !
    !       Subsaturation case with respect to rsi(T,P) (and case rv<0):
    !       Evaporation of rc>0 (while enough) to decrease the lack of vapor 
    !
            ZW5 (:,:,:)= MIN( PRCS , ZW4(:,:,:)/PTSTEP - PRVS(:,:,:) )                ! RVCNDC
            PRVS(:,:,:)= PRVS(:,:,:) + ZW5(:,:,:) 
            PRCS(:,:,:)= PRCS(:,:,:) - ZW5(:,:,:) 
            PTHS(:,:,:)= PTHS(:,:,:) - ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
    !
          END WHERE
    !  
          WHERE ( PRVS(:,:,:)*PTSTEP .GT. ZW3(:,:,:) )
    !
    !       Supersaturation case with respect to rsw(T,P):
    !       Condensation of the vapor that is left 
    !
            ZW5 (:,:,:)= PRVS(:,:,:) - ZW3(:,:,:)/PTSTEP 
            PRVS(:,:,:)= PRVS(:,:,:) - ZW5(:,:,:)                                  ! RVCNDC 
            PRCS(:,:,:)= PRCS(:,:,:) + ZW5(:,:,:) 
            PTHS(:,:,:)= PTHS(:,:,:) + ZW5(:,:,:) * ZLV(:,:,:) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
    !
          END WHERE
    !
          WHERE ( PRCS(:,:,:).GT.0. .AND. ZT(:,:,:).LT.ZT00 )
    !
    !       Treatment of rc>0 if T<T00:
    !
            PRIS(:,:,:)= PRIS(:,:,:) + PRCS(:,:,:) 
            PTHS(:,:,:)= PTHS(:,:,:) +                                                     &
                         PRCS(:,:,:) * (ZLS(:,:,:)-ZLV(:,:,:)) /(ZCPH(:,:,:)*PEXNREF(:,:,:))
            PRCS(:,:,:)= 0. 
    !
          END WHERE
    !
    !*       4.2    compute the intermediate temperature at t+1, T*
    !  
          ZT(:,:,:) = ( PTHS(:,:,:) * PTSTEP ) * ZEXNS(:,:,:)
    !
        END IF  !end PRETREATMENT 
    !
    !*       4.3    compute the saturation vapor pressures at t+1
    !
        ZW1(:,:,:) = EXP( XALPW - XBETAW/ZT(:,:,:) - XGAMW*ALOG(ZT(:,:,:)) ) ! e_sw
        ZW2(:,:,:) = EXP( XALPI - XBETAI/ZT(:,:,:) - XGAMI*ALOG(ZT(:,:,:)) ) ! e_si
        ZW1(:,:,:) =  MIN(PPABST(:,:,:)/2.,ZW1(:,:,:))   ! safety limitation
        ZW2(:,:,:) =  MIN(PPABST(:,:,:)/2.,ZW2(:,:,:))   ! safety limitation
    !
    !*       4.4    compute the saturation mixing ratios at t+1
    !
        ZW3(:,:,:) =  ZW1(:,:,:) * ZEPS     &
    
                   / ( PPABST(:,:,:) - ZW1(:,:,:) ) ! r_sw
    
        ZW4(:,:,:) =  ZW2(:,:,:) * ZEPS     &
    
                   / ( PPABST(:,:,:)  - ZW2(:,:,:) ) ! r_si
    
    !
    !*       4.5    compute the saturation mixing ratio derivatives (r'_vs)
    !
        ZW1(:,:,:) = (( XBETAW/ZT(:,:,:) - XGAMW ) / ZT(:,:,:)) & ! r'_sw
                        * ZW3(:,:,:) * ( 1. + ZW3(:,:,:)/ZEPS )
        ZW2(:,:,:) = (( XBETAI/ZT(:,:,:) - XGAMI ) / ZT(:,:,:)) & ! r'_si
                        * ZW4(:,:,:) * ( 1. + ZW4(:,:,:)/ZEPS )
    !
        IF (LNEW_ADJUST) THEN
          ZCND(:,:,:)= ( ZT(:,:,:) - ZT00 ) / (ZT0-ZT00)       ! Like Tao et al 89
          ZCND(:,:,:)= MAX ( MIN(ZCND(:,:,:),1.) , 0. )
        ELSE
          WHERE( (PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20 )      &
          ZCND(:,:,:)= PRCS(:,:,:) / (PRCS(:,:,:)+PRIS(:,:,:)) ! Like the original version
        END IF
    !
    !*       4.5    compute L_v CND + L_s DEP and F'(T)
    !
        WHERE( (PRCS(:,:,:)+PRIS(:,:,:)) .GT. 1.0E-20 )
    !
          ZW5(:,:,:) = ZLS(:,:,:) + ( ZLV(:,:,:)-ZLS(:,:,:) ) * ZCND(:,:,:)
          ZW6(:,:,:) = ZCPH(:,:,:) * ( PRCS(:,:,:) + PRIS(:,:,:) ) +               &
                       ZW5(:,:,:)  * ( PRCS(:,:,:)*ZW1(:,:,:)                      &
                                                   + PRIS(:,:,:)*ZW2(:,:,:) )
    !
    !*       4.6    compute Delta 2
    !
          ZW7(:,:,:) = (ZW5(:,:,:)/(ZW6(:,:,:)*ZT(:,:,:))) *                       &
                       ( PRCS(:,:,:)*ZW1(:,:,:) *                                  &
                ((-2.*XBETAW+XGAMW*ZT(:,:,:)) / (XBETAW-XGAMW*ZT(:,:,:))           &
                   + (XBETAW-XGAMW*ZT(:,:,:))*(1.0+2.0*ZW3(:,:,:)/ZEPS)/ZT(:,:,:)) &
                       + PRIS(:,:,:)*ZW2(:,:,:) *                                  &
                ((-2.*XBETAI+XGAMI*ZT(:,:,:)) / (XBETAI-XGAMI*ZT(:,:,:))           &
                   + (XBETAI-XGAMI*ZT(:,:,:))*(1.0+2.0*ZW4(:,:,:)/ZEPS)/ZT(:,:,:)) )
    !
    !*       4.7    compute Delta 1
    !
          ZW6(:,:,:) = ZW5(:,:,:)*( PRCS(:,:,:)*ZW3(:,:,:) + PRIS(:,:,:)*ZW4(:,:,:)&
                                  - PRVS(:,:,:)*PTSTEP*(PRCS(:,:,:) + PRIS(:,:,:))    &
                                  )/ZW6(:,:,:)
    !
    !*       4.8    compute the sources
    !
          ZW3(:,:,:) = (ZCPH(:,:,:)/ZW5(:,:,:)) *                               &
                       (-ZW6(:,:,:) * ( 1.0 + 0.5*ZW6(:,:,:)*ZW7(:,:,:) )) / PTSTEP
          ZW1(:,:,:) = ZW3(:,:,:)*ZCND(:,:,:)                               ! RVCNDC
          ZW2(:,:,:) = ZW3(:,:,:)*( 1.0 - ZCND(:,:,:) )                       ! RVDEPI
    !
        ELSEWHERE
    !
    !*       4.9    special case when both r_c and r_i are zero
    !
          ZW6(:,:,:) = ZCPH(:,:,:) + ZLV(:,:,:) * ZW1(:,:,:)               ! F'(T)
          ZW7(:,:,:) = (ZLV(:,:,:)/(ZW6(:,:,:)*ZT(:,:,:))) * &             ! Delta 2
                       ( ZW1(:,:,:) *                                              &
              ( (-2.*XBETAW+XGAMW*ZT(:,:,:)) / (XBETAW-XGAMW*ZT(:,:,:))            &
                  + (XBETAW-XGAMW*ZT(:,:,:))*(1.0+2.0*ZW3(:,:,:)/ZEPS)/ZT(:,:,:) ) )
          ZW6(:,:,:) = ZLV(:,:,:)*( ZW3(:,:,:)-PRVS(:,:,:)*PTSTEP )/ZW6(:,:,:)! Delta 1
          ZW1(:,:,:) = (ZCPH(:,:,:)/ZLV(:,:,:)) *                          &! RVCNDC
                       (-ZW6(:,:,:) * ( 1.0 + 0.5*ZW6(:,:,:)*ZW7(:,:,:) ))/PTSTEP
          ZW2(:,:,:) = 0.0                                                  ! RVDEPI
    !
        END WHERE
    !
      END IF 
    !
    !*       5.     COMPUTE THE SOURCES AND STORES THE CLOUD FRACTION
    !               -------------------------------------------------
    !
    !
    !*       5.1    compute the sources 
    ! 
      WHERE( ZW1(:,:,:) < 0.0 )
        ZW1(:,:,:) = MAX ( ZW1(:,:,:), -PRCS(:,:,:) )
      ELSEWHERE
        ZW1(:,:,:) = MIN ( ZW1(:,:,:),  PRVS(:,:,:) )
      END WHERE
    !
      PRVS(:,:,:) = PRVS(:,:,:) - ZW1(:,:,:)
      PRCS(:,:,:) = PRCS(:,:,:) + ZW1(:,:,:)  
      PTHS(:,:,:) = PTHS(:,:,:) +        &
                      ZW1(:,:,:) * ZLV(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
    !
      WHERE( ZW2(:,:,:) < 0.0 )
        ZW2(:,:,:) = MAX ( ZW2(:,:,:), -PRIS(:,:,:) )
      ELSEWHERE
        ZW2(:,:,:) = MIN ( ZW2(:,:,:),  PRVS(:,:,:) )
      END WHERE
    !
      PRVS(:,:,:) = PRVS(:,:,:) - ZW2(:,:,:)
      PRIS(:,:,:) = PRIS(:,:,:) + ZW2(:,:,:)
      PTHS(:,:,:) = PTHS(:,:,:) +        &
                    ZW2(:,:,:) * ZLS(:,:,:) / (ZCPH(:,:,:) * PEXNREF(:,:,:))
    !
    !
    END DO          ! end of the iterative loop
    !
    !
    !*       5.2    compute the cloud fraction PCLDFR
    !
    IF ( .NOT. OSUBG_COND ) THEN
      WHERE (PRCS(:,:,:) + PRIS(:,:,:) > 1.E-12 / PTSTEP)
        PCLDFR(:,:,:)  = 1.
      ELSEWHERE
        PCLDFR(:,:,:)  = 0. 
      ENDWHERE 
      IF ( SIZE(PSRCS,3) /= 0 ) THEN
        PSRCS(:,:,:) = PCLDFR(:,:,:) 
      END IF
    ELSE
      !We limit PRC_MF+PRI_MF to PRVS*PTSTEP to avoid negative humidity
      ZW1(:,:,:)=PRC_MF(:,:,:)
      ZW2(:,:,:)=PRI_MF(:,:,:)
      WHERE(ZW1(:,:,:)+ZW2(:,:,:)>PRVS(:,:,:)*PTSTEP)
        ZW1(:,:,:)=ZW1(:,:,:)*PRVS(:,:,:)*PTSTEP/(ZW1(:,:,:)+ZW2(:,:,:))
        ZW2(:,:,:)=PRVS(:,:,:)*PTSTEP-ZW1(:,:,:)
      ENDWHERE
      PCLDFR(:,:,:)=MIN(1.,PCLDFR(:,:,:)+PCF_MF(:,:,:))
      PRCS(:,:,:)=PRCS(:,:,:)+ZW1(:,:,:)/PTSTEP
      PRIS(:,:,:)=PRIS(:,:,:)+ZW2(:,:,:)/PTSTEP
      PRVS(:,:,:)=PRVS(:,:,:)-(ZW1(:,:,:)+ZW2(:,:,:))/PTSTEP
      PTHS(:,:,:) = PTHS(:,:,:) + &
                    (ZW1 * ZLV(:,:,:) + ZW2 * ZLS(:,:,:)) / ZCPH(:,:,:)     &
                    /  PEXNREF(:,:,:) /PTSTEP
    ENDIF
    !
    !
    !
    !*       6.  STORE THE BUDGET TERMS
    !            ----------------------
    !
    IF (LBUDGET_RV) CALL BUDGET (PRVS(:,:,:) * PRHODJ(:,:,:),6,'DEPI_BU_RRV')
    IF (LBUDGET_RC) CALL BUDGET (PRCS(:,:,:) * PRHODJ(:,:,:),7,'DEPI_BU_RRC')
    IF (LBUDGET_RI) CALL BUDGET (PRIS(:,:,:) * PRHODJ(:,:,:),9,'DEPI_BU_RRI')
    IF (LBUDGET_TH) CALL BUDGET (PTHS(:,:,:) * PRHODJ(:,:,:),4,'DEPI_BU_RTH')
    !
    !------------------------------------------------------------------------------
    !
    !
    END SUBROUTINE ICE_ADJUST