e_budget.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 e_budget(HISBA, HSNOW_ISBA, OFLOOD, OTEMP_ARP, HIMPLICIT_WIND,  &
                            ptstep, psodelx, puref, ppew_a_coef, ppew_b_coef,    &
                            ppet_a_coef, ppeq_a_coef, ppet_b_coef, ppeq_b_coef,  &
                            pvmod, pcd,                                          &
                            ptg, ptsm, pt2m, psnowalbm,                          &
                            psw_rad, plw_rad, pta, pqa, pps, prhoa,              &
                            pexns, pexna, pcps, plvtt, plstt,                    &
                            pveg, phug, phui, phv,                               &
                            pleg_delta, plegi_delta,                             &
                            pemis, palb, pra,                                    &
                            pct, pcg, ppsn, ppsnv, ppsng,                        &
                            pgrndflux, pflux_cor,                                &
                            pd_g, pdzg, pdzdif, psoilcondz, psoilhcapz,          &
                            palbt, pemist, pqsat, pdqsat,                        &
                            pfrozen1, ptdeep_a, ptdeep_b, pgammat,               &
                            pta_ic, pqa_ic, pustar2_ic,                          &
                            psnowfree_alb_veg, ppsnv_a,psnowfree_alb_soil,       &
                            pffg, pffv, pff, pffrozen, pfalb, pfemis, pdeep_flux,&
                            prestore                                             )  
!     ##########################################################################
!
!!****  *E_BUDGET*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the evolution of the surface and deep-soil temperature
!     (i.e., Ts and T2), as well as all the surface fluxes.
!         
!     
!!**  METHOD
!!    ------
!
!     1- find the grid-averaged albedo, emissivity, and roughness length
!     2- compute the za, zb, and zc terms involved in the numerical
!        resolution of the equations for Ts and T2.
!     3- find Ts(t) and T2(t).
!     4- derive the surface fluxes.
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!      
!!    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)  28/08/98 splits the routine in two for C02 computations
!!      (A.Boone)   15/03/99 Soil ice tendencies calculated here: heating/cooling
!!                           affects surface and deep soil temperatures.
!!      (A. Boone, V. Masson) 01/2003 Externalization
!!      (E. Martin)          07/05 implicit coupling (coeff ZA,ZB,ZC)
!!      (P. Le Moigne)       07/05 dependence on qs for cp
!!      (B. Decharme)        05/08 Add floodplains dependencies
!!      (B. Decharme)        01/09 optional deep soil temperature as in Arpege
!!      (R. Hamdi)           01/09 Cp and L are not constants (As in ALADIN)
!!      (B. Decharme)        09/09 When LCPL_ARP, do not calculate x2 each coef
!!      (A.Boone)            03/10 Add delta fnctions to force LEG ans LEGI=0
!!                                 when hug(i)Qsat < Qa and Qsat > Qa
!!      (B. Decharme)        09/12 new wind implicitation
!!      (V. Masson)          01/13 Deep soil flux implicitation
!!      (B. Decharme)        10/14 Bug in DIF composite budget
!!                                 Use harmonic mean to compute interfacial thermal conductivities
!!                                 "Restore" flux computed here
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,       ONLY : xlvtt, xlstt, xstefan, xcpd, xpi, xday, &
                              xtt, xcl, xcpv, xci  
USE modd_surf_par,   ONLY : xundef
USE modd_snow_par,   ONLY : xemissn, xemcrin
!
USE modd_surf_atm,   ONLY : lcpl_arp, lqvnplus
!
USE mode_thermos
!
USE modi_soil_heatdif 
USE modi_soil_temp_arp
!
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 2000)
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)   :: puref       ! reference height of the wind
REAL, DIMENSION(:), INTENT(IN)   :: psnowalbm
!                                     prognostic variables at time 't-dt'
!                                     PSNOWALBM = albedo of the snow
!
!
REAL, DIMENSION(:), INTENT (IN)  :: psw_rad, plw_rad, pps, prhoa, pta, pqa, pcd, pvmod
!                                     PSW_RAD = incoming solar radiation
!                                     PLW_RAD = atmospheric infrared radiation
!                                     PRHOA = near-ground air density
!                                     PPS = surface pressure
!                                     PTA = near-ground air temperature
!                                     PQA = near-ground air specific humidity
!                                     PCD = drag coefficient
!                                     PVMOD = wind speed
!
! implicit atmospheric coupling coefficients:
!
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)   :: pexns, pexna
REAL, DIMENSION(:), INTENT(IN)   :: pveg, phug, phui, phv
REAL, DIMENSION(:), INTENT(IN)   :: pemis, palb, pct, pcg
REAL, DIMENSION(:), INTENT(IN)   :: ppsnv, ppsng, ppsn
!                                     PVEG = fraction of vegetation
!                                     PHUG = relative humidity of the soil
!                                     PHV = Halstead coefficient
!                                     PEMIS = emissivity
!                                     PALB = albedo
!                                     PCT = area-averaged heat capacity
!                                     PCG = heat capacity of the ground
!                                     PPSN = grid fraction covered by snow
!                                     PPSNV = fraction of the vegetation covered by snow
!                                     PPSNG = fraction of the ground covered by snow 
!
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)      :: pgrndflux 
!                                      PGRNDFLUX = soil/snow interface flux (W/m2) using
!                                                  ISBA-SNOW3L option
!
REAL, DIMENSION(:,:), INTENT(IN)    :: pflux_cor
!                                      PFLUX_COR = correction flux to conserve energy (W/m2)
!
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)      :: 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(INOUT)   :: pleg_delta, plegi_delta
!                                      PLEG_DELTA = soil evaporation delta fn
!                                      PLEGI_DELTA = soil evaporation delta fn
!
REAL, DIMENSION(:), INTENT (OUT)   :: pqa_ic, pta_ic, pustar2_ic
!                                     PTA_IC = near-ground air temperature
!                                     PQA_IC = near-ground air specific humidity
!                                     PUSTAR2_IC = near-ground wind friction (m2/s2)
!                                           (modified if implicit coupling with
!                                            atmosphere used)
!
REAL, DIMENSION(:), INTENT(IN)   :: pra
!                                     PRA = aerodynamic surface resistance for
!                                           heat transfers
!
REAL, DIMENSION(:,:), INTENT(INOUT):: ptg
!                                     PTG    = soil temperature profile (K)
!
REAL, DIMENSION(:), INTENT (IN)     :: ptsm, pt2m
!                                      PTSM   = surface temperature at start 
!                                               of time step (K)
!                                      PT2M   = mean surface (or restore) temperature at start 
!                                               of time step (K)
!
REAL, DIMENSION(:), INTENT(INOUT)  :: pcps
!                                     PCPS   =  heat capacity at surface
!
REAL, DIMENSION(:), INTENT(OUT)  :: palbt, pemist, pdqsat
!                                     PALBT  = averaged albedo
!                                     PEMIST = averaged emissivity
!                                     PDQSAT = saturation vapor humidity derivative
REAL, DIMENSION(:), INTENT(IN)   :: pqsat
!                                     PQSAT  = saturation vapor humidity
!
REAL, DIMENSION(:), INTENT(IN)   :: pffv, pff, pffg, pfalb, pfemis, pffrozen
!                                   PFFG = Floodplain fraction over ground
!                                   PFFV = Floodplain fraction over vegetation
!                                   PFF  = Floodplain fraction at the surface
!                                   PFALB = Floodplain albedo
!                                   PFEMIS= Floodplain emis
!
REAL, DIMENSION(:), INTENT(INOUT)  :: plstt, plvtt
!
REAL, DIMENSION(:), INTENT(OUT)    :: pdeep_flux ! Heat flux at bottom of ISBA (W/m2)
!
REAL, DIMENSION(:), INTENT(OUT)    :: prestore
!                                     PRESTORE = surface restore flux (W m-2)
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PALB)) ::   zrora,                        &
!                                             rhoa / ra
!
                                  za,zb,zc  
!                                             terms for the calculation of Ts(t)
!
! ISBA-DF:
!
REAL, DIMENSION(SIZE(PALB)) ::  zcondavg, zcond1, zcond2, zterm2, zterm1
!
! implicit atmospheric coupling coefficients: (modified-form)
!
REAL, DIMENSION(SIZE(PALB)) :: zpet_a_coef, zpeq_a_coef, zpet_b_coef,      &
                               zpeq_b_coef, z_ccoef, zhums, zhuma, zlavg,  &
                               zhumsd, zhumad 
!                              ZPET_A_COEF = A-air temperature coefficient
!                              ZPET_B_COEF = B-air temperature coefficient
!                              ZPEQ_A_COEF = A-air specific humidity coefficient
!                              ZPEQ_B_COEF = B-air specific humidity coefficient
!                              Z_CCOEF     = C-working variable
!
REAL, DIMENSION(SIZE(PALB)) :: zustar2, zvmod
!                              ZUSTAR2 = friction     (m2/s2)
!                              ZVMOD   = wind modulus (m/s)
REAL, DIMENSION(SIZE(PALB)) :: zxcpv_xcl_avg
REAL, DIMENSION(SIZE(PALB)) :: zcnhuma, zpeqa2, zdpqb, zcdqsat, zincr, ztrad, &
                                zchums, zchuma, zpeta2, zpetb2,ztemp, zfgnfrz, &
                                zfgfrz, zfv, zfg, zfnfrz, zffrz, zfnsnow, zcps,&
                                zlvtt, zlstt
REAL                        :: zsnow
!
REAL(KIND=JPRB) :: zhook_handle
!
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('E_BUDGET',0,zhook_handle)
!
zcondavg(:)  = 0.0
zterm2(:)    = 0.0
zterm1(:)    = 0.0
zhumsd(:)    = 0.0
zhumad(:)    = 0.0
!
!-------------------------------------------------------------------------------
!
!*       1.     COEFFICIENTS FOR THE TIME INTEGRATION OF  TS 
!               --------------------------------------------
!
!
!                                              function dqsat(Ts,ps)
!
pdqsat(:) = dqsat(ptg(:,1),pps(:),pqsat(:))
!                                              function zrsra
!
! Modify flux-form implicit coupling coefficients:
! - wind components:
!
ztemp(:) = pcd(:)*pvmod(:)
!
IF(himplicit_wind=='OLD')THEN
! old implicitation (m2/s2)
  zustar2(:) = ztemp(:) * ppew_b_coef(:) / (1.0- ztemp(:)*prhoa(:)*ppew_a_coef(:)) 
ELSE
! new implicitation (m2/s2)
  zustar2(:) = ztemp(:) * (2.*ppew_b_coef(:)-pvmod(:)) / (1.0-2.0*ztemp(:)*prhoa(:)*ppew_a_coef(:)) 
ENDIF
!
!wind modulus at t+1 (m/s)
zvmod(:) = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
zvmod(:) = max(zvmod(:),0.)
!
WHERE(ppew_a_coef(:)/= 0.)
      zustar2(:) = max( ( zvmod(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0.)
ENDWHERE
!
zustar2(:) = max(zustar2(:),0.)
!
zrora(:)    = prhoa(:) / pra(:)
!
!                                              terms za, zb, and zc for the
!                                              calculation of ts(t)
!
! Modify flux-form implicit coupling coefficients:
! - air temperature:
!
ztemp(:) = ppet_a_coef(:)*zrora(:)
z_ccoef(:)     = (1.0 - ztemp(:))/pexna(:)
!
zpet_a_coef(:) = - ztemp(:)/pexns(:)/z_ccoef(:)
!
zpet_b_coef(:) = ppet_b_coef(:)/z_ccoef(:) 
!
!-------------------------------------------------------------------------------
!
!*       2.     AIR AND SOIL SPECIFIC HUMIDITIES 
!               --------------------------------
!
! - air specific humidity:
!
zfv(:) = pveg(:) * (1-ppsnv(:)-pffv(:))
zfg(:) = (1.-pveg(:))*(1.-ppsng(:)-pffg(:))
zfnfrz(:) = (1.-pffrozen(:))*pff(:) + zfv + zfg(:)*(1.-pfrozen1(:))
zffrz(:) = pffrozen(:)*pff(:) + zfg(:)*pfrozen1(:) + ppsn(:)
!
zsnow=1.
zfnsnow(:)=1.
zcps(:)=pcps(:)
!
IF (lcpl_arp) THEN

  ! currently this correction not applied for this option, but can be
  ! added later after testing...so delta fns set to 1 (turns OFF this correction)

  pleg_delta(:)  = 1.0
  plegi_delta(:) = 1.0

  zlavg(:)     = plvtt(:)*zfnfrz(:) + plstt(:)*zffrz(:)
  zxcpv_xcl_avg(:)= (xcpv-xcl)*zfnfrz(:) + (xcpv-xci)*zffrz(:) 

  zlvtt(:) = zlavg(:)
  zlstt(:) = zlavg(:)

ELSE

  IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
    zsnow = 0.
    zfnsnow(:) = 1. - ppsn(:)
    zcps(:)=xcpd
  ENDIF

  zlavg(:)     = xlvtt*zfnfrz(:) + xlstt*zffrz(:)

  zlvtt(:) = xlvtt
  zlstt(:) = xlstt

ENDIF
!
zfgnfrz(:) = zfg(:)*(1.-pfrozen1(:))*pleg_delta(:)
zfgfrz(:) = zfg(:)*pfrozen1(:)*plegi_delta(:)
!
zhuma(:)     = zlvtt(:)/zlavg(:) * ((1.-pffrozen(:))*pff(:) + zfv(:)*phv(:) + zfgnfrz(:))   +         &
                zlstt(:)/zlavg(:) * (pffrozen(:)*pff(:) + zfgfrz(:) + zsnow*ppsn(:) ) 
!
zhums(:)     = zlvtt(:)/zlavg(:) * ((1.-pffrozen(:))*pff(:) + zfv(:)*phv(:) + zfgnfrz(:)*phug(:)) +   &
                zlstt(:)/zlavg(:) * (pffrozen(:)*pff(:) + zfgfrz(:)*phui(:) + zsnow*ppsn(:) )  
!
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
!
! humidity considering no snow (done elsewhere) and flooded zones:
!
  zhumad(:) = pff(:) + zfv(:)*phv(:) + zfgnfrz(:) + zfgfrz(:)  
  zhumsd(:) = pff(:) + zfv(:)*phv(:) + zfgnfrz(:)*phug(:) + zfgfrz(:)*phui(:)          
ELSE
  zhumad(:) = zhuma(:)
  zhumsd(:) = zhums(:)
ENDIF
!
!-------------------------------------------------------------------------------
!
!*       3.     COEFFICIENTS FOR THE TIME INTEGRATION OF Q 
!               -------------------------------------------
!
! implicit q coefficients:
!
ztemp(:) = ppeq_a_coef(:)*zrora(:)
z_ccoef(:)     = 1.0 - ztemp(:)*zhumad(:)
!
zpeq_a_coef(:) = - ztemp(:)*pdqsat(:)*zhumsd(:)/z_ccoef(:)
!
zpeq_b_coef(:) = ( ppeq_b_coef(:) - ztemp(:)*zhumsd(:)* &
                  (pqsat(:) - pdqsat(:)*ptg(:,1)) )/z_ccoef(:)  
!
!-------------------------------------------------------------------------------
!
!*       4.     TOTAL ALBEDO AND EMISSIVITY 
!               ---------------------------
!
!
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN
!
! NON-SNOW covered Grid averaged albedo and emissivity for explicit
! snow scheme 
!
  IF(.NOT.oflood)THEN
!          
     palbt(:) = palb(:)
     pemist(:) = pemis(:)
!     
  ELSE
!
! Taking into account the floodplains with snow grid fractions :
!     PFF    1.-PFF-PPSN   PPSN
! |------------|----|---------------|
!
  WHERE(ppsn(:)<1.0)          
     palbt(:) = ((1.-pff(:)-ppsn(:))*palb(:)  + pff(:)*pfalb(:))/(1.-ppsn(:))
     pemist(:) = ((1.-pff(:)-ppsn(:))*pemis(:) + pff(:)*pfemis(:))/(1.-ppsn(:))
  ELSEWHERE
     palbt(:) = palb(:)
     pemist(:) = pemis(:)
  ENDWHERE
!
  ENDIF
!
!
ELSE
!
! Grid averaged albedo and emissivity for composite snow scheme:
!
   IF(hsnow_isba=='EBA') THEN
!
      palbt(:)  = (1-pveg(:))*(psnowfree_alb_soil(:)*(1-ppsng(:))+psnowalbm(:)*ppsng(:)) + &
                    pveg(:)*(psnowfree_alb_veg(:)*(1-ppsnv_a(:)) +             &
                    psnowalbm(:)*ppsnv_a(:))  
!
      pemist(:) = pemis(:)-ppsn(:)*(pemis(:)-xemcrin)
!      
   ELSE
!
      palbt(:) = ( 1.-ppsn(:)-pff(:))*palb(:) + ppsn(:)*psnowalbm(:) + pff(:)*pfalb(:)     
!
      pemist(:) = ( 1.-ppsn(:)-pff(:))*pemis(:) + ppsn(:)*xemissn + pff(:)*pfemis(:)
!
   ENDIF

ENDIF
!
!-------------------------------------------------------------------------------
!
!*       5. CALCULATION OF ZA, ZB, ZC
!       -----------------------------
!
!       5.1. Default
!       ------------
!
ztrad(:) =  pemist(:) * xstefan * (ptg(:,1)**3)
zchums(:) = zrora(:)*zlavg(:)*zhums(:)
zchuma(:) = zrora(:)*zlavg(:)*zhuma(:)
!
zpeta2(:) = 1./pexns(:) - zpet_a_coef(:)/pexna(:)
zpetb2(:) = zpet_b_coef(:)/pexna(:)
!
! Surface Energy Budget linearization coefficients for an explicit
! soil-flood-vegetation energy budget with an insulating fractional overlying
! layer of snow: fluxes partitioned between surface "felt" by atmosphere
! and surface in contact with base of snowpack (flux exchange between
! atmosphere and snow surface calculated in explicit snow routine)
! (Boone and Etchevers, 2001, J Hydromet.)
! NOTE for now, the meltwater advection term (heat source/sink)
! is OFF because the corresponding energy should be compensated for
! (but code is retained for possible future activation).
!
za(:) = 1. / ptstep + pct(:) *                            &
         ((zfnsnow(:) *                                    &
           ( 4.*ztrad(:) + zrora(:)*zcps(:)*zpeta2(:) ))   &
         + zchums(:)*pdqsat(:) - zchuma(:)*zpeq_a_coef(:)) &
         + 2. * xpi / xday        
!
zb(:) = 1. / ptstep + pct(:) * ( zfnsnow(:)* 3.*ztrad(:) + zchums(:)*pdqsat(:) ) 
!
zc(:) = 2. * xpi * ptg(:,2) / xday + pct(:) *              &
       ( zfnsnow(:) *                                       &
       ( zrora(:)*zcps(:)*zpetb2(:)                         &
       + psw_rad(:)*(1.-palbt(:)) + plw_rad(:)*pemist(:))   &
       - (zchums(:)*pqsat(:) - zchuma(:)*zpeq_b_coef(:)))           
!
IF(hsnow_isba == '3-L' .OR. hsnow_isba == 'CRO' .OR. hisba == 'DIF')THEN                                 
!
!       5.2. With CSNOW=SNOW3L or CSNOW=CRO or HISBA=DIF
!       -------------------------------------------------
!
   zc(:) = zc(:) + pct(:)*(ppsn(:)*pgrndflux(:)+pflux_cor(:,1))
!
ELSEIF (lcpl_arp) THEN
!
!       5.3. With Arpege
!       ----------------
!
zcdqsat(:) = (xcpv-xcpd)*zhums(:)*pdqsat(:)
zincr(:)= pct(:) * zrora(:) * &
            (zcdqsat(:) * ( zpeta2(:)*ptg(:,1) - zpetb2(:)) + &
            zxcpv_xcl_avg(:) * &
            (zhums(:)*pqsat(:) - zhuma(:) * (zpeq_b_coef(:) + zpeq_a_coef(:) * ptg(:,1)))) 

! Surface Energy Budget linearization coefficients for a composite 
! (soil-vegetation-flood-snow) energy budget: composite fluxes "felt" by
! atmosphere from a mixed soil,snow and vegetation surface:
! (Douville et al. 1995, J. Clim. Dyn.)
!

  za(:) = za(:) + zincr(:)

  zb(:) = zb(:) + zincr(:)            
           
  IF (lqvnplus) THEN
!
!       5.4. With  LQVNPLUS=TRUE
!       ------------------------
!
    zcnhuma(:)=(xcpv-xcpd)*(1.-zhuma(:))
    zpeqa2(:)=zcnhuma(:)*zpeq_a_coef(:)*zpeta2(:)*ptg(:,1)
    zdpqb(:)=zpeq_b_coef(:)-pqa(:)

    za(:) = za(:) + pct(:) * zrora(:) * &
             (2.* zpeqa2(:) + &
             zcnhuma(:) * (zdpqb(:)*zpeta2(:) - zpeq_a_coef(:)*zpetb2(:) )) 

    zb(:) = zb(:) + pct(:) * zrora(:) * zpeqa2(:)         
  
    zc(:) = zc(:) + pct(:)*zrora(:)*zcnhuma(:) *zdpqb(:)*zpetb2(:) 
           
  ENDIF
ENDIF
!

!-------------------------------------------------------------------------------
!
!*       6.     T AT TIME 'T+DT' (before snowmelt or soil ice evolution)
!               -----------------
!
IF(hisba == 'DIF')THEN                                                          
!
! First determine terms needed for implicit linearization of surface:
!
!  We use harmonic mean to compute the thermal conductivity at the layers interface
!
   zcond1(:) = pdzg(:,1)/((pdzg(:,1)+pdzg(:,2))*psoilcondz(:,1))
   zcond2(:) = pdzg(:,2)/((pdzg(:,1)+pdzg(:,2))*psoilcondz(:,2))
!
   zcondavg(:) = 1.0/(zcond1(:)+zcond2(:))
!   
   za(:)       = za(:) - (2. * xpi / xday) + zcondavg(:)*pcg(:)/pdzdif(:,1)
   zterm2(:)   = zcondavg(:)*pcg(:)/(za(:)*pdzdif(:,1))
   zterm1(:)   = (ptg(:,1)*zb(:) + (zc(:) - (2. * xpi * ptg(:,2) / xday)) )/za(:)  
!
! Determine the soil temperatures:
!
   CALL soil_heatdif(ptstep,pdzg,pdzdif,psoilcondz,      &
                     psoilhcapz,pcg,zterm1,zterm2,       &
                     ptdeep_a,ptdeep_b,ptg,pdeep_flux,   &
                     pflux_cor                           )
!
!
! "Restore" flux here is actually the heat flux between the surface
! and sub-surface layers (W m-2):
!
   prestore(:) = zcondavg(:)*(ptg(:,1)-ptg(:,2))/pdzdif(:,1)
!
ELSE
!
   IF(otemp_arp)THEN
!
      CALL soil_temp_arp(ptstep,za,zb,zc,pgammat,ptdeep_b,psodelx,ptg)
!
!     "Restore" flux between surface and deep layer(W m-2):
      prestore(:)=2.0*xpi*(ptg(:,1)-ptg(:,2))/(pct(:)*xday*psodelx(1)*(psodelx(1)+psodelx(2)))
!      
   ELSE
!
      ptg(:,1) = ( ptg(:,1)*zb(:) + zc(:) ) / za(:)
!
      WHERE(ptdeep_b(:) /= xundef .AND. pgammat(:) /= xundef)
            ptg(:,2) = (ptg(:,2) + (ptstep/xday)*(ptg(:,1) + pgammat(:)*ptdeep_b(:)))/ &
                         (1.+(ptstep/xday)*(1.0+pgammat(:)))  
      ELSEWHERE
            ptg(:,2) = (ptg(:,2) + (ptstep/xday)*ptg(:,1))/                          &
                         (1.+(ptstep/xday) )  
      END WHERE
!
!     "Restore" flux between surface and deep layer(W m-2):
      prestore(:) = 2.0*xpi*(ptg(:,1)-pt2m(:))/(pct(:)*xday)  
!
   ENDIF
!
ENDIF
!
!-------------------------------------------------------------------------------
!*       7.     TA and QA AT TIME 'T+DT' 
!               ------------------------
!               (QA and TA are only modified by these expressions
!                if the implicit atmospheric coupling is used)
!
pqa_ic(:) =  zpeq_a_coef(:)*ptg(:,1)   + zpeq_b_coef(:)
!
pta_ic(:) =  zpet_a_coef(:)*ptg(:,1)   + zpet_b_coef(:)
!
pustar2_ic(:) =  zustar2(:)
!
!--------------------------------------------------------------------------------------
!*       8.     Update of LSTT and LVTT for Arpege
!               ----------------------------------
!
IF (lcpl_arp) THEN

  IF (.NOT.lqvnplus) THEN
    pcps(:) =  pcps(:) + (xcpv-xcpd) *zhums(:)*pdqsat(:)*(ptg(:,1)-ptsm(:))
  ENDIF


  IF (lqvnplus) THEN
    pcps(:) =  pcps(:) + (xcpv-xcpd) *zhums(:)*pdqsat(:)*(ptg(:,1)-ptsm(:))  &
                       + (xcpv-xcpd) *(1-zhuma(:))*(pqa_ic(:)-pqa(:))  
  ENDIF

  plstt(:) = plstt(:) + (xcpv-xci)*(ptg(:,1)-ptsm(:))

  plvtt(:) = plvtt(:) + (xcpv-xcl)*(ptg(:,1)-ptsm(:))


ENDIF
!
IF (lhook) CALL dr_hook('E_BUDGET',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE e_budget