e_budget_meb.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_meb(IO, KK, PK, PEK, DK, DEK, DMK,     &
                              PTSTEP, PLTT, PPS, PCT, PTDEEP_A, PD_G, PSOILCONDZ,   &
                              PSOILHCAPZ, PSNOWRHO, PSNOWCONDZ, PSNOWHCAPZ, PTAU_N, &
                              PLWNET_V_DTV, PLWNET_V_DTG, PLWNET_V_DTN,    &
                              PLWNET_G_DTV, PLWNET_G_DTG, PLWNET_G_DTN,    &
                              PLWNET_N_DTV, PLWNET_N_DTG, PLWNET_N_DTN,    &
                              PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF,       &
                              PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF,       &
                              PTHRMA_TA, PTHRMB_TA, PTHRMA_TC, PTHRMB_TC,  &
                              PTHRMA_TG, PTHRMB_TG, PTHRMA_TV, PTHRMB_TV,  &
                              PTHRMA_TN, PTHRMB_TN, PQSAT_G, PQSAT_V,      &
                              PQSATI_N, PPSNA, PPSNCV, PCHEATV, PCHEATG,   &
                              PCHEATN, PLEG_DELTA, PLEGI_DELTA, PHUGI,     &
                              PHVG, PHVN, PFROZEN1, PFLXC_C_A, PFLXC_G_C,  &
                              PFLXC_VG_C, PFLXC_VN_C, PFLXC_N_C, PFLXC_N_A,&
                              PFLXC_MOM, PTG, PSNOWLIQ, PFLXC_V_C, PHVGS, PHVNS, & 
                              PDQSAT_G, PDQSAT_V, PDQSATI_N, PTA_IC,       &
                              PQA_IC, PUSTAR2_IC, PVMOD, PDELTAT_G,        &
                              PDELTAT_V, PDELTAT_N, PGRNDFLUX, PDEEP_FLUX, &
                              PDELHEATV_SFC, PDELHEATG_SFC, PDELHEATG     )
!     ##########################################################################
!
!!****  *E_BUDGET*  
!!
!!    PURPOSE
!!    -------
!
!     Calculates the implicit simulataneous solution of multiple surface
!     energy budgets (understory vegetation-soil composite, vegetation canopy
!     and explicit snowpack) and sub-surface soil temperatures. Canopy
!     air space temperature and specific humidities are also diagnosed.
!     Fully implicit coupling with the atmosphere allowed. 
!     Note that a 'test' snow temperature is computed herein which is just
!     use to compute the surface fluxes. These fluxes are then used within
!     the snow scheme as an upper boundary condition to compute the snow 
!     temperatures. In theory, they should be very close, but this method is 
!     done to ensure a high level of energy conservation.
!         
!     
!!**  METHOD
!!    ------
!
!     1- compute coefficients for each energy budget
!     2- solve the equations for snow, vegetation canopy and soil-understory
!        vegetation composite
!     3- solve soil T profile
!     4- diagnose canopy air space variables, and lowest level atmospheric
!        state variable values at t+dt
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!    Noilhan and Planton (1989)
!!    Belair (1995)
!!    * to be done * (2011)
!!      
!!    AUTHOR
!!    ------
!!
!!      A. Boone           * CNRM-GAME, Meteo-France *
!!      P. Samuelsson      * SMHI *
!!      S. Gollvik         * SMHI * 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    22/01/11 
!!                  10/10/14 (A. Boone) Removed understory vegetation
!!
!-------------------------------------------------------------------------------
!
!*       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_evap_isba_n, ONLY : diag_evap_isba_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t
!
USE modd_csts,                    ONLY : xlvtt, xlstt, xtt, xcpd, xcpv, xcl, &
                                         xday, xpi, xlmtt, xrholw
USE modd_surf_atm,                ONLY : lcpl_arp
USE modd_surf_par,                ONLY : xundef
USE modd_snow_metamo,             ONLY : xsnowdzmin
!
USE mode_thermos
USE mode_meb,                     ONLY : sfc_heatcap_veg
USE mode_snow3l,                  ONLY : snow3lhold
!
USE modi_tridiag_ground_rm_coefs
USE modi_tridiag_ground_rm_soln
!
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_evap_isba_t), INTENT(INOUT) :: DEK
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DMK
!
REAL,                 INTENT(IN)   :: PTSTEP
!                                     PTSTEP = timestep of the integration (s)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PCT
!                                     PCT     = surface understory or composite thermal inertia (m2 K J-1)
!
REAL, DIMENSION(:), INTENT(IN)     :: PTDEEP_A
!                                      PTDEEP_A = Deep soil temperature
!                                                 coefficient depending on flux
!                                      Tdeep = IP%XTDEEP + PTDEEP_A * PDEEP_FLUX
!                                              (with PDEEP_FLUX in W/m2)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PPS
!                                     PPS  = surface pressure (Pa)
!
REAL, DIMENSION(:,:), INTENT(IN)   :: PD_G, PSOILCONDZ, PSOILHCAPZ
!                                     PD_G       = soil layer depth      (m)
!                                     PSOILCONDZ = soil thermal conductivity (W m-1 K-1)
!                                     PSOILHCAPZ = soil heat capacity        (J m-3 K-1)
!
REAL, DIMENSION(:,:), INTENT(IN)   :: PSNOWCONDZ, PSNOWHCAPZ, PSNOWRHO
!                                     PSNOWCONDZ = snow thermal conductivity (W m-1 K-1)
!                                     PSNOWHCAPZ = snow heat capacity        (J m-3 K-1)
!                                     PSNOWRHO   = snow layer density        (kg/m3)
!
REAL, DIMENSION(:,:), INTENT(IN)   :: PTAU_N
!                                     PTAU_N   = shortwave radiation transmission through (at the base of)
!                                                each snow layer. Implicitly PTAU_N(0) = 1 (-)
!                                                It decreases with depth to zero somewhere in a deep snowpack,
!                                                but can be above 0 at snowpack base for shallow snow, then
!                                                it is presumed to go into the uppermost soil layer.
!
REAL, DIMENSION(:),   INTENT(IN)   :: PLWNET_V_DTV, PLWNET_V_DTG, PLWNET_V_DTN
!                                     PLWNET_V_DTV, PLWNET_V_DTG, PLWNET_V_DTN = Vegetation canopy net LW radiation 
!                                     derivatives w/r/t surface temperature(s) (W m-2 K-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PLWNET_G_DTV, PLWNET_G_DTG, PLWNET_G_DTN
!                                     PLWNET_G_DTV, PLWNET_G_DTG, PLWNET_G_DTN = Understory-ground net LW radiation 
!                                          derivatives w/r/t surface temperature(s) (W m-2 K-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PLWNET_N_DTV, PLWNET_N_DTG, PLWNET_N_DTN
!                                     PLWNET_N_DTV, PLWNET_N_DTG, PLWNET_N_DTN = Ground-based snow net LW radiation 
!                                          derivatives w/r/t surface temperature(s) (W m-2 K-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PTHRMA_TA, PTHRMB_TA, PTHRMA_TC, PTHRMB_TC,                     &
                                      PTHRMA_TG, PTHRMB_TG, PTHRMA_TV, PTHRMB_TV, PTHRMA_TN, PTHRMB_TN
!                                     PTHRMA_TA                                                    (J kg-1 K-1)
!                                     PTHRMB_TA = linear transform coefficinets for atmospheric
!                                                 thermal variable for lowest atmospheric level.   (J kg-1)
!                                                 Transform T to dry static energy or enthalpy.
!                                     PTHRMA_TC                                                    (J kg-1 K-1)
!                                     PTHRMB_TC = linear transform coefficinets for atmospheric
!                                                 thermal variable for canopy air                  (J kg-1)
!                                                 Transform T to dry static energy or enthalpy.
!                                     PTHRMA_TG,V,N                                                (J kg-1 K-1)
!                                     PTHRMB_TG,V,N = linear transform coefficinets for atmospheric
!                                                 thermal variable for surfaces (G, V, and N)      (J kg-1)
!                                                 Transform T to dry static energy or enthalpy.
!
REAL, DIMENSION(:),   INTENT(IN)   :: PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF,             &
                                      PPEW_A_COEF, PPEW_B_COEF
!                                     PPEW_A_COEF = wind atmospheric coupling coefficient 
!                                     PPEW_B_COEF = wind atmospheric coupling coefficient 
!                                     PPET_A_COEF = A-air temperature atmospheric coupling coefficient 
!                                     PPET_B_COEF = B-air temperature atmospheric coupling coefficient
!                                     PPEQ_A_COEF = A-air specific humidity atmospheric coupling coefficient
!                                     PPEQ_B_COEF = B-air specific humidity atmospheric coupling coefficient
!
REAL, DIMENSION(:),   INTENT(IN)   :: PQSAT_G, PQSAT_V, PQSATI_N
!                                     PQSAT_G  = saturation specific humidity for understory surface    (kg kg-1)
!                                     PQSAT_V  = saturation specific humidity for the vegetation canopy (kg kg-1)
!                                     PQSATI_N = saturation specific humidity over ice for the snowpack (kg kg-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PPSNA, PPSNCV
!                                     PPSN     = fraction of snow on ground and understory vegetation          (-)
!                                     PPSNA    = fraction of vegetation canopy buried by ground-based snowpack (-)
!                                     PPSNCV   = fraction of vegetation canopy covered by intercepted snow     (-)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PLEG_DELTA, PLEGI_DELTA, PHUGI, PHVG, PHVN, PFROZEN1
!                                     PHUGI       = relative humidity of the frozen surface soil                      (-)                         
!                                     PHVG        = Halstead coefficient of non-buried (snow) canopy vegetation       (-)                         
!                                     PHVN        = Halstead coefficient of paritally-buried (snow) canopy vegetation (-)  
!                                     PLEG_DELTA  = soil evaporation delta fn                                         (-)
!                                     PLEGI_DELTA = soil sublimation delta fn                                         (-)
!                                     PFROZEN1    = fraction of surface soil layer which is frozen                    (-)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PFLXC_C_A, PFLXC_G_C, PFLXC_VG_C, PFLXC_VN_C, PFLXC_N_C, PFLXC_N_A, PFLXC_MOM
!                                     PFLXC_C_A  = Flux form heat transfer coefficient: canopy air to atmosphere (kg m-2 s-1)
!                                     PFLXC_G_C  = As above, but for : ground-understory to canopy air           (kg m-2 s-1)
!                                     PFLXC_VG_C = As above, but for : non-snow buried canopy to canopy air      (kg m-2 s-1)
!                                     PFLXC_VN_C = As above, but for : partially snow-buried canopy air to canopy 
!                                                  air                                                           (kg m-2 s-1)
!                                     PFLXC_N_C  = As above, but for : ground-based snow to atmosphere           (kg m-2 s-1)
!                                     PFLXC_N_A  = As above, but for : ground-based snow to canopy air           (kg m-2 s-1)
!                                     PFLXC_MOM  = flux form drag transfer coefficient: canopy air to atmosphere (kg m-2 s-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: PLTT
!                                     PLTT    = latent heat normalization factor (J kg-1)
!
REAL, DIMENSION(:,:), INTENT(INOUT):: PTG
!                                     PTG    = Soil temperature profile (K)
!
REAL, DIMENSION(:,:), INTENT(INOUT):: PSNOWLIQ
!                                     PSNOWLIQ = Snow layer liquid water content       (m)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PCHEATV, PCHEATG, PCHEATN
!                                     PCHEATV = Vegetation canopy *effective surface* heat capacity          (J m-2 K-1)
!                                     PCHEATG = Understory-ground *effective surface* heat capacity (J m-2 K-1)
!                                     PCHEATN = Ground-based snow *effective surface* heat capacity          (J m-2 K-1)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PDQSAT_G, PDQSAT_V, PDQSATI_N
!                                     PQSAT_G  = saturation specific humidity derivative for understory 
!                                                surface               (kg kg-1 K-1)
!                                     PQSAT_V  = saturation specific humidity derivative for the vegetation 
!                                                canopy                (kg kg-1 K-1)
!                                     PQSATI_N = saturation specific humidity derivative over ice for the 
!                                                ground-based snowpack (kg kg-1 K-1)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PFLXC_V_C, PHVGS, PHVNS
!                                     PFLXC_V_C = Flux form heat transfer coefficient: total canopy (partially 
!                                                 snow-buried and non-buried parts) to canopy air               (kg m-2 s-1)
!                                     PHVGS     = Dimensionless pseudo humidity factor for computing vapor
!                                                 fluxes from the non-buried part of the canopy to the canopy air    (-)
!                                     PHVNS     = Dimensionless pseudo humidity factor for computing vapor
!                                                 fluxes from the partly-buried part of the canopy to the canopy air (-)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PTA_IC, PQA_IC, PUSTAR2_IC, PVMOD
!                                     PTA_IC      = Near-ground air temperature        (K)
!                                     PQA_IC      = Near-ground air specific humidity  (kg kg-1)
!                                     PUSTAR2_IC  = Surface friction velocity squared  (m2 s-2)
!                                                   (modified if implicit coupling with
!                                                   atmosphere used)
!                                     PVMOD       = lowest level wind speed            (m s-1)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PDELTAT_V, PDELTAT_N, PDELTAT_G
!                                     PDELTAT_V = Time change in vegetation canopy temperature (K)
!                                     PDELTAT_N = Time change in snowpack surface temperature  (K)
!                                     PDELTAT_G = Time change in soil surface temperature      (K)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: PGRNDFLUX 
!                                     PGRNDFLUX = Flux between snowpack base and ground surface (W m-2)
!
REAL, DIMENSION(:),  INTENT(OUT)   :: PDEEP_FLUX ! Heat flux at bottom of ISBA (W/m2)
!
REAL, DIMENSION(:),  INTENT(OUT)   :: PDELHEATV_SFC, PDELHEATG_SFC, PDELHEATG
!                                     PDELHEATV_SFC = change in heat storage of the vegetation canopy layer over the current time step (W m-2)
!                                     PDELHEATG_SFC = change in heat storage of the surface soil layer over the current time step (W m-2)
!                                     PDELHEATG     = change in heat storage of the entire soil column over the current time step (W m-2)
!
!*      0.2    declarations of local variables
!
!
INTEGER                                   :: JNSNOW, JNGRND, JNPTS, JJ, JK, JL
!
REAL, DIMENSION(SIZE(PPS))                :: ZHN, ZHS, ZHVS, ZHNS
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2))  :: ZTGO
!
REAL, DIMENSION(SIZE(DMK%XSNOWTEMP,1),SIZE(DMK%XSNOWTEMP,2))  :: ZTNO
!
REAL, DIMENSION(SIZE(PTG,1))              :: ZTVO
!
REAL, DIMENSION(SIZE(PPS))                :: ZPSNAG, ZWORK, ZFFF, ZGCOND1
!
REAL, DIMENSION(SIZE(PPS))                :: ZPET_A_COEF_P, ZPET_B_COEF_P, ZPET_C_COEF_P
!
REAL, DIMENSION(SIZE(PPS))                :: ZPEQ_A_COEF_P, ZPEQ_B_COEF_P, ZPEQ_C_COEF_P
!
REAL, DIMENSION(SIZE(PPS))                :: ZCOEFA_TC, ZCOEFB_TC, ZCOEFC_TC, ZCOEFD_TC,          &
                                             ZCOEFA_QC, ZCOEFB_QC, ZCOEFC_QC, ZCOEFD_QC
!
REAL, DIMENSION(SIZE(PPS))                :: ZBETA_V, ZALPHA_V, ZGAMMA_V,                         &
                                             ZBETA_G, ZALPHA_G, ZGAMMA_G,                         &
                                             ZBETA_N, ZALPHA_N, ZGAMMA_N,                         &
                                             ZBETA_P_V, ZALPHA_P_V, ZBETA_P_N, ZALPHA_P_N
!
REAL, DIMENSION(SIZE(PPS))                :: ZRNET_NN, ZRNET_N_DTNN, ZRNET_N_DTGN, ZRNET_N_DTVN
!
REAL, DIMENSION(SIZE(PPS))                :: ZVMOD, ZUSTAR2, ZPSNA
!
REAL, DIMENSION(SIZE(PPS))                :: ZTCONDA_DELZ_G, ZTCONDA_DELZ_N, ZTCONDA_DELZ_NG
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2))  :: ZSOIL_COEF_A, ZSOIL_COEF_B
!
REAL, DIMENSION(SIZE(DMK%XSNOWDZ,1),SIZE(DMK%XSNOWDZ,2)):: ZSNOW_COEF_A, ZSNOW_COEF_B
!
REAL, DIMENSION(SIZE(DMK%XSNOWDZ,1),SIZE(DMK%XSNOWDZ,2)):: ZWHOLDMAX
!
REAL, DIMENSION(SIZE(PD_G,1),SIZE(PD_G,2)+SIZE(DMK%XSNOWDZ,2)) :: ZD, ZT, ZHCAPZ, ZCONDZ,         &
                                             ZCOEF_A, ZCOEF_B, ZSOURCE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
REAL, PARAMETER                          :: ZERTOL       = 1.0e-6  ! -
REAL, PARAMETER                          :: ZERTOL_FLX_C = 1.0e-12 ! -
!
!-------------------------------------------------------------------------------
!
!*       0.     Initialization:
!               ---------------
!
IF (lhook) CALL dr_hook('E_BUDGET_MEB',0,zhook_handle)
!
ztgo(:,:) = ptg(:,:)
ztno(:,:) = dmk%XSNOWTEMP(:,:)
ztvo(:)   = pek%XTV(:)

! To prevent possible numerical problems in the limit as psna==>1, we
! limit the value for certain computations:

zpsna(:)  = min(1.0-zertol, ppsna(:))

jnsnow       = SIZE(dmk%XSNOWTEMP,2)
jngrnd       = SIZE(ptg,2)
jnpts        = SIZE(ptg,1)

! sub-surface / surface coupling coefficients:

zsoil_coef_a(:,:) = 0.0 
zsoil_coef_b(:,:) = 0.0
zsnow_coef_a(:,:) = 0.0
zsnow_coef_b(:,:) = 0.0

!-------------------------------------------------------------------------------
!
!*       1.     Some variables/coefficients needed for coupling or solution
!               -----------------------------------------------------------
! - Effective Surface heat capacities for each energy budget (J m-2 K-1)
!
pcheatg(:)    = 1/pct(:)                                               ! understory soil/floodplain
!
pcheatv(:)    = sfc_heatcap_veg(pek%XWRVN(:),pek%XWR(:),pek%XCV)         ! vegetation canopy heat capacity
!
pcheatn(:)    = psnowhcapz(:,1)*dmk%XSNOWDZ(:,1)                      ! snow surface layer
!
! - Specific humidity derivatives (kg kg-1 K-1)
!
pdqsat_g(:)   = dqsat(ztgo(:,1),  pps(:),pqsat_g(:)  )
pdqsat_v(:)   = dqsat(ztvo(:),    pps(:),pqsat_v(:)  )
pdqsati_n(:)  = dqsati(ztno(:,1), pps(:),pqsati_n(:) )
!
! - Precipitation heating term for snowpack:
!   Rainfall renders it's heat to the snow when it enters
!   the snowpack:
! NOTE: for now set to zero because we'd need to remove this heat from
!       the atmosphere to conserve energy. This can be done, but should
!       be tested within an atmos model first probably...
!       ALSO, water should include canopy drip etc....but again, conservation
!       would need to be carefully accounted for. So for now within MEB,
!       set to 0. But code within MEB and ES can accept this term, so
!       if this process is added at some point, it would be here. 
!       Also, it is treated explicitly...this might also have to be changed
!       to implicit depending how this term is expressed. That would require
!       modification to this routine and flux routine:  but retain it here for possible
!       future implementation.
!
dmk%XHPSNOW(:)   = 0.0   ! W m-2
!
! Snowmeltwater advection which can heat uppermost soil layer if below freezing.
! NOTE this term is currently OFF for MEB since it would mean keeping meltwater
! flux from previous timestep (since snow called *after* energy budget). Also,
! the heat would need to be removed from snowpack to conserve energy (and currently
! not done) so we turn off for now (set to 0), but retain it here for possible
! future implementation.
!
dek%XMELTADV(:)  = 0.0   ! W m-2
!
! - Implicit wind speed at lowest atmospheric level for this patch:
!
zustar2(:)   =    (pflxc_mom(:)*ppew_b_coef(:))/        &
              (1.0-pflxc_mom(:)*ppew_a_coef(:))
!
zvmod(:)     = ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
!
pvmod(:)     = max(zvmod,0.)
!
WHERE (ppew_a_coef(:) /= 0.) 
     zustar2(:) = max(0., ( pvmod(:) - ppew_b_coef(:) ) / ppew_a_coef(:) )
END WHERE
!
pustar2_ic(:)= zustar2(:) 
!
! NOTE: put in new option HIMPLICIT_WIND=='OLD' or 'NEW' ?
!
zsource(:,:) = 0. ! heat Eq source term initialization: c dz dT/dt = d/dz(k dT/dz) + Source
!
IF(io%CISBA == 'DIF')THEN
!
!*       2.a    Compute sub-surface soil coupling coefficients: upward sweep
!               ------------------------------------------------------------
! Upward sweep of the tridiagnal matrix using R&M method. Also compute
! interfacial thermal conductivity to layer thickness ratio (W m-2 K-1) 
! These coefficients are used to compute the temperature profile implicitly.
!
   CALL tridiag_ground_rm_coefs(ptstep, pd_g, ztgo, psoilhcapz, psoilcondz,   &
                                zsource(:,1:jngrnd), ptdeep_a, kk%XTDEEP,     &
                                ztconda_delz_g, zsoil_coef_a, zsoil_coef_b)
!
! Here we repeat this for the snowpack: but note, this is just
! to better estimate surface fluxes via the surface-sub surface heat
! Flux term (we do not actually solve
! the snow T profile here since fractional coverage not included). 
! We start from base of soil up through snow
! to surface. This gives coefficients which are fully implicit
! from the base of the soil to the snow surface, so numerically it is
! quite robust. Note, this only corresponds to the snow-covered part of grid box,
! so it is more accurate as the snow fraction approaches unity.
! Starting from snowpack surface downward to base of ground:
!
   jl                     = 1
   zd(:,jl)               = dmk%XSNOWDZ(:,1)
   zt(:,jl)               = ztno(:,1)
   zhcapz(:,jl)           = psnowhcapz(:,1)
   zcondz(:,jl)           = psnowcondz(:,1)
   zsource(:,jl)          = 0.               ! Included in energy budget
   DO jk=2,jnsnow
      jl                  = jl + 1
      DO jj=1,jnpts
         zd(jj,jl)        = zd(jj,jl-1) + dmk%XSNOWDZ(jj,jk)
         zt(jj,jl)        = ztno(jj,jk)
         zhcapz(jj,jl)    = psnowhcapz(jj,jk)
         zcondz(jj,jl)    = psnowcondz(jj,jk)
         zsource(jj,jl)   = dek%XSWNET_NS(jj)*(ptau_n(jj,jk-1)-ptau_n(jj,jk))
      ENDDO
   ENDDO
   jl                     = jl + 1
   zd(:,jl)               = pd_g(:,1)
   zt(:,jl)               = ztgo(:,1)
   zhcapz(:,jl)           = pcheatg(:)/pd_g(:,1)
   zcondz(:,jl)           = psoilcondz(:,1)
   zsource(:,jl)          = dek%XSWNET_NS(:)*ptau_n(:,jnsnow)
   DO jk=2,jngrnd
      jl                  = jl + 1
      DO jj=1,jnpts
         zd(jj,jl)        = pd_g(jj,jk)
         zt(jj,jl)        = ztgo(jj,jk)
         zhcapz(jj,jl)    = psoilhcapz(jj,jk)
         zcondz(jj,jl)    = psoilcondz(jj,jk)
         zsource(jj,jl)   = 0.
      ENDDO
   ENDDO
!
! Get coefficients from upward sweep (starting from soil base to snow surface):
!
   CALL tridiag_ground_rm_coefs(ptstep, zd, zt, zhcapz, zcondz, zsource, &
                                ptdeep_a, kk%XTDEEP, ztconda_delz_n, zcoef_a, zcoef_b)
!
   zsnow_coef_a(:,2)  = zcoef_a(:,2)
   zsnow_coef_b(:,2)  = zcoef_b(:,2)
!
   zgcond1(:)         = psoilcondz(:,1) ! save uppermost soil thermal condcutivity
!
ELSE
!
   zsoil_coef_a(:,2) = (ptstep/xday)/(1.0 + ptstep/xday)
   zsoil_coef_b(:,2) = ztgo(:,2)    /(1.0 + ptstep/xday)

   ztconda_delz_g(:) = (2*xpi/xday)/pct(:)

! - Soil thermal conductivity
!   is implicit in Force-Restore soil method, so it
!   must be backed-out of surface thermal coefficients
!   (Etchevers and Martin 1997):

   zgcond1(:)            = (4*xpi/xday)/( dmk%XCG(:)*dmk%XCG(:)/(pd_g(:,1)*pct(:)) )


   jl                    = 1
   zd(:,jl)              = dmk%XSNOWDZ(:,1)
   zt(:,jl)              = ztno(:,1)
   zhcapz(:,jl)          = psnowhcapz(:,1)
   zcondz(:,jl)          = psnowcondz(:,1)
   zsource(:,jl)         = 0.               ! Included in energy budget
   DO jk=2,jnsnow
      jl                 = jl + 1
      DO jj=1,jnpts
         zd(jj,jl)       = zd(jj,jl-1) + dmk%XSNOWDZ(jj,jk)
         zt(jj,jl)       = ztno(jj,jk)
         zhcapz(jj,jl)   = psnowhcapz(jj,jk)
         zcondz(jj,jl)   = psnowcondz(jj,jk)
         zsource(jj,jl)  = dek%XSWNET_NS(jj)*(ptau_n(jj,jk-1)-ptau_n(jj,jk))
      ENDDO
   ENDDO
   jl                    = jl + 1
   zd(:,jl)              = pd_g(:,1)
   zt(:,jl)              = ztgo(:,1)
   zhcapz(:,jl)          = 1/pct(:)
   zcondz(:,jl)          = zgcond1(:)
   zsource(:,jl)         = dek%XSWNET_NS(:)*ptau_n(:,jnsnow)

! Get coefficients from upward sweep (starting from soil base to snow surface):
!
   CALL tridiag_ground_rm_coefs(ptstep, zd(:,1:jl), zt(:,1:jl), zhcapz(:,1:jl), zcondz(:,1:jl), &
                                zsource(:,1:jl), ptdeep_a, kk%XTDEEP, ztconda_delz_n, &
                                zcoef_a(:,1:jl), zcoef_b(:,1:jl))

   zsnow_coef_a(:,2)     = zcoef_a(:,2)
   zsnow_coef_b(:,2)     = zcoef_b(:,2)
!
ENDIF
!
! Interfacial ground-snowbase thermal conductivity divided by interfacial dz:
!
ztconda_delz_ng(:) = 2/((dmk%XSNOWDZ(:,jnsnow)/psnowcondz(:,jnsnow))+(pd_g(:,1)/zgcond1(:) ))
!
!
!*       3.     Pseudo humidity factors (-) and effective latent heat (J kg-1)
!               ---------------------------------------------------------
!
! Compute the average flux heat exchange coefficient for the canopy: (kg m-2 s-1)
! Numerical: let get very small (fluxes can become essentially negligible), but not zero
!
pflxc_v_c(:)   = pflxc_vg_c(:)*(1.-pek%XPSN(:)) + pflxc_vn_c(:)*pek%XPSN(:)*(1.-zpsna(:))
pflxc_v_c(:)   = max(pflxc_v_c(:), zertol_flx_c)

! Understory vegetation and ground factors:

zfff(:)        = kk%XFF(:)*( (1.0 - kk%XFFROZEN(:))*(xlvtt/pltt(:)) +      &
                                    kk%XFFROZEN(:) *(xlstt/pltt(:)) )
 
zhn(:)         = (1.0-pek%XPSN(:)-kk%XFF(:))*(                                        &
                            pleg_delta(:) *(1.0-pfrozen1(:))*(xlvtt/pltt(:))          &
                 +          plegi_delta(:)*      pfrozen1(:) *(xlstt/pltt(:)) ) + zfff(:)         

zhs(:)         = (1.0-pek%XPSN(:)-kk%XFF(:))*(                                           &
                   dk%XHUG(:)  *pleg_delta(:) *(1.0-pfrozen1(:)) *(xlvtt/pltt(:))       &
                     + phugi(:)*plegi_delta(:)*      pfrozen1(:)  *(xlstt/pltt(:)) ) + zfff(:)

! adjust for local use in solution (since they are multiplied by 1-psn herein):

zhn(:)         = zhn(:)/max(1.0 - pek%XPSN(:) , zertol)

zhs(:)         = zhs(:)/max(1.0 - pek%XPSN(:) , zertol)

! Vegetation canopy factors:

phvgs(:)       = (1.-zpsna(:))*pek%XPSN(:) *phvn(:)*(pflxc_vn_c(:)/pflxc_v_c(:)) + &
                           (1.-pek%XPSN(:))*phvg(:)*(pflxc_vg_c(:)/pflxc_v_c(:))

phvns(:)       = (1.-zpsna(:))*pek%XPSN(:) *        (pflxc_vn_c(:)/pflxc_v_c(:)) + &
                           (1.-pek%XPSN(:))*        (pflxc_vg_c(:)/pflxc_v_c(:))

! - total canopy H factor (including intercepted snow)

zhvs(:)        = (1.-ppsncv(:))*(xlvtt/pltt(:))*phvgs(:) +  &
                     ppsncv(:) *(xlstt/pltt(:))*phvns(:)   

! Snow latent heating factor:
! NOTE, for now we consider only snow sublimation 
!       (not evaporation from snow liquid)
!
zhns           = (xlstt/pltt(:))
!
!
!*       4.     Transform atmospheric coupling coefficients for T and q
!               ---------------------------------------------------------
!
!  Transform coupling coefficients from flux form to scalar form for local
!  use in the energy budgets.
!
! - first, to shorten computations a bit, define a snow fraction product

zpsnag(:)         = 1.0 - pek%XPSN(:)*ppsna(:)

! T coefficients: where : TA = PET_B_COEF_P + PET_A_COEF_P*TC + PET_C_COEF_P*TN

zwork(:)          = pthrma_ta(:)*( 1.0 + ppet_a_coef(:)*(                                        &
                    pflxc_c_a(:)*zpsnag(:) + pflxc_n_a(:)*pek%XPSN(:)*ppsna(:)) )
zpet_a_coef_p(:)  =   ppet_a_coef(:)*pflxc_c_a(:)*zpsnag(:)*pthrma_tc(:)/zwork(:)
zpet_b_coef_p(:)  = ( ppet_b_coef(:) - pthrmb_ta(:) +                                            &
                      ppet_a_coef(:)*(pflxc_c_a(:)*zpsnag(:)*(pthrmb_tc(:)-pthrmb_ta(:)) +       &
                           pflxc_n_a(:)*pek%XPSN(:)*ppsna(:)*(pthrmb_tn(:)-pthrmb_ta(:)) ) ) /zwork(:)
zpet_c_coef_p(:)  =   ppet_a_coef(:)*pflxc_n_a(:)*pek%XPSN(:)*ppsna(:)*pthrma_tn(:)       /zwork(:)

! q coefficients:

zwork(:)          = 1.0 + ppeq_a_coef(:)*(pflxc_c_a(:)*zpsnag(:) +                               &
                                          pflxc_n_a(:)*pek%XPSN(:)*ppsna(:)*zhns)
zpeq_a_coef_p(:)  = ppeq_a_coef(:)*pflxc_c_a(:)*zpsnag(:)                /zwork(:)
zpeq_b_coef_p(:)  = ppeq_b_coef(:)                                       /zwork(:)
zpeq_c_coef_p(:)  = ppeq_a_coef(:)*pflxc_n_a(:)*pek%XPSN(:)*ppsna(:)*zhns/zwork(:)


!*       5.     Compute canopy air coefficients (T and q)
!               ---------------------------------------------------------


! - Canopy air T coefs, where : TC = COEFA_TC + COEFB_TC*TV + COEFC_TC*TG + COEFD_TC*TN

zwork(:)     = pflxc_c_a(:) *(pthrma_tc(:)-pthrma_ta(:)*zpet_a_coef_p(:))*zpsnag(:)                 &
                   + pflxc_v_c(:) *pthrma_tc(:)                                                     &
                   + pflxc_g_c(:) *pthrma_tc(:)*(1.0-pek%XPSN(:))                                       &
                   + pflxc_n_c(:) *pthrma_tc(:)*     pek%XPSN(:) *(1.0-ppsna(:))

zcoefa_tc(:) = (pflxc_c_a(:) * zpsnag(:) *(pthrma_ta(:)*zpet_b_coef_p(:)-pthrmb_tc(:)+pthrmb_ta(:)) &    
                   + pflxc_v_c(:) * (pthrmb_tv(:)-pthrmb_tc(:))                                     &
                   + pflxc_g_c(:) * (pthrmb_tg(:)-pthrmb_tc(:))*(1.0-pek%XPSN(:))                       &
                   + pflxc_n_c(:) * (pthrmb_tn(:)-pthrmb_tc(:))*     pek%XPSN(:) *(1.0-ppsna(:))        &
                                                                      )/zwork(:)

zcoefb_tc(:) = pflxc_v_c(:)*pthrma_tv(:)/zwork(:)

zcoefc_tc(:) = pflxc_g_c(:)*pthrma_tg(:)*(1.0-pek%XPSN(:))/zwork(:)

zcoefd_tc(:) =(pflxc_n_c(:)*pthrma_tn(:)*pek%XPSN(:)*(1.0-ppsna(:)) +                         &
               pflxc_c_a(:) *pthrma_ta(:)*zpet_c_coef_p(:)*zpsnag(:) ) /zwork(:)

!-  Canopy air q coefs, where : QC = COEFA_QC + COEFB_QC*TV + COEFC_QC*TG + COEFD_QC*TN

zwork(:)       = pflxc_c_a(:) *(1.-zpeq_a_coef_p(:))*zpsnag(:)                                      &
                   + pflxc_v_c(:)* zhvs(:)                                                          &
                   + pflxc_g_c(:) *zhn(:) *(1.0-pek%XPSN(:))                                       &
                   + pflxc_n_c(:) *zhns(:)*     pek%XPSN(:) *(1.0-ppsna(:))
zwork(:)       = max(zertol, zwork(:))

zcoefa_qc(:)   = ( pflxc_c_a(:) *zpeq_b_coef_p(:)*zpsnag(:)                                         &
                 + pflxc_v_c(:) *zhvs(:)*(pqsat_v(:)-pdqsat_v(:)*ztvo(:)  )                         &
                 + pflxc_g_c(:) *zhs(:) *(pqsat_g(:)-pdqsat_g(:)*ztgo(:,1))*(1.0-pek%XPSN(:))      &
                 + pflxc_n_c(:) *zhns(:)*(pqsati_n(:)-pdqsati_n(:)*ztno(:,1))*                      &
                                                                      pek%XPSN(:)*(1.0-ppsna(:))   &
                                                                                )/zwork(:) 

zcoefb_qc(:)   = pflxc_v_c(:) *zhvs(:)*pdqsat_v(:) /zwork(:)

zcoefc_qc(:)   = pflxc_g_c(:) *zhs(:) *pdqsat_g(:)*(1.0-pek%XPSN(:))/zwork(:)

zcoefd_qc(:)   = pflxc_n_c(:) *zhns(:)*pdqsati_n(:)*     pek%XPSN(:)*(1.0-ppsna(:))/zwork(:)

!*       6.     Surface Energy Budget(s) coefficients
!               -------------------------------------
! Each of the 'N' energy budgets is linearized, so we have 
! 'N' linear equations and 'N' unknowns. Here we set up the coefficients.
! For computations here, make snow radiative terms relative to snow surface:

zwork(:)        = 1/max(zertol,       pek%XPSN(:))

zrnet_nn(:)     = (dek%XSWNET_NS(:) + dek%XLWNET_N(:))*zwork(:)
zrnet_n_dtnn(:) = plwnet_n_dtn(:)            *zwork(:)
zrnet_n_dtgn(:) = plwnet_n_dtg(:)            *zwork(:)
zrnet_n_dtvn(:) = plwnet_n_dtv(:)            *zwork(:)

! Tv coefs, where TV = BETA_V + ALPHA_V*TG + GAMMA_V*TN

!
zwork(:)    =   (pcheatv(:)/ptstep) - plwnet_v_dtv(:)                                               &
              + pflxc_v_c(:) * ( pthrma_tv(:) - pthrma_tc(:)*zcoefb_tc(:)                           &
              + pltt(:)*zhvs(:)*(pdqsat_v(:) - zcoefb_qc(:)) )

zbeta_v(:)  = ( (pcheatv(:)/ptstep)*ztvo(:) + dek%XLWNET_V(:) + dek%XSWNET_V(:)                  &
              - plwnet_v_dtv(:)*ztvo(:) - plwnet_v_dtg(:)*ztgo(:,1) - plwnet_v_dtn(:)*ztno(:,1)  &
              - pflxc_v_c(:)*( pthrmb_tv(:)-pthrmb_tc(:)-pthrma_tc(:)*zcoefa_tc(:)               &
              + pltt(:)*zhvs(:)*(pqsat_v(:) - pdqsat_v(:)*ztvo(:)                                &
              - zcoefa_qc(:)) ) )/zwork(:)

zalpha_v(:) = (plwnet_v_dtg(:) + pflxc_v_c(:)*(pthrma_tc(:)*zcoefc_tc(:)                            &
              + pltt(:)*zhvs(:)*zcoefc_qc(:) ) )/zwork(:)

zgamma_v(:) = (plwnet_v_dtn(:) + pflxc_v_c(:)*(pthrma_tc(:)*zcoefd_tc(:)                            &
              + pltt(:)*zhvs(:)*zcoefd_qc(:) ) )/zwork(:)

! Tg coefs, where TG = BETA_G + ALPHA_G*TV + GAMMA_G*TN

zwork(:)    =   (pcheatg(:)/ptstep) - plwnet_g_dtg(:)                                               &
              + (1.0-pek%XPSN(:))*pflxc_g_c(:)*( (pthrma_tg(:) - pthrma_tc(:)*zcoefc_tc(:))        &
              + pltt(:)*(zhs(:)*pdqsat_g(:) - zhn(:)*zcoefc_qc(:)) )                                  &
              + ztconda_delz_g(:)*(1.0-zsoil_coef_a(:,2))                                           &
              + pek%XPSN(:)*ztconda_delz_ng(:)

zbeta_g(:)  = ( (pcheatg(:)/ptstep)*ztgo(:,1)  + dek%XLWNET_G(:) + dek%XSWNET_G(:)              &
              - plwnet_g_dtv(:)*ztvo(:) - plwnet_g_dtg(:)*ztgo(:,1) - plwnet_g_dtn(:)*ztno(:,1)     &
              - (1.0-pek%XPSN(:))*pflxc_g_c(:)*( pthrmb_tg(:)-pthrmb_tc(:)-pthrma_tc(:)*zcoefa_tc(:)  &
              + pltt(:)*(zhs(:)*(pqsat_g(:) - pdqsat_g(:)*ztgo(:,1))                                  &
              - zhn(:)*zcoefa_qc(:)) )                                                              &
              + ztconda_delz_g(:)*zsoil_coef_b(:,2)                                                 &
              + pek%XPSN(:)*ztconda_delz_ng(:)*ztno(:,jnsnow) )/zwork(:)

zalpha_g(:) =  (plwnet_g_dtv(:) + pflxc_g_c(:)*(1.0-pek%XPSN(:))*( pthrma_tc(:)*zcoefb_tc(:)       &
              + pltt(:)*zhn(:)*zcoefb_qc(:) ) )/zwork(:)

zgamma_g(:) =  (plwnet_g_dtn(:) + pflxc_g_c(:)*(1.0-pek%XPSN(:))*( pthrma_tc(:)*zcoefd_tc(:)       &
              + pltt(:)*zhn(:)*zcoefd_qc(:) ) )/zwork(:)

! Tn coefs, where TN = BETA_N + ALPHA_N*TV + GAMMA_N*TG

zwork(:)    =   (pcheatn(:)/ptstep) - zrnet_n_dtnn(:)                                                     &
              + pflxc_n_c(:)*(1.-ppsna(:))*( pthrma_tn(:) - pthrma_tc(:)*zcoefd_tc(:)                     &
              + (pltt(:)*zhns)*(pdqsati_n(:) - zcoefd_qc(:)) )                                              &
              + pflxc_n_a(:)*    ppsna(:)*(                                                               &
                pthrma_tn(:) - pthrma_ta(:)*(zpet_a_coef_p(:)*zcoefd_tc(:) + zpet_c_coef_p(:))            &
              + (pltt(:)*zhns)*(pdqsati_n(:)*(1.0-zpeq_c_coef_p(:)) - zpeq_a_coef_p(:)*zcoefd_qc(:)) )      &
              + ztconda_delz_n(:)*(1.0-zsnow_coef_a(:,2))

zbeta_n(:)  = ( (pcheatn(:)/ptstep)*ztno(:,1)  + zrnet_nn(:) + dmk%XHPSNOW(:) + dek%XMELTADV(:)        &
              - zrnet_n_dtvn(:)*ztvo(:) - zrnet_n_dtgn(:)*ztgo(:,1) - zrnet_n_dtnn(:)*ztno(:,1)           &
              - pflxc_n_c(:)*(1.-ppsna(:))*( pthrmb_tn(:) - pthrmb_tc(:) - pthrma_tc(:)*zcoefa_tc(:)      &
              + (pltt(:)*zhns)*(pqsati_n(:) - pdqsati_n(:)*ztno(:,1)                                        &
              - zcoefa_qc(:)) )                                                                           &
              - pflxc_n_a(:)*ppsna(:)*( pthrmb_tc(:) - pthrmb_ta(:) -                                     &
                                          pthrma_ta(:)*(zpet_b_coef_p(:) + zcoefa_tc(:)*zpet_a_coef_p(:)) &
              + (pltt(:)*zhns)*((pqsati_n(:) - pdqsati_n(:)*ztno(:,1))*(1.0-zpeq_c_coef_p(:))               &
              - zpeq_b_coef_p(:) - zpeq_a_coef_p(:)*zcoefa_qc(:)) )                                       &
              + ztconda_delz_n(:)*zsnow_coef_b(:,2) )/zwork(:)

zalpha_n(:) = ( zrnet_n_dtvn(:) + pflxc_n_c(:)*(1.-ppsna(:))*( pthrma_tc(:)*zcoefb_tc(:)                  &
              + (pltt(:)*zhns)*zcoefb_qc(:) )                                                               &
              + pflxc_n_a(:)*    ppsna(:) *( pthrma_ta(:)*zcoefb_tc(:)*zpet_a_coef_p(:)                   &
              + (pltt(:)*zhns)*zcoefb_qc(:)*zpeq_a_coef_p(:) ) )/zwork(:)

zgamma_n(:) = ( zrnet_n_dtgn(:) + pflxc_n_c(:)*(1.-ppsna(:))*( pthrma_tc(:)*zcoefc_tc(:)                  &
              + (pltt(:)*zhns)*zcoefc_qc(:))                                                                &
              + pflxc_n_a(:)*    ppsna(:) *( pthrma_ta(:)*zcoefc_tc(:)*zpet_a_coef_p(:)                   &
              + (pltt(:)*zhns)*zcoefc_qc(:)*zpeq_a_coef_p(:) ) )/zwork(:)

!*       7.     Solve multiple energy budgets simultaneously (using an implicit time scheme)
!               ----------------------------------------------------------------------------
! Also compute certain needed diagnostics, like the canopy air T and q, and T and q at
! the lowest model level: They are used within the current patch to compute implicit or explicit
! fluxes (depending upon the values of the atmospheric Implicit Coupling (IC) coefficients)


WHERE(pek%XPSN(:) > 0.0)

   zwork(:)      = 1.0 - zalpha_v(:)*zalpha_g(:)
   zbeta_p_v(:)  = (zbeta_v(:)  + zalpha_v(:)*zbeta_g(:) )/zwork(:)
   zalpha_p_v(:) = (zgamma_v(:) + zalpha_v(:)*zgamma_g(:))/zwork(:)

   zwork(:)      = 1.0 - zgamma_g(:)*zgamma_n(:)
   zbeta_p_n(:)  = (zbeta_n(:)  + zgamma_n(:)*zbeta_g(:) )/zwork(:)
   zalpha_p_n(:) = (zalpha_n(:) + zgamma_n(:)*zalpha_g(:))/zwork(:)

   dmk%XSNOWTEMP(:,1)      = (zbeta_p_n(:) + zalpha_p_n(:)*zbeta_p_v(:))/           &
                   (1.0 - zalpha_p_n(:)*zalpha_p_v(:)      )

! Since the fluxes are passed to the snow scheme, we can simply limit Tn here
! for flux computations (to make sure fluxes correspond to a snow sfc which
! doesn't exceed it's physical limit, Tf). The new real snow Tn consistent with these fluxes
! (and the T-profile within the snow) will be computed within the snow scheme.

   psnowliq(:,1) = psnowliq(:,1) + &
                   max(0., (dmk%XSNOWTEMP(:,1)-xtt)*psnowhcapz(:,1)*dmk%XSNOWDZ(:,1)/(xlmtt*xrholw)) ! m

   dmk%XSNOWTEMP(:,1)      = min(xtt, dmk%XSNOWTEMP(:,1))

   pek%XTV(:)        = zbeta_p_v(:) + zalpha_p_v(:)*dmk%XSNOWTEMP(:,1)

   ptg(:,1)      = zbeta_g(:) + zalpha_g(:)*pek%XTV(:) + zgamma_g(:)*dmk%XSNOWTEMP(:,1)

   pek%XTC(:)        = zcoefa_tc(:) + zcoefb_tc(:)*pek%XTV(:) + zcoefc_tc(:)*ptg(:,1) + zcoefd_tc(:)*dmk%XSNOWTEMP(:,1) 

   pek%XQC(:)        = zcoefa_qc(:) + zcoefb_qc(:)*pek%XTV(:) + zcoefc_qc(:)*ptg(:,1) + zcoefd_qc(:)*dmk%XSNOWTEMP(:,1) 

! Lowest atmospheric level air temperature (for this patch):

   pta_ic(:)     = zpet_b_coef_p(:) + zpet_a_coef_p(:) *pek%XTC(:) + zpet_c_coef_p(:) *dmk%XSNOWTEMP(:,1)

! Lowest atmospheric level specific humidity (for this patch):
! - First, compute the surface specific humidity of the snow:

   zwork(:)      = zhns(:)*( pqsati_n(:) + pdqsati_n(:)*(dmk%XSNOWTEMP(:,1) - ztno(:,1)) ) ! q_n+

   pqa_ic(:)     = zpeq_b_coef_p(:) + zpeq_a_coef_p(:) *pek%XQC(:) + zpeq_c_coef_p(:) *zwork(:)

! diagnostic snow thermal flux: surface to sub-surface thermal transfer

!   ZRESTOREN(:)  = PEK%XPSN(:)*ZTCONDA_DELZ_N(:)*(DMK%XSNOWTEMP(:)*(1.0-ZSNOW_COEF_A(:,2)) - ZSNOW_COEF_B(:,2))

ELSEWHERE ! snow free canopy-understory case:

   ptg(:,1)      = (zbeta_g(:) + zalpha_g(:)*zbeta_v(:))/           &
                   (1.0 - zalpha_g(:)*zalpha_v(:)      )

   pek%XTV(:)        = zbeta_v(:) + zalpha_v(:)*ptg(:,1)

   pek%XTC(:)        = zcoefa_tc(:) + zcoefb_tc(:)*pek%XTV(:) + zcoefc_tc(:)*ptg(:,1) 

   pek%XQC(:)        = zcoefa_qc(:) + zcoefb_qc(:)*pek%XTV(:) + zcoefc_qc(:)*ptg(:,1) 

! Lowest atmospheric level air temperature (for this patch):

   pta_ic(:)     = zpet_b_coef_p(:) + zpet_a_coef_p(:) *pek%XTC(:)

! Lowest atmospheric level specific humidity (for this patch):

   pqa_ic(:)     = zpeq_b_coef_p(:) + zpeq_a_coef_p(:) *pek%XQC(:)

! arbitrary: (as no snow mass present)

   dmk%XSNOWTEMP(:,1)      = xtt

!   ZRESTOREN(:)  = 0.0

END WHERE
!
!*       8.     Solve for sub-surface test snow temperature profile 
!               -----------------------------------------------------------------------------
! Compute test sub-surface snow temperatures: this improves time split estimate of 
! surface to sub-surface flux estimates. Note that the sub-surface
! snow temperatures are "test" temperatures, with final "true" values
! computed within the snow routine. 
! But we also include simple hydrology and refreezing since this can have
! a significant impact during melt events on the sub-surface test
! snow Temperature profile...
!
!
 CALL tridiag_ground_rm_soln(zt,zcoef_a,zcoef_b)
!
! Update Test T and Liquid content of snow:
!
DO jk=2,jnsnow
   DO jj=1,jnpts
      IF(pek%XPSN(jj) > 0.0)THEN
         psnowliq(jj,jk)     = psnowliq(jj,jk) + max(0., (zt(jj,jk)-xtt)* &
                               psnowhcapz(jj,jk)*dmk%XSNOWDZ(jj,jk)/(xlmtt*xrholw)) ! m
         dmk%XSNOWTEMP(jj,jk) = min(xtt,zt(jj,jk))
      ENDIF
   ENDDO
ENDDO
!
! Tipping bucket snow hydrology: update Liq
! NOTE this mimicks what is assumed to be done in the snow scheme. If a more sophisticated
! snow hydrology scheme is used, this code should be adapted.
!
zwholdmax(:,:)        = snow3lhold(psnowrho,dmk%XSNOWDZ) ! m
zwork(:)              = max(0., psnowliq(:,1)-zwholdmax(:,1))
psnowliq(:,1)         = psnowliq(:,1) - zwork(:)
DO jk=2,jnsnow
   DO jj=1,jnpts
      psnowliq(jj,jk) = psnowliq(jj,jk) + zwork(jj)
      zwork(jj)       = max(0., psnowliq(jj,jk)-zwholdmax(jj,jk))
      psnowliq(jj,jk) = psnowliq(jj,jk) - zwork(jj)
   ENDDO
ENDDO
psnowliq(:,:) = max(0.0, psnowliq(:,:)) ! numerical check
!
! Now, possibly refreeze liquid water which has flowed into a layer,
! thus reducing liquid water and warming the snowpack:
!
DO jk=2,jnsnow
   DO jj=1,jnpts
      IF(pek%XPSN(jj) > 0.0)THEN
         dmk%XSNOWTEMP(jj,jk) = dmk%XSNOWTEMP(jj,jk) + psnowliq(jj,jk)*(xlmtt*xrholw)/ &
                                (psnowhcapz(jj,jk)*max(1.e-6,dmk%XSNOWDZ(jj,jk)))     ! K
         zwork(jj)            = max(0., (xtt-dmk%XSNOWTEMP(jj,jk))*                    &
                                psnowhcapz(jj,jk)*dmk%XSNOWDZ(jj,jk)/(xlmtt*xrholw)) ! m 
         psnowliq(jj,jk)      = psnowliq(jj,jk) - zwork(jj)
         dmk%XSNOWTEMP(jj,jk) = min(xtt,dmk%XSNOWTEMP(jj,jk))
      ENDIF
   ENDDO
ENDDO
!
!*       9.     Solve soil temperature profile (implicitly coupled to surface energy budgets)
!               -----------------------------------------------------------------------------
! The back-substitution of sub-surface soil T profile
! (update all temperatures *below* the surface)
!
IF(io%CISBA == 'DIF')THEN
   CALL tridiag_ground_rm_soln(ptg,zsoil_coef_a,zsoil_coef_b)
ELSE
   ptg(:,2) = zsoil_coef_b(:,2) + zsoil_coef_a(:,2)*ptg(:,1)
ENDIF
!
! Diagnose (semi-implicit) flux across ground base:
! (semi-implicit if T imposed, explicit if flux imposed)
!
WHERE(kk%XTDEEP(:) == xundef)
   pdeep_flux(:) = 0.0
ELSEWHERE
   zwork(:)      = psoilcondz(:,jngrnd)*2/(pd_g(:,jngrnd)-pd_g(:,jngrnd-1)) ! (W/m2/K)
   pdeep_flux(:) = zwork(:)*(kk%XTDEEP(:) - ptg(:,jngrnd))/                 &
                              (1. - zwork(:)*ptdeep_a(:))                   ! (W/m2)
END WHERE
!
!*       10.     Temperature tendencies
!               ----------------------
!
! Compute energy budget tendencies (for flux computations) 
! and sub-sfc soil tendencies (K) (for phase changes)

pdeltat_g(:)   = ptg(:,1) - ztgo(:,1)
pdeltat_v(:)   = pek%XTV(:)  - ztvo(:)
pdeltat_n(:)   = dmk%XSNOWTEMP(:,1) - ztno(:,1)
!
!
!*      11.     Flux Diagnostics
!               ----------------
!
! Locally implicit (with respect to the ground T) snow-ground heat flux (W m-2).
! This is a first estimate, since the final fully implicit flux
! is computed in the snow routine: any differences
! between this flux estimate and the final one are
! transferred to the soil as heating/cooling, thus
! conserving energy.
!
pgrndflux(:)    = pek%XPSN(:)*ztconda_delz_ng(:)*( dmk%XSNOWTEMP(:,jnsnow) - ptg(:,1) )
!
!
!*      12.     Energy Storage Diagnostics (W m-2)
!               ----------------------------------
!
pdelheatg_sfc(:) = pcheatg(:)*pdeltat_g(:)/ptstep 
pdelheatv_sfc(:) = pcheatv(:)*pdeltat_v(:)/ptstep 
!
IF(io%CISBA == 'DIF')THEN

! Flux between surface and sub-surface (W m-2):

   dek%XRESTORE(:)  = (ptg(:,1) - ptg(:,2))* 2/( ((pd_g(:,2)-pd_g(:,1))/psoilcondz(:,2)) +     &
                                                   ( pd_g(:,1)           /psoilcondz(:,1)) )
!
!
!*      12.a    Energy Storage Diagnostics (W m-2): DIF
!               ---------------------------------------
!
! These terms are used for computing energy budget diagnostics.
!
!
   pdelheatg(:)        = pdelheatg_sfc(:) ! initialize
   DO jk=2,jngrnd
      DO jj=1,jnpts
         pdelheatg(jj) =  pdelheatg(jj) + psoilhcapz(jj,jk)*(pd_g(jj,jk)-pd_g(jj,jk-1))*    &
                          (ptg(jj,jk) - ztgo(jj,jk))/ptstep
      ENDDO
   ENDDO
!
ELSE

!*      12.b    Energy Storage Diagnostics (W m-2): Force Restore
!               -------------------------------------------------
!
! Flux between surface and sub-surface (W m-2):

   dek%XRESTORE(:)  = (2*xpi/xday)*(ptg(:,1) - ptg(:,2))/pct(:)

! These terms are used for computing energy budget diagnostics.
!
   pdelheatg(:) = pdelheatg_sfc(:) + dek%XRESTORE(:)
!
ENDIF

!
!*      13.     Additional Diagnostics
!               ----------------------
!
! Here compute the updated (time t+dt) effective reference level specific heat capacity:
! (just a diagnostic here):
!
IF (io%CCPSURF=='DRY') THEN
   pk%XCPS(:) = xcpd
ELSEIF(.NOT.lcpl_arp)THEN
   pk%XCPS(:) = xcpd + ( xcpv - xcpd ) * pek%XQC(:)
ENDIF
!
! These latent heats are used to convert between
! latent heat and water mass fluxes:
!
pk%XLVTT(:)   = pltt(:)
pk%XLSTT(:)   = pltt(:) 
!
!
IF (lhook) CALL dr_hook('E_BUDGET_MEB',1,zhook_handle)
!
END SUBROUTINE e_budget_meb