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 (IO, KK, PK, PEK)
!          ################################################
!
!
!!****  *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_surfex_mpi, ONLY : nrank, npio, ncomm, nproc
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_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_vegtype_to_patch
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
!
#ifdef SFX_MPI
include "mpif.h"
#endif
!
!*      0.1    declarations of arguments
!
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
!*      0.2    declarations of local parameter
!
REAL, PARAMETER :: ZRHOL1 = 150.
!
!*      0.3    declarations of local variables
!
INTEGER                             :: JL      ! 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
!
#ifdef SFX_MPI
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ISTATUS
#endif
INTEGER :: INFOMPI
INTEGER :: ISNOW          ! patch number where permanent snow is
!
REAL, DIMENSION(0:NPROC-1) :: ZPSN0
REAL :: ZSUM_PSN
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(pek%TSNOW%SCHEME=='3-L'.OR.pek%TSNOW%SCHEME=='CRO')THEN
  zrhosmax=xrhosmax_es
ENDIF
!
ALLOCATE(zpsn(SIZE(pek%XTG,1)))

zpsn(:) = min( pk%XVEGTYPE_PATCH(:,nvt_snow) , 0.9999 )
!
!* if no permanent snow present
!
zsum_psn = sum(zpsn(:))
IF (nproc>1) THEN
#ifdef SFX_MPI
  CALL mpi_allgather(zsum_psn,kind(zsum_psn)/4,mpi_real,&
                     zpsn0,kind(zpsn0)/4,mpi_real,ncomm,infompi)
#endif
ELSE
  zpsn0(:) = zsum_psn
ENDIF

IF (all(zpsn0(:)==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(pek%XTG,1)))
zwsnow_perm(:) = wsnow_from_snow_frac_ground(zpsn)
!
!* limitation of maximum snow amount
!
IF(io%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.io%LGLACIER)THEN
  WHERE(pek%XTG(:,SIZE(pek%XTG,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(pek%XTG,1),pek%TSNOW%NLAYER))
!
DO jl=1,pek%TSNOW%NLAYER
!
  gwork(:,jl)=.false.
!
  IF(io%LGLACIER)THEN
    WHERE(zwsnow_perm(:)>0.)gwork(:,jl)=.true.
  ELSE
    WHERE(zwsnow_perm(:)>0..AND.pek%TSNOW%WSNOW(:,jl)==0.)gwork(:,jl)=.true.
  ENDIF
!
!* rho
!
  WHERE(gwork(:,jl))
    pek%TSNOW%RHO(:,jl) = zrhosmax
  END WHERE
!
!* albedo
!
  IF(io%LGLACIER)THEN
    WHERE(gwork(:,jl))
      pek%TSNOW%ALB(:) = (xaglamax+xaglamin)/2.0
    END WHERE
  ELSE
    WHERE(gwork(:,jl))
      pek%TSNOW%ALB(:) = (xansmax+xansmin)/2.0
    END WHERE
  ENDIF
!
END DO
!
IF (pek%TSNOW%SCHEME=='3-L'.OR.pek%TSNOW%SCHEME=='CRO') THEN
!
! * optimized rho perm snow profile
!
  IF(io%LGLACIER.AND.pek%TSNOW%NLAYER>=6)THEN
    WHERE(gwork(:,1))
      pek%TSNOW%RHO(:,1) = zrhol1
    END WHERE 
    IF(pek%TSNOW%NLAYER>=6.AND.pek%TSNOW%NLAYER<12)THEN
      WHERE(gwork(:,2))
       pek%TSNOW%RHO(:,2) = zrhol1 + 100.
      END WHERE 
      WHERE(gwork(:,3))
       pek%TSNOW%RHO(:,3) = zrhol1 + 250.
      END WHERE 
    ELSE
      DO jl=2,pek%TSNOW%NLAYER
        WHERE(gwork(:,jl))
          pek%TSNOW%RHO(:,jl) = min(zrhosmax,pek%TSNOW%RHO(:,jl-1)+100.)
        END WHERE     
      ENDDO
    ENDIF
  ENDIF
!
! * Snow age profile
!
  DO jl=1,pek%TSNOW%NLAYER/4
    WHERE(gwork(:,jl))
      pek%TSNOW%AGE(:,jl) = 365.0*float(jl-1)/ float(pek%TSNOW%NLAYER)
    END WHERE
  END DO
  DO jl=1+pek%TSNOW%NLAYER/4,pek%TSNOW%NLAYER
    WHERE(gwork(:,jl))
      pek%TSNOW%AGE(:,jl) = 3650.*float(jl-1)/ float(pek%TSNOW%NLAYER) 
    END WHERE
  END DO
!
  IF(io%LGLACIER)THEN
    WHERE(gwork(:,:))pek%TSNOW%AGE(:,:) = 0.0
  ENDIF
!
END IF
!
IF (pek%TSNOW%SCHEME=='CRO') THEN
DO jl=1,pek%TSNOW%NLAYER/4
  WHERE(gwork(:,jl))
    pek%TSNOW%GRAN1(:,jl) = min(-1.,-99.* (1.-4*float(jl)/float(pek%TSNOW%NLAYER))) 
    pek%TSNOW%GRAN2(:,jl) = 50. 
    pek%TSNOW%HIST(:,jl) = 0 
  END WHERE
END DO
DO jl=1+pek%TSNOW%NLAYER/4,pek%TSNOW%NLAYER
  WHERE(gwork(:,jl))
    pek%TSNOW%GRAN1(:,jl) = 99. 
    pek%TSNOW%GRAN2(:,jl) = 0.0003 
    pek%TSNOW%HIST(:,jl) = 0 
  END WHERE
END DO
END IF
!
!-------------------------------------------------------------------------------------
!
!*       3.    Modification of snow reservoir profile
!              --------------------------------------
!
!* initial snow content
!
ALLOCATE(zwsnow(SIZE(pek%XTG,1)))
zwsnow(:) = 0.
DO jl=1,pek%TSNOW%NLAYER
  zwsnow(:) = zwsnow(:) + pek%TSNOW%WSNOW(:,jl) 
END DO
!
!* new total snow content
!
zwsnow_perm(:) = max(zwsnow_perm(:),zwsnow(:))
!
!* new total snow depth
!
ALLOCATE(zd(SIZE(pek%XTG,1)))
zd(:) = 0.
DO jl=1,pek%TSNOW%NLAYER
  zd(:) = zd(:) + pek%TSNOW%WSNOW(:,jl)/pek%TSNOW%RHO(:,jl)
END DO
zd(:) = zd(:) + (zwsnow_perm(:)-zwsnow(:))/zrhosmax
!
!* modified snow content profile
!
SELECT CASE(pek%TSNOW%SCHEME)
  CASE('D95','1-L','EBA')
    gwork(:,1)=.false.
    IF(io%LGLACIER)THEN
       WHERE(zwsnow(:)>=0..AND.pek%TSNOW%WSNOW(:,1)/=xundef)gwork(:,1)=.true.
    ELSE
       WHERE(zwsnow(:)==0..AND.pek%TSNOW%WSNOW(:,1)/=xundef)gwork(:,1)=.true.
    ENDIF
    WHERE(gwork(:,1))
      pek%TSNOW%WSNOW(:,1) = zwsnow_perm(:)
    END WHERE
  CASE('3-L','CRO')
    !* grid
    ALLOCATE(zdepth(SIZE(pek%XTG,1),pek%TSNOW%NLAYER))
    CALL snow3lgrid(zdepth,zd)
    DO jl=1,pek%TSNOW%NLAYER
      WHERE(zwsnow(:)>= 0. .AND. pek%TSNOW%WSNOW(:,jl)/=xundef)
        pek%TSNOW%WSNOW(:,jl) = zdepth(:,jl) * pek%TSNOW%RHO(:,jl)
      END WHERE
   END DO
   DEALLOCATE(zdepth)

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