Pb in calculation of pcapa, pcapa_en and pkappa

[jgipsl]

In the trunk today rev[2247] and since long-time (the same is found in tag 1_9_5), the variables pcapa, pcapa_en and pkappa are calculated only for the levels 1:ngrnd-2. Why? See here ​thermosoil rev 2247.

See here the corresponding lines for the calculation. The calculations are now in subroutine thermosoil_getdiff_old_thermix_with_snow since r[222] called from thermsoil_coef. Before the same lines where performed directly in thermosoil_coef.

    ! Computation of the soil thermal properties; snow properties are also accounted for

    DO ji = 1,kjpindex
      snow_h(ji) = snow(ji) / sn_dens

      IF ( snow_h(ji) .GT. zz_coef(1) ) THEN
          pcapa(ji,1) = sn_capa
          pcapa_en(ji,1) = sn_capa
          pkappa(ji,1) = sn_cond
      ELSE IF ( snow_h(ji) .GT. zero ) THEN
          pcapa_en(ji,1) = sn_capa
          zx1(ji) = snow_h(ji) / zz_coef(1)
          zx2(ji) = ( zz_coef(1) - snow_h(ji)) / zz_coef(1)
          pcapa(ji,1) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
          pkappa(ji,1) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
      ELSE
          pcapa(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
          pkappa(ji,1) = so_cond_dry + shum_ngrnd_perma(ji,1)*(so_cond_wet - so_cond_dry)
          pcapa_en(ji,1) = so_capa_dry + shum_ngrnd_perma(ji,1)*(so_capa_wet - so_capa_dry)
      ENDIF
      !
      DO jg = 2, ngrnd - 2
        IF ( snow_h(ji) .GT. zz_coef(jg) ) THEN
            pcapa(ji,jg) = sn_capa
            pkappa(ji,jg) = sn_cond
            pcapa_en(ji,jg) = sn_capa
        ELSE IF ( snow_h(ji) .GT. zz_coef(jg-1) ) THEN
            zx1(ji) = (snow_h(ji) - zz_coef(jg-1)) / (zz_coef(jg) - zz_coef(jg-1))
            zx2(ji) = ( zz_coef(jg) - snow_h(ji)) / (zz_coef(jg) - zz_coef(jg-1))
            pcapa(ji,jg) = zx1(ji) * sn_capa + zx2(ji) * so_capa_wet
            pkappa(ji,jg) = un / ( zx1(ji) / sn_cond + zx2(ji) / so_cond_wet )
            pcapa_en(ji,jg) = sn_capa
        ELSE
            pcapa(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
            pkappa(ji,jg) = so_cond_dry + shum_ngrnd_perma(ji,jg)*(so_cond_wet - so_cond_dry)
            pcapa_en(ji,jg) = so_capa_dry + shum_ngrnd_perma(ji,jg)*(so_capa_wet - so_capa_dry)
        ENDIF
      ENDDO
    
    ENDDO ! DO ji = 1,kjpindex

Note also that in the initialization part, in thermosoil_var_init, the same variables are calculated differently for all 1:ngrnd levels. This is the reason why the code don't stop with error of non initialization.

Proposed solution :

I don't see the reason why the second DO loop can not go up to 2:ngrnd. zz_coef is calcalated for all 1:ngrnd levels.

[jpolcher]

I think there is a need to simplify and clean-up the code here in order to eliminate this and probably as well other bugs.

• thermosoil_coef is the reference for computing pkappa, pcapa, cgrnd and dgrnd and should be the only routines where they are defined.

• The restart procedure should only archive and retrieve the variables needed as input into thermosoil_coef (ptn, ...).

• Once these variables are restored in the code then thermosoil_coef is called and the model can start.

I am convinced that among all the variables in the restart today, some are of the old timestep and thus destroying the 1+1=2 test.

[jgipsl]

I agree with all you said. I have integrated this in a local working version. There are several tickets opened for this purpose, see #92, #143. And that's why this ticket needs to be solved.

[jgipsl]

Jan Polcher said, 2 oct 2014:
The last 2 layers must not be filled with snow. Do the following:

DO jg = ngrnd - 1, ngrnd
       pcapa(ji,jg) = so_capa_dry
       pkappa(ji,jg) = so_cond_dry
       pcapa_en(ji,jg) = so_capa_dry
END DO

[jgipsl]

Done in trunk rev [2405]