isba_fluxes.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.
!     ######spl
      SUBROUTINE isba_fluxes(IO, KK, PK, PEK, DMK, PTSTEP, &
                             PSW_RAD, PLW_RAD, PTA, PQA, PRHOA, PEXNS, PEXNA, &
                             PHUG, PHUI, PLEG_DELTA, PLEGI_DELTA, PDELTA, PF5, PCS, PTSM, PT2M, &
                             PFROZEN1, PALBT, PEMIST, PQSAT, PDQSAT, PSNOW_THRUFAL, &
                             PRN, PH, PLE, PLEG, PLEGI, PLEV,  PLES, PLER, PLETR, PEVAP, &
                             PGFLUX, PMELTADV, PMELT, PSOILCONDZ, PLE_FLOOD, PLEI_FLOOD)
!     ##########################################################################
!
!!****  *ISBA_FLUXES*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the simple snowpack schemes melt and the surface fluxes.
!         
!     
!!**  METHOD
!!    ------
!
!     1- snow melt latent heat, liquid rate (DEF option)
!     2- derive the surface fluxes.
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!    Douville et al. (1995)
!!    Boone et al. (2000)
!!      
!!    AUTHOR
!!    ------
!!
!!      S. Belair           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    14/03/95 
!!      (J.Stein)   15/11/95  use the wind components in the flux computation
!!      (J.Noilhan) 15/03/96  use the potential temperature instead of the
!!                            temperature for the heat flux computation 
!!      (J.Stein)   27/03/96  use only H and LE in the soil scheme
!!      (A.Boone V.Masson) 05/10/98 splits e_budget in two for CO2
!!      (A.Boone)   03/10/99  explicit latent heat of sublimation calculated 
!!      (A.Boone)   08/11/00  snow melt changes herein
!!      (A.Boone)   06/05/02  Updates, ordering. 
!!                            Introduction of snow melt timescale to 'DEF' snow option
!!      (P.LeMoigne) 01/07/05 expression of latent heat flux as a function of
!!                            w'theta' instead of w'T' (divison by surface exner)
!!      (P.LeMoigne) 28/07/05 dependence on qs for cp
!!      (A. Dziedzic and PLM) 10/2006 EBA snow option
!!      (B. Decharme)01/2009  Floodplains
!!      (R. Hamdi)   01/09    Cp and L are not constants (As in ALADIN)
!!      (B. Decharme)09/2009  Close the energy budget with the D95 snow scheme
!!      (A.Boone)    03/2010  Add delta fnctions to force LEG ans LEGI=0
!!                            when hug(i)Qsat < Qa and Qsat > Qa
!!      (A.Boone)    11/2011  Add RS_max limit to Etv
!!      (B. Decharme)07/2012  Error in restore flux calculation (only for diag)
!!      (B. Decharme)10/2012  Melt rate with D95 computed using max(XTAU,PTSTEP)
!!      (A.Boone)    02/2013  Split soil phase changes into seperate routine
!!      (B. Decharme)04/2013  Pass soil phase changes routines in hydro.F90
!!      (B. Decharme)04/2013  Delete PTS_RAD because wrong diagnostic
!!      (B. Decharme)10/14    "Restore" flux computed in e_budget
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
USE modd_csts,       ONLY : xstefan, xcpd, xlstt, xlvtt, xcl, xtt, xpi, xday, &
                            xci, xrholi, xlmtt, xrholw, xg, xcl, xcondi
USE modd_surf_par,   ONLY : xundef
USE modd_isba_par,   ONLY : xwgmin, xsphsoil, xdrywght, xrs_max
USE modd_snow_par,   ONLY : xtau_smelt
!
USE mode_thermos
!
USE mode_surf_snow_frac
!
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) :: KK
TYPE(isba_p_t), INTENT(INOUT) :: PK
TYPE(isba_pe_t), INTENT(INOUT) :: PEK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
REAL, INTENT (IN)                   :: PTSTEP     ! model time step (s)
!
REAL, DIMENSION(:), INTENT (IN)     :: PSW_RAD, PLW_RAD, PTA, PQA, PRHOA
!                                      PSW_RAD = incoming solar radiation
!                                      PLW_RAD = atmospheric infrared radiation
!                                      PTA = near-ground air temperature
!                                      PQA = near-ground air specific humidity
!                                      PRHOA = near-ground air density
!
REAL, DIMENSION(:), INTENT(IN)      :: PEXNS, PEXNA
REAL, DIMENSION(:), INTENT(IN)      :: PHUG, PHUI, PDELTA, PF5
REAL, DIMENSION(:), INTENT(IN)      :: PFROZEN1
REAL, DIMENSION(:), INTENT(IN)      :: PALBT, PEMIST
REAL, DIMENSION(:), INTENT(IN)      :: PQSAT, PDQSAT
REAL, DIMENSION(:), INTENT(IN)      :: PLEG_DELTA, PLEGI_DELTA
!                                      PHUG = relative humidity of the soil
!                                      PF5 = water stress numerical correction factor (based on F2)
!                                      PDELTA = fraction of the foliage covered
!                                               by intercepted water
!                                      PFROZEN1 = fraction of ice in near-surface
!                                                 ground
!                                      PALBT = area averaged albedo
!                                      PEMIST = area averaged emissivity
!                                      PQSAT = stauration vapor humidity at 't'
!                                      PDQSAT= stauration vapor humidity derivative at 't'
!                                      PLEG_DELTA = soil evaporation delta fn
!                                      PLEGI_DELTA = soil evaporation delta fn
!
REAL, DIMENSION(:), INTENT (IN)     :: PCS, PT2M, PTSM
!                                      PCS    = heat capacity of the snow (K m2 J-1)
!                                      PT2M   = mean surface (or restore) temperature at start 
!                                               of time step (K)
!                                      PTSM   = surface temperature at start 
!                                               of time step (K)
REAL, DIMENSION(:), INTENT(IN)      :: PSNOW_THRUFAL
!                                      PSNOW_THRUFAL = rate that liquid water leaves snow pack: 
!                                                     ISBA-ES [kg/(m2 s)]
REAL, DIMENSION(:,:), INTENT(IN)    :: PSOILCONDZ
!                                      PSOILCONDZ= ISBA-DF Soil conductivity profile  [W/(m K)]
!
REAL, DIMENSION(:), INTENT(OUT)     :: PLE_FLOOD, PLEI_FLOOD !Floodplains latent heat flux [W/m²]
!
REAL, DIMENSION(:), INTENT(OUT)     :: PRN, PH, PLE, PLEG, PLEV, PLES
REAL, DIMENSION(:), INTENT(OUT)     :: PLER, PLETR, PEVAP, PGFLUX, PMELTADV, PMELT
!                                     PRN = net radiation at the surface
!                                     PH = sensible heat flux
!                                     PLE = latent heat flux
!                                     PLEG = latent heat flux from the soil surface
!                                     PLEV = latent heat flux from the vegetation
!                                     PLES = latent heat flux from the snow
!                                     PLER = direct evaporation from the fraction
!                                            delta of the foliage
!                                     PLETR = transpiration of the remaining
!                                             part of the leaves
!                                     PEVAP = total evaporative flux (kg/m2/s)
!                                     PGFLUX = ground flux
!                                     PMELTADV = heat advection by melting snow
!                                                (acts to restore temperature to
!                                                 melting point) (W/m2)
!                                     PMELT = melting rate of snow (kg m-2 s-1)
!
REAL, DIMENSION(:), INTENT(OUT)     :: PLEGI
!                                      PLEGI   = sublimation component of the 
!                                                latent heat flux from the soil surface
!
!*      0.2    declarations of local variables
!
REAL                        :: ZKSOIL     ! coefficient for soil freeze/thaw
!
REAL, DIMENSION(SIZE(PTA))  :: ZZHV, ZTN, ZDT
!                                      ZZHV = for the calculation of the latent
!                                             heat of evapotranspiration
!!                                     ZTN  = average temperature used in 
!                                             the calculation of the 
!                                             melting effect
!                                      ZDT  = temperature change (K)
!
REAL, DIMENSION(SIZE(PTA))  ::  ZPSN, ZPSNV, ZPSNG, ZFRAC
!                               ZPSN, ZPSNV, ZPSNG = snow fractions corresponding to
!                                                    dummy arguments PEK%XPSN(:), PEK%XPSNG(:), PEK%XPSNV(:)
!                                                    if PEK%TSNOW%SCHEME = 'DEF' (composite
!                                                    or Force-Restore snow scheme), else
!                                                    they are zero for explicit snow case
!                                                    as snow fluxes calculated outside of
!                                                    this routine using the 
!                                                    PEK%TSNOW%SCHEME = '3-L' option.
!
REAL, DIMENSION(SIZE(PTA))  ::  ZNEXTSNOW
!                               ZNEXTSNOW = Future snow reservoir to close the
!                                           energy budget (see hydro_snow.f90)
!
REAL, DIMENSION(SIZE(PTA))  ::  ZCONDAVG
!                               ZCONDAVG   = average thermal conductivity of surface
!                                            and sub-surface layers (W m-1 K-1)
!
!
!*      0.2    local arrays for EBA scheme
!
REAL                            :: ZEPS1
!
!*      0.3    declarations of local parameters
!
INTEGER         :: JJ
!
REAL, DIMENSION(SIZE(PTA))      :: ZWORK1, ZWORK2, ZWORK3
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
!
!*       0.     Initialization
!               --------------
IF (lhook) CALL dr_hook('ISBA_FLUXES',0,zhook_handle)
!
IF (pek%TSNOW%SCHEME == 'EBA') zeps1=1.0e-8
!
pmelt(:)        = 0.0
pler(:)         = 0.0 
!
ztn(:)          = 0.0
zdt(:)          = 0.0
!
! If ISBA-ES option in use, then snow covered surface
! fluxes calculated outside of this routine, so set
! the local snow fractions here to zero:
! 
IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN
   zpsn(:)      = 0.0
   zpsng(:)     = 0.0
   zpsnv(:)     = 0.0
ELSE
   zpsn(:)      = pek%XPSN(:)
   zpsng(:)     = pek%XPSNG(:)+kk%XFFG(:)
   zpsnv(:)     = pek%XPSNV(:)+kk%XFFV(:)
   zfrac(:)     = pek%XPSNG(:)
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       1.     FLUX CALCULATIONS
!               -----------------
!
DO jj=1,SIZE(pek%XTG,1)
!                                            temperature change
  zdt(jj) = pek%XTG(jj,1) - ptsm(jj)
!
!                                            net radiation
!
  prn(jj) = (1. - palbt(jj)) * psw_rad(jj) + pemist(jj) *      &
           (plw_rad(jj) - xstefan * (ptsm(jj)** 3)*(4.*pek%XTG(jj,1) - 3.*ptsm(jj)))
!
!                                            sensible heat flux
!
  ph(jj) = prhoa(jj) * pk%XCPS(jj) * (pek%XTG(jj,1) - pta(jj)*pexns(jj)/pexna(jj)) &
           / pek%XRESA(jj) / pexns(jj)
!
  zwork1(jj) = prhoa(jj) * (1.-pek%XVEG(jj))*(1.-zpsng(jj)) / pek%XRESA(jj)
  zwork2(jj) = pqsat(jj)+pdqsat(jj)*zdt(jj) 
!                                            latent heat of sublimation from
!                                            the ground
!
  plegi(jj) = zwork1(jj) * pk%XLSTT(jj) * ( phui(jj) * zwork2(jj) - pqa(jj)) * pfrozen1(jj) * plegi_delta(jj)
!
!                                            total latent heat of evaporation from
!                                            the ground
!
  pleg(jj) = zwork1(jj) * pk%XLVTT(jj) * ( phug(jj) * zwork2(jj) - pqa(jj)) * (1.-pfrozen1(jj)) * pleg_delta(jj)
!
  zwork2(jj) = prhoa(jj) * (zwork2(jj) - pqa(jj))
  zwork3(jj) = zwork2(jj) / pek%XRESA(jj)
!                                            latent heat of evaporation from 
!                                            the snow canopy
!
  ples(jj)     = pk%XLSTT(jj) * zpsn(jj) * zwork3(jj)
!
!                                            latent heat of evaporation from
!                                            evaporation
!
  plev(jj)     = pk%XLVTT(jj) * pek%XVEG(jj)*(1.-zpsnv(jj)) * dmk%XHV(jj) * zwork3(jj)
!
!                                            latent heat of evapotranspiration
!                                            
  zzhv(jj) = max(0., sign(1.,pqsat(jj) - pqa(jj)))
  pletr(jj) = zzhv(jj) * (1. - pdelta(jj)) * pk%XLVTT(jj) * pek%XVEG(jj)*(1-zpsnv(jj))          &
              * zwork2(jj) *( (1/(pek%XRESA(jj) + dmk%XRS(jj))) - ((1.-pf5(jj))/(pek%XRESA(jj) + xrs_max)) )
!               
!
  pler(jj)     = plev(jj) - pletr(jj)
!
!                                            latent heat of free water (floodplains)
!
  ple_flood(jj)  = pk%XLVTT(jj) * (1.-kk%XFFROZEN(jj)) * kk%XFF(jj) * zwork3(jj) 
!
  plei_flood(jj) = pk%XLSTT(jj) * kk%XFFROZEN(jj) * kk%XFF(jj) * zwork3(jj) 
!
!                                            total latent heat of evaporation
!                                            without flood
!
  ple(jj)      = pleg(jj) + plev(jj) + ples(jj) + plegi(jj)
!
!                                            heat flux into the ground
!                                            without flood
!
  pgflux(jj)   = prn(jj) - ph(jj) - ple(jj)
!
!                                            heat flux due to snow melt
!                                            (ISBA-ES/SNOW3L)
!
  pmeltadv(jj) = psnow_thrufal(jj)*xcl*(xtt - pek%XTG(jj,1))
!
!                                            restore heat flux in FR mode,
!                                            or surface to sub-surface heat
!                                            flux using the DIF mode.
!
!
  pevap(jj)    = ((plev(jj) + pleg(jj))/pk%XLVTT(jj)) + ((plegi(jj) + ples(jj))/pk%XLSTT(jj))
!                                            total evaporative flux (kg/m2/s)
!                                            without flood
!
ENDDO
!
!-------------------------------------------------------------------------------
!
IF(pek%TSNOW%SCHEME == 'D95')THEN
  DO jj=1,SIZE(pek%XTG,1)
    ple(jj)  = ple(jj) + ple_flood(jj) + plei_flood(jj)
    pgflux(jj)  = pgflux(jj) - ple_flood(jj) - plei_flood(jj)
    pevap(jj)  = pevap(jj) + ple_flood(jj)/pk%XLVTT(jj) + plei_flood(jj)/pk%XLSTT(jj)
  ENDDO
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       3.     SNOWMELT LATENT HEATING EFFECTS ('DEF' option)
!               ----------------------------------------------
!
IF( (pek%TSNOW%SCHEME == 'D95' .OR. pek%TSNOW%SCHEME == 'EBA') .AND. io%CISBA /= 'DIF' )THEN
!                                            temperature tn
!
    IF (pek%TSNOW%SCHEME == 'D95') THEN
!           
      ztn(:) = (1.-pek%XVEG(:))*pek%XTG(:,1) + pek%XVEG(:)*pt2m(:)
!
!     Only diag
      dmk%XSNOWTEMP(:,1) = ztn(:)
!
!
!                                            melting rate
!                                            there is melting only if T > T0 and
!                                            of course when SNOWSWE > 0.
!
      WHERE ( ztn(:) > xtt .AND. pek%TSNOW%WSNOW(:,1) > 0.0 )
        pmelt(:) = zpsn(:)*(ztn(:)-xtt) / (pcs(:)*xlmtt*max(xtau_smelt,ptstep))
      END WHERE
!
!                                            close the energy budget: cannot melt 
!                                            more than the futur available snow
!      
      znextsnow(:) = pek%TSNOW%WSNOW(:,1) + ptstep * (dmk%XSRSFC(:) - ples(:) / pk%XLSTT(:))
!
      WHERE ( pmelt(:) > 0.0 )
!              
              pmelt(:)=min(pmelt(:),znextsnow(:)/ptstep)      
              znextsnow(:) = znextsnow(:) - ptstep * pmelt
!              
!             removes very small fraction
              zfrac(:) = snow_frac_ground(znextsnow(:))
              WHERE(zfrac(:)<1.0e-4)
                    pmelt(:)     = pmelt(:) + znextsnow(:) / ptstep       
              ENDWHERE   
!       
      ENDWHERE   
!    
    ELSEIF (pek%TSNOW%SCHEME == 'EBA') THEN
!    
      pmelt(:)=min( pek%TSNOW%WSNOW(:,1)/ptstep + dmk%XSRSFC(:) - ples(:)/ pk%XLSTT(:) , &
                  max(0.0,(pek%XTG(:,1)-xtt))  / max(zeps1,dmk%XCT*ptstep) / xlmtt )
!
    ENDIF
!
!                                            new temperature Ts(t) after melting
!                                            (cooling due to the melting of the
!                                            snow)
!
  pek%XTG(:,1) = pek%XTG(:,1) - dmk%XCT(:)*xlmtt*pmelt(:)*ptstep
!
ENDIF
!
!
IF (lhook) CALL dr_hook('ISBA_FLUXES',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE isba_fluxes