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, ICP, IG, I, UG, U, USS, &
                      hprogram,hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!     #################################################################################
!
!!****  *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
!!------------------------------------------------------------------
!
!
!
!
!
!
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_isba_canopy_n, ONLY : isba_canopy_t
USE modd_isba_grid_n, ONLY : isba_grid_t
USE modd_isba_n, ONLY : isba_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modi_prep_hor_isba_field
USE modi_prep_ver_isba
USE modi_prep_output_grid
USE modi_get_luout
USE modi_prep_isba_canopy
!
USE modd_read_namelist,  ONLY : lnam_read
USE modd_surf_atm,       ONLY : lvertshift
USE modd_data_cover_par, ONLY : nvt_snow
!
!                           
USE modd_deepsoil,    ONLY : lphysdomc
USE modd_csts,        ONLY : xtt, xg, xlmtt
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_prep,        ONLY : xzs_ls

USE modd_prep_snow,   ONLY : lsnow_prep_perm
!
USE modn_prep_isba
USE modn_prep_isba_snow, ONLY : lswemax, xswemax 
!
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 yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_clean_prep_output_grid
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(isba_canopy_t), INTENT(INOUT) :: icp
TYPE(isba_grid_t), INTENT(INOUT) :: ig
TYPE(isba_t), INTENT(INOUT) :: i
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
TYPE(surf_atm_sso_t), INTENT(INOUT) :: uss
!
 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
!
INTEGER :: iluout, ini
INTEGER :: jp, jl, jj
INTEGER :: isnow          ! patch number where permanent snow is
REAL    :: zwork, zlog, zwtot, zmatpot, zwl
!
REAL,             DIMENSION(1)   :: zsw_bands ! middle wavelength of each band
REAL,             DIMENSION(SIZE(I%XLAI,1),SIZE(I%XLAI,2)) :: zdir_alb  ! direct albedo for each band
REAL,             DIMENSION(SIZE(I%XLAI,1),SIZE(I%XLAI,2)) :: zsca_alb  ! diffuse albedo for each band
REAL,             DIMENSION(SIZE(I%XLAI,1))   :: zemis     ! emissivity
REAL,             DIMENSION(SIZE(I%XLAI,1))   :: zzenith   ! solar zenithal angle
REAL,             DIMENSION(SIZE(I%XLAI,1))   :: ztsurf     ! surface effective temperature
!
LOGICAL         :: gpermsnow
LOGICAL         :: gtemp2wgi
LOGICAL         :: gwg
LOGICAL         :: gwgi
LOGICAL         :: gtg
!
REAL            :: smax
!
INTEGER         :: isize_lmeb_patch
!
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(i%LMEB_PATCH(:))
!
!*      1.1    Default
!
 CALL get_luout(hprogram,iluout)
!
 CALL prep_output_grid(ug, u, &
                       iluout,ig%CGRID,ig%XGRID_PAR,ig%XLAT,ig%XLON)
!
!-------------------------------------------------------------------------------------
!
!*      2.     Reading and horizontal interpolations
!
!
!*      2.0    Large scale orography
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'ZS     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
!*      2.1    Soil Water reservoirs
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WG     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,gwg)
!
!*      2.2    Soil ice reservoirs
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WGI    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,gwgi)
!
!*      2.3    Leaves interception water reservoir
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WR     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
!*      2.4    Temperature profile
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'TG     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,gtg)
!
!*      2.5    Snow variables
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'SN_VEG ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype,gpermsnow)
!
!*      2.6    LAI
!
 CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'LAI    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
!*      2.7    GLACIER
!
IF(i%LGLACIER)THEN
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'ICE_STO',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
ENDIF
!
!*      2.8    Canopy vegetation temperature and interception reservoirs and air variables
!
IF(isize_lmeb_patch>0)THEN
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'TV     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                        hprogram,'TL     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WRL    ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WRLI   ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'WRVN   ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'TC     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
  CALL prep_hor_isba_field(dtco, ig, i, ug, u, uss, &
                          hprogram,'QC     ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      3.    Physical limitation: 
!
! No ice for force restore third layer:
IF (i%CISBA == '3-L') THEN
   DO jp=1,i%NPATCH
      WHERE(i%XWG(:,3,jp) /= xundef)
        i%XWG(:,3,jp)  = min(i%XWG(:,3,jp)+i%XWGI(:,3,jp),i%XWSAT(:,3))
        i%XWGI(:,3,jp) = 0.
      END WHERE
   ENDDO
ENDIF
!
! Total water content should not exceed saturation:
DO jp=1,i%NPATCH
   WHERE(i%XWG(:,:,jp) /= xundef .AND. (i%XWG(:,:,jp) + i%XWGI(:,:,jp)) > i%XWSAT(:,:) )
      i%XWGI(:,:,jp) = i%XWSAT(:,:) - i%XWG(:,:,jp)
   END WHERE
ENDDO
!
!-------------------------------------------------------------------------------------
!
!*      3.     Vertical interpolations of all variables
!
IF(lvertshift)THEN
  CALL prep_ver_isba(i)
ENDIF
!
DEALLOCATE(xzs_ls)
!-------------------------------------------------------------------------------------
!
!*      4.     Treatment of permanent snow
!
IF (gpermsnow.AND.lsnow_prep_perm) THEN
  isnow = vegtype_to_patch(nvt_snow,i%NPATCH)
  CALL prep_perm_snow(i, &
                      i%TSNOW,i%XTG(:,:,isnow),i%XVEGTYPE_PATCH(:,:,isnow),isnow)
ENDIF
!
 CALL init_snow_lw(xemissn,i%TSNOW)
!
IF (lphysdomc) THEN
   i%TSNOW%WSNOW(:,:,:)=0.
ENDIF 
!------------------------------------------------------------------------------------- 
! 
!*      4.b     Possibility for setting an upper limit on the initial snow water equivalent field 
IF (lswemax) THEN 
  smax = maxval(i%TSNOW%WSNOW(:,:,:)) 
  WRITE(*,*) ' MAX(Snow content (kg/m2)): ', smax 
  WRITE(*,*) ' Set MAX to', xswemax, '(kg/m2)' 
  i%TSNOW%WSNOW(:,:,:) = min(i%TSNOW%WSNOW(:,:,:),xswemax) 
  smax = maxval(i%TSNOW%WSNOW(:,:,:)) 
  WRITE(*,*) ' MAX(Snow content (kg/m2)): ', smax 
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      5.     coherence between soil temperature and liquid/solid water
!
gtemp2wgi=(gwg.OR.gwgi.OR.gtg)
!
IF (i%CISBA == 'DIF'.AND.gtemp2wgi) THEN
   ini=SIZE(i%XWSAT,1)
   DO jp=1,i%NPATCH
      DO jl=1,i%NGROUND_LAYER
         DO jj=1,ini
            IF(i%XWG(jj,jl,jp)/=xundef)THEN
!     
!             total soil moisture
              zwtot = i%XWG(jj,jl,jp)+i%XWGI(jj,jl,jp)
              zwtot = min(zwtot,i%XWSAT(jj,jl))
!              
!             total matric potential
!             psi=mpotsat*(w/wsat)**(-bcoef)
              zwork   = zwtot/i%XWSAT(jj,jl)
              zlog    = i%XBCOEF(jj,jl)*log(zwork)
              zmatpot = i%XMPOTSAT(jj,jl)*exp(-zlog)
!
!             soil liquid water content computation
!             w=wsat*(psi/mpotsat)**(-1/bcoef)
              zmatpot       = min(i%XMPOTSAT(jj,jl),xlmtt*(i%XTG(jj,jl,jp)-xtt)/(xg*i%XTG(jj,jl,jp)))
              zwork         = max(1.0,zmatpot/i%XMPOTSAT(jj,jl))
              zlog          = log(zwork)
              zwl           = i%XWSAT(jj,jl)*exp(-zlog/i%XBCOEF(jj,jl))
              zwl           = max(zwl,xwgmin)
              i%XWG(jj,jl,jp) = min(zwl,zwtot )
!        
!             soil ice computation    
              i%XWGI(jj,jl,jp) = max(0.0,zwtot-i%XWG(jj,jl,jp))
! 
!             supress numerical artefact
              IF(i%XTG(jj,jl,jp)>=xtt)THEN
                i%XWG (jj,jl,jp) = min(i%XWG(jj,jl,jp)+i%XWGI(jj,jl,jp),i%XWSAT(jj,jl))
                i%XWGI(jj,jl,jp) = 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.
!
ALLOCATE(i%XRESA(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
i%XRESA = 100.
!
ALLOCATE(i%XTSRAD_NAT(SIZE(i%XLAI,1)))
zzenith(:)=0.
zsw_bands(:)=0.
!
ALLOCATE(i%XALBNIR(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XALBVIS(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XALBUV(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
i%XALBNIR = 0.0
i%XALBVIS = 0.0
i%XALBUV = 0.0
!
ALLOCATE(i%XALBNIR_SOIL(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XALBVIS_SOIL(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XALBUV_SOIL(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
 CALL soil_albedo(i%CALBEDO, i%XWSAT(:,1),i%XWG(:,1,:),     &
                    i%XALBVIS_DRY,i%XALBNIR_DRY,i%XALBUV_DRY, &
                    i%XALBVIS_WET,i%XALBNIR_WET,i%XALBUV_WET, &
                    i%XALBVIS_SOIL,i%XALBNIR_SOIL,i%XALBUV_SOIL )
!
ALLOCATE(i%XPSN   (SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XPSNG  (SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XPSNV  (SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
ALLOCATE(i%XPSNV_A(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
i%XPSN    = 0.0
i%XPSNG   = 0.0
i%XPSNV   = 0.0
i%XPSNV_A = 0.0
ALLOCATE(i%XDIR_ALB_WITH_SNOW(SIZE(i%XLAI,1),1,SIZE(i%XLAI,2)))
ALLOCATE(i%XSCA_ALB_WITH_SNOW(SIZE(i%XLAI,1),1,SIZE(i%XLAI,2)))
i%XDIR_ALB_WITH_SNOW = 0.0
i%XSCA_ALB_WITH_SNOW = 0.0
 CALL averaged_albedo_emis_isba(i, &
                               .false., i%CALBEDO, zzenith,                &
                                 i%XVEG,i%XZ0,i%XLAI,                          &
                                 i%LMEB_PATCH,i%XGNDLITTER,i%XZ0LITTER,i%XLAIGV, &
                                 i%XH_VEG, i%XTV,               &
                                 i%XTG(:,1,:),i%XPATCH, zsw_bands,           &
                                 i%XALBNIR_VEG,i%XALBVIS_VEG,i%XALBUV_VEG,     &
                                 i%XALBNIR_SOIL,i%XALBVIS_SOIL,i%XALBUV_SOIL,  &
                                 i%XEMIS,                                  &
                                 i%TSNOW,                                  &
                                 i%XALBNIR,i%XALBVIS,i%XALBUV,                 &
                                 zdir_alb, zsca_alb,                     &
                                 zemis,i%XTSRAD_NAT,ztsurf                 )
DEALLOCATE(i%XPSN)
DEALLOCATE(i%XPSNG)
DEALLOCATE(i%XPSNV)
DEALLOCATE(i%XPSNV_A)
DEALLOCATE(i%XDIR_ALB_WITH_SNOW)
DEALLOCATE(i%XSCA_ALB_WITH_SNOW)
!
!-------------------------------------------------------------------------------------
!
!*      7.     Isba-Ags prognostic fields
!
IF (i%CPHOTO /= 'NON') THEN
!
   ALLOCATE(i%XAN(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
   i%XAN = 0.
!
   ALLOCATE(i%XANDAY(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
   i%XANDAY = 0.
!
   ALLOCATE(i%XANFM(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
   i%XANFM = xanfminit
!
   ALLOCATE(i%XLE(SIZE(i%XLAI,1),SIZE(i%XLAI,2)))
   i%XLE = 0.
!
   ALLOCATE(i%XRESP_BIOMASS(SIZE(i%XLAI,1),i%NNBIOMASS,SIZE(i%XLAI,2)))
   i%XRESP_BIOMASS(:,:,:) = 0.
!
ENDIF
!
IF (i%CPHOTO == 'AGS' .OR. i%CPHOTO == 'AST') THEN
!
   ALLOCATE(i%XBIOMASS(SIZE(i%XLAI,1),i%NNBIOMASS,SIZE(i%XLAI,2)))
   i%XBIOMASS(:,:,:) = 0.
!
ELSEIF (i%CPHOTO == 'LAI' .OR. i%CPHOTO == 'LST') THEN
!
   ALLOCATE(i%XBIOMASS(SIZE(i%XLAI,1),i%NNBIOMASS,SIZE(i%XLAI,2)))
   WHERE(i%XLAI(:,:)/=xundef)
         i%XBIOMASS(:,1,:) = i%XLAI(:,:) * i%XBSLAI(:,:)
   ELSEWHERE
         i%XBIOMASS(:,1,:) = 0.0
   ENDWHERE
!
ELSEIF (i%CPHOTO == 'NIT' .OR. i%CPHOTO == 'NCB') THEN
!
   CALL prep_hor_isba_cc_field(dtco, u, &
                               ig, i, &
                               hprogram,'BIOMASS ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)   
!
ENDIF
!
!-------------------------------------------------------------------------------------
!
!*      8.     Isba-CC prognostic fields
!
IF (i%CRESPSL == 'CNT') THEN
!
!*      8.1    Litter
!
 CALL prep_hor_isba_cc_field(dtco, u, &
                               ig, i, &
                               hprogram,'LITTER  ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
!*      8.2    Soil carbon
!
 CALL prep_hor_isba_cc_field(dtco, u, &
                               ig, i, &
                               hprogram,'SOILCARB',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
!*      8.2    lignin
!
 CALL prep_hor_isba_cc_field(dtco, u, &
                               ig, i, &
                               hprogram,'LIGNIN  ',hatmfile,hatmfiletype,hpgdfile,hpgdfiletype)
!
ENDIF
!
!-------------------------------------------------------------------------------------
 CALL clean_prep_output_grid
!-------------------------------------------------------------------------------------
!
!*      10.     Preparation of canopy air variables
!
!
i%LCANOPY = lisba_canopy
IF (i%LCANOPY) CALL prep_isba_canopy(icp, ig)
!
IF (lhook) CALL dr_hook('PREP_ISBA',1,zhook_handle)
!
!-------------------------------------------------------------------------------------
!
END SUBROUTINE prep_isba