Ticket #412: hydrol_bugksat.f90

File hydrol_bugksat.f90, 361.7 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       CALL xios_orchidee_send_field("kk_moy",kk_moy)
921    END IF
922    CALL xios_orchidee_send_field("profil_froz_hydro_ns", profil_froz_hydro_ns)
923
924
925
926    !! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
927    humtot_lut(:,:)=0
928    humtot_top_lut(:,:)=0
929    mrro_lut(:,:)=0
930    DO jv=1,nvm
931       jst=pref_soil_veg(jv) ! soil tile index
932       IF (natural(jv)) THEN
933          humtot_lut(:,id_psl) = humtot_lut(:,id_psl) + tmc(:,jst)*veget_max(:,jv)
934          humtot_top_lut(:,id_psl) = humtot_top_lut(:,id_psl) + tmc_top(:,jst)*veget_max(:,jv)
935          mrro_lut(:,id_psl) = mrro_lut(:,id_psl) + (dr_ns(:,jst)+ru_ns(:,jst))*veget_max(:,jv)
936       ELSE
937          humtot_lut(:,id_crp) = humtot_lut(:,id_crp) + tmc(:,jst)*veget_max(:,jv)
938          humtot_top_lut(:,id_crp) = humtot_top_lut(:,id_crp) + tmc_top(:,jst)*veget_max(:,jv)
939          mrro_lut(:,id_crp) = mrro_lut(:,id_crp) + (dr_ns(:,jst)+ru_ns(:,jst))*veget_max(:,jv)
940       ENDIF
941    END DO
942
943    WHERE (fraclut(:,id_psl)>min_sechiba)
944       humtot_lut(:,id_psl) = humtot_lut(:,id_psl)/fraclut(:,id_psl)
945       humtot_top_lut(:,id_psl) = humtot_top_lut(:,id_psl)/fraclut(:,id_psl)
946       mrro_lut(:,id_psl) = mrro_lut(:,id_psl)/fraclut(:,id_psl)/dt_sechiba
947    ELSEWHERE
948       humtot_lut(:,id_psl) = val_exp
949       humtot_top_lut(:,id_psl) = val_exp
950       mrro_lut(:,id_psl) = val_exp
951    END WHERE
952    WHERE (fraclut(:,id_crp)>min_sechiba)
953       humtot_lut(:,id_crp) = humtot_lut(:,id_crp)/fraclut(:,id_crp)
954       humtot_top_lut(:,id_crp) = humtot_top_lut(:,id_crp)/fraclut(:,id_crp)
955       mrro_lut(:,id_crp) = mrro_lut(:,id_crp)/fraclut(:,id_crp)/dt_sechiba
956    ELSEWHERE
957       humtot_lut(:,id_crp) = val_exp
958       humtot_top_lut(:,id_crp) = val_exp
959       mrro_lut(:,id_crp) = val_exp
960    END WHERE
961
962    humtot_lut(:,id_pst) = val_exp
963    humtot_lut(:,id_urb) = val_exp
964    humtot_top_lut(:,id_pst) = val_exp
965    humtot_top_lut(:,id_urb) = val_exp
966    mrro_lut(:,id_pst) = val_exp
967    mrro_lut(:,id_urb) = val_exp
968
969    CALL xios_orchidee_send_field("humtot_lut",humtot_lut)
970    CALL xios_orchidee_send_field("humtot_top_lut",humtot_top_lut)
971    CALL xios_orchidee_send_field("mrro_lut",mrro_lut)
972
973
974    IF ( .NOT. almaoutput ) THEN
975       CALL histwrite_p(hist_id, 'frac_bare', kjit, frac_bare, kjpindex*nvm, indexveg)
976
977       DO jst=1,nstm
978          ! var_name= "mc_1" ... "mc_3"
979          WRITE (var_name,"('moistc_',i1)") jst
980          CALL histwrite_p(hist_id, TRIM(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
981
982          ! var_name= "kfactroot_1" ... "kfactroot_3"
983          WRITE (var_name,"('kfactroot_',i1)") jst
984          CALL histwrite_p(hist_id, TRIM(var_name), kjit, kfact_root(:,:,jst), kjpindex*nslm, indexlayer)
985
986          ! var_name= "vegetsoil_1" ... "vegetsoil_3"
987          WRITE (var_name,"('vegetsoil_',i1)") jst
988          CALL histwrite_p(hist_id, TRIM(var_name), kjit,vegetmax_soil(:,:,jst), kjpindex*nvm, indexveg)
989       ENDDO
990       CALL histwrite_p(hist_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
991       CALL histwrite_p(hist_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
992       CALL histwrite_p(hist_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
993       CALL histwrite_p(hist_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
994       CALL histwrite_p(hist_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
995       ! mrso is a perfect duplicate of humtot
996       CALL histwrite_p(hist_id, 'humtot', kjit, humtot, kjpindex, index)
997       CALL histwrite_p(hist_id, 'mrso', kjit, humtot, kjpindex, index)
998       CALL histwrite_p(hist_id, 'mrsos', kjit, humtot_top, kjpindex, index)
999       njsc_tmp(:)=njsc(:)
1000       CALL histwrite_p(hist_id, 'soilindex', kjit, njsc_tmp, kjpindex, index)
1001       CALL histwrite_p(hist_id, 'humrel',   kjit, humrel,   kjpindex*nvm, indexveg)
1002       CALL histwrite_p(hist_id, 'drainage', kjit, drainage, kjpindex, index)
1003       ! NB! According to histdef in intersurf, the variables 'runoff' and 'mrros' have different units
1004       CALL histwrite_p(hist_id, 'runoff', kjit, runoff, kjpindex, index)
1005       CALL histwrite_p(hist_id, 'mrros', kjit, runoff, kjpindex, index)
1006       histvar(:)=(runoff(:)+drainage(:))
1007       CALL histwrite_p(hist_id, 'mrro', kjit, histvar, kjpindex, index)
1008       CALL histwrite_p(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
1009       CALL histwrite_p(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
1010
1011       histvar(:)=(precip_rain(:)-SUM(throughfall(:,:),dim=2))
1012       CALL histwrite_p(hist_id, 'prveg', kjit, histvar, kjpindex, index)
1013
1014       CALL histwrite_p(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
1015       CALL histwrite_p(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
1016       CALL histwrite_p(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
1017       CALL histwrite_p(hist_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
1018       CALL histwrite_p(hist_id, 'shumdiag_perma',kjit,shumdiag_perma,kjpindex*nslm,indexnslm)
1019
1020       IF ( do_floodplains ) THEN
1021          CALL histwrite_p(hist_id, 'floodout', kjit, floodout, kjpindex, index)
1022       ENDIF
1023       !
1024       IF ( hist2_id > 0 ) THEN
1025          DO jst=1,nstm
1026             ! var_name= "mc_1" ... "mc_3"
1027             WRITE (var_name,"('moistc_',i1)") jst
1028             CALL histwrite_p(hist2_id, TRIM(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
1029
1030             ! var_name= "kfactroot_1" ... "kfactroot_3"
1031             WRITE (var_name,"('kfactroot_',i1)") jst
1032             CALL histwrite_p(hist2_id, TRIM(var_name), kjit, kfact_root(:,:,jst), kjpindex*nslm, indexlayer)
1033
1034             ! var_name= "vegetsoil_1" ... "vegetsoil_3"
1035             WRITE (var_name,"('vegetsoil_',i1)") jst
1036             CALL histwrite_p(hist2_id, TRIM(var_name), kjit,vegetmax_soil(:,:,jst), kjpindex*nvm, indexveg)
1037          ENDDO
1038          CALL histwrite_p(hist2_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
1039          CALL histwrite_p(hist2_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
1040          CALL histwrite_p(hist2_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
1041          CALL histwrite_p(hist2_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
1042          CALL histwrite_p(hist2_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
1043          ! mrso is a perfect duplicate of humtot
1044          CALL histwrite_p(hist2_id, 'humtot', kjit, humtot, kjpindex, index)
1045          CALL histwrite_p(hist2_id, 'mrso', kjit, humtot, kjpindex, index)
1046          CALL histwrite_p(hist2_id, 'mrsos', kjit, humtot_top, kjpindex, index)
1047          njsc_tmp(:)=njsc(:)
1048          CALL histwrite_p(hist2_id, 'soilindex', kjit, njsc_tmp, kjpindex, index)
1049          CALL histwrite_p(hist2_id, 'humrel',   kjit, humrel,   kjpindex*nvm, indexveg)
1050          CALL histwrite_p(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
1051          ! NB! According to histdef in intersurf, the variables 'runoff' and 'mrros' have different units
1052          CALL histwrite_p(hist2_id, 'runoff', kjit, runoff, kjpindex, index)
1053          CALL histwrite_p(hist2_id, 'mrros', kjit, runoff, kjpindex, index)
1054          histvar(:)=(runoff(:)+drainage(:))
1055          CALL histwrite_p(hist2_id, 'mrro', kjit, histvar, kjpindex, index)
1056
1057          IF ( do_floodplains ) THEN
1058             CALL histwrite_p(hist2_id, 'floodout', kjit, floodout, kjpindex, index)
1059          ENDIF
1060          CALL histwrite_p(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
1061          CALL histwrite_p(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
1062          CALL histwrite_p(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
1063          CALL histwrite_p(hist2_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
1064          CALL histwrite_p(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
1065          CALL histwrite_p(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
1066       ENDIF
1067    ELSE
1068       CALL histwrite_p(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
1069       CALL histwrite_p(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
1070       CALL histwrite_p(hist_id, 'Qs', kjit, runoff, kjpindex, index)
1071       CALL histwrite_p(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
1072       CALL histwrite_p(hist_id, 'Qsm', kjit, snowmelt, kjpindex, index)
1073       CALL histwrite_p(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
1074       CALL histwrite_p(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
1075       CALL histwrite_p(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
1076       !
1077       CALL histwrite_p(hist_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
1078       CALL histwrite_p(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
1079       !
1080       CALL histwrite_p(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
1081       CALL histwrite_p(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
1082       !
1083       IF (.NOT. ok_explicitsnow) CALL histwrite_p(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
1084       !
1085       IF ( hist2_id > 0 ) THEN
1086          CALL histwrite_p(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
1087          CALL histwrite_p(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
1088          CALL histwrite_p(hist2_id, 'Qs', kjit, runoff, kjpindex, index)
1089          CALL histwrite_p(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
1090          CALL histwrite_p(hist2_id, 'Qsm', kjit, snowmelt, kjpindex, index)
1091          CALL histwrite_p(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
1092          CALL histwrite_p(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
1093          CALL histwrite_p(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
1094          !
1095          CALL histwrite_p(hist2_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
1096          CALL histwrite_p(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
1097          !
1098          CALL histwrite_p(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
1099          CALL histwrite_p(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
1100          !
1101          IF (.NOT. ok_explicitsnow) CALL histwrite_p(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
1102       ENDIF
1103    ENDIF
1104
1105    IF (ok_freeze_cwrr) THEN
1106       CALL histwrite_p(hist_id, 'profil_froz_hydro', kjit,profil_froz_hydro , kjpindex*nslm, indexlayer)
1107       DO jst=1,nstm
1108          WRITE (var_name,"('profil_froz_hydro_',i1)") jst
1109          CALL histwrite_p(hist_id, TRIM(var_name), kjit, profil_froz_hydro_ns(:,:,jst), kjpindex*nslm, indexlayer)
1110       ENDDO
1111       CALL histwrite_p(hist_id, 'temp_hydro', kjit,temp_hydro , kjpindex*nslm, indexlayer)
1112       CALL histwrite_p(hist_id, 'kk_moy', kjit, kk_moy,kjpindex*nslm, indexlayer) ! averaged over soiltiles
1113    ENDIF
1114
1115
1116    ! Copy soilmoist into a local variable to be sent to thermosoil
1117    soilmoist_out(:,:) = soilmoist(:,:)
1118
1119    IF (printlev>=3) WRITE (numout,*) ' hydrol_main Done '
1120
1121  END SUBROUTINE hydrol_main
1122
1123
1124!! ================================================================================================================================
1125!! SUBROUTINE   : hydrol_finalize
1126!!
1127!>\BRIEF         
1128!!
1129!! DESCRIPTION : This subroutine writes the module variables and variables calculated in hydrol to restart file
1130!!
1131!! MAIN OUTPUT VARIABLE(S) :
1132!!
1133!! REFERENCE(S) :
1134!!
1135!! FLOWCHART    : None
1136!! \n
1137!_ ================================================================================================================================
1138
1139  SUBROUTINE hydrol_finalize( kjit,           kjpindex,   rest_id,  vegstress,  &
1140                              qsintveg,       humrel,     snow,     snow_age, snow_nobio, &
1141                              snow_nobio_age, snowrho,    snowtemp, snowdz,     &
1142                              snowheat,       snowgrain,  &
1143                              drysoil_frac, evap_bare_lim)
1144
1145    !! 0. Variable and parameter declaration
1146    !! 0.1 Input variables
1147    INTEGER(i_std), INTENT(in)                           :: kjit           !! Time step number
1148    INTEGER(i_std), INTENT(in)                           :: kjpindex       !! Domain size
1149    INTEGER(i_std),INTENT (in)                           :: rest_id        !! Restart file identifier
1150    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: vegstress      !! Veg. moisture stress (only for vegetation growth)
1151    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintveg       !! Water on vegetation due to interception
1152    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: humrel
1153    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow           !! Snow mass [Kg/m^2]
1154    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: snow_age       !! Snow age
1155    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio     !! Water balance on ice, lakes, .. [Kg/m^2]
1156    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio_age !! Snow age on ice, lakes, ...
1157    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowrho        !! Snow density
1158    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowtemp       !! Snow temperature
1159    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowdz         !! Snow layer thickness
1160    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowheat       !! Snow heat content
1161    REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in)  :: snowgrain      !! Snow grainsize
1162    REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: drysoil_frac   !! function of litter wetness
1163    REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: evap_bare_lim
1164
1165    !! 0.4 Local variables
1166    INTEGER(i_std)                                       :: jst, jsl
1167   
1168!_ ================================================================================================================================
1169
1170
1171    IF (printlev>=3) WRITE (numout,*) 'Write restart file with HYDROLOGIC variables '
1172
1173    DO jst=1,nstm
1174       ! var_name= "mc_1" ... "mc_3"
1175       WRITE (var_name,"('moistc_',i1)") jst
1176       CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mc(:,:,jst), 'scatter',  nbp_glo, index_g)
1177    END DO
1178
1179    DO jst=1,nstm
1180       ! var_name= "mcl_1" ... "mcl_3"
1181       WRITE (var_name,"('moistcl_',i1)") jst
1182       CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mcl(:,:,jst), 'scatter',  nbp_glo, index_g)
1183    END DO
1184   
1185    IF (ok_nudge_mc) THEN
1186       DO jst=1,nstm
1187          WRITE (var_name,"('mc_read_next_',i1)") jst
1188          CALL restput_p(rest_id, var_name, nbp_glo,  nslm, 1, kjit, mc_read_next(:,:,jst), 'scatter',  nbp_glo, index_g)
1189       END DO
1190    END IF
1191
1192    IF (ok_nudge_snow) THEN
1193       CALL restput_p(rest_id, 'snowdz_read_next', nbp_glo,  nsnow, 1, kjit, snowdz_read_next(:,:), &
1194            'scatter',  nbp_glo, index_g)
1195       CALL restput_p(rest_id, 'snowrho_read_next', nbp_glo,  nsnow, 1, kjit, snowrho_read_next(:,:), &
1196            'scatter',  nbp_glo, index_g)
1197       CALL restput_p(rest_id, 'snowtemp_read_next', nbp_glo,  nsnow, 1, kjit, snowtemp_read_next(:,:), &
1198            'scatter',  nbp_glo, index_g)
1199    END IF
1200
1201
1202           
1203    DO jst=1,nstm
1204       DO jsl=1,nslm
1205          ! var_name= "us_1_01" ... "us_3_11"
1206          WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
1207          CALL restput_p(rest_id, var_name, nbp_glo,nvm, 1,kjit,us(:,:,jst,jsl),'scatter',nbp_glo,index_g)
1208       END DO
1209    END DO
1210   
1211    CALL restput_p(rest_id, 'free_drain_coef', nbp_glo,   nstm, 1, kjit,  free_drain_coef, 'scatter',  nbp_glo, index_g)
1212    CALL restput_p(rest_id, 'zwt_force', nbp_glo,   nstm, 1, kjit,  zwt_force, 'scatter',  nbp_glo, index_g)
1213    CALL restput_p(rest_id, 'water2infilt', nbp_glo,   nstm, 1, kjit,  water2infilt, 'scatter',  nbp_glo, index_g)
1214    CALL restput_p(rest_id, 'ae_ns', nbp_glo,   nstm, 1, kjit,  ae_ns, 'scatter',  nbp_glo, index_g)
1215    CALL restput_p(rest_id, 'vegstress', nbp_glo,   nvm, 1, kjit,  vegstress, 'scatter',  nbp_glo, index_g)
1216    CALL restput_p(rest_id, 'snow', nbp_glo,   1, 1, kjit,  snow, 'scatter',  nbp_glo, index_g)
1217    CALL restput_p(rest_id, 'snow_age', nbp_glo,   1, 1, kjit,  snow_age, 'scatter',  nbp_glo, index_g)
1218    CALL restput_p(rest_id, 'snow_nobio', nbp_glo,   nnobio, 1, kjit,  snow_nobio, 'scatter', nbp_glo, index_g)
1219    CALL restput_p(rest_id, 'snow_nobio_age', nbp_glo,   nnobio, 1, kjit,  snow_nobio_age, 'scatter', nbp_glo, index_g)
1220    CALL restput_p(rest_id, 'qsintveg', nbp_glo, nvm, 1, kjit,  qsintveg, 'scatter',  nbp_glo, index_g)
1221    CALL restput_p(rest_id, 'evap_bare_lim_ns', nbp_glo, nstm, 1, kjit,  evap_bare_lim_ns, 'scatter',  nbp_glo, index_g)
1222    CALL restput_p(rest_id, 'evap_bare_lim', nbp_glo, 1, 1, kjit,  evap_bare_lim, 'scatter',  nbp_glo, index_g)
1223    CALL restput_p(rest_id, 'resdist', nbp_glo, nstm, 1, kjit,  resdist, 'scatter',  nbp_glo, index_g) 
1224    CALL restput_p(rest_id, 'vegtot_old', nbp_glo, 1, 1, kjit,  vegtot_old, 'scatter',  nbp_glo, index_g)           
1225    CALL restput_p(rest_id, 'drysoil_frac', nbp_glo,   1, 1, kjit, drysoil_frac, 'scatter', nbp_glo, index_g)
1226    CALL restput_p(rest_id, 'humrel', nbp_glo,   nvm, 1, kjit,  humrel, 'scatter',  nbp_glo, index_g)
1227
1228    CALL restput_p(rest_id, 'tot_watveg_beg', nbp_glo,  1, 1, kjit,  tot_watveg_beg, 'scatter',  nbp_glo, index_g)
1229    CALL restput_p(rest_id, 'tot_watsoil_beg', nbp_glo, 1, 1, kjit,  tot_watsoil_beg, 'scatter',  nbp_glo, index_g)
1230    CALL restput_p(rest_id, 'snow_beg', nbp_glo,        1, 1, kjit,  snow_beg, 'scatter',  nbp_glo, index_g)
1231   
1232   
1233    ! Write variables for explictsnow module to restart file
1234    IF (ok_explicitsnow) THEN
1235       CALL explicitsnow_finalize ( kjit,     kjpindex, rest_id,    snowrho,   &
1236                                    snowtemp, snowdz,   snowheat,   snowgrain)
1237    END IF
1238
1239  END SUBROUTINE hydrol_finalize
1240
1241
1242!! ================================================================================================================================
1243!! SUBROUTINE   : hydrol_init
1244!!
1245!>\BRIEF        Initializations and memory allocation   
1246!!
1247!! DESCRIPTION  :
1248!! - 1 Some initializations
1249!! - 2 make dynamic allocation with good dimension
1250!! - 2.1 array allocation for soil textur
1251!! - 2.2 Soil texture choice
1252!! - 3 Other array allocation
1253!! - 4 Open restart input file and read data for HYDROLOGIC process
1254!! - 5 get restart values if none were found in the restart file
1255!! - 6 Vegetation array     
1256!! - 7 set humrelv from us
1257!!
1258!! RECENT CHANGE(S) : None
1259!!
1260!! MAIN OUTPUT VARIABLE(S) :
1261!!
1262!! REFERENCE(S) :
1263!!
1264!! FLOWCHART    : None
1265!! \n
1266!_ ================================================================================================================================
1267!!_ hydrol_init
1268
1269  SUBROUTINE hydrol_init(kjit, kjpindex, index, rest_id, veget_max, soiltile, &
1270       humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
1271       snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, &
1272       drysoil_frac, evap_bare_lim)
1273   
1274
1275    !! 0. Variable and parameter declaration
1276
1277    !! 0.1 Input variables
1278
1279    INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number
1280    INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
1281    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: index              !! Indeces of the points on the map
1282    INTEGER(i_std), INTENT (in)                         :: rest_id            !! _Restart_ file identifier
1283    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max          !! Carte de vegetation max
1284    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in)  :: soiltile           !! Fraction of each soil tile within vegtot (0-1, unitless)
1285
1286    !! 0.2 Output variables
1287
1288    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: humrel             !! Stress hydrique, relative humidity
1289    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: vegstress          !! Veg. moisture stress (only for vegetation growth)
1290    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: snow               !! Snow mass [Kg/m^2]
1291    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: snow_age           !! Snow age
1292    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio       !! Snow on ice, lakes, ...
1293    REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age   !! Snow age on ice, lakes, ...
1294    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)    :: qsintveg         !! Water on vegetation due to interception
1295    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowdz           !! Snow depth
1296    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowgrain        !! Snow grain size
1297    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowheat         !! Snow heat content
1298    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowtemp         !! Snow temperature
1299    REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out)    :: snowrho          !! Snow density
1300    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: drysoil_frac     !! function of litter wetness
1301    REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: evap_bare_lim
1302
1303    !! 0.4 Local variables
1304
1305    INTEGER(i_std)                                     :: ier                   !! Error code
1306    INTEGER(i_std)                                     :: ji                    !! Index of land grid cells (1)
1307    INTEGER(i_std)                                     :: jv                    !! Index of PFTs (1)
1308    INTEGER(i_std)                                     :: jst                   !! Index of soil tiles (1)
1309    INTEGER(i_std)                                     :: jsl                   !! Index of soil layers (1)
1310    INTEGER(i_std)                                     :: jsc                   !! Index of soil texture (1)
1311    INTEGER(i_std), PARAMETER                          :: error_level = 3       !! Error level for consistency check
1312                                                                                !! Switch to 2 tu turn fatal errors into warnings 
1313    REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: free_drain_max        !! Temporary var for initialization of free_drain_coef
1314    REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: zwt_default           !! Temporary variable for initialization of zwt_force
1315    LOGICAL                                            :: zforce                !! To test if we force the WT in any of the soiltiles
1316
1317!_ ================================================================================================================================
1318
1319    !! 1 Some initializations
1320    !
1321    !Config Key   = DO_PONDS
1322    !Config Desc  = Should we include ponds
1323    !Config Def   = n
1324    !Config If    = HYDROL_CWRR
1325    !Config Help  = This parameters allows the user to ask the model
1326    !Config         to take into account the ponds and return
1327    !Config         the water into the soil moisture. If this is
1328    !Config         activated, then there is no reinfiltration
1329    !Config         computed inside the hydrol module.
1330    !Config Units = [FLAG]
1331    !
1332    doponds = .FALSE.
1333    CALL getin_p('DO_PONDS', doponds)
1334
1335    !Config Key   = FROZ_FRAC_CORR
1336    !Config Desc  = Coefficient for the frozen fraction correction
1337    !Config Def   = 1.0
1338    !Config If    = HYDROL_CWRR and OK_FREEZE
1339    !Config Help  =
1340    !Config Units = [-]
1341    froz_frac_corr = 1.0
1342    CALL getin_p("FROZ_FRAC_CORR", froz_frac_corr)
1343
1344    !Config Key   = MAX_FROZ_HYDRO
1345    !Config Desc  = Coefficient for the frozen fraction correction
1346    !Config Def   = 1.0
1347    !Config If    = HYDROL_CWRR and OK_FREEZE
1348    !Config Help  =
1349    !Config Units = [-]
1350    max_froz_hydro = 1.0
1351    CALL getin_p("MAX_FROZ_HYDRO", max_froz_hydro)
1352
1353    !Config Key   = SMTOT_CORR
1354    !Config Desc  = Coefficient for the frozen fraction correction
1355    !Config Def   = 2.0
1356    !Config If    = HYDROL_CWRR and OK_FREEZE
1357    !Config Help  =
1358    !Config Units = [-]
1359    smtot_corr = 2.0
1360    CALL getin_p("SMTOT_CORR", smtot_corr)
1361
1362    !Config Key   = DO_RSOIL
1363    !Config Desc  = Should we reduce soil evaporation with a soil resistance
1364    !Config Def   = n
1365    !Config If    = HYDROL_CWRR
1366    !Config Help  = This parameters allows the user to ask the model
1367    !Config         to calculate a soil resistance to reduce the soil evaporation
1368    !Config Units = [FLAG]
1369    !
1370    do_rsoil = .FALSE.
1371    CALL getin_p('DO_RSOIL', do_rsoil) 
1372
1373    !Config Key   = OK_DYNROOT
1374    !Config Desc  = Calculate dynamic root profile to optimize soil moisture usage 
1375    !Config Def   = n
1376    !Config If    = HYDROL_CWRR
1377    !Config Help  =
1378    !Config Units = [FLAG]
1379    ok_dynroot = .FALSE.
1380    CALL getin_p('OK_DYNROOT',ok_dynroot)
1381
1382    !! 2 make dynamic allocation with good dimension
1383
1384    !! 2.1 array allocation for soil texture
1385
1386    ALLOCATE (nvan(nscm),stat=ier)
1387    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan','','')
1388
1389    ALLOCATE (avan(nscm),stat=ier)
1390    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan','','')
1391
1392    ALLOCATE (mcr(nscm),stat=ier)
1393    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcr','','')
1394
1395    ALLOCATE (mcs(nscm),stat=ier)
1396    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcs','','')
1397
1398    ALLOCATE (ks(nscm),stat=ier)
1399    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ks','','')
1400
1401    ALLOCATE (pcent(nscm),stat=ier)
1402    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','')
1403
1404    ALLOCATE (mcfc(nscm),stat=ier)
1405    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcfc','','')
1406
1407    ALLOCATE (mcw(nscm),stat=ier)
1408    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcw','','')
1409   
1410    ALLOCATE (mc_awet(nscm),stat=ier)
1411    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','')
1412
1413    ALLOCATE (mc_adry(nscm),stat=ier)
1414    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','')
1415       
1416    !!__2.2 Soil texture choose
1417
1418    SELECTCASE (nscm)
1419    CASE (3)
1420             
1421       nvan(:) = nvan_fao(:)       
1422       avan(:) = avan_fao(:)
1423       mcr(:) = mcr_fao(:)
1424       mcs(:) = mcs_fao(:)
1425       ks(:) = ks_fao(:)
1426       pcent(:) = pcent_fao(:)
1427       mcfc(:) = mcf_fao(:)
1428       mcw(:) = mcw_fao(:)
1429       mc_awet(:) = mc_awet_fao(:)
1430       mc_adry(:) = mc_adry_fao(:)
1431    CASE (12)
1432       
1433       nvan(:) = nvan_usda(:)
1434       avan(:) = avan_usda(:)
1435       mcr(:) = mcr_usda(:)
1436       mcs(:) = mcs_usda(:)
1437       ks(:) = ks_usda(:)
1438       pcent(:) = pcent_usda(:)
1439       mcfc(:) = mcf_usda(:)
1440       mcw(:) = mcw_usda(:)
1441       mc_awet(:) = mc_awet_usda(:)
1442       mc_adry(:) = mc_adry_usda(:)
1443       
1444    CASE DEFAULT
1445       WRITE (numout,*) 'Unsupported soil type classification. Choose between zobler and usda according to the map'
1446       CALL ipslerr_p(3,'hydrol_init','Unsupported soil type classification. ',&
1447            'Choose between zobler and usda according to the map','')
1448    ENDSELECT
1449
1450
1451    !! 2.3 Read in the run.def the parameters values defined by the user
1452
1453    !Config Key   = CWRR_N_VANGENUCHTEN
1454    !Config Desc  = Van genuchten coefficient n
1455    !Config If    = HYDROL_CWRR
1456    !Config Def   = 1.89, 1.56, 1.31
1457    !Config Help  = This parameter will be constant over the entire
1458    !Config         simulated domain, thus independent from soil
1459    !Config         texture.   
1460    !Config Units = [-]
1461    CALL getin_p("CWRR_N_VANGENUCHTEN",nvan)
1462
1463    !! Check parameter value (correct range)
1464    IF ( ANY(nvan(:) <= zero) ) THEN
1465       CALL ipslerr_p(error_level, "hydrol_init.", &
1466            &     "Wrong parameter value for CWRR_N_VANGENUCHTEN.", &
1467            &     "This parameter should be positive. ", &
1468            &     "Please, check parameter value in run.def. ")
1469    END IF
1470
1471
1472    !Config Key   = CWRR_A_VANGENUCHTEN
1473    !Config Desc  = Van genuchten coefficient a
1474    !Config If    = HYDROL_CWRR
1475    !Config Def   = 0.0075, 0.0036, 0.0019
1476    !Config Help  = This parameter will be constant over the entire
1477    !Config         simulated domain, thus independent from soil
1478    !Config         texture.   
1479    !Config Units = [1/mm] 
1480    CALL getin_p("CWRR_A_VANGENUCHTEN",avan)
1481
1482    !! Check parameter value (correct range)
1483    IF ( ANY(avan(:) <= zero) ) THEN
1484       CALL ipslerr_p(error_level, "hydrol_init.", &
1485            &     "Wrong parameter value for CWRR_A_VANGENUCHTEN.", &
1486            &     "This parameter should be positive. ", &
1487            &     "Please, check parameter value in run.def. ")
1488    END IF
1489
1490
1491    !Config Key   = VWC_RESIDUAL
1492    !Config Desc  = Residual soil water content
1493    !Config If    = HYDROL_CWRR
1494    !Config Def   = 0.065, 0.078, 0.095
1495    !Config Help  = This parameter will be constant over the entire
1496    !Config         simulated domain, thus independent from soil
1497    !Config         texture.   
1498    !Config Units = [m3/m3] 
1499    CALL getin_p("VWC_RESIDUAL",mcr)
1500
1501    !! Check parameter value (correct range)
1502    IF ( ANY(mcr(:) < zero) .OR. ANY(mcr(:) > 1.)  ) THEN
1503       CALL ipslerr_p(error_level, "hydrol_init.", &
1504            &     "Wrong parameter value for VWC_RESIDUAL.", &
1505            &     "This parameter is ranged between 0 and 1. ", &
1506            &     "Please, check parameter value in run.def. ")
1507    END IF
1508
1509   
1510    !Config Key   = VWC_SAT
1511    !Config Desc  = Saturated soil water content
1512    !Config If    = HYDROL_CWRR
1513    !Config Def   = 0.41, 0.43, 0.41
1514    !Config Help  = This parameter will be constant over the entire
1515    !Config         simulated domain, thus independent from soil
1516    !Config         texture.   
1517    !Config Units = [m3/m3] 
1518    CALL getin_p("VWC_SAT",mcs)
1519
1520    !! Check parameter value (correct range)
1521    IF ( ANY(mcs(:) < zero) .OR. ANY(mcs(:) > 1.) .OR. ANY(mcs(:) <= mcr(:)) ) THEN
1522       CALL ipslerr_p(error_level, "hydrol_init.", &
1523            &     "Wrong parameter value for VWC_SAT.", &
1524            &     "This parameter should be greater than VWC_RESIDUAL and less than 1. ", &
1525            &     "Please, check parameter value in run.def. ")
1526    END IF
1527
1528
1529    !Config Key   = CWRR_KS
1530    !Config Desc  = Hydraulic conductivity Saturation
1531    !Config If    = HYDROL_CWRR
1532    !Config Def   = 1060.8, 249.6, 62.4
1533    !Config Help  = This parameter will be constant over the entire
1534    !Config         simulated domain, thus independent from soil
1535    !Config         texture.   
1536    !Config Units = [mm/d]   
1537    CALL getin_p("CWRR_KS",ks)
1538
1539    !! Check parameter value (correct range)
1540    IF ( ANY(ks(:) <= zero) ) THEN
1541       CALL ipslerr_p(error_level, "hydrol_init.", &
1542            &     "Wrong parameter value for CWRR_KS.", &
1543            &     "This parameter should be positive. ", &
1544            &     "Please, check parameter value in run.def. ")
1545    END IF
1546
1547
1548    !Config Key   = WETNESS_TRANSPIR_MAX
1549    !Config Desc  = Soil moisture above which transpir is max
1550    !Config If    = HYDROL_CWRR
1551    !Config Def   = 0.5, 0.5, 0.5
1552    !Config Help  = This parameter is independent from soil texture for
1553    !Config         the time being.
1554    !Config Units = [-]   
1555    CALL getin_p("WETNESS_TRANSPIR_MAX",pcent)
1556
1557    !! Check parameter value (correct range)
1558    IF ( ANY(pcent(:) <= zero) .OR. ANY(pcent(:) > 1.) ) THEN
1559       CALL ipslerr_p(error_level, "hydrol_init.", &
1560            &     "Wrong parameter value for WETNESS_TRANSPIR_MAX.", &
1561            &     "This parameter should be positive and less or equals than 1. ", &
1562            &     "Please, check parameter value in run.def. ")
1563    END IF
1564
1565
1566    !Config Key   = VWC_FC
1567    !Config Desc  = Volumetric water content field capacity
1568    !Config If    = HYDROL_CWRR
1569    !Config Def   = 0.32, 0.32, 0.32
1570    !Config Help  = This parameter is independent from soil texture for
1571    !Config         the time being.
1572    !Config Units = [m3/m3]   
1573    CALL getin_p("VWC_FC",mcfc)
1574
1575    !! Check parameter value (correct range)
1576    IF ( ANY(mcfc(:) > mcs(:)) ) THEN
1577       CALL ipslerr_p(error_level, "hydrol_init.", &
1578            &     "Wrong parameter value for VWC_FC.", &
1579            &     "This parameter should be less than VWC_SAT. ", &
1580            &     "Please, check parameter value in run.def. ")
1581    END IF
1582
1583
1584    !Config Key   = VWC_WP
1585    !Config Desc  = Volumetric water content Wilting pt
1586    !Config If    = HYDROL_CWRR
1587    !Config Def   = 0.10, 0.10, 0.10
1588    !Config Help  = This parameter is independent from soil texture for
1589    !Config         the time being.
1590    !Config Units = [m3/m3]   
1591    CALL getin_p("VWC_WP",mcw)
1592
1593    !! Check parameter value (correct range)
1594    IF ( ANY(mcw(:) > mcfc(:)) .OR. ANY(mcw(:) < mcr(:)) ) THEN
1595       CALL ipslerr_p(error_level, "hydrol_init.", &
1596            &     "Wrong parameter value for VWC_WP.", &
1597            &     "This parameter should be greater or equal than VWC_RESIDUAL and less or equal than VWC_SAT.", &
1598            &     "Please, check parameter value in run.def. ")
1599    END IF
1600
1601
1602    !Config Key   = VWC_MIN_FOR_WET_ALB
1603    !Config Desc  = Vol. wat. cont. above which albedo is cst
1604    !Config If    = HYDROL_CWRR
1605    !Config Def   = 0.25, 0.25, 0.25
1606    !Config Help  = This parameter is independent from soil texture for
1607    !Config         the time being.
1608    !Config Units = [m3/m3] 
1609    CALL getin_p("VWC_MIN_FOR_WET_ALB",mc_awet)
1610
1611    !! Check parameter value (correct range)
1612    IF ( ANY(mc_awet(:) < 0) ) THEN
1613       CALL ipslerr_p(error_level, "hydrol_init.", &
1614            &     "Wrong parameter value for VWC_MIN_FOR_WET_ALB.", &
1615            &     "This parameter should be positive. ", &
1616            &     "Please, check parameter value in run.def. ")
1617    END IF
1618
1619
1620    !Config Key   = VWC_MAX_FOR_DRY_ALB
1621    !Config Desc  = Vol. wat. cont. below which albedo is cst
1622    !Config If    = HYDROL_CWRR
1623    !Config Def   = 0.1, 0.1, 0.1
1624    !Config Help  = This parameter is independent from soil texture for
1625    !Config         the time being.
1626    !Config Units = [m3/m3]   
1627    CALL getin_p("VWC_MAX_FOR_DRY_ALB",mc_adry)
1628
1629    !! Check parameter value (correct range)
1630    IF ( ANY(mc_adry(:) < 0) .OR. ANY(mc_adry(:) > mc_awet(:)) ) THEN
1631       CALL ipslerr_p(error_level, "hydrol_init.", &
1632            &     "Wrong parameter value for VWC_MAX_FOR_DRY_ALB.", &
1633            &     "This parameter should be positive and not greater than VWC_MIN_FOR_WET_ALB.", &
1634            &     "Please, check parameter value in run.def. ")
1635    END IF
1636
1637
1638    !! 3 Other array allocation
1639
1640
1641    ALLOCATE (mask_veget(kjpindex,nvm),stat=ier)
1642    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_veget','','')
1643
1644    ALLOCATE (mask_soiltile(kjpindex,nstm),stat=ier)
1645    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_soiltile','','')
1646
1647    ALLOCATE (humrelv(kjpindex,nvm,nstm),stat=ier)
1648    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable humrelv','','')
1649
1650    ALLOCATE (vegstressv(kjpindex,nvm,nstm),stat=ier) 
1651    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegstressv','','')
1652
1653    ALLOCATE (us(kjpindex,nvm,nstm,nslm),stat=ier) 
1654    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable us','','')
1655
1656    ALLOCATE (precisol(kjpindex,nvm),stat=ier) 
1657    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable precisol','','')
1658
1659    ALLOCATE (throughfall(kjpindex,nvm),stat=ier) 
1660    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable throughfall','','')
1661
1662    ALLOCATE (precisol_ns(kjpindex,nstm),stat=ier) 
1663    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable precisol_nc','','')
1664
1665    ALLOCATE (free_drain_coef(kjpindex,nstm),stat=ier) 
1666    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable free_drain_coef','','')
1667
1668    ALLOCATE (zwt_force(kjpindex,nstm),stat=ier) 
1669    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zwt_force','','')
1670
1671    ALLOCATE (frac_bare_ns(kjpindex,nstm),stat=ier) 
1672    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable frac_bare_ns','','')
1673
1674    ALLOCATE (water2infilt(kjpindex,nstm),stat=ier)
1675    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable water2infilt','','')
1676
1677    ALLOCATE (ae_ns(kjpindex,nstm),stat=ier) 
1678    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ae_ns','','')
1679
1680    ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier) 
1681    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable evap_bare_lim_ns','','')
1682
1683    ALLOCATE (rootsink(kjpindex,nslm,nstm),stat=ier) 
1684    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable rootsink','','')
1685
1686    ALLOCATE (subsnowveg(kjpindex),stat=ier) 
1687    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsnowveg','','')
1688
1689    ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier) 
1690    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsnownobio','','')
1691
1692    ALLOCATE (icemelt(kjpindex),stat=ier) 
1693    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable icemelt','','')
1694
1695    ALLOCATE (subsinksoil(kjpindex),stat=ier) 
1696    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable subsinksoil','','')
1697
1698    ALLOCATE (mx_eau_var(kjpindex),stat=ier)
1699    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mx_eau_var','','')
1700
1701    ALLOCATE (vegtot(kjpindex),stat=ier) 
1702    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegtot','','')
1703
1704    ALLOCATE (vegtot_old(kjpindex),stat=ier) 
1705    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegtot_old','','')
1706
1707    ALLOCATE (resdist(kjpindex,nstm),stat=ier)
1708    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable resdist','','')
1709
1710    ALLOCATE (humtot(kjpindex),stat=ier)
1711    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable humtot','','')
1712
1713    ALLOCATE (resolv(kjpindex),stat=ier) 
1714    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable resolv','','')
1715
1716    ALLOCATE (k(kjpindex,nslm),stat=ier) 
1717    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k','','')
1718
1719    IF (ok_freeze_cwrr) THEN
1720       ALLOCATE (kk_moy(kjpindex,nslm),stat=ier) 
1721       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kk_moy','','')
1722       kk_moy(:,:) = 276.48
1723
1724       ALLOCATE (kk(kjpindex,nslm,nstm),stat=ier) 
1725       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kk','','')
1726       kk(:,:,:) = 276.48
1727    ENDIF
1728
1729    ALLOCATE (a(kjpindex,nslm),stat=ier) 
1730    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a','','')
1731
1732    ALLOCATE (b(kjpindex,nslm),stat=ier)
1733    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b','','')
1734
1735    ALLOCATE (d(kjpindex,nslm),stat=ier)
1736    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d','','')
1737
1738    ALLOCATE (e(kjpindex,nslm),stat=ier) 
1739    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable e','','')
1740
1741    ALLOCATE (f(kjpindex,nslm),stat=ier) 
1742    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable f','','')
1743
1744    ALLOCATE (g1(kjpindex,nslm),stat=ier) 
1745    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable g1','','')
1746
1747    ALLOCATE (ep(kjpindex,nslm),stat=ier)
1748    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ep','','')
1749
1750    ALLOCATE (fp(kjpindex,nslm),stat=ier)
1751    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable fp','','')
1752
1753    ALLOCATE (gp(kjpindex,nslm),stat=ier)
1754    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable gp','','')
1755
1756    ALLOCATE (rhs(kjpindex,nslm),stat=ier)
1757    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable rhs','','')
1758
1759    ALLOCATE (srhs(kjpindex,nslm),stat=ier)
1760    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable srhs','','')
1761
1762    ALLOCATE (tmc(kjpindex,nstm),stat=ier)
1763    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc','','')
1764
1765    ALLOCATE (tmcs(kjpindex,nstm),stat=ier)
1766    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcs','','')
1767
1768    ALLOCATE (tmcr(kjpindex,nstm),stat=ier)
1769    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcr','','')
1770
1771    ALLOCATE (tmcfc(kjpindex,nstm),stat=ier)
1772    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcfc','','')
1773
1774    ALLOCATE (tmcw(kjpindex,nstm),stat=ier)
1775    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmcw','','')
1776
1777    ALLOCATE (tmc_litter(kjpindex,nstm),stat=ier)
1778    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter','','')
1779
1780    ALLOCATE (tmc_litt_mea(kjpindex),stat=ier)
1781    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_mea','','')
1782
1783    ALLOCATE (tmc_litter_res(kjpindex,nstm),stat=ier)
1784    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_res','','')
1785
1786    ALLOCATE (tmc_litter_wilt(kjpindex,nstm),stat=ier)
1787    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_wilt','','')
1788
1789    ALLOCATE (tmc_litter_field(kjpindex,nstm),stat=ier)
1790    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_field','','')
1791
1792    ALLOCATE (tmc_litter_sat(kjpindex,nstm),stat=ier)
1793    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_sat','','')
1794
1795    ALLOCATE (tmc_litter_awet(kjpindex,nstm),stat=ier)
1796    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_awet','','')
1797
1798    ALLOCATE (tmc_litter_adry(kjpindex,nstm),stat=ier)
1799    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litter_adry','','')
1800
1801    ALLOCATE (tmc_litt_wet_mea(kjpindex),stat=ier)
1802    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_wet_mea','','')
1803
1804    ALLOCATE (tmc_litt_dry_mea(kjpindex),stat=ier)
1805    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmc_litt_dry_mea','','')
1806
1807    ALLOCATE (v1(kjpindex,nstm),stat=ier)
1808    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable v1','','')
1809
1810    ALLOCATE (ru_ns(kjpindex,nstm),stat=ier)
1811    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ru_ns','','')
1812    ru_ns(:,:) = zero
1813
1814    ALLOCATE (dr_ns(kjpindex,nstm),stat=ier)
1815    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dr_ns','','')
1816    dr_ns(:,:) = zero
1817
1818    ALLOCATE (tr_ns(kjpindex,nstm),stat=ier)
1819    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tr_ns','','')
1820
1821    ALLOCATE (vegetmax_soil(kjpindex,nvm,nstm),stat=ier)
1822    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable vegetmax_soil','','')
1823
1824    ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier)
1825    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc','','')
1826
1827
1828    ! Variables for nudging of soil moisture
1829    IF (ok_nudge_mc) THEN
1830       ALLOCATE (mc_read_prev(kjpindex,nslm,nstm),stat=ier)
1831       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_read_prev','','')
1832       ALLOCATE (mc_read_next(kjpindex,nslm,nstm),stat=ier)
1833       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_read_next','','')
1834       ALLOCATE (mask_mc_interp(kjpindex,nslm,nstm),stat=ier)
1835       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_mc_interp','','')
1836    END IF
1837
1838    ! Variables for nudging of snow variables
1839    IF (ok_nudge_snow) THEN
1840       ALLOCATE (snowdz_read_prev(kjpindex,nsnow),stat=ier)
1841       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowdz_read_prev','','')
1842       ALLOCATE (snowdz_read_next(kjpindex,nsnow),stat=ier)
1843       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowdz_read_next','','')
1844       
1845       ALLOCATE (snowrho_read_prev(kjpindex,nsnow),stat=ier)
1846       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowrho_read_prev','','')
1847       ALLOCATE (snowrho_read_next(kjpindex,nsnow),stat=ier)
1848       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowrho_read_next','','')
1849       
1850       ALLOCATE (snowtemp_read_prev(kjpindex,nsnow),stat=ier)
1851       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowtemp_read_prev','','')
1852       ALLOCATE (snowtemp_read_next(kjpindex,nsnow),stat=ier)
1853       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snowtemp_read_next','','')
1854       
1855       ALLOCATE (mask_snow_interp(kjpindex,nsnow),stat=ier)
1856       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mask_snow_interp','','')
1857    END IF
1858
1859    ALLOCATE (mcl(kjpindex, nslm, nstm),stat=ier)
1860    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcl','','')
1861
1862    IF (ok_freeze_cwrr) THEN
1863       ALLOCATE (profil_froz_hydro(kjpindex, nslm),stat=ier)
1864       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable profil_froz_hydrol','','')
1865       profil_froz_hydro(:,:) = zero
1866       
1867       ALLOCATE (temp_hydro(kjpindex, nslm),stat=ier)
1868       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable temp_hydro','','')
1869       temp_hydro(:,:) = 280.
1870    ENDIF
1871   
1872    ALLOCATE (profil_froz_hydro_ns(kjpindex, nslm, nstm),stat=ier)
1873    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable profil_froz_hydro_ns','','')
1874    profil_froz_hydro_ns(:,:,:) = zero
1875   
1876    ALLOCATE (soilmoist(kjpindex,nslm),stat=ier)
1877    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soilmoist','','')
1878
1879    ALLOCATE (soilmoist_liquid(kjpindex,nslm),stat=ier)
1880    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soilmoist_liquid','','')
1881
1882    ALLOCATE (soil_wet_ns(kjpindex,nslm,nstm),stat=ier)
1883    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_ns','','')
1884
1885    ALLOCATE (soil_wet_litter(kjpindex,nstm),stat=ier)
1886    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_litter','','')
1887
1888    ALLOCATE (qflux(kjpindex,nslm,nstm),stat=ier) 
1889    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable qflux','','')
1890
1891    ALLOCATE (tmat(kjpindex,nslm,3),stat=ier)
1892    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tmat','','')
1893
1894    ALLOCATE (stmat(kjpindex,nslm,3),stat=ier)
1895    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable stmat','','')
1896
1897    ALLOCATE (nroot(kjpindex,nvm, nslm),stat=ier)
1898    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nroot','','')
1899
1900    ALLOCATE (kfact_root(kjpindex, nslm, nstm), stat=ier)
1901    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','')
1902
1903    ALLOCATE (kfact(nslm, nscm),stat=ier)
1904    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','')
1905
1906    ALLOCATE (zz(nslm),stat=ier)
1907    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zz','','')
1908
1909    ALLOCATE (dz(nslm),stat=ier)
1910    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dz','','')
1911   
1912    ALLOCATE (dh(nslm),stat=ier)
1913    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','')
1914
1915    ALLOCATE (mc_lin(imin:imax, nscm),stat=ier)
1916    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','')
1917
1918    ALLOCATE (k_lin(imin:imax, nslm, nscm),stat=ier)
1919    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','')
1920
1921    ALLOCATE (d_lin(imin:imax, nslm, nscm),stat=ier)
1922    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','')
1923
1924    ALLOCATE (a_lin(imin:imax, nslm, nscm),stat=ier)
1925    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','')
1926
1927    ALLOCATE (b_lin(imin:imax, nslm, nscm),stat=ier)
1928    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','')
1929
1930    ALLOCATE (undermcr(kjpindex),stat=ier)
1931    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable undermcr','','')
1932
1933    ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
1934    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watveg_beg','','')
1935   
1936    ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
1937    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watvag_end','','')
1938   
1939    ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
1940    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watsoil_beg','','')
1941   
1942    ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
1943    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable tot_watsoil_end','','')
1944   
1945    ALLOCATE (delsoilmoist(kjpindex),stat=ier)
1946    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delsoilmoist','','')
1947   
1948    ALLOCATE (delintercept(kjpindex),stat=ier)
1949    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delintercept','','')
1950   
1951    ALLOCATE (delswe(kjpindex),stat=ier)
1952    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable delswe','','')
1953   
1954    ALLOCATE (snow_beg(kjpindex),stat=ier)
1955    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snow_beg','','')
1956   
1957    ALLOCATE (snow_end(kjpindex),stat=ier)
1958    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable snow_end','','')
1959   
1960    !! 4 Open restart input file and read data for HYDROLOGIC process
1961       IF (printlev>=3) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
1962
1963       IF (is_root_prc) CALL ioconf_setatt_p('UNITS', '-')
1964       !
1965       DO jst=1,nstm
1966          ! var_name= "mc_1" ... "mc_3"
1967           WRITE (var_name,"('moistc_',I1)") jst
1968           IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
1969           CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc(:,:,jst), "gather", nbp_glo, index_g)
1970       END DO
1971
1972       IF (ok_nudge_mc) THEN
1973          DO jst=1,nstm
1974             WRITE (var_name,"('mc_read_next_',I1)") jst
1975             IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME','Soil moisture read from nudging file')
1976             CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc_read_next(:,:,jst), &
1977                  "gather", nbp_glo, index_g)
1978          END DO
1979       END IF
1980
1981       IF (ok_nudge_snow) THEN
1982          IF (is_root_prc) THEN
1983             CALL ioconf_setatt_p('UNITS', 'm')
1984             CALL ioconf_setatt_p('LONG_NAME','Snow layer thickness read from nudging file')
1985          ENDIF
1986          CALL restget_p (rest_id, 'snowdz_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowdz_read_next, &
1987               "gather", nbp_glo, index_g)
1988
1989          IF (is_root_prc) THEN
1990             CALL ioconf_setatt_p('UNITS', 'kg/m^3')
1991             CALL ioconf_setatt_p('LONG_NAME','Snow density profile read from nudging file')
1992          ENDIF
1993          CALL restget_p (rest_id, 'snowrho_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowrho_read_next, &
1994               "gather", nbp_glo, index_g)
1995
1996          IF (is_root_prc) THEN
1997             CALL ioconf_setatt_p('UNITS', 'K')
1998             CALL ioconf_setatt_p('LONG_NAME','Snow temperature read from nudging file')
1999          ENDIF
2000          CALL restget_p (rest_id, 'snowtemp_read_next', nbp_glo, nsnow, 1, kjit, .TRUE., snowtemp_read_next, &
2001               "gather", nbp_glo, index_g)
2002       END IF
2003     
2004       DO jst=1,nstm
2005          ! var_name= "mcl_1" ... "mcl_3"
2006           WRITE (var_name,"('moistcl_',I1)") jst
2007           IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
2008           CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mcl(:,:,jst), "gather", nbp_glo, index_g)
2009       END DO
2010
2011       IF (is_root_prc) CALL ioconf_setatt_p('UNITS', '-')
2012       DO jst=1,nstm
2013          DO jsl=1,nslm
2014             ! var_name= "us_1_01" ... "us_3_11"
2015             WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
2016             IF (is_root_prc) CALL ioconf_setatt_p('LONG_NAME',var_name)
2017             CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., us(:,:,jst,jsl), "gather", nbp_glo, index_g)
2018          END DO
2019       END DO
2020       !
2021       var_name= 'free_drain_coef'
2022       IF (is_root_prc) THEN
2023          CALL ioconf_setatt_p('UNITS', '-')
2024          CALL ioconf_setatt_p('LONG_NAME','Coefficient for free drainage at bottom of soil')
2025       ENDIF
2026       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., free_drain_coef, "gather", nbp_glo, index_g)
2027       !
2028       var_name= 'zwt_force'
2029       IF (is_root_prc) THEN
2030          CALL ioconf_setatt_p('UNITS', 'm')
2031          CALL ioconf_setatt_p('LONG_NAME','Prescribed water table depth')
2032       ENDIF
2033       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., zwt_force, "gather", nbp_glo, index_g)
2034       !
2035       var_name= 'water2infilt'
2036       IF (is_root_prc) THEN
2037          CALL ioconf_setatt_p('UNITS', '-')
2038          CALL ioconf_setatt_p('LONG_NAME','Remaining water to be infiltrated on top of the soil')
2039       ENDIF
2040       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., water2infilt, "gather", nbp_glo, index_g)
2041       !
2042       var_name= 'ae_ns'
2043       IF (is_root_prc) THEN
2044          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
2045          CALL ioconf_setatt_p('LONG_NAME','Bare soil evap on each soil type')
2046       ENDIF
2047       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., ae_ns, "gather", nbp_glo, index_g)
2048       !
2049       var_name= 'snow'       
2050       IF (is_root_prc) THEN
2051          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
2052          CALL ioconf_setatt_p('LONG_NAME','Snow mass')
2053       ENDIF
2054       CALL restget_p (rest_id, var_name, nbp_glo, 1  , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
2055       !
2056       var_name= 'snow_age'
2057       IF (is_root_prc) THEN
2058          CALL ioconf_setatt_p('UNITS', 'd')
2059          CALL ioconf_setatt_p('LONG_NAME','Snow age')
2060       ENDIF
2061       CALL restget_p (rest_id, var_name, nbp_glo, 1  , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
2062       !
2063       var_name= 'snow_nobio'
2064       IF (is_root_prc) THEN
2065          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
2066          CALL ioconf_setatt_p('LONG_NAME','Snow on other surface types')
2067       ENDIF
2068       CALL restget_p (rest_id, var_name, nbp_glo, nnobio  , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
2069       !
2070       var_name= 'snow_nobio_age'
2071       IF (is_root_prc) THEN
2072          CALL ioconf_setatt_p('UNITS', 'd')
2073          CALL ioconf_setatt_p('LONG_NAME','Snow age on other surface types')
2074       ENDIF
2075       CALL restget_p (rest_id, var_name, nbp_glo, nnobio  , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
2076       !
2077       var_name= 'qsintveg'
2078       IF (is_root_prc) THEN
2079          CALL ioconf_setatt_p('UNITS', 'kg/m^2')
2080          CALL ioconf_setatt_p('LONG_NAME','Intercepted moisture')
2081       ENDIF
2082       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
2083
2084       var_name= 'evap_bare_lim_ns'
2085       IF (is_root_prc) THEN
2086          CALL ioconf_setatt_p('UNITS', '?')
2087          CALL ioconf_setatt_p('LONG_NAME','?')
2088       ENDIF
2089       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., evap_bare_lim_ns, "gather", nbp_glo, index_g)
2090       CALL setvar_p (evap_bare_lim_ns, val_exp, 'NO_KEYWORD', 0.0)
2091
2092       var_name= 'resdist'
2093       IF (is_root_prc) THEN
2094          CALL ioconf_setatt_p('UNITS', '-')
2095          CALL ioconf_setatt_p('LONG_NAME','soiltile values from previous time-step')
2096       ENDIF
2097       CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
2098
2099       var_name= 'vegtot_old'
2100       IF (is_root_prc) THEN
2101          CALL ioconf_setatt_p('UNITS', '-')
2102          CALL ioconf_setatt_p('LONG_NAME','vegtot from previous time-step')
2103       ENDIF
2104       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_old, "gather", nbp_glo, index_g)       
2105       
2106       ! Read drysoil_frac. It will be initalized later in hydrol_var_init if the varaible is not find in restart file.
2107       IF (is_root_prc) THEN
2108          CALL ioconf_setatt_p('UNITS', '')
2109          CALL ioconf_setatt_p('LONG_NAME','Function of litter wetness')
2110       ENDIF
2111       CALL restget_p (rest_id, 'drysoil_frac', nbp_glo, 1  , 1, kjit, .TRUE., drysoil_frac, "gather", nbp_glo, index_g)
2112
2113
2114    !! 5 get restart values if none were found in the restart file
2115       !
2116       !Config Key   = HYDROL_MOISTURE_CONTENT
2117       !Config Desc  = Soil moisture on each soil tile and levels
2118       !Config If    = HYDROL_CWRR       
2119       !Config Def   = 0.3
2120       !Config Help  = The initial value of mc if its value is not found
2121       !Config         in the restart file. This should only be used if the model is
2122       !Config         started without a restart file.
2123       !Config Units = [m3/m3]
2124       !
2125       CALL setvar_p (mc, val_exp, 'HYDROL_MOISTURE_CONTENT', 0.3_r_std)
2126
2127       ! Initialize mcl as mc if it is not found in the restart file
2128       IF ( ALL(mcl(:,:,:)==val_exp) ) THEN
2129          mcl(:,:,:) = mc(:,:,:)
2130       END IF
2131
2132
2133       
2134       !Config Key   = US_INIT
2135       !Config Desc  = US_NVM_NSTM_NSLM
2136       !Config If    = HYDROL_CWRR       
2137       !Config Def   = 0.0
2138       !Config Help  = The initial value of us (relative moisture) if its value is not found
2139       !Config         in the restart file. This should only be used if the model is
2140       !Config         started without a restart file.
2141       !Config Units = [-]
2142       !
2143       DO jsl=1,nslm
2144          CALL setvar_p (us(:,:,:,jsl), val_exp, 'US_INIT', zero)
2145       ENDDO
2146       !
2147       !Config Key   = ZWT_FORCE
2148       !Config Desc  = Prescribed water depth, dimension nstm
2149       !Config If    = HYDROL_CWRR       
2150       !Config Def   = undef undef undef
2151       !Config Help  = The initial value of zwt_force if its value is not found
2152       !Config         in the restart file. undef corresponds to a case whith no forced WT.
2153       !Config         This should only be used if the model is started without a restart file.
2154       !Config Units = [m]
2155       
2156       ALLOCATE (zwt_default(nstm),stat=ier)
2157       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable zwt_default','','')
2158       zwt_default(:) = undef_sechiba
2159       CALL setvar_p (zwt_force, val_exp, 'ZWT_FORCE', zwt_default )
2160
2161       zforce = .FALSE.
2162       DO jst=1,nstm
2163          IF (zwt_force(1,jst) <= zmaxh) zforce = .TRUE. ! AD16*** check if OK with vertical_soil
2164       ENDDO
2165       !
2166       !Config Key   = FREE_DRAIN_COEF
2167       !Config Desc  = Coefficient for free drainage at bottom, dimension nstm
2168       !Config If    = HYDROL_CWRR       
2169       !Config Def   = 1.0 1.0 1.0
2170       !Config Help  = The initial value of free drainage coefficient if its value is not found
2171       !Config         in the restart file. This should only be used if the model is
2172       !Config         started without a restart file.
2173       !Config Units = [-]
2174             
2175       ALLOCATE (free_drain_max(nstm),stat=ier)
2176       IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable free_drain_max','','')
2177       free_drain_max(:)=1.0
2178
2179       CALL setvar_p (free_drain_coef, val_exp, 'FREE_DRAIN_COEF', free_drain_max)
2180       IF (printlev>=2) WRITE (numout,*) ' hydrol_init => free_drain_coef = ',free_drain_coef(1,:)
2181       DEALLOCATE(free_drain_max)
2182
2183       !
2184       !Config Key   = WATER_TO_INFILT
2185       !Config Desc  = Water to be infiltrated on top of the soil
2186       !Config If    = HYDROL_CWRR   
2187       !Config Def   = 0.0
2188       !Config Help  = The initial value of free drainage if its value is not found
2189       !Config         in the restart file. This should only be used if the model is
2190       !Config         started without a restart file.
2191       !Config Units = [mm]
2192       !
2193       CALL setvar_p (water2infilt, val_exp, 'WATER_TO_INFILT', zero)
2194       !
2195       !Config Key   = EVAPNU_SOIL
2196       !Config Desc  = Bare soil evap on each soil if not found in restart
2197       !Config If    = HYDROL_CWRR 
2198       !Config Def   = 0.0
2199       !Config Help  = The initial value of bare soils evap if its value is not found
2200       !Config         in the restart file. This should only be used if the model is
2201       !Config         started without a restart file.
2202       !Config Units = [mm]
2203       !
2204       CALL setvar_p (ae_ns, val_exp, 'EVAPNU_SOIL', zero)
2205       !
2206       !Config Key  = HYDROL_SNOW
2207       !Config Desc  = Initial snow mass if not found in restart
2208       !Config If    = OK_SECHIBA
2209       !Config Def   = 0.0
2210       !Config Help  = The initial value of snow mass if its value is not found
2211       !Config         in the restart file. This should only be used if the model is
2212       !Config         started without a restart file.
2213       !Config Units =
2214       !
2215       CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', zero)
2216       !
2217       !Config Key   = HYDROL_SNOWAGE
2218       !Config Desc  = Initial snow age if not found in restart
2219       !Config If    = OK_SECHIBA
2220       !Config Def   = 0.0
2221       !Config Help  = The initial value of snow age if its value is not found
2222       !Config         in the restart file. This should only be used if the model is
2223       !Config         started without a restart file.
2224       !Config Units = ***
2225       !
2226       CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', zero)
2227       !
2228       !Config Key   = HYDROL_SNOW_NOBIO
2229       !Config Desc  = Initial snow amount on ice, lakes, etc. if not found in restart
2230       !Config If    = OK_SECHIBA
2231       !Config Def   = 0.0
2232       !Config Help  = The initial value of snow if its value is not found
2233       !Config         in the restart file. This should only be used if the model is
2234       !Config         started without a restart file.
2235       !Config Units = [mm]
2236       !
2237       CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero)
2238       !
2239       !Config Key   = HYDROL_SNOW_NOBIO_AGE
2240       !Config Desc  = Initial snow age on ice, lakes, etc. if not found in restart
2241       !Config If    = OK_SECHIBA
2242       !Config Def   = 0.0
2243       !Config Help  = The initial value of snow age if its value is not found
2244       !Config         in the restart file. This should only be used if the model is
2245       !Config         started without a restart file.
2246       !Config Units = ***
2247       !
2248       CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero)
2249       !
2250       !Config Key   = HYDROL_QSV
2251       !Config Desc  = Initial water on canopy if not found in restart
2252       !Config If    = OK_SECHIBA
2253       !Config Def   = 0.0
2254       !Config Help  = The initial value of moisture on canopy if its value
2255       !Config         is not found in the restart file. This should only be used if
2256       !Config         the model is started without a restart file.
2257       !Config Units = [mm]
2258       !
2259       CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', zero)
2260
2261    !! 6 Vegetation array     
2262       !
2263       ! If resdist is not in restart file, initialize with soiltile
2264       IF ( MINVAL(resdist) .EQ.  MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
2265          resdist(:,:) = soiltile(:,:)
2266       ENDIF
2267       
2268       !
2269       !  Remember that it is only frac_nobio + SUM(veget_max(,:)) that is equal to 1. Thus we need vegtot
2270       !
2271       IF ( ALL(vegtot_old(:) == val_exp) ) THEN
2272          ! vegtot_old was not found in restart file
2273          DO ji = 1, kjpindex
2274             vegtot_old(ji) = SUM(veget_max(ji,:))
2275          ENDDO
2276       ENDIF
2277       
2278       ! In the initialization phase, vegtot must take the value from previous time-step.
2279       ! This is because hydrol_main is done before veget_max is updated in the end of the time step.
2280       vegtot(:) = vegtot_old(:)
2281       
2282       !
2283       !
2284       ! compute the masks for veget
2285
2286       mask_veget(:,:) = 0
2287       mask_soiltile(:,:) = 0
2288
2289       DO jst=1,nstm
2290          DO ji = 1, kjpindex
2291             IF(soiltile(ji,jst) .GT. min_sechiba) THEN
2292                mask_soiltile(ji,jst) = 1
2293             ENDIF
2294          END DO
2295       ENDDO
2296         
2297       DO jv = 1, nvm
2298          DO ji = 1, kjpindex
2299             IF(veget_max(ji,jv) .GT. min_sechiba) THEN
2300                mask_veget(ji,jv) = 1
2301             ENDIF
2302          END DO
2303       END DO
2304
2305       humrelv(:,:,:) = SUM(us,dim=4)
2306
2307         
2308       !! 7a. Set vegstress
2309     
2310       var_name= 'vegstress'
2311       IF (is_root_prc) THEN
2312          CALL ioconf_setatt_p('UNITS', '-')
2313          CALL ioconf_setatt_p('LONG_NAME','Vegetation growth moisture stress')
2314       ENDIF
2315       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
2316
2317       vegstressv(:,:,:) = humrelv(:,:,:)
2318       ! Calculate vegstress if it is not found in restart file
2319       IF (ALL(vegstress(:,:)==val_exp)) THEN
2320          DO jv=1,nvm
2321             DO ji=1,kjpindex
2322                vegstress(ji,jv)=vegstress(ji,jv) + vegstressv(ji,jv,pref_soil_veg(jv))
2323             END DO
2324          END DO
2325       END IF
2326       !! 7b. Set humrel   
2327       ! Read humrel from restart file
2328       var_name= 'humrel'
2329       IF (is_root_prc) THEN
2330          CALL ioconf_setatt_p('UNITS', '')
2331          CALL ioconf_setatt_p('LONG_NAME','Relative humidity')
2332       ENDIF
2333       CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g)
2334
2335       ! Calculate humrel if it is not found in restart file
2336       IF (ALL(humrel(:,:)==val_exp)) THEN
2337          ! set humrel from humrelv, assuming equi-repartition for the first time step
2338          humrel(:,:) = zero
2339          DO jv=1,nvm
2340             DO ji=1,kjpindex
2341                humrel(ji,jv)=humrel(ji,jv) + humrelv(ji,jv,pref_soil_veg(jv))     
2342             END DO
2343          END DO
2344       END IF
2345
2346       ! Read evap_bare_lim from restart file
2347       var_name= 'evap_bare_lim'
2348       IF (is_root_prc) THEN
2349          CALL ioconf_setatt_p('UNITS', '')
2350          CALL ioconf_setatt_p('LONG_NAME','Limitation factor for bare soil evaporation')
2351       ENDIF
2352       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., evap_bare_lim, "gather", nbp_glo, index_g)
2353
2354       ! Calculate evap_bare_lim if it was not found in the restart file.
2355       IF ( ALL(evap_bare_lim(:) == val_exp) ) THEN
2356          DO ji = 1, kjpindex
2357             evap_bare_lim(ji) =  SUM(evap_bare_lim_ns(ji,:)*vegtot(ji)*soiltile(ji,:))
2358          ENDDO
2359       END IF
2360
2361
2362    ! Read from restart file       
2363    ! The variables tot_watsoil_beg, tot_watsoil_beg and snwo_beg will be initialized in the end of
2364    ! hydrol_initialize if they were not found in the restart file.
2365       
2366    var_name= 'tot_watveg_beg'
2367    IF (is_root_prc) THEN
2368       CALL ioconf_setatt_p('UNITS', '?')
2369       CALL ioconf_setatt_p('LONG_NAME','?')
2370    ENDIF
2371    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tot_watveg_beg, "gather", nbp_glo, index_g)
2372   
2373    var_name= 'tot_watsoil_beg'
2374    IF (is_root_prc) THEN
2375       CALL ioconf_setatt_p('UNITS', '?')
2376       CALL ioconf_setatt_p('LONG_NAME','?')
2377    ENDIF
2378    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., tot_watsoil_beg, "gather", nbp_glo, index_g)
2379   
2380    var_name= 'snow_beg'
2381    IF (is_root_prc) THEN
2382       CALL ioconf_setatt_p('UNITS', '?')
2383       CALL ioconf_setatt_p('LONG_NAME','?')
2384    ENDIF
2385    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., snow_beg, "gather", nbp_glo, index_g)
2386       
2387 
2388    ! Initialize variables for explictsnow module by reading restart file
2389    IF (ok_explicitsnow) THEN
2390       CALL explicitsnow_initialize( kjit,     kjpindex, rest_id,    snowrho,   &
2391                                     snowtemp, snowdz,   snowheat,   snowgrain)
2392    END IF
2393
2394
2395    ! Initialize soil moisture for nudging if not found in restart file
2396    IF (ok_nudge_mc) THEN
2397       IF ( ALL(mc_read_next(:,:,:)==val_exp) ) mc_read_next(:,:,:) = mc(:,:,:)
2398    END IF
2399   
2400    ! Initialize snow variables for nudging if not found in restart file
2401    IF (ok_nudge_snow) THEN
2402       IF ( ALL(snowdz_read_next(:,:)==val_exp) ) snowdz_read_next(:,:) = snowdz(:,:)
2403       IF ( ALL(snowrho_read_next(:,:)==val_exp) ) snowrho_read_next(:,:) = snowrho(:,:)
2404       IF ( ALL(snowtemp_read_next(:,:)==val_exp) ) snowtemp_read_next(:,:) = snowtemp(:,:)
2405    END IF
2406   
2407   
2408    IF (printlev>=3) WRITE (numout,*) ' hydrol_init done '
2409   
2410  END SUBROUTINE hydrol_init
2411
2412
2413!! ================================================================================================================================
2414!! SUBROUTINE   : hydrol_clear
2415!!
2416!>\BRIEF        Deallocate arrays
2417!!
2418!_ ================================================================================================================================
2419!_ hydrol_clear
2420
2421  SUBROUTINE hydrol_clear()
2422
2423    ! Allocation for soiltile related parameters
2424    IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
2425    IF ( ALLOCATED (avan)) DEALLOCATE (avan)
2426    IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
2427    IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
2428    IF ( ALLOCATED (ks)) DEALLOCATE (ks)
2429    IF ( ALLOCATED (pcent)) DEALLOCATE (pcent)
2430    IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
2431    IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
2432    IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet)
2433    IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry)
2434    ! Other arrays
2435    IF (ALLOCATED (mask_veget)) DEALLOCATE (mask_veget)
2436    IF (ALLOCATED (mask_soiltile)) DEALLOCATE (mask_soiltile)
2437    IF (ALLOCATED (humrelv)) DEALLOCATE (humrelv)
2438    IF (ALLOCATED (vegstressv)) DEALLOCATE (vegstressv)
2439    IF (ALLOCATED (us)) DEALLOCATE (us)
2440    IF (ALLOCATED  (precisol)) DEALLOCATE (precisol)
2441    IF (ALLOCATED  (throughfall)) DEALLOCATE (throughfall)
2442    IF (ALLOCATED  (precisol_ns)) DEALLOCATE (precisol_ns)
2443    IF (ALLOCATED  (free_drain_coef)) DEALLOCATE (free_drain_coef)
2444    IF (ALLOCATED  (frac_bare_ns)) DEALLOCATE (frac_bare_ns)
2445    IF (ALLOCATED  (water2infilt)) DEALLOCATE (water2infilt)
2446    IF (ALLOCATED  (ae_ns)) DEALLOCATE (ae_ns)
2447    IF (ALLOCATED  (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
2448    IF (ALLOCATED  (rootsink)) DEALLOCATE (rootsink)
2449    IF (ALLOCATED  (subsnowveg)) DEALLOCATE (subsnowveg)
2450    IF (ALLOCATED  (subsnownobio)) DEALLOCATE (subsnownobio)
2451    IF (ALLOCATED  (icemelt)) DEALLOCATE (icemelt)
2452    IF (ALLOCATED  (subsinksoil)) DEALLOCATE (subsinksoil)
2453    IF (ALLOCATED  (mx_eau_var)) DEALLOCATE (mx_eau_var)
2454    IF (ALLOCATED  (vegtot)) DEALLOCATE (vegtot)
2455    IF (ALLOCATED  (vegtot_old)) DEALLOCATE (vegtot_old)
2456    IF (ALLOCATED  (resdist)) DEALLOCATE (resdist)
2457    IF (ALLOCATED  (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
2458    IF (ALLOCATED  (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
2459    IF (ALLOCATED  (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
2460    IF (ALLOCATED  (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
2461    IF (ALLOCATED  (delsoilmoist)) DEALLOCATE (delsoilmoist)
2462    IF (ALLOCATED  (delintercept)) DEALLOCATE (delintercept)
2463    IF (ALLOCATED  (snow_beg)) DEALLOCATE (snow_beg)
2464    IF (ALLOCATED  (snow_end)) DEALLOCATE (snow_end)
2465    IF (ALLOCATED  (delswe)) DEALLOCATE (delswe)
2466    IF (ALLOCATED  (undermcr)) DEALLOCATE (undermcr)
2467    IF (ALLOCATED  (v1)) DEALLOCATE (v1)
2468    IF (ALLOCATED  (humtot)) DEALLOCATE (humtot)
2469    IF (ALLOCATED  (resolv)) DEALLOCATE (resolv)
2470    IF (ALLOCATED  (k)) DEALLOCATE (k)
2471    IF (ALLOCATED  (kk)) DEALLOCATE (kk)
2472    IF (ALLOCATED  (kk_moy)) DEALLOCATE (kk_moy)
2473    IF (ALLOCATED  (a)) DEALLOCATE (a)
2474    IF (ALLOCATED  (b)) DEALLOCATE (b)
2475    IF (ALLOCATED  (d)) DEALLOCATE (d)
2476    IF (ALLOCATED  (e)) DEALLOCATE (e)
2477    IF (ALLOCATED  (f)) DEALLOCATE (f)
2478    IF (ALLOCATED  (g1)) DEALLOCATE (g1)
2479    IF (ALLOCATED  (ep)) DEALLOCATE (ep)
2480    IF (ALLOCATED  (fp)) DEALLOCATE (fp)
2481    IF (ALLOCATED  (gp)) DEALLOCATE (gp)
2482    IF (ALLOCATED  (rhs)) DEALLOCATE (rhs)
2483    IF (ALLOCATED  (srhs)) DEALLOCATE (srhs)
2484    IF (ALLOCATED  (tmc)) DEALLOCATE (tmc)
2485    IF (ALLOCATED  (tmcs)) DEALLOCATE (tmcs)
2486    IF (ALLOCATED  (tmcr)) DEALLOCATE (tmcr)
2487    IF (ALLOCATED  (tmcfc)) DEALLOCATE (tmcfc)
2488    IF (ALLOCATED  (tmcw)) DEALLOCATE (tmcw)
2489    IF (ALLOCATED  (tmc_litter)) DEALLOCATE (tmc_litter)
2490    IF (ALLOCATED  (tmc_litt_mea)) DEALLOCATE (tmc_litt_mea)
2491    IF (ALLOCATED  (tmc_litter_res)) DEALLOCATE (tmc_litter_res)
2492    IF (ALLOCATED  (tmc_litter_wilt)) DEALLOCATE (tmc_litter_wilt)
2493    IF (ALLOCATED  (tmc_litter_field)) DEALLOCATE (tmc_litter_field)
2494    IF (ALLOCATED  (tmc_litter_sat)) DEALLOCATE (tmc_litter_sat)
2495    IF (ALLOCATED  (tmc_litter_awet)) DEALLOCATE (tmc_litter_awet)
2496    IF (ALLOCATED  (tmc_litter_adry)) DEALLOCATE (tmc_litter_adry)
2497    IF (ALLOCATED  (tmc_litt_wet_mea)) DEALLOCATE (tmc_litt_wet_mea)
2498    IF (ALLOCATED  (tmc_litt_dry_mea)) DEALLOCATE (tmc_litt_dry_mea)
2499    IF (ALLOCATED  (ru_ns)) DEALLOCATE (ru_ns)
2500    IF (ALLOCATED  (dr_ns)) DEALLOCATE (dr_ns)
2501    IF (ALLOCATED  (tr_ns)) DEALLOCATE (tr_ns)
2502    IF (ALLOCATED  (vegetmax_soil)) DEALLOCATE (vegetmax_soil)
2503    IF (ALLOCATED  (mc)) DEALLOCATE (mc)
2504    IF (ALLOCATED  (soilmoist)) DEALLOCATE (soilmoist)
2505    IF (ALLOCATED  (soilmoist_liquid)) DEALLOCATE (soilmoist_liquid)
2506    IF (ALLOCATED  (soil_wet_ns)) DEALLOCATE (soil_wet_ns)
2507    IF (ALLOCATED  (soil_wet_litter)) DEALLOCATE (soil_wet_litter)
2508    IF (ALLOCATED  (qflux)) DEALLOCATE (qflux)
2509    IF (ALLOCATED  (tmat)) DEALLOCATE (tmat)
2510    IF (ALLOCATED  (stmat)) DEALLOCATE (stmat)
2511    IF (ALLOCATED  (nroot)) DEALLOCATE (nroot)
2512    IF (ALLOCATED  (kfact_root)) DEALLOCATE (kfact_root)
2513    IF (ALLOCATED  (kfact)) DEALLOCATE (kfact)
2514    IF (ALLOCATED  (zz)) DEALLOCATE (zz)
2515    IF (ALLOCATED  (dz)) DEALLOCATE (dz)
2516    IF (ALLOCATED  (dh)) DEALLOCATE (dh)
2517    IF (ALLOCATED  (mc_lin)) DEALLOCATE (mc_lin)
2518    IF (ALLOCATED  (k_lin)) DEALLOCATE (k_lin)
2519    IF (ALLOCATED  (d_lin)) DEALLOCATE (d_lin)
2520    IF (ALLOCATED  (a_lin)) DEALLOCATE (a_lin)
2521    IF (ALLOCATED  (b_lin)) DEALLOCATE (b_lin)
2522
2523  END SUBROUTINE hydrol_clear
2524
2525!! ================================================================================================================================
2526!! SUBROUTINE   : hydrol_tmc_update
2527!!
2528!>\BRIEF        This routine updates the soil moisture profiles when the vegetation fraction have changed.
2529!!
2530!! DESCRIPTION  :
2531!!
2532!!    This routine update tmc and mc with variation of veget_max (LAND_USE or DGVM activated)
2533!!
2534!!
2535!!
2536!!
2537!! RECENT CHANGE(S) : Adaptation to excluding nobio from soiltile(1)
2538!!
2539!! MAIN OUTPUT VARIABLE(S) :
2540!!
2541!! REFERENCE(S) :
2542!!
2543!! FLOWCHART    : None
2544!! \n
2545!_ ================================================================================================================================
2546!_ hydrol_tmc_update
2547  SUBROUTINE hydrol_tmc_update ( kjpindex, veget_max, soiltile, qsintveg, drain_upd, runoff_upd)
2548
2549    !! 0.1 Input variables
2550    INTEGER(i_std), INTENT(in)                            :: kjpindex      !! domain size
2551    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)     :: veget_max     !! max fraction of vegetation type
2552    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)   :: soiltile      !! Fraction of each soil tile (0-1, unitless)
2553
2554    !! 0.2 Output variables
2555    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: drain_upd        !! Change in drainage due to decrease in vegtot
2556                                                                              !! on mc [kg/m2/dt]
2557    REAL(r_std),DIMENSION (kjpindex), INTENT(out)         :: runoff_upd       !! Change in runoff due to decrease in vegtot
2558                                                                              !! on water2infilt[kg/m2/dt]
2559   
2560    !! 0.3 Modified variables
2561    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)  :: qsintveg   !! Amount of water in the canopy interception
2562
2563    !! 0.4 Local variables
2564    INTEGER(i_std)                           :: ji, jv, jst,jsl
2565    LOGICAL                                  :: soil_upd        !! True if soiltile changed since last time step
2566    LOGICAL                                  :: vegtot_upd      !! True if vegtot changed since last time step
2567    LOGICAL                                  :: error=.FALSE.   !! If true, exit in the end of subroutine
2568    REAL(r_std), DIMENSION(kjpindex,nstm)    :: vmr             !! Change in soiltile (within vegtot)
2569    REAL(r_std), DIMENSION(kjpindex)         :: vmr_sum
2570    REAL(r_std), DIMENSION(kjpindex)         :: delvegtot   
2571    REAL(r_std), DIMENSION(kjpindex,nslm)    :: mc_dilu         !! Total loss of moisture content
2572    REAL(r_std), DIMENSION(kjpindex)         :: infil_dilu      !! Total loss for water2infilt
2573    REAL(r_std), DIMENSION(kjpindex,nstm)    :: tmc_old         !! tmc before calculations
2574    REAL(r_std), DIMENSION(kjpindex,nstm)    :: water2infilt_old!! water2infilt before calculations
2575    REAL(r_std), DIMENSION (kjpindex,nvm)    :: qsintveg_old    !! qsintveg before calculations
2576    REAL(r_std), DIMENSION(kjpindex)         :: test
2577    REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: mcaux        !! serves to hold the chnage in mc when vegtot decreases
2578
2579    !! 0. For checks
2580
2581    IF (check_cwrr) THEN
2582       ! Save soil moisture for later use
2583       tmc_old(:,:) = tmc(:,:) 
2584       water2infilt_old(:,:) = water2infilt(:,:)
2585       qsintveg_old(:,:) = qsintveg(:,:)
2586    ENDIF
2587   
2588    !! 1. If a PFT has disapperead as result from a veget_max change,
2589    !!    then add canopy water to surface water.
2590    !     Other adaptations of qsintveg are delt by the normal functioning of hydrol_canop
2591
2592    DO ji=1,kjpindex
2593       IF (vegtot_old(ji) .GT.min_sechiba) THEN
2594          DO jv=1,nvm
2595             IF ((veget_max(ji,jv).LT.min_sechiba).AND.(qsintveg(ji,jv).GT.0.)) THEN
2596                jst=pref_soil_veg(jv) ! soil tile index
2597                water2infilt(ji,jst) = water2infilt(ji,jst) + qsintveg(ji,jv)/(resdist(ji,jst)*vegtot_old(ji))
2598                qsintveg(ji,jv) = zero
2599             ENDIF
2600          ENDDO
2601       ENDIF
2602    ENDDO
2603   
2604    !! 2. We now deal with the changes of soiltile and corresponding soil moistures
2605    !!    Because sum(soiltile)=1 whatever vegtot, we need to distinguish two cases:
2606    !!    - when vegtot changes (meaning that the nobio fraction changes too),
2607    !!    - and when vegtot does not changes (a priori the most frequent case)
2608
2609    vegtot_upd = SUM(ABS((vegtot(:)-vegtot_old(:)))) .GT. zero ! True if at least one land point with a vegtot change
2610    runoff_upd(:) = zero
2611    drain_upd(:) = zero
2612    IF (vegtot_upd) THEN
2613       ! We find here the processing specific to the chnages of nobio fraction and vegtot
2614
2615       delvegtot(:) = vegtot(:) - vegtot_old(:)
2616
2617       DO jst=1,nstm
2618          DO ji=1,kjpindex
2619
2620             IF (delvegtot(ji) .GT. min_sechiba) THEN
2621
2622                !! 2.1. If vegtot increases (nobio decreases), then the mc in each soiltile is decreased
2623                !!      assuming the same proportions for each soiltile, and each soil layer
2624               
2625                mc(ji,:,jst) = mc(ji,:,jst) * vegtot_old(ji)/vegtot(ji) ! vegtot cannot be zero as > vegtot_old
2626                water2infilt(ji,jst) = water2infilt(ji,jst) * vegtot_old(ji)/vegtot(ji)
2627
2628             ELSE
2629
2630                !! 2.2 If vegtot decreases (nobio increases), then the mc in each soiltile should increase,
2631                !!     but should not exceed mcs
2632                !!     For simplicity, we choose to send the corresponding water volume to drainage
2633                !!     We do the same for water2infilt but send the excess to surface runoff
2634
2635                IF (vegtot(ji) .GT.min_sechiba) THEN
2636                   mcaux(ji,:,jst) =  mc(ji,:,jst) * (vegtot_old(ji)-vegtot(ji))/vegtot(ji) ! mcaux is the delta mc
2637                ELSE ! we just have nobio in the grid-cell
2638                   mcaux(ji,:,jst) =  mc(ji,:,jst)
2639                ENDIF
2640               
2641                drain_upd(ji) = drain_upd(ji) + dz(2) * ( trois*mcaux(ji,1,jst) + mcaux(ji,2,jst) )/huit
2642                DO jsl = 2,nslm-1
2643                   drain_upd(ji) = drain_upd(ji) + dz(jsl) * (trois*mcaux(ji,jsl,jst)+mcaux(ji,jsl-1,jst))/huit &
2644                        + dz(jsl+1) * (trois*mcaux(ji,jsl,jst)+mcaux(ji,jsl+1,jst))/huit
2645                ENDDO
2646                drain_upd(ji) = drain_upd(ji) + dz(nslm) * (trois*mcaux(ji,nslm,jst) + mcaux(ji,nslm-1,jst))/huit
2647
2648                IF (vegtot(ji) .GT.min_sechiba) THEN
2649                   runoff_upd(ji) = runoff_upd(ji) + water2infilt(ji,jst) * (vegtot_old(ji)-vegtot(ji))/vegtot(ji)
2650                ELSE ! we just have nobio in the grid-cell
2651                   runoff_upd(ji) = runoff_upd(ji) + water2infilt(ji,jst)
2652                ENDIF
2653
2654             ENDIF
2655             
2656          ENDDO
2657       ENDDO
2658       
2659    ENDIF
2660   
2661    !! 3. At the end of step 2, we are back to a case where vegtot changes are treated, so we can use soiltile
2662    !!    as a fraction of vegtot to process the mc transfers between soil tiles due to the changes of vegetation map
2663   
2664    !! 3.1 Check if soiltiles changed since last time step
2665    soil_upd=SUM(ABS(soiltile(:,:)-resdist(:,:))) .GT. zero
2666    IF (printlev>=3) WRITE (numout,*) 'soil_upd ', soil_upd
2667       
2668    IF (soil_upd) THEN
2669     
2670       !! 3.2 Define the change in soiltile
2671       vmr(:,:) = soiltile(:,:) - resdist(:,:)  ! resdist is the previous values of soiltiles, previous timestep, so before new map
2672
2673       ! Total area loss by the three soil tiles
2674       DO ji=1,kjpindex
2675          vmr_sum(ji)=SUM(vmr(ji,:),MASK=vmr(ji,:).LT.zero)
2676       ENDDO
2677
2678       !! 3.3 Shrinking soil tiles
2679       !! 3.3.1 Total loss of moisture content from the shrinking soil tiles, expressed by soil layer
2680       mc_dilu(:,:)=zero
2681       DO jst=1,nstm
2682          DO jsl = 1, nslm
2683             DO ji=1,kjpindex
2684                IF ( vmr(ji,jst) < -min_sechiba ) THEN
2685                   mc_dilu(ji,jsl) = mc_dilu(ji,jsl) + mc(ji,jsl,jst) * vmr(ji,jst) / vmr_sum(ji)
2686                ENDIF
2687             ENDDO
2688          ENDDO
2689       ENDDO
2690
2691       !! 3.3.2 Total loss of water2inft from the shrinking soil tiles
2692       infil_dilu(:)=zero
2693       DO jst=1,nstm
2694          DO ji=1,kjpindex
2695             IF ( vmr(ji,jst) < -min_sechiba ) THEN
2696                infil_dilu(ji) = infil_dilu(ji) + water2infilt(ji,jst) * vmr(ji,jst) / vmr_sum(ji)
2697             ENDIF
2698          ENDDO
2699       ENDDO
2700
2701       !! 3.4 Each gaining soil tile gets moisture proportionally to both the total loss and its areal increase
2702
2703       ! As the original mc from each soil tile are in [mcr,mcs] and we do weighted avrage, the new mc are in [mcr,mcs]
2704       ! The case where the soiltile is created (soiltile_old=0) works as the other cases
2705
2706       ! 3.4.1 Update mc(kjpindex,nslm,nstm) !m3/m3
2707       DO jst=1,nstm
2708          DO jsl = 1, nslm
2709             DO ji=1,kjpindex
2710                IF ( vmr(ji,jst) > min_sechiba ) THEN
2711                   mc(ji,jsl,jst) = ( mc(ji,jsl,jst) * resdist(ji,jst) + mc_dilu(ji,jsl) * vmr(ji,jst) ) / soiltile(ji,jst)
2712                   ! NB : soiltile can not be zero for case vmr > zero, see slowproc_veget
2713                ENDIF
2714             ENDDO
2715          ENDDO
2716       ENDDO
2717       
2718       ! 3.4.2 Update water2inft
2719       DO jst=1,nstm
2720          DO ji=1,kjpindex
2721             IF ( vmr(ji,jst) > min_sechiba ) THEN !donc soiltile>0     
2722                water2infilt(ji,jst) = ( water2infilt(ji,jst) * resdist(ji,jst) + infil_dilu(ji) * vmr(ji,jst) ) / soiltile(ji,jst)
2723             ENDIF !donc resdist>0
2724          ENDDO
2725       ENDDO
2726
2727       ! 3.4.3 Case where soiltile < min_sechiba
2728       DO jst=1,nstm
2729          DO ji=1,kjpindex
2730             IF ( soiltile(ji,jst) .LT. min_sechiba ) THEN
2731                water2infilt(ji,jst) = zero
2732                mc(ji,:,jst) = zero
2733             ENDIF
2734          ENDDO
2735       ENDDO
2736
2737    ENDIF ! soil_upd
2738
2739    !! 4. Update tmc and humtot
2740   
2741    DO jst=1,nstm
2742       DO ji=1,kjpindex
2743             tmc(ji,jst) = dz(2) * ( trois*mc(ji,1,jst) + mc(ji,2,jst) )/huit
2744             DO jsl = 2,nslm-1
2745                tmc(ji,jst) = tmc(ji,jst) + dz(jsl) * (trois*mc(ji,jsl,jst)+mc(ji,jsl-1,jst))/huit &
2746                     + dz(jsl+1) * (trois*mc(ji,jsl,jst)+mc(ji,jsl+1,jst))/huit
2747             ENDDO
2748             tmc(ji,jst) = tmc(ji,jst) + dz(nslm) * (trois*mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
2749             tmc(ji,jst) = tmc(ji,jst) + water2infilt(ji,jst)
2750             ! WARNING tmc is increased by includes water2infilt(ji,jst)
2751       ENDDO
2752    ENDDO
2753
2754    humtot(:) = zero
2755    DO jst=1,nstm
2756       DO ji=1,kjpindex
2757          humtot(ji) = humtot(ji) + vegtot(ji) * soiltile(ji,jst) * tmc(ji,jst) ! average over grid-cell (i.e. total land)
2758       ENDDO
2759    ENDDO
2760
2761    !! 5. Check
2762    IF (check_cwrr) THEN
2763       DO ji=1,kjpindex
2764          test(ji) = SUM(tmc(ji,:)*soiltile(ji,:)*vegtot(ji)) - SUM(tmc_old(ji,:)*resdist(ji,:)*vegtot_old(ji)) + &
2765               SUM(qsintveg(ji,:)) - SUM(qsintveg_old(ji,:)) + (drain_upd(ji) + runoff_upd(ji))   
2766          IF ( ABS(test(ji)) .GT.  10.*allowed_err ) THEN
2767             WRITE(numout,*) 'tmc update WRONG: ji',ji
2768             WRITE(numout,*) 'tot water avant:',SUM(tmc_old(ji,:)*resdist(ji,:)*vegtot_old(ji)) + SUM(qsintveg_old(ji,:))
2769             WRITE(numout,*) 'tot water apres:',SUM(tmc(ji,:)*soiltile(ji,:)*vegtot(ji)) + SUM(qsintveg(ji,:))
2770             WRITE(numout,*) 'err:',test(ji)
2771             WRITE(numout,*) 'allowed_err:',allowed_err
2772             WRITE(numout,*) 'tmc:',tmc(ji,:)
2773             WRITE(numout,*) 'tmc_old:',tmc_old(ji,:)
2774             WRITE(numout,*) 'qsintveg:',qsintveg(ji,:)
2775             WRITE(numout,*) 'qsintveg_old:',qsintveg_old(ji,:)
2776             WRITE(numout,*) 'SUMqsintveg:',SUM(qsintveg(ji,:))
2777             WRITE(numout,*) 'SUMqsintveg_old:',SUM(qsintveg_old(ji,:))
2778             WRITE(numout,*) 'veget_max:',veget_max(ji,:)
2779             WRITE(numout,*) 'soiltile:',soiltile(ji,:)
2780             WRITE(numout,*) 'resdist:',resdist(ji,:)
2781             WRITE(numout,*) 'vegtot:',vegtot(ji)
2782             WRITE(numout,*) 'vegtot_old:',vegtot_old(ji)
2783             WRITE(numout,*) 'drain_upd:',drain_upd(ji)
2784             WRITE(numout,*) 'runoff_upd:',runoff_upd(ji)
2785             WRITE(numout,*) 'vmr:',vmr(ji,:)
2786             WRITE(numout,*) 'vmr_sum:',vmr_sum(ji)
2787             DO jst=1,nstm
2788                WRITE(numout,*) 'mc(',jst,'):',mc(ji,:,jst)
2789             ENDDO
2790             WRITE(numout,*) 'water2infilt:',water2infilt(ji,:)
2791             WRITE(numout,*) 'water2infilt_old:',water2infilt_old(ji,:)
2792             WRITE(numout,*) 'infil_dilu:',infil_dilu(ji)
2793             WRITE(numout,*) 'mc_dilu:',mc_dilu(ji,:)
2794
2795             error=.TRUE.
2796             CALL ipslerr_p(2, 'hydrol_tmc_update', 'Error in water balance', 'We STOP in the end of this subroutine','')
2797          ENDIF
2798       ENDDO
2799    ENDIF
2800
2801    !! Now that the work is done, update resdist
2802    resdist(:,:) = soiltile(:,:)
2803
2804    !
2805    !!  Exit if error was found previously in this subroutine
2806    !
2807    IF ( error ) THEN
2808       WRITE(numout,*) 'One or more errors have been detected in hydrol_tmc_update. Model stops.'
2809       CALL ipslerr_p(3, 'hydrol_tmc_update', 'We will STOP now.',&
2810                  & 'One or several fatal errors were found previously.','')
2811    END IF
2812
2813    IF (printlev>=3) WRITE (numout,*) ' hydrol_tmc_update done '
2814
2815  END SUBROUTINE hydrol_tmc_update
2816
2817!! ================================================================================================================================
2818!! SUBROUTINE   : hydrol_var_init
2819!!
2820!>\BRIEF        This routine initializes hydrologic parameters to define K and D, and diagnostic hydrologic variables. 
2821!!
2822!! DESCRIPTION  :
2823!! - 1 compute the depths
2824!! - 2 compute the profile for roots
2825!! - 3 compute the profile for ksat, a and n Van Genuchten parameter
2826!! - 4 compute the linearized values of k, a, b and d for the resolution of Fokker Planck equation
2827!! - 5 water reservoirs initialisation
2828!!
2829!! RECENT CHANGE(S) : None
2830!!
2831!! MAIN OUTPUT VARIABLE(S) :
2832!!
2833!! REFERENCE(S) :
2834!!
2835!! FLOWCHART    : None
2836!! \n
2837!_ ================================================================================================================================
2838!_ hydrol_var_init
2839
2840  SUBROUTINE hydrol_var_init (kjpindex, veget, veget_max, soiltile, njsc, &
2841       mx_eau_var, shumdiag_perma, &
2842       drysoil_frac, qsintveg, mc_layh, mcl_layh) 
2843
2844    ! interface description
2845
2846    !! 0. Variable and parameter declaration
2847
2848    !! 0.1 Input variables
2849
2850    ! input scalar
2851    INTEGER(i_std), INTENT(in)                          :: kjpindex      !! Domain size (number of grid cells) (1)
2852    ! input fields
2853    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max     !! PFT fractions within grid-cells (1; 1)
2854    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget         !! Effective fraction of vegetation by PFT (1; 1)
2855    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc          !! Index of the dominant soil textural class
2856                                                                         !! in the grid cell (1-nscm, unitless)
2857    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltile      !! Fraction of each soil tile within vegtot (0-1, unitless)
2858
2859    !! 0.2 Output variables
2860
2861    REAL(r_std),DIMENSION (kjpindex), INTENT (out)      :: mx_eau_var    !! Maximum water content of the soil
2862                                                                         !! @tex $(kg m^{-2})$ @endtex
2863    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out) :: shumdiag_perma!! Percent of porosity filled with water (mc/mcs)
2864                                                                         !! used for the thermal computations
2865    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)    :: drysoil_frac  !! function of litter humidity
2866    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out):: mc_layh       !! Volumetric soil moisture content for each layer in hydrol(liquid+ice) [m3/m3]
2867    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out):: mcl_layh      !! Volumetric soil moisture content for each layer in hydrol(liquid) [m3/m3]
2868
2869    !! 0.3 Modified variables
2870    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout)  :: qsintveg    !! Water on vegetation due to interception
2871                                                                         !! @tex $(kg m^{-2})$ @endtex 
2872
2873    !! 0.4 Local variables
2874
2875    INTEGER(i_std)                                      :: ji, jv        !! Grid-cell and PFT indices (1)
2876    INTEGER(i_std)                                      :: jst, jsc, jsl !! Soiltile, Soil Texture, and Soil layer indices (1)
2877    INTEGER(i_std)                                      :: i             !! Index (1)
2878    REAL(r_std)                                         :: m             !! m=1-1/n (unitless)
2879    REAL(r_std)                                         :: frac          !! Relative linearized VWC (unitless)
2880    REAL(r_std)                                         :: avan_mod      !! VG parameter a modified from  exponantial profile
2881                                                                         !! @tex $(mm^{-1})$ @endtex
2882    REAL(r_std)                                         :: nvan_mod      !! VG parameter n  modified from  exponantial profile
2883                                                                         !! (unitless)
2884    REAL(r_std), DIMENSION(nslm,nscm)                   :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
2885                                                                         !! (unitless)
2886    ! parameters for "soil densification" with depth
2887    REAL(r_std)                                         :: dp_comp       !! Depth at which the 'compacted' value of ksat
2888                                                                         !! is reached (m)
2889    REAL(r_std)                                         :: f_ks          !! Exponential factor for decay of ksat with depth
2890                                                                         !! @tex $(m^{-1})$ @endtex
2891    ! Fixed parameters from fitted relationships
2892    REAL(r_std)                                         :: n0            !! fitted value for relation log((n-n0)/(n_ref-n0)) =
2893                                                                         !! nk_rel * log(k/k_ref)
2894                                                                         !! (unitless)
2895    REAL(r_std)                                         :: nk_rel        !! fitted value for relation log((n-n0)/(n_ref-n0)) =
2896                                                                         !! nk_rel * log(k/k_ref)
2897                                                                         !! (unitless)
2898    REAL(r_std)                                         :: a0            !! fitted value for relation log((a-a0)/(a_ref-a0)) =
2899                                                                         !! ak_rel * log(k/k_ref)
2900                                                                         !! @tex $(mm^{-1})$ @endtex
2901    REAL(r_std)                                         :: ak_rel        !! fitted value for relation log((a-a0)/(a_ref-a0)) =
2902                                                                         !! ak_rel * log(k/k_ref)
2903                                                                         !! (unitless)
2904    REAL(r_std)                                         :: kfact_max     !! Maximum factor for Ks decay with depth (unitless)
2905    REAL(r_std)                                         :: k_tmp, tmc_litter_ratio
2906    INTEGER(i_std), PARAMETER                           :: error_level = 3 !! Error level for consistency check
2907                                                                           !! Switch to 2 tu turn fatal errors into warnings
2908    REAL(r_std), DIMENSION (kjpindex,nslm)              :: ksat            !! Saturated hydraulic conductivity at each node (mm/d)
2909    INTEGER(i_std)                                      :: jiref           !! To identify the mc_lins where k_lin and d_lin
2910                                                                           !! need special treatment
2911
2912!_ ================================================================================================================================
2913
2914!!??Aurelien: Les 3 parametres qui suivent pourait peut-être mis dans hydrol_init?
2915    !
2916    !
2917    !Config Key   = CWRR_NKS_N0
2918    !Config Desc  = fitted value for relation log((n-n0)/(n_ref-n0)) = nk_rel * log(k/k_ref)
2919    !Config Def   = 0.95
2920    !Config If    = HYDROL_CWRR
2921    !Config Help  =
2922    !Config Units = [-]
2923    n0 = 0.95
2924    CALL getin_p("CWRR_NKS_N0",n0)
2925
2926    !! Check parameter value (correct range)
2927    IF ( n0 < zero ) THEN
2928       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2929            &     "Wrong parameter value for CWRR_NKS_N0.", &
2930            &     "This parameter should be non-negative. ", &
2931            &     "Please, check parameter value in run.def. ")
2932    END IF
2933
2934
2935    !Config Key   = CWRR_NKS_POWER
2936    !Config Desc  = fitted value for relation log((n-n0)/(n_ref-n0)) = nk_rel * log(k/k_ref)
2937    !Config Def   = 0.34
2938    !Config If    = HYDROL_CWRR
2939    !Config Help  =
2940    !Config Units = [-]
2941    nk_rel = 0.34
2942    CALL getin_p("CWRR_NKS_POWER",nk_rel)
2943
2944    !! Check parameter value (correct range)
2945    IF ( nk_rel < zero ) THEN
2946       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2947            &     "Wrong parameter value for CWRR_NKS_POWER.", &
2948            &     "This parameter should be non-negative. ", &
2949            &     "Please, check parameter value in run.def. ")
2950    END IF
2951
2952
2953    !Config Key   = CWRR_AKS_A0
2954    !Config Desc  = fitted value for relation log((a-a0)/(a_ref-a0)) = ak_rel * log(k/k_ref)
2955    !Config Def   = 0.00012
2956    !Config If    = HYDROL_CWRR
2957    !Config Help  =
2958    !Config Units = [1/mm]
2959    a0 = 0.00012
2960    CALL getin_p("CWRR_AKS_A0",a0)
2961
2962    !! Check parameter value (correct range)
2963    IF ( a0 < zero ) THEN
2964       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2965            &     "Wrong parameter value for CWRR_AKS_A0.", &
2966            &     "This parameter should be non-negative. ", &
2967            &     "Please, check parameter value in run.def. ")
2968    END IF
2969
2970
2971    !Config Key   = CWRR_AKS_POWER
2972    !Config Desc  = fitted value for relation log((a-a0)/(a_ref-a0)) = ak_rel * log(k/k_ref)
2973    !Config Def   = 0.53
2974    !Config If    = HYDROL_CWRR
2975    !Config Help  =
2976    !Config Units = [-]
2977    ak_rel = 0.53
2978    CALL getin_p("CWRR_AKS_POWER",ak_rel)
2979
2980    !! Check parameter value (correct range)
2981    IF ( nk_rel < zero ) THEN
2982       CALL ipslerr_p(error_level, "hydrol_var_init.", &
2983            &     "Wrong parameter value for CWRR_AKS_POWER.", &
2984            &     "This parameter should be non-negative. ", &
2985            &     "Please, check parameter value in run.def. ")
2986    END IF
2987
2988
2989    !Config Key   = KFACT_DECAY_RATE
2990    !Config Desc  = Factor for Ks decay with depth
2991    !Config Def   = 2.0
2992    !Config If    = HYDROL_CWRR
2993    !Config Help  = 
2994    !Config Units = [1/m]
2995    f_ks = 2.0
2996    CALL getin_p ("KFACT_DECAY_RATE", f_ks)
2997
2998    !! Check parameter value (correct range)
2999    IF ( f_ks < zero ) THEN
3000       CALL ipslerr_p(error_level, "hydrol_var_init.", &
3001            &     "Wrong parameter value for KFACT_DECAY_RATE.", &
3002            &     "This parameter should be positive. ", &
3003            &     "Please, check parameter value in run.def. ")
3004    END IF
3005
3006
3007    !Config Key   = KFACT_STARTING_DEPTH
3008    !Config Desc  = Depth for compacted value of Ks
3009    !Config Def   = 0.3
3010    !Config If    = HYDROL_CWRR
3011    !Config Help  = 
3012    !Config Units = [m]
3013    dp_comp = 0.3
3014    CALL getin_p ("KFACT_STARTING_DEPTH", dp_comp)
3015
3016    !! Check parameter value (correct range)
3017    IF ( dp_comp <= zero ) THEN
3018       CALL ipslerr_p(error_level, "hydrol_var_init.", &
3019            &     "Wrong parameter value for KFACT_STARTING_DEPTH.", &
3020            &     "This parameter should be positive. ", &
3021            &     "Please, check parameter value in run.def. ")
3022    END IF
3023
3024
3025    !Config Key   = KFACT_MAX
3026    !Config Desc  = Maximum Factor for Ks increase due to vegetation
3027    !Config Def   = 10.0
3028    !Config If    = HYDROL_CWRR
3029    !Config Help  =
3030    !Config Units = [-]
3031    kfact_max = 10.0
3032    CALL getin_p ("KFACT_MAX", kfact_max)
3033
3034    !! Check parameter value (correct range)
3035    IF ( kfact_max < 10. ) THEN
3036       CALL ipslerr_p(error_level, "hydrol_var_init.", &
3037            &     "Wrong parameter value for KFACT_MAX.", &
3038            &     "This parameter should be greater than 10. ", &
3039            &     "Please, check parameter value in run.def. ")
3040    END IF
3041
3042   
3043    !-
3044    !! 1 Create local variables in mm for the vertical depths
3045    !!   Vertical depth variables (znh, dnh, dlh) are stored in module vertical_soil_var in m.
3046    DO jsl=1,nslm
3047       zz(jsl) = znh(jsl)*mille
3048       dz(jsl) = dnh(jsl)*mille
3049       dh(jsl) = dlh(jsl)*mille
3050    ENDDO
3051
3052    !-
3053    !! 2 Compute the root density profile if not ok_dynroot
3054    !!   For the case with ok_dynroot, the calculations are done at each time step in hydrol_soil
3055    IF (.NOT. ok_dynroot) THEN
3056       DO ji=1, kjpindex
3057          !-
3058          !! The three following equations concerning nroot computation are derived from the integrals
3059          !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158).
3060          !! The occasional absence of minus sign before humcste parameter is correct.
3061          DO jv = 1,nvm
3062             DO jsl = 2, nslm-1
3063                nroot(ji,jv,jsl) = (EXP(-humcste(jv)*zz(jsl)/mille)) * &
3064                     & (EXP(humcste(jv)*dz(jsl)/mille/deux) - &
3065                     & EXP(-humcste(jv)*dz(jsl+1)/mille/deux))/ &
3066                     & (EXP(-humcste(jv)*dz(2)/mille/deux) &
3067                     & -EXP(-humcste(jv)*zz(nslm)/mille))
3068             ENDDO
3069             nroot(ji,jv,1) = zero
3070
3071             nroot(ji,jv,nslm) = (EXP(humcste(jv)*dz(nslm)/mille/deux) -un) * &
3072                  & EXP(-humcste(jv)*zz(nslm)/mille) / &
3073                  & (EXP(-humcste(jv)*dz(2)/mille/deux) &
3074                  & -EXP(-humcste(jv)*zz(nslm)/mille))
3075          ENDDO
3076       ENDDO
3077    END IF
3078
3079    !-
3080    !! 3 Compute the profile for ksat, a and n
3081    !-
3082
3083    ! For every soil texture
3084    DO jsc = 1, nscm 
3085       DO jsl=1,nslm
3086          ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2
3087          ! Calibrated against Hapex-Sahel measurements
3088          kfact(jsl,jsc) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
3089          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14
3090         
3091          nfact(jsl,jsc) = ( kfact(jsl,jsc) )**nk_rel
3092          afact(jsl,jsc) = ( kfact(jsl,jsc) )**ak_rel
3093       ENDDO
3094    ENDDO
3095
3096    ! Output of ksat at each node for CMIP6
3097    DO jsl = 1, nslm
3098       ksat(:,jsl) = ks(njsc(:)) * kfact(jsl,njsc(:))
3099    ENDDO
3100    CALL xios_orchidee_send_field("ksat",ksat) ! mm/d
3101
3102    ! For every soil texture
3103    DO jsc = 1, nscm
3104       !-
3105       !! 4 compute the linearized values of k, a, b and d
3106       !-
3107       ! Calculate the matrix coef for Dublin model (de Rosnay, 1999; p149)
3108       ! piece-wise linearised hydraulic conductivity k_lin=alin * mc_lin + b_lin
3109       ! and diffusivity d_lin in each interval of mc, called mc_lin,
3110       ! between imin, for residual mcr, and imax for saturation mcs.
3111
3112       ! We define 51 bounds for 50 bins of mc between mcr and mcs
3113       mc_lin(imin,jsc)=mcr(jsc)
3114       mc_lin(imax,jsc)=mcs(jsc)
3115       DO ji= imin+1, imax-1 ! ji=2,50
3116          mc_lin(ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin)
3117       ENDDO
3118
3119       DO jsl = 1, nslm
3120          ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42
3121          nvan_mod = n0 + (nvan(jsc)-n0) * nfact(jsl,jsc)
3122          avan_mod = a0 + (avan(jsc)-a0) * afact(jsl,jsc)
3123          m = un - un / nvan_mod
3124          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(jsc) * kfact(jsl,jsc)
3125          DO ji = imax,imin,-1 
3126             frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
3127             k_lin(ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
3128          ENDDO
3129
3130          ! k_lin should not be zero, nor too small
3131          ! We track jiref, the bin under which mc is too small and we may get zero k_lin     
3132          ji=imax-1
3133          DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0))
3134             jiref=ji
3135             ji=ji-1
3136          ENDDO
3137          DO ji=jiref-1,imin,-1
3138             k_lin(ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10.
3139          ENDDO
3140         
3141          DO ji = imin,imax-1 ! ji=1,50
3142             ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin
3143             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))
3144             b_lin(ji,jsl,jsc)  = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc)
3145
3146             ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta)
3147             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin
3148             IF (ji.NE.imin .AND. ji.NE.imax-1) THEN
3149                frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
3150                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) *  &
3151                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) * &
3152                     (  frac**(-un/m) -un ) ** (-m)
3153                frac=MIN(un,(mc_lin(ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
3154                d_lin(ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*&
3155                     ( (frac**(-un/m))/(mc_lin(ji+1,jsc)-mcr(jsc)) ) * &
3156                     (  frac**(-un/m) -un ) ** (-m)
3157                d_lin(ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc))
3158             ELSE IF(ji.EQ.imax-1) THEN
3159                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &
3160                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) *  &
3161                     (  frac**(-un/m) -un ) ** (-m)
3162             ENDIF
3163          ENDDO
3164
3165          ! Special case for ji=imin
3166          d_lin(imin,jsl,jsc) = d_lin(imin+1,jsl,jsc)/1000.
3167
3168          ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations
3169          ! We don't want d_lin = zero
3170          DO ji=jiref-1,imin,-1
3171             d_lin(ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10.
3172          ENDDO
3173
3174       ENDDO
3175    ENDDO
3176
3177    !! 5 Water reservoir initialisation
3178    !
3179!!$    DO jst = 1,nstm
3180!!$       DO ji = 1, kjpindex
3181!!$          mx_eau_var(ji) = mx_eau_var(ji) + soiltile(ji,jst)*&
3182!!$               &   zmaxh*mille*mcs(njsc(ji))
3183!!$       END DO
3184!!$    END DO
3185!!$    IF (check_CWRR) THEN
3186!!$       IF ( ANY ( ABS( mx_eau_var(:) - zmaxh*mille*mcs(njsc(:)) ) > min_sechiba ) ) THEN
3187!!$          ji=MAXLOC ( ABS( mx_eau_var(:) - zmaxh*mille*mcs(njsc(:)) ) , 1)
3188!!$          WRITE(numout, *) "Erreur formule simplifiée mx_eau_var ! ", mx_eau_var(ji), zmaxh*mille*mcs(njsc(ji))
3189!!$          WRITE(numout, *) "err = ",ABS(mx_eau_var(ji) - zmaxh*mille*mcs(njsc(ji)))
3190!!$          STOP 1
3191!!$       ENDIF
3192!!$    ENDIF
3193
3194    mx_eau_var(:) = zero
3195    mx_eau_var(:) = zmaxh*mille*mcs(njsc(:)) 
3196
3197    DO ji = 1,kjpindex 
3198       IF (vegtot(ji) .LE. zero) THEN
3199          mx_eau_var(ji) = mx_eau_nobio*zmaxh
3200          ! Aurelien: what does vegtot=0 mean? is it like frac_nobio=1? But if 0<frac_nobio<1 ???
3201       ENDIF
3202
3203    END DO
3204
3205    ! Compute the litter humidity, shumdiag and fry
3206    shumdiag_perma(:,:) = zero
3207    humtot(:) = zero
3208    tmc(:,:) = zero
3209
3210    ! Loop on soiltiles to compute the variables (ji,jst)
3211    DO jst=1,nstm 
3212       DO ji = 1, kjpindex
3213          tmcs(ji,jst)=zmaxh* mille*mcs(njsc(ji))
3214          tmcr(ji,jst)=zmaxh* mille*mcr(njsc(ji))
3215          tmcfc(ji,jst)=zmaxh* mille*mcfc(njsc(ji))
3216          tmcw(ji,jst)=zmaxh* mille*mcw(njsc(ji))
3217       ENDDO
3218    ENDDO
3219       
3220    ! The total soil moisture for each soiltile:
3221    DO jst=1,nstm
3222       DO ji=1,kjpindex
3223          tmc(ji,jst)= dz(2) * ( trois*mc(ji,1,jst)+ mc(ji,2,jst))/huit
3224       END DO
3225    ENDDO
3226
3227    DO jst=1,nstm 
3228       DO jsl=2,nslm-1
3229          DO ji=1,kjpindex
3230             tmc(ji,jst) = tmc(ji,jst) + dz(jsl) * ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
3231                  & + dz(jsl+1)*(trois*mc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
3232          END DO
3233       END DO
3234    ENDDO
3235
3236    DO jst=1,nstm 
3237       DO ji=1,kjpindex
3238          tmc(ji,jst) = tmc(ji,jst) +  dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
3239          tmc(ji,jst) = tmc(ji,jst) + water2infilt(ji,jst)
3240       ENDDO
3241    END DO
3242
3243!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.   
3244!    ! If veget has been updated before restart (with LAND USE or DGVM),
3245!    ! tmc and mc must be modified with respect to humtot conservation.
3246!   CALL hydrol_tmc_update ( kjpindex, veget_max, soiltile, qsintveg)
3247
3248    ! The litter variables:
3249    ! level 1
3250    DO jst=1,nstm 
3251       DO ji=1,kjpindex
3252          tmc_litter(ji,jst) = dz(2) * (trois*mc(ji,1,jst)+mc(ji,2,jst))/huit
3253          tmc_litter_wilt(ji,jst) = dz(2) * mcw(njsc(ji)) / deux
3254          tmc_litter_res(ji,jst) = dz(2) * mcr(njsc(ji)) / deux
3255          tmc_litter_field(ji,jst) = dz(2) * mcfc(njsc(ji)) / deux
3256          tmc_litter_sat(ji,jst) = dz(2) * mcs(njsc(ji)) / deux
3257          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(njsc(ji)) / deux
3258          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(njsc(ji)) / deux
3259       ENDDO
3260    END DO
3261    ! sum from level 2 to 4
3262    DO jst=1,nstm 
3263       DO jsl=2,4
3264          DO ji=1,kjpindex
3265             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 
3266                  & ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
3267                  & + dz(jsl+1)*(trois*mc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
3268             tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + &
3269                  &(dz(jsl)+ dz(jsl+1))*& 
3270                  & mcw(njsc(ji))/deux
3271             tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + &
3272                  &(dz(jsl)+ dz(jsl+1))*& 
3273                  & mcr(njsc(ji))/deux
3274             tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + &
3275                  &(dz(jsl)+ dz(jsl+1))* & 
3276                  & mcs(njsc(ji))/deux
3277             tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + &
3278                  & (dz(jsl)+ dz(jsl+1))* & 
3279                  & mcfc(njsc(ji))/deux
3280             tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + &
3281                  &(dz(jsl)+ dz(jsl+1))* & 
3282                  & mc_awet(njsc(ji))/deux
3283             tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + &
3284                  & (dz(jsl)+ dz(jsl+1))* & 
3285                  & mc_adry(njsc(ji))/deux
3286          END DO
3287       END DO
3288    END DO
3289
3290
3291    DO jst=1,nstm 
3292       DO ji=1,kjpindex
3293          ! here we set that humrelv=0 in PFT1
3294          humrelv(ji,1,jst) = zero
3295       ENDDO
3296    END DO
3297
3298
3299    ! Calculate shumdiag_perma for thermosoil
3300    ! Use resdist instead of soiltile because we here need to have
3301    ! shumdiag_perma at the value from previous time step.
3302    ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma
3303    ! (based on resdist=soiltile from previous timestep, but normally equal to soiltile)
3304    ! For consistency with hydrol_soil, we want to calculate a grid-cell average
3305    soilmoist(:,:) = zero
3306    DO jst=1,nstm
3307       DO ji=1,kjpindex
3308          soilmoist(ji,1) = soilmoist(ji,1) + resdist(ji,jst) * &
3309               dz(2) * ( trois*mc(ji,1,jst) + mc(ji,2,jst) )/huit
3310          DO jsl = 2,nslm-1
3311             soilmoist(ji,jsl) = soilmoist(ji,jsl) + resdist(ji,jst) * &
3312                  ( dz(jsl) * (trois*mc(ji,jsl,jst)+mc(ji,jsl-1,jst))/huit &
3313                  + dz(jsl+1) * (trois*mc(ji,jsl,jst)+mc(ji,jsl+1,jst))/huit )
3314          END DO
3315          soilmoist(ji,nslm) = soilmoist(ji,nslm) + resdist(ji,jst) * &
3316               dz(nslm) * (trois*mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
3317       ENDDO
3318    ENDDO
3319    DO ji=1,kjpindex
3320        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average
3321    ENDDO
3322   
3323    ! -- shumdiag_perma for restart
3324    !  For consistency with hydrol_soil, we want to calculate a grid-cell average
3325    DO jsl = 1, nslm
3326       DO ji=1,kjpindex       
3327          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji)))
3328          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
3329       ENDDO
3330    ENDDO
3331               
3332    ! Calculate drysoil_frac if it was not found in the restart file
3333    ! For simplicity, we set drysoil_frac to 0.5 in this case
3334    IF (ALL(drysoil_frac(:) == val_exp)) THEN
3335       DO ji=1,kjpindex
3336          drysoil_frac(ji) = 0.5
3337       END DO
3338    END IF
3339
3340    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in
3341    !! thermosoil for the thermal conductivity.
3342    ! These values are only used in thermosoil_init in absence of a restart file
3343    mc_layh(:,:) = zero
3344    mcl_layh(:,:) = zero
3345    DO jst=1,nstm
3346       DO jsl=1,nslm
3347          DO ji=1,kjpindex
3348            mc_layh(ji,jsl) = mc_layh(ji,jsl) + mc(ji,jsl,jst) * resdist(ji,jst)  * vegtot_old(ji)
3349            mcl_layh(ji,jsl) = mcl_layh(ji,jsl) + mcl(ji,jsl,jst) * resdist(ji,jst) * vegtot_old(ji)
3350         ENDDO
3351      END DO
3352    END DO
3353
3354    IF (printlev>=3) WRITE (numout,*) ' hydrol_var_init done '
3355
3356  END SUBROUTINE hydrol_var_init
3357
3358
3359!! ================================================================================================================================
3360!! SUBROUTINE   : hydrol_snow
3361!!
3362!>\BRIEF        This routine computes snow processes.
3363!!
3364!! DESCRIPTION  :
3365!! - 0 initialisation
3366!! - 1 On vegetation
3367!! - 1.1 Compute snow masse
3368!! - 1.2 Sublimation
3369!! - 1.2.1 Check that sublimation on the vegetated fraction is possible.
3370!! - 1.3. snow melt only if temperature positive
3371!! - 1.3.1 enough snow for melting or not
3372!! - 1.3.2 not enough snow
3373!! - 1.3.3 negative snow - now snow melt
3374!! - 1.4 Snow melts only on weight glaciers
3375!! - 2 On Land ice
3376!! - 2.1 Compute snow
3377!! - 2.2 Sublimation
3378!! - 2.3 Snow melt only for continental ice fraction
3379!! - 2.3.1 If there is snow on the ice-fraction it can melt
3380!! - 2.4 Snow melts only on weight glaciers
3381!! - 3 On other surface types - not done yet
3382!! - 4 computes total melt (snow and ice)
3383!! - 5 computes snow age on veg and ice (for albedo)
3384!! - 5.1 Snow age on vegetation
3385!! - 5.2 Snow age on ice
3386!! - 6 Diagnose the depth of the snow layer
3387!!
3388!! RECENT CHANGE(S) : None
3389!!
3390!! MAIN OUTPUT VARIABLE(S) :
3391!!
3392!! REFERENCE(S) :
3393!!
3394!! FLOWCHART    : None
3395!! \n
3396!_ ================================================================================================================================
3397!_ hydrol_snow
3398
3399  SUBROUTINE hydrol_snow (kjpindex, precip_rain, precip_snow , temp_sol_new, soilcap,&
3400       & frac_nobio, totfrac_nobio, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
3401       & tot_melt, snowdepth,snowmelt)
3402
3403    !
3404    ! interface description
3405
3406    !! 0. Variable and parameter declaration
3407
3408    !! 0.1 Input variables
3409
3410    ! input scalar
3411    INTEGER(i_std), INTENT(in)                               :: kjpindex      !! Domain size
3412    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_rain   !! Rainfall
3413    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_snow   !! Snow precipitation
3414    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: temp_sol_new  !! New soil temperature
3415    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: soilcap       !! Soil capacity
3416    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in)     :: frac_nobio    !! Fraction of continental ice, lakes, ...
3417    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
3418
3419    !! 0.2 Output variables
3420
3421    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: tot_melt      !! Total melt from snow and ice 
3422    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: snowmelt      !! Snow melt
3423    REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: snowdepth     !! Snow depth
3424
3425    !! 0.3 Modified variables
3426
3427    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: vevapsno      !! Snow evaporation
3428    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: snow          !! Snow mass [Kg/m^2]
3429    REAL(r_std), DIMENSION (kjpindex), INTENT(inout)         :: snow_age      !! Snow age
3430    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout)  :: snow_nobio    !! Ice water balance
3431    REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout)  :: snow_nobio_age!! Snow age on ice, lakes, ...
3432
3433    !! 0.4 Local variables
3434
3435    INTEGER(i_std)                               :: ji, jv
3436    REAL(r_std), DIMENSION (kjpindex)             :: d_age  !! Snow age change
3437    REAL(r_std), DIMENSION (kjpindex)             :: xx     !! temporary
3438    REAL(r_std)                                   :: snowmelt_tmp !! The name says it all !
3439    REAL(r_std)                                   :: snow_d1k !! The amount of snow that corresponds to a 1K cooling
3440
3441!_ ================================================================================================================================
3442
3443    !
3444    ! for continental points
3445    !
3446
3447    !
3448    !!_0 initialisation
3449    !
3450    DO jv = 1, nnobio
3451       DO ji=1,kjpindex
3452          subsnownobio(ji,jv) = zero
3453       ENDDO
3454    ENDDO
3455    DO ji=1,kjpindex
3456       subsnowveg(ji) = zero
3457       snowmelt(ji) = zero
3458       icemelt(ji) = zero
3459       subsinksoil(ji) = zero
3460       tot_melt(ji) = zero
3461    ENDDO
3462    !
3463    !! 1 On vegetation
3464    !
3465    DO ji=1,kjpindex
3466       !
3467    !! 1.1 Compute snow masse
3468       !
3469       snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji)
3470       !
3471       !
3472    !! 1.2 Sublimation
3473       !      Separate between vegetated and no-veget fractions
3474       !      Care has to be taken as we might have sublimation from the
3475       !      the frac_nobio while there is no snow on the rest of the grid.
3476       !
3477       IF ( snow(ji) > snowcri ) THEN
3478          subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
3479          subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
3480       ELSE
3481          ! Correction Nathalie - Juillet 2006.
3482          ! On doit d'abord tester s'il existe un frac_nobio!
3483          ! Pour le moment je ne regarde que le iice
3484          IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
3485             subsnownobio(ji,iice) = vevapsno(ji)
3486             subsnowveg(ji) = zero
3487          ELSE
3488             subsnownobio(ji,iice) = zero
3489             subsnowveg(ji) = vevapsno(ji)
3490          ENDIF
3491       ENDIF
3492       ! here vevapsno bas been separated into a bio and nobio fractions, without changing the total
3493       !
3494       !
3495    !! 1.2.1 Check that sublimation on the vegetated fraction is possible.
3496       !
3497       IF (subsnowveg(ji) .GT. snow(ji)) THEN
3498          ! What could not be sublimated goes into subsinksoil
3499          IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN
3500             subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji))
3501          END IF
3502          ! Sublimation is thus limited to what is available
3503          ! Then, evavpsnow is reduced, of subsinksoil
3504          subsnowveg(ji) = snow(ji)
3505          snow(ji) = zero
3506          vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
3507       ELSE
3508          snow(ji) = snow(ji) - subsnowveg(ji)
3509       ENDIF
3510       !
3511    !! 1.3. snow melt only if temperature positive
3512       !
3513       IF (temp_sol_new(ji).GT.tp_00) THEN
3514          !
3515          IF (snow(ji).GT.sneige) THEN
3516             !
3517             snowmelt(ji) = (un - frac_nobio(ji,iice))*(temp_sol_new(ji) - tp_00) * soilcap(ji) / chalfu0
3518             !
3519    !! 1.3.1 enough snow for melting or not
3520             !
3521             IF (snowmelt(ji).LT.snow(ji)) THEN
3522                snow(ji) = snow(ji) - snowmelt(ji)
3523             ELSE
3524                snowmelt(ji) = snow(ji)
3525                snow(ji) = zero
3526             END IF
3527             !
3528          ELSEIF (snow(ji).GE.zero) THEN
3529             !
3530    !! 1.3.2 not enough snow
3531             !
3532             snowmelt(ji) = snow(ji)
3533             snow(ji) = zero
3534          ELSE
3535             !
3536    !! 1.3.3 negative snow - now snow melt
3537             !
3538             snow(ji) = zero
3539             snowmelt(ji) = zero
3540             WRITE(numout,*) 'hydrol_snow: WARNING! snow was negative and was reset to zero. '
3541             !
3542          END IF
3543
3544       ENDIF
3545    !! 1.4 Snow melts above a threshold
3546       ! Ice melt only if there is more than a given mass : maxmass_snow,
3547       ! But the snow cannot melt more in one time step to what corresponds to
3548       ! a 1K cooling. This will lead to a progressive melting of snow above
3549       ! maxmass_snow but it is needed as a too strong cooling can destabilise the model.
3550       IF ( snow(ji) .GT. maxmass_snow ) THEN
3551          snow_d1k = un * soilcap(ji) / chalfu0
3552          snowmelt(ji) = snowmelt(ji) + MIN((snow(ji) - maxmass_snow),snow_d1k)
3553          snow(ji) = snow(ji) - snowmelt(ji)
3554          IF ( printlev >= 3 ) WRITE (numout,*) "Snow was above maxmass_snow (", maxmass_snow,") and we melted ", snowmelt(ji)
3555       ENDIF
3556       
3557    END DO
3558    !
3559    !! 2 On Land ice
3560    !
3561    DO ji=1,kjpindex
3562       !
3563    !! 2.1 Compute snow
3564       !
3565       !!??Aurelien: pkoi mettre precip_rain en dessous? We considere liquid precipitations becomes instantly snow? 
3566       snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
3567            & frac_nobio(ji,iice)*precip_rain(ji)
3568       !
3569    !! 2.2 Sublimation
3570       !      Was calculated before it can give us negative snow_nobio but that is OK
3571       !      Once it goes below a certain values (-maxmass_snow for instance) we should kill
3572       !      the frac_nobio(ji,iice) !
3573       !
3574       snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
3575       !
3576    !! 2.3 Snow melt only for continental ice fraction
3577       !
3578       snowmelt_tmp = zero
3579       IF (temp_sol_new(ji) .GT. tp_00) THEN
3580          !
3581    !! 2.3.1 If there is snow on the ice-fraction it can melt
3582          !
3583          snowmelt_tmp = frac_nobio(ji,iice)*(temp_sol_new(ji) - tp_00) * soilcap(ji) / chalfu0
3584          !
3585          IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN
3586             snowmelt_tmp = MAX( zero, snow_nobio(ji,iice))
3587          ENDIF
3588          snowmelt(ji) = snowmelt(ji) + snowmelt_tmp
3589          snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp
3590          !
3591       ENDIF
3592       !
3593    !! 2.4 Snow melts over a threshold
3594       !   Ice melt only if there is more than a given mass : maxmass_snow,
3595       !   But the snow cannot melt more in one time step to what corresponds to
3596       !   a 1K cooling. This will lead to a progressive melting of snow above
3597       !   maxmass_snow but it is needed as a too strong cooling can destabilise the model.
3598       !
3599       IF ( snow_nobio(ji,iice) .GT. maxmass_snow ) THEN
3600          snow_d1k = un * soilcap(ji) / chalfu0
3601          icemelt(ji) = MIN((snow_nobio(ji,iice) - maxmass_snow),snow_d1k)
3602          snow_nobio(ji,iice) = snow_nobio(ji,iice) - icemelt(ji)
3603
3604          IF ( printlev >= 3 ) WRITE (numout,*) "Snow was above maxmass_snow ON ICE (", maxmass_snow,") and we melted ", icemelt(ji)
3605       ENDIF
3606
3607    END DO
3608
3609    !
3610    !! 3 On other surface types - not done yet
3611    !
3612    IF ( nnobio .GT. 1 ) THEN
3613       WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
3614       WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
3615       CALL ipslerr_p(3,'hydrol_snow','nnobio > 1 not allowded','Cannot treat snow on these surface types.','')
3616    ENDIF
3617
3618    !
3619    !! 4 computes total melt (snow and ice)
3620    !
3621    DO ji = 1, kjpindex
3622       tot_melt(ji) = icemelt(ji) + snowmelt(ji)
3623    ENDDO
3624
3625    !
3626    !! 5 computes snow age on veg and ice (for albedo)
3627    !
3628    DO ji = 1, kjpindex
3629       !
3630    !! 5.1 Snow age on vegetation
3631       !
3632       IF (snow(ji) .LE. zero) THEN
3633          snow_age(ji) = zero
3634       ELSE
3635          snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
3636               & * EXP(-precip_snow(ji) / snow_trans)
3637       ENDIF
3638       !
3639    !! 5.2 Snow age on ice
3640       !
3641       ! age of snow on ice: a little bit different because in cold regions, we really
3642       ! cannot negect the effect of cold temperatures on snow metamorphism any more.
3643       !
3644       IF (snow_nobio(ji,iice) .LE. zero) THEN
3645          snow_nobio_age(ji,iice) = zero
3646       ELSE
3647          !
3648          d_age(ji) = ( snow_nobio_age(ji,iice) + &
3649               &  (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
3650               &  EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
3651          IF (d_age(ji) .GT. min_sechiba ) THEN
3652             xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
3653             xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
3654             d_age(ji) = d_age(ji) / (un+xx(ji))
3655          ENDIF
3656          snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
3657          !
3658       ENDIF
3659
3660    ENDDO
3661
3662    !
3663    !! 6 Diagnose the depth of the snow layer
3664    !
3665
3666    DO ji = 1, kjpindex
3667       snowdepth(ji) = snow(ji) /sn_dens
3668    ENDDO
3669
3670    IF (printlev>=3) WRITE (numout,*) ' hydrol_snow done '
3671
3672  END SUBROUTINE hydrol_snow
3673
3674   
3675!! ================================================================================================================================
3676!! SUBROUTINE   : hydrol_canop
3677!!
3678!>\BRIEF        This routine computes canopy processes.
3679!!
3680!! DESCRIPTION  :
3681!! - 1 evaporation off the continents
3682!! - 1.1 The interception loss is take off the canopy.
3683!! - 1.2 precip_rain is shared for each vegetation type
3684!! - 1.3 Limits the effect and sum what receives soil
3685!! - 1.4 swap qsintveg to the new value
3686!!
3687!! RECENT CHANGE(S) : None
3688!!
3689!! MAIN OUTPUT VARIABLE(S) :
3690!!
3691!! REFERENCE(S) :
3692!!
3693!! FLOWCHART    : None
3694!! \n
3695!_ ================================================================================================================================
3696!_ hydrol_canop
3697
3698  SUBROUTINE hydrol_canop (kjpindex, precip_rain, vevapwet, veget_max, veget, qsintmax, &
3699       & qsintveg,precisol,tot_melt)
3700
3701    !
3702    ! interface description
3703    !
3704
3705    !! 0. Variable and parameter declaration
3706
3707    !! 0.1 Input variables
3708
3709    INTEGER(i_std), INTENT(in)                               :: kjpindex    !! Domain size
3710    ! input fields
3711    REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: precip_rain !! Rain precipitation
3712    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: vevapwet    !! Interception loss
3713    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: veget_max   !! max fraction of vegetation type
3714    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: veget       !! Fraction of vegetation type
3715    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: qsintmax    !! Maximum water on vegetation for interception
3716    REAL(r_std), DIMENSION  (kjpindex), INTENT (in)          :: tot_melt    !! Total melt
3717
3718    !! 0.2 Output variables
3719
3720    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out)       :: precisol    !! Water fallen onto the ground (throughfall+Totmelt)
3721
3722    !! 0.3 Modified variables
3723
3724    REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout)     :: qsintveg    !! Water on vegetation due to interception
3725
3726    !! 0.4 Local variables
3727
3728    INTEGER(i_std)                                           :: ji, jv
3729    REAL(r_std), DIMENSION (kjpindex,nvm)                    :: zqsintvegnew
3730
3731!_ ================================================================================================================================
3732
3733    ! boucle sur les points continentaux
3734    ! calcul de qsintveg au pas de temps suivant
3735    ! par ajout du flux interception loss
3736    ! calcule par enerbil en fonction
3737    ! des calculs faits dans diffuco
3738    ! calcul de ce qui tombe sur le sol
3739    ! avec accumulation dans precisol
3740    ! essayer d'harmoniser le traitement du sol nu
3741    ! avec celui des differents types de vegetation
3742    ! fait si on impose qsintmax ( ,1) = 0.0
3743    !
3744    ! loop for continental subdomain
3745    !
3746    !
3747    !! 1 evaporation off the continents
3748    !
3749    !! 1.1 The interception loss is take off the canopy.
3750    DO jv=2,nvm
3751       qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv)
3752    END DO
3753
3754    !     It is raining :
3755    !! 1.2 precip_rain is shared for each vegetation type
3756    !
3757    qsintveg(:,1) = zero
3758    DO jv=2,nvm
3759       qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
3760    END DO
3761
3762    !
3763    !! 1.3 Limits the effect and sum what receives soil
3764    !
3765    precisol(:,1)=veget_max(:,1)*precip_rain(:)
3766    DO jv=2,nvm
3767       DO ji = 1, kjpindex
3768          zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv)) 
3769          precisol(ji,jv) = (veget(ji,jv)*throughfall_by_pft(jv)*precip_rain(ji)) + &
3770               qsintveg(ji,jv) - zqsintvegnew (ji,jv) + &
3771               (veget_max(ji,jv) - veget(ji,jv))*precip_rain(ji)
3772       ENDDO
3773    END DO
3774       
3775    ! Precisol is currently the same as throughfall, save it for diagnostics
3776    throughfall(:,:) = precisol(:,:)
3777
3778    DO jv=1,nvm
3779       DO ji = 1, kjpindex
3780          IF (vegtot(ji).GT.min_sechiba) THEN
3781             precisol(ji,jv) = precisol(ji,jv)+tot_melt(ji)*veget_max(ji,jv)/vegtot(ji)
3782          ENDIF
3783       ENDDO
3784    END DO
3785    !   
3786    !
3787    !! 1.4 swap qsintveg to the new value
3788    !
3789    DO jv=2,nvm
3790       qsintveg(:,jv) = zqsintvegnew (:,jv)
3791    END DO
3792
3793    IF (printlev>=3) WRITE (numout,*) ' hydrol_canop done '
3794
3795  END SUBROUTINE hydrol_canop
3796
3797
3798!! ================================================================================================================================
3799!! SUBROUTINE   : hydrol_vegupd
3800!!
3801!>\BRIEF        Vegetation update   
3802!!
3803!! DESCRIPTION  :
3804!!   The vegetation cover has changed and we need to adapt the reservoir distribution
3805!!   and the distribution of plants on different soil types.
3806!!   You may note that this occurs after evaporation and so on have been computed. It is
3807!!   not a problem as a new vegetation fraction will start with humrel=0 and thus will have no
3808!!   evaporation. If this is not the case it should have been caught above.
3809!!
3810!! - 1 Update of vegetation is it needed?
3811!! - 2 calculate water mass that we have to redistribute
3812!! - 3 put it into reservoir of plant whose surface area has grown
3813!! - 4 Soil tile gestion
3814!! - 5 update the corresponding masks
3815!!
3816!! RECENT CHANGE(S) : None
3817!!
3818!! MAIN OUTPUT VARIABLE(S) :
3819!!
3820!! REFERENCE(S) :
3821!!
3822!! FLOWCHART    : None
3823!! \n
3824!_ ================================================================================================================================
3825!_ hydrol_vegupd
3826
3827  SUBROUTINE hydrol_vegupd(kjpindex, veget, veget_max, soiltile, qsintveg, frac_bare, drain_upd, runoff_upd)
3828
3829
3830    !! 0. Variable and parameter declaration
3831
3832    !! 0.1 Input variables
3833
3834    ! input scalar
3835    INTEGER(i_std), INTENT(in)                            :: kjpindex 
3836    ! input fields
3837    REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in)    :: veget            !! New vegetation map
3838    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)     :: veget_max        !! Max. fraction of vegetation type