Pb when using 12 USDA textures with ok_freeze

[aducharne]

Reported by Jinfeng Chang on 28/11/2105:
I am doing global simulation with ORCHIDEE trunk rev3006, which is close to current reference run rev3013.
The simulation goes well with USDA soil classification map and without activating explicit snow and soil freezing.
However, when I try to activate the explicit snow and soil freezing and use USDA soil, the model keep running but gives empty values for some grids over the world, while the exactly the same setting except using FAO soil map instead, the model works well and give reasonable value.
Over 62482 grids over the world, 14 grids of type 8, 10 grids of type 11, and 1408 grids of type 12 shows the pb (undef values for runoff and soil moisture, but also other variables related to soil, such as HETERO_RESP.)

Proposed solution by Agnès Ducharne, 11/12/2015:
I've been focusing on the part of hydrol which calculates k_lin and d_lin, in hydrol_var_init.
To this end, I extracted the corresponding code for some tests.
With jsc=12 (clay), you get a lot of k_lin = 0 and d_lin=0, for low mc values.
Same thing with jsc =8 to 11
And when you activate freeze, you reduce the effective liquid moisture, which may explain why this revealed serious pbs.

I tested the removal of the modifications of avan and nvan introduced by Tristan d'Orgeval, which do not make much sense to me. These modifs are explained in section 3.3 in the note on hydrol.f90, where you can also find infos regarding the 12 USDA texture classes:
​http://forge.ipsl.jussieu.fr/orchidee/raw-attachment/wiki/Documentation/UserGuide/eqs_hydrol.pdf
To cancel these changes, you need to set n0 = 0., nk_rel = 0., a0 = 0., ak_rel = 0.
It does improve the values (less zero k_lin), but it's not enough to always get strictly positive k_lin and increasing d_lin with mc_lin, which is necessary for the numerical scheme dealing with water diffusion to work (thanks to Pascal Maugis for this information.

The modified code below (in hydrol_var_init) aims at satisfying the above constraint. It gives gives very small but positive k_lin (down to E-64 mm/d) at the bottom of the soil for textures 11 and 12, with the n0, nk_rel etc. of Tristan. And when you cancel the changes of Tristan, you get pretty decent values.

  DO jsl = 1, nslm
     nvan_mod = n0 + (nvan(jsc)-n0) * nfact(jsl,jsc)
     avan_mod = a0 + (avan(jsc)-a0) * afact(jsl,jsc)
     m = un - un / nvan_mod
     mc_lin(imin,jsc)=mcr(jsc)
     mc_lin(imax,jsc)=mcs(jsc)
     DO ji= imin+1, imax-1 
        mc_lin(ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin)
        frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
        k_lin(ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
     ENDDO

! we track jiref, the bin under which mc is too small and we may get zero k_lin     
     ji=imax-1
     DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0))
        !write(*,*) k_lin(ji,jsl,jsc)-1.e-32
        jiref=ji
        ji=ji-1
     ENDDO
     !write(*,*) 'jsl,jiref=',jsl,jiref
     DO ji=jiref-1,imin,-1
        k_lin(ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10.
     ENDDO    
     
     DO ji = imin,imax-1 
        a_lin(ji,jsl,jsc) = (k_lin(ji+1,jsl,jsc)-k_lin(ji,jsl,jsc)) / (mc_lin(ji+1,jsc)-mc_lin(ji,jsc))
        b_lin(ji,jsl,jsc)  = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc)
        
!        IF(ji.NE.imin.AND.ji.NE.imax-1) THEN
         IF(ji.NE.imax-1) THEN
              
           frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
           d_lin(ji,jsl,jsc) =MAX(1.e-32,(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) *  &
                &  ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) * &
                &  (  frac**(-un/m) -un ) ** (-m))
           
           frac=MIN(un,(mc_lin(ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
           d_lin(ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*&
                &  ( (frac**(-un/m))/(mc_lin(ji+1,jsc)-mcr(jsc)) ) * &
                &  (  frac**(-un/m) -un ) ** (-m)
           
           d_lin(ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc))    
                     
!        ELSEIF(ji.EQ.imin) THEN
!           d_lin(ji,jsl,jsc) = zero ! we don't want zero diffusivities
        ELSEIF(ji.EQ.imax-1) THEN
           frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
           d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &
                & ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) *  &
                & (  frac**(-un/m) -un ) ** (-m)         
        ENDIF
     ENDDO

! we adjust d_lin where k_lin was previously adjusted, or we get non-monotonous variations
     DO ji=jiref-1,imin,-1
        d_lin(ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10.
     ENDDO
        
  ENDDO

Verification in simulations:
Jinfeng performed 5 simulations A1-A5 with 2 degree CRU-NCEP forcing over 1961-2012 (results were analysed over 2001-2010). Soil freezing and explicit snow were activated in all simulations:
A1: USDA + OLD + NVAN
A2: USDA + NEW + NVAN
A3: USDA + NEW + NONVAN
A4: Zobler + OLD + NVAN
A5: Zobler + NEW + NVAN
USDA=USDA soil map; Zobler=Zobler soil map; OLD=original hydrol_var_init; NEW=new hydrol_var_init with k_lin, d_lin modified; NVAN=with Tristan's changes of avan/nvan; NONVAN=cancelling Tristan's changes of avan/nvan

The output are available on ​http://dods.extra.cea.fr/work/p529chan/OL2/PROD/rev3006cncep/isi.cncep.HIS.3006.snow.A2
You can also find it on CURIE: /ccc/store/cont003/dsm/p529chan/IGCM_OUT/OL2/PROD/rev3006cncep/
The folders 'isi.cncep.HIS.3006.snow.A1’ - 'isi.cncep.HIS.3006.snow.A5’ contain all the output files.

Main results:

• the correction works (no more undef values)
• evaporation is larger with FAO (3 classes, sim A5) than with USDA (12 textures, simA2), and cancelling the changes in avan/nvan (sim A3) further reduces evap; with logical related changes of drainage and soil moisture (which increase when evap decrease)

[jgipsl]

Correction done in the trunk rev [3082]. Note that some small differences where done compared to what is set above.

Agnès, Josefine