prep_isba.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_isba (DTCO, UG, U, USS, GCP, SB, IG, IO, S, NK, NP, NPE,  &
                      HPROGRAM,HATMFILE,HATMFILETYPE,HPGDFILE,HPGDFILETYPE,YDCTL)
!     #################################################################################
!
!!****  *PREP_ISBA* - Prepares ISBA fields
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      Modified by P. Le Moigne (11/2004): AGS fields
!!      Modified by B. Decharme   (2008)  : Floodplains
!!      Modified by B. Decharme  (01/2009): Consistency with Arpege deep soil
!!                                          temperature
!!      Modified by B. Decharme  (03/2009): Consistency with Arpege permanent
!!                                          snow/ice treatment
!!      A.L. Gibelin 04/2009 : BIOMASS and RESP_BIOMASS arrays 
!!      A.L. Gibelin 06/2009 : Soil carbon variables for CNT option
!!      Modified by S. Riette    (06/2009): PREP_ISBA_CANOPY has no more arg.
!!      Modified by S. Riette    (04/2010): ecmwf ice content is computed during
!!                                          grib reading (no longer here)
!!      B. Decharme  (10/2012): coherence between soil temp and liquid/solid water with DIF
!!                              bug in biomass prognostic fields calculation
!!      B. Decharme  (06/2013): XPSNV_A for EBA snow scheme not allocated
!!      M. Lafaysse (04/2014) : LSNOW_PREP_PERM
!!      B. Decharme  (04/2013): Good computation for coherence between soil temp and 
!!                              liquid/solid water with DIF (results don't change)
!!                              if lglacier in input file, do not initialize again
!!      P. Samuelsson            (10/2014): MEB
!!      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_canopy_n, ONLY : canopy_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_nk_t, isba_k_t, isba_np_t, isba_p_t, &
                        isba_npe_t, isba_pe_t
!
USE modn_prep_isba
USE modn_prep_isba_snow, ONLY : lswemax, xswemax 
!
USE modd_data_cover_par, ONLY : nvt_snow
USE modd_snow_par,    ONLY : xemissn
USE modd_isba_par,    ONLY : xwgmin
USE modd_co2v_par,    ONLY : xanfminit
USE modd_surf_par,    ONLY : xundef
!
USE modd_read_namelist,  ONLY : lnam_read
USE modd_surf_atm,       ONLY : lvertshift
!          
USE modd_deepsoil,    ONLY : lphysdomc
USE modd_csts,        ONLY : xtt, xg, xlmtt
USE modd_prep,        ONLY : xzs_ls

USE modd_prep_snow,   ONLY : lsnow_prep_perm
!
USE mode_prep_ctl, ONLY : prep_ctl
!
USE modi_prep_hor_isba_field
USE modi_prep_ver_isba
USE modi_prep_output_grid
USE modi_get_luout
USE modi_prep_sbl
USE modi_vegtype_to_patch
!
USE modi_prep_perm_snow
USE modi_init_snow_lw
USE modi_averaged_albedo_emis_isba
USE modi_prep_hor_isba_cc_field
USE modi_soil_albedo
!
USE modi_clean_prep_output_grid
!
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(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(canopy_t), INTENT(INOUT) :: SB
TYPE(grid_t), INTENT(INOUT) :: IG
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_nk_t), INTENT(INOUT) :: NK
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE
!
type(prep_ctl),    INTENT(INOUT) :: ydctl
!
 CHARACTER(LEN=6),   INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 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
!
!*      0.2    declarations of local variables
!
TYPE(isba_k_t), POINTER :: KK
TYPE(isba_p_t), POINTER :: PK
TYPE(isba_pe_t), POINTER :: PEK
!
INTEGER :: ILUOUT, INI
INTEGER :: JP, JL, JJ
REAL    :: ZWORK, ZLOG, ZWTOT, ZMATPOT, ZWL
!
REAL,             DIMENSION(1)   :: ZSW_BANDS ! middle wavelength of each band
REAL,             DIMENSION(U%NSIZE_NATURE,IO%NPATCH) :: ZDIR_ALB, ZTG1  ! direct albedo for each band
REAL,             DIMENSION(U%NSIZE_NATURE,IO%NPATCH) :: ZSCA_ALB  ! diffuse albedo for each band
REAL,             DIMENSION(U%NSIZE_NATURE)   :: ZEMIS     ! emissivity
REAL,             DIMENSION(U%NSIZE_NATURE)   :: ZZENITH   ! solar zenithal angle
REAL,             DIMENSION(U%NSIZE_NATURE)   :: ZTSURF     ! surface effective temperature
!
LOGICAL         :: GPERMSNOW
LOGICAL         :: GTEMP2WGI
LOGICAL         :: GWG
LOGICAL         :: GWGI
LOGICAL         :: GTG
!
REAL            :: SMAX
!
INTEGER         :: ISIZE_LMEB_PATCH, ISNOW
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('PREP_ISBA',0,zhook_handle)
!
!*      1.     Default of configuration
!
gpermsnow = .true.
gwg       = .true.
gwgi      = .true.
gtg       = .true.
!
isize_lmeb_patch=count(io%LMEB_PATCH(:))
!
!*      1.1    Default
!
 CALL get_luout(hprogram,iluout)
!
 CALL prep_output_grid(ug%G, ig, u%NSIZE_FULL, iluout)
!
!-------------------------------------------------------------------------------------
!
!*      2.     Reading and horizontal interpolations
!
!
!*      2.0    Large scale orography
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'ZS     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
!
!*      2.1    Soil Water reservoirs
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WG     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl,gwg)
!
!*      2.2    Soil ice reservoirs
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WGI    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl,gwgi)
!
!*      2.3    Leaves interception water reservoir
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WR     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
!
!*      2.4    Temperature profile
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'TG     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl,gtg)
!
!*      2.5    Snow variables
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'SN_VEG ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl,gpermsnow)
!
!*      2.6    LAI
!
 CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'LAI    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
!
!*      2.7    GLACIER
!
IF(io%LGLACIER)THEN
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'ICE_STO',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
ENDIF
!
!*      2.8    Canopy vegetation temperature and interception reservoirs and air variables
!
IF(isize_lmeb_patch>0)THEN
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'TV     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                        hprogram,'TL     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WRL    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WRLI   ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'WRVN   ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'TC     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  CALL prep_hor_isba_field(dtco, ug, u, uss, gcp, ig, io, s, nk, np, npe, s%TTIME, &
                          hprogram,'QC     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
ENDIF
!
!*      7.     Isba-Ags prognostic fields
!
IF (io%CPHOTO == 'NIT' .OR. io%CPHOTO == 'NCB') THEN
  CALL prep_hor_isba_cc_field(dtco, u, gcp, SIZE(ig%XLAT), io, s, nk, np, npe,   &
                             hprogram,'BIOMASS ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
ENDIF
!
!*      8.     Isba-CC prognostic fields
!
IF (io%CPHOTO/='NON' .AND. io%CRESPSL == 'CNT') THEN
  !
  !*      8.1    Litter
  !
    CALL prep_hor_isba_cc_field(dtco, u, gcp, SIZE(ig%XLAT), io, s, nk, np, npe,  &
                                 hprogram,'LITTER  ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  !
  !*      8.2    Soil carbon
  !
    CALL prep_hor_isba_cc_field(dtco, u, gcp, SIZE(ig%XLAT), io, s, nk, np, npe,  &
                               hprogram,'SOILCARB',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
  !
  !*      8.2    lignin
  !
   CALL prep_hor_isba_cc_field(dtco, u, gcp, SIZE(ig%XLAT), io, s, nk, np, npe,  &
                               hprogram,'LIGNIN  ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,ydctl)
ENDIF
!
 CALL clean_prep_output_grid

IF (ydctl%LPART6) THEN
!-------------------------------------------------------------------------------------
!
!*      3.    Physical limitation: 
!
  DO jp=1,io%NPATCH
    pek => npe%AL(jp)
    kk => nk%AL(jp)
  !
  ! No ice for force restore third layer:
    IF (io%CISBA == '3-L') THEN
    !
      WHERE(pek%XWG(:,3) /= xundef)
        pek%XWG(:,3) = min(pek%XWG(:,3)+ pek%XWGI(:,3),kk%XWSAT(:,3))
        pek%XWGI(:,3) = 0.
      END WHERE
    ENDIF
  !
  ! Total water content should not exceed saturation:
    WHERE(pek%XWG(:,:) /= xundef .AND. (pek%XWG(:,:) + pek%XWGI(:,:)) > kk%XWSAT(:,:) )
      pek%XWGI(:,:) = kk%XWSAT(:,:) - pek%XWG(:,:)
    END WHERE

  ENDDO
!
!-------------------------------------------------------------------------------------
!
!*      3.     Vertical interpolations of all variables
!
  IF(lvertshift)THEN
    CALL prep_ver_isba(io, npe, s%XZS, np)
  ENDIF
!
  DEALLOCATE(xzs_ls)
!-------------------------------------------------------------------------------------
!
!*      4.     Treatment of permanent snow
!
  IF (gpermsnow.AND.lsnow_prep_perm) THEN
    isnow = vegtype_to_patch(nvt_snow,io%NPATCH)
    CALL prep_perm_snow(io, nk%AL(isnow), np%AL(isnow), npe%AL(isnow))
  ENDIF
!
  DO jp = 1,io%NPATCH
    CALL init_snow_lw(xemissn,npe%AL(jp)%TSNOW)
  ENDDO
!
  IF (lphysdomc) THEN
    DO jp = 1,io%NPATCH
      npe%AL(jp)%TSNOW%WSNOW(:,:)=0.
    ENDDO
  ENDIF 
!------------------------------------------------------------------------------------- 
! 
!*      4.b     Possibility for setting an upper limit on the initial snow water equivalent field 
!
  IF (lswemax) THEN 
    DO jp = 1,io%NPATCH
      smax = maxval(npe%AL(jp)%TSNOW%WSNOW(:,:)) 
      WRITE(*,*) ' MAX(Snow content (kg/m2)): ', smax 
      WRITE(*,*) ' Set MAX to', xswemax, '(kg/m2)' 
      npe%AL(jp)%TSNOW%WSNOW(:,:) = min(npe%AL(jp)%TSNOW%WSNOW(:,:),xswemax) 
      smax = maxval(npe%AL(jp)%TSNOW%WSNOW(:,:)) 
      WRITE(*,*) ' MAX(Snow content (kg/m2)): ', smax 
    ENDDO
  ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      5.     coherence between soil temperature and liquid/solid water
!
  gtemp2wgi=(gwg.OR.gwgi.OR.gtg)
!
  IF (io%CISBA == 'DIF'.AND.gtemp2wgi) THEN
    DO jp=1,io%NPATCH
      pek => npe%AL(jp)
      pk => np%AL(jp)
      kk => nk%AL(jp)
      !
      DO jl=1,io%NGROUND_LAYER
        !
        DO jj=1,pk%NSIZE_P
          !
          IF(pek%XWG(jj,jl)/=xundef)THEN
!     
!         total soil moisture
            zwtot = pek%XWG(jj,jl)+pek%XWGI(jj,jl)
            zwtot = min(zwtot,kk%XWSAT(jj,jl))
!                
!           total matric potential
!           psi=mpotsat*(w/wsat)**(-bcoef)
            zwork   = zwtot/kk%XWSAT(jj,jl)
            zlog    = kk%XBCOEF(jj,jl)*log(zwork)
            zmatpot = kk%XMPOTSAT(jj,jl)*exp(-zlog)
!
!         soil liquid water content computation
!         w=wsat*(psi/mpotsat)**(-1/bcoef)
            zmatpot       = min(kk%XMPOTSAT(jj,jl),xlmtt*(pek%XTG(jj,jl)-xtt)/(xg*pek%XTG(jj,jl)))
            zwork         = max(1.0,zmatpot/kk%XMPOTSAT(jj,jl))
            zlog          = log(zwork)
            zwl           = kk%XWSAT(jj,jl)*exp(-zlog/kk%XBCOEF(jj,jl))
            zwl           = max(zwl,xwgmin)
            pek%XWG(jj,jl) = min(zwl,zwtot )
!        
!         soil ice computation    
            pek%XWGI(jj,jl) = max(0.0,zwtot-pek%XWG(jj,jl))
! 
!         supress numerical artefact
            IF(pek%XTG(jj,jl)>=xtt)THEN
              pek%XWG (jj,jl) = min(pek%XWG(jj,jl)+pek%XWGI(jj,jl),kk%XWSAT(jj,jl))
              pek%XWGI(jj,jl) = 0.0
            ENDIF
!
          ENDIF
        ENDDO        
      ENDDO        
    ENDDO
  ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      6.     Half prognostic fields
!              The only variable used from the AVERAGED_ALBEDO_EMIS_ISBA call
!              is XTSRAD_NAT. All other variables are treated as dummies.
!
  DO jp = 1,io%NPATCH
  !
    pek => npe%AL(jp)
    pk => np%AL(jp)
    kk => nk%AL(jp)
    !
    ALLOCATE(pek%XRESA(pk%NSIZE_P))
    pek%XRESA = 100.
    !
    ALLOCATE(pek%XALBNIR(pk%NSIZE_P))
    ALLOCATE(pek%XALBVIS(pk%NSIZE_P))
    ALLOCATE(pek%XALBUV(pk%NSIZE_P))
    pek%XALBNIR = 0.0
    pek%XALBVIS = 0.0
    pek%XALBUV = 0.0
    !
    ALLOCATE(pek%XALBNIR_SOIL(pk%NSIZE_P))
    ALLOCATE(pek%XALBVIS_SOIL(pk%NSIZE_P))
    ALLOCATE(pek%XALBUV_SOIL (pk%NSIZE_P))
    CALL soil_albedo (io%CALBEDO, kk%XWSAT(:,1),pek%XWG(:,1), kk, pek, "ALL" )
    !
    ALLOCATE(pek%XPSN   (pk%NSIZE_P))
    ALLOCATE(pek%XPSNG  (pk%NSIZE_P))
    ALLOCATE(pek%XPSNV  (pk%NSIZE_P))
    ALLOCATE(pek%XPSNV_A(pk%NSIZE_P))
    pek%XPSN    = 0.0
    pek%XPSNG   = 0.0
    pek%XPSNV   = 0.0
    pek%XPSNV_A = 0.0
    ALLOCATE(kk%XDIR_ALB_WITH_SNOW(pk%NSIZE_P,1))
    ALLOCATE(kk%XSCA_ALB_WITH_SNOW(pk%NSIZE_P,1))
    !
    ztg1(1:pk%NSIZE_P,jp) = pek%XTG(:,1)
    !
  ENDDO
  !
  ALLOCATE(s%XTSRAD_NAT(u%NSIZE_NATURE))
  zzenith(:)=0.
  zsw_bands(:)=0.
  !
  CALL averaged_albedo_emis_isba(io, s, nk, np, npe, &
                                 zzenith, ztg1, zsw_bands, zdir_alb, zsca_alb,   &
                                 zemis, s%XTSRAD_NAT, ztsurf              )
  !
  DO jp = 1,io%NPATCH
    DEALLOCATE(npe%AL(jp)%XPSN)
    DEALLOCATE(npe%AL(jp)%XPSNG)
    DEALLOCATE(npe%AL(jp)%XPSNV)
    DEALLOCATE(npe%AL(jp)%XPSNV_A)
    DEALLOCATE(nk%AL(jp)%XDIR_ALB_WITH_SNOW)
    DEALLOCATE(nk%AL(jp)%XSCA_ALB_WITH_SNOW)
  ENDDO
  !
  !-------------------------------------------------------------------------------------
  !
  !*      7.     Isba-Ags prognostic fields
  !
  IF (io%CPHOTO /= 'NON') THEN
  !
    DO jp = 1,io%NPATCH
      !
      pek => npe%AL(jp)
      pk => np%AL(jp)
      !
      ALLOCATE(pek%XAN(pk%NSIZE_P))
      pek%XAN = 0.
      !
      ALLOCATE(pek%XANDAY(pk%NSIZE_P))
      pek%XANDAY = 0.
      !
      ALLOCATE(pek%XANFM(pk%NSIZE_P))
      pek%XANFM = xanfminit
      !
      ALLOCATE(pek%XLE(pk%NSIZE_P))
      pek%XLE = 0.
      !
      ALLOCATE(pek%XRESP_BIOMASS(pk%NSIZE_P,io%NNBIOMASS))
      pek%XRESP_BIOMASS(:,:) = 0.
      !
    ENDDO
    !
  ENDIF
  !
  IF (io%CPHOTO == 'AST') THEN
  !
    DO jp = 1,io%NPATCH
      !
      pek => npe%AL(jp)
      pk => np%AL(jp)
      !
      ALLOCATE(pek%XBIOMASS(pk%NSIZE_P,io%NNBIOMASS))
      pek%XBIOMASS(:,:) = 0.
      !
    ENDDO
  !
  ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      10.     Preparation of canopy air variables
!
!
  io%LCANOPY = lisba_canopy
  IF (io%LCANOPY) CALL prep_sbl(ig%NDIM, sb)
!
ENDIF
!
IF (lhook) CALL dr_hook('PREP_ISBA',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_isba