= Martial notes for the merge =

== messages de TdO dans le module d'hydro ==
   * l1 
{{{
!tdo - enlever toute profondeur variable pour voir l'effet sur l'efficacite du code
}}}
   * l153
{{{
!!! A CHANGER DANS TOUT HYDROL: tmc_litter_res et sat ne devraient pas dependre de ji - tdo
}}}
   * hydrol_soil :
{{{
  SUBROUTINE hydrol_soil (kjpindex, dtradia, veget_max, soiltile, njsc, reinf_slope, &
       & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, returnflow, reinfiltration, irrigation, &
       & tot_melt, evap_bare_lim,  shumdiag, k_litt, litterhumdiag, humrel,vegstress, drysoil_frac)
    ! 
    ! interface description
    ! input scalar 
    INTEGER(i_std), INTENT(in)                               :: kjpindex         !!
    !To be removed later
}}}
     kjpindex doit-il être supprimé de cette fonction ?? ou est-ce un commentaire obsolète ?? à rechercher dans les versions antérieures.
 * deux commentaires dans la même fonction :
{{{
       !-
       !- step 6.5 : compute dr_ns with the bottom boundary condition 
       !- step 6.5 : initialize qflux at bottom of diffusion and avoid over saturated or under residual soil moisture 
       !-
}}}
   Le second à l'air obsolète ... à effacer ?
 * Toujours dans hydrol_soil, un commentaire important l3211 :
{{{
             !- Here we make the assumption that roots do not take water from the 1st layer. 
             !- Comment the us=0 if you want to change this.
}}}
  Il me semble utile de rajouter un getin pour cela, non ?
 * Concernant les erreurs de fermeture de bilan (cf variable waterbal_error), on indique plusieurs fois dans le code que l'on va l'arrêter :
{{{
             waterbal_error=.TRUE.
             CALL ipslerr(2, 'hydrol_split_soil', 'We will STOP after hydrol_split_soil.',&
                  & 'check_CWRR','PRECISOL SPLIT FALSE')
}}} 
   sans le faire ! Est-ce à rajouter ? 
== problèmes conversion (frac_bare/veget) versus (veget,veget_max) ==
 * dans hydrol_init : on passait veget (ie l'ancien veget_max) et on a pas changer le code par rapport
   à l'ancienne version. On obtient peut-être des incohérences : 
{{{
       IF ( MINVAL(resdist) .EQ.  MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
          resdist = veget
       ENDIF
       !
       !  Remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1. Thus we need vegtot
       !
       DO ji = 1, kjpindex
          vegtot(ji) = SUM(veget(ji,:))
       ENDDO
}}}
   Dois-je laisser le veget de la version standard ? D'après les discussions lors de la réunion du 20/10/2011, on met bien le veget_max dans ces définitions pour être cohérents avec le reste de hydrol.
 * dans hydrol_vegetupd on utilise vraiment frac_bare :
{{{
    DO jv = 1, nvm
       DO ji =1, kjpindex
          tot_frac_bare(ji) = tot_frac_bare(ji) + veget_max(ji,jv) * frac_bare(ji,jv)
       ENDDO
    END DO
}}}
   On doit donc le recalculer à partir de la formule dans slowproc :
{{{
    frac_bare(:,:) = zero
    frac_bare(:,1) = un
    IF (extcoef .LT. 100) THEN
       DO jv=2,nvm
          frac_bare(:,jv) = EXP(-extcoef * lai(:,jv))
       ENDDO
    ENDIF
}}}
   et la définition actuelle du veget :
{{{
    ! Ajout Nouveau calcul (stomate-like) 
    DO ji = 1, kjpindex
       SUMveg = 0.0
       veget(ji,1) = veget_max(ji,1)
       DO jv = 2, nvm
          veget(ji,jv) = veget_max(ji,jv) * ( 1. - exp( - lai(ji,jv) * ext_coef(jv) ) )
          veget(ji,1) = veget(ji,1) + (veget_max(ji,jv) - veget(ji,jv))
       ENDDO
 [...]
}}}
>> Commentaire de Jan : [[BR]]
>> Ne serait il pas mieux d'avoir des fonctions définies dans Slowproc qui donnent ces informations. Des genres de "slowproc_query veget". [[BR]]
>> On va s'amuser à répliquer des formules avec tous les risques d'erreures car on a oublié de changer dans une partie du code. [[BR]]
> Pour l'instant (merge de hydrol et de routing), il n'y a qu'à un seul endroit ou l'on utilise frac_bare. Le veget et le veget_max étant passé en paramètre de toutes les fonctions 
> de enerbil, condveg et diffuco qui pourraient l'utiliser, il est très de réécrire la formule suivante.

   On déduit alors en combinant ces deux définitions :
{{{
    frac_bare(:,1) = un
    frac_bare(:,2:nvm) = undef_sechiba
    DO jv = 2, nvm
       DO ji = 1, kjpindex
          IF ( veget_max(ji,jv) .GT. min_sechiba ) &
               frac_bare(ji,jv) = 1 - veget(ji,jv) / veget_max(ji,jv)
       ENDDO
    ENDDO
}}}
   L'utilisation du undef_sechiba me paraît cohérente lorsque l'on a pas de présence du PFT sur le pixel.[[BR]]
   On avait d'ailleurs une erreur avant avec le frac_bare car il était parfois définit (utilisé ?) lorsque veget_max(ji,jv) était nul. [[BR]]
   Un problème surgit alors lorsque l'on définit :
{{{
    tot_frac_bare(:) = zero
[...]
    DO jv = 1, nvm
       DO ji =1, kjpindex
          tot_frac_bare(ji) = tot_frac_bare(ji) + veget_max(ji,jv) * frac_bare(ji,jv)
       ENDDO
    ENDDO
}}}
   car on ne teste pas veget_max > min_sechiba. Je le rajoute.
 * '''hydrol_canop''' : je me pose la question du '''veget/veget_max''' ici car on définit '''qsintveg et precisol''' dans cette routine. 
   Si l'on passe veget_max, on ne tient plus compte du LAI pour calculer ces deux variables importantes. Je laisse donc veget.
   C'était bien un bogue de la version LMD. [[BR]]
   Voici les corrections que l'on a appliqué pour rendre le modèle conservatif :
{{{
    precisol(:,1)=veget_max(:,1)*precip_rain(:)
    DO jv=2,nvm
       DO ji = 1, kjpindex
          zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv)) 
          IF ( ok_throughfall_by_pft ) THEN
             ! correction throughfall Nathalie - Juin 2006
             precisol(ji,jv) = (veget(ji,jv)*throughfall_by_pft(jv)*precip_rain(ji)) + &
                  qsintveg(ji,jv) - zqsintvegnew (ji,jv) + &
                  (veget_max(ji,jv) - veget(ji,jv))*precip_rain(ji)
          ELSE
             precisol(ji,jv) = qsintveg(ji,jv) - zqsintvegnew (ji,jv) + &
                  (veget_max(ji,jv) - veget(ji,jv))*precip_rain(ji)
          ENDIF
       ENDDO
    END DO
[...]
    DO jv=2,nvm
       qsintveg(:,jv) = zqsintvegnew (:,jv)
    END DO
}}}
   On ne fait donc plus aucun calcul de l'interception pour le PFT 1. On prends bien en compte le precip_rain sur les surface fracnobio dans hydrol_snow.

== Formules de convolution/déconvolution ==
La convolution est définit par les corr_veg_soil ("(:,nvm,nstm) percentage of each veg. type on each soil of each grid point")
et les cvs_over_veg ("(:,nvm,nstm) old value of corr_veg_soil/veget_max kept from diag to next split" : commentaire sûrement faux d'après la définition).
{{{
       ! somme(corr_veg_soil / vegtot / veget_max ) = 1 - vegtot = frac_nobio
       DO ji=1,kjpindex
          tot_corr_veg_soil(ji)=zero
          DO jst = 1, nstm
             DO jv = 1,nvm
                tot_corr_veg_soil(ji)=tot_corr_veg_soil(ji)+corr_veg_soil(ji,jv,jst) / vegtot(ji) / veget_max(ji,jv)
             ENDDO
          ENDDO
       ENDDO

       DO ji=1,kjpindex
          IF ( ABS( tot_corr_veg_soil(ji) - 2 ) .GT. EPS1 ) THEN
             WRITE(numout,*) 'corr_veg_soil SPLIT FALSE:ji=',ji,tot_corr_veg_soil(ji)
             WRITE(numout,*) 'err',ABS( tot_corr_veg_soil(ji) - 2 )
             WRITE(numout,*) 'vegtot',vegtot(ji)
             DO jv=1,nvm
                WRITE(numout,*) 'jv,veget_max,corr_veg_soil',jv,veget_max(ji,jv),corr_veg_soil(ji,jv,:)
             END DO
          ENDIF
       ENDDO
}}}
et 
{{{
    cvs_over_veg(:,:,:) = zero
    DO jv=1,nvm
       DO ji=1,kjpindex
          IF(veget_max(ji,jv).GT.min_sechiba) THEN
             DO jst=1,nstm
                cvs_over_veg(ji,jv,jst) = corr_veg_soil(ji,jv,jst)/vegtot(ji) / veget_max(ji,jv)
             ENDDO
          ENDIF
       END DO
    END DO
}}}

Donc pour chaque point de terre :
{{{
#!formula
#nomode
#density=200
\[
\begin{array}{lcl}
\displaystyle \sum_{jv=1,nvm} \sum_{jst=1,nstm}\frac {corr\_veg\_soil(jv,jst)}{vegtot \times veget\_max(jv)} &=& \\
\displaystyle \sum_{jv=1,nvm} \sum_{jst=1,nstm} cvs\_over\_veg(jv,jst)  &=& 2
\end{array}
\]
}}}
soit 
{{{
#!formula
#nomode
#density=200
\[
\begin{array}{lcl}
\displaystyle \sum_{jv=1,nvm} \sum_{jst=1,nstm}\frac {corr\_veg\_soil(jv,jst)}{veget\_max(jv)} &=& 2\,vegtot \\
&=& 2 \times vegtot 
\end{array}
\]
}}}
Ce qui semble assez étrange, 

La formule de convolution est alors :
{{{
    !
    !
    ! split 2d variables into 3d variables, per soil type
    !
    precisol_ns(:,:)=zero
    DO jv=1,nvm
       DO jst=1,nstm
          DO ji=1,kjpindex
             IF(veget_max(ji,jv).GT.min_sechiba) THEN
                precisol_ns(ji,jst)=precisol_ns(ji,jst)+precisol(ji,jv)* &
                     & corr_veg_soil(ji,jv,jst) / vegtot(ji) / veget_max(ji,jv))
             ENDIF
          END DO
       END DO
    END DO
}}}
et la fomule de déconvolution est donnée par le test :
{{{
    !
    ! Now we check if the deconvolution is correct and conserves the fluxes:

    IF (check_cwrr) THEN


       tmp_check1(:)=zero
       tmp_check2(:)=zero  

       ! First we check the precisol and evapnu

       DO jst=1,nstm
          DO ji=1,kjpindex
             tmp_check1(ji)=tmp_check1(ji) + &
                  & precisol_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji)
          END DO
       END DO

       DO jv=1,nvm
          DO ji=1,kjpindex
             tmp_check2(ji)=tmp_check2(ji) + precisol(ji,jv)
          END DO
       END DO

       DO ji=1,kjpindex   

          IF(ABS(tmp_check1(ji)- tmp_check2(ji)).GT.allowed_err) THEN
             WRITE(numout,*) 'PRECISOL SPLIT FALSE:ji=',ji,tmp_check1(ji),tmp_check2(ji)
             WRITE(numout,*) 'err',ABS(tmp_check1(ji)- tmp_check2(ji))
             WRITE(numout,*) 'vegtot',vegtot(ji)

             DO jv=1,nvm
                WRITE(numout,*) 'jv,veget_max, precisol',jv,veget_max(ji,jv),precisol(ji,jv)
                DO jst=1,nstm
                   WRITE(numout,*) 'corr_veg_soil:jst',jst,corr_veg_soil(ji,jv,jst)
                END DO
             END DO

             DO jst=1,nstm
                WRITE(numout,*) 'jst,precisol_ns',jst,precisol_ns(ji,jst)
                WRITE(numout,*) 'soiltile', soiltile(ji,jst)
             END DO
             waterbal_error=.TRUE.
             CALL ipslerr(2, 'hydrol_split_soil', 'We will STOP after hydrol_split_soil.',&
                  & 'check_CWRR','PRECISOL SPLIT FALSE')
          ENDIF

       END DO

    ENDIF
}}}
On a aussi la formule d'intégration de humrel sur les tile qui semble un peu différente :
{{{
    DO jst=1,nstm
       DO jv=1,nvm
          DO ji=1,kjpindex
             humrel(ji,jv)=humrel(ji,jv)+humrelv(ji,jv,jst)*soiltile(ji,jst) &
                  & * cvs_over_veg(ji,jv,jst)*vegtot(ji)
             humrel(ji,jv)=MAX(humrel(ji,jv),zero)
          END DO
       END DO
    END DO
}}}

=== Première version conservative ===
Conservative version of new hydrology can be found and downloaded with
revision 606 of branches/Hydrology :
{{{
login> recup_my_ORCHIDEE your.login branches/Hydrology Martial.Mancip@ipsl.jussieu.fr 606
}}}


Your would be very kind to begin some tests with this first version and
give us some backup after compare your old reference versions.
There must be some changes ! we hope to understand all of them ;-)

=== Variables normalisées ===

La normalisation par vegtot (donc en mailles totales, sans tenir compte des surfaces en nobio) de certaines variables semble incorrecte:
 1. humrel et vegstress :
   a. dans hydrol_init :
{{{
                vegstress(ji,jv)=vegstress(ji,jv) + vegstressv(ji,jv,jst) * &
                     & soiltile(ji,jst) * cvs_over_veg(ji,jv,jst) * vegtot(ji)

                humrel(ji,jv)=humrel(ji,jv) + humrelv(ji,jv,jst) * & 
                     & soiltile(ji,jst) * cvs_over_veg(ji,jv,jst) * vegtot(ji)
                humrel(ji,jv)=MAX(humrel(ji,jv), zero)* mask_veget(ji,jv)           
}}}
   b. dans hydrol_diag_soil :
{{{
                vegstress(ji,jv)=vegstress(ji,jv)+vegstressv(ji,jv,jst)*soiltile(ji,jst) &
                     & * corr_veg_soil(ji,jv,jst) *vegtot(ji)/veget_max(ji,jv)
                vegstress(ji,jv)= MAX(vegstress(ji,jv),zero)
[...]
                 cvs_over_veg(ji,jv,jst) = corr_veg_soil(ji,jv,jst)/vegtot(ji) / veget_max(ji,jv)
[...]
             humrel(ji,jv)=humrel(ji,jv)+humrelv(ji,jv,jst)*soiltile(ji,jst) &
                  & * cvs_over_veg(ji,jv,jst)*vegtot(ji)
             humrel(ji,jv)=MAX(humrel(ji,jv),zero)
}}}
   Ce qui rend les définitions de vegstress incompatibles ! car on a pas : 
{{{
      cvs_over_veg(ji,jv,jst) * vegtot(ji) != corr_veg_soil(ji,jv,jst) *vegtot(ji)/veget_max(ji,jv)
}}}
   En fait, d'après Jan, on devrait toujours avoir vegstress = humrel dans ce code.
 1. définition à préciser de frac_bare_ns dans hydrol_vegupd
{{{
                frac_bare_ns(ji,jst) = frac_bare_ns(ji,jst) + corr_veg_soil(ji,jv,jst) * frac_bare(ji,jv) / vegtot(ji)
}}}
    et de precisol_ns dans hydrol_split_soil 
{{{
                precisol_ns(ji,jst)=precisol_ns(ji,jst)+precisol(ji,jv)* &
                     & corr_veg_soil(ji,jv,jst) /vegtot(ji) / veget_max(ji,jv)
}}}
    Pourquoi normaliser ces variables ?

== hydrol_tmc_update : LAND_USE or DGVM == 
The possibility to activate LAND_USE change or DGVM with STOMATE with the new Hydrology is now avaible.

This updates of veget_max are daily (for DGVM) or yearly (for LAND USE scheme) and conservation of tmc (humtot) and mc profile had to be developped. To give right new tmc, I have used the same algorithm as in stomate_lccchange that conserve carbon in stomate with LAND_USE. And the mc profile is updated in consequence. The main problem was to "create" a new soil (for example when natural grass fields are replaced by bare soil).

In [http://dods.ipsl.jussieu.fr/orchidee/HYDROLOGY/MONITORING_HYDRO2-CWRR-SECSTO-test-DGVM2.tgz test Hydro DGVM], you will find a MONITORING directory with my test run :
 * DGVM activated with STOMATE
 * IMPOSE_VEG as initial condition (with only PFT 10 at start)
 * global WG_cru 4° as forcing file
 * only 5 days
You may see the maxvegetfrac modifications (for PFT 1, 9, 10, 11) and humtot conservation (no jump after end of day). 

This version of new hydrology can be found and downloaded with
revision *610* of branches/Hydrology :
{{{
login> recup_my_ORCHIDEE your.login branches/Hydrology Martial.Mancip@ipsl.jussieu.fr 610
}}}

== traitement des corrections de Nathalie ==
=== Gestion du throughfall_by_pft ===
J'ai rajouté un booléen ok_throughfall_by_pft.
{{{
       IF ( ok_throughfall_by_pft ) THEN
          ! Correction Nathalie - Juin 2006 - une partie de la pluie arrivera toujours sur le sol
          ! sorte de throughfall supplementaire
          qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
       ELSE
          qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * precip_rain(:)
       ENDIF
}}}
pour toutes ces corrections (dans hydrol_canop).

J'ai globalisé le OFF_LINE_MODE qui depuis intersurf indiquait si l'on utilisait une interface "intersurf_main*" (utilisation off-line du modèle)
ou bien "intersurf_gathered*" (pour les utilisations on-line avec le GCM).

Comme convenu, j'ai donc imposé comme valeur par défaut de ce ok_throughfall_by_pft :
{{{
       ok_throughfall_by_pft = .FALSE.
       IF ( .NOT. OFF_LINE_MODE ) ok_throughfall_by_pft = .TRUE.
}}}  
== merge de routing ==
=== Messages de TdO dans le module de routage ===
Dans la procédure routing_truncate, un message et une modification de Tristan a attiré mon attention :
{{{
       ! tdo - To take into account rivers that do not flow to the oceans 
       IF ( route_tobasin(ff(1), ff(2)) .GT. nbasmax ) THEN
!       IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN
}}}
est-ce correcte ?

=== notes sur le merge ===
... et sur la version précédente !

 * beaucoup de modification. Ajout de la sauvegarde d'un fichier de diagnostique, soit en netcdf90 si le paramètre RIVER_DESC_FILE contient une extension ".nc" dans routing_diagncfile, soit en texte. 
   Il faudra reconstruire les appels nf90 pour ce fichier avec fliocom.f90.
 * des INTENT(in/inout/out) ont été rajouté par rapport aux versions pour diminuer les messages de warning des compilateurs. Il faudra vérifier leur validité avec lf95.
 * dans routing_simplify, todosz est bien initialisé, mais jamais utilisé... cela semble être un oubli : 
{{{
                todosz(itodo) = outsz(ip)
}}}
 * dans routing_irrigmap, deux tableaux diagnostiques, floodmap et irrigmap ne sont jamais testés (ni sauvés car ce sont directement floodplains et irrigated qui le sont).
 * Lors de la modification du fichier de paramètres sechiba.def, je viens de m'apercevoir d'un '''bogue dans la définition des paramètres de vitesse de routage''' de la nouvelle version :
{{{
  !  The time constants are in days.
  ! 
  REAL(r_std), PARAMETER                            :: slow_tcst_cwrr = 3.0 
  REAL(r_std), PARAMETER                            :: fast_tcst_cwrr = 3.0 
  REAL(r_std), PARAMETER                            :: stream_tcst_cwrr = 0.24 
  REAL(r_std), PARAMETER                            :: flood_tcst_cwrr = 4.0 
  REAL(r_std), PARAMETER                            :: swamp_cst_cwrr = 0.2
[...]
    !
    IF ( hydrol_cwrr ) THEN
       fast_tcst = slow_tcst_cwrr
    ELSE
       fast_tcst = fast_tcst_chois
    ENDIF
}}} 
   Il parait peut probable que la vitesse d'écoulement du répertoire rapide soit la même que la vitesse d'écoulement 
   du répertoire lent (3jours) alors que la vitesse d'écoulement du réservoir stream est de 0.24 jours.

> Commentaire de Matthieu G.:[[BR]]
> Ceci n'est pas un bug dans le code. Quand on se réfère à la thèse de Tristan (en haut de la page 105), on trouve la justification de mettre la constante du réservoir slow égale à celle du fast quand on utilise l'hydrologie 11 couches.[[BR]]
> Donc, ce n'est pas un bug mais en revanche, la justification de Tristan peut être rediscutée d'un point de vue scientifique...

   Je corrige ce bogue en remettant la valeur par défaut de choisnel de 25 j pour slow_tcst_cwrr 
   et en corrigeant la définition de fast_tcst qui doit valoir fast_tcst_cwrr dans le cas hydrol_cwrr.
== merge de sechiba ==
 * A la fin de l'appel de hydrol_main, on a dans la version LMD :
{{{
          rsol(:) = -un
}}}
   Pourquoi cette initialisation a disparut dans la version standard.
 * A revoir : on a encore des ntsm ?? ça ne devrait pas être nscm, comme dans hydrol ?
{{{
       IF ( control_in%hydrol_cwrr ) THEN
          CALL histwrite(hist_id, 'soiltile',  kjit, soiltile, kjpindex*nstm, indexsoil)
}}}

== merge de intersurf ==
 * nbdl devient nslm dans la définition de diaglev.
 * On définit diaglev dans intsurf_history :
{{{
    DO jv = 1, nslm-1
       diaglev(jv) = dpu_max/(2**(nslm-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / deux
    ENDDO
    diaglev(nslm) = dpu_max
    !
}}}
  On aurait pu le définir dans intsurf_config, non ?

== merge de slowproc ==
 * à vérifier avec les autres driver : déplacement/redéfinition de diaglev + validation stomate activé
 * ancien bug toujours présent dans la version LMD dans la définition des soiltiles :
{{{
          ! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
          ! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
          ! in the other modules to perform separated energy balances
          njsc(:) = 0
          soiltile(:,:) = zero
          DO ji = 1, kjpindex
             njsc(ji) = MAXLOC(soilclass(ji,:),1)
             soiltile(:,1) = SUM(frac_nobio(ji,:))
          ENDDO
          DO jv = 1, nvm
             jst = pref_soil_veg(jv)
             DO ji = 1, kjpindex
                soiltile(ji,jst) = soiltile(ji,jst) + veget(ji,jv)
             ENDDO
          ENDDO
}}}
   On prend alors ici toujours le frac_nobio du dernier point. Le code corrigé est :
{{{
          ! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
          ! on a changing vegetation (but then some adaptation should be made to hydrol) and be also used
          ! in the other modules to perform separated energy balances
          njsc(:) = 0
          soiltile(:,:) = zero
          DO ji = 1, kjpindex
             njsc(ji) = MAXLOC(soilclass(ji,:),1)
             soiltile(ji,1) = SUM(frac_nobio(ji,:))
          ENDDO
          DO jv = 1, nvm
             jst = pref_soil_veg(jv)
             DO ji = 1, kjpindex
                soiltile(ji,jst) = soiltile(ji,jst) + veget(ji,jv)
             ENDDO
          ENDDO
}}}
 * Définitions des soilclass, clayfraction, soiltile et njsc : On fait deux appel à slowproc_soilt
{{{
       IF ( MINVAL(soilclass) .EQ. MAXVAL(soilclass) .AND. MAXVAL(soilclass) .EQ. val_exp .OR.&
            & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN

          CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)

       ENDIF

       IF ( MINVAL(soiltile) .EQ. MAXVAL(soiltile) .AND. MAXVAL(soiltile) .EQ. val_exp .OR.&
            & MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int .OR.&
            & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN
          CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)

          ! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
}}}
  Le second appel me semble inutile. Si les soiltile et njsc ne sont pas dans le restart, il suffit simplement de les recalculer. Je simplifie donc :
{{{
       IF ( MINVAL(soilclass) .EQ. MAXVAL(soilclass) .AND. MAXVAL(soilclass) .EQ. val_exp .OR.&
            & MINVAL(clayfraction) .EQ. MAXVAL(clayfraction) .AND. MAXVAL(clayfraction) .EQ. val_exp) THEN

          CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)

       ENDIF

       IF ( MINVAL(soiltile) .EQ. MAXVAL(soiltile) .AND. MAXVAL(soiltile) .EQ. val_exp .OR.&
            & MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int) THEN

          ! Soiltiles are only used in hydrol, but we fix them in here because some time it might depend
}}}
 * définition de ext_coef dans slowproc_lai : incohérence manifeste avec l'utilisation de stomate dans la version LMD. [[BR]]
   Lorsque active stomate avec la version LMD, on n'utilise slowproc_lai uniquement si on a pas de restart (dans slowproc_init).
   Or on définit un extcoef (constant pour tous les PFTs) dans cette routine par un getin_p qui n'est pas a prioiri cohérent avec
   le '''ext_coef(j par PFT)''' utilisé dans diffuco et slowproc_veget et définit dans constantes_co2
{{{
! extinction coefficient of the Monsi&Seaki relationship (1953)
 &  ext_coef = (/.5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5, .5/)
}}}
   '''et''' le '''ext_coeff(j par PFT)''' utilisé dans lpj_establish et lpj_light [[BR]]
   (définit dans stomate_constants, mais qui est initialisé avec le ext_coef de constantes_co2 dans stomate_data)
{{{
! extinction coefficient of the Monsi&Seaki (53) relationship
  REAL(r_std), SAVE, DIMENSION(nvm)                     :: ext_coeff
}}}
   Comme convenu, on ne définira plus le frac_bare dans la routine, mais on utilisera partout le ext_coef de constantes_co2 (paramétré par un getin dans slowproc_init) pour simplifié. 

> Commentaire de Matthieu G.:[[BR]]
> 2 remarques au sujet de ce coeff d'extinction:[[BR]]
>   * Il y a une grande incertitude sur sa valeur qui nécessitera, à mon avis, des tests de sensibilité en temps voulu.[[BR]]
>       - ext_coeff=0.5 par défaut dans le modèle[[BR]]
>       - ext_coeff=3.0 d'après Perrier & Guimberteau (on avait regardé l'évaporation du sol nu principalement sur les USA)[[BR]]
>       - ext_coeff=4.0 d'après la thèse de D'Orgeval (simulation CWRR2 en bas de la page 87). En haut de la page 92, ce coefficient aurait été fixé (donc =4) "à une valeur élevée principalement pour éviter une réactivité trop forte via l'évaporation dès que le LAI diminue et éviter des incohérences avec les mesures d'Hapex-Sahel". Il préconise même de faire d'autres tests pour mieux caler ce paramètre.[[BR]]
>   * Actuellement, le coef ne dépend pas du PFT. L'attribution de différentes valeurs de extcoef pour les différentes PFTs est par conséquent souhaitable.

 * Dans slowproc_lai de la version tagguée, le commentaire "Test Nathalie. On impose LAI PFT 1 a 0" est obsolète car on doit toujours bouclé sur 2,nvm pour le calcul du lai. Le lai du PFT 1 vaut
   toujours zéro par convention. 
 * Enfin dans le cas read_lai=.TRUE., on conserve (pour l'instant ?) la suppression de l'interpolation de la carte de lai lorsque type_of_lai vaut "mean" pour un PFT du fait du bogue de neige dans l'ancienne carte de lai.
   La nouvelle carte de LAI corrige ce problème et on a donc type_of_lai toujours égale à "inter" dans les versions suivants la orchidee_1_9_5.
 * il me semble que le commentaire de la routine slowproc_soilt devrait être revu :
{{{
  SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction)
    !
    !
    !   This subroutine should read the Zobler map and interpolate to the model grid. The method
    !   is to get fraction of the three main soiltypes for each grid box.
    !   The soil fraction are going to be put into the array soiltype in the following order :
    !   coarse, medium and fine.
    !
}}}

>
> soiltile = fraction of PFT on each soil-hydrology tile [[BR]]
>
> njsc = indexing of PFT to soiltile
>

== Notes sur les tests ==

 1. Choisnel sans routage
 1. CWRR sans routage ni irrigation
 1. CWRR avec routage mais sans irrigation ou plaines d'inondation
 1. CWRR avec routage et irrigation+plaines
 1. Choisnel avec routage
 1. Choisnel avec routage et irrigation. 
 1. Activation de stomate avec la même liste !
    1. Choisnel sans routage
    1. CWRR sans routage ni irrigation
    1. CWRR avec routage mais sans irrigation ou plaines d'inondation
    1. CWRR avec routage et irrigation+plaines
    1. Choisnel avec routage
    1. Choisnel avec routage et irrigation. 

== merge de thermosoil ==
voir le paragraphe suivant [wiki:Branches/MergeHydro/Martial_notes_on_merge#NotesdeIsabelleGoutevinsurlaversionhydro].

On a modifier l'appel et le code de thermosoil_humlev pour tenir compte de l'épaisseur de neige dans le calcul du wetdiag. Voir le commit [606].

= Notes de Isabelle Goutevin sur la version hydro =

Explications dans le fichier [attachment:qques_bugs.odt] : 
 1. prise en compte de la hauteur de neige dans le shumdiag utilisé par thermosoil.
 1. indices de l'échelle diag (nbdl) vs indices de l'échelle hydro (nslm).
Les fichiers [attachment:correction_neige] et [attachment:correction_indices] contiennent les corrections proposées.