prep_perm_snow.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 prep_perm_snow (I, &
                           tpsnow,ptg,pperm_snow_frac,ksnow)
!          ################################################
!
!
!!****  *PREP_PERM_SNOW* - takes into account permanent snow into prognostic snow
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      B. Decharme 03/2009: Consistency with Arpege permanent
!!                                          snow/ice treatment
!!      B. Decharme 07/2012: 3-L or Crocus adjustments
!!      M. Lafaysse 09/2012: adaptation with new snow age in Crocus
!!------------------------------------------------------------------
!

!
USE modd_isba_n, ONLY : isba_t
!
USE modd_type_snow
USE modd_csts,           ONLY : xtt
USE modd_data_cover_par, ONLY : nvt_snow
USE modd_snow_par,       ONLY : xrhosmax, xansmax, xansmin, &
                                xaglamax, xaglamin, xhgla,  &
                                xrhosmax_es
USE modd_surf_par,       ONLY : xundef
!
USE modd_isba_par,       ONLY : xwgmin
!
USE modi_snow_heat_to_t_wliq
USE modi_snow_t_wliq_to_heat
USE modi_mkflag_snow
USE mode_surf_snow_frac
USE mode_snow3l
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(isba_t), INTENT(INOUT) :: i
!
TYPE(surf_snow), INTENT(INOUT) :: tpsnow            ! snow mantel characteristics
REAL, DIMENSION(:,:),  INTENT(IN):: ptg             ! soil temperature for patch KSNOW
REAL, DIMENSION(:,:),  INTENT(IN):: pperm_snow_frac ! fraction of permanent snow for patch KSNOW
INTEGER,               INTENT(IN):: ksnow           ! patch number where permanent snow is
!
!*      0.2    declarations of local parameter
!
REAL, PARAMETER :: zrhol1 = 150.
!
!*      0.3    declarations of local variables
!
INTEGER                             :: jlayer      ! loop counter on snow layers
REAL, DIMENSION(:),   ALLOCATABLE   :: zwsnow_perm ! snow total reservoir due to perm. snow
REAL, DIMENSION(:),   ALLOCATABLE   :: zwsnow      ! initial snow total reservoir
REAL, DIMENSION(:),   ALLOCATABLE   :: zd          ! new snow total depth
REAL, DIMENSION(:,:), ALLOCATABLE   :: zdepth      ! depth of each layer
REAL, DIMENSION(:,:,:), ALLOCATABLE :: zt          ! new snow temperature profile
REAL, DIMENSION(:),   ALLOCATABLE   :: zpsn        ! permanent snow fraction
REAL, DIMENSION(:,:), ALLOCATABLE   :: zwat        ! 
!
LOGICAL, DIMENSION(:,:), ALLOCATABLE :: gwork
INTEGER                              :: iwork
!
REAL              ::zrhosmax
REAL              ::zage_now
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------------
!
!*       1.    Snow where permanent snow is
!              ----------------------------
!
!* snow fraction must be at least equal to permanent snow fraction
!  The snow fraction is computed as Wsnow/(Wsnow+XWCRN)
!
!
IF (lhook) CALL dr_hook('PREP_PERM_SNOW',0,zhook_handle)
!
zrhosmax=xrhosmax
IF(tpsnow%SCHEME=='3-L'.OR.tpsnow%SCHEME=='CRO')THEN
  zrhosmax=xrhosmax_es
ENDIF
!
ALLOCATE(zpsn(SIZE(ptg,1)))
zpsn(:) = min( pperm_snow_frac(:,nvt_snow) , 0.9999 )
!
!* if no permanent snow present
!
IF (all(zpsn(:)==0.)) THEN
  DEALLOCATE(zpsn) 
  IF (lhook) CALL dr_hook('PREP_PERM_SNOW',1,zhook_handle)
  RETURN
END IF
!
!* total snow amount due to permanent snow
!
ALLOCATE(zwsnow_perm(SIZE(ptg,1)))
zwsnow_perm(:) = wsnow_from_snow_frac_ground(zpsn)
!
!* limitation of maximum snow amount
!
IF(i%LGLACIER)THEN
!  limited to 33.3 meters of aged snow
   zwsnow_perm(:) = min(zwsnow_perm(:),xhgla * zrhosmax )
ELSE
!  limited to 2. meters of aged snow
   zwsnow_perm(:) = min(zwsnow_perm(:),2.0 * zrhosmax )
ENDIF
!
!* permanent snow can be added only if deep soil temperature is below 5 C
!  (glaciers on subgrid mountains tops that are contained in the grid mesh are neglected)
!
IF(.NOT.i%LGLACIER)THEN
  WHERE(ptg(:,SIZE(ptg,2))>xtt+5.) zwsnow_perm(:) = 0.
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*       2.    Other parameters of new snow, except temperature
!              ------------------------------------------------
!
!* rho must be defined for snow 3-L before temperature and heat computations
!
ALLOCATE(gwork(SIZE(ptg,1),tpsnow%NLAYER))
!
DO jlayer=1,tpsnow%NLAYER
!
  gwork(:,jlayer)=.false.
!
  IF(i%LGLACIER)THEN
      WHERE(zwsnow_perm(:)>0.)gwork(:,jlayer)=.true.
  ELSE
      WHERE(zwsnow_perm(:)>0..AND.tpsnow%WSNOW(:,jlayer,ksnow)==0.)gwork(:,jlayer)=.true.
  ENDIF
!
!* rho
!
  WHERE(gwork(:,jlayer))
    tpsnow%RHO(:,jlayer,ksnow) = zrhosmax
  END WHERE
!
!* albedo
!
  IF(i%LGLACIER)THEN
    WHERE(gwork(:,jlayer))
         tpsnow%ALB(:,ksnow) = (xaglamax+xaglamin)/2.0
    END WHERE
  ELSE
    WHERE(gwork(:,jlayer))
         tpsnow%ALB(:,ksnow) = (xansmax+xansmin)/2.0
    END WHERE
  ENDIF
!
END DO
!
IF (tpsnow%SCHEME=='3-L'.OR.tpsnow%SCHEME=='CRO') THEN
!
! * optimized rho perm snow profile
!
  IF(i%LGLACIER.AND.tpsnow%NLAYER>=6)THEN
    WHERE(gwork(:,1))
      tpsnow%RHO(:,1,ksnow) = zrhol1
    END WHERE 
    IF(tpsnow%NLAYER>=6.AND.tpsnow%NLAYER<12)THEN
      WHERE(gwork(:,2))
       tpsnow%RHO(:,2,ksnow) = zrhol1 + 100.
      END WHERE 
      WHERE(gwork(:,3))
       tpsnow%RHO(:,3,ksnow) = zrhol1 + 250.
      END WHERE 
    ELSE
      DO jlayer=2,tpsnow%NLAYER
         WHERE(gwork(:,jlayer))
              tpsnow%RHO(:,jlayer,ksnow) = min(zrhosmax,tpsnow%RHO(:,jlayer-1,ksnow)+100.)
         END WHERE     
      ENDDO
    ENDIF
  ENDIF
!
! * Snow age profile
!
  DO jlayer=1,tpsnow%NLAYER/4
    WHERE(gwork(:,jlayer))
             tpsnow%AGE(:,jlayer,ksnow) = 365.0*float(jlayer-1)/ &
                                          float(tpsnow%NLAYER)
    END WHERE
  END DO
  DO jlayer=1+tpsnow%NLAYER/4,tpsnow%NLAYER
    WHERE(gwork(:,jlayer))
             tpsnow%AGE(:,jlayer,ksnow) = 3650.*float(jlayer-1)/ &
                                          float(tpsnow%NLAYER) 
    END WHERE
  END DO
!
  IF(i%LGLACIER)THEN
    WHERE(gwork(:,:))tpsnow%AGE(:,:,ksnow) = 0.0
  ENDIF
!
END IF
!
IF (tpsnow%SCHEME=='CRO') THEN
DO jlayer=1,tpsnow%NLAYER/4
  WHERE(gwork(:,jlayer))
           tpsnow%GRAN1(:,jlayer,ksnow) = min(-1.,-99.*     &
                  (1.-4*float(jlayer)/float(tpsnow%NLAYER))) 
           tpsnow%GRAN2(:,jlayer,ksnow) = 50. 
           tpsnow%HIST(:,jlayer,ksnow) = 0 
  END WHERE
END DO
DO jlayer=1+tpsnow%NLAYER/4,tpsnow%NLAYER
  WHERE(gwork(:,jlayer))
           tpsnow%GRAN1(:,jlayer,ksnow) = 99. 
           tpsnow%GRAN2(:,jlayer,ksnow) = 0.0003 
           tpsnow%HIST(:,jlayer,ksnow) = 0 
  END WHERE
END DO
END IF
!
!-------------------------------------------------------------------------------------
!
!*       3.    Modification of snow reservoir profile
!              --------------------------------------
!
!* initial snow content
!
ALLOCATE(zwsnow(SIZE(ptg,1)))
zwsnow(:) = 0.
DO jlayer=1,tpsnow%NLAYER
  zwsnow(:) = zwsnow(:) + tpsnow%WSNOW(:,jlayer,ksnow) 
END DO
!
!* new total snow content
!
zwsnow_perm(:) = max(zwsnow_perm(:),zwsnow(:))
!
!* new total snow depth
!
ALLOCATE(zd(SIZE(ptg,1)))
zd(:) = 0.
DO jlayer=1,tpsnow%NLAYER
  zd(:) = zd(:) + tpsnow%WSNOW(:,jlayer,ksnow)/tpsnow%RHO(:,jlayer,ksnow)
END DO
zd(:) = zd(:) + (zwsnow_perm(:)-zwsnow(:))/zrhosmax
!
!* modified snow content profile
!
SELECT CASE(tpsnow%SCHEME)
  CASE('D95','1-L','EBA')
    gwork(:,1)=.false.
    IF(i%LGLACIER)THEN
       WHERE(zwsnow(:)>=0..AND.tpsnow%WSNOW(:,1,ksnow)/=xundef)gwork(:,1)=.true.
    ELSE
       WHERE(zwsnow(:)==0..AND.tpsnow%WSNOW(:,1,ksnow)/=xundef)gwork(:,1)=.true.
    ENDIF
    WHERE(gwork(:,1))
      tpsnow%WSNOW(:,1,ksnow) = zwsnow_perm(:)
    END WHERE
  CASE('3-L','CRO')
    !* grid
    ALLOCATE(zdepth(SIZE(ptg,1),tpsnow%NLAYER))
    CALL snow3lgrid(zdepth,zd)
    DO jlayer=1,tpsnow%NLAYER
      WHERE(zwsnow(:)>= 0. .AND. tpsnow%WSNOW(:,jlayer,ksnow)/=xundef)
        tpsnow%WSNOW(:,jlayer,ksnow) = zdepth(:,jlayer) * tpsnow%RHO(:,jlayer,ksnow)
      END WHERE
   END DO
   DEALLOCATE(zdepth)

END SELECT
!
DEALLOCATE(zd)
!-------------------------------------------------------------------------------------
!
!*       4.    Temperature of new snow
!              -----------------------
!
ALLOCATE(zt(SIZE(tpsnow%WSNOW,1),SIZE(tpsnow%WSNOW,2),SIZE(tpsnow%WSNOW,3)))
!       
SELECT CASE(tpsnow%SCHEME)
  CASE('1-L')
    zt(:,:,:) = tpsnow%T (:,:,:)
  CASE('3-L','CRO')
    CALL snow_heat_to_t_wliq(tpsnow%HEAT,tpsnow%RHO,zt)
END SELECT
!
!* new snow is set to deep ground temperature
!
DO jlayer=1,tpsnow%NLAYER
!
  gwork(:,jlayer)=.false.
!
  IF(i%LGLACIER)THEN
      WHERE(zwsnow_perm(:)>0.)gwork(:,jlayer)=.true.
  ELSE
      WHERE(zwsnow_perm(:)>0. .AND. zwsnow(:)==0)gwork(:,jlayer)=.true.
  ENDIF
!  
  WHERE(gwork(:,jlayer))
      zt(:,jlayer,ksnow) = min(ptg(:,SIZE(ptg,2)),xtt)
  END WHERE
!
END DO
!
!
SELECT CASE(tpsnow%SCHEME)
  CASE('1-L')
    tpsnow%T (:,:,:) = zt(:,:,:)
  CASE('3-L','CRO')
    CALL snow_t_wliq_to_heat(tpsnow%HEAT,tpsnow%RHO,zt)
END SELECT
!
DEALLOCATE(zt   )
DEALLOCATE(gwork)
!
!
!-------------------------------------------------------------------------------------
!
!*       5.    Soil ice initialization for LGLACIER
!              -----------------------
!
ALLOCATE(zwat(SIZE(ptg,1),SIZE(ptg,2)))
!
IF(i%LGLACIER)THEN
!
  IF (i%CISBA == 'DIF') THEN
      iwork=i%NGROUND_LAYER
      zwat(:,:)=i%XWFC(:,:)
  ELSE
      iwork=2
      zwat(:,:)=i%XWSAT(:,:)
  ENDIF
!
  DO jlayer=1,iwork
     WHERE(pperm_snow_frac(:,nvt_snow)>0.0)
           i%XWGI(:,jlayer,ksnow) = max(i%XWGI(:,jlayer,ksnow),zwat(:,jlayer)*zpsn(:))
           i%XWG (:,jlayer,ksnow) = min(i%XWG (:,jlayer,ksnow),max(i%XWSAT(:,jlayer)-i%XWGI(:,jlayer,ksnow),xwgmin))
     END WHERE
     WHERE(i%XWG(:,jlayer,ksnow) /= xundef .AND. (i%XWG(:,jlayer,ksnow) + i%XWGI(:,jlayer,ksnow)) > i%XWSAT(:,jlayer) )
           i%XWGI(:,jlayer,ksnow) = i%XWSAT(:,jlayer)-i%XWG (:,jlayer,ksnow) !WGT<=WSAT
     END WHERE
  ENDDO
!
ENDIF
!
DEALLOCATE(zwat)
DEALLOCATE(zpsn)
!
!-------------------------------------------------------------------------------------
!
!*       6.    Masking where there is no snow
!              ------------------------------
!
 CALL mkflag_snow(tpsnow)
IF (lhook) CALL dr_hook('PREP_PERM_SNOW',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_perm_snow