prep_hor_snow_field.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_hor_snow_field (DTCO, G, U, GCP, HPROGRAM,      &
                                HFILE,HFILETYPE,                &
                                HFILEPGD,HFILEPGDTYPE,          &
                                KLUOUT,OUNIF,HSNSURF,KPATCH,    &
                                KTEB_PATCH,                     &
                                KL,TNPSNOW, TPTIME,             &
                                PUNIF_WSNOW, PUNIF_RSNOW,       &
                                PUNIF_TSNOW, PUNIF_LWCSNOW,     &
                                PUNIF_ASNOW, OSNOW_IDEAL,       &
                                PUNIF_SG1SNOW, PUNIF_SG2SNOW,   &
                                PUNIF_HISTSNOW,PUNIF_AGESNOW, YDCTL,   &
                                PVEGTYPE_PATCH, PPATCH,         &
                                KSIZE_P, KR_P, PDEPTH  )
!     #######################################################
!
!!****  *PREP_HOR_SNOW_FIELD* - reads, interpolates and prepares a snow field
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      P. Le Moigne 10/2005, Phasage Arome
!!      B. Decharme  10/2013, Phasage Arpège-Climat
!!      M. Lafaysse 11/2012, snow liquid water content
!!      B. Decharme  04/2014, external init with FA files
!!                            new init for ES
!!      P. Marguinaud10/2014, Support for a 2-part PREP
!!------------------------------------------------------------------
!
USE modd_data_cover_n, ONLY : data_cover_t
!
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_type_snow
USE modd_type_date_surf, ONLY : date_time
!
USE mode_prep_ctl, ONLY : prep_ctl, prep_ctl_int_part2, prep_ctl_int_part4
!
USE modd_grid_grib, ONLY : cinmodel
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc
!
USE modd_csts,           ONLY : xtt
USE modd_prep_snow,      ONLY : xgrid_snow
USE modd_surf_par,       ONLY : xundef
USE modd_data_cover_par, ONLY : nvegtype, nvt_snow
USE modd_prep,           ONLY : linterp,cinterp_type,cingrid_type, cmask
!
USE modd_snow_par, ONLY : xansmax
!
USE modi_prep_grib_grid
USE modi_prep_snow_grib
USE modi_prep_snow_unif
USE modi_prep_snow_extern
USE modi_prep_snow_buffer
USE modi_hor_interpol
USE modi_vegtype_grid_to_patch_grid
USE modi_snow_t_wliq_to_heat
USE modi_vegtype_to_patch
USE modi_pack_same_rank
USE modi_get_prep_interp
USE modi_put_on_all_vegtypes
!
USE mode_snow3l, ONLY : snow3lgrid
!
USE modi_abor1_sfx
!
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(data_cover_t), INTENT(INOUT) :: DTCO
!
TYPE(grid_t), INTENT(INOUT) :: G
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
type(prep_ctl),    INTENT (INOUT) :: ydctl
!
 CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=28),  INTENT(IN)  :: HFILE     ! file name
 CHARACTER(LEN=6),   INTENT(IN)  :: HFILETYPE ! file type
 CHARACTER(LEN=28),  INTENT(IN)  :: HFILEPGD     ! file name
 CHARACTER(LEN=6),   INTENT(IN)  :: HFILEPGDTYPE ! file type
INTEGER,            INTENT(IN)  :: KLUOUT    ! logical unit of output listing
LOGICAL,            INTENT(IN)  :: OUNIF     ! flag for prescribed uniform field
 CHARACTER(LEN=10)               :: HSNSURF   ! type of field
INTEGER,            INTENT(IN)  :: KPATCH    ! patch number for output scheme
INTEGER,            INTENT(IN) :: KTEB_PATCH
TYPE(nsurf_snow), INTENT(INOUT) :: TNPSNOW    ! snow fields
INTEGER,            INTENT(IN)  :: KL        ! number of points
TYPE(date_time),    INTENT(IN)  :: TPTIME    ! date and time
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_WSNOW ! prescribed snow content (kg/m2)
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_RSNOW ! prescribed density (kg/m3)
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_TSNOW ! prescribed temperature (K)
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_LWCSNOW ! prescribed snow liquid water content (kg/m3)
REAL,               INTENT(IN)  :: PUNIF_ASNOW ! prescribed albedo (-)
LOGICAL,            INTENT(INOUT)  :: OSNOW_IDEAL
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_SG1SNOW ! 
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_SG2SNOW ! 
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_HISTSNOW ! 
REAL, DIMENSION(:), INTENT(IN)  :: PUNIF_AGESNOW ! 
!
REAL,DIMENSION(:,:,:),  INTENT(IN) :: PVEGTYPE_PATCH ! fraction of each vegtype per patch
REAL,DIMENSION(:,:),  INTENT(IN) :: PPATCH ! fraction of each patch
INTEGER, DIMENSION(:), INTENT(IN) :: KSIZE_P
INTEGER,DIMENSION(:,:),  INTENT(IN) :: KR_P
!
REAL,DIMENSION(:,:,:),INTENT(IN), OPTIONAL :: PDEPTH ! thickness of each snow layer
!
!*      0.2    declarations of local variables
!
TYPE fout
  REAL, DIMENSION(:,:), ALLOCATABLE :: ZOUT
END TYPE fout
TYPE nfout
  TYPE(fout), DIMENSION(:), ALLOCATABLE :: AL
END TYPE nfout
type(nfout) :: zw
TYPE(surf_snow), POINTER :: SK
!
REAL, POINTER, DIMENSION(:,:,:)     :: ZFIELDIN   ! field to interpolate horizontally
REAL, POINTER, DIMENSION(:,:,:) :: ZFIELDOUTP ! field interpolated   horizontally
REAL, POINTER, DIMENSION(:,:,:) :: ZFIELDOUTV !
REAL, ALLOCATABLE, DIMENSION(:)   :: ZD        ! snow depth (x, kpatch)
REAL, ALLOCATABLE, DIMENSION(:,:) :: ZTEMP
REAL, ALLOCATABLE, DIMENSION(:,:) :: ZWLIQ     ! liquid water snow pack content
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZGRID     ! grid array (x, output snow grid, kpatch)
!
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDEPTH, ZPATCH
!
type(date_time)              :: tztime_grib    ! current date and time
INTEGER                       :: JP, IP    ! loop on patches
INTEGER                       :: JL    ! loop on layers
INTEGER :: INFOMPI, INL, INP, ISNOW_NLAYER, IMASK, JI
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!----------------------------------------------------------------------------
!
!*      1.     Does the field exist?
!
!
IF (lhook) CALL dr_hook('PREP_HOR_SNOW_FIELD',0,zhook_handle)
!
isnow_nlayer = tnpsnow%AL(1)%NLAYER
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      2.     Reading of input  configuration (Grid and interpolation type)
!
NULLIFY (zfieldin, zfieldoutp, zfieldoutv)
!
IF (ydctl%LPART1) THEN
  IF (ounif) THEN
    CALL prep_snow_unif(kluout,hsnsurf,zfieldin, tptime, osnow_ideal,       &
                        punif_wsnow, punif_rsnow, punif_tsnow,              &
                        punif_lwcsnow, punif_asnow, punif_sg1snow,          &
                        punif_sg2snow, punif_histsnow, punif_agesnow,       &
                        isnow_nlayer                                     )
  ELSE IF (hfiletype=='GRIB  ') THEN
    CALL prep_grib_grid(hfile,kluout,cinmodel,cingrid_type,cinterp_type,tztime_grib)            
    IF (nrank==npio) CALL prep_snow_grib(hprogram,hsnsurf,hfile,kluout,isnow_nlayer,zfieldin)        
  ELSE IF (hfiletype=='MESONH' .OR. hfiletype=='ASCII ' .OR. hfiletype=='LFI   '&
          .OR. hfiletype=='FA    '.OR. hfiletype=='AROME '.OR.hfiletype=='NC    ') THEN
    CALL prep_snow_extern(gcp,hprogram,hsnsurf,hfile,hfiletype,hfilepgd,hfilepgdtype, &
                          kluout,zfieldin,osnow_ideal,isnow_nlayer,kteb_patch)
  ELSE IF (hfiletype=='BUFFER') THEN
    CALL prep_snow_buffer(g, u, hprogram,hsnsurf,kluout,isnow_nlayer,zfieldin)
  ELSE
    CALL abor1_sfx('PREP_HOR_SNOW_FIELD: data file type not supported : '//hfiletype)
  END IF
ENDIF
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      3.     Horizontal interpolation
!
 CALL prep_ctl_int_part2 (ydctl, hsnsurf, 'SNOW', cmask, zfieldin)
!
IF (ydctl%LPART3) THEN
!
  IF (nrank==npio) THEN
    inl = SIZE(zfieldin,2)
    inp = SIZE(zfieldin,3)
  ELSEIF (.NOT.ASSOCIATED(zfieldin)) THEN
   ALLOCATE(zfieldin(0,0,0))
  ENDIF
!
  IF (nproc>1) THEN
#ifdef SFX_MPI
    CALL mpi_bcast(inl,kind(inl)/4,mpi_integer,npio,ncomm,infompi)
    CALL mpi_bcast(inp,kind(inp)/4,mpi_integer,npio,ncomm,infompi)
#endif
  ENDIF
!    
  ALLOCATE(zfieldoutp(kl,inl,inp))
!
! ZPATCH is the array of output patches put on the input patches
  ALLOCATE(zpatch(kl,inp))
  zpatch(:,:) = 0.
!
! if the number of input patches is NVEGTYPE
  IF (inp==nvegtype) THEN
    DO jp = 1,nvegtype
    ! each vegtype takes the output contribution of the patch it is in
      ip = vegtype_to_patch(jp,kpatch)
      zpatch(:,jp) = pvegtype_patch(:,jp,ip)
    ENDDO
  ENDIF
!
  CALL get_prep_interp(inp,kpatch,zpatch,ppatch,zpatch,kr_p)
!
! the same for depth that is defined on the output patches
  IF (PRESENT(pdepth)) THEN
  !
    ALLOCATE(zdepth(kl,inp))
    zdepth(:,:) = 0.
  !
    IF (inp==nvegtype) THEN
      DO jp = 1,nvegtype
        ip = vegtype_to_patch(jp,kpatch)
        zdepth(:,jp) = pdepth(:,1,ip)
      ENDDO
    ENDIF
    ! 
    CALL get_prep_interp(inp,kpatch,zdepth,pdepth(:,1,:),zdepth,kr_p)
    !
  ENDIF
!
  DO jp = 1, inp
    ! ZDEPTH and ZPATCH are defined on the size on the patch for the snow: use of
    ! the mask
    IF (PRESENT(pdepth)) THEN
      linterp(:) = ( zdepth(:,jp) /= 0. .AND. zdepth(:,jp) /= xundef )
    ENDIF
    linterp(:) = (linterp(:) .AND. zpatch(:,jp)>0.)
    !* horizontal interpolation
    CALL hor_interpol(dtco, u, gcp, kluout,zfieldin(:,:,jp),zfieldoutp(:,:,jp))
    !
    linterp(:) = .true.
  END DO
  !
  DEALLOCATE(zfieldin, zpatch )
  IF (PRESENT(pdepth)) DEALLOCATE(zdepth)
  !
ENDIF
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      4.     Transformation from vegtype grid to patch grid, if any
!
 CALL prep_ctl_int_part4 (ydctl, hsnsurf, 'SNOW', cmask, zfieldin, zfieldoutp)

IF (ydctl%LPART5) THEN
!
  ALLOCATE(zw%AL(kpatch))
  !
  IF (kpatch/=inp.and.inp/=1) THEN
    !
    ALLOCATE(zfieldoutv(kl,inl,nvegtype))
    CALL put_on_all_vegtypes(kl,inl,inp,nvegtype,zfieldoutp,zfieldoutv)
    !
    !*      6.     Transformation from vegtype grid to patch grid
    !
    DEALLOCATE(zfieldoutp)
    !
    DO jp = 1,kpatch
      !
      ALLOCATE(zw%AL(jp)%ZOUT(ksize_p(jp),inl))
      !
      CALL vegtype_grid_to_patch_grid(jp, kpatch, pvegtype_patch(1:ksize_p(jp),:,jp),  &
                                      ppatch(1:ksize_p(jp),jp), kr_p(1:ksize_p(jp),jp), &
                                      zfieldoutv, zw%AL(jp)%ZOUT)
    ENDDO  
    !
    DEALLOCATE(zfieldoutv)
    !
  ELSEIF (inp==1) THEN
    !
    DO jp = 1,kpatch
      !
      ALLOCATE(zw%AL(jp)%ZOUT(ksize_p(jp),inl))
      !
      CALL pack_same_rank(kr_p(1:ksize_p(jp),jp),zfieldoutp(:,:,1),zw%AL(jp)%ZOUT)
      !
    ENDDO    
    !
  ELSE
    !
    DO jp = 1,kpatch
      !
      ALLOCATE(zw%AL(jp)%ZOUT(ksize_p(jp),inl))
      !
      CALL pack_same_rank(kr_p(1:ksize_p(jp),jp),zfieldoutp(:,:,jp),zw%AL(jp)%ZOUT)
      !
    ENDDO
    !
    DEALLOCATE(zfieldoutp)  
    !
  ENDIF
  !

  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  !*      5.     Defines normalized output grid, if depths of snow layers are present
  !
  IF ( PRESENT(pdepth) .AND. .NOT.osnow_ideal ) THEN
    !
    ALLOCATE(zd(SIZE(pdepth,1)))
    !
    ALLOCATE(zgrid(SIZE(pdepth,1),isnow_nlayer,kpatch))
    DO jp = 1,kpatch
      !* total snow depth
      !
      zd(:)=0.
      DO jl = 1,isnow_nlayer
        WHERE (pdepth(1:ksize_p(jp),jl,jp)/=xundef) zd(1:ksize_p(jp)) = zd(1:ksize_p(jp)) + pdepth(1:ksize_p(jp),jl,jp)
      END DO
      !
      !* grid at center of layers
      !
      zgrid(1:ksize_p(jp),1,jp) = pdepth(1:ksize_p(jp),1,jp)
      IF(isnow_nlayer>1) THEN
        DO jl = 2,isnow_nlayer
          zgrid(1:ksize_p(jp),jl,jp) = zgrid(1:ksize_p(jp),jl-1,jp) + pdepth(1:ksize_p(jp),jl,jp)
        ENDDO
      ENDIF
      !
      ! * normalized grid
      !
      DO jl=1,isnow_nlayer
        WHERE (zd(1:ksize_p(jp))/=0.)
          zgrid(1:ksize_p(jp),jl,jp) = zgrid(1:ksize_p(jp),jl,jp) / zd(1:ksize_p(jp))
        ELSEWHERE
          zgrid(1:ksize_p(jp),jl,jp) = 1.0
        END WHERE
      END DO
      !
    ENDDO
    !
    DEALLOCATE(zd)
    !
  ELSEIF (.NOT.osnow_ideal) THEN
    IF (hsnsurf(1:3)=='RHO' .OR. hsnsurf(1:3)=='HEA') THEN
      WRITE(kluout,*) 'when interpolation profiles of snow pack quantities,'
      WRITE(kluout,*) 'depth of snow layers must be given'
      CALL abor1_sfx('PREP_HOR_SNOW_FIELD: DEPTH OF SNOW LAYERS NEEDED')
    END IF
  END IF
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  !*      6.     Return to historical variable
  !
  DO jp = 1,kpatch
    !
    sk => tnpsnow%AL(jp)
    !
    SELECT CASE (hsnsurf(1:3))
      !
      CASE('WWW')  ! total snow content
        !
        IF (osnow_ideal) THEN
          sk%WSNOW(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSE
          DO jl=1,SIZE(sk%WSNOW,2)
            sk%WSNOW(:,jl) = zw%AL(jp)%ZOUT(:,1)
          ENDDO
        ENDIF
        !
        DO jl = 1,isnow_nlayer
          WHERE(ppatch(1:ksize_p(jp),jp)==0.)
            sk%WSNOW(:,jl) = xundef
          END WHERE
        ENDDO
        !
      CASE('DEP')  ! snow thickness
        !
        IF (osnow_ideal) THEN
          sk%DEPTH(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSE
          CALL snow3lgrid(sk%DEPTH(:,:),zw%AL(jp)%ZOUT(:,1))                
          !DO JL=1,SIZE(SK%DEPTH,2)
          !  SK%DEPTH(:,JL) = ZW%AL(JP)%ZOUT(:,1)
          !ENDDO
        ENDIF
        !
        DO jl = 1,isnow_nlayer
          WHERE(ppatch(1:ksize_p(jp),jp)==0.)
            sk%DEPTH(:,jl) = xundef
          END WHERE
        ENDDO
      !
      !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      !
      CASE('RHO') 
        !
        IF (osnow_ideal) THEN
          sk%RHO(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSEIF(inl==1) THEN
          DO jl = 1,isnow_nlayer
            sk%RHO(:,jl) = zw%AL(jp)%ZOUT(:,1)
          ENDDO      
        ELSE
          !* interpolation on snow levels
          CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%RHO)
        ENDIF
        !
        !* mask for areas where there is no snow
        DO jl=1,isnow_nlayer
          WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%RHO(:,jl) = xundef
        END DO
        !
        !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        !
      CASE('ALB')
        !
        sk%ALB(:) = zw%AL(jp)%ZOUT(:,1)
        !
        !* mask for areas where there is no snow
        WHERE(pdepth(1:ksize_p(jp),1,jp)==0. .OR. pdepth(1:ksize_p(jp),1,jp)==xundef)  sk%ALB(:) = xundef
      !
      !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      !
      CASE('HEA') 
        !
        IF (sk%SCHEME=='3-L' .OR. sk%SCHEME=='CRO') THEN
          !
          ALLOCATE(ztemp(ksize_p(jp),isnow_nlayer))
          ALLOCATE(zwliq(ksize_p(jp),isnow_nlayer))
          zwliq(:,:) = 0.0
          !
          IF (osnow_ideal) THEN
            ztemp(:,:) = zw%AL(jp)%ZOUT(:,:)
          ELSEIF (inl==1) THEN
            DO jl = 1,isnow_nlayer 
              ztemp(:,jl) = zw%AL(jp)%ZOUT(:,1)
            ENDDO
          ELSE
            !* interpolation of heat on snow levels
            CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),ztemp)
          ENDIF
          !
          CALL snow_t_wliq_to_heat(sk%HEAT,sk%RHO,ztemp,zwliq)
          DEALLOCATE(ztemp)
          DEALLOCATE(zwliq)
          !
          !* mask for areas where there is no snow
          DO jl=1,isnow_nlayer
            WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%HEAT(:,jl) = xundef
          END DO
          !
        ELSE IF (sk%SCHEME=='1-L') THEN
          !* interpolation of heat on snow levels
          !
          IF (osnow_ideal) THEN
            sk%T(:,:) = zw%AL(jp)%ZOUT(:,:)
          ELSEIF(inl==1) THEN
            DO jl = 1,isnow_nlayer
              sk%T(:,jl) = zw%AL(jp)%ZOUT(:,1)
            ENDDO        
          ELSE
            CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%T)
          ENDIF
          !
          !* transformation from heat to temperature
          WHERE (sk%T>xtt) sk%T = xtt
          !
          !* mask for areas where there is no snow
          DO jl=1,isnow_nlayer
            WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%T(:,jl) = xundef
          END DO
          !
        END IF
        !
        !
      CASE('SG1')
        !
        IF (osnow_ideal) THEN
          sk%GRAN1(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSEIF(inl==1) THEN
         DO jl = 1,isnow_nlayer
             sk%GRAN1(:,jl) = zw%AL(jp)%ZOUT(:,1)
          ENDDO      
        ELSE
          !* interpolation of heat on snow levels
          CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%GRAN1)
        ENDIF
        !
        !* mask for areas where there is no snow
        DO jl=1,isnow_nlayer
          WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%GRAN1(:,jl) = xundef
        END DO
        !
      CASE('SG2')
        !
        IF (osnow_ideal) THEN
          sk%GRAN2(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSEIF(SIZE(zw%AL(jp)%ZOUT,2)==1) THEN
          DO jl = 1,isnow_nlayer
            sk%GRAN2(:,jl) = zw%AL(jp)%ZOUT(:,1)
          ENDDO      
        ELSE
          !* interpolation of heat on snow levels
          CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%GRAN2)
        ENDIF
        !
        !* mask for areas where there is no snow
        DO jl=1,isnow_nlayer
          WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%GRAN2(:,jl) = xundef
        END DO
        !
      CASE('HIS')
        !
        IF (osnow_ideal) THEN
          sk%HIST(:,:) = zw%AL(jp)%ZOUT(:,:)
        ELSEIF(inl==1) THEN
          DO jl = 1,isnow_nlayer
            sk%HIST(:,jl) = zw%AL(jp)%ZOUT(:,1)
          ENDDO      
        ELSE
          !* interpolation of heat on snow levels
          CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%HIST)
        ENDIF
        !
        !* mask for areas where there is no snow
        DO jl=1,isnow_nlayer
          WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%HIST(:,jl) = xundef
        END DO
        !
      CASE('AGE')
        !
        IF (sk%SCHEME=='3-L'.AND.(.NOT.osnow_ideal).AND.(.NOT.ounif))THEN
          sk%AGE(:,:) = 0.0
        ELSE
          IF (osnow_ideal) THEN
            sk%AGE(:,:) = zw%AL(jp)%ZOUT(:,:)
          ELSEIF(inl==1) THEN
            DO jl = 1,isnow_nlayer
              sk%AGE(:,jl) = zw%AL(jp)%ZOUT(:,1)
            ENDDO        
          ELSE
            !* interpolation of heat on snow levels
            CALL init_from_ref_grid(xgrid_snow,zw%AL(jp)%ZOUT,zgrid(1:ksize_p(jp),:,jp),sk%AGE)
          ENDIF
        ENDIF
        !
        !* mask for areas where there is no snow
        DO jl=1,isnow_nlayer
          WHERE(pdepth(1:ksize_p(jp),jl,jp)==0. .OR. pdepth(1:ksize_p(jp),jl,jp)==xundef) sk%AGE(:,jl) = xundef
        END DO
        !
    END SELECT
    !
  ENDDO
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      7.     Deallocations
!
  IF (PRESENT(pdepth) .AND. .NOT.osnow_ideal) DEALLOCATE(zgrid    )
  DO jp =1,kpatch
    DEALLOCATE(zw%AL(jp)%ZOUT)
  ENDDO
  DEALLOCATE(zw%AL)
!
ENDIF
!
IF (lhook) CALL dr_hook('PREP_HOR_SNOW_FIELD',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
CONTAINS
!
!-------------------------------------------------------------------------------------
!
SUBROUTINE init_from_ref_grid(PGRID1,PT1,PD2,PT2)
!
USE modi_interp_grid_nat
!
REAL, DIMENSION(:,:), INTENT(IN)  :: PT1    ! variable profile
REAL, DIMENSION(:),     INTENT(IN)  :: PGRID1 ! normalized grid
REAL, DIMENSION(:,:), INTENT(IN)  :: PD2    ! output layer thickness
REAL, DIMENSION(:,:), INTENT(OUT) :: PT2    ! variable profile
!
INTEGER                                  :: JL, JL1  ! loop counter
REAL, DIMENSION(SIZE(PT1,1),SIZE(PGRID1)) :: ZT1
REAL, DIMENSION(SIZE(PT1,1),SIZE(PGRID1)) :: ZD1 ! input grid
REAL, DIMENSION(SIZE(PD2,1),SIZE(PD2,2)) :: ZD2 ! output grid
INTEGER                       :: JP    ! loop on patches
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',0,zhook_handle)
!
zd2(:,:) = 0.
DO jl=1,SIZE(zd2,2)
  zd2(:,jl) = pd2(:,jl)
END DO
!
DO jl=1,SIZE(pgrid1)
  jl1 = min(jl,SIZE(pt1,2))
  zt1(:,jl) = pt1(:,jl1)
  zd1(:,jl) = pgrid1(jl)
END DO
!
 CALL interp_grid_nat(zd1,zt1(:,:),zd2,pt2(:,:))
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',1,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END SUBROUTINE init_from_ref_grid
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_snow_field