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 (GCP,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_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_surfex_mpi, ONLY : nrank, npio
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 yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_town_presence
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_dealloc_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_snow_heat_to_t_wliq
USE modi_read_teb_patch
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
 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(INOUT)  :: 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 :: ZHEAT
REAL, DIMENSION(:),   ALLOCATABLE :: ZD           ! total snow depth
REAL, DIMENSION(:,:), ALLOCATABLE :: ZDEPTH       ! thickness of each layer (m)
REAL, DIMENSION(:,:), ALLOCATABLE :: ZGRID        ! normalized input grid
!
INTEGER, DIMENSION(:), ALLOCATABLE :: IMASK_P
!
LOGICAL                           :: GTOWN          ! town variables written in the file
 CHARACTER(LEN=12)                 :: YRECFM         ! record name
INTEGER                           :: IRESP          ! error return code
INTEGER                           :: IVERSION_PGD, IVERSION_PREP       ! SURFEX version
LOGICAL                           :: GOLD_NAME      ! old name flag 
INTEGER                           :: IBUGFIX_PGD, IBUGFIX_PREP        ! SURFEX bug version
INTEGER                           :: IVEGTYPE       ! actual number of vegtypes
INTEGER                           :: JL         ! 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                           :: JP         ! 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_pgd,iresp,hdir='-')
 CALL read_surf(hfilepgdtype,'BUG',ibugfix_pgd,iresp,hdir='-')
gold_name=(iversion_pgd<7 .OR. (iversion_pgd==7 .AND. ibugfix_pgd<3))
!
 CALL prep_grid_extern(gcp,hfilepgdtype,kluout,cingrid_type,cinterp_type,ini)
!
 CALL town_presence(hfilepgdtype,gtown,hdir='-')
!
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
!
 CALL open_aux_io_surf(hfilepgd,hfilepgdtype,ymask)
!
IF (yarea(1:4)=='VEG ') THEN
  yrecfm = 'PATCH_NUMBER'
  CALL read_surf(hfilepgdtype,yrecfm,ipatch,iresp,hdir='-')
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
  IF (gtown) THEN   
    CALL read_teb_patch(hfilepgd,hfilepgdtype,iversion_pgd,ibugfix_pgd,iteb_patch,hdir='-')
  ELSE
    iteb_patch = 1
  ENDIF    
  IF (iteb_patch>1) THEN
    WRITE(yprefix,fmt='(A,I1,A)') 'T',min(kteb_patch,iteb_patch),'_'
  END IF  
END IF
!
 CALL close_aux_io_surf(hfilepgd,hfilepgdtype)
!
 CALL open_aux_io_surf(hfile,hfiletype,'FULL  ')
 CALL read_surf(hfiletype,'VERSION',iversion_prep,iresp,hdir='-')
 CALL read_surf(hfiletype,'BUG',ibugfix_prep,iresp,hdir='-')
 CALL close_aux_io_surf(hfile,hfiletype)
!
!-------------------------------------------------------------------------------------
!
!*      2.     Reading of grid
!              ---------------
!
IF (nrank/=npio) ini = 0
!
!-------------------------------------------------------------------------------------
!
!*      4.     Reading of snow data
!              ---------------------
!
ALLOCATE(imask_p(ini))
DO ji = 1,ini
  imask_p(ji) = ji
ENDDO
!
DO jp = 1,ipatch
  !
  IF (yarea(1:2)=='RO' .OR. yarea(1:2)=='GA' .OR. yarea(1:2)=='RF' .OR. yarea(1:2)=='RD') THEN
    IF (.NOT. gtown) THEN
      tzsnow%SCHEME='1-L'
      tzsnow%NLAYER=1
      CALL allocate_gr_snow(tzsnow,ini)
    ELSE
      CALL open_aux_io_surf(hfile,hfiletype,ymask)
      CALL read_gr_snow(hfiletype,trim(yarea),yprefix,ini,ini,imask_p,0,tzsnow, &
                        hdir='E',kversion=iversion_prep,kbugfix=ibugfix_prep)
      CALL close_aux_io_surf(hfile,hfiletype)
    ENDIF
  ELSE
    CALL open_aux_io_surf(hfile,hfiletype,ymask)
    CALL read_gr_snow(hfiletype,trim(yarea),yprefix,ini,ini,imask_p,jp,tzsnow, &
                      hdir='E',kversion=iversion_prep,kbugfix=ibugfix_prep,knpatch=ipatch)
    CALL close_aux_io_surf(hfile,hfiletype)
  ENDIF
  !
  !-------------------------------------------------------------------------------------
  !
  !*      5.     Total snow content
  !              ------------------
  !
  IF (nrank==npio) THEN
    !
    SELECT CASE (hsurf(1:3))
      CASE ('WWW')
        IF (osnow_ideal) THEN
          IF (jp<=1) ALLOCATE(pfield(ini,klayer,ipatch))
          pfield(:,:,jp) = tzsnow%WSNOW(:,1:klayer)
        ELSE
          IF (jp<=1) ALLOCATE(pfield(ini,1,ipatch))
          pfield(:,:,jp) = 0.
          DO jl=1,tzsnow%NLAYER
            pfield(:,1,jp) = pfield(:,1,jp) + tzsnow%WSNOW(:,jl)
          END DO 
          WHERE ( pfield(:,1,jp)>xundef ) pfield(:,1,jp)=xundef
        ENDIF
  !
  !-------------------------------------------------------------------------------------
  !
  !*      6.     Albedo
  !              ------
  !
    CASE ('ALB')
      IF (jp<=1) ALLOCATE(pfield(ini,1,ipatch))
      pfield(:,1,jp) = tzsnow%ALB(:)
  !
  !-------------------------------------------------------------------------------------
  !
  !*      7.     Total depth to snow grid
  !              ------------------------
  !
    CASE ('DEP')
      IF (osnow_ideal) THEN    
        IF (jp<=1) ALLOCATE(pfield(ini,klayer,ipatch))  
        pfield(:,:,jp) = tzsnow%WSNOW(:,1:klayer)/tzsnow%RHO(:,1:klayer)
        WHERE(tzsnow%WSNOW(:,1:klayer)==xundef) pfield(:,:,jp)=xundef
      ELSE     
        ALLOCATE(zd(ini))
        zd(:) = 0.0
        DO jl=1,tzsnow%NLAYER
          WHERE (tzsnow%WSNOW(:,jl)/=xundef)
            zd(:) = zd(:) + tzsnow%WSNOW(:,jl)/tzsnow%RHO(:,jl)
          ENDWHERE
        END DO
        IF (jp<=1) ALLOCATE(pfield(ini,1,ipatch))
        pfield(:,1,jp) = zd(:)
        DEALLOCATE(zd)
      ENDIF
  !
  !-------------------------------------------------------------------------------------
  !
  !*      8.     Density or heat profile
  !              -----------------------
  !
    CASE ('RHO','HEA','SG1','SG2','HIS','AGE')
  !
      SELECT CASE (tzsnow%SCHEME)
        CASE ('D95','1-L','EBA')
          IF (jp<=1) ALLOCATE(pfield(ini,1,ipatch))   
          !* computes output physical variable
          IF (hsurf(1:3)=='RHO') pfield(:,1,jp) = tzsnow%RHO(:,1)
          IF (hsurf(1:3)=='HEA') THEN
            IF (tzsnow%SCHEME=='D95'.OR.tzsnow%SCHEME=='EBA') pfield(:,1,jp) = xtt-2.
            IF (tzsnow%SCHEME=='1-L') pfield(:,1,jp) = tzsnow%T(:,1)
          END IF
          IF (hsurf(1:3)=='SG1') pfield(:,1,jp) = -20.0
          IF (hsurf(1:3)=='SG2') pfield(:,1,jp) = 80.0
          IF (hsurf(1:3)=='HIS') pfield(:,1,jp) = 0.0
          IF (hsurf(1:3)=='AGE') pfield(:,1,jp) = 3.0

        CASE ('3-L','CRO')
          ALLOCATE(zfield(ini,tzsnow%NLAYER))
          !* input physical variable
          IF (hsurf(1:3)=='RHO') zfield(:,:) = tzsnow%RHO (:,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 ( hsurf(1:3)=='HEA') THEN
            ALLOCATE(zheat(ini,tzsnow%NLAYER))
            zheat(:,:) = tzsnow%HEAT(:,1:tzsnow%NLAYER)
            CALL snow_heat_to_t_wliq(zheat,tzsnow%RHO,zfield)
            WHERE (zfield>xtt.AND.zfield/=xundef) zfield = xtt
            DEALLOCATE(zheat)
          ENDIF
          !
          IF (osnow_ideal) THEN 
            IF (jp<=1) ALLOCATE(pfield(ini,klayer,ipatch))                  
            pfield(:,:,jp) = zfield(:,:)
          ELSE
            !
            IF (jp<=1) ALLOCATE(pfield(ini,ngrid_level,ipatch))
            !* input snow layer thickness
            ALLOCATE(zdepth(ini,tzsnow%NLAYER))
            zdepth(:,:) = tzsnow%WSNOW(:,:)/tzsnow%RHO(:,:)
            !
            !* total depth
            ALLOCATE(zd(ini))
            zd(:) = 0.
            DO jl=1,tzsnow%NLAYER
              zd(:) = zd(:) + zdepth(:,jl)
            ENDDO
            !
            !* input normalized grid
            ALLOCATE(zgrid(ini,tzsnow%NLAYER))
            DO ji=1,ini
              IF(zd(ji)==0.0)THEN
                DO jl = 1,tzsnow%NLAYER
                  zgrid(ji,jl)=REAL(jl)/REAL(tzsnow%nlayer)
                ENDDO
              ELSE
                DO jl = 1,tzsnow%NLAYER
                  IF(jl==1)THEN
                    zgrid(ji,jl)=zdepth(ji,jl)/ zd(ji)
                  ELSE
                    zgrid(ji,jl) = zgrid(ji,jl-1) + zdepth(ji,jl)/zd(ji)
                  ENDIF
                ENDDO
              ENDIF
            ENDDO
            DEALLOCATE(zdepth)
            DEALLOCATE(zd)
            !    
            ! * interpolation of profile onto fine normalized snow grid
            CALL interp_grid_nat(zgrid(:,:),zfield(:,:),    &
                               xgrid_snow(:), pfield(:,:,jp))
            DEALLOCATE(zgrid)
          ENDIF
          DEALLOCATE(zfield)

        END SELECT
        !* put field form patch to all vegtypes    
    END SELECT
  !
  ENDIF
  !
  CALL dealloc_gr_snow(tzsnow)
  !
ENDDO
!
DEALLOCATE(imask_p)
!
!-------------------------------------------------------------------------------------
!
!*      9.     End of IO
!              ---------
!
IF (lhook) CALL dr_hook('PREP_SNOW_EXTERN',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
END SUBROUTINE prep_snow_extern