prep_hor_teb_garden_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_teb_garden_field (DTCO, UG, U, USS, GCP, IO, S, K, P, PEK, TG, TOP,  &
                                      HPROGRAM,HSURF,HATMFILE,HATMFILETYPE,HPGDFILE,HPGDFILETYPE,KPATCH,YDCTL)
!     #################################################################################
!
!!****  *PREP_HOR_TEB_GARDEN_FIELD* - reads, interpolates and prepares an ISBA field
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      P. Le Moigne 10/2005, Phasage Arome
!!      P. Le Moigne 03/2007, Ajout initialisation par ascllv
!!      B. Decharme  01/2009, Optional Arpege deep soil temperature initialization
!!      B. Decharme  03/2014, external init with FA files
!!                            new vertical interpol
!!      P. Marguinaud10/2014, Support for a 2-part PREP
!!------------------------------------------------------------------
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_sso_n, ONLY : sso_t
USE modd_grid_conf_proj_n, ONLY : grid_conf_proj_t
!
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_p_t, isba_pe_t
USE modd_isba_options_n, ONLY : isba_options_t
!
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_teb_option_n, ONLY : teb_options_t
!
USE modd_type_snow
USE modd_type_date_surf, ONLY : date_time
USE modd_grid_grib, ONLY : cinmodel
USE modd_surfex_mpi, ONLY : nrank, npio, ncomm, nproc
USE modd_prep,            ONLY : cingrid_type, cinterp_type, xzs_ls, linterp, cmask

USE modd_prep_teb_garden, ONLY : xgrid_soil, ngrid_level,                  &
                                 xwsnow_gd, xrsnow_gd, xtsnow_gd, xlwcsnow_gd, &
                                 xagesnow_gd, xasnow_gd, lsnow_ideal_gd
!
USE modd_isba_par,        ONLY : xwgmin
USE modd_data_cover_par,  ONLY : nvegtype
USE modd_surf_par,        ONLY : xundef
!
USE mode_prep_ctl, ONLY : prep_ctl, prep_ctl_can
!
USE modi_prep_grib_grid
USE modi_read_prep_teb_garden_conf
USE modi_read_prep_garden_snow
USE modi_prep_teb_garden_ascllv
USE modi_prep_teb_garden_grib
USE modi_prep_teb_garden_unif
USE modi_prep_teb_garden_buffer
USE modi_hor_interpol
USE modi_put_on_all_vegtypes
USE modi_vegtype_grid_to_patch_grid
USE modi_prep_hor_snow_fields
USE modi_get_luout
USE modi_prep_teb_garden_extern
USE modi_abor1_sfx
USE modi_allocate_gr_snow
!
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(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
TYPE(grid_conf_proj_t),INTENT(INOUT) :: GCP
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_p_t), INTENT(INOUT) :: P
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(grid_t), INTENT(INOUT) :: TG
TYPE(teb_options_t), INTENT(INOUT) :: TOP
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)  :: HATMFILE    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: HATMFILETYPE! type of the Atmospheric file
 CHARACTER(LEN=28),  INTENT(IN)  :: HPGDFILE    ! name of the Atmospheric file
 CHARACTER(LEN=6),   INTENT(IN)  :: HPGDFILETYPE! type of the Atmospheric file
!
INTEGER,            INTENT(IN)  :: KPATCH
!
!*      0.2    declarations of local variables
!
 CHARACTER(LEN=6)              :: YFILETYPE ! type of input file
 CHARACTER(LEN=28)             :: YFILE     ! name of file
 CHARACTER(LEN=6)              :: YFILEPGDTYPE ! type of input file
 CHARACTER(LEN=28)             :: YFILEPGD     ! name of file
 CHARACTER(LEN=6)              :: YFILETYPE_SNOW ! type of input file
 CHARACTER(LEN=28)             :: YFILE_SNOW     ! name of file
 CHARACTER(LEN=6)              :: YFILEPGDTYPE_SNOW ! type of input file
 CHARACTER(LEN=28)             :: YFILEPGD_SNOW     ! name of file 
REAL, POINTER,     DIMENSION(:,:,:) :: ZFIELDIN=>null()  ! field to interpolate horizontally
!
TYPE(nsurf_snow) :: TNPSNOW
!
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZFIELDOUTP ! field interpolated   horizontally
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZFIELDOUTV !
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: ZVEGTYPE_PATCH ! vegtype for each patch
REAL, ALLOCATABLE, DIMENSION(:,:)   :: ZW        ! work array (x, fine   soil grid)
REAL, ALLOCATABLE, DIMENSION(:)     :: ZSUM
REAL, ALLOCATABLE, DIMENSION(:,:)   :: ZF        ! work array (x, output soil grid)
REAL, ALLOCATABLE, DIMENSION(:,:)   :: ZDG       ! out T grid (x, output soil grid)
REAL, ALLOCATABLE, DIMENSION(:,:)   :: ZPATCH    ! work array for patches
REAL, ALLOCATABLE, DIMENSION(:)     :: ZSG1SNOW, ZSG2SNOW, ZHISTSNOW
INTEGER                             :: ILUOUT    ! output listing logical unit
!
type(date_time)                :: tztime_grib    ! current date and time
LOGICAL                             :: GUNIF     ! flag for prescribed uniform field
LOGICAL                             :: GUNIF_SNOW! flag for prescribed uniform field
INTEGER                             :: JVEGTYPE, JPATCH  ! loop on vegtypes
INTEGER                             :: JLAYER    ! loop on layers
INTEGER                             :: JI, INP, INL, INI
INTEGER                             :: IWORK     ! Work integer
INTEGER :: INFOMPI
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
!
!
!*      1.     Reading of input file name and type
!
IF (lhook) CALL dr_hook('PREP_HOR_TEB_GARDEN_FIELD',0,zhook_handle)
!
IF (.NOT. prep_ctl_can(ydctl)) THEN
  CALL abor1_sfx('PREP_HOR_TEB_GARDEN_FIELD: TWO STEP PREP NOT IMPLEMENTED')
ENDIF
!
 CALL get_luout(hprogram,iluout)
!
 CALL read_prep_teb_garden_conf(hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,&
                               hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,iluout,gunif)
!
cmask = 'TOWN  '
!
ini=SIZE(tg%XLAT)
!-------------------------------------------------------------------------------------
!
!*      2.     Snow variables case?
!
IF (hsurf=='SN_VEG ') THEN
  CALL read_prep_garden_snow(hprogram,pek%TSNOW%SCHEME,pek%TSNOW%NLAYER,yfile_snow,&
        yfiletype_snow,yfilepgd_snow,yfilepgdtype_snow,gunif_snow)
  !
  IF(.NOT.gunif_snow.AND.len_trim(yfile_snow)==0.AND.len_trim(yfiletype_snow)==0)THEN
    !IF(LEN_TRIM(YFILE)/=0.AND.LEN_TRIM(YFILETYPE)/=0)THEN
    IF (yfiletype=='GRIB') THEN
      yfile_snow        = yfile
      yfiletype_snow    = yfiletype
      yfilepgd_snow     = yfilepgd
      yfilepgdtype_snow = yfilepgdtype       
    ELSE          
      gunif_snow = .true.
      IF(all(xwsnow_gd==xundef)) xwsnow_gd = 0.0 
    ENDIF 
  ENDIF
  !
  ALLOCATE(zsg1snow(SIZE(xwsnow_gd)))
  ALLOCATE(zsg2snow(SIZE(xwsnow_gd)))
  ALLOCATE(zhistsnow(SIZE(xwsnow_gd)))
  !
  ALLOCATE(tnpsnow%AL(1))
  tnpsnow%AL(1)%SCHEME = pek%TSNOW%SCHEME
  tnpsnow%AL(1)%NLAYER = pek%TSNOW%NLAYER
  !
  CALL prep_hor_snow_fields(dtco, tg, u, gcp, hprogram,hsurf, &
                            yfile,yfiletype,                &
                            yfilepgd, yfilepgdtype,         &
                            iluout,gunif_snow,1, kpatch,    &
                            ini, tnpsnow, top%TTIME,        &
                            xwsnow_gd, xrsnow_gd, xtsnow_gd,&
                            xlwcsnow_gd, xasnow_gd,         &
                            lsnow_ideal_gd, zsg1snow,       &
                            zsg2snow, zhistsnow, xagesnow_gd, ydctl, &
                            pvegtype_patch=s%XVEGTYPE_PATCH, ppatch=s%XPATCH )
  !
  CALL allocate_gr_snow(pek%TSNOW,ini)
  pek%TSNOW%WSNOW = tnpsnow%AL(1)%WSNOW
  pek%TSNOW%RHO   = tnpsnow%AL(1)%RHO
  pek%TSNOW%ALB   = tnpsnow%AL(1)%ALB
  IF (pek%TSNOW%SCHEME/='D95') pek%TSNOW%HEAT = tnpsnow%AL(1)%HEAT
  IF (pek%TSNOW%SCHEME=='CRO'.OR.pek%TSNOW%SCHEME=='3-L') &
    pek%TSNOW%AGE = tnpsnow%AL(1)%AGE
  IF (pek%TSNOW%SCHEME=='CRO') THEN
    pek%TSNOW%GRAN1 = tnpsnow%AL(1)%GRAN1
    pek%TSNOW%GRAN2 = tnpsnow%AL(1)%GRAN2
    pek%TSNOW%HIST = tnpsnow%AL(1)%HIST
  ENDIF
  !
  CALL type_snow_init(tnpsnow%AL(1))
  DEALLOCATE(tnpsnow%AL)
  !
  DEALLOCATE(zsg1snow)
  DEALLOCATE(zsg2snow)
  DEALLOCATE(zhistsnow)  
  IF (lhook) CALL dr_hook('PREP_HOR_TEB_GARDEN_FIELD',1,zhook_handle)
  RETURN
END IF
!
!-------------------------------------------------------------------------------------
!
!*      3.     Reading of input  configuration (Grid and interpolation type)
!
IF (gunif) THEN
  CALL prep_teb_garden_unif(iluout,hsurf,zfieldin)
ELSE IF (yfiletype=='ASCLLV') THEN
  CALL prep_teb_garden_ascllv(dtco, ug, u, uss, hprogram,hsurf,iluout,zfieldin)
ELSE IF (yfiletype=='GRIB  ') THEN
  CALL prep_grib_grid(yfile,iluout,cinmodel,cingrid_type,cinterp_type,tztime_grib)            
   IF (nrank==npio) CALL prep_teb_garden_grib(hprogram,hsurf,yfile,iluout,zfieldin)        
ELSE IF (yfiletype=='MESONH' .OR. yfiletype=='ASCII ' .OR. yfiletype=='LFI   '&
        .OR.yfiletype=='FA    '.OR. yfiletype=='AROME '.OR.yfiletype=='NC    ') THEN
   CALL prep_teb_garden_extern(dtco, io, u, gcp, &
                               hprogram,hsurf,yfile,yfiletype,yfilepgd,yfilepgdtype,iluout,kpatch,zfieldin)
ELSE IF (yfiletype=='BUFFER') THEN
   CALL prep_teb_garden_buffer(hprogram,hsurf,iluout,zfieldin)
ELSE
   CALL abor1_sfx('PREP_HOR_TEB_GARDEN_FIELD: data file type not supported : '//yfiletype)
END IF
!
!-------------------------------------------------------------------------------------
!
!*      5.     Horizontal interpolation
!
IF (nrank==npio) THEN
  inl = SIZE(zfieldin,2)
  inp = SIZE(zfieldin,3)
ELSE
 IF (.NOT.ASSOCIATED(zfieldin)) 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(ini,inl,inp))
!
DO jpatch = 1, inp
  IF (inp==nvegtype) linterp = (s%XVEGTYPE(:,jpatch) > 0.)
  CALL hor_interpol(dtco, u, gcp, iluout,zfieldin(:,:,jpatch),zfieldoutp(:,:,jpatch))
  linterp = .true.
END DO
!
DEALLOCATE(zfieldin )
!
ALLOCATE(zw(ini,inl))
zw = 0.
!
IF (1/=inp) THEN
  !
  ALLOCATE(zfieldoutv(ini,inl,nvegtype))
  CALL put_on_all_vegtypes(ini,inl,inp,nvegtype,zfieldoutp,zfieldoutv)
  !
  ALLOCATE(zsum(ini))
  DO jlayer=1,SIZE(zw,2)
    zsum(:) = sum(s%XVEGTYPE(:,:),2,zfieldoutv(:,jlayer,:)/=xundef)
    DO jvegtype=1,nvegtype
      WHERE (zfieldoutv(:,jlayer,jvegtype)/=xundef) 
        zw(:,jlayer) = zw(:,jlayer) + s%XVEGTYPE(:,jvegtype) * zfieldoutv(:,jlayer,jvegtype) / zsum(:)
      END WHERE
    END DO
    DO ji=1,SIZE(zw,1)
      IF (all(zfieldoutv(ji,jlayer,:)==xundef)) zw(ji,jlayer) = xundef
    ENDDO
  END DO
  DEALLOCATE(zfieldoutv)
  DEALLOCATE(zsum)
  !
ELSE
  !
  zw(:,:) = zfieldoutp(:,:,1)
  !
ENDIF
!
DEALLOCATE(zfieldoutp)
!
!-------------------------------------------------------------------------------------
!
!*      7.     Return to historical variable
!
!
SELECT CASE (hsurf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WG     ') 
  ALLOCATE(zf(ini,io%NGROUND_LAYER))
  !
  !* interpolates on output levels
  CALL init_from_ref_grid(xgrid_soil,zw,p%XDG(:,:),zf)
  !
  !* retrieves soil water content from soil relative humidity
  ALLOCATE(pek%XWG(ini,io%NGROUND_LAYER))
  pek%XWG(:,:) = k%XWWILT + zf(:,:) * (k%XWFC-k%XWWILT)
  pek%XWG(:,:) = max(min(pek%XWG(:,:),k%XWSAT),xwgmin)
  !
  WHERE(zf(:,:)==xundef)pek%XWG(:,:)=xundef
  !
  DEALLOCATE(zf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WGI    ')
  ALLOCATE(zf(ini,io%NGROUND_LAYER))
  !
  !* interpolates on output levels
  CALL init_from_ref_grid(xgrid_soil,zw,p%XDG(:,:),zf) 
  !
  !* retrieves soil ice content from soil relative humidity
  ALLOCATE(pek%XWGI(ini,io%NGROUND_LAYER))
  pek%XWGI(:,:) = zf(:,:) * k%XWSAT
  pek%XWGI(:,:) = max(min(pek%XWGI(:,:),k%XWSAT),0.)
  !
  WHERE(zf(:,:)==xundef)pek%XWGI(:,:)=xundef
  !
  DEALLOCATE(zf)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('TG     ') 
  iwork=io%NGROUND_LAYER
  ALLOCATE(pek%XTG(ini,iwork))
  ALLOCATE(zdg(SIZE(p%XDG,1),iwork))
  IF (io%CISBA=='2-L'.OR.io%CISBA=='3-L') THEN
    zdg(:,1) = 0.01
    zdg(:,2) = 0.40   ! deep temperature for force-restore taken at 20cm
    IF(io%CISBA=='3-L') zdg(:,3) = 5.00   ! climatological temperature, usually not used
  ELSE
    !* diffusion method, the soil grid is the same as for humidity
    zdg(:,:) = p%XDG(:,:)
  END IF
  CALL init_from_ref_grid(xgrid_soil,zw,zdg,pek%XTG(:,:))
  DEALLOCATE(zdg)
  !
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('WR     ') 
  ALLOCATE(pek%XWR(ini))
  pek%XWR(:) = zw(:,1)
  !
  !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  !
 CASE('LAI    ') 
  !* LAI is updated only if present and pertinent (evolutive LAI) in input file

   WHERE (zw(:,1)/=xundef) pek%XLAI(:) = zw(:,1)
  !
END SELECT
!
DEALLOCATE(zw)
!-------------------------------------------------------------------------------------
!
!*      8.     Deallocations
!
!
IF (lhook) CALL dr_hook('PREP_HOR_TEB_GARDEN_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                                  :: JI, JL  ! loop counter
REAL, DIMENSION(SIZE(PT1),SIZE(PT1,2)) :: ZD1 ! input grid
!
INTEGER :: ILAYER1, ILAYER2
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',0,zhook_handle)
!
IF (SIZE(pt1,2)==3) THEN
!
!* 1. case with only 3 input levels (typically coming from 'UNIF')
!     -----------------------------
!
  IF (io%CISBA=='2-L' .OR. io%CISBA=='3-L') THEN
    !* Possible LTEMP_ARP case
    IF(SIZE(pt2,2)>3)THEN
       ilayer1=3
       ilayer2=SIZE(pt2,2)
    ELSE
       ilayer1=SIZE(pt2,2)
       ilayer2=0
    ENDIF
    !* historical 2L or 3L ISBA version
    pt2(:,1:ilayer1) = pt1(:,1:ilayer1) 
    !* Possible LTEMP_ARP case
    IF(ilayer2>0)THEN
       DO jl=ilayer1+1,ilayer2
         pt2(:,jl) = pt2(:,ilayer1)
       ENDDO
    ENDIF
!    
  ELSEIF(io%CISBA=='DIF')THEN
       !surface layer (generally 0.01m imposed)
       pt2(:,1) = pt1(:,1) 
       !deep layers
       DO jl=2,io%NGROUND_LAYER
          pt2(:,jl) = pt1(:,3)
       END DO
       !if root layers
       DO ji=1,SIZE(pt1,1)
          DO jl=2,io%NGROUND_LAYER
             IF(p%XROOTFRAC(ji,jl)<=1.0)THEN 
                pt2(ji,jl) = pt1(ji,2)
                EXIT
             ENDIF
          END DO
       END DO 
  END IF    
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ELSE
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 2. case with fine grid as input (general case)
!     ----------------------------
!
  DO jl=1,SIZE(pt1,2)
    zd1(:,jl) = pgrid1(jl)
  END DO
!
  CALL interp_grid_nat(zd1,pt1(:,:),pd2,pt2(:,:))
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ENDIF
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
IF (lhook) CALL dr_hook('INIT_FROM_REF_GRID',1,zhook_handle)
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END SUBROUTINE init_from_ref_grid
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_hor_teb_garden_field