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, &
                                 ig, u, &
                                hprogram,                       &
                                hfile,hfiletype,                &
                                hfilepgd,hfilepgdtype,          &
                                kluout,ounif,hsnsurf,kpatch,    &
                                kteb_patch, &
                                kl,tpsnow, tptime,              &
                                punif_wsnow, punif_rsnow,       &
                                punif_tsnow, punif_lwcsnow,     &
                                punif_asnow, osnow_ideal,       &
                                punif_sg1snow, punif_sg2snow,   &
                                punif_histsnow,punif_agesnow,   &                                
                                pf,pdepth,pvegtype,             &
                                pvegtype_patch,ppatch           )
!     #######################################################
!
!!****  *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
!!------------------------------------------------------------------
!
!
!
!
!
!
USE modd_data_cover_n, ONLY : data_cover_t
!
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
!
USE modd_type_snow
USE modd_type_date_surf, ONLY : date_time
!
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
!
USE modd_snow_par, ONLY : xansmax
!
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_heat_to_t_wliq
USE modi_vegtype_to_patch
!
USE modi_abor1_sfx
!
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(isba_grid_t), INTENT(INOUT) :: ig
TYPE(surf_atm_t), INTENT(INOUT) :: u
!
 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
INTEGER,            INTENT(IN)  :: kl        ! number of points
TYPE(surf_snow)                 :: tpsnow    ! 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(IN)  :: 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(OUT),OPTIONAL :: pf     ! output field (x,kpatch)
REAL,DIMENSION(:,:,:),INTENT(IN), OPTIONAL :: pdepth ! thickness of each snow layer
REAL,DIMENSION(:,:),  INTENT(IN), OPTIONAL :: pvegtype ! fraction of each vegtype
REAL,DIMENSION(:,:,:),  INTENT(IN), OPTIONAL :: pvegtype_patch ! fraction of each vegtype per patch
REAL,DIMENSION(:,:),  INTENT(IN), OPTIONAL :: ppatch ! fraction of each patch
!
!
!*      0.2    declarations of local variables
!
REAL, POINTER, DIMENSION(:,:,:)     :: zfieldin  ! field to interpolate horizontally
REAL, POINTER, DIMENSION(:,:)       :: zfield ! field to interpolate horizontally
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zfieldout ! field interpolated   horizontally
REAL, ALLOCATABLE, DIMENSION(:,:)   :: zd        ! snow depth (x, kpatch)
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zw        ! work array (x, fine   snow grid, kpatch)
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zheat     ! work array (x, output snow grid, kpatch)
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zgrid     ! grid array (x, output snow grid, kpatch)
!
LOGICAL                       :: gsnow_ideal
INTEGER                       :: jpatch    ! loop on patches
INTEGER                       :: jvegtype  ! loop on vegtypes
INTEGER                       :: jlayer    ! loop on layers
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
!*      1.     Does the field exist?
!
!
IF (lhook) CALL dr_hook('PREP_HOR_SNOW_FIELD',0,zhook_handle)
!
gsnow_ideal = .false.
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      2.     Reading of input  configuration (Grid and interpolation type)
!
IF (ounif) THEN
  gsnow_ideal = osnow_ideal
  CALL prep_snow_unif(kluout,hsnsurf,zfieldin, tptime, gsnow_ideal,       &
                      punif_wsnow, punif_rsnow, punif_tsnow,              &
                      punif_lwcsnow, punif_asnow, punif_sg1snow,          &
                      punif_sg2snow, punif_histsnow, punif_agesnow,       &
                      tpsnow%NLAYER                                       )
ELSE IF (hfiletype=='GRIB  ') THEN
  CALL prep_snow_grib(hprogram,hsnsurf,hfile,kluout,tpsnow%NLAYER,zfieldin)
ELSE IF (hfiletype=='MESONH' .OR. hfiletype=='ASCII ' .OR. hfiletype=='LFI   '.OR. hfiletype=='FA    ') THEN
  gsnow_ideal = osnow_ideal
  CALL prep_snow_extern(&
                        hprogram,hsnsurf,hfile,hfiletype,hfilepgd,hfilepgdtype,&
                        kluout,zfieldin,gsnow_ideal,tpsnow%NLAYER,kteb_patch)
ELSE IF (hfiletype=='BUFFER') THEN
  CALL prep_snow_buffer(ig, u, &
                        hprogram,hsnsurf,kluout,tpsnow%NLAYER,zfieldin)
ELSE
  CALL abor1_sfx('PREP_HOR_SNOW_FIELD: data file type not supported : '//hfiletype)
END IF
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      3.     Horizontal interpolation
!
ALLOCATE(zfieldout(kl,SIZE(zfieldin,2),SIZE(zfieldin,3)))
ALLOCATE(zfield(SIZE(zfieldin,1),SIZE(zfieldin,2)))
!
DO jvegtype = 1, SIZE(zfieldin,3)
  !* horizontal interpolation
  zfield=zfieldin(:,:,jvegtype)
  IF(present(pvegtype).AND.SIZE(zfieldin,3)==nvegtype) linterp(:) = (pvegtype(:,jvegtype) > 0.)
  IF(present(pdepth)) THEN
     jpatch = 1
     IF (kpatch>1) jpatch = vegtype_to_patch(jvegtype,kpatch)
     linterp(:) = (linterp(:).AND.pdepth(:,1,jpatch)>0..AND.pdepth(:,1,jpatch)<xundef)
  ENDIF
  !
  CALL hor_interpol(dtco, u, &
                    kluout,zfield,zfieldout(:,:,jvegtype))
  !
  linterp(:) = .true.
END DO
!
DEALLOCATE(zfield)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      4.     Transformation from vegtype grid to patch grid, if any
!
ALLOCATE(zw(SIZE(zfieldout,1),SIZE(zfieldout,2),kpatch))
!
zw = 0.
IF (SIZE(zfieldout,3)==nvegtype.AND.SIZE(pvegtype_patch,2)==nvegtype) THEN
  CALL vegtype_grid_to_patch_grid(kpatch,pvegtype_patch,ppatch,zfieldout,zw)
ELSE
  DO jpatch=1,kpatch
    zw(:,:,jpatch) = zfieldout(:,:,1)
  END DO
END IF
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      5.     Defines normalized output grid, if depths of snow layers are present
!
IF (present(pdepth) .AND. .NOT.gsnow_ideal) THEN
!
!* total snow depth
!
  ALLOCATE(zd(SIZE(tpsnow%WSNOW,1),kpatch))
  zd(:,:)=0.
  DO jpatch=1,kpatch
    DO jlayer=1,tpsnow%NLAYER
      WHERE (pdepth(:,jlayer,jpatch)/=xundef) zd(:,jpatch) = zd(:,jpatch) + pdepth(:,jlayer,jpatch)
    END DO
  END DO
!
!* grid at center of layers
!
  ALLOCATE(zgrid(SIZE(zw,1),tpsnow%NLAYER,kpatch))
  zgrid(:,1,:) = pdepth(:,1,:)
  IF(tpsnow%NLAYER>1)THEN
    DO jpatch=1,kpatch
      DO jlayer=2,tpsnow%NLAYER
        zgrid(:,jlayer,jpatch) = zgrid(:,jlayer-1,jpatch) + pdepth(:,jlayer,jpatch)
      ENDDO
    ENDDO
  ENDIF
!
!* normalized grid
!
  DO jpatch=1,kpatch
    DO jlayer=1,tpsnow%NLAYER
      WHERE (zd(:,jpatch)/=0.)
        zgrid(:,jlayer,jpatch) = zgrid(:,jlayer,jpatch) / zd(:,jpatch)
      ELSEWHERE
        zgrid(:,jlayer,jpatch) = 1.0
      END WHERE
    END DO
  END DO
!
  DEALLOCATE(zd)
!
ELSEIF (.NOT.gsnow_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
!
SELECT CASE (hsnsurf(1:3))
  !
  CASE('WWW')  ! total snow content
    !
    IF (gsnow_ideal) THEN
      pf(:,:,:) = zw(:,:,:)
    ELSE
      DO jlayer=1,SIZE(pf,2)
        pf(:,jlayer,:) = zw(:,1,:)
      ENDDO
    ENDIF
    !
    IF (present(ppatch)) THEN
      DO jlayer = 1,tpsnow%NLAYER
        WHERE(ppatch(:,:)==0.)
          pf(:,jlayer,:) = xundef
        END WHERE
      ENDDO
    ENDIF
  !
  CASE('DEP')  ! snow thickness
    !
    IF (gsnow_ideal) THEN
      pf(:,:,:) = zw(:,:,:)
    ELSE
      DO jlayer=1,SIZE(pf,2)
        pf(:,jlayer,:) = zw(:,jlayer,:)
      ENDDO
    ENDIF
    !
    IF (present(ppatch)) THEN
      DO jlayer = 1,tpsnow%NLAYER
        WHERE(ppatch(:,:)==0.)
          pf(:,jlayer,:) = xundef
        END WHERE
      ENDDO
    ENDIF
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CASE('RHO') 
    !
    IF (gsnow_ideal) THEN
      tpsnow%RHO(:,:,:) = zw(:,:,:)
    ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
      tpsnow%RHO(:,:,:) = zw(:,:,:)      
    ELSE
      !* interpolation on snow levels
      CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%RHO)
    ENDIF
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      DO jlayer=1,tpsnow%NLAYER
        WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%RHO(:,jlayer,jpatch) = xundef
      END DO
    END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CASE('ALB')
    !
    DO jpatch=1,kpatch
      tpsnow%ALB(:,jpatch) = zw(:,1,jpatch)
    END DO
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      WHERE(pdepth(:,1,jpatch)==0. .OR. pdepth(:,1,jpatch)==xundef)  tpsnow%ALB(:,jpatch) = xundef
    END DO
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
  CASE('HEA') 
    !
    IF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN
      !
      IF (gsnow_ideal) THEN
        tpsnow%HEAT(:,:,:) = zw(:,:,:)
      ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
        tpsnow%HEAT(:,:,:) = zw(:,:,:)
      ELSE
        !* interpolation of heat on snow levels
        CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%HEAT)
      ENDIF
      !
      !* mask for areas where there is no snow
      DO jpatch=1,kpatch
        DO jlayer=1,tpsnow%NLAYER
          WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%HEAT(:,jlayer,jpatch) = xundef
        END DO
      END DO
      !
    ELSE IF (tpsnow%SCHEME=='1-L') THEN
      !* interpolation of heat on snow levels
      ALLOCATE(zheat(SIZE(zfieldout,1),tpsnow%NLAYER,kpatch))
      !
      IF (gsnow_ideal) THEN
        zheat(:,:,:) = zw(:,:,:)
      ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
        zheat(:,:,:) = zw(:,:,:)
      ELSE
        CALL init_from_ref_grid(xgrid_snow,zw,zgrid,zheat)
      ENDIF
      !
      !* transformation from heat to temperature
      CALL snow_heat_to_t_wliq(zheat,tpsnow%RHO,tpsnow%T)
      WHERE (tpsnow%T>xtt) tpsnow%T = xtt
      DEALLOCATE(zheat)
      !
      !* mask for areas where there is no snow
      DO jpatch=1,kpatch
        DO jlayer=1,tpsnow%NLAYER
          WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%T(:,jlayer,jpatch) = xundef
        END DO
      END DO
      !
    END IF
  !
  !
  CASE('SG1')
    !
    IF (gsnow_ideal) THEN
      tpsnow%GRAN1(:,:,:) = zw(:,:,:)
    ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
      tpsnow%GRAN1(:,:,:) = zw(:,:,:)
    ELSE
      !* interpolation of heat on snow levels
      CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%GRAN1)
    ENDIF
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      DO jlayer=1,tpsnow%NLAYER
        WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%GRAN1(:,jlayer,jpatch) = xundef
      END DO
    END DO
    !
  CASE('SG2')
    !
    IF (gsnow_ideal) THEN
      tpsnow%GRAN2(:,:,:) = zw(:,:,:)
    ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
      tpsnow%GRAN2(:,:,:) = zw(:,:,:)
    ELSE
      !* interpolation of heat on snow levels
      CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%GRAN2)
    ENDIF
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      DO jlayer=1,tpsnow%NLAYER
        WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%GRAN2(:,jlayer,jpatch) = xundef
      END DO
    END DO
    !
  CASE('HIS')
    !
    IF (gsnow_ideal) THEN
      tpsnow%HIST(:,:,:) = zw(:,:,:)
    ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
      tpsnow%HIST(:,:,:) = zw(:,:,:)
    ELSE
      !* interpolation of heat on snow levels
      CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%HIST)
    ENDIF
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      DO jlayer=1,tpsnow%NLAYER
        WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%HIST(:,jlayer,jpatch) = xundef
      END DO
    END DO
    !
  CASE('AGE')
    !
    IF (tpsnow%SCHEME=='3-L'.AND.(.NOT.gsnow_ideal).AND.(.NOT.ounif))THEN
       tpsnow%AGE(:,:,:) = 0.0
    ELSE
      IF (gsnow_ideal) THEN
        tpsnow%AGE(:,:,:) = zw(:,:,:)
      ELSEIF(SIZE(zw,2)==tpsnow%NLAYER)THEN
        tpsnow%AGE(:,:,:) = zw(:,:,:)
      ELSE
        !* interpolation of heat on snow levels
        CALL init_from_ref_grid(xgrid_snow,zw,zgrid,tpsnow%AGE)
      ENDIF
    ENDIF
    !
    !* mask for areas where there is no snow
    DO jpatch=1,kpatch
      DO jlayer=1,tpsnow%NLAYER
        WHERE(pdepth(:,jlayer,jpatch)==0. .OR. pdepth(:,jlayer,jpatch)==xundef) tpsnow%AGE(:,jlayer,jpatch) = xundef
      END DO
    END DO
    !
END SELECT
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!*      7.     Deallocations
!
DEALLOCATE(zfieldin )
DEALLOCATE(zfieldout)
IF (present(pdepth) .AND. .NOT.gsnow_ideal) DEALLOCATE(zgrid    )
DEALLOCATE(zw       )
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  ! loop counter
REAL, DIMENSION(SIZE(PT1,1),SIZE(PT1,2)) :: zd1 ! input grid
REAL, DIMENSION(SIZE(PD2,1),SIZE(PD2,2)) :: zd2 ! output grid
INTEGER                       :: jpatch    ! loop on patches
REAL(KIND=JPRB) :: zhook_handle
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',0,zhook_handle)
DO jpatch=1,kpatch
  zd2(:,:) = 0.
  DO jl=1,SIZE(zd2,2)
    zd2(:,jl) = pd2(:,jl,jpatch)
  END DO
  !
  DO jl=1,SIZE(pt1,2)
    zd1(:,jl) = pgrid1(jl)
  END DO
  !
  CALL interp_grid_nat(zd1,pt1(:,:,jpatch),zd2,pt2(:,:,jpatch))
END DO
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',1,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END SUBROUTINE init_from_ref_grid
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_snow_field