init_isba_sbl.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 init_isba_sbl(IO, K, NP, NPE, SB, PTSTEP, PPA, PPS, PTA, PQA, PRHOA, PU, PV, &
                             PDIR_SW, PSCA_SW, PSW_BANDS, PRAIN, PSNOW, PZREF, PUREF, PSSO_SLOPE )  
!     #################################################################################
!
!!****  *INIT_WATER_SBL* - inits water SBL profiles
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     S. Riette
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    03/2010
!!------------------------------------------------------------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_np_t, isba_npe_t, isba_p_t, isba_pe_t
USE modd_canopy_n, ONLY : canopy_t
!
USE modd_type_snow
!
USE modd_surf_par, ONLY : xundef
!
USE modd_csts,             ONLY : xcpd, xrd, xp00, xg, xlvtt
USE modd_surf_atm,         ONLY : lnosof
USE modd_canopy_turb,      ONLY : xalpsbl
!
USE modi_cls_tq
USE modi_isba_snow_frac
USE modi_wet_leaves_frac
USE modi_veg
USE modi_drag
USE modi_cls_wind
!
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_k_t), INTENT(INOUT) :: K
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(isba_npe_t), INTENT(INOUT) :: NPE
TYPE(canopy_t), INTENT(INOUT) :: SB
!
REAL,               INTENT(IN)   :: PTSTEP   ! timestep of the integration
REAL, DIMENSION(:), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
REAL, DIMENSION(:), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(:), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
REAL, DIMENSION(:), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(:), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
REAL, DIMENSION(:), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
REAL, DIMENSION(:), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
REAL, DIMENSION(:,:),INTENT(IN) :: PDIR_SW   ! direct  solar radiation (on horizontal surf.)
!                                            !                                       (W/m2)
REAL, DIMENSION(:,:),INTENT(IN) :: PSCA_SW   ! diffuse solar radiation (on horizontal surf.)
!                                            !                                        (W/m2)
REAL, DIMENSION(:), INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(:), INTENT(IN)  :: PSNOW     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(:), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
REAL, DIMENSION(:), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
REAL, DIMENSION(:), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
REAL, DIMENSION(:), INTENT(IN)  :: PSSO_SLOPE! slope of S.S.O.                         (-)
!                                            ! canopy       
!
!*      0.2    declarations of local variables
!
TYPE(isba_pe_t), POINTER :: PEK
TYPE(isba_p_t), POINTER :: PK
!
!* forcing variables
!
REAL, DIMENSION(SIZE(PTA))   :: ZWIND    ! lowest atmospheric level wind speed           (m/s)
REAL, DIMENSION(SIZE(PTA))   :: ZEXNA    ! Exner function at lowest SBL scheme level     (-)
REAL, DIMENSION(SIZE(PTA))   :: ZQA      ! specific humidity                             (kg/m3)
!
! SBL turbulence scheme
!
REAL, DIMENSION(SIZE(PTA))   ::ZRI
REAL, DIMENSION(SIZE(PTA))   ::ZCD
REAL, DIMENSION(SIZE(PTA))   ::ZCDN
REAL, DIMENSION(SIZE(PTA))   ::ZCH
REAL, DIMENSION(SIZE(PTA))   ::ZTNM
REAL, DIMENSION(SIZE(PTA))   ::ZQNM
REAL, DIMENSION(SIZE(PTA))   ::ZHUNM
REAL, DIMENSION(SIZE(PTA))   ::ZP_SLOPE_COS
REAL, DIMENSION(SIZE(PTA))   ::ZZ0
REAL, DIMENSION(SIZE(PTA))   ::ZZ0H
REAL, DIMENSION(SIZE(PTA))   ::ZEXNS
REAL, DIMENSION(SIZE(PTA))   ::ZTS
REAL, DIMENSION(SIZE(PTA))   ::ZHU
REAL, DIMENSION(SIZE(PTA))   ::ZQS
REAL, DIMENSION(SIZE(PTA))   ::ZZ0EFF
REAL, DIMENSION(SIZE(PTA))   ::ZWG
REAL, DIMENSION(SIZE(PTA))   ::ZWGI
REAL, DIMENSION(SIZE(PTA))   ::ZVEG
REAL, DIMENSION(SIZE(PTA))   ::ZRESA
REAL, DIMENSION(SIZE(PTA))   ::ZHUG
REAL, DIMENSION(SIZE(PTA))   ::ZHUGI
REAL, DIMENSION(SIZE(PTA))   ::ZHV
REAL, DIMENSION(SIZE(PTA))   ::ZCPS
REAL, DIMENSION(SIZE(PTA))   ::ZWRMAX_CF
REAL, DIMENSION(SIZE(PTA))   ::ZWR
REAL, DIMENSION(SIZE(PTA))   ::ZZ0_WITH_SNOW
REAL, DIMENSION(SIZE(PTA))   ::ZPSNG
REAL, DIMENSION(SIZE(PTA))   ::ZPSNV
REAL, DIMENSION(SIZE(PTA))   ::ZPSNV_A
REAL, DIMENSION(SIZE(PTA))   ::ZPSN
REAL, DIMENSION(SIZE(PTA))   ::ZSNOWALB
REAL, DIMENSION(SIZE(PTA))   ::ZFFG
REAL, DIMENSION(SIZE(PTA))   ::ZFFGNOS
REAL, DIMENSION(SIZE(PTA))   ::ZFFV
REAL, DIMENSION(SIZE(PTA))   ::ZFFVNOS
REAL, DIMENSION(SIZE(PTA))   ::ZFF
REAL, DIMENSION(SIZE(PTA))   ::ZRS
REAL, DIMENSION(SIZE(PTA))   ::ZP_GLOBAL_SW
REAL, DIMENSION(SIZE(PTA))   ::ZF2
REAL, DIMENSION(SIZE(PTA))   ::ZF5
REAL, DIMENSION(SIZE(PTA))   ::ZLAI
REAL, DIMENSION(SIZE(PTA))   ::ZGAMMA
REAL, DIMENSION(SIZE(PTA))   ::ZRGL
REAL, DIMENSION(SIZE(PTA))   ::ZRSMIN
REAL, DIMENSION(SIZE(PTA))   ::ZDELTA
REAL, DIMENSION(SIZE(PTA))   ::ZWRMAX
REAL, DIMENSION(SIZE(PTA))   ::ZCLS_WIND_ZON
REAL, DIMENSION(SIZE(PTA))   ::ZCLS_WIND_MER
REAL, DIMENSION(SIZE(PTA))   ::ZSUM
REAL, DIMENSION(SIZE(PTA))   :: ZLEG_DELTA  ! soil evaporation delta fn
REAL, DIMENSION(SIZE(PTA))   :: ZLEGI_DELTA ! soil sublimation delta fn
REAL, DIMENSION(SIZE(PTA))   :: ZLVTT
!
REAL, DIMENSION(:,:), ALLOCATABLE ::ZSNOWSWE
REAL, DIMENSION(:,:), ALLOCATABLE ::ZSNOWRHO
REAL, DIMENSION(:,:), ALLOCATABLE ::ZSUM_LAYER
!
INTEGER                     :: JSWB
INTEGER                     :: JL, JI, JP, IMASK
INTEGER :: ISNOW_LAYER
!
REAL, DIMENSION(SIZE(PTA),IO%NPATCH) ::ZWSNOW
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INIT_ISBA_SBL',0,zhook_handle)
!
zts(:) = 0.
zwg(:) = 0.
zwgi(:) = 0.
zz0(:) = 0.
!
zz0h(:) = 0.
zveg(:) = 0.
!
zresa(:) = 0.
!
zrgl(:) = 0.
zrsmin(:) = 0.
zgamma(:) = 0.
!
!Means over patches
DO jp = 1,io%NPATCH
  pk => np%AL(jp)
  pek => npe%AL(jp)

  DO ji=1,pk%NSIZE_P
    imask = pk%NR_P(ji)
    !
    zts(imask) = zts(imask) + pek%XTG (ji,1) * pk%XPATCH(ji)
    zwg(imask) = zwg(imask) + pek%XWG (ji,1) * pk%XPATCH(ji)
    zwgi(imask) = zwgi(imask) + pek%XWGI(ji,1) * pk%XPATCH(ji)
    zz0(imask) = zz0(imask) + pek%XZ0(ji)    * pk%XPATCH(ji)
    !
    zz0h(imask) = zz0h(imask) + pk%XPATCH(ji) * pek%XZ0 (ji) / pk%XZ0_O_Z0H(ji)
    zveg(imask) = zveg(imask) + pk%XPATCH(ji) * pek%XVEG(ji) 
    !    
    zresa(imask) = zresa(imask) + pk%XPATCH(ji) * pek%XRESA(ji)
    !
    zrgl(imask) = zrgl(imask) + pk%XPATCH(ji) * pek%XRGL  (ji)
    zrsmin(imask) = zrsmin(imask) + pk%XPATCH(ji) * pek%XRSMIN(ji)
    zgamma(imask) = zgamma(imask) + pk%XPATCH(ji) * pek%XGAMMA(ji)
    !  
  ENDDO
ENDDO
!
!We choose to set ZZ0EFF and ZZ0_WITH_SNOW equal to ZZ0
zz0eff(:)        = zz0(:)
zz0_with_snow(:) = zz0(:)
!
!
zlai(:) = 0.
zwrmax_cf(:) = 0.
zwr(:) = 0.
!
DO jp = 1,io%NPATCH
  pk => np%AL(jp)
  pek => npe%AL(jp)

  DO ji=1,pk%NSIZE_P
    imask = pk%NR_P(ji)
    !
    IF (zveg(ji)>0.) THEN
      zlai(imask) = zlai(imask) + pk%XPATCH(ji) * pek%XVEG(ji) *  pek%XLAI(ji)
      zwrmax_cf(imask) = zwrmax_cf(imask) + pk%XPATCH(ji) * pek%XVEG(ji) *  pek%XWRMAX_CF(ji)
      zwr(imask) = zwr(imask) + pk%XPATCH(ji) * pek%XVEG(ji) *  pek%XWR(ji)
    ELSEIF (jp==1) THEN
      zlai(imask) = pek%XLAI     (ji)
      zwrmax_cf(imask) = pek%XWRMAX_CF(ji)
      zwr(imask) = pek%XWR      (ji)
    ENDIF
    !
  ENDDO
ENDDO
!
WHERE (zveg(:)>0)
  zlai(:)= zlai(:) / zveg(:)
  zwrmax_cf(:)= zwrmax_cf(:) / zveg(:)
  zwr(:)= zwr(:) / zveg(:)
ENDWHERE
!
!
isnow_layer = npe%AL(1)%TSNOW%NLAYER
ALLOCATE(zsnowswe(SIZE(pta),isnow_layer))
ALLOCATE(zsum_layer(SIZE(pta),isnow_layer))
ALLOCATE(zsnowrho(SIZE(pta),isnow_layer))
zsnowswe(:,:) = 0.
zsum_layer(:,:) = 0.
zsnowrho(:,:) = 0.
!
DO jl=1,isnow_layer
  !
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)

    DO ji=1,pk%NSIZE_P
      imask = pk%NR_P(ji)
      !
      IF (pek%TSNOW%WSNOW(ji,jl)>0.) THEN
        zsnowswe(imask,jl) = zsnowswe(imask,jl) + pk%XPATCH(ji) * pek%TSNOW%WSNOW(ji,jl)
        zsum_layer(imask,jl) = zsum_layer(imask,jl) + pk%XPATCH(ji)
      ENDIF
      !
    ENDDO
  ENDDO
  !
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)

    DO ji=1,pk%NSIZE_P
      imask = pk%NR_P(ji)
      !
      IF (zsum_layer(imask,jl)>0.) THEN  
        zsnowrho(imask,jl) = zsnowrho(imask,jl) + pk%XPATCH(ji) * pek%TSNOW%RHO(ji,jl)
      ELSEIF (jp==1) THEN
        zsnowrho(imask,jl) = pek%TSNOW%RHO(ji,jl)
      ENDIF
      !
    ENDDO
  ENDDO
  !
END DO
!
WHERE (zsnowswe(:,:)==0.) zsnowrho(:,:) = xundef
WHERE (zsum_layer(:,:)>0.) 
  zsnowrho(:,:) = zsnowrho(:,:) / zsum_layer(:,:)
END WHERE
!
zsum(:)=sum(zsum_layer(:,:),dim=2)
DEALLOCATE(zsum_layer)
! 
zwsnow(:,:) = 0.
DO jl = 1,isnow_layer
  DO jp = 1,io%NPATCH
    pk => np%AL(jp)
    pek => npe%AL(jp)

    DO ji=1,pk%NSIZE_P
      imask = pk%NR_P(ji)
      !
      zwsnow(imask,jp) = zwsnow(imask,jp) + pek%TSNOW%WSNOW(ji,jl)
      !
    ENDDO
  ENDDO
ENDDO    
!
zsnowalb(:) = 0.
DO jp = 1,io%NPATCH
  pk => np%AL(jp)
  pek => npe%AL(jp)

  DO ji=1,pk%NSIZE_P
    imask = pk%NR_P(ji)
    !

    IF (zsum(imask)>0.) THEN
      IF (zwsnow(imask,jp)>0.) THEN
        zsnowalb(imask) = zsnowalb(imask) + pk%XPATCH(ji) * pek%TSNOW%ALB(ji)
      ENDIF
    ELSEIF (jp==1) THEN
      zsnowalb(imask) = pek%TSNOW%ALB(ji)
    ENDIF
    !
  ENDDO
ENDDO
!
WHERE(zsum(:)>0)         
  zsnowalb(:) = zsnowalb(:) / zsum(:)      
ENDWHERE
!
!
zexna(:) = (ppa(:)/xp00)**(xrd/xcpd)
zexns(:) = (pps(:)/xp00)**(xrd/xcpd)
zqa(:) = pqa(:) / prhoa(:)
zwind(:) = sqrt(pu**2+pv**2)
!
!We compute the snow fractions
 CALL isba_snow_frac(pek%TSNOW%SCHEME, zsnowswe, zsnowrho, zsnowalb,   &
                     zveg, zlai, zz0, zpsn, zpsnv_a, zpsng, zpsnv   )  
!
DEALLOCATE(zsnowswe, zsnowrho)
!
!We compute total shortwave incoming radiation needed by veg
zp_global_sw(:) = 0.
DO jswb=1,SIZE(psw_bands)
  zp_global_sw(:)   = zp_global_sw(:) + (pdir_sw(:,jswb) + psca_sw(:,jswb))
END DO
!
!
!We choose the case HPHOTO=='NON' and a humid soil (ZF2=1) to compute ZRS
zf2(:)=1.0
 CALL veg(zp_global_sw, pta, zqa, pps, zrgl, zlai, zrsmin, zgamma, zf2, zrs)
!Calculation of ZDELTA
 CALL wet_leaves_frac(zwr, zveg, zwrmax_cf, zz0_with_snow, zlai, zwrmax, zdelta)
!
!We choose the case LFLOOD=false
zffg(:) = 0.0
zffgnos(:) = 0.0
zffv(:) = 0.0
zffvnos(:) = 0.0
zff(:) = 0.0
!
zf5(:) = 1.0
zlvtt(:) = xlvtt
!
zp_slope_cos(:) = 1./sqrt(1.+psso_slope(:)**2)
IF (lnosof) zp_slope_cos(:) = 1.0
!
!We compute ZCD, ZCH and ZRI
 CALL drag(io%CISBA, pek%TSNOW%SCHEME, io%CCPSURF,  ptstep, zts, zwg, zwgi, &
           zexns, zexna, pta, zwind, zqa, prain, psnow, pps, zrs, zveg,    &
           zz0, zz0eff, zz0h, k%XWFC(:,1), k%XWSAT(:,1), zpsng, zpsnv,   &
           pzref, puref, zp_slope_cos, zdelta, zf5, zresa, zch, zcd, zcdn, &
           zri, zhug, zhugi, zhv, zhu, zcps, zqs, zffg, zffv, zff, zffgnos,&
           zffvnos, zleg_delta, zlegi_delta, zwr, prhoa, zlvtt            )  
!
!Initialisation of T, Q, Wind and TKE on all canopy levels
DO jl=1,sb%NLVL
  !
  CALL cls_tq(pta, zqa, ppa, pps, pzref, zcd, zch, zri, zts, zhu, zz0h, &
              sb%XZ(:,jl), ztnm, zqnm, zhunm           ) 
  ! 
  sb%XT(:,jl)=ztnm
  sb%XQ(:,jl)=zqnm
  !
  CALL cls_wind(pu, pv, puref, zcd, zcdn, zri, sb%XZ(:,jl), &
                zcls_wind_zon, zcls_wind_mer                 )
  !
  sb%XU   (:,jl) = sqrt( zcls_wind_zon(:)**2 + zcls_wind_mer(:)**2 )
  sb%XTKE (:,jl) = xalpsbl * zcd(:) * ( pu(:)**2 + pv(:)**2 )
  sb%XP   (:,jl) = ppa(:) + xg * prhoa(:) * (sb%XZ(:,sb%NLVL) - sb%XZ(:,jl))
  !
ENDDO
!
IF (lhook) CALL dr_hook('INIT_ISBA_SBL',1,zhook_handle) 
!
END SUBROUTINE init_isba_sbl