Ticket #412: hydrol_bugksat2.f90

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