z0eff.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 z0eff (I, &
                      hrough, omeb, palfa, pzref, puref, pz0, pz0rel, ppsn,      &
                      ppalphan,pz0litter, pwsnow,                               &
                      pz0effip,pz0effim,pz0effjp,pz0effjm,pff,pz0_flood,        &
                      paosip,paosim,paosjp,paosjm,pho2ip,pho2im,pho2jp,pho2jm,  &
                      pz0_o_z0h, pz0_with_snow, pz0h_with_snow,pz0eff,          &
                      pz0g_without_snow,                                        &
                      pz0_mebv,pz0h_mebv,pz0eff_mebv,                           &
                      pz0_mebn,pz0h_mebn,pz0eff_mebn                            )

!   ############################################################################
!
!!****  *Z0EFF*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the z0eff for momentum fluxes according to wind direction.
!         
!     
!!**  METHOD
!!    ------
!
!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Mascart et al. (1995)
!!    Belair (1995)
!!      
!!    AUTHOR
!!    ------
!!
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    13/03/95 
!!      (J.Stein)   15/11/95  use the potential temperature to compute Ri
!!                            and PVMOD instead of ZVMOD
!!      (P.Lacarrere)15/03/96 replace * PEXNS by / PEXNS
!!      (V.Masson)   22/12/97 computation of z0eff after snow treatment
!!      (V.Masson)   05/10/98 clear routine
!!      (A.Boone)    11/26/98 Option for PDELTA: forested vs default surface
!!      (V Masson)   12/07/01 new formulation for aggregation with snow z0
!!      (P.LeMoigne) 09/02/06 computation of z0h in presence of snow
!!      (B. Decharme)    2008 floodplains
!!      (P. Samuelsson) 10/2014 MEB
!!      P. LeMoigne  12/2014 EBA scheme update
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_isba_n, ONLY : isba_t
!
USE modd_csts,     ONLY : xpi, xg
USE modd_snow_par, ONLY : xz0sn, xwcrn, xz0hsn
!
USE modi_subscale_z0eff
USE modd_surf_atm, ONLY : laldz0h
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
TYPE(isba_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=*),   INTENT(IN)  :: hrough         ! type of roughness length
LOGICAL, INTENT(IN)             :: omeb           ! True = patch with multi-energy balance 
!                                                 ! False = patch with classical ISBA
REAL, DIMENSION(:), INTENT(IN)  :: palfa          ! wind direction from J axis (clockwise)
REAL, DIMENSION(:), INTENT(IN)  :: pzref          ! height of atmospheric level
REAL, DIMENSION(:), INTENT(IN)  :: puref          ! reference height for wind
REAL, DIMENSION(:), INTENT(IN)  :: pz0            ! vegetation roughness length
REAL, DIMENSION(:), INTENT(IN)  :: pz0rel         ! 1d orographic roughness length
REAL, DIMENSION(:), INTENT(IN)  :: ppsn           ! fraction of snow
REAL, DIMENSION(:), INTENT(IN)  :: ppalphan       ! snow/canopy transition coefficient
REAL, DIMENSION(:), INTENT(IN)  :: pz0effip       ! z0eff for increasing x
REAL, DIMENSION(:), INTENT(IN)  :: pz0effim       ! z0eff for decreasing x
REAL, DIMENSION(:), INTENT(IN)  :: pz0effjp       ! z0eff for increasing y
REAL, DIMENSION(:), INTENT(IN)  :: pz0effjm       ! z0eff for decreasing y
REAL, DIMENSION(:), INTENT(IN)  :: paosip         ! A/S for increasing x
REAL, DIMENSION(:), INTENT(IN)  :: paosim         ! A/S for decreasing x
REAL, DIMENSION(:), INTENT(IN)  :: paosjp         ! A/S for increasing y
REAL, DIMENSION(:), INTENT(IN)  :: paosjm         ! A/S for decreasing y
REAL, DIMENSION(:), INTENT(IN)  :: pho2ip         ! h/2 for increasing x
REAL, DIMENSION(:), INTENT(IN)  :: pho2im         ! h/2 for decreasing x
REAL, DIMENSION(:), INTENT(IN)  :: pho2jp         ! h/2 for increasing y
REAL, DIMENSION(:), INTENT(IN)  :: pho2jm         ! h/2 for decreasing y
REAL, DIMENSION(:), INTENT(IN)  :: pz0_o_z0h      ! ratio between heat and momentum z0
!
REAL, DIMENSION(:), INTENT(IN)  :: pff            ! fraction of flood
REAL, DIMENSION(:), INTENT(IN)  :: pz0_flood      ! floodplains roughness length
!
! For multi-energy balance
REAL, DIMENSION(:), INTENT(IN)  :: pz0litter      ! ground litter roughness length for MEB
!
REAL, DIMENSION(:), INTENT(IN)  :: pwsnow         ! equivalent snow water content
!
REAL, DIMENSION(:), INTENT(OUT) :: pz0_with_snow  ! vegetation z0 modified by snow
REAL, DIMENSION(:), INTENT(OUT) :: pz0h_with_snow ! vegetation z0h modified by snow
REAL, DIMENSION(:), INTENT(OUT) :: pz0eff         ! effective z0
!
! For multi-energy balance
REAL, DIMENSION(:), INTENT(OUT) :: pz0g_without_snow  ! roughness length for momentum at snow-free canopy floor
!
REAL, DIMENSION(:), INTENT(OUT) :: pz0_mebv           ! roughness length for momentum over MEB vegetation part of patch
REAL, DIMENSION(:), INTENT(OUT) :: pz0h_mebv          ! roughness length for heat over MEB vegetation part of path
REAL, DIMENSION(:), INTENT(OUT) :: pz0eff_mebv        ! roughness length for momentum over MEB vegetation part of patch
!                                                     ! eventually including orograhic roughness
REAL, DIMENSION(:), INTENT(OUT) :: pz0_mebn           ! roughness length for momentum over MEB snow part of patch
REAL, DIMENSION(:), INTENT(OUT) :: pz0h_mebn          ! roughness length for heat over MEB snow part of path
REAL, DIMENSION(:), INTENT(OUT) :: pz0eff_mebn        ! roughness length for momentum over MEB snow part of patch
!
!
!
!*      0.2    declarations of local variables
!
!
!
REAL, DIMENSION(SIZE(PZ0EFF)) :: zwork, zalfa,       &
                                   zz0effip, zz0effim, &
                                   zz0effjp, zz0effjm, &
                                   zpff
!                                              effective roughness length in 4
!                                              directions
REAL                          :: z0cr, zuz0cn, zalrcn1, zalrcn2
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('Z0EFF',0,zhook_handle)
zalrcn1=1.e-02
zalrcn2=2.5e-03
z0cr = zalrcn1
zuz0cn=1./zalrcn2
zalfa(:) = palfa(:)
WHERE(zalfa(:)<=-xpi) zalfa = zalfa + 2.*xpi
WHERE(zalfa(:)>  xpi) zalfa = zalfa - 2.*xpi
!
! Initialisation of MEB roughness lengths
pz0g_without_snow=0.
pz0_mebv=0.
pz0h_mebv=0.
pz0eff_mebv=0.
pz0_mebn=0.
pz0h_mebn=0.
pz0eff_mebn=0.
!
!*       1.     GRID-AVERAGED ROUGHNESS LENGTHS
!               -------------------------------
!       (considering the effect of snow-flood-covered surfaces and orography)
!
!*       1.1    for heat
!               --------
!
pz0_with_snow(:)  = pz0(:)
pz0h_with_snow(:) = pz0(:) / pz0_o_z0h(:)
!
IF(i%TSNOW%SCHEME=='EBA') THEN
!        
   WHERE (ppsn(:)>0.)
!
!!!!!Flooding scheme not implemented with this option 
      pz0_with_snow(:) = pz0_with_snow(:) + ( z0cr - pz0(:))* &
        pwsnow(:)/(pwsnow(:) + xwcrn*(1.0+zuz0cn*pz0(:)))  
!        
   END WHERE 


   IF (laldz0h) THEN  
     WHERE (ppsn(:)>0.)
         pz0h_with_snow(:) = pz0h_with_snow(:) + ( z0cr - pz0h_with_snow(:))* &
           pwsnow(:)/(pwsnow(:) + xwcrn*(1.0+zuz0cn*pz0h_with_snow(:)))   
     END WHERE  
  END IF   
    
!        
ELSE
!        
   WHERE (ppsn(:)>0..OR.pff(:)>0.)
!        
      zwork(:) =  (            ppsn(:) /(log(puref(:)/xz0sn       ))**2 ) &
                + (            pff(:) /(log(puref(:)/pz0_flood(:)))**2 ) &
                + ((1.-ppsn(:)-pff(:))/(log(puref(:)/pz0(:)      ))**2 )  
!
      pz0_with_snow(:) = puref(:) /exp( sqrt( 1./zwork(:) ) )
!
      zwork(:) =  (            ppsn(:) /(log(pzref(:)/xz0hsn                      ))**2 ) &
                + (            pff(:) /(log(pzref(:)/(pz0_flood(:)/ pz0_o_z0h(:))))**2 ) &
                + ((1.-ppsn(:)-pff(:))/(log(pzref(:)/(pz0(:)/pz0_o_z0h(:))       ))**2 )  
!
      pz0h_with_snow(:) = pzref(:) /exp( sqrt( 1./zwork(:) ) )
!
   END WHERE
!
ENDIF
!
! For multi-energy balance
IF(omeb)THEN

! roughness length for momentum at snow-free canopy floor
  pz0g_without_snow(:) = pz0litter
  WHERE (pff(:)>0.)
    zpff(:)=pff(:)/(1-ppsn(:)+1.e-6)
    zwork(:) =  (    zpff(:)  * log(pz0_flood(:)) )   &
              + ( (1.-zpff(:)) * log(pz0litter(:)) )
    pz0g_without_snow(:) = exp( zwork(:) )
  END WHERE
!
! roughness length for momentum over MEB vegetation part of patch
  pz0_mebv(:)  = pz0(:)
!
! roughness length for momentum over MEB snow part of patch
  zwork(:) =  (    ppalphan(:) /(log(puref(:)/xz0sn       ))**2 ) &
            + ((1.-ppalphan(:))/(log(puref(:)/pz0_mebv(:)      ))**2 )
  pz0_mebn(:) = puref(:) /exp( sqrt( 1./zwork(:) ) )
!

! roughness length for momentum over MEB total patch
  zwork(:) =  (    ppsn(:) /(log(puref(:)/pz0_mebn(:)  ))**2 ) &
            + ((1.-ppsn(:))/(log(puref(:)/pz0_mebv(:)  ))**2 )
  pz0_with_snow(:) = puref(:) /exp( sqrt( 1./zwork(:) ) )
!
! roughness length for heat over MEB vegetation part of path
  pz0h_mebv(:) = pz0_mebv(:)/pz0_o_z0h(:)
! for nordic forest, z0h=z0m according to M&#195;&#182;lder (tested in Hirlam):
!
! PZ0H_MEBV(:) = PZ0_MEBV(:)   
!
! roughness length for heat over MEB snow part of path
  zwork(:) =  (     ppalphan(:) /(log(pzref(:)/xz0hsn          ))**2 ) &
            + ( (1.-ppalphan(:))/(log(pzref(:)/pz0h_mebv(:)    ))**2 )
  pz0h_mebn(:) = pzref(:) /exp( sqrt( 1./zwork(:) ) )
!

! roughness length for heat over MEB total path
  zwork(:) =  (     ppsn(:) /(log(pzref(:)/pz0h_mebn(:)    ))**2 ) &
            + ( (1.-ppsn(:))/(log(pzref(:)/pz0h_mebv(:)    ))**2 )
  pz0h_with_snow(:) = pzref(:) /exp( sqrt( 1./zwork(:) ) )
!
ENDIF
!
!*       1.2    for momentum
!               ------------
!
!
!                                     In this particular case, we now use
!                                     the roughness length due to the coupled
!                                     effect of vegetation and topography
!                                     Snow and Flood effects are yet taken
!                                     into account through ZZ0EFF
!
IF (hrough=='Z04D') THEN
!
! For multi-energy balance (MEB): the HROUGH=='Z04D' option is not considered yet!
  !
  zz0effip(:) = pz0effip(:)
  zz0effim(:) = pz0effim(:)
  zz0effjp(:) = pz0effjp(:)
  zz0effjm(:) = pz0effjm(:)
  !
  CALL subscale_z0eff(paosip,paosim,paosjp,paosjm,               &
                        pho2ip,pho2im,pho2jp,pho2jm,pz0_with_snow, &
                        zz0effip,zz0effim,zz0effjp,zz0effjm,       &
                        omask=(ppsn>0..OR.pff(:)>0.)               )  
  !
  WHERE(zalfa(:)>=0. .AND. zalfa(:)<xpi/2.)
    pz0eff(:)=zz0effip(:)*sin(zalfa(:))**2 + zz0effjp(:)*cos(zalfa(:))**2
  END WHERE
  WHERE(zalfa(:)>=xpi/2. .AND. zalfa(:)<=xpi)
    pz0eff(:)=zz0effip(:)*sin(zalfa(:))**2 + zz0effjm(:)*cos(zalfa(:))**2
  END WHERE
  WHERE (zalfa(:)>=-xpi/2 .AND. zalfa(:)<0.)
    pz0eff(:)=zz0effim(:)*sin(zalfa(:))**2 + zz0effjp(:)*cos(zalfa(:))**2
  END WHERE
  WHERE (zalfa(:)>=-xpi .AND. zalfa(:)<-xpi/2.)
    pz0eff(:)=zz0effim(:)*sin(zalfa(:))**2 + zz0effjm(:)*cos(zalfa(:))**2
  END WHERE
!
ELSE IF (hrough=='Z01D') THEN
  pz0eff(:) = pz0_with_snow(:) + pz0rel(:)
  IF(omeb)THEN
    pz0eff_mebv(:) = pz0_mebv(:) + pz0rel(:)
    pz0eff_mebn(:) = pz0_mebn(:) + pz0rel(:)
  ENDIF
  IF (laldz0h) THEN
     ! Aladin dynamic z0 contains already orographic component
     pz0eff(:) = pz0eff(:) - pz0rel(:)
     IF(omeb)THEN
       pz0eff_mebv(:) = pz0eff_mebv(:) - pz0rel(:)
       pz0eff_mebn(:) = pz0eff_mebn(:) - pz0rel(:)
     ENDIF
!     PZ0H_WITH_SNOW(:) = PZ0EFF(:) / PZ0_O_Z0H(:)   ! it is aleardy corrected under IF statement of TSNOW%SCHEME  
  ENDIF
ELSE
  pz0eff(:) = pz0_with_snow(:)
  IF(omeb)THEN
    pz0eff_mebv(:) = pz0_mebv(:)
    pz0eff_mebn(:) = pz0_mebn(:)
  ENDIF
END IF
IF (lhook) CALL dr_hook('Z0EFF',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE z0eff