diaglev : trunk needs to be upated for stomate and ok_freeze after change in vertical soil discretization

[jgipsl]

Part 1) Subroutine hydrol_calculate_temp_hydro calculates temp_hydro using stempdiag and diaglev.
We think:

• in thermosoil_main: the only output variable should be ptn, so we can remove stempdiag and ptnlev1.
• In hydrol: ptn can be used directly instead of temp_hydro. Stempdiag is not needed any more and not the subroutine hydrol_calculate_temp_hydro.

Part 2) Subroutine hydrol_calculate_frac_hydro_diag also uses diaglev.
We think:

• If stomate can deal with the prognostic vertical axis instead of diaglev, then subroutine hydrol_calculate_frac_hydro_diag can be removed because the variable frac_hydro_diag becomes one. It is not needed any more in the calculations of shumdiag and shumdiag_perma.

Part 3) Currently diaglev is used in stomate to define z_soil. Before the commits of the vertical discretization the diaglev were corresponding to the depths of the bottom of the hydrological layers. After the first commit of the vertical discretization, diaglev is set by the subroutine vertical_getdiag to tdepth : the depths of the nodes of the thermal layers.

• We think this might be wrong. The stomate group needs to tell which variable(s) is needed in stomate instead of diaglev.

See also ticket for vertical discretization #190

Fuxing, Agnes, Josefine

[jpolcher]

Yes, this needs simplification now.

But PTN should not leave thermosoil as it is a prognostic variable of this module. Please keep still keep the variable stempdiag into which you copy ptn.

It is written in the coding guidelines that prognostic variables should not leave the module else they might get changed outside and thus not be restarted correctly or other stupid mistake happen.

[jgipsl]

Replying to jpolcher:

But PTN should not leave thermosoil as it is a prognostic variable of this module. Please keep still keep the variable stempdiag into which you copy ptn.

ptn is only "leaving" as an INTENT(OUT) variable. In hydrol, it is an INTENT(IN). It can not be changed outside the module thermosoil.

[jpolcher]

To me this is not enough protection. There will be somebody who will change the INTENT(OUT) into INTENT(INOUT) as he is changing PTN outside and the compiler is complaining.

The prognostic variables do not leave their module, this is a fundamental rule and makes total sense.
Jan

[aducharne]

This has been discussed during the technical meeting of Tuesday 22, 2015.
The conclusion is to overlook this comment, as the complexity it induces overcomes the security, and many prognostic variables in stomate are not treated like this.

[maignan]

Answer to point 1:

---

• Be careful regarding the snow model:
  • If the old snow model is used, stempdiag is still needed because the soil levels are modified in hydrol.f90 when snow is present;
  • otherwise with the new snow model, since the snow pack is above the soil surface, ptn is the only needed variable.

For the same reason, in subroutine hydrol_calculate_temp_hydro, we have to modify the present code:

       IF (ok_explicitsnow) THEN 
          snow_h=SUM(snowdz(ji,:))
       ELSE  
          snow_h=snow(ji)/sn_dens
       ENDIF

into:

       IF (ok_explicitsnow) THEN 
          snow_h=0
       ELSE  
          snow_h=snow(ji)/sn_dens
       ENDIF

to account for the fact that the snow pack is above the surface soil in the new snow model.

• Secondly, profil_froz_hydro should be kept for the moment, if we want to use the thermodynamical option ; we just have to change temp_hydro into ptn. Besides,this thermodynamical option needs to be more thoroughly tested because some preliminary tests have shown instabilities in particular temperature et soil moisture conditions (see #197).

• Thirdly, the parameter fr_dT which defines the freezing window in thermosoil_getdiff is not correctly used in hydrol_soil_coef. The present code:

               IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(njsc(ji))+min_sechiba))) THEN
                  ! Linear soil freezing or soil moisture below residual
                  IF (temp_hydro(ji, jsl).GE.273._r_std) THEN
                     x=1._r_std
                  ELSE IF (273._r_std.GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.271._r_std) THEN 
                     x=(temp_hydro(ji, jsl)-271._r_std)/fr_dT
                  ELSE 
                     x=0._r_std
                  ENDIF
               ELSE IF (ok_thermodynamical_freezing) THEN
                  ! Thermodynamical soil freezing
                  IF (temp_hydro(ji, jsl).GE.273._r_std) THEN
                     x=1._r_std
                  ELSE IF (273._r_std.GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.271._r_std) THEN 
                     x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
                          *((1000.*avan(njsc(ji))*(ZeroCelsius-temp_hydro(ji, jsl)) &
                          *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m))/(mc(ji,jsl, ins) &
                          -mcr(njsc(ji))),1._r_std)                
                  ELSE
                     x=0._r_std
                  ENDIF
               ENDIF
  

should be changed into:

             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(njsc(ji))+min_sechiba))) THEN
                ! Linear soil freezing or soil moisture below residual
                IF (temp_hydro(ji, jsl).GE. (ZeroCelsius+fr_dT/2.)) THEN
                   x=1._r_std
                ELSE IF ((ZeroCelsius+fr_dT/2.).GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.(ZeroCelsius-fr_dT/2.)) THEN 
                   x=(temp_hydro(ji, jsl)-(ZeroCelsius-fr_dT/2.))/fr_dT
                ELSE 
                   x=0._r_std
                ENDIF
             ELSE IF (ok_thermodynamical_freezing) THEN
                ! Thermodynamical soil freezing
                IF (temp_hydro(ji, jsl).GE.(ZeroCelsius+fr_dT/2.)) THEN
                   x=1._r_std
                ELSE IF ((ZeroCelsius+fr_dT/2.).GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.(ZeroCelsius-fr_dT/2.)) THEN 
                   x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
                        *((1000.*avan(njsc(ji))*(ZeroCelsius-temp_hydro(ji, jsl)) &
                        *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m))/(mc(ji,jsl, ins) &
                        -mcr(njsc(ji))),1._r_std)                
                ELSE
                   x=0._r_std
                ENDIF
             ENDIF

Catherine, Sarah and Fabienne.

[maignan]

Comment 6 committed in r3077.

[jgipsl]

Il y a un erreur sur la ligne 5540 dans rev [3077] car (ZeroCelsius-temp_hydro(ji, jsl)) peut devenir negative :

 ELSE IF ((ZeroCelsius+fr_dT/2.).GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.(ZeroCelsius-fr_dT/2.)) THEN
    x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
        *((1000.*avan(njsc(ji))*(ZeroCelsius-temp_hydro(ji, jsl)) &
        *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m))/(mc(ji,jsl, ins) &
        -mcr(njsc(ji))),1._r_std)      

Le problème se passe quand temp_hydro est dans l'intervalle [ZeroCelsius, ZeroCelsius+fr_dT/2.[, car on prend une puissance -m de (ZeroCelsius-temp_hydro(ji, jsl)), avec m qui est dans ]0,1[. Ca interdit les valeurs LE 0 pour (ZeroCelsius-temp_hydro(ji, jsl)).

Cela a été corrigé dans rev [3080] en changeant la température de référence pour ce calcul, puisque les bornes de l'intervalle ont changé:

  ELSE IF ((ZeroCelsius+fr_dT/2.).GT.temp_hydro(ji, jsl).AND.temp_hydro(ji, jsl).GE.(ZeroCelsius-fr_dT/2.)) THEN
     x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
          *((1000.*avan(njsc(ji))*(ZeroCelsius + fr_dT/2. - temp_hydro(ji, jsl)) &
          *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m))/(mc(ji,jsl, ins) &
          -mcr(njsc(ji))),1._r_std)

Agnès, Josefine

[fwang]

Fuxing is checking with Bertrand about the variables required by stomate module.

[aducharne]

Après vérifications nombreuses par Fuxing, Bertrand, et Agnès, on peut fermer ce ticket en s'appuyant sur les recommandations du document ​Ticket_195_end_FW2.doc

[jgipsl]

Done in trunk :

• [4631] : remove nbdl dimension which is only a duplication of nslm since new vertical discretization. This is ok also for Choisnel. No change in results.
• [4633] : Change diaglev to take the values at lower-layer interface (zlt) instead of values at nodes(znt). Change in results.

Still to be done:

• Change calculation of shumdiag, shumdiag_perma to be done without frac_hydro_diag. Remove subroutine hydrol_calculate_frac_hydro_diag and variable frac_hydro_diag which should not be needed anymore.
• Remove calculation of tmc_layh and send soilmoist instead to thermosoil

Agnès, Josefine

[aducharne]

The last two above points have been committed in [4637], so that we don't use frac_hydro_diag and subroutine hydrol_calculate_frac_hydro_diag anymore, and we don't calculate tmc_lay but transfer soilmoist to thermosoil.

In theory, these changes should not change the results. To check this, tests were run over 5 years with and without the commit. They reveal very slight changes in both the restart and history files. Logically, the most affected variables are related to thermal properties, snow, and some stomate output (cim, gpp, lai, etc.) but the changes are very very small and can likely be attributed to the effect of machine precision when doing differently the same operations.

Remaining actions to fully benefit from using the same vertical discretisation for heat and soil moisture are reported to new tickets:
#397 Can we get rid of temp_hydro and hydrol_calculate_temp_hydro in hydrol?
#398 We can remove the diaglevs when we get rid of Choisnel