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(HISBA,HCPSURF,PTSTEP,                                           &
           pps,pcg,pct,pcv,pwrvn,pwr,                                                         &
           ptdeep_a,ptdeep_b,pd_g,psoilcondz,psoilhcapz,                                      &  
           psnowdz,psnowcondz,psnowhcapz,                                                     &  
           pswnet_v,pswnet_g,pswnet_n,ptau_n,                                                 &
           plwnet_v,plwnet_g,plwnet_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,                                                          &
           pff,pffrozen,ppsn,ppsna,ppsncv,                                                    &
           pcheatv,pcheatg,pcheatn,                                                           &
           pleg_delta,plegi_delta,phug,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,ptv,ptn,                                                                       &
           pflxc_v_c,phvgs,phvns,                                                             & 
           pdqsat_g,pdqsat_v,pdqsati_n,                                                       & 
           ptc,pqc,pta_ic,pqa_ic,pustar2_ic,pvmod,                                            &
           pdeltat_g,pdeltat_v,pdeltat_n,pgrndflux,pcps,plvtt,plstt,                          &
           phpsnow,pmeltadv,prestore,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_csts,                    ONLY : xlvtt, xlstt, xtt, xcpd, xcpv, xcl, &
                                         xday, xpi
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 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
!
 CHARACTER(LEN=*),      INTENT(IN)  :: hisba     ! type of ISBA version:
!                                               ! '2-L' (default)
!                                               ! '3-L'
!                                               ! 'DIF'
!
 CHARACTER(LEN=*),      INTENT(IN)  :: hcpsurf   ! Specific heat
!                                               ! 'DRY' = dry Cp
!                                               ! 'HUM' = humid Cp fct of qc for MEB
!
REAL,                 INTENT(IN)   :: ptstep
!                                     PTSTEP = timestep of the integration (s)
!
REAL, DIMENSION(:),   INTENT(IN)   :: pct, pcv, pwrvn, pwr, pcg
!                                     PCG     = area-averaged soil heat capacity (m2 K J-1)
!                                     PCT     = surface understory or composite thermal inertia (m2 K J-1)
!                                     PCV     = vegetation canopy thermal inertia (m2 K J-1)
!                                     PWRVN   = liquid water equivalent snow  intercepted on the canopy (kg m-2)
!                                     PWR     = liquid water intercepted on the canopy (kg m-2)
!
REAL, DIMENSION(:), INTENT(IN)     :: ptdeep_a, ptdeep_b
!                                      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)
!
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)   :: psnowdz, psnowcondz, psnowhcapz
!                                     PSNOWDZ    = snow layer thickness      (m)
!                                     PSNOWCONDZ = snow thermal conductivity (W m-1 K-1)
!                                     PSNOWHCAPZ = snow heat capacity        (J m-3 K-1)
!
REAL, DIMENSION(:),   INTENT(IN)   :: pswnet_v, pswnet_g, pswnet_n
!                                     PSWNET_G = Understory-ground net SW radiation explicit term (W m-2)
!                                     PSWNET_V = Vegetation canopy net SW radiation explicit term (W m-2)
!                                     PSWNET_N = Ground-based snow net SW radiation explicit term (W m-2)
!
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, plwnet_g, plwnet_n
!                                     PLWNET_G = Understory-ground net LW radiation explicit term (W m-2)
!                                     PLWNET_V = Vegetation canopy net LW radiation explicit term (W m-2)
!                                     PLWNET_N = Ground-based snow net LW radiation explicit term (W m-2)
!
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)   :: pff, pffrozen, ppsn, ppsna, ppsncv
!                                     PFF      = total flooded fraction                                        (-) 
!                                     PFFROZEN = total frozen flooded fraction                                 (-) 
!                                     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, phug, phugi, phvg, phvn, pfrozen1
!                                     PHUG        = relative humidity of the surface soil                             (-)  
!                                     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(INOUT):: ptn
!                                     PTN    = Ground-based snow temperature profile (K)
!
REAL, DIMENSION(:),   INTENT(INOUT):: ptv
!                                     PTV    = Vegetation canopy temperature (K)
!
REAL, DIMENSION(:,:), INTENT(INOUT):: ptg
!                                     PTG    = Soil temperature profile (K)
!
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)  :: ptc, pqc, pta_ic, pqa_ic, pustar2_ic, pvmod
!                                     PTC         = Canopy air space temperature       (K)
!                                     PQC         = Canopy air space specific humidity (kg kg-1)
!                                     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)  :: pcps, plvtt, plstt
!                                     PCPS = specific heat capacity of the surface (J kg-1 K-1) 
!                                     PLVTT= latent heat of vaporization (J kg-1)
!                                     PLSTT= latent heat of sublimation  (J kg-1)
!
REAL, DIMENSION(:),   INTENT(OUT)  :: phpsnow, pmeltadv, prestore
!                                     PHPSNOW = Precipitation heating term for ground-based snowpack (W m-2)
!                                     PMELTADV= Heating resulting when meltwater enters soil at sub-freezing
!                                               temperatures (W m-2)
!                                     PRESTORE= flux between surface and sub-surface in soil (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                                      :: zhns
!
REAL, DIMENSION(SIZE(PTG,1),SIZE(PTG,2))  :: ztgo
!
REAL, DIMENSION(SIZE(PTN,1),SIZE(PTN,2))  :: ztno
!
REAL, DIMENSION(SIZE(PTG,1))              :: ztvo
!
REAL, DIMENSION(SIZE(PPS))                :: zhn, zhs, zhvs
!
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(PSNOWDZ,1),SIZE(PSNOWDZ,2)):: zsnow_coef_a, zsnow_coef_b, zsnowdz
!
REAL, DIMENSION(SIZE(PD_G,1),SIZE(PD_G,2)+SIZE(PSNOWDZ,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(:,:) = ptn(:,:)
ztvo(:)   = ptv(:)

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

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

zsnowdz(:,:) = 0.0 ! local working variable

jnsnow       = SIZE(ptn,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(pwrvn,pwr,pcv)                    ! vegetation canopy heat capacity
!
pcheatn(:)    = psnowhcapz(:,1)*psnowdz(:,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.
!
phpsnow(:)   = 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.
!
pmeltadv(:)  = 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(hisba == '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,ptdeep_b,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:
!
   zsnowdz(:,:)           = max(xsnowdzmin, psnowdz(:,:))
!

   jl                     = 1
   zd(:,jl)               = zsnowdz(:,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)        = zsnowdz(jj,jk)
         zt(jj,jl)        = ztno(jj,jk)
         zhcapz(jj,jl)    = psnowhcapz(jj,jk)
         zcondz(jj,jl)    = psnowcondz(jj,jk)
         zsource(jj,jl)   = pswnet_n(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)          = pswnet_n(:)*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,ptdeep_b,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)/( pcg(:)*pcg(:)/(pd_g(:,1)*pct(:)) )


   jl                    = 1
   zd(:,jl)              = zsnowdz(:,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)       = zsnowdz(jj,jk)
         zt(jj,jl)       = ztno(jj,jk)
         zhcapz(jj,jl)   = psnowhcapz(jj,jk)
         zcondz(jj,jl)   = psnowcondz(jj,jk)
         zsource(jj,jl)  = pswnet_n(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)         = pswnet_n(:)*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,ptdeep_b,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/((zsnowdz(:,jnsnow)/psnowcondz(:,jnsnow))+(pd_g(:,1)/zgcond1(:) ))
!
!
!*       3.     Pseudo humidity factors (-)
!               ---------------------------------------------------------
!
! First, 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.-ppsn(:)) + pflxc_vn_c(:)*ppsn(:)*(1.-zpsna(:))
pflxc_v_c(:)   = max(pflxc_v_c(:), zertol_flx_c)

! Understory vegetation and ground factors:

zfff(:)        = pff(:)*( 1.0 - pffrozen(:)*(1.0 - (xlstt/xlvtt)) )
 
zhn(:)         = (1.0-ppsn(:)-zfff(:))*(                                             &
                            pleg_delta(:) *(1.0-pfrozen1(:))                         &
                 +          plegi_delta(:)*     pfrozen1(:)*(xlstt/xlvtt) ) + zfff(:)         

zhs(:)         = (1.0-ppsn(:)-zfff(:))*(                                             &
                   phug(:) *pleg_delta(:) *(1.0-pfrozen1(:))                         &
                 + phugi(:)*plegi_delta(:)*     pfrozen1(:)*(xlstt/xlvtt) ) + zfff(:)

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

zhn(:)         = zhn(:)/max(1.0 - ppsn(:) , zertol)

zhs(:)         = zhs(:)/max(1.0 - ppsn(:) , zertol)

! Vegetation canopy factor:

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

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

! - total canopy H factor (including intercepted snow)

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

! Snow latent heating factor:
! NOTE, for now we consider only snow sublimation 
!       (not evaporation from snow liquid)
!
zhns           = (xlstt/xlvtt)
!
!
!*       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 - ppsn(:)*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(:)*ppsn(:)*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(:)*ppsn(:)*ppsna(:)*(pthrmb_tn(:)-pthrmb_ta(:)) ) )   /zwork(:)
zpet_c_coef_p(:)  =   ppet_a_coef(:)*pflxc_n_a(:)*ppsn(:)*ppsna(:)*pthrma_tn(:)       /zwork(:)

! q coefficients:

zwork(:)          = 1.0 + ppeq_a_coef(:)*(pflxc_c_a(:)*zpsnag(:) +                               &
                                          pflxc_n_a(:)*ppsn(:)*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(:)*ppsn(:)*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-ppsn(:))                                       &
                   + pflxc_n_c(:) *pthrma_tc(:)*     ppsn(:) *(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-ppsn(:))                       &
                   + pflxc_n_c(:) * (pthrmb_tn(:)-pthrmb_tc(:))*     ppsn(:) *(1.0-ppsna(:))        &
                                                                      )/zwork(:)

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

zcoefc_tc(:) = pflxc_g_c(:)*pthrma_tg(:)*(1.0-ppsn(:))                 /zwork(:)

zcoefd_tc(:) =(pflxc_n_c(:)*pthrma_tn(:)*          ppsn(:)*(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-ppsn(:))                                            &
                   + pflxc_n_c(:) *zhns   *     ppsn(:) *(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-ppsn(:))           &
                 + pflxc_n_c(:) *zhns   *(pqsati_n(:)-pdqsati_n(:)*ztno(:,1))*                      &
                                                                      ppsn(:)*(1.0-ppsna(:))        &
                                                                                )/zwork(:) 

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

zcoefc_qc(:)   = pflxc_g_c(:) *zhs(:) *pdqsat_g(:)*(1.0-ppsn(:))              /zwork(:)

zcoefd_qc(:)   = pflxc_n_c(:) *zhns   *pdqsati_n(:)*     ppsn(:)*(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,       ppsn(:))

zrnet_nn(:)     = (pswnet_n(:) + plwnet_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(:)                              &
              + xlvtt*zhvs(:)*(pdqsat_v(:) - zcoefb_qc(:)) )           

zbeta_v(:)  = ( (pcheatv(:)/ptstep)*ztvo(:) + plwnet_v(:) + pswnet_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(:)                  &
              + xlvtt*zhvs(:)*(pqsat_v(:) - pdqsat_v(:)*ztvo(:)                                     &
              - zcoefa_qc(:)) ) )/zwork(:)

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

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

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

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

zbeta_g(:)  = ( (pcheatg(:)/ptstep)*ztgo(:,1)  + plwnet_g(:) + pswnet_g(:)                          &
              - plwnet_g_dtv(:)*ztvo(:) - plwnet_g_dtg(:)*ztgo(:,1) - plwnet_g_dtn(:)*ztno(:,1)     &
              - (1.0-ppsn(:))*pflxc_g_c(:)*( pthrmb_tg(:)-pthrmb_tc(:)-pthrma_tc(:)*zcoefa_tc(:)    &
              + xlvtt*(zhs(:)*(pqsat_g(:) - pdqsat_g(:)*ztgo(:,1))                                  &
              - zhn(:)*zcoefa_qc(:)) )                                                              &
              + ztconda_delz_g(:)*zsoil_coef_b(:,2)                                                 &
              + ppsn(:)*ztconda_delz_ng(:)*ztno(:,jnsnow) )/zwork(:)

zalpha_g(:) =  (plwnet_g_dtv(:) + pflxc_g_c(:)*(1.0-ppsn(:))*( pthrma_tc(:)*zcoefb_tc(:)             &
              + xlvtt*zhn(:)*zcoefb_qc(:) ) )/zwork(:)

zgamma_g(:) =  (plwnet_g_dtn(:) + pflxc_g_c(:)*(1.0-ppsn(:))*( pthrma_tc(:)*zcoefd_tc(:)             &
              + xlvtt*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(:)                     &
              + (xlvtt*zhns)*(pdqsati_n(:) - zcoefd_qc(:)) )                                              &
              + pflxc_n_a(:)*    ppsna(:)*(                                                               &
                pthrma_tn(:) - pthrma_ta(:)*(zpet_a_coef_p(:)*zcoefd_tc(:) + zpet_c_coef_p(:))            &
              + (xlvtt*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(:) + phpsnow(:) + pmeltadv(:)                   &
              - 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(:)      &
              + (xlvtt*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(:)) &
              + (xlvtt*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(:)                  &
              + (xlvtt*zhns)*zcoefb_qc(:) )                                                               &
              + pflxc_n_a(:)*    ppsna(:) *( pthrma_ta(:)*zcoefb_tc(:)*zpet_a_coef_p(:)                   &
              + (xlvtt*zhns)*zcoefb_qc(:)*zpeq_a_coef_p(:) ) )/zwork(:)

zgamma_n(:) = ( zrnet_n_dtgn(:) + pflxc_n_c(:)*(1.-ppsna(:))*( pthrma_tc(:)*zcoefc_tc(:)                  &
              + (xlvtt*zhns)*zcoefc_qc(:))                                                                &
              + pflxc_n_a(:)*    ppsna(:) *( pthrma_ta(:)*zcoefc_tc(:)*zpet_a_coef_p(:)                   &
              + (xlvtt*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(ppsn(:) > 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(:)

   ptn(:,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.

   ptn(:,1)      = min(xtt, ptn(:,1))

   ptv(:)        = zbeta_p_v(:) + zalpha_p_v(:)*ptn(:,1)

   ptg(:,1)      = zbeta_g(:) + zalpha_g(:)*ptv(:) + zgamma_g(:)*ptn(:,1)

   ptc(:)        = zcoefa_tc(:) + zcoefb_tc(:)*ptv(:) + zcoefc_tc(:)*ptg(:,1) + zcoefd_tc(:)*ptn(:,1) 

   pqc(:)        = zcoefa_qc(:) + zcoefb_qc(:)*ptv(:) + zcoefc_qc(:)*ptg(:,1) + zcoefd_qc(:)*ptn(:,1) 

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

   pta_ic(:)     = zpet_b_coef_p(:) + zpet_a_coef_p(:) *ptc(:) + zpet_c_coef_p(:) *ptn(:,1)

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

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

   pqa_ic(:)     = zpeq_b_coef_p(:) + zpeq_a_coef_p(:) *pqc(:) + zpeq_c_coef_p(:) *zwork(:)

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

!   ZRESTOREN(:)  = PPSN(:)*ZTCONDA_DELZ_N(:)*(PTN(:,1)*(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(:)      )

   ptv(:)        = zbeta_v(:) + zalpha_v(:)*ptg(:,1)

   ptc(:)        = zcoefa_tc(:) + zcoefb_tc(:)*ptv(:) + zcoefc_tc(:)*ptg(:,1) 

   pqc(:)        = zcoefa_qc(:) + zcoefb_qc(:)*ptv(:) + zcoefc_qc(:)*ptg(:,1) 

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

   pta_ic(:)     = zpet_b_coef_p(:) + zpet_a_coef_p(:) *ptc(:) 

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

   pqa_ic(:)     = zpeq_b_coef_p(:) + zpeq_a_coef_p(:) *pqc(:) 

! arbitrary: (as no snow mass present)

   ptn(:,1)      = xtt

!   ZRESTOREN(:)  = 0.0

END WHERE
!
!
! 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.
!

 CALL tridiag_ground_rm_soln(zt,zcoef_a,zcoef_b)
DO jk=2,jnsnow
   WHERE(ppsn(:) > 0.0)
      ptn(:,jk) = min(xtt,zt(:,jk))
   ENDWHERE
ENDDO
!
!
!*       8.     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(hisba == '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(ptdeep_b(:) == xundef)
   pdeep_flux(:) = 0.0
ELSEWHERE
   zwork(:)      = psoilcondz(:,jngrnd)*2/(pd_g(:,jngrnd)-pd_g(:,jngrnd-1)) ! (W/m2/K)
   pdeep_flux(:) = zwork(:)*(ptdeep_b(:) - ptg(:,jngrnd))/                 &
                              (1. - zwork(:)*ptdeep_a(:))                   ! (W/m2)
END WHERE
!
!*       9.     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(:)   = ptv(:)   - ztvo(:)
pdeltat_n(:)   = ptn(:,1) - ztno(:,1)
!
!
!*      10.     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(:)    = ppsn(:)*ztconda_delz_ng(:)*( ptn(:,jnsnow) - ptg(:,1) )
!
!
!*      10.     Energy Storage Diagnostics (W m-2)
!               ----------------------------------
!
pdelheatg_sfc(:) = pcheatg(:)*pdeltat_g(:)/ptstep 
pdelheatv_sfc(:) = pcheatv(:)*pdeltat_v(:)/ptstep 
!
IF(hisba == 'DIF')THEN

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

   prestore(:)  = (ptg(:,1) - ptg(:,2))* 2/( ((pd_g(:,2)-pd_g(:,1))/psoilcondz(:,2)) +     &
                                             ( pd_g(:,1)           /psoilcondz(:,1)) )
!
!
!*      10.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

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

   prestore(:)  = (2*xpi/xday)*(ptg(:,1) - ptg(:,2))/pct(:)

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

!
!*      11.     Additional Diagnostics
!               ----------------------
!
! Here compute the updated (time t+dt) effective reference level specific heat capacity:
! (just a diagnostic here):
!
IF (hcpsurf=='DRY') THEN
   pcps(:) = xcpd
ELSEIF(.NOT.lcpl_arp)THEN
   pcps(:) = xcpd + ( xcpv - xcpd ) * pqc(:)   
ENDIF
!
plvtt(:)   = xlvtt
plstt(:)   = xlvtt ! latent heat of sublimation already accounted for
!
!
IF (lhook) CALL dr_hook('E_BUDGET_MEB',1,zhook_handle)
!
END SUBROUTINE e_budget_meb