snow3L_isba.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 snow3l_isba(HISBA, HSNOW_ISBA, HSNOWRES, OMEB, OGLACIER, HIMPLICIT_WIND,          &
                         tptime, ptstep, pvegtype,                                           &
                         psnowswe, psnowheat, psnowrho, psnowalb,                            &
                         psnowgran1, psnowgran2, psnowhist,psnowage,                         &
                         ptg, pcg, pct, psoilhcapz, psoilcondz,                              &
                         pps, pta, psw_rad, pqa, pvmod, plw_rad, prr, psr,                   &
                         prhoa, puref, pexns, pexna, pdircoszw, plvtt, plstt,                &
                         pzref, pz0nat, pz0eff, pz0hnat, palb, pd_g, pdzg,                   &
                         ppew_a_coef, ppew_b_coef,                                           &
                         ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef,                 &
                         pthrufal, pgrndflux, pflsn_cor, pgsfcsnow, pevapcor,                &
                         pswnetsnow, pswnetsnows, plwnetsnow,                                &
                         prnsnow, phsnow, pgfluxsnow, phpsnow, ples3l, plel3l, pevap,        &
                         psndrift, pustarsnow, ppsn, psrsfc, prrsfc, psnowsfch,              &
                         pdelheatn, pdelheatn_sfc,                                           &
                         pemisnow, pcdsnow, pchsnow, psnowtemp, psnowliq, psnowdz,           &
                         psnowhmass, pri, pzenith, pdelheatg, pdelheatg_sfc, plat, plon, pqs,&
                         osnowdrift,osnowdrift_sublim,osnow_abs_zenith,                      &
                         hsnowmetamo, hsnowrad                                               )                               
!     ######################################################################################
!
!!****  *SNOW3L_ISBA*  
!!
!!    PURPOSE
!!    -------
!
!     3-Layer snow scheme option (Boone and Etchevers 1999)
!     This routine is NOT called as snow depth goes below
!     a critical threshold which is vanishingly small.
!     This routine acts as an interface between SNOW3L and ISBA.
!     
!!**  METHOD
!!    ------
!
!     Direct calculation
!
!!    EXTERNAL
!!    --------
!
!     None
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Boone and Etchevers (1999)
!!    Belair (1995)
!!    Noilhan and Planton (1989)
!!    Noilhan and Mahfouf (1996)
!!      
!!    AUTHOR
!!    ------
!!      A. Boone           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original        7/99  Boone
!!      Packing added   4/00  Masson & Boone
!!      z0h and snow    2/06  LeMoigne
!!
!!      Modified by B. Decharme  (03/2009): Consistency with Arpege permanent
!!                                          snow/ice treatment
!!      Modified by A. Boone     (04/2010): Implicit coupling with atmosphere permitted.
!!
!!      Modified by B. Decharme  (04/2010): check suspicious low temperature for ES and CROCUS
!!      Modified by B. Decharme  (08/2013): Qsat as argument (needed for coupling with atm)
!!      Modified by A. Boone     (10/2014): MEB: pass in fluxes when using MEB
!!      Modified by B. Decharme  (03/2016): No snowdrift under forest
!!
!-------------------------------------------------------------------------------
!
USE modd_csts,       ONLY : xtt, xpi, xday, xlmtt, xlstt
USE modd_snow_par,   ONLY : xrhosmax_es, xsnowdmin, xrhosmin_es, xemissn
USE modd_surf_par,   ONLY : xundef
USE modd_type_date_surf, ONLY: date_time
!
USE modd_data_cover_par, ONLY : nvt_snow,                       &
                                nvt_tebd, nvt_trbe, nvt_bone,   &
                                nvt_trbd, nvt_tebe, nvt_tene,   &
                                nvt_bobd, nvt_bond, nvt_shrb
!
USE modi_snow3l
USE modi_snowcro
!
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!
!
!*      0.1    declarations of arguments
!
REAL, INTENT(IN)                    :: ptstep
!                                      PTSTEP    = time step of the integration
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pvegtype ! fraction of each vegetation
!
 CHARACTER(LEN=*),     INTENT(IN)    :: hisba
!                                      HISBA     = FLAG to use Force-Restore or DIFfusion
!                                      soil heat and mass transfer method
!
 CHARACTER(LEN=*),     INTENT(IN)    :: hsnow_isba
!                                      HSNOW_ISBA = FLAG to use SNOW3L or not 
!                                      (or default FR method)
!
 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)                 :: omeb       ! True = coupled to MEB. This means surface fluxes ae IMPOSED
!                                                 ! as an upper boundary condition to the explicit snow schemes. 
!                                                 ! If = False, then energy
!                                                 ! budget and fluxes are computed herein.
!
LOGICAL, INTENT(IN)                 :: oglacier   ! True = Over permanent snow and ice, 
!                                                     initialise WGI=WSAT,
!                                                     Hsnow>=10m and allow 0.8<SNOALB<0.85
!
 CHARACTER(LEN=*),     INTENT(IN)    :: himplicit_wind   ! wind implicitation option
!                                                       ! 'OLD' = direct
!                                                       ! 'NEW' = Taylor serie, order 1
!
TYPE(date_time), INTENT(IN)         :: tptime     ! current date and time
!
!
REAL, DIMENSION(:,:), INTENT(IN)    :: psoilhcapz, pd_g, pdzg
REAL, DIMENSION(:),   INTENT(IN)    :: pcg, pct, psoilcondz  
!                                      PD_G      = Depth to bottom of each soil layer (m)
!                                      PDZG      = Soil layer thicknesses (m)
!                                      PCG       = area-averaged soil heat capacity [(K m2)/J]
!                                      PCT       = area-averaged surface heat capacity [(K m2)/J]
!                                      PSOILCONDZ= soil thermal conductivity (W m-1 K-1)
!                                      PSOILHCAPZ= soil heat capacity (J m-3 K-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)      :: pzref, puref, pexns, pexna, pdircoszw, prhoa, pz0nat, pz0eff, pz0hnat, palb, &
                                       plvtt, plstt
!                                      PZ0EFF    = roughness length for momentum 
!                                      PZ0NAT    = grid box average roughness length
!                                      PZ0HNAT   = grid box average roughness length
!                                      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
!                                      PLVTT     = latent heat of vaporization-hydrology (J/kg)
!                                      PLSTT     = latent heat of sublimation-hydrology  (J/kg)
!
REAL, DIMENSION(:), INTENT(IN)      :: ppsn
!                                      PPSN     = Snow cover fraction (total) 
!
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
!                                      PPEW_B_COEF = wind coefficient
!                                      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) :: ptg
!                                      PTG       = Soil temperature profile (K)
!
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/m3)
!                                      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 layer(s) grain parameter 1
!                                      PSNOWGRAN2 = Snow layer(s) grain parameter 2
!                                      PSNOWHIST  = Snow layer(s) grain historical parameter
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowage  ! Snow grain age
!
!
REAL, DIMENSION(:), INTENT(INOUT)   :: prnsnow, phsnow, ples3l, plel3l,     &
                                       phpsnow, pemisnow, pevap, pgrndflux, pswnetsnow, &
                                       plwnetsnow, pswnetsnows, pdelheatg, pdelheatg_sfc
!                                      PLEL3L        = evaporation heat flux from snow (W/m2)
!                                      PLES3L        = 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)
!                                      PEMISNOW      = snow surface emissivity
!                                      PEVAP         = total evaporative flux from snow (kg/m2/s)
!                                      PGRNDFLUX     = soil/snow interface heat flux (W/m2)
!                                      PSWNETSNOW    = net shortwave radiation entering top of snowpack (W/m2)
!                                      PSWNETSNOWS   = net shortwave radiation in uppermost layer of snowpack (W/m2)
!                                                      (for surface energy budget closure diagnostics)
!                                      PLWNETSNOW    = net longwave radiation entering top of snowpack (W/m2)
!                                      PDELHEATG     = ground heat content change (diagnostic) (W/m2)
!                                                      note, modified if ground-snow flux adjusted
!                                      PDELHEATG_SFC = ground heat content change in sfc only (diagnostic) (W/m2)
!                                                      note, modified if ground-snow flux adjusted
!
!
REAL, DIMENSION(:), INTENT(OUT)     :: pgfluxsnow
!                                      PGFLUXSNOW    = net heat flux from snow (W/m2)
!
REAL, DIMENSION(:), INTENT(INOUT)   :: pustarsnow, pcdsnow, pchsnow, pri
!                                      PCDSNOW    = drag coefficient for momentum over snow (-)
!                                      PUSTARSNOW = friction velocity over snow (m/s)
!                                      PCHSNOW    = drag coefficient for heat over snow (-)
!                                      PRI        = Richardson number (-)
!
REAL, DIMENSION(:), INTENT(OUT)     :: pthrufal, pflsn_cor, pevapcor, psnowhmass, pgsfcsnow
!                                      PTHRUFAL  = rate that liquid water leaves snow pack: 
!                                                  paritioned into soil infiltration/runoff 
!                                                  by ISBA [kg/(m2 s)]
!                                      PFLSN_COR = soil/snow correction heat flux (W/m2) (not MEB)
!                                      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)]
!                                      PSNOWHMASS = heat content change due to mass
!                                                   changes in snowpack (J/m2): for budget
!                                                   calculations only.
!                                      PGSFCSNOW  = heat flux between the surface and sub-surface 
!                                                   snow layers (for energy budget diagnostics) (W/m2)
!
REAL, DIMENSION(:), INTENT(OUT)     :: psndrift
!                                      PSNDRIFT    = blowing snow sublimation (kg/m2/s)
!
REAL, DIMENSION(:), INTENT(OUT)     :: psrsfc, prrsfc, psnowsfch, pdelheatn, pdelheatn_sfc
!                                      PSRSFC = snow rate on soil/veg surface when SNOW3L in use
!                                      PRRSFC = rain rate on soil/veg surface when SNOW3L in use
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: psnowtemp
REAL, DIMENSION(:,:), INTENT(OUT)   :: psnowliq, 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)     :: pqs
!                                      PQS = surface humidity (kg/kg)
!
! ajout_EB pour prendre en compte angle zenithal du soleil dans LRAD
! puis plus tard dans LALB
REAL, DIMENSION(:), INTENT(IN)      :: pzenith    ! solar zenith angle
REAL, DIMENSION(:), INTENT(IN)      :: plat
REAL, DIMENSION(:), INTENT(IN)      :: plon
!
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
                                         !-----------------------
                                         ! Crocus 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
                                         !-----------------------
                                         ! Crocus 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, PARAMETER                     :: zcheck_temp = 50.0 
!                                      Limit to check suspicious low temperature (K)
!
INTEGER                             :: jwrk, jj ! Loop control
!
INTEGER                             :: inlvls   ! maximum number of snow layers
INTEGER                             :: inlvlg   ! number of ground layers
!
REAL, DIMENSION(SIZE(PTA))          :: zrrsnow, zsoilcond, zsnow, zsnowfall,  &
                                       zsnowablat_delta, zsnowswe_1d, zsnowd, & 
                                       zsnowh, zsnowh1, zgrndfluxn, zpsn,     &
                                       zsoilcor, zsnowswe_out, zthrufal,      &
                                       zsnow_mass_budget
!                                      ZSOILCOND    = soil thermal conductivity [W/(m K)]
!                                      ZRRSNOW      = rain rate over snow [kg/(m2 s)]
!                                      ZSNOW        = snow depth (m) 
!                                      ZSNOWFALL    = minimum equivalent snow depth
!                                                     for snow falling during the
!                                                     current time step (m)
!                                      ZSNOWABLAT_DELTA = FLAG =1 if snow ablates completely
!                                                     during current time step, else=0
!                                      ZSNOWSWE_1D  = TOTAL snowpack SWE (kg m-2)
!                                      ZSNOWD       = snow depth
!                                      ZSNOWH       = snow total heat content (J m-2)
!                                      ZSNOWH1      = snow surface layer heat content (J m-2)
!                                      ZGRNDFLUXN   = corrected snow-ground flux (if snow fully ablated during timestep)
!                                      ZPSN         = snow fraction working array
!                                      ZSOILCOR = for vanishingy thin snow cover,
!                                                 allow any excess evaporation
!                                                 to be extracted from the soil
!                                                 to maintain an accurate water
!                                                 balance [kg/(m2 s)]
!                                      ZSNOW_MASS_BUDGET = snow water equivalent budget (kg/m2/s)
!
!*      0.3    declarations of packed  variables
!
INTEGER                            :: isize_snow ! number of points where computations are done
INTEGER, DIMENSION(SIZE(PTA))      :: nmask      ! indices correspondance between arrays
!
LOGICAL, DIMENSION(SIZE(PTA))      :: lremove_snow
!
REAL(KIND=JPRB) :: zhook_handle
!
! - - ---------------------------------------------------
!
IF (lhook) CALL dr_hook('SNOW3L_ISBA',0,zhook_handle)
!
!*       0.     Initialize variables:
!               ---------------------
!
pflsn_cor(:)   = 0.0
pthrufal(:)    = 0.0
pevapcor(:)    = 0.0
psndrift(:)    = 0.0
psnowhmass(:)  = 0.0
psrsfc(:)      = psr(:)         ! these are snow and rain rates passed to ISBA,
prrsfc(:)      = prr(:)         ! so initialize here if SNOW3L not used:
pqs(:)         = xundef
!
zsnow(:)       = 0.0
zsnowd(:)      = 0.0
zgrndfluxn(:)  = 0.0
zsnowh(:)      = 0.0
zsnowh1(:)     = 0.0
zsnowswe_1d(:) = 0.0
zsnowswe_out(:)= 0.0
zsoilcond(:)   = 0.0
zrrsnow(:)     = 0.0
zsnowfall(:)   = 0.0
zsnowablat_delta(:) = 0.0
psnowliq(:,:)  = 0.0
psnowdz(:,:)   = 0.0
!
inlvls          = SIZE(psnowswe(:,:),2)                         
inlvlg          = min(SIZE(pd_g(:,:),2),SIZE(ptg(:,:),2))                         
!
IF(.NOT.omeb)THEN 
!
! If MEB activated, these values are input, else initialize here:
!
   pgrndflux(:)   = 0.0
   ples3l(:)      = 0.0
   plel3l(:)      = 0.0
   pevap(:)       = 0.0
   prnsnow(:)     = 0.0
   phsnow(:)      = 0.0
   pgfluxsnow(:)  = 0.0
   phpsnow(:)     = 0.0
   pswnetsnow(:)  = 0.0
   pswnetsnows(:) = 0.0
   plwnetsnow(:)  = 0.0
   pemisnow(:)    = xemissn
   pustarsnow(:)  = 0.0
   pcdsnow(:)     = 0.0
   pchsnow(:)     = 0.0
   pri(:)         = xundef
END IF
!
!
! Use ISBA-SNOW3L or NOT: NOTE that if explicit soil diffusion method in use,
! then *must* use explicit snow model:
!
IF (hsnow_isba=='3-L' .OR. hisba == 'DIF' .OR. hsnow_isba == 'CRO') THEN
!
! - Snow and rain falling onto the 3-L grid space:
!
   psrsfc(:)=0.0
!
   DO jj=1,SIZE(psr)
      zrrsnow(jj)        = ppsn(jj)*prr(jj)
      prrsfc(jj)         = prr(jj) - zrrsnow(jj)
      zsnowfall(jj)      = psr(jj)*ptstep/xrhosmax_es    ! maximum possible snowfall depth (m)
   ENDDO
!
! Calculate preliminary snow depth (m)

   zsnow(:)=0.
   zsnowh(:)=0.
   zsnowswe_1d(:)=0.
   zsnowh1(:)              = psnowheat(:,1)*psnowswe(:,1)/psnowrho(:,1) ! sfc layer only
!
   DO jwrk=1,SIZE(psnowswe,2)
      DO jj=1,SIZE(psnowswe,1)
         zsnowswe_1d(jj)     = zsnowswe_1d(jj) + psnowswe(jj,jwrk)
         zsnow(jj)           = zsnow(jj)       + psnowswe(jj,jwrk)/psnowrho(jj,jwrk)
         zsnowh(jj)          = zsnowh(jj)      + psnowheat(jj,jwrk)*psnowswe(jj,jwrk)/psnowrho(jj,jwrk)
      END DO
   ENDDO
!
   IF(hisba == 'DIF')THEN
      zsoilcond(:)   = psoilcondz(:)
   ELSE
!
! - Soil thermal conductivity
!   is implicit in Force-Restore soil method, so it
!   must be backed-out of surface thermal coefficients
!   (Etchevers and Martin 1997):
!
      zsoilcond(:)    = 4.*xpi/( pcg(:)*pcg(:)*xday/(pd_g(:,1)*pct(:)) )
!
   ENDIF
!
! ===============================================================
! === Packing: Only call snow model when there is snow on the surface
!              exceeding a minimum threshold OR if the equivalent
!              snow depth falling during the current time step exceeds 
!              this limit.
!
! counts the number of points where the computations will be made
!
!
   isize_snow = 0
   nmask(:) = 0
!
   DO jj=1,SIZE(zsnow)
      IF (zsnow(jj) >= xsnowdmin .OR. zsnowfall(jj) >= xsnowdmin) THEN
         isize_snow = isize_snow + 1
         nmask(isize_snow) = jj
      ENDIF
   ENDDO
!  
   IF (isize_snow>0) CALL call_model(isize_snow,inlvls,inlvlg,nmask)
!
! ===============================================================
!
! Remove trace amounts of snow and reinitialize snow prognostic variables
! if snow cover is ablated.
! If MEB used, soil T already computed, therefore correct heating/cooling
! effect of updated snow-soil flux
!
   zsnowd(:) = 0.
   zsnowswe_out(:) = 0.
   DO jwrk=1,SIZE(psnowswe,2)
      DO jj=1,SIZE(psnowswe,1)
         zsnowd(jj) = zsnowd(jj) + psnowswe(jj,jwrk)/psnowrho(jj,jwrk)
         zsnowswe_out(jj) = zsnowswe_out(jj) + psnowswe(jj,jwrk)
      ENDDO
   END DO
!
   lremove_snow(:)=(zsnowd(:)<xsnowdmin*1.1)
!
!
   IF(omeb)THEN
     zpsn(:)=1.0
   ELSE
!    To Conserve mass in ISBA without MEB, 
!    EVAP must be weignted by the snow fraction
!    in the calulation of THRUFAL
     zpsn(:)=ppsn(:)
   ENDIF
!
   zsnowablat_delta(:) = 0.0
   zthrufal(:) = pthrufal(:)
!
   WHERE(lremove_snow(:))
      zsnowswe_out(:)     = 0.0
      ples3l(:)           = min(ples3l(:), xlstt*(zsnowswe_1d(:)/ptstep + psr(:)))
      plel3l(:)           = 0.0
      pevap(:)            = ples3l(:)/plstt(:)
      pthrufal(:)         = max(0.0, zsnowswe_1d(:)/ptstep + psr(:) - pevap(:)*zpsn(:) + zrrsnow(:)) ! kg m-2 s-1
      zthrufal(:)         = max(0.0, zsnowswe_1d(:)/ptstep + psr(:) - pevap(:)         + zrrsnow(:)) ! kg m-2 s-1
      psrsfc(:)           = 0.0
      prrsfc(:)           = prrsfc(:)
      zsnowablat_delta(:) = 1.0
      psnowalb(:)         = xundef
      pevapcor(:)         = 0.0
      zsoilcor(:)         = 0.0
      pgfluxsnow(:)       = prnsnow(:) - phsnow(:) - ples3l(:) - plel3l(:)
      psnowhmass(:)       = -psr(:)*(xlmtt*ptstep)
      pgsfcsnow(:)        = 0.0
      pdelheatn(:)        = -zsnowh(:) /ptstep
      pdelheatn_sfc(:)    = -zsnowh1(:)/ptstep
      psnowsfch(:)        = pdelheatn_sfc(:) - (pswnetsnows(:) + plwnetsnow(:)    &
                          - phsnow(:) - ples3l(:) - plel3l(:)) + pgsfcsnow(:)     &
                          - psnowhmass(:)/ptstep 
      zgrndfluxn(:)       = (zsnowh(:)+psnowhmass(:))/ptstep + pgfluxsnow(:)
      ptg(:,1)            = ptg(:,1) + ptstep*pct(:)*zpsn(:)*(zgrndfluxn(:) - pgrndflux(:) - pflsn_cor(:))
      pdelheatg(:)        = pdelheatg(:)     + zpsn(:)*(zgrndfluxn(:) - pgrndflux(:) - pflsn_cor(:))
      pdelheatg_sfc(:)    = pdelheatg_sfc(:) + zpsn(:)*(zgrndfluxn(:) - pgrndflux(:) - pflsn_cor(:))
      pgrndflux(:)        = zgrndfluxn(:)
      pflsn_cor(:)        = 0.0
     END WHERE
!
!
   DO jwrk=1,inlvls
      DO jj=1,SIZE(psnowswe,1)
         psnowswe(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowswe(jj,jwrk)
         psnowheat(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowheat(jj,jwrk)
         psnowrho(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowrho(jj,jwrk)  + &
                                    zsnowablat_delta(jj)*xrhosmin_es  
         psnowtemp(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowtemp(jj,jwrk) + &
                                    zsnowablat_delta(jj)*xtt  
         psnowliq(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowliq(jj,jwrk)        
         psnowdz(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowdz(jj,jwrk)
         psnowage(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowage(jj,jwrk)
      ENDDO
   ENDDO
!  
   IF (hsnow_isba=='CRO') THEN
      DO jwrk=1,inlvls
         DO jj=1,SIZE(psnowgran1,1)
            psnowgran1(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowgran1(jj,jwrk) 
            psnowgran2(jj,jwrk)  = (1.0-zsnowablat_delta(jj))*psnowgran2(jj,jwrk)
            psnowhist(jj,jwrk)   = (1.0-zsnowablat_delta(jj))*psnowhist(jj,jwrk)
         ENDDO
      ENDDO
   ENDIF
!
!  ===============================================================
!
!  Compute snow mass budget 
!
   zsnow_mass_budget(:) = (zsnowswe_1d(:)-zsnowswe_out(:))/ptstep + psr(:)+zrrsnow(:) &
                                                                  - pevap(:)-zthrufal(:) &
                                                                  + pevapcor(:)+zsoilcor(:)
!
!
!  ===============================================================
!
!  To Conserve mass in ISBA, the latent heat flux part of 
!  the EVAPCOR term must be weignted by the snow fraction 
!
   pevapcor(:) = pevapcor(:)*zpsn(:) + zsoilcor(:)
!
! ===============================================================
!
! check suspicious low temperature
!
   DO jwrk=1,inlvls
      DO jj=1,SIZE(psnowswe,1)
         IF(psnowswe(jj,jwrk)>0.0.AND.psnowtemp(jj,jwrk)<zcheck_temp)THEN
            WRITE(*,*) 'Suspicious low temperature :',psnowtemp(jj,jwrk)
            WRITE(*,*) 'At point and location      :',jj,'LAT=',plat(jj),'LON=',plon(jj)
            WRITE(*,*) 'At snow level / total layer:',jwrk,'/',inlvls
            WRITE(*,*) 'SNOW MASS BUDGET (kg/m2/s) :',zsnow_mass_budget(jj)
            WRITE(*,*) 'SWE BY LAYER      (kg/m2)  :',psnowswe(jj,1:inlvls)
            WRITE(*,*) 'DEPTH BY LAYER      (m)    :',psnowdz(jj,1:inlvls)
            WRITE(*,*) 'DENSITY BY LAYER   (kg/m3) :',psnowrho(jj,1:inlvls)
            WRITE(*,*) 'TEMPERATURE BY LAYER (K)   :',psnowtemp(jj,1:inlvls)
            CALL abor1_sfx('SNOW3L_ISBA: Suspicious low temperature')                
         ENDIF
      ENDDO
   ENDDO
!
! ===============================================================
!
ENDIF
!
IF (lhook) CALL dr_hook('SNOW3L_ISBA',1,zhook_handle)
!
 CONTAINS
!
!================================================================
SUBROUTINE call_model(KSIZE1,KSIZE2,KSIZE3,KMASK)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: ksize1
INTEGER, INTENT(IN) :: ksize2
INTEGER, INTENT(IN) :: ksize3
INTEGER, DIMENSION(:), INTENT(IN) :: kmask
!
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowswe
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowdz
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowrho
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowheat
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowtemp
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowliq
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowgran1
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowgran2
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowhist
REAL, DIMENSION(KSIZE1,KSIZE2) :: zp_snowage
REAL, DIMENSION(KSIZE1)        :: zp_snowalb
REAL, DIMENSION(KSIZE1)        :: zp_swnetsnow
REAL, DIMENSION(KSIZE1)        :: zp_swnetsnows
REAL, DIMENSION(KSIZE1)        :: zp_lwnetsnow
REAL, DIMENSION(KSIZE1)        :: zp_ps
REAL, DIMENSION(KSIZE1)        :: zp_srsnow
REAL, DIMENSION(KSIZE1)        :: zp_rrsnow
REAL, DIMENSION(KSIZE1)        :: zp_psn3l
REAL, DIMENSION(KSIZE1)        :: zp_ta
REAL, DIMENSION(KSIZE1)        :: zp_ct
REAL, DIMENSION(KSIZE1,KSIZE3) :: zp_tg
REAL, DIMENSION(KSIZE1,KSIZE3) :: zp_d_g
REAL, DIMENSION(KSIZE1,KSIZE3) :: zp_dzg
REAL, DIMENSION(KSIZE1,KSIZE3) :: zp_soilhcapz
REAL, DIMENSION(KSIZE1)        :: zp_soild
REAL, DIMENSION(KSIZE1)        :: zp_delheatg
REAL, DIMENSION(KSIZE1)        :: zp_delheatg_sfc
REAL, DIMENSION(KSIZE1)        :: zp_sw_rad
REAL, DIMENSION(KSIZE1)        :: zp_qa
REAL, DIMENSION(KSIZE1)        :: zp_lvtt
REAL, DIMENSION(KSIZE1)        :: zp_lstt
REAL, DIMENSION(KSIZE1)        :: zp_vmod
REAL, DIMENSION(KSIZE1)        :: zp_lw_rad
REAL, DIMENSION(KSIZE1)        :: zp_rhoa
REAL, DIMENSION(KSIZE1)        :: zp_uref
REAL, DIMENSION(KSIZE1)        :: zp_exns
REAL, DIMENSION(KSIZE1)        :: zp_exna
REAL, DIMENSION(KSIZE1)        :: zp_dircoszw
REAL, DIMENSION(KSIZE1)        :: zp_zref
REAL, DIMENSION(KSIZE1)        :: zp_z0nat
REAL, DIMENSION(KSIZE1)        :: zp_z0hnat
REAL, DIMENSION(KSIZE1)        :: zp_z0eff
REAL, DIMENSION(KSIZE1)        :: zp_alb
REAL, DIMENSION(KSIZE1)        :: zp_soilcond
REAL, DIMENSION(KSIZE1)        :: zp_thrufal
REAL, DIMENSION(KSIZE1)        :: zp_grndflux
REAL, DIMENSION(KSIZE1)        :: zp_flsn_cor
REAL, DIMENSION(KSIZE1)        :: zp_gsfcsnow
REAL, DIMENSION(KSIZE1)        :: zp_evapcor
REAL, DIMENSION(KSIZE1)        :: zp_soilcor
REAL, DIMENSION(KSIZE1)        :: zp_gflxcor
REAL, DIMENSION(KSIZE1)        :: zp_rnsnow
REAL, DIMENSION(KSIZE1)        :: zp_hsnow
REAL, DIMENSION(KSIZE1)        :: zp_gfluxsnow
REAL, DIMENSION(KSIZE1)        :: zp_delheatn
REAL, DIMENSION(KSIZE1)        :: zp_delheatn_sfc
REAL, DIMENSION(KSIZE1)        :: zp_snowsfch
REAL, DIMENSION(KSIZE1)        :: zp_hpsnow
REAL, DIMENSION(KSIZE1)        :: zp_les3l
REAL, DIMENSION(KSIZE1)        :: zp_lel3l
REAL, DIMENSION(KSIZE1)        :: zp_evap
REAL, DIMENSION(KSIZE1)        :: zp_sndrift
REAL, DIMENSION(KSIZE1)        :: zp_ri
REAL, DIMENSION(KSIZE1)        :: zp_qs
REAL, DIMENSION(KSIZE1)        :: zp_emisnow
REAL, DIMENSION(KSIZE1)        :: zp_cdsnow
REAL, DIMENSION(KSIZE1)        :: zp_ustarsnow
REAL, DIMENSION(KSIZE1)        :: zp_chsnow
REAL, DIMENSION(KSIZE1)        :: zp_snowhmass
REAL, DIMENSION(KSIZE1)        :: zp_vegtype
REAL, DIMENSION(KSIZE1)        :: zp_forest
REAL, DIMENSION(KSIZE1)        :: zp_pew_a_coef
REAL, DIMENSION(KSIZE1)        :: zp_pew_b_coef
REAL, DIMENSION(KSIZE1)        :: zp_pet_a_coef
REAL, DIMENSION(KSIZE1)        :: zp_pet_b_coef
REAL, DIMENSION(KSIZE1)        :: zp_peq_a_coef
REAL, DIMENSION(KSIZE1)        :: zp_peq_b_coef
REAL, DIMENSION(KSIZE1)        :: zp_zenith
REAL, DIMENSION(KSIZE1)        :: zp_lat,zp_lon
REAL, DIMENSION(KSIZE1)        :: zp_psn_inv
REAL, DIMENSION(KSIZE1)        :: zp_psn
REAL, DIMENSION(KSIZE1)        :: zp_psn_gflxcor
REAL, DIMENSION(KSIZE1)        :: zp_work
!
REAL, PARAMETER :: zdepthabs = 0.60 ! m
!
INTEGER :: jwrk, jj, ji
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('SNOW3L_ISBA:CALL_MODEL',0,zhook_handle)
!
! Initialize:
!
zp_psn_gflxcor(:)  = 0.
zp_work(:)         = 0.
zp_soild(:)        = 0.
!
! pack the variables
!
DO jwrk=1,ksize2
   DO jj=1,ksize1
      ji = kmask(jj)
      zp_snowswe(jj,jwrk) = psnowswe(ji,jwrk)
      zp_snowrho(jj,jwrk) = psnowrho(ji,jwrk)
      zp_snowheat(jj,jwrk) = psnowheat(ji,jwrk)
      zp_snowtemp(jj,jwrk) = psnowtemp(ji,jwrk)
      zp_snowliq(jj,jwrk) = psnowliq(ji,jwrk)
      zp_snowdz(jj,jwrk) = psnowdz(ji,jwrk)
      zp_snowage(jj,jwrk) = psnowage(ji,jwrk)
   ENDDO
ENDDO
!
IF (hsnow_isba=='CRO') THEN
   DO jwrk=1,ksize2
      DO jj=1,ksize1
         ji = kmask(jj)
         zp_snowgran1(jj,jwrk) = psnowgran1(ji,jwrk)
         zp_snowgran2(jj,jwrk) = psnowgran2(ji,jwrk)
         zp_snowhist(jj,jwrk) = psnowhist(ji,jwrk)
      ENDDO
   ENDDO
ELSE
   DO jwrk=1,ksize2
      DO jj=1,ksize1
         zp_snowgran1(jj,jwrk) = xundef
         zp_snowgran2(jj,jwrk) = xundef
         zp_snowhist(jj,jwrk) = xundef
      ENDDO
   ENDDO
ENDIF
!  
DO jwrk=1,ksize3
   DO jj=1,ksize1
      ji                    = kmask(jj)
      zp_tg(jj,jwrk) = ptg(ji,jwrk)
      zp_d_g(jj,jwrk) = pd_g(ji,jwrk)
      zp_soilhcapz(jj,jwrk) = psoilhcapz(ji,jwrk)
   ENDDO
ENDDO
!
IF (omeb) THEN
  DO jwrk=1,ksize3
    DO jj=1,ksize1
      ji                    = kmask(jj)
      zp_dzg(jj,jwrk) = pdzg(ji,jwrk)
    ENDDO
  ENDDO
ENDIF
!
DO jj=1,ksize1
   ji = kmask(jj)
   zp_snowalb(jj) = psnowalb(ji)    
   zp_ps(jj) = pps(ji)
   zp_srsnow(jj) = psr(ji)
   zp_rrsnow(jj) = zrrsnow(ji)
   zp_psn3l(jj) = ppsn(ji)
   zp_ct(jj) = pct(ji)
   zp_ta(jj) = pta(ji)
   zp_delheatg(jj) = pdelheatg(ji)
   zp_delheatg_sfc(jj) = pdelheatg_sfc(ji)
   zp_sw_rad(jj) = psw_rad(ji)
   zp_qa(jj) = pqa(ji)
   zp_vmod(jj) = pvmod(ji)
   zp_lw_rad(jj) = plw_rad(ji)
   zp_rhoa(jj) = prhoa(ji)
   zp_uref(jj) = puref(ji)
   zp_exns(jj) = pexns(ji)
   zp_exna(jj) = pexna(ji)
   zp_lvtt(jj) = plvtt(ji)
   zp_lstt(jj) = plstt(ji)
   zp_dircoszw(jj) = pdircoszw(ji)
   zp_zref(jj) = pzref(ji)
   zp_z0nat(jj) = pz0nat(ji)
   zp_z0hnat(jj) = pz0hnat(ji)
   zp_z0eff(jj) = pz0eff(ji)
   zp_alb(jj) = palb(ji)
   zp_soilcond(jj) = zsoilcond(ji)
   !  
   zp_pew_a_coef(jj) = ppew_a_coef(ji)
   zp_pew_b_coef(jj) = ppew_b_coef(ji)
   zp_pet_a_coef(jj) = ppet_a_coef(ji)
   zp_peq_a_coef(jj) = ppeq_a_coef(ji)      
   zp_pet_b_coef(jj) = ppet_b_coef(ji)
   zp_peq_b_coef(jj) = ppeq_b_coef(ji)
   !
   zp_lat(jj)      = plat(ji)
   zp_lon(jj)      = plon(ji)
   zp_zenith(jj)     = pzenith(ji)
!
   zp_grndflux(jj) = pgrndflux(ji)
   zp_rnsnow(jj) = prnsnow(ji)
   zp_hsnow(jj) = phsnow(ji)
   zp_delheatn(jj) = pdelheatn(ji)
   zp_delheatn_sfc(jj) = pdelheatn_sfc(ji)
   zp_snowsfch(jj) = psnowsfch(ji)
   zp_hpsnow(jj) = phpsnow(ji)
   zp_les3l(jj) = ples3l(ji) 
   zp_lel3l(jj) = plel3l(ji)  
   zp_evap(jj) = pevap(ji)
   zp_emisnow(jj) = pemisnow(ji) 
   zp_swnetsnow(jj) = pswnetsnow(ji) 
   zp_swnetsnows(jj) = pswnetsnows(ji) 
   zp_lwnetsnow(jj) = plwnetsnow(ji) 
ENDDO
!
DO jj=1,ksize1
   ji = kmask(jj)
   zp_vegtype(jj) = pvegtype(ji,nvt_snow)
   zp_forest(jj) = pvegtype(ji,nvt_tebd) + pvegtype(ji,nvt_trbe) + pvegtype(ji,nvt_bone)   &
                   + pvegtype(ji,nvt_trbd) + pvegtype(ji,nvt_tebe) + pvegtype(ji,nvt_tene)   & 
                   + pvegtype(ji,nvt_bobd) + pvegtype(ji,nvt_bond) + pvegtype(ji,nvt_shrb)   
ENDDO
!
!
! ===============================================================
! conversion of snow heat from J/m3 into J/m2
WHERE(zp_snowswe(:,:)>0.) &
  zp_snowheat(:,:) = zp_snowheat(:,:) / zp_snowrho(:,:) * zp_snowswe(:,:)  
! ===============================================================
!
zp_psn_inv(:)       = 0.0
zp_psn(:)           = zp_psn3l(:)
!
IF(omeb)THEN
!
!   MEB (case of imposed surface fluxes)
!   - Prepare inputs for explicit snow scheme(s):
!     If using MEB, these are INPUTs ONLY:
!     divide fluxes by snow fraction to make "snow-relative"
!
   zp_psn(:)         = max(1.e-4, zp_psn3l(:))
   zp_psn_inv(:)     = 1.0/zp_psn(:)
!
   zp_rnsnow(:)      = zp_rnsnow(:)      *zp_psn_inv(:)
   zp_swnetsnow(:)   = zp_swnetsnow(:)   *zp_psn_inv(:)
   zp_swnetsnows(:)  = zp_swnetsnows(:)  *zp_psn_inv(:)
   zp_lwnetsnow(:)   = zp_lwnetsnow(:)   *zp_psn_inv(:)
   zp_hsnow(:)       = zp_hsnow(:)       *zp_psn_inv(:)
   zp_gfluxsnow(:)   = zp_gfluxsnow(:)   *zp_psn_inv(:) 
   zp_gsfcsnow(:)    = zp_gsfcsnow(:)    *zp_psn_inv(:) 
   zp_snowhmass(:)   = zp_snowhmass(:)   *zp_psn_inv(:)
   zp_les3l(:)       = zp_les3l(:)       *zp_psn_inv(:)
   zp_lel3l(:)       = zp_lel3l(:)       *zp_psn_inv(:)
   zp_grndflux(:)    = zp_grndflux(:)    *zp_psn_inv(:)
   zp_evap(:)        = zp_evap(:)        *zp_psn_inv(:)
   zp_hpsnow(:)      = zp_hpsnow(:)      *zp_psn_inv(:)
   zp_delheatn(:)    = zp_delheatn(:)    *zp_psn_inv(:)
   zp_delheatn_sfc(:)= zp_delheatn_sfc(:)*zp_psn_inv(:)
   zp_snowsfch(:)    = zp_snowsfch(:)    *zp_psn_inv(:)

   zp_srsnow(:)      = zp_srsnow(:)      *zp_psn_inv(:)
   zp_rrsnow(:)      = zp_rrsnow(:)      *zp_psn_inv(:)

   DO jj=1,ksize2
      DO ji=1,ksize1
         zp_snowswe(ji,jj)  = zp_snowswe(ji,jj) *zp_psn_inv(ji)
         zp_snowheat(ji,jj) = zp_snowheat(ji,jj)*zp_psn_inv(ji)
         zp_snowdz(ji,jj)   = zp_snowdz(ji,jj)  *zp_psn_inv(ji)
      ENDDO
   ENDDO
   !
ENDIF
!
! Call ISBA-SNOW3L model:  
!  
IF (hsnow_isba=='CRO') THEN 

   CALL snowcro(hsnowres, tptime, oglacier, himplicit_wind,                &
             zp_pew_a_coef, zp_pew_b_coef,                                 &
             zp_pet_a_coef, zp_peq_a_coef, zp_pet_b_coef, zp_peq_b_coef,   &
             zp_snowswe,zp_snowrho, zp_snowheat, zp_snowalb,               &
             zp_snowgran1, zp_snowgran2, zp_snowhist, zp_snowage, ptstep,  &
             zp_ps, zp_srsnow, zp_rrsnow ,zp_psn3l, zp_ta, zp_tg(:,1),     &
             zp_sw_rad, zp_qa, zp_vmod, zp_lw_rad, zp_rhoa, zp_uref,       &
             zp_exns, zp_exna, zp_dircoszw, zp_zref, zp_z0nat, zp_z0eff,   &
             zp_z0hnat, zp_alb, zp_soilcond, zp_d_g(:,1), zp_snowliq,      &
             zp_snowtemp, zp_snowdz, zp_thrufal, zp_grndflux, zp_evapcor,  &
             zp_rnsnow, zp_hsnow, zp_gfluxsnow, zp_hpsnow, zp_les3l,       &
             zp_lel3l, zp_evap, zp_sndrift, zp_ri,                         &
             zp_emisnow, zp_cdsnow, zp_ustarsnow,                          &
             zp_chsnow, zp_snowhmass, zp_qs, zp_vegtype, zp_zenith,        &
             zp_lat, zp_lon, osnowdrift,osnowdrift_sublim,                 &
             osnow_abs_zenith, hsnowmetamo,hsnowrad                        )
!
  zp_gflxcor(:) = 0.0
  zp_flsn_cor(:) = 0.0
  zp_soilcor(:) = 0.0
!
ELSE 
!
  CALL snow3l(hsnowres, tptime, omeb, himplicit_wind,                      &
             zp_pew_a_coef, zp_pew_b_coef,                                 &
             zp_pet_a_coef, zp_peq_a_coef,zp_pet_b_coef, zp_peq_b_coef,    &
             zp_snowswe, zp_snowrho, zp_snowheat, zp_snowalb,              &
             zp_snowgran1, zp_snowgran2, zp_snowhist, zp_snowage, ptstep,  &
             zp_ps, zp_srsnow, zp_rrsnow, zp_psn3l, zp_ta, zp_tg(:,1),     &
             zp_sw_rad, zp_qa, zp_vmod, zp_lw_rad, zp_rhoa, zp_uref,       &
             zp_exns, zp_exna, zp_dircoszw, zp_zref, zp_z0nat, zp_z0eff,   &
             zp_z0hnat, zp_alb, zp_soilcond, zp_d_g(:,1),                  &
             zp_lvtt, zp_lstt, zp_snowliq,                                 &
             zp_snowtemp, zp_snowdz, zp_thrufal, zp_grndflux ,             &
             zp_evapcor, zp_soilcor, zp_gflxcor, zp_snowsfch,              &
             zp_delheatn, zp_delheatn_sfc,                                 &
             zp_swnetsnow, zp_swnetsnows, zp_lwnetsnow, zp_gsfcsnow,       &
             zp_rnsnow, zp_hsnow, zp_gfluxsnow, zp_hpsnow, zp_les3l,       &
             zp_lel3l, zp_evap, zp_sndrift, zp_ri,                         &
             zp_emisnow, zp_cdsnow, zp_ustarsnow,                          &
             zp_chsnow, zp_snowhmass, zp_qs, zp_vegtype, zp_forest,        &
             zp_zenith, zp_lat, zp_lon, osnowdrift, osnowdrift_sublim     )
!
  IF(omeb)THEN
!
! - reverse transform: back to surface-relative
!
     zp_rnsnow(:)      = zp_rnsnow(:)      /zp_psn_inv(:)
     zp_swnetsnow(:)   = zp_swnetsnow(:)   /zp_psn_inv(:)
     zp_swnetsnows(:)  = zp_swnetsnows(:)  /zp_psn_inv(:)
     zp_lwnetsnow(:)   = zp_lwnetsnow(:)   /zp_psn_inv(:)
     zp_hsnow(:)       = zp_hsnow(:)       /zp_psn_inv(:)
     zp_les3l(:)       = zp_les3l(:)       /zp_psn_inv(:)
     zp_lel3l(:)       = zp_lel3l(:)       /zp_psn_inv(:)
     zp_grndflux(:)    = zp_grndflux(:)    /zp_psn_inv(:)
     zp_evap(:)        = zp_evap(:)        /zp_psn_inv(:)
     zp_hpsnow(:)      = zp_hpsnow(:)      /zp_psn_inv(:)
     zp_gfluxsnow(:)   = zp_gfluxsnow(:)   /zp_psn_inv(:) 
     zp_delheatn(:)    = zp_delheatn(:)    /zp_psn_inv(:) 
     zp_delheatn_sfc(:)= zp_delheatn_sfc(:)/zp_psn_inv(:) 
     zp_snowsfch(:)    = zp_snowsfch(:)    /zp_psn_inv(:) 
     zp_gsfcsnow(:)    = zp_gsfcsnow(:)    /zp_psn_inv(:) 

     zp_srsnow(:)      = zp_srsnow(:)      /zp_psn_inv(:)
     zp_rrsnow(:)      = zp_rrsnow(:)      /zp_psn_inv(:)
     DO jj=1,ksize2
        DO ji=1,ksize1
           zp_snowswe(ji,jj)  = zp_snowswe(ji,jj) /zp_psn_inv(ji)
           zp_snowheat(ji,jj) = zp_snowheat(ji,jj)/zp_psn_inv(ji)
           zp_snowdz(ji,jj)   = zp_snowdz(ji,jj)  /zp_psn_inv(ji)
        ENDDO
     ENDDO
     
     zp_snowhmass(:)  = zp_snowhmass(:)/zp_psn_inv(:)
     zp_thrufal(:)    = zp_thrufal(:)  /zp_psn_inv(:)
!
!    Final Adjustments:
!    ------------------
!    Add cooling/heating flux correction to underlying soil.
!    This term is usually active for vanishingly thin snowpacks..
!    it is put outside of the snow scheme owing to it's dependence on
!    snow fraction. It is related to a possible correction to the ground-snow
!    heat flux when it is imposed (using MEB).
!    Also, it is added as a heat sink/source here since
!    fluxes have already be computed and should not be adjusted at this point:
!    applying it to the soil has the same impact as soil freeze-thaw, in the
!    sense it is computed after the fluxes have been updated.
!    (and update heat storage diagnostic in a consistent manner)
!
!    Energy is thickness weighted, thus thicker layers receive more energy and energy
!    is evenly distributed to depth ZDEPTHABS. An
!    alternate method is to weight near surface layers more and diminish weights
!    (thus eenrgy received by each layer) with depth. Both methods conserve energy as 
!    long as vertical weights are normalized.

!    i) Determine soil depth for energy absorption:

     zp_soild(:) = zp_dzg(:,1)
     DO jj=2,ksize3
        DO ji=1,ksize1
           IF(zp_dzg(ji,jj) <= zdepthabs)THEN
              zp_soild(ji) = zp_dzg(ji,jj)
           ENDIF
        ENDDO
     ENDDO

!    ii) Distribute (possible) energy to absorb vertically over some layer (defined above):

     zp_psn_gflxcor(:)  = zp_psn(:)*zp_gflxcor(:)                                ! (W/m2)
     zp_work(:)         = zp_psn_gflxcor(:)*ptstep/zp_soild(:)

     zp_tg(:,1)         = zp_tg(:,1)         + zp_work(:)*zp_ct(:)*zp_d_g(:,1)   ! (K)
     DO jj=2,ksize3
        DO ji=1,ksize1
           IF(zp_soild(ji) <= zdepthabs)THEN
              zp_tg(ji,jj) = zp_tg(ji,jj)    + zp_work(ji)/zp_soilhcapz(ji,jj)   ! K
           ENDIF
        ENDDO
     ENDDO

     zp_delheatg(:)     = zp_delheatg(:)     + zp_psn_gflxcor(:)                 ! (W/m2)
     zp_delheatg_sfc(:) = zp_delheatg_sfc(:) + zp_psn_gflxcor(:)                 ! (W/m2)
!
     zp_flsn_cor(:)     = 0.0
!
  ELSE
!
!    To conserve energy in ISBA, the correction flux must be distributed at least
!    over the first 60cm depth. This method prevent numerical oscillations
!    especially when explicit snow vanishes. Final Adjustments are done in ISBA_CEB
!
     zp_flsn_cor(:) = zp_gflxcor(:) ! (W/m2)
!
  ENDIF
!
ENDIF
!
!===============================================================
!conversion of snow heat from J/m2 into J/m3
WHERE(zp_snowswe(:,:)>0.)
      zp_snowheat(:,:)=zp_snowheat(:,:)*zp_snowrho(:,:)/zp_snowswe(:,:)  
ENDWHERE
!===============================================================
!
! === Packing:
!
! unpack variables
!
DO jwrk=1,ksize2
  DO jj=1,ksize1
    ji = kmask(jj)
    psnowswe(ji,jwrk) = zp_snowswe(jj,jwrk)
    psnowrho(ji,jwrk) = zp_snowrho(jj,jwrk)
    psnowheat(ji,jwrk) = zp_snowheat(jj,jwrk)
    psnowtemp(ji,jwrk) = zp_snowtemp(jj,jwrk)
    psnowliq(ji,jwrk) = zp_snowliq(jj,jwrk)
    psnowdz(ji,jwrk) = zp_snowdz(jj,jwrk)
    psnowage(ji,jwrk) = zp_snowage(jj,jwrk)
  ENDDO
ENDDO
!
IF (hsnow_isba=='CRO') THEN
  DO jwrk=1,ksize2
    DO jj=1,ksize1
      ji = kmask(jj)
      psnowgran1(ji,jwrk) = zp_snowgran1(jj,jwrk)
      psnowgran2(ji,jwrk) = zp_snowgran2(jj,jwrk)
      psnowhist(ji,jwrk) = zp_snowhist(jj,jwrk)
    ENDDO
  ENDDO
ENDIF
!
DO jwrk=1,ksize3
   DO jj=1,ksize1
      ji              = kmask(jj)
      ptg(ji,jwrk)= zp_tg(jj,jwrk)
   ENDDO
ENDDO
!
DO jj=1,ksize1
  ji                  = kmask(jj)
  pdelheatg(ji)   = zp_delheatg(jj)
  pdelheatg_sfc(ji)   = zp_delheatg_sfc(jj)
  psnowalb(ji)   = zp_snowalb(jj)
  pthrufal(ji)   = zp_thrufal(jj)
  pevapcor(ji)   = zp_evapcor(jj)
  zsoilcor(ji)   = zp_soilcor(jj)
  pri(ji)   = zp_ri(jj)
  pqs(ji)   = zp_qs(jj)
  pcdsnow(ji)   = zp_cdsnow(jj)
  pustarsnow(ji)   = zp_ustarsnow(jj)
  pchsnow(ji)   = zp_chsnow(jj)
  psnowhmass(ji)   = zp_snowhmass(jj)
  pgrndflux(ji)   = zp_grndflux(jj)
  pflsn_cor(ji)   = zp_flsn_cor(jj)
  prnsnow(ji)   = zp_rnsnow(jj)
  phsnow(ji)   = zp_hsnow(jj)
  pgfluxsnow(ji)   = zp_gfluxsnow(jj)
  pdelheatn(ji)   = zp_delheatn(jj)
  pdelheatn_sfc(ji)   = zp_delheatn_sfc(jj)
  psnowsfch(ji)   = zp_snowsfch(jj)
  pgsfcsnow(ji)   = zp_gsfcsnow(jj)
  phpsnow(ji)   = zp_hpsnow(jj)
  ples3l(ji)   = zp_les3l(jj)
  plel3l(ji)   = zp_lel3l(jj)
  pevap(ji)   = zp_evap(jj)
  pemisnow(ji)   = zp_emisnow(jj)
  pswnetsnow(ji)   = zp_swnetsnow(jj)
  pswnetsnows(ji)   = zp_swnetsnows(jj)
  plwnetsnow(ji)   = zp_lwnetsnow(jj)
ENDDO
!
IF (lhook) CALL dr_hook('SNOW3L_ISBA:CALL_MODEL',1,zhook_handle)
!
END SUBROUTINE call_model
!
END SUBROUTINE snow3l_isba