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(IO, KK, PK, PEK, DK, DMK, HIMPLICIT_WIND,  &
                          PTSTEP, PUREF, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, &
                          PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PVMOD, PTSM, PT2M, &
                          PSW_RAD, PLW_RAD, PTA, PQA, PPS, PRHOA, PEXNS, PEXNA, &
                          PHUI, PLEG_DELTA, PLEGI_DELTA, PGRNDFLUX, PFLUX_COR, &
                          PSOILCONDZ, PSOILHCAPZ, PALBT, PEMIST, PQSAT, PDQSAT, &
                          PFROZEN1, PTDEEP_A,PTA_IC, PQA_IC, PUSTAR2_IC, 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 (199:)
!!      
!!    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_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_k_t, isba_p_t, isba_pe_t
USE modd_diag_n, ONLY : diag_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
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
!
!
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_t), INTENT(INOUT) :: DK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
 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)   :: PUREF       ! reference height of the wind
REAL, DIMENSION(:), INTENT (IN)  :: PSW_RAD, PLW_RAD, PPS, PRHOA, PTA, PQA, 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
!                                     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)   :: PHUI
!
REAL, DIMENSION(:), INTENT(IN)     :: PFROZEN1
!                                     PFROZEN1 = ice fraction in supurficial soil
REAL, DIMENSION(:), INTENT(IN)     :: PTDEEP_A
!                                      PTDEEP_A = Deep soil temperature
!                                                 coefficient depending on flux
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)    :: PSOILCONDZ, PSOILHCAPZ
!                                      PSOILCONDZ= ISBA-DF Soil conductivity profile  [W/(m K)]
!                                      PSOILHCAPZ=ISBA-DF Soil heat capacity profile [J/(m3 K)]
!
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)     :: 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(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(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(PEK%XSNOWFREE_ALB)) :: ZRORA,                        &
!                                             rhoa / ra
!
                                               za,zb,zc  
!                                             terms for the calculation of Ts(t)
!
! ISBA-DF:
!
REAL, DIMENSION(SIZE(PEK%XSNOWFREE_ALB)) ::  ZCONDAVG, ZCOND1, ZCOND2, ZTERM2, ZTERM1
!
! implicit atmospheric coupling coefficients: (modified-form)
!
REAL, DIMENSION(SIZE(PEK%XSNOWFREE_ALB)) :: 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(PEK%XSNOWFREE_ALB)) :: ZUSTAR2, ZVMOD
!                              ZUSTAR2 = friction     (m2/s2)
!                              ZVMOD   = wind modulus (m/s)
REAL, DIMENSION(SIZE(PEK%XSNOWFREE_ALB)) :: ZXCPV_XCL_AVG
REAL, DIMENSION(SIZE(PEK%XSNOWFREE_ALB)) :: ZCNHUMA, ZPEQA2, ZDKQB, 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(pek%XTG(:,1),pps(:),pqsat(:))
!                                              function zrsra
!
! Modify flux-form implicit coupling coefficients:
! - wind components:
!
ztemp(:) = dk%XCD(:)*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(:) / pek%XRESA(:)
!
!                                              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(:) = pek%XVEG(:)     *(1-pek%XPSNV(:)-kk%XFFV(:))
zfg(:) = (1.-pek%XVEG(:))*(1.-pek%XPSNG(:)-kk%XFFG(:))
zfnfrz(:) = (1.-kk%XFFROZEN(:))*kk%XFF(:) + zfv + zfg(:)*(1.-pfrozen1(:))
zffrz(:)  = kk%XFFROZEN(:)     *kk%XFF(:) +       zfg(:)*pfrozen1(:)      + pek%XPSN(:)
!
zsnow      = 1.
zfnsnow(:) = 1.
zcps(:)    = pk%XCPS(:)
!
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(:)        = pk%XLVTT(:)*zfnfrz(:) + pk%XLSTT(:)*zffrz(:)
  zxcpv_xcl_avg(:)=   (xcpv-xcl)*zfnfrz(:) + (xcpv-xci)  *zffrz(:) 

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

ELSE

  IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN
    zsnow = 0.
    zfnsnow(:) = 1. - pek%XPSN(:)
    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.-kk%XFFROZEN(:))*kk%XFF(:) + zfv(:)*dmk%XHV(:) + zfgnfrz(:))   +  &
           zlstt(:)/zlavg(:) * (kk%XFFROZEN(:)     *kk%XFF(:) + zfgfrz(:) + zsnow*pek%XPSN(:) ) 
!
zhums(:) = zlvtt(:)/zlavg(:) * ((1.-kk%XFFROZEN(:))*kk%XFF(:) + zfv(:)*dmk%XHV(:) + zfgnfrz(:)*dk%XHUG(:)) +   &
           zlstt(:)/zlavg(:) * (kk%XFFROZEN(:)     *kk%XFF(:) + zfgfrz(:)*phui(:)  + zsnow*pek%XPSN(:) )  
!
IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN
!
! humidity considering no snow (done elsewhere) and flooded zones:
!
  zhumad(:) = kk%XFF(:) + zfv(:)*dmk%XHV(:) + zfgnfrz(:) + zfgfrz(:)  
  zhumsd(:) = kk%XFF(:) + zfv(:)*dmk%XHV(:) + zfgnfrz(:)*dk%XHUG(:) + 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(:)*pek%XTG(:,1)) )/z_ccoef(:)  
!
!-------------------------------------------------------------------------------
!
!*       4.     TOTAL ALBEDO AND EMISSIVITY 
!               ---------------------------
!
!
IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN
!
! NON-SNOW covered Grid averaged albedo and emissivity for explicit
! snow scheme 
!
  IF(.NOT.io%LFLOOD)THEN
!          
     palbt (:) = pek%XSNOWFREE_ALB(:)
     pemist(:) = pek%XEMIS(:)
!     
  ELSE
!
! Taking into account the floodplains with snow grid fractions :
!     PFF    1.-PFF-PEK%XPSN(:)   PEK%XPSN(:)
! |------------|----|---------------|
!
  WHERE(pek%XPSN(:)<1.0)          
     palbt (:) = ((1.-kk%XFF(:)-pek%XPSN(:))*pek%XSNOWFREE_ALB(:) + kk%XFF(:)*kk%XALBF(:) )/(1.-pek%XPSN(:))
     pemist(:) = ((1.-kk%XFF(:)-pek%XPSN(:))*pek%XEMIS(:)         + kk%XFF(:)*kk%XEMISF(:))/(1.-pek%XPSN(:))
  ELSEWHERE
     palbt (:) = pek%XSNOWFREE_ALB(:)
     pemist(:) = pek%XEMIS(:)
  ENDWHERE
!
  ENDIF
!
!
ELSE
!
! Grid averaged albedo and emissivity for composite snow scheme:
!
   IF(pek%TSNOW%SCHEME=='EBA') THEN
!
      palbt(:)  = (1-pek%XVEG(:)) * (pek%XSNOWFREE_ALB_SOIL(:)*(1-pek%XPSNG(:)) + &
                                                 pek%TSNOW%ALB(:)*pek%XPSNG(:)) + &
                     pek%XVEG(:)  * (pek%XSNOWFREE_ALB_VEG (:)*(1-pek%XPSNV_A(:)) +      &
                                                pek%TSNOW%ALB(:)*pek%XPSNV_A(:))  
!
      pemist(:) = pek%XEMIS(:)-pek%XPSN(:)*(pek%XEMIS(:)-xemcrin)
!      
   ELSE
!
      palbt(:) = ( 1.-pek%XPSN(:)-kk%XFF(:))* pek%XSNOWFREE_ALB(:) + &
                       pek%XPSN(:)           * pek%TSNOW%ALB(:)     + kk%XFF(:)*kk%XALBF(:)     
!
      pemist(:) = ( 1.-pek%XPSN(:)-kk%XFF(:))* pek%XEMIS(:) + &
                       pek%XPSN(:)           * xemissn      + kk%XFF(:)*kk%XEMISF(:)
!
   ENDIF

ENDIF
!
!-------------------------------------------------------------------------------
!
!*       5. CALCULATION OF ZA, ZB, ZC
!       -----------------------------
!
!       5.1. Default
!       ------------
!
ztrad(:) =  pemist(:) * xstefan * (pek%XTG(:,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 + dmk%XCT(:) *                         &
         ((zfnsnow(:) *                                    &
           ( 4.*ztrad(:) + zrora(:)*zcps(:)*zpeta2(:) ))   &
         + zchums(:)*pdqsat(:) - zchuma(:)*zpeq_a_coef(:)) &
         + 2. * xpi / xday        
!
zb(:) = 1. / ptstep + dmk%XCT(:) * ( zfnsnow(:)* 3.*ztrad(:) + zchums(:)*pdqsat(:) ) 
!
zc(:) = 2. * xpi * pek%XTG(:,2) / xday + dmk%XCT(:) *       &
       ( zfnsnow(:) *                                       &
       ( zrora(:)*zcps(:)*zpetb2(:)                         &
       + psw_rad(:)*(1.-palbt(:)) + plw_rad(:)*pemist(:))   &
       - (zchums(:)*pqsat(:) - zchuma(:)*zpeq_b_coef(:)))           
!
IF(pek%TSNOW%SCHEME == '3-L' .OR. pek%TSNOW%SCHEME == 'CRO' .OR. io%CISBA == 'DIF')THEN                                 
!
!       5.2. With CSNOW=SNOW3L or CSNOW=CRO or IO%CISBA=DIF
!       -------------------------------------------------
!
   zc(:) = zc(:) + dmk%XCT(:)*(pek%XPSN(:)*pgrndflux(:)+pflux_cor(:,1))
!
ELSEIF (lcpl_arp) THEN
!
!       5.3. With Arpege
!       ----------------
!
zcdqsat(:) = (xcpv-xcpd)*zhums(:)*pdqsat(:)
zincr(:)= dmk%XCT(:) * zrora(:) * &
          (zcdqsat(:) * ( zpeta2(:)*pek%XTG(:,1) - zpetb2(:)) + &
          zxcpv_xcl_avg(:) * &
          (zhums(:)*pqsat(:) - zhuma(:) * (zpeq_b_coef(:) + zpeq_a_coef(:) * pek%XTG(:,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(:)*pek%XTG(:,1)
    zdkqb(:)=zpeq_b_coef(:)-pqa(:)

    za(:) = za(:) + dmk%XCT(:) * zrora(:) * &
             (2.* zpeqa2(:) + &
             zcnhuma(:) * (zdkqb(:)*zpeta2(:) - zpeq_a_coef(:)*zpetb2(:) )) 

    zb(:) = zb(:) + dmk%XCT(:) * zrora(:) * zpeqa2(:)         
  
    zc(:) = zc(:) + dmk%XCT(:)*zrora(:)*zcnhuma(:) *zdkqb(:)*zpetb2(:) 
           
  ENDIF
ENDIF
!

!-------------------------------------------------------------------------------
!
!*       6.     T AT TIME 'T+DT' (before snowmelt or soil ice evolution)
!               -----------------
!
IF(io%CISBA == '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(:) = pk%XDZG(:,1)/((pk%XDZG(:,1)+pk%XDZG(:,2))*psoilcondz(:,1))
   zcond2(:) = pk%XDZG(:,2)/((pk%XDZG(:,1)+pk%XDZG(:,2))*psoilcondz(:,2))
!
   zcondavg(:) = 1.0/(zcond1(:)+zcond2(:))
!   
   za(:)       = za(:) - (2. * xpi / xday) + zcondavg(:)*dmk%XCG(:)/pk%XDZDIF(:,1)
   zterm2(:)   = zcondavg(:)*dmk%XCG(:)/(za(:)*pk%XDZDIF(:,1))
   zterm1(:)   = (pek%XTG(:,1)*zb(:) + (zc(:) - (2. * xpi * pek%XTG(:,2) / xday)) )/za(:)  
!
! Determine the soil temperatures:
!
   CALL soil_heatdif(ptstep,pk%XDZG(:,:),pk%XDZDIF(:,:),psoilcondz,  &
                     psoilhcapz,dmk%XCG,zterm1,zterm2, ptdeep_a,     &
                     kk%XTDEEP, pek%XTG(:,:), pdeep_flux, pflux_cor   )
!
!
! "Restore" flux here is actually the heat flux between the surface
! and sub-surface layers (W m-2):
!
   prestore(:) = zcondavg(:)*(pek%XTG(:,1)-pek%XTG(:,2))/pk%XDZDIF(:,1)
!
ELSE
!
  IF(io%LTEMP_ARP)THEN
!
    CALL soil_temp_arp(ptstep,za,zb,zc,kk%XGAMMAT,kk%XTDEEP,io%XSODELX,pek%XTG(:,:))
!
!     "Restore" flux between surface and deep layer(W m-2):
    prestore(:) = 2.0 * xpi * (pek%XTG(:,1)-pek%XTG(:,2)) / &
                         ( dmk%XCT(:)*xday*io%XSODELX(1)*(io%XSODELX(1)+io%XSODELX(2)) )
!      
  ELSE
!
    pek%XTG(:,1) = ( pek%XTG(:,1)*zb(:) + zc(:) ) / za(:)
!
    WHERE(kk%XTDEEP(:) /= xundef .AND. kk%XGAMMAT(:) /= xundef)
      pek%XTG(:,2) = (pek%XTG(:,2) + (ptstep/xday)*(pek%XTG(:,1) + kk%XGAMMAT(:)*kk%XTDEEP(:))) / &
                         (1.+(ptstep/xday)*(1.0+kk%XGAMMAT(:)))  
    ELSEWHERE
      pek%XTG(:,2) = (pek%XTG(:,2) + (ptstep/xday)*pek%XTG(:,1)) /        &
                         (1.+(ptstep/xday) )  
    END WHERE
!
!     "Restore" flux between surface and deep layer(W m-2):
    prestore(:) = 2.0*xpi*(pek%XTG(:,1)-pt2m(:))/(dmk%XCT(:)*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(:)*pek%XTG(:,1)   + zpeq_b_coef(:)
!
pta_ic(:) =  zpet_a_coef(:)*pek%XTG(:,1)   + zpet_b_coef(:)
!
pustar2_ic(:) =  zustar2(:)
!
!--------------------------------------------------------------------------------------
!*       8.     Update of LSTT and LVTT for Arpege
!               ----------------------------------
!
IF (lcpl_arp) THEN

  IF (.NOT.lqvnplus) THEN
    pk%XCPS(:) =  pk%XCPS(:) + (xcpv-xcpd) *zhums(:)*pdqsat(:)*(pek%XTG(:,1)-ptsm(:))
  ENDIF


  IF (lqvnplus) THEN
    pk%XCPS(:) =  pk%XCPS(:) + (xcpv-xcpd) *zhums(:)*pdqsat(:)*(pek%XTG(:,1)-ptsm(:))  &
                       + (xcpv-xcpd) *(1-zhuma(:))*(pqa_ic(:)-pqa(:))  
  ENDIF

  pk%XLSTT(:) = pk%XLSTT(:) + (xcpv-xci)*(pek%XTG(:,1)-ptsm(:))

  pk%XLVTT(:) = pk%XLVTT(:) + (xcpv-xcl)*(pek%XTG(:,1)-ptsm(:))


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