!MNH_LIC Copyright 2018-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.
!-----------------------------------------------------------------
! Modifications:
!  P. Wautelet 28/05/2018: corrected truncated integer division (1*10**(-6) -> 1E-6)
!  P. Wautelet 28/05/2019: move COUNTJV function to tools.f90
!  P. Wautelet 01/03/2021: bugfix: correct intent of PSVT in module interface
!-----------------------------------------------------------------
!      ####################
       MODULE MODI_SUBL_BLOWSNOW
!      ####################
!
INTERFACE
      SUBROUTINE SUBL_BLOWSNOW(PZZ, PRHODJ , PRHODREF, PEXNREF , PPABST,   &
                         PTHT, PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PSVT,    &
                         PTHS, PRVS, PSVS,PSNWSUBL3D,PVGK)

REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ     ! Height (z)
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ  ! Dry density * Jacobian
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODREF! Reference density
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF ! Reference Exner function
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PPABST  ! absolute pressure at t

REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT    ! Theta at time t
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(IN)    :: PRRT    ! Rain water m.r. at t
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRIT    ! Pristine ice m.r. at t
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRST    ! Snow/aggregate m.r. at t
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRGT    ! Graupel/hail m.r. at t
REAL, DIMENSION(:,:,:,:),   INTENT(INOUT)  :: PSVT    ! Blowing snow concentration
REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PVGK    ! Mass averaged settling velocity



REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
REAL, DIMENSION(:,:,:,:),   INTENT(INOUT) :: PSVS  ! Blowing snow source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PSNWSUBL3D ! Blowing snow sublimation flux (kg/m3/s)



END SUBROUTINE SUBL_BLOWSNOW
END INTERFACE
END MODULE MODI_SUBL_BLOWSNOW

      SUBROUTINE SUBL_BLOWSNOW(PZZ, PRHODJ , PRHODREF, PEXNREF , PPABST,    &
                         PTHT, PRVT, PRCT, PRRT, PRIT, PRST, PRGT, PSVT,     &
                         PTHS, PRVS, PSVS,PSNWSUBL3D,PVGK)

USE MODD_BLOWSNOW
USE MODD_CST
USE MODD_CSTS_BLOWSNOW
USE MODD_PARAMETERS

USE MODE_BLOWSNOW_PSD
use mode_tools,           only: Countjv

USE MODI_GAMMA
USE MODI_GAMMA_INC
USE MODI_GAMMA_INC_LOW

IMPLICIT NONE
!
!*       0.1   Declarations of dummy arguments :
!
!

REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PZZ     ! Height (z)
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODJ  ! Dry density * Jacobian
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRHODREF! Reference density
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PEXNREF ! Reference Exner function
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PPABST  ! absolute pressure at t

REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PTHT    ! Theta at time t
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(IN)    :: PRRT    ! Rain water m.r. at t 
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRIT    ! Pristine ice m.r. at t
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRST    ! Snow/aggregate m.r. at t
REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRGT    ! Graupel/hail m.r. at t
REAL, DIMENSION(:,:,:,:),   INTENT(INOUT)  :: PSVT    ! Drifting snow concentration at t
REAL, DIMENSION(:,:,:),   INTENT(IN)  :: PVGK    ! ! Mass averaged settling velocity


REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PTHS    ! Theta source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PRVS    ! Water vapor m.r. source
REAL, DIMENSION(:,:,:,:),   INTENT(INOUT) :: PSVS  ! Drifting snow source
REAL, DIMENSION(:,:,:),   INTENT(INOUT) :: PSNWSUBL3D ! Drifting snow sublimation flux (kg/m3/s)
!
!*       0.2   Declarations of local variables :
!
!
INTEGER :: JN            ! Loop index for numerical integration 
INTEGER :: IIB           !  Define the domain where is 
INTEGER :: IIE           !  the microphysical sources have to be computed
INTEGER :: IJB           !
INTEGER :: IJE           !
INTEGER :: IKB           !
INTEGER :: IKE           !
!
REAL,  DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3))   ::  ZBETA
REAL,  DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3))   ::  ZT
REAL,    DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3))   &
                                  :: ZW ! work array
LOGICAL,    DIMENSION(SIZE(PEXNREF,1),SIZE(PEXNREF,2),SIZE(PEXNREF,3))   &
                                  :: GSUBL ! Test where to compute sublimation

REAL, DIMENSION(:), ALLOCATABLE :: ZRVT    ! Water vapor m.r. at t
REAL, DIMENSION(:), ALLOCATABLE :: ZRCT    ! Cloud water m.r. at t 
REAL, DIMENSION(:), ALLOCATABLE :: ZRRT    ! Rain water m.r. at t 
REAL, DIMENSION(:), ALLOCATABLE :: ZRIT    ! Pristine ice m.r. at t
REAL, DIMENSION(:), ALLOCATABLE :: ZRST    ! Snow/aggregate m.r. at t
REAL, DIMENSION(:), ALLOCATABLE :: ZRGT    ! Graupel m.r. at t

REAL, DIMENSION(:), ALLOCATABLE :: ZRVS    ! Water vapor m.r. source

REAL, DIMENSION(:,:), ALLOCATABLE :: ZSVT    ! Drifting snow m.r. at t
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSVS    ! Drifting snow m.r. source


REAL, DIMENSION(:), ALLOCATABLE :: ZTHS    ! Theta source

INTEGER, DIMENSION(:), ALLOCATABLE :: NMAX    ! Maximum index for numerical integration


REAL, DIMENSION(:), ALLOCATABLE &
         ::  ZZT,      & ! Temperature
             ZRHODREF, & ! RHO Dry REFerence
             ZPRES,    & ! Pressure
             ZKA,      & ! Thermal conductivity of the air
             ZDV,      & ! Diffusivity of water vapor in the air
             ZUSI,     & ! Undersaturation over ice
             ZZW,      & ! Work array  
             ZAI,      & ! Denominator in Thorpe and Masson (66) formulation 
             ZAA,      & ! Constant in Carrier's equation for settling velocity
             ZBB,      & ! Constant in Carrier's equation for settling velocity
             ZR1,      & ! 1st limit radius in Mitchell's formulation 
             ZR2,      & ! 2nd limit radius in Mitchell's formulation 
             ZAM1,     & ! Constant in Mitchell's fall speed : v = a * R^b 
             ZAM2,     & ! Constant in Mitchell's fall speed : v = a * R^b 
             ZAM3,     & ! Constant in Mitchell's fall speed : v = a * R^b 
             ZZBETA,   & ! Scale parameter
             ZEXNREF,  & ! EXNer Pressure REFerence
             ZMU,      & ! Air kinematic viscosity
             ZZZ,      & ! Height
             ZLSFACT,  & ! L_s/(Pi_ref*C_ph)     
             ZSNWSUBL, & ! Snow sublimation rate kg.m{-3}.s{-1}
             ZVGK        ! Mass averaged settling velocity


REAL                                                       :: ZGAM,ZVEL_CARRIER,ZR,ZVEL_VENT
REAL                                                       :: ZW_M0, ZNU , ZMASS
REAL                                                       :: ZSUM_SUBL,ZNUM,ZMOB,ZTEMP
REAL                                                       :: ZDELTAR
REAL                                                       :: ZGAM_BM1,ZGAM_BM2,ZGAM_BM3
REAL                                                       :: ZR_EFF

INTEGER                           :: IMICRO
INTEGER , DIMENSION(SIZE(GSUBL)) :: I1,I2,I3 ! Used to replace the COUNT
INTEGER                           :: JL, JSV       ! and PACK intrinsics
LOGICAL :: LNONEFFICIENT
LOGICAL :: LSUBL_PIEKTUK
LOGICAL :: LSUBL_ALPINE3D
!
!              Initialize variables
ZDELTAR = 1e-6                   ! Bin size (m)
ZGAM     = GAMMA(XALPHA_SNOW)
ZGAM_BM1 = GAMMA(1.5*XBM1+XALPHA_SNOW+1)
ZGAM_BM2 = GAMMA(1.5*XBM2+XALPHA_SNOW+1)
ZGAM_BM3 = GAMMA(1.5*XBM3+XALPHA_SNOW+1)
!
LSUBL_PIEKTUK = .TRUE.  ! Compute sublimation according to PIEKTUK (Dery and Yau, 1999)
                        ! Use mass-averaged settling velocity as ventilation
                        ! velocity 
                        ! Save computational time compared to numerical
                        ! integration of Carrier's or Mitchell's formulation
!
LSUBL_ALPINE3D    = .FALSE.  ! Compute sublimation using the method of reprsentative
                        ! radius implemented in Alpine 3D (Groot and al, 2011)



!              Air Temperature
ZT(:,:,:) = PTHT(:,:,:) * ( PPABST(:,:,:) / XP00 ) ** (XRD/XCPD)

!
!-------------------------------------------------------------------------------
!
!*       1.     COMPUTE THE LOOP BOUNDS
!   	        -----------------------
!
IIB=1+JPHEXT
IIE=SIZE(PZZ,1) - JPHEXT
IJB=1+JPHEXT
IJE=SIZE(PZZ,2) - JPHEXT
IKB=1+JPVEXT
IKE=SIZE(PZZ,3) - JPVEXT
!

!
!-------------------------------------------------------------------------------
!
!*       2. COMPUTE THE BLOWINGG SNOW SUBLIMATION
!
! Optimization by looking for locations where
! the blowing snow fields are larger than a minimal value only !!!
!
!     compute parameters of the snow particle distribution
!
CALL PPP2SNOW(PSVT, PRHODREF, PBET3D=ZBETA)
!

GSUBL(:,:,:) = .FALSE.

GSUBL(IIB:IIE,IJB:IJE,IKB:IKE) =                               &
                PSVT(IIB:IIE,IJB:IJE,IKB:IKE,1)>10 .AND. &
                PSVT(IIB:IIE,IJB:IJE,IKB:IKE,2)>1e-20

!GSUBL(IIB:IIE,IJB:IJE,IKB:IKE) =                               &
!                PSVT(IIB:IIE,IJB:IJE,IKB:IKE,1)>0. .OR. &
!                PSVT(IIB:IIE,IJB:IJE,IKB:IKE,2)>0.

IMICRO = COUNTJV( GSUBL(:,:,:),I1(:),I2(:),I3(:))
IF( IMICRO >= 0 ) THEN
  ALLOCATE(ZRVT(IMICRO))
  ALLOCATE(ZRCT(IMICRO))
  ALLOCATE(ZRRT(IMICRO))
  ALLOCATE(ZRIT(IMICRO))
  ALLOCATE(ZRST(IMICRO))
  ALLOCATE(ZRGT(IMICRO)) 
  ALLOCATE(ZRVS(IMICRO))
  ALLOCATE(ZSVT(IMICRO, NBLOWSNOW3D )) 
  ALLOCATE(ZSVS(IMICRO, NBLOWSNOW3D ))
  ALLOCATE(ZTHS(IMICRO))
  ALLOCATE(ZZT(IMICRO))
  ALLOCATE(ZRHODREF(IMICRO))
  ALLOCATE(ZPRES(IMICRO))
  ALLOCATE(ZZBETA(IMICRO))
  ALLOCATE(ZEXNREF(IMICRO))
  ALLOCATE(ZZZ(IMICRO))
  ALLOCATE(ZVGK(IMICRO))
  ALLOCATE(ZSNWSUBL(IMICRO))


  DO JL=1,IMICRO
        ZRVT(JL)     = PRVT(I1(JL),I2(JL),I3(JL))
        ZRCT(JL)     = PRCT(I1(JL),I2(JL),I3(JL))
        ZRRT(JL)     = PRRT(I1(JL),I2(JL),I3(JL))
        ZRIT(JL)     = PRIT(I1(JL),I2(JL),I3(JL))
        ZRST(JL)     = PRST(I1(JL),I2(JL),I3(JL))
        ZRGT(JL)     = PRGT(I1(JL),I2(JL),I3(JL))
        ZRVS(JL)     = PRVS(I1(JL),I2(JL),I3(JL))
        ZSVT(JL,:)   = PSVT(I1(JL),I2(JL),I3(JL),:)
        ZSVS(JL,:)   = PSVS(I1(JL),I2(JL),I3(JL),:)
        ZTHS(JL)     = PTHS(I1(JL),I2(JL),I3(JL))
        ZZT(JL)      = ZT(I1(JL),I2(JL),I3(JL))
        ZRHODREF(JL) = PRHODREF(I1(JL),I2(JL),I3(JL))
        ZPRES(JL)    = PPABST(I1(JL),I2(JL),I3(JL))
        ZZBETA(JL)   = ZBETA(I1(JL),I2(JL),I3(JL))
        ZEXNREF(JL)  = PEXNREF(I1(JL),I2(JL),I3(JL))
        ZZZ(JL)      = PZZ(I1(JL),I2(JL),I3(JL))
        ZVGK(JL)   = PVGK(I1(JL),I2(JL),I3(JL))
        ZSNWSUBL(JL)   = PSNWSUBL3D(I1(JL),I2(JL),I3(JL))
   END DO
  ALLOCATE(ZZW(IMICRO))
  ALLOCATE(ZUSI(IMICRO))
  ZZW(:) = EXP( XALPI - XBETAI/ZZT(:) - XGAMI*ALOG(ZZT(:) ) )
  ZUSI(:) = ZRVT(:)*( ZPRES(:)-ZZW(:) ) / ( (XMV/XMD) * ZZW(:) ) - 1.0
                                                      ! Undersaturation over ice
 
  ALLOCATE(ZLSFACT(IMICRO))
    ZZW(:)  = ZEXNREF(:)*( XCPD+XCPV*ZRVT(:)+XCL*(ZRCT(:)+ZRRT(:)) &
                                    +XCI*(ZRIT(:)+ZRST(:)+ZRGT(:)) )
    ZLSFACT(:) = (XLSTT+(XCPV-XCI)*(ZZT(:)-XTT))/ZZW(:) ! L_s/(Pi_ref*C_ph)                                                     
  ALLOCATE(ZKA(IMICRO))
  ALLOCATE(ZDV(IMICRO))
  ALLOCATE(ZMU(IMICRO))
  ALLOCATE(ZAI(IMICRO))
  ALLOCATE(ZAA(IMICRO))
  ALLOCATE(ZBB(IMICRO))
  ALLOCATE(ZR1(IMICRO))
  ALLOCATE(ZR2(IMICRO))
  ALLOCATE(ZAM1(IMICRO))
  ALLOCATE(ZAM2(IMICRO))
  ALLOCATE(ZAM3(IMICRO))
  ALLOCATE(NMAX(IMICRO))


CALL SNOW_SUBL

  ZW(:,:,:) = PRVS(:,:,:)
  PRVS(:,:,:) = UNPACK( ZRVS(:),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )
  ZW(:,:,:) = PTHS(:,:,:)
  PTHS(:,:,:) = UNPACK( ZTHS(:),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )
  ZW(:,:,:) = PSVS(:,:,:,1)
  PSVS(:,:,:,1) = UNPACK( ZSVS(:,1),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )
  ZW(:,:,:) = PSVS(:,:,:,2)
  PSVS(:,:,:,2) = UNPACK( ZSVS(:,2),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )
!  ZW(:,:,:) = PSVS(:,:,:,3)
!  PSVS(:,:,:,3) = UNPACK( ZSVS(:,3),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )
  ZW(:,:,:) = PSNWSUBL3D(:,:,:)
  PSNWSUBL3D(:,:,:) = UNPACK( ZSNWSUBL(:),MASK=GSUBL(:,:,:),FIELD=ZW(:,:,:) )


  DEALLOCATE(ZRVT)
  DEALLOCATE(ZRCT)
  DEALLOCATE(ZRRT)
  DEALLOCATE(ZRIT)
  DEALLOCATE(ZRST)
  DEALLOCATE(ZRGT) 
  DEALLOCATE(ZRVS)
  DEALLOCATE(ZSVT) 
  DEALLOCATE(ZSVS)
  DEALLOCATE(ZTHS)
  DEALLOCATE(ZZT)
  DEALLOCATE(ZRHODREF)
  DEALLOCATE(ZPRES)
  DEALLOCATE(ZKA)
  DEALLOCATE(ZDV)
  DEALLOCATE(ZUSI)
  DEALLOCATE(ZZW)
  DEALLOCATE(ZAI)
  DEALLOCATE(ZAA)
  DEALLOCATE(ZBB)
  DEALLOCATE(ZR1)
  DEALLOCATE(ZR2)
  DEALLOCATE(ZAM1)
  DEALLOCATE(ZAM2)
  DEALLOCATE(ZAM3)
  DEALLOCATE(ZZBETA)
  DEALLOCATE(NMAX)
  DEALLOCATE(ZEXNREF)
  DEALLOCATE(ZLSFACT)
  DEALLOCATE(ZZZ)
  DEALLOCATE(ZMU)
  DEALLOCATE(ZSNWSUBL)
  DEALLOCATE(ZVGK)


END IF
!
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
!
CONTAINS

!
!-------------------------------------------------------------------------------
!

SUBROUTINE SNOW_SUBL

IMPLICIT NONE

! Sutherland's equation for kinematic viscosity
ZMU(:)=1.8325d-5*416.16/(ZZT(:)+120)*(ZZT(:)/296.16)*SQRT(ZZT(:)/296.16)/ZRHODREF(:)
! Thermal conductivity of the air
  ZKA(:) = 2.38E-2 + 0.0071E-2 * ( ZZT(:) - XTT )          ! k_a
! Diffusivity of water vapor in the air. 
  ZDV(:) = 0.211E-4 * (ZZT(:)/XTT)**1.94 * (XP00/ZPRES(:)) ! D_v
!
!*       Compute the denominator in the Thorpe and Masson (66) equation
!  
ZAI(:) = EXP( XALPI - XBETAI/ZZT(:) - XGAMI*ALOG(ZZT(:) ) ) ! es_i
ZAI(:) = ( XLSTT + (XCPV-XCI)*(ZZT(:)-XTT) ) / (ZKA(:)*ZZT(:))                  &
             *( ( XLSTT + (XCPV-XCI)*(ZZT(:)-XTT) ) / (XRV*ZZT(:)) - 1.)  &
             + (XRV*ZZT(:)) / (ZDV(:)*ZAI(:))

IF(LSUBL_ALPINE3D) THEN
         ZR_EFF = 73.5e-6  ! Effective radisus computed following the Swiss
                          ! method. This effective radius give the same total 
                          ! sublimation for a equal concentration an ensemble of
                          ! gamma distributed particles with rm = 35e-6 m and
                          ! alpha=3 
! Compute coefficient for settling velocity following Carrier (1953)
        ZAA(:) = 6.203*ZMU(:)/2.
        ZBB(:) = 5.516*XRHOLI/(4.*ZRHODREF(:))*XG
       DO JL=1,IMICRO
          ZSUM_SUBL = 0.
          ZUSI(JL) = MIN(ZUSI(JL), 0.)   !Only the undersaturation over ice is considered.
          ! Ventilation velocity taken as settling velocity of particle of mean
          ! radius
          ZVEL_VENT = - ZAA(JL)/ZR_EFF+((ZAA(JL)/ZR_EFF)**2+ZBB(JL)*ZR_EFF)**0.5
!       Nusselt Number                                                 
            ZNU    =    NUSSELT(ZR_EFF,ZMU(JL),ZVEL_VENT)
!       Rate of change of mass for a subliming ice sphere of radius ZR_EFF            
            ZMASS  =    2*XPI*ZR_EFF*ZNU*ZUSI(JL)/ZAI(JL)
!       Integration over the radius spectrum 
        ZSUM_SUBL = ZMASS*ZSVT(JL,2)/(4./3.*XPI*XRHOLI*ZR_EFF**2)

    ZSUM_SUBL =  MIN( ZRVS(JL),ZSUM_SUBL)*(0.5+SIGN(0.5,ZSUM_SUBL)) &
                   - MIN(ZSVS(JL,2),ABS(ZSUM_SUBL))*(0.5-SIGN(0.5,ZSUM_SUBL))
    ZSUM_SUBL=MIN(0.,ZSUM_SUBL) ! Sink of snow               
!   Change in concentration rate  Sn = Sb*N/qb (Dery and Yau,2000)
    ZNUM      = ZSUM_SUBL*ZSVT(JL,1)/ZSVT(JL,2)
!   Change in mobility index : value mob*M3 is reduced according to reduction of
!   M3 due to sublimation so that mob is constant due to sublimation
!    ZMOB      = ZSUM_SUBL*ZSVT(JL,3)/ZSVT(JL,2)
!   Update tendencies for snow particles, water vapour and potential temperature 
    ZSVS(JL,2) = ZSVS(JL,2) + ZSUM_SUBL             ! Particle mixing ratio
    ZSVS(JL,1) = ZSVS(JL,1) + ZNUM                  ! Particles number
!    ZSVS(JL,3) = ZSVS(JL,3) + ZMOB
    ZRVS(JL) = ZRVS(JL) - ZSUM_SUBL                 ! Water vapour
    ZTHS(JL) = ZTHS(JL) + ZSUM_SUBL*ZLSFACT(JL)     ! Potential temperature
    ZSNWSUBL(JL) = ZSNWSUBL(JL)+ZSUM_SUBL*ZRHODREF(JL)  ! Sublimation rate kg/m3/s

END DO
ELSE IF(LSUBL_PIEKTUK) THEN 
DO JL=1,IMICRO
   ZSUM_SUBL=0.
   ZUSI(JL) = MIN(ZUSI(JL), 0.)   !Only the undersaturation over ice is considered.
!       Ventilation velocity as mass averaged settling velocity            
   ZVEL_VENT = ZVGK(JL)      
!       Nusselt Number using mean radius of particle size distribution  
   ZNU    =    NUSSELT(XALPHA_SNOW*ZZBETA(JL),ZMU(JL),ZVEL_VENT)
   ! mass averaged sublimation rate follows Dery and Yan (1999) and avoids
   ! numerical integration over the particle spectrum
   ZSUM_SUBL = ZSVT(JL,2)*ZNU*ZUSI(JL)/   &
                (ZAI(JL)*2*XRHOLI*(XALPHA_SNOW*ZZBETA(JL))**2)
!   Restriction of ZSUM_SUBL   
    ZTEMP=ZSUM_SUBL
    ZSUM_SUBL =  MIN( ZRVS(JL),ZSUM_SUBL)*(0.5+SIGN(0.5,ZSUM_SUBL)) &
                   - MIN(ZSVS(JL,2),ABS(ZSUM_SUBL))*(0.5-SIGN(0.5,ZSUM_SUBL))
    ZSUM_SUBL=MIN(0.,ZSUM_SUBL) ! Sink of snow               
    IF(ZSUM_SUBL>0) THEN  
       write(*,*) 'Warning Subl',JL,'Subl',ZSUM_SUBL,'TEMP',ZTEMP
       write(*,*) 'Warning Subl ZSVT',ZSVT(JL,1),ZSVT(JL,2)
       write(*,*) 'Warning vap',ZRVS(JL),'ZSVs',ZSVS(JL,2)
    END IF
!   Change in concentration rate  Sn = Sb*N/qb (Dery and Yau,2000)
    ZNUM      = ZSUM_SUBL*ZSVT(JL,1)/ZSVT(JL,2)
!   Change in mobility index : value mob*M3 is reduced according to reduction of
!   M3 due to sublimation so that mob is constant due to sublimation
!    ZMOB      = ZSUM_SUBL*ZSVT(JL,3)/ZSVT(JL,2)
!   Update tendencies for snow particles, water vapour and potential temperature 
    ZSVS(JL,2) = ZSVS(JL,2) + ZSUM_SUBL             ! Particle mixing ratio
    ZSVS(JL,1) = ZSVS(JL,1) + ZNUM                  ! Particles number
!    ZSVS(JL,3) = ZSVS(JL,3) + ZMOB
    ZRVS(JL) = ZRVS(JL) - ZSUM_SUBL                 ! Water vapour
    ZTHS(JL) = ZTHS(JL) + ZSUM_SUBL*ZLSFACT(JL)     ! Potential temperature
    ZSNWSUBL(JL) = ZSNWSUBL(JL)+ZSUM_SUBL*ZRHODREF(JL)  ! Sublimation rate kg/m3/s

END DO

ELSE        
 !
!*       Compute the constants in Carrier equation
!
IF(CSNOWSEDIM=='CARR') THEN
ZAA(:) = 6.203*ZMU(:)/2.
ZBB(:) = 5.516*XRHOLI/(4.*ZRHODREF(:))*XG
NMAX=GET_INDEX(ZZBETA(:),ZDELTAR)

DO JL=1,IMICRO
   ZSUM_SUBL=0.
   ZUSI(JL) = MIN(ZUSI(JL), 0.)   !Only the undersaturation over ice is considered.  
   DO JN=1,NMAX(JL)
            ZR = 1E-6+(JN-0.5)*ZDELTAR
!       Carrier settling velocity
            ZVEL_CARRIER = - ZAA(JL)/ZR+((ZAA(JL)/ZR)**2+ZBB(JL)*ZR)**0.5
!       Weight of the corresponding bin following the gamma distribution            
            ZW_M0=ZSVT(JL,1)*ZR**(XALPHA_SNOW-1)*exp(-ZR/ZZBETA(JL))/(ZZBETA(JL))**XALPHA_SNOW*ZGAM
!       Ventilation velocity as a sum of settling velocity and relative
!       turbulent velocity fluctuations            
            ZVEL_VENT = ZVEL_CARRIER!+TURB_FLUC(ZR,ZMU(JL),ZVEL_CARRIER,ZRHODREF(JL),    & 
!                                                 ZZZ(JL),ZVMOD(JL))
!       Nusselt Number                                                 
            ZNU    =    NUSSELT(ZR,ZMU(JL),ZVEL_VENT)
!       Rate of change of mass for a subliming ice sphere            
            ZMASS  =    2*XPI*ZR*ZNU*ZUSI(JL)/ZAI(JL)
!       Integration over the radius spectrum            
            ZSUM_SUBL = ZSUM_SUBL+ZW_M0*ZMASS*ZDELTAR 
    END DO
!   Restriction of ZSUM_SUBL    
    ZSUM_SUBL =  MIN( ZRVS(JL),ZSUM_SUBL)*(0.5+SIGN(0.5,ZSUM_SUBL)) &
                   - MIN( ZSVS(JL,2),ABS(ZSUM_SUBL) )*(0.5-SIGN(0.5,ZSUM_SUBL))
!   Change in concentration rate  Sn = Sb*N/qb (Dery and Yau,2000)
    ZNUM      = ZSUM_SUBL*ZSVT(JL,1)/ZSVT(JL,2)
!   Change in mobility index : value mob*M3 is reduced according to reduction of
!   M3 due to sublimation so that mob is constant due to sublimation
!    ZMOB      = ZSUM_SUBL*ZSVT(JL,3)/ZSVT(JL,2)
!   Update tendencies for snow particles, water vapour and potential temperature 
    ZSVS(JL,2) = ZSVS(JL,2) + ZSUM_SUBL             ! Particle mixing ratio
    ZSVS(JL,1) = ZSVS(JL,1) + ZNUM                  ! Particles number
!    ZSVS(JL,3) = ZSVS(JL,3) + ZMOB
    ZRVS(JL) = ZRVS(JL) - ZSUM_SUBL                 ! Water vapour
    ZTHS(JL) = ZTHS(JL) + ZSUM_SUBL*ZLSFACT(JL)     ! Potential temperature
    ZSNWSUBL(JL) = ZSUM_SUBL*ZRHODREF(JL)           ! Sublimation rate kg/m3/s
END DO
END IF

IF(CSNOWSEDIM=='MITC') THEN
LNONEFFICIENT = .FALSE.
!        write(*,*) 'MITC'
  ! Compute limit radius for integration of Mitchell's formulation
ZR1(:)=RLIM(ZMU,ZRHODREF,XBESTL_1)
ZR2(:)=RLIM(ZMU,ZRHODREF,XBESTL_2)
  ! Compute parameter avr for integration of Mitchell's formulation
ZAM1(:)=AVR(XAM1,XBM1,ZRHODREF,ZMU)
ZAM2(:)=AVR(XAM2,XBM2,ZRHODREF,ZMU)
ZAM3(:)=AVR(XAM3,XBM3,ZRHODREF,ZMU)

DO JL=1,IMICRO
    ZUSI(JL) = MIN(ZUSI(JL), 0.)  !Only the undersaturation over ice is considered.
!                               no water deposition on blown snow particles
IF(LNONEFFICIENT) THEN
    ZSUM_SUBL = 2*XPI*ZUSI(JL)*ZSVT(JL,1)/ZAI(JL)*(XANU*ZZBETA(JL)*XALPHA_SNOW  +      &
                 XBNU/ZGAM*(2/ZMU(JL))**0.5*(                                          & 
                 ZZBETA(JL)**(1.5*XBM1+1)*ZAM1(JL)**0.5*ZGAM_BM1*                        &
                       GAMMA_INC(1.5*XBM1+XALPHA_SNOW+1,ZR1(JL)/ZZBETA(JL)) +            &
                 ZZBETA(JL)**(1.5*XBM2+1)*ZAM2(JL)**0.5*ZGAM_BM2*                        &
                       (GAMMA_INC(1.5*XBM2+XALPHA_SNOW+1,ZR2(JL)/ZZBETA(JL))-            &
                       GAMMA_INC(1.5*XBM2+XALPHA_SNOW+1,ZR1(JL)/ZZBETA(JL)))+            & 
                 ZZBETA(JL)**(1.5*XBM3+1)*ZAM3(JL)**0.5*ZGAM_BM3*                        & 
                   (1.-GAMMA_INC(1.5*XBM3+XALPHA_SNOW+1,ZR2(JL)/ZZBETA(JL)))))
ELSE
    ZSUM_SUBL = 2*XPI*ZUSI(JL)*ZSVT(JL,1)/ZAI(JL)*(XANU*ZZBETA(JL)*XALPHA_SNOW  +      &
                 XBNU/ZGAM*(2/ZMU(JL))**0.5*(                                          & 
                 ZZBETA(JL)**(1.5*XBM1+1)*ZAM1(JL)**0.5*                        &
                       GAMMA_INC_LOW(1.5*XBM1+XALPHA_SNOW+1,ZR1(JL)/ZZBETA(JL)) +            &
                 ZZBETA(JL)**(1.5*XBM2+1)*ZAM2(JL)**0.5*                        &
                       (GAMMA_INC_LOW(1.5*XBM2+XALPHA_SNOW+1,ZR2(JL)/ZZBETA(JL))-            &
                       GAMMA_INC_LOW(1.5*XBM2+XALPHA_SNOW+1,ZR1(JL)/ZZBETA(JL)))+            & 
                 ZZBETA(JL)**(1.5*XBM3+1)*ZAM3(JL)**0.5*                        & 
                   (ZGAM_BM3-GAMMA_INC_LOW(1.5*XBM3+XALPHA_SNOW+1,ZR2(JL)/ZZBETA(JL)))))
END IF
!   Restriction of ZSUM_SUBL   
    ZTEMP=ZSUM_SUBL
    ZSUM_SUBL =  MIN( ZRVS(JL),ZSUM_SUBL)*(0.5+SIGN(0.5,ZSUM_SUBL)) &
                   - MIN(ZSVS(JL,2),ABS(ZSUM_SUBL))*(0.5-SIGN(0.5,ZSUM_SUBL))
    ZSUM_SUBL=MIN(0.,ZSUM_SUBL) ! Sink of snow               
    IF(ZSUM_SUBL>0) THEN  
       write(*,*) 'Warning Subl',JL,'Subl',ZSUM_SUBL,'TEMP',ZTEMP
       write(*,*) 'Warning Subl ZSVT',ZSVT(JL,1),ZSVT(JL,2)
       write(*,*) 'Warning vap',ZRVS(JL),'ZSVs',ZSVS(JL,2)
    END IF
!   Change in concentration rate  Sn = Sb*N/qb (Dery and Yau,2000)
    ZNUM      = ZSUM_SUBL*ZSVT(JL,1)/ZSVT(JL,2)
!   Change in mobility index : value mob*M3 is reduced according to reduction of
!   M3 due to sublimation so that mob is constant due to sublimation
!    ZMOB      = ZSUM_SUBL*ZSVT(JL,3)/ZSVT(JL,2)
!   Update tendencies for snow particles, water vapour and potential temperature 
    ZSVS(JL,2) = ZSVS(JL,2) + ZSUM_SUBL             ! Particle mixing ratio
    ZSVS(JL,1) = ZSVS(JL,1) + ZNUM                  ! Particles number
!    ZSVS(JL,3) = ZSVS(JL,3) + ZMOB
    ZRVS(JL) = ZRVS(JL) - ZSUM_SUBL                 ! Water vapour
    ZTHS(JL) = ZTHS(JL) + ZSUM_SUBL*ZLSFACT(JL)     ! Potential temperature
    ZSNWSUBL(JL) = ZSUM_SUBL*ZRHODREF(JL)           ! Sublimation rate kg/m3/s
END DO


END IF
END IF

END SUBROUTINE SNOW_SUBL
!
!-------------------------------------------------------------------------------
!
FUNCTION GET_INDEX(PBETA,PDELTAR) RESULT(KMAX)
!
!!    PURPOSE
!!    -------
!     Calculate the upper index in numerical integration of Carrier's formulation 
!     Index equals to 5* mean radius
!
!
USE MODD_BLOWSNOW,     ONLY : XALPHA_SNOW


!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)                                  :: PDELTAR      ! (-)
REAL, DIMENSION(:), INTENT(IN)                                  :: PBETA ! (kg/m3)

!
INTEGER, DIMENSION(SIZE(PBETA,1)) :: KMAX ! (-)
!

KMAX(:)=int(PBETA(:)*XALPHA_SNOW*5/PDELTAR)


END FUNCTION GET_INDEX
!
!-------------------------------------------------------------------------------
!
FUNCTION RLIM(PMU,PRHODREF,PBEST_LIM) RESULT(PRLIM)
!
!!    PURPOSE
!!    -------
!     Calculate the radius of a sperical particle for a given Best Number 
!
!
USE MODD_CSTS_BLOWSNOW,     ONLY : XRHOLI,XG
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)                                  :: PRHODREF ! (kg/m3)
REAL, DIMENSION(:), INTENT(IN)                                  :: PMU      ! (m2/s)
REAL,               INTENT(IN)                                  :: PBEST_LIM! (-)

!
REAL, DIMENSION(SIZE(PMU,1)) :: PRLIM ! (m)
!
PRLIM(:)=(3./32.*PRHODREF(:)/(XRHOLI*XG)*PMU(:)**2.*PBEST_LIM)**0.333333333

END FUNCTION RLIM

FUNCTION AVR(PARE,PBRE,PRHODREF,PMU) RESULT(PAVR)
!
!!    PURPOSE
!!    -------
!     Calculate the parameter av_r in KC02 formulation (Eq. 3.1)
!
!
USE MODD_CSTS_BLOWSNOW,     ONLY : XRHOLI,XG


!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL,               INTENT(IN)                                  :: PARE      ! (-)
REAL,               INTENT(IN)                                  :: PBRE      ! (-)
REAL, DIMENSION(:), INTENT(IN)                                  :: PRHODREF ! (kg/m3)
REAL, DIMENSION(:), INTENT(IN)                                  :: PMU      ! (m2/s)

!
REAL, DIMENSION(SIZE(PMU,1)) :: PAVR ! (-)
!


PAVR(:)=2.**(3.*PBRE-1.)*PARE*PMU(:)**(1.-2.*PBRE)*(4./3.*XRHOLI/PRHODREF(:)*XG)**PBRE

END FUNCTION AVR
!
!-------------------------------------------------------------------------------
!
FUNCTION TURB_FLUC(PR,PMU,PCARRIER,PRHODREF,PZZ,PVMOD) RESULT(PSIG)
!
!!    PURPOSE
!!    -------
!     Calculate the relative turbulent velocity fluctuations for a given radius. 
!     Used to compute the ventilation velocity.
!     Formulation based on Dover (1993)
!
USE MODD_CSTS
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                                  :: PR ! (m)
REAL, INTENT(IN)                                  :: PMU ! (m2/s)
REAL, INTENT(IN)                                  :: PCARRIER ! (m/s)
REAL, INTENT(IN)                                  :: PRHODREF ! (kg/m3)
REAL, INTENT(IN)                                  :: PZZ ! (m)
REAL, INTENT(IN)                                  :: PVMOD ! (m/s)
!
REAL              :: PSIG ! (m/s)
!
!
!*      0.2    declaration of local variables
!
REAL                  :: ZFCRI1,ZFCRI2,ZFCRI
REAL                  :: ZS0,ZSIGU,ZSIGV,ZSIGW,ZUSTAR
!
!
!*      1    Calculate critical frequency
!
ZFCRI1 = 9*PRHODREF*PMU/(4*XPI*PR**2*XRHOLI)
ZFCRI2 = 0.363*PRHODREF*PCARRIER/(XPI*PR*XRHOLI)
ZFCRI = MAX(ZFCRI1,ZFCRI2) 
!
!*      2    Calculate variances of the horizontal and vertical velocity components
!
ZS0   = ZFCRI*PZZ/PVMOD
ZSIGU = 4.77 *ZUSTAR**2/ (1+33*ZS0)**0.66666
ZSIGV = 2.76 *ZUSTAR**2/ (1+9.5*ZS0)**0.66666
ZSIGW = 1.31 *ZUSTAR**2/ (1+3.12*ZS0)**0.66666

PSIG = (ZSIGU+ZSIGV+ZSIGW)**0.5
END FUNCTION TURB_FLUC
!
FUNCTION NUSSELT(PR,PMU,PVEL_VENT) RESULT(PNU)
!
!!    PURPOSE
!!    -------
!     Calculate the Nusselt number for a given particle radius 
!     Formulation based on Lee (1975)
!
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                                  :: PR ! (m)
REAL, INTENT(IN)                                  :: PMU ! (m2/s)
REAL, INTENT(IN)                                  :: PVEL_VENT ! (m/s)
!
REAL              :: PNU ! (m/s)
!
!
!*      0.2    declaration of local variables
!
REAL                  :: ZRE
!
!
!*      1    Calculate Reynolds number 
!
ZRE = 2*PR*PVEL_VENT/(PMU)
!
!*      2    Calculate Nusselt number
!
IF(ZRE<10) THEN
        PNU = 1.79+0.606*ZRE**0.5
ELSE
        PNU = 1.88+0.580*ZRE**0.5
END IF

END FUNCTION NUSSELT
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE SUBL_BLOWSNOW



