Context Navigation

hydrol.f90 @ 7444

Last change on this file since 7444 was 7444, checked in by agnes.ducharne, 20 months ago
Finishes r7443 for WETNESS_TRANSPIR_MAX.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 348.7 KB

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