source: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90 @ 7337

Last change on this file since 7337 was 7337, checked in by agnes.ducharne, 8 months ago

Simplification of soil texture processing, cf ticket #416: when using the Zobler map, the soil parameters are no more taken from 3-value "FAO" vectors in constantes_soil_var.f90, but from 13-value USDA vectors, owing to a pointer fao2usda. A 5-day running test with the Zobler map shows no changes, but over Greenland (class 6=ice in Zobler map).

  • Property svn:keywords set to Revision Date HeadURL Date Author Revision
File size: 348.6 KB
Line 
1! ===================================================================================================\n
2! MODULE        : hydrol
3!
4! CONTACT       : orchidee-help _at_ listes.ipsl.fr
5!
6! LICENCE       : IPSL (2006)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF        This module computes the soil moisture processes on continental points.
10!!
11!!\n DESCRIPTION : contains hydrol_main, hydrol_initialize, hydrol_finalise, hydrol_init,
12!!                 hydrol_var_init, hydrol_waterbal, hydrol_alma,
13!!                 hydrol_vegupd, hydrol_canop, hydrol_flood, hydrol_soil.
14!!                 The assumption in this module is that very high vertical resolution is
15!!                 needed in order to properly resolve the vertical diffusion of water in
16!!                 the soils. Furthermore we have taken into account the sub-grid variability
17!!                 of soil properties and vegetation cover by allowing the co-existence of
18!!                 different soil moisture columns in the same grid box.
19!!                 This routine was originaly developed by Patricia deRosnay.
20!!
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!!                 
27!! REFERENCE(S) :
28!! - de Rosnay, P., J. Polcher, M. Bruen, and K. Laval, Impact of a physically based soil
29!! water flow and soil-plant interaction representation for modeling large-scale land surface
30!! processes, J. Geophys. Res, 107 (10.1029), 2002. \n
31!! - de Rosnay, P. and Polcher J. (1998) Modeling root water uptake in a complex land surface scheme coupled
32!! to a GCM. Hydrology and Earth System Sciences, 2(2-3):239-256. \n
33!! - de Rosnay, P., M. Bruen, and J. Polcher, Sensitivity of surface fluxes to the number of layers in the soil
34!! model used in GCMs, Geophysical research letters, 27 (20), 3329 - 3332, 2000. \n
35!! - d’Orgeval, T., J. Polcher, and P. De Rosnay, Sensitivity of the West African hydrological
36!! cycle in ORCHIDEE to infiltration processes, Hydrol. Earth Syst. Sci. Discuss, 5, 2251 - 2292, 2008. \n
37!! - Carsel, R., and R. Parrish, Developing joint probability distributions of soil water retention
38!! characteristics, Water Resources Research, 24 (5), 755 - 769, 1988. \n
39!! - Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated porous
40!! media, Water Resources Research, 12 (3), 513 - 522, 1976. \n
41!! - Van Genuchten, M., A closed-form equation for predicting the hydraulic conductivity of
42!! unsaturated soils, Soil Science Society of America Journal, 44 (5), 892 - 898, 1980. \n
43!! - Campoy, A., Ducharne, A., Cheruy, F., Hourdin, F., Polcher, J., and Dupont, J.-C., Response
44!! of land surface fluxes and precipitation to different soil bottom hydrological conditions in a
45!! general circulation model,  J. Geophys. Res, in press, 2013. \n
46!! - Gouttevin, I., Krinner, G., Ciais, P., Polcher, J., and Legout, C. , 2012. Multi-scale validation
47!! of a new soil freezing scheme for a land-surface model with physically-based hydrology.
48!! 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
51!!
52!! SVN          :
53!! $HeadURL$
54!! $Date$
55!! $Revision$
56!! \n
57!_ ===============================================================================================\n
58MODULE hydrol
59
60  USE ioipsl
61  USE xios_orchidee
62  USE constantes
63  USE time, ONLY : one_day, dt_sechiba, julian_diff
64  USE constantes_soil
65  USE pft_parameters
66  USE sechiba_io_p
67  USE grid
68  USE explicitsnow
69
70  IMPLICIT NONE
71
72  PRIVATE
73  PUBLIC :: hydrol_main, hydrol_initialize, hydrol_finalize, hydrol_clear
74
75  !
76  ! variables used inside hydrol module : declaration and initialisation
77  !
78  LOGICAL, SAVE                                   :: doponds=.FALSE.           !! Reinfiltration flag (true/false)
79!$OMP THREADPRIVATE(doponds)
80  REAL(r_std), SAVE                               :: froz_frac_corr            !! Coefficient for water frozen fraction correction
81!$OMP THREADPRIVATE(froz_frac_corr)
82  REAL(r_std), SAVE                               :: max_froz_hydro            !! Coefficient for water frozen fraction correction
83!$OMP THREADPRIVATE(max_froz_hydro)
84  REAL(r_std), SAVE                               :: smtot_corr                !! Coefficient for water frozen fraction correction
85!$OMP THREADPRIVATE(smtot_corr)
86  LOGICAL, SAVE                                   :: do_rsoil=.FALSE.          !! Flag to calculate rsoil for bare soile evap
87                                                                               !! (true/false)
88!$OMP THREADPRIVATE(do_rsoil)
89  LOGICAL, SAVE                                   :: ok_dynroot                !! Flag to activate dynamic root profile to optimize soil 
90                                                                               !! moisture usage, similar to Beer et al.2007
91!$OMP THREADPRIVATE(ok_dynroot)
92  CHARACTER(LEN=80) , SAVE                        :: var_name                  !! To store variables names for I/O
93!$OMP THREADPRIVATE(var_name)
94  !
95  REAL(r_std), PARAMETER                          :: allowed_err =  2.0E-8_r_std
96  REAL(r_std), PARAMETER                          :: EPS1 = EPSILON(un)      !! A small number
97 
98  ! one dimension array allocated, computed, saved and got in hydrol module
99  ! Values per soil type
100  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pcent               !! Fraction of saturated volumetric soil moisture above
101                                                                         !! which transpir is max (0-1, unitless)
102!$OMP THREADPRIVATE(pcent)                                                               
103  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mc_awet             !! Vol. wat. cont. above which albedo is cst
104                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
105!$OMP THREADPRIVATE(mc_awet)                                             
106  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mc_adry             !! Vol. wat. cont. below which albedo is cst
107                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
108!$OMP THREADPRIVATE(mc_adry)                                             
109  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tot_watveg_beg   !! Total amount of water on vegetation at start of time
110                                                                         !! step @tex $(kg m^{-2})$ @endtex
111!$OMP THREADPRIVATE(tot_watveg_beg)                                     
112  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tot_watveg_end   !! Total amount of water on vegetation at end of time step
113                                                                         !!  @tex $(kg m^{-2})$ @endtex
114!$OMP THREADPRIVATE(tot_watveg_end)                                     
115  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tot_watsoil_beg  !! Total amount of water in the soil at start of time step
116                                                                         !!  @tex $(kg m^{-2})$ @endtex
117!$OMP THREADPRIVATE(tot_watsoil_beg)                                     
118  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tot_watsoil_end  !! Total amount of water in the soil at end of time step
119                                                                         !!  @tex $(kg m^{-2})$ @endtex
120!$OMP THREADPRIVATE(tot_watsoil_end)                                     
121  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: snow_beg         !! Total amount of snow at start of time step
122                                                                         !!  @tex $(kg m^{-2})$ @endtex
123!$OMP THREADPRIVATE(snow_beg)                                           
124  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: snow_end         !! Total amount of snow at end of time step
125                                                                         !!  @tex $(kg m^{-2})$ @endtex
126!$OMP THREADPRIVATE(snow_end)                                           
127  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: delsoilmoist     !! Change in soil moisture @tex $(kg m^{-2})$ @endtex
128!$OMP THREADPRIVATE(delsoilmoist)                                         
129  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: delintercept     !! Change in interception storage
130                                                                         !!  @tex $(kg m^{-2})$ @endtex
131!$OMP THREADPRIVATE(delintercept)                                       
132  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: delswe           !! Change in SWE @tex $(kg m^{-2})$ @endtex
133!$OMP THREADPRIVATE(delswe)
134  REAL(r_std),ALLOCATABLE, SAVE, DIMENSION (:)       :: undermcr         !! Nb of tiles under mcr for a given time step
135!$OMP THREADPRIVATE(undermcr)
136  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_veget       !! zero/one when veget fraction is zero/higher (1)
137!$OMP THREADPRIVATE(mask_veget)
138  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_soiltile    !! zero/one where soil tile is zero/higher (1)
139!$OMP THREADPRIVATE(mask_soiltile)
140  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: humrelv          !! Water stress index for transpiration
141                                                                         !! for each soiltile x PFT couple (0-1, unitless)
142!$OMP THREADPRIVATE(humrelv)
143  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: vegstressv       !! Water stress index for vegetation growth
144                                                                         !! for each soiltile x PFT couple (0-1, unitless)
145!$OMP THREADPRIVATE(vegstressv)
146  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:):: us               !! Water stress index for transpiration
147                                                                         !! (by soil layer and PFT) (0-1, unitless)
148!$OMP THREADPRIVATE(us)
149  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: precisol         !! Throughfall+Totmelt per PFT
150                                                                         !!  @tex $(kg m^{-2})$ @endtex
151!$OMP THREADPRIVATE(precisol)
152  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: throughfall      !! Throughfall per PFT
153                                                                         !!  @tex $(kg m^{-2})$ @endtex
154!$OMP THREADPRIVATE(throughfall)
155  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: precisol_ns      !! Throughfall per soiltile
156                                                                         !!  @tex $(kg m^{-2})$ @endtex
157!$OMP THREADPRIVATE(precisol_ns)
158  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: ae_ns            !! Bare soil evaporation per soiltile
159                                                                         !!  @tex $(kg m^{-2})$ @endtex
160!$OMP THREADPRIVATE(ae_ns)
161  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: free_drain_coef  !! Coefficient for free drainage at bottom
162                                                                         !!  (0-1, unitless)
163!$OMP THREADPRIVATE(free_drain_coef)
164  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: zwt_force        !! Prescribed water table depth (m)
165!$OMP THREADPRIVATE(zwt_force)
166  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: frac_bare_ns     !! Evaporating bare soil fraction per soiltile
167                                                                         !!  (0-1, unitless)
168!$OMP THREADPRIVATE(frac_bare_ns)
169  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: rootsink         !! Transpiration sink by soil layer and soiltile
170                                                                         !! @tex $(kg m^{-2})$ @endtex
171!$OMP THREADPRIVATE(rootsink)
172  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: subsnowveg       !! Sublimation of snow on vegetation
173                                                                         !!  @tex $(kg m^{-2})$ @endtex
174!$OMP THREADPRIVATE(subsnowveg)
175  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: subsnownobio     !! Sublimation of snow on other surface types 
176                                                                         !! (ice, lakes,...) @tex $(kg m^{-2})$ @endtex
177!$OMP THREADPRIVATE(subsnownobio)
178  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: icemelt          !! Ice melt @tex $(kg m^{-2})$ @endtex
179!$OMP THREADPRIVATE(icemelt)
180  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: subsinksoil      !! Excess of sublimation as a sink for the soil
181                                                                         !! @tex $(kg m^{-2})$ @endtex
182!$OMP THREADPRIVATE(subsinksoil)
183  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: vegtot           !! Total Total fraction of grid-cell covered by PFTs
184                                                                         !! (bare soil + vegetation) (1; 1)
185!$OMP THREADPRIVATE(vegtot)
186  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: resdist          !! Soiltile values from previous time-step (1; 1)
187!$OMP THREADPRIVATE(resdist)
188  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: vegtot_old       !! Total Total fraction of grid-cell covered by PFTs
189                                                                         !! from previous time-step (1; 1)
190!$OMP THREADPRIVATE(vegtot_old)
191  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: mx_eau_var       !! Maximum water content of the soil @tex $(kg m^{-2})$ @endtex
192!$OMP THREADPRIVATE(mx_eau_var)
193
194  ! arrays used by cwrr scheme
195  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: nroot            !! Normalized root length fraction in each soil layer
196                                                                         !! (0-1, unitless)
197                                                                         !! DIM = kjpindex * nvm * nslm
198!$OMP THREADPRIVATE(nroot)
199  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: kfact_root       !! Factor to increase Ks towards the surface
200                                                                         !! (unitless)
201                                                                         !! DIM = kjpindex * nslm * nstm
202!$OMP THREADPRIVATE(kfact_root)
203  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: kfact            !! Factor to reduce Ks with depth (unitless)
204                                                                         !! DIM = nslm * kjpindex
205!$OMP THREADPRIVATE(kfact)
206  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: zz               !! Depth of nodes [znh in vertical_soil] transformed into (mm)
207!$OMP THREADPRIVATE(zz)
208  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: dz               !! Internode thickness [dnh in vertical_soil] transformed into (mm)
209!$OMP THREADPRIVATE(dz)
210  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: dh               !! Layer thickness [dlh in vertical_soil] transformed into (mm)
211!$OMP THREADPRIVATE(dh)
212  INTEGER(i_std), SAVE                               :: itopmax          !! Number of layers where the node is above 0.1m depth
213!$OMP THREADPRIVATE(itopmax)
214  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: mc_lin   !! 50 Vol. Wat. Contents to linearize K and D, for each texture
215                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex
216                                                                 !! DIM = imin:imax * kjpindex
217!$OMP THREADPRIVATE(mc_lin)
218  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: k_lin    !! 50 values of unsaturated K, for each soil layer and texture
219                                                                 !!  @tex $(mm d^{-1})$ @endtex
220                                                                 !! DIM = imin:imax * nslm * kjpindex
221!$OMP THREADPRIVATE(k_lin)
222  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: d_lin    !! 50 values of diffusivity D, for each soil layer and texture
223                                                                 !!  @tex $(mm^2 d^{-1})$ @endtex
224                                                                 !! DIM = imin:imax * nslm * kjpindex
225!$OMP THREADPRIVATE(d_lin)
226  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: a_lin    !! 50 values of the slope in K=a*mc+b, for each soil layer and texture
227                                                                 !!  @tex $(mm d^{-1})$ @endtex
228                                                                 !! DIM = imin:imax * nslm * kjpindex
229!$OMP THREADPRIVATE(a_lin)
230  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: b_lin    !! 50 values of y-intercept in K=a*mc+b, for each soil layer and texture
231                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex
232                                                                 !! DIM = imin:imax * nslm * kjpindex
233!$OMP THREADPRIVATE(b_lin)
234
235  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: humtot   !! Total Soil Moisture @tex $(kg m^{-2})$ @endtex
236!$OMP THREADPRIVATE(humtot)
237  LOGICAL, ALLOCATABLE, SAVE, DIMENSION (:)          :: resolv   !! Mask of land points where to solve the diffusion equation
238                                                                 !! (true/false)
239!$OMP THREADPRIVATE(resolv)
240
241!! for output
242  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: kk_moy   !! Mean hydraulic conductivity over soiltiles (mm/d)
243!$OMP THREADPRIVATE(kk_moy)
244  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: kk       !! Hydraulic conductivity for each soiltiles (mm/d)
245!$OMP THREADPRIVATE(kk)
246  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: avan_mod_tab  !! VG parameter a modified from  exponantial profile
247                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,kjpindex)
248!$OMP THREADPRIVATE(avan_mod_tab) 
249  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: nvan_mod_tab  !! VG parameter n  modified from  exponantial profile
250                                                                      !! (unitless) !! DIMENSION (nslm,kjpindex) 
251!$OMP THREADPRIVATE(nvan_mod_tab)
252 
253!! linarization coefficients of hydraulic conductivity K (hydrol_soil_coef)
254  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: k        !! Hydraulic conductivity K for each soil layer
255                                                                 !!  @tex $(mm d^{-1})$ @endtex
256                                                                 !! DIM = (:,nslm)
257!$OMP THREADPRIVATE(k)
258  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: a        !! Slope in K=a*mc+b(:,nslm)
259                                                                 !!  @tex $(mm d^{-1})$ @endtex
260                                                                 !! DIM = (:,nslm)
261!$OMP THREADPRIVATE(a)
262  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: b        !! y-intercept in K=a*mc+b
263                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex
264                                                                 !! DIM = (:,nslm)
265!$OMP THREADPRIVATE(b)
266!! linarization coefficients of hydraulic diffusivity D (hydrol_soil_coef)
267  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: d        !! Diffusivity D for each soil layer
268                                                                 !!  @tex $(mm^2 d^{-1})$ @endtex
269                                                                 !! DIM = (:,nslm)
270!$OMP THREADPRIVATE(d)
271!! matrix coefficients (hydrol_soil_tridiag and hydrol_soil_setup), see De Rosnay (1999), p155-157
272  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: e        !! Left-hand tridiagonal matrix coefficients
273!$OMP THREADPRIVATE(e)
274  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: f        !! Left-hand tridiagonal matrix coefficients
275!$OMP THREADPRIVATE(f)
276  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: g1       !! Left-hand tridiagonal matrix coefficients
277!$OMP THREADPRIVATE(g1)
278
279  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: ep       !! Right-hand matrix coefficients
280!$OMP THREADPRIVATE(ep)
281  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: fp       !! Right-hand atrix coefficients
282!$OMP THREADPRIVATE(fp)
283  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: gp       !! Right-hand atrix coefficients
284!$OMP THREADPRIVATE(gp)
285  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: rhs      !! Right-hand system
286!$OMP THREADPRIVATE(rhs)
287  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: srhs     !! Temporarily stored rhs
288!$OMP THREADPRIVATE(srhs)
289  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: tmat             !! Left-hand tridiagonal matrix
290!$OMP THREADPRIVATE(tmat)
291  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: stmat            !! Temporarily stored tmat
292  !$OMP THREADPRIVATE(stmat)
293  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: water2infilt     !! Water to be infiltrated
294                                                                         !! @tex $(kg m^{-2})$ @endtex
295!$OMP THREADPRIVATE(water2infilt)
296  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc              !! Total moisture content per soiltile
297                                                                         !!  @tex $(kg m^{-2})$ @endtex
298!$OMP THREADPRIVATE(tmc)
299  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmcr             !! Total moisture content at residual per soiltile
300                                                                         !!  @tex $(kg m^{-2})$ @endtex
301!$OMP THREADPRIVATE(tmcr)
302  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmcs             !! Total moisture content at saturation per soiltile
303                                                                         !!  @tex $(kg m^{-2})$ @endtex
304!$OMP THREADPRIVATE(tmcs)
305  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmcfc            !! Total moisture content at field capacity per soiltile
306                                                                         !!  @tex $(kg m^{-2})$ @endtex
307!$OMP THREADPRIVATE(tmcfc)
308  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmcw             !! Total moisture content at wilting point per soiltile
309                                                                         !!  @tex $(kg m^{-2})$ @endtex
310!$OMP THREADPRIVATE(tmcw)
311  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter       !! Total moisture in the litter per soiltile
312                                                                         !!  @tex $(kg m^{-2})$ @endtex
313!$OMP THREADPRIVATE(tmc_litter)
314  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tmc_litt_mea     !! Total moisture in the litter over the grid
315                                                                         !!  @tex $(kg m^{-2})$ @endtex
316!$OMP THREADPRIVATE(tmc_litt_mea)
317  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_wilt  !! Total moisture of litter at wilt point per soiltile
318                                                                         !!  @tex $(kg m^{-2})$ @endtex
319!$OMP THREADPRIVATE(tmc_litter_wilt)
320  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_field !! Total moisture of litter at field cap. per soiltile
321                                                                         !!  @tex $(kg m^{-2})$ @endtex
322!$OMP THREADPRIVATE(tmc_litter_field)
323!!! A CHANGER DANS TOUT HYDROL: tmc_litter_res et sat ne devraient pas dependre de ji - tdo
324  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_res   !! Total moisture of litter at residual moisture per soiltile
325                                                                         !!  @tex $(kg m^{-2})$ @endtex
326!$OMP THREADPRIVATE(tmc_litter_res)
327  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_sat   !! Total moisture of litter at saturation per soiltile
328                                                                         !!  @tex $(kg m^{-2})$ @endtex
329!$OMP THREADPRIVATE(tmc_litter_sat)
330  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_awet  !! Total moisture of litter at mc_awet per soiltile
331                                                                         !!  @tex $(kg m^{-2})$ @endtex
332!$OMP THREADPRIVATE(tmc_litter_awet)
333  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tmc_litter_adry  !! Total moisture of litter at mc_adry per soiltile
334                                                                         !!  @tex $(kg m^{-2})$ @endtex
335!$OMP THREADPRIVATE(tmc_litter_adry)
336  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tmc_litt_wet_mea !! Total moisture in the litter over the grid below which
337                                                                         !! albedo is fixed constant
338                                                                         !!  @tex $(kg m^{-2})$ @endtex
339!$OMP THREADPRIVATE(tmc_litt_wet_mea)
340  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: tmc_litt_dry_mea !! Total moisture in the litter over the grid above which
341                                                                         !! albedo is constant
342                                                                         !!  @tex $(kg m^{-2})$ @endtex
343!$OMP THREADPRIVATE(tmc_litt_dry_mea)
344  LOGICAL, SAVE                                      :: tmc_init_updated = .FALSE. !! Flag allowing to determine if tmc is initialized.
345!$OMP THREADPRIVATE(tmc_init_updated)
346
347  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: v1               !! Temporary variable (:)
348!$OMP THREADPRIVATE(v1)
349
350  !! par type de sol :
351  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: ru_ns            !! Surface runoff per soiltile
352                                                                         !!  @tex $(kg m^{-2})$ @endtex
353!$OMP THREADPRIVATE(ru_ns)
354  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: dr_ns            !! Drainage per soiltile
355                                                                         !!  @tex $(kg m^{-2})$ @endtex
356!$OMP THREADPRIVATE(dr_ns)
357  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: tr_ns            !! Transpiration per soiltile
358!$OMP THREADPRIVATE(tr_ns)
359  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: vegetmax_soil    !! (:,nvm,nstm) percentage of each veg. type on each soil
360                                                                         !! of each grid point
361!$OMP THREADPRIVATE(vegetmax_soil)
362  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: mc               !! Total volumetric water content at the calculation nodes
363                                                                         !! (eg : liquid + frozen)
364                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
365!$OMP THREADPRIVATE(mc)
366
367   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_prev       !! Soil moisture from file at previous timestep in the file
368!$OMP THREADPRIVATE(mc_read_prev)
369   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_next       !! Soil moisture from file at next time step in the file
370!$OMP THREADPRIVATE(mc_read_next)
371   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_current    !! For nudging, linear time interpolation bewteen mc_read_prev and mc_read_next
372!$OMP THREADPRIVATE(mc_read_current)
373   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mask_mc_interp     !! Mask of valid data in soil moisture nudging file
374!$OMP THREADPRIVATE(mask_mc_interp)
375   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: tmc_aux            !! Temporary variable needed for the calculation of diag nudgincsm for nudging
376!$OMP THREADPRIVATE(tmc_aux)
377   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowdz_read_prev   !! snowdz read from file at previous timestep in the file
378!$OMP THREADPRIVATE(snowdz_read_prev)
379   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowdz_read_next   !! snowdz read from file at next time step in the file
380!$OMP THREADPRIVATE(snowdz_read_next)
381   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowrho_read_prev  !! snowrho read from file at previous timestep in the file
382!$OMP THREADPRIVATE(snowrho_read_prev)
383   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowrho_read_next  !! snowrho read from file at next time step in the file
384!$OMP THREADPRIVATE(snowrho_read_next)
385   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowtemp_read_prev !! snowtemp read from file at previous timestep in the file
386!$OMP THREADPRIVATE(snowtemp_read_prev)
387   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: snowtemp_read_next !! snowtemp read from file at next time step in the file
388!$OMP THREADPRIVATE(snowtemp_read_next)
389   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)   :: mask_snow_interp   !! Mask of valid data in snow nudging file
390!$OMP THREADPRIVATE(mask_snow_interp)
391
392   REAL(r_std),ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: mcl              !! Liquid water content
393                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
394!$OMP THREADPRIVATE(mcl)
395  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soilmoist        !! (:,nslm) Mean of each soil layer's moisture
396                                                                         !! across soiltiles
397                                                                         !!  @tex $(kg m^{-2})$ @endtex
398!$OMP THREADPRIVATE(soilmoist)
399  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soilmoist_liquid !! (:,nslm) Mean of each soil layer's liquid moisture
400                                                                         !! across soiltiles
401                                                                         !!  @tex $(kg m^{-2})$ @endtex
402!$OMP THREADPRIVATE(soilmoist_liquid)
403  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: soil_wet_ns      !! Soil wetness above mcw (0-1, unitless)
404!$OMP THREADPRIVATE(soil_wet_ns)
405  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soil_wet_litter  !! Soil wetness aove mvw in the litter (0-1, unitless)
406!$OMP THREADPRIVATE(soil_wet_litter)
407  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: qflux_ns         !! Diffusive water fluxes between soil layers
408                                                                         !! (at lower interface)
409!$OMP THREADPRIVATE(qflux_ns)
410  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: check_top_ns     !! Diagnostic calculated in hydrol_diag_soil_flux
411                                                                         !! (water balance residu of top soil layer)
412!$OMP THREADPRIVATE(check_top_ns)
413  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: profil_froz_hydro     !! Frozen fraction for each hydrological soil layer
414!$OMP THREADPRIVATE(profil_froz_hydro)
415  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: profil_froz_hydro_ns  !! As  profil_froz_hydro per soiltile
416!$OMP THREADPRIVATE(profil_froz_hydro_ns)
417
418
419CONTAINS
420
421!! ================================================================================================================================
422!! SUBROUTINE   : hydrol_initialize
423!!
424!>\BRIEF         Allocate module variables, read from restart file or initialize with default values
425!!
426!! DESCRIPTION :
427!!
428!! MAIN OUTPUT VARIABLE(S) :
429!!
430!! REFERENCE(S) :
431!!
432!! FLOWCHART    : None
433!! \n
434!_ ================================================================================================================================
435
436  SUBROUTINE hydrol_initialize ( ks,             nvan,      avan,          mcr,              &
437                                 mcs,            mcfc,      mcw,           kjit,             &
438                                 kjpindex,       index,     rest_id,                         &
439                                 njsc,           soiltile,  veget,         veget_max,        &
440                                 humrel,    vegstress,  drysoil_frac,        &
441                                 shumdiag_perma,    qsintveg,                        &
442                                 evap_bare_lim,  evap_bare_lim_ns,  snow,      snow_age,      snow_nobio,       &
443                                 snow_nobio_age, snowrho,   snowtemp,      snowgrain,        &
444                                 snowdz,         snowheat,  &
445                                 mc_layh,        mcl_layh,  soilmoist_out)
446
447    !! 0. Variable and parameter declaration
448    !! 0.1 Input variables
449
450    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number
451    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
452    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: index            !! Indeces of the points on the map
453    INTEGER(i_std),INTENT (in)                         :: rest_id          !! Restart file identifier
454    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: njsc             !! Index of the dominant soil textural class in the
455                                                                           !! grid cell (1-nscm, unitless) 
456    ! 2D soil parameters
457    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
458    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
459    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
460    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
461    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
462    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
463    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
464   
465    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltile         !! Fraction of each soil tile within vegtot (0-1, unitless)
466    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget            !! Fraction of vegetation type           
467    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget_max        !! Max. fraction of vegetation type (LAI -> infty)
468
469   
470    !! 0.2 Output variables
471    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: humrel         !! Relative humidity
472    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: vegstress      !! Veg. moisture stress (only for vegetation growth)
473    REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: drysoil_frac   !! function of litter wetness
474    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)  :: shumdiag_perma !! Percent of porosity filled with water (mc/mcs) used for the thermal computations
475    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: qsintveg       !! Water on vegetation due to interception
476    REAL(r_std),DIMENSION (kjpindex), INTENT(out)        :: evap_bare_lim  !! Limitation factor for bare soil evaporation
477    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT(out)   :: evap_bare_lim_ns !! Limitation factor for bare soil evaporation
478    REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: snow           !! Snow mass [Kg/m^2]
479    REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: snow_age       !! Snow age
480    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out):: snow_nobio     !! Water balance on ice, lakes, .. [Kg/m^2]
481    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out):: snow_nobio_age !! Snow age on ice, lakes, ...
482    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho        !! Snow density
483    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp       !! Snow temperature
484    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain      !! Snow grainsize
485    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz         !! Snow layer thickness
486    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat       !! Snow heat content
487    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)  :: mc_layh        !! Volumetric moisture content for each layer in hydrol (liquid+ice) m3/m3
488    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)  :: mcl_layh       !! Volumetric moisture content for each layer in hydrol (liquid) m3/m3
489    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)  :: soilmoist_out  !! Total soil moisture content for each layer in hydrol (liquid+ice), mm
490    REAL(r_std),DIMENSION (kjpindex)                     :: soilwetdummy   !! Temporary variable never used
491
492    !! 0.4 Local variables
493    INTEGER(i_std)                                       :: jsl
494   
495!_ ================================================================================================================================
496
497    CALL hydrol_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc, kjit, kjpindex, index, rest_id, veget_max, soiltile, &
498         humrel, vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
499         snowdz, snowgrain, snowrho,    snowtemp,   snowheat, &
500         drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
501   
502    CALL hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget, veget_max, &
503         soiltile, njsc, mx_eau_var, shumdiag_perma, &
504         drysoil_frac, qsintveg, mc_layh, mcl_layh) 
505
506    !! Initialize hydrol_alma routine if the variables were not found in the restart file. This is done in the end of
507    !! hydrol_initialize so that all variables(humtot,..) that will be used are initialized.
508    IF (ALL(tot_watveg_beg(:)==val_exp) .OR.  ALL(tot_watsoil_beg(:)==val_exp) .OR. ALL(snow_beg(:)==val_exp)) THEN
509       ! The output variable soilwetdummy is not calculated at first call to hydrol_alma.
510       CALL hydrol_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwetdummy)
511    END IF
512   
513    !! Calculate itopmax indicating the number of layers where the node is above 0.1m depth
514    itopmax=1
515    DO jsl = 1, nslm
516       ! znh : depth of nodes
517       IF (znh(jsl) <= 0.1) THEN
518          itopmax=jsl
519       END IF
520    END DO
521    IF (printlev>=3) WRITE(numout,*) "Number of layers where the node is above 0.1m depth: itopmax=",itopmax
522
523    ! Copy soilmoist into a local variable to be sent to thermosoil
524    soilmoist_out(:,:) = soilmoist(:,:)
525
526  END SUBROUTINE hydrol_initialize
527
528
529!! ================================================================================================================================
530!! SUBROUTINE   : hydrol_main
531!!
532!>\BRIEF         
533!!
534!! DESCRIPTION :
535!! - called every time step
536!! - initialization and finalization part are not done in here
537!!
538!! - 1 computes snow  ==> explicitsnow
539!! - 2 computes vegetations reservoirs  ==> hydrol_vegupd
540!! - 3 computes canopy  ==> hydrol_canop
541!! - 4 computes surface reservoir  ==> hydrol_flood
542!! - 5 computes soil hydrology ==> hydrol_soil
543!!
544!! IMPORTANT NOTICE : The water fluxes are used in their integrated form, over the time step
545!! dt_sechiba, with a unit of kg m^{-2}.
546!!
547!! RECENT CHANGE(S) : None
548!!
549!! MAIN OUTPUT VARIABLE(S) :
550!!
551!! REFERENCE(S) :
552!!
553!! FLOWCHART    : None
554!! \n
555!_ ================================================================================================================================
556
557  SUBROUTINE hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw,  &
558       & kjit, kjpindex, &
559       & index, indexveg, indexsoil, indexlayer, indexnslm, &
560       & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
561       & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
562       & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
563       & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, evap_bare_lim_ns, &
564       & flood_frac, flood_res, &
565       & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope, rest_id, hist_id, hist2_id,&
566       & contfrac, stempdiag, &
567       & temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
568       & snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
569       & grndflux,gtemp,tot_bare_soil, &
570       & lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
571       & mc_layh, mcl_layh, soilmoist_out )
572
573    !! 0. Variable and parameter declaration
574
575    !! 0.1 Input variables
576 
577    INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number
578    INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
579    INTEGER(i_std),INTENT (in)                         :: rest_id,hist_id  !! _Restart_ file and _history_ file identifier
580    INTEGER(i_std),INTENT (in)                         :: hist2_id         !! _history_ file 2 identifier
581    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: index            !! Indeces of the points on the map
582    INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg        !! Indeces of the points on the 3D map for veg
583    INTEGER(i_std),DIMENSION (kjpindex*nstm), INTENT (in):: indexsoil      !! Indeces of the points on the 3D map for soil
584    INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in):: indexlayer     !! Indeces of the points on the 3D map for soil layers
585    INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in):: indexnslm      !! Indeces of the points on the 3D map for of diagnostic soil layers
586
587    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: precip_rain      !! Rain precipitation
588    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: precip_snow      !! Snow precipitation
589    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: returnflow       !! Routed water which comes back into the soil (from the
590                                                                           !! bottom)
591    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinfiltration   !! Routed water which comes back into the soil (at the
592                                                                           !! top)
593    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: irrigation       !! Water from irrigation returning to soil moisture 
594    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: temp_sol_new     !! New soil temperature
595
596    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)   :: njsc             !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
597    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio     !! Fraction of ice, lakes, ...
598    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: totfrac_nobio    !! Total fraction of ice+lakes+...
599    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: soilcap          !! Soil capacity
600    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltile         !! Fraction of each soil tile within vegtot (0-1, unitless)
601    REAL(r_std),DIMENSION (kjpindex,nlut), INTENT (in) :: fraclut          !! Fraction of each landuse tile (0-1, unitless)
602    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vevapwet         !! Interception loss
603    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget            !! Fraction of vegetation type           
604    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: veget_max        !! Max. fraction of vegetation type (LAI -> infty)
605    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: qsintmax         !! Maximum water on vegetation for interception
606    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration
607    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef
608
609    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
610    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
611    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
612    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
613    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
614    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
615    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
616 
617    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation
618    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot_penm      !! Soil Potential Evaporation Correction
619    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: flood_frac       !! flood fraction
620    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: contfrac         !! Fraction of continent in the grid
621    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in) :: stempdiag        !! Diagnostic temp profile from thermosoil
622    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: temp_air         !! Air temperature
623    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: u,v              !! Horizontal wind speed
624    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: tq_cdrag         !! Surface drag coefficient (-)
625    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: pb               !! Surface pressure
626    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: swnet            !! Net shortwave radiation
627    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: pgflux           !! Net energy into snowpack
628    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: gtemp            !! First soil layer temperature
629    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: tot_bare_soil    !! Total evaporating bare soil fraction
630    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: lambda_snow      !! Coefficient of the linear extrapolation of surface temperature
631    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (in):: cgrnd_snow       !! Integration coefficient for snow numerical scheme
632    REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (in):: dgrnd_snow       !! Integration coefficient for snow numerical scheme
633    REAL(r_std),DIMENSION (kjpindex), INTENT(in)       :: frac_snow_veg    !! Snow cover fraction on vegetation   
634
635    !! 0.2 Output variables
636
637    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress        !! Veg. moisture stress (only for vegetation growth)
638    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: drysoil_frac     !! function of litter wetness
639    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out):: shumdiag         !! Relative soil moisture in each soil layer
640                                                                           !! with respect to (mcfc-mcw)
641                                                                           !! (unitless; can be out of 0-1)
642    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out):: shumdiag_perma   !! Percent of porosity filled with water (mc/mcs) used for the thermal computations
643    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: k_litt           !! litter approximate conductivity
644    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: litterhumdiag    !! litter humidity
645    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: tot_melt         !! Total melt   
646    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: floodout         !! Flux out of floodplains
647   
648    !! 0.3 Modified variables
649
650    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(inout):: qsintveg         !! Water on vegetation due to interception
651    REAL(r_std),DIMENSION (kjpindex), INTENT(inout)    :: evap_bare_lim    !! Limitation factor (beta) for bare soil evaporation
652    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(inout):: evap_bare_lim_ns !! Limitation factor (beta) for bare soil evaporation   
653    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(inout):: humrel           !! Relative humidity
654    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: vevapnu          !! Bare soil evaporation
655    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: vevapsno         !! Snow evaporation
656    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: vevapflo         !! Floodplain evaporation
657    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: flood_res        !! flood reservoir estimate
658    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: snow             !! Snow mass [kg/m^2]
659    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: snow_age         !! Snow age
660    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio  !! Water balance on ice, lakes, .. [Kg/m^2]
661    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
662    !! We consider that any water on the ice is snow and we only peforme a water balance to have consistency.
663    !! The water balance is limite to + or - 10^6 so that accumulation is not endless
664
665    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: runoff       !! Complete surface runoff
666    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: drainage     !! Drainage
667    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho      !! Snow density
668    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp     !! Snow temperature
669    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain    !! Snow grainsize
670    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz       !! Snow layer thickness
671    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat     !! Snow heat content
672    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out)   :: snowliq      !! Snow liquid content (m)
673    REAL(r_std), DIMENSION (kjpindex), INTENT(out)         :: grndflux     !! Net flux into soil W/m2
674    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT(out)     :: mc_layh      !! Volumetric moisture content for each layer in hydrol(liquid + ice) [m3/m3)]
675    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT(out)     :: mcl_layh     !! Volumetric moisture content for each layer in hydrol(liquid) [m3/m3]
676    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT(out)     :: soilmoist_out!! Total soil moisture content for each layer in hydrol(liquid + ice) [mm]
677    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)       :: temp_sol_add !! additional surface temperature due to the melt of first layer
678                                                                           !! at the present time-step @tex ($K$) @endtex
679
680    !! 0.4 Local variables
681    INTEGER(i_std)                                     :: jst              !! Index of soil tiles (unitless, 1-3)
682    INTEGER(i_std)                                     :: jsl              !! Index of soil layers (unitless)
683    INTEGER(i_std)                                     :: ji, jv
684    CHARACTER(LEN=80)                                  :: kfact_root_type  !! read from run.def: when equal to 'cons', it indicates that
685                                                                           !! ks does not increase in the rootzone, ie, kfact_root=1;
686                                                                           !! else, kfact_root defined as usual
687    REAL(r_std),DIMENSION (kjpindex)                   :: soilwet          !! A temporary diagnostic of soil wetness
688    REAL(r_std),DIMENSION (kjpindex)                   :: snowdepth_diag   !! Depth of snow layer containing default values, only for diagnostics
689    REAL(r_std),DIMENSION (kjpindex, nsnow)            :: snowdz_diag      !! Depth of snow layer on all layers containing default values,
690                                                                           !! only for diagnostics
691    REAL(r_std),DIMENSION (kjpindex)                   :: njsc_tmp         !! Temporary REAL value for njsc to write it
692    REAL(r_std), DIMENSION (kjpindex)                  :: snowmelt         !! Snow melt [mm/dt_sechiba]
693    REAL(r_std), DIMENSION (kjpindex,nstm)             :: tmc_top          !! Moisture content in the itopmax upper layers, per tile
694    REAL(r_std), DIMENSION (kjpindex)                  :: humtot_top       !! Moisture content in the itopmax upper layers, for diagnistics
695    REAL(r_std), DIMENSION(kjpindex)                   :: histvar          !! Temporary variable when computations are needed
696    REAL(r_std), DIMENSION (kjpindex,nvm)              :: frac_bare        !! Fraction(of veget_max) of bare soil in each vegetation type
697    INTEGER(i_std), DIMENSION(kjpindex*imax)           :: mc_lin_axis_index
698    REAL(r_std), DIMENSION(kjpindex)                   :: twbr             !! Grid-cell mean of TWBR Total Water Budget Residu[kg/m2/dt]
699    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_nroot       !! To ouput the grid-cell mean of nroot
700    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_dlh         !! To ouput the soil layer thickness on all grid points [m]
701    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_mcs         !! To ouput the mean of mcs
702    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_mcfc        !! To ouput the mean of mcfc
703    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_mcw         !! To ouput the mean of mcw
704    REAL(r_std),DIMENSION (kjpindex,nslm)              :: land_mcr         !! To ouput the mean of mcr
705    REAL(r_std),DIMENSION (kjpindex)                   :: land_tmcs        !! To ouput the grid-cell mean of tmcs
706    REAL(r_std),DIMENSION (kjpindex)                   :: land_tmcfc       !! To ouput the grid-cell mean of tmcfc
707    REAL(r_std),DIMENSION (kjpindex)                   :: drain_upd        !! Change in drainage due to decrease in vegtot
708                                                                           !! on mc [kg/m2/dt]
709    REAL(r_std),DIMENSION (kjpindex)                   :: runoff_upd       !! Change in runoff due to decrease in vegtot
710                                                                           !! on water2infilt[kg/m2/dt]
711    REAL(r_std),DIMENSION (kjpindex)                   :: mrsow            !! Soil wetness above wilting point for CMIP6 (humtot-WP)/(SAT-WP)
712    REAL(r_std), DIMENSION (kjpindex,nlut)             :: humtot_lut       !! Moisture content on landuse tiles, for diagnostics
713    REAL(r_std), DIMENSION (kjpindex,nlut)             :: humtot_top_lut   !! Moisture content in upper layers on landuse tiles, for diagnostics
714    REAL(r_std), DIMENSION (kjpindex,nlut)             :: mrro_lut         !! Total runoff from landuse tiles, for diagnostics
715
716!_ ================================================================================================================================
717    !! 1. Update vegtot_old and recalculate vegtot
718    vegtot_old(:) = vegtot(:)
719
720    DO ji = 1, kjpindex
721       vegtot(ji) = SUM(veget_max(ji,:))
722    ENDDO
723
724
725    !! 2. Applay nudging for soil moisture and/or snow variables
726
727    ! For soil moisture, here only read and interpolate the soil moisture from file to current time step.
728    ! The values will be applayed in hydrol_soil after the soil moisture has been updated.
729    IF (ok_nudge_mc) THEN
730       CALL hydrol_nudge_mc_read(kjit)
731    END IF
732
733    ! Read, interpolate and applay nudging of snow variables
734    IF ( ok_nudge_snow) THEN
735     CALL hydrol_nudge_snow(kjit, kjpindex, snowdz, snowrho, snowtemp )
736    END IF
737
738
739    !! 3. Shared time step
740    IF (printlev>=3) WRITE (numout,*) 'hydrol pas de temps = ',dt_sechiba
741
742    !
743    !! 3.1 Calculate snow processes with explicit snow model
744    CALL explicitsnow_main(kjpindex,    precip_rain,  precip_snow,   temp_air,    pb,       &
745         u,           v,            temp_sol_new,  soilcap,     pgflux,   &
746         frac_nobio,  totfrac_nobio,gtemp,                                &
747         lambda_snow, cgrnd_snow,   dgrnd_snow,    contfrac,              & 
748         vevapsno,    snow_age,     snow_nobio_age,snow_nobio,  snowrho,  &
749         snowgrain,   snowdz,       snowtemp,      snowheat,    snow,     &
750         temp_sol_add,                                                         &
751         snowliq,     subsnownobio, grndflux,      snowmelt,    tot_melt, &
752         subsinksoil)           
753       
754    !
755    !! 3.2 computes vegetations reservoirs  ==>hydrol_vegupd
756! Modif temp vuichard
757    CALL hydrol_vegupd(kjpindex, veget, veget_max, soiltile, qsintveg, frac_bare, drain_upd, runoff_upd)
758
759    !! Calculate kfact_root
760    !! An exponential factor is used to increase ks near the surface depending on the amount of roots in the soil
761    !! through a geometric average over the vegets
762    !! This comes from the PhD thesis of d'Orgeval, 2006, p82; d'Orgeval et al. 2008, Eqs. 3-4
763    !! (Calibrated against Hapex-Sahel measurements)
764    !! Since rev 2916: veget_max/2 is used instead of veget
765    kfact_root(:,:,:) = un
766    DO jsl = 1, nslm
767       DO jv = 2, nvm
768          jst = pref_soil_veg(jv)
769          DO ji = 1, kjpindex
770             IF(soiltile(ji,jst) .GT. min_sechiba) THEN
771                kfact_root(ji,jsl,jst) = kfact_root(ji,jsl,jst) * &
772                     & MAX((MAXVAL(ks_usda)/ks(ji))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &
773                     un) 
774             ENDIF
775          ENDDO
776       ENDDO
777    ENDDO
778
779    !! KFACT_ROOT_TYPE = cons is used to impose that kfact_root = 1 in every soil layer,
780    !! so that ks does not increase in the rootzone; else, kfact_root defined as usual
781    !Config Key   = KFACT_ROOT_TYPE
782    !Config Desc  = keyword added for spmip exp1 and exp4 to get a constant ks over soil depth and rootzone
783    !Config If    = spmip exp1 or exp4
784    !Config Def   = var
785    !Config Help  = can have two values: 'cons' or 'var'. If var then no changes, if cons then kfact_root=un
786    !Config Units = [mm/d]
787    kfact_root_type = 'var'
788    CALL getin_p("KFACT_ROOT_TYPE",kfact_root_type)
789
790    IF (kfact_root_type=='cons') THEN
791       kfact_root(:,:,:) = un
792    ENDIF
793
794    !
795    !! 3.3 computes canopy  ==>hydrol_canop
796    CALL hydrol_canop(kjpindex, precip_rain, vevapwet, veget_max, veget, qsintmax, qsintveg,precisol,tot_melt)
797
798    !
799    !! 3.4 computes surface reservoir  ==>hydrol_flood
800    CALL hydrol_flood(kjpindex,  vevapflo, flood_frac, flood_res, floodout)
801
802    !
803    !! 3.5 computes soil hydrology ==>hydrol_soil
804
805    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
806         transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
807         returnflow, reinfiltration, irrigation, &
808         tot_melt,evap_bare_lim,evap_bare_lim_ns, shumdiag, shumdiag_perma, &
809         k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,&
810         stempdiag,snow,snowdz, tot_bare_soil,  u, v, tq_cdrag, &
811         mc_layh, mcl_layh)
812
813    ! The update fluxes come from hydrol_vegupd
814    drainage(:) =  drainage(:) +  drain_upd(:)
815    runoff(:) =  runoff(:) +  runoff_upd(:)
816
817
818    !! 4 write out file  ==> hydrol_alma/histwrite(*)
819    !
820    ! If we use the ALMA standards
821    CALL hydrol_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet)
822   
823
824    ! Calculate the moisture in the upper itopmax layers corresponding to 0.1m (humtot_top):
825    ! For ORCHIDEE with nslm=11 and zmaxh=2, itopmax=6.
826    ! We compute tmc_top as tmc but only for the first itopmax layers. Then we compute a humtot with this variable.
827    DO jst=1,nstm
828       DO ji=1,kjpindex
829          tmc_top(ji,jst) = dz(2) * ( trois*mc(ji,1,jst) + mc(ji,2,jst) )/huit
830          DO jsl = 2, itopmax
831             tmc_top(ji,jst) = tmc_top(ji,jst) + dz(jsl) * (trois*mc(ji,jsl,jst)+mc(ji,jsl-1,jst))/huit &
832                  + dz(jsl+1) * (trois*mc(ji,jsl,jst)+mc(ji,jsl+1,jst))/huit
833          ENDDO
834       ENDDO
835    ENDDO
836 
837    ! We average the values of each soiltile and multiply by vegtot to transform to a grid-cell mean
838    humtot_top(:) = zero
839    DO jst=1,nstm
840       DO ji=1,kjpindex
841          humtot_top(ji) = humtot_top(ji) + soiltile(ji,jst) * tmc_top(ji,jst) * vegtot(ji)
842       ENDDO
843    ENDDO
844
845    ! Calculate the Total Water Budget Residu (in kg/m2 over dt_sechiba)
846    ! All the delstocks and fluxes below are averaged over the mesh
847    ! snow_nobio included in delswe
848    ! Does not include the routing reservoirs, although the flux to/from routing are integrated
849    DO ji=1,kjpindex
850       twbr(ji) = (delsoilmoist(ji) + delintercept(ji) + delswe(ji)) &
851            - ( precip_rain(ji) + precip_snow(ji) + irrigation(ji) + floodout(ji) &
852            + returnflow(ji) + reinfiltration(ji) ) &
853            + ( runoff(ji) + drainage(ji) + SUM(vevapwet(ji,:)) &
854            + SUM(transpir(ji,:)) + vevapnu(ji) + vevapsno(ji) + vevapflo(ji) ) 
855    ENDDO
856    ! Transform unit from kg/m2/dt to kg/m2/s (or mm/s)
857    CALL xios_orchidee_send_field("twbr",twbr/dt_sechiba)
858    CALL xios_orchidee_send_field("undermcr",undermcr) ! nb of tiles undermcr at end of timestep
859
860    ! Calculate land_nroot : grid-cell mean of nroot
861    ! Do not treat PFT1 because it has no roots
862    land_nroot(:,:) = zero
863    DO jsl=1,nslm
864       DO jv=2,nvm
865          DO ji=1,kjpindex
866               IF ( vegtot(ji) > min_sechiba ) THEN
867               land_nroot(ji,jsl) = land_nroot(ji,jsl) + veget_max(ji,jv) * nroot(ji,jv,jsl) / vegtot(ji) 
868            END IF
869          END DO
870       ENDDO
871    ENDDO
872    CALL xios_orchidee_send_field("nroot",land_nroot)   
873
874    DO jsl=1,nslm
875       land_dlh(:,jsl)=dlh(jsl)
876    ENDDO
877    CALL xios_orchidee_send_field("dlh",land_dlh)
878
879    ! Particular soil moisture values, spatially averaged over the grid-cell
880    ! (a) total SM in kg/m2
881    !     we average the total values of each soiltile and multiply by vegtot to transform to a grid-cell mean (over total land)
882    land_tmcs(:) = zero
883    land_tmcfc(:) = zero
884    DO jst=1,nstm
885       DO ji=1,kjpindex
886          land_tmcs(ji) = land_tmcs(ji) + soiltile(ji,jst) * tmcs(ji,jst) * vegtot(ji)
887          land_tmcfc(ji) = land_tmcfc(ji) + soiltile(ji,jst) * tmcfc(ji,jst) * vegtot(ji)
888       ENDDO
889    ENDDO
890    CALL xios_orchidee_send_field("tmcs",land_tmcs) ! in kg/m2
891    CALL xios_orchidee_send_field("tmcfc",land_tmcfc) ! in kg/m2
892
893    ! (b) volumetric moisture content by layers in m3/m3
894    !     mcs etc are identical in all layers (no normalization by vegtot to be comparable to mc)
895    DO jsl=1,nslm
896       land_mcs(:,jsl) = mcs(:)
897       land_mcfc(:,jsl) = mcfc(:)
898       land_mcw(:,jsl) = mcw(:)
899       land_mcr(:,jsl) = mcr(:)
900    ENDDO
901    CALL xios_orchidee_send_field("mcs",land_mcs) ! in m3/m3
902    CALL xios_orchidee_send_field("mcfc",land_mcfc) ! in m3/m3
903    CALL xios_orchidee_send_field("mcw",land_mcw) ! in m3/m3
904    CALL xios_orchidee_send_field("mcr",land_mcr) ! in m3/m3
905
906     
907    CALL xios_orchidee_send_field("water2infilt",water2infilt)   
908    CALL xios_orchidee_send_field("mc",mc)
909    CALL xios_orchidee_send_field("kfact_root",kfact_root)
910    CALL xios_orchidee_send_field("vegetmax_soil",vegetmax_soil)
911    CALL xios_orchidee_send_field("evapnu_soil",ae_ns/dt_sechiba)
912    CALL xios_orchidee_send_field("drainage_soil",dr_ns/dt_sechiba)
913    CALL xios_orchidee_send_field("transpir_soil",tr_ns/dt_sechiba)
914    CALL xios_orchidee_send_field("runoff_soil",ru_ns/dt_sechiba)
915    CALL xios_orchidee_send_field("humrel",humrel)     
916    CALL xios_orchidee_send_field("drainage",drainage/dt_sechiba) ! [kg m-2 s-1]
917    CALL xios_orchidee_send_field("runoff",runoff/dt_sechiba) ! [kg m-2 s-1]
918    CALL xios_orchidee_send_field("precisol",precisol/dt_sechiba)
919    CALL xios_orchidee_send_field("throughfall",throughfall/dt_sechiba)
920    CALL xios_orchidee_send_field("precip_rain",precip_rain/dt_sechiba)
921    CALL xios_orchidee_send_field("precip_snow",precip_snow/dt_sechiba)
922    CALL xios_orchidee_send_field("qsintmax",qsintmax)
923    CALL xios_orchidee_send_field("qsintveg",qsintveg)
924    CALL xios_orchidee_send_field("qsintveg_tot",SUM(qsintveg(:,:),dim=2))
925    histvar(:)=(precip_rain(:)-SUM(throughfall(:,:),dim=2))
926    CALL xios_orchidee_send_field("prveg",histvar/dt_sechiba)
927
928    IF ( do_floodplains ) THEN
929       CALL xios_orchidee_send_field("floodout",floodout/dt_sechiba)
930    END IF
931
932    CALL xios_orchidee_send_field("snowmelt",snowmelt/dt_sechiba)
933    CALL xios_orchidee_send_field("tot_melt",tot_melt/dt_sechiba)
934
935    CALL xios_orchidee_send_field("soilmoist",soilmoist)
936    CALL xios_orchidee_send_field("soilmoist_liquid",soilmoist_liquid)
937    CALL xios_orchidee_send_field("humtot_frozen",SUM(soilmoist(:,:),2)-SUM(soilmoist_liquid(:,:),2))
938    CALL xios_orchidee_send_field("tmc",tmc)
939    CALL xios_orchidee_send_field("humtot",humtot)
940    CALL xios_orchidee_send_field("humtot_top",humtot_top)
941
942    ! For the soil wetness above wilting point for CMIP6 (mrsow)
943    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(:) ) &
944         / ( zmaxh*mille*( mcs(:) - mcw(:) ) )
945    CALL xios_orchidee_send_field("mrsow",mrsow)
946
947
948   
949    ! Prepare diagnostic snow variables
950    !  Add XIOS default value where no snow
951    DO ji=1,kjpindex
952       IF (snow(ji) > 0) THEN
953          snowdz_diag(ji,:) = snowdz(ji,:)
954          snowdepth_diag(ji) = SUM(snowdz(ji,:))*(1-totfrac_nobio(ji))*frac_snow_veg(ji)
955       ELSE
956          snowdz_diag(ji,:) = xios_default_val
957          snowdepth_diag(ji) = xios_default_val             
958       END IF
959    END DO
960    CALL xios_orchidee_send_field("snowdz",snowdz_diag)
961    CALL xios_orchidee_send_field("snowdepth",snowdepth_diag)
962
963    CALL xios_orchidee_send_field("frac_bare",frac_bare)
964    CALL xios_orchidee_send_field("soilwet",soilwet)
965    CALL xios_orchidee_send_field("delsoilmoist",delsoilmoist)
966    CALL xios_orchidee_send_field("delswe",delswe)
967    CALL xios_orchidee_send_field("delintercept",delintercept) 
968
969    IF (ok_freeze_cwrr) THEN
970       CALL xios_orchidee_send_field("profil_froz_hydro",profil_froz_hydro)
971    END IF
972    CALL xios_orchidee_send_field("profil_froz_hydro_ns", profil_froz_hydro_ns)
973    CALL xios_orchidee_send_field("kk_moy",kk_moy) ! in mm/d
974
975    !! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
976    humtot_lut(:,:)=0
977    humtot_top_lut(:,:)=0
978    mrro_lut(:,:)=0
979    DO jv=1,nvm
980       jst=pref_soil_veg(jv) ! soil tile index
981       IF (natural(jv)) THEN
982          humtot_lut(:,id_psl) = humtot_lut(:,id_psl) + tmc(:,jst)*veget_max(:,jv)
983          humtot_top_lut(:,id_psl) = humtot_top_lut(:,id_psl) + tmc_top(:,jst)*veget_max(:,jv)
984          mrro_lut(:,id_psl) = mrro_lut(:,id_psl) + (dr_ns(:,jst)+ru_ns(:,jst))*veget_max(:,jv)
985       ELSE
986          humtot_lut(:,id_crp) = humtot_lut(:,id_crp) + tmc(:,jst)*veget_max(:,jv)
987          humtot_top_lut(:,id_crp) = humtot_top_lut(:,id_crp) + tmc_top(:,jst)*veget_max(:,jv)
988          mrro_lut(:,id_crp) = mrro_lut(:,id_crp) + (dr_ns(:,jst)+ru_ns(:,jst))*veget_max(:,jv)
989       ENDIF
990    END DO
991
992    WHERE (fraclut(:,id_psl)>min_sechiba)
993       humtot_lut(:,id_psl) = humtot_lut(:,id_psl)/fraclut(:,id_psl)
994       humtot_top_lut(:,id_psl) = humtot_top_lut(:,id_psl)/fraclut(:,id_psl)
995       mrro_lut(:,id_psl) = mrro_lut(:,id_psl)/fraclut(:,id_psl)/dt_sechiba
996    ELSEWHERE
997       humtot_lut(:,id_psl) = val_exp
998       humtot_top_lut(:,id_psl) = val_exp
999       mrro_lut(:,id_psl) = val_exp
1000    END WHERE
1001    WHERE (fraclut(:,id_crp)>min_sechiba)
1002       humtot_lut(:,id_crp) = humtot_lut(:,id_crp)/fraclut(:,id_crp)
1003       humtot_top_lut(:,id_crp) = humtot_top_lut(:,id_crp)/fraclut(:,id_crp)
1004       mrro_lut(:,id_crp) = mrro_lut(:,id_crp)/fraclut(:,id_crp)/dt_sechiba
1005    ELSEWHERE
1006       humtot_lut(:,id_crp) = val_exp
1007       humtot_top_lut(:,id_crp) = val_exp
1008       mrro_lut(:,id_crp) = val_exp
1009    END WHERE
1010
1011    humtot_lut(:,id_pst) = val_exp
1012    humtot_lut(:,id_urb) = val_exp
1013    humtot_top_lut(:,id_pst) = val_exp
1014    humtot_top_lut(:,id_urb) = val_exp
1015    mrro_lut(:,id_pst) = val_exp
1016    mrro_lut(:,id_urb) = val_exp
1017
1018    CALL xios_orchidee_send_field("humtot_lut",humtot_lut)
1019    CALL xios_orchidee_send_field("humtot_top_lut",humtot_top_lut)
1020    CALL xios_orchidee_send_field("mrro_lut",mrro_lut)
1021
1022    ! Write diagnistic for soil moisture nudging
1023    IF (ok_nudge_mc) CALL hydrol_nudge_mc_diag(kjpindex, soiltile)
1024
1025
1026    IF ( .NOT. almaoutput ) THEN
1027       CALL histwrite_p(hist_id, 'frac_bare', kjit, frac_bare, kjpindex*nvm, indexveg)
1028
1029       DO jst=1,nstm
1030          ! var_name= "mc_1" ... "mc_3"
1031          WRITE (var_name,"('moistc_',i1)") jst
1032          CALL histwrite_p(hist_id, TRIM(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
1033
1034          ! var_name= "kfactroot_1" ... "kfactroot_3"
1035          WRITE (var_name,"('kfactroot_',i1)") jst
1036          CALL histwrite_p(hist_id, TRIM(var_name), kjit, kfact_root(:,:,jst), kjpindex*nslm, indexlayer)
1037
1038          ! var_name= "vegetsoil_1" ... "vegetsoil_3"
1039          WRITE (var_name,"('vegetsoil_',i1)") jst
1040          CALL histwrite_p(hist_id, TRIM(var_name), kjit,vegetmax_soil(:,:,jst), kjpindex*nvm, indexveg)
1041       ENDDO
1042       CALL histwrite_p(hist_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
1043       CALL histwrite_p(hist_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
1044       CALL histwrite_p(hist_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
1045       CALL histwrite_p(hist_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
1046       CALL histwrite_p(hist_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
1047       ! mrso is a perfect duplicate of humtot
1048       CALL histwrite_p(hist_id, 'humtot', kjit, humtot, kjpindex, index)
1049       CALL histwrite_p(hist_id, 'mrso', kjit, humtot, kjpindex, index)
1050       CALL histwrite_p(hist_id, 'mrsos', kjit, humtot_top, kjpindex, index)
1051       njsc_tmp(:)=njsc(:)
1052       CALL histwrite_p(hist_id, 'soilindex', kjit, njsc_tmp, kjpindex, index)
1053       CALL histwrite_p(hist_id, 'humrel',   kjit, humrel,   kjpindex*nvm, indexveg)
1054       CALL histwrite_p(hist_id, 'drainage', kjit, drainage, kjpindex, index)
1055       ! NB! According to histdef in intersurf, the variables 'runoff' and 'mrros' have different units
1056       CALL histwrite_p(hist_id, 'runoff', kjit, runoff, kjpindex, index)
1057       CALL histwrite_p(hist_id, 'mrros', kjit, runoff, kjpindex, index)
1058       histvar(:)=(runoff(:)+drainage(:))
1059       CALL histwrite_p(hist_id, 'mrro', kjit, histvar, kjpindex, index)
1060       CALL histwrite_p(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
1061       CALL histwrite_p(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
1062
1063       histvar(:)=(precip_rain(:)-SUM(throughfall(:,:),dim=2))
1064       CALL histwrite_p(hist_id, 'prveg', kjit, histvar, kjpindex, index)
1065
1066       CALL histwrite_p(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
1067       CALL histwrite_p(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
1068       CALL histwrite_p(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
1069       CALL histwrite_p(hist_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
1070       CALL histwrite_p(hist_id, 'shumdiag_perma',kjit,shumdiag_perma,kjpindex*nslm,indexnslm)
1071
1072       IF ( do_floodplains ) THEN
1073          CALL histwrite_p(hist_id, 'floodout', kjit, floodout, kjpindex, index)
1074       ENDIF
1075       !
1076       IF ( hist2_id > 0 ) THEN
1077          DO jst=1,nstm
1078             ! var_name= "mc_1" ... "mc_3"
1079             WRITE (var_name,"('moistc_',i1)") jst
1080             CALL histwrite_p(hist2_id, TRIM(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
1081
1082             ! var_name= "kfactroot_1" ... "kfactroot_3"
1083             WRITE (var_name,"('kfactroot_',i1)") jst
1084             CALL histwrite_p(hist2_id, TRIM(var_name), kjit, kfact_root(:,:,jst), kjpindex*nslm, indexlayer)
1085
1086             ! var_name= "vegetsoil_1" ... "vegetsoil_3"
1087             WRITE (var_name,"('vegetsoil_',i1)") jst
1088             CALL histwrite_p(hist2_id, TRIM(var_name), kjit,vegetmax_soil(:,:,jst), kjpindex*nvm, indexveg)
1089          ENDDO
1090          CALL histwrite_p(hist2_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
1091          CALL histwrite_p(hist2_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
1092          CALL histwrite_p(hist2_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
1093          CALL histwrite_p(hist2_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
1094          CALL histwrite_p(hist2_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
1095          ! mrso is a perfect duplicate of humtot
1096          CALL histwrite_p(hist2_id, 'humtot', kjit, humtot, kjpindex, index)
1097          CALL histwrite_p(hist2_id, 'mrso', kjit, humtot, kjpindex, index)
1098          CALL histwrite_p(hist2_id, 'mrsos', kjit, humtot_top, kjpindex, index)
1099          njsc_tmp(:)=njsc(:)
1100          CALL histwrite_p(hist2_id, 'soilindex', kjit, njsc_tmp, kjpindex, index)
1101          CALL histwrite_p(hist2_id, 'humrel',   kjit, humrel,   kjpindex*nvm, indexveg)
1102          CALL histwrite_p(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
1103          ! NB! According to histdef in intersurf, the variables 'runoff' and 'mrros' have different units
1104          CALL histwrite_p(hist2_id, 'runoff', kjit, runoff, kjpindex, index)
1105          CALL histwrite_p(hist2_id, 'mrros', kjit, runoff, kjpindex, index)
1106          histvar(:)=(runoff(:)+drainage(:))
1107          CALL histwrite_p(hist2_id, 'mrro', kjit, histvar, kjpindex, index)
1108
1109          IF ( do_floodplains ) THEN
1110             CALL histwrite_p(hist2_id, 'floodout', kjit, floodout, kjpindex, index)
1111          ENDIF
1112          CALL histwrite_p(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
1113          CALL histwrite_p(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
1114          CALL histwrite_p(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
1115          CALL histwrite_p(hist2_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
1116          CALL histwrite_p(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
1117          CALL histwrite_p(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
1118       ENDIF
1119    ELSE
1120       CALL histwrite_p(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
1121       CALL histwrite_p(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
1122       CALL histwrite_p(hist_id, 'Qs', kjit, runoff, kjpindex, index)
1123       CALL histwrite_p(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
1124       CALL histwrite_p(hist_id, 'Qsm', kjit, snowmelt, kjpindex, index)
1125       CALL histwrite_p(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
1126       CALL histwrite_p(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
1127       CALL histwrite_p(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
1128       !
1129       CALL histwrite_p(hist_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
1130       CALL histwrite_p(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
1131       !
1132       CALL histwrite_p(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
1133       CALL histwrite_p(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
1134
1135       IF ( hist2_id > 0 ) THEN
1136          CALL histwrite_p(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
1137          CALL histwrite_p(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
1138          CALL histwrite_p(hist2_id, 'Qs', kjit, runoff, kjpindex, index)
1139          CALL histwrite_p(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
1140          CALL histwrite_p(hist2_id, 'Qsm', kjit, snowmelt, kjpindex, index)
1141          CALL histwrite_p(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
1142          CALL histwrite_p(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
1143          CALL histwrite_p(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
1144          !
1145          CALL histwrite_p(hist2_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
1146          CALL histwrite_p(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
1147          !
1148          CALL histwrite_p(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
1149          CALL histwrite_p(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
1150       ENDIF
1151    ENDIF
1152
1153    IF (ok_freeze_cwrr) THEN
1154       CALL histwrite_p(hist_id, 'profil_froz_hydro', kjit,profil_froz_hydro , kjpindex*nslm, indexlayer)
1155    ENDIF
1156    CALL histwrite_p(hist_id, 'kk_moy', kjit, kk_moy,kjpindex*nslm, indexlayer) ! averaged over soiltiles
1157    DO jst=1,nstm
1158       WRITE (var_name,"('profil_froz_hydro_',i1)") jst
1159       CALL histwrite_p(hist_id, TRIM(var_name), kjit, profil_froz_hydro_ns(:,:,jst), kjpindex*nslm, indexlayer)
1160    ENDDO
1161
1162    ! Copy soilmoist into a local variable to be sent to thermosoil
1163    soilmoist_out(:,:) = soilmoist(:,:)
1164
1165    IF (printlev>=3) WRITE (numout,*) ' hydrol_main Done '
1166
1167  END SUBROUTINE hydrol_main
1168
1169
1170!! ================================================================================================================================
1171!! SUBROUTINE   : hydrol_finalize
1172!!
1173!>\BRIEF         
1174!!
1175!! DESCRIPTION : This subroutine writes the module variables and variables calculated in hydrol to restart file
1176!!
1177!! MAIN OUTPUT VARIABLE(S) :
1178!!
1179!! REFERENCE(S) :
1180!!
1181!! FLOWCHART    : None
1182!! \n
1183!_ ================================================================================================================================
1184
1185  SUBROUTINE hydrol_finalize( kjit,           kjpindex,   rest_id,  vegstress,  &
1186                              qsintveg,       humrel,     snow,     snow_age, snow_nobio, &
1187                              snow_nobio_age, snowrho,    snowtemp, snowdz,     &
1188                              snowheat,       snowgrain,  &
1189                              drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
1190
1191    !! 0. Variable and parameter declaration
1192    !! 0.1 Input variables
1193    INTEGER(i_std), INTENT(in)                           :: kjit           !! Time step number
1194    INTEGER(i_std), INTENT(in)                           :: kjpindex       !! Domain size
1195    INTEGER(i_std),INTENT (in)                           :: rest_id        !! Restart file identifier
1196    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: vegstress      !! Veg. moisture stress (only for vegetation growth)
1197    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintveg       !! Water on vegetation due to interception
1198    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: humrel
1199    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow           !! Snow mass [Kg/m^2]
1200    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow_age       !! Snow age
1201    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio     !! Water balance on ice, lakes, .. [Kg/m^2]
1202    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio_age !! Snow age on ice, lakes, ...
1203    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowrho        !! Snow density
1204    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowtemp       !! Snow temperature
1205    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowdz         !! Snow layer thickness
1206    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowheat       !! Snow heat content
1207    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowgrain      !! Snow grainsize
1208    REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: drysoil_frac   !! function of litter wetness
1209    REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: evap_bare_lim
1210    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(in)     :: evap_bare_lim_ns
1211
1212    !! 0.4 Local variables
1213    INTEGER(i_std)                                       :: jst, jsl
1214   
1215!_ ================================================================================================================================
1216
1217
1218    IF (printlev>=3) WRITE (numout,*) 'Write restart file with HYDROLOGIC variables '
1219
1220    DO jst=1,nstm
1221       ! var_name= "mc_1" ... "mc_3"
1222       WRITE (var_name,"('moistc_',i1)") jst
1223       CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mc(:,:,jst), 'scatter',  nbp_glo, index_g)
1224    END DO
1225
1226    DO jst=1,nstm
1227       ! var_name= "mcl_1" ... "mcl_3"
1228       WRITE (var_name,"('moistcl_',i1)") jst
1229       CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mcl(:,:,jst), 'scatter',  nbp_glo, index_g)
1230    END DO
1231   
1232    IF (ok_nudge_mc) THEN
1233       DO jst=1,nstm
1234          WRITE (var_name,"('mc_read_next_',i1)") jst
1235          CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mc_read_next(:,:,jst), 'scatter',  nbp_glo, index_g)
1236       END DO
1237    END IF
1238
1239    IF (ok_nudge_snow) THEN
1240       CALL restput_p(rest_id, 'snowdz_read_next', nbp_glo,  nsnow, 1, kjit, snowdz_read_next(:,:), &
1241            'scatter',  nbp_glo, index_g)
1242       CALL restput_p(rest_id, 'snowrho_read_next', nbp_glo,  nsnow, 1, kjit, snowrho_read_next(:,:), &
1243            'scatter',  nbp_glo, index_g)
1244       CALL restput_p(rest_id, 'snowtemp_read_next', nbp_glo,  nsnow, 1, kjit, snowtemp_read_next(:,:), &
1245            'scatter',  nbp_glo, index_g)
1246    END IF
1247
1248
1249           
1250    DO jst=1,nstm
1251       DO jsl=1,nslm
1252          ! var_name= "us_1_01" ... "us_3_11"
1253          WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
1254          CALL restput_p(rest_id, var_name, nbp_glo,nvm, 1,kjit,us(:,:,jst,jsl),'scatter',nbp_glo,index_g)
1255       END DO
1256    END DO
1257   
1258    CALL restput_p(rest_id, 'free_drain_coef', nbp_glo,   nstm, 1, kjit,  free_drain_coef, 'scatter',  nbp_glo, index_g)
1259    CALL restput_p(rest_id, 'zwt_force', nbp_glo,   nstm, 1, kjit,  zwt_force, 'scatter',  nbp_glo, index_g)
1260    CALL restput_p(rest_id, 'water2infilt', nbp_glo,   nstm, 1, kjit,  water2infilt, 'scatter',  nbp_glo, index_g)
1261    CALL restput_p(rest_id, 'ae_ns', nbp_glo,   nstm, 1, kjit,  ae_ns, 'scatter',  nbp_glo, index_g)
1262    CALL restput_p(rest_id, 'vegstress', nbp_glo,   nvm, 1, kjit,  vegstress, 'scatter',  nbp_glo, index_g)
1263    CALL restput_p(rest_id, 'snow', nbp_glo,   1, 1, kjit,  snow, 'scatter',  nbp_glo, index_g)
1264    CALL restput_p(rest_id, 'snow_age', nbp_glo,   1, 1, kjit,  snow_age, 'scatter',  nbp_glo, index_g)
1265    CALL restput_p(rest_id, 'snow_nobio', nbp_glo,   nnobio, 1, kjit,  snow_nobio, 'scatter', nbp_glo, index_g)
1266    CALL restput_p(rest_id, 'snow_nobio_age', nbp_glo,   nnobio, 1, kjit,  snow_nobio_age, 'scatter', nbp_glo, index_g)
1267    CALL restput_p(rest_id, 'qsintveg', nbp_glo, nvm, 1, kjit,  qsintveg, 'scatter',  nbp_glo, index_g)
1268    CALL restput_p(rest_id, 'evap_bare_lim_ns', nbp_glo, nstm, 1, kjit,  evap_bare_lim_ns, 'scatter',  nbp_glo, index_g)
1269    CALL restput_p(rest_id, 'evap_bare_lim', nbp_glo, 1, 1, kjit,  evap_bare_lim, 'scatter',  nbp_glo, index_g)
1270    CALL restput_p(rest_id, 'resdist', nbp_glo, nstm, 1, kjit,  resdist, 'scatter',  nbp_glo, index_g) 
1271    CALL restput_p(rest_id, 'vegtot_old', nbp_glo, 1, 1, kjit,  vegtot_old, 'scatter',  nbp_glo, index_g)           
1272    CALL restput_p(rest_id, 'drysoil_frac', nbp_glo,   1, 1, kjit, drysoil_frac, 'scatter', nbp_glo, index_g)
1273    CALL restput_p(rest_id, 'humrel', nbp_glo,   nvm, 1, kjit,  humrel, 'scatter',  nbp_glo, index_g)
1274
1275    CALL restput_p(rest_id, 'tot_watveg_beg', nbp_glo,  1, 1, kjit,  tot_watveg_beg, 'scatter',  nbp_glo, index_g)
1276    CALL restput_p(rest_id, 'tot_watsoil_beg', nbp_glo, 1, 1, kjit,  tot_watsoil_beg, 'scatter',  nbp_glo, index_g)
1277    CALL restput_p(rest_id, 'snow_beg', nbp_glo,        1, 1, kjit,  snow_beg, 'scatter',  nbp_glo, index_g)
1278   
1279   
1280    ! Write variables for explictsnow module to restart file
1281    CALL explicitsnow_finalize ( kjit,     kjpindex, rest_id,    snowrho,   &
1282         snowtemp, snowdz,   snowheat,   snowgrain)
1283
1284  END SUBROUTINE hydrol_finalize
1285
1286
1287!! ================================================================================================================================
1288!! SUBROUTINE   : hydrol_init
1289!!
1290!>\BRIEF        Initializations and memory allocation   
1291!!
1292!! DESCRIPTION  :
1293!! - 1 Some initializations
1294!! - 2 make dynamic allocation with good dimension
1295!! - 2.1 array allocation for soil textur
1296!! - 2.2 Soil texture choice
1297!! - 3 Other array allocation
1298!! - 4 Open restart input file and read data for HYDROLOGIC process
1299!! - 5 get restart values if none were found in the restart file
1300!! - 6 Vegetation array     
1301!! - 7 set humrelv from us
1302!!
1303!! RECENT CHANGE(S) : None
1304!!
1305!! MAIN OUTPUT VARIABLE(S) :
1306!!
1307!! REFERENCE(S) :
1308!!
1309!! FLOWCHART    : None
1310!! \n
1311!_ ================================================================================================================================
1312!!_ hydrol_init
1313
1314  SUBROUTINE hydrol_init(ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc,&
1315       kjit, kjpindex, index, rest_id, veget_max, soiltile, &
1316       humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
1317       snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, &
1318       drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
1319   
1320
1321    !! 0. Variable and parameter declaration
1322
1323    !! 0.1 Input variables
1324    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc               !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
1325    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number
1326    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
1327    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map
1328    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier
1329    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max          !! Carte de vegetation max
1330    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in)  :: soiltile           !! Fraction of each soil tile within vegtot (0-1, unitless)
1331   
1332    !! 0.2 Output variables
1333
1334    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: ks               !! Hydraulic conductivity at saturation (mm {-1})
1335    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: nvan             !! Van Genuchten coeficients n (unitless)
1336    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: avan             !! Van Genuchten coeficients a (mm-1})
1337    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
1338    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
1339    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
1340    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
1341
1342    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: humrel             !! Stress hydrique, relative humidity
1343    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: vegstress          !! Veg. moisture stress (only for vegetation growth)
1344    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: snow               !! Snow mass [Kg/m^2]
1345    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: snow_age           !! Snow age
1346    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio       !! Snow on ice, lakes, ...
1347    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age   !! Snow age on ice, lakes, ...
1348    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)    :: qsintveg         !! Water on vegetation due to interception
1349    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowdz           !! Snow depth
1350    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowgrain        !! Snow grain size
1351    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowheat         !! Snow heat content
1352    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowtemp         !! Snow temperature
1353    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowrho          !! Snow density
1354    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: drysoil_frac     !! function of litter wetness
1355    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: evap_bare_lim
1356    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(out)     :: evap_bare_lim_ns
1357
1358    !! 0.4 Local variables
1359
1360    INTEGER(i_std)                                     :: ier                   !! Error code
1361    INTEGER(i_std)                                     :: ji                    !! Index of land grid cells (1)
1362    INTEGER(i_std)                                     :: jv                    !! Index of PFTs (1)
1363    INTEGER(i_std)                                     :: jst                   !! Index of soil tiles (1)
1364    INTEGER(i_std)                                     :: jsl                   !! Index of soil layers (1)
1365    INTEGER(i_std)                                     :: jsc                   !! Index of soil texture (1)
1366    INTEGER(i_std), PARAMETER                          :: error_level = 3       !! Error level for consistency check
1367    !! Switch to 2 tu turn fatal errors into warnings
1368    REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: free_drain_max        !! Temporary var for initialization of free_drain_coef
1369    REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: zwt_default           !! Temporary variable for initialization of zwt_force
1370    LOGICAL                                            :: zforce                !! To test if we force the WT in any of the soiltiles
1371   
1372
1373!_ ================================================================================================================================
1374
1375    !! 1 Some initializations
1376    !
1377    !Config Key   = DO_PONDS
1378    !Config Desc  = Should we include ponds
1379    !Config Def   = n
1380    !Config If    =
1381    !Config Help  = This parameters allows the user to ask the model
1382    !Config         to take into account the ponds and return
1383    !Config         the water into the soil moisture. If this is
1384    !Config         activated, then there is no reinfiltration
1385    !Config         computed inside the hydrol module.
1386    !Config Units = [FLAG]
1387    !
1388    doponds = .FALSE.
1389    CALL getin_p('DO_PONDS', doponds)
1390
1391    !Config Key   = FROZ_FRAC_CORR
1392    !Config Desc  = Coefficient for the frozen fraction correction
1393    !Config Def   = 1.0
1394    !Config If    = OK_FREEZE
1395    !Config Help  =
1396    !Config Units = [-]
1397    froz_frac_corr = 1.0
1398    CALL getin_p("FROZ_FRAC_CORR", froz_frac_corr)
1399
1400    !Config Key   = MAX_FROZ_HYDRO
1401    !Config Desc  = Coefficient for the frozen fraction correction
1402    !Config Def   = 1.0
1403    !Config If    = OK_FREEZE
1404    !Config Help  =
1405    !Config Units = [-]
1406    max_froz_hydro = 1.0
1407    CALL getin_p("MAX_FROZ_HYDRO", max_froz_hydro)
1408
1409    !Config Key   = SMTOT_CORR
1410    !Config Desc  = Coefficient for the frozen fraction correction
1411    !Config Def   = 2.0
1412    !Config If    = OK_FREEZE
1413    !Config Help  =
1414    !Config Units = [-]
1415    smtot_corr = 2.0
1416    CALL getin_p("SMTOT_CORR", smtot_corr)
1417
1418    !Config Key   = DO_RSOIL
1419    !Config Desc  = Should we reduce soil evaporation with a soil resistance
1420    !Config Def   = n
1421    !Config If    =
1422    !Config Help  = This parameters allows the user to ask the model
1423    !Config         to calculate a soil resistance to reduce the soil evaporation
1424    !Config Units = [FLAG]
1425    !
1426    do_rsoil = .FALSE.
1427    CALL getin_p('DO_RSOIL', do_rsoil) 
1428
1429    !Config Key   = OK_DYNROOT
1430    !Config Desc  = Calculate dynamic root profile to optimize soil moisture usage 
1431    !Config Def   = n
1432    !Config If    =
1433    !Config Help  =
1434    !Config Units = [FLAG]
1435    ok_dynroot = .FALSE.
1436    CALL getin_p('OK_DYNROOT',ok_dynroot)
1437
1438    !! 2 make dynamic allocation with good dimension
1439
1440    !! 2.1 array allocation for soil texture
1441
1442    ALLOCATE (pcent(nscm),stat=ier)
1443    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','')
1444   
1445    ALLOCATE (mc_awet(nscm),stat=ier)
1446    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','')
1447
1448    ALLOCATE (mc_adry(nscm),stat=ier)
1449    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','')
1450       
1451    !! 2.2 Soil texture parameters
1452         
1453    pcent(:) = pcent_usda(:) 
1454    mc_awet(:) = mc_awet_usda(:)
1455    mc_adry(:) = mc_adry_usda(:) 
1456
1457    !! 2.3 Read in the run.def the parameters values defined by the user
1458
1459    !Config Key   = WETNESS_TRANSPIR_MAX
1460    !Config Desc  = Soil moisture above which transpir is max
1461    !Config If    =
1462    !Config Def   = 0.5, 0.5, 0.5
1463    !Config Help  = This parameter is independent from soil texture for
1464    !Config         the time being.
1465    !Config Units = [-]   
1466    CALL getin_p("WETNESS_TRANSPIR_MAX",pcent)
1467
1468    !! Check parameter value (correct range)
1469    IF ( ANY(pcent(:) <= zero) .OR. ANY(pcent(:) > 1.) ) THEN
1470       CALL ipslerr_p(error_level, "hydrol_init.", &
1471            &     "Wrong parameter value for WETNESS_TRANSPIR_MAX.", &
1472            &     "This parameter should be positive and less or equals than 1. ", &
1473            &     "Please, check parameter value in run.def. ")
1474    END IF
1475   
1476
1477    !Config Key   = VWC_MIN_FOR_WET_ALB
1478    !Config Desc  = Vol. wat. cont. above which albedo is cst
1479    !Config If    =
1480    !Config Def   = 0.25, 0.25, 0.25
1481    !Config Help  = This parameter is independent from soil texture for
1482    !Config         the time being.
1483    !Config Units = [m3/m3] 
1484    CALL getin_p("VWC_MIN_FOR_WET_ALB",mc_awet)
1485
1486    !! Check parameter value (correct range)
1487    IF ( ANY(mc_awet(:) < 0) ) THEN
1488       CALL ipslerr_p(error_level, "hydrol_init.", &
1489            &     "Wrong parameter value for VWC_MIN_FOR_WET_ALB.", &
1490            &     "This parameter should be positive. ", &
1491            &     "Please, check parameter value in run.def. ")
1492    END IF
1493
1494
1495    !Config Key   = VWC_MAX_FOR_DRY_ALB
1496    !Config Desc  = Vol. wat. cont. below which albedo is cst
1497    !Config If    =
1498    !Config Def   = 0.1, 0.1, 0.1
1499    !Config Help  = This parameter is independent from soil texture for
1500    !Config         the time being.
1501    !Config Units = [m3/m3]   
1502    CALL getin_p("VWC_MAX_FOR_DRY_ALB",mc_adry)
1503
1504    !! Check parameter value (correct range)
1505    IF ( ANY(mc_adry(:) < 0) .OR. ANY(mc_adry(:) > mc_awet(:)) ) THEN
1506       CALL ipslerr_p(error_level, "hydrol_init.", &
1507            &     "Wrong parameter value for VWC_MAX_FOR_DRY_ALB.", &
1508            &     "This parameter should be positive and not greater than VWC_MIN_FOR_WET_ALB.", &
1509            &     "Please, check parameter value in run.def. ")
1510    END IF
1511
1512
1513    !! 3 Other array allocation
1514
1515
1516    ALLOCATE (mask_veget(kjpindex,nvm),stat=ier)
1517    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_veget','','')
1518
1519    ALLOCATE (mask_soiltile(kjpindex,nstm),stat=ier)
1520    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_soiltile','','')
1521
1522    ALLOCATE (humrelv(kjpindex,nvm,nstm),stat=ier)
1523    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable humrelv','','')
1524
1525    ALLOCATE (vegstressv(kjpindex,nvm,nstm),stat=ier) 
1526    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegstressv','','')
1527
1528    ALLOCATE (us(kjpindex,nvm,nstm,nslm),stat=ier) 
1529    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable us','','')
1530
1531    ALLOCATE (precisol(kjpindex,nvm),stat=ier) 
1532    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable precisol','','')
1533
1534    ALLOCATE (throughfall(kjpindex,nvm),stat=ier) 
1535    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable throughfall','','')
1536
1537    ALLOCATE (precisol_ns(kjpindex,nstm),stat=ier) 
1538    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable precisol_nc','','')
1539
1540    ALLOCATE (free_drain_coef(kjpindex,nstm),stat=ier) 
1541    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable free_drain_coef','','')
1542
1543    ALLOCATE (zwt_force(kjpindex,nstm),stat=ier) 
1544    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zwt_force','','')
1545
1546    ALLOCATE (frac_bare_ns(kjpindex,nstm),stat=ier) 
1547    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable frac_bare_ns','','')
1548
1549    ALLOCATE (water2infilt(kjpindex,nstm),stat=ier)
1550    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable water2infilt','','')
1551
1552    ALLOCATE (ae_ns(kjpindex,nstm),stat=ier) 
1553    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ae_ns','','')
1554
1555    ALLOCATE (rootsink(kjpindex,nslm,nstm),stat=ier) 
1556    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable rootsink','','')
1557
1558    ALLOCATE (subsnowveg(kjpindex),stat=ier) 
1559    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsnowveg','','')
1560
1561    ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier) 
1562    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsnownobio','','')
1563
1564    ALLOCATE (icemelt(kjpindex),stat=ier) 
1565    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable icemelt','','')
1566
1567    ALLOCATE (subsinksoil(kjpindex),stat=ier) 
1568    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsinksoil','','')
1569
1570    ALLOCATE (mx_eau_var(kjpindex),stat=ier)
1571    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mx_eau_var','','')
1572
1573    ALLOCATE (vegtot(kjpindex),stat=ier) 
1574    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegtot','','')
1575
1576    ALLOCATE (vegtot_old(kjpindex),stat=ier) 
1577    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegtot_old','','')
1578
1579    ALLOCATE (resdist(kjpindex,nstm),stat=ier)
1580    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable resdist','','')
1581
1582    ALLOCATE (humtot(kjpindex),stat=ier)
1583    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable humtot','','')
1584
1585    ALLOCATE (resolv(kjpindex),stat=ier) 
1586    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable resolv','','')
1587
1588    ALLOCATE (k(kjpindex,nslm),stat=ier) 
1589    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k','','')
1590
1591    ALLOCATE (kk_moy(kjpindex,nslm),stat=ier) 
1592    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kk_moy','','')
1593    kk_moy(:,:) = 276.48
1594   
1595    ALLOCATE (kk(kjpindex,nslm,nstm),stat=ier) 
1596    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kk','','')
1597    kk(:,:,:) = 276.48
1598   
1599    ALLOCATE (avan_mod_tab(nslm,kjpindex),stat=ier) 
1600    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan_mod_tab','','')
1601   
1602    ALLOCATE (nvan_mod_tab(nslm,kjpindex),stat=ier) 
1603    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan_mod_tab','','')
1604
1605    ALLOCATE (a(kjpindex,nslm),stat=ier) 
1606    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a','','')
1607
1608    ALLOCATE (b(kjpindex,nslm),stat=ier)
1609    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b','','')
1610
1611    ALLOCATE (d(kjpindex,nslm),stat=ier)
1612    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d','','')
1613
1614    ALLOCATE (e(kjpindex,nslm),stat=ier) 
1615    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable e','','')
1616
1617    ALLOCATE (f(kjpindex,nslm),stat=ier) 
1618    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable f','','')
1619
1620    ALLOCATE (g1(kjpindex,nslm),stat=ier) 
1621    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable g1','','')
1622
1623    ALLOCATE (ep(kjpindex,nslm),stat=ier)
1624    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ep','','')
1625
1626    ALLOCATE (fp(kjpindex,nslm),stat=ier)
1627    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable fp','','')
1628
1629    ALLOCATE (gp(kjpindex,nslm),stat=ier)
1630    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable gp','','')
1631
1632    ALLOCATE (rhs(kjpindex,nslm),stat=ier)
1633    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable rhs','','')
1634
1635    ALLOCATE (srhs(kjpindex,nslm),stat=ier)
1636    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable srhs','','')
1637
1638    ALLOCATE (tmc(kjpindex,nstm),stat=ier)
1639    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc','','')
1640
1641    ALLOCATE (tmcs(kjpindex,nstm),stat=ier)
1642    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcs','','')
1643
1644    ALLOCATE (tmcr(kjpindex,nstm),stat=ier)
1645    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcr','','')
1646
1647    ALLOCATE (tmcfc(kjpindex,nstm),stat=ier)
1648    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcfc','','')
1649
1650    ALLOCATE (tmcw(kjpindex,nstm),stat=ier)
1651    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcw','','')
1652
1653    ALLOCATE (tmc_litter(kjpindex,nstm),stat=ier)
1654    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter','','')
1655
1656    ALLOCATE (tmc_litt_mea(kjpindex),stat=ier)
1657    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_mea','','')
1658
1659    ALLOCATE (tmc_litter_res(kjpindex,nstm),stat=ier)
1660    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_res','','')
1661
1662    ALLOCATE (tmc_litter_wilt(kjpindex,nstm),stat=ier)
1663    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_wilt','','')
1664
1665    ALLOCATE (tmc_litter_field(kjpindex,nstm),stat=ier)
1666    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_field','','')
1667
1668    ALLOCATE (tmc_litter_sat(kjpindex,nstm),stat=ier)
1669    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_sat','','')
1670
1671    ALLOCATE (tmc_litter_awet(kjpindex,nstm),stat=ier)
1672    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_awet','','')
1673
1674    ALLOCATE (tmc_litter_adry(kjpindex,nstm),stat=ier)
1675    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_adry','','')
1676
1677    ALLOCATE (tmc_litt_wet_mea(kjpindex),stat=ier)
1678    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_wet_mea','','')
1679
1680    ALLOCATE (tmc_litt_dry_mea(kjpindex),stat=ier)
1681    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_dry_mea','','')
1682
1683    ALLOCATE (v1(kjpindex,nstm),stat=ier)
1684    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable v1','','')
1685
1686    ALLOCATE (ru_ns(kjpindex,nstm),stat=ier)
1687    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ru_ns','','')
1688    ru_ns(:,:) = zero
1689
1690    ALLOCATE (dr_ns(kjpindex,nstm),stat=ier)
1691    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dr_ns','','')
1692    dr_ns(:,:) = zero
1693
1694    ALLOCATE (tr_ns(kjpindex,nstm),stat=ier)
1695    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tr_ns','','')
1696
1697    ALLOCATE (vegetmax_soil(kjpindex,nvm,nstm),stat=ier)
1698    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegetmax_soil','','')
1699
1700    ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier)
1701    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc','','')
1702
1703
1704    ! Variables for nudging of soil moisture
1705    IF (ok_nudge_mc) THEN
1706       ALLOCATE (mc_read_prev(kjpindex,nslm,nstm),stat=ier)
1707       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_read_prev','','')
1708       ALLOCATE (mc_read_next(kjpindex,nslm,nstm),stat=ier)
1709       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_read_next','','')
1710       ALLOCATE (mc_read_current(kjpindex,nslm,nstm),stat=ier)
1711       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_read_current','','')
1712       ALLOCATE (mask_mc_interp(kjpindex,nslm,nstm),stat=ier)
1713       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_mc_interp','','')
1714       ALLOCATE (tmc_aux(kjpindex,nstm),stat=ier)
1715       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_aux','','')
1716    END IF
1717
1718    ! Variables for nudging of snow variables
1719    IF (ok_nudge_snow) THEN
1720       ALLOCATE (snowdz_read_prev(kjpindex,nsnow),stat=ier)
1721       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowdz_read_prev','','')
1722       ALLOCATE (snowdz_read_next(kjpindex,nsnow),stat=ier)
1723       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowdz_read_next','','')
1724       
1725       ALLOCATE (snowrho_read_prev(kjpindex,nsnow),stat=ier)
1726       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowrho_read_prev','','')
1727       ALLOCATE (snowrho_read_next(kjpindex,nsnow),stat=ier)
1728       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowrho_read_next','','')
1729       
1730       ALLOCATE (snowtemp_read_prev(kjpindex,nsnow),stat=ier)
1731       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowtemp_read_prev','','')
1732       ALLOCATE (snowtemp_read_next(kjpindex,nsnow),stat=ier)
1733       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowtemp_read_next','','')
1734       
1735       ALLOCATE (mask_snow_interp(kjpindex,nsnow),stat=ier)
1736       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_snow_interp','','')
1737    END IF
1738
1739    ALLOCATE (mcl(kjpindex, nslm, nstm),stat=ier)
1740    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcl','','')
1741
1742    IF (ok_freeze_cwrr) THEN
1743       ALLOCATE (profil_froz_hydro(kjpindex, nslm),stat=ier)
1744       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable profil_froz_hydrol','','')
1745       profil_froz_hydro(:,:) = zero
1746    ENDIF
1747   
1748    ALLOCATE (profil_froz_hydro_ns(kjpindex, nslm, nstm),stat=ier)
1749    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable profil_froz_hydro_ns','','')
1750    profil_froz_hydro_ns(:,:,:) = zero
1751   
1752    ALLOCATE (soilmoist(kjpindex,nslm),stat=ier)
1753    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soilmoist','','')
1754
1755    ALLOCATE (soilmoist_liquid(kjpindex,nslm),stat=ier)
1756    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soilmoist_liquid','','')
1757
1758    ALLOCATE (soil_wet_ns(kjpindex,nslm,nstm),stat=ier)
1759    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_ns','','')
1760
1761    ALLOCATE (soil_wet_litter(kjpindex,nstm),stat=ier)
1762    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_litter','','')
1763
1764    ALLOCATE (qflux_ns(kjpindex,nslm,nstm),stat=ier) 
1765    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable qflux_ns','','')
1766
1767    ALLOCATE (check_top_ns(kjpindex,nstm),stat=ier) 
1768    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable check_top_ns','','')
1769
1770    ALLOCATE (tmat(kjpindex,nslm,3),stat=ier)
1771    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmat','','')
1772
1773    ALLOCATE (stmat(kjpindex,nslm,3),stat=ier)
1774    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable stmat','','')
1775
1776    ALLOCATE (nroot(kjpindex,nvm, nslm),stat=ier)
1777    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nroot','','')
1778
1779    ALLOCATE (kfact_root(kjpindex, nslm, nstm), stat=ier)
1780    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','')
1781
1782    ALLOCATE (kfact(nslm, kjpindex),stat=ier)
1783    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','')
1784
1785    ALLOCATE (zz(nslm),stat=ier)
1786    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zz','','')
1787
1788    ALLOCATE (dz(nslm),stat=ier)
1789    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dz','','')
1790   
1791    ALLOCATE (dh(nslm),stat=ier)
1792    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','')
1793
1794    ALLOCATE (mc_lin(imin:imax, kjpindex),stat=ier)
1795    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','')
1796
1797    ALLOCATE (k_lin(imin:imax, nslm, kjpindex),stat=ier)
1798    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','')
1799
1800    ALLOCATE (d_lin(imin:imax, nslm, kjpindex),stat=ier)
1801    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','')
1802
1803    ALLOCATE (a_lin(imin:imax, nslm, kjpindex),stat=ier)
1804    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','')
1805
1806    ALLOCATE (b_lin(imin:imax, nslm, kjpindex),stat=ier)
1807    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','')
1808
1809    ALLOCATE (undermcr(kjpindex),stat=ier)
1810    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable undermcr','','')
1811
1812    ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
1813    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watveg_beg','','')
1814   
1815    ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
1816    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watvag_end','','')
1817   
1818    ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
1819    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watsoil_beg','','')
1820   
1821    ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
1822    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watsoil_end','','')
1823   
1824    ALLOCATE (delsoilmoist(kjpindex),stat=ier)
1825    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delsoilmoist','','')
1826   
1827    ALLOCATE (delintercept(kjpindex),stat=ier)
1828    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delintercept','','')
1829   
1830    ALLOCATE (delswe(kjpindex),stat=ier)
1831    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delswe','','')
1832   
1833    ALLOCATE (snow_beg(kjpindex),stat=ier)
1834    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snow_beg','','')
1835   
1836    ALLOCATE (snow_end(kjpindex),stat=ier)
1837    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snow_end','','')
1838   
1839    !! 4 Open restart input file and read data for HYDROLOGIC process
1840       IF (printlev>=3) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
1841
1842       IF (is_root_prc) CALL ioconf_setatt_p('UNITS', '-')
1843       !
1844       DO jst=1,nstm
1845          ! var_name= "mc_1" ... "mc_3"
1846           WRITE (var_name,"('moistc_',I1)") jst
1847           IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
1848           CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc(:,:,jst), "gather", nbp_glo, index_g)
1849       END DO
1850
1851       IF (ok_nudge_mc) THEN
1852          DO jst=1,nstm
1853             WRITE (var_name,"('mc_read_next_',I1)") jst
1854             IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME','Soil moisture read from nudging file')
1855             CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc_read_next(:,:,jst), &
1856                  "gather", nbp_glo, index_g)
1857          END DO
1858       END IF
1859
1860       IF (ok_nudge_snow) THEN
1861          IF (is_root_prc) THEN
1862             CALL ioconf_setatt_p('UNITS', 'm')
1863             CALL ioconf_setatt_p('LONG_NAME','Snow layer thickness read from nudging file')
1864          ENDIF
1865          CALL restget_p (rest_id, 'snowdz_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowdz_read_next, &
1866               "gather", nbp_glo, index_g)
1867
1868          IF (is_root_prc) THEN
1869             CALL ioconf_setatt_p('UNITS', 'kg/m^3')
1870             CALL ioconf_setatt_p('LONG_NAME','Snow density profile read from nudging file')
1871          ENDIF
1872          CALL restget_p (rest_id, 'snowrho_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowrho_read_next, &
1873               "gather", nbp_glo, index_g)
1874
1875          IF (is_root_prc) THEN
1876             CALL ioconf_setatt_p('UNITS', 'K')
1877             CALL ioconf_setatt_p('LONG_NAME','Snow temperature read from nudging file')
1878          ENDIF
1879          CALL restget_p (rest_id, 'snowtemp_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowtemp_read_next, &
1880               "gather", nbp_glo, index_g)
1881       END IF
1882     
1883       DO jst=1,nstm
1884          ! var_name= "mcl_1" ... "mcl_3"
1885           WRITE (var_name,"('moistcl_',I1)") jst
1886           IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
1887           CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mcl(:,:,jst), "gather", nbp_glo, index_g)
1888       END DO
1889
1890       IF (is_root_prc) CALL ioconf_setatt_p('UNITS', '-')
1891       DO jst=1,nstm
1892          DO jsl=1,nslm
1893             ! var_name= "us_1_01" ... "us_3_11"
1894             WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
1895             IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
1896             CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., us(:,:,jst,jsl), "gather", nbp_glo, index_g)
1897          END DO
1898       END DO
1899       !
1900       var_name= 'free_drain_coef'
1901       IF (is_root_prc) THEN
1902          CALL ioconf_setatt_p('UNITS', '-')
1903          CALL ioconf_setatt_p('LONG_NAME','Coefficient for free drainage at bottom of soil')
1904       ENDIF
1905       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., free_drain_coef, "gather", nbp_glo, index_g)
1906       !
1907       var_name= 'zwt_force'
1908       IF (is_root_prc) THEN
1909          CALL ioconf_setatt_p('UNITS', 'm')
1910          CALL ioconf_setatt_p('LONG_NAME','Prescribed water table depth')
1911       ENDIF
1912       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., zwt_force, "gather", nbp_glo, index_g)
1913       !
1914       var_name= 'water2infilt'
1915       IF (is_root_prc) THEN
1916          CALL ioconf_setatt_p('UNITS', '-')
1917          CALL ioconf_setatt_p('LONG_NAME','Remaining water to be infiltrated on top of the soil')
1918       ENDIF
1919       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., water2infilt, "gather", nbp_glo, index_g)
1920       !
1921       var_name= 'ae_ns'
1922       IF (is_root_prc) THEN
1923          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1924          CALL ioconf_setatt_p('LONG_NAME','Bare soil evap on each soil type')
1925       ENDIF
1926       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., ae_ns, "gather", nbp_glo, index_g)
1927       !
1928       var_name= 'snow'       
1929       IF (is_root_prc) THEN
1930          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1931          CALL ioconf_setatt_p('LONG_NAME','Snow mass')
1932       ENDIF
1933       CALL restget_p (rest_id, var_name, nbp_glo, 1  , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
1934       !
1935       var_name= 'snow_age'
1936       IF (is_root_prc) THEN
1937          CALL ioconf_setatt_p('UNITS', 'd')
1938          CALL ioconf_setatt_p('LONG_NAME','Snow age')
1939       ENDIF
1940       CALL restget_p (rest_id, var_name, nbp_glo, 1  , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
1941       !
1942       var_name= 'snow_nobio'
1943       IF (is_root_prc) THEN
1944          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1945          CALL ioconf_setatt_p('LONG_NAME','Snow on other surface types')
1946       ENDIF
1947       CALL restget_p (rest_id, var_name, nbp_glo, nnobio  , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
1948       !
1949       var_name= 'snow_nobio_age'
1950       IF (is_root_prc) THEN
1951          CALL ioconf_setatt_p('UNITS', 'd')
1952          CALL ioconf_setatt_p('LONG_NAME','Snow age on other surface types')
1953       ENDIF
1954       CALL restget_p (rest_id, var_name, nbp_glo, nnobio  , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
1955       !
1956       var_name= 'qsintveg'
1957       IF (is_root_prc) THEN
1958          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1959          CALL ioconf_setatt_p('LONG_NAME','Intercepted moisture')
1960       ENDIF
1961       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
1962
1963       var_name= 'evap_bare_lim_ns'
1964       IF (is_root_prc) THEN
1965          CALL ioconf_setatt_p('UNITS', '?')
1966          CALL ioconf_setatt_p('LONG_NAME','?')
1967       ENDIF
1968       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., evap_bare_lim_ns, "gather", nbp_glo, index_g)
1969       CALL setvar_p (evap_bare_lim_ns, val_exp, 'NO_KEYWORD', 0.0)
1970
1971       var_name= 'resdist'
1972       IF (is_root_prc) THEN
1973          CALL ioconf_setatt_p('UNITS', '-')
1974          CALL ioconf_setatt_p('LONG_NAME','soiltile values from previous time-step')
1975       ENDIF
1976       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
1977
1978       var_name= 'vegtot_old'
1979       IF (is_root_prc) THEN
1980          CALL ioconf_setatt_p('UNITS', '-')
1981          CALL ioconf_setatt_p('LONG_NAME','vegtot from previous time-step')
1982       ENDIF
1983       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_old, "gather", nbp_glo, index_g)       
1984       
1985       ! Read drysoil_frac. It will be initalized later in hydrol_var_init if the varaible is not find in restart file.
1986       IF (is_root_prc) THEN
1987          CALL ioconf_setatt_p('UNITS', '')
1988          CALL ioconf_setatt_p('LONG_NAME','Function of litter wetness')
1989       ENDIF
1990       CALL restget_p (rest_id, 'drysoil_frac', nbp_glo, 1  , 1, kjit, .TRUE., drysoil_frac, "gather", nbp_glo, index_g)
1991
1992
1993    !! 5 get restart values if none were found in the restart file
1994       !
1995       !Config Key   = HYDROL_MOISTURE_CONTENT
1996       !Config Desc  = Soil moisture on each soil tile and levels
1997       !Config If    =
1998       !Config Def   = 0.3
1999       !Config Help  = The initial value of mc if its value is not found
2000       !Config         in the restart file. This should only be used if the model is
2001       !Config         started without a restart file.
2002       !Config Units = [m3/m3]
2003       !
2004       CALL setvar_p (mc, val_exp, 'HYDROL_MOISTURE_CONTENT', 0.3_r_std)
2005
2006       ! Initialize mcl as mc if it is not found in the restart file
2007       IF ( ALL(mcl(:,:,:)==val_exp) ) THEN
2008          mcl(:,:,:) = mc(:,:,:)
2009       END IF
2010
2011
2012       
2013       !Config Key   = US_INIT
2014       !Config Desc  = US_NVM_NSTM_NSLM
2015       !Config If    =
2016       !Config Def   = 0.0
2017       !Config Help  = The initial value of us (relative moisture) if its value is not found
2018       !Config         in the restart file. This should only be used if the model is
2019       !Config         started without a restart file.
2020       !Config Units = [-]
2021       !
2022       DO jsl=1,nslm
2023          CALL setvar_p (us(:,:,:,jsl), val_exp, 'US_INIT', zero)
2024       ENDDO
2025       !
2026       !Config Key   = ZWT_FORCE
2027       !Config Desc  = Prescribed water depth, dimension nstm
2028       !Config If    =
2029       !Config Def   = undef undef undef
2030       !Config Help  = The initial value of zwt_force if its value is not found
2031       !Config         in the restart file. undef corresponds to a case whith no forced WT.
2032       !Config         This should only be used if the model is started without a restart file.
2033       !Config Units = [m]
2034       
2035       ALLOCATE (zwt_default(nstm),stat=ier)
2036       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zwt_default','','')
2037       zwt_default(:) = undef_sechiba
2038       CALL setvar_p (zwt_force, val_exp, 'ZWT_FORCE', zwt_default )
2039
2040       zforce = .FALSE.
2041       DO jst=1,nstm
2042          IF (zwt_force(1,jst) <= zmaxh) zforce = .TRUE. ! AD16*** check if OK with vertical_soil
2043       ENDDO
2044       !
2045       !Config Key   = FREE_DRAIN_COEF
2046       !Config Desc  = Coefficient for free drainage at bottom, dimension nstm
2047       !Config If    =
2048       !Config Def   = 1.0 1.0 1.0
2049       !Config Help  = The initial value of free drainage coefficient if its value is not found
2050       !Config         in the restart file. This should only be used if the model is
2051       !Config         started without a restart file.
2052       !Config Units = [-]
2053             
2054       ALLOCATE (free_drain_max(nstm),stat=ier)
2055       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable free_drain_max','','')
2056       free_drain_max(:)=1.0
2057
2058       CALL setvar_p (free_drain_coef, val_exp, 'FREE_DRAIN_COEF', free_drain_max)
2059       IF (printlev>=2) WRITE (numout,*) ' hydrol_init => free_drain_coef = ',free_drain_coef(1,:)
2060       DEALLOCATE(free_drain_max)
2061
2062       !
2063       !Config Key   = WATER_TO_INFILT
2064       !Config Desc  = Water to be infiltrated on top of the soil
2065       !Config If    =
2066       !Config Def   = 0.0
2067       !Config Help  = The initial value of free drainage if its value is not found
2068       !Config         in the restart file. This should only be used if the model is
2069       !Config         started without a restart file.
2070       !Config Units = [mm]
2071       !
2072       CALL setvar_p (water2infilt, val_exp, 'WATER_TO_INFILT', zero)
2073       !
2074       !Config Key   = EVAPNU_SOIL
2075       !Config Desc  = Bare soil evap on each soil if not found in restart
2076       !Config If    =
2077       !Config Def   = 0.0
2078       !Config Help  = The initial value of bare soils evap if its value is not found
2079       !Config         in the restart file. This should only be used if the model is
2080       !Config         started without a restart file.
2081       !Config Units = [mm]
2082       !
2083       CALL setvar_p (ae_ns, val_exp, 'EVAPNU_SOIL', zero)
2084       !
2085       !Config Key  = HYDROL_SNOW
2086       !Config Desc  = Initial snow mass if not found in restart
2087       !Config If    = OK_SECHIBA
2088       !Config Def   = 0.0
2089       !Config Help  = The initial value of snow mass if its value is not found
2090       !Config         in the restart file. This should only be used if the model is
2091       !Config         started without a restart file.
2092       !Config Units =
2093       !
2094       CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', zero)
2095       !
2096       !Config Key   = HYDROL_SNOWAGE
2097       !Config Desc  = Initial snow age if not found in restart
2098       !Config If    = OK_SECHIBA
2099       !Config Def   = 0.0
2100       !Config Help  = The initial value of snow age if its value is not found
2101       !Config         in the restart file. This should only be used if the model is
2102       !Config         started without a restart file.
2103       !Config Units = ***
2104       !
2105       CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', zero)
2106       !
2107       !Config Key   = HYDROL_SNOW_NOBIO
2108       !Config Desc  = Initial snow amount on ice, lakes, etc. if not found in restart
2109       !Config If    = OK_SECHIBA
2110       !Config Def   = 0.0
2111       !Config Help  = The initial value of snow if its value is not found
2112       !Config         in the restart file. This should only be used if the model is
2113       !Config         started without a restart file.
2114       !Config Units = [mm]
2115       !
2116       CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero)
2117       !
2118       !Config Key   = HYDROL_SNOW_NOBIO_AGE
2119       !Config Desc  = Initial snow age on ice, lakes, etc. if not found in restart
2120       !Config If    = OK_SECHIBA
2121       !Config Def   = 0.0
2122       !Config Help  = The initial value of snow age if its value is not found
2123       !Config         in the restart file. This should only be used if the model is
2124       !Config         started without a restart file.
2125       !Config Units = ***
2126       !
2127       CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero)
2128       !
2129       !Config Key   = HYDROL_QSV
2130       !Config Desc  = Initial water on canopy if not found in restart
2131       !Config If    = OK_SECHIBA
2132       !Config Def   = 0.0
2133       !Config Help  = The initial value of moisture on canopy if its value
2134       !Config         is not found in the restart file. This should only be used if
2135       !Config         the model is started without a restart file.
2136       !Config Units = [mm]
2137       !
2138       CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', zero)
2139
2140    !! 6 Vegetation array     
2141       !
2142       ! If resdist is not in restart file, initialize with soiltile
2143       IF ( MINVAL(resdist) .EQ.  MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
2144          resdist(:,:) = soiltile(:,:)
2145       ENDIF
2146       
2147       !
2148       !  Remember that it is only frac_nobio + SUM(veget_max(,:)) that is equal to 1. Thus we need vegtot
2149       !
2150       IF ( ALL(vegtot_old(:) == val_exp) ) THEN
2151          ! vegtot_old was not found in restart file
2152          DO ji = 1, kjpindex
2153             vegtot_old(ji) = SUM(veget_max(ji,:))
2154          ENDDO
2155       ENDIF
2156       
2157       ! In the initialization phase, vegtot must take the value from previous time-step.
2158       ! This is because hydrol_main is done before veget_max is updated in the end of the time step.
2159       vegtot(:) = vegtot_old(:)
2160       
2161       !
2162       !
2163       ! compute the masks for veget
2164
2165       mask_veget(:,:) = 0
2166       mask_soiltile(:,:) = 0
2167
2168       DO jst=1,nstm
2169          DO ji = 1, kjpindex
2170             IF(soiltile(ji,jst) .GT. min_sechiba) THEN
2171                mask_soiltile(ji,jst) = 1
2172             ENDIF
2173          END DO
2174       ENDDO
2175         
2176       DO jv = 1, nvm
2177          DO ji = 1, kjpindex
2178             IF(veget_max(ji,jv) .GT. min_sechiba) THEN
2179                mask_veget(ji,jv) = 1
2180             ENDIF
2181          END DO
2182       END DO
2183
2184       humrelv(:,:,:) = SUM(us,dim=4)
2185
2186         
2187       !! 7a. Set vegstress
2188     
2189       var_name= 'vegstress'
2190       IF (is_root_prc) THEN
2191          CALL ioconf_setatt_p('UNITS', '-')
2192          CALL ioconf_setatt_p('LONG_NAME','Vegetation growth moisture stress')
2193       ENDIF
2194       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
2195
2196       vegstressv(:,:,:) = humrelv(:,:,:)
2197       ! Calculate vegstress if it is not found in restart file
2198       IF (ALL(vegstress(:,:)==val_exp)) THEN
2199          DO jv=1,nvm
2200             DO ji=1,kjpindex
2201                vegstress(ji,jv)=vegstress(ji,jv) + vegstressv(ji,jv,pref_soil_veg(jv))
2202             END DO
2203          END DO
2204       END IF
2205       !! 7b. Set humrel   
2206       ! Read humrel from restart file
2207       var_name= 'humrel'
2208       IF (is_root_prc) THEN
2209          CALL ioconf_setatt_p('UNITS', '')
2210          CALL ioconf_setatt_p('LONG_NAME','Relative humidity')
2211       ENDIF
2212       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g)
2213
2214       ! Calculate humrel if it is not found in restart file
2215       IF (ALL(humrel(:,:)==val_exp)) THEN
2216          ! set humrel from humrelv, assuming equi-repartition for the first time step
2217          humrel(:,:) = zero
2218          DO jv=1,nvm
2219             DO ji=1,kjpindex
2220                humrel(ji,jv)=humrel(ji,jv) + humrelv(ji,jv,pref_soil_veg(jv))     
2221             END DO
2222          END DO
2223       END IF
2224
2225       ! Read evap_bare_lim from restart file
2226       var_name= 'evap_bare_lim'
2227       IF (is_root_prc) THEN
2228          CALL ioconf_setatt_p('UNITS', '')
2229          CALL ioconf_setatt_p('LONG_NAME','Limitation factor for bare soil evaporation')
2230       ENDIF
2231       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evap_bare_lim, "gather", nbp_glo, index_g)
2232
2233       ! Calculate evap_bare_lim if it was not found in the restart file.
2234       IF ( ALL(evap_bare_lim(:) == val_exp) ) THEN
2235          DO ji = 1, kjpindex
2236             evap_bare_lim(ji) =  SUM(evap_bare_lim_ns(ji,:)*vegtot(ji)*soiltile(ji,:))
2237          ENDDO
2238       END IF
2239
2240
2241    ! Read from restart file       
2242    ! The variables tot_watsoil_beg, tot_watsoil_beg and snwo_beg will be initialized in the end of
2243    ! hydrol_initialize if they were not found in the restart file.
2244       
2245    var_name= 'tot_watveg_beg'
2246    IF (is_root_prc) THEN
2247       CALL ioconf_setatt_p('UNITS', '?')
2248       CALL ioconf_setatt_p('LONG_NAME','?')
2249    ENDIF
2250    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tot_watveg_beg, "gather", nbp_glo, index_g)
2251   
2252    var_name= 'tot_watsoil_beg'
2253    IF (is_root_prc) THEN
2254       CALL ioconf_setatt_p('UNITS', '?')
2255       CALL ioconf_setatt_p('LONG_NAME','?')
2256    ENDIF
2257    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tot_watsoil_beg, "gather", nbp_glo, index_g)
2258   
2259    var_name= 'snow_beg'
2260    IF (is_root_prc) THEN
2261       CALL ioconf_setatt_p('UNITS', '?')
2262       CALL ioconf_setatt_p('LONG_NAME','?')
2263    ENDIF
2264    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., snow_beg, "gather", nbp_glo, index_g)
2265       
2266 
2267    ! Initialize variables for explictsnow module by reading restart file
2268    CALL explicitsnow_initialize( kjit,     kjpindex, rest_id,    snowrho,   &
2269         snowtemp, snowdz,   snowheat,   snowgrain)
2270
2271
2272    ! Initialize soil moisture for nudging if not found in restart file
2273    IF (ok_nudge_mc) THEN
2274       IF ( ALL(mc_read_next(:,:,:)==val_exp) ) mc_read_next(:,:,:) = mc(:,:,:)
2275    END IF
2276   
2277    ! Initialize snow variables for nudging if not found in restart file
2278    IF (ok_nudge_snow) THEN
2279       IF ( ALL(snowdz_read_next(:,:)==val_exp) ) snowdz_read_next(:,:) = snowdz(:,:)
2280       IF ( ALL(snowrho_read_next(:,:)==val_exp) ) snowrho_read_next(:,:) = snowrho(:,:)
2281       IF ( ALL(snowtemp_read_next(:,:)==val_exp) ) snowtemp_read_next(:,:) = snowtemp(:,:)
2282    END IF
2283   
2284   
2285    IF (printlev>=3) WRITE (numout,*) ' hydrol_init done '
2286   
2287  END SUBROUTINE hydrol_init
2288
2289
2290!! ================================================================================================================================
2291!! SUBROUTINE   : hydrol_clear
2292!!
2293!>\BRIEF        Deallocate arrays
2294!!
2295!_ ================================================================================================================================
2296!_ hydrol_clear
2297
2298  SUBROUTINE hydrol_clear()
2299
2300    ! Allocation for soiltile related parameters
2301   
2302    IF ( ALLOCATED (pcent)) DEALLOCATE (pcent)
2303    IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet)
2304    IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry)
2305    ! Other arrays
2306    IF (ALLOCATED (mask_veget)) DEALLOCATE (mask_veget)
2307    IF (ALLOCATED (mask_soiltile)) DEALLOCATE (mask_soiltile)
2308    IF (ALLOCATED (humrelv)) DEALLOCATE (humrelv)
2309    IF (ALLOCATED (vegstressv)) DEALLOCATE (vegstressv)
2310    IF (ALLOCATED (us)) DEALLOCATE (us)
2311    IF (ALLOCATED  (precisol)) DEALLOCATE (precisol)
2312    IF (ALLOCATED  (throughfall)) DEALLOCATE (throughfall)
2313    IF (ALLOCATED  (precisol_ns)) DEALLOCATE (precisol_ns)
2314    IF (ALLOCATED  (free_drain_coef)) DEALLOCATE (free_drain_coef)
2315    IF (ALLOCATED  (frac_bare_ns)) DEALLOCATE (frac_bare_ns)
2316    IF (ALLOCATED  (water2infilt)) DEALLOCATE (water2infilt)
2317    IF (ALLOCATED  (ae_ns)) DEALLOCATE (ae_ns)
2318    IF (ALLOCATED  (rootsink)) DEALLOCATE (rootsink)
2319    IF (ALLOCATED  (subsnowveg)) DEALLOCATE (subsnowveg)
2320    IF (ALLOCATED  (subsnownobio)) DEALLOCATE (subsnownobio)
2321    IF (ALLOCATED  (icemelt)) DEALLOCATE (icemelt)
2322    IF (ALLOCATED  (subsinksoil)) DEALLOCATE (subsinksoil)
2323    IF (ALLOCATED  (mx_eau_var)) DEALLOCATE (mx_eau_var)
2324    IF (ALLOCATED  (vegtot)) DEALLOCATE (vegtot)
2325    IF (ALLOCATED  (vegtot_old)) DEALLOCATE (vegtot_old)
2326    IF (ALLOCATED  (resdist)) DEALLOCATE (resdist)
2327    IF (ALLOCATED  (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
2328    IF (ALLOCATED  (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
2329    IF (ALLOCATED  (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
2330    IF (ALLOCATED  (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
2331    IF (ALLOCATED  (delsoilmoist)) DEALLOCATE (delsoilmoist)
2332    IF (ALLOCATED  (delintercept)) DEALLOCATE (delintercept)
2333    IF (ALLOCATED  (snow_beg)) DEALLOCATE (snow_beg)
2334    IF (ALLOCATED  (snow_end)) DEALLOCATE (snow_end)
2335    IF (ALLOCATED  (delswe)) DEALLOCATE (delswe)
2336    IF (ALLOCATED  (undermcr)) DEALLOCATE (undermcr)
2337    IF (ALLOCATED  (v1)) DEALLOCATE (v1)
2338    IF (ALLOCATED  (humtot)) DEALLOCATE (humtot)
2339    IF (ALLOCATED  (resolv)) DEALLOCATE (resolv)
2340    IF (ALLOCATED  (k)) DEALLOCATE (k)
2341    IF (ALLOCATED  (kk)) DEALLOCATE (kk)
2342    IF (ALLOCATED  (kk_moy)) DEALLOCATE (kk_moy)
2343    IF (ALLOCATED  (avan_mod_tab)) DEALLOCATE (avan_mod_tab)
2344    IF (ALLOCATED  (nvan_mod_tab)) DEALLOCATE (nvan_mod_tab)
2345    IF (ALLOCATED  (a)) DEALLOCATE (a)
2346    IF (ALLOCATED  (b)) DEALLOCATE (b)
2347    IF (ALLOCATED  (d)) DEALLOCATE (d)
2348    IF (ALLOCATED  (e)) DEALLOCATE (e)
2349    IF (ALLOCATED  (f)) DEALLOCATE (f)
2350    IF (ALLOCATED  (g1)) DEALLOCATE (g1)
2351    IF (ALLOCATED  (ep)) DEALLOCATE (ep)
2352    IF (ALLOCATED  (fp)) DEALLOCATE (fp)
2353    IF (ALLOCATED  (gp)) DEALLOCATE (gp)
2354    IF (ALLOCATED  (rhs)) DEALLOCATE (rhs)
2355    IF (ALLOCATED  (srhs)) DEALLOCATE (srhs)
2356    IF (ALLOCATED  (tmc)) DEALLOCATE (tmc)
2357    IF (ALLOCATED  (tmcs)) DEALLOCATE (tmcs)
2358    IF (ALLOCATED  (tmcr)) DEALLOCATE (tmcr)
2359    IF (ALLOCATED  (tmcfc)) DEALLOCATE (tmcfc)
2360    IF (ALLOCATED  (tmcw)) DEALLOCATE (tmcw)
2361    IF (ALLOCATED  (tmc_litter)) DEALLOCATE (tmc_litter)
2362    IF (ALLOCATED  (tmc_litt_mea)) DEALLOCATE (tmc_litt_mea)
2363    IF (ALLOCATED  (tmc_litter_res)) DEALLOCATE (tmc_litter_res)
2364    IF (ALLOCATED  (tmc_litter_wilt)) DEALLOCATE (tmc_litter_wilt)
2365    IF (ALLOCATED  (tmc_litter_field)) DEALLOCATE (tmc_litter_field)
2366    IF (ALLOCATED  (tmc_litter_sat)) DEALLOCATE (tmc_litter_sat)
2367    IF (ALLOCATED  (tmc_litter_awet)) DEALLOCATE (tmc_litter_awet)
2368    IF (ALLOCATED  (tmc_litter_adry)) DEALLOCATE (tmc_litter_adry)
2369    IF (ALLOCATED  (tmc_litt_wet_mea)) DEALLOCATE (tmc_litt_wet_mea)
2370    IF (ALLOCATED  (tmc_litt_dry_mea)) DEALLOCATE (tmc_litt_dry_mea)
2371    IF (ALLOCATED  (ru_ns)) DEALLOCATE (ru_ns)
2372    IF (ALLOCATED  (dr_ns)) DEALLOCATE (dr_ns)
2373    IF (ALLOCATED  (tr_ns)) DEALLOCATE (tr_ns)
2374    IF (ALLOCATED  (vegetmax_soil)) DEALLOCATE (vegetmax_soil)
2375    IF (ALLOCATED  (mc)) DEALLOCATE (mc)
2376    IF (ALLOCATED  (soilmoist)) DEALLOCATE (soilmoist)
2377    IF (ALLOCATED  (soilmoist_liquid)) DEALLOCATE (soilmoist_liquid)
2378    IF (ALLOCATED  (soil_wet_ns)) DEALLOCATE (soil_wet_ns)
2379    IF (ALLOCATED  (soil_wet_litter)) DEALLOCATE (soil_wet_litter)
2380    IF (ALLOCATED  (qflux_ns)) DEALLOCATE (qflux_ns)
2381    IF (ALLOCATED  (tmat)) DEALLOCATE (tmat)
2382    IF (ALLOCATED  (stmat)) DEALLOCATE (stmat)
2383    IF (ALLOCATED  (nroot)) DEALLOCATE (nroot)
2384    IF (ALLOCATED  (kfact_root)) DEALLOCATE (kfact_root)
2385    IF (ALLOCATED  (kfact)) DEALLOCATE (kfact)
2386    IF (ALLOCATED  (zz)) DEALLOCATE (zz)
2387    IF (ALLOCATED  (dz)) DEALLOCATE (dz)
2388    IF (ALLOCATED  (dh)) DEALLOCATE (dh)
2389    IF (ALLOCATED  (mc_lin)) DEALLOCATE (mc_lin)
2390    IF (ALLOCATED  (k_lin)) DEALLOCATE (k_lin)
2391    IF (ALLOCATED  (d_lin)) DEALLOCATE (d_lin)
2392    IF (ALLOCATED  (a_lin)) DEALLOCATE (a_lin)
2393    IF (ALLOCATED  (b_lin)) DEALLOCATE (b_lin)
2394
2395  END SUBROUTINE hydrol_clear
2396
2397!! ================================================================================================================================
2398!! SUBROUTINE   : hydrol_tmc_update
2399!!
2400!>\BRIEF        This routine updates the soil moisture profiles when the vegetation fraction have changed.
2401!!
2402!! DESCRIPTION  :
2403!!
2404!!    This routine update tmc and mc with variation of veget_max (LAND_USE or DGVM activated)
2405!!
2406!!
2407!!
2408!!
2409!! RECENT CHANGE(S) : Adaptation to excluding nobio from soiltile(1)
2410!!
2411!! MAIN OUTPUT VARIABLE(S) :
2412!!
2413!! REFERENCE(S) :
2414!!
2415!! FLOWCHART    : None
2416!! \n
2417!_ ================================================================================================================================
2418!_ hydrol_tmc_update
2419  SUBROUTINE hydrol_tmc_update ( kjpindex, veget_max, soiltile, qsintveg, drain_upd, runoff_upd)
2420
2421    !! 0.1 Input variables
2422    INTEGER(i_std), INTENT(in)                            :: kjpindex      !! domain size
2423    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)     :: veget_max     !! max fraction of vegetation type
2424    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)   :: soiltile      !! Fraction of each soil tile (0-1, unitless)
2425
2426    !! 0.2 Output variables
2427    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: drain_upd        !! Change in drainage due to decrease in vegtot
2428                                                                              !! on mc [kg/m2/dt]
2429    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: runoff_upd       !! Change in runoff due to decrease in vegtot
2430                                                                              !! on water2infilt[kg/m2/dt]
2431   
2432    !! 0.3 Modified variables
2433    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)  :: qsintveg   !! Amount of water in the canopy interception
2434
2435    !! 0.4 Local variables
2436    INTEGER(i_std)                           :: ji, jv, jst,jsl
2437    LOGICAL                                  :: soil_upd        !! True if soiltile changed since last time step
2438    LOGICAL                                  :: vegtot_upd      !! True if vegtot changed since last time step
2439    REAL(r_std), DIMENSION(kjpindex,nstm)    :: vmr             !! Change in soiltile (within vegtot)
2440    REAL(r_std), DIMENSION(kjpindex)         :: vmr_sum
2441    REAL(r_std), DIMENSION(kjpindex)         :: delvegtot   
2442    REAL(r_std), DIMENSION(kjpindex,nslm)    :: mc_dilu         !! Total loss of moisture content
2443    REAL(r_std), DIMENSION(kjpindex)         :: infil_dilu      !! Total loss for water2infilt
2444    REAL(r_std), DIMENSION(kjpindex,nstm)    :: tmc_old         !! tmc before calculations
2445    REAL(r_std), DIMENSION(kjpindex,nstm)    :: water2infilt_old!! water2infilt before calculations
2446    REAL(r_std), DIMENSION (kjpindex,nvm)    :: qsintveg_old    !! qsintveg before calculations
2447    REAL(r_std), DIMENSION(kjpindex)         :: test
2448    REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: mcaux        !! serves to hold the chnage in mc when vegtot decreases
2449
2450   
2451    !! 1. If a PFT has disapperead as result from a veget_max change,
2452    !!    then add canopy water to surface water.
2453    !     Other adaptations of qsintveg are delt by the normal functioning of hydrol_canop
2454
2455    DO ji=1,kjpindex
2456       IF (vegtot_old(ji) .GT.min_sechiba) THEN
2457          DO jv=1,nvm
2458             IF ((veget_max(ji,jv).LT.min_sechiba).AND.(qsintveg(ji,jv).GT.0.)) THEN
2459                jst=pref_soil_veg(jv) ! soil tile index
2460                water2infilt(ji,jst) = water2infilt(ji,jst) + qsintveg(ji,jv)/(resdist(ji,jst)*vegtot_old(ji))
2461                qsintveg(ji,jv) = zero
2462             ENDIF
2463          ENDDO
2464       ENDIF
2465    ENDDO
2466   
2467    !! 2. We now deal with the changes of soiltile and corresponding soil moistures
2468    !!    Because sum(soiltile)=1 whatever vegtot, we need to distinguish two cases:
2469    !!    - when vegtot changes (meaning that the nobio fraction changes too),
2470    !!    - and when vegtot does not changes (a priori the most frequent case)
2471
2472    vegtot_upd = SUM(ABS((vegtot(:)-vegtot_old(:)))) .GT. zero ! True if at least one land point with a vegtot change
2473    runoff_upd(:) = zero
2474    drain_upd(:) = zero
2475    IF (vegtot_upd) THEN
2476       ! We find here the processing specific to the chnages of nobio fraction and vegtot
2477
2478       delvegtot(:) = vegtot(:) - vegtot_old(:)
2479
2480       DO jst=1,nstm
2481          DO ji=1,kjpindex
2482
2483             IF (delvegtot(ji) .GT. min_sechiba) THEN
2484
2485                !! 2.1. If vegtot increases (nobio decreases), then the mc in each soiltile is decreased
2486                !!      assuming the same proportions for each soiltile, and each soil layer
2487               
2488                mc(ji,:,jst) = mc(ji,:,jst) * vegtot_old(ji)/vegtot(ji) ! vegtot cannot be zero as > vegtot_old
2489                water2infilt(ji,jst) = water2infilt(ji,jst) * vegtot_old(ji)/vegtot(ji)
2490
2491             ELSE
2492
2493                !! 2.2 If vegtot decreases (nobio increases), then the mc in each soiltile should increase,
2494                !!     but should not exceed mcs
2495                !!     For simplicity, we choose to send the corresponding water volume to drainage
2496                !!     We do the same for water2infilt but send the excess to surface runoff
2497
2498                IF (vegtot(ji) .GT.min_sechiba) THEN
2499                   mcaux(ji,:,jst) =  mc(ji,:,jst) * (vegtot_old(ji)-vegtot(ji))/vegtot(ji) ! mcaux is the delta mc
2500                ELSE ! we just have nobio in the grid-cell
2501                   mcaux(ji,:,jst) =  mc(ji,:,jst)
2502                ENDIF
2503               
2504                drain_upd(ji) = drain_upd(ji) + dz(2) * ( trois*mcaux(ji,1,jst) + mcaux(ji,2,jst) )/huit
2505                DO jsl = 2,nslm-1
2506                   drain_upd(ji) = drain_upd(ji) + dz(jsl) * (trois*mcaux(ji,jsl,jst)+mcaux(ji,jsl-1,jst))/huit &
2507                        + dz(jsl+1) * (trois*mcaux(ji,jsl,jst)+mcaux(ji,jsl+1,jst))/huit
2508                ENDDO
2509                drain_upd(ji) = drain_upd(ji) + dz(nslm) * (trois*mcaux(ji,nslm,jst) + mcaux(ji,nslm-1,jst))/huit
2510
2511                IF (vegtot(ji) .GT.min_sechiba) THEN
2512                   runoff_upd(ji) = runoff_upd(ji) + water2infilt(ji,jst) * (vegtot_old(ji)-vegtot(ji))/vegtot(ji)
2513                ELSE ! we just have nobio in the grid-cell
2514                   runoff_upd(ji) = runoff_upd(ji) + water2infilt(ji,jst)
2515                ENDIF
2516
2517             ENDIF
2518             
2519          ENDDO
2520       ENDDO
2521       
2522    ENDIF
2523   
2524    !! 3. At the end of step 2, we are back to a case where vegtot changes are treated, so we can use soiltile
2525    !!    as a fraction of vegtot to process the mc transfers between soil tiles due to the changes of vegetation map
2526   
2527    !! 3.1 Check if soiltiles changed since last time step
2528    soil_upd=SUM(ABS(soiltile(:,:)-resdist(:,:))) .GT. zero
2529    IF (printlev>=3) WRITE (numout,*) 'soil_upd ', soil_upd
2530       
2531    IF (soil_upd) THEN
2532     
2533       !! 3.2 Define the change in soiltile
2534       vmr(:,:) = soiltile(:,:) - resdist(:,:)  ! resdist is the previous values of soiltiles, previous timestep, so before new map
2535
2536       ! Total area loss by the three soil tiles
2537       DO ji=1,kjpindex
2538          vmr_sum(ji)=SUM(vmr(ji,:),MASK=vmr(ji,:).LT.zero)
2539       ENDDO
2540
2541       !! 3.3 Shrinking soil tiles
2542       !! 3.3.1 Total loss of moisture content from the shrinking soil tiles, expressed by soil layer
2543       mc_dilu(:,:)=zero
2544       DO jst=1,nstm
2545          DO jsl = 1, nslm
2546             DO ji=1,kjpindex
2547                IF ( vmr(ji,jst) < -min_sechiba ) THEN
2548                   mc_dilu(ji,jsl) = mc_dilu(ji,jsl) + mc(ji,jsl,jst) * vmr(ji,jst) / vmr_sum(ji)
2549                ENDIF
2550             ENDDO
2551          ENDDO
2552       ENDDO
2553
2554       !! 3.3.2 Total loss of water2inft from the shrinking soil tiles
2555       infil_dilu(:)=zero
2556       DO jst=1,nstm
2557          DO ji=1,kjpindex
2558             IF ( vmr(ji,jst) < -min_sechiba ) THEN
2559                infil_dilu(ji) = infil_dilu(ji) + water2infilt(ji,jst) * vmr(ji,jst) / vmr_sum(ji)
2560             ENDIF
2561          ENDDO
2562       ENDDO
2563
2564       !! 3.4 Each gaining soil tile gets moisture proportionally to both the total loss and its areal increase
2565
2566       ! As the original mc from each soil tile are in [mcr,mcs] and we do weighted avrage, the new mc are in [mcr,mcs]
2567       ! The case where the soiltile is created (soiltile_old=0) works as the other cases
2568
2569       ! 3.4.1 Update mc(kjpindex,nslm,nstm) !m3/m3
2570       DO jst=1,nstm
2571          DO jsl = 1, nslm
2572             DO ji=1,kjpindex
2573                IF ( vmr(ji,jst) > min_sechiba ) THEN
2574                   mc(ji,jsl,jst) = ( mc(ji,jsl,jst) * resdist(ji,jst) + mc_dilu(ji,jsl) * vmr(ji,jst) ) / soiltile(ji,jst)
2575                   ! NB : soiltile can not be zero for case vmr > zero, see slowproc_veget
2576                ENDIF
2577             ENDDO
2578          ENDDO
2579       ENDDO
2580       
2581       ! 3.4.2 Update water2inft
2582       DO jst=1,nstm
2583          DO ji=1,kjpindex
2584             IF ( vmr(ji,jst) > min_sechiba ) THEN !donc soiltile>0     
2585                water2infilt(ji,jst) = ( water2infilt(ji,jst) * resdist(ji,jst) + infil_dilu(ji) * vmr(ji,jst) ) / soiltile(ji,jst)
2586             ENDIF !donc resdist>0
2587          ENDDO
2588       ENDDO
2589
2590       ! 3.4.3 Case where soiltile < min_sechiba
2591       DO jst=1,nstm
2592          DO ji=1,kjpindex
2593             IF ( soiltile(ji,jst) .LT. min_sechiba ) THEN
2594                water2infilt(ji,jst) = zero
2595                mc(ji,:,jst) = zero
2596             ENDIF
2597          ENDDO
2598       ENDDO
2599
2600    ENDIF ! soil_upd
2601
2602    !! 4. Update tmc and humtot
2603   
2604    DO jst=1,nstm
2605       DO ji=1,kjpindex
2606             tmc(ji,jst) = dz(2) * ( trois*mc(ji,1,jst) + mc(ji,2,jst) )/huit
2607             DO jsl = 2,nslm-1
2608                tmc(ji,jst) = tmc(ji,jst) + dz(jsl) * (trois*mc(ji,jsl,jst)+mc(ji,jsl-1,jst))/huit &
2609                     + dz(jsl+1) * (trois*mc(ji,jsl,jst)+mc(ji,jsl+1,jst))/huit
2610             ENDDO
2611             tmc(ji,jst) = tmc(ji,jst) + dz(nslm) * (trois*mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
2612             tmc(ji,jst) = tmc(ji,jst) + water2infilt(ji,jst)
2613             ! WARNING tmc is increased by includes water2infilt(ji,jst)
2614       ENDDO
2615    ENDDO
2616
2617    humtot(:) = zero
2618    DO jst=1,nstm
2619       DO ji=1,kjpindex
2620          humtot(ji) = humtot(ji) + vegtot(ji) * soiltile(ji,jst) * tmc(ji,jst) ! average over grid-cell (i.e. total land)
2621       ENDDO
2622    ENDDO
2623
2624
2625    !! Now that the work is done, update resdist
2626    resdist(:,:) = soiltile(:,:)
2627
2628    IF (printlev>=3) WRITE (numout,*) ' hydrol_tmc_update done '
2629
2630  END SUBROUTINE hydrol_tmc_update
2631
2632!! ================================================================================================================================
2633!! SUBROUTINE   : hydrol_var_init
2634!!
2635!>\BRIEF        This routine initializes hydrologic parameters to define K and D, and diagnostic hydrologic variables. 
2636!!
2637!! DESCRIPTION  :
2638!! - 1 compute the depths
2639!! - 2 compute the profile for roots
2640!! - 3 compute the profile for a and n Van Genuchten parameter
2641!! - 4 compute the linearized values of k, a, b and d for the resolution of Fokker Planck equation
2642!! - 5 water reservoirs initialisation
2643!!
2644!! RECENT CHANGE(S) : None
2645!!
2646!! MAIN OUTPUT VARIABLE(S) :
2647!!
2648!! REFERENCE(S) :
2649!!
2650!! FLOWCHART    : None
2651!! \n
2652!_ ================================================================================================================================
2653!_ hydrol_var_init
2654
2655  SUBROUTINE hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
2656       kjpindex, veget, veget_max, soiltile, njsc, &
2657       mx_eau_var, shumdiag_perma, &
2658       drysoil_frac, qsintveg, mc_layh, mcl_layh) 
2659
2660    ! interface description
2661
2662    !! 0. Variable and parameter declaration
2663
2664    !! 0.1 Input variables
2665    ! input scalar
2666    INTEGER(i_std), INTENT(in)                          :: kjpindex      !! Domain size (number of grid cells) (1)
2667    ! input fields
2668    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max     !! PFT fractions within grid-cells (1; 1)
2669    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget         !! Effective fraction of vegetation by PFT (1; 1)
2670    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc          !! Index of the dominant soil textural class
2671                                                                         !! in the grid cell (1-nscm, unitless)
2672    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltile      !! Fraction of each soil tile within vegtot (0-1, unitless)
2673   
2674    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
2675    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
2676    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
2677    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
2678    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
2679    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
2680    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
2681 
2682    !! 0.2 Output variables
2683
2684    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: mx_eau_var    !! Maximum water content of the soil
2685                                                                         !! @tex $(kg m^{-2})$ @endtex
2686    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: shumdiag_perma!! Percent of porosity filled with water (mc/mcs)
2687                                                                         !! used for the thermal computations
2688    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)    :: drysoil_frac  !! function of litter humidity
2689    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out):: mc_layh       !! Volumetric soil moisture content for each layer in hydrol(liquid+ice) [m3/m3]
2690    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out):: mcl_layh      !! Volumetric soil moisture content for each layer in hydrol(liquid) [m3/m3]
2691
2692    !! 0.3 Modified variables
2693    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)  :: qsintveg    !! Water on vegetation due to interception
2694                                                                         !! @tex $(kg m^{-2})$ @endtex 
2695
2696    !! 0.4 Local variables
2697
2698    INTEGER(i_std)                                      :: ji, jv        !! Grid-cell and PFT indices (1)
2699    INTEGER(i_std)                                      :: jst, jsc, jsl !! Soiltile, Soil Texture, and Soil layer indices (1)
2700    INTEGER(i_std)                                      :: i             !! Index (1)
2701    REAL(r_std)                                         :: m             !! m=1-1/n (unitless)
2702    REAL(r_std)                                         :: frac          !! Relative linearized VWC (unitless)
2703    REAL(r_std)                                         :: avan_mod      !! VG parameter a modified from  exponantial profile
2704                                                                         !! @tex $(mm^{-1})$ @endtex
2705    REAL(r_std)                                         :: nvan_mod      !! VG parameter n  modified from  exponantial profile
2706                                                                         !! (unitless)
2707    REAL(r_std), DIMENSION(nslm,kjpindex)               :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
2708                                                                         !! (unitless)
2709    ! parameters for "soil densification" with depth
2710    REAL(r_std)                                         :: dp_comp       !! Depth at which the 'compacted' value of ksat
2711                                                                         !! is reached (m)
2712    REAL(r_std)                                         :: f_ks          !! Exponential factor for decay of ksat with depth
2713                                                                         !! @tex $(m^{-1})$ @endtex
2714    ! Fixed parameters from fitted relationships
2715    REAL(r_std)                                         :: n0            !! fitted value for relation log((n-n0)/(n_ref-n0)) =
2716                                                                         !! nk_rel * log(k/k_ref)
2717                                                                         !! (unitless)
2718    REAL(r_std)                                         :: nk_rel        !! fitted value for relation log((n-n0)/(n_ref-n0)) =
2719                                                                         !! nk_rel * log(k/k_ref)
2720                                                                         !! (unitless)
2721    REAL(r_std)                                         :: a0            !! fitted value for relation log((a-a0)/(a_ref-a0)) =
2722                                                                         !! ak_rel * log(k/k_ref)
2723                                                                         !! @tex $(mm^{-1})$ @endtex
2724    REAL(r_std)                                         :: ak_rel        !! fitted value for relation log((a-a0)/(a_ref-a0)) =
2725                                                                         !! ak_rel * log(k/k_ref)
2726                                                                         !! (unitless)
2727    REAL(r_std)                                         :: kfact_max     !! Maximum factor for Ks decay with depth (unitless)
2728    REAL(r_std)                                         :: k_tmp, tmc_litter_ratio
2729    INTEGER(i_std), PARAMETER                           :: error_level = 3 !! Error level for consistency check
2730                                                                           !! Switch to 2 tu turn fatal errors into warnings
2731    REAL(r_std), DIMENSION (kjpindex,nslm)              :: alphavg         !! VG param a modified with depth at each node
2732                                                                           !! @tex $(mm^{-1})$ @endtexe
2733    REAL(r_std), DIMENSION (kjpindex,nslm)              :: nvg             !! VG param n modified with depth at each node
2734                                                                           !! (unitless)
2735                                                                           !! need special treatment
2736    INTEGER(i_std)                                      :: ii
2737    INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin
2738                                                                           !! need special treatment
2739
2740!_ ================================================================================================================================
2741
2742    !Config Key   = CWRR_NKS_N0
2743    !Config Desc  = fitted value for relation log((n-n0)/(n_ref-n0)) = nk_rel * log(k/k_ref)
2744    !Config Def   = 0.0
2745    !Config If    =
2746    !Config Help  =
2747    !Config Units = [-]
2748    n0 = 0.0
2749    CALL getin_p("CWRR_NKS_N0",n0)
2750
2751    !! Check parameter value (correct range)
2752    IF ( n0 < zero ) THEN
2753       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2754            &     "Wrong parameter value for CWRR_NKS_N0.", &
2755            &     "This parameter should be non-negative. ", &
2756            &     "Please, check parameter value in run.def. ")
2757    END IF
2758
2759
2760    !Config Key   = CWRR_NKS_POWER
2761    !Config Desc  = fitted value for relation log((n-n0)/(n_ref-n0)) = nk_rel * log(k/k_ref)
2762    !Config Def   = 0.0
2763    !Config If    =
2764    !Config Help  =
2765    !Config Units = [-]
2766    nk_rel = 0.0
2767    CALL getin_p("CWRR_NKS_POWER",nk_rel)
2768
2769    !! Check parameter value (correct range)
2770    IF ( nk_rel < zero ) THEN
2771       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2772            &     "Wrong parameter value for CWRR_NKS_POWER.", &
2773            &     "This parameter should be non-negative. ", &
2774            &     "Please, check parameter value in run.def. ")
2775    END IF
2776
2777
2778    !Config Key   = CWRR_AKS_A0
2779    !Config Desc  = fitted value for relation log((a-a0)/(a_ref-a0)) = ak_rel * log(k/k_ref)
2780    !Config Def   = 0.0
2781    !Config If    =
2782    !Config Help  =
2783    !Config Units = [1/mm]
2784    a0 = 0.0
2785    CALL getin_p("CWRR_AKS_A0",a0)
2786
2787    !! Check parameter value (correct range)
2788    IF ( a0 < zero ) THEN
2789       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2790            &     "Wrong parameter value for CWRR_AKS_A0.", &
2791            &     "This parameter should be non-negative. ", &
2792            &     "Please, check parameter value in run.def. ")
2793    END IF
2794
2795
2796    !Config Key   = CWRR_AKS_POWER
2797    !Config Desc  = fitted value for relation log((a-a0)/(a_ref-a0)) = ak_rel * log(k/k_ref)
2798    !Config Def   = 0.0
2799    !Config If    =
2800    !Config Help  =
2801    !Config Units = [-]
2802    ak_rel = 0.0
2803    CALL getin_p("CWRR_AKS_POWER",ak_rel)
2804
2805    !! Check parameter value (correct range)
2806    IF ( nk_rel < zero ) THEN
2807       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2808            &     "Wrong parameter value for CWRR_AKS_POWER.", &
2809            &     "This parameter should be non-negative. ", &
2810            &     "Please, check parameter value in run.def. ")
2811    END IF
2812
2813
2814    !Config Key   = KFACT_DECAY_RATE
2815    !Config Desc  = Factor for Ks decay with depth
2816    !Config Def   = 2.0
2817    !Config If    =
2818    !Config Help  = 
2819    !Config Units = [1/m]
2820    f_ks = 2.0
2821    CALL getin_p ("KFACT_DECAY_RATE", f_ks)
2822
2823    !! Check parameter value (correct range)
2824    IF ( f_ks < zero ) THEN
2825       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2826            &     "Wrong parameter value for KFACT_DECAY_RATE.", &
2827            &     "This parameter should be positive. ", &
2828            &     "Please, check parameter value in run.def. ")
2829    END IF
2830
2831
2832    !Config Key   = KFACT_STARTING_DEPTH
2833    !Config Desc  = Depth for compacted value of Ks
2834    !Config Def   = 0.3
2835    !Config If    =
2836    !Config Help  = 
2837    !Config Units = [m]
2838    dp_comp = 0.3
2839    CALL getin_p ("KFACT_STARTING_DEPTH", dp_comp)
2840
2841    !! Check parameter value (correct range)
2842    IF ( dp_comp <= zero ) THEN
2843       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2844            &     "Wrong parameter value for KFACT_STARTING_DEPTH.", &
2845            &     "This parameter should be positive. ", &
2846            &     "Please, check parameter value in run.def. ")
2847    END IF
2848
2849
2850    !Config Key   = KFACT_MAX
2851    !Config Desc  = Maximum Factor for Ks increase due to vegetation
2852    !Config Def   = 10.0
2853    !Config If    =
2854    !Config Help  =
2855    !Config Units = [-]
2856    kfact_max = 10.0
2857    CALL getin_p ("KFACT_MAX", kfact_max)
2858
2859    !! Check parameter value (correct range)
2860    IF ( kfact_max < 10. ) THEN
2861       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2862            &     "Wrong parameter value for KFACT_MAX.", &
2863            &     "This parameter should be greater than 10. ", &
2864            &     "Please, check parameter value in run.def. ")
2865    END IF
2866
2867   
2868    !-
2869    !! 1 Create local variables in mm for the vertical depths
2870    !!   Vertical depth variables (znh, dnh, dlh) are stored in module vertical_soil_var in m.
2871    DO jsl=1,nslm
2872       zz(jsl) = znh(jsl)*mille
2873       dz(jsl) = dnh(jsl)*mille
2874       dh(jsl) = dlh(jsl)*mille
2875    ENDDO
2876
2877    !-
2878    !! 2 Compute the root density profile if not ok_dynroot
2879    !!   For the case with ok_dynroot, the calculations are done at each time step in hydrol_soil
2880    IF (.NOT. ok_dynroot) THEN
2881       DO ji=1, kjpindex
2882          !-
2883          !! The three following equations concerning nroot computation are derived from the integrals
2884          !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158).
2885          !! The occasional absence of minus sign before humcste parameter is correct.
2886          DO jv = 1,nvm
2887             DO jsl = 2, nslm-1
2888                nroot(ji,jv,jsl) = (EXP(-humcste(jv)*zz(jsl)/mille)) * &
2889                     & (EXP(humcste(jv)*dz(jsl)/mille/deux) - &
2890                     & EXP(-humcste(jv)*dz(jsl+1)/mille/deux))/ &
2891                     & (EXP(-humcste(jv)*dz(2)/mille/deux) &
2892                     & -EXP(-humcste(jv)*zz(nslm)/mille))
2893             ENDDO
2894             nroot(ji,jv,1) = zero
2895
2896             nroot(ji,jv,nslm) = (EXP(humcste(jv)*dz(nslm)/mille/deux) -un) * &
2897                  & EXP(-humcste(jv)*zz(nslm)/mille) / &
2898                  & (EXP(-humcste(jv)*dz(2)/mille/deux) &
2899                  & -EXP(-humcste(jv)*zz(nslm)/mille))
2900          ENDDO
2901       ENDDO
2902    END IF
2903
2904 
2905
2906    !-
2907    !! 3 Compute the profile for a and n
2908    !-
2909    DO ji = 1, kjpindex
2910       DO jsl=1,nslm
2911          ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2
2912          ! Calibrated against Hapex-Sahel measurements
2913          kfact(jsl,ji) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
2914          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14
2915
2916          nfact(jsl,ji) = ( kfact(jsl,ji) )**nk_rel
2917          afact(jsl,ji) = ( kfact(jsl,ji) )**ak_rel
2918       ENDDO
2919    ENDDO
2920   
2921    ! For every grid cell
2922     DO ji = 1, kjpindex
2923       !-
2924       !! 4 Compute the linearized values of k, a, b and d
2925       !!   The effect of kfact_root on ks thus on k, a, n and d, is taken into account further in the code,
2926       !!   in hydrol_soil_coef.
2927       !-
2928       ! Calculate the matrix coef for Dublin model (de Rosnay, 1999; p149)
2929       ! piece-wise linearised hydraulic conductivity k_lin=alin * mc_lin + b_lin
2930       ! and diffusivity d_lin in each interval of mc, called mc_lin,
2931       ! between imin, for residual mcr, and imax for saturation mcs.
2932
2933       ! We define 51 bounds for 50 bins of mc between mcr and mcs
2934       mc_lin(imin,ji)=mcr(ji)
2935       mc_lin(imax,ji)=mcs(ji)
2936       DO ii= imin+1, imax-1 ! ii=2,50
2937          mc_lin(ii,ji) = mcr(ji) + (ii-imin)*(mcs(ji)-mcr(ji))/(imax-imin)
2938       ENDDO
2939
2940       DO jsl = 1, nslm
2941          ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42
2942          nvan_mod = n0 + (nvan(ji)-n0) * nfact(jsl,ji)
2943          avan_mod = a0 + (avan(ji)-a0) * afact(jsl,ji)
2944          m = un - un / nvan_mod
2945          ! Creation of arrays for SP-MIP output by landpoint
2946          nvan_mod_tab(jsl,ji) = nvan_mod
2947          avan_mod_tab(jsl,ji) = avan_mod
2948          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(ji) * kfact(jsl,ji)
2949          DO ii = imax,imin,-1
2950             frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
2951             k_lin(ii,jsl,ji) = ks(ji) * kfact(jsl,ji) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
2952          ENDDO
2953
2954          ! k_lin should not be zero, nor too small
2955          ! We track iiref, the bin under which mc is too small and we may get zero k_lin
2956          !salma: ji replaced with ii and jiref replaced with iiref and jsc with ji
2957          ii=imax-1
2958          DO WHILE ((k_lin(ii,jsl,ji) > 1.e-32) .and. (ii>0))
2959             iiref=ii
2960             ii=ii-1
2961          ENDDO
2962          DO ii=iiref-1,imin,-1
2963             k_lin(ii,jsl,ji)=k_lin(ii+1,jsl,ji)/10.
2964          ENDDO
2965
2966          DO ii = imin,imax-1 ! ii=1,50
2967             ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin
2968             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))
2969             b_lin(ii,jsl,ji)  = k_lin(ii,jsl,ji) - a_lin(ii,jsl,ji)*mc_lin(ii,ji)
2970
2971             ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta)
2972             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin
2973             IF (ii.NE.imin .AND. ii.NE.imax-1) THEN
2974                frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
2975                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) *  &
2976                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * &
2977                     (  frac**(-un/m) -un ) ** (-m)
2978                frac=MIN(un,(mc_lin(ii+1,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
2979                d_lin(ii+1,jsl,ji) =(k_lin(ii+1,jsl,ji) / (avan_mod*m*nvan_mod))*&
2980                     ( (frac**(-un/m))/(mc_lin(ii+1,ji)-mcr(ji)) ) * &
2981                     (  frac**(-un/m) -un ) ** (-m)
2982                d_lin(ii,jsl,ji) = undemi * (d_lin(ii,jsl,ji)+d_lin(ii+1,jsl,ji))
2983             ELSE IF(ii.EQ.imax-1) THEN
2984                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * &
2985                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) *  &
2986                     (  frac**(-un/m) -un ) ** (-m)
2987             ENDIF
2988          ENDDO
2989
2990          ! Special case for ii=imin
2991          d_lin(imin,jsl,ji) = d_lin(imin+1,jsl,ji)/1000.
2992
2993          ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations
2994          ! We don't want d_lin = zero
2995          DO ii=iiref-1,imin,-1
2996             d_lin(ii,jsl,ji)=d_lin(ii+1,jsl,ji)/10.
2997          ENDDO
2998
2999       ENDDO
3000    ENDDO
3001
3002
3003    ! Output of alphavg and nvg at each node for SP-MIP
3004    DO jsl = 1, nslm
3005       alphavg(:,jsl) = avan_mod_tab(jsl,:)*1000. ! from mm-1 to m-1
3006       nvg(:,jsl) = nvan_mod_tab(jsl,:)
3007    ENDDO
3008    CALL xios_orchidee_send_field("alphavg",alphavg) ! in m-1
3009    CALL xios_orchidee_send_field("nvg",nvg) ! unitless
3010
3011    !! 5 Water reservoir initialisation
3012    !
3013!!$    DO jst = 1,nstm
3014!!$       DO ji = 1, kjpindex
3015!!$          mx_eau_var(ji) = mx_eau_var(ji) + soiltile(ji,jst)*&
3016!!$               &   zmaxh*mille*mcs(njsc(ji))
3017!!$       END DO
3018!!$    END DO
3019
3020    mx_eau_var(:) = zero
3021    mx_eau_var(:) = zmaxh*mille*mcs(:)
3022
3023    DO ji = 1,kjpindex
3024       IF (vegtot(ji) .LE. zero) THEN
3025          mx_eau_var(ji) = mx_eau_nobio*zmaxh
3026          ! Aurelien: what does vegtot=0 mean? is it like frac_nobio=1? But if 0<frac_nobio<1 ???
3027       ENDIF
3028
3029    END DO
3030
3031    ! Compute the litter humidity, shumdiag and fry
3032    shumdiag_perma(:,:) = zero
3033    humtot(:) = zero
3034    tmc(:,:) = zero
3035
3036    ! Loop on soiltiles to compute the variables (ji,jst)
3037    DO jst=1,nstm
3038       DO ji = 1, kjpindex
3039          tmcs(ji,jst)=zmaxh* mille*mcs(ji)
3040          tmcr(ji,jst)=zmaxh* mille*mcr(ji)
3041          tmcfc(ji,jst)=zmaxh* mille*mcfc(ji)
3042          tmcw(ji,jst)=zmaxh* mille*mcw(ji)
3043       ENDDO
3044    ENDDO
3045
3046    ! The total soil moisture for each soiltile:
3047    DO jst=1,nstm
3048       DO ji=1,kjpindex
3049          tmc(ji,jst)= dz(2) * ( trois*mc(ji,1,jst)+ mc(ji,2,jst))/huit
3050       END DO
3051    ENDDO
3052
3053    DO jst=1,nstm
3054       DO jsl=2,nslm-1
3055          DO ji=1,kjpindex
3056             tmc(ji,jst) = tmc(ji,jst) + dz(jsl) * ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
3057                  & + dz(jsl+1)*(trois*mc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
3058          END DO
3059       END DO
3060    ENDDO
3061
3062    DO jst=1,nstm
3063       DO ji=1,kjpindex
3064          tmc(ji,jst) = tmc(ji,jst) +  dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
3065          tmc(ji,jst) = tmc(ji,jst) + water2infilt(ji,jst)
3066       ENDDO
3067    END DO
3068
3069!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.
3070!    ! If veget has been updated before restart (with LAND USE or DGVM),
3071!    ! tmc and mc must be modified with respect to humtot conservation.
3072!   CALL hydrol_tmc_update ( kjpindex, veget_max, soiltile, qsintveg)
3073
3074    ! The litter variables:
3075    ! level 1
3076    DO jst=1,nstm
3077       DO ji=1,kjpindex
3078          tmc_litter(ji,jst) = dz(2) * (trois*mcl(ji,1,jst)+mcl(ji,2,jst))/huit
3079          tmc_litter_wilt(ji,jst) = dz(2) * mcw(ji) / deux
3080          tmc_litter_res(ji,jst) = dz(2) * mcr(ji) / deux
3081          tmc_litter_field(ji,jst) = dz(2) * mcfc(ji) / deux
3082          tmc_litter_sat(ji,jst) = dz(2) * mcs(ji) / deux
3083          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(njsc(ji)) / deux
3084          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(njsc(ji)) / deux
3085       ENDDO
3086    END DO
3087    ! sum from level 2 to 4
3088    DO jst=1,nstm
3089       DO jsl=2,4
3090          DO ji=1,kjpindex
3091             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * &
3092                  & ( trois*mcl(ji,jsl,jst) + mcl(ji,jsl-1,jst))/huit &
3093                  & + dz(jsl+1)*(trois*mcl(ji,jsl,jst) + mcl(ji,jsl+1,jst))/huit
3094             tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + &
3095                  &(dz(jsl)+ dz(jsl+1))*&
3096                  & mcw(ji)/deux
3097             tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + &
3098                  &(dz(jsl)+ dz(jsl+1))*&
3099                  & mcr(ji)/deux
3100             tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + &
3101                  &(dz(jsl)+ dz(jsl+1))* &
3102                  & mcs(ji)/deux
3103             tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + &
3104                  & (dz(jsl)+ dz(jsl+1))* &
3105                  & mcfc(ji)/deux
3106             tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + &
3107                  &(dz(jsl)+ dz(jsl+1))* &
3108                  & mc_awet(njsc(ji))/deux
3109             tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + &
3110                  & (dz(jsl)+ dz(jsl+1))* &
3111                  & mc_adry(njsc(ji))/deux
3112          END DO
3113       END DO
3114    END DO
3115
3116
3117    DO jst=1,nstm
3118       DO ji=1,kjpindex
3119          ! here we set that humrelv=0 in PFT1
3120         humrelv(ji,1,jst) = zero
3121       ENDDO
3122    END DO
3123
3124
3125    ! Calculate shumdiag_perma for thermosoil
3126    ! Use resdist instead of soiltile because we here need to have
3127    ! shumdiag_perma at the value from previous time step.
3128    ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma
3129    ! (based on resdist=soiltile from previous timestep, but normally equal to soiltile)
3130    ! For consistency with hydrol_soil, we want to calculate a grid-cell average
3131    soilmoist(:,:) = zero
3132    DO jst=1,nstm
3133       DO ji=1,kjpindex
3134          soilmoist(ji,1) = soilmoist(ji,1) + resdist(ji,jst) * &
3135               dz(2) * ( trois*mc(ji,1,jst) + mc(ji,2,jst) )/huit
3136          DO jsl = 2,nslm-1
3137             soilmoist(ji,jsl) = soilmoist(ji,jsl) + resdist(ji,jst) * &
3138                  ( dz(jsl) * (trois*mc(ji,jsl,jst)+mc(ji,jsl-1,jst))/huit &
3139                  + dz(jsl+1) * (trois*mc(ji,jsl,jst)+mc(ji,jsl+1,jst))/huit )
3140          END DO
3141          soilmoist(ji,nslm) = soilmoist(ji,nslm) + resdist(ji,jst) * &
3142               dz(nslm) * (trois*mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
3143       ENDDO
3144    ENDDO
3145    DO ji=1,kjpindex
3146        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average
3147    ENDDO
3148
3149    ! -- shumdiag_perma for restart
3150   !  For consistency with hydrol_soil, we want to calculate a grid-cell average
3151    DO jsl = 1, nslm
3152       DO ji=1,kjpindex
3153          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji))
3154          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero)
3155       ENDDO
3156    ENDDO
3157
3158    ! Calculate drysoil_frac if it was not found in the restart file
3159    ! For simplicity, we set drysoil_frac to 0.5 in this case
3160    IF (ALL(drysoil_frac(:) == val_exp)) THEN
3161       DO ji=1,kjpindex
3162          drysoil_frac(ji) = 0.5
3163       END DO
3164    END IF
3165
3166    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in
3167    !! thermosoil for the thermal conductivity.
3168    ! These values are only used in thermosoil_init in absence of a restart file
3169
3170    mc_layh(:,:) = zero
3171    mcl_layh(:,:) = zero
3172     
3173    DO jst=1,nstm
3174       DO jsl=1,nslm
3175          DO ji=1,kjpindex
3176            mc_layh(ji,jsl) = mc_layh(ji,jsl) + mc(ji,jsl,jst) * resdist(ji,jst)  * vegtot_old(ji)
3177            mcl_layh(ji,jsl) = mcl_layh(ji,jsl) + mcl(ji,jsl,jst) * resdist(ji,jst) * vegtot_old(ji)
3178         ENDDO
3179      END DO
3180    END DO
3181
3182    IF (printlev>=3) WRITE (numout,*) ' hydrol_var_init done '
3183
3184  END SUBROUTINE hydrol_var_init
3185
3186
3187
3188   
3189!! ================================================================================================================================
3190!! SUBROUTINE   : hydrol_canop
3191!!
3192!>\BRIEF        This routine computes canopy processes.
3193!!
3194!! DESCRIPTION  :
3195!! - 1 evaporation off the continents
3196!! - 1.1 The interception loss is take off the canopy.
3197!! - 1.2 precip_rain is shared for each vegetation type
3198!! - 1.3 Limits the effect and sum what receives soil
3199!! - 1.4 swap qsintveg to the new value
3200!!
3201!! RECENT CHANGE(S) : None
3202!!
3203!! MAIN OUTPUT VARIABLE(S) :
3204!!
3205!! REFERENCE(S) :
3206!!
3207!! FLOWCHART    : None
3208!! \n
3209!_ ================================================================================================================================
3210!_ hydrol_canop
3211
3212  SUBROUTINE hydrol_canop (kjpindex, precip_rain, vevapwet, veget_max, veget, qsintmax, &
3213       & qsintveg,precisol,tot_melt)
3214
3215    !
3216    ! interface description
3217    !
3218
3219    !! 0. Variable and parameter declaration
3220
3221    !! 0.1 Input variables
3222
3223    INTEGER(i_std), INTENT(in)                               :: kjpindex    !! Domain size
3224    ! input fields
3225    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_rain !! Rain precipitation
3226    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: vevapwet    !! Interception loss
3227    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: veget_max   !! max fraction of vegetation type
3228    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: veget       !! Fraction of vegetation type
3229    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: qsintmax    !! Maximum water on vegetation for interception
3230    REAL(r_std), DIMENSION  (kjpindex), INTENT (in)          :: tot_melt    !! Total melt
3231
3232    !! 0.2 Output variables
3233
3234    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out)       :: precisol    !! Water fallen onto the ground (throughfall+Totmelt)
3235
3236    !! 0.3 Modified variables
3237
3238    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout)     :: qsintveg    !! Water on vegetation due to interception
3239
3240    !! 0.4 Local variables
3241
3242    INTEGER(i_std)                                           :: ji, jv
3243    REAL(r_std), DIMENSION (kjpindex,nvm)                    :: zqsintvegnew
3244
3245!_ ================================================================================================================================
3246
3247    ! boucle sur les points continentaux
3248    ! calcul de qsintveg au pas de temps suivant
3249    ! par ajout du flux interception loss
3250    ! calcule par enerbil en fonction
3251    ! des calculs faits dans diffuco
3252    ! calcul de ce qui tombe sur le sol
3253    ! avec accumulation dans precisol
3254    ! essayer d'harmoniser le traitement du sol nu
3255    ! avec celui des differents types de vegetation
3256    ! fait si on impose qsintmax ( ,1) = 0.0
3257    !
3258    ! loop for continental subdomain
3259    !
3260    !
3261    !! 1 evaporation off the continents
3262    !
3263    !! 1.1 The interception loss is take off the canopy.
3264    DO jv=2,nvm
3265       qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv)
3266    END DO
3267
3268    !     It is raining :
3269    !! 1.2 precip_rain is shared for each vegetation type
3270    !
3271    qsintveg(:,1) = zero
3272    DO jv=2,nvm
3273       qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
3274    END DO
3275
3276    !
3277    !! 1.3 Limits the effect and sum what receives soil
3278    !
3279    precisol(:,1)=veget_max(:,1)*precip_rain(:)
3280    DO jv=2,nvm
3281       DO ji = 1, kjpindex
3282          zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv)) 
3283          precisol(ji,jv) = (veget(ji,jv)*throughfall_by_pft(jv)*precip_rain(ji)) + &
3284               qsintveg(ji,jv) - zqsintvegnew (ji,jv) + &
3285               (veget_max(ji,jv) - veget(ji,jv))*precip_rain(ji)
3286       ENDDO
3287    END DO
3288       
3289    ! Precisol is currently the same as throughfall, save it for diagnostics
3290    throughfall(:,:) = precisol(:,:)
3291
3292    DO jv=1,nvm
3293       DO ji = 1, kjpindex
3294          IF (vegtot(ji).GT.min_sechiba) THEN
3295             precisol(ji,jv) = precisol(ji,jv)+tot_melt(ji)*veget_max(ji,jv)/vegtot(ji)
3296          ENDIF
3297       ENDDO
3298    END DO
3299    !   
3300    !
3301    !! 1.4 swap qsintveg to the new value
3302    !
3303    DO jv=2,nvm
3304       qsintveg(:,jv) = zqsintvegnew (:,jv)
3305    END DO
3306
3307    IF (printlev>=3) WRITE (numout,*) ' hydrol_canop done '
3308
3309  END SUBROUTINE hydrol_canop
3310
3311
3312!! ================================================================================================================================
3313!! SUBROUTINE   : hydrol_vegupd
3314!!
3315!>\BRIEF        Vegetation update   
3316!!
3317!! DESCRIPTION  :
3318!!   The vegetation cover has changed and we need to adapt the reservoir distribution
3319!!   and the distribution of plants on different soil types.
3320!!   You may note that this occurs after evaporation and so on have been computed. It is
3321!!   not a problem as a new vegetation fraction will start with humrel=0 and thus will have no
3322!!   evaporation. If this is not the case it should have been caught above.
3323!!
3324!! - 1 Update of vegetation is it needed?
3325!! - 2 calculate water mass that we have to redistribute
3326!! - 3 put it into reservoir of plant whose surface area has grown
3327!! - 4 Soil tile gestion
3328!! - 5 update the corresponding masks
3329!!
3330!! RECENT CHANGE(S) : None
3331!!
3332!! MAIN OUTPUT VARIABLE(S) :
3333!!
3334!! REFERENCE(S) :
3335!!
3336!! FLOWCHART    : None
3337!! \n
3338!_ ================================================================================================================================
3339!_ hydrol_vegupd
3340
3341  SUBROUTINE hydrol_vegupd(kjpindex, veget, veget_max, soiltile, qsintveg, frac_bare, drain_upd, runoff_upd)
3342
3343
3344    !! 0. Variable and parameter declaration
3345
3346    !! 0.1 Input variables
3347
3348    ! input scalar
3349    INTEGER(i_std), INTENT(in)                            :: kjpindex 
3350    ! input fields
3351    REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in)    :: veget            !! New vegetation map
3352    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)     :: veget_max        !! Max. fraction of vegetation type
3353    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)   :: soiltile         !! Fraction of each soil tile within vegtot (0-1, unitless)
3354
3355    !! 0.2 Output variables
3356    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out)    :: frac_bare        !! Fraction(of veget_max) of bare soil
3357                                                                              !! in each vegetation type
3358    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: drain_upd        !! Change in drainage due to decrease in vegtot
3359                                                                              !! on mc [kg/m2/dt]
3360    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: runoff_upd       !! Change in runoff due to decrease in vegtot
3361                                                                              !! on water2infilt[kg/m2/dt]
3362   
3363
3364    !! 0.3 Modified variables
3365
3366    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)  :: qsintveg         !! Water on old vegetation
3367
3368    !! 0.4 Local variables
3369
3370    INTEGER(i_std)                                 :: ji,jv,jst
3371
3372!_ ================================================================================================================================
3373
3374    !! 1 If veget has been updated at last time step (with LAND USE or DGVM),
3375    !! tmc and mc must be modified with respect to humtot conservation.
3376    CALL hydrol_tmc_update ( kjpindex, veget_max, soiltile, qsintveg, drain_upd, runoff_upd)
3377
3378
3379    ! Compute the masks for veget
3380   
3381    mask_veget(:,:) = 0
3382    mask_soiltile(:,:) = 0
3383   
3384    DO jst=1,nstm
3385       DO ji = 1, kjpindex
3386          IF(soiltile(ji,jst) .GT. min_sechiba) THEN
3387             mask_soiltile(ji,jst) = 1
3388          ENDIF
3389       END DO
3390    ENDDO
3391         
3392    DO jv = 1, nvm
3393       DO ji = 1, kjpindex
3394          IF(veget_max(ji,jv) .GT. min_sechiba) THEN
3395             mask_veget(ji,jv) = 1
3396          ENDIF
3397       END DO
3398    END DO
3399
3400    ! Compute vegetmax_soil
3401    vegetmax_soil(:,:,:) = zero
3402    DO jv = 1, nvm
3403       jst = pref_soil_veg(jv)
3404       DO ji=1,kjpindex
3405          ! for veget distribution used in sechiba via humrel
3406          IF (mask_soiltile(ji,jst).GT.0 .AND. vegtot(ji) > min_sechiba) THEN
3407             vegetmax_soil(ji,jv,jst)=veget_max(ji,jv)/soiltile(ji,jst)
3408          ENDIF
3409       ENDDO
3410    ENDDO
3411
3412    ! Calculate frac_bare (previosly done in slowproc_veget)
3413    DO ji =1, kjpindex
3414       IF( veget_max(ji,1) .GT. min_sechiba ) THEN
3415          frac_bare(ji,1) = un
3416       ELSE
3417          frac_bare(ji,1) = zero
3418       ENDIF
3419    ENDDO
3420    DO jv = 2, nvm
3421       DO ji =1, kjpindex
3422          IF( veget_max(ji,jv) .GT. min_sechiba ) THEN
3423             frac_bare(ji,jv) = un - veget(ji,jv)/veget_max(ji,jv)
3424          ELSE
3425             frac_bare(ji,jv) = zero
3426          ENDIF
3427       ENDDO
3428    ENDDO
3429
3430    ! Tout dans cette routine est maintenant certainement obsolete (veget_max etant constant) en dehors des lignes
3431    ! suivantes et le calcul de frac_bare:
3432    frac_bare_ns(:,:) = zero
3433    DO jst = 1, nstm
3434       DO jv = 1, nvm
3435          DO ji =1, kjpindex
3436             IF(vegtot(ji) .GT. min_sechiba) THEN
3437                frac_bare_ns(ji,jst) = frac_bare_ns(ji,jst) + vegetmax_soil(ji,jv,jst) * frac_bare(ji,jv) / vegtot(ji)
3438             ENDIF
3439          END DO
3440       ENDDO
3441    END DO
3442   
3443    IF (printlev>=3) WRITE (numout,*) ' hydrol_vegupd done '
3444
3445  END SUBROUTINE hydrol_vegupd
3446
3447
3448!! ================================================================================================================================
3449!! SUBROUTINE   : hydrol_flood
3450!!
3451!>\BRIEF        This routine computes the evolution of the surface reservoir (floodplain). 
3452!!
3453!! DESCRIPTION  :
3454!! - 1 Take out vevapflo from the reservoir and transfer the remaining to subsinksoil
3455!! - 2 Compute the total flux from floodplain floodout (transfered to routing)
3456!! - 3 Discriminate between precip over land and over floodplain
3457!!
3458!! RECENT CHANGE(S) : None
3459!!
3460!! MAIN OUTPUT VARIABLE(S) :
3461!!
3462!! REFERENCE(S) :
3463!!
3464!! FLOWCHART    : None
3465!! \n
3466!_ ================================================================================================================================
3467!_ hydrol_flood
3468
3469  SUBROUTINE hydrol_flood (kjpindex, vevapflo, flood_frac, flood_res, floodout)
3470
3471    !! 0. Variable and parameter declaration
3472
3473    !! 0.1 Input variables
3474
3475    ! input scalar
3476    INTEGER(i_std), INTENT(in)                               :: kjpindex         !!
3477    ! input fields
3478    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: flood_frac       !! Fraction of floodplains in grid box
3479
3480    !! 0.2 Output variables
3481
3482    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: floodout         !! Flux to take out from floodplains
3483
3484    !! 0.3 Modified variables
3485
3486    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: flood_res        !! Floodplains reservoir estimate
3487    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: vevapflo         !! Evaporation over floodplains
3488
3489    !! 0.4 Local variables
3490
3491    INTEGER(i_std)                                           :: ji, jv           !! Indices
3492    REAL(r_std), DIMENSION (kjpindex)                        :: temp             !!
3493
3494!_ ================================================================================================================================
3495    !-
3496    !! 1 Take out vevapflo from the reservoir and transfer the remaining to subsinksoil
3497    !-
3498    DO ji = 1,kjpindex
3499       temp(ji) = MIN(flood_res(ji), vevapflo(ji))
3500    ENDDO
3501    DO ji = 1,kjpindex
3502       flood_res(ji) = flood_res(ji) - temp(ji)
3503       subsinksoil(ji) = subsinksoil(ji) + vevapflo(ji) - temp(ji)
3504       vevapflo(ji) = temp(ji)
3505    ENDDO
3506
3507    !-
3508    !! 2 Compute the total flux from floodplain floodout (transfered to routing)
3509    !-
3510    DO ji = 1,kjpindex
3511       floodout(ji) = vevapflo(ji) - flood_frac(ji) * SUM(precisol(ji,:))
3512    ENDDO
3513
3514    !-
3515    !! 3 Discriminate between precip over land and over floodplain
3516    !-
3517    DO jv=1, nvm
3518       DO ji = 1,kjpindex
3519          precisol(ji,jv) = precisol(ji,jv) * (1 - flood_frac(ji))
3520       ENDDO
3521    ENDDO 
3522
3523    IF (printlev>=3) WRITE (numout,*) ' hydrol_flood done'
3524
3525  END SUBROUTINE hydrol_flood
3526
3527!! ================================================================================================================================
3528!! SUBROUTINE   : hydrol_soil
3529!!
3530!>\BRIEF        This routine computes soil processes with CWRR scheme (Richards equation solved by finite differences).
3531!! Note that the water fluxes are in kg/m2/dt_sechiba.
3532!!
3533!! DESCRIPTION  :
3534!! 0. Initialisation, and split 2d variables to 3d variables, per soil tile
3535!! -- START MAIN LOOP (prognostic loop to update mc and mcl) OVER SOILTILES
3536!! 1. FIRSTLY, WE CHANGE MC BASED ON EXTERNAL FLUXES, ALL APPLIED AT THE SOIL SURFACE
3537!! 1.1 Reduces water2infilt and water2extract to their difference
3538!! 1.2 To remove water2extract (including bare soilevaporation) from top layer
3539!! 1.3 Infiltration
3540!! 1.4 Reinfiltration of surface runoff : compute temporary surface water and extract from runoff
3541!! 2. SECONDLY, WE UPDATE MC FROM DIFFUSION, INCLUDING DRAINAGE AND ROOTSINK
3542!!    This will act on mcl (liquid water content) only
3543!! 2.1 K and D are recomputed after infiltration
3544!! 2.2 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme
3545!! 2.3 We define mcl (liquid water content) based on mc and profil_froz_hydro_ns
3546!! 2.4 We calculate the total SM at the beginning of the routine tridiag for water conservation check
3547!! 2.5 Defining where diffusion is solved : everywhere
3548!! 2.6 We define the system of linear equations for mcl redistribution
3549!! 2.7 Solves diffusion equations
3550!! 2.8 Computes drainage = bottom boundary condition, consistent with rhs(ji,jsl=nslm)
3551!! 2.9 For water conservation check during redistribution, we calculate the total liquid SM
3552!!     at the end of the routine tridiag, and we compare the difference with the flux...
3553!! 3. AFTER DIFFUSION/REDISTRIBUTION
3554!! 3.1 Updating mc, as all the following checks against saturation will compare mc to mcs
3555!! 3.2 Correct here the possible over-saturation values (subroutine hydrol_soil_smooth_over_mcs2 acts on mc)
3556!!     Here hydrol_soil_smooth_over_mcs2 discard all excess as ru_corr_ns, oriented to either ru_ns or dr_ns
3557!! 3.3 Negative runoff is reported to drainage
3558!! 3.4 Optional block to force saturation below zwt_force
3559!! 3.5 Diagnosing the effective water table depth
3560!! 3.6 Diagnose under_mcr to adapt water stress calculation below
3561!! 4. At the end of the prognostic calculations, we recompute important moisture variables
3562!! 4.1 Total soil moisture content (water2infilt added below)
3563!! 4.2 mcl is a module variable; we update it here for calculating bare soil evaporation,
3564!! 5. Optional check of the water balance of soil column (if check_cwrr)
3565!! 5.1 Computation of the vertical water fluxes
3566!! 6. SM DIAGNOSTICS FOR OTHER ROUTINES, MODULES, OR NEXT STEP
3567!! 6.1 Total soil moisture, soil moisture at litter levels, soil wetness, us, humrelv, vesgtressv
3568!! 6.2 We need to turn off evaporation when is_under_mcr
3569!! 6.3 Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in thermosoil
3570!! 6.4 The hydraulic conductivities exported here are the ones used in the diffusion/redistribution
3571!! -- ENDING THE MAIN LOOP ON SOILTILES
3572!! 7. Summing 3d variables into 2d variables
3573!! 8. XIOS export of local variables, including water conservation checks
3574!! 9. COMPUTING EVAP_BARE_LIM_NS FOR NEXT TIME STEP, WITH A LOOP ON SOILTILES
3575!!    The principle is to run a dummy integration of the water redistribution scheme
3576!!    to check if the SM profile can sustain a potential evaporation.
3577!!    If not, the dummy integration is redone from the SM profile of the end of the normal integration,
3578!!    with a boundary condition leading to a very severe water limitation: mc(1)=mcr
3579!! 10. evap_bar_lim is the grid-cell scale beta
3580!!
3581!! RECENT CHANGE(S) : 2016 by A. Ducharne
3582!!
3583!! MAIN OUTPUT VARIABLE(S) :
3584!!
3585!! REFERENCE(S) :
3586!!
3587!! FLOWCHART    : None
3588!! \n
3589!_ ================================================================================================================================
3590!_ hydrol_soil
3591  SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
3592       kjpindex, veget_max, soiltile, njsc, reinf_slope, &
3593       & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
3594       & returnflow, reinfiltration, irrigation, &
3595       & tot_melt, evap_bare_lim, evap_bare_lim_ns, shumdiag, shumdiag_perma,&
3596       & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, &
3597       & stempdiag,snow, &
3598       & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh)
3599    !
3600    ! interface description
3601
3602    !! 0. Variable and parameter declaration
3603
3604    !! 0.1 Input variables
3605   
3606    INTEGER(i_std), INTENT(in)                               :: kjpindex 
3607    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in)       :: veget_max        !! Map of max vegetation types [-]
3608    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: njsc             !! Index of the dominant soil textural class
3609                                                                                 !!   in the grid cell (1-nscm, unitless)
3610   
3611    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1})
3612    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless)
3613    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1})
3614    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
3615    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
3616    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
3617    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
3618   
3619    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)      :: soiltile         !! Fraction of each soil tile within vegtot (0-1, unitless)
3620    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: transpir         !! Transpiration 
3621                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3622    REAL(r_std), DIMENSION (kjpindex), INTENT (in)           :: reinf_slope      !! Fraction of surface runoff that reinfiltrates
3623                                                                                 !!  (unitless, [0-1])
3624    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: returnflow       !! Water returning to the soil from the bottom
3625                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3626    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: reinfiltration   !! Water returning to the top of the soil
3627                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3628    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: irrigation       !! Irrigation
3629                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3630    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: evapot           !! Potential evaporation
3631                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3632    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: evapot_penm      !! Potential evaporation "Penman" (Milly's correction)
3633                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3634    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: tot_melt         !! Total melt from snow and ice
3635                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3636    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)       :: stempdiag        !! Diagnostic temp profile from thermosoil
3637    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: snow             !! Snow mass
3638                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3639    REAL(r_std), DIMENSION (kjpindex,nsnow),INTENT(in)       :: snowdz           !! Snow depth (m)
3640    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: tot_bare_soil    !! Total evaporating bare soil fraction
3641                                                                                 !!  (unitless, [0-1])
3642    REAL(r_std),DIMENSION (kjpindex), INTENT(in)             :: u,v              !! Horizontal wind speed
3643    REAL(r_std),DIMENSION (kjpindex), INTENT(in)             :: tq_cdrag         !! Surface drag coefficient
3644
3645    !! 0.2 Output variables
3646
3647    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: runoff           !! Surface runoff
3648                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3649    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: drainage         !! Drainage
3650                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3651    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: evap_bare_lim    !! Limitation factor (beta) for bare soil evaporation 
3652                                                                                 !! on each soil column (unitless, [0-1])
3653    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out)      :: evap_bare_lim_ns !! Limitation factor (beta) for bare soil evaporation 
3654                                                                                 !! on each soil column (unitless, [0-1])
3655    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiag         !! Relative soil moisture in each diag soil layer
3656                                                                                 !! with respect to (mcfc-mcw) (unitless, [0-1])
3657    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiag_perma   !! Percent of porosity filled with water (mc/mcs)
3658                                                                                 !! in each diag soil layer (for the thermal computations)
3659                                                                                 !! (unitless, [0-1])
3660    REAL(r_std), DIMENSION (kjpindex), INTENT (out)          :: k_litt           !! Litter approximated hydraulic conductivity
3661                                                                                 !!  @tex $(mm d^{-1})$ @endtex
3662    REAL(r_std), DIMENSION (kjpindex), INTENT (out)          :: litterhumdiag    !! Mean of soil_wet_litter across soil tiles
3663                                                                                 !! (unitless, [0-1])
3664    REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out)      :: vegstress        !! Veg. moisture stress (only for vegetation
3665                                                                                 !! growth) (unitless, [0-1])
3666    REAL(r_std), DIMENSION (kjpindex), INTENT (out)          :: drysoil_frac     !! Function of the litter humidity
3667    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: mc_layh          !! Volumetric water content (liquid + ice) for each soil layer
3668                                                                                 !! averaged over the mesh (for thermosoil)
3669                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex
3670    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: mcl_layh         !! Volumetric liquid water content for each soil layer
3671                                                                                 !! averaged over the mesh (for thermosoil)
3672                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex
3673    !! 0.3 Modified variables
3674
3675    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: vevapnu          !! Bare soil evaporation
3676                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3677    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (inout)    :: humrel           !! Relative humidity (0-1, dimensionless)
3678
3679    !! 0.4 Local variables
3680
3681    INTEGER(i_std)                                 :: ji, jv, jsl, jst           !! Indices
3682    REAL(r_std), PARAMETER                         :: frac_mcs = 0.66            !! Temporary depth
3683    REAL(r_std), DIMENSION(kjpindex)               :: temp                       !! Temporary value for fluxes
3684    REAL(r_std), DIMENSION(kjpindex)               :: tmcold                     !! Total SM at beginning of hydrol_soil (kg/m2)
3685    REAL(r_std), DIMENSION(kjpindex)               :: tmcint                     !! Ancillary total SM (kg/m2)
3686    REAL(r_std), DIMENSION(kjpindex,nslm)          :: mcint                      !! To save mc values for future use
3687    REAL(r_std), DIMENSION(kjpindex,nslm)          :: mclint                     !! To save mcl values for future use
3688    LOGICAL, DIMENSION(kjpindex,nstm)              :: is_under_mcr               !! Identifies under residual soil moisture points
3689    LOGICAL, DIMENSION(kjpindex)                   :: is_over_mcs                !! Identifies over saturated soil moisture points
3690    REAL(r_std), DIMENSION(kjpindex)               :: deltahum,diff              !!
3691    LOGICAL(r_std), DIMENSION(kjpindex)            :: test                       !!
3692    REAL(r_std), DIMENSION(kjpindex)               :: water2extract              !! Water flux to be extracted at the soil surface
3693                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3694    REAL(r_std), DIMENSION(kjpindex)               :: returnflow_soil            !! Water from the routing back to the bottom of
3695                                                                                 !! the soil @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3696    REAL(r_std), DIMENSION(kjpindex)               :: reinfiltration_soil        !! Water from the routing back to the top of the
3697                                                                                 !! soil @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3698    REAL(r_std), DIMENSION(kjpindex)               :: irrigation_soil            !! Water from irrigation returning to soil moisture
3699                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
3700    REAL(r_std), DIMENSION(kjpindex)               :: flux_infilt                !! Water to infiltrate
3701                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3702    REAL(r_std), DIMENSION(kjpindex)               :: flux_bottom                !! Flux at bottom of the soil column
3703                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3704    REAL(r_std), DIMENSION(kjpindex)               :: flux_top                   !! Flux at top of the soil column (for bare soil evap)
3705                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3706    REAL(r_std), DIMENSION (kjpindex,nstm)         :: qinfilt_ns                 !! Effective infiltration flux per soil tile
3707                                                                                 !!  @tex $(kg m^{-2})$ @endtex   
3708    REAL(r_std), DIMENSION (kjpindex)              :: qinfilt                    !! Effective infiltration flux 
3709                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3710    REAL(r_std), DIMENSION (kjpindex,nstm)         :: ru_infilt_ns               !! Surface runoff from hydrol_soil_infilt per soil tile
3711                                                                                 !!  @tex $(kg m^{-2})$ @endtex   
3712    REAL(r_std), DIMENSION (kjpindex)              :: ru_infilt                  !! Surface runoff from hydrol_soil_infilt
3713                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3714    REAL(r_std), DIMENSION (kjpindex,nstm)         :: ru_corr_ns                 !! Surface runoff produced to correct excess per soil tile
3715                                                                                 !!  @tex $(kg m^{-2})$ @endtex
3716    REAL(r_std), DIMENSION (kjpindex)              :: ru_corr                    !! Surface runoff produced to correct excess
3717                                                                                 !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex 
3718    REAL(r_std), DIMENSION (kjpindex,nstm)         :: ru_corr2_ns                !! Correction of negative surface runoff per soil tile
3719