isba_ceb.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 isba_ceb(HISBA, HSNOW_ISBA, HCPSURF, OFLOOD, OTEMP_ARP, HIMPLICIT_WIND, &
                    ptstep, psodelx, ppew_a_coef, ppew_b_coef, ppet_a_coef,        &
                    ppeq_a_coef, ppet_b_coef, ppeq_b_coef, psnowalb,               &
                    psw_rad, plw_rad, pwg, pwgi, pexns, pexna, pta, pvmod,         &
                    pqa, prr, psr, pps, prs, pveg, pz0_with_snow, pz0eff,          &
                    pz0h_with_snow, pwfc, pwsat, ppsn, ppsng, ppsnv, pzref,        &
                    puref, pdircoszw, pf5, pffg, pffv, pff, pffg_nosnow,           &
                    pffv_nosnow, pwr, prhoa, pemis, palb, pct, pcs, pcg,           &
                    pd_g, pdzg, pdzdif, psoilcondz, psoilhcapz,  pfrozen1,         &
                    ptdeep_a, ptdeep_b, pgammat,  ppsnv_a, psnowfree_alb_veg,      &
                    psnowfree_alb_soil, pgrndflux, pflsn_cor,                      &
                    psnow_thrufal, pffrozen, pfalb, pfemis, psnowswe, psrsfc,      &
                    ptg, presa, plvtt, plstt, pcps, pdelta, pch, pcd, pcdn,        &
                    pri, phug, phugi, phv, phu, pqs, palbt, pemist, pdeep_flux,    &
                    prn, ph, ple, pleg, plegi, plev, ples, pler, pletr, pevap,     &
                    pgflux, pmeltadv, pmelt, prestore, pustar, ple_flood,          &
                    plei_flood, psnowtemp, pac_agg, phu_agg                        )
!     ##########################################################################
!
!
!!****  *ISBA_CEB*  
!!
!!    PURPOSE
!!    -------
!       Monitor for the calculation of the surface composit energy budget
!!      call of drag, e_budget, and isba_fluxes
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!      
!!    AUTHOR
!!    ------
!!      B. Decharme           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    10/03/95 
!!
!!      (B. Decharme)   03/16 Bug : limitation of Er for Interception reservoir
!!                                  PTSTEP insted of ZTSTEP in drag.F90
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
USE modd_surf_par,   ONLY : xundef
!
USE modd_surf_atm,   ONLY : lcpl_arp
!
USE modi_drag
USE modi_e_budget
USE modi_isba_fluxes
!
USE mode_thermos
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!              -------------------------
!
! in
!
!
 CHARACTER(LEN=*),     INTENT(IN) :: hisba      ! type of ISBA version:
!                                              ! '2-L' (default)
!                                              ! '3-L'
!                                              ! 'DIF'
 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 2000)
 CHARACTER(LEN=*),     INTENT(IN) :: hcpsurf    ! option for specific heat Cp:
!                                              ! 'DRY' = dry Cp
!                                              ! 'HUM' = Cp as a function of qs
!
LOGICAL, INTENT(IN)              :: oflood     ! Activation of the flooding scheme
LOGICAL, INTENT(IN)              :: otemp_arp  ! True  = time-varying force-restore soil temperature (as in ARPEGE)
                                               ! False = No time-varying force-restore soil temperature (Default)
!
 CHARACTER(LEN=*),    INTENT(IN) :: himplicit_wind   ! wind implicitation option
!                                                    ! 'OLD' = direct
!                                                    ! 'NEW' = Taylor serie, order 1
!
REAL,                 INTENT(IN) :: ptstep    ! timestep of the integration
!
REAL, DIMENSION(:),   INTENT(IN) ::  psodelx  ! Pulsation for each layer (Only used if LTEMP_ARP=True)
!
REAL, DIMENSION(:),  INTENT(IN)  :: ppew_a_coef, ppew_b_coef,                   &
                                    ppet_a_coef, ppeq_a_coef, ppet_b_coef,      &
                                    ppeq_b_coef  
!                                  PPEW_A_COEF = A-wind coefficient (m2s/kg)
!                                  PPEW_B_COEF = B-wind coefficient (m/s)
!                                  PPET_A_COEF = A-air temperature coefficient
!                                  PPET_B_COEF = B-air temperature coefficient
!                                  PPEQ_A_COEF = A-air specific humidity coefficient
!                                  PPEQ_B_COEF = B-air specific humidity coefficient
!
REAL, DIMENSION(:), INTENT(IN)   :: psnowalb
!                                     prognostic variables at time 't-dt'
!                                     PSNOWALB = albedo of the snow
!
REAL, DIMENSION(:,:), INTENT(IN) :: pwg, pwgi, pwfc, pwsat
!                                     PWG     = near-surface volumetric water content
!                                     PWGI    = near-surface frozen volumetric water content
!                                     PWFC = field capacity volumetric water content
!                                     PWSAT = volumetric water content at saturation
!
REAL, DIMENSION(:), INTENT(IN)   :: psw_rad, plw_rad
!                                     PSW_RAD = incoming solar radiation
!                                     PLW_RAD = atmospheric infrared radiation

REAL, DIMENSION(:), INTENT(IN)   :: pexna, pexns, pta, pvmod, pqa, prr, psr, pps
!                                     PEXNA= Exner function near surface atmospheric variables
!                                     PEXNS   = Exner function at the surface
!                                     PTA = 2m temperature
!                                     PVMOD = module of the horizontal wind
!                                             NOTE it should be input as >= 1. (m)
!                                     PQA = specific humidity
!                                     PPS = surface pressure
!                                     PRR = rain rate    
!                                     PSR = snow rate             
!
REAL, DIMENSION(:), INTENT(IN)   :: prs, pveg
!                                     PRS = stomatal resistance
!                                     PVEG = vegetation fraction
REAL, DIMENSION(:), INTENT(IN)   :: ppsn, ppsng, ppsnv
!                                     PPSN = grid fraction covered by snow
!                                     PPSNG = fraction of the bare ground covered
!                                             by snow
!                                     PPSNV = fraction of the vegetation covered
!                                             by snow
REAL, DIMENSION(:), INTENT(IN)   :: pzref, puref
!                                     PZREF = reference height of the first
!                                             atmospheric level 
!                                     PUREF = reference height of the wind
!                                             NOTE this is different from ZZREF
!                                             ONLY in stand-alone/forced mode,
!                                             NOT when coupled to a model (MesoNH)
!
REAL, DIMENSION(:), INTENT(IN)    :: pz0eff         ! roughness length for momentum
REAL, DIMENSION(:), INTENT(IN)    :: pz0_with_snow  ! roughness length for momentum
!                                                   ! (with snow taken into account)
REAL, DIMENSION(:), INTENT(IN)    :: pz0h_with_snow ! roughness length for heat
!                                                   ! (with snow taken into account)
!
REAL, DIMENSION(:), INTENT(IN)    :: pf5
!                                     PF5 = water stress function for Hv
REAL, DIMENSION(:), INTENT(IN)    :: pdircoszw 
!                                     PDIRCOSZW = Cosinus of the angle between the normal to the surface and the vertical
!
REAL, DIMENSION(:), INTENT(IN)   :: pffv, pff, pffg, pffg_nosnow, pffv_nosnow
!                                   PFFG = Floodplain fraction over ground
!                                   PFFV = Floodplain fraction over vegetation
!                                   PFF  = Floodplain fraction at the surface
!
REAL, DIMENSION(:), INTENT(IN)   :: pwr, prhoa
!                                   PWR = liquid water retained on the foliage
!                                   PRHOA = near-ground air density
!
REAL, DIMENSION(:), INTENT(IN)   :: pemis, palb, pct, pcs, pcg
!                                     PEMIS = emissivity
!                                     PALB = albedo
!                                     PCT = area-averaged heat capacity
!                                     PCS    = heat capacity of the snow (K m2 J-1)
!                                     PCG    = heat capacity of the soil (K m2 J-1)
REAL, DIMENSION(:), INTENT(IN)   :: pgrndflux, pflsn_cor, psnow_thrufal 
!                                     PGRNDFLUX = soil/snow interface flux (W/m2) using
!                                                 ISBA-SNOW3L option
!                                     PFLSN_COR = soil/snow interface correction flux to conserve energy (W/m2)
!                                     PSNOW_THRUFAL  = snow runoff/melt leaving pack and available
!                                                  at the surface for runoff or infiltration
!                                                  [kg/(m2 s)]
REAL, DIMENSION(:,:), INTENT(IN)  :: pd_g,  psoilcondz, psoilhcapz
!                                     PD_G      = Depth of bottom of Soil layers (m)
!                                     PSOILCONDZ= ISBA-DF Soil conductivity profile  [W/(m K)]
!                                     PSOILHCAPZ=ISBA-DF Soil heat capacity profile [J/(m3 K)]
!
REAL, DIMENSION(:,:), INTENT(IN)  :: pdzg       ! soil layers thicknesses (DIF option) (m)
REAL, DIMENSION(:,:), INTENT(IN)  :: pdzdif     ! distance between consecuative layer mid-points (DIF option) (m)
!
REAL, DIMENSION(:), INTENT(IN)    :: pfrozen1
!                                     PFROZEN1 = ice fraction in supurficial soil
!
REAL, DIMENSION(:), INTENT(IN)     :: ptdeep_a, ptdeep_b, pgammat
!                                      PTDEEP_A = Deep soil temperature
!                                                 coefficient depending on flux
!                                      PTDEEP_B = Deep soil temperature (prescribed)
!                                               which models heating/cooling from
!                                               below the diurnal wave penetration
!                                               (surface temperature) depth. If it
!                                               is FLAGGED as undefined, then the zero
!                                               flux lower BC is applied.
!                                      Tdeep = PTDEEP_B + PTDEEP_A * PDEEP_FLUX
!                                              (with PDEEP_FLUX in W/m2)
!                                     PGAMMAT  = Deep soil heat transfer coefficient:
!                                                assuming homogeneous soil so that
!                                                this can be prescribed in units of 
!                                                (1/days): associated time scale with
!                                                PTDEEP.
!
REAL, DIMENSION(:), INTENT(IN)      :: psnowfree_alb_veg  !snow free albedo of vegetation for EBA
REAL, DIMENSION(:), INTENT(IN)      :: psnowfree_alb_soil !snow free albedo of soil for EBA option
REAL, DIMENSION(:), INTENT(IN)      :: ppsnv_a !fraction of the the vegetation covered by snow for EBA scheme
!
REAL, DIMENSION(:), INTENT(IN)   :: pfalb, pfemis, pffrozen
!                                   PFALB = Floodplain albedo
!                                   PFEMIS= Floodplain emis
!
REAL, DIMENSION(:), INTENT(IN)   :: psnowswe
!                                   PSNOWSWE = equivalent water content of
!                                              the D95 snow reservoir (kg m-2)
!
REAL, DIMENSION(:), INTENT(IN)   :: psrsfc  ! Snow rate falling outside of snow
!                                             covered grid area [kg/(m2 s)]
!
!inout
!
REAL, DIMENSION(:,:), INTENT(INOUT) :: ptg
!                                      surface temperature
REAL, DIMENSION(:),   INTENT(INOUT) :: presa
!                                      aerodynamic surface resistance
REAL, DIMENSION(:),   INTENT(INOUT) :: plvtt
!                                      PLVTT=Vaporization heat
REAL, DIMENSION(:),   INTENT(INOUT) :: plstt
!                                      PLSTT=Sublimation heat
REAL, DIMENSION(:),   INTENT(INOUT) :: pcps
!                                      PCPS=specific heat at surface
REAL, DIMENSION(:),   INTENT(INOUT) :: pdelta
!                                      PDELTA = fraction of the foliage covered
!                                              by intercepted water
!
!out
!
REAL, DIMENSION(:), INTENT(OUT)  :: pch, pcd, pcdn, pri
!                                     PCH = drag coefficient for heat
!                                     PCD = drag coefficient for momentum
!                                     PCDN= neutral drag coefficient for momentum
!                                     PRI = Richardson number
!
REAL, DIMENSION(:), INTENT(OUT)  :: phug, phugi, phv, phu, pqs
!                                     PHUG  = ground relative humidity
!                                     PHUGI = ground (ice) relative humidity
!                                     PHV = Halstead coefficient
!                                     PHU = relative humidity at the surface
!                                     PQS = humidity at surface
!
REAL, DIMENSION(:), INTENT(OUT)  :: palbt, pemist
!                                     PALBT  = averaged albedo
!                                     PEMIST = averaged emissivity
!
REAL, DIMENSION(:), INTENT(OUT)   :: pdeep_flux ! Heat flux at bottom of ISBA (W/m2)
!
REAL, DIMENSION(:), INTENT(OUT)     :: prn, ph, ple, pleg, plev, ples
REAL, DIMENSION(:), INTENT(OUT)     :: pler, pletr, pevap, pgflux, pmeltadv, pmelt, prestore
!                                     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)
!                                     PRESTORE = surface restore flux (W m-2)
!
REAL, DIMENSION(:), INTENT(OUT)     :: plegi, pustar
!                                      PLEGI   = sublimation component of the 
!                                                latent heat flux from the soil surface
!                                      PUSTAR  = friction velocity
!
REAL, DIMENSION(:), INTENT(OUT)     :: ple_flood, plei_flood !Floodplains latent heat flux [W/m2]
REAL, DIMENSION(:), INTENT(OUT)     :: psnowtemp  ! snow layer temperatures (K)
!
REAL, DIMENSION(:),  INTENT(OUT) :: pac_agg  ! aggregated aerodynamic conductance
                                     ! for evaporative flux calculations
REAL, DIMENSION(:),  INTENT(OUT) :: phu_agg  ! aggregated relative humidity
                                    ! for evaporative flux calculations
!
!
!*      0.2    declarations of local parameters
!
REAL, PARAMETER            :: zdepth_cor = 0.6
!                             ZDEPTH_COR = depth over which the correction flux is applied
!
REAL, PARAMETER            :: zdtg1_cor = 10.0 
!                             ZDTG1_COR = Delta temperature limit to comput the correction flux (K)
!
!
!*      0.3    declarations of local variables
!
REAL, DIMENSION(SIZE(PTA)) :: zqsat     ! expression for the saturation specific humidity 
!
REAL, DIMENSION(SIZE(PTA)) :: zdqsat    ! expression for the saturation specific humidity derivative
!
REAL, DIMENSION(SIZE(PTA)) :: zta_ic, zqa_ic, zustar2_ic ! TA, QA and friction updated values
!                                                        ! if implicit coupling with atmosphere used.
!
REAL, DIMENSION(SIZE(PTA)) :: zleg_delta  ! soil evaporation delta fn
REAL, DIMENSION(SIZE(PTA)) :: zlegi_delta ! soil sublimation delta fn
!
REAL, DIMENSION(SIZE(PTA)) :: zt2m     ! restore temperature before time integration (K)
REAL, DIMENSION(SIZE(PTA)) :: ztsm     ! surface temperature before time integration (K)
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2)) :: zflux_cor, zlayerhcap
!                                           ZFLUX_COR = correction flux by layer to conserve energy (W/m2)
!
REAL, DIMENSION(SIZE(PTA)) :: zgrndflux, ztotalhcap, zwork
!
INTEGER                    :: ini, inl, ji, jl
LOGICAL                    :: lexplicit_snow
!
!*      0.4    declarations of local time spliting variables
!
! Working arrays for flux averaging over time split
!
REAL, DIMENSION(SIZE(PTA)) :: zdeep_flux, zle_flood, zlei_flood, &
                              zrn, zh, zle, zleg, zlev,  &
                              zles, zler, zletr, zevap,      &
                              zgflux, zmeltadv, zmelt,           &
                              zrestore, zlegi, zustar2,          &
                              zac_agg, zhu_agg
!
REAL, DIMENSION(SIZE(PTA)) :: zdeep_flux_sum, zle_flood_sum, zlei_flood_sum, &
                              zrn_sum, zh_sum, zle_sum, zleg_sum, zlev_sum,  &
                              zles_sum, zler_sum, zletr_sum, zevap_sum,      &
                              zgflux_sum, zmeltadv_sum, zmelt_sum,           &
                              zrestore_sum, zlegi_sum, zustar2_sum,          &
                              zac_agg_sum, zhu_agg_sum
!
REAL, PARAMETER            :: ztsplit  = 300. ! s Minimum time tstep required to time-split energy budget
INTEGER                    :: itsplit, jsplit
REAL                       :: ztstep, znsplit
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ISBA_CEB',0,zhook_handle)
!
!-------------------------------------------------------------------------------
!
!*      1.0    Preliminaries
!              -------------
!
ini=SIZE(ptg,1)
inl=SIZE(ptg,2)
!
!local init
!
zqsat(:) = xundef
zta_ic(:) = xundef
zqa_ic(:) = xundef
zustar2_ic(:) = xundef
zleg_delta(:) = xundef
zlegi_delta(:) = xundef
zdqsat(:) = xundef
!
zdeep_flux(:) = xundef
zle_flood(:) = xundef
zlei_flood(:) = xundef
zrn(:) = xundef
zh(:) = xundef
zle(:) = xundef
zleg(:) = xundef
zlev(:) = xundef
zles(:) = xundef
zler(:) = xundef
zletr(:) = xundef
zevap(:) = xundef
zgflux(:) = xundef
zmeltadv(:) = xundef
zmelt(:) = xundef
zrestore(:) = xundef
zlegi(:) = xundef
zac_agg(:) = xundef
zhu_agg(:) = xundef
!
zustar2_sum(:) = 0.0
zevap_sum(:) = 0.0
!
zrn_sum(:) = 0.0
zh_sum(:) = 0.0
zgflux_sum(:) = 0.0
zle_sum(:) = 0.0
!
zleg_sum(:) = 0.0
zlegi_sum(:) = 0.0
zlev_sum(:) = 0.0
zles_sum(:) = 0.0
zler_sum(:) = 0.0
zletr_sum(:) = 0.0
zle_flood_sum(:) = 0.0
zlei_flood_sum(:) = 0.0
!
zdeep_flux_sum(:) = 0.0
zmeltadv_sum(:) = 0.0
zmelt_sum(:) = 0.0
zrestore_sum(:) = 0.0
zac_agg_sum(:) = 0.0
zhu_agg_sum(:) = 0.0
!
!-------------------------------------------------------------------------------
!
zgrndflux(:  ) = pgrndflux(:)
!
zflux_cor(:,:) = 0.0
!
!-------------------------------------------------------------------------------
!
!*      2.0    Correction flux to conserv energy budget
!              ----------------------------------------
!
lexplicit_snow=(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO')
!
IF(lexplicit_snow.AND.hisba/='DIF')THEN
!
  zflux_cor(:,1)=ppsn(:)*pflsn_cor(:)
!
ELSEIF(lexplicit_snow.AND.hisba=='DIF')THEN
!
  zlayerhcap(:,:) = 0.0
  ztotalhcap(:) = 0.0
!
! To conserv energy, the correction flux is distributed at least
! over the first layers of the soil, ZDEPTH_COR. This method prevent 
! numerical oscillations especially when explicit snow vanishes
!
  zwork(:)=min(pd_g(:,inl),zdepth_cor)
!
  zlayerhcap(:,1)= 1.0/pct(:)
  ztotalhcap(:  )= 1.0/pct(:)
  DO jl=2,inl
     DO ji=1,ini
        zlayerhcap(ji,jl)=psoilhcapz(ji,jl)*min(pdzg(ji,jl),max(0.0,zwork(ji)-pd_g(ji,jl)+pdzg(ji,jl)))
        ztotalhcap(ji   )=ztotalhcap(ji)+zlayerhcap(ji,jl)
     ENDDO
  ENDDO
!
  DO jl=1,inl
     DO ji=1,ini
        IF(ztotalhcap(ji)>0.0)THEN
          zflux_cor(ji,jl)=ppsn(ji)*pflsn_cor(ji)*zlayerhcap(ji,jl)/ztotalhcap(ji)
        ENDIF
    ENDDO
  ENDDO
!
! The second correction is computed if the delta temperature
! due to snow/soil ground flux is superior to ZDTG1_COR (K)
! Especially relevant when PPSN ~ 1 over vegetated area
!
  zwork(:)=ptstep*pct(:)*ppsn(:)*abs(pgrndflux(:))
!
  WHERE(ztotalhcap(:)>0.0.AND.zwork(:)>=zdtg1_cor)
       zgrndflux(:) = pgrndflux(:)*zlayerhcap(:,1)/ztotalhcap(:)
  ENDWHERE
!
  DO jl=2,inl
     DO ji=1,ini
        IF(ztotalhcap(ji)>0.0.AND.zwork(ji)>=zdtg1_cor)THEN
          zflux_cor(ji,jl)=zflux_cor(ji,jl)+ppsn(ji)*pgrndflux(ji) &
                                           *zlayerhcap(ji,jl)/ztotalhcap(ji)
        ENDIF
     ENDDO
  ENDDO
!
ENDIF
!
!-------------------------------------------------------------------------------
!
IF(lcpl_arp)THEN
  itsplit = 1
ELSE
  itsplit  = max(1,nint(ptstep/ztsplit))  ! number of split-time steps
ENDIF
!
znsplit  = REAL(itsplit)
!
ztstep  = ptstep/znsplit             ! split time step
!
!-------------------------------------------------------------------------------
!
DO jsplit=1,itsplit
!
!  Save surface and sub-surface temperature values at beginning of time step for 
!  budget and flux calculations:
!
   ztsm(:) = ptg(:,1)
   zt2m(:) = ptg(:,2)
!
!
!*      3.0    Aerodynamic drag and heat transfer coefficients
!              -----------------------------------------------
!
!  In DRAG, we use the timestep of ISBA (PTSTEP) and not the split time step (ZTSTEP)
!  because diagnostic canopy evaporation (Er) must be consistent with PWR water
!  mass to limit negative dripping in hydro_veg
!
   CALL drag(hisba, hsnow_isba, hcpsurf, ptstep,                                         &
             ptg(:,1), pwg(:,1), pwgi(:,1), pexns, pexna, pta, pvmod, pqa, prr, psr,     &
             pps, prs, pveg, pz0_with_snow, pz0eff, pz0h_with_snow,                      &
             pwfc(:,1), pwsat(:,1), ppsng, ppsnv, pzref, puref,                          &
             pdircoszw, pdelta, pf5, presa,  pch, pcd, pcdn, pri, phug, phugi,           &
             phv, phu, pcps, pqs, pffg, pffv, pff, pffg_nosnow, pffv_nosnow,             &
             zleg_delta, zlegi_delta, pwr, prhoa, plvtt, pqsat=zqsat                     )  
!
!
!
!*      4.0    Resolution of the surface and soil energy budget
!              ------------------------------------------------
!
   CALL e_budget(hisba, hsnow_isba, oflood, otemp_arp, himplicit_wind,                   &
                 ztstep, psodelx, puref,                                                 &
                 ppew_a_coef, ppew_b_coef, ppet_a_coef, ppeq_a_coef, ppet_b_coef,        &
                 ppeq_b_coef, pvmod, pcd, ptg, ztsm, zt2m, psnowalb, psw_rad, plw_rad,   &
                 pta, pqa, pps, prhoa, pexns, pexna, pcps, plvtt, plstt,  pveg,          &
                 phug, phugi, phv, zleg_delta, zlegi_delta, pemis, palb, presa,          &
                 pct, pcg, ppsn, ppsnv, ppsng, zgrndflux, zflux_cor,                     &
                 pd_g, pdzg, pdzdif, psoilcondz, psoilhcapz,  palbt, pemist,             &
                 zqsat, zdqsat, pfrozen1, ptdeep_a, ptdeep_b, pgammat,                   &
                 zta_ic, zqa_ic, zustar2_ic,                                             &
                 psnowfree_alb_veg, ppsnv_a, psnowfree_alb_soil,                         &
                 pffg, pffv, pff, pffrozen, pfalb, pfemis, zdeep_flux, zrestore          )
!
!
!*      5.0    Energy and momentum fluxes
!              --------------------------
!
!*******************************************************************************
! WARNING: at this stage, ZALBT and ZEMIST have two different meanings according
!          to the ISBA snow-scheme option:
!  'D95' : they represent aggregated (snow + flood + snow-flood-free) albedo and emissivity
!  '3-L' : they represent                    flood + snow-flood-free  albedo and emissivity
!*******************************************************************************
!
   CALL isba_fluxes(hisba, hsnow_isba, otemp_arp, ztstep, psodelx,                       &
                    psw_rad, plw_rad, zta_ic, zqa_ic,                                    &
                    prhoa, pexns, pexna, pcps, plvtt, plstt,                             &
                    pveg, phug, phugi, phv, zleg_delta, zlegi_delta, pdelta, presa,      &
                    pf5, prs, pcs, pcg, pct, psnowswe, ztsm, zt2m,                       &
                    ppsn, ppsnv, ppsng, pfrozen1,                                        &
                    palbt, pemist, zqsat, zdqsat, psnow_thrufal,                         &
                    zrn, zh, zle, zleg, zlegi, zlev, zles, zler, zletr, zevap, zgflux,   &
                    zmeltadv, zmelt,                                                     &
                    psoilcondz,  pd_g, pdzg, ptg,                                        &
                    psrsfc, ppsnv_a, pffg, pffv, pff, pffrozen,                          &
                    zle_flood, zlei_flood, psnowtemp                                     ) 
!
!
!*      6.0    Aggregated coefficients
!              -----------------------
!
!  Compute aggregated coefficients for evaporation
!  Sum(LEV+LEG+LEGI+LES) = ACagg * Lv * RHOA * (HUagg.Qsat - Qa)
!
   zac_agg(:) =   1. / presa(:) / plvtt(:)                                &
        * ( plvtt(:)*    pveg(:) *(1.-ppsnv(:))                 *phv(:)   &
          + plvtt(:)*(1.-pveg(:))*(1.-ppsng(:))*(1.-pfrozen1(:))          &
          + plstt(:)*(1.-pveg(:))*(1.-ppsng(:))*    pfrozen1(:)           &
          + plstt(:)*                 ppsn(:)                            )  
!
   WHERE(zac_agg(:)>0.0)
         zhu_agg(:) =   1. / (presa(:) * zac_agg(:)) / plvtt(:)                 &
              * ( plvtt(:)*    pveg(:) *(1.-ppsnv(:))                 *phv(:)   &
                + plvtt(:)*(1.-pveg(:))*(1.-ppsng(:))*(1.-pfrozen1(:))*phug(:)  &
                + plstt(:)*(1.-pveg(:))*(1.-ppsng(:))*    pfrozen1(:) *phugi(:) &
                + plstt(:)*                 ppsn(:)                            )  
   ENDWHERE
!
   zustar2_sum(:) = zustar2_sum(:) + zustar2_ic(:)
!
   zevap_sum(:) = zevap_sum(:) + zevap(:)
!
   zrn_sum(:) = zrn_sum(:) + zrn(:)
   zh_sum(:) = zh_sum(:) + zh(:)
   zgflux_sum(:) = zgflux_sum(:) + zgflux(:)
   zle_sum(:) = zle_sum(:) + zle(:)
!
   zleg_sum(:) = zleg_sum(:) + zleg(:)
   zlegi_sum(:) = zlegi_sum(:) + zlegi(:)   
   zlev_sum(:) = zlev_sum(:) + zlev(:)
   zles_sum(:) = zles_sum(:) + zles(:)
   zler_sum(:) = zler_sum(:) + zler(:)
   zletr_sum(:) = zletr_sum(:) + zletr(:)
   zle_flood_sum(:) = zle_flood_sum(:) + zle_flood(:)
   zlei_flood_sum(:) = zlei_flood_sum(:) + zlei_flood(:)   
!
   zdeep_flux_sum(:) = zdeep_flux_sum(:) + zdeep_flux(:)
   zmeltadv_sum(:) = zmeltadv_sum(:) + zmeltadv(:)
   zmelt_sum(:) = zmelt_sum(:) + zmelt(:)
   zrestore_sum(:) = zrestore_sum(:) + zrestore(:)
   zac_agg_sum(:) = zac_agg_sum(:) + zac_agg(:)
   zhu_agg_sum(:) = zhu_agg_sum(:) + zhu_agg(:)
!
!-------------------------------------------------------------------------------
!
ENDDO
!
pustar(:) = sqrt(zustar2_sum(:)/znsplit)
!
pevap(:) = zevap_sum(:) / znsplit
!
prn(:) = zrn_sum(:) / znsplit
ph(:) = zh_sum(:) / znsplit
pgflux(:) = zgflux_sum(:) / znsplit
ple(:) = zle_sum(:) / znsplit
!
pleg(:) = zleg_sum(:) / znsplit
plegi(:) = zlegi_sum(:) / znsplit
plev(:) = zlev_sum(:) / znsplit
ples(:) = zles_sum(:) / znsplit
pler(:) = zler_sum(:) / znsplit
pletr(:) = zletr_sum(:) / znsplit
ple_flood(:) = zle_flood_sum(:) / znsplit
plei_flood(:) = zlei_flood_sum(:) / znsplit
!
pdeep_flux(:) = zdeep_flux_sum(:) / znsplit
pmeltadv(:) = zmeltadv_sum(:) / znsplit
pmelt(:) = zmelt_sum(:) / znsplit
prestore(:) = zrestore_sum(:) / znsplit
pac_agg(:) = zac_agg_sum(:) / znsplit
phu_agg(:) = zhu_agg_sum(:) / znsplit
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('ISBA_CEB',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE isba_ceb