prep_hor_snow_fields.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_fields (DTCO, G, U, GCP, HPROGRAM,HSURF,&
                                HFILE,HFILETYPE,             &
                                HFILEPGD,HFILEPGDTYPE,       &
                                KLUOUT,OUNIF,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, KSIZE_P, KR_P,&
                                PPATCH, OKEY                 )  
!     #######################################################
!
!
!!****  *PREP_HOR_SNOW_FIELDS* - prepares all snow fields for one surface scheme.
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      B. Decharme 10/2013, Phasage Arpege-Climat
!!      B. Decharme 04/2014, Init permsnow
!!      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 modd_surf_par,       ONLY : xundef
USE modd_snow_par,       ONLY : xaglamin, xaglamax
!
USE modd_data_cover_par, ONLY : nvegtype
!
USE mode_prep_ctl, ONLY : prep_ctl, prep_ctl_int_part2, prep_ctl_int_part4
!
USE modi_allocate_gr_snow
USE modi_prep_hor_snow_field
USE mode_snow3l
USE modi_open_aux_io_surf
USE modi_read_surf
USE modi_close_aux_io_surf
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      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=7),   INTENT(IN)  :: HSURF     ! type of field
 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
INTEGER,            INTENT(IN)  :: KPATCH    ! patch number for output scheme
INTEGER,            INTENT(IN) :: KTEB_PATCH
INTEGER,            INTENT(IN)  :: KL        ! number of points
TYPE(nsurf_snow), INTENT(INOUT) :: TNPSNOW   ! snow fields
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 ), OPTIONAL :: PVEGTYPE_PATCH ! fraction of each vegtype per patch
INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: KSIZE_P
INTEGER, DIMENSION(:,:), INTENT(IN), OPTIONAL :: KR_P
REAL,DIMENSION(:,:),    INTENT(IN ), OPTIONAL :: PPATCH         ! fraction of each patch
LOGICAL,                INTENT(OUT), OPTIONAL :: OKEY
!
!
!*      0.2    declarations of local variables
!
TYPE(surf_snow), POINTER :: SK
 CHARACTER(LEN=10)                   :: YSNSURF   ! type of field
REAL,ALLOCATABLE,DIMENSION(:)     :: ZWRHO     ! total snow content from rho profile alone
REAL,ALLOCATABLE,DIMENSION(:,:,:) :: ZDEPTH    ! snow depth of each layer
REAL,ALLOCATABLE,DIMENSION(:)     :: ZDTOT     ! total snow depth
INTEGER, DIMENSION(KPATCH) :: ISIZE_P
INTEGER,DIMENSION(KL,KPATCH) :: IR_P    ! fraction of each patch
REAL,DIMENSION(KL,KPATCH)         :: ZPATCH    ! fraction of each patch
REAL,DIMENSION(:,:,:), ALLOCATABLE  :: ZVEGTYPE_PATCH    ! fraction of each vegtype per patch
!
INTEGER                             :: JP, ISNOW_NLAYER    ! loop counter on patches
INTEGER                             :: JL, JI, ISIZE   ! loop counter on layers
INTEGER                             :: IVERSION  ! surface version
 CHARACTER(LEN=16)                   :: YRECFM    ! record name
INTEGER                             :: IRESP     ! error return code
LOGICAL                             :: GGLACIER
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!---------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PREP_HOR_SNOW_FIELDS',0,zhook_handle)
!
isnow_nlayer = tnpsnow%AL(1)%NLAYER
!
IF (PRESENT(ksize_p)) THEN
  isize_p(:) = ksize_p(:)
ELSE
  isize_p(:) = kl
ENDIF
!
IF (PRESENT(kr_p)) THEN
  ir_p(:,:) = kr_p(:,:)
ELSE
  DO ji = 1,kl
    ir_p(ji,:) = ji
  ENDDO
ENDIF
!
IF (PRESENT(ppatch)) THEN
   zpatch = ppatch
ELSE
   zpatch = 1.
ENDIF
IF (PRESENT(pvegtype_patch)) THEN
  ALLOCATE(zvegtype_patch(kl,SIZE(pvegtype_patch,2),kpatch))
  zvegtype_patch = pvegtype_patch
ELSE
  ALLOCATE(zvegtype_patch(kl,1,kpatch))
  zvegtype_patch = 1.
ENDIF
!
!*      1.     Allocation of output field
!
DO jp = 1,kpatch
  CALL allocate_gr_snow(tnpsnow%AL(jp),isize_p(jp))
ENDDO
!
!---------------------------------------------------------------------------
!
!*      2.     Find if PERMSNOW must be done
!
IF(PRESENT(okey))THEN
  !  
  IF ( (hfiletype=='MESONH' .OR. hfiletype=='ASCII ' .OR. hfiletype=='LFI   '.OR. hfiletype=='FA    ') &
       .AND. (hsurf=='SN_VEG ')  ) THEN
    !       
    CALL open_aux_io_surf(hfile,hfiletype,'FULL  ')
    yrecfm='VERSION'
    CALL read_surf(hfiletype,yrecfm,iversion,iresp)    
    CALL close_aux_io_surf(hfile,hfiletype)
    !  
    IF(iversion>7)THEN       
      CALL open_aux_io_surf(hfile,hfiletype,'NATURE')
      yrecfm='GLACIER'
      CALL read_surf(hfiletype,yrecfm,gglacier,iresp)    
      CALL close_aux_io_surf(hfile,hfiletype)  
      IF(gglacier)okey=.false.
    ENDIF
    !
  ENDIF
  !  
  IF(osnow_ideal)okey=.false.
  !
ENDIF
!
!---------------------------------------------------------------------------
!
!*      3.     Treatment of total snow content (kg/m2)
!
ysnsurf='WWW'//hsurf
 CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram, hfile, hfiletype, hfilepgd, hfilepgdtype,  &
                         kluout, ounif, ysnsurf, 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,  &                      
                         zvegtype_patch, zpatch, isize_p, ir_p         )  
!
!----------------------------------------------------------------------------
!
!*      4.     Treatment of total snow depth
!
DO jp = 1,kpatch
  ALLOCATE(tnpsnow%AL(jp)%DEPTH(isize_p(jp),isnow_nlayer))
ENDDO
!
ysnsurf='DEP'//hsurf
 CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram, hfile, hfiletype, hfilepgd, hfilepgdtype,  &
                         kluout, ounif, ysnsurf, 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,  &
                         zvegtype_patch, zpatch, isize_p, ir_p        )
!
!* snow layer thickness definition
!
ALLOCATE(zdepth(kl,isnow_nlayer,kpatch))
!
DO jp = 1,kpatch
  !
  IF (osnow_ideal .OR. isnow_nlayer==1 .OR. tnpsnow%AL(jp)%SCHEME=='3-L') THEN
    zdepth(1:isize_p(jp),:,jp) = tnpsnow%AL(jp)%DEPTH(:,:)
  ELSEIF (tnpsnow%AL(1)%SCHEME=='CRO') THEN
    ALLOCATE(zdtot(isize_p(jp)))
    zdtot(:) = 0.0
    DO jl=1,isnow_nlayer
       zdtot(:) = zdtot(:) + tnpsnow%AL(jp)%DEPTH(:,jl)
    END DO
    CALL snow3lgrid(zdepth(1:isize_p(jp),:,jp),zdtot(:))
    DEALLOCATE(zdtot)
  ENDIF
  !
ENDDO
!
DO jp = 1,kpatch
  DEALLOCATE(tnpsnow%AL(jp)%DEPTH)
ENDDO
!
!----------------------------------------------------------------------------
!
!*      4.     Snow density profile
!              --------------------
!
!* density profile
ysnsurf='RHO'//hsurf
 CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,           &
                         kluout,ounif,ysnsurf, 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,       &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth  )  
!
!----------------------------------------------------------------------------
!
!*      5.     Snow water content profile
!              --------------------------
!
IF (.NOT.osnow_ideal) THEN
  !
  ALLOCATE(zwrho(kl))
  ALLOCATE(zdtot(kl))
  !
  DO jp = 1,kpatch
    !
    sk => tnpsnow%AL(jp)
    !
    zwrho(:) = 0.0
    zdtot(:) = 0.0
    !
    isize = isize_p(jp)
    !
    !* snow depth estimated from rho profile
    DO jl=1,isnow_nlayer
      WHERE (zpatch(1:isize,jp)>0. .AND. sk%RHO(:,jl)/=xundef)
        zwrho(1:isize) = zwrho(1:isize) + sk%RHO(:,jl) * zdepth(1:isize,jl,jp)
      ELSEWHERE
        zwrho(1:isize) = xundef
      END WHERE
    ENDDO
    !
    DO jl = 1,isnow_nlayer
      !* modification of snow depth: coherence between rho profile, total snow and total depth
      WHERE(zpatch(1:isize,jp)>0. .AND. zwrho(1:isize)/=0. &
           .AND. zwrho(1:isize)/=xundef .AND. sk%WSNOW(:,1)>0.0)
        zdtot(1:isize) = zdtot(1:isize) + zdepth(1:isize,jl,jp) * sk%WSNOW(:,1) / zwrho(1:isize)
      ENDWHERE
    END DO
    !
    IF (isnow_nlayer > 1) THEN
      CALL snow3lgrid(zdepth(1:isize,:,jp),zdtot(1:isize))
    ELSE
      zdepth(1:isize,1,jp) = zdtot(1:isize)
    ENDIF
    !
    !* snow content profile for each grid level
    DO jl=1,isnow_nlayer
      WHERE(zpatch(1:isize,jp)>0..AND.sk%RHO(:,jl)/=xundef.AND.zdtot(1:isize)>0.)
        sk%WSNOW(:,jl) = sk%RHO(:,jl) * zdepth(1:isize,jl,jp)
      ELSEWHERE(zpatch(1:isize,jp)>0..AND.(sk%RHO(:,jl)==xundef.OR.zdtot(1:isize)==0.0))
        sk%WSNOW(:,jl) = 0.0
      ELSEWHERE
        sk%WSNOW(:,jl) = xundef
      END WHERE
    END DO
    !
  END DO
  !
  DEALLOCATE(zwrho)
  DEALLOCATE(zdtot)
  !
ENDIF
!
!----------------------------------------------------------------------------
!
!*      6.     Albedo, snow heat content, and age
!              ----------------------------------
!
!* albedo
ysnsurf='ALB'//hsurf
 CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,         &
                         kluout,ounif,ysnsurf, 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,      &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth  ) 
!
IF (tnpsnow%AL(1)%SCHEME/='D95') THEN
  !
  !* heat in snowpack profile
  ysnsurf='HEA'//hsurf
  CALL prep_hor_snow_field(dtco, g, u, gcp, &
                           hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,          &
                           kluout,ounif,ysnsurf, 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,      &
                           zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth    )
  !
ENDIF
!
IF (tnpsnow%AL(1)%SCHEME=='CRO'.OR. tnpsnow%AL(1)%SCHEME=='3-L') THEN
  !
  !* age in snowpack profile
  ysnsurf='AGE'//hsurf
  CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,        &
                         kluout,ounif,ysnsurf, 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,      &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth  )   
  !
  DO jp = 1,kpatch
    WHERE(tnpsnow%AL(jp)%WSNOW(:,1)>0.0 .AND. tnpsnow%AL(jp)%WSNOW(:,1)/=xundef .AND. &
          tnpsnow%AL(jp)%AGE(:,1)==0.0  .AND. tnpsnow%AL(jp)%ALB(:)<xaglamin)
      tnpsnow%AL(jp)%ALB(:)=(xaglamin+xaglamax)/2.0
    ENDWHERE
  ENDDO
  !
ENDIF
!
!----------------------------------------------------------------------------
!
!*      7.     Crocus specific parameters
!              --------------------------
!
IF (tnpsnow%AL(1)%SCHEME=='CRO') THEN
  !
  ysnsurf='SG1'//hsurf
  CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,        &
                         kluout,ounif,ysnsurf, 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,      &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth )   
  !
  ysnsurf='SG2'//hsurf
  CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,        &
                         kluout,ounif,ysnsurf, 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,      &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth  )   
  !
  ysnsurf='HIS'//hsurf
  CALL prep_hor_snow_field(dtco, g, u, gcp, &
                         hprogram,hfile,hfiletype,hfilepgd,hfilepgdtype,        &
                         kluout,ounif,ysnsurf, 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,      &
                         zvegtype_patch, zpatch, isize_p, ir_p, pdepth=zdepth  )   
  !
ENDIF
!
!*      8.     Deallocations
!
DEALLOCATE(zdepth  )
DEALLOCATE(zvegtype_patch)
!
IF (lhook) CALL dr_hook('PREP_HOR_SNOW_FIELDS',1,zhook_handle)
!
!----------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_snow_fields