snowcro.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ##########################################################################
      SUBROUTINE snowcro(HSNOWRES, TPTIME, OGLACIER, HIMPLICIT_WIND,     &
               ppew_a_coef, ppew_b_coef,                                 &
               ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef,       &
               psnowswe,psnowrho,psnowheat,psnowalb,                     &
               psnowgran1,psnowgran2,psnowhist,psnowage,                 &
               ptstep,pps,psr,prr,ppsn3l,                                &
               pta,ptg,psw_rad,pqa,pvmod,plw_rad, prhoa,                 &
               puref,pexns,pexna,pdircoszw,                              &
               pzref,pz0,pz0eff,pz0h,palb,                               &
               psoilcond,pd_g,                                           &
               psnowliq,psnowtemp,psnowdz,                               &
               pthrufal,pgrndflux,pevapcor,prnsnow,phsnow,pgfluxsnow,    &
               phpsnow,ples3l,plel3l,pevap,psndrift,pri,                 &
               pemisnow,pcdsnow,pustar,pchsnow,psnowhmass,pqs,           &
               ppermsnowfrac,pzenith,pxlat,pxlon,                        &
               osnowdrift,osnowdrift_sublim,osnow_abs_zenith,            &
               hsnowmetamo,hsnowrad) 
!     ##########################################################################
!
!!****  *SNOWCRO*
!!
!!    PURPOSE
!!    -------
!
!     Detailed snowpack scheme Crocus, computationnally based on the
!     3-Layer snow scheme option (Boone and Etchevers 1999)
!     For shallow snow cover, Default method of Douville et al. (1995)
!     used with this option: Model "turns on" when snow sufficiently
!     deep/above some preset critical snow depth.
!
!
!
!
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!
!!
!!    REFERENCE
!!    ---------
!!
!!    ISBA: Belair (1995)
!!    ISBA: Noilhan and Planton (1989)
!!    ISBA: Noilhan and Mahfouf (1996)
!!    ISBA-ES: Boone and Etchevers (2001)
!!    Crocus : Brun et al., 1989 (J. Glaciol.)
!!    Crocus : Brun et al., 1992 (J. Glaciol.)
!!    Crocus : Vionnet et al., in prep (Geosci. Mod. Devel. Discuss.)
!!
!!
!!    AUTHOR
!!    ------
!!      A. Boone           * Meteo-France *
!!      V. Vionnet         * Meteo-France *
!!      E. Brun            * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    7/99
!!      Modified by A.Boone 05/02 (code, not physics)
!!      Modified by A.Boone 11/04 i) maximum density limit imposed (although
!!                                rarely if ever reached), ii) check to
!!                                see if upermost layer completely sublimates
!!                                during a timestep (as snowpack becomes vanishly
!!                                thin), iii) impose maximum grain size limit
!!                                in radiation transmission computation.
!!
!!      Modified by B. Decharme  (03/2009): Consistency with Arpege permanent
!!                                          snow/ice treatment (LGLACIER for alb)
!!      Modified by A. Boone     (04/2010): Implicit coupling and replace Qsat and DQsat
!!                                          by Qsati and DQsati, respectively.
!!      Modified by E. Brun, V. Vionnet, S. Morin (05/2011): 
!!                                          Addition of Crocus processes and 
!!                                          parametrizations to 
!!                                          the SNOW-3L code. This includes the dynamic handling
!!                                          of snow layers and the inclusion of snow metamorphism
!!                                          rules similar to the original Crocus implementation.
!!      Modified by B. Decharme  (09/2012): New wind implicitation
!!
!!      Modified by M. Lafaysse (07/2012) : 
!!                                          * Albedo and roughness parametrizations
!!                                            for surface ice over glaciers
!!                                                     MODIF 2012-10-03 : don't modify roughness if implicit coupling
!!                                                                 (test PPEW_A_COEF == 0. )
!!                                          * SNOWCROALB is now called by SNOWCRORAD to remove duplicated code
!!                                          * Parameters for albedo are moved to modd_snow_par
!!                                          * PSNOWAGE is stored as an age
!!                                            (days since snowfall) and not as a date
!!                                            to allow spinup simulations
!!                                          * New rules for optimal discretization of very thick snowpacks
!!                                          * Optional outputs for debugging
!!
!!       Modified by E. Brun and M. Lafaysse (07/2012) :
!!                                          * Implement sublimation in SNOWDRIFT
!!                                          * Flag in namelist to activate SNOWDRIFT and SNOWDRIFT_SUBLIM
!!       Modified by E. Brun and M. Lafaysse (08/2012) :
!!                                          * XUEPSI replaced by 0 in the if statement of case 1.3.3.2 (SNOWCROMETAMO)
!!                                          * If SNOWDRIFT is activated the wind do not modify grain types during snowfall
!!                                            (redundant with snowdrift)
!!       Modified by E. Brun (24/09/2012) :
!!                                          * Correction coupling coefficient for specific humidity in SNOWCROEBUD
!!                                          * PSFCFRZ(:)  = 1.0 for systematic solid/vapor latent fluxes in SNOWCROEBUD
!!       Modified by C. Carmagnola (3/2013): 
!!                                          * Dendricity and size replaced by the optical diameter
!!                                          * Test of different evolution laws for the optical diameter
!!
!!       Modified by B. Decharme  (08/2013): Qsat as argument (needed for coupling with atm)
!!                                           add PSNDRIFT
!!
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_type_date_surf, ONLY: date_time
!
USE modd_csts, ONLY : xtt, xrholw, xlmtt,xlstt,xlvtt, xcl, xci, xpi, xrholi
USE modd_snow_par, ONLY : xz0icez0snow, xrhothreshold_ice
USE modd_snow_metamo
USE modd_const_tartes, ONLY: npnimp, xpsnowg0, xpsnowy0, xpsnoww0, xpsnowb0
!
USE mode_snow3l
USE mode_tartes, ONLY : snowcro_tartes
!
USE mode_crodebug
!
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
! this module is not used anymore
! USE MODI_GREGODSTRATI
!
IMPLICIT NONE
!
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!                                      PTSTEP    = time step of the integration
TYPE(date_time), INTENT(IN)         :: tptime      ! current date and time
!
 CHARACTER(LEN=*), INTENT(IN)       :: hsnowres
!                                      HSNOWRES  = ISBA-SNOW3L turbulant exchange option
!                                      'DEF' = Default: Louis (ISBA: Noilhan and Mahfouf 1996)
!                                      'RIL' = Limit Richarson number under very stable
!                                              conditions (currently testing)
LOGICAL, INTENT(IN)                 :: oglacier   ! True = Over permanent snow and ice, 
!                                                     initialise WGI=WSAT,
!                                                     Hsnow>=10m and allow 0.8<SNOALB<0.85
                                                  ! False = No specific treatment
!
 CHARACTER(LEN=*), INTENT(IN)       :: himplicit_wind   ! wind implicitation option
!                                                     ! 'OLD' = direct
!                                                     ! 'NEW' = Taylor serie, order 1
!
REAL, DIMENSION(:), INTENT(IN)      :: pps, pta, psw_rad, pqa, pvmod, plw_rad, psr, prr 
!                                      PSW_RAD = incoming solar radiation (W/m2)
!                                      PLW_RAD = atmospheric infrared radiation (W/m2)
!                                      PRR     = rain rate [kg/(m2 s)]
!                                      PSR     = snow rate (SWE) [kg/(m2 s)]
!                                      PTA     = atmospheric temperature at level za (K)
!                                      PVMOD   = modulus of the wind parallel to the orography (m/s)
!                                      PPS     = surface pressure
!                                      PQA     = atmospheric specific humidity
!                                                at level za
!
REAL, DIMENSION(:), INTENT(IN)      :: ptg, psoilcond, pd_g, ppsn3l
!                                      PTG       = Surface soil temperature (effective
!                                                  temperature the of layer lying below snow)
!                                      PSOILCOND = soil thermal conductivity [W/(m K)]
!                                      PD_G      = Assumed first soil layer thickness (m)
!                                                  Used to calculate ground/snow heat flux
!                                      PPSN3L    = snow fraction
!
REAL, DIMENSION(:), INTENT(IN)      :: pzref, puref, pexns, pexna, pdircoszw, prhoa, pz0, pz0eff, &
                                       palb, pz0h, ppermsnowfrac 
!                                      PZ0EFF    = roughness length for momentum
!                                      PZ0       = grid box average roughness length
!                                      PZ0H      = grid box average roughness length for heat
!                                      PZREF     = reference height of the first
!                                                  atmospheric level
!                                      PUREF     = reference height of the wind
!                                      PRHOA     = air density
!                                      PEXNS     = Exner function at surface
!                                      PEXNA     = Exner function at lowest atmos level
!                                      PDIRCOSZW = Cosinus of the angle between the
!                                                  normal to the surface and the vertical
!                                      PALB      = soil/vegetation albedo
!                                      PPERMSNOWFRAC  = fraction of permanet snow/ice
!
REAL, DIMENSION(:), INTENT(IN)      :: ppew_a_coef, ppew_b_coef,                   &
                                        ppet_a_coef, ppeq_a_coef, ppet_b_coef,      &
                                        ppeq_b_coef  
!                                      PPEW_A_COEF = wind coefficient (m2s/kg)
!                                      PPEW_B_COEF = wind coefficient (m/s)
!                                      PPET_A_COEF = A-air temperature coefficient
!                                      PPET_B_COEF = B-air temperature coefficient
!                                      PPEQ_A_COEF = A-air specific humidity coefficient
!                                      PPEQ_B_COEF = B-air specific humidity coefficient
!
REAL, DIMENSION(:), INTENT(INOUT)   :: psnowalb
!                                      PSNOWALB = Prognostic surface snow albedo
!                                                 (does not include anything but
!                                                 the actual snow cover)
!
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowheat, psnowrho, psnowswe
!                                      PSNOWHEAT = Snow layer(s) heat content (J/m2)
!                                      PSNOWRHO  = Snow layer(s) averaged density (kg/m3)
!                                      PSNOWSWE  = Snow layer(s) liquid Water Equivalent (SWE:kg m-2)
!
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowgran1, psnowgran2, psnowhist
!                                      PSNOWGRAN1 = Snow layers grain feature 1
!                                      PSNOWGRAN2 = Snow layer grain feature 2
!                                      PSNOWHIST  = Snow layer grain historical
!                                                   parameter (only for non
!                                                   dendritic snow)
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowage  ! Snow grain age
!
REAL, DIMENSION(:,:), INTENT(OUT)    :: psnowliq, psnowtemp, psnowdz
!                                      PSNOWLIQ  = Snow layer(s) liquid water content (m)
!                                      PSNOWTEMP = Snow layer(s) temperature (m)
!                                      PSNOWDZ   = Snow layer(s) thickness (m)
!
REAL, DIMENSION(:), INTENT(OUT)      :: pthrufal, pgrndflux, pevapcor
!                                      PTHRUFAL  = rate that liquid water leaves snow pack:
!                                                  paritioned into soil infiltration/runoff
!                                                  by ISBA [kg/(m2 s)]
!                                      PGRNDFLUX = soil/snow interface heat flux (W/m2)
!                                      PEVAPCOR  = evaporation/sublimation correction term:
!                                                  extract any evaporation exceeding the
!                                                  actual snow cover (as snow vanishes)
!                                                  and apply it as a surface soil water
!                                                  sink. [kg/(m2 s)]
!
REAL, DIMENSION(:), INTENT(OUT)      :: prnsnow, phsnow, pgfluxsnow, ples3l, plel3l, &
                                        phpsnow, pcdsnow, pustar, pevap, psndrift
!                                      PLES3L      = evaporation heat flux from snow (W/m2)
!                                      PLEL3L      = sublimation (W/m2)
!                                      PHPSNOW     = heat release from rainfall (W/m2)
!                                      PRNSNOW     = net radiative flux from snow (W/m2)
!                                      PHSNOW      = sensible heat flux from snow (W/m2)
!                                      PGFLUXSNOW  = net heat flux from snow (W/m2)
!                                      PCDSNOW     = drag coefficient for momentum over snow
!                                      PUSTAR      = friction velocity over snow (m/s)
!                                      PEVAP       = total evaporative flux (kg/m2/s)
!                                      PSNDRIFT    = blowing snow sublimation (kg/m2/s)
!
REAL, DIMENSION(:), INTENT(OUT)      :: pchsnow, pemisnow, psnowhmass
!                                      PEMISNOW    = snow surface emissivity
!                                      PCHSNOW     = drag coefficient for heat over snow
!                                      PSNOWHMASS  = heat content change due to mass
!                                                    changes in snowpack (J/m2): for budget
!                                                    calculations only.
!
REAL, DIMENSION(:), INTENT(OUT)      :: pri, pqs
!                                      PRI = Ridcharson number
!                                      PQS = surface humidity
!
REAL, DIMENSION(:), INTENT(IN)        :: pzenith ! solar zenith angle
REAL, DIMENSION(:), INTENT(IN)        :: pxlat,pxlon ! LAT/LON after packing
!
LOGICAL, INTENT(IN)                   :: osnowdrift, osnowdrift_sublim ! activate snowdrift, sublimation during drift
LOGICAL, INTENT(IN)                   :: osnow_abs_zenith ! activate parametrization of solar absorption for polar regions
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo, hsnowrad
                                         !-----------------------
                                         ! Metamorphism scheme
                                         ! HSNOWMETAMO=B92 Brun et al 1992
                                         ! HSNOWMETAMO=C13 Carmagnola et al 2014
                                         ! HSNOWMETAMO=T07 Taillandier et al 2007
                                         ! HSNOWMETAMO=F06 Flanner et al 2006
                                         !-----------------------
                                         ! Radiative transfer scheme
                                         ! HSNOWMETAMO=B92 Brun et al 1992
                                         ! HSNOWMETAMO=TAR TARTES (Libois et al 2013)
                                         ! HSNOWMETAMO=TA1 TARTES with constant impurities
                                         ! HSNOWMETAMO=TA2 TARTES with constant impurities as function of ageing
                                         !-----------------------                                         
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),NPNIMP) :: zsnowimp_density !impurities density (kg/m^3) (npoints,nlayer,ntypes_impurities)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2),NPNIMP) :: zsnowimp_content !impurities content (g/g) (npoints,nlayer,ntypes_impurities)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowtemp, zscap, zsnowdzn, zscond, zradsink 
!                                      ZSNOWTEMP  = Snow layer(s) averaged temperature (K)
!                                      ZSCAP      = Snow layer(s) heat capacity [J/(K m3)]
!                                      ZSNOWDZN   = Updated snow layer thicknesses (m)
!                                      ZSCOND     = Snow layer(s) thermal conducivity [W/(m K)]
!                                      ZRADSINK   = Snow solar Radiation source terms (W/m2)
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zwholdmax 
!
!For now these values are constant
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowg0 ! asymmetry parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit) (npoints,nlayer)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowy0 ! Value of y of snow grains at nr=1.3 (no unit
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnoww0 ! Value of W of snow grains at nr=1.3 (no unit)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowb0 ! absorption enhancement parameter of snow grains at nr=1.3 and at non absorbing wavelengths (no unit)
!
!spectral albedo (3 bands for now) :: ready to output if necessary
REAL, DIMENSION(SIZE(PSNOWRHO,1),3) :: zspectralalbedo
!  
REAL, DIMENSION(SIZE(PTA))          :: zsnowbis
!                                     ZSNOWBIS      = Total snow depth after snowfall  
!
REAL, DIMENSION(SIZE(PTA))          :: zsnow, zsfcfrz, ztsterm1, ztsterm2, zct, zra, zsnowflux, zsnowtempo1 
!                                      ZSNOW      = Total snow depth (m)
!                                      ZCT        = inverse of the product of snow heat capacity
!                                                   and layer thickness [(m2 K)/J]
!                                      ZRA        = Surface aerodynamic resistance
!                                      ZTSTERM1,ZTSTERM2 = Surface energy budget coefficients
!                                      ZSNOWFLUX  = heat flux between 1st and 2nd snow layers:
!                                                   used during surface melting (W/m2)
!                                      ZSNOWTEMPO1= value of uppermost snow temperature
!                                                   before time integration (K)
!
REAL, DIMENSION(SIZE(PTA))          :: zrsra, zdqsat, zqsat, zradxs, zliqheatxs, zlwupsnow 
!                                      ZRSRA    = air density over aerodynamic resistance
!                                      ZDQSAT   = derrivative of saturation specific humidity
!                                      ZQSAT    = saturation specific humidity
!                                      ZRADXS   = shortwave radiation absorbed by soil surface
!                                                 (for thin snow sover) (W m-2)
!                                      ZLIQHEATXS = excess snowpack heating for vanishingly thin
!                                                 snow cover: add energy to snow/ground heat
!                                                 flux (W m-2)
!                                      ZLWUPSNOW = upwelling longwave raaditive flux (W m-2)
!
REAL, DIMENSION(SIZE(PTA))          :: zustar2_ic, zta_ic, zqa_ic, &
                                       zpet_a_coef_t, zpeq_a_coef_t, zpet_b_coef_t, zpeq_b_coef_t  
!                                      ZUSTAR2_IC    = implicit lowest atmospheric level friction (m2/s2)
!                                      ZTA_IC        = implicit lowest atmospheric level air temperature
!                                      ZQA_IC        = implicit lowest atmospheric level specific humidity
!                                      ZPET_A_COEF_T = transformed A-air temperature coefficient
!                                      ZPET_B_COEF_T = transformed B-air temperature coefficient
!                                      ZPEQ_A_COEF_T = transformed A-air specific humidity coefficient
!                                      ZPEQ_B_COEF_T = transformed B-air specific humidity coefficient
!
REAL, DIMENSION(SIZE(PTA))          :: zsnowrhof, zsnowdzf, zsnowgran1f, zsnowgran2f, zsnowhistf 
REAL, DIMENSION(SIZE(PTA))          :: zsnowagef

! New roughness lengths in case of glaciers without snow.
REAL, DIMENSION(SIZE(PTA))          :: zz0_snowice, zz0h_snowice, zz0eff_snowice
!
!To control and print eneregy balance
REAL , DIMENSION(SIZE(PTA))         :: zsummass_ini,zsumheat_ini,zsummass_fin,zsumheat_fin
!
REAL, DIMENSION(SIZE(PTA))          :: zmassbalance, zenergybalance, zevapcor2
!
INTEGER, DIMENSION(SIZE(PTA))       :: inlvls_use ! varying number of effective layers 
!
LOGICAL, DIMENSION(SIZE(PTA))       :: gsnowfall,gmodif_maillage
!                                      GSNOWFALL  = FLAG if snowfall exceed PSNOW/10, used for 
!                                                   grid updating. 
!
REAL :: ztstepdays ! time step in days
!
LOGICAL :: gcond_grain, gcond_yen
!
LOGICAL :: gcroinfoprint ! print daily informations
LOGICAL :: gcrodebugprint, gcrodebugdetailsprint, gcrodebugprintatm ! print diagnostics for debugging
LOGICAL :: gcrodebugprintbalance
!
INTEGER :: jj,jst  ! looping indexes
INTEGER :: iprint  ! gridpoint number to be printed
INTEGER :: idebug
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO',0,zhook_handle)
!
!***************************************PRINT IN**********************************************
! Look if we have to print snow profiles for debugging
gcroinfoprint = lcrodailyinfo .AND. (tptime%TIME ==0.0)
!***************************************PRINT OUT*********************************************
!***************************************DEBUG IN**********************************************
gcrodebugprintbalance = ( tptime%TDATE%YEAR*10000 + tptime%TDATE%MONTH*100 + tptime%TDATE%DAY &
                        >= ntimecrodebug ) .AND. &
                        ( tptime%TIME/3600. >= nhourcrodebug ) .AND. &
                        ( tptime%TDATE%YEAR*10000 + tptime%TDATE%MONTH*100 + tptime%TDATE%DAY &
                        < nendcrodebug )
!
IF (lcrodebug) THEN
  gcrodebugprint        = gcrodebugprintbalance
  gcrodebugdetailsprint = lcrodebugdetails .AND. gcrodebugprint 
  gcrodebugprintatm     = lcrodebugatm .AND. gcrodebugprint
ELSE
  gcrodebugprint        = .false.
  gcrodebugdetailsprint = .false.
  gcrodebugprintatm     = .false.
END IF

!
! Look if we have to compute and print energy balance control
gcrodebugprintbalance = lcontrolbalance .AND. gcrodebugprintbalance 
!
IF ( lcrodebug .OR. gcroinfoprint .OR. gcrodebugprintbalance ) THEN
  !
  IF ( xlatcrodebug >= -90 .AND. xloncrodebug >= -180. ) THEN
    CALL getpoint_crodebug(pxlat,pxlon,idebug)
  ELSE
    idebug = npointcrodebug
  END IF
  !
  ! For parallel runs : we just want to do this for the thread where there is this point.
  IF ( xlatcrodebug >= -90 .AND. xloncrodebug >= -180. ) THEN
    IF ( abs( pxlat(idebug)-xlatcrodebug ) + abs(pxlon(idebug) - xloncrodebug) > 0.1 ) THEN
      gcrodebugprint        = .false.
      gcrodebugdetailsprint = .false.
      gcrodebugprintatm     = .false.
    END IF
  END IF
  !
END IF
!***************************************DEBUG OUT*********************************************
!
IF ( hsnowrad=="TAR" .OR. hsnowrad=="TA1" .OR. hsnowrad=="TA2" ) THEN
  !For now fix constant values
  zsnowg0 = xpsnowg0
  zsnowy0 = xpsnowy0
  zsnoww0 = xpsnoww0  
  zsnowb0 = xpsnowb0
  !
  zsnowimp_density = 1500.
  ! ZSNOWIMP_CONTENT=25.0E-9
  !
END IF
!
zustar2_ic = 0.0
zta_ic     = 0.0
zqa_ic     = 0.0
!
gcond_grain = .true.
gcond_yen   = .true.!FALSE. !(if TRUE : use of the Yen (1981) thermal conductivity paramztrization ; 
!                    otherwise, use the default ISBA-ES thermal conductivity parametrization)
!
pgrndflux  = 0.
psnowhmass = 0.
phsnow     = 0.
prnsnow    = 0.
ples3l     = 0.
plel3l     = 0.
phpsnow    = 0.
pevapcor   = 0.
pthrufal   = 0.
!
! pour imprimer des diagnostics sur un des points
iprint = 1 
!
! - - ---------------------------------------------------
!
!       0.     Initialization
!               --------------
! NOTE that snow layer thickness is used throughout this code: SWE
! is only used to diagnose the thickness at the beginning of this routine
! and it is updated at the end of this routine.
!
! Initialization of the actual number of snow layers, total snow depth 
!  and layer thicknesses
!
zsnowtemp(:,:) = 0.
!
gsnowfall(:) = .false.
inlvls_use(:) = 0
DO jst = 1,SIZE(psnowswe(:,:),2)
  DO jj = 1,SIZE(zsnow)
    IF ( psnowswe(jj,jst)>0. ) THEN 
      psnowdz(jj,jst) = psnowswe(jj,jst) / psnowrho(jj,jst)
      inlvls_use(jj) = jst
    ELSE
      psnowdz(jj,jst) = 0.
    ENDIF
  ENDDO  !  end loop snow layers
ENDDO    ! end loop grid points
! Incrementation of snow layers age
ztstepdays = ptstep/86400. ! time step in days
WHERE ( psnowswe>0 ) psnowage = psnowage + ztstepdays
!
!***************************************PRINT IN**********************************************
!
!Compute total SWE and heat for energy control 
IF ( gcrodebugprintbalance ) THEN
  DO jj = 1,SIZE(zsnow)
    zsummass_ini(jj) = sum(psnowswe(jj,1:inlvls_use(jj)))
    zsumheat_ini(jj) = sum(psnowheat(jj,1:inlvls_use(jj)))
  ENDDO    ! end loop grid points
ENDIF
!
!  Print of some simulation characteristics
IF(gcroinfoprint) THEN
  CALL snowcroprintdate()
  WRITE(*,fmt="(A12,I3,A12,I4)") 'nlayer:',inlvls_use(idebug), ' nbpoints:', SIZE(zsnow) 
!   WRITE(*,*) 'PZ0H: ', PZ0H(IDEBUG)
  WRITE(*,*) 'Snow fraction =',ppsn3l(idebug)                   
ENDIF
!
!***************************************PRINT OUT*********************************************
!***************************************DEBUG IN**********************************************
IF (gcrodebugprint) THEN
  CALL snowcroprintdate()
  CALL snowcroprintprofile("crocus initialization",inlvls_use(idebug),lprintgran,      &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo)
END IF
!
IF (gcrodebugprintatm) THEN
  CALL snowcroprintatm("forcing data :",pta(idebug),pqa(idebug),pvmod(idebug),   &
                       prr(idebug),psr(idebug),psw_rad(idebug),plw_rad(idebug),  &
                       ptg(idebug),psoilcond(idebug),pd_g(idebug),ppsn3l(idebug) )
END IF
!***************************************DEBUG OUT********************************************
!
!*       1.     Snow total depth
!               ----------------
!
zsnow(:) = 0.
DO jj = 1,SIZE(zsnow)
  zsnow(jj) = sum(psnowdz(jj,1:inlvls_use(jj)))
ENDDO
!
zsnowbis(:) = zsnow(:)
!
!*       2.     Snowfall
!               --------
! Calculate uppermost density and thickness changes due to snowfall,
! and add heat content of falling snow
!
 CALL snownlfall_upgrid(tptime, oglacier,                                       &
                        ptstep,psr,pta,pvmod,zsnowbis,psnowrho,psnowdz,         &
                        psnowheat,psnowhmass,psnowalb,ppermsnowfrac,            &
                        psnowgran1,psnowgran2,gsnowfall,zsnowdzn,               &
                        zsnowrhof,zsnowdzf,zsnowgran1f,zsnowgran2f, zsnowhistf, &
                        zsnowagef,gmodif_maillage,inlvls_use,osnowdrift,pz0eff,puref,&
                        hsnowmetamo) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWFALL_UPGRID",inlvls_use(idebug),lprintgran,      &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************   
!
zsnow(:) = zsnowbis(:)
!
!*       3.     Update grid/discretization
!               --------------------------
! Reset grid to conform to model specifications:
!
DO jj=1,SIZE(zsnow)
  !
  IF ( gmodif_maillage(jj) ) THEN    
    CALL snownlgridfresh_1d(jj,zsnow(jj),psnowdz(jj,:),zsnowdzn(jj,:),psnowrho(jj,:),    &
                            psnowheat(jj,:),psnowgran1(jj,:),psnowgran2(jj,:),           &
                            psnowhist(jj,:),psnowage(jj,:),gsnowfall(jj),zsnowrhof(jj),  &
                            zsnowdzf(jj),psnowhmass(jj),zsnowgran1f(jj),zsnowgran2f(jj), &
                            zsnowhistf(jj),zsnowagef(jj),inlvls_use(jj),hsnowmetamo      ) 
  ENDIF
  !
ENDDO 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWNLGRIDFRESH_1D",inlvls_use(idebug),lprintgran,   &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************
!
!*       4.     Liquid water content and snow temperature
!               -----------------------------------------
!
! First diagnose snow temperatures and liquid
! water portion of the snow from snow heat content:
! update some snow layers parameters after new discretization 
!
DO jj = 1,SIZE(zsnow)
  !
  ! active layers
  DO jst=1,inlvls_use(jj)
    psnowswe(jj,jst) = psnowdz(jj,jst) * psnowrho(jj,jst)
!
    zscap(jj,jst) = psnowrho(jj,jst) * xci
!
    zsnowtemp(jj,jst) = xtt + &
                        ( ( psnowheat(jj,jst)/psnowdz(jj,jst) + xlmtt*psnowrho(jj,jst) )/zscap(jj,jst) ) 
!
    psnowliq(jj,jst) = max( 0.0, zsnowtemp(jj,jst)-xtt ) * zscap(jj,jst) * &
                        psnowdz(jj,jst) / (xlmtt*xrholw) 
!
    zsnowtemp(jj,jst) = min( xtt, zsnowtemp(jj,jst) )
  ENDDO  !  end loop active snow layers
  !
  ! unactive layers  
  DO jst = inlvls_use(jj)+1,SIZE(psnowswe,2)
    psnowswe(jj,jst) = 0.0
    psnowrho(jj,jst) = 999.
    psnowdz(jj,jst) = 0.
    psnowgran1(jj,jst) = 0.
    psnowgran2(jj,jst) = 0.
    psnowhist(jj,jst) = 0.
    psnowage(jj,jst) = 0.
    psnowheat(jj,jst) = 0.
    zsnowtemp(jj,jst) = xtt
    psnowliq(jj,jst) = 0.
  ENDDO  !  end loop unactive snow layers
  !
ENDDO    ! end loop grid points
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after liquid water/temperature diagnostic",                &
                           inlvls_use(idebug),lprintgran,                              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************
!
!        4.BIS   Snow metamorphism
!                ----------------- 
!  
 CALL snowcrometamo(psnowdz,psnowgran1,psnowgran2,psnowhist,zsnowtemp, &
                    psnowliq,ptstep,psnowswe,inlvls_use,psnowage,hsnowmetamo       ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROMETAMO", inlvls_use(idebug),lprintgran,       &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************
!
!*       5.     Snow Compaction
!               ---------------
! Calculate snow density: compaction/aging: density increases
!
 CALL snowcrocompactn(ptstep,psnowrho,psnowdz,zsnowtemp,zsnow,                      &
                      psnowgran1,psnowgran2,psnowhist,psnowliq,inlvls_use,pdircoszw,&
                      hsnowmetamo)
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROCOMPACTN", inlvls_use(idebug),lprintgran,     &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************
!       
!*       5.1    Snow Compaction and Metamorphism due to snow drift
!               ---------------
psndrift(:) = 0.0
IF (osnowdrift) THEN
  CALL snowdrift(ptstep, pvmod, psnowrho,psnowdz, zsnow,                      &
                 psnowgran1,psnowgran2,psnowhist,inlvls_use,pta,pqa,pps,prhoa,&
                 pz0eff,puref,osnowdrift_sublim,hsnowmetamo,psndrift)
ENDIF
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWDRIFT", inlvls_use(idebug),lprintgran,           &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),   &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),&
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),&
                           hsnowmetamo )
ENDIF
!***************************************DEBUG OUT**********************************************
!
! Update snow heat content (J/m2) using dry density instead of total density:
!
DO jj = 1,SIZE(zsnow)
  DO jst = 1,inlvls_use(jj)
    zscap(jj,jst) = ( psnowrho(jj,jst) - &
                      psnowliq(jj,jst) * xrholw / max( psnowdz(jj,jst),xsnowdzmin) ) * xci
    psnowheat(jj,jst) = psnowdz(jj,jst) * &
                        ( zscap(jj,jst)*(zsnowtemp(jj,jst)-xtt) - xlmtt*psnowrho(jj,jst) ) + &
                        xlmtt * xrholw * psnowliq(jj,jst) 
  ENDDO  !  end loop snow layers
ENDDO    ! end loop grid points
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after  update snow heat content", inlvls_use(idebug),lprintgran,&
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo      )
ENDIF
!***************************************DEBUG OUT********************************************
!
!*       6.     Solar radiation transmission
!               -----------------------------
!
! Heat source (-sink) term due to shortwave
! radiation transmission within the snowpack:
!
SELECT CASE (hsnowrad)
  CASE ("TA1")
    zsnowimp_content(:,:,1) = 0.0
  CASE ("TA2")
    zsnowimp_content(:,:,1) = 100.0e-9
  CASE ("TAR")
    zsnowimp_content(:,:,1) = 2. * psnowage(:,:) * 1e-9
  CASE default
END SELECT
!
SELECT CASE (hsnowrad)
  CASE ("B92")
    CALL snowcrorad(tptime,oglacier,                      &
                    psw_rad,psnowalb,psnowdz,psnowrho,    &
                    palb,zradsink,zradxs,                 &
                    psnowgran1, psnowgran2, psnowage,pps, &
                    pzenith, ppermsnowfrac,inlvls_use,    &
                    osnow_abs_zenith,hsnowmetamo) 
  !
  CASE ("TAR","TA1","TA2")
    CALL snowcro_tartes(psnowgran1,psnowgran2,psnowrho,psnowdz,zsnowg0,zsnowy0,zsnoww0, &
                        zsnowb0,zsnowimp_density,zsnowimp_content,palb,psw_rad,pzenith, &
                        inlvls_use,psnowalb,zradsink,zradxs,gcrodebugdetailsprint,hsnowmetamo)
  !
  CASE default
    CALL abor1_sfx("UNKNOWN CSNOWRAD OPTION")
  !
END SELECT
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCRORAD", inlvls_use(idebug),lprintgran,               &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
!*       7.     Heat transfer and surface energy budget
!               ---------------------------------------
! Snow thermal conductivity:
!
 CALL snowcrothrm(psnowrho,zscond,zsnowtemp,pps,psnowliq,gcond_grain,gcond_yen)
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROTHRM", inlvls_use(idebug),lprintgran,              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
! Precipitation heating term:
! Rainfall renders it's heat to the snow when it enters
! the snowpack:
!
phpsnow(:) = prr(:) * xcl * ( max( xtt,pta(:) ) - xtt )    ! (W/m2)
!
! Surface Energy Budget calculations using ISBA linearized form
! and standard ISBA turbulent transfer formulation
!
IF ( all(ppew_a_coef==0.) ) THEN
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Modif Matthieu Lafaysse for glaciers
  ! For surface ice, modify roughness lengths
  ! Only if not implicit coupling
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  WHERE( psnowrho(:,1)>xrhothreshold_ice )
    zz0_snowice    = pz0    * xz0icez0snow
    zz0h_snowice   = pz0h   * xz0icez0snow
    zz0eff_snowice = pz0eff * xz0icez0snow
  ELSEWHERE
    zz0_snowice    = pz0
    zz0h_snowice   = pz0h
    zz0eff_snowice = pz0eff
  ENDWHERE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ELSE
  zz0_snowice    = pz0
  zz0h_snowice   = pz0h
  zz0eff_snowice = pz0eff
END IF

 CALL snowcroebud(hsnowres, himplicit_wind,                                    &
                  ppew_a_coef, ppew_b_coef,                                    &
                  ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef,          &
                  xsnowdzmin,                                                  &
                  pzref,zsnowtemp(:,1),psnowrho(:,1),psnowliq(:,1),zscap(:,1), &
                  zscond(:,1),zscond(:,2),                                     &
                  puref,pexns,pexna,pdircoszw,pvmod,                           &
                  plw_rad,psw_rad,pta,pqa,pps,ptstep,                          &
                  psnowdz(:,1),psnowdz(:,2),psnowalb,zz0_snowice,              &
                  zz0eff_snowice,zz0h_snowice,                                 &
                  zsfcfrz,zradsink(:,1),phpsnow,                               &
                  zct,pemisnow,prhoa,ztsterm1,ztsterm2,zra,pcdsnow,pchsnow,    &
                  zqsat, zdqsat, zrsra, zustar2_ic, pri,                       &
                  zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t      ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROEBUD", inlvls_use(idebug),lprintgran,              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
! Heat transfer: simple diffusion along the thermal gradient
!
zsnowtempo1(:) = zsnowtemp(:,1) ! save surface snow temperature before update
!
 CALL snowcrosolvt(ptstep,xsnowdzmin,psnowdz,zscond,zscap,ptg,                &
                   psoilcond,pd_g,zradsink,zct,ztsterm1,ztsterm2,             &
                   zpet_a_coef_t,zpeq_a_coef_t,zpet_b_coef_t,zpeq_b_coef_t,   &
                   zta_ic,zqa_ic,pgrndflux, zsnowtemp ,zsnowflux,             &
                   inlvls_use                                                 ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROSOLVT", inlvls_use(idebug),lprintgran,             &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
!*       8.     Surface fluxes
!               --------------
!
 CALL snowcroflux(zsnowtemp(:,1),psnowdz(:,1),pexns,pexna,            &
                  zustar2_ic,                                         &
                  ptstep,psnowalb,psw_rad,pemisnow,zlwupsnow,plw_rad, &
                  zta_ic,zsfcfrz,zqa_ic,phpsnow,                      &
                  zsnowtempo1,zsnowflux,zct,zradsink(:,1),            &
                  zqsat,zdqsat,zrsra,                                 &
                  prnsnow,phsnow,pgfluxsnow,ples3l,plel3l,pevap,      &
                  pustar                                              ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROFLUX", inlvls_use(idebug),lprintgran,              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
!*       9.     Snow melt
!               ---------
!
! First Test to see if snow pack vanishes during this time step:
!
 CALL snowcrogone(ptstep,plel3l,ples3l,psnowrho,                  &
                  psnowheat,zradsink,pevapcor,pthrufal,pgrndflux, &
                  pgfluxsnow,psnowdz,psnowliq,zsnowtemp,zradxs,   &
                  prr,inlvls_use                                  ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROGONE", inlvls_use(idebug),lprintgran,              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
! Add radiation not absorbed by snow to soil/vegetation interface flux
! (for thin snowpacks):
!
pgrndflux(:) = pgrndflux(:) + zradxs(:)
!
! Second Test to see if one or several snow layers vanishe during this time
! step. In such a case, the concerned snow layers are agregated to neighbours

 CALL snowcrolayer_gone(ptstep,zscap,zsnowtemp,psnowdz,          &
                        psnowrho,psnowliq,psnowgran1,psnowgran2, &
                        psnowhist,psnowage,ples3l, inlvls_use    ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROLAYER_GONE", inlvls_use(idebug),lprintgran,        &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
! For partial melt: transform excess heat content into snow liquid:
!
 CALL snowcromelt(zscap,zsnowtemp,psnowdz,psnowrho,psnowliq,inlvls_use)
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROMELT", inlvls_use(idebug),lprintgran,              &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************             
!
!*      10.     Snow water flow and refreezing
!               ------------------------------
! Liquid water vertical transfer and possible snowpack runoff
! And refreezing/freezing of meltwater/rainfall (ripening of the snow)
!
 CALL snowcrorefrz(ptstep,prr,psnowrho,zsnowtemp,psnowdz,psnowliq,pthrufal, &
                   zscap,plel3l,inlvls_use                                  ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROREFRZ", inlvls_use(idebug),lprintgran,             &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
!*      11.     Snow Evaporation/Sublimation mass updates:
!               ------------------------------------------
!
 CALL snowcroevapn(ples3l,ptstep,zsnowtemp(:,1),psnowrho(:,1), &
                   psnowdz(:,1),pevapcor,psnowhmass            ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROEVAPN", inlvls_use(idebug),lprintgran,             &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
! If all snow in uppermost layer evaporates/sublimates, re-distribute
! grid (below could be evoked for vanishingly thin snowpacks):
!
 CALL snowcroevapgone(psnowheat,psnowdz,psnowrho,zsnowtemp,psnowliq,psnowgran1, &
                      psnowgran2,psnowhist,psnowage,inlvls_use,hsnowmetamo      ) 
!
!***************************************DEBUG IN**********************************************
IF (gcrodebugdetailsprint) THEN
  CALL snowcroprintprofile("after SNOWCROEVAPGONE", inlvls_use(idebug),lprintgran,          &
                           psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),        &
                           psnowliq(idebug,:),psnowheat(idebug,:),psnowgran1(idebug,:),     &
                           psnowgran2(idebug,:),psnowhist(idebug,:),psnowage(idebug,:),     &
                           hsnowmetamo)
ENDIF
!***************************************DEBUG OUT********************************************
!
!*      12.     Update surface albedo:
!               ----------------------
! Snow clear sky albedo:
!
IF ( hsnowrad=='B92' ) THEN
  CALL snowcroalb(tptime,oglacier,                                             &
                  psnowalb,zspectralalbedo,psnowdz(:,1),psnowrho(:,1:2),       &
                  ppermsnowfrac,psnowgran1(:,1),psnowgran2(:,1),               &
                  psnowage(:,1),psnowgran1(:,2),psnowgran2(:,2),psnowage(:,2), &
                  pps, pzenith, inlvls_use, hsnowmetamo) 
  !the albedo is not updated in the case of TARTES scheme                 
ENDIF
!
!*      13.     Update snow heat content:
!               -------------------------
! Update the heat content (variable stored each time step)
! using current snow temperature and liquid water content:
!
! First, make check to make sure heat content not too large
! (this can result due to signifigant heating of thin snowpacks):
! add any excess heat to ground flux:
!
DO jj = 1,SIZE(zsnow)
! active layers
  DO jst = 1,inlvls_use(jj)
    zwholdmax(jj,jst) = snowcrohold( psnowrho(jj,jst),psnowliq(jj,jst),psnowdz(jj,jst) )
    zliqheatxs(jj)     = max( 0.0, (psnowliq(jj,jst) - zwholdmax(jj,jst)) * xrholw ) * xlmtt/ptstep 
    psnowliq(jj,jst) = psnowliq(jj,jst) - zliqheatxs(jj)*ptstep/(xrholw*xlmtt)
    psnowliq(jj,jst) = max( 0.0, psnowliq(jj,jst) )
    pgrndflux(jj)     = pgrndflux(jj) + zliqheatxs(jj)
    psnowtemp(jj,jst) = zsnowtemp(jj,jst)
!   Heat content using total density
    zscap(jj,jst) = psnowrho(jj,jst) * xci
    psnowheat(jj,jst) = psnowdz(jj,jst) * &
                        ( zscap(jj,jst)*(psnowtemp(jj,jst)-xtt) - xlmtt*psnowrho(jj,jst) ) + &
                        xlmtt * xrholw * psnowliq(jj,jst) 
!
    psnowswe(jj,jst)  = psnowdz(jj,jst) * psnowrho(jj,jst)
  ENDDO  !  end loop active snow layers
!  
! unactive layers
  DO jst = inlvls_use(jj)+1,SIZE(psnowswe,2)
    psnowswe(jj,jst)  = 0.
    psnowheat(jj,jst) = 0.
    psnowrho(jj,jst)  = 999.
    psnowtemp(jj,jst) = 0.
    psnowdz(jj,jst)   = 0.
  ENDDO  !  end loop unactive snow layers
!  
ENDDO    ! end loop grid points
!
! print some final diagnostics
! ! ! IF (INLVLS_USE(I)>0) THEN
! ! !  WRITE(*,FMT="(I4,2I4,F4.0,A7,F5.2,A10,F7.1,A11,F6.2,A13,F6.2)")  &
! ! !      TPTIME%TDATE%YEAR,TPTIME%TDATE%MONTH,TPTIME%TDATE%DAY,TPTIME%TIME/3600.,&
! ! !           'HTN= ',SUM(PSNOWDZ(I,1:INLVLS_USE(I))), 'FLUX Sol=', PGRNDFLUX(I),&
! ! !   'Tsurf_sol=',PTG(I)-273.16, 'Tbase_neige=',PSNOWTEMP(I,INLVLS_USE(I))-273.16 
! ! ! ENDIF
!
!***************************************DEBUG IN*********************************************
IF (gcrodebugprint) THEN
  CALL snowcroprintdate()
  CALL snowcroprintprofile("CROCUS : end of time step",inlvls_use(idebug),lprintgran, &
        psnowdz(idebug,:),psnowrho(idebug,:),psnowtemp(idebug,:),psnowliq(idebug,:),  &
        psnowheat(idebug,:),psnowgran1(idebug,:),psnowgran2(idebug,:),                &
        psnowhist(idebug,:),psnowage(idebug,:), hsnowmetamo)
END IF
!***************************************DEBUG OUT********************************************
!***************************************PRINT IN*********************************************
! check suspect low temperature
DO jj = 1,SIZE(zsnow)
!IF(INLVLS_USE(JJ)>0) WRITE(*,*) 'SOL:',JJ,INLVLS_USE(JJ),PGRNDFLUX(JJ),PTG(JJ),&
! PSNOWTEMP(jj,INLVLS_USE(JJ)),PSNOWTEMP(jj,1),PZENITH(JJ)           
  DO jst = 1,inlvls_use(jj)        
    IF ( psnowtemp(jj,jst) < 100. ) THEN
      WRITE(*,*) 'pb tempe snow :',psnowtemp(jj,jst)
      WRITE(*,fmt='("DATE:",2(I2.2,"/"),I4.4,F3.0)')          &
        tptime%TDATE%DAY,tptime%TDATE%MONTH,tptime%TDATE%YEAR,tptime%TIME/3600.
      WRITE(*,*) 'point',jj,"/",SIZE(zsnow)
      WRITE(*,*) 'layer',jst
      WRITE(*,*) 'pressure',pps(jj)
      WRITE(*,*) 'slope',acos(pdircoszw(jj))*(180./xpi),"deg"
      WRITE(*,*) 'XLAT=',pxlat(jj),'XLON=',pxlon(jj)
      WRITE(*,*) 'solar radiation=',psw_rad(jj)
      WRITE(*,*) 'INLVLS_USE(JJ):',inlvls_use(jj)
      WRITE(*,*) psnowdz(jj,1:inlvls_use(jj))
      WRITE(*,*) psnowrho(jj,1:inlvls_use(jj))
      WRITE(*,*) psnowtemp(jj,1:inlvls_use(jj))
      CALL abor1_sfx('SNOWCRO: erreur tempe snow')
    ENDIF
  ENDDO
ENDDO
!***************************************PRINT OUT*********************************************
!***************************************DEBUG IN*********************************************
!Control and print energy balance
IF (gcrodebugprintbalance) THEN
  !
  zsummass_fin(idebug) = sum( psnowswe(idebug,1:inlvls_use(idebug)) )
  zsumheat_fin(idebug) = sum( psnowheat(idebug,1:inlvls_use(idebug)) )
  !
  CALL get_balance(zsummass_ini(idebug),zsumheat_ini(idebug),zsummass_fin(idebug), &
                   zsumheat_fin(idebug),psr(idebug),prr(idebug),pthrufal(idebug),  &
                   pevap(idebug),pevapcor(idebug),pgrndflux(idebug),phsnow(idebug),&
                   prnsnow(idebug),plel3l(idebug),ples3l(idebug),phpsnow(idebug),  &
                   psnowhmass(idebug),psnowdz(idebug,1),ptstep,                    &
                   zmassbalance(idebug),zenergybalance(idebug),zevapcor2(idebug)   )
  !
  CALL snowcroprintbalance(zsummass_ini(idebug),zsumheat_ini(idebug),zsummass_fin(idebug), &
                           zsumheat_fin(idebug),psr(idebug),prr(idebug),pthrufal(idebug),  &
                           pevap(idebug),pevapcor(idebug),pgrndflux(idebug),phsnow(idebug),&
                           prnsnow(idebug),plel3l(idebug),ples3l(idebug),phpsnow(idebug),  &
                           psnowhmass(idebug),psnowdz(idebug,1),ptstep,                    &
                           zmassbalance(idebug),zenergybalance(idebug),zevapcor2(idebug))
  !
ENDIF
!
IF (lpstopbalance) THEN
  !
  ! bilan pour tous points pour test eventuel sur depassement seuil des residus
  DO jj=1, SIZE(zsnow)
    !
    zsummass_fin(jj) = sum( psnowswe(jj,1:inlvls_use(jj)) )
    zsumheat_fin(jj) = sum( psnowheat(jj,1:inlvls_use(jj)) )
    !
    CALL get_balance(zsummass_ini(jj),zsumheat_ini(jj),zsummass_fin(jj), &
                     zsumheat_fin(jj),psr(jj),prr(jj),pthrufal(jj),      &
                     pevap(jj),pevapcor(jj),pgrndflux(jj),phsnow(jj),    &
                     prnsnow(jj),plel3l(jj),ples3l(jj),phpsnow(jj),      &
                     psnowhmass(jj),psnowdz(jj,1),ptstep,                &
                     zmassbalance(jj),zenergybalance(jj),zevapcor2(jj)   )
    !
  ENDDO    ! end loop grid points
  !
  CALL snowcrostopbalance(zmassbalance(:),zenergybalance(:))
  !
END IF
!***************************************DEBUG OUT********************************************
!
pqs(:) = zqsat(:)
!
IF (lhook) CALL dr_hook('SNOWCRO',1,zhook_handle)
!
 CONTAINS
!
!####################################################################
!####################################################################
!####################################################################
        SUBROUTINE snowcrocompactn(PTSTEP,PSNOWRHO,PSNOWDZ,                         &
                                   psnowtemp,psnow,psnowgran1,psnowgran2,psnowhist, &
                                   psnowliq,inlvls_use,pdircoszw,hsnowmetamo        ) 
!
!!    PURPOSE
!!    -------
!     Snow compaction due to overburden and settling.
!     Mass is unchanged: layer thickness is reduced
!     in proportion to density increases. Method
!     directly inherited from Crocus v2.4 and
!     coarsely described in Brun et al., J. Glac 1989 and 1992
!
!     de/e = -sigma/eta * dt
!
!     where e is layer thickness, sigma is the vertical stress, dt is the
!     time step and eta is the snow viscosity
!     * sigma is currently calculated taking into account only the overburden
!     (a term representing "metamorphism stress" in fresh snow may be added
!      in the future)
!     * eta is computed as a function of snowtype, density and temperature
!
!     The local slope is taken into account, through the variable PDIRCOSZW
!     which is directly the cosine of the local slope
!
!
!     HISTORY:
!     Basic structure from ISBA-ES model (Boone and Etchevers, 2001)
!     Implementation of Crocus laws : E. Brun, S. Morin, J.-M. Willemet July 2010.
!     Implementation of slope effect on settling : V. Vionnet, S. Morin May 2011
!
!
USE modd_csts,     ONLY : xtt, xg
USE modd_snow_par, ONLY : xrhosmax_es
USE modd_snow_metamo
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep       ! Time step UNIT : s
REAL, DIMENSION(:), INTENT(IN)      :: pdircoszw    ! cosine of local slope
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowtemp    ! Snow temperature UNIT : K
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz   ! Density UNIT : kg m-3, Layer thickness UNIT : m
!
REAL, DIMENSION(:), INTENT(OUT)     :: psnow        ! Snowheight UNIT : m
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowgran1, psnowgran2, psnowhist, &!Snowtype variables
                                        psnowliq     ! Snow liquid water content UNIT ??? 
INTEGER, DIMENSION(:), INTENT(IN)   :: inlvls_use   ! Number of snow layers used
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo ! metamorphism scheme
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho2,    &! work snow density UNIT : kg m-3
                                                      zviscosity,   &! Snow viscosity UNIT : N m-2 s (= Pa s)
                                                      zsmass        !, &  ! overburden mass for a given layer UNIT : kg m-2 
!                                                      ZWSNOWDZ       ! mass of each snow layer UNIT : kg m-2
!
INTEGER   :: jj,jst   ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
! 1. Cumulative snow mass (kg/m2):
! --------------------------------
!
IF (lhook) CALL dr_hook('SNOWCROCOMPACTN',0,zhook_handle)
!
DO jj = 1,SIZE(psnow) 
  zsmass(jj,1) = 0.0
  DO jst = 2,inlvls_use(jj)
    zsmass(jj,jst) = zsmass(jj,jst-1) + psnowdz(jj,jst-1) * psnowrho(jj,jst-1)
  ENDDO
  zsmass(jj,1) = 0.5 * psnowdz(jj,1) * psnowrho(jj,1)  ! overburden of half the mass of the uppermost layer applied to itself 
ENDDO

!
! 2. Compaction/Settling
! ----------------------
!
DO jj = 1,SIZE(psnow)
  !
  DO jst = 1,inlvls_use(jj)
    !
    ! Snow viscosity basic equation (depend on temperature and density only):
    zviscosity(jj,jst) = xvvisc1 * &
                         exp( xvvisc3*psnowrho(jj,jst) + xvvisc4*abs(xtt-psnowtemp(jj,jst)) ) * &
                         psnowrho(jj,jst) / xvro11 
    !
    ! Equations below apply changes to the basic viscosity value, based on snow microstructure properties
    IF ( psnowliq(jj,jst)>0.0 ) THEN
      zviscosity(jj,jst) = zviscosity(jj,jst) / &
                           ( xvvisc5 + xvvisc6*psnowliq(jj,jst)/psnowdz(jj,jst) ) 
    ENDIF
    !
    IF( psnowliq(jj,jst)/psnowdz(jj,jst)<=0.01 .AND. psnowhist(jj,jst)>=nvhis2 ) THEN   
      zviscosity(jj,jst) = zviscosity(jj,jst) * xvvisc7
    ENDIF
    !
    IF ( psnowhist(jj,jst)==nvhis1 ) THEN
      !
      IF ( hsnowmetamo=="B92" ) THEN
        !
        IF ( psnowgran1(jj,jst)>=0. .AND. psnowgran1(jj,jst)<xvgran6 ) THEN   
          zviscosity(jj,jst) = zviscosity(jj,jst) * &
                               min( 4., exp( min( xvdiam1, &
                                                  psnowgran2(jj,jst)             -xvdiam4 ) / xvdiam6 ) )
        ENDIF
        !
      ELSEIF ( psnowgran1(jj,jst)>=xvdiam6*(4.-psnowgran2(jj,jst)) .AND. psnowgran2(jj,jst)<xvgran6/xvgran1 ) THEN 
        zviscosity(jj,jst) = zviscosity(jj,jst) * &
                             min( 4., exp( min( xvdiam1, &
                                               (xvdiam6*(4.-psnowgran2(jj,jst)))-xvdiam4 ) / xvdiam6 ) )
      ENDIF
      !
    ENDIF
    !
    ! Calculate new snow snow density: compaction from weight/over-burden
    zsnowrho2(jj,jst) = psnowrho(jj,jst) + psnowrho(jj,jst) * ptstep * &
                                           ( xg*pdircoszw(jj)*zsmass(jj,jst)/zviscosity(jj,jst) )
    !    
    ! Calculate new grid thickness in response to density change
    psnowdz(jj,jst) = psnowdz(jj,jst) * ( psnowrho(jj,jst)/zsnowrho2(jj,jst) )
    !
    !  Update density (kg m-3):
    psnowrho(jj,jst) = zsnowrho2(jj,jst)
    !
  ENDDO    ! end loop snow layers
  !
ENDDO      ! end loop grid points
!
!
! 3. Update total snow depth:
! -----------------------------------------------
!
DO jj = 1,SIZE(psnowdz,1)
  psnow(jj) = sum( psnowdz(jj,1:inlvls_use(jj)) )
ENDDO
!
IF (lhook) CALL dr_hook('SNOWCROCOMPACTN',1,zhook_handle)

!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snowcrocompactn


!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrometamo(PSNOWDZ,PSNOWGRAN1, PSNOWGRAN2,         &
                         psnowhist, psnowtemp, psnowliq, ptstep, &
                         psnowswe,inlvls_use, psnowage, hsnowmetamo)  
!                               

!**** *METAMO* - METAMORPHOSE DES GRAINS
!              - SNOW METAMORPHISM
!     OBJET.
!     ------
!     METAMORPHOSE DU MANTEAU NEIGEUX.
!     EVOLUTION DU TYPE DE GRAINS 
!     MISE A JOUR DES VARIABLES HISTORIQUES.
!     METAMORPHISM OF THE SNOW GRAINS,
!     HISTORICAL VARIABLES

!**   INTERFACE.
!     ----------
!     FORMALISME ADOPTE POUR LA REPRESENTATION DES GRAINS :
!     FORMALISM FOR THE REPRESENTATION OF GRAINS
!     -----------------------------------------------------


!                    1       - -1                 NEIGE FRAICHE
!                   / \      |                    -------------
!                  /   \     |  DENDRICITE        DECRITE EN TERME
!                 /     \    |  DENDRICITY        DE DENDRICITE ET
!                /       \   |                    SPHERICITE
!               2---------3  -  0                 DESCRIBED WITH
!                                                 SPHERICITY AND
!               |---------|                       DENDRICITY
!               0         1
!               SPHERICITE
!               SPHERICITY

!               4---------5  -
!               |         |  |
!               |         |  | DIAMETRE  (OU TAILLE)
!               |         |  | DIAMETER  (OR SIZE  )
!               |         |  |                 
!               |         |  |                   NEIGE NON DENDRITIQUE
!               6---------7  -                   ---------------------

!                                                SPHERICITE ET TAILLE
!                                                SPHERICITY AND SIZE 

!              LES VARIABLES DU MODELE : 
!              -------------------------
!              CAS DENDRITIQUE             CAS NON DENDRITIQUE
!  
!            SGRAN1(JST) : DENDRICITE      SGRAN1(JST) : SPHERICITE
!            SGRAN2(JST) : SPHERICITE      SGRAN2(JST) : TAILLE (EN METRE)
!                                                        SIZE

! 
!    CAS DENDRITIQUE/ DENDRITIC CASE
!    -------------------------------
!    SGRAN1(JST) VARIE DE -XVGRAN1 (-99 PAR DEFAUT) (ETOILE) A 0
!  (DENDRICITY)  >D OU LA DIVISION PAR -XVGRAN1 POUR OBTENIR DES VALEURS
!                 ENTRE 1 ET 0
!                 VARIES FROM -XVGRAN1 (DEFAULT -99) (FRESH SNOW) TO 0
!                 DIVISION BY -XVGRAN1 TO OBTAIN VALUES BETWEEN 0 AND 1

!    SGRAN2(JST) VARIE DE 0 (CAS COMPLETEMENT ANGULEUX) A XVGRAN1
!  (SPHERICITY)  (99 PAR DEFAUT)
!                >D OU LA DIVISION PAR XVGRAN1 POUR OBTENIR DES VALEURS
!                 ENTRE 0 ET 1
!                 VARIES FROM 0 (SPHERICITY=0) TO XVGRAN1


!    CAS NON DENDRITIQUE / NON DENDRITIC CASE
!    ---------------------------------------

!    SGRAN1(JST) VARIE DE 0 (CAS COMPLETEMENT ANGULEUX) A XVGRAN1
!  (SPHERICITY)  (99 PAR DEFAUT) (CAS SPHERIQUE)
!                >D OU LA DIVISION PAR XVGRAN1 POUR OBTENIR DES VALEURS
!                 ENTRE 0 ET 1
!                 VARIES FROM 0 TO 99

!    SGRAN2(JST) EST SUPERIEUR A XVDIAM1-SPHERICITE (3.E-4 M) ET NE FAIT QUE CROITRE
!     (SIZE)     IS GREATER THAN XVDIAM1-SPHERICITE (3.E-4 M) ALWAYS INCREASE


!    EXEMPLES : POINTS CARACTERISTIQUES DE LA FIGURE
!    --------

!                 SGRAN1     SGRAN2    DENDRICITE  SPHERICITE  TAILLE
!                                      DENDRICITY  SPHERICITY  SIZE
!      --------------------------------------------------------------
!                                                               (M)
!        1        -XVGRAN1    VNSPH3        1           0.5
!        2           0         0           0            0
!        3           0       XVGRAN1        0            1
!        4           0       XVDIAM1                     0        4.E-4 
!        5         XVGRAN1    XVDIAM1-XVSPHE1              1        3.E-4
!        6           0         --                       0        --
!        7         XVGRAN1      --                       1        --

!     PAR DEFAUT : XVGRAN1 =99   VNSPH3=50 XVSPHE1=1. XVDIAM1=4.E-4


!     METHODE.
!     --------
!     EVOLUTION DES TYPES DE GRAINS : SELON LES LOIS DECRITES 
!     DANS BRUN ET AL (1992)
!     PLUSIEURS CAS SONT A DISTINGUER
!      1.2 NEIGE HUMIDE
!      1.3 METAMORPHOSE NEIGE SECHE
!        1.3.1 FAIBLE GRADIENT
!        1.3.2 GRADIENT MOYEN
!        1.3.3 FORT GRADIENT
!     DANS CHAQUE CAS ON SEPARE NEIGE DENDRITIQUE ET NON DENDRITIQUE
!     LE PASSAGE DENDRITIQUE => NON DENDRITIQUE SE FAIT LORSQUE 
!     SGRAN1 DEVIENT > 0
         
!     TASSEMENT : LOIS DE VISCOSITE ADAPTEE SELON LE TYPE DE GRAINS

!     VARIABLES HISTORIQUES (CAS NON DENDRITIQUE SEULEMENT)

!     MSHIST DEFAUT
!        0           CAS NORMAL
!     NVHIS1   1     GRAINS ANGULEUX
!     NVHIS2   2     GRAINS AYANT ETE EN PRESENCE D EAU LIQUIDE 
!                    MAIS N'AYANT PAS EU DE CARATERE ANGULEUX 
!     NVHIS3   3     GRAINS AYANT ETE EN PRESENCE D EAU LIQUIDE
!                    AYANT EU AUPARAVANT UN CARACTERE ANGULEUX 

!     GRAIN METAMORPHISM ACCORDING TO BRUN ET AL (1992)
!     THE DIFFERENT CASES ARE :
!     1.2 WET SNOW
!     1.3 DRY SNOW
!       1.3.1. LOW      TEMPERATURE GRADIENT
!       1.3.2. MODERATE TEMPERATURE GRADIENT
!       1.3.3. HIGH     TEMPERATURE GRADIENTi
!     THE CASE OF DENTRITIC OR NON DENDRITIC SNOW IS TREATED SEPARATELY
!     THE LIMIT DENTRITIC ==> NON DENDRITIC IS REACHED WHEN SGRAN1>0

!     SNOW SETTLING : VISCOSITY DEPENDS ON THE GRAIN TYPES

!     HISTORICAL VARIABLES (NON DENDRITIC CASE)
!     MSHIST DEFAUT
!        0           CAS NORMAL
!     NVHIS1   1     FACETED CRISTAL
!     NVHIS2   2     LIQUID WATER AND NO FACETED CRISTALS BEFORE
!     NVHIS3   3     LIQUID WATER AND FACETED CRISTALS BEFORE

!     EXTERNES.
!     ---------

!     REFERENCES.
!     -----------

!     AUTEURS.
!     --------
!        ERIC BRUN ET AL. - JOURNAL OF GLACIOLOGY 1989/1992.

!     MODIFICATIONS.
!     --------------
!        08/95: YANNICK DANIELOU - CODAGE A LA NORME DOCTOR.
!        09/96: ERIC MARTIN      - CORRECTION COMMENTAIRES
!        03/06: JM Willemet      - F90 and SI units
!        08/06: JM Willemet      - new formulation for TEL (Mwat/(Mice+Mwat) instead of Mwat/Mice.
!                                  Threshold on the diameter increasing of the wet grains.
!        01/07 : JM Willemet     - CORRECTION DES COUCHES SATUREES SUBISSANT DU TASSEMENT
!                                  CORRECTION ON THE SATURATED LAYERS WHICH ARE SETTLED
!        12/12: CM Carmagnola    - Dendricity and size replaced by the optical diameter
!                                - Test of different evolution laws for the optical diameter
!        08/13: M Lafaysse       - Simplification of historical parameter computation (logicals GNONDENDRITIC, GFACETED, GSPHE_LW)
   !
USE modd_snow_metamo
USE modd_csts, ONLY : xtt, xpi, xrholw, xrholi
USE modd_surf_par, ONLY : xundef
!
USE mode_snow3l
!
IMPLICIT NONE 
!
!     0.1 declarations of arguments  
!      
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowdz, psnowtemp, psnowliq, psnowswe 
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowgran1, psnowgran2, psnowhist              
!   
REAL, INTENT(IN)                    :: ptstep    
! 
INTEGER, DIMENSION(:), INTENT(IN)   :: inlvls_use      
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowage
!
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo ! metamorphism scheme
!
!     0.2 declaration of local variables      
!     
REAL :: zgradt, ztelm, zvdent, zdent, zsphe, zvap, zdangl, &
        zsize, zssa, zssa0, zssa_t, zssa_t_dt, za, zb, zc, &
        za2, zb2, zc2, zoptd, zoptr, zoptr0, zdrdt
REAL :: zvdent1, zvdent2, zvsphe, zcoef_sph
REAL :: zdenom1, zdenom2, zfact1, zfact2                          
INTEGER :: inlvls
INTEGER :: jst,jj                                !Loop controls 
INTEGER :: idrho, idgrad, idtemp           !Indices for values from Flanner 2006
LOGICAL :: gnondendritic ,gfaceted, gsphe_lw
LOGICAL :: gcond_b92, gcond_c13, gcond_sph
!
REAL(KIND=JPRB) :: zhook_handle
!
!     INITIALISATION
!     --------------
!
IF (lhook) CALL dr_hook('SNOWCROMETAMO',0,zhook_handle)
!
inlvls = SIZE(psnowgran1(:,:),2)  ! total snow layers
!
!*    1. METAMORPHOSES DANS LES STRATES. / METAMORPHISM
!        -----------------------------------------------
DO jj = 1,SIZE(psnowrho,1)
  !
  DO jst = 1,inlvls_use(jj)
    !
    ! 1.1 INITIALISATION: GRADIENT DE TEMPERATURE / TEMPERATURE GRADIENT
    IF ( jst==inlvls_use(jj) ) THEN
      zgradt = abs(psnowtemp(jj,jst)   - psnowtemp(jj,jst-1))*2. / (psnowdz(jj,jst-1) + psnowdz(jj,jst)) 
    ELSEIF ( jst==1 ) THEN
      zgradt = abs(psnowtemp(jj,jst+1) - psnowtemp(jj,jst)  )*2. / (psnowdz(jj,jst) + psnowdz(jj,jst+1))
    ELSE
      zgradt = abs(psnowtemp(jj,jst+1) - psnowtemp(jj,jst-1))*2. / &
                                            (psnowdz(jj,jst-1) + psnowdz(jj,jst)*2. + psnowdz(jj,jst+1)) 
    ENDIF
    !
    IF ( psnowliq(jj,jst)>xuepsi ) THEN
      ! 1.2 METAMORPHOSE HUMIDE. / WET SNOW METAMORPHISM
      !
      ! TENEUR EN EAU LIQUIDE / LIQUID WATER CONTENT
      ztelm  = xupourc * psnowliq(jj,jst) * xrholw / psnowswe(jj,jst)
      !
      ! VITESSES DE DIMINUTION DE LA DENDRICITE / RATE OF THE DENDRICITY DECREASE
      zvdent1 = max( xvdent2 * ztelm**nvdent1, xvdent1 * exp(xvvap1/xtt) )
      zvdent2 = zvdent1
      ! CONDITION POUR LE CAS NON DENDRITIQUE NON SPHERIQUE
      gcond_b92  = ( psnowgran1(jj,jst)<xvgran1-xuepsi )
      gcond_c13 = .true. ! CONDITION POUR LE CALCUL DE SNOWGRAN1
      ! X COEF
      zvsphe = xvsphe1
      ! FOR C13
      zcoef_sph = 2.      
      !
    ELSEIF ( zgradt<xvgrat1 ) THEN
      ! 1.3.1 METAMORPHOSE SECHE FAIBLE/ DRY LOW GRADIENT (0-5 DEG/M).
      !
      zvap = exp( xvvap1/psnowtemp(jj,jst) )
      !
      ! VITESSES DE DIMINUTION DE LA DENDRICITE / RATE OF THE DENDRICITY DECREASE
      zvdent1 = xvdent1 * zvap
      zvdent2 = xvsphe2 * zvap
      ! CONDITION POUR LE CAS NON DENDRITIQUE SPHERICITE NON LIMITEE
      gcond_b92  = ( psnowhist(jj,jst)/=nvhis1 .OR. psnowgran2(jj,jst)<xvdiam2 )
      gcond_c13 = ( hsnowmetamo=='C13' )  ! CONDITION POUR LE CALCUL DE SNOWGRAN1
      ! X COEF
      zvsphe = xvsphe1
      ! FOR C13
      zcoef_sph = 2.
      !
    ELSE
      ! 1.3.2 METAMORPHOSE SECHE GRADIENT MOYEN / DRY MODERATE (5-15).
      ! 1.3.3 METAMORPHOSE SECHE FORT / DRY HIGH GRADIENT
      !
      zvap = exp( xvvap1/psnowtemp(jj,jst) ) * (zgradt)**xvvap2
      !
      ! VITESSES DE DIMINUTION DE LA DENDRICITE / RATE OF THE DENDRICITY DECREASE
      zvdent1 = xvdent1 * zvap
      zvdent2 = - xvdent1 * zvap
      ! CONDITION POUR LE CAS NON DENDRITIQUE NON COMPLETEMENT ANGULEUX
      gcond_b92  = ( zgradt<xvgrat2 .OR. psnowgran1(jj,jst)>0. )
      gcond_c13 = ( hsnowmetamo=='C13' ) ! CONDITION POUR LE CALCUL DE SNOWGRAN1
      ! X COEF
      zvsphe = xundef 
      ! FOR C13
      zcoef_sph = 3.
      !
    ENDIF       
    !
    IF ( hsnowmetamo=="B92" ) THEN
      !
      !------------------------------------------------
      !    BRUN et al. 1992 (B92)
      !
      ! -> Wet snow and dry snow 
      ! -> Evolution of dendricity, sphericity and size 
      !------------------------------------------------
      !
      IF ( psnowgran1(jj,jst)<-xuepsi ) THEN
        ! 1.2.1 CAS DENDRITIQUE/DENDRITIC CASE. 
        !
        ! / CALCUL NOUVELLE DENDRICITE ET SPHERICITE.
        zdent = - psnowgran1(jj,jst)/xvgran1 - zvdent1 * ptstep
        zsphe =   psnowgran2(jj,jst)/xvgran1 + zvdent2 * ptstep
        CALL set_thresh(zgradt,psnowliq(jj,jst),zsphe)
        IF( zdent<=xuepsi ) THEN
          ! EVOLUTION DE SGRAN1 ET SGRAN2 ET TEST PASSAGE DENDRITIQUE > NON DENDRITIQUE.
          psnowgran1(jj,jst) =  zsphe * xvgran1
          psnowgran2(jj,jst) = xvdiam1 - xvdiam5 * min( zsphe, zvsphe )
        ELSE
          psnowgran1(jj,jst) = -zdent * xvgran1
          psnowgran2(jj,jst) =  zsphe * xvgran1
        ENDIF          
        !
      ELSEIF ( gcond_b92 ) THEN
        ! 1.2.2 CAS NON DENDRITIQUE ET
        !             NON COMPLETEMENT SPHERIQUE / NON DENDRITIC AND NOT COMPLETELY SPHERIC CASE
        ! OU          SPHERICITE NON LIMITEE
        ! OU          NON COMPLETEMENT ANGULEUX
        !
        ! . EVOLUTION DE LA SPHERICITE SEULEMENT / EVOLUTION OF SPHERICITY ONLY (NOT SIZE)
        zsphe = psnowgran1(jj,jst)/xvgran1 + zvdent2 * ptstep
        CALL set_thresh(zgradt,psnowliq(jj,jst),zsphe)
        psnowgran1(jj,jst) = zsphe * xvgran1
        !          
      ELSEIF ( psnowliq(jj,jst)>xuepsi ) THEN
        ! 1.2.3 CAS NON DENDRITIQUE ET SPHERIQUE/NON DENDRITIC AND SPHERIC EN METAMORPHOSE HUMIDE
        !
        ! EVOLUTION DE LA TAILLE SEULEMENT/EVOLUTION OF SIZE ONLY                  
        CALL get_gran(ptstep,ztelm,psnowgran2(jj,jst))        
        !
      ELSEIF ( zgradt<xvgrat1 ) THEN
        ! 1.2.4. CAS HISTORIQUE=2 OU 3 ET GROS GRAINS SPHERICITE LIMITEE / CASE HISTORY=2 OR 3 AND BIG GRAINS LIMITED SPHERICITY 
        !
        zsphe =  psnowgran1(jj,jst)/xvgran1 + &
                 zvdent2 * ptstep * exp( min( 0., xvdiam3-psnowgran2(jj,jst) ) / xvdiam6 ) 
        zsphe = min( zsphe, xvsphe3 )
        CALL set_thresh(zgradt,psnowliq(jj,jst),zsphe)          
        psnowgran1(jj,jst) = zsphe * xvgran1
        !
      ELSE
        ! 1.2.5. CAS NON DENDRITIQUE ET ANGULEUX/DENDRITIC AND SPERICITY=0.
        !
        zdangl = snow3l_marbouty(psnowrho(jj,jst),psnowtemp(jj,jst),zgradt) 
        psnowgran2(jj,jst) = psnowgran2(jj,jst) + zdangl * xvfi * ptstep
        !
      ENDIF
      !
    ELSE
      !
      !------------------------------------------------
      !    CARMAGNOLA et al. 2013 (C13)
      !
      ! -> Wet snow
      ! -> Evolution of optical diameter and sphericity
      !------------------------------------------------
      !
      ! SPHERICITY
      zsphe = psnowgran2(jj,jst) + zvdent2 * ptstep
      CALL set_thresh(zgradt,psnowliq(jj,jst),zsphe)      
      IF ( psnowliq(jj,jst)>xuepsi .OR. zgradt<xvgrat1 ) THEN
        gcond_sph = ( zsphe < 1.-xuepsi )
      ELSE 
        gcond_sph = ( zsphe > xuepsi )
      ENDIF
      !      
      IF ( gcond_c13 .AND. psnowgran1(jj,jst)<xvdiam6*(4.-zsphe)-xuepsi ) THEN
        ! 1.1.1 CAS DENDRITIQUE/DENDRITIC CASE.
        !
        IF ( gcond_sph ) THEN
          psnowgran1(jj,jst) = psnowgran1(jj,jst) + xvdiam6 * ptstep * &
                               ( zvdent2*(psnowgran1(jj,jst)/xvdiam6-1.)/(zsphe-3.) - &
                                 zvdent1*(zsphe-3.) )
        ELSE
           psnowgran1(jj,jst) =  psnowgran1(jj,jst) + xvdiam6 * ptstep * zvdent1 * zcoef_sph
        ENDIF          
        !
      ELSEIF ( gcond_c13 .AND. gcond_sph ) THEN
        ! 1.2.2 CAS NON DENDRITIQUE ET
        !             NON COMPLETEMENT SPHERIQUE / NON DENDRITIC AND NOT COMPLETELY SPHERIC CASE
        ! OU          NON COMPLETEMENT ANGULEUX
        !        
        psnowgran1(jj,jst) = psnowgran1(jj,jst) - xvdiam6 * ptstep * zvdent2 * 2.* zsphe
        !
      ELSEIF ( psnowliq(jj,jst)>xuepsi ) THEN 
        ! 1.2.3 CAS NON DENDRITIQUE ET SPHERIQUE/NON DENDRITIC AND SPHERIC EN METAMORPHOSE HUMIDE      
        !
        ! NON DENDRITIC AND SPHERIC: EVOLUTION OF SIZE ONLY 
        CALL get_gran(ptstep,ztelm,psnowgran1(jj,jst)) 
        !
      ELSEIF ( gcond_c13 .AND. zgradt>=xvgrat2 ) THEN
        !
        zdangl = snow3l_marbouty(psnowrho(jj,jst),psnowtemp(jj,jst),zgradt)   
        psnowgran1(jj,jst) = psnowgran1(jj,jst) + 0.5 * zdangl * xvfi * ptstep
        !
      ENDIF          
      !
      psnowgran2(jj,jst) = zsphe
      !
      !---------------------------------
      !    TAILLANDIER et al. 2007 (T07)
      !
      ! -> Dry snow
      ! -> Evolution of optical diameter
      !---------------------------------
      !
      IF ( psnowliq(jj,jst)<=xuepsi .AND. hsnowmetamo=='T07' ) THEN   
        !    
     !   WRITE(*,*) csnowmetamo,': you are using T07 formulation!!'
        !
        ! Coefficients from Taillander et al. 2007
        zssa0 = 6./( xrholi*xvdiam6 ) * 10.
        !
        za  =  0.659*zssa0 - 27.2 * ( psnowtemp(jj,jst)-273.15-2.03 )      ! TG conditions
        zb  = 0.0961*zssa0 - 3.44 * ( psnowtemp(jj,jst)-273.15+1.90 )
        zc  = -0.341*zssa0 - 27.2 * ( psnowtemp(jj,jst)-273.15-2.03 )
        za2 =  0.629*zssa0 - 15.0 * ( psnowtemp(jj,jst)-273.15-11.2 )     ! ET conditions
        zb2 = 0.0760*zssa0 - 1.76 * ( psnowtemp(jj,jst)-273.15-2.96 )
        zc2 = -0.371*zssa0 - 15.0 * ( psnowtemp(jj,jst)-273.15-11.2 )
        !
        ! Compute SSA (method from Jacobi, 2010)
!   ZSSA = 6./(XRHOLI*PSNOWGRAN1(JJ,JST))*10.
!    ZSSA_t = (0.5+0.5*TANH(0.5*(ZGRADT-10.)))*(ZA-ZB*LOG(PSNOWAGE(JJ,JST)*24+EXP(ZC/ZB))) + &
!         (0.5-0.5*TANH(0.5*(ZGRADT-10.)))*(ZA2-ZB2*LOG(PSNOWAGE(JJ,JST)*24+EXP(ZC2/ZB2)))
!
!   ZSSA_t_dt = (0.5+0.5*TANH(0.5*(ZGRADT-10.)))*(ZA-ZB*LOG(PSNOWAGE(JJ,JST)*24+PTSTEP/3600.+EXP(ZC/ZB))) + &
!             (0.5-0.5*TANH(0.5*(ZGRADT-10.)))*(ZA2-ZB2*LOG(PSNOWAGE(JJ,JST)*24+PTSTEP/3600.+EXP(ZC2/ZB2)))
! 
!   ZSSA = ZSSA + (ZSSA_t_dt-ZSSA_t)
!
!   ZSSA = MAX(ZSSA,8.*10.)
!
!   PSNOWGRAN1(JJ,JST) = 6./(XRHOLI*ZSSA)*10.
        !
        ! Compute SSA (rate equation with Taylor series)
        zssa = 6./( xrholi*psnowgran1(jj,jst) ) * 10.
        !
        zdenom1 = (psnowage(jj,jst)*24.) + exp(zc/zb)
        zdenom2 = (psnowage(jj,jst)*24.) + exp(zc2/zb2)
        zfact1 = 0.5 + 0.5*tanh( 0.5*(zgradt-10.) )
        zfact2 = 0.5 - 0.5*tanh( 0.5*(zgradt-10.) )
        zssa = zssa + (ptstep/3600.) * (                   zfact1 * (-zb/zdenom1)    + zfact2 * (-zb2/zdenom2)   + &
                                        (ptstep/3600.) * ( zfact1 * (zb/zdenom1**2.) + zfact2 * (zb2/zdenom2**2.) ) * 1./2. )
        !ZSSA = ZSSA + (PTSTEP/3600.) * (  ZFACT1 * ZB /ZDENOM1 * ( 1./ZDENOM1 * (PTSTEP/3600.) * 1./2. - 1. ) + &
        !                                  ZFACT2 * ZB2/ZDENOM2 * ( 1./ZDENOM2 * (PTSTEP/3600.) * 1./2. - 1. ) )
         !                                      
        zssa = max( zssa, 8.*10. )
        !
        psnowgran1(jj,jst) = 6./( xrholi*zssa ) * 10.
        !
      !---------------------------------
      !    FLANNER et al. 2006 (F06)
      !
      ! -> Dry snow
      ! -> Evolution of optical diameter
      !---------------------------------
      ELSEIF ( psnowliq(jj,jst)<=xuepsi .AND. hsnowmetamo=='F06' )THEN   
        !
      !  WRITE(*,*) CSNOWMETAMO,': you are using F06 formulation!!'
        !
        ! XDRDT0(dens,gradT,T), XTAU(dens,gradT,T), XKAPPA(dens,gradT,T)
        ! dens: [1-8 <-> 50.-400. kg/m3]
        ! gradT: [1-31 <-> 0.-300. K/m]
        ! T: [1-11 <-> 223.15-273.15 K]
        !
        !  Select indices of density, temperature gradient and temperature
        idrho  = min( abs( int( (psnowrho(jj,jst)-25.)/50.        ) + 1 ), 8  )
        idgrad = min( abs( int( (zgradt-5.)/10.+2.                )     ), 31 )
        idtemp = min( abs( int( (psnowtemp(jj,jst)-225.65 )/5.+2. )     ), 11 )
        IF ( psnowtemp(jj,jst)<221. ) idtemp = 1
        !
        ! Compute SSA
        zoptr0 = xvdiam6/2. * 10.**6.
        zoptr  = psnowgran1(jj,jst)/2. * 10.**6.
        zdrdt  = xdrdt0(idrho,idgrad,idtemp) * &
                 ( xtau(idrho,idgrad,idtemp) / &
                   ( zoptr - zoptr0 + xtau(idrho,idgrad,idtemp) ) )**(1./xkappa(idrho,idgrad,idtemp))
        zoptr  = zoptr + zdrdt * ptstep/3600.
        zoptr  = min( zoptr, 3./(xrholi*2.) * 10.**6.)
        !
        psnowgran1(jj,jst) = zoptr * 2./10.**6.
        !
      ENDIF
      !
    ENDIF
    !
  ENDDO
  !
ENDDO      

!*    2. MISE A JOUR VARIABLES HISTORIQUES (CAS NON DENDRITIQUE).
!        UPDATE OF THE HISTORICAL VARIABLES
!        --------------------------------------------------------
DO jj = 1,SIZE(psnowrho,1)
  !
  DO jst = 1,inlvls_use(jj)
    !
    IF ( hsnowmetamo=='B92' ) THEN   
      !
      !non dendritic
      gnondendritic = ( psnowgran1(jj,jst)>=0. )
      IF ( gnondendritic ) THEN
        !faceted crystals
        gfaceted = ( psnowgran1(jj,jst)<xvsphe4 .AND. psnowhist(jj,jst)==0. )
        !spheric and liquid water
        gsphe_lw = ( xvgran1-psnowgran1(jj,jst)<xvsphe4 .AND. psnowliq(jj,jst)/psnowdz(jj,jst)>xvtelv1 )
      END IF
      !  
    ELSE
      !
      !non dendritic
      gnondendritic = ( psnowgran1(jj,jst)>=xvdiam6*(4.-psnowgran2(jj,jst))-xuepsi )
      IF ( gnondendritic ) THEN
        !faceted crystals
        gfaceted = ( psnowgran2(jj,jst)<xvsphe4/xvgran1 .AND. psnowhist(jj,jst)==0. )
        !spheric and liquid water
        gsphe_lw = ( xvsphe1-psnowgran2(jj,jst)<xvsphe4/xvgran1 .AND. psnowliq(jj,jst)/psnowdz(jj,jst)>xvtelv1 )
      END IF    
      !
    ENDIF
    !    
    IF ( gnondendritic ) THEN
      !
      IF ( gfaceted ) THEN
        !
        psnowhist(jj,jst) = nvhis1
        !
      ELSEIF ( gsphe_lw ) THEN 
        !
        IF (psnowhist(jj,jst)==0.)     psnowhist(jj,jst) = nvhis2
        IF (psnowhist(jj,jst)==nvhis1) psnowhist(jj,jst) = nvhis3
        !
      ELSEIF ( psnowtemp(jj,jst) < xtt ) THEN
        !
        IF(psnowhist(jj,jst)==nvhis2) psnowhist(jj,jst) = nvhis4
        IF(psnowhist(jj,jst)==nvhis3) psnowhist(jj,jst) = nvhis5
        !
      ENDIF
      !
    ENDIF
    !
  ENDDO
  !
ENDDO
!
IF (lhook) CALL dr_hook('SNOWCROMETAMO',1,zhook_handle)
!
END SUBROUTINE snowcrometamo  
!
!####################################################################
!####################################################################
SUBROUTINE set_thresh(PGRADT,PSNOWLIQ,PSPHE)
!
USE modd_snow_metamo, ONLY : xuepsi, xvgrat1
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: pgradt
REAL, INTENT(IN) :: psnowliq
REAL, INTENT(INOUT) :: psphe
!
IF ( psnowliq>xuepsi .OR. pgradt<xvgrat1 ) THEN
  psphe = min(1.,psphe)
ELSE 
  psphe = max(0.,psphe)
ENDIF
!
END SUBROUTINE set_thresh
!####################################################################
!####################################################################
SUBROUTINE get_gran(PTSTEP,PTELM,PGRAN)
!
USE modd_csts, ONLY : xpi
USE modd_snow_metamo, ONLY : xvtail1, xvtail2, nvdent1
!
IMPLICIT NONE
!
REAL, INTENT(IN)  :: ptstep, ptelm
REAL, INTENT(INOUT) :: pgran
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_GRAN',0,zhook_handle)
!
pgran = 2. * ( 3./(4.*xpi) * &
                  ( 4. * xpi/3. * (pgran/2.)**3 + &
                   ( xvtail1 + xvtail2 * ptelm**nvdent1 ) * ptstep ) )**(1./3.)
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_GRAN',1,zhook_handle)
!
END SUBROUTINE get_gran
!  
!####################################################################
!####################################################################
!####################################################################
!
SUBROUTINE snowcroalb(TPTIME,OGLACIER,                            &
                      palbedosc,pspectralalbedo,psnowdz,          &
                      psnowrho,ppermsnowfrac,                     &
                      psnowgran1_top,psnowgran2_top,psnowage_top, &
                      psnowgran1_bot,psnowgran2_bot,psnowage_bot, &
                      pps, pzenith, knlvls_use ,hsnowmetamo       ) 
!
!!    PURPOSE
!!    -------
!     Calculate the snow surface albedo. Use the method of original
!     Crocus which considers a specified spectral distribution of solar 
!     solar radiation (to be replaced by an input forcing when available)
!     In addition to original crocus, the top 2 surface snow layers are
!     considered in the calculation, using an arbitrary weighting, in order
!     to avoid time discontinuities due to layers agregation
!     Ageing depends on the presence of permanent snow cover 
!
USE modd_snow_par, ONLY : xansmax, xansmin,xaglamin, xaglamax, &
                          xvrpre1,xvrpre2,xvaging_noglacier,   &
                          xvaging_glacier, xvspec1,xvspec2,    &
                          xvspec3, xvw1,xvw2,xvd1,xvd2
 
USE modd_type_date_surf, ONLY : date_time                                                    
!
USE mode_snow3l
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(date_time), INTENT(IN)       :: tptime      ! current date and time
LOGICAL, INTENT(IN)               :: oglacier    ! True = Over permanent snow and ice, 
!                                                   initialise WGI=WSAT,
!                                                   Hsnow>=10m and allow 0.8<SNOALB<0.85
                                                 ! False = No specific treatment
REAL, DIMENSION(:), INTENT(IN)    :: psnowdz,ppermsnowfrac
!
REAL,DIMENSION(:,:), INTENT(IN)   :: psnowrho ! For now only the 2 first layers are required
!
REAL, DIMENSION(:), INTENT(INOUT) :: palbedosc
!  
REAL, DIMENSION(:,:), INTENT(OUT) :: pspectralalbedo   ! Albedo in the different spectral bands
!
REAL, DIMENSION(:), INTENT(IN)    :: psnowgran1_top,psnowgran2_top,psnowage_top, &
                                     psnowgran1_bot,psnowgran2_bot,psnowage_bot, pps 
INTEGER, DIMENSION(:), INTENT(IN) :: knlvls_use                
! 
REAL, DIMENSION(:), INTENT(IN)    :: pzenith ! solar zenith angle for future use
!
 CHARACTER(3),INTENT(IN)           :: hsnowmetamo ! metamorphism scheme
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),3) :: zalb_top, zalb_bot
!
REAL, DIMENSION(SIZE(PSNOWRHO,1))   :: zansmin, zansmax, zmin, zmax
REAL, DIMENSION(SIZE(PSNOWRHO,1))   :: zfac_top, zfac_bot
!
REAL, DIMENSION(SIZE(PALBEDOSC))  :: zvage1   
!
!REAL            :: ZAGE_NOW
!
INTEGER         :: jj   ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('SNOWCROALB',0,zhook_handle)
!
! 0. Initialize:
! ------------------
!
! PRINT*,XVAGING_NOGLACIER,XVAGING_GLACIER
IF ( oglacier ) THEN
  zansmin(:) = xaglamin * ppermsnowfrac(:) + xansmin * (1.0-ppermsnowfrac(:))
  zansmax(:) = xaglamax * ppermsnowfrac(:) + xansmax * (1.0-ppermsnowfrac(:))  
  zvage1(:)  = xvaging_glacier * ppermsnowfrac(:) + xvaging_noglacier * (1.0-ppermsnowfrac(:)) 
ELSE
  zansmin(:) = xansmin
  zansmax(:) = xansmax   
  zvage1(:)  = xvaging_noglacier
ENDIF
!
! ! ! ! ! !  date computation for ageing effects
! ! ! ! !   CALL GREGODSTRATI(TPTIME%TDATE%YEAR,TPTIME%TDATE%MONTH,TPTIME%TDATE%DAY,   &
! ! ! ! !                       TPTIME%TIME,ZAGE_NOW) 
!
! coherence control
! to remove when initialization routines will be updated
IF ( minval(psnowage_bot)<0. ) THEN
  CALL abor1_sfx('FATAL ERROR in SNOWCRO: Snow layer age inconsistent : check initialization routine. !')
END IF
!
! ! ! ! ! ! should be moved with other time controls to not compute MAXVAL(PSNOWAGE_TOP) at each time step
! ! ! ! ! IF ((ZAGE_NOW - MAXVAL(PSNOWAGE_TOP))<-0.001) THEN
! ! ! ! !       WRITE(*,*),"ZAGE_NOW=",ZAGE_NOW
! ! ! ! !       WRITE(*,*),"MAXVAL(PSNOWAGE_TOP)=",MAXVAL(PSNOWAGE_TOP)
! ! ! ! !       CALL ABOR1_SFX(&
! ! ! ! !       'FATAL ERROR in SNOWCRO: Snow layer date inconsistent with the current day !')
! ! ! ! ! END IF
!
DO jj=1, SIZE(palbedosc)
  !
  IF ( knlvls_use(jj)==0 ) THEN
    ! case with no snow on the ground 
    palbedosc(jj) = zansmin(jj)
  ELSE
    !
    CALL get_alb(jj,psnowrho(jj,1),pps(jj),zvage1(jj),psnowgran1_top(jj),&
                 psnowgran2_top(jj),psnowage_top(jj),zalb_top(jj,:),hsnowmetamo)
    !                  
!      IF (KNLVLS_USE(JJ)>=1) THEN
    IF ( knlvls_use(jj)>=2 ) THEN !modif ML
      ! second surface layer when it exists 
      !
      CALL get_alb(jj,psnowrho(jj,2),pps(jj),zvage1(jj),psnowgran1_bot(jj),&
                   psnowgran2_bot(jj),min(365.,psnowage_bot(jj)),zalb_bot(jj,:),hsnowmetamo)
      !
    ELSE
      ! when it does not exist, the second surface layer gets top layer albedo   
      zalb_bot(jj,1) = zalb_top(jj,1)
      zalb_bot(jj,2) = zalb_top(jj,2)
      zalb_bot(jj,3) = zalb_top(jj,3)
    ENDIF
    ! 
    ! computation of spectral albedo over 3 bands taking into account the respective
    ! depths of top layers 
    zmin(jj) = min( 1., psnowdz(jj)/xvd1 )
    zmax(jj) = max( 0., (psnowdz(jj)-xvd1)/xvd2 )
    zfac_top(jj) = xvw1 * zmin(jj) + xvw2 * min( 1., zmax(jj) )
    zfac_bot(jj) = xvw1 * ( 1. - zmin(jj) ) + xvw2 * ( 1. - min( 1., zmax(jj) ) )
    pspectralalbedo(jj,1) = zfac_top(jj) * zalb_top(jj,1) + zfac_bot(jj) * zalb_bot(jj,1) 
    pspectralalbedo(jj,2) = zfac_top(jj) * zalb_top(jj,2) + zfac_bot(jj) * zalb_bot(jj,2)
    pspectralalbedo(jj,3) = zfac_top(jj) * zalb_top(jj,3) + zfac_bot(jj) * zalb_bot(jj,3)
    !
    ! arbitrarily specified spectral distribution  
    ! to be changed when solar radiation distribution is an input variable 
    palbedosc(jj) = xvspec1 * pspectralalbedo(jj,1) + &
                    xvspec2 * pspectralalbedo(jj,2) + &
                    xvspec3 * pspectralalbedo(jj,3) 
    !    
  ENDIF ! end case with snow on the ground
  !
ENDDO ! end loop grid points
!
IF (lhook) CALL dr_hook('SNOWCROALB',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE snowcroalb
!####################################################################
SUBROUTINE get_alb(KJ,PSNOWRHO_IN,PPS_IN,PVAGE1,PSNOWGRAN1,PSNOWGRAN2,PSNOWAGE,PALB,&
                   hsnowmetamo)
!
USE modd_snow_par, ONLY : xalbice1, xalbice2, xalbice3,   &
                          xrhothreshold_ice,              &
                          xvalb2, xvalb3, xvalb4, xvalb5, &
                          xvalb6, xvalb7, xvalb8, xvalb9, &
                          xvalb10, xvalb11, xvdiop1,      &
                          xvrpre1, xvrpre2, xvpres1
!
USE mode_snow3l, ONLY : get_diam
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: kj
REAL, INTENT(IN) :: psnowrho_in, pps_in
REAL, INTENT(IN) :: pvage1
REAL, INTENT(IN) :: psnowgran1, psnowgran2, psnowage
REAL, DIMENSION(3), INTENT(OUT) :: palb
!
 CHARACTER(3),INTENT(IN)::hsnowmetamo
!
REAL :: zdiam, zdiam_sqrt
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_ALB',0,zhook_handle)
!
IF ( psnowrho_in<xrhothreshold_ice ) THEN
  ! Normal case (snow)
  CALL get_diam(psnowgran1,psnowgran2,zdiam,hsnowmetamo)
  zdiam_sqrt = sqrt(zdiam)
  palb(1) = min( xvalb2 - xvalb3*zdiam_sqrt, xvalb4 )
  palb(2) = max( 0.3, xvalb5 - xvalb6*zdiam_sqrt )
  zdiam   = min( zdiam, xvdiop1 )
  zdiam_sqrt = sqrt(zdiam)  
  palb(3) = max( 0., xvalb7*zdiam - xvalb8*zdiam_sqrt + xvalb9 ) 
 ! AGE CORRECTION ONLY FOR VISIBLE BAND
 
! ! ! ! !               PALB(1)=MAX(XVALB11,PALB(1)-MIN(MAX(PPS_IN/XVPRES1,XVRPRE1), &
! ! ! ! !                       XVRPRE2)*XVALB10*MIN(365.,ZAGE_NOW-PSNOWAGE)/PVAGE1)

  palb(1) = max( xvalb11, palb(1) - min( max(pps_in/xvpres1,xvrpre1), xvrpre2 ) * &
                   xvalb10 * psnowage / pvage1 )
ELSE
  ! Prescribed spectral albedo for surface ice
  palb(1) = xalbice1
  palb(2) = xalbice2
  palb(3) = xalbice3
ENDIF
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_ALB',1,zhook_handle)
!
END SUBROUTINE get_alb
!
!####################################################################
!####################################################################
SUBROUTINE snowcrorad(TPTIME, OGLACIER,                        &
                      psw_rad, psnowalb, psnowdz,              &
                      psnowrho, palb, pradsink, pradxs,        &
                      psnowgran1, psnowgran2, psnowage,pps,    &
                      pzenith, ppermsnowfrac,knlvls_use,       &
                      osnow_abs_zenith,hsnowmetamo) 
!
!!    PURPOSE
!!    -------
!     Calculate the transmission of shortwave (solar) radiation
!     through the snowpack (using a form of Beer's Law: exponential
!     decay of radiation with increasing snow depth).
!     Needs a first calculation of the albedo to stay coherent with
!     ISBA-ES ==> make sure to keep SNOWCRORAD coherent with SNOWCROALB
!
USE modd_snow_par, ONLY : xwcrn, xansmax, xansmin, xans_todry,          &
                          xsnowdmin, xans_t, xaglamin, xaglamax,        &
                          xd1, xd2, xd3, xx, xvspec1, xvspec2, xvspec3, &
                          xvbeta1, xvbeta2, xvbeta3, xvbeta4, xvbeta5
USE modd_type_date_surf, ONLY : date_time                      
!
USE mode_snow3l, ONLY : get_diam
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(date_time), INTENT(IN)       :: tptime      ! current date and time
LOGICAL, INTENT(IN)               :: oglacier   ! True = Over permanent snow and ice, 
!                                                   initialise WGI=WSAT,
!                                                   Hsnow>=10m and allow 0.8<SNOALB<0.85
                                                ! False = No specific treatment
!
REAL, DIMENSION(:), INTENT(IN)      :: psw_rad, psnowalb, palb,ppermsnowfrac
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowrho, psnowdz
!
LOGICAL, INTENT(IN)                 :: osnow_abs_zenith ! parametrization for polar regions (not physic but better results)
!                                                       ! default FALSE
 CHARACTER(3), INTENT(IN)            :: hsnowmetamo
!
REAL, DIMENSION(:), INTENT(OUT)     :: pradxs
!
REAL, DIMENSION(:,:), INTENT(OUT)   :: pradsink
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowgran1, psnowgran2, psnowage 
REAL, DIMENSION(:), INTENT(IN)      :: pps 
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use
REAL, DIMENSION(:), INTENT(IN)      :: pzenith 
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1))    :: zradtot
REAL, DIMENSION(SIZE(PSNOWRHO,1))    :: zalb_new
REAL, DIMENSION(SIZE(PSNOWRHO,1),3)  :: zalb !albedo 3 bands
REAL, DIMENSION(SIZE(PSNOWRHO,2))    :: zdiam
REAL, DIMENSION(SIZE(PSNOWRHO,2),3)  :: zbeta
REAL, DIMENSION(3) :: zopticalpath, zfact
REAL :: zprojlat
!
INTEGER :: jj,jst,jb   ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCRORAD',0,zhook_handle)
!
! 0. Initialization:
! ------------------
!
pradsink(:,:) = 0.
!  
! 1. Computation of the new albedo (see SNOWCROALB):
! -----------------------------------
!
 CALL snowcroalb(tptime,oglacier,                                             &
                 zalb_new,zalb,psnowdz(:,1),psnowrho(:,1:2),                  &
                 ppermsnowfrac,psnowgran1(:,1),psnowgran2(:,1),               &
                 psnowage(:,1),psnowgran1(:,2),psnowgran2(:,2),psnowage(:,2), &
                 pps, pzenith, knlvls_use, hsnowmetamo                        ) 
!
DO jj = 1,SIZE(psw_rad)
  !
  DO jst = 1,knlvls_use(jj)
    CALL get_diam(psnowgran1(jj,jst),psnowgran2(jj,jst),zdiam(jst),hsnowmetamo)
  ENDDO    ! end loop snow layers
  !
  ! 2. Extinction of net shortwave radiation
  ! ----------------------------------------
  ! First calculates extinction coefficients fn of grain size and density
  ! then calculates exctinction in the layer and increases optical path length
  !
  ! Coefficient for taking into account the increase of path length of rays
  ! in snow due to zenithal angle
  zprojlat = 1. / max( xuepsi, cos(pzenith(jj)) )
  !
  pradsink(jj,:) = -psw_rad(jj) * ( 1.-psnowalb(jj) ) / ( 1.-zalb_new(jj) )
  !
  !   Initialize optical depth  
  zopticalpath(1) = 0.
  zopticalpath(2) = 0.
  zopticalpath(3) = 0.
  !
  DO jst = 1,knlvls_use(jj)
    !
    zbeta(jst,1) = max( xvbeta1 * psnowrho(jj,jst) / sqrt(zdiam(jst)), xvbeta2 )
    zbeta(jst,2) = max( xvbeta3 * psnowrho(jj,jst) / sqrt(zdiam(jst)), xvbeta4 )
    zbeta(jst,3) = xvbeta5
    !
    zfact(:) = 0.
    DO jb = 1,3
      zopticalpath(jb) = zopticalpath(jb) + zbeta(jst,jb) * psnowdz(jj,jst)
      IF (osnow_abs_zenith) THEN
        !This formulation is incorrect but it compensate partly the fact that the albedo formulation does not account for zenithal angle
        zfact(jb) = (1.-zalb(jj,jb)) * exp( -zopticalpath(jb)*zprojlat)
      ELSE
        zfact(jb) = (1.-zalb(jj,jb)) * exp( -zopticalpath(jb) )
      ENDIF
    ENDDO
    !
    pradsink(jj,jst) = pradsink(jj,jst) * &
                       ( xvspec1*zfact(1) + xvspec2*zfact(2) + xvspec3*zfact(3) )
    !
  ENDDO    ! end loop snow layers
  !
  ! For thin snow packs, radiation might reach base of
  ! snowpack and the reflected energy can be absorbed by the bottom of snow layer:
  ! THIS PROCESS IS NOT SIMULATED
  !
  ! 4. Excess radiation not absorbed by snowpack (W/m2)JJ
  ! ----------------------------------------------------
  !
  pradxs(jj) = -pradsink( jj,knlvls_use(jj) )
  !
ENDDO    !end loop grid points
!
IF (lhook) CALL dr_hook('SNOWCRORAD',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE snowcrorad
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrothrm(PSNOWRHO,PSCOND,PSNOWTEMP,PPS,PSNOWLIQ, &
                       ocond_grain,ocond_yen                   ) 
!
!!    PURPOSE
!!    -------
!     Calculate snow thermal conductivity from
!     Sun et al. 1999, J. of Geophys. Res., 104, 19587-19579
!     (vapor) and Anderson, 1976, NOAA Tech. Rep. NWS 19 (snow).
!
!     Upon activation of flag OCOND_YEN, use the Yen (1981) formula for thermal conductivity
!     This formula was originally used in Crocus.
!
USE modd_csts, ONLY : xp00, xcondi, xrholw
USE modd_snow_par, ONLY : xsnowthrmcond1, xsnowthrmcond2, xsnowthrmcond_avap, &
                          xsnowthrmcond_bvap, xsnowthrmcond_cvap, xvrkz6
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)      :: pps
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowtemp, psnowrho, psnowliq
REAL, DIMENSION(:,:), INTENT(OUT)   :: pscond
LOGICAL, INTENT(IN)                 :: ocond_grain, ocond_yen
!
!*      0.2    declarations of local variables
!
INTEGER :: jj, jst ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROTHRM',0,zhook_handle)
!
! 1. Snow thermal conductivity
! ----------------------------
!
DO jst = 1,SIZE(psnowrho(:,:),2)
  !
  DO jj = 1,SIZE(psnowrho(:,:),1)
    !
    IF ( ocond_yen ) THEN
      pscond(jj,jst) = xcondi * exp( xvrkz6 * log( psnowrho(jj,jst)/xrholw ) )
    ELSE
      pscond(jj,jst) = ( xsnowthrmcond1 + &
                         xsnowthrmcond2 * psnowrho(jj,jst) * psnowrho(jj,jst) ) + &
                         max( 0.0, ( xsnowthrmcond_avap + &
                                    ( xsnowthrmcond_bvap/(psnowtemp(jj,jst) + xsnowthrmcond_cvap) ) ) &
                                   * (xp00/pps(jj)) ) 
    ENDIF
    !
    ! Snow thermal conductivity is set to be above 0.04 W m-1 K-1
    IF ( ocond_grain ) THEN                  
      pscond(jj,jst) = max( 0.04, pscond(jj,jst) )
      ! Snow thermal conductivity is annihilated in presence of liquid water            
      IF( psnowliq(jj,jst)>xuepsi ) pscond(jj,jst) = 0.01 * pscond(jj,jst)  
    ENDIF
     !
   ENDDO ! end loop JST
   !
ENDDO ! end loop JST
!
IF (lhook) CALL dr_hook('SNOWCROTHRM',1,zhook_handle)
!
END SUBROUTINE snowcrothrm
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcroebud(HSNOWRES, HIMPLICIT_WIND,                                   &
                       ppew_a_coef, ppew_b_coef,                                   &
                       ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef,         &
                       psnowdzmin,                                                 &
                       pzref,pts,psnowrho,psnowliq,pscap,pscond1,pscond2,          &
                       puref,pexns,pexna,pdircoszw,pvmod,                          &
                       plw_rad,psw_rad,pta,pqa,pps,ptstep,                         &
                       psnowdz1,psnowdz2,palbt,pz0,pz0eff,pz0h,                    &
                       psfcfrz,pradsink,phpsnow,                                   &
                       pct,pemist,prhoa,ptsterm1,ptsterm2,pra,pcdsnow,pchsnow,     &
                       pqsat,pdqsat,prsra,pustar2_ic, pri,                         &
                       ppet_a_coef_t,ppeq_a_coef_t,ppet_b_coef_t,ppeq_b_coef_t     ) 
!
!!    PURPOSE
!!    -------
!     Calculate surface energy budget linearization (terms) and turbulent
!     exchange coefficients/resistance between surface and atmosphere.
!     (Noilhan and Planton 1989; Giordani 1993; Noilhan and Mahfouf 1996)
!
!!    MODIFICATIONS
!!    -------------
!!      Original A. Boone 
!!      Modified by E. Brun (24/09/2012) :
!!      * Correction coupling coefficient for specific humidity in SNOWCROEBUD
!!      * PSFCFRZ(:)  = 1.0 for systematic solid/vapor latent fluxes in SNOWCROEBUD
!!      Modified by B. Decharme 09/12  new wind implicitation
!
USE modd_surf_par, ONLY : xundef
USE modd_csts,     ONLY : xcpd, xrholw, xstefan, xlvtt, xlstt, xrholw
USE modd_snow_par, ONLY : x_ri_max, xemissn
!
USE mode_thermos
!
USE modi_surface_ri
USE modi_surface_aero_cond
USE modi_surface_cd
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                   :: ptstep, psnowdzmin
!
 CHARACTER(LEN=*),     INTENT(IN)  :: hsnowres ! type of sfc resistance
!                                      DEFAULT=Louis (1979), standard ISBA
!                                      method. Option to limit Ri number
!                                      for very stable conditions
!
 CHARACTER(LEN=*),     INTENT(IN)  :: himplicit_wind   ! wind implicitation option
!                                                     ! 'OLD' = direct
!                                                     ! 'NEW' = Taylor serie, order 1
!
REAL, DIMENSION(:), INTENT(IN)     :: ppew_a_coef, ppew_b_coef,                   &
                                      ppet_a_coef, ppeq_a_coef, ppet_b_coef,      &
                                      ppeq_b_coef  
!                                      PPEW_A_COEF = wind coefficient (m2s/kg)
!                                      PPEW_B_COEF = wind coefficient (m/s)
!                                      PPET_A_COEF = A-air temperature coefficient
!                                      PPET_B_COEF = B-air temperature coefficient
!                                      PPEQ_A_COEF = A-air specific humidity coefficient
!                                      PPEQ_B_COEF = B-air specific humidity coefficient                                    
!
REAL, DIMENSION(:), INTENT(IN)     :: pzref, pts, psnowdz1, psnowdz2,        &
                                      pradsink, psnowrho, psnowliq, pscap,   &
                                      pscond1, pscond2,                      &
                                      pz0, phpsnow,                          &
                                      palbt, pz0eff, pz0h 
!
REAL, DIMENSION(:), INTENT(IN)     :: psw_rad, plw_rad, pta, pqa, pps, prhoa
!
REAL, DIMENSION(:), INTENT(IN)     :: puref, pexns, pexna, pdircoszw, pvmod
!
REAL, DIMENSION(:), INTENT(OUT)    :: ptsterm1, ptsterm2, pemist, pra,         &
                                      pct, psfcfrz, pcdsnow, pchsnow,          &
                                      pqsat, pdqsat, prsra 
!
REAL, DIMENSION(:), INTENT(OUT)    :: pustar2_ic,                      &
                                      ppet_a_coef_t, ppeq_a_coef_t,    &
                                      ppet_b_coef_t, ppeq_b_coef_t 
!
REAL, DIMENSION(:), INTENT(OUT)    :: pri
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PTS))        :: zac, zri,                        &
                                     zsconda, za, zb, zc,             &
                                     zcdn, zsnowdzm1, zsnowdzm2,      &
                                     zvmod, zustar2, zts3, zlvt,      &
                                     z_ccoef 
REAL, DIMENSION(SIZE(PTS))        :: zsnowevapx, zdenom, znumer   
!
INTEGER :: jj   ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
! 1. New saturated specific humidity and derrivative:
! ---------------------------------------------------
!
IF (lhook) CALL dr_hook('SNOWCROEBUD',0,zhook_handle)
!
zri(:) = xundef
!
pqsat(:) = qsati(pts(:),pps(:))
pdqsat(:) = dqsati(pts(:),pps(:),pqsat(:))
!
! 2. Surface properties:
! ----------------------
! It might be of interest to use snow-specific roughness
! or a temperature dependence on emissivity:
! but for now, use ISBA defaults.
!
pemist(:) = xemissn
!
! 2. Computation of resistance and drag coefficient
! -------------------------------------------------
!
 CALL surface_ri(pts, pqsat, pexns, pexna, pta, pqa,  &
                 pzref, puref, pdircoszw, pvmod, zri  ) 
!
! Simple adaptation of method by Martin and Lejeune (1998)
! to apply a lower limit to turbulent transfer coef
! by defining a maximum Richarson number for stable
! conditions:
!
IF ( hsnowres=='RIL' ) THEN
  DO jj=1,SIZE(zri)
    zri(jj) = min( x_ri_max, zri(jj) )
  ENDDO
ENDIF
!
pri(:) = zri(:)
!
! Surface aerodynamic resistance for heat transfers
!
 CALL surface_aero_cond(zri, pzref, puref, pvmod, pz0, pz0h, zac, pra, pchsnow)
!
prsra(:) = prhoa(:) / pra(:)
!
! For atmospheric model coupling:
!
 CALL surface_cd(zri, pzref, puref, pz0eff, pz0h, pcdsnow, zcdn)
!
!
! Modify flux-form implicit coupling coefficients:
! - wind components:
!
znumer(:) = pcdsnow(:)*pvmod(:)
zdenom(:) = prhoa(:) * pcdsnow(:) * pvmod(:) * ppew_a_coef(:)
IF(himplicit_wind=='OLD')THEN
! old implicitation
  zustar2(:) = ( znumer(:) * ppew_b_coef(:) ) / ( 1.0 - zdenom(:) ) 
ELSE                 
! new implicitation
  zustar2(:) = ( znumer(:) * ( 2.*ppew_b_coef(:) - pvmod(:) ) ) / ( 1.0 - 2.0*zdenom(:) )
ENDIF               
!
zvmod(:) = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
zvmod(:) = max( zvmod(:),0. )
!
WHERE ( ppew_a_coef(:)/= 0. )
  zustar2(:) = max( ( zvmod(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0. )
ENDWHERE
!
! implicit wind friction
zustar2(:) = max( zustar2(:),0. )
!
pustar2_ic(:) = zustar2(:)
!
! 3. Calculate linearized surface energy budget components:
! ---------------------------------------------------------
! To prevent numerical difficulties for very thin snow
! layers, limit the grid "thinness": this is important as
! layers become vanishing thin:
!
zsnowdzm1(:) = max( psnowdz1(:), psnowdzmin )
zsnowdzm2(:) = max( psnowdz2(:), psnowdzmin )
!
! Surface thermal inertia:
!
pct(:) = 1.0 / ( pscap(:)*zsnowdzm1(:) )
!
! Fraction of surface frozen (sublimation) with the remaining
! fraction being liquid (evaporation):
!
psfcfrz(:) = 1.0 
!
! Thermal conductivity between uppermost and lower snow layers:
!
zsconda(:) = ( zsnowdzm1(:)*pscond1(:) + zsnowdzm2(:)*pscond2(:) ) / &
             ( zsnowdzm1(:)            + zsnowdzm2(:)            ) 
!
! Transform implicit coupling coefficients: 
! Note, surface humidity is 100% over snow.
!
! - specific humidity:
!
z_ccoef(:) = 1.0 - ppeq_a_coef(:) * prsra(:)
!
ppeq_a_coef_t(:) = - ppeq_a_coef(:) * prsra(:) * pdqsat(:) / z_ccoef(:)
!
ppeq_b_coef_t(:) = ( ppeq_b_coef(:) &
                   - ppeq_a_coef(:) * prsra(:) * (pqsat(:) - pdqsat(:)*pts(:)) ) / z_ccoef(:)  
!
! - air temperature:
!   (assumes A and B correspond to potential T):
!
z_ccoef(:) = ( 1.0 - ppet_a_coef(:) * prsra(:) ) / pexna(:)
!
ppet_a_coef_t(:) = - ppet_a_coef(:) * prsra(:) / ( pexns(:) * z_ccoef(:) )
!
ppet_b_coef_t(:) = ppet_b_coef(:) / z_ccoef(:)
!
!
! Energy budget solution terms:
!
zts3(:) = pemist(:) * xstefan * pts(:)**3
zlvt(:) = (1.-psfcfrz(:))*xlvtt + psfcfrz(:)*xlstt
!
za(:) = 1./ptstep + pct(:) * &
         ( 4. * zts3(:) + prsra(:) * zlvt(:) * ( pdqsat(:) - ppeq_a_coef_t(:) )            &
                        + prsra(:) * xcpd * ( (1./pexns(:))-(ppet_a_coef_t(:)/pexna(:)) )  &
                        + ( 2*zsconda(:) / ( zsnowdzm2(:)+zsnowdzm1(:) ) ) )  
!
zb(:) = 1./ptstep + pct(:) * &
         ( 3. * zts3(:) + prsra(:) * zlvt(:) * pdqsat(:) )  
!
zc(:) = pct(:) * ( - prsra(:) * zlvt(:) * ( pqsat(:) - ppeq_b_coef_t(:) )  & 
                   + prsra(:) * xcpd * ppet_b_coef_t(:) / pexna(:)         &
                   + psw_rad(:) * (1. - palbt(:)) + pemist(:) * plw_rad(:) &
                   + phpsnow(:) + pradsink(:) )  
!
!
! Coefficients needed for implicit solution
! of linearized surface energy budget:
!
ptsterm2(:) = 2. * zsconda(:) * pct(:) / ( za(:) * (zsnowdzm2(:)+zsnowdzm1(:) ) )
!
ptsterm1(:) = ( pts(:)*zb(:) + zc(:) ) / za(:)
!
IF (lhook) CALL dr_hook('SNOWCROEBUD',1,zhook_handle)
!
END SUBROUTINE snowcroebud
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrosolvt(PTSTEP,PSNOWDZMIN,                     &
                        psnowdz,pscond,pscap,ptg,              &
                        psoilcond,pd_g,                        &
                        pradsink,pct,pterm1,pterm2,            &
                        ppet_a_coef_t,ppeq_a_coef_t,           &
                        ppet_b_coef_t,ppeq_b_coef_t,           &
                        pta_ic, pqa_ic,                        &
                        pgbas,psnowtemp,psnowflux,             &
                        knlvls_use                             ) 
!
!!    PURPOSE
!!    -------
!     This subroutine solves the 1-d diffusion of 'ZSNOWTEMP' using a
!     layer averaged set of equations which are time differenced
!     using the backward-difference scheme (implicit).
!     The eqs are solved rapidly by taking advantage of the
!     fact that the matrix is tridiagonal. This is a very
!     general routine and can be used for the 1-d diffusion of any
!     quantity as long as the diffusity is not a function of the
!     quantity being diffused. Aaron Boone 8-98. Soln to the eq:
!
!                 c  dQ    d  k dQ    dS
!                    -- =  --   -- -  --
!                    dt    dx   dx    dx
!
!     where k = k(x) (thermal conductivity), c = c(x) (heat capacity)
!     as an eg. for temperature/heat eq. S is a sink (-source) term.
!     Diffusivity is k/c
!
!!     MODIFICATIONS
!!     -------------
!!      Original A. Boone 
!!       05/2011: Brun  Special treatment to tackle the variable number
!!                      of snow layers
!
USE modd_csts, ONLY : xtt
!
USE modi_tridiag_ground_snowcro
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep, psnowdzmin
!
REAL, DIMENSION(:), INTENT(IN)      :: ptg, psoilcond, pd_g,        &
                                        pct, pterm1, pterm2 

!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowdz, pscond, pscap,      &
                                        pradsink 
!
REAL, DIMENSION(:), INTENT(IN)      :: ppet_a_coef_t, ppeq_a_coef_t, &
                                        ppet_b_coef_t, ppeq_b_coef_t 
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowtemp
!
REAL, DIMENSION(:), INTENT(OUT)     :: pgbas, psnowflux, pta_ic, pqa_ic
!
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: zsnowtemp, zdterm, zcterm, &
                                                    zfrcv, zamtrx, zbmtrx,     &
                                                    zcmtrx 
!
REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)) :: zwork1, zwork2, zdzdif,    &
                                                    zsnowdzm 
!
REAL, DIMENSION(SIZE(PSNOWDZ,1),SIZE(PSNOWDZ,2)-1) :: zsnowtemp_m,             &
                                                      zfrcv_m, zamtrx_m,       &
                                                      zbmtrx_m, zcmtrx_m 
!
REAL, DIMENSION(SIZE(PTG)) :: zgbas, zsnowtemp_delta
!
INTEGER :: jj, jst   ! looping indexes
INTEGER :: inlvls
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROSOLVT',0,zhook_handle)
!
! 0. Initialize:
! ------------------
!
zsnowtemp(:,:) = psnowtemp(:,:)
inlvls         = SIZE(psnowdz(:,:),2)
!
! 1. Calculate tri-diagnoal matrix coefficients:
! ----------------------------------------------
! For heat transfer, assume a minimum grid
! thickness (to prevent numerical
! problems for very thin snow cover):
!
DO jj=1,SIZE(ptg)
  !
  DO jst = knlvls_use(jj),1,-1
    !
    zsnowdzm(jj,jst) = max( psnowdz(jj,jst), psnowdzmin )
    !
    zwork1(jj,jst)   = zsnowdzm(jj,jst) * pscond(jj,jst) 
    !
    IF ( jst<knlvls_use(jj) ) THEN
      !
      zdzdif(jj,jst) = zsnowdzm(jj,jst) + zsnowdzm(jj,jst+1) 
      !
      zwork2(jj,jst) = zsnowdzm(jj,jst+1) * pscond(jj,jst+1)
      !
    ELSE
      !
      zdzdif(jj,jst) = zsnowdzm(jj,jst) + pd_g(jj)
      !
      zwork2(jj,jst) = pd_g(jj) * psoilcond(jj)
      !
    ENDIF
    !
    zdterm(jj,jst) = 2.0 * ( zwork1(jj,jst) + zwork2(jj,jst) ) / zdzdif(jj,jst)**2 
    !
    zcterm(jj,jst) = pscap(jj,jst) * zsnowdzm(jj,jst) / ptstep 
    !
  ENDDO
  !
ENDDO
!
! 2. Set up tri-diagonal matrix
! -----------------------------
!
zamtrx(:,:) = 0.
zbmtrx(:,:) = 0.
zcmtrx(:,:) = 0.
zfrcv(:,:) = 0.
! Upper BC
!
zamtrx(:,1) =  0.0
zbmtrx(:,1) =  1. / ( pct(:)*ptstep )
zcmtrx(:,1) = - pterm2(:) * zbmtrx(:,1)
zfrcv(:,1) =   pterm1(:) * zbmtrx(:,1)
!
DO jj = 1,SIZE(ptg)
  !
  DO jst = 2,knlvls_use(jj)
    !
    ! Interior Grid & Lower BC
    zamtrx(jj,jst) = -zdterm(jj,jst-1)
    zbmtrx(jj,jst) =  zcterm(jj,jst) + zdterm(jj,jst-1) + zdterm(jj,jst) 
    zfrcv(jj,jst) =  zcterm(jj,jst)*psnowtemp(jj,jst) - (pradsink(jj,jst-1)-pradsink(jj,jst))     
    !
    IF ( jst<knlvls_use(jj) ) THEN
      zcmtrx(jj,jst) = -zdterm(jj,jst)
    ELSE
      zcmtrx(jj,jst) = 0.0
      zfrcv(jj,jst) = zfrcv(jj,jst) + zdterm(jj,jst)*ptg(jj)
    ENDIF
    !
  ENDDO
  !
ENDDO
!
! 4. Compute solution vector
! --------------------------
!
 CALL tridiag_ground_snowcro(zamtrx,zbmtrx,zcmtrx,zfrcv,zsnowtemp,knlvls_use,0) 
!
! Heat flux between surface and 2nd snow layers: (W/m2)
!
psnowflux(:) = zdterm(:,1) * ( zsnowtemp(:,1) - zsnowtemp(:,2) )
!
! 5. Snow melt case
! -----------------
! If melting in uppermost layer, assume surface layer
! temperature at freezing point and re-evaluate lower
! snowpack temperatures. This is done as it is most likely
! most signigant melting will occur within a time step in surface layer.
! Surface energy budget (and fluxes) will
! be re-calculated (outside of this routine):
!
zamtrx_m(:,1) =  0.0
zbmtrx_m(:,1) =  zcterm(:,2) + zdterm(:,1) + zdterm(:,2)
zcmtrx_m(:,1) = -zdterm(:,2)
zfrcv_m(:,1) =  zcterm(:,2)*psnowtemp(:,2) - (pradsink(:,1)-pradsink(:,2)) + zdterm(:,1)*xtt 
!
DO jj = 1,SIZE(ptg)
  DO jst = 2,knlvls_use(jj)-1
    zamtrx_m(jj,jst) = zamtrx(jj,jst+1)
    zbmtrx_m(jj,jst) = zbmtrx(jj,jst+1)
    zcmtrx_m(jj,jst) = zcmtrx(jj,jst+1)
    zfrcv_m(jj,jst) = zfrcv(jj,jst+1)
    zsnowtemp_m(jj,jst) = psnowtemp(jj,jst+1)
  ENDDO
ENDDO
!
 CALL tridiag_ground_snowcro(zamtrx_m,zbmtrx_m,zcmtrx_m,zfrcv_m,zsnowtemp_m,knlvls_use,1) 
!
! If melting for 2 consecuative time steps, then replace current T-profile
! with one assuming T=Tf in surface layer:
!
zsnowtemp_delta(:) = 0.0
!
WHERE( zsnowtemp(:,1)>xtt .AND. psnowtemp(:,1)>=xtt )
  psnowflux(:) = zdterm(:,1) * ( xtt-zsnowtemp_m(:,1) )
  zsnowtemp_delta(:) = 1.0
END WHERE
!
DO jj = 1,SIZE(ptg)
  DO jst = 2,knlvls_use(jj)
    zsnowtemp(jj,jst) = zsnowtemp_delta(jj)       * zsnowtemp_m(jj,jst-1) + &
                        (1.0-zsnowtemp_delta(jj)) * zsnowtemp(jj,jst)
  ENDDO
ENDDO
!
! 6. Lower boundary flux:
! -----------------------
! NOTE: evaluate this term assuming the snow layer
! can't exceed the freezing point as this adjustment
! is made in melting routine. Then must adjust temperature
! to conserve energy:
!
DO jj=1, SIZE(ptg)
 zgbas(jj) = zdterm(jj,knlvls_use(jj)) * ( zsnowtemp(jj,knlvls_use(jj))             - ptg(jj) ) 
 pgbas(jj) = zdterm(jj,knlvls_use(jj)) * ( min( xtt, zsnowtemp(jj,knlvls_use(jj)) ) - ptg(jj) )
 zsnowtemp(jj,knlvls_use(jj)) = zsnowtemp(jj,knlvls_use(jj)) + &
                                ( zgbas(jj)-pgbas(jj) ) / zcterm(jj,knlvls_use(jj))      
ENDDO
!
! 7. Update temperatute profile in time:
! --------------------------------------
!
DO jj=1, SIZE(ptg)
  psnowtemp(jj,1:knlvls_use(jj)) = zsnowtemp(jj,1:knlvls_use(jj))
ENDDO
!
!
! 8. Compute new (implicit) air T and specific humidity
! -----------------------------------------------------
!
pta_ic(:) = ppet_b_coef_t(:) + ppet_a_coef_t(:) * psnowtemp(:,1)
!
pqa_ic(:) = ppeq_b_coef_t(:) + ppeq_a_coef_t(:) * psnowtemp(:,1)
!
IF (lhook) CALL dr_hook('SNOWCROSOLVT',1,zhook_handle)
!
END SUBROUTINE snowcrosolvt
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcromelt(PSCAP,PSNOWTEMP,PSNOWDZ,         &
                       psnowrho,psnowliq,knlvls_use     ) 
!
!!    PURPOSE
!!    -------
!     Calculate snow melt (resulting from surface fluxes, ground fluxes,
!     or internal shortwave radiation absorbtion). It is used to
!     augment liquid water content, maintain temperatures
!     at or below freezing, and possibly reduce the mass
!     or compact the layer(s).
!
USE modd_csts,ONLY : xtt, xlmtt, xrholw, xrholi
!
USE mode_snow3l
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pscap
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz, psnowtemp, psnowrho,   &
                                          psnowliq 
!
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use 
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zphase, zcmprsfact,   &
                                                      zsnowlwe,             &
                                                      zsnowmelt, zsnowtemp 
!
INTEGER :: jj, jst ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROMELT',0,zhook_handle)
!
! 0. Initialize:
! ---------------------------
!
DO jj = 1,SIZE(psnowdz,1)
   DO jst = 1,knlvls_use(jj)
     zphase(jj,jst) = 0.0
     zcmprsfact(jj,jst) = 0.0
     zsnowlwe(jj,jst) = 0.0
     zsnowmelt(jj,jst) = 0.0
     zsnowtemp(jj,jst) = 0.0
   ENDDO
ENDDO
!
! 1. Determine amount of melt in each layer:
! ------------------------------------------
!
DO jj = 1,SIZE(psnowdz,1)
  !
   DO jst = 1,knlvls_use(jj)
     !
     ! total liquid equivalent water content of snow(m):
     zsnowlwe(jj,jst) = psnowrho(jj,jst) * psnowdz(jj,jst) / xrholw
     !
     ! melt snow if excess energy and snow available:
     ! phase change(j/m2)
     zphase(jj,jst) = min( pscap(jj,jst) * max(0.0,psnowtemp(jj,jst)-xtt) * psnowdz(jj,jst), &
                           max(0.0,zsnowlwe(jj,jst)-psnowliq(jj,jst)) * xlmtt * xrholw ) 
     !
     ! update snow liq water content and temperature if melting:
     ! liquid water available for next layer from melting of snow
     ! which is assumed to be leaving the current layer(m):
     zsnowmelt(jj,jst) = zphase(jj,jst) / (xlmtt*xrholw)
     !
     ! cool off snow layer temperature due to melt:
     zsnowtemp(jj,jst) = psnowtemp(jj,jst) - zphase(jj,jst) / (pscap(jj,jst)*psnowdz(jj,jst))
     !
     ! difference with isba_es: zmeltxs should never be different of 0.
     ! because of the introduction of the tests in llayergone
     psnowtemp(jj,jst) =  zsnowtemp(jj,jst)
     !
     ! the control below should be suppressed after further tests
     IF (psnowtemp(jj,jst)-xtt > xuepsi) THEN
       WRITE(*,*) 'pb dans MELT PSNOWTEMP(JJ,JST) >XTT:', jj,jst, psnowtemp(jj,jst)
       CALL abor1_sfx('SNOWCRO: pb dans MELT')
     ENDIF
    !
    ! Loss of snowpack depth: (m) and liquid equiv (m):
    ! Compression factor for melt loss: this decreases
    ! layer thickness and increases density thereby leaving
    ! total SWE constant. 
    !
    ! Difference with ISBA_ES: All melt is considered to decrease the depth
    ! without consideration to the irreducible content
    !
    zcmprsfact(jj,jst) = ( zsnowlwe(jj,jst) - (psnowliq(jj,jst)+zsnowmelt(jj,jst)) ) &
                       / ( zsnowlwe(jj,jst) - psnowliq(jj,jst) )
    psnowdz(jj,jst) = psnowdz(jj,jst) * zcmprsfact(jj,jst)
    psnowrho(jj,jst) = zsnowlwe(jj,jst) * xrholw / psnowdz(jj,jst)
    !
    ! 2. Add snow melt to current snow liquid water content:
    ! ------------------------------------------------------
    !
    psnowliq(jj,jst) = psnowliq(jj,jst) + zsnowmelt(jj,jst)
    !
  ENDDO   ! loop JST active snow layers
ENDDO   ! loop JJ grid points 
!
IF (lhook) CALL dr_hook('SNOWCROMELT',1,zhook_handle)
!
END SUBROUTINE snowcromelt
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrorefrz(PTSTEP,PRR,                            &
                        psnowrho,psnowtemp,psnowdz,psnowliq,   &
                        pthrufal, pscap, plel3l,knlvls_use     ) 
!
!!    PURPOSE
!!    -------
!     Calculate any freezing/refreezing of liquid water in the snowpack.
!     Also, calculate liquid water transmission and snow runoff.
!     Refreezing causes densification of a layer.
!
USE modd_csts,     ONLY : xtt, xlmtt, xrholw, xci,xrholi
USE modd_snow_par, ONLY : xsnowdmin
!
USE mode_snow3l
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                      :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)        :: prr
!
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowdz, psnowtemp, psnowliq, psnowrho
!
REAL, DIMENSION(:), INTENT(INOUT)     :: pthrufal
!
! modifs_EB layers
INTEGER, DIMENSION(:), INTENT(IN)      :: knlvls_use
REAL, DIMENSION(:,:), INTENT(IN)       :: pscap
REAL, DIMENSION(:), INTENT(IN)         :: plel3l
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zphase,              &
                                                      zsnowliq, zsnowrho,  &
                                                      zwholdmax, zsnowdz,  &
                                                      zsnowtemp 
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),0:SIZE(PSNOWRHO,2)) :: zflowliq
!
REAL :: zdenom, znumer
!
INTEGER :: jj, jst   ! looping indexes
INTEGER :: inlvls     ! maximum snow layers number
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROREFRZ',0,zhook_handle)
!
! 0. Initialize:
! --------------
!
inlvls = SIZE(psnowdz,2)
!
DO jj=1,SIZE(psnowdz,1)        
  DO jst=1,knlvls_use(jj)
    zsnowrho(jj,jst) = psnowrho(jj,jst)
    zsnowtemp(jj,jst) = psnowtemp(jj,jst)
    zwholdmax(jj,jst) = snowcrohold( psnowrho(jj,jst),psnowliq(jj,jst),psnowdz(jj,jst) )
  ENDDO
ENDDO
!
DO jj = 1,SIZE(psnowdz,1)  ! loop JJ grid points        
  !
  ! 1. Increases Liquid Water of top layer from rain
  !    ---------------------------------------------
  !
  !  Rainfall (m) initialises the liquid flow whih feeds the top layer 
  !  and evaporation/condensation are taken into account
  ! 
  IF ( knlvls_use(jj)>0. ) THEN
    zflowliq(jj,0) = prr(jj) * ptstep / xrholw 
    zflowliq(jj,0) = max(0., zflowliq(jj,0) - plel3l(jj)*ptstep/(xlvtt*xrholw)) 
  ELSE
    zflowliq(jj,0) = 0
  ENDIF
  !
  DO jst=1,knlvls_use(jj) ! loop JST active snow layers
    !
    ! 2. Increases Liquid Water from the upper layers flow (or rain for top layer) 
    !    -----------------------------
    psnowliq(jj,jst) = psnowliq(jj,jst) + zflowliq(jj,jst-1)
    !                        
    ! 3. Freezes liquid water in any cold layers
    !    ---------------------------------------
    !                    
    ! Calculate the maximum possible refreezing 
    zphase(jj,jst) = min( pscap(jj,jst)* max(0.0, xtt - zsnowtemp(jj,jst)) * psnowdz(jj,jst), &
                          psnowliq(jj,jst) * xlmtt * xrholw ) 
    !
    ! Reduce liquid content if freezing occurs:
    zsnowliq(jj,jst) = psnowliq(jj,jst) - zphase(jj,jst)/(xlmtt*xrholw)       
    !
    ! Warm layer and reduce liquid if freezing occurs:
    zsnowdz(jj,jst) = max(xsnowdmin/inlvls, psnowdz(jj,jst))
    !
    !
    ! Difference with ISBA-ES: a possible cooling of current refreezing water
    !                          is taken into account to calculate temperature change
    znumer = ( zsnowrho(jj,jst) * zsnowdz(jj,jst) - ( psnowliq(jj,jst) - zflowliq(jj,jst-1) ) * xrholw )
    zdenom = ( zsnowrho(jj,jst) * zsnowdz(jj,jst) - ( zsnowliq(jj,jst) - zflowliq(jj,jst-1) ) * xrholw )
    !
    psnowtemp(jj,jst) = xtt + ( zsnowtemp(jj,jst)-xtt )*znumer/zdenom + zphase(jj,jst)/( xci*zdenom ) 
    !
    ! 4. Calculate flow from the excess of holding capacity
    !    --------------------------------------------------------------
    !
    ! Any water in excess of the maximum holding space for liquid water
    ! amount is drained into next layer down.
    zflowliq(jj,jst) = max( 0., zsnowliq(jj,jst)-zwholdmax(jj,jst) )
    !
    zsnowliq(jj,jst) = zsnowliq(jj,jst) - zflowliq(jj,jst)
    !
    ! 5. Density is adjusted to conserve the mass
    !    --------------------------------------------------------------
    znumer = ( zsnowrho(jj,jst) * psnowdz(jj,jst) - ( zflowliq(jj,jst) - zflowliq(jj,jst-1) ) * xrholw )
    !
    zsnowrho(jj,jst) = znumer / zsnowdz(jj,jst) 
    !
    ! keeps snow denisty below ice density
    IF ( zsnowrho(jj,jst)>xrholi ) THEN
      psnowdz(jj,jst) = psnowdz(jj,jst) * zsnowrho(jj,jst) / xrholi
      zsnowrho(jj,jst) = xrholi
    ENDIF
    !
    ! 6. Update thickness and density and any freezing:
    !    ----------------------------------------------
    psnowrho(jj,jst) = zsnowrho(jj,jst)
    psnowliq(jj,jst) = zsnowliq(jj,jst)
    !
  ENDDO ! loop JST active snow layers
  !
  ! Any remaining throughflow after freezing is available to
  ! the soil for infiltration or surface runoff (m).
  ! I.E. This is the amount of water leaving the snowpack:
  ! Rate water leaves the snowpack [kg/(m2 s)]:
  !
  pthrufal(jj)  = pthrufal(jj) + zflowliq(jj,knlvls_use(jj)) * xrholw / ptstep
  !
ENDDO ! loop JJ grid points
!
IF (lhook) CALL dr_hook('SNOWCROREFRZ',1,zhook_handle)
!
END SUBROUTINE snowcrorefrz
!####################################################################
SUBROUTINE get_rho(PRHO_IN,PDZ,PSNOWLIQ,PFLOWLIQ,PRHO_OUT)
!
USE modd_csts,     ONLY : xrholw
!
IMPLICIT NONE
!
REAL, INTENT(IN)  :: prho_in, pdz, psnowliq,pflowliq
REAL, INTENT(OUT) :: prho_out 
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_RHO',0,zhook_handle)
!
prho_out =  ( prho_in * pdz - ( psnowliq - pflowliq ) * xrholw )
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_RHO',1,zhook_handle)
!
END SUBROUTINE get_rho
!####################################################################
!####################################################################
SUBROUTINE snowcroflux(PSNOWTEMP,PSNOWDZ,PEXNS,PEXNA,          &
                       pustar2_ic,                             &
                       ptstep,palbt,psw_rad,pemist,plwupsnow,  &
                       plw_rad,pta,psfcfrz,pqa,phpsnow,        &
                       psnowtempo1,psnowflux,pct,pradsink,     &
                       pqsat,pdqsat,prsra,                     &
                       prn,ph,pgflux,ples3l,plel3l,pevap,      &
                       pustar                                  ) 
!
!!    PURPOSE
!!    -------
!     Calculate the surface fluxes (atmospheric/surface).
!     (Noilhan and Planton 1989; Noilhan and Mahfouf 1996)
!
USE modd_csts,ONLY : xtt
!
USE mode_thermos
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)      :: psnowdz, psnowtempo1, psnowflux, pct, &
                                        pradsink, pexns, pexna 
!
REAL, DIMENSION(:), INTENT(IN)      :: palbt, psw_rad, pemist, plw_rad,      &
                                        pta, psfcfrz, pqa,                    &
                                        phpsnow, pqsat, pdqsat, prsra,        &
                                        pustar2_ic 
!
REAL, DIMENSION(:), INTENT(INOUT)   :: psnowtemp
!
REAL, DIMENSION(:), INTENT(OUT)     :: prn, ph, pgflux, ples3l, plel3l,      &
                                        pevap, plwupsnow, pustar 
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWDZ))      :: zevapc, zsnowtemp
REAL :: zsmsnow, zgflux
!
INTEGER :: jj
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROFLUX',0,zhook_handle)
!
! 0. Initialize:
! --------------
!
! 1. Flux calculations when melt not occuring at surface (W/m2):
! --------------------------------------------------------------
!
DO jj = 1,SIZE(palbt)
  !
  CALL get_flux(palbt(jj),pemist(jj),psw_rad(jj),plw_rad(jj),   &
                pexns(jj),pexna(jj),pta(jj),pqa(jj),prsra(jj), &
                pqsat(jj),pdqsat(jj),psfcfrz(jj),phpsnow(jj),   &
                psnowtemp(jj),psnowtempo1(jj),                  &
                prn(jj),ph(jj),zevapc(jj),                      &
                ples3l(jj),plel3l(jj),zgflux                    )
  !
  IF ( psnowtemp(jj)>xtt ) THEN
    !
    IF ( psnowtempo1(jj)<xtt ) THEN
      !
      ! 2. Initial melt adjustment
      ! --------------------------
      ! If energy avalabile to melt snow, then recalculate fluxes
      ! at the freezing point and add residual heat to layer
      ! average heat.
      !
      ! A) If temperature change is > 0 and passes freezing point this timestep,
      !    then recalculate fluxes at freezing point and excess energy
      !    will be used outside of this routine to change snow heat content:
      !
      ! WRITE (*,*) 'attention test LFLUX traitement XTT supprime!'   
      !  
      CALL get_flux(palbt(jj),pemist(jj),psw_rad(jj),plw_rad(jj),   &
                    pexns(jj),pexna(jj), pta(jj),pqa(jj),prsra(jj), &
                    pqsat(jj),pdqsat(jj),psfcfrz(jj),phpsnow(jj),   &
                    xtt,psnowtempo1(jj),                            &
                    prn(jj),ph(jj),zevapc(jj),                      &
                    ples3l(jj),plel3l(jj),pgflux(jj)                )
      !
      zsmsnow = zgflux - pgflux(jj)
      !
      ! This will be used to change heat content of snow:
      zsnowtemp(jj) = psnowtemp(jj) - zsmsnow * ptstep * pct(jj)
      !
    ELSE
      !
      ! 3. Ongoing melt adjustment: explicit solution
      ! ---------------------------------------------
      !    If temperature change is 0 and at freezing point this timestep,
      !    then recalculate fluxes and surface temperature *explicitly*
      !    as this is *exact* for snow at freezing point (Brun, Martin)
      !
      CALL get_flux(palbt(jj),pemist(jj),psw_rad(jj),plw_rad(jj),   &
                    pexns(jj),pexna(jj), pta(jj),pqa(jj),prsra(jj), &
                    pqsat(jj),pdqsat(jj),psfcfrz(jj),phpsnow(jj),   &
                    xtt,xtt,                                        &
                    prn(jj),ph(jj),zevapc(jj),                      &
                    ples3l(jj),plel3l(jj),pgflux(jj)                )
      !
      zsnowtemp(jj) = xtt + ptstep * pct(jj) * ( pgflux(jj) + pradsink(jj) - psnowflux(jj) )
      !
    ENDIF
    !
  ELSE
    !
    zsnowtemp(jj) = psnowtemp(jj)
    !
    pgflux(jj) = zgflux
    !
  ENDIF
  !
ENDDO
!
! 4. Update surface temperature:
! ------------------------------
!
psnowtemp(:) = zsnowtemp(:)
!
! 5. Final evaporative flux (kg/m2/s)
!
pevap(:) = zevapc(:)
 !WRITE(*,*) 'Flux surface:',PGFLUX(1),PRN(1),PH(1), ZLE(1), PHPSNOW(1)
!
! 5. Friction velocity
! --------------------
!
pustar(:) = sqrt(pustar2_ic(:)) 
!
IF (lhook) CALL dr_hook('SNOWCROFLUX',1,zhook_handle)
!
END SUBROUTINE snowcroflux
!####################################################################
SUBROUTINE get_flux(PALBT,PEMIST,PSW_RAD,PLW_RAD,PEXNS,PEXNA,   &
                    pta,pqa,prsra,pqsat,pdqsat,psfcfrz,phpsnow, &
                    psnowtemp,psnowtempo1,                      &
                    prn,ph,pevapc,ples3l,plel3l,pgflux          )
!
USE modd_csts,ONLY : xstefan, xcpd, xlstt, xlvtt
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: palbt, pemist
REAL, INTENT(IN) :: psw_rad, plw_rad
REAL, INTENT(IN) :: pexns, pexna
REAL, INTENT(IN) :: pta, pqa, prsra, pqsat, pdqsat, psfcfrz, phpsnow
REAL, INTENT(IN) :: psnowtemp,psnowtempo1
REAL, INTENT(OUT):: prn, ph, pevapc, ples3l, plel3l, pgflux
!
REAL :: zle, zdeltat, zlwupsnow, zsnowto3
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_FLUX',0,zhook_handle)
!
zsnowto3  = psnowtempo1**3  ! to save some CPU time, store this
!
zdeltat   = psnowtemp - psnowtempo1   ! surface T time change:
!
zlwupsnow = pemist * xstefan * zsnowto3 * ( psnowtempo1 + 4.*zdeltat )
!
prn       = ( 1.-palbt )*psw_rad + pemist*plw_rad - zlwupsnow
!
ph        = prsra * xcpd * ( psnowtemp/pexns - pta/pexna )
!
pevapc    = prsra * ( (pqsat - pqa) + pdqsat*zdeltat )
!
ples3l    = psfcfrz      * xlstt * pevapc
!
plel3l    = (1.-psfcfrz) * xlvtt * pevapc
!
zle       = ples3l + plel3l
!
pgflux    = prn - ph - zle + phpsnow
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_FLUX',1,zhook_handle)
!
END SUBROUTINE get_flux
!
!####################################################################
!####################################################################
SUBROUTINE snowcroevapn(PLES3L,PTSTEP,PSNOWTEMP,PSNOWRHO, &
                       psnowdz,pevapcor,psnowhmass        ) 
!
!!    PURPOSE
!!    -------
!     Remove mass from uppermost snow layer in response to
!     evaporation (liquid) and sublimation.
!
!!     MODIFICATIONS
!!     -------------
!!      Original A. Boone 
!!      05/2011: E. Brun  Takes only into account sublimation and solid
!!                         condensation. Evaporation and liquid condensation
!!                         are taken into account in SNOWCROREFRZ
!
USE modd_csts,     ONLY : xlstt, xlmtt, xci, xtt
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)      :: psnowtemp
!
REAL, DIMENSION(:), INTENT(IN)      :: ples3l   ! (W/m2)
!
REAL, DIMENSION(:), INTENT(INOUT)   :: psnowrho, psnowdz, psnowhmass, &
                                        pevapcor 
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PLES3L))       :: zsnowevaps, zsnowevap, zsnowevapx,          &
                                       zsnowdz, zevapcor 
!                                      ZEVAPCOR = for vanishingy thin snow cover,
!                                                 allow any excess evaporation
!                                                 to be extracted from the soil
!                                                 to maintain an accurate water
!                                                 balance [kg/(m2 s)]
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROEVAPN',0,zhook_handle)
!
! 0. Initialize:
! --------------
!
zevapcor(:) = 0.0
zsnowevaps(:) = 0.0
zsnowevap(:) = 0.0
zsnowevapx(:) = 0.0
zsnowdz(:) = 0.0
!
WHERE ( psnowdz>0.0 )
  !
  ! 1. Sublimation/condensation of snow ice
  ! ----------------------------------------
  ! Reduce layer thickness and total snow depth
  ! if sublimation: add to correction term if potential
  ! sublimation exceeds available snow cover.
  !
  zsnowevaps(:) = ples3l(:) * ptstep / ( xlstt*psnowrho(:) )
  zsnowdz(:)    = psnowdz(:) - zsnowevaps(:)
  psnowdz(:)    = max( 0.0, zsnowdz(:) )
  zevapcor(:)   = zevapcor(:) + max(0.0,-zsnowdz(:)) * psnowrho(:) / ptstep
  !
  ! Total heat content change due to snowfall and sublimation (added here):
  ! (for budget calculations):
  !
  psnowhmass(:) = psnowhmass(:) &
                  - ples3l(:) * (ptstep/xlstt) * ( xci * (psnowtemp(:)-xtt) - xlmtt )
  ! 
END WHERE
!
! 3. Update evaporation correction term:
! --------------------------------------
!
pevapcor(:) = pevapcor(:) + zevapcor(:)
!
IF (lhook) CALL dr_hook('SNOWCROEVAPN',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE snowcroevapn
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrogone(PTSTEP,PLEL3L,PLES3L,PSNOWRHO,                     &
                       psnowheat,pradsink_2d,pevapcor,pthrufal,pgrndflux, &
                       pgfluxsnow,psnowdz,psnowliq,psnowtemp,pradxs,      &
                       prr,knlvls_use                                     )
!
!!    PURPOSE
!!    -------
!     Account for the case when the last trace of snow melts
!     during a time step: ensure mass and heat balance of
!     snow AND underlying surface.
!     Original A. Boone 
!     05/2011: E. Brun  Takes into account sublimation and PGRNDFLUX 
!                       Adds rain and evaporation/liquid condensation
!                       in PTHRUFAL
!
USE modd_csts,ONLY : xtt, xlstt, xlvtt
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)      :: plel3l, ples3l, pgfluxsnow,prr
!  
REAL, DIMENSION(:,:), INTENT(IN)    :: pradsink_2d
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psnowrho, psnowheat
!
REAL, DIMENSION(:), INTENT(INOUT)   :: pgrndflux, pradxs
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowdz, psnowliq, psnowtemp
!
REAL, DIMENSION(:), INTENT(OUT)     :: pthrufal   ! melt water [kg/(m2 s)]
!
REAL, DIMENSION(:), INTENT(OUT)     :: pevapcor   ! [kg/(m2 s)]
!                                      PEVAPCOR = for vanishingy thin snow cover,
!                                                 allow any excess evaporation
!                                                 to be extracted from the soil
!                                                 to maintain an accurate water
!                                                 balance.
!
INTEGER, DIMENSION(:), INTENT(INOUT)      :: knlvls_use
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PLES3L))       :: zradsink
REAL, DIMENSION(SIZE(PLES3L))       :: zsnowheatc
INTEGER, DIMENSION(SIZE(PLES3L))    :: isnowgone_delta
!
INTEGER  :: jj
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROGONE',0,zhook_handle)
!
! 0. Initialize:
! --------------
!
pevapcor(:) = 0.0
pthrufal(:) = 0.0
!
DO jj = 1,SIZE(zradsink)
  zradsink(jj) = pradsink_2d(jj,inlvls_use(jj)) 
  zsnowheatc(jj) = sum(psnowheat(jj,1:inlvls_use(jj))) !total heat content (J m-2)
END DO
!
isnowgone_delta(:) = 1
!
! 1. Simple test to see if snow vanishes:
! ---------------------------------------
! If so, set thicknesses (and therefore mass and heat) and liquid content
! to zero, and adjust fluxes of water, evaporation and heat into underlying
! surface.
!
! takes into account the heat content corresponding to the occasional
! sublimation  and then PGRNDFLUX
!
zsnowheatc(:) = zsnowheatc(:) + max( 0., ples3l(:)*ptstep/xlstt ) * xlmtt 
!
WHERE ( pgfluxsnow(:)+zradsink(:)-pgrndflux(:) >= (-zsnowheatc(:)/ptstep) )
  pgrndflux(:)       = pgfluxsnow(:) + (zsnowheatc(:)/ptstep)
  pevapcor(:)       = ples3l(:)/xlstt
  pradxs(:)       = 0.0
  isnowgone_delta(:) = 0          ! FLAG...if=0 then snow vanishes, else=1
END WHERE
!
! 2. Final update of snow state and computation of corresponding flow
!    Only if snow vanishes
! -----------------------------
!
pthrufal(:) = 0.
!
DO jj=1, SIZE(zradsink)
  !
  IF(isnowgone_delta(jj) == 0 ) THEN
    pthrufal(jj) = pthrufal(jj) + &
                   sum( psnowrho(jj,1:inlvls_use(jj))*psnowdz(jj,1:inlvls_use(jj)) ) / ptstep 
! takes into account rain and condensation/evaporation
    pthrufal(jj) = pthrufal(jj) + prr(jj) - plel3l(jj)/xlvtt
    psnowtemp(jj,:) = xtt
    psnowdz(jj,:) = 0.
    psnowliq(jj,:) = 0.
    inlvls_use(jj) = 0
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('SNOWCROGONE',1,zhook_handle)
!
END SUBROUTINE snowcrogone
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcroevapgone(PSNOWHEAT,PSNOWDZ,PSNOWRHO,PSNOWTEMP,PSNOWLIQ,      &
                           psnowgran1,psnowgran2,psnowhist,psnowage,knlvls_use,&
                           hsnowmetamo) 
!
!!    PURPOSE
!!    -------
!
!     If all snow in uppermost layer evaporates/sublimates, re-distribute
!     grid (below assumes very thin snowpacks so layer-thicknesses are
!     constant).
!     Original A. Boone 
!     05/2011: E. Brun  Takes into account previous changes in the energy
!                       content
!
!
USE modd_csts,     ONLY : xtt, xrholw, xlmtt, xci
USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, xrhosmax_es
USE mode_snow3l
USE modd_snow_metamo
USE modd_type_date_surf
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowrho   ! snow density profile                (kg/m3)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowdz    ! snow layer thickness profile        (m)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowheat  ! snow heat content/enthalpy          (J/m2)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowgran1 ! snow grain parameter 1              (-)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowgran2 ! snow grain parameter 2              (-)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowhist  ! snow grain historical variable      (-)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowage  ! Snow grain age
! 
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowtemp  ! snow temperature profile            (K)
REAL, DIMENSION(:,:), INTENT(INOUT)   :: psnowliq   ! snow liquid water profile           (m)
! 
INTEGER, DIMENSION(:), INTENT(IN)      :: knlvls_use
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo ! metamorphism scheme
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnowheat_1d ! total heat content                (J/m2)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnowrho_1d  ! total snowpack average density    (kg/m3)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnow        ! total snow depth                  (m)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zscap        ! Snow layer heat capacity          (J/K/m3)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zndent       ! Number of dendritic layers        (-)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: znvieu       ! Number of non dendritic layers    (-)
REAL, DIMENSION(SIZE(PSNOWDZ,1)) :: zsnowage_1d  ! total snowpack average 
!age (days)
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowgran1n,           &
                                                      zsnowgran2n,zsnowhistn 
!
LOGICAL :: gdendritic
!
INTEGER :: jj, jst          ! loop control
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWCROEVAPGONE',0,zhook_handle)
!
! Initialize:
!
zsnowheat_1d(:) = 0.
zsnow(:)        = 0.
zsnowrho_1d(:)  = 0.
zndent(:)       = 0.
znvieu(:)       = 0.  
zsnowage_1d(:)  = 0.
zscap(:)        = 0.
!
! First, determine where uppermost snow layer has completely
! evaporated/sublimated (as it becomes thin):
DO jj = 1,SIZE(psnowrho,1)
   !
   IF ( psnowdz(jj,1)==0.0 ) THEN
     !
     DO jst = 2,knlvls_use(jj)
       !
       zsnowheat_1d(jj) = zsnowheat_1d(jj) + psnowdz(jj,jst) * &
                          ( psnowrho(jj,jst)*xci * (zsnowtemp(jj,jst)-xtt) &
                            - xlmtt * psnowrho(jj,jst) ) &
                          + xlmtt * xrholw * psnowliq(jj,jst) 
       zsnow(jj) = zsnow(jj) + psnowdz(jj,jst)
       zsnowrho_1d(jj) = zsnowrho_1d(jj) + psnowdz(jj,jst) * psnowrho(jj,jst)     
       zsnowage_1d(jj) = zsnowage_1d(jj) + psnowdz(jj,jst) * psnowrho(jj,jst) * psnowage(jj,jst) 
       !  
       ! snow grains
       IF ( hsnowmetamo=='B92' ) THEN 
         gdendritic = ( psnowgran1(jj,jst)<-xepsi )
       ELSE
         gdendritic = ( psnowgran1(jj,jst)<xvdiam6*(4.-psnowgran2(jj,jst))-xuepsi )
       ENDIF
       !
       IF ( gdendritic ) THEN   ! Dendritic snow
         zndent(jj) = zndent(jj) + 1.0
       ELSE                                    ! Non dendritic snow
         znvieu(jj) = znvieu(jj) + 1.0
       ENDIF
       !
     ENDDO
     !
   ENDIF
   !
END DO
!
zsnowrho_1d(:) = zsnowrho_1d(:) / max( xsnowdmin, zsnow(:) )
zsnowage_1d(:) = zsnowage_1d(:) / max( xsnowdmin, zsnow(:) * zsnowrho_1d(:) )
zsnowrho_1d(:) = max( xrhosmin_es, min( xrhosmax_es, zsnowrho_1d(:) ) )
!
! Where uppermost snow layer has vanished, redistribute vertical
! snow mass and heat profiles (and associated quantities):
!
 CALL snow3lavgrain(psnowgran1,psnowgran2,psnowhist,                 &
                    zsnowgran1n,zsnowgran2n,zsnowhistn,zndent,znvieu,&
                    hsnowmetamo)         
!
DO jj=1,SIZE(psnowrho,1)
  !
  IF( zsnow(jj)/=0.0 ) THEN
    !
    psnowdz(jj,1:knlvls_use(jj)) = zsnow(jj) / knlvls_use(jj)
    psnowheat(jj,1:knlvls_use(jj)) = zsnowheat_1d(jj) / knlvls_use(jj)
    psnowrho(jj,1:knlvls_use(jj)) = zsnowrho_1d(jj)
    !
    zscap(jj) = zsnowrho_1d(jj) * xci
    !
    DO jst = 1,knlvls_use(jj)
      !
      psnowtemp(jj,jst) = xtt + ( ( (psnowheat(jj,jst)/psnowdz(jj,jst))   &
                                    + xlmtt*psnowrho(jj,jst) ) / zscap(jj) ) 
      psnowtemp(jj,jst) = min( xtt, psnowtemp(jj,jst) )  
      !      
      psnowliq(jj,jst) = max( 0.0, psnowtemp(jj,jst)-xtt ) * zscap(jj) * &
                                    psnowdz(jj,jst) / (xlmtt*xrholw) 
      !
    ENDDO
    !
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('SNOWCROEVAPGONE',1,zhook_handle)
!
END SUBROUTINE snowcroevapgone
!
!####################################################################
!####################################################################
!####################################################################
SUBROUTINE snownlfall_upgrid(TPTIME, OGLACIER,PTSTEP,PSR,PTA,PVMOD,        &
                             psnow,psnowrho,psnowdz,psnowheat,psnowhmass,  &
                             psnowalb,ppermsnowfrac,psnowgran1,psnowgran2, &
                             gsnowfall,psnowdzn,psnowrhof,psnowdzf,        &
                             psnowgran1f,psnowgran2f,psnowhistf,psnowagef, &
                             omodif_grid,knlvls_use,osnowdrift,pz0eff,puref,&
                             hsnowmetamo) 
!
!!    PURPOSE
!!    -------
! Adds new snowfall and updates the vertical grid in order to keep an
! optimal discertisation
!
!!    AUTHOR
!!    ------
!!      E. Brun           * Meteo-France *
!!
!
!!
!!     MODIFICATIONS
!!    ------
!!    
!!     2014-02-05 V. Vionnet: wind speed in the parameterization for new snow
!!                            density and characteristic of grains of new snow 
!!                            are taken at a reference height
!!     2014-06-03 M. Lafaysse : threshold on PZ0EFF
!!
USE modd_type_date_surf,  ONLY: date_time
USE modd_csts,     ONLY : xlmtt, xtt, xci
USE modd_snow_metamo, ONLY : xnden1, xnden2, xnden3, xgran, &
                             xnsph1, xnsph2, xnsph3, xnsph4
!
USE modd_snow_par, ONLY : xrhosmin_es, xsnowdmin, xansmax, xaglamax, xsnowcritd,   &
                          xdzmin_top, xdzmin_top_bis, xdzmin_bot, xsplit_coef,     &
                          xagreg_coef_1, xagreg_coef_2, xdz1, xdz2, xdz3, xdz3_bis,&
                          xdz4, xdz5, xdz_base, xdz_internal, xscale_cm,           &
                          xdzmax_internal, xdzmin_top_extrem, xsnowfall_threshold, &
                          xratio_newlayer, xdepth_threshold1, xdepth_threshold2,   &
                          xdepth_surface, xdiff_1, xdiff_max, xscale_diff,         &
                          xsnowfall_a_sn, xsnowfall_b_sn, xsnowfall_c_sn
!
USE mode_snow3l
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(date_time), INTENT(IN)          :: tptime      ! current date and time
LOGICAL, INTENT(IN)                  :: oglacier    ! True = Over permanent snow and ice, 
!                                                     initialise WGI=WSAT,
!                                                     Hsnow>=10m and allow 0.8<SNOALB<0.85
                                                    ! False = No specific treatment
!
REAL, INTENT(IN)                     :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)       :: psr, pta, pvmod, ppermsnowfrac
!
REAL, DIMENSION(:),INTENT(IN)       :: pz0eff,puref
!
REAL, DIMENSION(:), INTENT(INOUT)   :: psnow, psnowalb
! 
REAL, DIMENSION(:,:), INTENT(IN)     :: psnowrho, psnowdz, psnowheat
!
REAL, DIMENSION(:), INTENT(OUT)      :: psnowhmass
!
REAL, DIMENSION(:,:), INTENT(IN)     :: psnowgran1, psnowgran2
!
LOGICAL, DIMENSION(:), INTENT(INOUT) :: gsnowfall   
!
! Fresh snow characteristics
REAL, DIMENSION(:), INTENT(OUT)      :: psnowrhof, psnowdzf 
REAL, DIMENSION(:), INTENT(OUT)      :: psnowgran1f, psnowgran2f, psnowhistf 
REAL, DIMENSION(:), INTENT(OUT)      :: psnowagef  
! New vertical grid
REAL, DIMENSION(:,:), INTENT(OUT)    :: psnowdzn
!
LOGICAL, DIMENSION(:), INTENT(OUT)   :: omodif_grid
! 
INTEGER, DIMENSION(:), INTENT(INOUT) :: knlvls_use

LOGICAL,INTENT(IN) :: osnowdrift ! if snowdrift then grain types are not modified by wind
 CHARACTER(3), INTENT(IN)              :: hsnowmetamo ! metamorphism scheme
!*      0.2    declarations of local variables
!
!
LOGICAL, DIMENSION(SIZE(PTA))       :: gagreg_surf 
!
REAL, DIMENSION(SIZE(PTA))          :: zsnowfall, zsnowtemp, zscap, zansmax
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zdzopt 
!
REAL :: zz0eff
!
REAL :: zage_now
REAL :: zsnow_upper, zsnow_upper2 ! snow depth treatednormally (<= XDEPTH_SURFACE)
REAL :: zcoef_depth !coefficient for repartition of deep snow above 3 meters
REAL :: zthickness_intermediate, zthickness2
REAL :: zpenalty, zdiftype_inf, zdiftype_sup, zcritsize, zcritsize_inf, zcritsize_sup
REAL :: zsnow2l, zcoef                
!
INTEGER :: inb_deep_layer, inb_upper_layer !separation between deep and upper layers
                                           ! if snow depth below XDEPTH_SURFACE then INB_DEEP_LAYER=0
INTEGER :: inb_min_layers    ! why this test ?
INTEGER :: inb_intermediate  ! number of intermediate layers (constant optimal gridding)
INTEGER :: iend_intermediate ! layer indice for bottom of intermediate layers
INTEGER :: jstdeep, jstend    
INTEGER :: jst_1, jj_a_agreg_sup, jj_a_agreg_inf, jj_a_dedoub
INTEGER :: inlvls, inlvlsmin, inlvlsmax, jj, jst
!
! Coefficient to adjust wind speed at the height used in the parameterization
! for:
!         - density of new snow
!         - sphericity and dendricity of new snow
! Default values : 10 m for new snow (Pahaut, 1976) and 5 m for characteristics
! of snow grains (Guyomarc'h et Merindol, 1998)
REAL, PARAMETER                    :: pphref_wind_rho   = 10.
REAL, PARAMETER                    :: pphref_wind_grain = 5.
REAL, PARAMETER                    :: pphref_wind_min = min(pphref_wind_rho,pphref_wind_grain)*0.5
REAL, DIMENSION(SIZE(PTA))         :: zwind_rho
REAL, DIMENSION(SIZE(PTA))         :: zwind_grain
!
REAL(KIND=JPRB) :: zhook_handle
!
!*      1.0   Initialization and snowage calculation for the present date 
!
IF (lhook) CALL dr_hook('SNOWNLFALL_UPGRID',0,zhook_handle)
!
inlvls    = SIZE (psnowrho(:,:),2)
inlvlsmax = SIZE (psnowrho(:,:),2)
inlvlsmin = 3 
!
zsnowtemp(:) = xtt
zsnowfall(:) = 0.0 !Matthieu Lafaysse 21/09/2012
!
gsnowfall(:) =.false.
gagreg_surf(:) =.false.
!
psnowhmass(:) = 0.0 
psnowrhof(:) = 0.0 
psnowdzf(:) = 0.0 
psnowgran1f(:) = 0.0 
psnowgran2f(:) = 0.0      
psnowhistf(:) = 0.0
psnowdzn(:,:) = psnowdz(:,:)
!
omodif_grid(:) = .false.
!
!************************************************************************************
!*      1.1   Calculation of the optimal vertical grid size ZDZOPT
!             as a function of maximum number of layers and of current 
!             snow depth (modified 05/06/2012 by Matthieu Lafaysse)
!
! KNLVLS_USE(JJ) > INB_MIN_LAYERS =>
! KNLVLS_USE(JJ) > 2 + INLVLSMAX/3 =>
! ( KNLVLS_USE(JJ) + INLVLSMAX ) / 6 > (2 + INLVLSMAX/3 + INLVLSMAX) / 6 =>
! INB_DEEP_LAYER > (2 + 4*INLVLSMAX/3 ) / 6 >= 1
inb_min_layers = 2 + inlvlsmax/3
!
DO jj = 1,SIZE(psnow(:))
  !
  IF ( psnow(jj)>xdepth_threshold2 .AND. knlvls_use(jj)>inb_min_layers ) THEN
    ! for very thick snowpack with enough snow layers
    ! special treatment
    ! we put the highest thickness in the lowest layers
    ! about 1/3 of layers for all snow except XDEPTH_SURFACE=3 first meters
    !
    !number of "deep layers"
    inb_deep_layer  = ( knlvls_use(jj) + inlvlsmax ) / 6
    !
    !number of "upper layers"
    inb_upper_layer = knlvls_use(jj) - inb_deep_layer
    !
    !thickness of "upper layers"
    zsnow_upper = xdepth_surface
    !
    !Arithmetic serie : 1+2+3+...+INB_DEEP_LAYER=INB_DEEP_LAYER*(INB_DEEP_LAYER+1)/2
    zcoef_depth = ( psnow(jj) - xdepth_surface ) * 2. / ( (inb_deep_layer+1) * inb_deep_layer )
    !
    ! deep layers optimal thickness :
    ! increasing thickness with depth
    DO jstdeep = 1,inb_deep_layer
      jst = inb_upper_layer + jstdeep
      zdzopt(jj,jst) = zcoef_depth * jstdeep
      !This sum is equal to PSNOW(JJ)-XDEPTH_SURFACE
    ENDDO
    !
  ELSE
    !        
    inb_upper_layer = knlvls_use(jj)
    !
    zsnow_upper = psnow(jj)
    !
  END IF
  !
  !on force le ZDZOPT des 3 premières couches à ZSNOW_UPPER/3 maximum, chacune. 
  ! => si on n'a qu'une couche, ZDZOPT(1) = ZSNOW_UPPER/3
  ! quel que soit INB_UPPER_LAYER
  !
  zsnow_upper2 = zsnow_upper / max( inlvlsmin, inb_upper_layer )
  !
  zdzopt(jj,1) = min( xdz1, zsnow_upper2 )
  IF ( knlvls_use(jj)>=2 ) zdzopt(jj,2) = min( xdz2, zsnow_upper2 )
  IF ( knlvls_use(jj)>=3 ) zdzopt(jj,3) = min( xdz3, zsnow_upper2 )
  !
  IF ( inb_upper_layer>0 ) THEN
    !
    zsnow_upper2 = zsnow_upper / inb_upper_layer
    !
    ! dans ce cas, à partir de la 3ème couche, on prend la fraction du nombre de
    ! couches supérieures total, pour les couches jusqu'à 5
    !
    !ML : replace > by >= on 12-12-20 because the last layer was not initialised in case of thick snowpacks
    IF ( inb_upper_layer>=3 ) zdzopt(jj,3) = min( xdz3_bis, zsnow_upper2 )
    IF ( inb_upper_layer>=4 ) zdzopt(jj,4) = min( xdz4    , zsnow_upper2 )
    IF ( inb_upper_layer>=5 ) zdzopt(jj,5) = min( xdz5    , zsnow_upper2 )
    !
    IF ( inb_upper_layer==knlvls_use(jj) ) THEN
      ! si on n'a pas de couches profondes
      !
      ! dans ce cas, on reprend ZSNOW_UPPER/3 maximum pour la dernière couche
      !
      ! last layer of' upper layers' : normal case : thin layer
      zdzopt(jj,inb_upper_layer) = min( xdz_base, zsnow_upper/max(inlvlsmin,inb_upper_layer) )
      !
      ! ZTHICKNESS_INTERMEDIATE contient ce qu'il reste d'épaisseur disponible
      ! dans les couches supérieures
      !remaining snow for remaining layers
      zthickness_intermediate = zsnow_upper - sum(zdzopt(jj,1:5)) - zdzopt(jj,inb_upper_layer)

      IF ( zsnow_upper<=xdepth_threshold1 .OR. inb_upper_layer<8 ) THEN             
        inb_intermediate  = inb_upper_layer - 6
        iend_intermediate = inb_upper_layer - 1
      ELSE
        ! si INB_UPPER_LAYER>=8, les avant et avant-dernière couches ne sont pas
        ! considérées commes intermédiaires
        inb_intermediate  = inb_upper_layer - 8
        iend_intermediate = inb_upper_layer - 3
        ! dans ce cas, on garde un peu d'épaisseur pour les deux couches restantes
        IF ( inb_intermediate>0 ) THEN
          zthickness_intermediate = zthickness_intermediate * inb_intermediate / float(inb_intermediate+1)
        END IF
      END IF
      !
    ELSE
      ! si on a des couches profondes, les couches intermédiaires sont celles
      ! qui restent quand on a enlevé les 5 premières des couches supérieures
      ! 
      ! case with very thick snowpacks :
      ! the last layer of upper layers is not an exception
      zthickness_intermediate = zsnow_upper - sum(zdzopt(jj,1:5))
      inb_intermediate  = inb_upper_layer - 5
      iend_intermediate = inb_upper_layer 
      !  
    END IF
    !
    ! For thick snowpack : add maximum value of optimal thickness to avoid too 
    ! large differencies between layers
    IF ( inb_intermediate>0 ) THEN
      !
      zthickness2 = max( xdz_internal, zthickness_intermediate/inb_intermediate )
      !
      jstend = min( iend_intermediate,10 )
      DO jst = 6,jstend
        zdzopt(jj,jst) = min( xdzmax_internal(jst-5), zthickness2 ) 
      END DO
      !
      IF ( iend_intermediate>10 ) THEN
        DO jst = 11,iend_intermediate
          zdzopt(jj,jst) = zthickness2
        END DO
      END IF
      !
    END IF
    !
    IF ( zsnow_upper>=xdepth_threshold1 .AND. inb_upper_layer>=8 ) THEN
      !Linear interpolation of optimal thickness between layers N-3 and N :
      zdzopt(jj,inb_upper_layer-2) = 0.34*zdzopt(jj,inb_upper_layer) + &
                                     0.66*zdzopt(jj,inb_upper_layer-3)
      zdzopt(jj,inb_upper_layer-1) = 0.66*zdzopt(jj,inb_upper_layer) + &
                                     0.34*zdzopt(jj,inb_upper_layer-3)
    ENDIF
    !
  END IF
  !
END DO
!
!************************************************************************************
!This was the initial code for optimal layers until may 2012
!
! ! ! ! ! 
! ! ! ! ! !*      1.1   Calculation of the optimal vertical grid size 
! ! ! ! ! !             as a function of maximum number of layers and of current 
! ! ! ! ! !             snow depth
! ! ! ! ! !
! ! ! ! ! DO JJ=1, SIZE(PSNOW(:))
! ! ! ! !    ZDZOPT(JJ,1) = MIN(XDZ1,PSNOW(JJ)/MAX(INLVLSMIN,KNLVLS_USE(JJ))) 
! ! ! ! !    ZDZOPT(JJ,2) = MIN(XDZ2,PSNOW(JJ)/MAX(INLVLSMIN,KNLVLS_USE(JJ))) 
! ! ! ! !    ZDZOPT(JJ,3) = MIN(XDZ3,PSNOW(JJ)/MAX(INLVLSMIN,KNLVLS_USE(JJ))) 
! ! ! ! !    IF (KNLVLS_USE(JJ)>3)  ZDZOPT(JJ,3) = MIN(XDZ3_BIS,PSNOW(JJ)/KNLVLS_USE(JJ))
! ! ! ! !    IF (KNLVLS_USE(JJ)>4)  ZDZOPT(JJ,4) = MIN(XDZ4,PSNOW(JJ)/KNLVLS_USE(JJ))
! ! ! ! !    IF (KNLVLS_USE(JJ)>5)  ZDZOPT(JJ,5) = MIN(XDZ5,PSNOW(JJ)/KNLVLS_USE(JJ))
! ! ! ! !    IF (KNLVLS_USE(JJ)>0) ZDZOPT(JJ,KNLVLS_USE(JJ))=   &
! ! ! ! !               MIN(XDZ_BASE,PSNOW(JJ)/MAX(INLVLSMIN,KNLVLS_USE(JJ)))  
! ! ! ! !    DO JST=6,KNLVLS_USE(JJ)-1,1
! ! ! ! !        ZDZOPT(JJ,JST) = MAX(XDZ_INTERNAL,(PSNOW(JJ) - SUM(ZDZOPT(JJ,1:5))-    &
! ! ! ! !                          ZDZOPT(JJ,KNLVLS_USE(JJ))) /(KNLVLS_USE(JJ)-6)) 
! ! ! ! !    END DO
! ! ! ! ! END DO
! ! ! ! ! !
! ! ! ! ! 
!
!************************************************************************************
!
!*      2.0   Fresh snow characteristics
!
!
!
! Heat content of newly fallen snow (J/m2):
! NOTE for now we assume the snowfall has
! the temperature of the snow surface upon reaching the snow.
! This is done as opposed to using the air temperature since
! this flux is quite small and has little to no impact
! on the time scales of interest. If we use the above assumption
! then, then the snowfall advective heat flux is zero.
!!
DO jj = 1,SIZE(psnow(:))
  !
  IF ( psr(jj)>0.0 ) THEN  
    !    
    ! newly fallen snow characteristics:
    IF ( knlvls_use(jj)>0 ) THEN !Case of new snowfall on a previously snow-free surface 
      zscap(jj) = xci*psnowrho(jj,1)
      zsnowtemp(jj) = xtt + ( psnowheat(jj,1) + xlmtt*psnowrho(jj,1)*psnowdz(jj,1) ) / &
                            ( zscap(jj) * max( xsnowdmin/inlvls, psnowdz(jj,1) ) ) 
    ELSE  ! case with bare ground
      zsnowtemp(jj) = pta(jj)
    ENDIF
    zsnowtemp(jj) = min( xtt, zsnowtemp(jj) )
    !
    !
    ! Wind speeds at reference heights for new snow density and charactristics of
    ! grains of new snow 
    ! Computed from PVMOD at PUREF (m) assuming a log profile in the SBL 
    ! and a roughness length equal to PZ0EFF
    !
    zz0eff=min(pz0eff(jj),puref(jj)*0.5,pphref_wind_min)
    
    zwind_rho(jj)   = pvmod(jj)*log(pphref_wind_rho/zz0eff)/          &
                               log(puref(jj)/zz0eff)
    zwind_grain(jj) = pvmod(jj)*log(pphref_wind_grain/zz0eff)/        &
                               log(puref(jj)/zz0eff)    
    
    psnowhmass(jj) = psr(jj) * ( xci * ( zsnowtemp(jj)-xtt ) - xlmtt ) * ptstep
    !
    psnowrhof(jj) = max( xrhosmin_es, xsnowfall_a_sn + &
                                       xsnowfall_b_sn * ( pta(jj)-xtt ) + &
                                       xsnowfall_c_sn * min( pvmod(jj), sqrt(zwind_rho(jj) ) ) ) 
    zsnowfall(jj) = psr(jj) * ptstep / psnowrhof(jj)    ! snowfall thickness (m)
    psnow(jj) = psnow(jj) + zsnowfall(jj)
    psnowdzf(jj) = zsnowfall(jj)
    !
    IF ( hsnowmetamo=='B92' ) THEN
      !
      IF ( osnowdrift ) THEN
        psnowgran1f(jj) = -xgran
        psnowgran2f(jj) = xnsph3
      ELSE
        psnowgran1f(jj) = max( min( xnden1*zwind_grain(jj)-xnden2, xnden3 ), -xgran )
        psnowgran2f(jj) = min( max( xnsph1*zwind_grain(jj)+xnsph2, xnsph3 ), xnsph4 )     
      END IF
      !
    ELSE
      !
      IF ( osnowdrift ) THEN
        psnowgran1f(jj) = xvdiam6
        psnowgran2f(jj) = xnsph3/xgran
      ELSE
        psnowgran2f(jj) = min( max( xnsph1*zwind_grain(jj)+xnsph2, xnsph3 ), xnsph4 ) / xgran
        zcoef = max( min( xnden1*zwind_grain(jj)-xnden2, xnden3 ), -xgran ) / ( -xgran )
        psnowgran1f(jj) = xvdiam6 * &
                        ( zcoef + ( 1.- zcoef ) * &
                                  ( 3.*psnowgran2f(jj) + 4.*(1.-psnowgran2f(jj)) ) )
      END IF      
      !
    ENDIF
    !
    psnowhistf(jj) = 0.0
    psnowagef(jj) = 0.0
    gsnowfall(jj) = .true. 
    omodif_grid(jj) = .true.
    !
  ENDIF
  !
ENDDO
!                   
!  intialize the albedo:
!  penser a changer 0.000001 par XUEPSI
IF(oglacier)THEN
  zansmax(:) = xaglamax * ppermsnowfrac(:) + xansmax * (1.0-ppermsnowfrac(:))
ELSE
  zansmax(:) = xansmax
ENDIF
!
WHERE( gsnowfall(:) .AND. abs(psnow(:)-zsnowfall(:))< 0.000001 )
  psnowalb(:) = zansmax(:)
END WHERE
!
! Computation of the new grid size
! It starts with successive exclusive cases
! Each case is described inside the corresponding condition
! 
! cases with fresh snow
!
DO jj=1,SIZE(psnow(:)) ! grid point loop
  !
  IF( .NOT.gsnowfall(jj) .AND. psnow(jj)>=xsnowcritd .AND. knlvls_use(jj)>=inlvlsmin ) THEN 
    !
    ! no fresh snow + deep enough snowpack + enough snow layers ==> no change
    !
  ELSEIF( psnow(jj)<xsnowcritd .OR. knlvls_use(jj)<inlvlsmin .OR. psnow(jj)==zsnowfall(jj) ) THEN 
    !
    ! too shallow snowpack or too few layers or only fresh snow
    ! ==> uniform grid and identical snow layers / number depends on snow depth 
    omodif_grid(jj) = .true.
    knlvls_use(jj) = max( inlvlsmin, min( inlvlsmax, int(psnow(jj)*xscale_cm) ) )
    psnowdzn(jj,1:knlvls_use(jj)) = psnow(jj) / knlvls_use(jj)
    !
  ELSE
    !
    ! fresh snow over snow covered ground + enough snow layers 
    omodif_grid(jj) = .true.
    zdiftype_sup = snow3ldiftyp( psnowgran1(jj,1),psnowgran1f(jj), &
                                 psnowgran2(jj,1),psnowgran2f(jj),hsnowmetamo ) 
    !
    IF ( ( zdiftype_sup<xdiff_1        .AND. psnowdz(jj,1)<   zdzopt(jj,1) ) .OR. &
         ( psr(jj)<xsnowfall_threshold .AND. psnowdz(jj,1)<2.*zdzopt(jj,1) ) .OR. &
                                             psnowdz(jj,1)<xdzmin_top_extrem ) THEN
      !
      ! Fresh snow is similar to a shallow surface layer (< ZDZOPT)
      ! or snowfall is very low and the surface layer not too deep (< 2*ZDZOPT) [NEW CONDITION 11/2012]
      ! or the surface layer is extremely thin (< XDZMIN_TOP_EXTREM) [NEW CONDITION 11/2012]
      ! The two new conditions are necessary for forcings with very low precipitation
      ! (e.g. ERA interim reanalyses, or climate models)
      ! ==> fresh snow is agregated to the surface layer
      psnowdzn(jj,1) = psnowdz(jj,1) + psnowdzf(jj)
      DO jst = knlvls_use(jj),2,-1
        psnowdzn(jj,jst) = psnowdz(jj,jst)
      ENDDO
      !
    ELSEIF ( knlvls_use(jj)<inlvlsmax ) THEN
      !
      ! fresh snow is too different from the surface or the surface is too deep
      ! and there is room for extra layers ==> we create a new layer
      knlvls_use(jj)=knlvls_use(jj)+1
      !
      IF ( psnowdzf(jj)>xratio_newlayer*psnowdz(jj,2) ) THEN
        !       
        ! Snowfall is sufficient to create a new layer not lower than 1/10 of the second layer
        psnowdzn(jj,1) = psnowdzf(jj)
        DO jst = knlvls_use(jj),2,-1
          psnowdzn(jj, jst) = psnowdz(jj,jst-1)
        ENDDO
        !
      ELSE 
        ! The ratio would be lower than 1/10 : [NEW : 11/2012]
        ! aggregate a part of the old layer with fresh snow to limit the ratio to 1/10.
        zsnow2l = psnowdzf(jj) + psnowdz(jj,1)
        psnowdzn(jj,1) = xratio_newlayer      * zsnow2l
        psnowdzn(jj,2) = (1.-xratio_newlayer) * zsnow2l
        DO jst = knlvls_use(jj),3,-1
          psnowdzn(jj,jst) = psnowdz(jj,jst-1)
        ENDDO
        !
      ENDIF     
      !
    ELSE
      !
      ! fresh snow is too different from the surface or the surface is too deep
      ! and there is no room for extra layers
      ! ==> we agregate internal most similar snowlayers and create a new surface layer
      jj_a_agreg_sup = 1
      jj_a_agreg_inf = 2
      !
      DO jst = 1,knlvls_use(jj)
        !
        IF ( jst>1 ) THEN
          !
          zcritsize_sup = xscale_diff * ( psnowdz(jj,jst)  /zdzopt(jj,jst)    + &
                                          psnowdz(jj,jst-1)/zdzopt(jj,jst-1) )
          zdiftype_sup  = snow3ldiftyp( psnowgran1(jj,jst-1),psnowgran1(jj,jst), &
                                        psnowgran2(jj,jst-1),psnowgran2(jj,jst), &
                                        hsnowmetamo )
          !
          IF ( zdiftype_sup+zcritsize_sup<zpenalty ) THEN
            zpenalty = zdiftype_sup + zcritsize_sup
            jj_a_agreg_sup = jst - 1
            jj_a_agreg_inf = jst
          ENDIF
          !
        ENDIF
        !
        IF ( jst<knlvls_use(jj) ) THEN
          !
          zcritsize_inf = xscale_diff * ( psnowdz(jj,jst)  /zdzopt(jj,jst)    + &
                                          psnowdz(jj,jst+1)/zdzopt(jj,jst+1) )
          !
          IF ( jst==1 ) THEN
            zdiftype_inf  = snow3ldiftyp( psnowgran1(jj,1),psnowgran1f(jj), &
                                          psnowgran2(jj,1),psnowgran2f(jj), &
                                          hsnowmetamo)
            !
            zpenalty = zdiftype_inf + zcritsize_inf
          ELSE
            zdiftype_inf  = snow3ldiftyp( psnowgran1(jj,jst+1),psnowgran1(jj,jst), &
                                          psnowgran2(jj,jst+1),psnowgran2(jj,jst), &
                                          hsnowmetamo)
            !
            IF ( zdiftype_inf+zcritsize_inf<zpenalty ) THEN
              zpenalty = zdiftype_inf + zcritsize_inf
              jj_a_agreg_sup = jst
              jj_a_agreg_inf = jst + 1
            ENDIF
          ENDIF
          !
        ENDIF
        !
      ENDDO
      !
      ! agregation of the similar layers and shift of upper layers
      psnowdzn(jj,jj_a_agreg_inf) = psnowdz(jj,jj_a_agreg_inf) + psnowdz(jj,jj_a_agreg_sup) 
      DO jst = jj_a_agreg_sup,2,-1
        psnowdzn(jj,jst) = psnowdz(jj,jst-1)
      ENDDO
      psnowdzn(jj,1) = psnowdzf(jj)
      !
      ! Limit the ratio between the new layer and the one beneath (ratio 1/10)
      ! [NEW : 11/2012]
      IF( psnowdzn(jj,1)<xratio_newlayer*psnowdzn(jj,2) ) THEN 
        zsnow2l = psnowdzn(jj,1) + psnowdzn(jj,2)
        psnowdzn(jj,1) = xratio_newlayer      * zsnow2l
        psnowdzn(jj,2) = (1.-xratio_newlayer) * zsnow2l
      ENDIF    
      !
    ENDIF
    !
  ENDIF ! end of the case with fresh snow
  !
ENDDO ! end loop grid points
!
! cases with no fresh snow and no previous grid resize
!
IF ( inlvlsmin==inlvlsmax ) THEN ! specific case with INLSVSMIN = INLVLSMAX  (INLVLS)
  !
  ! check if surface layer depth is too small
  ! in such a case looks for an other layer to be split
  DO jj = 1,SIZE(psnow(:)) ! loop grid points
    !
    IF ( .NOT.omodif_grid(jj) .AND. psnowdz(jj,1)<xdzmin_top ) THEN
      omodif_grid(jj) = .true.
      CALL get_snowdzn_deb(inlvls,psnowdz(jj,:),zdzopt(jj,:),psnowdzn(jj,:))
      gagreg_surf(jj) = .true.
    ENDIF
    ! check if bottom layer depth is too small
    ! in such a case agregation with upper layer and
    ! looks for an other layer to be splitted
    IF( .NOT.omodif_grid(jj) .AND. psnowdz(jj,inlvls)<xdzmin_top ) THEN
      omodif_grid(jj) = .true.
      CALL get_snowdzn_end(inlvls,psnowdz(jj,:),zdzopt(jj,:),psnowdzn(jj,:))
    ENDIF
    !
  ENDDO ! end grid points loop
  !
ENDIF  ! end specific case INLSVSMIN = INLVLSMAX
!
! case without new snowfall and INVLS > INLVLSMIN
!
DO jj=1,SIZE(psnow(:))
  !
  ! check if surface layer depth is too small
  ! in such a case agregation with layer beneath
  ! in case of reaching INLVLSMIN, looks for an other layer to be splitted
  IF( .NOT.gsnowfall(jj) .AND. psnow(jj)>xsnowcritd .AND. &
      .NOT.omodif_grid(jj) .AND. psnowdz(jj,1)<xdzmin_top_bis ) THEN ! case shallow surface layer 
    !
    omodif_grid(jj) = .true.
    !
    IF( knlvls_use(jj)>inlvlsmin ) THEN ! case minimum not reached
      knlvls_use(jj) = knlvls_use(jj) - 1
      psnowdzn(jj,1) = psnowdz(jj,1) + psnowdz(jj,2)
      DO jst = 2,knlvls_use(jj)
        psnowdzn(jj,jst) = psnowdz(jj,jst+1)
      ENDDO
    ELSE ! case minimum reached
      CALL get_snowdzn_deb(knlvls_use(jj),psnowdz(jj,:),zdzopt(jj,:),psnowdzn(jj,:))
    ENDIF ! end case minimum reached end case shallow surface layer
    !
    gagreg_surf(jj) = .true.
    !
  ENDIF
  !
  ! check if bottom layer depth is too small
  ! in such a case agregation with above layer 
  ! in case of reaching INLVLSMIN, looks for an other layer to be splitted
  ! case shallow bottom layer
  IF( .NOT.gsnowfall(jj) .AND. psnow(jj)> xsnowcritd                   .AND. &
      .NOT.omodif_grid(jj) .AND. psnowdz(jj,knlvls_use(jj))<xdzmin_top .AND. &
      .NOT.gagreg_surf(jj) ) THEN
    !
    omodif_grid(jj) = .true.
    !
    IF ( knlvls_use(jj)>inlvlsmin ) THEN ! case minimum not reached
      knlvls_use(jj) = knlvls_use(jj) - 1
      psnowdzn(jj,knlvls_use(jj)) = psnowdz(jj,knlvls_use(jj)) + psnowdz(jj,knlvls_use(jj)+1)
    ELSE  ! case minimum reached
      CALL get_snowdzn_end(knlvls_use(jj),psnowdz(jj,:),zdzopt(jj,:),psnowdzn(jj,:))
    ENDIF ! end case minimum reached end case shallow surface layer
    !
  ENDIF
  !
ENDDO ! end grid points loop  
!
! case whithout new snow fall and without a previous grid resize 
! looks for a shallow layer to be splitted according to its depth and to
! the optimal grid size
DO jj = 1,SIZE(psnow(:))
  !
  IF ( .NOT.omodif_grid(jj) .AND. inlvls_use(jj)<inlvls-3 )THEN
    !  
    DO jst = 1,inlvls-4
      !
      IF ( jst<=knlvls_use(jj) .AND. .NOT.omodif_grid(jj) ) THEN
        !
        IF( psnowdz(jj,jst) > &
            ( xsplit_coef - float( inlvls-knlvls_use(jj) )/max( 1., float( inlvls-inlvlsmin ) ) ) &
              * zdzopt(jj,jst) ) THEN
          !
          DO jst_1 = knlvls_use(jj)+1,jst+2,-1
            psnowdzn(jj,jst_1) = psnowdz(jj,jst_1-1)
            zdzopt(jj,jst_1) = zdzopt(jj,jst_1-1)
          ENDDO
          !
          ! generale case : old layer divided in two equal layers
          IF ( jst/=1 .OR. psnowdz(jj,jst)<3.*zdzopt(jj,1) ) THEN
            psnowdzn(jj,jst+1) = 0.5*psnowdz(jj,jst)
            psnowdzn(jj,jst)   = psnowdzn(jj,jst+1)
          ELSE
            ! if thick surface layer : force the surface layer to this value to avoid successive resizing
            ! [NEW : 11/2012]
            psnowdzn(jj,1) = 1.5 * zdzopt(jj,1)
            psnowdzn(jj,2) = psnowdz(jj,jst) - psnowdzn(jj,1) 
          ENDIF
          !
          knlvls_use(jj) = knlvls_use(jj) + 1
          omodif_grid(jj) = .true. 
          !
        ENDIF
        !
      ENDIF
      !
    ENDDO
    !
  ENDIF
  !
ENDDO
!
! case whithout new snow fall and without a previous grid resize 
! looks for a deep layer to be agregated to the layer beneath if similar
! according to its depth and to the optimal grid size
!
!NB : allow these changes for 5 layers and more [NEW] (before : 6 layers)
!
DO jj = 1,SIZE(psnow(:))
  !
  IF ( .NOT.omodif_grid(jj) ) THEN
    !
    DO jst = 2,inlvls
      !
      IF ( jst<=knlvls_use(jj)-1 .AND. knlvls_use(jj)>inlvlsmin+1 .AND. .NOT.omodif_grid(jj) ) THEN
        !
        zdiftype_inf = snow3ldiftyp( psnowgran1(jj,jst+1),psnowgran1(jj, jst), &
                                     psnowgran2(jj,jst+1),psnowgran2(jj, jst), &
                                     hsnowmetamo) 
        zdiftype_inf = max( xdiff_1, min( xdiff_max, zdiftype_inf ) )
        !
        IF( psnowdz(jj,jst) < zdzopt(jj,jst) * xagreg_coef_1 / zdiftype_inf .AND. &
            psnowdz(jj,jst) + psnowdz(jj,jst+1) < &
                xagreg_coef_2 * max( zdzopt(jj,jst),zdzopt(jj,jst+1) ) ) THEN
          !
          psnowdzn(jj,jst) = psnowdz(jj,jst) + psnowdz(jj,jst+1)
          zdzopt(jj,jst) = zdzopt(jj,jst+1)
          DO jst_1 = jst+1,knlvls_use(jj)-1
            psnowdzn(jj,jst_1) = psnowdz(jj,jst_1+1)
            zdzopt(jj,jst_1) = zdzopt(jj,jst_1+1)
          ENDDO
          knlvls_use(jj) = knlvls_use(jj)-1
          omodif_grid(jj)=.true.
          !
        ENDIF
        !
      ENDIF
      !
    ENDDO
    !
  ENDIF
  ! 
ENDDO
!
! [NEW : 11/2012]
! In case of very low snow fall checks if a new internal snow layer is too shallow
! even if a the grid has already been resized in this time step
! starts from bottom to INLVS_USE-3 until old and new grid differ
DO jj = 1,SIZE(psnow(:))
  !
  IF ( .NOT.gsnowfall(jj) .OR. knlvls_use(jj)<inlvlsmin+3 ) cycle ! go to next point 
  !
  IF( abs( psnowdzn(jj,knlvls_use(jj)) - psnowdz(jj,knlvls_use(jj)) ) > xuepsi ) cycle ! go to next point 
  !
  ! bottom layer
  IF( psnowdzn(jj,knlvls_use(jj))<xdzmin_top ) THEN ! case shallow bottom layer
    !
    knlvls_use(jj) = knlvls_use(jj)-1
    psnowdzn(jj,knlvls_use(jj)) = psnowdzn(jj,knlvls_use(jj)) + psnowdzn(jj,knlvls_use(jj)+1)
    psnowdzn(jj,knlvls_use(jj)+1) = 0.
    !
  ELSE
    !
    ! internal layer
    DO jst = knlvls_use(jj)-1,4,-1
      !
      IF ( abs( psnowdzn(jj,jst) - psnowdz(jj,jst) ) > xuepsi ) EXIT ! old/new grid differ ==> go to next grid point
      !
      IF ( psnowdzn(jj,jst)> 0.001 ) cycle
      !
      ! If an internal layer is too shallow, it is merged with the upper layer              
      psnowdzn(jj,jst-1) = psnowdzn(jj,jst) + psnowdzn(jj,jst-1)
      knlvls_use(jj)   = knlvls_use(jj) - 1
      !
      ! shifts the lower layers       
      DO jst_1 = jst,inlvls_use(jj)
        psnowdzn(jj,jst_1) = psnowdz(jj,jst_1+1)
        zdzopt(jj,jst_1) = zdzopt(jj,jst_1+1)
      ENDDO 
      psnowdzn(jj,inlvls_use(jj)+1) = 0.
      !      
      EXIT ! goto to next grid point
      !
    ENDDO ! end loop internal layers     
    !    
  ENDIF
  !
ENDDO ! end grid loops for checking shallow layers   
!
!final check of the consistensy of the new grid size
!
DO jj = 1,SIZE(psnow(:))
  !
  IF ( abs( sum( psnowdzn(jj,1:knlvls_use(jj)) ) - psnow(jj) ) > xuepsi ) THEN
    !
    WRITE(*,*) 'error in grid resizing', jj, knlvls_use(jj), sum( psnowdzn(jj,1:knlvls_use(jj)) ),  &
                                         psnow(jj), sum( psnowdzn(jj,1:inlvls_use(jj)) )-psnow(jj), &
                                         zsnowfall(jj)
    WRITE( *,*) 'JJ , PSNOWDZ(JJ):',jj ,  psnowdz(jj,:)
    WRITE( *,*) 'JJ , PSNOWDZN(JJ):',jj , psnowdzn(jj,:)
    !
    CALL abor1_sfx("SNOWCRO: error in grid resizing")
    !
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('SNOWNLFALL_UPGRID',1,zhook_handle)
!
END SUBROUTINE snownlfall_upgrid

!###############################################################################
SUBROUTINE get_snowdzn_deb(KNLVLS,PSNOWDZ,PDZOPT,PSNOWDZN)
!
USE modd_snow_par, ONLY : xdzmin_top, xdzmin_bot
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: knlvls
REAL, DIMENSION(:), INTENT(IN)  :: psnowdz, pdzopt
REAL, DIMENSION(:), INTENT(OUT) :: psnowdzn
!
REAL :: zpenalty, zcritsize
INTEGER :: jj_a_dedoub, jst
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_SNOWDZN_DEB',0,zhook_handle)
!
zpenalty = psnowdz(2) / pdzopt(2)
IF( psnowdz(2)<xdzmin_top ) zpenalty = 0.
jj_a_dedoub = 2
!
DO jst = 3,knlvls
  zcritsize = psnowdz(jst) / pdzopt(jst)
  IF ( jst==knlvls .AND. psnowdz(jst)<xdzmin_bot ) zcritsize = 0.
  IF ( zcritsize>zpenalty ) THEN
    zpenalty    = zcritsize
    jj_a_dedoub = jst
  ENDIF
ENDDO
!
IF ( jj_a_dedoub==2 ) THEN ! case splitted layer == 2
  psnowdzn(1) = 0.5 * ( psnowdz(1) + psnowdz(2) )
  psnowdzn(2) = psnowdzn(1)
ELSE ! case splitted layer =/ 2
  psnowdzn(1) = psnowdz(1) + psnowdz(2)
  DO jst = 2,jj_a_dedoub-2
    psnowdzn(jst) = psnowdz(jst+1)
  ENDDO
  psnowdzn(jj_a_dedoub-1) = 0.5 * psnowdz(jj_a_dedoub)  
  psnowdzn(jj_a_dedoub)   = psnowdzn(jj_a_dedoub-1)
ENDIF ! end case splitted layer =/ 2
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_SNOWDZN_DEB',1,zhook_handle)
!
END SUBROUTINE get_snowdzn_deb
!
!###############################################################################
SUBROUTINE get_snowdzn_end(KNLVLS,PSNOWDZ,PDZOPT,PSNOWDZN)
!
USE modd_snow_par, ONLY : xdzmin_top, xdzmin_bot
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: knlvls
REAL, DIMENSION(:), INTENT(IN)  :: psnowdz, pdzopt
REAL, DIMENSION(:), INTENT(OUT) :: psnowdzn
!
REAL :: zpenalty, zcritsize
INTEGER :: jj_a_dedoub, jst
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_SNOWDZN_END',0,zhook_handle)
!
zpenalty = psnowdz(knlvls-2) / pdzopt(knlvls-2)
jj_a_dedoub = knlvls - 2
!
DO jst = max(1,knlvls-3),1,-1
  zcritsize = psnowdz(jst) / pdzopt(jst)
  IF ( jst==1 .AND. psnowdz(jst)<xdzmin_bot ) zcritsize = 0.
  IF ( zcritsize>zpenalty ) THEN
    zpenalty    = zcritsize
    jj_a_dedoub = jst
  ENDIF
ENDDO
!
IF ( jj_a_dedoub==knlvls-1 ) THEN ! case splitted layer == 2
  psnowdzn(knlvls)   = 0.5 * (psnowdz(knlvls-1)+psnowdz(knlvls))
  psnowdzn(knlvls-1) = psnowdzn(knlvls)
ELSE ! case splitted layer =/ 2
  psnowdzn(knlvls) = psnowdz(knlvls-1) + psnowdz(knlvls)
  DO jst = knlvls-1,jj_a_dedoub+2,-1
    psnowdzn(jst) = psnowdz(jst-1)
  ENDDO
  psnowdzn(jj_a_dedoub+1) = 0.5 * psnowdz(jj_a_dedoub)
  psnowdzn(jj_a_dedoub  ) = psnowdzn(jj_a_dedoub+1)
ENDIF ! end case splitted layer =/ 2
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_SNOWDZN_END',1,zhook_handle)
!
END SUBROUTINE get_snowdzn_end
!
!###############################################################################
!################################################################################
!################################################################################
!
SUBROUTINE snownlgridfresh_1d (KJ,PSNOW,PSNOWDZ,PSNOWDZN,                  &
                               psnowrho,psnowheat,psnowgran1,psnowgran2,   &
                               psnowhist,psnowage,gsnowfall,               &
                               psnowrhof, psnowdzf,psnowheatf,psnowgran1f, &
                               psnowgran2f, psnowhistf,psnowagef,          &
                               knlvls_use, hsnowmetamo                     ) 
!
!!    PURPOSE
!!    -------
!     Snow mass,heat and characteristics redistibution in case of
!     grid resizing. Total mass and heat content of the overall snowpack
!     unchanged/conserved within this routine.
!     Grain size and type of mixed layers is deduced from the conservation
!     of the average optical size
!
!!    AUTHOR
!!    ------
!!      E. Brun           * Meteo-France *
!!
!
USE modd_snow_par, ONLY : xd1,xd2,xd3,xx,xvalb5,xvalb6
USE mode_snow3l, ONLY : get_mass_heat
!
IMPLICIT NONE
!
!
!*      0.1    declarations of arguments
!
INTEGER, INTENT(IN)               :: kj
REAL, INTENT(IN)                  :: psnow
!
REAL, DIMENSION(:), INTENT(INOUT) :: psnowheat, psnowrho, psnowdz,     &
                                     psnowdzn, psnowgran1, psnowgran2, &
                                     psnowhist 
REAL, DIMENSION(:), INTENT(INOUT) :: psnowage  
REAL,  INTENT(IN)                 :: psnowrhof, psnowdzf,psnowheatf,   &
                                     psnowgran1f,psnowgran2f, psnowhistf 
REAL, INTENT(IN)                  :: psnowagef                                 
!
INTEGER, INTENT(IN)               :: knlvls_use
!
LOGICAL, INTENT(IN)               :: gsnowfall
!
 CHARACTER(3),INTENT(IN)           :: hsnowmetamo
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1)+1) :: zsnowrhoo,zsnowgran1o,zsnowgran2o, &
                                       zsnowheato,zsnowhisto,zsnowdzo,    &
                                       zsnowztop_old,zsnowzbot_old
REAL,DIMENSION(SIZE(PSNOWRHO,1)+1)  :: zsnowageo  
!
REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: zsnowrhon,zsnowgran1n,zsnowgran2n,   &
                                     zsnowheatn,zsnowhistn,               &
                                     zsnowztop_new,zsnowzbot_new
REAL,DIMENSION(SIZE(PSNOWRHO,1)) ::zsnowagen
!
REAL :: zmastotn, zmastoto, zsnowhean, zsnowheao 
REAL :: zpsnow_old, zpsnow_new
!
INTEGER :: inlvls_old, inlvls_new
INTEGER :: jst 
!
LOGICAL :: gdiam
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWNLGRIDFRESH_1D',0,zhook_handle)
!
! 0. Initialization:
! ------------------
!
! starts by checking the consistency between both vertical grid sizes
inlvls_new = knlvls_use 
inlvls_old = -1
!
zpsnow_new = 0.
zpsnow_old = 0.
!
DO jst = 1,inlvls_new
  zpsnow_new = zpsnow_new + psnowdzn(jst)
ENDDO
!
IF ( abs( zpsnow_new - psnowdzf )<xuepsi ) THEN
  inlvls_old = 0
ELSE
  DO jst = 1,SIZE(psnowrho)
    IF ( psnowdz(jst)>=xuepsi ) THEN
      zpsnow_old = zpsnow_old + psnowdz(jst)
      IF ( abs( zpsnow_new - psnowdzf - zpsnow_old )<xuepsi ) THEN
        inlvls_old = jst
      ENDIF
    ENDIF
  ENDDO
  IF ( inlvls_old==-1 ) THEN
    WRITE(*,*)'pb INLVLS_OLD INLVLS_NEW=',inlvls_new
    WRITE(*,*)'pb INLVLS_OLD',psnowdzf
    WRITE(*,*)'pb INLVLS_OLD',psnowdzn
    WRITE(*,*)'pb INLVLS_OLD',psnowdz
    CALL abor1_sfx('SNOWCRO: Error INLVLS_OLD')
  ENDIF
ENDIF
!
IF ( gsnowfall ) inlvls_old = inlvls_old + 1
!
zpsnow_old = psnow
zpsnow_new = zpsnow_old
!
! initialization of variables describing the initial snowpack + new snowfall
!
IF ( gsnowfall ) THEN
  DO jst = 2,inlvls_old
    zsnowdzo(jst) = psnowdz(jst-1)
    zsnowrhoo(jst) = psnowrho(jst-1)
    zsnowheato(jst) = psnowheat(jst-1)
    zsnowgran1o(jst) = psnowgran1(jst-1)
    zsnowgran2o(jst) = psnowgran2(jst-1)
    zsnowhisto(jst) = psnowhist(jst-1)
    zsnowageo(jst) = psnowage(jst-1)
  ENDDO
  zsnowdzo(1) = psnowdzf
  zsnowrhoo(1) = psnowrhof
  zsnowheato(1) = psnowheatf
  zsnowgran1o(1) = psnowgran1f
  zsnowgran2o(1) = psnowgran2f
  zsnowhisto(1) = psnowhistf
  zsnowageo(1) = psnowagef
ELSE
  DO jst = 1,inlvls_old
    zsnowdzo(jst) = psnowdz(jst)
    zsnowrhoo(jst) = psnowrho(jst)
    zsnowheato(jst) = psnowheat(jst)
    zsnowgran1o(jst) = psnowgran1(jst)
    zsnowgran2o(jst) = psnowgran2(jst)
    zsnowhisto(jst) = psnowhist(jst)
    zsnowageo(jst) = psnowage(jst)
  ENDDO
ENDIF    
!
! 1. Calculate vertical grid limits (m):
! --------------------------------------
!
zsnowztop_old(1) = zpsnow_old
zsnowztop_new(1) = zpsnow_new
!
DO jst = 1,inlvls_old
  IF ( jst>1 ) zsnowztop_old(jst) = zsnowzbot_old(jst-1)
  zsnowzbot_old(jst) = zsnowztop_old(jst) - zsnowdzo(jst)
ENDDO
!
DO jst = 1,inlvls_new
  IF ( jst>1 ) zsnowztop_new(jst) = zsnowzbot_new(jst-1)
  zsnowzbot_new(jst) = zsnowztop_new(jst) - psnowdzn(jst)
ENDDO
!
! Check consistency
IF ( abs(zsnowzbot_old(inlvls_old)) > 0.00001 ) WRITE (*,*) 'Error bottom OLD'
!
zsnowzbot_old(inlvls_old) = 0.
!
! Check consistency
if ( abs(zsnowzbot_new(inlvls_new)) > 0.00001 ) WRITE (*,*) 'Error bottom NEW'
!
zsnowzbot_new(inlvls_new) = 0.
!
! 3. Calculate mass, heat, charcateristics mixing due to vertical grid resizing:
! --------------------------------------------------------------------
!
! loop over the new snow layers
! Summ or avergage of the constituting quantities of the old snow layers
! which are totally or partially inserted in the new snow layer
 CALL get_mass_heat(kj,inlvls_new,inlvls_old,                                &
                    zsnowztop_old,zsnowztop_new,zsnowzbot_old,zsnowzbot_new, &
                    zsnowrhoo,zsnowdzo,zsnowgran1o,zsnowgran2o,zsnowhisto,   &
                    zsnowageo,zsnowheato,                                    &
                    zsnowrhon,psnowdzn,zsnowgran1n,zsnowgran2n,zsnowhistn,   &
                    zsnowagen,zsnowheatn,hsnowmetamo                         )     
!
! check of consistency between new and old snowpacks
zsnowhean  = 0.
zmastotn   = 0.
zsnowheao  = 0.
zmastoto   = 0.    
zpsnow_new = 0.
zpsnow_old = 0.
!
DO jst = 1,inlvls_new
  zsnowhean  = zsnowhean  + zsnowheatn(jst)
  zmastotn   = zmastotn   + zsnowrhon(jst) * psnowdzn(jst)
  zpsnow_new = zpsnow_new + psnowdzn(jst)
ENDDO
!
DO jst = 1,inlvls_old
  zsnowheao  = zsnowheao  + zsnowheato(jst)
  zmastoto   = zmastoto   + zsnowrhoo(jst) * zsnowdzo(jst)
  zpsnow_old = zpsnow_old + zsnowdzo(jst)
ENDDO
!
IF ( abs( zsnowhean-zsnowheao )>0.0001 .OR. abs( zmastotn-zmastoto )>0.0001 .OR. &
     abs( zpsnow_new-zpsnow_old )> 0.0001 ) THEN 
  WRITE(*,*) 'Warning diff', zsnowhean-zsnowheao,zmastotn-zmastoto,zpsnow_new-zpsnow_old
ENDIF
!
! 5. Update mass (density and thickness) and heat:
! ------------------------------------------------
!
psnowdz(:) = psnowdzn(:)
!
psnowrho(:) = zsnowrhon(:)
psnowheat(:) = zsnowheatn(:)
psnowgran1(:) = zsnowgran1n(:)
psnowgran2(:) = zsnowgran2n(:)
psnowhist(:) =  zsnowhistn(:)
!
psnowage(:) =  zsnowagen(:)
!
IF (lhook) CALL dr_hook('SNOWNLGRIDFRESH_1D',1,zhook_handle)
!
END SUBROUTINE snownlgridfresh_1d
!####################################################################
!####################################################################
!###################################################################
SUBROUTINE snowdrift(PTSTEP,PVMOD,PSNOWRHO,PSNOWDZ,PSNOW,        &
                     psnowgran1,psnowgran2,psnowhist,knlvls_use, &
                     pta,pqa,pps,prhoa,pz0eff,puref,             &
                     osnowdrift_sublim,hsnowmetamo,psndrift      )
!
!!    PURPOSE
!!    -------
!     Snow compaction  and metamorphism due to drift
!     Mass is unchanged: layer thickness is reduced
!     in proportion to density increases. Method inspired from
!     Brun et al. (1997) and Guyomarch
!     
!     - computes a mobility index of each snow layer from its grains, density
!                 and history
!     - computes a drift index of each layer from its mobility and wind speed
!     - computes a transport index with an exponential decay taking into
!                 account its depth and the mobility of upper layers
!     - increases density and changes grains in case of transport
!
!     HISTORY:
!     Basic parameterization from Crocus/ARPEGE Coupling (1997) 
!     Implementation in V5 
!     Insertion in V6 of grains type evolution in case of dendritic snow (V.
!     Vionnet) 
!     07/2012 (for V7.3): E. Brun, M. Lafaysse : optional sublimation of drifted snow
!     2012-09-20 : bug correction : ZFF was not computed if LSNOWDRIFT_SUBLIM=FALSE.
!
!     2014-02-05 V. Vionnet: systematic use of 5m wind speed to compute drift index
!     2014-06-03 M. Lafaysse: threshold on PZ0EFF

USE modd_csts,ONLY : xtt
USE mode_thermos

USE modd_snow_par, ONLY : xvtime, xvromax, xvromin, xvmob1,  &
                          xvmob2, xvmob3, xvmob4, xvdrift1, xvdrift2, xvdrift3, &
                          xvsizemin, xcoef_ff, xcoef_effect, xqs_ref
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!
REAL, DIMENSION(:), INTENT(IN)      :: pta, pqa, pps, prhoa
!
REAL, DIMENSION(:), INTENT(IN)      :: pvmod
!
INTEGER, DIMENSION(:), INTENT(IN)   :: knlvls_use
!
REAL, DIMENSION(:),INTENT(IN)       :: pz0eff,puref
!
LOGICAL,INTENT(IN)                  :: osnowdrift_sublim
!
 CHARACTER(3), INTENT(IN)            :: hsnowmetamo ! metamorphism scheme
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowrho, psnowdz,psnowgran1, &
                                       psnowgran2,psnowhist
REAL, DIMENSION(:), INTENT(OUT)     :: psnow
REAL, DIMENSION(:), INTENT(OUT)     :: psndrift !blowing snow sublimation (kg/m2/s)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: zsnowrho2
REAL, DIMENSION(SIZE(PSNOWRHO,1)                 ) :: zsnowdz1
!
REAL, DIMENSION(SIZE(PSNOWRHO,1))   :: zqsati, zff ! QS wrt ice, gust speed
!
REAL     :: zz0eff 
!
REAL     :: zprofequ, zrmob, zrdrift, zrt, zdro, zdgr1, zdgr2
REAL     :: zvt ! 5m wind speed threshold for surface
!transport
REAL     :: zqs_effect ! effect of QS on snow
REAL     :: zwind_effect ! effect of wind on snow
REAL     :: zdrift_effect ! effect of QS and wind on snow
! transformation
REAL     :: zqs !Blowing snow sublimation (kg/m2/s) 
REAL     :: zrhi, zfact
!
INTEGER  :: jj,jst   ! looping indexes
!
REAL(KIND=JPRB) :: zhook_handle
!
! Reference height for wind speed used to dertermine the occurrence of blowing snow
REAL, PARAMETER :: pphref_wind=5.
REAL, PARAMETER :: pphref_min=pphref_wind/2.
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('SNOWDRIFT',0,zhook_handle)
!
! 0. Initialization:
! ------------------
!
zsnowdz1(:) = psnowdz(:,1)
!
DO jj = 1,SIZE(psnow)
  DO jst = 1,knlvls_use(jj)
    zsnowrho2(jj,jst) = psnowrho(jj,jst)
  ENDDO
ENDDO
!
IF ( osnowdrift_sublim ) THEN
  zqsati(:) = qsati( pta(:),pps(:) )
END IF
!
! 1. Computation of drift and induced settling and metamorphism
! ------------------
!
DO jj=1, SIZE(psnow)
  !
  ! gust speed at 5m above the snowpack
  ! Computed from PVMOD at PUREF (m) assuming a log profile in the SBL 
  ! and a roughness length equal to PZ0EFF
  zz0eff=min(pz0eff(jj),puref(jj)*0.5,pphref_min)
  zff(jj) = xcoef_ff*pvmod(jj)*log(pphref_wind/zz0eff)/log(puref(jj)/zz0eff)
  !
  ! initialization decay coeff
  zprofequ = 0.
  !
  DO jst = 1,knlvls_use(jj)
    !
    zfact = 1.25 - 1.25 * ( max( psnowrho(jj,jst), xvromin ) - xvromin )/1000./xvmob1
    !
    IF ( hsnowmetamo=='B92' ) THEN
      !
      !  mobility index computation of a layer as a function of its properties
      IF( psnowgran1(jj,jst)<0. ) THEN        
        ! dendritic case              
        zrmob = 0.34 * ( 0.5  - ( 0.75*psnowgran1(jj,jst) + 0.5*psnowgran2(jj,jst) )/99. ) + &
                0.66 * zfact
      ELSE
        ! non dendritic case         
        zrmob = 0.34 * ( xvmob2 - xvmob2*psnowgran1(jj,jst)/99. - xvmob3*psnowgran2(jj,jst)*1000. ) + &
                0.66 * zfact
      ENDIF
      !
    ELSE
      !
      IF ( psnowgran1(jj,jst)<xvdiam6*(4.-psnowgran2(jj,jst))-xuepsi ) THEN  
        ! dendritic case
        zrmob = 0.34 * ( 0.5 + 0.75 * &
               ( psnowgran1(jj,jst)/xvdiam6-4.+psnowgran2(jj,jst) )/( psnowgran2(jj,jst)-3. ) &
                             - 0.5 * psnowgran2(jj,jst) ) + &
                0.66 * zfact  
      ELSE
        ! non dendritic case
        zrmob = 0.34 * ( xvmob2 - xvmob2 * psnowgran2(jj,jst) &
                                - xvmob3 * (4.-psnowgran2(jj,jst))*xvdiam6*1000. ) + &
                0.66 * zfact
      ENDIF      
      !
    ENDIF
    !
    ! correction in case of former wet snow
    IF ( psnowhist(jj,jst) >= 2. ) zrmob = min(zrmob, xvmob4)
    !      
    ! computation of drift index supposing no overburden snow
    zrdrift = zrmob - ( xvdrift1 * exp( -xvdrift2*zff(jj) ) - 1.)
    ! modif_EB exit loop if there is no drift
    IF ( zrdrift<=0. ) EXIT
    !    
    ! update the decay coeff by half the current layer
    zprofequ = zprofequ + 0.5 * psnowdz(jj,jst) * 0.1 * ( xvdrift3 - zrdrift )
    ! computation of the drift index inclunding the decay by overburden snow 
    zrt = max( 0., zrdrift * exp( -zprofequ*100 ) )
    !
    IF ( osnowdrift_sublim .AND. jst==1 ) THEN
      !Specific case for blowing snow sublimation
      ! computation of wind speed threshold QSATI and RH withe respect to ice
      zvt  = -log( (zrmob+1.)/xvdrift1 ) / xvdrift2
      zrhi = pqa(jj) / zqsati(jj)
      ! computation of sublimation rate according to Gordon's PhD
      zqs = 0.0018 * (xtt/pta(jj))**4 * zvt * prhoa(jj) * zqsati(jj) * &
            (1.-zrhi) * (zff(jj)/zvt)**3.6
      ! WRITE(*,*) 'surface Vt vent*coef  ZRDRIFT ZRMOB :',ZVT,&
      ! ZFF(JJ),ZRDRIFT,ZRMOB 
      ! WRITE(*,*) 'V>Vt ZQS   :',ZQS
      ! surface depth decrease in case of blowing snow sublimation
      ! WRITE(*,*) 'V>Vt DSWE DZ Z:',- MAX(0.,ZQS)*PTSTEP/COEF_FF,
      ! - MAX(0.,ZQS)*PTSTEP/COEF_FF/PSNOWRHO(JJ,JST),PSNOWDZ(JJ,JST)
      ! 2 lignes ci-dessous a valider pour avoir sublim drift
      psnowdz(jj,jst) = max( 0.5*psnowdz(jj,jst), &
                             psnowdz(jj,jst) - max(0.,zqs) * ptstep/xcoef_ff/psnowrho(jj,jst) )
      psndrift(jj) = (zsnowdz1(jj)-psnowdz(jj,jst))*psnowrho(jj,jst)/ptstep
    ELSE
      zqs = 0.
    END IF
    !
    zqs_effect    = min( 3., max( 0.,zqs )/xqs_ref ) * zrt
    zwind_effect  = xcoef_effect * zrt
    zdrift_effect = ( zqs_effect + zwind_effect ) * ptstep / xcoef_ff / xvtime
    ! WRITE(*,*) 'ZQS_EFFECT,ZWIND_EFFECT,ZDRIFT_EFFECT:',ZQS_EFFECT,ZWIND_EFFECT,ZDRIFT_EFFECT
    !
    ! settling by wind transport only in case of not too dense snow
    IF( psnowrho(jj,jst) < xvromax ) THEN
      zdro = zdrift_effect * ( xvromax - psnowrho(jj,jst) )
      psnowrho(jj,jst) = min( xvromax , psnowrho(jj,jst) + zdro )
      psnowdz(jj,jst) = psnowdz(jj,jst) * zsnowrho2(jj,jst) / psnowrho(jj,jst)
    ENDIF
    !
    IF ( hsnowmetamo=='B92' ) THEN
      !  
      ! metamorphism induced by snow drift
      IF ( psnowgran1(jj,jst)<0. ) THEN        
        ! dendritic case              
        zdgr1 = zdrift_effect * ( -psnowgran1(jj,jst) ) * 0.5
        psnowgran1(jj,jst) = psnowgran1(jj,jst) + min( zdgr1, -0.99 * psnowgran1(jj,jst) )
        ! modif_VV_140910          
        zdgr2 = zdrift_effect * ( 99. - psnowgran2(jj,jst) )
        psnowgran2(jj,jst) = min( 99., psnowgran2(jj,jst) + zdgr2 )
        ! fin modif_VV_140910          
      ELSE
        ! non dendritic case         
        zdgr1 = zdrift_effect * ( 99. - psnowgran1(jj,jst) )
        zdgr2 = zdrift_effect * 5. / 10000.
        psnowgran1(jj,jst) = min( 99., psnowgran1(jj,jst) + zdgr1 )
        psnowgran2(jj,jst) = max( xvsizemin, psnowgran2(jj,jst) - zdgr2 )
      ENDIF
      !
    ELSE
      !
      ! dendritic case        
      IF ( psnowgran1(jj,jst)<xvdiam6*(4.-psnowgran2(jj,jst))-xuepsi ) THEN  
        !
        zdgr1 = min( zdrift_effect * ( ( psnowgran1(jj,jst)/xvdiam6-4.+psnowgran2(jj,jst) )/ &
                                        (psnowgran2(jj,jst)-3.) ) * 0.5, &
                              0.99 * ( ( psnowgran1(jj,jst)/xvdiam6-4.+ psnowgran2(jj,jst))/ &
                                        (psnowgran2(jj,jst)-3.) ) )
        zdgr2 = zdrift_effect * ( 1.-psnowgran2(jj,jst) )
        !
        psnowgran1(jj,jst) = psnowgran1(jj,jst) + xvdiam6 * &
                             ( zdgr2 * ( (psnowgran1(jj,jst)/xvdiam6-1.)/(psnowgran2(jj,jst)-3.) ) - &
                               zdgr1 * ( psnowgran2(jj,jst)-3. ) )   
          psnowgran2(jj,jst) = min(1.,psnowgran2(jj,jst)+zdgr2)
      ! non dendritic case       
      ELSE
        !
        zdgr1 = zdrift_effect * 5./10000.
        zdgr2 = zdrift_effect * (1.-psnowgran2(jj,jst))
        !
        psnowgran1(jj,jst) = psnowgran1(jj,jst) - 2. * xvdiam6 * psnowgran2(jj,jst) * zdgr2
        psnowgran2(jj,jst) = min( 1., psnowgran2(jj,jst)+zdgr2 )
        !
      ENDIF        
      !
    ENDIF
    !
    ! update the decay coeff by half the current layer
    zprofequ = zprofequ + 0.5 * psnowdz(jj,jst) * 0.1 * ( xvdrift3 - zrdrift )
    !      
  ENDDO  ! snow layers loop
  !
ENDDO    ! grid points loop
!
! 2. Update total snow depth:
! -----------------------------------------------
!
! Compaction of total snowpack depth
!
DO jj = 1,SIZE(psnowdz,1)
  psnow(jj) = sum( psnowdz(jj,1:knlvls_use(jj)) )
ENDDO
!
IF (lhook) CALL dr_hook('SNOWDRIFT',1,zhook_handle)
!
END SUBROUTINE snowdrift
!####################################################################
!###################################################################
!####################################################################
!####################################################################
SUBROUTINE snowcrolayer_gone(PTSTEP,PSCAP,PSNOWTEMP,PSNOWDZ,          &
                             psnowrho,psnowliq,psnowgran1,psnowgran2, &
                             psnowhist,psnowage,ples3l,knlvls_use     ) 
!
!
!!    PURPOSE
!     Account for the case when one or several snow layers melt
!     during a time step:
!     in that case, merge these layers with the underlying layer
!     except for the bottom layer which is merged to the abovelying layer
!     energy and mass are conserved
!     a new merged layer keeps the grain, histo and age properties of the 
!     non-melted layer
!
USE modd_csts,ONLY : xtt, xlmtt, xrholw, xrholi, xlvtt, xci
!
USE mode_snow3l
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                     :: ptstep
!
REAL, DIMENSION(:,:), INTENT(INOUT)  :: pscap
!
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowdz, psnowtemp, psnowrho, psnowliq 
REAL, DIMENSION(:,:), INTENT(INOUT)  :: psnowgran1,psnowgran2,psnowhist,psnowage
!
INTEGER, DIMENSION(:), INTENT(INOUT) :: knlvls_use ! 
!
REAL, DIMENSION(:), INTENT(IN) :: ples3l
!
!*      0.2    declarations of local variables
!
REAL :: zheat, zmass, zdz, zliq, zsnowlwe
!
INTEGER :: jj,jst,jst_1, jst_2, jst_max, idiff_layer ! loop counter
INTEGER :: id_1, id_2
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('SNOWCROLAYER_GONE',0,zhook_handle)
!
DO jj=1,SIZE(psnowrho,1)  ! loop on gridpoints
  !
  jst_max = knlvls_use(jj)
  !
  idiff_layer = 0  ! used as shift counter of previously melted layers
  !     
  DO jst_1 = jst_max,1-jst_max,-1 ! loop on 2 x layers in case of multi melt
    !
    jst = jst_1 + idiff_layer
    !
    ! Merge is possible only in case of 2 active layers or more
    IF ( jst>=1 .AND. knlvls_use(jj)>1 ) THEN
      !
      ! Total Liquid equivalent water content of snow (m):
      zsnowlwe = psnowrho(jj,jst) * psnowdz(jj,jst) / xrholw
      !
      ! Consideration of sublimation if any       
      IF ( jst==1 ) zsnowlwe = zsnowlwe - max( 0., ples3l(jj)*ptstep/(xlstt*xrholw) )
      !
      ! Test if avalaible energy exceeds total latent heat
      IF ( pscap(jj,jst) * max( 0.0, psnowtemp(jj,jst)-xtt ) * psnowdz(jj,jst) >=  &
           ( ( zsnowlwe-psnowliq(jj,jst) ) * xlmtt * xrholw ) - xuepsi ) THEN
        !
        IF ( jst==knlvls_use(jj) ) THEN
          id_1 = jst-1
          id_2 = jst
        ELSE
          id_1 = jst
          id_2 = jst + 1
        ENDIF
        !
        ! Case of a total melt of the bottom layer: merge with above layer 
        !        which keeps its grain, histo and age properties
        zheat = 0.
        zmass = 0.
        zdz   = 0.
        zliq  = 0.
        DO jst_2 = id_1,id_2
          zheat = zheat + &
                  psnowdz(jj,jst_2) * &
                  ( pscap(jj,jst_2)*( psnowtemp(jj,jst_2)-xtt ) - xlmtt*psnowrho(jj,jst_2) ) + &
                  xlmtt * xrholw * psnowliq(jj,jst_2)
          zmass = zmass + psnowdz(jj,jst_2) * psnowrho(jj,jst_2)
          zdz   = zdz   + psnowdz(jj,jst_2)
          zliq  = zliq  + psnowliq(jj,jst_2)
        ENDDO
        !
        psnowdz(jj,id_1) = zdz
        psnowrho(jj,id_1) = zmass / zdz
        psnowliq(jj,id_1) = zliq
        !
        ! Temperature of the merged layer is deduced from the heat content 
        pscap(jj,id_1) = ( psnowrho(jj,id_1) - &
                               psnowliq(jj,id_1) * xrholw / &
                               max( psnowdz(jj,id_1),xsnowdzmin ) ) * xci
        psnowtemp(jj,id_1) = xtt + &
          ( ( ( ( zheat - xlmtt*xrholw*psnowliq(jj,id_1) ) / psnowdz(jj,id_1) ) + &
              xlmtt*psnowrho(jj,id_1) ) &
            / pscap(jj,id_1) )
        !
        IF( jst/=knlvls_use(jj) ) THEN
          !
          psnowgran1(jj,jst) = psnowgran1(jj,jst+1)
          psnowgran2(jj,jst) = psnowgran2(jj,jst+1)
          psnowhist(jj,jst) = psnowhist(jj,jst+1)
          psnowage(jj,jst) = psnowage(jj,jst+1)
          !
          ! Shift the above layers
          DO jst_2 = jst+1,knlvls_use(jj)-1
            psnowtemp(jj,jst_2) = psnowtemp(jj,jst_2+1)
            pscap(jj,jst_2) = pscap(jj,jst_2+1)
            psnowdz(jj,jst_2) = psnowdz(jj,jst_2+1)
            psnowrho(jj,jst_2) = psnowrho(jj,jst_2+1)
            psnowliq(jj,jst_2) = psnowliq(jj,jst_2+1)
            psnowgran1(jj,jst_2) = psnowgran1(jj,jst_2+1)
            psnowgran2(jj,jst_2) = psnowgran2(jj,jst_2+1)
            psnowhist(jj,jst_2) = psnowhist(jj,jst_2+1)
            psnowage(jj,jst_2) = psnowage(jj,jst_2+1)
          ENDDO !  loop JST_2
          !        
          ! Update the shift counter IDIFF_LAYER
          idiff_layer = idiff_layer + 1
          !
        ENDIF ! end test of bottom layer
        !
        ! Decrease the number of active snow layers              
        knlvls_use(jj) = knlvls_use(jj) - 1
        !
      ENDIF ! end test on availibility of energy
      !
    ENDIF ! end test on the number of remaining active layers
    !
  ENDDO ! end loop on the snow layers
  !
ENDDO ! end loop gridpoints
!
IF (lhook) CALL dr_hook('SNOWCROLAYER_GONE',1,zhook_handle)
!
END SUBROUTINE snowcrolayer_gone
!####################################################################
!###################################################################
!####################################################################
!###################################################################
SUBROUTINE snowcroprintprofile(HINFO,KLAYERS,OPRINTGRAN,PSNOWDZ,PSNOWRHO, &
                               psnowtemp,psnowliq,psnowheat,psnowgran1,   &
                               psnowgran2,psnowhist,psnowage,hsnowmetamo  )
!
! Matthieu Lafaysse 08/06/2012
! This routine prints the snow profile of a given point for debugging
!
!to compute SSA
USE modd_csts, ONLY : xrholi
USE modd_snow_par, ONLY : xd1, xd2, xd3, xx
!
IMPLICIT NONE
!
 CHARACTER(*), INTENT(IN) :: hinfo
LOGICAL,       INTENT(IN) :: oprintgran
INTEGER,       INTENT(IN) :: klayers
REAL, DIMENSION(:), INTENT(IN) :: psnowdz,psnowrho,psnowtemp,psnowliq, &
                                  psnowheat,psnowgran1,psnowgran2,     &
                                  psnowhist,psnowage
 CHARACTER(3), INTENT(IN)       :: hsnowmetamo
!
REAL, DIMENSION(KLAYERS) :: zsnowssa
REAL :: zdiam
!
INTEGER :: jst
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCROPRINTPROFILE',0,zhook_handle)
!
WRITE(*,*)
WRITE(*,*)trim(hinfo)
!
IF (oprintgran) THEN
  !
  ! Compute SSA from SNOWGRAN1 and SNOWGRAN2
  IF ( hsnowmetamo=='B92' ) THEN
    !
    DO jst = 1,klayers
      !
      IF ( psnowgran1(jst)<0. ) THEN
        zdiam =  -psnowgran1(jst)*xd1/xx + (1.+psnowgran1(jst)/xx) * &
                ( psnowgran2(jst)*xd2/xx + (1.-psnowgran2(jst)/xx) * xd3 ) 
        zdiam = zdiam/10000.
      ELSE
        zdiam = psnowgran2(jst)*psnowgran1(jst)/xx + &
                max( 0.0004, 0.5*psnowgran2(jst) ) * ( 1.-psnowgran1(jst)/xx )                
      ENDIF
      zsnowssa(jst) = 6. / (xrholi*zdiam)
      !
    END DO
    !
  ELSE
    !
    zsnowssa = 6. / (xrholi*psnowgran1)
    !
  ENDIF 
  !   
  WRITE(*,'(9(A12,"|"))')"-------------","-------------","-------------",&
        "-------------","-------------","-------------","-------------",&
        "-------------","-------------"
  WRITE(*,'(9(A12,"|"))')"PSNOWDZ","PSNOWRHO","PSNOWTEMP","PSNOWLIQ","PSNOWHEAT",&
        "PSNOWGRAN1","PSNOWGRAN2","PSNOWHIST","PSNOWAGE"
  WRITE(*,'(9(A12,"|"))')"-------------","-------------","-------------",&
        "-------------","-------------","-------------","-------------",&
        "-------------","-------------"
  DO jst = 1,klayers
    WRITE(*,'(9(ES12.3,"|")," L",I2.2)') psnowdz(jst),psnowrho(jst),psnowtemp(jst),    &
                                          psnowliq(jst),psnowheat(jst),psnowgran1(jst), &
                                          psnowgran2(jst),psnowhist(jst),psnowage(jst),jst
  ENDDO
  WRITE(*,'(9(A12,"|"))')"-------------","-------------","-------------",&
        "-------------","-------------","-------------","-------------",&
        "-------------","-------------"
  !
ELSE
  !
  WRITE(*,'(5(A12,"|"))')"------------","------------","------------",&
        "------------","------------"
  WRITE(*,'(5(A12,"|"))')"PSNOWDZ","PSNOWRHO","PSNOWTEMP","PSNOWLIQ","PSNOWHEAT"
  WRITE(*,'(5(A12,"|"))')"------------","------------","------------",&
        "------------","------------"
  DO jst = 1,klayers
    WRITE(*,'(5(ES12.3,"|")," L",I2.2)') psnowdz(jst),psnowrho(jst),psnowtemp(jst),&
                                         psnowliq(jst),psnowheat(jst),jst
  ENDDO
  WRITE(*,'(5(A12,"|"))')"------------","------------","------------",&
        "------------","------------"
  !
END IF
!
WRITE(*,*)
!
IF (lhook) CALL dr_hook('SNOWCROPRINTPROFILE',1,zhook_handle)
!
END SUBROUTINE snowcroprintprofile
!####################################################################
!###################################################################
SUBROUTINE snowcroprintatm(CINFO,PTA,PQA,PVMOD,PRR,PSR,PSW_RAD,PLW_RAD, &
                           ptg, psoilcond,pd_g,ppsn3l                 )

! Matthieu Lafaysse 08/06/2012
! This routine prints the atmospheric forcing of a given point for debugging
! and ground data

IMPLICIT NONE

 CHARACTER(*), INTENT(IN) :: cinfo
REAL,          INTENT(IN) :: pta,pqa,pvmod,prr,psr,psw_rad,plw_rad
REAL,          INTENT(IN) :: ptg, psoilcond, pd_g, ppsn3l
!
INTEGER :: jst
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCROPRINTATM',0,zhook_handle)
!
 CALL snowcroprintdate()
!
WRITE(*,*)
WRITE(*,*)trim(cinfo)
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
WRITE(*,'(4(A12,"|"))')"PTA","PQA","PRR","PSR"
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
WRITE(*,'(4(ES12.3,"|")," meteo1")')pta,pqa,prr,psr
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
WRITE(*,'(3(A12,"|"))')"------------","------------","------------"
WRITE(*,'(3(A12,"|"))')"PSW_RAD","PLW_RAD","PVMOD"
WRITE(*,'(3(A12,"|"))')"------------","------------","------------"
WRITE(*,'(3(ES12.3,"|")," meteo2")')psw_rad,plw_rad,pvmod
WRITE(*,'(3(A12,"|"))')"------------","------------","------------"
WRITE(*,*)
WRITE(*,*)"Ground :"
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
WRITE(*,'(4(A12,"|"))')"PTG","PSOILCOND","PD_G","PPSN3L"
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
WRITE(*,'(4(ES12.3,"|")," soil")')ptg,psoilcond,pd_g,ppsn3l
WRITE(*,'(4(A12,"|"))')"------------","------------","------------",&
"------------"
!
IF (lhook) CALL dr_hook('SNOWCROPRINTATM',1,zhook_handle)
!
END SUBROUTINE snowcroprintatm
!
!####################################################################
SUBROUTINE snowcrostopbalance(PMASSBALANCE,PENERGYBALANCE)
!
USE mode_crodebug, ONLY : xwarning_massbalance, xwarning_energybalance
!
USE modi_abor1_sfx
!
! stop if energy and mass balances are not closed
!
IMPLICIT NONE
!
REAL , DIMENSION(:), INTENT(IN) :: pmassbalance, penergybalance
!
REAL,DIMENSION(SIZE(PSR)) :: zmassbalance,zenergybalance
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCROSTOPBALANCE',0,zhook_handle)
!
IF ( any( pmassbalance   > xwarning_massbalance   ) ) &
        CALL abor1_sfx("SNOWCRO: WARNING MASS BALANCE !")
IF ( any( penergybalance > xwarning_energybalance ) ) &
        CALL abor1_sfx("SNOWCRO: WARNING ENERGY BALANCE !")
!
IF (lhook) CALL dr_hook('SNOWCROSTOPBALANCE',1,zhook_handle)
!
END SUBROUTINE snowcrostopbalance
!
!###################################################################
SUBROUTINE snowcroprintbalance(PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN, &
                               psr,prr,pthrufal,pevap,pevapcor,pgrndflux,phsnow,    &
                               prnsnow,plel3l,ples3l,phpsnow,psnowhmass,psnowdz,    &
                               ptstep,pmassbalance,penergybalance,pevapcor2         )
!
! Matthieu Lafaysse / Eric Brun 03/10/2012
! Print energy and mass balances.
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: psummass_ini,psumheat_ini,psummass_fin,psumheat_fin
REAL, INTENT(IN) :: psr,prr,pthrufal,pevap,pevapcor
REAL, INTENT(IN) :: pgrndflux,phsnow,prnsnow,plel3l,ples3l,phpsnow,psnowhmass
REAL, INTENT(IN) :: psnowdz !first layer
REAL, INTENT(IN) :: ptstep !time step
REAL, INTENT(IN) :: pmassbalance, penergybalance, pevapcor2
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCROPRINTBALANCE',0,zhook_handle)
!
WRITE(*,*) ' '
WRITE(*,fmt='(A1,67("+"),A1)') "+","+"
!
 CALL snowcroprintdate()
!
WRITE(*,*) ' '
!
! print des residus de bilan et des differents termes pour le point
WRITE (*,fmt="(A25,1x,E17.10)") 'final mass (kg/m2) =' , psummass_fin
WRITE (*,fmt="(A25,1x,E17.10)") 'final energy (J/m2) =', zsumheat_fin
WRITE(*,*) ' '
!
WRITE(*,fmt="(A25,1x,E17.10)") 'mass balance (kg/m2) =', pmassbalance
!
WRITE(*,*) ' '
WRITE(*,fmt="(A35)") 'mass balance contribution (kg/m2) ' 
WRITE(*,fmt="(A51,1x,E17.10)") 'delta mass:', (psummass_fin-psummass_ini)
WRITE(*,fmt="(A51,1x,E17.10)") 'hoar or condensation (>0 towards snow):', -pevap * ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'rain:',    prr * ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'snow:',    psr * ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'run-off:', pthrufal * ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'evapcor:', pevapcor * ptstep
!
WRITE(*,fmt='(A1,55("-"),A1)')"+","+"
WRITE(*,*) ' '
!
WRITE(*,fmt="(A25,4(1x,E17.10))") 'energy balance (W/m2)=',penergybalance
!
WRITE(*,*) ' '
WRITE(*,fmt="(A55)") 'energy balance contribution (W/m2) >0 towards snow :' 
WRITE(*,fmt="(A51,1x,E17.10)") 'delta heat:', (zsumheat_fin-zsumheat_ini)/ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'radiation (LW + SW):', prnsnow
WRITE(*,fmt="(A51,1x,E17.10)") 'sensible flux :',      -phsnow
WRITE(*,fmt="(A51,1x,E17.10)") 'ground heat flux :',   -pgrndflux
WRITE(*,fmt="(A51,1x,E17.10)") 'liquid latent flux:',  -plel3l
WRITE(*,fmt="(A51,1x,E17.10)") 'solid latent flux:',   -ples3l
WRITE(*,fmt="(A51,1x,E17.10)") 'rain sensible heat:',  phpsnow
WRITE(*,fmt="(A51,1x,E17.10)") 'snowfall/hoar heat (sensible + melt heat):', psnowhmass/ptstep
WRITE(*,fmt="(A51,1x,E17.10)") 'evapcor:', pevapcor2
WRITE(*,fmt='(A1,67("+"),A1)')"+","+"
!
IF (lhook) CALL dr_hook('SNOWCROPRINTBALANCE',1,zhook_handle)
!
END SUBROUTINE snowcroprintbalance
!
!####################################################################
SUBROUTINE get_balance(PSUMMASS_INI,PSUMHEAT_INI,PSUMMASS_FIN,PSUMHEAT_FIN, &
                       psr,prr,pthrufal,pevap,pevapcor,pgrndflux,phsnow,    &
                       prnsnow,plel3l,ples3l,phpsnow,psnowhmass,psnowdz,    &
                       ptstep,pmassbalance,penergybalance,pevapcor2         )
!
IMPLICIT NONE
!
REAL, INTENT(IN) :: psummass_ini,psumheat_ini,psummass_fin,psumheat_fin
REAL, INTENT(IN) :: psr,prr,pthrufal,pevap,pevapcor
REAL, INTENT(IN) :: pgrndflux,phsnow,prnsnow,plel3l,ples3l,phpsnow,psnowhmass
REAL, INTENT(IN) :: psnowdz !first layer
REAL, INTENT(IN) :: ptstep !time step
!
REAL, INTENT(OUT)  :: pmassbalance, penergybalance, pevapcor2
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_BALANCE',0,zhook_handle)
!
pmassbalance = psummass_fin - psummass_ini - &
               ( psr + prr - pthrufal - pevap + pevapcor ) * ptstep
!
pevapcor2 = pevapcor * psnowdz / max( xuepsi,psnowdz ) *  &
           ( abs(plel3l) * xlvtt / max( xuepsi,abs(plel3l) ) + &
             abs(ples3l) * xlstt / max( xuepsi,abs(ples3l) ) ) 
!
penergybalance = ( psumheat_fin-psumheat_ini ) / ptstep - &
                 ( -pgrndflux - phsnow + prnsnow - plel3l - ples3l + phpsnow ) - &
                 psnowhmass / ptstep - pevapcor2
!
IF (lhook) CALL dr_hook('SNOWCRO:GET_BALANCE',1,zhook_handle)
!
END SUBROUTINE get_balance
!
!###################################################################
SUBROUTINE snowcroprintdate()
!
IMPLICIT NONE
!
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOWCROPRINTDATE',0,zhook_handle)
!
WRITE(*,fmt='(I4.4,2("-",I2.2)," Hour=",F5.2)') &
  tptime%TDATE%YEAR, tptime%TDATE%MONTH, tptime%TDATE%DAY, tptime%TIME/3600.
!
IF (lhook) CALL dr_hook('SNOWCROPRINTDATE',1,zhook_handle)
!
END SUBROUTINE snowcroprintdate
!####################################################################
!###################################################################
!
END SUBROUTINE snowcro