!MNH_LIC Copyright 1994-2021 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 MODE_ICE4_FAST_RH
IMPLICIT NONE
CONTAINS
SUBROUTINE ICE4_FAST_RH(KPROMA,KSIZE, LDSOFT, PCOMPUTE, PWETG, &
                       &PRHODREF, PLVFACT, PLSFACT, PPRES, &
                       &PDV, PKA, PCJ, &
                       &PLBDAS, PLBDAG, PLBDAR, PLBDAH, &
                       &PT,  PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PRHT, &
                       &PRCWETH, PRIWETH, PRSWETH, PRGWETH, PRRWETH, &
                       &PRCDRYH, PRIDRYH, PRSDRYH, PRRDRYH, PRGDRYH, PRDRYHG, PRHMLTR, &
                       &PRH_TEND, &
                       &PA_TH, PA_RC, PA_RR, PA_RI, PA_RS, PA_RG, PA_RH)
!!
!!**  PURPOSE
!!    -------
!!      Computes the fast rh process
!!
!!    AUTHOR
!!    ------
!!      S. Riette from the splitting of rain_ice source code (nov. 2014)
!!
!!    MODIFICATIONS
!!    -------------
!!
!  P. Wautelet 26/04/2019: replace non-standard FLOAT function by REAL function
!  P. Wautelet 29/05/2019: remove PACK/UNPACK intrinsics (to get more performance and better OpenACC support)
!!     R. El Khatib 24-Aug-2021 Optimizations
!
!
!*      0. DECLARATIONS
!          ------------
!
USE MODD_CST,            ONLY: XALPI, XALPW, XBETAI, XBETAW, XGAMW, XCI, XCL, XCPV, XESTT, XGAMI, XLMTT, &
                             & XLVTT, XMD, XMV, XRV, XTT, XEPSILO
USE MODD_PARAM_ICE,      ONLY: LCONVHG, LEVLIMIT, LNULLWETH, LWETHPOST
USE MODD_RAIN_ICE_DESCR, ONLY: XBG, XBS, XCEXVT, XCXG, XCXH, XCXS, XDH, XRTMIN
USE MODD_RAIN_ICE_PARAM, ONLY: NWETLBDAG, NWETLBDAH, NWETLBDAR, NWETLBDAS, X0DEPH, X1DEPH, XCOLEXGH, XCOLEXIH, &
                             & XCOLGH, XCOLIH, XCOLEXSH, XCOLSH, XEX0DEPH, XEX1DEPH, XFGWETH, XFRWETH, &
                             & XFSWETH, XFWETH, XKER_GWETH, XKER_RWETH, XKER_SWETH, XLBGWETH1, XLBGWETH2, &
                             & XLBGWETH3, XLBRWETH1, XLBRWETH2, XLBRWETH3, XLBSWETH1, XLBSWETH2, XLBSWETH3, &
                             & XWETINTP1G, XWETINTP1H, XWETINTP1R, XWETINTP1S, XWETINTP2G, XWETINTP2H, &
                             & XWETINTP2R, XWETINTP2S
!
USE PARKIND1, ONLY : JPRB
USE YOMHOOK , ONLY : LHOOK, DR_HOOK
!
IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
INTEGER,                      INTENT(IN)    :: KPROMA,KSIZE
LOGICAL,                      INTENT(IN)    :: LDSOFT
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PCOMPUTE
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PWETG    ! 1. where graupel grows in wet mode, 0. elsewhere
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRHODREF ! Reference density
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLVFACT
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLSFACT
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PPRES    ! absolute pressure at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PDV      ! Diffusivity of water vapor in the air
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PKA      ! Thermal conductivity of the air
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PCJ      ! Function to compute the ventilation coefficient
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLBDAS   ! Slope parameter of the aggregate distribution
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLBDAG   ! Slope parameter of the graupel   distribution
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLBDAR   ! Slope parameter of the rain      distribution
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PLBDAH   ! Slope parameter of the hail      distribution
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PT       ! Temperature
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRVT     ! Water vapor m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRCT     ! Cloud water m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRRT     ! Rain m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRIT     ! Pristine ice m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRST     ! Snow/aggregate m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRGT     ! Graupel m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(IN)    :: PRHT     ! Hail m.r. at t
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRCWETH  ! Dry growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRIWETH  ! Dry growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRSWETH  ! Dry growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRGWETH  ! Dry growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRRWETH  ! Dry growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRCDRYH  ! Wet growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRIDRYH  ! Wet growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRSDRYH  ! Wet growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRRDRYH  ! Wet growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRGDRYH  ! Wet growth of hailstone
REAL, DIMENSION(KSIZE),       INTENT(OUT)   :: PRDRYHG  ! Conversion of hailstone into graupel
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PRHMLTR  ! Melting of the hailstones
REAL, DIMENSION(KPROMA, 10),  INTENT(INOUT) :: PRH_TEND ! Individual tendencies
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_TH
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RC
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RR
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RI
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RS
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RG
REAL, DIMENSION(KSIZE),       INTENT(INOUT) :: PA_RH
!
!*       0.2  declaration of local variables
!
INTEGER, PARAMETER :: IRCWETH=1, IRRWETH=2, IRIDRYH=3, IRIWETH=4, IRSDRYH=5, IRSWETH=6, IRGDRYH=7, IRGWETH=8, &
                    & IFREEZ1=9, IFREEZ2=10
LOGICAL, DIMENSION(KSIZE) :: GWET
REAL, DIMENSION(KSIZE) :: ZHAIL, ZWET, ZMASK, ZWETH, ZDRYH
INTEGER :: IHAIL, IGWET
INTEGER, DIMENSION(KSIZE) :: I1
REAL, DIMENSION(KSIZE) :: ZVEC1, ZVEC2, ZVEC3
INTEGER, DIMENSION(KSIZE) :: IVEC1, IVEC2
REAL, DIMENSION(KSIZE) :: ZZW, &
                                   ZRDRYH_INIT, ZRWETH_INIT, &
                                   ZRDRYHG
INTEGER :: JJ, JL
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('ICE4_FAST_RH',0,ZHOOK_HANDLE)
!
!
!*       7.2    compute the Wet and Dry growth of hail
!
DO JL=1, KSIZE
  ZMASK(JL)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JL))) * & ! WHERE(PRHT(:)>XRTMIN(7))
           &MAX(0., -SIGN(1., XRTMIN(2)-PRCT(JL))) * & ! WHERE(PRCT(:)>XRTMIN(2))
           &PCOMPUTE(JL)
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IRCWETH)=ZMASK(JL) * PRH_TEND(JL, IRCWETH)
  ENDDO
ELSE
  PRH_TEND(:, IRCWETH)=0.
  WHERE(ZMASK(1:KSIZE)==1.)
    ZZW(1:KSIZE) = PLBDAH(1:KSIZE)**(XCXH-XDH-2.0) * PRHODREF(1:KSIZE)**(-XCEXVT)
    PRH_TEND(1:KSIZE, IRCWETH)=XFWETH * PRCT(1:KSIZE) * ZZW(1:KSIZE)    ! RCWETH
  END WHERE
ENDIF
DO JL=1, KSIZE
  ZMASK(JL)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JL))) * & ! WHERE(PRHT(:)>XRTMIN(7))
           &MAX(0., -SIGN(1., XRTMIN(4)-PRIT(JL))) * & ! WHERE(PRIT(:)>XRTMIN(4))
           &PCOMPUTE(JL)
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IRIWETH)=ZMASK(JL) * PRH_TEND(JL, IRIWETH)
    PRH_TEND(JL, IRIDRYH)=ZMASK(JL) * PRH_TEND(JL, IRIDRYH)
  ENDDO
ELSE
  PRH_TEND(:, IRIWETH)=0.
  PRH_TEND(:, IRIDRYH)=0.
  WHERE(ZMASK(1:KSIZE)==1.)
    ZZW(1:KSIZE) = PLBDAH(1:KSIZE)**(XCXH-XDH-2.0) * PRHODREF(1:KSIZE)**(-XCEXVT)
    PRH_TEND(1:KSIZE, IRIWETH)=XFWETH * PRIT(1:KSIZE) * ZZW(1:KSIZE)   ! RIWETH
    PRH_TEND(1:KSIZE, IRIDRYH)=PRH_TEND(1:KSIZE, IRIWETH)*(XCOLIH*EXP(XCOLEXIH*(PT(1:KSIZE)-XTT)))   ! RIDRYH
  END WHERE
ENDIF

!
!*       7.2.1  accretion of aggregates on the hailstones
!
IGWET = 0
DO JJ = 1, KSIZE
  ZWET(JJ) = MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JJ))) * & ! WHERE(PRHT(:)>XRTMIN(7))
            &MAX(0., -SIGN(1., XRTMIN(5)-PRST(JJ))) * & ! WHERE(PRST(:)>XRTMIN(5))
            &PCOMPUTE(JJ)
  IF (ZWET(JJ)>0) THEN
    IGWET = IGWET + 1
    I1(IGWET) = JJ
    GWET(JJ) = .TRUE.
  ELSE
    GWET(JJ) = .FALSE.
  ENDIF
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IRSWETH)=ZWET(JL) * PRH_TEND(JL, IRSWETH)
    PRH_TEND(JL, IRSDRYH)=ZWET(JL) * PRH_TEND(JL, IRSDRYH)
  ENDDO
ELSE
  PRH_TEND(:, IRSWETH)=0.
  PRH_TEND(:, IRSDRYH)=0.
  IF(IGWET>0)THEN
    !
    !*       7.2.3  select the (PLBDAH,PLBDAS) couplet
    !
    DO JJ = 1, IGWET
      ZVEC1(JJ) = PLBDAH(I1(JJ))
      ZVEC2(JJ) = PLBDAS(I1(JJ))
    ENDDO
    !
    !*       7.2.4  find the next lower indice for the PLBDAG and for the PLBDAS
    !               in the geometrical set of (Lbda_h,Lbda_s) couplet use to
    !               tabulate the SWETH-kernel
    !
    ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAH)-0.00001,           &
                          XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + XWETINTP2H ) )
    IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) )
    ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) )
    !
    ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAS)-0.00001,           &
                          XWETINTP1S * LOG( ZVEC2(1:IGWET) ) + XWETINTP2S ) )
    IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) )
    ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) )
    !
    !*       7.2.5  perform the bilinear interpolation of the normalized
    !               SWETH-kernel
    !
    DO JJ = 1,IGWET
      ZVEC3(JJ) = (  XKER_SWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                   - XKER_SWETH(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                                 * ZVEC1(JJ) &
                  - ( XKER_SWETH(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                  - XKER_SWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                         * (ZVEC1(JJ) - 1.0)
    END DO
    ZZW(:) = 0.
    DO JJ = 1, IGWET
      ZZW(I1(JJ)) = ZVEC3(JJ)
    END DO
    !
    WHERE(GWET(1:KSIZE))
      PRH_TEND(1:KSIZE, IRSWETH)=XFSWETH*ZZW(1:KSIZE)                       & ! RSWETH
                    *( PLBDAS(1:KSIZE)**(XCXS-XBS) )*( PLBDAH(1:KSIZE)**XCXH )  &
                       *( PRHODREF(1:KSIZE)**(-XCEXVT-1.) )               &
                       *( XLBSWETH1/( PLBDAH(1:KSIZE)**2              ) + &
                          XLBSWETH2/( PLBDAH(1:KSIZE)   * PLBDAS(1:KSIZE)   ) + &
                          XLBSWETH3/(               PLBDAS(1:KSIZE)**2) )
      PRH_TEND(1:KSIZE, IRSDRYH)=PRH_TEND(1:KSIZE, IRSWETH)*(XCOLSH*EXP(XCOLEXSH*(PT(1:KSIZE)-XTT)))
    END WHERE
  ENDIF
ENDIF
!
!*       7.2.6  accretion of graupeln on the hailstones
!
IGWET = 0
DO JJ = 1, KSIZE
  ZWET(JJ)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JJ))) * & ! WHERE(PRHT(:)>XRTMIN(7))
          &MAX(0., -SIGN(1., XRTMIN(6)-PRGT(JJ))) * & ! WHERE(PRGT(:)>XRTMIN(6))
          &PCOMPUTE(JJ)
  IF (ZWET(JJ)>0) THEN
    IGWET = IGWET + 1
    I1(IGWET) = JJ
    GWET(JJ) = .TRUE.
  ELSE
    GWET(JJ) = .FALSE.
  END IF
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IRGWETH)=ZWET(JL) * PRH_TEND(JL, IRGWETH)
    PRH_TEND(JL, IRGDRYH)=ZWET(JL) * PRH_TEND(JL, IRGDRYH)
  ENDDO
ELSE
  PRH_TEND(:, IRGWETH)=0.
  PRH_TEND(:, IRGDRYH)=0.
  IF(IGWET>0)THEN
    !
    !*       7.2.8  select the (PLBDAH,PLBDAG) couplet
    !
    DO JJ = 1, IGWET
      ZVEC1(JJ) = PLBDAH(I1(JJ))
      ZVEC2(JJ) = PLBDAG(I1(JJ))
    END DO
    !
    !*       7.2.9  find the next lower indice for the PLBDAH and for the PLBDAG
    !               in the geometrical set of (Lbda_h,Lbda_g) couplet use to
    !               tabulate the GWETH-kernel
    !
    ZVEC1(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAG)-0.00001,           &
                          XWETINTP1H * LOG( ZVEC1(1:IGWET) ) + XWETINTP2H ) )
    IVEC1(1:IGWET) = INT( ZVEC1(1:IGWET) )
    ZVEC1(1:IGWET) = ZVEC1(1:IGWET) - REAL( IVEC1(1:IGWET) )
    !
    ZVEC2(1:IGWET) = MAX( 1.00001, MIN( REAL(NWETLBDAG)-0.00001,           &
                          XWETINTP1G * LOG( ZVEC2(1:IGWET) ) + XWETINTP2G ) )
    IVEC2(1:IGWET) = INT( ZVEC2(1:IGWET) )
    ZVEC2(1:IGWET) = ZVEC2(1:IGWET) - REAL( IVEC2(1:IGWET) )
    !
    !*       7.2.10 perform the bilinear interpolation of the normalized
    !               GWETH-kernel
    !
    DO JJ = 1,IGWET
      ZVEC3(JJ) = (  XKER_GWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                   - XKER_GWETH(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                                 * ZVEC1(JJ) &
                - (  XKER_GWETH(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                   - XKER_GWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                        * (ZVEC1(JJ) - 1.0)
    END DO
    ZZW(:) = 0.
    DO JJ = 1, IGWET
      ZZW(I1(JJ)) = ZVEC3(JJ)
    END DO
    !
    WHERE(GWET(1:KSIZE))
      PRH_TEND(1:KSIZE, IRGWETH)=XFGWETH*ZZW(1:KSIZE)                       & ! RGWETH
                    *( PLBDAG(1:KSIZE)**(XCXG-XBG) )*( PLBDAH(1:KSIZE)**XCXH )  &
                       *( PRHODREF(1:KSIZE)**(-XCEXVT-1.) )               &
                       *( XLBGWETH1/( PLBDAH(1:KSIZE)**2              ) + &
                          XLBGWETH2/( PLBDAH(1:KSIZE)   * PLBDAG(1:KSIZE)   ) + &
                          XLBGWETH3/(               PLBDAG(1:KSIZE)**2) )
      PRH_TEND(1:KSIZE, IRGDRYH)=PRH_TEND(1:KSIZE, IRGWETH)
    END WHERE
    !When graupel grows in wet mode, graupel is wet (!) and collection efficiency must remain the same
    WHERE(GWET(1:KSIZE) .AND. .NOT. PWETG(1:KSIZE)==1.)
      PRH_TEND(1:KSIZE, IRGDRYH)=PRH_TEND(1:KSIZE, IRGDRYH)*(XCOLGH*EXP(XCOLEXGH*(PT(1:KSIZE)-XTT)))
    END WHERE
  END IF
ENDIF
!
!*       7.2.11  accretion of raindrops on the hailstones
!
IGWET = 0
DO JJ = 1, KSIZE
  ZWET(JJ)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JJ))) * & ! WHERE(PRHT(:)>XRTMIN(7))
          &MAX(0., -SIGN(1., XRTMIN(3)-PRRT(JJ))) * & ! WHERE(PRRT(:)>XRTMIN(3))
          &PCOMPUTE(JJ)
  IF (ZWET(JJ)>0) THEN
    IGWET = IGWET + 1
    I1(IGWET) = JJ
    GWET(JJ) = .TRUE.
  ELSE
    GWET(JJ) = .FALSE.
  ENDIF
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IRRWETH)=ZWET(JL) * PRH_TEND(JL, IRRWETH)
  ENDDO
ELSE
  PRH_TEND(:, IRRWETH)=0.
  IF(IGWET>0)THEN
    !
    !*       7.2.12  select the (PLBDAH,PLBDAR) couplet
    !
    DO JJ = 1, IGWET
      ZVEC1(JJ) = PLBDAH(I1(JJ))
      ZVEC2(JJ) = PLBDAR(I1(JJ))
    ENDDO
    !
    !*       7.2.13 find the next lower indice for the PLBDAH and for the PLBDAR
    !               in the geometrical set of (Lbda_h,Lbda_r) couplet use to
    !               tabulate the RWETH-kernel
    !
    ZVEC1(1:IGWET)=MAX(1.00001, MIN( REAL(NWETLBDAH)-0.00001,           &
                          XWETINTP1H*LOG(ZVEC1(1:IGWET))+XWETINTP2H))
    IVEC1(1:IGWET)=INT(ZVEC1(1:IGWET))
    ZVEC1(1:IGWET)=ZVEC1(1:IGWET)-REAL(IVEC1(1:IGWET))
    !
    ZVEC2(1:IGWET)=MAX(1.00001, MIN( REAL(NWETLBDAR)-0.00001,           &
                          XWETINTP1R*LOG(ZVEC2(1:IGWET))+XWETINTP2R))
    IVEC2(1:IGWET)=INT(ZVEC2(1:IGWET))
    ZVEC2(1:IGWET)=ZVEC2(1:IGWET)-REAL(IVEC2(1:IGWET))
    !
    !*       7.2.14 perform the bilinear interpolation of the normalized
    !               RWETH-kernel
    !
    DO JJ=1, IGWET
      ZVEC3(JJ)= (  XKER_RWETH(IVEC1(JJ)+1,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                    - XKER_RWETH(IVEC1(JJ)+1,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                                  * ZVEC1(JJ) &
                 - (  XKER_RWETH(IVEC1(JJ)  ,IVEC2(JJ)+1)* ZVEC2(JJ)          &
                    - XKER_RWETH(IVEC1(JJ)  ,IVEC2(JJ)  )*(ZVEC2(JJ) - 1.0) ) &
                                                         *(ZVEC1(JJ) - 1.0)
    END DO
    ZZW(:) = 0.
    DO JJ = 1, IGWET
      ZZW(I1(JJ)) = ZVEC3(JJ)
    END DO
    !
    WHERE(GWET(1:KSIZE))
      PRH_TEND(1:KSIZE, IRRWETH) = XFRWETH*ZZW(1:KSIZE)                    & ! RRWETH
                        *( PLBDAR(1:KSIZE)**(-4) )*( PLBDAH(1:KSIZE)**XCXH ) &
                               *( PRHODREF(1:KSIZE)**(-XCEXVT-1.) )   &
                    *( XLBRWETH1/( PLBDAH(1:KSIZE)**2              ) + &
                       XLBRWETH2/( PLBDAH(1:KSIZE)   * PLBDAR(1:KSIZE)   ) + &
                       XLBRWETH3/(               PLBDAR(1:KSIZE)**2) )
    END WHERE
  ENDIF
ENDIF
!
DO JL=1, KSIZE
  ZRDRYH_INIT(JL)=PRH_TEND(JL, IRCWETH)+PRH_TEND(JL, IRIDRYH)+ &
                 &PRH_TEND(JL, IRSDRYH)+PRH_TEND(JL, IRRWETH)+PRH_TEND(JL, IRGDRYH)
ENDDO
!
!*       7.3    compute the Wet growth of hail
!
DO JL=1, KSIZE
  ZHAIL(JL)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JL))) * & ! WHERE(PRHT(:)>XRTMIN(7))
           &PCOMPUTE(JL)
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRH_TEND(JL, IFREEZ1)=ZHAIL(JL) * PRH_TEND(JL, IFREEZ1)
    PRH_TEND(JL, IFREEZ2)=ZHAIL(JL) * PRH_TEND(JL, IFREEZ2)
  ENDDO
ELSE
  DO JL=1, KSIZE
    PRH_TEND(JL, IFREEZ1)=PRVT(JL)*PPRES(JL)/(XEPSILO+PRVT(JL)) ! Vapor pressure
  ENDDO
  IF(LEVLIMIT) THEN
    WHERE(ZHAIL(1:KSIZE)==1.)
      PRH_TEND(1:KSIZE, IFREEZ1)=MIN(PRH_TEND(1:KSIZE, IFREEZ1), EXP(XALPI-XBETAI/PT(1:KSIZE)-XGAMI*ALOG(PT(1:KSIZE)))) ! min(ev, es_i(T))
    END WHERE
  ENDIF
  PRH_TEND(:, IFREEZ2)=0.
  WHERE(ZHAIL(1:KSIZE)==1.)
    PRH_TEND(1:KSIZE, IFREEZ1)=PKA(1:KSIZE)*(XTT-PT(1:KSIZE)) +                              &
             (PDV(1:KSIZE)*(XLVTT+(XCPV-XCL)*(PT(1:KSIZE)-XTT)) &
                           *(XESTT-PRH_TEND(1:KSIZE, IFREEZ1))/(XRV*PT(1:KSIZE))           )
    PRH_TEND(1:KSIZE, IFREEZ1)=PRH_TEND(1:KSIZE, IFREEZ1)* ( X0DEPH*       PLBDAH(1:KSIZE)**XEX0DEPH +     &
                           X1DEPH*PCJ(1:KSIZE)*PLBDAH(1:KSIZE)**XEX1DEPH )/ &
                          ( PRHODREF(1:KSIZE)*(XLMTT-XCL*(XTT-PT(1:KSIZE))) )
    PRH_TEND(1:KSIZE, IFREEZ2)=(PRHODREF(1:KSIZE)*(XLMTT+(XCI-XCL)*(XTT-PT(1:KSIZE)))   ) / &
                          ( PRHODREF(1:KSIZE)*(XLMTT-XCL*(XTT-PT(1:KSIZE))) )
  END WHERE
ENDIF
DO JL=1, KSIZE
  !We must agregate, at least, the cold species
  ZRWETH_INIT(JL)=ZHAIL(JL) * MAX(PRH_TEND(JL, IRIWETH)+PRH_TEND(JL, IRSWETH)+PRH_TEND(JL, IRGWETH), &
                                 &MAX(0., PRH_TEND(JL, IFREEZ1) + &
                                         &PRH_TEND(JL, IFREEZ2) * ( &
                     &PRH_TEND(JL, IRIWETH)+PRH_TEND(JL, IRSWETH)+PRH_TEND(JL, IRGWETH) )))
ENDDO
!
!*       7.4    Select Wet or Dry case
!
!Wet case
DO JL=1, KSIZE
  ZWETH(JL) = ZHAIL(JL) * &
            & MAX(0., SIGN(1., MAX(0., ZRDRYH_INIT(JL)-PRH_TEND(JL, IRIDRYH)-PRH_TEND(JL, IRSDRYH)-PRH_TEND(JL, IRGDRYH)) - &
                              &MAX(0., ZRWETH_INIT(JL)-PRH_TEND(JL, IRIWETH)-PRH_TEND(JL, IRSWETH)-PRH_TEND(JL, IRGWETH))))
ENDDO
IF(LNULLWETH) THEN
  DO JL=1, KSIZE
    ZWETH(JL) = ZWETH(JL) * MAX(0., -SIGN(1., -ZRDRYH_INIT(JL))) ! WHERE(ZRDRYH_INIT(:)>0.)
  ENDDO
ELSE
  DO JL=1, KSIZE
    ZWETH(JL) = ZWETH(JL) * MAX(0., -SIGN(1., -ZRWETH_INIT(JL))) ! WHERE(ZRWETH_INIT(:)>0.)
  ENDDO
ENDIF
IF(.NOT. LWETHPOST) THEN
  DO JL=1, KSIZE
    ZWETH(JL) = ZWETH(JL) * MAX(0., -SIGN(1., PT(JL)-XTT)) ! WHERE(PT(:)<XTT)
  ENDDO
ENDIF
DO JL=1, KSIZE
  ZDRYH(JL) = ZHAIL(JL) * &
            & MAX(0., -SIGN(1., PT(JL)-XTT)) * & ! WHERE(PT(:)<XTT)
            & MAX(0., -SIGN(1., 1.E-20-ZRDRYH_INIT(JL))) * & !WHERE(ZRDRYH_INIT(:)>0.)
            & MAX(0., -SIGN(1., MAX(0., ZRDRYH_INIT(JL)-PRH_TEND(JL, IRIDRYH)-PRH_TEND(JL, IRSDRYH)) - &
                               &MAX(0., ZRWETH_INIT(JL)-PRH_TEND(JL, IRIWETH)-PRH_TEND(JL, IRSWETH))))
ENDDO
!
ZRDRYHG(:)=0.
IF(LCONVHG)THEN
  WHERE(ZDRYH(:)==1.)
    ZRDRYHG(:)=ZRDRYH_INIT(:)*ZRWETH_INIT(:)/(ZRDRYH_INIT(:)+ZRWETH_INIT(:))
  END WHERE
ENDIF
DO JL=1, KSIZE
  PRCWETH(JL) = ZWETH(JL) * PRH_TEND(JL, IRCWETH)
  PRIWETH(JL) = ZWETH(JL) * PRH_TEND(JL, IRIWETH)
  PRSWETH(JL) = ZWETH(JL) * PRH_TEND(JL, IRSWETH)
  PRGWETH(JL) = ZWETH(JL) * PRH_TEND(JL, IRGWETH)
  !Collected minus aggregated
  PRRWETH(JL) = ZWETH(JL) * (ZRWETH_INIT(JL) - PRH_TEND(JL, IRIWETH) - &
                             PRH_TEND(JL, IRSWETH) - PRH_TEND(JL, IRGWETH) - &
                             PRH_TEND(JL, IRCWETH))

  PRCDRYH(JL) = ZDRYH(JL) * PRH_TEND(JL, IRCWETH)
  PRIDRYH(JL) = ZDRYH(JL) * PRH_TEND(JL, IRIDRYH)
  PRSDRYH(JL) = ZDRYH(JL) * PRH_TEND(JL, IRSDRYH)
  PRRDRYH(JL) = ZDRYH(JL) * PRH_TEND(JL, IRRWETH)
  PRGDRYH(JL) = ZDRYH(JL) * PRH_TEND(JL, IRGDRYH)
  PRDRYHG(JL) = ZDRYH(JL) * ZRDRYHG(JL)

  PA_RC(JL) = PA_RC(JL) - PRCWETH(JL)
  PA_RI(JL) = PA_RI(JL) - PRIWETH(JL)
  PA_RS(JL) = PA_RS(JL) - PRSWETH(JL)
  PA_RG(JL) = PA_RG(JL) - PRGWETH(JL)
  PA_RH(JL) = PA_RH(JL) + PRCWETH(JL)+PRIWETH(JL)+PRSWETH(JL)+PRGWETH(JL)+PRRWETH(JL)
  PA_RR(JL) = PA_RR(JL) - PRRWETH(JL)
  PA_TH(JL) = PA_TH(JL) + (PRRWETH(JL)+PRCWETH(JL))*(PLSFACT(JL)-PLVFACT(JL))
  PA_RC(JL) = PA_RC(JL) - PRCDRYH(JL)
  PA_RI(JL) = PA_RI(JL) - PRIDRYH(JL)
  PA_RS(JL) = PA_RS(JL) - PRSDRYH(JL)
  PA_RR(JL) = PA_RR(JL) - PRRDRYH(JL)
  PA_RG(JL) = PA_RG(JL) - PRGDRYH(JL) + PRDRYHG(JL)
  PA_RH(JL) = PA_RH(JL) + PRCDRYH(JL)+PRIDRYH(JL)+PRSDRYH(JL)+&
                         &PRRDRYH(JL)+PRGDRYH(JL) - PRDRYHG(JL)
  PA_TH(JL) = PA_TH(JL) + (PRCDRYH(JL)+PRRDRYH(JL))*(PLSFACT(JL)-PLVFACT(JL))
ENDDO
!
!*       7.5    Melting of the hailstones
!
DO JL=1, KSIZE
  ZMASK(JL)=MAX(0., -SIGN(1., XRTMIN(7)-PRHT(JL))) * & ! WHERE(PRHT(:)>XRTMIN(7))
           &MAX(0., -SIGN(1., XTT-PT(JL))) * & ! WHERE(PT(:)>XTT)
           &PCOMPUTE(JL)
ENDDO
IF(LDSOFT) THEN
  DO JL=1, KSIZE
    PRHMLTR(JL)=ZMASK(JL)*PRHMLTR(JL)
  ENDDO
ELSE
  DO JL=1, KSIZE
    PRHMLTR(JL) = ZMASK(JL)* PRVT(JL)*PPRES(JL)/(XEPSILO+PRVT(JL)) ! Vapor pressure
  ENDDO
  IF(LEVLIMIT) THEN
    WHERE(ZMASK(:)==1.)
      PRHMLTR(:)=MIN(PRHMLTR(:), EXP(XALPW-XBETAW/PT(:)-XGAMW*ALOG(PT(:)))) ! min(ev, es_w(T))
    END WHERE
  ENDIF
  DO JL=1, KSIZE
    PRHMLTR(JL) = ZMASK(JL)* (PKA(JL)*(XTT-PT(JL)) +                              &
           ( PDV(JL)*(XLVTT + ( XCPV - XCL ) * ( PT(JL) - XTT )) &
                           *(XESTT-PRHMLTR(JL))/(XRV*PT(JL))         ))
  ENDDO
  WHERE(ZMASK(1:KSIZE)==1.)
    !
    ! compute RHMLTR
    !
    PRHMLTR(1:KSIZE)  = MAX( 0.0,( -PRHMLTR(1:KSIZE) *                     &
                           ( X0DEPH*       PLBDAH(1:KSIZE)**XEX0DEPH +     &
                             X1DEPH*PCJ(1:KSIZE)*PLBDAH(1:KSIZE)**XEX1DEPH ) -   &
                         ( PRH_TEND(1:KSIZE, IRCWETH)+PRH_TEND(1:KSIZE, IRRWETH) )*        &
                               ( PRHODREF(1:KSIZE)*XCL*(XTT-PT(1:KSIZE))) ) /    &
                                             ( PRHODREF(1:KSIZE)*XLMTT ) )
  END WHERE
END IF
DO JL=1, KSIZE
  PA_RR(JL) = PA_RR(JL) + PRHMLTR(JL)
  PA_RH(JL) = PA_RH(JL) - PRHMLTR(JL)
  PA_TH(JL) = PA_TH(JL) - PRHMLTR(JL)*(PLSFACT(JL)-PLVFACT(JL))
ENDDO
!
IF (LHOOK) CALL DR_HOOK('ICE4_FAST_RH', 1, ZHOOK_HANDLE)
!
END SUBROUTINE ICE4_FAST_RH
END MODULE MODE_ICE4_FAST_RH
