Energy Balance

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 dozen of W/m2 still remains. Some energy seems to be lost in explicitsnow_melt_refrz, we are still investigating a precise case study on a pixel.