Energy Balance

Update 27/03/2018

Energy balance in the soil has to be checked too.

Update 05/09/2017

Working group: FM, CO, PP
The energy budget is not closed only when the amount of snow and the corresponding fraction are small, but it is difficult to pinpoint exactly what the problem is.
There are indeed several suspicious facts:

• A minimal threshold psnowdzmin is used in thermosoil_coef for the computation of the snow coefficients.
• The computations made in the explicit_snow_melt_refrz subroutine use the amount of snow and the snow temperature, which includes the additional temperature temp_sol_add computed in enerbil_flux from a surface energy budget equation, which does not take into account the amount of snow.

Meanwhile, to have outputs available for the group to look at, we could first commit the implementation of the diagnostics:

• in sechiba.90: TEBR_enerbil and TEBR_snow
• in explicit_snow.f90:

A TEBR can be computed for every subroutine to better understand why/where the budget is not closed. This implies additional computations of the total snow heat; the snowheat variable represents for each snow layer the energy required to bring the temperature to 0°C and to melt the non-liquid snow fraction. For example, we could compute and output the following variables:

    CALL xios_orchidee_send_field("snowheattot_begin",snowheattot_begin)
    CALL xios_orchidee_send_field("snowheattot_snowfall",snowheattot_snowfall)
    CALL xios_orchidee_send_field("snowheattot_redist",snowheattot_redist)
    CALL xios_orchidee_send_field("snowheattot_compact",snowheattot_compact)
    CALL xios_orchidee_send_field("snowheattot_profile",snowheattot_profile)
    CALL xios_orchidee_send_field("snowheattot_gone",snowheattot_gone)
    CALL xios_orchidee_send_field("snowheattot_melt_refrz",snowheattot_melt_refrz)
    CALL xios_orchidee_send_field("snowheattot_subli",snowheattot_subli)
    CALL xios_orchidee_send_field("del_snowheattot_snowfall",del_snowheattot_snowfall)
    CALL xios_orchidee_send_field("del_snowheattot_redist",del_snowheattot_redist)
    CALL xios_orchidee_send_field("del_snowheattot_compact",del_snowheattot_compact)
    CALL xios_orchidee_send_field("del_snowheattot_profile",del_snowheattot_profile)
    CALL xios_orchidee_send_field("del_snowheattot_gone",del_snowheattot_gone)
    CALL xios_orchidee_send_field("del_snowheattot_melt_refrz",del_snowheattot_melt_refrz)
    CALL xios_orchidee_send_field("del_snowheattot_subli",del_snowheattot_subli)
    CALL xios_orchidee_send_field("TEBR_explicitsnow_gone",TEBR_explicitsnow_gone) 
    CALL xios_orchidee_send_field("TEBR_explicitsnow_melt_refrz",TEBR_explicitsnow_melt_refrz)
    CALL xios_orchidee_send_field("TEBR_explicitsnow_subli",TEBR_explicitsnow_subli)
    CALL xios_orchidee_send_field("TEBR_explicitsnow_total",TEBR_explicitsnow_total)

Special cases are to be handled properly, like for example: when there is a snowfall that melts completely during the time step.

Update 16/05/2017

Working group: FM, CO, PP
A diagnostic analog to twbr is computed in sechiba.f90 after the call to enerbil_main:

    tebr(:)=swnet(:) + lwnet(:) - fluxsens(:) - fluxlat(:) + soilflx(:) - soilcap(:)*(temp_sol_new(:) - temp_sol(:))/dt_sechiba
    CALL xios_orchidee_send_field("TEBR",tebr)

The budget is closed here with a mean value on the order of 1e-15 W/m2 and min/max values on the order of -/+1e-13 W/m2 on a global run at 2 degrees with monthly outputs.

A second diagnostic is computed after the call to explicit_snow and the possible modification of temp_sol_new in thermosoil_main. This one is not closed with a mean bias of 0.9 W/m2.

A first modification is proposed at the end of thermosoil_main. Instead of a callback to tp_00 if snow is present, we do a callback to a weighted temperature, similar to the computation of temp_sol_eff:

    !! Surface temperature is forced to zero celcius if its value is larger than melting point, only for explicit snow scheme
    IF  (ok_explicitsnow) THEN
       DO ji=1,kjpindex
          IF  (SUM(snowdz(ji,:)) .GT. 0.0) THEN
             IF (temp_sol_new(ji) .GE. tp_00) THEN
                !temp_sol_new(ji) = tp_00
                temp_sol_new(ji) = snowtemp(ji,1)*frac_snow_veg(ji)*(1-totfrac_nobio(ji))+ &      ! weights related to first layer of snow cover fraction on vegetation
                       temp_sol_new(ji)*SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji)+ &           ! weights related to SCF on nobio
                       temp_sol_new(ji)*(1-(frac_snow_veg(ji)*(1-totfrac_nobio(ji))+SUM(frac_snow_nobio(ji,:))*totfrac_nobio(ji))) ! weights related to non snow fraction
             ENDIF
          END IF
       END DO
    ENDIF

The mean bias is reduced to 0.2 W/m2.

Other modifications are proposed in enerbil.f90 for the computation of pgflux, zgflux and temp_sol_add:
PHPSNOW is removed because the rain does not presently fall on the snow and, if considered, should also be taken into account in the general energy balance scheme.
Basically, pgflux is now equal to -soilflx (not sure why we're not dealing with the same sign convention).

      !! To calculate net energy flux into the snowpack
      IF (ok_explicitsnow) THEN
          !PHPSNOW(ji) = precip_rain(ji)*(4.218E+3)*(MAX(tp_00,temp_air(ji))-tp_00)/dt_sechiba ! (w/m2)
          !pgflux(ji)  = netrad(ji) - fluxsens(ji) - fluxlat(ji) + PHPSNOW(ji)
          pgflux(ji)  = netrad(ji) - fluxsens(ji) - fluxlat(ji) - soilcap(ji)*(temp_sol_new(ji) - temp_sol(ji))/dt_sechiba
      ENDIF

zgflux and temp_sol_add are computed considering a snowpack at 0°C and the heat capacity of the snow:

         !zgflux(ji)  = netrad_tmp(ji) - fluxsens_tmp(ji) - fluxlat_tmp(ji) - soilcap(ji)*(tp_00 - temp_sol(ji))/dt_sechiba
         zgflux(ji)  = netrad_tmp(ji) - fluxsens_tmp(ji) - fluxlat_tmp(ji) - snowcap(ji)*(tp_00 - temp_sol(ji))/dt_sechiba

         !temp_sol_add(ji) = -(pgflux(ji) - zgflux(ji))*dt_sechiba/soilcap(ji)
         temp_sol_add(ji) = -(pgflux(ji) - zgflux(ji))*dt_sechiba/snowcap(ji)

The mean bias is reduced to -0.04 W/m2, however min/max values on the order of several dozens of W/m2 still remain. Some energy seems to be lost in explicitsnow_melt_refrz, we are still investigating a precise case study on a pixel.