enerbil_fusion transforms snowmet into temperature change even if it is not true melt

[jpolcher]

When snow is melted in hydrol or hydrolc, enerbil_fusion will transform this tot_melt water flux into a temperature reduction.

The problem is that this temperature change is also done in enerbil_fusion also when the snowmelt comes because we have reached the maximum snow amount (maxmass_glacier). This will induce dramatic temperature drops as large tot_melt can be generated and which are unrelated to any temperature above the freezing point.

This can be avoided in enerbil_fusion by something like

IF (temp_sol_new(ji).GT.tp_00) THEN

!! \latexonly
!! \input{enerbilfusion1.tex}
!! \endlatexonly
fusion(ji) = tot_melt(ji) * chalfu0 / dtradia

!! \latexonly
!! \input{enerbilfusion2.tex}
!! \endlatexonly
temp_sol_new(ji) = temp_sol_new(ji) - MIN(((tot_melt(ji) * chalfu0) / soilcap(ji)), temp_sol_new(ji)-tp_00)

ENDIF

Obviously in this case we do not conserve energy ... but in any case it is impossible to melt snow above maxmass_glacier and conserve energy at the same time.

Another solution, and probably a cleaner one, would be to separate the flux which indeed comes from temperature above the freezing level from the other melt fluxes (the numerically related ones).

[jgipsl]

Note : The variable maxmass_glacier changed name into maxmass_snow [2053].

[jpolcher]

One option here would be to separate true melt from the melting to avoid excess snow (i.e. snow above maxmass_glacier).

In this case the energy associated to the "excess snow" melting could be added to the temperature of the lower layers of the soil in thermosoil. This would ensure energy conservation in the model and not perturb too much the surface temperature calculations.

[maignan]

A more sophisticated solution proposed by J. Polcher

diff hydrol.f90:

--- src_sechiba/hydrol.f90    (revision 2165)
+++ src_sechiba/hydrol.f90    (working copy)
@@ -2978,6 +2978,7 @@
     REAL(r_std), DIMENSION (kjpindex)             :: d_age  !! Snow age change
     REAL(r_std), DIMENSION (kjpindex)             :: xx     !! temporary
     REAL(r_std)                                   :: snowmelt_tmp !! The name says it all !
+    REAL(r_std)                                   :: snow_d1k !! The amount of snow that corresponds to a 1K cooling

 !_ ================================================================================================================================

@@ -3084,11 +3085,16 @@
        ENDIF
     !! 1.4 Snow melts only on weight glaciers
        ! Ice melt only if there is more than a given mass : maxmass_snow,
-       ! Ajouts Edouard Davin / Nathalie de Noblet add extra to melting
+       ! But the snow cannot melt more in one time step to what corresponds to
+       ! a 1K cooling. This will lead to a progressive melting of snow above
+       ! maxmass_snow but it is needed as a too strong cooling can destabilise the model.
        !
        IF ( snow(ji) .GT. maxmass_snow ) THEN
-          snowmelt(ji) = snowmelt(ji) + (snow(ji) - maxmass_snow)
-          snow(ji) = maxmass_snow
+          snow_d1k = un * soilcap(ji) / chalfu0
+          snowmelt(ji) = snowmelt(ji) + MIN((snow(ji) - maxmass_snow),snow_d1k)
+          snow(ji) = snow(ji) - snowmelt(ji)
+!!!          IF (long_print)
+          WRITE (numout,*) "Snow was above maxmass_snow (", maxmass_snow,") and we melted ", snowmelt(ji)
        ENDIF
        !
     END DO

There is a remaining print for testing.

[maignan]

A similar correction is to be added for the ice case:

    !! 2.4 Snow melts only on weight glaciers 
       !      Ice melt only if there is more than a given mass : maxmass_snow, 
       ! But the snow cannot melt more in one time step to what corresponds to
       ! a 1K cooling. This will lead to a progressive melting of snow above
       ! maxmass_snow but it is needed as a too strong cooling can destabilise the model.
       !
       IF ( snow_nobio(ji,iice) .GT. maxmass_snow ) THEN
          !icemelt(ji) = snow_nobio(ji,iice) - maxmass_snow
          !snow_nobio(ji,iice) = maxmass_snow
          snow_d1k = un * soilcap(ji) / chalfu0
          icemelt(ji) = MIN((snow_nobio(ji,iice) - maxmass_snow),snow_d1k)
          snow_nobio(ji,iice) = snow_nobio(ji,iice) - icemelt(ji)
!!!          IF (long_print)
          WRITE (numout,*) "Snow was above maxmass_snow ON ICE (", maxmass_snow,") and we melted ", icemelt(ji)
       ENDIF

With these two modifications on sections 1.4 and 2.4, no more problems are seen on global evapotranspiration on a run at 2 degrees (1901-2012).

[jgipsl]

Done in rev [2360] in hydrol + hydrolc.

Question for hydrolc_init. Is following section needed as snow and snow_nobio are previously read from restart file?

  !! 6. We limit the initial value of snow water equivalent 
    DO ik=1, kjpindex
       IF (snow(ik).gt.maxmass_snow) THEN
          WRITE(numout,*)' The stock of snow from restart file will be limited', &
               ' as it is higher than the limit maxmass_snow'
          snow(ik)=maxmass_snow
       END IF
       IF (snow_nobio(ik,iice).gt.maxmass_snow) THEN
          WRITE(numout,*)' The stock of snow_nobio from restart file will be limited', &
               ' as it is higher than the limit maxmass_snow'
          snow_nobio(ik,iice)=maxmass_snow
       END IF
    ENDDO

[jgipsl]

Jan agree.
Done rev [2361]