update_rad_isban.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 update_rad_isba_n (IO, S, KK, PK, PEK, KPATCH, PZENITH, PSW_BANDS, &
                              PDIR_ALB_WITH_SNOW,PSCA_ALB_WITH_SNOW, PEMIST,  &
                              PDIR_SW, PSCA_SW     )
!     ####################################################################
!
!!****  *UPDATE_RAD_ISBA_n * - Calculate snow/flood fraction, dir/dif albedo
!!                             and emissivity at t+1 in order to close the 
!!                             energy budget between the atmospheric model 
!!                             and surfex  
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!,ZEMIST,PEMIST,ZPUT0)!
!!    AUTHOR
!!    ------
!!     B. Decharme 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    09/2009
!!      P. Samuelsson 02/2012 MEB
!!      A. Boone      03/2015 MEB-use TR_ML scheme for SW radiation
!!------------------------------------------------------------------
!
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY: isba_s_t, isba_k_t, isba_p_t, isba_pe_t
!
USE modd_type_snow
!
USE modd_csts,      ONLY : xtt
USE modd_surf_par,  ONLY : xundef
USE modd_snow_par,  ONLY : xrhosmin_es,xrhosmax_es,xsnowdmin,xemissn
USE modd_water_par, ONLY : xalbsca_wat, xemiswat, xalbwatice, xemiswatice 
USE modd_meb_par,   ONLY : xsw_wght_vis, xsw_wght_nir
!
USE mode_surf_flood_frac
USE mode_surf_snow_frac      
USE mode_meb,       ONLY : meb_shield_factor, mebpalphan
!
USE modi_albedo_ta96
USE modi_albedo_from_nir_vis
USE modi_pack_same_rank
USE modi_unpack_same_rank
USE modi_isba_snow_frac
USE modi_isba_emis_meb
USE modi_radiative_transfert
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
!
INTEGER, INTENT(IN) :: KPATCH
!
REAL, DIMENSION(:),     INTENT(IN)   :: PZENITH   ! Zenithal angle at t+1
REAL, DIMENSION(:),     INTENT(IN)   :: PSW_BANDS ! mean wavelength of each shortwave band (m)
!
REAL, DIMENSION(:,:), INTENT(OUT)  :: PDIR_ALB_WITH_SNOW ! Total direct albedo at t+1
REAL, DIMENSION(:,:), INTENT(OUT)  :: PSCA_ALB_WITH_SNOW ! Total diffuse albedo at t+1
REAL, DIMENSION(:),   INTENT(OUT)  :: PEMIST             ! Total emissivity at t+1
!
REAL, DIMENSION(:,:),   INTENT(IN), OPTIONAL   :: PDIR_SW   ! direct  solar radiation (on horizontal surf.)
REAL, DIMENSION(:,:),   INTENT(IN), OPTIONAL   :: PSCA_SW   ! diffuse solar radiation (on horizontal surf.)
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(PK%NSIZE_P) :: ZVEG
REAL, DIMENSION(PK%NSIZE_P) :: ZPSNV_A
REAL, DIMENSION(PK%NSIZE_P) :: ZALBF_DIR
REAL, DIMENSION(PK%NSIZE_P) :: ZALBF_SCA
REAL, DIMENSION(PK%NSIZE_P) :: ZEMISF   
REAL, DIMENSION(PK%NSIZE_P) :: ZFF
REAL, DIMENSION(PK%NSIZE_P) :: ZALBNIR_WITH_SNOW
REAL, DIMENSION(PK%NSIZE_P) :: ZALBVIS_WITH_SNOW
REAL, DIMENSION(PK%NSIZE_P) :: ZALBUV_WITH_SNOW
REAL, DIMENSION(PK%NSIZE_P) :: ZZENITH
REAL, DIMENSION(PK%NSIZE_P) :: ZSNOWDEPTH, ZPALPHAN
REAL, DIMENSION(PK%NSIZE_P) :: ZSWUP
REAL, DIMENSION(PK%NSIZE_P) :: ZGLOBAL_SW
REAL, DIMENSION(PK%NSIZE_P) :: ZALBT, ZEMIST
REAL, DIMENSION(PK%NSIZE_P) :: ZPSNA, ZSIGMA_F, ZSIGMA_FN, ZEMISSN
REAL, DIMENSION(PK%NSIZE_P,SIZE(PSW_BANDS)) :: ZDIR_SW, ZSCA_SW
REAL, DIMENSION(PK%NSIZE_P)            :: ZLAIN, ZALBVIS_TSOIL, ZALBNIR_TSOIL
REAL, DIMENSION(PK%NSIZE_P)            :: ZFAPIR, ZFAPAR, ZFAPIR_BS, ZFAPAR_BS
REAL, DIMENSION(PK%NSIZE_P,SIZE(S%XABC)) :: ZIACAN_SUNLIT, ZIACAN_SHADE, ZFRAC_SUN, ZIACAN
LOGICAL, DIMENSION(PK%NSIZE_P)         :: GSHADE
!
REAL, PARAMETER :: ZPUT0 = 0.0
INTEGER :: ISWB
INTEGER :: JSWB
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------------
! Initialization
!-------------------------------------------------------------------------------------

IF (lhook) CALL dr_hook('UPDATE_RAD_ISBA_N',0,zhook_handle)
iswb   = SIZE(psw_bands)
!
!-------------------------------------------------------------------------------------
!
zveg(:) = pek%XVEG(:)
!
IF(io%LMEB_PATCH(kpatch).OR.io%LFLOOD)THEN
  !
  CALL pack_same_rank(pk%NR_P,pzenith(:),zzenith(:))
  !
  IF(io%LMEB_PATCH(kpatch))THEN
    !
    zveg(:)=0.   ! Set veg=0 for MEB to get bare soil conditions for snow cover and
    !            ! flood fraction
    !
    IF(PRESENT(pdir_sw))THEN
      !
      CALL pack_same_rank(pk%NR_P,pdir_sw(:,:), zdir_sw(:,:))
      CALL pack_same_rank(pk%NR_P,psca_sw(:,:), zsca_sw(:,:))  
      !
    ENDIF
    !
  ENDIF
  !
ENDIF
!   
!-------------------------------------------------------------------------------
!
 CALL isba_snow_frac(pek%TSNOW%SCHEME, pek%TSNOW%WSNOW, pek%TSNOW%RHO, pek%TSNOW%ALB, &
                     zveg, pek%XLAI, pek%XZ0, pek%XPSN, zpsnv_a, pek%XPSNG, pek%XPSNV )  
!
IF ( pek%TSNOW%SCHEME=='EBA' ) pek%XPSNV_A(:) = zpsnv_a(:)
!
!-------------------------------------------------------------------------------
!
! Flood fractions and properties
!
IF(io%LFLOOD)THEN  
  !
  kk%XFFG(:) = flood_frac_ground(pek%XPSNG,kk%XFFLOOD)
  kk%XFFV(:) = flood_frac_veg(pek%XLAI,pek%XPSNV,kk%XFFLOOD)
  kk%XFF (:) = flood_frac_nat(pek%XVEG,kk%XFFG,kk%XFFV,kk%XFFLOOD)
  !
  WHERE(kk%XFFLOOD(:)==0.0)
    zalbf_dir(:) = xundef
    zalbf_sca(:) = xundef
    kk%XALBF   (:) = xundef
    kk%XEMISF  (:) = xundef
    kk%XFFROZEN(:) = 0.0
  ELSEWHERE
    WHERE(pek%XTG(:,1)>=xtt)
      zalbf_dir(:) = albedo_ta96(zzenith(:))
      zalbf_sca(:) = xalbsca_wat
      kk%XEMISF  (:) = xemiswat
      kk%XFFROZEN(:) = 0.0
    ELSEWHERE
      zalbf_dir(:) = xalbwatice
      zalbf_sca(:) = xalbwatice
      kk%XEMISF  (:) = xemiswatice
      kk%XFFROZEN(:) = 1.0
    END WHERE
    kk%XALBF(:)=0.5*(zalbf_dir(:)+zalbf_sca(:))
  ENDWHERE
  !
  zemisf(:) = kk%XEMISF(:)
  zff(:) = kk%XFF(:)
  !        
ELSE
  zalbf_dir(:)=0.0
  zalbf_sca(:)=0.0
  zemisf(:)=0.0
  zff(:)=0.0
ENDIF        
!-------------------------------------------------------------------------------
!
IF(io%LMEB_PATCH(kpatch))THEN
  !
  zsnowdepth(:) = sum(pek%TSNOW%WSNOW(:,:)/pek%TSNOW%RHO(:,:),2)
  zpalphan(:) = mebpalphan(zsnowdepth,pek%XH_VEG)
  !
  kk%XDIR_ALB_WITH_SNOW(:,:) = xundef
  kk%XSCA_ALB_WITH_SNOW(:,:) = xundef
  !
  IF(PRESENT(pdir_sw))THEN
    !
    ! Albedo
    !
    ! - just extract some parameters for call, but no need to update 
    !   the cummulative variables in this routine:
    !
    DO jswb=1,iswb
      zglobal_sw(:) = zdir_sw(:,jswb) + zsca_sw(:,jswb)
      !
      WHERE(pek%TSNOW%ALB(:)/=xundef .AND. pek%TSNOW%ALBVIS(:)/=xundef .AND. pek%TSNOW%ALBNIR(:)/=xundef)
        zlain(:)         = pek%XLAI(:)*(1.0-zpalphan(:))
        zalbvis_tsoil(:) = pek%XALBVIS_SOIL(:)*(1.-pek%XPSN(:)) + pek%XPSN(:)*pek%TSNOW%ALBVIS(:)
        zalbnir_tsoil(:) = pek%XALBNIR_SOIL(:)*(1.-pek%XPSN(:)) + pek%XPSN(:)*pek%TSNOW%ALBNIR(:)
      ELSEWHERE
        zlain(:)         = pek%XLAI(:)
        zalbvis_tsoil(:) = pek%XALBVIS_SOIL(:)
        zalbnir_tsoil(:) = pek%XALBNIR_SOIL(:)
      END WHERE
      !
      CALL radiative_transfert(io%LAGRI_TO_GRASS, kk%XVEGTYPE,                   &
              pek%XALBVIS_VEG, zalbvis_tsoil, pek%XALBNIR_VEG, zalbnir_tsoil,    &
              zglobal_sw, zlain, zzenith, s%XABC,                                &
              pek%XFAPARC, pek%XFAPIRC, pek%XMUS, pek%XLAI_EFFC, gshade, ziacan, &              
              ziacan_sunlit, ziacan_shade, zfrac_sun,                            &
              zfapar, zfapir, zfapar_bs, zfapir_bs                               )    

      ! Total effective surface (canopy, ground/flooded zone, snow) all-wavelength
      ! albedo: diagnosed from shortwave energy budget closure.
      ! Final note: purely diagnostic - apply limits for night time

      zalbt(:)      = 1. - (xsw_wght_vis*(zfapar(:)+zfapar_bs(:)) + xsw_wght_nir*(zfapir(:)+zfapir_bs(:)))
      zswup(:)      = zglobal_sw(:)*zalbt(:)
      zalbt(:)      = zswup(:)/max(1.e-5, zglobal_sw(:))
      !
      kk%XDIR_ALB_WITH_SNOW(:,jswb)=zalbt(:)
      kk%XSCA_ALB_WITH_SNOW(:,jswb)=zalbt(:) 
      !
    END DO
    !
  ENDIF
  !
  ! Emissivity
  !
  zemissn(:)   = xemissn
  zpsna(:)   = 0.
  zsigma_f(:)  = 1.0 - meb_shield_factor(pek%XLAI,zpsna)
  zsigma_fn(:) = 1.0 - meb_shield_factor(pek%XLAI,zpalphan)
  !
  CALL isba_emis_meb(pek%XPSN, zpalphan, zsigma_f, zsigma_fn, zemissn, zemist  )
  !
ELSE
  !        
  !  * albedo for near-infra-red and visible over snow-covered and snow-flood-free surface
  !
  zalbnir_with_snow(:) = pek%XALBNIR(:) * (1.-pek%XPSN(:)-zff(:)) + pek%TSNOW%ALB (:) * pek%XPSN(:)   
  zalbvis_with_snow(:) = pek%XALBVIS(:) * (1.-pek%XPSN(:)-zff(:)) + pek%TSNOW%ALB (:) * pek%XPSN(:)  
  zalbuv_with_snow(:) = pek%XALBUV (:) * (1.-pek%XPSN(:)-zff(:)) + pek%TSNOW%ALB (:) * pek%XPSN(:)  
  !
  !  * snow-flood-covered surface albedo for each wavelength (needed for outputs)
  !
  CALL albedo_from_nir_vis(psw_bands,                                              &
                           zalbnir_with_snow,  zalbvis_with_snow, zalbuv_with_snow,&
                           kk%XDIR_ALB_WITH_SNOW, kk%XSCA_ALB_WITH_SNOW          )  
  !
  DO jswb=1,iswb
    kk%XDIR_ALB_WITH_SNOW(:,jswb)=kk%XDIR_ALB_WITH_SNOW(:,jswb) + zff(:)*zalbf_dir(:)
    kk%XSCA_ALB_WITH_SNOW(:,jswb)=kk%XSCA_ALB_WITH_SNOW(:,jswb) + zff(:)*zalbf_sca(:)
  ENDDO
  !
  !-------------------------------------------------------------------------------
  !
  ! longwave computations for outputs (emissivity for radiative scheme)
  !
  zemist(:) = (1.-pek%XPSN(:)-zff(:))*pek%XEMIS(:) + pek%XPSN(:) * xemissn + zff(:)*zemisf(:)
  !
ENDIF
!
!Update albedo with snow for the next time step
!
 CALL unpack_same_rank(pk%NR_P,kk%XDIR_ALB_WITH_SNOW, pdir_alb_with_snow,zput0)
 CALL unpack_same_rank(pk%NR_P,kk%XSCA_ALB_WITH_SNOW, psca_alb_with_snow,zput0)
 CALL unpack_same_rank(pk%NR_P,zemist,pemist,zput0)
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('UPDATE_RAD_ISBA_N',1,zhook_handle)
!
END SUBROUTINE update_rad_isba_n