prep_snow_extern.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_snow_extern (&
                             hprogram,hsurf,hfile,hfiletype,hfilepgd,hfilepgdtype,&
                            kluout,pfield,osnow_ideal,klayer,kteb_patch)
!     #################################################################################
!
!
!!****  *PREP_SNOW_EXTERN*  
!!
!!    PURPOSE
!!    -------
!       Read and prepare initial snow fields from external files
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    
!!    AUTHOR
!!    ------
!!         * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    ?
!!       02/2014 E. Martin : cor. for passing from from multilayer to a single layer
!!      B. Decharme  04/2014, external init with FA files
!!                            improve vertical interpolation
!-------------------------------------------------------------------------------
!
!
!
!
!
USE modd_type_snow
USE modd_prep,           ONLY : cingrid_type, cinterp_type
USE modd_prep_snow,      ONLY : xgrid_snow, ngrid_level
USE modd_data_cover_par, ONLY : nvegtype
USE modd_surf_par,       ONLY : xundef
USE modd_csts,           ONLY : xtt
!
USE mode_snow3l
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_town_presence
USE modi_put_on_all_vegtypes
USE modi_abor1_sfx
USE modi_prep_grid_extern
USE modi_open_aux_io_surf
USE modi_close_aux_io_surf
USE modi_allocate_gr_snow
USE modi_interp_grid_nat
USE modi_read_gr_snow
USE modi_read_surf
USE modi_snow_t_wliq_to_heat
USE modi_read_teb_patch
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
 CHARACTER(LEN=6),   INTENT(IN)  :: hprogram  ! program calling surf. schemes
 CHARACTER(LEN=10),  INTENT(IN)  :: hsurf     ! type of field
 CHARACTER(LEN=28),  INTENT(IN)  :: hfile     ! name of file
 CHARACTER(LEN=6),   INTENT(IN)  :: hfiletype ! type of file
 CHARACTER(LEN=28),  INTENT(IN)  :: hfilepgd     ! name of file
 CHARACTER(LEN=6),   INTENT(IN)  :: hfilepgdtype ! type of file
INTEGER,            INTENT(IN)  :: kluout    ! logical unit of output listing
REAL,DIMENSION(:,:,:), POINTER  :: pfield    ! field to interpolate horizontally
LOGICAL,            INTENT(IN)  :: osnow_ideal
INTEGER,            INTENT(IN)  :: klayer    ! Number of layer of output snow scheme
INTEGER,            INTENT(IN) :: kteb_patch
!
!*      0.2    declarations of local variables
!
TYPE(surf_snow)                    :: tzsnow ! snow characteristics

REAL, DIMENSION(:,:,:), ALLOCATABLE :: zfield       ! work field on input snow grid
REAL, DIMENSION(:,:,:), ALLOCATABLE :: zfield_fine  ! work field on fine snow grid
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ztemp        ! snow temperature
REAL, DIMENSION(:,:,:), ALLOCATABLE :: zwliq        ! liquid water snow pack content
REAL, DIMENSION(:,:),   ALLOCATABLE :: zd           ! total snow depth
REAL, DIMENSION(:,:,:), ALLOCATABLE :: zdepth       ! thickness of each layer (m)
REAL, DIMENSION(:,:,:), ALLOCATABLE :: zgrid        ! normalized input grid
!
LOGICAL                           :: gtown          ! town variables written in the file
 CHARACTER(LEN=12)                 :: yrecfm         ! record name
INTEGER                           :: iresp          ! error return code
INTEGER                           :: iversion       ! SURFEX version
LOGICAL                           :: gold_name      ! old name flag 
INTEGER                           :: ibugfix        ! SURFEX bug version
INTEGER                           :: ivegtype       ! actual number of vegtypes
INTEGER                           :: jlayer         ! loop on snow vertical grids
INTEGER                           :: ji             ! loop on pts
INTEGER                           :: ini
 CHARACTER(LEN=8)                  :: yarea          ! area treated ('ROOF','ROAD','VEG ')
 CHARACTER(LEN=3)                  :: yprefix        ! prefix to identify patch
INTEGER                           :: ipatch         ! number of input patch
INTEGER                           :: iteb_patch     ! number of input patch for TEB
INTEGER                           :: jpatch         ! loop on patch
 CHARACTER(LEN=6)                  :: ymask          ! type of tile mask
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------------
!
!*      3.     Area being treated
!              ------------------
!
IF (lhook) CALL dr_hook('PREP_SNOW_EXTERN',0,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
yarea='        '
yarea(1:4) = hsurf(7:10)
!
IF (yarea(1:4)=='VEG ') THEN
  ivegtype = nvegtype
  ymask = 'NATURE'
  yprefix = '   '  
ELSE
  ivegtype = 1
  ymask    = 'TOWN  '
  ipatch   = 1
  yprefix = '   '  
END IF
!
!*      1.     Preparation of IO for reading in the file
!              -----------------------------------------
!
!* Note that all points are read, even those without physical meaning.
!  These points will not be used during the horizontal interpolation step.
!  Their value must be defined as XUNDEF.
!
!* reading of version of the file being read
 CALL open_aux_io_surf(&
                      hfilepgd,hfilepgdtype,'FULL  ')
 CALL read_surf(&
               hfilepgdtype,'VERSION',iversion,iresp)
 CALL read_surf(&
               hfilepgdtype,'BUG',ibugfix,iresp)
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
!
gold_name=(iversion<7 .OR. (iversion==7 .AND. ibugfix<3))
!
 CALL open_aux_io_surf(&
                      hfilepgd,hfilepgdtype,ymask)
!
IF (yarea(1:4)=='VEG ') THEN
  yrecfm = 'PATCH_NUMBER'
  CALL read_surf(&
               hfilepgdtype,yrecfm,ipatch,iresp)
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
ELSE
  IF (.NOT.gold_name) THEN
     IF (yarea(1:4)=='ROOF') yarea(1:4) = 'RF  '
     IF (yarea(1:4)=='ROAD') yarea(1:4) = 'RD  '
   ENDIF
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
  CALL read_teb_patch(&
                      hfilepgd,hfilepgdtype,iteb_patch)
  IF (iteb_patch>1) THEN
    WRITE(yprefix,fmt='(A,I1,A)') 'T',min(kteb_patch,iteb_patch),'_'
  END IF  
END IF
!
!
!-------------------------------------------------------------------------------------
!
!*      2.     Reading of grid
!              ---------------
!
 CALL open_aux_io_surf(&
                      hfilepgd,hfilepgdtype,'FULL  ')
!
 CALL prep_grid_extern(&
                      hfilepgdtype,kluout,cingrid_type,cinterp_type,ini)
!
!-------------------------------------------------------------------------------------
!
!*      4.     Reading of snow data
!              ---------------------
!
IF (yarea(1:2)=='RO' .OR. yarea(1:2)=='GA' .OR. yarea(1:2)=='RF' .OR. yarea(1:2)=='RD') THEN
  CALL town_presence(&
                     hfilepgdtype,gtown)
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
  IF (.NOT. gtown) THEN
    tzsnow%SCHEME='1-L'
    tzsnow%NLAYER=1
    CALL allocate_gr_snow(tzsnow,ini,ipatch)
  ELSE
    CALL open_aux_io_surf(&
                      hfile,hfiletype,ymask)
    CALL read_gr_snow(&
                      hfiletype,trim(yarea),yprefix,ini,ipatch,tzsnow, &
                      hdir='A',kversion=iversion,kbugfix=ibugfix)
    CALL close_aux_io_surf(hfile,hfiletype)
  ENDIF
ELSE
  CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
  CALL open_aux_io_surf(&
                      hfile,hfiletype,ymask)
  CALL read_gr_snow(&
                      hfiletype,trim(yarea),yprefix,ini,ipatch,tzsnow, &
                    hdir='A',kversion=iversion,kbugfix=ibugfix)
  CALL close_aux_io_surf(hfile,hfiletype)
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      5.     Total snow content
!              ------------------
!
SELECT CASE (hsurf(1:3))
  CASE ('WWW')
    IF (osnow_ideal) THEN
      ALLOCATE(zfield(ini,klayer,ipatch))
      zfield(:,:,:) = tzsnow%WSNOW(:,1:klayer,:)
      ALLOCATE(pfield(ini,klayer,ivegtype))
      CALL put_on_all_vegtypes(ini,klayer,ipatch,ivegtype,zfield,pfield)
    ELSE
      ALLOCATE(zfield(ini,1,ipatch))
      zfield(:,1,:) = 0.0
      DO jlayer=1,tzsnow%NLAYER
        zfield(:,1,:) = zfield(:,1,:) + tzsnow%WSNOW(:,jlayer,:)
      END DO 
      WHERE ( zfield(:,1,:)>xundef ) zfield(:,1,:)=xundef
      ALLOCATE(pfield(ini,1,ivegtype))
      CALL put_on_all_vegtypes(ini,1,ipatch,ivegtype,zfield,pfield)
    ENDIF
    DEALLOCATE(zfield)
!
!-------------------------------------------------------------------------------------
!
!*      6.     Albedo
!              ------
!
  CASE ('ALB')
    ALLOCATE(zfield(ini,1,ipatch))
    zfield(:,1,:) = tzsnow%ALB(:,:)
    !
    ALLOCATE(pfield(ini,1,ivegtype))
    CALL put_on_all_vegtypes(ini,1,ipatch,ivegtype,zfield,pfield)
    DEALLOCATE(zfield)
!
!-------------------------------------------------------------------------------------
!
!*      7.     Total depth to snow grid
!              ------------------------
!
  CASE ('DEP')
    ALLOCATE(zfield_fine(ini,klayer,ipatch))
    IF (osnow_ideal) THEN
      zfield_fine(:,:,:) = tzsnow%WSNOW(:,1:klayer,:)/tzsnow%RHO(:,1:klayer,:)
      WHERE(tzsnow%WSNOW(:,1:klayer,:)==xundef) zfield_fine(:,:,:)=xundef
    ELSE          
      ALLOCATE(zdepth(ini,tzsnow%NLAYER,ipatch))
      zdepth(:,:,:) = tzsnow%WSNOW(:,:,:)/tzsnow%RHO(:,:,:)
      WHERE(tzsnow%WSNOW(:,:,:)==xundef) zdepth(:,:,:)=xundef
      IF(tzsnow%NLAYER/=klayer)THEN
        !* total depth
        ALLOCATE(zd(ini,ipatch))
        zd(:,:) = 0.0
        DO jpatch=1,ipatch
           DO jlayer=1,tzsnow%NLAYER
              DO ji=1,ini
                 IF(zdepth(ji,jlayer,jpatch)/=xundef)THEN
                   zd(ji,jpatch) = zd(ji,jpatch) + zdepth(ji,jlayer,jpatch)
                 ENDIF
              ENDDO
           ENDDO
        ENDDO
        !* fine grid
        DO jpatch=1,ipatch
           CALL snow3lgrid(zfield_fine(:,:,jpatch),zd(:,jpatch))
        ENDDO  
        DEALLOCATE(zd)
      ELSE
        zfield_fine(:,:,:)=zdepth(:,:,:)
      ENDIF
      DEALLOCATE(zdepth)
    ENDIF
    ALLOCATE(pfield(ini,klayer,ivegtype))
    CALL put_on_all_vegtypes(ini,klayer,ipatch,ivegtype,zfield_fine,pfield)
    DEALLOCATE(zfield_fine)
!
!-------------------------------------------------------------------------------------
!
!*      8.     Density or heat profile
!              -----------------------
!
  CASE ('RHO','HEA','SG1','SG2','HIS','AGE')
    ALLOCATE(zfield(ini,tzsnow%NLAYER,ipatch))
!
    SELECT CASE (tzsnow%SCHEME)
      CASE ('D95','1-L','EBA')
        ALLOCATE(zfield_fine(ini,ngrid_level,ipatch))      
        !* computes output physical variable
        IF (hsurf(1:3)=='RHO') zfield(:,1,:) = tzsnow%RHO(:,1,:)
        IF (hsurf(1:3)=='HEA') THEN
          ALLOCATE(ztemp(ini,tzsnow%NLAYER,ipatch))
          ALLOCATE(zwliq(ini,tzsnow%NLAYER,ipatch))
          IF (tzsnow%SCHEME=='D95'.OR.tzsnow%SCHEME=='EBA') ztemp(:,1,:) = xtt-2.
          IF (tzsnow%SCHEME=='1-L') ztemp(:,1,:) = tzsnow%T(:,1,:)
          zwliq(:,:,:) = 0.0
          CALL snow_t_wliq_to_heat(zfield,tzsnow%RHO,ztemp,zwliq)
          DEALLOCATE(ztemp)
          DEALLOCATE(zwliq)
        END IF
        IF (hsurf(1:3)=='SG1') zfield(:,1,:) = -20.0
        IF (hsurf(1:3)=='SG2') zfield(:,1,:) = 80.0
        IF (hsurf(1:3)=='HIS') zfield(:,1,:) = 0.0
        IF (hsurf(1:3)=='AGE') zfield(:,1,:) = 3.0
        !* put profile on fine snow grid
        DO jlayer=1,ngrid_level
          zfield_fine(:,jlayer,:) = zfield(:,1,:)
        END DO
        ALLOCATE(pfield(ini,ngrid_level,ivegtype))
        CALL put_on_all_vegtypes(ini,ngrid_level,ipatch,ivegtype,zfield_fine,pfield)

      CASE ('3-L','CRO')
        !* input physical variable
        IF (hsurf(1:3)=='RHO') zfield(:,:,:) = tzsnow%RHO (:,1:tzsnow%NLAYER,:)
        IF (hsurf(1:3)=='HEA') zfield(:,:,:) = tzsnow%HEAT(:,1:tzsnow%NLAYER,:)
        IF (hsurf(1:3)=='AGE') zfield(:,:,:) = tzsnow%AGE (:,1:tzsnow%NLAYER,:)
        IF (tzsnow%SCHEME=='CRO')THEN
           IF (hsurf(1:3)=='SG1') zfield(:,:,:) = tzsnow%GRAN1(:,1:tzsnow%NLAYER,:)
           IF (hsurf(1:3)=='SG2') zfield(:,:,:) = tzsnow%GRAN2(:,1:tzsnow%NLAYER,:)
           IF (hsurf(1:3)=='HIS') zfield(:,:,:) = tzsnow%HIST (:,1:tzsnow%NLAYER,:)
        ELSE
           IF (hsurf(1:3)=='SG1') zfield(:,:,:) = -20.0
           IF (hsurf(1:3)=='SG2') zfield(:,:,:) = 80.0
           IF (hsurf(1:3)=='HIS') zfield(:,:,:) = 0.0
        ENDIF
        !
        IF (osnow_ideal) THEN
          ALLOCATE(zfield_fine(ini,klayer,ipatch))
          zfield_fine(:,:,:) = zfield(:,:,:)
          ALLOCATE(pfield(ini,klayer,ivegtype))
          CALL put_on_all_vegtypes(ini,klayer,ipatch,ivegtype,zfield_fine,pfield)
        ELSE
          ALLOCATE(zfield_fine(ini,ngrid_level,ipatch))  
          !
          !* input snow layer thickness
          ALLOCATE(zdepth(ini,tzsnow%NLAYER,ipatch))
          DO jpatch=1,ipatch
              zdepth(:,:,jpatch) = tzsnow%WSNOW(:,:,jpatch)/tzsnow%RHO(:,:,jpatch)
          END DO
          !
          !* total depth
          ALLOCATE(zd(ini,ipatch))
          zd(:,:) = 0.
          DO jlayer=1,tzsnow%NLAYER
            zd(:,:) = zd(:,:) + zdepth(:,jlayer,:)
          ENDDO
          !
          !* input normalized grid
          ALLOCATE(zgrid(ini,tzsnow%NLAYER,ipatch))
          DO jpatch=1,ipatch
             DO ji=1,ini
                IF(zd(ji,jpatch)==0.0)THEN
                   DO jlayer = 1,tzsnow%NLAYER
                      zgrid(ji,jlayer,jpatch)=REAL(jlayer)/REAL(tzsnow%nlayer)
                   ENDDO
                ELSE
                   DO jlayer = 1,tzsnow%NLAYER
                      IF(jlayer==1)THEN
                        zgrid(ji,jlayer,jpatch)=zdepth(ji,jlayer,jpatch)/ zd(ji,jpatch)
                      ELSE
                        zgrid(ji,jlayer,jpatch) = zgrid(ji,jlayer-1,jpatch) + zdepth(ji,jlayer,jpatch)/zd(ji,jpatch)
                      ENDIF
                   ENDDO
                ENDIF
             ENDDO
          ENDDO
          DEALLOCATE(zdepth)
          DEALLOCATE(zd)
          !    
          !* interpolation of profile onto fine normalized snow grid
          DO jpatch=1,ipatch
            CALL interp_grid_nat(zgrid(:,:,jpatch),zfield(:,:,jpatch),    &
                             xgrid_snow(:), zfield_fine(:,:,jpatch))
          END DO
          DEALLOCATE(zgrid)
          ALLOCATE(pfield(ini,ngrid_level,ivegtype))
          CALL put_on_all_vegtypes(ini,ngrid_level,ipatch,ivegtype,zfield_fine,pfield)
        ENDIF
      END SELECT
    !
    !* put field form patch to all vegtypes
    DEALLOCATE(zfield)
    DEALLOCATE(zfield_fine)
!
END SELECT
!
!-------------------------------------------------------------------------------------
!
!*      9.     End of IO
!              ---------
!
IF (lhook) CALL dr_hook('PREP_SNOW_EXTERN',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
END SUBROUTINE prep_snow_extern