Changeset 6954


Ignore:
Timestamp:
12/03/20 18:16:43 (10 months ago)
Author:
agnes.ducharne
Message:

Introducing the changes of branch SP-MIP r6950.
These developments were made by Salma Tafasca durinh her PhD and allow ORCHIDEE to:
(1) read soil hydraulic parameters instead of soil texture maps (with dimension changes) or impose a uniform soil texture all over land, using keywords SPMIPEXP and EXP4;
(2) cancel the effect of roots on Ks with new keyword KFACT_ROOT_TYPE;
(3) introduce a 13th soil texture class for clay oxisols in addition to the 12 USDA classes (with thermal parameters equal to the ones of clay).
These changes to branch 2.2 have been tested in Salma's PhD simulations, and are fully compatible with former ways to control soil texture and hydraulic parameters.

Location:
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90

    r6189 r6954  
    1919!!                 This routine was originaly developed by Patricia deRosnay. 
    2020!! 
    21 !! RECENT CHANGE(S) : None 
    22 !! 
     21!! RECENT CHANGE(S) : November 2020: It is possible to define soil hydraulic parameters from maps, 
     22!!                    as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 
     23!!                    Here, it leads to change dimensions and indices.  
     24!!                    We can also impose kfact_root=1 in all soil layers to cancel the effect of 
     25!!                    roots on ks profile (keyword KFACT_ROOT_TYPE). 
     26!!                  
    2327!! REFERENCE(S) : 
    2428!! - de Rosnay, P., J. Polcher, M. Bruen, and K. Laval, Impact of a physically based soil 
     
    4347!! of a new soil freezing scheme for a land-surface model with physically-based hydrology. 
    4448!! The Cryosphere, 6, 407-430, doi: 10.5194/tc-6-407-2012. \n 
     49!! - Tafasca S. (2020). Evaluation de l’impact des propriétés du sol sur l’hydrologie simulee dans le 
     50!! modÚle ORCHIDEE, PhD thesis, Sorbonne Universite. \n 
    4551!! 
    4652!! SVN          : 
     
    9197  ! one dimension array allocated, computed, saved and got in hydrol module 
    9298  ! Values per soil type 
    93   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: nvan                !! Van Genuchten coeficients n (unitless) 
    94                                                                           ! RK: 1/n=1-m 
    95 !$OMP THREADPRIVATE(nvan)                                                  
    96   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: avan                !! Van Genuchten coeficients a 
    97                                                                          !!  @tex $(mm^{-1})$ @endtex 
    98 !$OMP THREADPRIVATE(avan)                                                 
    99   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcr                 !! Residual volumetric water content 
    100                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex 
    101 !$OMP THREADPRIVATE(mcr)                                                  
    102   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                 !! Saturated volumetric water content 
    103                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex 
    104 !$OMP THREADPRIVATE(mcs)                                                   
    105   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: ks                  !! Hydraulic conductivity at saturation 
    106                                                                          !!  @tex $(mm d^{-1})$ @endtex 
    107 !$OMP THREADPRIVATE(ks)                                                   
     99                                                               
    108100  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pcent               !! Fraction of saturated volumetric soil moisture above  
    109101                                                                         !! which transpir is max (0-1, unitless) 
    110 !$OMP THREADPRIVATE(pcent) 
    111   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcfc                !! Volumetric water content at field capacity 
    112                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex  
    113 !$OMP THREADPRIVATE(mcfc)                                                  
    114   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcw                 !! Volumetric water content at wilting point 
    115                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex  
    116 !$OMP THREADPRIVATE(mcw)                                                  
     102!$OMP THREADPRIVATE(pcent)                                                       
    117103  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mc_awet             !! Vol. wat. cont. above which albedo is cst 
    118104                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex  
     
    216202!$OMP THREADPRIVATE(kfact_root) 
    217203  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: kfact            !! Factor to reduce Ks with depth (unitless) 
    218                                                                          !! DIM = nslm * nscm 
     204                                                                         !! DIM = nslm * kjpindex 
    219205!$OMP THREADPRIVATE(kfact) 
    220206  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: zz               !! Depth of nodes [znh in vertical_soil] transformed into (mm)  
     
    228214  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: mc_lin   !! 50 Vol. Wat. Contents to linearize K and D, for each texture  
    229215                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex 
    230                                                                  !! DIM = imin:imax * nscm 
     216                                                                 !! DIM = imin:imax * kjpindex 
    231217!$OMP THREADPRIVATE(mc_lin) 
    232218  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: k_lin    !! 50 values of unsaturated K, for each soil layer and texture 
    233219                                                                 !!  @tex $(mm d^{-1})$ @endtex 
    234                                                                  !! DIM = imin:imax * nslm * nscm 
     220                                                                 !! DIM = imin:imax * nslm * kjpindex 
    235221!$OMP THREADPRIVATE(k_lin) 
    236222  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: d_lin    !! 50 values of diffusivity D, for each soil layer and texture 
    237223                                                                 !!  @tex $(mm^2 d^{-1})$ @endtex 
    238                                                                  !! DIM = imin:imax * nslm * nscm 
     224                                                                 !! DIM = imin:imax * nslm * kjpindex 
    239225!$OMP THREADPRIVATE(d_lin) 
    240226  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: a_lin    !! 50 values of the slope in K=a*mc+b, for each soil layer and texture 
    241227                                                                 !!  @tex $(mm d^{-1})$ @endtex 
    242                                                                  !! DIM = imin:imax * nslm * nscm 
     228                                                                 !! DIM = imin:imax * nslm * kjpindex 
    243229!$OMP THREADPRIVATE(a_lin) 
    244230  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: b_lin    !! 50 values of y-intercept in K=a*mc+b, for each soil layer and texture 
    245231                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex 
    246                                                                  !! DIM = imin:imax * nslm * nscm 
     232                                                                 !! DIM = imin:imax * nslm * kjpindex 
    247233!$OMP THREADPRIVATE(b_lin) 
    248234 
     
    259245!$OMP THREADPRIVATE(kk) 
    260246  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: avan_mod_tab  !! VG parameter a modified from  exponantial profile 
    261                                                                       !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,nscm) 
     247                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,kjpindex) 
    262248!$OMP THREADPRIVATE(avan_mod_tab)   
    263249  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: nvan_mod_tab  !! VG parameter n  modified from  exponantial profile 
    264                                                                       !! (unitless) !! DIMENSION (nslm,nscm 
     250                                                                      !! (unitless) !! DIMENSION (nslm,kjpindex 
    265251!$OMP THREADPRIVATE(nvan_mod_tab) 
    266252   
     
    450436!_ ================================================================================================================================ 
    451437 
    452   SUBROUTINE hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          & 
     438  SUBROUTINE hydrol_initialize ( ks,             nvan,      avan,          mcr,              & 
     439                                 mcs,            mcfc,      mcw,           kjit,             & 
     440                                 kjpindex,       index,     rest_id,                         & 
    453441                                 njsc,           soiltile,  veget,         veget_max,        & 
    454442                                 humrel,         vegstress, drysoil_frac,                    & 
     
    461449    !! 0. Variable and parameter declaration 
    462450    !! 0.1 Input variables 
     451 
     452 
     453    !salma: added soil params as input variables 
     454    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     455    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless) 
     456    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1}) 
     457    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     458    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     459    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     460    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     461    !end salma: added soil params as input variables 
     462 
    463463    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number  
    464464    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size 
     
    495495    !! 0.4 Local variables 
    496496    INTEGER(i_std)                                       :: jsl 
     497     
     498    !salma: added soil params which depend on soil texture 
     499    REAL(r_std), DIMENSION (nscm)                        :: nvan_text 
     500    REAL(r_std), DIMENSION (nscm)                        :: avan_text 
     501    REAL(r_std), DIMENSION (nscm)                        :: mcr_text 
     502    REAL(r_std), DIMENSION (nscm)                        :: mcs_text 
     503    REAL(r_std), DIMENSION (nscm)                        :: ks_text 
     504    REAL(r_std), DIMENSION (nscm)                        :: pcent_text 
     505    REAL(r_std), DIMENSION (nscm)                        :: mcfc_text 
     506    REAL(r_std), DIMENSION (nscm)                        :: mcw_text 
     507    REAL(r_std), DIMENSION (nscm)                        :: mc_awet_text 
     508    REAL(r_std), DIMENSION (nscm)                        :: mc_adry_text 
     509     !end salma: added soil params which depend on soil texture 
    497510!_ ================================================================================================================================ 
    498511 
    499     CALL hydrol_init (kjit, kjpindex, index, rest_id, veget_max, soiltile, & 
     512    CALL hydrol_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc, kjit, kjpindex, index, rest_id, veget_max, soiltile, & 
    500513         humrel, vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, & 
    501514         snowdz, snowgrain, snowrho,    snowtemp,   snowheat, & 
    502515         drysoil_frac, evap_bare_lim, evap_bare_lim_ns) 
    503516     
    504     CALL hydrol_var_init (kjpindex, veget, veget_max, & 
     517    CALL hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget, veget_max, & 
    505518         soiltile, njsc, mx_eau_var, shumdiag_perma, & 
    506519         drysoil_frac, qsintveg, mc_layh, mcl_layh)  
     
    557570!_ ================================================================================================================================ 
    558571 
    559   SUBROUTINE hydrol_main (kjit, kjpindex, & 
     572  SUBROUTINE hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw,  & 
     573       & kjit, kjpindex, & 
    560574       & index, indexveg, indexsoil, indexlayer, indexnslm, & 
    561575       & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, & 
     
    607621    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration 
    608622    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef 
     623 
     624    !salma: added input soil params: 
     625    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     626    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless) 
     627    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1}) 
     628    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     629    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     630    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     631    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     632    !end salma: added input soil params 
     633  
    609634    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation 
    610635    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot_penm      !! Soil Potential Evaporation Correction 
     
    672697 
    673698    !! 0.4 Local variables 
    674  
     699    !salma: added variable kfact_root_type 
    675700    INTEGER(i_std)                                     :: jst              !! Index of soil tiles (unitless, 1-3) 
    676701    INTEGER(i_std)                                     :: jsl              !! Index of soil layers (unitless) 
    677702    INTEGER(i_std)                                     :: ji, jv 
     703    CHARACTER(LEN=80)                                  :: kfact_root_type  !! read from run.def: when equal to 'cons', it indicates that ks does not increase 
     704                                                                           !!   in the rootzone, ie, kfact_root=1; else, kfact_root defined as usual 
    678705    REAL(r_std),DIMENSION (kjpindex)                   :: soilwet          !! A temporary diagnostic of soil wetness 
    679706    REAL(r_std),DIMENSION (kjpindex)                   :: snowdepth_diag   !! Depth of snow layer containing default values, only for diagnostics 
     
    760787             IF(soiltile(ji,jst) .GT. min_sechiba) THEN 
    761788                kfact_root(ji,jsl,jst) = kfact_root(ji,jsl,jst) * & 
    762                      & MAX((MAXVAL(ks_usda)/ks(njsc(ji)))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), & 
     789                     & MAX((MAXVAL(ks_usda)/ks(ji))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), & 
    763790                     un)  
    764791             ENDIF 
     
    767794    ENDDO 
    768795 
     796    !salma: added config key of KFACT_ROOT_TYPE 
     797    !! KFACT_ROOT_TYPE = cons is used to impose that kfact_root = 1 in every soil layer, 
     798    !! so that ks does not increase in the rootzone; else, kfact_root defined as usual 
     799    !Config Key   = KFACT_ROOT_TYPE 
     800    !Config Desc  = keyword added for spmip exp1 and exp4 to get a constant ks over soil depth and rootzone 
     801    !Config If    = spmip exp1 or exp4 
     802    !Config Def   = var 
     803    !Config Help  = can have two values: 'cons' or 'var'. If var then no changes, if cons then kfact_root=un 
     804    !Config Units = [mm/d] 
     805    CALL getin_p("KFACT_ROOT_TYPE",kfact_root_type) 
     806 
     807    IF (kfact_root_type=='cons') THEN 
     808    kfact_root(:,:,:) = un 
     809    ENDIF 
     810    !end salma: added config key 
     811 
    769812    ! 
    770813    !! 3.3 computes canopy  ==>hydrol_canop 
     
    778821    !! 3.5 computes soil hydrology ==>hydrol_soil 
    779822 
    780     CALL hydrol_soil(kjpindex, veget_max, soiltile, njsc, reinf_slope,  & 
     823    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  & 
    781824         transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &  
    782825         returnflow, reinfiltration, irrigation, & 
     
    869912    !     mcs etc are identical in all layers (no normalization by vegtot to be comparable to mc) 
    870913    DO jsl=1,nslm 
    871        land_mcs(:,jsl) = mcs(njsc(:)) 
    872        land_mcfc(:,jsl) = mcfc(njsc(:)) 
    873        land_mcw(:,jsl) = mcw(njsc(:)) 
    874        land_mcr(:,jsl) = mcr(njsc(:)) 
     914       !salma: replaced mcs(njsc(:)) by mcs(:) and same for other variables 
     915       land_mcs(:,jsl) = mcs(:) 
     916       land_mcfc(:,jsl) = mcfc(:) 
     917       land_mcw(:,jsl) = mcw(:) 
     918       land_mcr(:,jsl) = mcr(:) 
    875919    ENDDO 
    876920    CALL xios_orchidee_send_field("mcs",land_mcs) ! in m3/m3 
     
    878922    CALL xios_orchidee_send_field("mcw",land_mcw) ! in m3/m3  
    879923    CALL xios_orchidee_send_field("mcr",land_mcr) ! in m3/m3  
    880            
     924 
     925       
    881926    CALL xios_orchidee_send_field("water2infilt",water2infilt)    
    882927    CALL xios_orchidee_send_field("mc",mc) 
     
    915960 
    916961    ! For the soil wetness above wilting point for CMIP6 (mrsow) 
    917     mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(njsc(:)) ) & 
    918          / ( zmaxh*mille*( mcs(njsc(:)) - mcw(njsc(:)) ) ) 
     962    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(:) ) & 
     963         / ( zmaxh*mille*( mcs(:) - mcw(:) ) ) 
    919964    CALL xios_orchidee_send_field("mrsow",mrsow) 
    920965 
     
    12881333!_ ================================================================================================================================ 
    12891334!!_ hydrol_init 
    1290  
    1291   SUBROUTINE hydrol_init(kjit, kjpindex, index, rest_id, veget_max, soiltile, & 
     1335!salma: added variables and njsc in arguments and input variables 
     1336 
     1337  SUBROUTINE hydrol_init(ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc,& 
     1338       kjit, kjpindex, index, rest_id, veget_max, soiltile, & 
    12921339       humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, & 
    12931340       snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, & 
     
    12981345 
    12991346    !! 0.1 Input variables 
    1300  
     1347    !salma introduced njsc variable 
     1348    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc               !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 
    13011349    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number  
    13021350    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size 
     
    13051353    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max          !! Carte de vegetation max 
    13061354    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in)  :: soiltile           !! Fraction of each soil tile within vegtot (0-1, unitless) 
    1307  
     1355    
    13081356    !! 0.2 Output variables 
     1357 
     1358    !salma: added input variables 
     1359    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     1360    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: nvan             !! Van Genuchten coeficients n (unitless) 
     1361    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: avan             !! Van Genuchten coeficients a (mm-1}) 
     1362    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     1363    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     1364    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     1365    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
    13091366 
    13101367    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: humrel             !! Stress hydrique, relative humidity 
     
    14071464    !! 2.1 array allocation for soil texture 
    14081465 
    1409     ALLOCATE (nvan(nscm),stat=ier) 
    1410     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan','','') 
    1411  
    1412     ALLOCATE (avan(nscm),stat=ier) 
    1413     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan','','') 
    1414  
    1415     ALLOCATE (mcr(nscm),stat=ier) 
    1416     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcr','','') 
    1417  
    1418     ALLOCATE (mcs(nscm),stat=ier) 
    1419     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcs','','') 
    1420  
    1421     ALLOCATE (ks(nscm),stat=ier) 
    1422     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ks','','') 
    1423  
    1424     ALLOCATE (pcent(nscm),stat=ier) 
     1466    
     1467    ALLOCATE (pcent(kjpindex),stat=ier) 
    14251468    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','') 
    1426  
    1427     ALLOCATE (mcfc(nscm),stat=ier) 
    1428     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcfc','','') 
    1429  
    1430     ALLOCATE (mcw(nscm),stat=ier) 
    1431     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcw','','') 
    1432      
    1433     ALLOCATE (mc_awet(nscm),stat=ier) 
     1469     
     1470    ALLOCATE (mc_awet(kjpindex),stat=ier) 
    14341471    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','') 
    14351472 
    1436     ALLOCATE (mc_adry(nscm),stat=ier) 
     1473    ALLOCATE (mc_adry(kjpindex),stat=ier) 
    14371474    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','') 
    14381475        
     
    14421479    CASE (3) 
    14431480               
    1444        nvan(:) = nvan_fao(:)        
    1445        avan(:) = avan_fao(:) 
    1446        mcr(:) = mcr_fao(:) 
    1447        mcs(:) = mcs_fao(:) 
    1448        ks(:) = ks_fao(:) 
    1449        pcent(:) = pcent_fao(:) 
    1450        mcfc(:) = mcf_fao(:) 
    1451        mcw(:) = mcw_fao(:) 
    1452        mc_awet(:) = mc_awet_fao(:) 
    1453        mc_adry(:) = mc_adry_fao(:) 
    1454     CASE (12) 
    1455         
    1456        nvan(:) = nvan_usda(:) 
    1457        avan(:) = avan_usda(:) 
    1458        mcr(:) = mcr_usda(:) 
    1459        mcs(:) = mcs_usda(:) 
    1460        ks(:) = ks_usda(:) 
    1461        pcent(:) = pcent_usda(:) 
    1462        mcfc(:) = mcf_usda(:) 
    1463        mcw(:) = mcw_usda(:) 
    1464        mc_awet(:) = mc_awet_usda(:) 
    1465        mc_adry(:) = mc_adry_usda(:) 
     1481       pcent(:) = pcent_fao(njsc(:)) 
     1482       mc_awet(:) = mc_awet_fao(njsc(:)) 
     1483       mc_adry(:) = mc_adry_fao(njsc(:)) 
     1484    CASE (13) 
     1485            
     1486       pcent(:) = pcent_usda(njsc(:)) 
     1487       mc_awet(:) = mc_awet_usda(njsc(:)) 
     1488       mc_adry(:) = mc_adry_usda(njsc(:)) 
    14661489        
    14671490    CASE DEFAULT 
     
    14741497    !! 2.3 Read in the run.def the parameters values defined by the user 
    14751498 
    1476     !Config Key   = CWRR_N_VANGENUCHTEN 
    1477     !Config Desc  = Van genuchten coefficient n 
    1478     !Config If    =  
    1479     !Config Def   = 1.89, 1.56, 1.31 
    1480     !Config Help  = This parameter will be constant over the entire  
    1481     !Config         simulated domain, thus independent from soil 
    1482     !Config         texture.    
    1483     !Config Units = [-] 
    1484     CALL getin_p("CWRR_N_VANGENUCHTEN",nvan) 
    1485  
    1486     !! Check parameter value (correct range) 
    1487     IF ( ANY(nvan(:) <= zero) ) THEN 
    1488        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1489             &     "Wrong parameter value for CWRR_N_VANGENUCHTEN.", & 
    1490             &     "This parameter should be positive. ", & 
    1491             &     "Please, check parameter value in run.def. ") 
    1492     END IF 
    1493  
    1494  
    1495     !Config Key   = CWRR_A_VANGENUCHTEN 
    1496     !Config Desc  = Van genuchten coefficient a 
    1497     !Config If    =  
    1498     !Config Def   = 0.0075, 0.0036, 0.0019 
    1499     !Config Help  = This parameter will be constant over the entire  
    1500     !Config         simulated domain, thus independent from soil 
    1501     !Config         texture.    
    1502     !Config Units = [1/mm]   
    1503     CALL getin_p("CWRR_A_VANGENUCHTEN",avan) 
    1504  
    1505     !! Check parameter value (correct range) 
    1506     IF ( ANY(avan(:) <= zero) ) THEN 
    1507        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1508             &     "Wrong parameter value for CWRR_A_VANGENUCHTEN.", & 
    1509             &     "This parameter should be positive. ", & 
    1510             &     "Please, check parameter value in run.def. ") 
    1511     END IF 
    1512  
    1513  
    1514     !Config Key   = VWC_RESIDUAL 
    1515     !Config Desc  = Residual soil water content 
    1516     !Config If    =  
    1517     !Config Def   = 0.065, 0.078, 0.095 
    1518     !Config Help  = This parameter will be constant over the entire  
    1519     !Config         simulated domain, thus independent from soil 
    1520     !Config         texture.    
    1521     !Config Units = [m3/m3]   
    1522     CALL getin_p("VWC_RESIDUAL",mcr) 
    1523  
    1524     !! Check parameter value (correct range) 
    1525     IF ( ANY(mcr(:) < zero) .OR. ANY(mcr(:) > 1.)  ) THEN 
    1526        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1527             &     "Wrong parameter value for VWC_RESIDUAL.", & 
    1528             &     "This parameter is ranged between 0 and 1. ", & 
    1529             &     "Please, check parameter value in run.def. ") 
    1530     END IF 
    1531  
    1532      
    1533     !Config Key   = VWC_SAT 
    1534     !Config Desc  = Saturated soil water content 
    1535     !Config If    =  
    1536     !Config Def   = 0.41, 0.43, 0.41 
    1537     !Config Help  = This parameter will be constant over the entire  
    1538     !Config         simulated domain, thus independent from soil 
    1539     !Config         texture.    
    1540     !Config Units = [m3/m3]   
    1541     CALL getin_p("VWC_SAT",mcs) 
    1542  
    1543     !! Check parameter value (correct range) 
    1544     IF ( ANY(mcs(:) < zero) .OR. ANY(mcs(:) > 1.) .OR. ANY(mcs(:) <= mcr(:)) ) THEN 
    1545        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1546             &     "Wrong parameter value for VWC_SAT.", & 
    1547             &     "This parameter should be greater than VWC_RESIDUAL and less than 1. ", & 
    1548             &     "Please, check parameter value in run.def. ") 
    1549     END IF 
    1550  
    1551  
    1552     !Config Key   = CWRR_KS  
    1553     !Config Desc  = Hydraulic conductivity Saturation 
    1554     !Config If    =  
    1555     !Config Def   = 1060.8, 249.6, 62.4 
    1556     !Config Help  = This parameter will be constant over the entire  
    1557     !Config         simulated domain, thus independent from soil 
    1558     !Config         texture.    
    1559     !Config Units = [mm/d]    
    1560     CALL getin_p("CWRR_KS",ks) 
    1561  
    1562     !! Check parameter value (correct range) 
    1563     IF ( ANY(ks(:) <= zero) ) THEN 
    1564        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1565             &     "Wrong parameter value for CWRR_KS.", & 
    1566             &     "This parameter should be positive. ", & 
    1567             &     "Please, check parameter value in run.def. ") 
    1568     END IF 
     1499 
     1500 
     1501 
     1502 
     1503 
     1504 
    15691505 
    15701506 
     
    15871523 
    15881524 
    1589     !Config Key   = VWC_FC  
    1590     !Config Desc  = Volumetric water content field capacity 
    1591     !Config If    =  
    1592     !Config Def   = 0.32, 0.32, 0.32 
    1593     !Config Help  = This parameter is independent from soil texture for 
    1594     !Config         the time being. 
    1595     !Config Units = [m3/m3]    
    1596     CALL getin_p("VWC_FC",mcfc) 
    1597  
    1598     !! Check parameter value (correct range) 
    1599     IF ( ANY(mcfc(:) > mcs(:)) ) THEN 
    1600        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1601             &     "Wrong parameter value for VWC_FC.", & 
    1602             &     "This parameter should be less than VWC_SAT. ", & 
    1603             &     "Please, check parameter value in run.def. ") 
    1604     END IF 
    1605  
    1606  
    1607     !Config Key   = VWC_WP 
    1608     !Config Desc  = Volumetric water content Wilting pt 
    1609     !Config If    =  
    1610     !Config Def   = 0.10, 0.10, 0.10  
    1611     !Config Help  = This parameter is independent from soil texture for 
    1612     !Config         the time being. 
    1613     !Config Units = [m3/m3]    
    1614     CALL getin_p("VWC_WP",mcw) 
    1615  
    1616     !! Check parameter value (correct range) 
    1617     IF ( ANY(mcw(:) > mcfc(:)) .OR. ANY(mcw(:) < mcr(:)) ) THEN 
    1618        CALL ipslerr_p(error_level, "hydrol_init.", & 
    1619             &     "Wrong parameter value for VWC_WP.", & 
    1620             &     "This parameter should be greater or equal than VWC_RESIDUAL and less or equal than VWC_SAT.", & 
    1621             &     "Please, check parameter value in run.def. ") 
    1622     END IF 
    1623  
     1525    
    16241526 
    16251527    !Config Key   = VWC_MIN_FOR_WET_ALB 
     
    17451647    kk(:,:,:) = 276.48 
    17461648     
    1747     ALLOCATE (avan_mod_tab(nslm,nscm),stat=ier)  
     1649    ALLOCATE (avan_mod_tab(nslm,kjpindex),stat=ier)  
    17481650    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan_mod_tab','','') 
    17491651     
    1750     ALLOCATE (nvan_mod_tab(nslm,nscm),stat=ier)  
     1652    ALLOCATE (nvan_mod_tab(nslm,kjpindex),stat=ier)  
    17511653    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan_mod_tab','','') 
    17521654 
     
    19321834    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','') 
    19331835 
    1934     ALLOCATE (kfact(nslm, nscm),stat=ier) 
     1836    ALLOCATE (kfact(nslm, kjpindex),stat=ier) 
    19351837    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','') 
    19361838 
     
    19441846    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','') 
    19451847 
    1946     ALLOCATE (mc_lin(imin:imax, nscm),stat=ier) 
     1848    ALLOCATE (mc_lin(imin:imax, kjpindex),stat=ier) 
    19471849    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','') 
    19481850 
    1949     ALLOCATE (k_lin(imin:imax, nslm, nscm),stat=ier) 
     1851    ALLOCATE (k_lin(imin:imax, nslm, kjpindex),stat=ier) 
    19501852    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','') 
    19511853 
    1952     ALLOCATE (d_lin(imin:imax, nslm, nscm),stat=ier) 
     1854    ALLOCATE (d_lin(imin:imax, nslm, kjpindex),stat=ier) 
    19531855    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','') 
    19541856 
    1955     ALLOCATE (a_lin(imin:imax, nslm, nscm),stat=ier) 
     1857    ALLOCATE (a_lin(imin:imax, nslm, kjpindex),stat=ier) 
    19561858    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','') 
    19571859 
    1958     ALLOCATE (b_lin(imin:imax, nslm, nscm),stat=ier) 
     1860    ALLOCATE (b_lin(imin:imax, nslm, kjpindex),stat=ier) 
    19591861    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','') 
    19601862 
     
    24512353 
    24522354    ! Allocation for soiltile related parameters 
    2453     IF ( ALLOCATED (nvan)) DEALLOCATE (nvan) 
    2454     IF ( ALLOCATED (avan)) DEALLOCATE (avan) 
    2455     IF ( ALLOCATED (mcr)) DEALLOCATE (mcr) 
    2456     IF ( ALLOCATED (mcs)) DEALLOCATE (mcs) 
    2457     IF ( ALLOCATED (ks)) DEALLOCATE (ks) 
     2355    
    24582356    IF ( ALLOCATED (pcent)) DEALLOCATE (pcent) 
    2459     IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc) 
    2460     IF ( ALLOCATED (mcw)) DEALLOCATE (mcw) 
    24612357    IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet) 
    24622358    IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry) 
     
    28112707!_ hydrol_var_init 
    28122708 
    2813   SUBROUTINE hydrol_var_init (kjpindex, veget, veget_max, soiltile, njsc, & 
     2709  SUBROUTINE hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 
     2710       kjpindex, veget, veget_max, soiltile, njsc, & 
    28142711       mx_eau_var, shumdiag_perma, & 
    28152712       drysoil_frac, qsintveg, mc_layh, mcl_layh)  
     
    28202717 
    28212718    !! 0.1 Input variables 
     2719    !salma: added the following soil parameters 
     2720    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     2721    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless) 
     2722    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1}) 
     2723    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     2724    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     2725    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     2726    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
    28222727 
    28232728    ! input scalar  
     
    28552760    REAL(r_std)                                         :: nvan_mod      !! VG parameter n  modified from  exponantial profile 
    28562761                                                                         !! (unitless) 
    2857     REAL(r_std), DIMENSION(nslm,nscm)                   :: afact, nfact  !! Multiplicative factor for decay of a and n with depth 
     2762    REAL(r_std), DIMENSION(nslm,kjpindex)               :: afact, nfact  !! Multiplicative factor for decay of a and n with depth 
    28582763                                                                         !! (unitless) 
    28592764    ! parameters for "soil densification" with depth 
     
    28832788    REAL(r_std), DIMENSION (kjpindex,nslm)              :: nvg             !! VG param n modified with depth at each node 
    28842789                                                                           !! (unitless) 
    2885     INTEGER(i_std)                                      :: jiref           !! To identify the mc_lins where k_lin and d_lin  
     2790                                                                           !! need special treatment 
     2791    !salma: added local variable ii and jiref replaced with iiref 
     2792    INTEGER(i_std)                                      :: ii 
     2793    INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin 
    28862794                                                                           !! need special treatment 
    28872795 
     
    30502958    END IF 
    30512959 
     2960  
     2961 
    30522962    !- 
    30532963    !! 3 Compute the profile for a and n 
    30542964    !- 
    3055  
    3056     ! For every soil texture 
    3057     DO jsc = 1, nscm  
     2965    !salma: for every grid cell 
     2966    DO ji = 1, kjpindex 
    30582967       DO jsl=1,nslm 
    30592968          ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2 
    30602969          ! Calibrated against Hapex-Sahel measurements 
    3061           kfact(jsl,jsc) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un) 
    3062           ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14  
    3063            
    3064           nfact(jsl,jsc) = ( kfact(jsl,jsc) )**nk_rel 
    3065           afact(jsl,jsc) = ( kfact(jsl,jsc) )**ak_rel 
    3066        ENDDO 
    3067     ENDDO 
    3068  
    3069     ! For every soil texture 
    3070     DO jsc = 1, nscm 
     2970          kfact(jsl,ji) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un) 
     2971          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14 
     2972 
     2973          nfact(jsl,ji) = ( kfact(jsl,ji) )**nk_rel 
     2974          afact(jsl,ji) = ( kfact(jsl,ji) )**ak_rel 
     2975       ENDDO 
     2976    ENDDO 
     2977    !end salma 
     2978     
     2979    ! For every grid cell 
     2980     DO ji = 1, kjpindex 
    30712981       !- 
    30722982       !! 4 Compute the linearized values of k, a, b and d 
     
    30802990 
    30812991       ! We define 51 bounds for 50 bins of mc between mcr and mcs 
    3082        mc_lin(imin,jsc)=mcr(jsc) 
    3083        mc_lin(imax,jsc)=mcs(jsc) 
    3084        DO ji= imin+1, imax-1 ! ji=2,50 
    3085           mc_lin(ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin) 
     2992       mc_lin(imin,ji)=mcr(ji) 
     2993       mc_lin(imax,ji)=mcs(ji) 
     2994       DO ii= imin+1, imax-1 ! ii=2,50 
     2995          mc_lin(ii,ji) = mcr(ji) + (ii-imin)*(mcs(ji)-mcr(ji))/(imax-imin) 
    30862996       ENDDO 
    30872997 
    30882998       DO jsl = 1, nslm 
    30892999          ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42 
    3090           nvan_mod = n0 + (nvan(jsc)-n0) * nfact(jsl,jsc) 
    3091           avan_mod = a0 + (avan(jsc)-a0) * afact(jsl,jsc) 
     3000          nvan_mod = n0 + (nvan(ji)-n0) * nfact(jsl,ji) 
     3001          avan_mod = a0 + (avan(ji)-a0) * afact(jsl,ji) 
    30923002          m = un - un / nvan_mod 
    30933003          ! Creation of arrays for SP-MIP output by landpoint 
    3094           nvan_mod_tab(jsl,jsc) = nvan_mod 
    3095           avan_mod_tab(jsl,jsc) = avan_mod 
    3096           ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(jsc) * kfact(jsl,jsc) 
    3097           DO ji = imax,imin,-1  
    3098              frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc))) 
    3099              k_lin(ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2 
     3004          nvan_mod_tab(jsl,ji) = nvan_mod 
     3005          avan_mod_tab(jsl,ji) = avan_mod 
     3006          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(ji) * kfact(jsl,ji) 
     3007          DO ii = imax,imin,-1 
     3008             frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 
     3009             k_lin(ii,jsl,ji) = ks(ji) * kfact(jsl,ji) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2 
    31003010          ENDDO 
    31013011 
    31023012          ! k_lin should not be zero, nor too small 
    3103           ! We track jiref, the bin under which mc is too small and we may get zero k_lin      
    3104           ji=imax-1 
    3105           DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0)) 
    3106              jiref=ji 
    3107              ji=ji-1 
     3013          ! We track iiref, the bin under which mc is too small and we may get zero k_lin 
     3014          !salma: ji replaced with ii and jiref replaced with iiref and jsc with ji 
     3015          ii=imax-1 
     3016          DO WHILE ((k_lin(ii,jsl,ji) > 1.e-32) .and. (ii>0)) 
     3017             iiref=ii 
     3018             ii=ii-1 
    31083019          ENDDO 
    3109           DO ji=jiref-1,imin,-1 
    3110              k_lin(ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10. 
     3020          DO ii=iiref-1,imin,-1 
     3021             k_lin(ii,jsl,ji)=k_lin(ii+1,jsl,ji)/10. 
    31113022          ENDDO 
    3112           
    3113           DO ji = imin,imax-1 ! ji=1,50 
     3023 
     3024          DO ii = imin,imax-1 ! ii=1,50 
    31143025             ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin 
    3115              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)) 
    3116              b_lin(ji,jsl,jsc)  = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc) 
     3026             a_lin(ii,jsl,ji) = (k_lin(ii+1,jsl,ji)-k_lin(ii,jsl,ji)) / (mc_lin(ii+1,ji)-mc_lin(ii,ji)) 
     3027             b_lin(ii,jsl,ji)  = k_lin(ii,jsl,ji) - a_lin(ii,jsl,ji)*mc_lin(ii,ji) 
    31173028 
    31183029             ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta) 
    3119              ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin  
    3120              IF (ji.NE.imin .AND. ji.NE.imax-1) THEN 
    3121                 frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc))) 
    3122                 d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) *  & 
    3123                      ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) * & 
     3030             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin 
     3031             IF (ii.NE.imin .AND. ii.NE.imax-1) THEN 
     3032                frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 
     3033                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) *  & 
     3034                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * & 
    31243035                     (  frac**(-un/m) -un ) ** (-m) 
    3125                 frac=MIN(un,(mc_lin(ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc))) 
    3126                 d_lin(ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*& 
    3127                      ( (frac**(-un/m))/(mc_lin(ji+1,jsc)-mcr(jsc)) ) * & 
     3036                frac=MIN(un,(mc_lin(ii+1,ji)-mcr(ji))/(mcs(ji)-mcr(ji))) 
     3037                d_lin(ii+1,jsl,ji) =(k_lin(ii+1,jsl,ji) / (avan_mod*m*nvan_mod))*& 
     3038                     ( (frac**(-un/m))/(mc_lin(ii+1,ji)-mcr(ji)) ) * & 
    31283039                     (  frac**(-un/m) -un ) ** (-m) 
    3129                 d_lin(ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc)) 
    3130              ELSE IF(ji.EQ.imax-1) THEN 
    3131                 d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * & 
    3132                      ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) *  & 
     3040                d_lin(ii,jsl,ji) = undemi * (d_lin(ii,jsl,ji)+d_lin(ii+1,jsl,ji)) 
     3041             ELSE IF(ii.EQ.imax-1) THEN 
     3042                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * & 
     3043                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) *  & 
    31333044                     (  frac**(-un/m) -un ) ** (-m) 
    31343045             ENDIF 
    3135           ENDDO 
    3136  
    3137           ! Special case for ji=imin 
    3138           d_lin(imin,jsl,jsc) = d_lin(imin+1,jsl,jsc)/1000. 
     3046          ENDDO !Salma end loop over landpoints 
     3047 
     3048          ! Special case for ii=imin 
     3049          d_lin(imin,jsl,ji) = d_lin(imin+1,jsl,ji)/1000. 
    31393050 
    31403051          ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations 
    31413052          ! We don't want d_lin = zero 
    3142           DO ji=jiref-1,imin,-1 
    3143              d_lin(ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10. 
     3053          DO ii=iiref-1,imin,-1 
     3054             d_lin(ii,jsl,ji)=d_lin(ii+1,jsl,ji)/10. 
    31443055          ENDDO 
    31453056 
    31463057       ENDDO 
    31473058    ENDDO 
     3059 
    31483060 
    31493061    ! Output of alphavg and nvg at each node for SP-MIP 
    31503062    DO jsl = 1, nslm 
    3151        alphavg(:,jsl) = avan_mod_tab(jsl,njsc(:))*1000. ! from mm-1 to m-1 
    3152        nvg(:,jsl) = nvan_mod_tab(jsl,njsc(:)) 
     3063       alphavg(:,jsl) = avan_mod_tab(jsl,:)*1000. ! from mm-1 to m-1 
     3064       nvg(:,jsl) = nvan_mod_tab(jsl,:) 
    31533065    ENDDO 
    31543066    CALL xios_orchidee_send_field("alphavg",alphavg) ! in m-1 
     
    31653077 
    31663078    mx_eau_var(:) = zero 
    3167     mx_eau_var(:) = zmaxh*mille*mcs(njsc(:))  
    3168  
    3169     DO ji = 1,kjpindex  
     3079    mx_eau_var(:) = zmaxh*mille*mcs(:) 
     3080 
     3081    DO ji = 1,kjpindex 
    31703082       IF (vegtot(ji) .LE. zero) THEN 
    31713083          mx_eau_var(ji) = mx_eau_nobio*zmaxh 
     
    31813093 
    31823094    ! Loop on soiltiles to compute the variables (ji,jst) 
    3183     DO jst=1,nstm  
     3095    DO jst=1,nstm 
    31843096       DO ji = 1, kjpindex 
    3185           tmcs(ji,jst)=zmaxh* mille*mcs(njsc(ji)) 
    3186           tmcr(ji,jst)=zmaxh* mille*mcr(njsc(ji)) 
    3187           tmcfc(ji,jst)=zmaxh* mille*mcfc(njsc(ji)) 
    3188           tmcw(ji,jst)=zmaxh* mille*mcw(njsc(ji)) 
    3189        ENDDO 
    3190     ENDDO 
    3191         
     3097          tmcs(ji,jst)=zmaxh* mille*mcs(ji) 
     3098          tmcr(ji,jst)=zmaxh* mille*mcr(ji) 
     3099          tmcfc(ji,jst)=zmaxh* mille*mcfc(ji) 
     3100          tmcw(ji,jst)=zmaxh* mille*mcw(ji) 
     3101       ENDDO 
     3102    ENDDO 
     3103 
    31923104    ! The total soil moisture for each soiltile: 
    31933105    DO jst=1,nstm 
     
    31973109    ENDDO 
    31983110 
    3199     DO jst=1,nstm  
     3111    DO jst=1,nstm 
    32003112       DO jsl=2,nslm-1 
    32013113          DO ji=1,kjpindex 
     
    32063118    ENDDO 
    32073119 
    3208     DO jst=1,nstm  
     3120    DO jst=1,nstm 
    32093121       DO ji=1,kjpindex 
    32103122          tmc(ji,jst) = tmc(ji,jst) +  dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit 
     
    32133125    END DO 
    32143126 
    3215 !JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.     
     3127!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty. 
    32163128!    ! If veget has been updated before restart (with LAND USE or DGVM), 
    32173129!    ! tmc and mc must be modified with respect to humtot conservation. 
     
    32203132    ! The litter variables: 
    32213133    ! level 1 
    3222     DO jst=1,nstm  
     3134    DO jst=1,nstm 
    32233135       DO ji=1,kjpindex 
    32243136          tmc_litter(ji,jst) = dz(2) * (trois*mcl(ji,1,jst)+mcl(ji,2,jst))/huit 
    3225           tmc_litter_wilt(ji,jst) = dz(2) * mcw(njsc(ji)) / deux 
    3226           tmc_litter_res(ji,jst) = dz(2) * mcr(njsc(ji)) / deux 
    3227           tmc_litter_field(ji,jst) = dz(2) * mcfc(njsc(ji)) / deux 
    3228           tmc_litter_sat(ji,jst) = dz(2) * mcs(njsc(ji)) / deux 
    3229           tmc_litter_awet(ji,jst) = dz(2) * mc_awet(njsc(ji)) / deux 
    3230           tmc_litter_adry(ji,jst) = dz(2) * mc_adry(njsc(ji)) / deux 
     3137          tmc_litter_wilt(ji,jst) = dz(2) * mcw(ji) / deux 
     3138          tmc_litter_res(ji,jst) = dz(2) * mcr(ji) / deux 
     3139          tmc_litter_field(ji,jst) = dz(2) * mcfc(ji) / deux 
     3140          tmc_litter_sat(ji,jst) = dz(2) * mcs(ji) / deux 
     3141          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(ji) / deux 
     3142          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(ji) / deux 
    32313143       ENDDO 
    32323144    END DO 
    32333145    ! sum from level 2 to 4 
    3234     DO jst=1,nstm  
     3146    DO jst=1,nstm 
    32353147       DO jsl=2,4 
    32363148          DO ji=1,kjpindex 
    3237              tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * &  
     3149             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 
    32383150                  & ( trois*mcl(ji,jsl,jst) + mcl(ji,jsl-1,jst))/huit & 
    32393151                  & + dz(jsl+1)*(trois*mcl(ji,jsl,jst) + mcl(ji,jsl+1,jst))/huit 
    32403152             tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + & 
    3241                   &(dz(jsl)+ dz(jsl+1))*&  
    3242                   & mcw(njsc(ji))/deux 
     3153                  &(dz(jsl)+ dz(jsl+1))*& 
     3154                  & mcw(ji)/deux 
    32433155             tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + & 
    3244                   &(dz(jsl)+ dz(jsl+1))*&  
    3245                   & mcr(njsc(ji))/deux 
     3156                  &(dz(jsl)+ dz(jsl+1))*& 
     3157                  & mcr(ji)/deux 
    32463158             tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + & 
    3247                   &(dz(jsl)+ dz(jsl+1))* &  
    3248                   & mcs(njsc(ji))/deux 
     3159                  &(dz(jsl)+ dz(jsl+1))* & 
     3160                  & mcs(ji)/deux 
    32493161             tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + & 
    3250                   & (dz(jsl)+ dz(jsl+1))* &  
    3251                   & mcfc(njsc(ji))/deux 
     3162                  & (dz(jsl)+ dz(jsl+1))* & 
     3163                  & mcfc(ji)/deux 
    32523164             tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + & 
    3253                   &(dz(jsl)+ dz(jsl+1))* &  
    3254                   & mc_awet(njsc(ji))/deux 
     3165                  &(dz(jsl)+ dz(jsl+1))* & 
     3166                  & mc_awet(ji)/deux 
    32553167             tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + & 
    3256                   & (dz(jsl)+ dz(jsl+1))* &  
    3257                   & mc_adry(njsc(ji))/deux 
     3168                  & (dz(jsl)+ dz(jsl+1))* & 
     3169                  & mc_adry(ji)/deux 
    32583170          END DO 
    32593171       END DO 
     
    32613173 
    32623174 
    3263     DO jst=1,nstm  
     3175    DO jst=1,nstm 
    32643176       DO ji=1,kjpindex 
    32653177          ! here we set that humrelv=0 in PFT1 
    3266           humrelv(ji,1,jst) = zero 
     3178         humrelv(ji,1,jst) = zero 
    32673179       ENDDO 
    32683180    END DO 
     
    32703182 
    32713183    ! Calculate shumdiag_perma for thermosoil 
    3272     ! Use resdist instead of soiltile because we here need to have  
     3184    ! Use resdist instead of soiltile because we here need to have 
    32733185    ! shumdiag_perma at the value from previous time step. 
    32743186    ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma 
     
    32903202    ENDDO 
    32913203    DO ji=1,kjpindex 
    3292         soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average  
    3293     ENDDO 
    3294      
     3204        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average 
     3205    ENDDO 
     3206 
    32953207    ! -- shumdiag_perma for restart 
    3296     !  For consistency with hydrol_soil, we want to calculate a grid-cell average 
     3208   !  For consistency with hydrol_soil, we want to calculate a grid-cell average 
    32973209    DO jsl = 1, nslm 
    3298        DO ji=1,kjpindex         
    3299           shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji))) 
    3300           shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero)  
    3301        ENDDO 
    3302     ENDDO 
    3303                 
     3210       DO ji=1,kjpindex 
     3211          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji)) 
     3212          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
     3213       ENDDO 
     3214    ENDDO 
     3215 
    33043216    ! Calculate drysoil_frac if it was not found in the restart file 
    33053217    ! For simplicity, we set drysoil_frac to 0.5 in this case 
     
    33103222    END IF 
    33113223 
    3312     !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in  
    3313     !! thermosoil for the thermal conductivity.  
     3224    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in 
     3225    !! thermosoil for the thermal conductivity. 
    33143226    ! These values are only used in thermosoil_init in absence of a restart file 
     3227 
    33153228    mc_layh(:,:) = zero 
    33163229    mcl_layh(:,:) = zero 
     3230       
    33173231    DO jst=1,nstm 
    33183232       DO jsl=1,nslm 
     
    36683582 
    36693583  END SUBROUTINE hydrol_flood 
    3670  
    36713584 
    36723585!! ================================================================================================================================ 
     
    37343647!_ ================================================================================================================================ 
    37353648!_ hydrol_soil 
    3736  
    3737   SUBROUTINE hydrol_soil (kjpindex, veget_max, soiltile, njsc, reinf_slope, & 
     3649  SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, & 
     3650       kjpindex, veget_max, soiltile, njsc, reinf_slope, & 
    37383651       & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
    37393652       & returnflow, reinfiltration, irrigation, & 
     
    37483661 
    37493662    !! 0.1 Input variables 
     3663    !salma added soil params to arguments and input variables 
     3664    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     3665    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless) 
     3666    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1}) 
     3667    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     3668    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     3669    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     3670    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     3671 
    37503672 
    37513673    INTEGER(i_std), INTENT(in)                               :: kjpindex  
     
    38073729                                                                                 !! averaged over the mesh (for thermosoil) 
    38083730                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex  
    3809  
    38103731    !! 0.3 Modified variables 
    38113732 
     
    39533874       ! These values will be kept till the end of the prognostic loop 
    39543875       DO jst=1,nstm 
    3955           CALL hydrol_soil_froz(kjpindex,jst,njsc) 
     3876          CALL hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,jst,njsc) 
    39563877       ENDDO 
    39573878 
     
    40804001       !! K and D are computed for the profile of mc before infiltration 
    40814002       !! They depend on the fraction of soil ice, given by profil_froz_hydro_ns 
    4082        CALL hydrol_soil_coef(kjpindex,jst,njsc) 
     4003       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc) 
    40834004 
    40844005       !! Infiltration and surface runoff are computed 
     
    40864007       !! The conductivity comes from hydrol_soil_coef and relates to the liquid phase only 
    40874008       !  This seems consistent with ok_freeze 
    4088        CALL hydrol_soil_infilt(kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, & 
     4009       CALL hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, & 
    40894010            check_infilt_ns) 
    40904011       ru_ns(:,jst) = ru_infilt_ns(:,jst)  
     
    41124033       !! 2.1 K and D are recomputed after infiltration 
    41134034       !! They depend on the fraction of soil ice, still given by profil_froz_hydro_ns  
    4114        CALL hydrol_soil_coef(kjpindex,jst,njsc) 
     4035       CALL hydrol_soil_coef(mcr, mcs,kjpindex,jst,njsc) 
    41154036  
    41164037       !! 2.2 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme 
     
    41214042       DO jsl = 1, nslm 
    41224043          DO ji =1, kjpindex 
    4123              mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + & 
    4124                   (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) ) 
     4044             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 
     4045                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 
    41254046             ! we always have mcl<=mc 
    41264047             ! if mc>mcr, then mcl>mcr; if mc=mcr,mcl=mcr; if mc<mcr, then mcl<mcr 
     
    41434064          DO ji = 1, kjpindex 
    41444065             IF (soiltile(ji,jst) .GT. zero) THEN 
    4145                 mvg = un - un / nvan_mod_tab(jsl,njsc(ji)) 
    4146                 avg = avan_mod_tab(jsl,njsc(ji))*1000. ! to convert in m-1 
    4147                 mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(njsc(ji)))/(mcs(njsc(ji)) - mcr(njsc(ji))) )**(-un/mvg) 
    4148                 psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,njsc(ji))) / avg ! in m  
     4066                mvg = un - un / nvan_mod_tab(jsl,ji) 
     4067                avg = avan_mod_tab(jsl,ji)*1000. ! to convert in m-1 
     4068                mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(ji))/(mcs(ji) - mcr(ji)) )**(-un/mvg) 
     4069                psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,ji)) / avg ! in m  
    41494070                psi_moy(ji,jsl) = psi_moy(ji,jsl) + soiltile(ji,jst) * psi ! average across soil tiles 
    41504071             ENDIF 
     
    42744195          DO ji =1, kjpindex 
    42754196             mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + & 
    4276                   profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) ) 
     4197                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) ) 
    42774198             ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl 
    42784199          ENDDO 
     
    42864207       dr_corr_ns(:,jst) = zero 
    42874208       ru_corr_ns(:,jst) = zero 
    4288        call hydrol_soil_smooth_over_mcs2(kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns) 
     4209       call hydrol_soil_smooth_over_mcs2(mcs, kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns) 
    42894210        
    42904211       ! In absence of freezing, if F is large enough, the correction of oversaturation is sent to drainage        
     
    43304251                dmc(ji,jsl) = zero 
    43314252                IF ( ( zz(jsl) >= zwt_force(ji,jst)*mille ) ) THEN 
    4332                    dmc(ji,jsl) = mcs(njsc(ji)) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value 
    4333                    mc(ji,jsl,jst) = mcs(njsc(ji)) 
     4253                   dmc(ji,jsl) = mcs(ji) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value 
     4254                   mc(ji,jsl,jst) = mcs(ji) 
    43344255                ENDIF 
    43354256             ENDDO 
     
    43664287          wtd_ns(ji,jst) = undef_sechiba ! in meters 
    43674288          jsl=nslm 
    4368           DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(njsc(ji))) .AND. (jsl > 1) ) 
     4289          DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(ji)) .AND. (jsl > 1) ) 
    43694290             wtd_ns(ji,jst) = zz(jsl)/mille ! in meters 
    43704291             jsl=jsl-1 
     
    43754296       !      This routine does not change tmc but decides where we should turn off ET to prevent further mc decrease 
    43764297       !      Like above, the tests are made on total mc, compared to mcr 
    4377        CALL hydrol_soil_smooth_under_mcr(kjpindex, jst, njsc, is_under_mcr, check_under_ns) 
     4298       CALL hydrol_soil_smooth_under_mcr(mcr, kjpindex, jst, njsc, is_under_mcr, check_under_ns) 
    43784299  
    43794300       !! 4. At the end of the prognostic calculations, we recompute important moisture variables 
     
    43994320       DO jsl = 1, nslm 
    44004321          DO ji =1, kjpindex 
    4401              mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + & 
    4402                   (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) ) 
     4322             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 
     4323                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 
    44034324             ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc 
    44044325          ENDDO 
     
    44514372       sm(:,1)  = dz(2) * (trois*mcl(:,1,jst) + mcl(:,2,jst))/huit 
    44524373       smt(:,1) = dz(2) * (trois*mc(:,1,jst) + mc(:,2,jst))/huit 
    4453        smw(:,1) = dz(2) * (quatre*mcw(njsc(:)))/huit 
    4454        smf(:,1) = dz(2) * (quatre*mcfc(njsc(:)))/huit 
    4455        sms(:,1) = dz(2) * (quatre*mcs(njsc(:)))/huit 
     4374       smw(:,1) = dz(2) * (quatre*mcw(:))/huit 
     4375       smf(:,1) = dz(2) * (quatre*mcfc(:))/huit 
     4376       sms(:,1) = dz(2) * (quatre*mcs(:))/huit 
    44564377       DO jsl = 2,nslm-1 
    44574378          sm(:,jsl)  = dz(jsl) * (trois*mcl(:,jsl,jst)+mcl(:,jsl-1,jst))/huit & 
     
    44594380          smt(:,jsl) = dz(jsl) * (trois*mc(:,jsl,jst)+mc(:,jsl-1,jst))/huit & 
    44604381               + dz(jsl+1) * (trois*mc(:,jsl,jst)+mc(:,jsl+1,jst))/huit 
    4461           smw(:,jsl) = dz(jsl) * ( quatre*mcw(njsc(:)) )/huit & 
    4462                + dz(jsl+1) * ( quatre*mcw(njsc(:)) )/huit 
    4463           smf(:,jsl) = dz(jsl) * ( quatre*mcfc(njsc(:)) )/huit & 
    4464                + dz(jsl+1) * ( quatre*mcfc(njsc(:)) )/huit 
    4465           sms(:,jsl) = dz(jsl) * ( quatre*mcs(njsc(:)) )/huit & 
    4466                + dz(jsl+1) * ( quatre*mcs(njsc(:)) )/huit 
     4382          smw(:,jsl) = dz(jsl) * ( quatre*mcw(:) )/huit & 
     4383               + dz(jsl+1) * ( quatre*mcw(:) )/huit 
     4384          smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit & 
     4385               + dz(jsl+1) * ( quatre*mcfc(:) )/huit 
     4386          sms(:,jsl) = dz(jsl) * ( quatre*mcs(:) )/huit & 
     4387               + dz(jsl+1) * ( quatre*mcs(:) )/huit 
    44674388       ENDDO 
    44684389       sm(:,nslm)  = dz(nslm) * (trois*mcl(:,nslm,jst) + mcl(:,nslm-1,jst))/huit   
    44694390       smt(:,nslm) = dz(nslm) * (trois*mc(:,nslm,jst) + mc(:,nslm-1,jst))/huit       
    4470        smw(:,nslm) = dz(nslm) * (quatre*mcw(njsc(:)))/huit 
    4471        smf(:,nslm) = dz(nslm) * (quatre*mcfc(njsc(:)))/huit 
    4472        sms(:,nslm) = dz(nslm) * (quatre*mcs(njsc(:)))/huit 
     4391       smw(:,nslm) = dz(nslm) * (quatre*mcw(:))/huit 
     4392       smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit 
     4393       sms(:,nslm) = dz(nslm) * (quatre*mcs(:))/huit 
    44734394       ! sm_nostress = soil moisture of each layer at which us reaches 1, here at the middle of [smw,smf] 
    44744395       DO jsl = 1,nslm 
    4475           sm_nostress(:,jsl) = smw(:,jsl) + pcent(njsc(:)) * (smf(:,jsl)-smw(:,jsl)) 
     4396          sm_nostress(:,jsl) = smw(:,jsl) + pcent(:) * (smf(:,jsl)-smw(:,jsl)) 
    44764397       END DO 
    44774398 
     
    46424563       !! The multiplication by vegtot creates grid-cell average values 
    46434564       ! *** To be checked for consistency with the use of nobio properties in thermosoil 
     4565            
    46444566       DO jsl=1,nslm 
    46454567          DO ji=1,kjpindex 
     
    46614583       DO jsl = 1, nslm 
    46624584          ksat(:,jsl) = ksat(:,jsl) + soiltile(:,jst) * & 
    4663                ( ks(njsc(:)) * kfact(jsl,njsc(:)) * kfact_root(:,jsl,jst) )  
     4585               ( ks(:) * kfact(jsl,:) * kfact_root(:,jsl,jst) )  
    46644586       ENDDO 
    46654587               
     
    46714593 
    46724594    !! 7. Summing 3d variables into 2d variables  
    4673     CALL hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 
     4595    CALL hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 
    46744596         & evapot, vevapnu, returnflow, reinfiltration, irrigation, & 
    46754597         & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt) 
     
    47624684       !! 8.2 Since we estimate bare soile evap for the next time step, we update profil_froz_hydro and mcl 
    47634685       !     (effect of mc only, the change in temp_hydro is neglected) 
    4764        IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(kjpindex,jst,njsc) 
     4686       IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(nvan, avan, mcr, mcs,kjpindex,jst,njsc) 
    47654687        DO jsl = 1, nslm 
    47664688          DO ji =1, kjpindex 
    4767              mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + & 
    4768                   (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) ) 
     4689             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + & 
     4690                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) ) 
    47694691             ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc 
    47704692          ENDDO 
     
    47724694 
    47734695       !! 8.3 K and D are recomputed for the updated profile of mc/mcl 
    4774        CALL hydrol_soil_coef(kjpindex,jst,njsc) 
     4696       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc) 
    47754697 
    47764698       !! 8.4 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme 
     
    48114733              
    48124734             flux_top(ji) = zero 
     4735             r_soil_ns(ji,jst) = zero 
    48134736              
    48144737          ENDIF 
     
    48634786       DO ji = 1, kjpindex 
    48644787          ! by construction, mc and mcl are always on the same side of mcr, so we can use mcl here 
    4865           resolv(ji) = (mcl(ji,1,jst).LT.(mcr(njsc(ji))).AND.flux_top(ji).GT.min_sechiba) 
     4788          resolv(ji) = (mcl(ji,1,jst).LT.(mcr(ji)).AND.flux_top(ji).GT.min_sechiba) 
    48664789       ENDDO 
    48674790       !! Reset the coefficient for diffusion (tridiag is only solved if resolv(ji) = .TRUE.)O 
     
    48794802          tmat(ji,1,2) = un 
    48804803          tmat(ji,1,3) = zero 
    4881           rhs(ji,1) = mcr(njsc(ji)) 
     4804          rhs(ji,1) = mcr(ji) 
    48824805       ENDDO 
    48834806        
     
    48984821          DO ji =1, kjpindex 
    48994822             mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + & 
    4900                   profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) ) 
     4823                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) ) 
    49014824             ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl 
    49024825          ENDDO 
     
    50504973!_ hydrol_soil_infilt 
    50514974 
    5052   SUBROUTINE hydrol_soil_infilt(kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check) 
     4975  SUBROUTINE hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check) 
    50534976 
    50544977    !! 0. Variable and parameter declaration 
     
    50574980 
    50584981    ! GLOBAL (in or inout) 
     4982    ! salma added input soil hydraulic parameters 
     4983    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     4984    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless) 
     4985    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1}) 
     4986    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     4987    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     4988    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     4989    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     4990 
     4991 
    50594992    INTEGER(i_std), INTENT(in)                        :: kjpindex        !! Domain size 
    50604993    INTEGER(i_std), INTENT(in)                        :: ins 
     
    51085041    DO ji = 1, kjpindex 
    51095042       !! First we fill up the first layer (about 1mm) without any resistance and quasi-immediately 
    5110        wat_inf_pot(ji) = MAX((mcs(njsc(ji))-mc(ji,1,ins)) * dz(2) / deux, zero) 
     5043       wat_inf_pot(ji) = MAX((mcs(ji)-mc(ji,1,ins)) * dz(2) / deux, zero) 
    51115044       wat_inf(ji) = MIN(wat_inf_pot(ji), flux_infilt(ji)) 
    51125045       mc(ji,1,ins) = mc(ji,1,ins) + wat_inf(ji) * deux / dz(2) 
     
    51275060          ! This is computed by an simple arithmetic average because  
    51285061          ! the time step (30min) is not appropriate for a geometric average (advised by Haverkamp and Vauclin) 
    5129           k_m = (k(ji,jsl) + ks(njsc(ji))*kfact(jsl-1,njsc(ji))*kfact_root(ji,jsl,ins)) / deux  
     5062          k_m = (k(ji,jsl) + ks(ji)*kfact(jsl-1,ji)*kfact_root(ji,jsl,ins)) / deux  
    51305063 
    51315064          IF (ok_freeze_cwrr) THEN 
     
    51415074 
    51425075          !! From which we deduce the time it takes to fill up the layer or to end the time step... 
    5143           wat_inf_pot(ji) =  MAX((mcs(njsc(ji))-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero) 
     5076          wat_inf_pot(ji) =  MAX((mcs(ji)-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero) 
    51445077          IF ( infilt_tmp(ji) > min_sechiba) THEN 
    51455078             dt_inf(ji) =  MIN(wat_inf_pot(ji)/infilt_tmp(ji), dt_tmp(ji)) 
     
    52235156!_ hydrol_soil_smooth_under_mcr 
    52245157 
    5225   SUBROUTINE hydrol_soil_smooth_under_mcr(kjpindex, ins, njsc, is_under_mcr, check) 
     5158  SUBROUTINE hydrol_soil_smooth_under_mcr(mcr, kjpindex, ins, njsc, is_under_mcr, check) 
    52265159 
    52275160    !- arguments 
     
    52315164    !! 0.1 Input variables 
    52325165 
     5166    !salma: added parameter mcr 
     5167    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr          !! Residual volumetric water content (m^{3} m^{-3}) 
    52335168    INTEGER(i_std), INTENT(in)                         :: kjpindex     !! Domain size 
    52345169    INTEGER(i_std), INTENT(in)                         :: ins          !! Soiltile index (1-nstm, unitless) 
     
    52685203    DO jsl = 1,nslm-2 
    52695204       DO ji=1, kjpindex 
    5270           excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero) 
     5205          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 
    52715206          mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 
    52725207          mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * & 
     
    52775212    jsl = nslm-1 
    52785213    DO ji=1, kjpindex 
    5279        excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero) 
     5214       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 
    52805215       mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 
    52815216       mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * & 
     
    52855220    jsl = nslm 
    52865221    DO ji=1, kjpindex 
    5287        excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero) 
     5222       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 
    52885223       mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 
    52895224       mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * & 
     
    52945229    DO jsl = nslm-1,2,-1 
    52955230       DO ji=1, kjpindex 
    5296           excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero) 
     5231          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero) 
    52975232          mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess 
    52985233          mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * & 
     
    53045239    ! excess > 0 
    53055240    DO ji=1, kjpindex 
    5306        excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(njsc(ji))-mc(ji,1,ins),zero) 
     5241       excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(ji)-mc(ji,1,ins),zero) 
    53075242    ENDDO 
    53085243    DO ji=1, kjpindex 
     
    53805315!_ hydrol_soil_smooth_over_mcs 
    53815316 
    5382   SUBROUTINE hydrol_soil_smooth_over_mcs(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 
     5317  SUBROUTINE hydrol_soil_smooth_over_mcs(mcs ,kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 
    53835318 
    53845319    !- arguments 
     
    53865321    !! 0. Variable and parameter declaration 
    53875322 
    5388     !! 0.1 Input variables 
    5389  
     5323    !! 0.1 Input variables  
     5324    !salma: added mcs 
     5325    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3}) 
    53905326    INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size 
    53915327    INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless) 
     
    54255361    DO jsl = 1, nslm-2 
    54265362       DO ji=1, kjpindex 
    5427           excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) 
     5363          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 
    54285364          mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 
    54295365          mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * & 
     
    54345370    jsl = nslm-1 
    54355371    DO ji=1, kjpindex 
    5436        excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) 
     5372       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 
    54375373       mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 
    54385374       mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * & 
     
    54425378    jsl = nslm 
    54435379    DO ji=1, kjpindex 
    5444        excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) 
     5380       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 
    54455381       mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 
    54465382       mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * & 
     
    54515387    DO jsl = nslm-1,2,-1 
    54525388       DO ji=1, kjpindex 
    5453           excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) 
     5389          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero) 
    54545390          mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess 
    54555391          mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * & 
     
    54615397 
    54625398    DO ji=1, kjpindex 
    5463        excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(njsc(ji)),zero) 
     5399       excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(ji),zero) 
    54645400       mc(ji,1,ins) = mc(ji,1,ins) - excess ! then mc(1)=mcs 
    54655401       rudr_corr(ji,ins) = rudr_corr(ji,ins) + excess * dz(2) / deux  
     
    55125448!_ hydrol_soil_smooth_over_mcs2 
    55135449 
    5514   SUBROUTINE hydrol_soil_smooth_over_mcs2(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 
     5450  SUBROUTINE hydrol_soil_smooth_over_mcs2(mcs, kjpindex, ins, njsc, is_over_mcs, rudr_corr, check) 
    55155451 
    55165452    !- arguments 
     
    55195455 
    55205456    !! 0.1 Input variables 
    5521  
     5457    !salma: added mcs 
     5458    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3}) 
    55225459    INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size 
    55235460    INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless) 
     
    55625499    DO jsl = 1, nslm 
    55635500       DO ji=1, kjpindex 
    5564           excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) ! >=0 
     5501          excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(ji),zero) ! >=0 
    55655502          mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess(ji,jsl) ! here mc either does not change or decreases 
    55665503       ENDDO 
     
    57965733!_ ================================================================================================================================ 
    57975734!_ hydrol_soil_coef 
    5798   
    5799   SUBROUTINE hydrol_soil_coef(kjpindex,ins,njsc) 
     5735 
     5736  SUBROUTINE hydrol_soil_coef(mcr, mcs, kjpindex,ins,njsc) 
    58005737 
    58015738    IMPLICIT NONE 
     
    58045741 
    58055742    !! 0.1 Input variables 
    5806  
     5743    !salma: added mcr and mcs 
     5744    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     5745    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
    58075746    INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size 
    58085747    INTEGER(i_std), INTENT(in)                        :: ins              !! Index of soil type 
     
    58365775             ! It corresponds to a liquid mc, but the expression is different from mcl in hydrol_soil, 
    58375776             ! to ensure that we get the a, b, d of the first bin when mcl<mcr 
    5838              mc_used = mcr(njsc(ji))+x*MAX((mc(ji,jsl, ins)-mcr(njsc(ji))),zero)  
     5777             mc_used = mcr(ji)+x*MAX((mc(ji,jsl, ins)-mcr(ji)),zero)  
    58395778             ! 
    58405779             ! calcul de k based on mc_liq 
    58415780             ! 
    5842              i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(njsc(ji)))/(mcs(njsc(ji))-mcr(njsc(ji)))))) 
    5843              a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d 
    5844              b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d 
    5845              d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d 
    5846              k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), & 
    5847                   a_lin(i,jsl,njsc(ji)) * mc_used + b_lin(i,jsl,njsc(ji))) ! in mm/d 
     5781             i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(ji))/(mcs(ji)-mcr(ji))))) 
     5782             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 
     5783             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 
     5784             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d 
     5785             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), & 
     5786                  a_lin(i,jsl,ji) * mc_used + b_lin(i,jsl,ji)) ! in mm/d 
    58485787          ENDDO ! loop on grid 
    58495788       ENDDO 
     
    58555794              
    58565795             ! it is impossible to consider a mc<mcr for the binning 
    5857              mc_ratio = MAX(mc(ji,jsl,ins)-mcr(njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji))) 
     5796             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(ji), zero)/(mcs(ji)-mcr(ji)) 
    58585797              
    58595798             i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin) 
    5860              a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d 
    5861              b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d 
    5862              d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d 
    5863              k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), & 
    5864                   a_lin(i,jsl,njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji)))  ! in mm/d 
     5799             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 
     5800             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d 
     5801             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d 
     5802             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), & 
     5803                  a_lin(i,jsl,ji) * mc(ji,jsl,ins) + b_lin(i,jsl,ji))  ! in mm/d 
    58655804          END DO  
    58665805       END DO 
     
    58875826!_ hydrol_soil_froz 
    58885827  
    5889   SUBROUTINE hydrol_soil_froz(kjpindex,ins,njsc) 
     5828  SUBROUTINE hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,ins,njsc) 
    58905829 
    58915830    IMPLICIT NONE 
     
    58945833 
    58955834    !! 0.1 Input variables 
     5835    !salma: added the following soil variables 
     5836    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: nvan             !! Van Genuchten coeficients n (unitless) 
     5837    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: avan             !! Van Genuchten coeficients a (mm-1}) 
     5838    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     5839    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
    58965840 
    58975841    INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size 
     
    59255869          DO ji=1,kjpindex 
    59265870             ! Van Genuchten parameter for thermodynamical calculation 
    5927              m = 1. -1./nvan(njsc(ji)) 
     5871             m = 1. -1./nvan(ji) 
    59285872            
    5929              IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(njsc(ji))+min_sechiba))) THEN 
     5873             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(ji)+min_sechiba))) THEN 
    59305874                ! Linear soil freezing or soil moisture below residual 
    59315875                IF (temp_hydro(ji, jsl).GE.(ZeroCelsius+fr_dT/2.)) THEN 
     
    59445888                     (temp_hydro(ji,jsl) .LT. (ZeroCelsius+fr_dT/2.)) ) THEN 
    59455889                   ! Factor 2.2 from the PhD of Isabelle Gouttevin 
    5946                    x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) & 
    5947                         *((2.2*1000.*avan(njsc(ji))*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) & 
    5948                         *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m)) / & 
    5949                         (mc(ji,jsl, ins)-mcr(njsc(ji))),1._r_std)                 
     5890                   x=MIN(((mcs(ji)-mcr(ji)) & 
     5891                        *((2.2*1000.*avan(ji)*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) & 
     5892                        *lhf/ZeroCelsius/10.)**nvan(ji)+1.)**(-m)) / & 
     5893                        (mc(ji,jsl, ins)-mcr(ji)),1._r_std)                 
    59505894                ELSE 
    59515895                   x=0._r_std  
     
    59555899             profil_froz_hydro_ns(ji,jsl,ins) = 1._r_std-x 
    59565900              
    5957              mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(njsc(ji)) 
     5901             mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(ji) 
    59585902 
    59595903          ENDDO ! loop on grid 
     
    63976341!_ hydrol_diag_soil 
    63986342 
    6399   SUBROUTINE hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 
     6343  SUBROUTINE hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, & 
    64006344       & evapot, vevapnu, returnflow, reinfiltration, irrigation, & 
    64016345       & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt) 
     
    64066350 
    64076351    !! 0.1 Input variables 
     6352    !salma: added the following soil variables 
     6353    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1}) 
     6354    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless) 
     6355    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1}) 
     6356    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3}) 
     6357    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3}) 
     6358    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     6359    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3}) 
    64086360 
    64096361    ! input scalar  
     
    65586510             i= MAX(MIN(INT((imax-imin)*tmc_litter_ratio)+imin, imax-1), imin) 
    65596511          ENDIF        
    6560           k_tmp = MAX(k_lin(i,1,njsc(ji))*ks(njsc(ji)), zero) 
     6512          k_tmp = MAX(k_lin(i,1,ji)*ks(ji), zero) 
    65616513          k_litt(ji) = k_litt(ji) + vegtot(ji)*soiltile(ji,jst) * SQRT(k_tmp) ! grid-cell average 
    65626514       ENDDO       
     
    66306582          DO ji=1,kjpindex 
    66316583             shumdiag(ji,jsl) = shumdiag(ji,jsl) + soil_wet_ns(ji,jsl,jst) * soiltile(ji,jst) * & 
    6632                                ((mcs(njsc(ji))-mcw(njsc(ji)))/(mcfc(njsc(ji))-mcw(njsc(ji)))) 
     6584                               ((mcs(ji)-mcw(ji))/(mcfc(ji)-mcw(ji))) 
    66336585             shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero)  
    66346586          ENDDO 
     
    66416593    DO jsl=1,nslm               
    66426594       DO ji=1,kjpindex 
    6643           shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji))) 
     6595          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji)) 
    66446596          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero)  
    66456597       ENDDO 
  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90

    r6319 r6954  
    1818!! which is called before the modules which calculate them. 
    1919!!  
    20 !! RECENT CHANGE(S): None  
     20!! RECENT CHANGE(S): November 2020: It is possible to define soil hydraulic parameters from maps, 
     21!!                    as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 
     22!!                    Here, it leads to declare and allocate global variables.  
    2123!!  
    2224!! REFERENCE(S) : None 
     
    113115  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration)  
    114116!$OMP THREADPRIVATE(reinf_slope) 
     117 
     118 
     119!salma: adding soil hydraulic params 
     120  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: ks              !! Saturated soil conductivity (mm d^{-1}) 
     121!$OMP THREADPRIVATE(ks) 
     122  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: nvan            !! Van Genushten n parameter (unitless) 
     123!$OMP THREADPRIVATE(nvan) 
     124  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: avan            !! Van Genushten alpha parameter (mm ^{-1}) 
     125!$OMP THREADPRIVATE(avan) 
     126  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcr             !! Residual soil moisture (m^{3} m^{-3}) 
     127!$OMP THREADPRIVATE(mcr) 
     128  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcs             !! Saturated soil moisture (m^{3} m^{-3}) 
     129!$OMP THREADPRIVATE(mcs) 
     130  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcfc            !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     131!$OMP THREADPRIVATE(mcfc) 
     132  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcw             !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     133!$OMP THREADPRIVATE(mcw) 
     134!end salma:adding soil hydraulic params 
     135 
     136 
    115137  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag       !! Mean relative soil moisture in the different levels used  
    116138                                                                     !! by thermosoil.f90 (unitless, 0-1) 
     
    431453                              index,         indexveg,     lalo,           neighbours,        & 
    432454                              resolution,    contfrac,     temp_air,                          & 
    433                               soiltile,      reinf_slope,  deadleaf_cover, assim_param,       & 
     455                              soiltile,      reinf_slope,  ks,             nvan,              & 
     456                              avan,          mcr,          mcs,            mcfc,              & 
     457                              mcw,           deadleaf_cover,               assim_param,       & 
    434458                              lai,           frac_age,     height,         veget,             & 
    435459                              frac_nobio,    njsc,         veget_max,      fraclut,           & 
     
    451475     
    452476    !! 1.7 Initialize remaining hydrological variables 
    453     CALL hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          & 
     477    CALL hydrol_initialize (ks,  nvan, avan, mcr, mcs, mcfc, mcw,    & 
     478         kjit,           kjpindex,  index,         rest_id,          & 
    454479         njsc,           soiltile,  veget,         veget_max,        & 
    455480         humrel,         vegstress, drysoil_frac,                    & 
     
    718743    !! 4. Compute hydrology 
    719744    !! 4.1 Water balance from CWRR module (11 soil layers) 
    720     CALL hydrol_main (kjit, kjpindex, & 
     745    CALL hydrol_main (ks,  nvan, avan, mcr, mcs, mcfc, mcw, kjit, kjpindex, & 
    721746         & index, indexveg, indexsoil, indexlayer, indexnslm, & 
    722747         & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, & 
     
    13731398                            njsc,       lai,       height,   veget,      & 
    13741399                            frac_nobio, veget_max, reinf_slope,          & 
     1400                            ks,         nvan,      avan,     mcr,        & 
     1401                            mcs,        mcfc,      mcw,                  & 
    13751402                            co2_to_bm,  assim_param, frac_age) 
    13761403     
     
    15301557    ALLOCATE (reinf_slope(kjpindex),stat=ier) 
    15311558    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','') 
     1559 
     1560    !salma: adding soil params 
     1561    ALLOCATE (ks(kjpindex),stat=ier) 
     1562    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ks','','') 
     1563 
     1564    ALLOCATE (nvan(kjpindex),stat=ier) 
     1565    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nvan ','','') 
     1566 
     1567    ALLOCATE (avan(kjpindex),stat=ier) 
     1568    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for avan','','') 
     1569 
     1570    ALLOCATE (mcr(kjpindex),stat=ier) 
     1571    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcr','','') 
     1572 
     1573    ALLOCATE (mcs(kjpindex),stat=ier) 
     1574    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcs','','') 
     1575 
     1576    ALLOCATE (mcfc(kjpindex),stat=ier) 
     1577    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcfc','','') 
     1578     
     1579    ALLOCATE (mcw(kjpindex),stat=ier) 
     1580    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcw','','') 
     1581    !end salma: adding soil params 
     1582 
     1583 
     1584 
    15321585 
    15331586    ALLOCATE (vbeta1(kjpindex),stat=ier) 
     
    18771930    IF ( ALLOCATED (njsc)) DEALLOCATE (njsc) 
    18781931    IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope) 
     1932 
     1933    !salma: adding soil hydraulic params 
     1934    IF ( ALLOCATED (ks)) DEALLOCATE (ks) 
     1935    IF ( ALLOCATED (nvan)) DEALLOCATE (nvan) 
     1936    IF ( ALLOCATED (avan)) DEALLOCATE (avan) 
     1937    IF ( ALLOCATED (mcr)) DEALLOCATE (mcr) 
     1938    IF ( ALLOCATED (mcs)) DEALLOCATE (mcs) 
     1939    IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc) 
     1940    IF ( ALLOCATED (mcw)) DEALLOCATE (mcw) 
     1941    !end salma: adding soil hydraulic params 
     1942 
    18791943    IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1) 
    18801944    IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4) 
  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

    r6319 r6954  
    1717!! 
    1818!! RECENT CHANGE(S): Allowed reading of USDA map, Nov 2014, ADucharne 
     19!!                   November 2020: It is possible to define soil hydraulic parameters from maps, 
     20!!                   as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 
     21!!                   Changes in slowproc_xios_initialize, slowproc_soilt and slowproc_finalize 
    1922!! 
    2023!! REFERENCE(S) : 
    21 !! 
     24!!- Tafasca S. (2020). Evaluation de l’impact des propriétés du sol sur l’hydrologie simulee dans le 
     25!! modÚle ORCHIDEE, PhD thesis, Sorbonne Universite. \n 
    2226!! SVN          : 
    2327!! $HeadURL$ 
     
    9296!! 
    9397!! DESCRIPTION:   Initialize xios dependant defintion before closing context defintion 
     98!!   
     99!! RECENT CHANGE(S): Initialization of XIOS to read soil hydraulic parameters from maps, 
     100!!                   as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes). 
    94101!! 
    95102!! \n 
     
    219226    ENDIF 
    220227 
     228    !Salma 
     229    !! 6. Prepare for reading of soil parameter files 
     230 
     231    ! Get the file name from run.def file and set file attributes accordingly 
     232    filename = 'params_sp_mip.nc' 
     233    CALL getin_p('PARAM_FILE',filename) 
     234    name = filename(1:LEN_TRIM(FILENAME)-3) 
     235    CALL xios_orchidee_set_file_attr("soilparam_file",name=name) 
     236    ! Determine if the file will be read by XIOS. If not, deactivate reading of the file. 
     237    IF (xios_interpolation .AND. restname_in=='NONE' .AND. .NOT. impsoilt) THEN 
     238       ! Reading will be done with XIOS later 
     239       IF (printlev>=2) WRITE(numout,*) 'Reading of soil hydraulic parameters file will be done later using XIOS. The filename is ', filename 
     240    ELSE 
     241       ! No reading by XIOS, deactivate soilparam_file and related variables declared in context_input_orchidee.xml. 
     242       ! If this is not done, the model will crash if the file is not available in the run directory. 
     243       IF (printlev>=2) WRITE(numout,*) 'Reading of soil parameter file will not be done with XIOS.' 
     244       CALL xios_orchidee_set_file_attr("soilparam_file",enabled=.FALSE.) 
     245       CALL xios_orchidee_set_field_attr("soilks",enabled=.FALSE.) 
     246       CALL xios_orchidee_set_field_attr("soilnvan",enabled=.FALSE.) 
     247       CALL xios_orchidee_set_field_attr("soilavan",enabled=.FALSE.) 
     248       CALL xios_orchidee_set_field_attr("soilmcr",enabled=.FALSE.) 
     249       CALL xios_orchidee_set_field_attr("soilmcs",enabled=.FALSE.) 
     250       CALL xios_orchidee_set_field_attr("soilmcfc",enabled=.FALSE.) 
     251       CALL xios_orchidee_set_field_attr("soilmcw",enabled=.FALSE.) 
     252    ENDIF 
     253 
    221254    IF (printlev_loc>=3) WRITE(numout,*) 'End slowproc_xios_intialize' 
    222255    
     
    244277                                  IndexLand,     indexveg,     lalo,           neighbours,        & 
    245278                                  resolution,    contfrac,     temp_air,                          & 
    246                                   soiltile,      reinf_slope,  deadleaf_cover, assim_param,       & 
     279                                  soiltile,      reinf_slope,  ks,             nvan,              & 
     280                                  avan,          mcr,          mcs,            mcfc,              & 
     281                                  mcw,           deadleaf_cover,               assim_param,       & 
    247282                                  lai,           frac_age,     height,         veget,             & 
    248283                                  frac_nobio,    njsc,         veget_max,      fraclut,           & 
     
    282317    REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out)    :: nwdFraclut     !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless) 
    283318    REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: reinf_slope    !! slope coef for reinfiltration 
     319    !Salma: adding soil params 
     320    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1}) 
     321    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless) 
     322    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1}) 
     323    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3}) 
     324    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3}) 
     325    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     326    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     327 
    284328    REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (out):: assim_param    !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1}) 
    285329    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless) 
     
    296340    CALL slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & 
    297341         rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, & 
     342         ks,  nvan, avan, mcr, mcs, mcfc, mcw, & 
    298343         veget_max, tot_bare_soil, njsc, & 
    299344         height, lcanop, veget_update, veget_year) 
     
    684729!>\BRIEF         Write to restart file variables for slowproc module and call finalization of stomate module 
    685730!! 
    686 !! DESCRIPTION :  
     731!! DESCRIPTION : 
     732!! 
     733!! RECENT CHANGE(S): Add arrays of soil hydraulic parameters to the restart file. 
     734!!                   Linked to SP-MIP project (Tafasca Salma and Ducharne Agnes). 
    687735!! 
    688736!! MAIN OUTPUT VARIABLE(S) :  
     
    697745                                njsc,       lai,       height,   veget,      & 
    698746                                frac_nobio, veget_max, reinf_slope,          & 
     747                                ks,  nvan, avan, mcr, mcs, mcfc, mcw,        & 
    699748                                co2_to_bm,  assim_param, frac_age ) 
    700749 
     
    711760    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget_max      !! Maximum fraction of vegetation type including none biological fraction (unitless) 
    712761    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: reinf_slope    !! slope coef for reinfiltration 
     762    !salma: added following soil params 
     763    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: ks             !! Hydraulic conductivity at saturation (mm {-1}) 
     764    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: nvan           !! Van Genuchten coeficients n (unitless) 
     765    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: avan           !! Van Genuchten coeficients a (mm-1}) 
     766    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3}) 
     767    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3}) 
     768    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     769    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     770 
    713771    REAL(r_std),DIMENSION (kjpindex,nvm),INTENT(in)      :: co2_to_bm      !! virtual gpp flux between atmosphere and biosphere 
    714772    REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param   !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1}) 
     
    759817    CALL restput_p (rest_id, 'clay_frac', nbp_glo, 1, 1, kjit, clayfraction, 'scatter',  nbp_glo, index_g) 
    760818    CALL restput_p (rest_id, 'sand_frac', nbp_glo, 1, 1, kjit, sandfraction, 'scatter',  nbp_glo, index_g) 
    761     ! 
     819    !salma: added the following lines for restput of the soil parameters 
     820    CALL restput_p (rest_id, 'ks', nbp_glo, 1, 1, kjit, ks, 'scatter',  nbp_glo, index_g) 
     821    CALL restput_p (rest_id, 'mcs', nbp_glo, 1, 1, kjit, mcs, 'scatter',  nbp_glo, index_g) 
     822    CALL restput_p (rest_id, 'mcr', nbp_glo, 1, 1, kjit, mcr, 'scatter',  nbp_glo, index_g) 
     823    CALL restput_p (rest_id, 'mcw', nbp_glo, 1, 1, kjit, mcw, 'scatter',  nbp_glo, index_g) 
     824    CALL restput_p (rest_id, 'mcfc', nbp_glo, 1, 1, kjit, mcfc, 'scatter',  nbp_glo, index_g) 
     825    CALL restput_p (rest_id, 'nvan', nbp_glo, 1, 1, kjit, nvan, 'scatter',  nbp_glo, index_g) 
     826    CALL restput_p (rest_id, 'avan', nbp_glo, 1, 1, kjit, avan, 'scatter',  nbp_glo, index_g) 
     827 
    762828    ! The height of the vegetation could in principle be recalculated at the beginning of the run. 
    763829    ! However, this is very tedious, as many special cases have to be taken into account. This variable 
     
    803869  SUBROUTINE slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, & 
    804870       rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, & 
     871       ks,  nvan, avan, mcr, mcs, mcfc, mcw, & 
    805872       veget_max, tot_bare_soil, njsc, & 
    806873       height, lcanop, veget_update, veget_year) 
     
    837904    REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out)   :: nwdfraclut     !! Fraction of non woody vegetation in each landuse tile 
    838905    REAL(r_std), DIMENSION (kjpindex), INTENT(out)        :: reinf_slope    !! slope coef for reinfiltration 
     906 
     907    !salma: added soil parameters     
     908    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1}) 
     909    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless) 
     910    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1}) 
     911    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3}) 
     912    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3}) 
     913    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     914    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     915 
    839916    INTEGER(i_std), DIMENSION(kjpindex), INTENT(out)      :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 
    840917     
     
    10151092       njsc = NINT(tmp_real) 
    10161093    END IF 
    1017      
     1094 
     1095    ! Salma: we define restget for the 7 soil variables following clay_frac 
     1096 
     1097    var_name= 'ks' 
     1098    CALL ioconf_setatt_p('UNITS', 'mm/d') 
     1099    CALL ioconf_setatt_p('LONG_NAME','Soil saturated water content') 
     1100    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., ks, "gather", nbp_glo, index_g) 
     1101 
     1102    var_name= 'mcs' 
     1103    CALL ioconf_setatt_p('UNITS', 'm3/m3') 
     1104    CALL ioconf_setatt_p('LONG_NAME','') 
     1105    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcs, "gather", nbp_glo, index_g) 
     1106 
     1107    var_name= 'mcr' 
     1108    CALL ioconf_setatt_p('UNITS', 'm3/m3') 
     1109    CALL ioconf_setatt_p('LONG_NAME','') 
     1110    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcr, "gather", nbp_glo, index_g) 
     1111 
     1112    var_name= 'mcfc' 
     1113    CALL ioconf_setatt_p('UNITS', 'm3/m3') 
     1114    CALL ioconf_setatt_p('LONG_NAME','') 
     1115    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcfc, "gather", nbp_glo, index_g) 
     1116 
     1117    var_name= 'mcw' 
     1118    CALL ioconf_setatt_p('UNITS', 'm3/m3') 
     1119    CALL ioconf_setatt_p('LONG_NAME','') 
     1120    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcw, "gather", nbp_glo, index_g) 
     1121 
     1122    var_name= 'nvan' 
     1123    CALL ioconf_setatt_p('UNITS', '-') 
     1124    CALL ioconf_setatt_p('LONG_NAME','') 
     1125    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., nvan, "gather", nbp_glo, index_g) 
     1126 
     1127    var_name= 'avan' 
     1128    CALL ioconf_setatt_p('UNITS', 'm-1') 
     1129    CALL ioconf_setatt_p('LONG_NAME','') 
     1130    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., avan, "gather", nbp_glo, index_g) 
     1131 
     1132    !end salma 
     1133 
    10181134    var_name= 'clay_frac' 
    10191135    CALL ioconf_setatt_p('UNITS', '-') 
     
    12061322               MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN 
    12071323              
    1208              CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, & 
     1324             CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, lalo, neighbours, resolution, contfrac, soilclass, & 
    12091325                  clayfraction, sandfraction, siltfraction) 
    12101326             njsc(:) = 0 
     
    13431459        
    13441460       IF (printlev_loc>=4) WRITE (numout,*) 'clayfraction or njcs were not in restart file, call slowproc_soilt' 
    1345        CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, & 
     1461       CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, lalo, neighbours, resolution, contfrac, soilclass, & 
    13461462            clayfraction, sandfraction, siltfraction) 
    13471463       IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_soilt' 
     
    26142730!! 
    26152731!! RECENT CHANGE(S): Nov 2014, ADucharne 
     2732!!                   Nov 2020, Salma Tafasca and Agnes Ducharne: adding a choice for spmipexp/SPMIPEXP, 
     2733!!                             and everything needed to read all maps and assign parameter values.   
    26162734!! 
    26172735!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction 
     
    26242742!! \n 
    26252743!_ ================================================================================================================================ 
    2626   SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction) 
     2744  SUBROUTINE slowproc_soilt(njsc,  ks,  nvan, avan, mcr, mcs, mcfc, mcw, nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction) 
    26272745 
    26282746    USE interpweight 
     
    26492767    !  0.2 OUTPUT 
    26502768    ! 
     2769    !salma: added soil params and njsc because needed in the calculation of the soil params 
     2770    INTEGER(i_std),DIMENSION (nbpt), INTENT (out)      :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless) 
     2771    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1}) 
     2772    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless) 
     2773    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1}) 
     2774    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3}) 
     2775    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3}) 
     2776    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3}) 
     2777    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3}) 
     2778 
    26512779    REAL(r_std), INTENT(out)      :: soilclass(nbpt, nscm)  !! Soil type map to be created from the Zobler map 
    26522780                                                            !! or a map defining the 12 USDA classes (e.g. Reynolds) 
     
    26602788    !  0.3 LOCAL 
    26612789    ! 
     2790    !salma: added the following local variable to be used for all the soil hydraulic parameters 
     2791    REAL(r_std), DIMENSION(nbpt)        :: param            !! to be introduced in function: interpweight 
     2792 
    26622793    CHARACTER(LEN=80) :: filename 
    26632794    INTEGER(i_std) :: ib, ilf, nbexp, i 
     
    26762807    REAL(r_std), DIMENSION(:,:), ALLOCATABLE             :: textrefrac       !! text fractions re-dimensioned 
    26772808    REAL(r_std), DIMENSION(nbpt)                         :: atext            !! Availability of the texture interpolation 
     2809    !salma added the following 3 variables to control the SP-MIP simulations 
     2810    REAL(r_std), DIMENSION(nbpt)                         :: aparam            !! Availability of the parameter interpolation 
     2811    CHARACTER(LEN=80)                                    :: spmipexp          !! designing the number of sp-mip experiment 
     2812    CHARACTER(LEN=80)                                    :: exp               !! designing the model of experiment 4 (sp_mip) 
     2813 
    26782814    REAL(r_std)                                          :: vmin, vmax       !! min/max values to use for the  
    26792815 
     
    27052841 
    27062842    IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt' 
    2707     ! 
    2708     !  Needs to be a configurable variable 
    2709     ! 
    2710     ! 
    2711     !Config Key   = SOILCLASS_FILE 
    2712     !Config Desc  = Name of file from which soil types are read 
    2713     !Config Def   = soils_param.nc 
    2714     !Config If    = NOT(IMPOSE_VEG) 
    2715     !Config Help  = The name of the file to be opened to read the soil types.  
    2716     !Config         The data from this file is then interpolated to the grid of 
    2717     !Config         of the model. The aim is to get fractions for sand loam and 
    2718     !Config         clay in each grid box. This information is used for soil hydrology 
    2719     !Config         and respiration. 
    2720     !Config Units = [FILE] 
    2721     ! 
    2722     ! soils_param.nc file is 1deg soil texture file (Zobler) 
    2723     ! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution) 
    2724  
    2725     filename = 'soils_param.nc' 
    2726     CALL getin_p('SOILCLASS_FILE',filename) 
    2727  
    2728     variablename = 'soiltext' 
    2729  
    2730     !! Variables for interpweight 
    2731     ! Type of calculation of cell fractions 
    2732     fractype = 'default' 
    2733  
    2734     IF (xios_interpolation) THEN 
    2735        IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " & 
    2736             // TRIM(filename) // " for variable " // TRIM(variablename) 
     2843 
     2844    ! Salma: sp_mip_experiment: select experiment:exp1, exp2, exp3 or exp4 in run.def 
     2845    !   exp1: Reading the soil parameter maps of SPMIP 
     2846    !   exp2: Reading a usda soil texture map, with 12 USDA classes + Oxisols 
     2847    !         Works for SP-MIP textur emap based on SoilGrids, Reynolds, and any map with USDA classes 
     2848    !   exp3: Reading the Zobler soil texture map 
     2849    !   exp4: Imposing uniform soil texture over the globe (4 texture options, with parameter values imposed by SP-MIP) 
     2850    ! We read a soil texture map in all experiments but exp4. 
     2851    ! Even with exp1, some parameters are defined based on texture. 
     2852     
     2853    ! IMPORTANT: if no exp is defined in run.def, the model works as before, by deriving the soil parameters  
     2854    ! from a soil texture map, itself defined by the SOILTYPE_CLASSIF keyword, and soil_classif variable: 
     2855    !   if soil_classif = zobler, spmipexp = exp3 
     2856    !   if soil_classif = usda, spmipexp = exp2 
     2857    ! For the SP-MIP experiments made by Salma Tafasca, exp1 was made by reading the SPMIP soil texture map (for clayfraction, etc.) 
     2858    ! But to get a uniform texture (exp 4), you need to select a soil texture map using soil_classif, even if it's not read 
     2859    ! WARNING: if you choose the zobler map with usda/exp2, the simulation will run, but not as you probably wish. 
     2860 
     2861    SELECTCASE (soil_classif) ! copied from src-parameters/control.f90 
     2862    CASE ('zobler','none') 
     2863       spmipexp='exp3' 
     2864    CASE ('usda') 
     2865       spmipexp='exp2' 
     2866    CASE DEFAULT 
     2867       WRITE(numout,*) "Unsupported soil type classification: soil_classif=",soil_classif 
     2868       WRITE(numout,*) "Choose between zobler, usda and none according to the map" 
     2869       CALL ipslerr_p(3,'control_initialize','Bad choice of soil_classif','Choose between zobler, usda and none','') 
     2870    ENDSELECT 
     2871 
     2872    !Config Key   = spmipexp 
     2873    !Config Desc  = number of sp_mip experiment 
     2874    !Config Help  = possible values: exp1, exp2, exp3 and exp4 
     2875    !Config Units = [-] 
     2876    CALL getin_p("SPMIPEXP",spmipexp) 
     2877 
     2878    IF (spmipexp == 'exp4') THEN 
     2879       ! case where exp4 is selected: uniform soil parameters 
     2880       ! the values of the hydraulic parameters below come from SP-MIP, 
     2881       ! and correspond to the Rosetta PTF (Schaap et al., 2001) 
     2882 
     2883       ! sp_mip_experiment_4: select another level of experiment: a, b, c or d in run.def 
     2884 
     2885       !Config Key   = EXP4 
     2886       !Config Desc  = number of sp_mip experiment 4 
     2887       !Config Help  = possible values: a, b, c and d 
     2888       !Config Units = [-] 
     2889       CALL getin_p("EXP4",exp) 
     2890 
     2891       SELECTCASE (exp) 
     2892 
     2893       CASE ('a') ! loamy sand 
     2894          clayfraction=0.06 
     2895          sandfraction=0.81 
     2896          siltfraction=0.13 
     2897         ! njsc=2 
     2898          DO ib=1 , nbpt 
     2899             njsc=2 
     2900             mcr(ib) = 0.049 
     2901             mcs(ib) = 0.39 
     2902             ks(ib) = (1.41e-5)*1000*24*3600 
     2903             avan(ib) = 3.475*(1e-3) 
     2904             nvan(ib) = 1.746 
     2905             mcfc(ib) = 0.1039 
     2906             mcw(ib) = 0.05221 
     2907          ENDDO 
     2908          ! definir clayfraction à partir du barycentre du domaine loamy sand 
     2909          ! definir njsc pour cette classe => pcent dans hydrol 
     2910 
     2911       CASE ('b') !loam 
     2912          clayfraction=0.2 
     2913          sandfraction=0.4 
     2914          siltfraction=0.4 
     2915          njsc=6 
     2916          DO ib=1, nbpt 
     2917             mcr(ib) = 0.061 
     2918             mcs(ib) = 0.399 
     2919             ks(ib) = (3.38e-6)*1000*24*3600 
     2920             avan(ib) = 1.112*(1e-3) 
     2921             nvan(ib) = 1.472 
     2922             mcfc(ib) = 0.236 
     2923             mcw(ib) = 0.09115 
     2924          ENDDO 
     2925          ! definir clayfraction à partir du barycentre du domaine 
     2926          ! definir njsc pour cette classe 
     2927 
     2928       CASE ('c') !silt 
     2929          clayfraction=0.1 
     2930          sandfraction=0.06 
     2931          siltfraction=0.84 
     2932          njsc=5 
     2933          DO ib=1, nbpt 
     2934             mcr(ib) = 0.05 
     2935             mcs(ib) = 0.489 
     2936             ks(ib) = (2.81e-6)*1000*24*3600 
     2937             avan(ib) = 0.6577*(1e-3) 
     2938             nvan(ib) = 1.679 
     2939             mcfc(ib) = 0.2854 
     2940             mcw(ib) = 0.06944 
     2941          ENDDO 
     2942 
     2943          ! definir clayfraction à partir du barycentre du domaine 
     2944          ! definir njsc pour cette classe 
     2945       CASE ('d')!clay 
     2946          clayfraction=0.55 
     2947          sandfraction=0.15 
     2948          siltfraction=0.3 
     2949          njsc=12 
     2950          DO ib=1, nbpt 
     2951             mcr(ib) = 0.098 
     2952             mcs(ib) = 0.459 
     2953             ks(ib) = (9.74e-7)*1000*24*3600 
     2954             avan(ib) = 1.496*(1e-3) 
     2955             nvan(ib) = 1.253 
     2956             mcfc(ib) = 0.3329 
     2957             mcw(ib) = 0.1897 
     2958          ENDDO 
     2959          ! definir clayfraction à partir du barycentre du domaine 
     2960          ! definir njsc pour cette classe 
     2961 
     2962       CASE DEFAULT 
     2963 
     2964          WRITE (numout,*) 'Unsupported experiment number. Choose between a, b, c or d according to sp_mip_experiment_4 number' 
     2965          CALL ipslerr_p(3,'hydrol_init','Unsupported experiment number. ',& 
     2966               'Choose between a,b,c or d','') 
     2967       ENDSELECT 
     2968 
     2969    ELSE !if we are here then it's either exp1, exp2 or exp3 
    27372970        
    2738        SELECT CASE(soil_classif) 
    2739  
    2740        CASE('none') 
    2741           ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
    2742           IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    2743           DO ib=1, nbpt 
    2744              soilclass(ib,:) = soilclass_default_fao 
    2745              clayfraction(ib) = clayfraction_default 
    2746           ENDDO 
    2747  
    2748  
    2749        CASE('zobler') 
    2750  
    2751          ! 
    2752           soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes 
    2753           ! 
    2754           IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification" 
    2755           ! 
    2756           ALLOCATE(textrefrac(nbpt,nzobler)) 
    2757           ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR) 
    2758           IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    2759           CALL get_soilcorr_zobler (nzobler, textfrac_table) 
    2760         
    2761           CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1)) 
    2762           CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2)) 
    2763           CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3)) 
    2764           CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4)) 
    2765           CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5)) 
    2766           CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6)) 
    2767           CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7)) 
    2768  
    2769  
    2770         
    2771           CALL get_soilcorr_zobler (nzobler, textfrac_table) 
    2772         ! 
    2773         ! 
    2774           DO ib =1, nbpt 
    2775             soilclass(ib,1)=textrefrac(ib,1) 
    2776             soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7) 
    2777             soilclass(ib,3)=textrefrac(ib,5) 
    2778              
    2779             ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture) 
    2780             ! over the zobler pixels composing the ORCHIDEE grid-cell 
    2781             clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + & 
    2782                                textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + & 
    2783                                textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7) 
    2784  
    2785             sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + & 
    2786                                textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + & 
    2787                                textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7) 
    2788  
    2789             siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + & 
    2790                                textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + & 
    2791                                textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7) 
    2792  
    2793             sgn=SUM(soilclass(ib,1:3)) 
    2794              
    2795             IF (sgn < min_sechiba) THEN 
    2796               soilclass(ib,:) = soilclass_default(:) 
    2797               clayfraction(ib) = clayfraction_default 
    2798               sandfraction(ib) = sandfraction_default 
    2799               siltfraction(ib) = siltfraction_default 
    2800               atext(ib)=0. 
    2801             ELSE 
    2802               atext(ib)=sgn 
    2803               clayfraction(ib) = clayfraction(ib) / sgn 
    2804               sandfraction(ib) = sandfraction(ib) / sgn 
    2805               siltfraction(ib) = siltfraction(ib) / sgn 
    2806               soilclass(ib,1:3) = soilclass(ib,1:3) / sgn 
    2807             ENDIF 
    2808              
    2809           ENDDO 
    2810            
    2811            
    2812         
    2813        CASE('usda') 
    2814  
    2815            IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda' 
    2816   
    2817            soilclass_default=soilclass_default_usda 
    2818            ! 
    2819            WRITE(numout,*) "Using a soilclass map with usda classification" 
    2820            ! 
    2821            ALLOCATE(textrefrac(nbpt,nscm)) 
    2822            ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
    2823            IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    2824      
    2825            CALL get_soilcorr_usda (nscm, textfrac_table) 
    2826      
    2827            IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda' 
    2828  
    2829           CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1)) 
    2830           CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2)) 
    2831           CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3)) 
    2832           CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4)) 
    2833           CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5)) 
    2834           CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6)) 
    2835           CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7)) 
    2836           CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8)) 
    2837           CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9)) 
    2838           CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10)) 
    2839           CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11)) 
    2840           CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12)) 
    2841  
    2842  
    2843         
    2844           CALL get_soilcorr_usda (nscm, textfrac_table) 
    2845           IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'       
    2846  
    2847           DO ib =1, nbpt 
    2848             clayfraction(ib) = 0.0 
    2849             DO ilf = 1,nscm 
    2850               soilclass(ib,ilf)=textrefrac(ib,ilf)       
    2851               clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf) 
    2852               sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf) 
    2853               siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf) 
    2854             ENDDO 
    2855  
    2856  
    2857             sgn=SUM(soilclass(ib,:)) 
    2858              
    2859             IF (sgn < min_sechiba) THEN 
    2860               soilclass(ib,:) = soilclass_default(:) 
    2861               clayfraction(ib) = clayfraction_default 
    2862               sandfraction(ib) = sandfraction_default 
    2863               siltfraction(ib) = siltfraction_default 
    2864               atext(ib)=0 
    2865             ELSE 
    2866               soilclass(ib,:) = soilclass(ib,:) / sgn 
    2867               clayfraction(ib) = clayfraction(ib) / sgn 
    2868               sandfraction(ib) = sandfraction(ib) / sgn 
    2869               siltfraction(ib) = siltfraction(ib) / sgn 
    2870               atext(ib)=sgn 
    2871             ENDIF 
    2872           ENDDO 
    2873  
    2874         CASE DEFAULT 
    2875              WRITE(numout,*) 'slowproc_soilt:' 
    2876              WRITE(numout,*) '  A non supported soil type classification has been chosen' 
    2877              CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','') 
    2878         END SELECT 
    2879  
    2880  
    2881  
    2882     ELSE              !    xios_interpolation  
    2883        ! Read and interpolate using stardard method with IOIPSL and aggregate 
    2884      
    2885        IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " & 
    2886             // TRIM(filename) // " for variable " // TRIM(variablename) 
    2887  
    2888  
    2889        ! Name of the longitude and latitude in the input file 
    2890        lonname = 'nav_lon' 
    2891        latname = 'nav_lat' 
    2892  
    2893        IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " & 
    2894             // TRIM(filename) // " for variable " // TRIM(variablename) 
    2895  
    2896        IF ( TRIM(soil_classif) /= 'none' ) THEN 
    2897  
    2898           ! Define a variable for the number of soil textures in the input file 
    2899           SELECTCASE(soil_classif) 
     2971       !  Needs to be a configurable variable 
     2972       ! 
     2973       !Config Key   = SOILCLASS_FILE 
     2974       !Config Desc  = Name of file from which soil types are read 
     2975       !Config Def   = soils_param.nc 
     2976       !Config If    = NOT(IMPOSE_VEG) 
     2977       !Config Help  = The name of the file to be opened to read the soil types. 
     2978       !Config         The data from this file is then interpolated to the grid of 
     2979       !Config         of the model. The aim is to get fractions for sand loam and 
     2980       !Config         clay in each grid box. This information is used for soil hydrology 
     2981       !Config         and respiration. 
     2982       !Config Units = [FILE] 
     2983       ! 
     2984       ! soils_param.nc file is 1deg soil texture file (Zobler) 
     2985       ! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution) 
     2986 
     2987       filename = 'soils_param.nc' 
     2988       CALL getin_p('SOILCLASS_FILE',filename) 
     2989 
     2990       variablename = 'soiltext' 
     2991 
     2992       !! Variables for interpweight 
     2993       ! Type of calculation of cell fractions 
     2994       fractype = 'default' 
     2995 
     2996       IF (xios_interpolation) THEN 
     2997          IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " & 
     2998               // TRIM(filename) // " for variable " // TRIM(variablename) 
     2999 
     3000          SELECT CASE(soil_classif) 
     3001 
     3002          CASE('none') 
     3003             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
     3004             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3005             DO ib=1, nbpt 
     3006                soilclass(ib,:) = soilclass_default_fao 
     3007                clayfraction(ib) = clayfraction_default 
     3008             ENDDO 
     3009 
    29003010          CASE('zobler') 
    2901              ntextinfile=nzobler 
     3011             ! 
     3012             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes 
     3013             ! 
     3014             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification, to be read using XIOS" 
     3015             ! 
     3016             ALLOCATE(textrefrac(nbpt,nzobler)) 
     3017             ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR) 
     3018             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3019             CALL get_soilcorr_zobler (nzobler, textfrac_table) 
     3020             CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1)) 
     3021             CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2)) 
     3022             CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3)) 
     3023             CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4)) 
     3024             CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5)) 
     3025             CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6)) 
     3026             CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7)) 
     3027 
     3028             CALL get_soilcorr_zobler (nzobler, textfrac_table) 
     3029             !             ! 
     3030             DO ib =1, nbpt 
     3031                soilclass(ib,1)=textrefrac(ib,1) 
     3032                soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7) 
     3033                soilclass(ib,3)=textrefrac(ib,5) 
     3034 
     3035                ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture) 
     3036                ! over the zobler pixels composing the ORCHIDEE grid-cell 
     3037                clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + & 
     3038                     textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + & 
     3039                     textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7) 
     3040 
     3041                sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + & 
     3042                     textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + & 
     3043                     textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7) 
     3044 
     3045                siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + & 
     3046                     textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + & 
     3047                     textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7) 
     3048 
     3049                sgn=SUM(soilclass(ib,1:3)) 
     3050 
     3051                IF (sgn < min_sechiba) THEN 
     3052                   soilclass(ib,:) = soilclass_default(:) 
     3053                   clayfraction(ib) = clayfraction_default 
     3054                   sandfraction(ib) = sandfraction_default 
     3055                   siltfraction(ib) = siltfraction_default 
     3056                   atext(ib)=0. 
     3057                ELSE 
     3058                   atext(ib)=sgn 
     3059                   clayfraction(ib) = clayfraction(ib) / sgn 
     3060                   sandfraction(ib) = sandfraction(ib) / sgn 
     3061                   siltfraction(ib) = siltfraction(ib) / sgn 
     3062                   soilclass(ib,1:3) = soilclass(ib,1:3) / sgn 
     3063                ENDIF 
     3064 
     3065             ENDDO 
     3066 
    29023067          CASE('usda') 
    2903              ntextinfile=nscm 
     3068 
     3069             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda' 
     3070 
     3071             soilclass_default=soilclass_default_usda 
     3072             ! 
     3073             WRITE(numout,*) "Using a soilclass map with usda classification, to be read using XIOS" 
     3074             ! 
     3075             ALLOCATE(textrefrac(nbpt,nscm)) 
     3076             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
     3077             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3078 
     3079             CALL get_soilcorr_usda (nscm, textfrac_table) 
     3080 
     3081             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda' 
     3082 
     3083             CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1)) 
     3084             CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2)) 
     3085             CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3)) 
     3086             CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4)) 
     3087             CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5)) 
     3088             CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6)) 
     3089             CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7)) 
     3090             CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8)) 
     3091             CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9)) 
     3092             CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10)) 
     3093             CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11)) 
     3094             CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))         
     3095             CALL xios_orchidee_recv_field('soiltext13',textrefrac(:,13)) 
     3096 
     3097             CALL get_soilcorr_usda (nscm, textfrac_table) 
     3098             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda' 
     3099 
     3100             DO ib =1, nbpt 
     3101                clayfraction(ib) = 0.0 
     3102                DO ilf = 1,nscm 
     3103                   soilclass(ib,ilf)=textrefrac(ib,ilf) 
     3104                   clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf) 
     3105                   sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf) 
     3106                   siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf) 
     3107                ENDDO 
     3108 
     3109                sgn=SUM(soilclass(ib,:)) 
     3110 
     3111                IF (sgn < min_sechiba) THEN 
     3112                   soilclass(ib,:) = soilclass_default(:) 
     3113                   clayfraction(ib) = clayfraction_default 
     3114                   sandfraction(ib) = sandfraction_default 
     3115                   siltfraction(ib) = siltfraction_default 
     3116                   atext(ib)=0 
     3117                ELSE 
     3118                   soilclass(ib,:) = soilclass(ib,:) / sgn 
     3119                   clayfraction(ib) = clayfraction(ib) / sgn 
     3120                   sandfraction(ib) = sandfraction(ib) / sgn 
     3121                   siltfraction(ib) = siltfraction(ib) / sgn 
     3122                   atext(ib)=sgn 
     3123                ENDIF 
     3124             ENDDO 
     3125 
    29043126          CASE DEFAULT 
    29053127             WRITE(numout,*) 'slowproc_soilt:' 
    29063128             WRITE(numout,*) '  A non supported soil type classification has been chosen' 
    29073129             CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','') 
    2908           ENDSELECT 
    2909  
    2910           ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR) 
    2911           IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',& 
    2912             '','') 
    2913  
    2914           ! Assigning values to vmin, vmax 
    2915           vmin = un 
    2916           vmax = ntextinfile*un 
    2917            
    2918           ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR) 
    2919           IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','') 
    2920           variabletypevals = -un 
    2921            
    2922           !! Variables for interpweight 
    2923           ! Should negative values be set to zero from input file? 
    2924           nonegative = .FALSE. 
    2925           ! Type of mask to apply to the input data (see header for more details) 
    2926           maskingtype = 'mabove' 
    2927           ! Values to use for the masking 
    2928           maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /) 
    2929           ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used) 
    2930           namemaskvar = '' 
    2931            
    2932           CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours,        & 
    2933              contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,    &  
    2934              maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext) 
    2935  
    2936           ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR) 
    2937           IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','') 
    2938           ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR) 
    2939           IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','') 
    2940            
    2941           IF (printlev_loc >= 5) THEN 
    2942              WRITE(numout,*)'  slowproc_soilt after interpweight_2D' 
    2943              WRITE(numout,*)'  slowproc_soilt before starting loop nbpt:', nbpt 
    2944              WRITE(numout,*)"  slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..." 
     3130          END SELECT 
     3131 
     3132       ELSE              !    xios_interpolation 
     3133          ! Read and interpolate using stardard method with IOIPSL and aggregate 
     3134 
     3135          IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " & 
     3136               // TRIM(filename) // " for variable " // TRIM(variablename) 
     3137 
     3138          ! Name of the longitude and latitude in the input file 
     3139          lonname = 'nav_lon' 
     3140          latname = 'nav_lat' 
     3141 
     3142          IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " & 
     3143               // TRIM(filename) // " for variable " // TRIM(variablename) 
     3144 
     3145          IF ( TRIM(soil_classif) /= 'none' ) THEN 
     3146 
     3147             ! Define a variable for the number of soil textures in the input file 
     3148             SELECTCASE(soil_classif) 
     3149             CASE('zobler') 
     3150                ntextinfile=nzobler 
     3151             CASE('usda') 
     3152                ntextinfile=nscm 
     3153             CASE DEFAULT 
     3154                WRITE(numout,*) 'slowproc_soilt:' 
     3155                WRITE(numout,*) '  A non supported soil type classification has been chosen' 
     3156                CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','') 
     3157             ENDSELECT 
     3158 
     3159             ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR) 
     3160             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',& 
     3161                  '','') 
     3162 
     3163             ! Assigning values to vmin, vmax 
     3164             vmin = un 
     3165             vmax = ntextinfile*un 
     3166 
     3167             ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR) 
     3168             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','') 
     3169             variabletypevals = -un 
     3170 
     3171             !! Variables for interpweight 
     3172             ! Should negative values be set to zero from input file? 
     3173             nonegative = .FALSE. 
     3174             ! Type of mask to apply to the input data (see header for more details) 
     3175             maskingtype = 'mabove' 
     3176             ! Values to use for the masking 
     3177             maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /) 
     3178             ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used) 
     3179 
     3180             namemaskvar = '' 
     3181 
     3182             CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours,        & 
     3183                  contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,    & 
     3184                  maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext) 
     3185 
     3186             ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR) 
     3187             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','') 
     3188             ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR) 
     3189             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','') 
     3190 
     3191             IF (printlev_loc >= 5) THEN 
     3192                WRITE(numout,*)'  slowproc_soilt after interpweight_2D' 
     3193                WRITE(numout,*)'  slowproc_soilt before starting loop nbpt:', nbpt 
     3194                WRITE(numout,*)"  slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..." 
     3195             END IF 
     3196          ELSE 
     3197             IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt using default values all points are propertly ' // & 
     3198                  'interpolated atext = 1. everywhere!' 
     3199             atext = 1. 
    29453200          END IF 
    2946        ELSE 
    2947          IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt using default values all points are propertly ' // & 
    2948            'interpolated atext = 1. everywhere!' 
    2949          atext = 1. 
    2950        END IF 
    2951  
    2952     nbexp = 0 
    2953     SELECTCASE(soil_classif) 
    2954     CASE('none') 
    2955        ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
    2956        IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    2957        DO ib=1, nbpt 
    2958           soilclass(ib,:) = soilclass_default_fao 
    2959           clayfraction(ib) = clayfraction_default 
    2960           sandfraction(ib) = sandfraction_default 
    2961           siltfraction(ib) = siltfraction_default 
    2962        ENDDO 
    2963     CASE('zobler') 
    2964        ! 
    2965        soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes 
    2966        ! 
    2967        IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification" 
    2968        ! 
    2969        ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR) 
    2970        IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    2971        CALL get_soilcorr_zobler (nzobler, textfrac_table) 
    2972        ! 
    2973        ! 
    2974        IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt after getting table of textures' 
    2975        DO ib =1, nbpt 
    2976           soilclass(ib,:) = zero 
    2977           clayfraction(ib) = zero 
    2978           sandfraction(ib) = zero 
    2979           siltfraction(ib) = zero 
    2980           ! 
    2981           ! vecpos: List of positions where textures were not zero 
    2982           ! vecpos(1): number of not null textures found 
    2983           vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq') 
    2984           fopt = vecpos(1) 
    2985  
    2986           IF ( fopt .EQ. 0 ) THEN 
    2987              ! No points were found for current grid box, use default values 
    2988              nbexp = nbexp + 1 
    2989              soilclass(ib,:) = soilclass_default(:) 
    2990              clayfraction(ib) = clayfraction_default 
    2991              sandfraction(ib) = sandfraction_default 
    2992              siltfraction(ib) = siltfraction_default 
    2993           ELSE 
    2994              IF (fopt == nzobler) THEN 
    2995                 ! All textures are not zero 
    2996                 solt=(/(i,i=1,nzobler)/) 
    2997              ELSE 
    2998                DO ilf = 1,fopt 
    2999                  solt(ilf) = vecpos(ilf+1) 
    3000                END DO 
    3001              END IF 
    3002              ! 
    3003              !   Compute the fraction of each textural class 
    3004              ! 
    3005              sgn = 0. 
    3006              DO ilf = 1,fopt 
    3007                    ! 
    3008                    ! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE 
    3009                    ! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine) 
    3010                    ! For type 6 = glacier, default values are set and it is also taken into account during the normalization  
    3011                    ! of the fractions (done in interpweight_2D) 
    3012                    ! Note that type 0 corresponds to ocean but it is already removed using the mask above. 
    3013                    ! 
    3014                 IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &  
    3015                      (solt(ilf) .NE. 6) ) THEN 
    3016                    SELECT CASE(solt(ilf)) 
    3017                      CASE(1) 
    3018                         soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf)) 
    3019                      CASE(2) 
    3020                         soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
    3021                      CASE(3) 
    3022                         soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
    3023                      CASE(4) 
    3024                         soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
    3025                      CASE(5) 
    3026                         soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf)) 
    3027                      CASE(7) 
    3028                         soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
    3029                      CASE DEFAULT 
    3030                         WRITE(numout,*) 'We should not be here, an impossible case appeared' 
    3031                         CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','') 
    3032                    END SELECT 
    3033                    ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture) 
    3034                    ! over the zobler pixels composing the ORCHIDEE grid-cell 
    3035                    clayfraction(ib) = clayfraction(ib) + & 
    3036                         & textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf)) 
    3037                    sandfraction(ib) = sandfraction(ib) + & 
    3038                         & textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf)) 
    3039                    siltfraction(ib) = siltfraction(ib) + & 
    3040                         & textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf)) 
    3041                    ! Sum the fractions which are not glaciers nor ocean 
    3042                    sgn = sgn + textrefrac(ib,solt(ilf)) 
    3043                 ELSE 
    3044                    IF (solt(ilf) .GT. nzobler) THEN 
    3045                       WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program' 
    3046                       CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','') 
    3047                    ENDIF 
    3048                 END IF 
    3049              ENDDO 
    3050  
    3051              IF ( sgn .LT. min_sechiba) THEN 
    3052                 ! Set default values if grid cells were only covered by glaciers or ocean  
    3053                 ! or if now information on the source grid was found. 
    3054                 nbexp = nbexp + 1 
    3055                 soilclass(ib,:) = soilclass_default(:) 
     3201 
     3202          nbexp = 0 
     3203          SELECTCASE(soil_classif) 
     3204          CASE('none') 
     3205             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
     3206             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3207             DO ib=1, nbpt 
     3208                soilclass(ib,:) = soilclass_default_fao 
    30563209                clayfraction(ib) = clayfraction_default 
    30573210                sandfraction(ib) = sandfraction_default 
    30583211                siltfraction(ib) = siltfraction_default 
    3059              ELSE 
    3060                 ! Normalize using the fraction of surface not including glaciers and ocean 
    3061                 soilclass(ib,:) = soilclass(ib,:)/sgn 
    3062                 clayfraction(ib) = clayfraction(ib)/sgn 
    3063                 sandfraction(ib) = sandfraction(ib)/sgn 
    3064                 siltfraction(ib) = siltfraction(ib)/sgn 
    3065              ENDIF 
    3066           ENDIF 
    3067        ENDDO 
    3068        
    3069        ! The "USDA" case reads a map of the 12 USDA texture classes,  
    3070        ! such as to assign the corresponding soil properties 
    3071        CASE("usda") 
    3072           IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification" 
    3073  
    3074           soilclass_default=soilclass_default_usda 
    3075  
    3076           ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
    3077           IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
    3078  
    3079           CALL get_soilcorr_usda (nscm, textfrac_table) 
    3080  
    3081           IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda' 
    3082           ! 
    3083           DO ib =1, nbpt 
    3084           ! GO through the point we have found 
    3085           ! 
    3086           ! 
    3087           ! Provide which textures were found 
    3088           ! vecpos: List of positions where textures were not zero 
    3089           !   vecpos(1): number of not null textures found 
    3090           vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq') 
    3091           fopt = vecpos(1) 
    3092           
    3093           ! 
    3094           !    Check that we found some points 
    3095           ! 
    3096           soilclass(ib,:) = 0.0 
    3097           clayfraction(ib) = 0.0 
    3098            
    3099           IF ( fopt .EQ. 0) THEN 
    3100              ! No points were found for current grid box, use default values 
    3101              IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib 
    3102              nbexp = nbexp + 1 
    3103              soilclass(ib,:) = soilclass_default 
    3104              clayfraction(ib) = clayfraction_default 
    3105              sandfraction(ib) = sandfraction_default 
    3106              siltfraction(ib) = siltfraction_default 
    3107           ELSE 
    3108              IF (fopt == nscm) THEN 
    3109                 ! All textures are not zero 
    3110                 solt(:) = (/(i,i=1,nscm)/) 
    3111              ELSE 
    3112                DO ilf = 1,fopt 
    3113                  solt(ilf) = vecpos(ilf+1)  
    3114                END DO 
    3115              END IF 
    3116              
     3212             ENDDO 
     3213          CASE('zobler') 
     3214             ! 
     3215             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes 
     3216             ! 
     3217             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification" 
     3218             ! 
     3219             ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR) 
     3220             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3221             CALL get_soilcorr_zobler (nzobler, textfrac_table) 
     3222                        
     3223             IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt after getting table of textures' 
     3224             DO ib =1, nbpt 
     3225                soilclass(ib,:) = zero 
     3226                clayfraction(ib) = zero 
     3227                sandfraction(ib) = zero 
     3228                siltfraction(ib) = zero 
    31173229                ! 
    3118                 ! 
    3119                 !   Compute the fraction of each textural class   
    3120                 ! 
    3121                 ! 
    3122                 DO ilf = 1,fopt 
    3123                    IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN 
    3124                       soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf)) 
    3125                       clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) *                & 
    3126                            textrefrac(ib,solt(ilf)) 
    3127                       sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * & 
    3128                            textrefrac(ib,solt(ilf)) 
    3129                       siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * & 
    3130                         textrefrac(ib,solt(ilf)) 
    3131                    ELSE 
    3132                       IF (solt(ilf) .GT. nscm) THEN 
    3133                          WRITE(*,*) 'The file contains a soil color class which is incompatible with this program' 
    3134                          CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','') 
    3135                       ENDIF 
    3136                    ENDIF 
    3137                    ! 
    3138                 ENDDO 
    3139  
    3140                 ! Set default values if the surface in source file is too small 
    3141                 IF ( atext(ib) .LT. min_sechiba) THEN 
     3230                ! vecpos: List of positions where textures were not zero 
     3231                ! vecpos(1): number of not null textures found 
     3232                vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq') 
     3233                fopt = vecpos(1) 
     3234 
     3235                IF ( fopt .EQ. 0 ) THEN 
     3236                   ! No points were found for current grid box, use default values 
    31423237                   nbexp = nbexp + 1 
    31433238                   soilclass(ib,:) = soilclass_default(:) 
     
    31453240                   sandfraction(ib) = sandfraction_default 
    31463241                   siltfraction(ib) = siltfraction_default 
     3242 
     3243                ELSE 
     3244                   IF (fopt == nzobler) THEN 
     3245                      ! All textures are not zero 
     3246                      solt=(/(i,i=1,nzobler)/) 
     3247                   ELSE 
     3248                      DO ilf = 1,fopt 
     3249                         solt(ilf) = vecpos(ilf+1) 
     3250                      END DO 
     3251                   END IF 
     3252                   ! 
     3253                   !   Compute the fraction of each textural class 
     3254                   ! 
     3255                   sgn = 0. 
     3256                   DO ilf = 1,fopt 
     3257                      ! 
     3258                      ! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE 
     3259                      ! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine) 
     3260                      ! For type 6 = glacier, default values are set and it is also taken into account during the normalization 
     3261                      ! of the fractions (done in interpweight_2D) 
     3262                      ! Note that type 0 corresponds to ocean but it is already removed using the mask above. 
     3263                     ! 
     3264                      IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. & 
     3265                           (solt(ilf) .NE. 6) ) THEN 
     3266                         SELECT CASE(solt(ilf)) 
     3267                         CASE(1) 
     3268                            soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf)) 
     3269                         CASE(2) 
     3270                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
     3271                         CASE(3) 
     3272                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
     3273                         CASE(4) 
     3274                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
     3275                         CASE(5) 
     3276                            soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf)) 
     3277                         CASE(7) 
     3278                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf)) 
     3279                         CASE DEFAULT 
     3280                            WRITE(numout,*) 'We should not be here, an impossible case appeared' 
     3281                            CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','') 
     3282                         END SELECT 
     3283                         ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture) 
     3284                         ! over the zobler pixels composing the ORCHIDEE grid-cell 
     3285                         clayfraction(ib) = clayfraction(ib) + & 
     3286                              & textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf)) 
     3287                         sandfraction(ib) = sandfraction(ib) + & 
     3288                              & textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf)) 
     3289                         siltfraction(ib) = siltfraction(ib) + & 
     3290                              & textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf)) 
     3291                         ! Sum the fractions which are not glaciers nor ocean 
     3292                         sgn = sgn + textrefrac(ib,solt(ilf)) 
     3293                      ELSE 
     3294                         IF (solt(ilf) .GT. nzobler) THEN 
     3295                            WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program' 
     3296                            CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','') 
     3297                         ENDIF 
     3298                      END IF 
     3299                   ENDDO 
     3300 
     3301                   IF ( sgn .LT. min_sechiba) THEN 
     3302                      ! Set default values if grid cells were only covered by glaciers or ocean 
     3303                      ! or if now information on the source grid was found. 
     3304                      nbexp = nbexp + 1 
     3305                      soilclass(ib,:) = soilclass_default(:) 
     3306                      clayfraction(ib) = clayfraction_default 
     3307                      sandfraction(ib) = sandfraction_default 
     3308                      siltfraction(ib) = siltfraction_default 
     3309                   ELSE 
     3310                      ! Normalize using the fraction of surface not including glaciers and ocean 
     3311                      soilclass(ib,:) = soilclass(ib,:)/sgn 
     3312                      clayfraction(ib) = clayfraction(ib)/sgn 
     3313                      sandfraction(ib) = sandfraction(ib)/sgn 
     3314                      siltfraction(ib) = siltfraction(ib)/sgn 
     3315                   ENDIF 
    31473316                ENDIF 
    3148              ENDIF 
    3149  
    3150           ENDDO 
    3151            
    3152           IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda' 
    3153            
    3154        CASE DEFAULT 
    3155           WRITE(numout,*) 'slowproc_soilt _______' 
    3156           WRITE(numout,*) '  A non supported soil type classification has been chosen' 
    3157           CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','') 
    3158        ENDSELECT 
    3159        IF (printlev_loc >= 5 ) WRITE(numout,*)'  slowproc_soilt end of type classification' 
    3160  
    3161        IF ( nbexp .GT. 0 ) THEN 
    3162           WRITE(numout,*) 'slowproc_soilt:' 
    3163           WRITE(numout,*) '  The interpolation of the bare soil albedo had ', nbexp 
    3164           WRITE(numout,*) '  points without data. This are either coastal points or ice covered land.' 
    3165           WRITE(numout,*) '  The problem was solved by using the default soil types.' 
    3166        ENDIF 
    3167  
    3168        IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals) 
    3169        IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac) 
    3170        IF (ALLOCATED(solt)) DEALLOCATE (solt) 
    3171        IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table) 
    3172      
    3173     ENDIF        !      xios_interpolation  
    3174      
     3317             ENDDO 
     3318 
     3319             ! The "USDA" case reads a map of the 12 USDA texture classes, 
     3320             ! such as to assign the corresponding soil properties 
     3321          CASE("usda") 
     3322             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification" 
     3323 
     3324             soilclass_default=soilclass_default_usda 
     3325 
     3326 
     3327             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR) 
     3328             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','') 
     3329 
     3330             CALL get_soilcorr_usda (nscm, textfrac_table) 
     3331 
     3332             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda' 
     3333             ! 
     3334             DO ib =1, nbpt 
     3335                ! GO through the point we have found 
     3336                ! 
     3337                ! Provide which textures were found 
     3338                ! vecpos: List of positions where textures were not zero 
     3339                !   vecpos(1): number of not null textures found 
     3340                vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq') 
     3341                fopt = vecpos(1) 
     3342                ! 
     3343                !    Check that we found some points 
     3344                ! 
     3345                soilclass(ib,:) = 0.0 
     3346                clayfraction(ib) = 0.0 
     3347                sandfraction(ib) = 0.0 
     3348                siltfraction(ib) = 0.0 
     3349 
     3350                IF ( fopt .EQ. 0) THEN 
     3351                   ! No points were found for current grid box, use default values 
     3352                   IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib 
     3353                   nbexp = nbexp + 1 
     3354                   soilclass(ib,:) = soilclass_default 
     3355                   clayfraction(ib) = clayfraction_default 
     3356                   sandfraction(ib) = sandfraction_default 
     3357                   siltfraction(ib) = siltfraction_default 
     3358                ELSE 
     3359                   IF (fopt == nscm) THEN 
     3360                      ! All textures are not zero 
     3361                      solt(:) = (/(i,i=1,nscm)/) 
     3362                   ELSE 
     3363                      DO ilf = 1,fopt 
     3364                         solt(ilf) = vecpos(ilf+1) 
     3365                      END DO 
     3366                   END IF 
     3367 
     3368                   !   Compute the fraction of each textural class 
     3369                   DO ilf = 1,fopt 
     3370                      IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN 
     3371                         soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf)) 
     3372                         clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) *                & 
     3373                              textrefrac(ib,solt(ilf)) 
     3374                         sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * & 
     3375                              textrefrac(ib,solt(ilf)) 
     3376                         siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * & 
     3377                              textrefrac(ib,solt(ilf)) 
     3378                      ELSE 
     3379                         IF (solt(ilf) .GT. nscm) THEN 
     3380                            WRITE(*,*) 'The file contains a soil color class which is incompatible with this program' 
     3381                            CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','') 
     3382                         ENDIF 
     3383                      ENDIF 
     3384                      ! 
     3385                   ENDDO 
     3386 
     3387                   ! Set default values if the surface in source file is too small 
     3388                   IF ( atext(ib) .LT. min_sechiba) THEN 
     3389                      nbexp = nbexp + 1 
     3390                      soilclass(ib,:) = soilclass_default(:) 
     3391                      clayfraction(ib) = clayfraction_default 
     3392                      sandfraction(ib) = sandfraction_default 
     3393                      siltfraction(ib) = siltfraction_default 
     3394                   ENDIF 
     3395                ENDIF 
     3396 
     3397             ENDDO 
     3398 
     3399             IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda' 
     3400 
     3401          CASE DEFAULT 
     3402             WRITE(numout,*) 'slowproc_soilt _______' 
     3403             WRITE(numout,*) '  A non supported soil type classification has been chosen' 
     3404             CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','') 
     3405          ENDSELECT 
     3406          IF (printlev_loc >= 5 ) WRITE(numout,*)'  slowproc_soilt end of type classification' 
     3407 
     3408          IF ( nbexp .GT. 0 ) THEN 
     3409             WRITE(numout,*) 'slowproc_soilt:' 
     3410             WRITE(numout,*) '  The interpolation of variable soiltext had ', nbexp 
     3411             WRITE(numout,*) '  points without data. This are either coastal points or ice covered land.' 
     3412             WRITE(numout,*) '  The problem was solved by using the default soil types.' 
     3413          ENDIF 
     3414 
     3415          IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals) 
     3416          IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac) 
     3417          IF (ALLOCATED(solt)) DEALLOCATE (solt) 
     3418          IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table) 
     3419 
     3420       ENDIF        !      xios_interpolation 
     3421 
     3422       !salma: calculate njsc 
     3423       njsc(:) = 0 
     3424       DO ib = 1, nbpt 
     3425          njsc(ib) = MAXLOC(soilclass(ib,:),1) 
     3426       ENDDO 
     3427 
     3428       IF (spmipexp == 'exp1') THEN 
     3429              IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt: Read soil hydraulic parameters with IOIPSL' 
     3430 
     3431              ! Read using IOIPSL and interpolate using aggregate tool in ORCHIDEE 
     3432 
     3433              !Config Key   = PARAM_FILE 
     3434              !Config Desc  = Name of file from which soil parameter  values are read 
     3435              !Config Def   = params_sp_mip.nc 
     3436              !Config Help  = The name of the file to be opened to read values of parameters. 
     3437              !Config         The data from this file is then interpolated to the grid of 
     3438              !Config         of the model. 
     3439              !Config Units = [FILE] 
     3440              ! 
     3441              ! params_sp_mip.nc file is 0.5 deg soil hydraulic parameters file provided by sp_mip 
     3442 
     3443              filename = 'params_sp_mip.nc' 
     3444              CALL getin_p('PARAM_FILE',filename) 
     3445 
     3446              !! Variables for interpweight 
     3447              ! Type of calculation of cell fractions 
     3448              fractype = 'default' 
     3449              ! Name of the longitude and latitude in the input file 
     3450              lonname = 'nav_lon' 
     3451              latname = 'nav_lat' 
     3452              ! Assigning values to vmin, vmax (there are not types/categories 
     3453              vmin =0. 
     3454              vmax = 99999. 
     3455              !! Variables for interpweight 
     3456              ! Should negative values be set to zero from input file? 
     3457              nonegative = .FALSE. 
     3458              ! Type of mask to apply to the input data (see header for more details) 
     3459              maskingtype = 'mabove' 
     3460              ! Values to use for the masking 
     3461              maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /) 
     3462              ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used) 
     3463              namemaskvar = '' 
     3464 
     3465              variablename = 'ks' 
     3466              IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " & 
     3467                   // TRIM(filename) // " for variable " // TRIM(variablename) 
     3468              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3469                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3470                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3471                   ks, aparam) 
     3472              WRITE(numout,*) 'ks map is read _______' 
     3473 
     3474              variablename = 'alpha' 
     3475              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3476                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3477                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3478                   avan, aparam) 
     3479              WRITE(numout,*) 'avan map read _______' 
     3480 
     3481              variablename = 'thetar' 
     3482              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3483                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3484                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3485                   mcr, aparam) 
     3486              WRITE(numout,*) 'thetar map read _______' 
     3487 
     3488              variablename = 'thetas' 
     3489              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3490                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3491                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3492                   mcs, aparam) 
     3493              WRITE(numout,*) 'thetas map read _______' 
     3494 
     3495              variablename = 'thetapwpvg' 
     3496              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3497                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3498                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3499                   mcw, aparam) 
     3500              WRITE(numout,*) 'thetapwpvg map read _______' 
     3501 
     3502              variablename = 'thetafcvg' 
     3503              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3504                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3505                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3506                   mcfc, aparam) 
     3507              WRITE(numout,*) 'thetafcvg map read _______' 
     3508 
     3509              variablename = 'nvg' 
     3510              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                & 
     3511                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     & 
     3512                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              & 
     3513                   nvan, aparam) 
     3514              WRITE(numout,*) 'nvan map read _______' 
     3515 
     3516      ELSE 
     3517          IF (spmipexp == 'exp2') THEN 
     3518            ! Texture map from SP-MIP, thus Soilgrids modified 
     3519            nvan(:) = nvan_usda(njsc(:)) 
     3520            avan(:) = avan_usda(njsc(:)) 
     3521            mcr(:) = mcr_usda(njsc(:)) 
     3522            mcs(:) = mcs_usda(njsc(:)) 
     3523            ks(:) = ks_usda(njsc(:)) 
     3524            mcfc(:) = mcf_usda(njsc(:)) 
     3525            mcw(:) = mcw_usda(njsc(:)) 
     3526            ! on aura pcent(:) = pcent(njsc(:)) dans hydrol 
     3527         ELSE 
     3528            ! salma: here we are in exp3 -- Zobler map 
     3529            nvan(:) = nvan_fao(njsc(:)) 
     3530            avan(:) = avan_fao(njsc(:)) 
     3531            mcr(:) = mcr_fao(njsc(:)) 
     3532            mcs(:) = mcs_fao(njsc(:)) 
     3533            ks(:) = ks_fao(njsc(:)) 
     3534            mcfc(:) = mcf_fao(njsc(:)) 
     3535            mcw(:) = mcw_fao(njsc(:)) 
     3536 
     3537          ENDIF !if EXP2 
     3538       ENDIF !if EXP1 
     3539     ENDIF !if EXP4 
     3540              
    31753541    ! Write diagnostics 
    31763542    CALL xios_orchidee_send_field("atext",atext) 
    3177  
    31783543    CALL xios_orchidee_send_field("interp_diag_atext",atext) 
    31793544    CALL xios_orchidee_send_field("interp_diag_soilclass",soilclass) 
     
    35853950    !! 0.4 Local variables 
    35863951 
    3587     INTEGER(i_std),PARAMETER :: nbtypes_usda = 12                    !! Number of USDA texture classes (unitless) 
     3952    INTEGER(i_std),PARAMETER :: nbtypes_usda = 13                    !! Number of USDA texture classes (unitless) 
    35883953    INTEGER(i_std) :: n                                              !! Index (unitless) 
    35893954     
     
    36193984    textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay 
    36203985    textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay 
     3986    textfrac_table(13,2:3) = (/ 0.15, 0.55 /) ! Clay 
    36213987 
    36223988    ! Fraction of silt 
Note: See TracChangeset for help on using the changeset viewer.