snow_cover_1layer.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 snow_cover_1layer(PTSTEP, PANSMIN, PANSMAX, PTODRY, PRHOSMIN, PRHOSMAX,   &
                                 PRHOFOLD, OALL_MELT, PDRAIN_TIME, PWCRN, PZ0SN, PZ0HSN, &
                                 TPSNOW, PTG, PTG_COEFA, PTG_COEFB, PABS_SW, PLW1, PLW2, &
                                 PTA, PQA, PVMOD, PPS, PRHOA, PSR, PZREF, PUREF, PRNSNOW,&
                                 PHSNOW, PLESNOW, PGSNOW, PMELT, PDQS_SNOW, PABS_LW  )  
!   ##########################################################################
!
!!****  *SNOW_COVER_1LAYER*  
!!
!!    PURPOSE
!!    -------
!
!     One layer snow mantel scheme
!         
!     
!!**  METHOD
!     ------
!
!
! The temperature equation is written as:
!
!              b T+ = y
!
!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    MODD_CST
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    08/09/98 
!!      J. Escobar 24/10/2012 : BUF PGI10.X , rewrite some 1 line WHERE statement
!!      V. Masson  13/09/2013 : implicitation of coupling with roof below
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_type_snow, ONLY : surf_snow
!
USE modd_csts,       ONLY : xtt, xci, xrholi, xrholw, xcpd, xlstt, xlmtt, xday, xcondi
USE modd_snow_par,   ONLY : xemissn
USE modd_surf_par,   ONLY : xundef
!
USE mode_thermos
!
USE modi_surface_ri
USE modi_surface_aero_cond
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
REAL,                 INTENT(IN)    :: PTSTEP   ! time step
REAL,                 INTENT(IN)    :: PANSMIN  ! minimum snow albedo
REAL,                 INTENT(IN)    :: PANSMAX  ! maximum snow albedo
REAL,                 INTENT(IN)    :: PTODRY   ! snow albedo decreasing constant
REAL,                 INTENT(IN)    :: PRHOSMIN ! minimum snow density
REAL,                 INTENT(IN)    :: PRHOSMAX ! maximum snow density
REAL,                 INTENT(IN)    :: PRHOFOLD ! snow density increasing constant
LOGICAL,              INTENT(IN)    :: OALL_MELT! T --> all snow runs off if
                                                ! lower surf. temperature is
                                                ! positive
REAL,                 INTENT(IN)    :: PDRAIN_TIME ! drainage folding time (days)
REAL,                 INTENT(IN)    :: PWCRN    ! critical snow amount necessary
                                                ! to cover the considered surface
REAL,                 INTENT(IN)    :: PZ0SN    ! snow roughness length for momentum
REAL,                 INTENT(IN)    :: PZ0HSN   ! snow roughness length for heat
TYPE(surf_snow), INTENT(INOUT) :: TPSNOW
REAL, DIMENSION(:), INTENT(IN)    :: PTG      ! underlying ground temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTG_COEFA! underlying ground temperature
REAL, DIMENSION(:), INTENT(IN)    :: PTG_COEFB! implicit terms
REAL, DIMENSION(:), INTENT(IN)    :: PABS_SW  ! absorbed SW energy (Wm-2)
REAL, DIMENSION(:), INTENT(IN)    :: PLW1     ! LW coef independant of TSNOW
                                              ! (Wm-2)     usually equal to:
                                              !      emis_snow * LW_down
                                              !
REAL, DIMENSION(:), INTENT(IN)    :: PLW2     ! LW coef dependant   of TSNOW
                                              ! (Wm-2 K-4) usually equal to:
                                              ! -1 * emis_snow * stefan_constant
                                              !
REAL, DIMENSION(:), INTENT(IN)    :: PTA      ! temperature at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PQA      ! specific humidity
                                                ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PVMOD    ! module of the horizontal wind
REAL, DIMENSION(:), INTENT(IN)    :: PPS      ! pressure at the surface
REAL, DIMENSION(:), INTENT(IN)    :: PRHOA    ! air density
                                                ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: PSR      ! snow rate
REAL, DIMENSION(:), INTENT(IN)    :: PZREF    ! reference height of the first
                                              ! atmospheric level (temperature)
REAL, DIMENSION(:), INTENT(IN)    :: PUREF    ! reference height of the first
                                              ! atmospheric level (wind)
REAL, DIMENSION(:), INTENT(OUT)   :: PRNSNOW  ! net radiation over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PHSNOW   ! sensible heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PLESNOW  ! latent heat flux over snow
REAL, DIMENSION(:), INTENT(OUT)   :: PGSNOW   ! flux under the snow
REAL, DIMENSION(:), INTENT(OUT)   :: PMELT    ! snow melting rate (kg/m2/s)
REAL, DIMENSION(:), INTENT(OUT)   :: PDQS_SNOW! heat storage inside snow
REAL, DIMENSION(:), INTENT(OUT)   :: PABS_LW ! absorbed LW rad by snow (W/m2)
!
!
!*      0.2    declarations of local variables
!
REAL :: ZEXPL = 0.
REAL :: ZIMPL = 1.
!
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZEXNS, ZEXNA, ZDIRCOSZW
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZZ0      ! roughness length for momentum
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZZ0H     ! roughness length forheat
!
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZRI      ! Richardson number
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZAC      ! aerodynamical conductance
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZRA      ! aerodynamical resistance
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZCH      ! drag coefficient for heat
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZB, ZY   ! coefficients in Ts eq.
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZWSNOW   ! snow before evolution
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZSNOW_HC ! snow heat capacity
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZSNOW_TC ! snow thermal conductivity
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZSNOW_D  ! snow depth
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZMELT    ! snow melting rate (kg/m3/s)
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZTS_SNOW ! snow surface temperature
                                          ! at previous time-step
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZQSAT    ! specific humidity
!                                         ! for ice
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZDQSAT   ! d(specific humidity)/dT
!                                         ! for ice
!
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZSR1, ZSR2   ! norm. snow precip.
!
LOGICAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: GSNOWMASK ! where snow is
!                                             ! at previuos time-step
LOGICAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: GFLUXMASK ! where fluxes can
!                                             ! be computed at
!                                             ! new time-step
!                                             ! i.e. snow occurence
!                                             ! at previous time-step
!                                             ! OR snow fall
INTEGER, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: JSNOWMASK1, JSNOWMASK2, JSNOWMASK3 ! where snow is or not
!                                                                      ! at previuos time-step
INTEGER, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: JFLUXMASK ! where fluxes can
!                                             ! be computed at
!                                             ! new time-step
!                                             ! i.e. snow occurence
!                                             ! at previous time-step
!                                             ! OR snow fall
!
REAL :: ZWSNOW_MIN = 0.1 ! minimum value of snow content (kg/m2) for prognostic
!                        ! computations
!
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZEI_SNOW  ! internal energy of snow
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZPEI_SNOW ! internal energy of snow at t+
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZWORK1
REAL, DIMENSION(SIZE(TPSNOW%WSNOW,1)) :: ZDQSATI, ZQSATI
!
INTEGER                         :: JJ, JI, JCOMPT_SNOW1, JCOMPT_SNOW2, JCOMPT_SNOW3, JCOMPT_FLUX
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
!
IF (lhook) CALL dr_hook('SNOW_COVER_1LAYER',0,zhook_handle)
zb(:)=0.
zy(:)=0.
zmelt(:) = 0.
pmelt(:) = 0.
prnsnow(:) = 0.
phsnow(:) = 0.
plesnow(:) = 0.
pgsnow(:) = 0.
!RJ: workaround to prevent decomposition unstable xundef masks for Tx_LWA_SN_RD fields
!RJ: in TEB_DIAGNOSTICS.nc during TEB_GARDEN_GREENROOF_BEM_3L_IRRIG_* tests
!RJ: problem with decomposition handling somewhere else
#ifdef RJ_PFIX
pabs_lw(:) = 0.0
#else
pabs_lw(:) = xundef
#endif
!
!* snow reservoir before evolution
!
zwsnow(:) = tpsnow%WSNOW(:,1)
zts_snow(:) = min(xtt,ptg(:))
!
zsnow_d(:) = 0.
zsnow_tc(:) = 0.
zsnow_hc(:) = 0.
!
!-------------------------------------------------------------------------------
!
!*      1.1    most useful masks
!              -----------------
!
gsnowmask(:)=.false.
gfluxmask(:)=.false.
jsnowmask1(:)=0.
jsnowmask2(:)=0.
jsnowmask3(:)=0.
jfluxmask(:)=0.

  !*      1.2    drag
!              ----
!
!*      1.2.1  defaults
!              --------
!
!* variation of temperature with altitude is neglected
!
zexns(:) = 1.
zexna(:) = 1.
!
!* slope is neglected in drag computation
!
zdircoszw(:) = 1.
!
!* roughness length are imposed:
!
zz0(:) = pz0sn
zz0h(:) = pz0hsn

!
!
!*      1.1    most useful masks
!              -----------------
!* snow occurence at previous time-step
!
!* snow occurence during the time-step for fluxes computation
!
jcompt_snow1=0
jcompt_snow2=0
jcompt_snow3=0
jcompt_flux=0
DO jj=1,SIZE(zwsnow)
  IF (zwsnow(jj)>0.) THEN
    gsnowmask(jj)=.true.
    !* surface temperature
    zts_snow(jj) = tpsnow%TS(jj)
    gfluxmask(jj)=.true.
    !gsnowmask=t
    jcompt_snow1=jcompt_snow1+1
    jsnowmask1(jcompt_snow1) = jj
    !gfluxmask=t
    jcompt_flux=jcompt_flux+1
    jfluxmask(jcompt_flux) = jj
    IF (zwsnow(jj)>=zwsnow_min) THEN
      !second snow mask
      jcompt_snow3=jcompt_snow3+1
      jsnowmask3(jcompt_snow3)=jj
    ELSE
      !lower limit of snow cover for prognostic computations
      !0.1 kg/m2 of snow water content
      tpsnow%T(jj,1) = min(ptg(jj),xtt)
    ENDIF
  ELSE
    tpsnow%T(jj,1) = min(ptg(jj),xtt)
    !gsnowmask=false
    jcompt_snow2=jcompt_snow2+1
    jsnowmask2(jcompt_snow2) = jj
    IF (psr(jj)>0.) THEN
      gfluxmask(jj)=.true.
      jcompt_flux=jcompt_flux+1
      jfluxmask(jcompt_flux) = jj
    ENDIF
  ENDIF
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      1.2    drag
!              ----
!
!*      1.2.2   qsat (Tsnow)
!              ------------
!
zqsat(:) = qsati(zts_snow(:), pps(:) )
!
!*      1.2.3  Richardson number
!              -----------------
!
!* snow is present on all the considered surface.
!* computation occurs where snow is and/or falls.
!
 CALL surface_ri(zts_snow, zqsat, zexns, zexna, pta, pqa, &
                  pzref, puref, zdircoszw, pvmod, zri      )  
!
!*      1.2.4  Aerodynamical conductance
!              -------------------------
!
 CALL surface_aero_cond(zri, pzref, puref, pvmod, zz0, zz0h, zac, zra, zch)
!
!-------------------------------------------------------------------------------
!
!*      2.     snow thermal characteristics
!              ----------------------------
!cdir nodep
DO jj=1,jcompt_snow1
  !
  ji = jsnowmask1(jj)
  !
  !*      2.1    snow heat capacity
  zsnow_hc(ji) = tpsnow%RHO(ji,1) * xci * xrholi / xrholw
!*      2.2    snow depth
  zsnow_d(ji) = zwsnow(ji) / tpsnow%RHO(ji,1)
!*      2.3    snow thermal conductivity
  zsnow_tc(ji) = xcondi * (tpsnow%RHO(ji,1)/xrholw)**1.885
!*      2.4    internal energy of snow
  zei_snow(ji) = zsnow_hc(ji)*zsnow_d(ji)*tpsnow%T(ji,1)
  !
ENDDO
!
!cdir nodep
DO jj=1,jcompt_snow2
  !
  ji = jsnowmask2(jj)
  !
  !*      2.1    snow heat capacity
  zsnow_hc(ji) = prhosmin * xci * xrholi / xrholw
!*      2.2    snow depth
  zsnow_d(ji) = ptstep * psr(ji) / prhosmin
!*      2.3    snow thermal conductivity
  zsnow_tc(ji) = xcondi * (prhosmin /xrholw)**1.885
!*      2.4    internal energy of snow
  zei_snow(ji) = 0.
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      3.     Snow temperature evolution
!              --------------------------
!
!*      3.5    dqsat/ dT (Tsnow)
!              -----------------
!
zdqsati = dqsati(zts_snow(:),pps(:),zqsat(:))
WHERE (gsnowmask(:) .AND. zwsnow(:)>=zwsnow_min)
  zdqsat(:) = zdqsati(:)
END WHERE
!
!*      3.1    coefficients from Temperature tendency
!              --------------------------------------
!
!cdir nodep
DO jj=1,jcompt_snow3
!
  ji=jsnowmask3(jj)

  zwork1(ji) = zsnow_d(ji) * zsnow_hc(ji) / ptstep
!
  zb(ji) = zb(ji) + zwork1(ji)
!
!*      3.2    coefficients from solar radiation
!          ---------------------------------
!    
  zy(ji) = zy(ji) + zwork1(ji) * tpsnow%T(ji,1) + pabs_sw(ji)
!
!
!*      3.3    coefficients from infra-red radiation
!              -------------------------------------
!
  zwork1(ji) = plw2(ji) * tpsnow%T(ji,1)**3
!
  zb(ji) = zb(ji) - 4 * zimpl * zwork1(ji)
!
  zy(ji) = zy(ji) + plw1(ji) + zwork1(ji) * (zexpl-3.*zimpl) * tpsnow%T(ji,1)
!
!
!*      3.4    coefficients from sensible heat flux
!              ------------------------------------
!
 zwork1(ji) = xcpd * prhoa(ji) * zac(ji)
! 
  zb(ji) = zb(ji) + zwork1(ji) *   zimpl
!
  zy(ji) = zy(ji) - zwork1(ji) * ( zexpl * tpsnow%T(ji,1) - pta(ji) )
!
!
!*      3.6    coefficients from latent heat flux
!              ----------------------------------
!
  zwork1(ji) =  xlstt * prhoa(ji) * zac(ji)
!
  zb(ji) = zb(ji) + zwork1(ji) *  zimpl * zdqsat(ji)
!
  zy(ji) = zy(ji) - zwork1(ji) * (  zqsat(ji) - pqa(ji) - zimpl * zdqsat(ji)*tpsnow%T(ji,1) )
!
!*      3.7    coefficients from conduction flux at snow base
!              ----------------------------------------------
!
  zwork1(ji) = zsnow_tc(ji)/(0.5*zsnow_d(ji))
!
  zb(ji) = zb(ji) + zwork1(ji) *  zimpl / ( 1. + zwork1(ji)*ptg_coefa(ji) )
!
  zy(ji) = zy(ji) - zwork1(ji) * (zexpl * tpsnow%T(ji,1) - ptg_coefb(ji)) &
                   / ( 1. + zwork1(ji)*ptg_coefa(ji) )
!
!*      3.8    guess of snow temperature before accumulation and melting
!              ---------------------------------------------------------
!
  tpsnow%T(ji,1) = zy(ji) / zb(ji)
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      4.     Snow melt
!              ---------
!
!*      4.1    melting
!              -------
!
!cdir nodep
DO jj=1,jcompt_snow1
!
  ji = jsnowmask1(jj)
!
  zmelt(ji)  = max( tpsnow%T(ji,1) - xtt , 0. ) * zsnow_hc(ji) /  xlmtt / ptstep
!
  zmelt(ji)  = min( zmelt(ji) , zwsnow(ji) / zsnow_d(ji) / ptstep )
!
  tpsnow%T(ji,1) = min( tpsnow%T(ji,1) , xtt )
!
ENDDO
!
!*      4.2    run-off of all snow if lower surface temperature is positive
!              ------------------------------------------------------------
!
!* this option is used when snow is located on sloping roofs for example.
!
IF (oall_melt) THEN
  WHERE ( gsnowmask(:) .AND. ptg(:)>xtt .AND. zwsnow(:)>=zwsnow_min )
    pmelt(:) = pmelt(:) + zwsnow(:) / ptstep
  END WHERE
END IF
!
!*      4.3    output melting in kg/m2/s
!              -------------------------
!
pmelt(:) = zmelt(:) * zsnow_d(:)
!
!-------------------------------------------------------------------------------
!
!*      5.     fluxes
!              ------
!
!*      5.3    qsat (Tsnow)
!              ------------
!
zqsati = qsati(tpsnow%T(:,1),pps(:))
WHERE (gfluxmask(:)) 
   zqsat(:) = zqsati(:)
END WHERE
!
!*      5.1    net radiation (with Ts lin. extrapolation)
!              -------------
!
!cdir nodep
DO jj = 1, jcompt_flux
!
  ji = jfluxmask(jj)
!
  pabs_lw(ji) =  plw1(ji) + plw2(ji) * tpsnow%T(ji,1)**4
!
  prnsnow(ji) = pabs_sw(ji) + pabs_lw(ji)
!
!
!*      5.2    sensible heat flux
!              ------------------
!
  phsnow(ji) = xcpd * prhoa(ji) * zac(ji) * ( tpsnow%T(ji,1) - pta(ji) )
!
!
!*      5.4    latent heat flux
!              ----------------
!
  plesnow(ji) = xlstt * prhoa(ji) * zac(ji) * ( zqsat(ji) - pqa(ji) )
  !
!
!*      5.5    Conduction heat flux
!              --------------------
!
  !PGSNOW(JI) = ZSNOW_TC(JI)/(0.5*ZSNOW_D(JI)) * ( TPSNOW%T(JI,1) - PTG(JI) )
  pgsnow(ji) = zsnow_tc(ji)/(0.5*zsnow_d(ji)) * ( tpsnow%T(ji,1) - ptg_coefb(ji) ) &
             / ( 1. + zsnow_tc(ji)/(0.5*zsnow_d(ji))*ptg_coefa(ji) )
!
!
!*      5.6    If ground T>0 C, Melting is estimated from conduction heat flux
!              ---------------------------------------------------------------
!
  IF (ptg(ji)>xtt)  pmelt(ji) = max(pmelt(ji), -pgsnow(ji)/xlmtt)
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      6.     reservoir evolution
!              -------------------
!
!cdir nodep
DO jj = 1, SIZE(tpsnow%WSNOW,1)
!
!*      6.1    snow fall
!              ---------
!
  tpsnow%WSNOW(jj,1) = tpsnow%WSNOW(jj,1) + ptstep * psr(jj)
!
!
!*      6.2    sublimation
!              -----------
!
  plesnow(jj) = min( plesnow(jj), xlstt*tpsnow%WSNOW(jj,1)/ptstep )
!
  tpsnow%WSNOW(jj,1)  = max( tpsnow%WSNOW(jj,1) - ptstep * plesnow(jj)/xlstt , 0.)
!
  IF ( tpsnow%WSNOW(jj,1)<1.e-8 * ptstep ) tpsnow%WSNOW(jj,1) = 0.
!
!*      6.3    melting
!              -------
!
  pmelt(jj) = min( pmelt(jj), tpsnow%WSNOW(jj,1)/ptstep )
!
  tpsnow%WSNOW(jj,1)= max( tpsnow%WSNOW(jj,1) - ptstep * pmelt(jj) , 0.)
!
  IF ( tpsnow%WSNOW(jj,1)<1.e-8 * ptstep ) tpsnow%WSNOW(jj,1) = 0.
!
  IF (tpsnow%WSNOW(jj,1)==0.) pgsnow(jj) = max( pgsnow(jj), - pmelt(jj)*xlmtt )
!
ENDDO
!
!*      6.4    time dependent drainage
!              -----------------------
!
IF (pdrain_time>0.) THEN
  WHERE ( tpsnow%WSNOW(:,1)>0.)
    tpsnow%WSNOW(:,1) = tpsnow%WSNOW(:,1) * exp(-ptstep/pdrain_time/xday)
  END WHERE
END IF
!
!*      6.5    melting of last 1mm of snow depth
!              ---------------------------------
!
WHERE ( tpsnow%WSNOW(:,1)<zwsnow_min .AND. pmelt(:)>0. .AND. psr(:)==0. )
  pmelt(:) = pmelt(:) + tpsnow%WSNOW(:,1) / ptstep
  tpsnow%WSNOW(:,1)=0.
END WHERE
!
WHERE ( tpsnow%WSNOW(:,1)<1.e-8 * ptstep ) 
   tpsnow%WSNOW(:,1) = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*      7.     albedo evolution
!              ----------------
!
!*      7.1    If melting occurs or not
!              -----------------------
!
!
!cdir nodep
DO jj=1,jcompt_snow1
!
  ji = jsnowmask1(jj)
!
  IF (pmelt(ji) > 0. ) THEN
!
    tpsnow%ALB(ji) = (tpsnow%ALB(ji)-pansmin)*exp(-prhofold*ptstep/xday) + pansmin   &
                     + psr(ji)*ptstep/pwcrn*pansmax  
!
  ELSEIF (pmelt(ji)==0.) THEN 
    tpsnow%ALB(ji) = tpsnow%ALB(ji) - ptodry*ptstep/xday  + psr(ji)*ptstep/pwcrn*pansmax  
!
  ENDIF
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      8.     density evolution
!              -----------------
!
!*      8.1    old snow
!              --------
!
!cdir nodep
DO jj = 1, jcompt_snow1
! 
  ji = jsnowmask1(jj)
!
  IF (tpsnow%WSNOW(ji,1)>0. ) THEN
!
    zsr1(ji) = max( tpsnow%WSNOW(ji,1) , psr(ji) * ptstep )
!
    tpsnow%RHO(ji,1) = (tpsnow%RHO(ji,1)-prhosmax)*exp(-prhofold*ptstep/xday) + prhosmax
    tpsnow%RHO(ji,1) = ( (zsr1(ji)-psr(ji)*ptstep) * tpsnow%RHO(ji,1)    &
                  + (psr(ji)*ptstep) * prhosmin ) / zsr1(ji)  
  ENDIF
!
ENDDO
!
!*      8.2    fresh snow
!              ----------
!
!cdir nodep
DO jj=1,SIZE(tpsnow%WSNOW,1)
  IF (  tpsnow%WSNOW(jj,1)>0. ) THEN
    tpsnow%ALB(jj) = max(tpsnow%ALB(jj),pansmin)
    tpsnow%ALB(jj) = min(tpsnow%ALB(jj),pansmax)
    IF (zwsnow(jj)==0.) THEN
      tpsnow%ALB  (jj) = pansmax
      tpsnow%EMIS (jj) = xemissn
      tpsnow%RHO(jj,1) = prhosmin
    ENDIF
  ENDIF
    ENDDO
!
!-------------------------------------------------------------------------------
!
!*      9.     fresh snow accumulation (if more than 1mm of snow depth)
!              -----------------------
!
!cdir nodep
DO jj=1,jcompt_snow3
!
  ji = jsnowmask3(jj)
!
  IF (psr(ji)>0. .AND. tpsnow%WSNOW(ji,1)>0.) THEN
!
    zsr2(ji) = min( tpsnow%WSNOW(ji,1) , psr(ji) * ptstep )
!
    tpsnow%T(ji,1) =( ( tpsnow%WSNOW(ji,1) - zsr2(ji) ) *  tpsnow%T(ji,1)        &
              +   zsr2(ji)  * min( pta(ji) ,xtt )) / ( tpsnow%WSNOW(ji,1) )  
  ENDIF
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*     10.     Surface temperature
!              -------------------
!
!* note that if the relation between snow pack temperature and its
!  surface temperature is modified, think to modify it also in
!  subroutine init_snow_lw.f90
!
WHERE (gsnowmask(:) )
  tpsnow%TS(:) = tpsnow%T(:,1)
END WHERE
!
!-------------------------------------------------------------------------------
!
!*     11.     bogus values
!              ------------
!
!*     11.1    snow characteristics where snow IS present at current time-step
!              ---------------------------------------------------------------
!
WHERE (tpsnow%WSNOW(:,1)==0.)
  tpsnow%T(:,1)   = xundef
  tpsnow%RHO(:,1) = xundef
  tpsnow%ALB(:)   = xundef
  tpsnow%TS(:)    = xundef
  tpsnow%EMIS(:)  = xundef
END WHERE
!
!
!-------------------------------------------------------------------------------
!
!*     12.     Heat storage inside snow pack
!
WHERE (gsnowmask(:))
  zpei_snow(:) = zsnow_hc(:)*zsnow_d(:)*tpsnow%T(:,1)
ELSEWHERE
  zpei_snow(:) = 0.
END WHERE
pdqs_snow(:) = (zpei_snow(:)-zei_snow(:))/ptstep
!
IF (lhook) CALL dr_hook('SNOW_COVER_1LAYER',1,zhook_handle)

!------------------------------------------------------------------------------- !
END SUBROUTINE snow_cover_1layer