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(HISBA, HSNOW_ISBA, OTEMP_ARP,                     &
                          ptstep, psodelx,                                     &
                          psw_rad, plw_rad, pta, pqa,                          &
                          prhoa, pexns, pexna, pcps, plvtt, plstt,             &
                          pveg, phug, phui, phv,                               &
                          pleg_delta, plegi_delta, pdelta, pra,                &
                          pf5, prs, pcs, pcg, pct, psnowswe, ptsm, pt2m,       &
                          ppsn, ppsnv, ppsng, pfrozen1,                        &
                          palbt, pemist, pqsat, pdqsat, psnow_thrufal,         &
                          prn, ph, ple, pleg, plegi, plev,                     &
                          ples, pler, pletr, pevap,                            &
                          pgflux, pmeltadv, pmelt,                             &
                          psoilcondz, pd_g, pdzg, ptg,                         &
                          psr, ppsnv_a,                                        &
                          pffg, pffv, pff, pffrozen,                           &
                          ple_flood, plei_flood, psnowtemp                     )
!     ##########################################################################
!
!!****  *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_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
!
 CHARACTER(LEN=*),   INTENT(IN)      :: hisba      ! type of soil (Force-Restore OR Diffusion)
!                                                 ! '2-L'
!                                                 ! '3-L'
!                                                 ! 'DIF'   ISBA-DF
!
 CHARACTER(LEN=*), INTENT(IN)        :: hsnow_isba ! 'DEF' = Default F-R snow scheme
!                                                 !         (Douville et al. 1995)
!                                                 ! '3-L' = 3-L snow scheme (option)
!                                                 !         (Boone and Etchevers 2001)
!
LOGICAL, INTENT(IN)                 :: otemp_arp  ! True  = time-varying force-restore soil temperature (as in ARPEGE)
                                                  ! False = No time-varying force-restore soil temperature (Default)
!
!
REAL, INTENT (IN)                   :: ptstep     ! model time step (s)
!
!
REAL, DIMENSION(:), INTENT(IN)      :: psodelx    ! Pulsation for each layer (Only used if LTEMP_ARP=True)
!
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)      :: pveg
REAL, DIMENSION(:), INTENT(IN)      :: phug, phui, phv, pdelta, pra, prs, pf5
REAL, DIMENSION(:), INTENT(IN)      :: ppsn, ppsnv, ppsng, pfrozen1
REAL, DIMENSION(:), INTENT(IN)      :: palbt, pemist
REAL, DIMENSION(:), INTENT(IN)      :: pqsat, pdqsat
REAL, DIMENSION(:), INTENT(IN)      :: pleg_delta, plegi_delta
!                                      PVEG = fraction of vegetation
!                                      PHUG = relative humidity of the soil
!                                      PHV = Halstead coefficient
!                                      PF5 = water stress numerical correction factor (based on F2)
!                                      PDELTA = fraction of the foliage covered
!                                               by intercepted water
!                                      PRA = aerodynamic surface resistance for
!                                            heat transfers
!                                      PRS = surface stomatal resistance
!                                      PPSN = grid fraction covered by snow
!                                      PPSNV = fraction of the vegetation covered
!                                              by snow
!                                      PPSNG = fraction of the ground covered by
!                                              snow 
!                                      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, pcg, pct, pt2m, ptsm, psnowswe
!                                      PCT    = area-averaged heat capacity (K m2 J-1)
!                                      PCS    = heat capacity of the snow (K m2 J-1)
!                                      PCG    = heat capacity of the soil (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)
!                                      PSNOWSWE = equivalent water content of
!                                               the snow reservoir (kg m-2)
!
REAL, DIMENSION(:), INTENT(IN)      :: psnow_thrufal
!                                      PSNOW_THRUFAL = rate that liquid water leaves snow pack: 
!                                                     ISBA-ES [kg/(m2 s)]
REAL, DIMENSION(:), INTENT(IN)      :: psr 
!                                      PSR = snow precipitation rate [kg/(m2 s)]
REAL, DIMENSION(:), INTENT(IN)      :: ppsnv_a
!                                     PPSNV_A = vegetation covered by snow EBA scheme
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pd_g, psoilcondz
!                                      PD_G   = Depth of bottom of Soil layers (m)
!                                      PSOILCONDZ= ISBA-DF Soil conductivity profile  [W/(m K)]
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pdzg
!                                      PDZG   = Layer thickness (DIF option)
!
REAL, DIMENSION(:), INTENT(IN)      :: pcps, plvtt, plstt
!                                      PCPS  = heat capacity at surface
!
REAL, DIMENSION(:), INTENT(IN)      :: pffv      !Floodplain fraction over vegetation
REAL, DIMENSION(:), INTENT(IN)      :: pff       !Floodplain fraction at the surface
REAL, DIMENSION(:), INTENT(IN)      :: pffg   !Efective floodplain fraction
REAL, DIMENSION(:), INTENT(IN)      :: pffrozen  !fraction of frozen flood
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: ptg 
!                                      PTG    = soil temperature profile (K)
!
REAL, DIMENSION(:), INTENT(OUT)     :: ple_flood, plei_flood !Floodplains latent heat flux [W/m²]
REAL, DIMENSION(:), INTENT(OUT)     :: psnowtemp  ! snow layer temperatures (K)
!
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 PPSN, PPSNG, PPSNV
!                                                    if HSNOW_ISBA = '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 
!                                                    HSNOW_ISBA = '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 (hsnow_isba == '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(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
   zpsn(:)      = 0.0
   zpsng(:)     = 0.0
   zpsnv(:)     = 0.0
ELSE
   zpsn(:)      = ppsn(:)
   zpsng(:)     = ppsng(:)+pffg(:)
   zpsnv(:)     = ppsnv(:)+pffv(:)
   zfrac(:)     = ppsng(:)
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       1.     FLUX CALCULATIONS
!               -----------------
!
DO jj=1,SIZE(ptg,1)
!                                            temperature change
  zdt(jj)      = ptg(jj,1) - ptsm(jj)
!
!                                            net radiation
!
  prn(jj)      = (1. - palbt(jj)) * psw_rad(jj) + pemist(jj) *      &
              (plw_rad(jj) - xstefan * (ptsm(jj)** 3)*(4.*ptg(jj,1) - 3.*ptsm(jj)))
!
!                                            sensible heat flux
!
  ph(jj)       = prhoa(jj) * pcps(jj) * (ptg(jj,1) - pta(jj)*pexns(jj)/pexna(jj)) &
              / pra(jj) / pexns(jj)
!
  zwork1(jj) = prhoa(jj) * (1.-pveg(jj))*(1.-zpsng(jj)) / pra(jj)
  zwork2(jj) = pqsat(jj)+pdqsat(jj)*zdt(jj) 
!                                            latent heat of sublimation from
!                                            the ground
!
  plegi(jj)    = zwork1(jj) * plstt(jj) * ( phui(jj) * zwork2(jj) - pqa(jj)) * pfrozen1(jj) * plegi_delta(jj)
!
!                                            total latent heat of evaporation from
!                                            the ground
!
  pleg(jj)     = zwork1(jj) * plvtt(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) / pra(jj)
!                                            latent heat of evaporation from 
!                                            the snow canopy
!
  ples(jj)     = plstt(jj) * zpsn(jj) * zwork3(jj)
!
!                                            latent heat of evaporation from
!                                            evaporation
!
  plev(jj)     = plvtt(jj) * pveg(jj)*(1.-zpsnv(jj)) * phv(jj) * zwork3(jj)
!
!                                            latent heat of evapotranspiration
!                                            
  zzhv(jj)     = max(0., sign(1.,pqsat(jj) - pqa(jj)))
  pletr(jj)    = zzhv(jj) * (1. - pdelta(jj)) * plvtt(jj) * pveg(jj)*(1-zpsnv(jj))          &
               * zwork2(jj) *( (1/(pra(jj) + prs(jj))) - ((1.-pf5(jj))/(pra(jj) + xrs_max)) )
!               
!
  pler(jj)     = plev(jj) - pletr(jj)
!
!                                            latent heat of free water (floodplains)
!
  ple_flood(jj)  = plvtt(jj) * (1.-pffrozen(jj)) * pff(jj) * zwork3(jj) 
!
  plei_flood(jj) = plstt(jj) * pffrozen(jj) * pff(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 - ptg(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))/plvtt(jj)) + ((plegi(jj) + ples(jj))/plstt(jj))
!                                            total evaporative flux (kg/m2/s)
!                                            without flood
!
ENDDO
!
!-------------------------------------------------------------------------------
!
IF(hsnow_isba == 'D95')THEN
  DO jj=1,SIZE(ptg,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)/plvtt(jj) + plei_flood(jj)/plstt(jj)
  ENDDO
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       3.     SNOWMELT LATENT HEATING EFFECTS ('DEF' option)
!               ----------------------------------------------
!
IF( (hsnow_isba == 'D95' .OR. hsnow_isba == 'EBA') .AND. hisba /= 'DIF' )THEN
!                                            temperature tn
!
    IF (hsnow_isba == 'D95') THEN
!           
      ztn(:) = (1.-pveg(:))*ptg(:,1) + pveg(:)*pt2m(:)
!
!     Only diag
      psnowtemp(:) = ztn(:)
!
!
!                                            melting rate
!                                            there is melting only if T > T0 and
!                                            of course when SNOWSWE > 0.
!
      WHERE ( ztn(:) > xtt .AND. psnowswe(:) > 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(:) = psnowswe(:) + ptstep * (psr(:) - ples(:) / plstt(:))
!
      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 (hsnow_isba == 'EBA') THEN
!    
      pmelt(:)=min( psnowswe(:)/ptstep + psr(:) - ples(:)/plstt(:) , &
                    max(0.0,(ptg(:,1)-xtt))  / max(zeps1,pct*ptstep) / xlmtt )
!
    ENDIF
!
!                                            new temperature Ts(t) after melting
!                                            (cooling due to the melting of the
!                                            snow)
!
  ptg(:,1) = ptg(:,1) - pct(:)*xlmtt*pmelt(:)*ptstep
!
ENDIF
!
!
IF (lhook) CALL dr_hook('ISBA_FLUXES',1,zhook_handle)
!-------------------------------------------------------------------------------
END SUBROUTINE isba_fluxes