Context Navigation

source: branches/publications/ORCHIDEE-Clateral/src_sechiba/hydrolc.f90 @ 7329

Last change on this file since 7329 was 3269, checked in by fabienne.maignan, 8 years ago
Correction of a bug on snowmelt with explicitsnow introduced in revision [3059]
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 197.8 KB

Line
1	! =================================================================================================================================
2	! MODULE : hydrolc
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.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 water budget of
10	!! land surfaces. It distinguishes three main reservoirs with separate water budgets :
11	!! the canopy interception reservoir, the snow pack, and the soil, described here using two
12	!! soil layers, following the so-called Choisnel scheme. The module contains the
13	!! following subroutines: hydrolc_main, hydrolc_init, hydrolc_clear, hydrolc_var_init,
14	!! hydrolc_snow, hydrolc_vegupd, hydrolc_canop, hydrolc_soil, hydrolc_waterbal, hydrolc_alma,
15	!! hydrolc_hdiff.
16	!!
17	!! \n DESCRIPTION : The aim of this module is to update the different water prognostic variables
18	!! (within canopy, snow and soil) and to compute the related liquid water fluxes (e.g.
19	!! throughfall, runoff, snow melt) and diagnostic variables used in other modules (e.g.
20	!! humrel and rsol to compute latent heat fluxes ; shumdiag for soil thermodynamics; snow mass
21	!! and snow age for snow albedo).
22	!!
23	!! The main input variables are precipitation (liquid and solid) and the different terms of
24	!! evapotranspiration (transpiration, bare soil evaporation, interception loss, sublimation).
25	!! The module organizes the required initialisation of water prognostic variables, their
26	!! integration at each timestep given the above forcings, and the required output (writing
27	!! of restart file, updated prognostic variables, diagnostic variables).
28	!!
29	!! The specificity of hydrolc in ORCHIDEE is to describe the soil using two soil layers,
30	!! following the so-called Choisnel scheme. The upper layer has a variable depth
31	!! and can disappear after dry spells (DucoudrÃ© et al., 1993). Soil moisture stress on
32	!! transpiration and bare soil evaporation acts via the height of dry soil at the top of soil.
33	!!
34	!! Runoff is produced when the entire soil column is saturated, as in the bucket scheme
35	!! of Manabe (1969). Drainage and surface runoff are only diagnosed, mostly for use in
36	!! the routing module, using a constant 95% - 5% redistribution (Ngo-Duc et al., 2005).
37	!! Irrigation can interact with the soil columns (de Rosnay et al., 2003), as well as
38	!! river discharge in floodplains (Ngo-Duc et al., 2005).
39	!!
40	!! The dynamic of the one-layer snow pack results from accumulation, sublimation and
41	!! snowmelt. Snow ageing is accounted for but it does not influence snow density.
42	!!
43	!! The subgrid variability of soil follows the one of vegetation, by introducing as many
44	!! soil columns as PFTs (de Rosnay & Polcher, 1998). The water budget is performed
45	!! independently in each PFT/soil column, but water can be exchanged laterally.
46	!! The soil moisture of the bottom layer is homogenized between the soil columns at
47	!! each timestep, and a diffusion between the top soil layers can be allowed if the
48	!! flag ok_hdiff is true (the default value is false).
49	!!
50	!! RECENT CHANGE(S) : None
51	!!
52	!! REFERENCE(S) :
53	!! - DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
54	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
55	!! Circulation Model. Journal of Climate, 6, pp. 248-273.
56	!! - Ducharne, A. Laval, K. and Polcher, J. (1998) Sensitivity of the hydrological cycle to the
57	!! parameterization of soil hydrology in a GCM. Climate Dynamics, 14:307-327.
58	!! - de Rosnay, P. and Polcher J. (1998) Modeling root water uptake in a complex land surface scheme
59	!! coupled to a GCM. Hydrology and Earth System Sciences, 2(2-3):239-256.
60	!! - de Rosnay, P., Polcher, J., Laval, K. et Sabre, M. (2003). Integrated parameterization of
61	!! irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula. Geophys. Res.
62	!! Lett, 30(19):HLS2-1 -
63	!! HLS2-4.
64	!! - Krinner, G., Viovy, N., de Noblet-DucoudrÃ©, N., Ogee, J., Polcher, J., Friedlingstein, P.,
65	!! Ciais, P., Sitch, S. et Prentice, I. (2005). A dynamic global vegetation model for studies of the
66	!! coupled atmosphere-biosphere system. Global Biogeochem. Cycles, 19(1).
67	!! - Ngo-Duc, T., Laval, K., Ramillien, G., Polcher, J. et Cazenave, A. (2007). Validation of the land
68	!! water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with
69	!! Gravity Recovery and Climate Experiment (GRACE) data. Water Resour. Res, 43:W04427.
70	!! - ALMA : http://www.lmd.jussieu.fr/~polcher/ALMA/
71	!! - DucoudrÃ©, N. (1990). Sensibilite du climat simule a la parametrisation des echanges de vapeur
72	!! d'eau entre la biosphere et l'atmosphere. These de Doctorat, UniversitÃ© Paris 6.
73	!! - Ducharne A (1997). Le cycle de l'eau : modelisation de l'hydrologie continentale, etude de ses
74	!! interactions avec le climat, These de Doctorat, UniversitÃ© Paris 6.
75	!! - de Rosnay, P. (1999). Representation des interactions sol-plante-atmosphere dans le modele de
76	!! circulation generale du LMD. These de doctorat, UniversitÃ© Paris 6.
77	!! - Guimberteau, M. (2010). Modelisation de l'hydrologie continentale et influences de l'irrigation
78	!! sur le cycle de l'eau, these de doctorat, UniversitÃ© Paris 6.
79	!!
80	!! SVN :
81	!! $HeadURL$
82	!! $Date$
83	!! $Revision$
84	!! \n
85	!_ ================================================================================================================================
86
87	MODULE hydrolc
88
89	! modules used
90	USE ioipsl
91	USE xios_orchidee
92	USE ioipsl_para
93	USE constantes ! contains technical constants
94	USE constantes_soil ! soil properties and geometry (underground)
95	USE pft_parameters ! surface and vegetation properties
96	USE sechiba_io ! for inpout/output
97	USE grid ! for grid, masks and calendar
98	USE mod_orchidee_para ! for variable and subroutines realted to the parallelism of orchidee
99	USE explicitsnow
100
101	IMPLICIT NONE
102
103	PRIVATE
104	PUBLIC :: hydrolc_main, hydrolc_initialize, hydrolc_finalize, hydrolc_clear
105
106
107	LOGICAL, SAVE :: ok_hdiff !! To do horizontal diffusion between soil columns
108	!$OMP THREADPRIVATE(ok_hdiff)
109	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: bqsb !! Water content in the bottom layer
110	!! @tex ($kg m^{-2}$) @endtex
111	!$OMP THREADPRIVATE(bqsb)
112	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gqsb !! Water content in the top layer @tex ($kg m^{-2}$) @endtex
113	!$OMP THREADPRIVATE(gqsb)
114	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsg !! Depth of the top layer (m)
115	!$OMP THREADPRIVATE(dsg)
116	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsp !! Depth of dry soil in the bottom layer plus depth of top
117	!! layer (m)
118	!$OMP THREADPRIVATE(dsp)
119	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_bqsb !! Mean water content in the bottom layer, averaged over the
120	!! soil columns @tex ($kg m^{-2}$) @endtex
121	!$OMP THREADPRIVATE(mean_bqsb)
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_gqsb !! Mean water content in the top layer, averaged over the
123	!! soil columns @tex ($kg m^{-2}$) @endtex
124	!$OMP THREADPRIVATE(mean_gqsb)
125	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_flux !! Total water flux
126	!$OMP THREADPRIVATE(tot_flux)
127	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_beg !! Total amount of water at start of time step
128	!! @tex ($kg m^{-2}$) @endtex
129	!$OMP THREADPRIVATE(tot_water_beg)
130	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_end !! Total amount of water at end of time step
131	!! @tex ($kg m^{-2}$) @endtex
132	!$OMP THREADPRIVATE(tot_water_end)
133	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_beg !! Total amount of intercepted water on vegetation at start of
134	!! time step @tex ($kg m^{-2}$) @endtex
135	!$OMP THREADPRIVATE(tot_watveg_beg)
136	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_end !! Total amount of intercepted water on vegetation at end of
137	!! time step @tex ($kg m^{-2}$) @endtex
138	!$OMP THREADPRIVATE(tot_watveg_end)
139	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_beg !! Total amount of water in the soil at start of time step
140	!! @tex ($kg m^{-2}$) @endtex
141	!$OMP THREADPRIVATE(tot_watsoil_beg)
142	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_end !! Total amount of water in the soil at end of time step
143	!! @tex ($kg m^{-2}$) @endtex
144	!$OMP THREADPRIVATE(tot_watsoil_end)
145	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_beg !! Total snow water equivalent at start of time step
146	!! @tex ($kg m^{-2}$) @endtex
147	!$OMP THREADPRIVATE(snow_beg)
148	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_end !! Total snow water equivalent at end of time step
149	!! @tex ($kg m^{-2}$) @endtex
150	!$OMP THREADPRIVATE(snow_end)
151	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delsoilmoist !! Change in soil moisture during time step
152	!! @tex ($kg m^{-2}$) @endtex
153	!$OMP THREADPRIVATE(delsoilmoist)
154	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delintercept !! Change in interception storage during time step
155	!! @tex ($kg m^{-2}$) @endtex
156	!$OMP THREADPRIVATE(delintercept)
157	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delswe !! Change in snow water equivalent during time step
158	!! @tex ($kg m^{-2}$) @endtex
159	!$OMP THREADPRIVATE(delswe)
160	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dss !! Depth of dry soil at the top, whether in the top or bottom
161	!! layer (m)
162	!$OMP THREADPRIVATE(dss)
163	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: hdry !! Mean top dry soil height (m) version beton
164	!$OMP THREADPRIVATE(hdry)
165	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Throughfall @tex ($kg m^{-2}$) @endtex
166	!$OMP THREADPRIVATE(precisol)
167	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsnowveg !! Sublimation of snow on vegetation
168	!! @tex ($kg m^{-2}$) @endtex
169	!$OMP THREADPRIVATE(subsnowveg)
170	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes,...)
171	!! @tex ($kg m^{-2}$) @endtex
172	!$OMP THREADPRIVATE(subsnownobio)
173	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snowmelt !! Snow melt @tex ($kg m^{-2}$) @endtex
174	!$OMP THREADPRIVATE(snowmelt)
175	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: icemelt !! Ice melt @tex ($kg m^{-2}$) @endtex
176	!$OMP THREADPRIVATE(icemelt)
177	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gdrainage !! Drainage between the two soil layers
178	!! @tex ($kg m^{-2}$) @endtex
179	!$OMP THREADPRIVATE(gdrainage)
180	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vegtot !! Total fraction of grid-cell covered by PFTs (bare soil +
181	!! vegetation) (0-1, unitless)
182	!$OMP THREADPRIVATE(vegtot)
183	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: resdist !! Previous values of "veget" (0-1, unitless) ! Last map of
184	!! PFT fractions used to define the soil columns
185	!$OMP THREADPRIVATE(resdist)
186	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var !! Maximum water content of the soil
187	!! @tex ($kg m^{-2}$) @endtex
188	!$OMP THREADPRIVATE(mx_eau_var)
189	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ruu_ch !! Volumetric soil water capacity @tex ($kg m^{-3}$) @endtex
190	!$OMP THREADPRIVATE(ruu_ch)
191	REAL(r_std), PARAMETER :: dsg_min = 0.001 !! Reference depth to define subgrid variability of saturated
192	!! top soil layer (m)
193	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsinksoil !! Excess of sublimation as a sink for the soil
194	!$OMP THREADPRIVATE(subsinksoil)
195
196	!_ ================================================================================================================================
197
198	CONTAINS
199
200	!! ================================================================================================================================
201	!! SUBROUTINE : hydrolc_initialize
202	!!
203	!>\BRIEF Allocate module variables, read from restart file or initialize with default values
204	!!
205	!! DESCRIPTION :
206	!!
207	!! MAIN OUTPUT VARIABLE(S) :
208	!!
209	!! REFERENCE(S) :
210	!!
211	!! FLOWCHART : None
212	!! \n
213	!_ ================================================================================================================================
214
215	SUBROUTINE hydrolc_initialize( kjit, kjpindex, index, rest_id, &
216	veget, veget_max, tot_bare_soil, &
217	rsol, drysoil_frac, snow, &
218	snow_age, snow_nobio, snow_nobio_age, humrel, &
219	vegstress, qsintveg, shumdiag, snowrho, &
220	snowtemp, snowgrain, snowdz, snowheat )
221
222	!! 0. Variable and parameter declaration
223
224	!! 0.1 Input variables
225	INTEGER(i_std), INTENT(in) :: kjit !! Current time step number (unitless)
226	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
227	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the land grid points on the map
228	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
229	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by
230	!! vegetation for each PFT, except for soil
231	!! (0-1, unitless)
232	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! PFT fractions within grid-cells: max value of
233	!! veget for vegetation PFTs and min value for
234	!! bare soil (0-1, unitless)
235	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
236
237	!! 0.2 Output variables
238	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistance to bare soil evaporation
239	!! @tex ($s m^{-1}$) @endtex
240	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visibly dry soil, for bare soil
241	!! albedo calculation in condveg.f90
242	!! (0-1, unitless)
243	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
244	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age (days)
245	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio !! Snow water equivalent on nobio areas
246	!! @tex ($kg m^{-2}$) @endtex
247	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age !! Snow age on ice, lakes, ... (days)
248	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Soil moisture stress factor on transpiration and
249
250	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Vegetation moisture stress (only for vegetation
251	!! growth) (0-1, unitless)
252	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Amount of water in the canopy interception
253	!! reservoir @tex ($kg m^{-2}$) @endtex
254	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! Mean relative soil moisture in the different
255	!! levels used by thermosoil.f90 (0-1, unitless)
256	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowrho !! Snow density
257	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowtemp !! Snow temperature
258	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowgrain !! Snow grainsize
259	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowdz !! Snow layer thickness
260	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowheat !! Snow heat content
261
262	!! 0.4 Local variables
263
264	!_ ================================================================================================================================
265
266	!! 1. Initialisation
267	CALL hydrolc_init (kjit, kjpindex, index, rest_id, &
268	veget, tot_bare_soil, &
269	humrel, vegstress, shumdiag, &
270	snow, snow_age, snow_nobio, snow_nobio_age, qsintveg, &
271	snowdz, snowgrain, snowrho, snowtemp,snowheat, &
272	drysoil_frac, rsol)
273
274	CALL hydrolc_var_init (kjpindex, veget, veget_max, tot_bare_soil, &
275	rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag)
276
277	! If we check the water balance, we first save the total amount of water
278	IF (check_waterbal) CALL hydrolc_waterbal_init (kjpindex, qsintveg, snow, snow_nobio)
279
280	IF (almaoutput) THEN
281	CALL hydrolc_alma_init(kjpindex, index, qsintveg, snow, snow_nobio)
282	ENDIF
283
284	END SUBROUTINE hydrolc_initialize
285
286
287	!! ================================================================================================================================
288	!! SUBROUTINE : hydrolc_main
289	!!
290	!>\BRIEF Main routine for water budget calculation on land surfaces.
291	!!
292	!! DESCRIPTION : This routine drives all the calls pertaining
293	!! to the land-surface water budget. It is called once during the initialisation of
294	!! ORCHIDEE, and then once for each time step, to drive the water budget calculations,
295	!! and the link to history files processing (histwrite). It is called one more time
296	!! at the last time step, for the writing of the restart file.
297	!!
298	!! The initialisation step calls hydrolc_init and hydrolc_var_init, plus hydrolc_waterbal
299	!! if we choose to check the water balance. Writing is performed to the restart file,
300	!! before the water budget calculations, at a timestep that is controlled by the flag
301	!! "ldrestart_write".
302	!!
303	!! The water budget calculations are separated into four subroutines related to the
304	!! main three water reservoirs over land grid-cells, and called at each time step : \n
305	!! - hydrolc_snow for snow process (including age of snow)
306	!! - hydrolc_vegupd to manage the required redistribution of water between reservoirs
307	!! when the vegetation fraction of PFTs has changed
308	!! - hydrolc_canop for canopy process (interception)
309	!! - hydrolc_soil for soil hydrological process (soil moisture and runoff), followed by
310	!! hydrolc_hdiff if we choose to account for horizontal diffusion between the soil columns.
311	!! If we choose to check the water balance, this is done over each timestep dt_sechiba,
312	!! by hydrolc_waterbal.
313	!!
314	!! The link to "hist" files processing has two kinds of controls:
315	!! - the standard output ncdf file corresponds to hist_id ; if hist2_id positive, a second
316	!! ncdf file is output, with a different frequency
317	!! - almaoutput defines the nature/name of the written variables, whether following the
318	!! ALMA norm or not
319	!! Note that histwrite does different actions depending on the current time step :
320	!! writing on "hist" files at a selected frequency ; summing of the output variables
321	!! in between in view of averaging.
322	!! The subroutine also handles some "CMIP5" output (mrro, mrros and prveg).
323	!!
324	!! We consider that any water on ice (nobio) is snow and we only peform a water balance to have consistency.
325	!! The water balance is limited to + or - \f$10^6\f$ so that accumulation is not endless
326	!!
327	!! IMPORTANT NOTE : The water fluxes are used in their integrated form, over the time step
328	!! dt_sechiba, with a unit of \f$kg.m^{-2}\f$.
329	!!
330	!! RECENT CHANGE(S) : None
331	!!
332	!! MAIN OUTPUT VARIABLE(S) : snow, snow_age, snow_nobio, snow_nobio_age, tot_melt,
333	!! returnflow, irrigation, humrel, vegstress, rsol, drysoil_frac, shumdiag, litterhumdiag
334	!! + variables declared in the module
335	!! + variables sent for ouput by histwrite
336	!!
337	!! REFERENCE(S) : Same as for module hydrolc
338	!!
339	!! FLOWCHART : None
340	!! \n
341	!_ ================================================================================================================================
342
343	SUBROUTINE hydrolc_main (kjit, kjpindex, index, indexveg, &
344	& temp_sol_new, floodout, run_off_tot, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
345	& qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, &
346	& precip_rain, precip_snow, returnflow, reinfiltration, irrigation, humrel, vegstress, rsol, drysoil_frac, &
347	& evapot, evapot_corr, flood_frac, flood_res, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id, &
348	& temp_air, pb, u, v, swnet, pgflux, &
349	& snowrho,snowtemp,snowgrain,snowdz,snowheat,snowliq, &
350	& grndflux,gtemp,tot_bare_soil, &
351	& lambda_snow,cgrnd_snow,dgrnd_snow,temp_sol_add)
352
353	!! 0. Variable and parameter declaration
354
355	!! 0.1 Input variables
356
357	INTEGER(i_std), INTENT(in) :: kjit !! Current time step number (unitless)
358	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
359	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
360	!! (unitless)
361	INTEGER(i_std),INTENT (in) :: hist2_id !! Second _history_ file identifier (unitless)
362	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the land grid points on the map
363	!! (unitless)
364	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indices of the PFT tiles / soil columns on
365	!! the 3D map (unitless)
366	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rainfall @tex ($kg m^{-2}$) @endtex
367	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snowfall @tex ($kg m^{-2}$) @endtex
368	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil
369	!! @tex ($kg m^{-2}$) @endtex
370	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: reinfiltration !! Routed water which comes back into the soil
371	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Irrigation water applied to soils
372	!! @tex ($kg m^{-2}$) @endtex
373	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature (K)
374	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of terrestrial ice, lakes, ...
375	!! (0-1, unitless)
376	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of terrestrial ice+lakes+...
377	!! (0-1, unitless)
378	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil heat capacity @tex ($J K^{-1}$) @endtex
379	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss over each PFT
380	!! @tex ($kg m^{-2}$) @endtex
381	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by
382	!! vegetation for each PFT, except for soil
383	!! (0-1, unitless)
384	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! PFT fractions within grid-cells: max value of
385	!! veget for vegetation PFTs and min value for bare soil (0-1, unitless)
386	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum amount of water in the canopy
387	!! interception reservoir
388	!! @tex ($kg m^{-2}$) @endtex
389	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration over each PFT
390	!! @tex ($kg m^{-2}$) @endtex
391	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! [DISPENSABLE] Soil Potential Evaporation
392	!! @tex ($kg m^{-2}$) @endtex
393	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! [DISPENSABLE] Soil Potential Evaporation
394	!! Correction
395	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: flood_frac !! Flooded fraction
396	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: temp_air !! Air temperature
397	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: u,v !! Horizontal wind speed
398	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: pb !! Surface pressure
399	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: swnet !! Net short wave radiation
400	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: pgflux !! Net energy into snowpack
401	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: gtemp !! First soil layer temperature
402	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
403	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: lambda_snow !! Coefficient of the linear extrapolation of surface temperature
404	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (in) :: cgrnd_snow !! Integration coefficient for snow numerical scheme
405	REAL(r_std),DIMENSION (kjpindex,nsnow), INTENT (in) :: dgrnd_snow !! Integration coefficient for snow numerical scheme
406
407	!! 0.2 Output variables
408
409	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: run_off_tot !! Diagnosed surface runoff
410	!! @tex ($kg m^{-2}$) @endtex
411	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: flood_res !! Flood reservoir estimate
412	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drainage !! Diagnosed rainage at the bottom of soil
413	!! @tex ($kg m^{-2}$) @endtex
414	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistance to bare soil evaporation
415	!! @tex ($s m^{-1}$) @endtex
416	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visibly dry soil, for bare soil
417	!! albedo calculation in condveg.f90
418	!! (0-1, unitless)
419	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! Litter humidity factor (0-1, unitless), used
420	!! in stomate
421	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt @tex ($kg m^{-2}$) @endtex
422
423
424	!! 0.3 Modified variables
425
426	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapnu !! Bare soil evaporation @tex ($kg m^{-2}$) @endtex
427	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
428	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapflo !! Floodplains evaporation
429	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
430	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow_age !! Snow age (days)
431	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio !! Snow water equivalent on nobio areas
432	!! @tex ($kg m^{-2}$) @endtex
433	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ... (days)
434
435	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Soil moisture stress factor on transpiration and
436	!! bare soil evaporation (0-1, unitless)
437	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vegstress !! Vegetation moisture stress (only for vegetation
438	!! growth) (0-1, unitless)
439	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Amount of water in the canopy interception
440	!! reservoir @tex ($kg m^{-2}$) @endtex
441	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: floodout !! flux out of floodplains
442	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout) :: shumdiag !! Mean relative soil moisture in the different
443	!! levels used by thermosoil.f90 (0-1, unitless)
444	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowrho !! Snow density
445	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowtemp !! Snow temperature
446	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowgrain !! Snow grainsize
447	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowdz !! Snow layer thickness
448	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(inout) :: snowheat !! Snow heat content
449	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(out) :: snowliq !! Liquid water content (m)
450	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: grndflux !! Net flux into soil W/m2
451	REAL(r_std), DIMENSION (kjpindex), INTENT (inout) :: temp_sol_add !! additional surface temperature due to the melt of first layer
452	!! at the present time-step @tex ($K$) @endtex
453
454
455	!! 0.4 Local variables
456
457	REAL(r_std),DIMENSION (kjpindex) :: soilwet !! Temporary diagnostic of relative humidity of
458	!! total soil (0-1, unitless)
459	REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Snow Depth (m)
460	INTEGER(i_std) :: ji,jv !! Grid-cell and PFT indices (unitless)
461	REAL(r_std), DIMENSION(kjpindex) :: histvar !! Ancillary variable when computations are needed
462	!! before writing to history files
463	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
464
465	!_ ================================================================================================================================
466
467	!! 1. Water budget calculations:
468
469	!! 1.a snow processes
470	IF (ok_explicitsnow) THEN
471
472	IF (printlev>=3) WRITE (numout,*) ' ok_explicitsnow : use multi-snow layer '
473	CALL explicitsnow_main(kjpindex, precip_rain, precip_snow, temp_air, pb, &
474	u, v, temp_sol_new, soilcap, pgflux, &
475	frac_nobio, totfrac_nobio,gtemp, &
476	lambda_snow, cgrnd_snow, dgrnd_snow, &
477	vevapsno, snow_age, snow_nobio_age,snow_nobio, snowrho, &
478	snowgrain, snowdz, snowtemp, snowheat, snow, &
479	temp_sol_add, &
480	snowliq, subsnownobio, grndflux, snowmelt, tot_melt, &
481	subsinksoil)
482	ELSE
483	CALL hydrolc_snow(kjpindex, precip_rain, precip_snow, temp_sol_new, soilcap, &
484	frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
485	tot_melt, snowdepth)
486	END IF
487
488	!! 1.b canopy processes
489	CALL hydrolc_vegupd(kjpindex, veget, tot_bare_soil, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss,dsp, resdist)
490
491	CALL hydrolc_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, tot_bare_soil, qsintveg, precisol)
492	!
493	! computes surface reservoir
494	!
495	CALL hydrolc_flood(kjpindex, vevapnu, vevapflo, flood_frac, flood_res, floodout)
496
497	!! 1.c soil hydrology processes
498	CALL hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, reinfiltration, irrigation, tot_melt, mx_eau_var, &
499	& veget, tot_bare_soil, ruu_ch, transpir,&
500	& gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, run_off_tot, drainage, humrel, &
501	& vegstress, shumdiag, litterhumdiag)
502
503	! computes horizontal diffusion between the water reservoirs
504	IF ( ok_hdiff ) THEN
505	CALL hydrolc_hdiff(kjpindex, veget, tot_bare_soil, ruu_ch, gqsb, bqsb, dsg, dss, dsp)
506	ENDIF
507
508	! If we check the water balance, we end with the comparison of total water change and fluxes
509	IF (check_waterbal) THEN
510	CALL hydrolc_waterbal(kjpindex, index, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
511	& precip_rain, precip_snow, returnflow, reinfiltration, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,&
512	& vevapflo, floodout, run_off_tot, drainage )
513	ENDIF
514
515	!! 2. Output
516
517	! If we use the ALMA standards (then 'humrel' is not written)
518	IF (almaoutput) THEN
519	CALL hydrolc_alma(kjpindex, index, qsintveg, snow, snow_nobio, soilwet)
520	ENDIF
521
522	CALL xios_orchidee_send_field("runoff",run_off_tot*one_day/dt_sechiba)
523	CALL xios_orchidee_send_field("drainage",drainage*one_day/dt_sechiba)
524	CALL xios_orchidee_send_field("rain",precip_rain*one_day/dt_sechiba)
525	CALL xios_orchidee_send_field("snowf",precip_snow*one_day/dt_sechiba)
526	CALL xios_orchidee_send_field("rain_alma",precip_rain/dt_sechiba)
527	CALL xios_orchidee_send_field("snowf_alma",precip_snow/dt_sechiba)
528	CALL xios_orchidee_send_field("dss",dss)
529	CALL xios_orchidee_send_field("bqsb",mean_bqsb)
530	CALL xios_orchidee_send_field("gqsb",mean_gqsb)
531	CALL xios_orchidee_send_field("qsintveg",qsintveg)
532	CALL xios_orchidee_send_field("precisol",precisol*one_day/dt_sechiba)
533	CALL xios_orchidee_send_field("floodout",floodout*one_day/dt_sechiba)
534	CALL xios_orchidee_send_field("humrel",humrel)
535	CALL xios_orchidee_send_field("qsintmax",qsintmax)
536	CALL xios_orchidee_send_field("mrros",run_off_tot/dt_sechiba)
537
538	histvar(:)=run_off_tot(:)+drainage(:)
539	CALL xios_orchidee_send_field("mrro",histvar/dt_sechiba)
540
541	histvar(:)=(precip_rain(:)-SUM(precisol(:,:),dim=2))
542	CALL xios_orchidee_send_field("prveg",histvar/dt_sechiba)
543
544	histvar(:)=zero
545	DO jv = 1, nvm
546	DO ji = 1, kjpindex
547	IF ( vegtot(ji) .GT. zero ) THEN
548	!MM veget(:,1) BUG ??!!!!!!!!!!!
549	IF (jv .EQ. 1) THEN
550	histvar(ji)=histvar(ji)+tot_bare_soil(ji)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
551	ELSE
552	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
553	ENDIF
554	ENDIF
555	ENDDO
556	ENDDO
557	CALL xios_orchidee_send_field("mrsos",histvar)
558
559	! Total soil moisture
560	histvar(:)=mean_bqsb(:)+mean_gqsb(:)
561	CALL xios_orchidee_send_field("mrso",histvar)
562	CALL xios_orchidee_send_field("SubSnow",vevapsno/dt_sechiba)
563	CALL xios_orchidee_send_field("SnowDepth",snowdepth)
564	CALL xios_orchidee_send_field("Qs",run_off_tot/dt_sechiba)
565	CALL xios_orchidee_send_field("Qsb",drainage/dt_sechiba)
566	CALL xios_orchidee_send_field("Qsm",tot_melt/dt_sechiba)
567
568	IF (almaoutput) THEN
569	CALL xios_orchidee_send_field("SoilWet",soilwet)
570	CALL xios_orchidee_send_field("DelSoilMoist",delsoilmoist)
571	CALL xios_orchidee_send_field("DelSWE",delswe)
572	CALL xios_orchidee_send_field("DelIntercept",delintercept)
573	CALL xios_orchidee_send_field("RootMoist",tot_watsoil_end)
574	END IF
575
576
577	IF ( .NOT. almaoutput ) THEN
578	CALL histwrite_p(hist_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
579	CALL histwrite_p(hist_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
580	CALL histwrite_p(hist_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
581	CALL histwrite_p(hist_id, 'runoff', kjit, run_off_tot, kjpindex, index) ! this is surface runoff = 5% of total runoff
582	CALL histwrite_p(hist_id, 'drainage', kjit, drainage, kjpindex, index) ! this 95% of total runoff
583	IF ( do_floodplains ) THEN
584	CALL histwrite_p(hist_id, 'floodout', kjit, floodout, kjpindex, index)
585	ENDIF
586	CALL histwrite_p(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
587	CALL histwrite_p(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
588	CALL histwrite_p(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
589	CALL histwrite_p(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
590	CALL histwrite_p(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
591	CALL histwrite_p(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg) ! this the transpiration stress factor
592
593	!! The output for "CMIP5" is handled here, assuming that fluxes in kg m^{-2} s^{-1}
594	!! But the transformation below on mrro, mrros and prveg lead to fluxes that are 48 times too small !***
595
596	histvar(:)=zero
597	DO jv = 1, nvm
598	DO ji = 1, kjpindex
599	IF ( vegtot(ji) .GT. zero ) THEN
600	!MM veget(:,1) BUG ??!!!!!!!!!!!
601	IF (jv .EQ. 1) THEN
602	histvar(ji)=histvar(ji)+tot_bare_soil(ji)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
603	ELSE
604	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
605	ENDIF
606	ENDIF
607	ENDDO
608	ENDDO
609	! this is soil moisture in the top 10 cms
610	CALL histwrite_p(hist_id, 'mrsos', kjit, histvar, kjpindex, index)
611
612	! this is total soil moisture
613	histvar(:)=mean_bqsb(:)+mean_gqsb(:)
614	CALL histwrite_p(hist_id, 'mrso', kjit, histvar, kjpindex, index)
615
616	CALL histwrite_p(hist_id, 'mrros', kjit, run_off_tot, kjpindex, index)
617
618	histvar(:)=run_off_tot(:)+drainage(:)
619	CALL histwrite_p(hist_id, 'mrro', kjit, histvar, kjpindex, index)
620
621	histvar(:)=(precip_rain(:)-SUM(precisol(:,:),dim=2))
622	CALL histwrite_p(hist_id, 'prveg', kjit, histvar, kjpindex, index)
623	CALL histwrite_p(hist_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
624
625	ELSE
626
627	! almaoutput=true
628	CALL histwrite_p(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
629	CALL histwrite_p(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
630
631	! surface runoff = 5% of total runoff
632	CALL histwrite_p(hist_id, 'Qs', kjit, run_off_tot, kjpindex, index)
633
634	! drainage = 95% of total runoff
635	CALL histwrite_p(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
636	CALL histwrite_p(hist_id, 'Qsm', kjit, snowmelt, kjpindex, index)
637	CALL histwrite_p(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
638	CALL histwrite_p(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
639	CALL histwrite_p(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
640	CALL histwrite_p(hist_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
641	CALL histwrite_p(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
642
643	! this the transpiration stress factor
644	CALL histwrite_p(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
645
646	CALL histwrite_p(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
647	CALL histwrite_p(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
648	CALL histwrite_p(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
649
650	ENDIF
651	IF ( hist2_id > 0 ) THEN
652	IF ( .NOT. almaoutput ) THEN
653	CALL histwrite_p(hist2_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
654	CALL histwrite_p(hist2_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
655	CALL histwrite_p(hist2_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
656	CALL histwrite_p(hist2_id, 'runoff', kjit, run_off_tot, kjpindex, index)
657	CALL histwrite_p(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
658	IF ( do_floodplains ) THEN
659	CALL histwrite_p(hist2_id, 'floodout', kjit, floodout, kjpindex, index)
660	ENDIF
661	CALL histwrite_p(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
662	CALL histwrite_p(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
663	CALL histwrite_p(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
664	CALL histwrite_p(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
665	CALL histwrite_p(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
666	CALL histwrite_p(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
667	CALL histwrite_p(hist2_id, 'snowmelt',kjit,snowmelt,kjpindex,index)
668
669	IF (check_waterbal) THEN
670	CALL histwrite_p(hist2_id, 'TotWater', kjit, tot_water_end, kjpindex, index)
671	CALL histwrite_p(hist2_id, 'TotWaterFlux', kjit, tot_flux, kjpindex, index)
672	ENDIF
673
674	histvar(:)=zero
675	DO jv = 1, nvm
676	DO ji = 1, kjpindex
677	IF ( vegtot(ji) .GT. zero ) THEN
678	!MM veget(:,1) BUG ??!!!!!!!!!!!
679	IF (jv .EQ. 1) THEN
680	histvar(ji)=histvar(ji)+tot_bare_soil(ji)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
681	ELSE
682	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))wmax_veg(jv), zero)
683	ENDIF
684	ENDIF
685	ENDDO
686	ENDDO
687	CALL histwrite_p(hist2_id, 'mrsos', kjit, histvar, kjpindex, index)
688
689	histvar(:)=run_off_tot(:)+drainage(:)
690	CALL histwrite_p(hist2_id, 'mrro', kjit, histvar, kjpindex, index)
691
692
693	! this is total soil moisture
694	histvar(:)=mean_bqsb(:)+mean_gqsb(:)
695	CALL histwrite_p(hist2_id, 'mrso', kjit, histvar, kjpindex, index)
696	CALL histwrite_p(hist2_id, 'mrros', kjit, run_off_tot, kjpindex, index)
697
698	ELSE
699	CALL histwrite_p(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
700	CALL histwrite_p(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
701	CALL histwrite_p(hist2_id, 'Qs', kjit, run_off_tot, kjpindex, index)
702	CALL histwrite_p(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
703	CALL histwrite_p(hist2_id, 'Qsm', kjit, snowmelt, kjpindex, index)
704	CALL histwrite_p(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
705	CALL histwrite_p(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
706	CALL histwrite_p(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
707	CALL histwrite_p(hist2_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
708	CALL histwrite_p(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
709
710	CALL histwrite_p(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
711	CALL histwrite_p(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
712
713	CALL histwrite_p(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
714
715	ENDIF
716	ENDIF
717
718	IF (printlev>=3) WRITE (numout,*) ' hydrolc_main Done '
719
720	END SUBROUTINE hydrolc_main
721
722
723	!! ================================================================================================================================
724	!! SUBROUTINE : hydrolc_finalize
725	!!
726	!>\BRIEF
727	!!
728	!! DESCRIPTION : This subroutine writes the module variables and variables calculated in hydrolc to restart file
729	!!
730	!! MAIN OUTPUT VARIABLE(S) :
731	!!
732	!! REFERENCE(S) :
733	!!
734	!! FLOWCHART : None
735	!! \n
736	!_ ================================================================================================================================
737
738	SUBROUTINE hydrolc_finalize( kjit, kjpindex, rest_id, snow, &
739	snow_age, snow_nobio, snow_nobio_age, humrel, &
740	vegstress, qsintveg, snowrho, snowtemp, &
741	snowdz, snowheat, snowgrain, &
742	drysoil_frac, rsol, shumdiag)
743
744
745	!! 0. Variable and parameter declaration
746	!! 0.1 Input variables
747	INTEGER(i_std), INTENT(in) :: kjit !! Current time step number (unitless)
748	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
749	INTEGER(i_std),INTENT (in) :: rest_id !! Restart file identifier
750	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent
751	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow_age !! Snow age (days)
752	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Snow water equivalent on nobio areas
753	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio_age !! Snow age on ice, lakes, ... (days)
754	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Soil moisture stress factor on transpiration and
755	!! bare soil evaporation (0-1, unitless)
756	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vegstress !! Vegetation moisture stress (only for vegetation
757	!! growth) (0-1, unitless)
758	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Amount of water in the canopy interception
759	!! reservoir @tex ($kg m^{-2}$) @endtex
760	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowrho !! Snow density
761	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowtemp !! Snow temperature
762	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowdz !! Snow layer thickness
763	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowheat !! Snow heat content
764	REAL(r_std), DIMENSION (kjpindex,nsnow), INTENT(in) :: snowgrain !! Snow grain size
765	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: drysoil_frac
766	REAL(r_std),DIMENSION (kjpindex),INTENT(in) :: rsol
767	REAL(r_std),DIMENSION (kjpindex,nbdl),INTENT(in) :: shumdiag
768	!_ ================================================================================================================================
769
770	IF (printlev>=3) WRITE (numout,*) ' we have to complete restart file with HYDROLOGIC variables '
771
772	CALL restput_p(rest_id, 'humrel', nbp_glo, nvm, 1, kjit, humrel, 'scatter', nbp_glo, index_g)
773	CALL restput_p(rest_id, 'vegstress', nbp_glo, nvm, 1, kjit, vegstress, 'scatter', nbp_glo, index_g)
774	CALL restput_p(rest_id, 'snow', nbp_glo, 1, 1, kjit, snow, 'scatter', nbp_glo, index_g)
775	CALL restput_p(rest_id, 'snow_age', nbp_glo, 1, 1, kjit, snow_age, 'scatter', nbp_glo, index_g)
776	CALL restput_p(rest_id, 'snow_nobio', nbp_glo, nnobio, 1, kjit, snow_nobio, 'scatter', nbp_glo, index_g)
777	CALL restput_p(rest_id, 'snow_nobio_age', nbp_glo, nnobio, 1, kjit, snow_nobio_age, 'scatter', nbp_glo, index_g)
778	CALL restput_p(rest_id, 'bqsb', nbp_glo, nvm, 1, kjit, bqsb, 'scatter', nbp_glo, index_g)
779	CALL restput_p(rest_id, 'gqsb', nbp_glo, nvm, 1, kjit, gqsb, 'scatter', nbp_glo, index_g)
780	CALL restput_p(rest_id, 'dsg', nbp_glo, nvm, 1, kjit, dsg, 'scatter', nbp_glo, index_g)
781	CALL restput_p(rest_id, 'dsp', nbp_glo, nvm, 1, kjit, dsp, 'scatter', nbp_glo, index_g)
782	CALL restput_p(rest_id, 'dss', nbp_glo, nvm, 1, kjit, dss, 'scatter', nbp_glo, index_g)
783	CALL restput_p(rest_id, 'hdry', nbp_glo, 1, 1, kjit, hdry, 'scatter', nbp_glo, index_g)
784	CALL restput_p(rest_id, 'qsintveg', nbp_glo, nvm, 1, kjit, qsintveg, 'scatter', nbp_glo, index_g)
785	CALL restput_p(rest_id, 'resdist', nbp_glo, nvm, 1, kjit, resdist, 'scatter', nbp_glo, index_g)
786	CALL restput_p(rest_id, 'drysoil_frac', nbp_glo, 1, 1, kjit, drysoil_frac, 'scatter', nbp_glo, index_g)
787	CALL restput_p(rest_id, 'rsol', nbp_glo, 1, 1, kjit, rsol, 'scatter', nbp_glo, index_g)
788	CALL restput_p(rest_id, 'shumdiag', nbp_glo, nbdl, 1, kjit, shumdiag, 'scatter', nbp_glo, index_g)
789	CALL restput_p(rest_id, 'mean_gqsb', nbp_glo, 1, 1, kjit, mean_gqsb, 'scatter', nbp_glo, index_g)
790	CALL restput_p(rest_id, 'mean_bqsb', nbp_glo, 1, 1, kjit, mean_bqsb, 'scatter', nbp_glo, index_g)
791	CALL restput_p(rest_id, 'mx_eau_var', nbp_glo, 1, 1, kjit, mx_eau_var, 'scatter', nbp_glo, index_g)
792	CALL restput_p(rest_id, 'vegtot', nbp_glo, 1, 1, kjit, vegtot, 'scatter', nbp_glo, index_g)
793
794
795	! Write variables for explictsnow module to restart file
796	IF (ok_explicitsnow) THEN
797	CALL explicitsnow_finalize ( kjit, kjpindex, rest_id, snowrho, &
798	snowtemp, snowdz, snowheat, snowgrain)
799	END IF
800
801	END SUBROUTINE hydrolc_finalize
802
803
804	!! ================================================================================================================================
805	!! SUBROUTINE : hydrolc_init
806	!!
807	!>\BRIEF This routine drives the initialisation of the water budget calculations.
808	!!
809	!! DESCRIPTION : The following sequences are performed :
810	!! - Setting ok_hdiff (default is false) for horizontal diffusion between soil columns
811	!! - Dynamic allocation of arrays that are local to module hydrolc
812	!! - Initializing prognostic variables from input restart file
813	!! - Default attribution is performed in case the model is started without a restart file
814	!!
815	!! RECENT CHANGE(S) : None
816	!!
817	!! MAIN OUTPUT VARIABLE(S) : humrel, vegstress,snow, snow_age, snow_nobio, snow_nobio_age, qsintveg
818	!!
819	!! REFERENCE(S) : None
820	!!
821	!! FLOWCHART11 : None
822	!! \n
823	!_ ================================================================================================================================
824
825	SUBROUTINE hydrolc_init(kjit, kjpindex, index, rest_id, &
826	veget, tot_bare_soil, &
827	humrel, vegstress, shumdiag, &
828	snow, snow_age, snow_nobio, snow_nobio_age, qsintveg, &
829	snowdz, snowgrain, snowrho, snowtemp, snowheat, &
830	drysoil_frac, rsol)
831
832	!! 0. Variable and parameter declaration
833
834	!! 0.1 Input variables
835
836	INTEGER(i_std), INTENT (in) :: kjit !! Current time step number (unitless)
837	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
838	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the land grid points on the map
839	!! (unitless)
840	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier (unitless)
841	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by
842	!! vegetation for each PFT, except for soil
843	!! (0-1, unitless)
844	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
845
846	!! 0.2 Output variables
847
848	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Soil moisture stress factor on transpiration
849	!! and bare soil evaporation (0-1, unitless)
850	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Vegetation moisture stress (only for
851	!! vegetation growth) (0-1, unitless)
852	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! Mean relative soil moisture in the different
853	!! levels used by thermosoil.f90 (0-1, unitless)
854	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
855	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age (days)
856	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out):: snow_nobio !! Snow water equivalent on nobio areas
857	!! @tex ($kg m^{-2}$) @endtex
858	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out):: snow_nobio_age !! Snow age on ice, lakes, ... (days)
859	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Amount of water in the canopy interception
860	!! reservoir @tex ($kg m^{-2}$) @endtex
861	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out) :: snowdz !! Snow depth
862	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out) :: snowgrain !! Snow grain size
863	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out) :: snowheat !! Snow heat content
864	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out) :: snowtemp !! Snow temperature
865	REAL(r_std),DIMENSION (kjpindex,nsnow),INTENT(out) :: snowrho !! Snow density
866	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: drysoil_frac
867	REAL(r_std),DIMENSION (kjpindex),INTENT(out) :: rsol !! Resistance to bare soil evaporation
868
869	!! 0.3 Modified variables
870
871	!! 0.4 Local variables
872
873	INTEGER(i_std) :: ier !! To receive error flag during allocation
874	INTEGER(i_std) :: ji,jv,ik !! Indices for grid-cells, PFTs, and grid-cells
875	!! (unitless)
876	REAL(r_std), DIMENSION (kjpindex,nvm) :: zdsp, tmpdss !! Ancillary variables for initializing dsp
877	!! and dss (m)
878	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: dsp_g !! Ancillary variable for initializing dsp (m)
879	!! dss (m)
880	REAL(r_std), DIMENSION(kjpindex) :: a_subgrd !! Diagnosed subgrid fraction of saturated soil in
881	!! the top layer, to calculate hdry (0-1, unitless)
882	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
883	!_ ================================================================================================================================
884
885	!! 1. Setting ok_hdiff for horizontal diffusion between soil columns
886
887	!Config Key = HYDROL_OK_HDIFF
888	!Config Desc = do horizontal diffusion?
889	!Config If = OK_SECHIBA and .NOT.(HYDROL_CWRR)
890	!Config Def = n
891	!Config Help = If TRUE, then water can diffuse horizontally between
892	!Config the PFTs' water reservoirs.
893	!Config Units = [FLAG]
894	ok_hdiff = .FALSE.
895	CALL getin_p('HYDROL_OK_HDIFF',ok_hdiff)
896
897	!! 2. Make dynamic allocation with the good dimensions
898
899	! one dimension array allocation with possible restart value
900	ALLOCATE (bqsb(kjpindex,nvm),stat=ier)
901	IF (ier.NE.0) THEN
902	WRITE (numout,*) ' error in bqsb allocation. We stop. We need kjpindex words = ',kjpindex
903	STOP 'hydrolc_init'
904	END IF
905	bqsb(:,:) = zero
906
907	ALLOCATE (gqsb(kjpindex,nvm),stat=ier)
908	IF (ier.NE.0) THEN
909	WRITE (numout,*) ' error in gqsb allocation. We stop. We need kjpindex words = ',kjpindex
910	STOP 'hydrolc_init'
911	END IF
912	gqsb(:,:) = zero
913
914	ALLOCATE (dsg(kjpindex,nvm),stat=ier)
915	IF (ier.NE.0) THEN
916	WRITE (numout,*) ' error in dsg allocation. We stop. We need kjpindex words = ',kjpindex
917	STOP 'hydrolc_init'
918	END IF
919	dsg(:,:) = zero
920
921	ALLOCATE (dsp(kjpindex,nvm),stat=ier)
922	IF (ier.NE.0) THEN
923	WRITE (numout,*) ' error in dsp allocation. We stop. We need kjpindex words = ',kjpindex
924	STOP 'hydrolc_init'
925	END IF
926	dsp(:,:) = zero
927
928	! one dimension array allocation
929	ALLOCATE (mean_bqsb(kjpindex),stat=ier)
930	IF (ier.NE.0) THEN
931	WRITE (numout,*) ' error in mean_bqsb allocation. We stop. We need kjpindex words = ',kjpindex
932	STOP 'hydrolc_init'
933	END IF
934	mean_bqsb(:) = zero
935
936	ALLOCATE (mean_gqsb(kjpindex),stat=ier)
937	IF (ier.NE.0) THEN
938	WRITE (numout,*) ' error in mean_gqsb allocation. We stop. We need kjpindex words = ',kjpindex
939	STOP 'hydrolc_init'
940	END IF
941	mean_gqsb(:) = zero
942
943	ALLOCATE (dss(kjpindex,nvm),stat=ier)
944	IF (ier.NE.0) THEN
945	WRITE (numout,*) ' error in dss allocation. We stop. We need kjpindex words = ',kjpindex
946	STOP 'hydrolc_init'
947	END IF
948	dss(:,:) = zero
949
950	ALLOCATE (hdry(kjpindex),stat=ier)
951	IF (ier.NE.0) THEN
952	WRITE (numout,*) ' error in hdry allocation. We stop. We need kjpindex words = ',kjpindex
953	STOP 'hydrolc_init'
954	END IF
955	hdry(:) = zero
956
957	ALLOCATE (precisol(kjpindex,nvm),stat=ier)
958	IF (ier.NE.0) THEN
959	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
960	STOP 'hydrolc_init'
961	END IF
962	precisol(:,:) = zero
963
964	ALLOCATE (gdrainage(kjpindex,nvm),stat=ier)
965	IF (ier.NE.0) THEN
966	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
967	STOP 'hydrolc_init'
968	END IF
969	gdrainage(:,:) = zero
970
971	ALLOCATE (subsnowveg(kjpindex),stat=ier)
972	IF (ier.NE.0) THEN
973	WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex
974	STOP 'hydrolc_init'
975	END IF
976	subsnowveg(:) = zero
977
978	ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier)
979	IF (ier.NE.0) THEN
980	WRITE (numout,) ' error in subsnownobio allocation. We stop. We need kjpindexnnobio words = ', &
981	kjpindex*nnobio
982	STOP 'hydrolc_init'
983	END IF
984	subsnownobio(:,:) = zero
985
986	ALLOCATE (snowmelt(kjpindex),stat=ier)
987	IF (ier.NE.0) THEN
988	WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex
989	STOP 'hydrolc_init'
990	END IF
991	snowmelt(:) = zero
992
993	ALLOCATE (icemelt(kjpindex),stat=ier)
994	IF (ier.NE.0) THEN
995	WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex
996	STOP 'hydrolc_init'
997	END IF
998	icemelt(:) = zero
999
1000	ALLOCATE (subsinksoil(kjpindex),stat=ier)
1001	IF (ier.NE.0) THEN
1002	WRITE (numout,*) ' error in subsinksoil allocation. We stop. We need kjpindex words = ',kjpindex
1003	STOP 'hydrolc_init'
1004	END IF
1005
1006
1007	ALLOCATE (mx_eau_var(kjpindex),stat=ier)
1008	IF (ier.NE.0) THEN
1009	WRITE (numout,*) ' error in mx_eau_var allocation. We stop. We need kjpindex words = ',kjpindex
1010	STOP 'hydrolc_init'
1011	END IF
1012	mx_eau_var(:) = zero
1013
1014	ALLOCATE (ruu_ch(kjpindex),stat=ier)
1015	IF (ier.NE.0) THEN
1016	WRITE (numout,*) ' error in ruu_ch allocation. We stop. We need kjpindex words = ',kjpindex
1017	STOP 'hydrolc_init'
1018	END IF
1019	ruu_ch(:) = zero
1020
1021	ALLOCATE (vegtot(kjpindex),stat=ier)
1022	IF (ier.NE.0) THEN
1023	WRITE (numout,) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindexnvm
1024	STOP 'hydrolc_init'
1025	END IF
1026	vegtot(:) = zero
1027
1028	ALLOCATE (resdist(kjpindex,nvm),stat=ier)
1029	IF (ier.NE.0) THEN
1030	WRITE (numout,) ' error in resdist allocation. We stop. We need kjpindex words = ',kjpindexnvm
1031	STOP 'hydrolc_init'
1032	END IF
1033	resdist(:,:) = zero
1034
1035	! If we check the water balance we need two more variables
1036	IF ( check_waterbal ) THEN
1037
1038	ALLOCATE (tot_water_beg(kjpindex),stat=ier)
1039	IF (ier.NE.0) THEN
1040	WRITE (numout,*) ' error in tot_water_beg allocation. We stop. We need kjpindex words = ',kjpindex
1041	STOP 'hydrolc_init'
1042	END IF
1043
1044	ALLOCATE (tot_water_end(kjpindex),stat=ier)
1045	IF (ier.NE.0) THEN
1046	WRITE (numout,*) ' error in tot_water_end allocation. We stop. We need kjpindex words = ',kjpindex
1047	STOP 'hydrolc_init'
1048	END IF
1049
1050	ALLOCATE (tot_flux(kjpindex),stat=ier)
1051	IF (ier.NE.0) THEN
1052	WRITE (numout,*) ' error in tot_flux allocation. We stop. We need kjpindex words = ',kjpindex
1053	STOP 'hydrol_init'
1054	END IF
1055
1056	ENDIF
1057
1058	! If we use the almaoutputs we need four more variables
1059	IF ( almaoutput ) THEN
1060
1061	ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
1062	IF (ier.NE.0) THEN
1063	WRITE (numout,*) ' error in tot_watveg_beg allocation. We stop. We need kjpindex words = ',kjpindex
1064	STOP 'hydrolc_init'
1065	END IF
1066
1067	ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
1068	IF (ier.NE.0) THEN
1069	WRITE (numout,*) ' error in tot_watveg_end allocation. We stop. We need kjpindex words = ',kjpindex
1070	STOP 'hydrolc_init'
1071	END IF
1072
1073	ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
1074	IF (ier.NE.0) THEN
1075	WRITE (numout,*) ' error in tot_watsoil_beg allocation. We stop. We need kjpindex words = ',kjpindex
1076	STOP 'hydrolc_init'
1077	END IF
1078
1079	ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
1080	IF (ier.NE.0) THEN
1081	WRITE (numout,*) ' error in tot_watsoil_end allocation. We stop. We need kjpindex words = ',kjpindex
1082	STOP 'hydrolc_init'
1083	END IF
1084
1085	ALLOCATE (delsoilmoist(kjpindex),stat=ier)
1086	IF (ier.NE.0) THEN
1087	WRITE (numout,*) ' error in delsoilmoist allocation. We stop. We need kjpindex words = ',kjpindex
1088	STOP 'hydrolc_init'
1089	END IF
1090
1091	ALLOCATE (delintercept(kjpindex),stat=ier)
1092	IF (ier.NE.0) THEN
1093	WRITE (numout,*) ' error in delintercept. We stop. We need kjpindex words = ',kjpindex
1094	STOP 'hydrolc_init'
1095	END IF
1096
1097	ALLOCATE (delswe(kjpindex),stat=ier)
1098	IF (ier.NE.0) THEN
1099	WRITE (numout,*) ' error in delswe. We stop. We need kjpindex words = ',kjpindex
1100	STOP 'hydrolc_init'
1101	ENDIF
1102
1103	ALLOCATE (snow_beg(kjpindex),stat=ier)
1104	IF (ier.NE.0) THEN
1105	WRITE (numout,*) ' error in snow_beg allocation. We stop. We need kjpindex words =',kjpindex
1106	STOP 'hydrolc_init'
1107	END IF
1108
1109	ALLOCATE (snow_end(kjpindex),stat=ier)
1110	IF (ier.NE.0) THEN
1111	WRITE (numout,*) ' error in snow_end allocation. We stop. We need kjpindex words =',kjpindex
1112	STOP 'hydrolc_init'
1113	END IF
1114
1115	ENDIF
1116
1117	!! 3. Initialization of hydrologic variables:
1118
1119	!! 3.a Read data from restart input file (opened by sechiba_init)
1120	!! for HYDROLOGIC processes
1121	IF (printlev>=3) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
1122
1123	var_name= 'snow'
1124	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1125	CALL ioconf_setatt_p('LONG_NAME','Snow mass')
1126	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
1127
1128	var_name= 'snow_age'
1129	CALL ioconf_setatt_p('UNITS', 'd')
1130	CALL ioconf_setatt_p('LONG_NAME','Snow age')
1131	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
1132
1133	var_name= 'snow_nobio'
1134	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1135	CALL ioconf_setatt_p('LONG_NAME','Snow on other surface types')
1136	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
1137
1138	var_name= 'snow_nobio_age'
1139	CALL ioconf_setatt_p('UNITS', 'd')
1140	CALL ioconf_setatt_p('LONG_NAME','Snow age on other surface types')
1141	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
1142
1143	var_name= 'humrel'
1144	CALL ioconf_setatt_p('UNITS', '-')
1145	CALL ioconf_setatt_p('LONG_NAME','Soil moisture stress')
1146	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g)
1147
1148	var_name= 'vegstress'
1149	CALL ioconf_setatt_p('UNITS', '-')
1150	CALL ioconf_setatt_p('LONG_NAME','Vegetation growth moisture stress')
1151	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
1152
1153	var_name= 'bqsb'
1154	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1155	CALL ioconf_setatt_p('LONG_NAME','Deep soil moisture')
1156	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., bqsb, "gather", nbp_glo, index_g)
1157
1158	var_name= 'gqsb'
1159	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1160	CALL ioconf_setatt_p('LONG_NAME','Surface soil moisture')
1161	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., gqsb, "gather", nbp_glo, index_g)
1162
1163	var_name= 'dsg'
1164	CALL ioconf_setatt_p('UNITS', 'm')
1165	CALL ioconf_setatt_p('LONG_NAME','Depth of upper reservoir')
1166	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsg, "gather", nbp_glo, index_g)
1167
1168	var_name= 'dsp'
1169	CALL ioconf_setatt_p('UNITS', 'm')
1170	CALL ioconf_setatt_p('LONG_NAME','Depth to lower reservoir')
1171	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsp, "gather", nbp_glo, index_g)
1172
1173	var_name= 'dss'
1174	CALL ioconf_setatt_p('UNITS', 'm')
1175	CALL ioconf_setatt_p('LONG_NAME','Hauteur au dessus du reservoir de surface')
1176	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dss, "gather", nbp_glo, index_g)
1177
1178	var_name= 'hdry'
1179	CALL ioconf_setatt_p('UNITS', 'm')
1180	CALL ioconf_setatt_p('LONG_NAME','Dry soil heigth in meters')
1181	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., hdry, "gather", nbp_glo, index_g)
1182
1183	var_name= 'qsintveg'
1184	CALL ioconf_setatt_p('UNITS', 'kg/m^2')
1185	CALL ioconf_setatt_p('LONG_NAME','Intercepted moisture')
1186	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
1187
1188	var_name= 'resdist'
1189	CALL ioconf_setatt_p('UNITS', '-')
1190	CALL ioconf_setatt_p('LONG_NAME','Distribution of reservoirs')
1191	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
1192
1193	! Read drysoil_frac. It will be initalized later in hydrolc_var_init if the varaible is not find in restart file.
1194	CALL restget_p (rest_id, 'drysoil_frac', nbp_glo, 1 , 1, kjit, .TRUE., drysoil_frac, "gather", nbp_glo, index_g)
1195
1196	! Read rsol. It will be initalized later in hydrolc_var_init if the varaible is not find in restart file.
1197	CALL restget_p (rest_id, 'rsol', nbp_glo, 1 , 1, kjit, .TRUE., rsol, "gather", nbp_glo, index_g)
1198
1199	! shumdiag : initialization if not in restart file will be done in hydrolc_var_init
1200	CALL restget_p (rest_id, 'shumdiag', nbp_glo, nbdl , 1, kjit, .TRUE., shumdiag, "gather", nbp_glo, index_g)
1201
1202	CALL restget_p (rest_id, 'mean_bqsb', nbp_glo, 1 , 1, kjit, .TRUE., mean_bqsb, "gather", nbp_glo, index_g)
1203	CALL restget_p (rest_id, 'mean_gqsb', nbp_glo, 1 , 1, kjit, .TRUE., mean_gqsb, "gather", nbp_glo, index_g)
1204
1205	var_name= 'mx_eau_var'
1206	CALL ioconf_setatt_p('UNITS', '')
1207	CALL ioconf_setatt_p('LONG_NAME','')
1208	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., mx_eau_var, "gather", nbp_glo, index_g)
1209
1210	var_name= 'vegtot'
1211	CALL ioconf_setatt_p('UNITS', '')
1212	CALL ioconf_setatt_p('LONG_NAME','')
1213	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., vegtot, "gather", nbp_glo, index_g)
1214
1215
1216	!! 3.b Assign default values if non were found in the restart file
1217	!Config Key = HYDROL_SNOW
1218	!Config Desc = Initial snow mass if not found in restart
1219	!Config If = OK_SECHIBA
1220	!Config Def = 0.0
1221	!Config Help = The initial value of snow mass if its value is not found
1222	!Config in the restart file. This should only be used if the model is
1223	!Config started without a restart file.
1224	!Config Units = [kg/m^2]
1225	CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', zero)
1226
1227	!Config Key = HYDROL_SNOWAGE
1228	!Config Desc = Initial snow age if not found in restart
1229	!Config If = OK_SECHIBA
1230	!Config Def = 0.0
1231	!Config Help = The initial value of snow age if its value is not found
1232	!Config in the restart file. This should only be used if the model is
1233	!Config started without a restart file.
1234	!Config Units = [days]
1235	CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', zero)
1236
1237	!Config Key = HYDROL_SNOW_NOBIO
1238	!Config Desc = Initial snow amount on ice, lakes, etc. if not found in restart
1239	!Config If = OK_SECHIBA
1240	!Config Def = 0.0
1241	!Config Help = The initial value of snow if its value is not found
1242	!Config in the restart file. This should only be used if the model is
1243	!Config started without a restart file.
1244	!Config Units = [m]
1245	CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero)
1246
1247	!Config Key = HYDROL_SNOW_NOBIO_AGE
1248	!Config Desc = Initial snow age on ice, lakes, etc. if not found in restart
1249	!Config If = OK_SECHIBA
1250	!Config Def = 0.0
1251	!Config Help = The initial value of snow age if its value is not found
1252	!Config in the restart file. This should only be used if the model is
1253	!Config started without a restart file.
1254	!Config Units = [days]
1255	CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero)
1256
1257	!Config Key = HYDROL_HUMR
1258	!Config Desc = Initial soil moisture stress if not found in restart
1259	!Config If = OK_SECHIBA
1260	!Config Def = 1.0
1261	!Config Help = The initial value of soil moisture stress if its value is not found
1262	!Config in the restart file. This should only be used if the model is
1263	!Config started without a restart file.
1264	!Config Units = [-]
1265	CALL setvar_p (humrel, val_exp,'HYDROL_HUMR', un)
1266	CALL setvar_p (vegstress, val_exp,'HYDROL_HUMR', un)
1267
1268	!Config Key = HYDROL_BQSB
1269	!Config Desc = Initial restart deep soil moisture if not found in restart
1270	!Config If = OK_SECHIBA
1271	!Config Def = 999999.
1272	!Config Help = The initial value of deep soil moisture if its value is not found
1273	!Config in the restart file. This should only be used if the model is
1274	!Config started without a restart file. Default behaviour is a saturated soil.
1275	!Config Units = [kg/m^2]
1276	CALL setvar_p (bqsb, val_exp, 'HYDROL_BQSB', wmax_veg*zmaxh)
1277
1278	!Config Key = HYDROL_GQSB
1279	!Config Desc = Initial upper soil moisture if not found in restart
1280	!Config If = OK_SECHIBA
1281	!Config Def = 0.0
1282	!Config Help = The initial value of upper soil moisture if its value is not found
1283	!Config in the restart file. This should only be used if the model is
1284	!Config started without a restart file.
1285	!Config Units = [kg/m^2]
1286	CALL setvar_p (gqsb, val_exp, 'HYDROL_GQSB', zero)
1287
1288	!Config Key = HYDROL_DSG
1289	!Config Desc = Initial upper reservoir depth if not found in restart
1290	!Config If = OK_SECHIBA
1291	!Config Def = 0.0
1292	!Config Help = The initial value of upper reservoir depth if its value is not found
1293	!Config in the restart file. This should only be used if the model is
1294	!Config started without a restart file.
1295	!Config Units = [m]
1296	CALL setvar_p (dsg, val_exp, 'HYDROL_DSG', zero)
1297
1298	!Config Key = HYDROL_QSV
1299	!Config Desc = Initial water on canopy if not found in restart
1300	!Config If = OK_SECHIBA
1301	!Config Def = 0.0
1302	!Config Help = The initial value of moisture on canopy if its value
1303	!Config is not found in the restart file. This should only be used if
1304	!Config the model is started without a restart file.
1305	!Config Units = [kg/m^2]
1306	CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', zero)
1307
1308	!! 3.c Specific calculations to define default values for dry soil depths : dsp, dss, and hdry
1309	! set inital value for dsp if needed
1310	!Config Key = HYDROL_DSP
1311	!Config Desc = Initial dry soil above upper reservoir if not found in restart
1312	!Config If = OK_SECHIBA
1313	!Config Def = 999999.
1314	!Config Help = The initial value of dry soil above upper reservoir if its value
1315	!Config is not found in the restart file. This should only be used if
1316	!Config the model is started without a restart file. The default behaviour
1317	!Config is to compute it from the variables above. Should be OK most of
1318	!Config the time.
1319	!Config Units = [m]
1320	!
1321	! Calculate temporary variable to use for initialiaion of dsp.
1322	! This variable will only be used if no value is set in run.def
1323	DO jv = 1,nvm
1324	zdsp(:,jv) = zmaxh - bqsb(:,jv) / wmax_veg(jv)
1325	END DO
1326
1327	CALL setvar_p(dsp, val_exp, 'HYDROL_DSP', zdsp)
1328
1329	! set inital value for dss
1330	DO jv = 1,nvm
1331	tmpdss(:,jv) = dsg(:,jv) - gqsb(:,jv) / wmax_veg(jv)
1332	END DO
1333
1334	! Initialize dss if it is not found in restart file
1335	IF ( ALL( dss(:,:) == val_exp ) ) THEN
1336	dss(:,:)=tmpdss(:,:)
1337	END IF
1338
1339	! set inital value for hdry if not found in restart file
1340	! see hydrolc_soil, step 8.4
1341	IF ( ALL( hdry(:) == val_exp ) ) THEN
1342	a_subgrd(:) = zero
1343	DO ji=1,kjpindex
1344	IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN
1345
1346	IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN
1347
1348	! Ajouts Nathalie - Fred - le 28 Mars 2006
1349	a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),zero)/dsg_min,un)
1350
1351	ENDIF
1352	ENDIF
1353	ENDDO
1354
1355	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
1356	! revu 22 novembre 2007
1357	hdry(:) = a_subgrd(:)dss(:,1) + (un-a_subgrd(:))dsp(:,1)
1358	ENDIF
1359
1360	! There is no need to configure the initialisation of resdist. If not available it is the vegetation map
1361	IF ( MINVAL(resdist) .EQ. MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
1362	resdist(:,1) = tot_bare_soil(:)
1363	resdist(:,2:nvm) = veget(:,2:nvm)
1364	ENDIF
1365
1366	!! 3.d Compute vegtot (remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1)
1367	IF (ALL(vegtot(:)==val_exp)) THEN
1368	DO ji = 1, kjpindex
1369	vegtot(ji) = SUM(veget(ji,2:nvm)) + tot_bare_soil(ji)
1370	ENDDO
1371	END IF
1372
1373	! Initialize variables for explictsnow module by reading restart file
1374	IF (ok_explicitsnow) THEN
1375	CALL explicitsnow_initialize( kjit, kjpindex, rest_id, snowrho, &
1376	snowtemp, snowdz, snowheat, snowgrain )
1377	END IF
1378
1379	IF (printlev>=3) WRITE (numout,*) ' hydrolc_init done '
1380
1381	END SUBROUTINE hydrolc_init
1382
1383
1384	!! ================================================================================================================================
1385	!! SUBROUTINE : hydrolc_clear
1386	!!
1387	!>\BRIEF Deallocates arrays that were allocated in hydrolc_init and hydrolc_var_init.
1388	!!
1389	!! DESCRIPTION : None
1390	!!
1391	!! RECENT CHANGES : None
1392	!!
1393	!! MAIN OUTPUT VARIABLE(S) : None
1394	!!
1395	!! REFERENCE(S) : None
1396	!!
1397	!! FLOWCHART : None
1398	!!\n
1399	!_ ================================================================================================================================
1400
1401	SUBROUTINE hydrolc_clear()
1402
1403	IF (ALLOCATED (bqsb)) DEALLOCATE (bqsb)
1404	IF (ALLOCATED (gqsb)) DEALLOCATE (gqsb)
1405	IF (ALLOCATED (dsg)) DEALLOCATE (dsg)
1406	IF (ALLOCATED (dsp)) DEALLOCATE (dsp)
1407	IF (ALLOCATED (mean_bqsb)) DEALLOCATE (mean_bqsb)
1408	IF (ALLOCATED (mean_gqsb)) DEALLOCATE (mean_gqsb)
1409	IF (ALLOCATED (dss)) DEALLOCATE (dss)
1410	IF (ALLOCATED (hdry)) DEALLOCATE (hdry)
1411	IF (ALLOCATED (precisol)) DEALLOCATE (precisol)
1412	IF (ALLOCATED (gdrainage)) DEALLOCATE (gdrainage)
1413	IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg)
1414	IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio)
1415	IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt)
1416	IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt)
1417	IF (ALLOCATED (subsinksoil)) DEALLOCATE (subsinksoil)
1418	IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var)
1419	IF (ALLOCATED (ruu_ch)) DEALLOCATE (ruu_ch)
1420	IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot)
1421	IF (ALLOCATED (resdist)) DEALLOCATE (resdist)
1422	IF (ALLOCATED (tot_water_beg)) DEALLOCATE (tot_water_beg)
1423	IF (ALLOCATED (tot_water_end)) DEALLOCATE (tot_water_end)
1424	IF (ALLOCATED (tot_flux)) DEALLOCATE (tot_flux)
1425	IF (ALLOCATED (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
1426	IF (ALLOCATED (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
1427	IF (ALLOCATED (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
1428	IF (ALLOCATED (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
1429	IF (ALLOCATED (delsoilmoist)) DEALLOCATE (delsoilmoist)
1430	IF (ALLOCATED (delintercept)) DEALLOCATE (delintercept)
1431	IF (ALLOCATED (snow_beg)) DEALLOCATE (snow_beg)
1432	IF (ALLOCATED (snow_end)) DEALLOCATE (snow_end)
1433	IF (ALLOCATED (delswe)) DEALLOCATE (delswe)
1434
1435	END SUBROUTINE hydrolc_clear
1436
1437
1438	!! ================================================================================================================================
1439	!! SUBROUTINE : hydrolc_var_init
1440	!!
1441	!>\BRIEF This routine initializes diagnostic hydrologic variables.
1442	!!
1443	!! DESCRIPTION : The following variables are assigned :
1444	!! - Soil characteristics : soil water capacities mx_eau_var and ruu_ch \f$(kg m^{-2})\f$ and \f$(kg m^{-3})\f$
1445	!! - Initial values for diagnostic variables : mean_bqsb, mean_gqsb, mean_dsg, shumdiag, drysoil_frac, rsol, litterhumdiag
1446	!!
1447	!! RECENT CHANGE(S) : None
1448	!!
1449	!! MAIN OUTPUT VARIABLE(S) : rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag
1450	!!
1451	!! REFERENCE(S) : None
1452	!!
1453	!! FLOWCHART : None
1454	!! \n
1455	!_ ================================================================================================================================
1456
1457	SUBROUTINE hydrolc_var_init (kjpindex, veget, veget_max, tot_bare_soil, &
1458	rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag)
1459
1460	!! 0. Variable and parameter declaration
1461
1462	!! 0.1 Input variables
1463
1464	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
1465	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by vegetation for
1466	!! each PFT, except for soil (0-1, unitless)
1467	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! PFT fractions within grid-cells: max value of veget for
1468	!! vegetation PFTs and min value for bare soil (0-1,
1469	!! unitless)
1470	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
1471
1472	!! 0.2 Output variables
1473
1474	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: rsol !! Resistance to bare soil evaporation
1475	!! @tex ($s m^{-1}$) @endtex
1476	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: drysoil_frac !! Fraction of visible dry soil (0-1, unitless)
1477	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: mx_eau_var !! Maximum water content of the soil
1478	!! @tex ($kg m^{-2}$) @endtex
1479	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ruu_ch !! Volumetric soil water capacity
1480	!! @tex ($kg m^{-3}$) @endtex
1481	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: shumdiag !! Mean relative soil moisture, diagnosed for
1482	!! thermosoil.f90 (0-1, unitless)
1483
1484	!! 0.3 Modified variables
1485
1486	!! 0.4 Local variables
1487
1488	INTEGER(i_std) :: ji,jv, jd !! Indices for grid-cells, PFTs and diagnostic levels in
1489	!! the soil (unitless)
1490	REAL(r_std), DIMENSION(kjpindex) :: mean_dsg !! Mean Depth of the top layer, averaged over the soil
1491	!! columns (m)
1492	REAL(r_std) :: gtr, btr !! Ancillary variables to compute shumdiag (0-1, unitless)
1493	REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl !! Temporary diagnostic levels in the soil (m)
1494	!_ ================================================================================================================================
1495
1496	!! 1. Vertical diagnostic levels to compute shumdiag (relative moisture for thermosoil)
1497
1498	tmp_dl(1) = 0
1499	tmp_dl(2:nbdl+1) = diaglev(1:nbdl)
1500
1501	!! 2. Calculation of mx_eau_var and ruu_ch (soil water capacities)
1502
1503	IF (ALL(mx_eau_var(:)==val_exp)) THEN
1504	mx_eau_var(:) = zero
1505
1506	DO ji = 1,kjpindex
1507	DO jv = 1,nvm
1508	!MM veget(:,1) BUG ??!!!!!!!!!!!
1509	IF (jv .EQ. 1) THEN
1510	mx_eau_var(ji) = mx_eau_var(ji) + tot_bare_soil(ji)wmax_veg(jv)zmaxh
1511	ELSE
1512	mx_eau_var(ji) = mx_eau_var(ji) + veget(ji,jv)wmax_veg(jv)zmaxh
1513	ENDIF
1514	END DO
1515	IF (vegtot(ji) .GT. zero) THEN
1516	mx_eau_var(ji) = mx_eau_var(ji)/vegtot(ji)
1517	ELSE
1518	! lakes, ice, cities...
1519	mx_eau_var(ji) = mx_eau_nobio*zmaxh
1520	ENDIF
1521	END DO
1522	END IF
1523	!! Initialize ruu_ch
1524	ruu_ch(:) = mx_eau_var(:) / zmaxh
1525
1526	!! 3.-4. Initial values of the mean soil layer depths and water contents and shumdiag
1527
1528	IF (ALL(mean_bqsb(:)==val_exp) .OR. ALL(mean_gqsb(:)==val_exp) .OR. ALL(shumdiag(:,:)==val_exp)) THEN
1529	!! 3. Initial values of the mean soil layer depths and water contents
1530	! could be done with SUM instruction but this kills vectorization
1531	mean_bqsb(:) = zero
1532	mean_gqsb(:) = zero
1533	mean_dsg(:) = zero
1534	DO jv = 1, nvm
1535	DO ji = 1, kjpindex
1536	mean_bqsb(ji) = mean_bqsb(ji) + resdist(ji,jv)*bqsb(ji,jv)
1537	mean_gqsb(ji) = mean_gqsb(ji) + resdist(ji,jv)*gqsb(ji,jv)
1538	mean_dsg(ji) = mean_dsg(ji) + resdist(ji,jv)*dsg(ji,jv)
1539	ENDDO
1540	ENDDO
1541	mean_dsg(:) = MAX( mean_dsg(:), mean_gqsb(:)/ruu_ch(:) )
1542
1543	DO ji = 1, kjpindex
1544	IF (vegtot(ji) .GT. zero) THEN
1545	mean_bqsb(ji) = mean_bqsb(ji)/vegtot(ji)
1546	mean_gqsb(ji) = mean_gqsb(ji)/vegtot(ji)
1547	mean_dsg(ji) = mean_dsg(ji)/vegtot(ji)
1548	ENDIF
1549	ENDDO
1550
1551
1552	!! 4. Initial value of shumdiag (see hydrolc_soil, 8.2 for details)
1553	DO jd = 1,nbdl
1554	DO ji = 1,kjpindex
1555	IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN
1556	shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji)
1557	ELSE
1558	IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN
1559	gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd))
1560	btr = 1 - gtr
1561	shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + &
1562	& btr*mean_bqsb(ji)/mx_eau_var(ji)
1563	ELSE
1564	shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji)
1565	ENDIF
1566	ENDIF
1567	shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero)
1568	ENDDO
1569	ENDDO
1570	END IF
1571
1572	!! 5. Initialize drysoil_frac if it was not found in the restart file
1573	IF (ALL(drysoil_frac(:) == val_exp)) THEN
1574	drysoil_frac(:) = MIN(MAX(dss(:,1),zero)*10._r_std, un)
1575	END IF
1576
1577	!! 6. Initial value of the resistance to bare soil evaporation (see hydrolc_soil, 8.4 for details)
1578
1579	IF (ALL(rsol(:)==val_exp)) THEN
1580	rsol(:) = -un
1581	DO ji = 1, kjpindex
1582	IF (tot_bare_soil(ji) .GE. min_sechiba) THEN
1583
1584	! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
1585	! on modifie le rsol pour que la resistance croisse subitement si on s'approche
1586	! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70
1587	!rsol(ji) = dss(ji,1) * rsol_cste
1588	!rsol(ji) = ( drysoil_frac(ji) + un/(10.(zmaxh - drysoil_frac(ji))+1.e-10)2 ) rsol_cste
1589	rsol(ji) = ( hdry(ji) + un/(10.(zmaxh - hdry(ji))+1.e-10)2 ) rsol_cste
1590	ENDIF
1591	ENDDO
1592	END IF
1593
1594	IF (printlev>=3) WRITE (numout,*) ' hydrolc_var_init done '
1595
1596	END SUBROUTINE hydrolc_var_init
1597
1598
1599	!! ================================================================================================================================
1600	!! SUBROUTINE : hydrolc_snow
1601	!!
1602	!>\BRIEF This routine computes accumulation, sublimation, snowmelt and snow ageing
1603	!! over vegetated and nobio (eg land-ice) areas. Snowdepth is then updated.
1604	!!
1605	!! DESCRIPTION : In this routine the snow-related processes accumulation, sublimation,
1606	!! melting and ageing are computed separately on vegetated and nobio areas. This subroutine
1607	!! is the first one among the ones that compute the physical processes. The water budgets
1608	!! of the interception reservoir and the soil are performed afterwards.\n
1609	!!
1610	!! - Values of the main parameters used in this routine are :\n
1611	!! tp_00=273.15K : freezing point\n
1612	!! snowcri=1.5\f$kg m^{-2}\f$ : critical snowmass for sublimation\n
1613	!! sneige=snowcri/1000._r_std : critical snowmass for snow melt\n
1614	!! maxmass_snow=3000 \f$kg m^{-2}\f$ (~ 10m snow) : critical snowmass for snow stock\n
1615	!! snow_trans=0.3 (\f$kg m^{-2}\f$) : critical constant for snow ageing\n
1616	!! max_snow_age=50 (days) : maximum snow age\n
1617	!! sn_dens=330 (\f$kg m^{-3}\f$) : snow density\n
1618	!! one_day=86400 (s) : one day, expressed in seconds...
1619	!!
1620	!! - Accumulation\n
1621	!! On the vegetated areas, the snow mass increases due to area-weighted snow
1622	!! precipitations precip_snow; on nobio areas, rain additionnaly converts into snow.\n
1623	!!
1624	!! - Sublimation\n
1625	!! Sublimation on vegetated and nobio areas is calculated as the area-weighed fraction
1626	!! of vevapsno (from enerbil). In case snow on vegetated areas is limited (less than snowcri)
1627	!! and a significant nobio area exists, the nobio area accomodates the whole sublimation
1628	!! vevapsno. The nobio area also accomodates the possible excess sublimation from vegetated
1629	!! areas, which otherwise goes into bare soil evaporation. In this case vevapsno is updated.\n
1630	!!
1631	!! - Melting\n
1632	!! Over vegetated and nobio areas, snowmelt occurs if the calculated soil temperature
1633	!! temp_soil_new(ji) exceeds the freezing point tp_00. The energy required to warm up the soil
1634	!! surface above tp_00 is converted into melting, according to the following equation :
1635	!! \latexonly
1636	!! \input{hydrolcsnow1.tex}
1637	!! \endlatexonly
1638	!! \n
1639	!! with soilcap the soil heat capacity @tex ($J K^{-1}$) @endtex and chalfu0 the latent heat of
1640	!! fusion.\n
1641	!! A special case occurs in areas with snowmass exceeding maxmass_snow, where
1642	!! all the snow in excess of maxmass_snow adds to the total melting tot_melt.
1643	!! This excess melting is considered as additionnal snowmelt for vegetated areas
1644	!! and as ice_melt for nobio areas.\n
1645	!!
1646	!! - Ageing\n
1647	!! Over vegetated areas, during the time step dt_radia (usually 1800 s), the snow age is
1648	!! incremented by a d_age calculated as follow:
1649	!! \latexonly
1650	!! \input{hydrolcsnow2.tex}
1651	!! \endlatexonly
1652	!! \n
1653	!! Over nobio areas and at subfreezing temperatures, snow ageing is limited by very
1654	!! slow snow metamorphism, hence d_age is reduced as follow:
1655	!! \latexonly
1656	!! \input{hydrolcsnow3.tex}
1657	!! \endlatexonly
1658	!! \n
1659	!! Scientific reference for this parametrization choice is obviously missing !
1660	!! At the end of the routine the snowheight is diagnosed using the fixed
1661	!! snow density sn_dens.
1662	!!
1663	!! RECENT CHANGE(S) : None
1664	!!
1665	!! MAIN OUTPUT VARIABLES:
1666	!! for vegetated and nobio areas respectively:\n
1667	!! snow(kjpindex) and snow_nobio(kjpindex)\n
1668	!! snow_age(kjpindex) and snow_nobio_age(kjpindex)\n
1669	!! for the whole grid-cell:\n
1670	!! vevapnu(kjpindex)
1671	!! vevapsno(kjpindex)
1672	!! tot_melt(kjpindex)
1673	!! snowdepth(kjpindex)
1674	!!
1675	!! REFERENCES : None
1676	!!
1677	!! FLOWCHART : None
1678	!! \n
1679	!_ ================================================================================================================================
1680
1681	SUBROUTINE hydrolc_snow (kjpindex, precip_rain, precip_snow , temp_sol_new, soilcap,&
1682	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
1683	& tot_melt, snowdepth)
1684
1685	!! 0. Variable and parameter declaration
1686
1687	!! 0.1 Input variabless
1688
1689	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (unitless)
1690	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall @tex ($kg m^{-2}$) @endtex
1691	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation @tex ($kg m^{-2}$) @endtex
1692	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! New soil temperature (K)
1693	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil heat capacity @tex ($J K^{-1}$) @endtex
1694	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
1695	!! (unitless)
1696	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
1697	!! (unitless)
1698
1699	!! 0.2 Output variables
1700
1701	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt @tex ($kg m^{-2}$) @endtex
1702	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowdepth !! Snow depth (m)
1703
1704	!! 0.3 Modified variables
1705
1706	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation @tex ($kg m^{-2}$) @endtex
1707	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
1708	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass @tex ($kg m^{-2}$) @endtex
1709	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age (days)
1710	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Snomass on nobio areas
1711	!! @tex ($kg m^{-2}$) @endtex
1712	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age !! Snow age on ice, lakes, ...(days)
1713
1714	!! 0.4 Local variables
1715
1716	INTEGER(i_std) :: ji, jv !! indices (unitless)
1717	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change (days)
1718	REAL(r_std), DIMENSION (kjpindex) :: xx !! temporary
1719	REAL(r_std) :: snowmelt_tmp !! temporary @tex ($kg m^{-2}$) @endtex
1720	REAL(r_std) :: snow_d1k !! The amount of snow that corresponds to a 1K cooling
1721	LOGICAL, DIMENSION (kjpindex) :: warnings
1722	LOGICAL :: any_warning
1723	!_ ================================================================================================================================
1724
1725	!! 1. initialisation
1726
1727	DO jv = 1, nnobio
1728	DO ji=1,kjpindex
1729	subsnownobio(ji,jv) = zero
1730	ENDDO
1731	ENDDO
1732	DO ji=1,kjpindex
1733	subsnowveg(ji) = zero
1734	snowmelt(ji) = zero
1735	icemelt(ji) = zero
1736	tot_melt(ji) = zero
1737	ENDDO
1738
1739	!! 2. On vegetation
1740
1741	! 2.1. It is snowing
1742	DO ji=1,kjpindex
1743	snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji)
1744	ENDDO
1745
1746	DO ji=1,kjpindex
1747
1748	! 2.2. Sublimation - separate between vegetated and no-veget fractions
1749	! Care has to be taken as we might have sublimation from the
1750	! the frac_nobio while there is no snow on the rest of the grid.
1751	IF ( snow(ji) > snowcri ) THEN
1752	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
1753	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
1754	ELSE
1755
1756	! Correction Nathalie - Juillet 2006.
1757	! On doit d'abord tester s'il existe un frac_nobio!
1758	! Pour le moment je ne regarde que le iice
1759	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1760	subsnownobio(ji,iice) = vevapsno(ji)
1761	subsnowveg(ji) = zero
1762	ELSE
1763	subsnownobio(ji,iice) = zero
1764	subsnowveg(ji) = vevapsno(ji)
1765	ENDIF
1766	ENDIF
1767	ENDDO
1768
1769	warnings(:) = .FALSE.
1770	any_warning = .FALSE.
1771	DO ji=1,kjpindex
1772
1773	! 2.2.1 Check that sublimation on the vegetated fraction is possible.
1774	IF (subsnowveg(ji) .GT. snow(ji)) THEN
1775
1776	! What could not be sublimated goes into bare soil evaporation
1777	! Nathalie - Juillet 2006 - il faut avant tout tester s'il existe du
1778	! frac_nobio sur ce pixel pour eviter de puiser dans le sol!
1779	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1780	subsnownobio(ji,iice) = subsnownobio(ji,iice) + (subsnowveg(ji) - snow(ji))
1781	ELSE
1782	vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji))
1783	warnings(ji) = .TRUE.
1784	any_warning = .TRUE.
1785	ENDIF
1786
1787	! Sublimation is thus limited to what is available
1788	subsnowveg(ji) = snow(ji)
1789	snow(ji) = zero
1790	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
1791	ELSE
1792	snow(ji) = snow(ji) - subsnowveg(ji)
1793	ENDIF
1794	ENDDO
1795	IF ( any_warning ) THEN
1796	WRITE(numout,*)' ATTENTION on prend de l eau au sol nu sur au moins un point car evapsno est trop fort!'
1797	!!$ DO ji=1,kjpindex
1798	!!$ IF ( warnings(ji) ) THEN
1799	!!$ WRITE(numout,*)' ATTENTION on prend de l eau au sol nu car evapsno est trop fort!'
1800	!!$ WRITE(numout,)' ',ji,' vevapnu (en mm/jour) = ',vevapnu(ji)one_day/dt_sechiba
1801	!!$ ENDIF
1802	!!$ ENDDO
1803	ENDIF
1804
1805	warnings(:) = .FALSE.
1806	any_warning = .FALSE.
1807	DO ji=1,kjpindex
1808
1809	! 2.3. snow melt only if temperature positive
1810	IF (temp_sol_new(ji).GT.tp_00) THEN
1811
1812	IF (snow(ji).GT.sneige) THEN
1813
1814	snowmelt(ji) = (un - frac_nobio(ji,iice))(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1815
1816	! 2.3.1 enough snow for melting or not
1817	IF (snowmelt(ji).LT.snow(ji)) THEN
1818	snow(ji) = snow(ji) - snowmelt(ji)
1819	ELSE
1820	snowmelt(ji) = snow(ji)
1821	snow(ji) = zero
1822	END IF
1823
1824	ELSEIF (snow(ji).GE.zero) THEN
1825
1826	! 2.3.2 not enough snow
1827	snowmelt(ji) = snow(ji)
1828	snow(ji) = zero
1829	ELSE
1830
1831	! 2.3.3 negative snow - now snow melt
1832	snow(ji) = zero
1833	snowmelt(ji) = zero
1834	warnings(ji) = .TRUE.
1835	any_warning = .TRUE.
1836
1837	END IF
1838
1839	ENDIF
1840	ENDDO
1841	IF ( any_warning ) THEN
1842	DO ji=1,kjpindex
1843	IF ( warnings(ji) ) THEN
1844	WRITE(numout,*) 'hydrolc_snow: WARNING! snow was negative and was reset to zero for point ',ji,'. '
1845	ENDIF
1846	ENDDO
1847	ENDIF
1848
1849
1850
1851	!! 2.4 Snow melts above a threshold
1852	! Ice melt only if there is more than a given mass : maxmass_snow,
1853	! But the snow cannot melt more in one time step to what corresponds to
1854	! a 1K cooling. This will lead to a progressive melting of snow above
1855	! maxmass_snow but it is needed as a too strong cooling can destabilise the model.
1856	DO ji=1,kjpindex
1857	IF ( snow(ji) .GT. maxmass_snow ) THEN
1858	snow_d1k = un * soilcap(ji) / chalfu0
1859	snowmelt(ji) = snowmelt(ji) + MIN((snow(ji) - maxmass_snow),snow_d1k)
1860	snow(ji) = snow(ji) - snowmelt(ji)
1861	IF ( printlev >= 3 ) WRITE (numout,*) "Snow was above maxmass_snow (", maxmass_snow,") and we melted ", snowmelt(ji)
1862	ENDIF
1863	END DO
1864
1865	!! 3. On Land ice
1866
1867	DO ji=1,kjpindex
1868
1869	!! 3.1. It is snowing
1870	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
1871	& frac_nobio(ji,iice)*precip_rain(ji)
1872
1873	!! 3.2. Sublimation - was calculated before it can give us negative snow_nobio but that is OK
1874	! Once it goes below a certain values (-maxmass_snow for instance) we should kill
1875	! the frac_nobio(ji,iice) !
1876	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
1877
1878	!! 3.3. snow melt only for continental ice fraction
1879	snowmelt_tmp = zero
1880	IF (temp_sol_new(ji) .GT. tp_00) THEN
1881
1882	! 3.3.1 If there is snow on the ice-fraction it can melt
1883	snowmelt_tmp = frac_nobio(ji,iice)(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1884
1885	IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN
1886	snowmelt_tmp = MAX( zero, snow_nobio(ji,iice))
1887	ENDIF
1888	snowmelt(ji) = snowmelt(ji) + snowmelt_tmp
1889	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp
1890
1891	ENDIF
1892
1893	!! 3.4 Snow melts over a threshold
1894	! Ice melt only if there is more than a given mass : maxmass_snow,
1895	! But the snow cannot melt more in one time step to what corresponds to
1896	! a 1K cooling. This will lead to a progressive melting of snow above
1897	! maxmass_snow but it is needed as a too strong cooling can destabilise the model.
1898	!
1899	IF ( snow_nobio(ji,iice) .GT. maxmass_snow ) THEN
1900	snow_d1k = un * soilcap(ji) / chalfu0
1901	icemelt(ji) = MIN((snow_nobio(ji,iice) - maxmass_snow),snow_d1k)
1902	snow_nobio(ji,iice) = snow_nobio(ji,iice) - icemelt(ji)
1903
1904	IF ( printlev >= 3 ) WRITE (numout,*) "Snow was above maxmass_snow ON ICE (", maxmass_snow,") and we melted ", icemelt(ji)
1905	ENDIF
1906
1907	END DO
1908
1909
1910	!! 4. On other surface types - not done yet
1911
1912	IF ( nnobio .GT. 1 ) THEN
1913	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
1914	CALL ipslerr_p (3,'hydrolc_snow', '', &
1915	& 'CANNOT TREAT SNOW ON THESE SURFACE TYPES', '')
1916	ENDIF
1917
1918
1919	!! 5. computes total melt (snow and ice)
1920
1921
1922	DO ji = 1, kjpindex
1923	tot_melt(ji) = icemelt(ji) + snowmelt(ji)
1924	ENDDO
1925
1926	!! 6. computes snow age on veg and ice (for albedo)
1927
1928	DO ji = 1, kjpindex
1929
1930	!! 6.1 Snow age on vegetation
1931	IF (snow(ji) .LE. zero) THEN
1932	snow_age(ji) = zero
1933	ELSE
1934	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dt_sechiba/one_day) &
1935	& * EXP(-precip_snow(ji) / snow_trans)
1936	ENDIF
1937
1938	!! 6.2 Snow age on ice
1939	! age of snow on ice: a little bit different because in cold regions, we really
1940	! cannot negect the effect of cold temperatures on snow metamorphism any more.
1941	IF (snow_nobio(ji,iice) .LE. zero) THEN
1942	snow_nobio_age(ji,iice) = zero
1943	ELSE
1944
1945	d_age(ji) = ( snow_nobio_age(ji,iice) + &
1946	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dt_sechiba/one_day ) * &
1947	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
1948	IF (d_age(ji) .GT. zero ) THEN
1949	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
1950	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
1951	d_age(ji) = d_age(ji) / (un+xx(ji))
1952	ENDIF
1953	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
1954
1955	ENDIF
1956
1957	ENDDO
1958
1959	!! 7. Diagnose the depth of the snow layer
1960
1961	DO ji = 1, kjpindex
1962	snowdepth(ji) = snow(ji) /sn_dens
1963	ENDDO
1964
1965	IF (printlev>=3) WRITE (numout,*) ' hydrolc_snow done '
1966
1967	END SUBROUTINE hydrolc_snow
1968
1969
1970
1971	!! ================================================================================================================================
1972	!! SUBROUTINE : hydrolc_canop
1973	!!
1974	!>\BRIEF This routine computes the water budget of the canopy interception reservoir.
1975	!!
1976	!! DESCRIPTION : The first sequence is to withdraw interception loss from the interception
1977	!! reservoir, with independent values for each PFT. The reservoir loading happens afterwards.
1978	!! Rain falls uniformly over the PFTs. It uniformly loads the canopy interception
1979	!! reservoir of each PFT, up to the interception capacity, but a fraction of rainfall always falls to the ground.
1980	!! Throughfall is thus comprised of the fraction that always falls through, plus the remining fraction of rainfall that
1981	!! exceeds the interception capacity, the interception loss having been removed.
1982	!! \latexonly
1983	!! \input{interception.tex}
1984	!! \endlatexonly
1985	!!
1986	!! IMPORTANT NOTE : Bare soil treatment is complex in hydrolc :
1987	!! - veget_max(:,1) is the fraction of bare soil from the "annual" vegetation map
1988	!! - veget(:,1) is not the fraction of vegetation over bare soil but the total fraction of bare soil in the grid-cell
1989	!! (including the bare soil fractions over the vegetation PFTs).
1990	!! *** A diagram should be added in the spatial entry of the "umbrella"
1991	!! - For interception to be zero on bare soil, we thus need to impose qsintmax(:,1)=0
1992	!!
1993	!! RECENT CHANGE(S) : None
1994	!!
1995	!! MAIN OUTPUT VARIABLE(S) : precisol (throughfall), qsintveg (amount of water in the canopy interception reservoir)
1996	!!
1997	!! REFERENCE(S) : None
1998	!!
1999	!! FLOWCHART : None
2000	!! \n
2001	!_ ================================================================================================================================
2002
2003	SUBROUTINE hydrolc_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, tot_bare_soil, qsintveg, precisol)
2004
2005	!! 0. Variable and parameter declaration
2006
2007	!! 0.1 Input variables
2008
2009	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
2010	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall @tex ($kg m^{-2}$) @endtex
2011	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: vevapwet !! Interception loss over each PFT
2012	!! @tex ($kg m^{-2}$) @endtex
2013	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Grid-cell fraction effectively covered by
2014	!! vegetation for each PFT, except for soil
2015	!! (0-1, unitless)
2016	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum amount of water in the canopy
2017	!! interception reservoir @tex ($kg m^{-2}$) @endtex
2018	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
2019
2020	!! 0.2 Output variables
2021
2022	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Throughfall @tex ($kg m^{-2}$) @endtex
2023
2024	!! 0.3 Modified variables
2025
2026	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Amount of water in the canopy interception
2027	!! reservoir @tex ($kg m^{-2}$) @endtex
2028
2029	!! 0.4 Local variabless
2030
2031	INTEGER(i_std) :: ji, jv !! Grid-cell and PFT indices (unitless)
2032	REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew !! Temporary variable for intercepted water
2033	!! amount @tex ($kg m^{-2}$) @endtex
2034
2035	!_ ================================================================================================================================
2036
2037	! calcul de qsintmax a prevoir a chaque pas de temps
2038	! dans ini_sechiba
2039	! boucle sur les points continentaux
2040	! calcul de qsintveg au pas de temps suivant
2041	! par ajout du flux interception loss
2042	! calcule par enerbil en fonction
2043	! des calculs faits dans diffuco
2044	! calcul de ce qui tombe sur le sol
2045	! avec accumulation dans precisol
2046	! essayer d'harmoniser le traitement du sol nu
2047	! avec celui des differents types de vegetation
2048	! fait si on impose qsintmax ( ,1) = 0.0
2049	!
2050	! loop for continental subdomain
2051	!
2052
2053	!
2054	! 1. evaporation off the continents
2055	!
2056	! 1.1 The interception loss is take off the canopy.
2057
2058	DO jv=1,nvm
2059	qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv)
2060	END DO
2061
2062	! 1.2 It is raining : precip_rain is shared for each vegetation
2063	! type
2064	! sum (veget (1,nvm)) must be egal to 1-totfrac_nobio.
2065	! iniveget computes veget each day
2066	!
2067	DO jv=1,nvm
2068	! Correction Nathalie - Juin 2006 - une partie de la pluie arrivera toujours sur le sol
2069	! sorte de throughfall supplementaire
2070	!qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * precip_rain(:)
2071	!MM veget(:,1) BUG ??!!!!!!!!!!!
2072	IF (jv .EQ. 1) THEN
2073	qsintveg(:,jv) = qsintveg(:,jv) + tot_bare_soil(:) * ((1-throughfall_by_pft(jv))*precip_rain(:))
2074	ELSE
2075	qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
2076	ENDIF
2077	END DO
2078
2079	!
2080	! 1.3 Limits the effect and sum what receives soil
2081	!
2082	precisol(:,:) = zero
2083	DO jv=1,nvm
2084	DO ji = 1, kjpindex
2085	zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv))
2086	! correction throughfall Nathalie - Juin 2006
2087	!precisol(ji,jv) = qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
2088	!MM veget(:,1) BUG ??!!!!!!!!!!!
2089	IF (jv .EQ. 1) THEN
2090	precisol(ji,jv) = (tot_bare_soil(ji)throughfall_by_pft(jv)precip_rain(ji)) + qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
2091	ELSE
2092	precisol(ji,jv) = (veget(ji,jv)throughfall_by_pft(jv)precip_rain(ji)) + qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
2093	ENDIF
2094	ENDDO
2095	ENDDO
2096	!
2097	! 1.4 swap qsintveg to the new value
2098	!
2099
2100	DO jv=1,nvm
2101	qsintveg(:,jv) = zqsintvegnew (:,jv)
2102	END DO
2103
2104	IF (printlev>=3) WRITE (numout,*) ' hydrolc_canop done '
2105
2106	END SUBROUTINE hydrolc_canop
2107
2108
2109
2110	!! ================================================================================================================================
2111	!! SUBROUTINE : hydrolc_vegupd
2112	!!
2113	!>\BRIEF This subroutines works at adapting the distribution of water
2114	!! in the soil and interception reservoirs between the different soil columns when the vegetation
2115	!! fraction of the PFTs (veget) has changed in slowproc.
2116	!!
2117	!! DESCRIPTION : Different vegetation changes are allowed in ORCHIDEE:\n
2118	!! - veget_max can be updated annually in slowproc (dynamic vegetation or prescribed vegetation change)\n
2119	!! - veget can be updated daily in slowproc (because of changes in veget_max or depending on LAI)
2120	!!
2121	!! This subroutine aims at adapting the distribution of water among the different liquid water reservoirs
2122	!! (interception reservoir and soil moisture) after these changes of veget :
2123	!! - the variable veget holds the "new" vegetation map
2124	!! - the variable resdist holds the previous vegetation map
2125	!! *** I see no flag or time step dependance to control the call to vegupd : is it called at each time step ?
2126	!!
2127	!! The principle is that PFTs where "veget" has shrunk keep the same water contents in kg.m^{-2},
2128	!! whereas the water contents are increased in PFTs where "veget" has grown.
2129	!! *** I still have questions about redistribution to interception reservoirs which have shrunk
2130	!! *** and about thresholding to the reservoirs' capacities (is it ensured that the capacities are not exceeded ?)
2131	!! You may note that this occurs after evaporation and so on have been computed. It is not a
2132	!! problem as a new vegetation fraction will start with humrel=0 and thus will have no evaporation.
2133	!! If this is not the case it should have been caught above.
2134	!!
2135	!! IMPORTANT NOTE : the definition of veget is not simple :
2136	!! - for the non-bare soil PFTs (iv > 2), veget is the fraction effectively covered by
2137	!! vegetation : veget \f$\le\f$ veget_max\n
2138	!! - for the bare soil PFT (iv=1), veget holds the fraction of bare soil, from all
2139	!! PFTs : veget(:,1) \f$\ge\f$ veget_max(:,1)\n
2140	!!
2141	!! RECENT CHANGE(S) : None
2142	!!
2143	!! MAIN OUTPUT VARIABLE(S) : qsintveg, gqsb, bqsb, dsg, dss, dsp, resdist
2144	!!
2145	!! REFERENCE(S) : None
2146	!!
2147	!! FLOWCHART : None
2148	!! \n
2149	!_ ================================================================================================================================
2150
2151	SUBROUTINE hydrolc_vegupd(kjpindex, veget, tot_bare_soil, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss, dsp, resdist)
2152
2153	!! 0. Variable and parameter declaration
2154
2155	!! 0.1 Input variables
2156
2157	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
2158	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! Grid-cell fraction effectively covered by vegetation
2159	!! for each PFT, except for soil (0-1, unitless)
2160	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
2161	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Volumetric soil water capacity
2162	!! @tex ($kg m^{-3}$) @endtex
2163
2164	!! 0.2 Output variables
2165
2166	!! 0.3 Modified variables
2167
2168	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation due to interception
2169	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Water content in the top layer
2170	!! @tex ($kg m^{-2}$) @endtex
2171	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Water content in the bottom layer
2172	!! @tex ($kg m^{-2}$) @endtex
2173	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Depth of the top layer (m)
2174	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Depth of dry soil at the top, whether in the top or
2175	!! bottom layer (m)
2176	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Depth of dry soil in the bottom layer plus depth of
2177	!! top layer (m)
2178	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(inout) :: resdist !! Previous values of "veget" (0-1, unitless)
2179
2180	!! 0.4 Local variables
2181
2182	INTEGER(i_std) :: ji,jv !! Grid-cell and PFT indices (unitless)
2183	REAL(r_std),DIMENSION (kjpindex,nvm) :: qsintveg2 !! Ancillary qsintveg @tex ($kg m^{-2}$) @endtex
2184	REAL(r_std), DIMENSION (kjpindex,nvm) :: bdq, gdq !! Changes in surface and bottom soil layer water
2185	!! content @tex ($kg m^{-2}; positive$) @endtex
2186	REAL(r_std), DIMENSION (kjpindex,nvm) :: qsdq !! Changes in interception reservoir water content
2187	!! @tex ($kg m^{-2}; positive$) @endtex
2188	REAL(r_std), DIMENSION (kjpindex,nvm) :: vmr !! Variation of "veget" since previous values
2189	!! (-1,1; unitless)
2190	REAL(r_std), DIMENSION(kjpindex) :: gtr !! Total water mass to redistribute between the top
2191	!! layers of a grid-cell @tex ($kg m^{-2}; positive)
2192	!! @ endtex
2193	REAL(r_std), DIMENSION(kjpindex) :: btr !! Total water mass to redistribute between the bottom
2194	!! layers of a grid-cell @tex ($kg m^{-2}; positive)
2195	!! @ endtex
2196	REAL(r_std), DIMENSION(kjpindex) :: vtr !! Total fraction over which "veget" has decreased
2197	!! (dimensionless; positive)
2198	REAL(r_std), DIMENSION(kjpindex) :: qstr !! Total water mass to redistribute between the
2199	!! interception reservoirs of a grid-cell
2200	!! @tex ($kg m^{-2}; positive$) @endtex
2201	REAL(r_std), DIMENSION(kjpindex) :: fra !! Weight for redistribution (dimensionless; negative)
2202	REAL(r_std), DIMENSION(kjpindex) :: vegchtot !! Total change of "veget" in the grid-point
2203	!! @tex ($kg m^{-2}$) @endtex
2204	REAL(r_std), PARAMETER :: EPS1 = EPSILON(un) !! Very small (scalar) *** why not use min_sechiba ?
2205	!_ ================================================================================================================================
2206
2207	!! 1. vmr is the change in veget in the different PFTs
2208
2209	! vmr is the change in veget in the different PFTs
2210	! (dimensionless; negative if the new "veget" is smaller, i.e. if the PFT has shrunk)
2211	! By construction, the algebraic sum of vmr in a grid-cell is zero
2212	DO jv = 1, nvm
2213	DO ji = 1, kjpindex
2214	!MM veget(:,1) BUG ??!!!!!!!!!!!
2215	IF (jv .EQ. 1) THEN
2216	IF ( ABS(tot_bare_soil(ji)-resdist(ji,jv)) .GT. EPS1 ) THEN
2217	vmr(ji,jv) = tot_bare_soil(ji)-resdist(ji,jv)
2218	ELSE
2219	vmr(ji,jv) = zero
2220	ENDIF
2221	ELSE
2222	IF ( ABS(veget(ji,jv)-resdist(ji,jv)) .GT. EPS1 ) THEN
2223	vmr(ji,jv) = veget(ji,jv)-resdist(ji,jv)
2224	ELSE
2225	vmr(ji,jv) = zero
2226	ENDIF
2227	ENDIF
2228	!! 2. qsintveg2 is the intercepted water in mm
2229
2230	! qsintveg2 is the intercepted water in mmif the total volume
2231	! was redistributed over the previous "veget" fractions
2232	! This the distribution of intercepted water that needs to be
2233	! changed if "veget" changes
2234
2235	IF (resdist(ji,jv) .GT. zero) THEN
2236	qsintveg2(ji,jv) = qsintveg(ji,jv)/resdist(ji,jv)
2237	ELSE
2238	qsintveg2(ji,jv) = zero
2239	ENDIF
2240	ENDDO
2241	ENDDO
2242
2243
2244	!! 3. vegchtot is the total change of "veget" in the grid-points
2245
2246	! vegchtot is the total change of "veget" in the grid-points,
2247	! integrated over the PFTs (vegchtot in kg m^{-2})
2248	! It is the sum of the absolute values of changes, it may thus
2249	! be larger than 1 : it serves as a flag of change in each grid-point
2250	vegchtot(:) = zero
2251	DO jv = 1, nvm
2252	DO ji = 1, kjpindex
2253	vegchtot(ji) = vegchtot(ji) + ABS( vmr(ji,jv) )
2254	ENDDO
2255	ENDDO
2256
2257
2258	!! 4. In the grid-points with "veget" change, we define changes in water content
2259
2260	DO jv = 1, nvm
2261	DO ji = 1, kjpindex
2262	IF ( vegchtot(ji) .GT. zero ) THEN
2263
2264	! change in surface soil layer water content @tex ($kg m^{-2}$) @endtex
2265	gdq(ji,jv) = ABS(vmr(ji,jv)) * gqsb(ji,jv)
2266
2267	! change in bottom soil layer water content @tex ($kg m^{-2}$) @endtex
2268	bdq(ji,jv) = ABS(vmr(ji,jv)) * bqsb(ji,jv)
2269
2270	! change in interception reservoir water content @tex ($kg m^{-2}$) @endtex
2271	qsdq(ji,jv) = ABS(vmr(ji,jv)) * qsintveg2(ji,jv)
2272	ENDIF
2273	ENDDO
2274	ENDDO
2275
2276	!! 5. Total water mass to redistribute in a grid-point = sum of water changes from PFTs where "veget" decreases
2277
2278	! Calculate the total water mass that we need to redistribute by grid-point, for each water reservoir
2279	! We sum up the water changes from PFTs where "veget" decreases : this is the water that needs to be redistributed
2280	! vtr is the total fraction over which "veget" has decreased (> 0)
2281	! By construction, it is balanced by the total fraction where "veget" has increased
2282
2283	gtr(:) = zero
2284	btr(:) = zero
2285	qstr(:) = zero
2286	vtr(:) = zero
2287	!
2288	!
2289	DO jv = 1, nvm
2290	DO ji = 1, kjpindex
2291	IF ( ( vegchtot(ji) .GT. zero ) .AND. ( vmr(ji,jv) .LT. zero ) ) THEN
2292	gtr(ji) = gtr(ji) + gdq(ji,jv)
2293	btr(ji) = btr(ji) + bdq(ji,jv)
2294	qstr(ji) = qstr(ji) + qsdq(ji,jv)
2295
2296	! vtr is necessarily le 0 since we enter here only if vmr <0
2297	vtr(ji) = vtr(ji) - vmr(ji,jv)
2298	ENDIF
2299	ENDDO
2300	ENDDO
2301
2302
2303	!! 6. Put the water to redistribute from the PFTs that "shrank" into the PFTs that "growed"
2304
2305	! In doing so, water contents are kept constant in PFTs that shrank, and they increase in PFTs that growed
2306	! fra is the weight for that redistribution, scaled to vtr which is the total amount to redistribute
2307	! *** Feasability of the redistribution : it doesn't seem to be checked ???
2308	! *** In the soil, if the water capacities are constant between PFTs, it's OK ;this is the default for wmax_veg in
2309	! *** constantes_veg.f90
2310	! *** But qsintmax is different between PFTsand also evolves in time : how is this managed ???
2311	DO jv = 1, nvm
2312	DO ji = 1, kjpindex
2313
2314	! "veget" changed
2315	IF ( vegchtot(ji) .GT. zero .AND. ABS(vtr(ji)) .GT. EPS1) THEN
2316
2317	! negative when vmr positive, thus in the condition below
2318	fra(ji) = vmr(ji,jv) / vtr(ji)
2319
2320	! the PFT growed, thus its water contents must be updated
2321	IF ( vmr(ji,jv) .GT. zero) THEN
2322	!MM veget(:,1) BUG ??!!!!!!!!!!!
2323	IF (jv .EQ. 1) THEN
2324	IF (tot_bare_soil(ji) .GT. zero) THEN
2325	gqsb(ji,jv) = (resdist(ji,jv)gqsb(ji,jv) + fra(ji)gtr(ji))/tot_bare_soil(ji)
2326	bqsb(ji,jv) = (resdist(ji,jv)bqsb(ji,jv) + fra(ji)btr(ji))/tot_bare_soil(ji)
2327	ENDIF
2328	ELSE
2329	IF (veget(ji,jv) .GT. zero) THEN
2330	gqsb(ji,jv) = (resdist(ji,jv)gqsb(ji,jv) + fra(ji)gtr(ji))/veget(ji,jv)
2331	bqsb(ji,jv) = (resdist(ji,jv)bqsb(ji,jv) + fra(ji)btr(ji))/veget(ji,jv)
2332	ENDIF
2333	ENDIF
2334	qsintveg(ji,jv) = qsintveg(ji,jv) + fra(ji)* qstr(ji)
2335	ELSE
2336
2337	! vmr negative *** I don't understand why we remove water from the interception reservoir in PFTs which shrank
2338	qsintveg(ji,jv) = qsintveg(ji,jv) - qsdq(ji,jv)
2339	ENDIF
2340
2341	! Then we update the soil layers depths.
2342	! But we do not change dss, so that the redistribution does directly affect transpiration
2343	! constantes : min_sechiba = 1.E-8_r_std
2344	IF (gqsb(ji,jv) .LT. min_sechiba) THEN
2345	dsg(ji,jv) = zero
2346	ELSE
2347	dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) &
2348	/ ruu_ch(ji)
2349	ENDIF
2350	dsp(ji,jv) = zmaxh - bqsb(ji,jv) / ruu_ch(ji)
2351	ENDIF
2352	ENDDO
2353	ENDDO
2354
2355
2356	!! 7. We update resdist for the next "veget" change
2357
2358	!MM veget(:,1) BUG ??!!!!!!!!!!!
2359	resdist(:,1) = tot_bare_soil(:)
2360	DO jv = 2, nvm
2361	resdist(:,jv) = veget(:,jv)
2362	ENDDO
2363
2364
2365	!! 8. Where vegetation fraction is zero, set water to that of bare soil.
2366
2367	! This does not create any additional water.
2368	DO jv = 2, nvm
2369	DO ji = 1, kjpindex
2370	IF ( veget(ji,jv) .LT. EPS1 ) THEN
2371	gqsb(ji,jv) = gqsb(ji,1)
2372	bqsb(ji,jv) = bqsb(ji,1)
2373	dsg(ji,jv) = dsg(ji,1)
2374	dss(ji,jv) = dss(ji,1)
2375	dsp(ji,jv) = dsp(ji,1)
2376	ENDIF
2377	ENDDO
2378	ENDDO
2379
2380	END SUBROUTINE hydrolc_vegupd
2381
2382	!! ================================================================================================================================
2383	!! SUBROUTINE : hydrolc_flood
2384	!!
2385	!>\BRIEF this routine computes the evolution of the surface reservoir (floodplain)
2386	!!
2387	!! DESCRIPTION :
2388	!!
2389	!! RECENT CHANGE(S) : None
2390	!!
2391	!! MAIN OUTPUT VARIABLE(S) : floodout
2392	!!
2393	!! REFERENCE(S) : Same as for module hydrolc
2394	!!
2395	!! FLOWCHART : None
2396	!! \n
2397	!_ ================================================================================================================================
2398
2399	SUBROUTINE hydrolc_flood (kjpindex, vevapnu, vevapflo, flood_frac, flood_res, floodout)
2400	! input scalar
2401	INTEGER(i_std), INTENT(in) :: kjpindex
2402	! input fields
2403	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: flood_frac !! Fraction of floodplains in grid box
2404	! modified fields
2405	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: floodout !! Flux to take out from floodplains
2406	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: flood_res !! Floodplains reservoir estimate
2407	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation
2408	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapflo !! Floodplains evaporation
2409	! local declaration
2410	INTEGER(i_std) :: ji, jst, jv !! indices
2411	REAL(r_std) :: k_m !! conductivity in the soil
2412	REAL(r_std) :: temp !!
2413
2414	!_ ================================================================================================================================
2415
2416	!-
2417	!- 1. Take out vevapflo from the reservoir and transfer the remaining to vevapnu
2418	!-
2419	DO ji = 1,kjpindex
2420	temp = MIN(flood_res(ji), vevapflo(ji))
2421	flood_res(ji) = flood_res(ji) - temp
2422	vevapnu(ji) = vevapnu(ji) + vevapflo(ji) - temp
2423	vevapflo(ji) = temp
2424	ENDDO
2425
2426	!-
2427	!- 2. Compute the total flux from floodplain floodout (transfered to routing)
2428	!-
2429	DO ji = 1,kjpindex
2430	floodout(ji) = vevapflo(ji) - flood_frac(ji) * SUM(precisol(ji,:))
2431	ENDDO
2432
2433	!-
2434	!- 3. Discriminate between precip over land and over floodplain
2435	!-
2436	DO jv=1, nvm
2437	DO ji = 1,kjpindex
2438	precisol(ji,jv) = precisol(ji,jv) * (1 - flood_frac(ji))
2439	ENDDO
2440	ENDDO
2441
2442	IF (printlev>=3) WRITE (numout,*) ' hydrolc_flood done'
2443
2444	END SUBROUTINE hydrolc_flood
2445
2446	!! ================================================================================================================================
2447	!! SUBROUTINE : hydrolc_soil
2448	!!
2449	!>\BRIEF This routines computes the soil water budget and the related soil water
2450	!! fluxes and soil moisture diagnostics using the two-layer Choisnel scheme.
2451	!!
2452	!! DESCRIPTION :
2453	!!
2454	!! 1. Main processes: The Choisnel scheme relies on two soil layers.
2455	!! As show in the figure below, the upper one has a variable depth
2456	!! and can disappear after dry spells. It is created by infiltration (of throughfall or snow melt),
2457	!! in which case the layer is saturated. If this top layer is already present, infiltration can either
2458	!! fill it or make it deeper (DucoudrÃ© et al., 1993).
2459	!! \latexonly
2460	!! \includegraphics[scale = 0.5]{choisnelvariables.pdf}
2461	!! \endlatexonly
2462	!!
2463	!! In this framework, most water fluxes updating soil moisture act from top:\n
2464	!! - throughfall, melted snow and ice, and irrigation increase soil moisture (of the top layer if present);
2465	!! - transpiration and bare soil evaporation reduce soil moisture (of the top layer if present);
2466	!! - return flow from rivers increases bottom soil moisture. \n
2467	!!
2468	!! Soil moisture stress on bare soil evaporation and transpiration (for diffuco), vegetation growth and
2469	!! litter processes (for stomate), proceed respectively via a soil resistance (rsol), and three soil
2470	!! moisture stress factor (humrel, vegstress, and litterhumdiag), which are all controlled by the
2471	!! height of dry soil at the top of soil: \n
2472	!! - Soil moisture stress factor on transpiration and bare soil evaporation: humrel = \f$U_s\f$ \n
2473	!! \latexonly
2474	!! \input{humrel.tex}
2475	!! \endlatexonly \n
2476	!! - Resistance to bare soil evaporation: rsol = \f$r_\mathrm{soil}\f$ \n
2477	!! \latexonly
2478	!! \input{rsol.tex}
2479	!! \endlatexonly
2480	!!
2481	!! Runoff is only produced when the entire soil column is saturated, as in the bucket
2482	!! scheme of Manabe (1969). Drainage at the bottom of soil and surface runoff are only diagnosed,
2483	!! mostly for use in the routing module, using a constant 95% - 5% redistribution (Ngo-Duc et al., 2005).
2484	!! Internal drainage is allowed from the top to the bottom layer, following Ducharne et al. (1998).
2485	!! \latexonly
2486	!! \input{gdrainage.tex}
2487	!! \endlatexonly
2488	!!
2489	!! Irrigation (de Rosnay et al., 2003) and return flow are optional inflows, calculated by the
2490	!! routing scheme (Ngo-Duc et al., 2005; Guimberteau, 2010).
2491	!!
2492	!! 2. Subgrid variability:
2493	!! The subgrid variability of soil is described by introducing as many soil columns as PFTs
2494	!! (de Rosnay & Polcher, 1998). The water budget is performed independently in each PFT/soil
2495	!! column, but the bottom soil moisture is merged between all PFTs at the end of the calculations.
2496	!! Therefore, vegetation withdraws from a shared bottom moisture (grid-cell weighted average).
2497	!! There can also be horizontal diffusion between the top layers, if the flag ok_hdiff is true (the default value
2498	!! is false). This is performed in the subroutine hydrolc_hdiff, call after the present subroutine.
2499	!!
2500	!! The areas of each soil column are defined by veget (the vegetation fraction of the PFTs),
2501	!! and not by veget_max (the maximum area of the PFT defined annually in slowproc).
2502	!! As veget can change at a daily time step, a redistribution of water between the soil
2503	!! columns is required : this is done in hydrolc_vegupd.
2504	!!
2505	!! 3. Water budget issues :
2506	!! Different issues can arise when solving the above processes :
2507	!! - negative total soil moisture (tracked using flag warning, but no action taken, cf. 4)
2508	!! - dsg > dsp after merging the bottom layers, solved iteratively (6.2)
2509	!! - we may also have a pathological behavior of rsol if hdry > dpu_cste (8.4)
2510	!!
2511	!! 4. Diagnostics for other subroutines: humrel for diffuco (7.1), vegstress for stomate (8.1),
2512	!! shumdiag for stomate (8.2), drysoil_frac for condveg (8.3), rsol for diffuco (8.4),
2513	!! litterhumdiag for stomate (8.5).
2514	!!
2515	!! MAIN OUTPUT VARIABLE(S) : rsol, drysoil_frac, hdry,
2516	!! run_off_tot (surface runoff), drainage, humrel, vegstress, shumdiag, litterhumdiag
2517	!! gqsb, bqsb, dsg, dss, dsp
2518	!!
2519	!! REFERENCE(S) :
2520	!! - DucoudrÃ©, N, Laval, K & Perrier, A, 1993. SECHIBA, a new set of parameterisations
2521	!! of the hydrologic exchanges at the land-atmosphere interface within the LMD Atmospheric General
2522	!! Circulation Model. Journal of Climate, 6, pp. 248-273.
2523	!! - Ducharne, A. Laval, K. and Polcher, J. (1998) Sensitivity of the hydrological cycle to the parameterization
2524	!! of soil hydrology in a GCM. Climate Dynamics, 14:307-327.
2525	!! - de Rosnay, P. and Polcher J. (1998) Modeling root water uptake in a complex land surface scheme coupled
2526	!! to a GCM. Hydrology and Earth System Sciences, 2(2-3):239-256.
2527	!! - de Rosnay, P., Polcher, J., Laval, K. et Sabre, M. (2003). Integrated parameterization of irrigation in
2528	!! the land surface model ORCHIDEE. Validation over Indian Peninsula. Geophys. Res. Lett, 30(19):HLS2-1 -
2529	!! HLS2-4.
2530	!! - Ngo-Duc, T., J. Polcher, and K. Laval (2005), A 53-year forcing data set for land surface models,
2531	!! J. Geophys. Res., 110, D06116.
2532	!! - ALMA : http://www.lmd.jussieu.fr/~polcher/ALMA/
2533	!! - DucoudrÃ©, N. (1990). Sensibilite du climat simule a la parametrisation des echanges de vapeur d'eau entre
2534	!! la biosphere et l'atmosphere. These de Doctorat, UniversitÃ© Paris 6.
2535	!! - Ducharne A (1997). Le cycle de l'eau : modelisation de l'hydrologie continentale, etude de ses interactions
2536	!! avec le climat, These de Doctorat, UniversitÃ© Paris 6.
2537	!! - de Rosnay, P. (1999). Representation des interactions sol-plante-atmosphere dans le modele de circulation generale
2538	!! du LMD. These de doctorat, UniversitÃ© Paris 6.
2539	!! - Guimberteau, M, 2010. Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau,
2540	!! These de doctorat, UniversitÃ© Paris 6.
2541	!!
2542	!! FLOWCHART : None
2543	!! \n
2544	!_ ================================================================================================================================
2545
2546
2547	SUBROUTINE hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, reinfiltration, irrigation, tot_melt, mx_eau_var, &
2548	& veget, tot_bare_soil, ruu_ch, transpir,&
2549	& gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, run_off_tot, drainage, &
2550	& humrel, vegstress, shumdiag, litterhumdiag)
2551
2552	!! 0. Variable and parameter declaration
2553
2554	!! 0.1 Input variables
2555
2556	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size (number of grid cells) (unitless)
2557	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: vevapnu !! Bare soil evaporation @tex ($kg m^{-2}$) @endtex
2558	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: precisol !! Throughfall @tex ($kg m^{-2}$) @endtex
2559	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Routed water which comes back into the soil
2560	!! @tex ($kg m^{-2}$) @endtex
2561	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: reinfiltration !! Water returning to the top reservoir
2562	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: irrigation !! Irrigation water applied to soils
2563	!! @tex ($kg m^{-2}$) @endtex
2564	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt !! Total melt @tex ($kg m^{-2}$) @endtex
2565	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: mx_eau_var !! Maximum water content of the soil
2566	!! @tex ($kg m^{-2}$) @endtex
2567	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! Grid-cell fraction effectively covered by vegetation
2568	!! for each PFT, except for soil (0-1, unitless)
2569	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
2570	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Volumetric soil water capacity
2571	!! @tex ($kg m^{-3}$) @endtex
2572	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! Transpiration over each PFT @tex ($kg m^{-2}$) @endtex
2573
2574	!! 0.2 Output variables
2575
2576	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: run_off_tot !! Diagnosed surface runoff @tex ($kg m^{-2}$) @endtex
2577	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! Diagnosed drainage at the bottom of soil
2578	!! @tex ($kg m^{-2}$) @endtex
2579	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: humrel !! Soil moisture stress factor on transpiration and bare
2580	!! soil evaporation (0-1, unitless) ! Relative humidity
2581	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: vegstress !! Vegetation moisture stress (only for vegetation
2582	!! growth) (0-1, unitless)
2583	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! Mean relative soil moisture in the different levels
2584	!! used by thermosoil.f90 (0-1, unitless)
2585	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! Litter humidity factor (0-1, unitless), used in
2586	!! stomate
2587	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: rsol !! Resistance to bare soil evaporation (s m^{-1})
2588	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drysoil_frac !! Fraction of visible dry soil (0-1, unitless)
2589	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: hdry !! Mean top dry soil height (m) version beton
2590
2591	!! 0.3 Modified variables
2592
2593	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout):: gqsb !! Water content in the top layer
2594	!! @tex ($kg m^{-2}$) @endtex
2595	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout):: bqsb !! Water content in the bottom layer
2596	!! @tex ($kg m^{-2}$) @endtex
2597	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout):: dsg !! Depth of the top layer (m)
2598	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout):: dss !! Depth of dry soil at the top, whether in the top or
2599	!! bottom layer (m)
2600	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout):: dsp !! Depth of dry soil in the bottom layer plus depth of
2601	!! top layer (m)
2602
2603	!! 0.4 Local variables
2604
2605	INTEGER(i_std) :: ji,jv, jd !! Indices for grid-cells, PFTs and diagnostic levels in
2606	!! the soil (unitless)
2607	REAL(r_std), DIMENSION(kjpindex,nvm) :: runoff !! Total runoff @tex ($kg m^{-2}$) @endtex
2608	REAL(r_std), DIMENSION(kjpindex) :: zhumrel_lo !! Soil moisture stress factors on transpiration for the
2609	!! bottom and top layers (0-1, unitless)
2610	REAL(r_std), DIMENSION(kjpindex) :: zhumrel_up !! Soil moisture stress factors on transpiration for the
2611	!! bottom and top layers (0-1, unitless)
2612	REAL(r_std), DIMENSION(kjpindex,nvm) :: zeflux !! Evaporation to be withdrawn from soil
2613	!! @tex ($kg m^{-2}) @endtex ! *** why divided by veget ?
2614	REAL(r_std), DIMENSION(kjpindex,nvm) :: zpreci !! Throughfall @tex ($kg m^{-2}$) @endtex
2615	!! ! *** why divided by veget ?
2616	LOGICAL, DIMENSION(kjpindex,nvm) :: warning !! To track negative total soil moisture cases in soil
2617	!! columns (true/false)
2618	REAL(r_std) :: gtr, btr !! Fractions of top and botttom soil layers to the
2619	!! different diagnostic soil layers (0-1, unitless)
2620	REAL(r_std), DIMENSION(kjpindex) :: mean_dsg !! Mean depth of water in top soil layer (m)
2621	LOGICAL :: OnceMore !! To repeat the correction to solve dsg > dsp after
2622	!! diffusion between the bottom layers (true/false)
2623	INTEGER(i_std), PARAMETER :: nitermax = 100 !! Maximum number of iterations to dsg > dsp after
2624	!! diffusion between the bottom layers (unitless)
2625	INTEGER(i_std) :: niter !! Counter of iterations to solve dsg > dsp after
2626	!! diffusion between the bottom layers (unitless)
2627	INTEGER(i_std) :: nbad !! Number of negative total soil moisture cases
2628	!! (unitless); this is before diffusion between the
2629	!! bottom layers, cf. "warning" ***
2630	LOGICAL, DIMENSION(kjpindex,nvm) :: lbad !! Tags "bad" PFTs where dsg > dsp after diffusion
2631	!! between the bottom layers (true/false)
2632	LOGICAL, DIMENSION(kjpindex) :: lbad_ij !! Tags "bad" grid-points where at least one PFT exhibits
2633	!! dsg > dsp after diffusion between the bottom layers
2634	!! (true/false)
2635	REAL(r_std) :: gqseuil !! Ancillary variables to compute drainage between the
2636	!! two soil layers @tex ($kg m^{-2}$) @endtex
2637	REAL(r_std) :: eausup !! Ancillary variables to compute drainage between the
2638	!! two soil layers @tex ($kg m^{-2}$) @endtex
2639	REAL(r_std) :: wd1 !! Ancillary variables to compute drainage between the
2640	!! two soil layers @tex ($kg m^{-2}$) @endtex
2641	REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl !! Temporary diagnostic levels in the soil (m)
2642	REAL(r_std), DIMENSION(kjpindex,nvm) :: a_subgrd !! Diagnosed subgrid fraction of saturated soil in the
2643	!! top layer, to calculate hdry (0-1, unitless)
2644	!_ ================================================================================================================================
2645
2646	!! 1. Preliminary calculations
2647
2648	!! We define input and output fluxes at the soil surface, for each PFT
2649	!! The input flux is throughfall.
2650	!! The ouput flux is the total evaporation withdrawn from the soil (transpiration and bare soil evaporation, BSE)
2651	! * I don't understand why the fluxes are divided by veget if they are in kg.m^{-2} within each PFT *
2652	DO jv=1,nvm
2653	DO ji = 1, kjpindex
2654	!MM veget(:,1) BUG ??!!!!!!!!!!!
2655	IF (jv .EQ. 1) THEN
2656	IF ( tot_bare_soil(ji) .GT. zero ) THEN
2657	zeflux(ji,jv) = transpir(ji,jv)/tot_bare_soil(ji)
2658	zpreci(ji,jv) = precisol(ji,jv)/tot_bare_soil(ji)
2659	ELSE
2660	zeflux(ji,jv) = zero
2661	zpreci(ji,jv) = zero
2662	ENDIF
2663	ELSE
2664	IF ( veget(ji,jv) .GT. zero ) THEN
2665	zeflux(ji,jv) = transpir(ji,jv)/veget(ji,jv)
2666	zpreci(ji,jv) = precisol(ji,jv)/veget(ji,jv)
2667	ELSE
2668	zeflux(ji,jv) = zero
2669	zpreci(ji,jv) = zero
2670	ENDIF
2671	ENDIF
2672	ENDDO
2673	ENDDO
2674
2675	!! For bare soil evaporation, we need to distinguish two cases.
2676	!! This should only apply if there is vegetation but we do not test this case.
2677	!! Case 1 - if bare soil is present, BSE is extracted from the bare soil fraction
2678	DO ji = 1, kjpindex
2679	IF ( (vegtot(ji) .GT. zero) .AND. (tot_bare_soil(ji) .GT. min_sechiba) ) THEN
2680	zeflux(ji,1) = vevapnu(ji)/tot_bare_soil(ji)
2681	ENDIF
2682	ENDDO
2683
2684	!! Case 2 - if bare soil is not present, BSE is uniformly redistributed among the vegetation fractions
2685	!! This case is possible because of transfers (snow for instance).
2686	!MM veget(:,1) BUG ??!!!!!!!!!!!
2687	DO jv = 2, nvm
2688	DO ji = 1, kjpindex
2689	IF ( (vegtot(ji) .GT. zero) .AND. (tot_bare_soil(ji) .LE. min_sechiba)&
2690	& .AND. (veget(ji,jv) .GT. min_sechiba)) THEN
2691	zeflux(ji,jv) = zeflux(ji,jv) + vevapnu(ji)/vegtot(ji)
2692	ENDIF
2693	ENDDO
2694	ENDDO
2695
2696	! Temporary diagnostic levels in the soil to diagnose shumdiag
2697	tmp_dl(1) = 0
2698	tmp_dl(2:nbdl+1) = diaglev(1:nbdl)
2699
2700	!! 2. Updating soil moisture
2701
2702	!! We update soil moisture:
2703	!! The top soil moisture reacts to evaporation, throughfall, snow and ice melt, and inflow from irrigation
2704	!! The bottom soil moisture can receive returnflow from routing.f90
2705	DO jv=1,nvm
2706	DO ji=1,kjpindex
2707
2708	! Evaporated water is taken out of the ground
2709	gqsb(ji,jv) = gqsb(ji,jv) - zeflux(ji,jv)
2710	!
2711	! 1.2 Add snow and ice melt, troughfall from vegetation, reinfiltration and irrigation.
2712	!
2713	IF(vegtot(ji) .NE. zero) THEN
2714
2715	! snow and ice melt, reinfiltration and troughfall from vegetation
2716	gqsb(ji,jv) = gqsb(ji,jv) + zpreci(ji,jv) + (tot_melt(ji)+reinfiltration(ji))/vegtot(ji)
2717
2718	! We take care to add the irrigation only to the vegetated part if possible
2719	!
2720	IF (ABS(vegtot(ji)-tot_bare_soil(ji)) .LE. min_sechiba) THEN
2721	gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/vegtot(ji)
2722	ELSE
2723
2724	! no vegetation, : we only add irrigation in the PFTs where vegetation can grow
2725	IF ( jv > 1 ) THEN
2726
2727	! Only add the irrigation to the upper soil if there is a reservoir.
2728	! Without this the water evaporates right away.
2729	IF ( gqsb(ji,jv) > zero ) THEN
2730	gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-tot_bare_soil(ji))
2731	ELSE
2732	bqsb(ji,jv) = bqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-tot_bare_soil(ji))
2733	ENDIF
2734	ENDIF
2735	ENDIF
2736
2737	! We add the water returning from rivers to the lower reservoir.
2738	bqsb(ji,jv) = bqsb(ji,jv) + returnflow(ji)/vegtot(ji)
2739
2740	ENDIF
2741
2742	END DO
2743	ENDDO
2744
2745	!! 3. We compute runoff and adjust the soil layers' depth and water content
2746
2747	!! The depth of top dry soil, dss, is of particular importance as it controls the soil moisture stresses
2748	!! to transpiration and BSE
2749	runoff(:,:) = zero
2750
2751	warning(:,:) = .FALSE.
2752	DO jv=1,nvm
2753	DO ji = 1, kjpindex
2754
2755	!! 3.1 Soil moisture in excess of total water holding capacity runs off
2756	runoff(ji,jv) = MAX(gqsb(ji,jv) + bqsb(ji,jv) - mx_eau_var(ji), zero)
2757
2758	!! 3.2 If the soil is saturated (runoff is generated): no top layer; dss = dsp
2759	IF (mx_eau_var(ji) .LE. (gqsb(ji,jv) + bqsb(ji,jv))) THEN
2760	!
2761	gqsb(ji,jv) = zero
2762	dsg(ji,jv) = zero
2763	bqsb(ji,jv) = mx_eau_var(ji)
2764	dsp(ji,jv) = zero
2765	dss(ji,jv) = dsp (ji,jv)
2766
2767	ELSEIF ((gqsb(ji,jv) + bqsb(ji,jv)).GE.zero) THEN
2768
2769	!! 3.3 Else, if the top layer holds more water than allowed by its depth, the top layer deepens
2770	!! The top layer is saturated, so that dss=0.
2771	!! In this case, dsp is useless, and is not updated,
2772	!! even if bqsb has changed because of irrigation and return flow.
2773	IF (gqsb(ji,jv) .GT. dsg(ji,jv) * ruu_ch(ji)) THEN
2774	!
2775	dsg(ji,jv) = gqsb(ji,jv) / ruu_ch(ji)
2776	dss(ji,jv) = zero
2777	ELSEIF (gqsb(ji,jv) .GT. zero ) THEN
2778
2779	!! 3.4 If the top layer is not saturated, its total depth does not change and we only update its dry soil depth
2780	!! In this case, dsp is useless, and is not updated,
2781	!! even if bqsb has changed because of irrigation and return flow.
2782	!
2783	dss(ji,jv) = ((dsg(ji,jv) * ruu_ch(ji)) - gqsb(ji,jv)) / ruu_ch(ji)
2784	ELSE
2785
2786	!! 3.5 If the top layer's moisture is negative, it vanishes and the required moisture is withdrawn from the bottom layer
2787	!! dsp is updated accordingly, and defines the depth of top dray soil dss.
2788	!
2789	bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv)
2790	dsp(ji,jv) = zmaxh - bqsb(ji,jv) / ruu_ch(ji) ! dsp>dpu if bqsb<0 *** we can be here with bsqb<0 and bqsb+gqsb>0
2791	gqsb(ji,jv) = zero
2792	dsg(ji,jv) = zero
2793	dss(ji,jv) = dsp(ji,jv)
2794	END IF
2795
2796	ELSE
2797
2798	!! 3.6 If the total soil moisture is negative: we keep track of the problem for further warning.
2799	!! In this extreme case of dry soil, the top layer is absent.
2800	!! The top dry soil depth, dsp, is equal to the soil depth.
2801	!! But we keep the negative moisture for water budget closure, and assign it to the bottom layer.
2802	! Ceci ne devrait jamais arriver plus d'une fois par point. C-a-d une fois la valeur negative
2803	! atteinte les flux doivent etre nuls. On ne signale que ce cas donc.
2804	! * But this is obviously not the case in CMIP5 runs *
2805	! *** Also, I don't see the use of conditionning upon zeflux(ji,jv) .GT. zero : negative moisture is always a problem !
2806
2807	IF ( ( zeflux(ji,jv) .GT. zero ) .AND. &
2808	( gqsb(ji,jv) + bqsb(ji,jv) .LT. -1.e-10 ) ) THEN
2809	warning(ji,jv) = .TRUE.
2810
2811	! WRITE (numout,*) 'WARNING! Soil Moisture will be negative'
2812	! WRITE (numout,*) 'ji, jv = ', ji,jv
2813	! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji)
2814	! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv)
2815	! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv)
2816	! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv)
2817	! WRITE (numout,*) 'dss = ', dss(ji,jv)
2818	! WRITE (numout,*) 'dsg = ', dsg(ji,jv)
2819	! WRITE (numout,*) 'dsp = ', dsp(ji,jv)
2820	! WRITE (numout,*) 'humrel = ', humrel(ji, jv)
2821	! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv)
2822	! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji)
2823	! WRITE (numout,*) '============================'
2824	ENDIF
2825
2826
2827	bqsb(ji,jv) = gqsb(ji,jv) + bqsb(ji,jv) ! this will be negative
2828	dsp(ji,jv) = zmaxh
2829	gqsb(ji,jv) = zero
2830	dsg(ji,jv) = zero
2831	dss(ji,jv) = dsp(ji,jv)
2832
2833	ENDIF
2834
2835	ENDDO
2836	ENDDO
2837
2838	!! 4. If there are some PFTs with negative moisture, it is written in the run log
2839
2840	nbad = COUNT( warning(:,:) .EQV. .TRUE. )
2841	IF ( nbad .GT. 0 ) THEN
2842	WRITE(numout,*) 'hydrolc_soil: WARNING! Soil moisture was negative at', &
2843	nbad, ' points:'
2844	!DO jv = 1, nvm
2845	! DO ji = 1, kjpindex
2846	! IF ( warning(ji,jv) ) THEN
2847	! WRITE(numout,*) ' ji,jv = ', ji,jv
2848	! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji)
2849	! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv)
2850	! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv)
2851	! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv)
2852	! WRITE (numout,*) 'runoff = ',runoff(ji,jv)
2853	! WRITE (numout,*) 'dss = ', dss(ji,jv)
2854	! WRITE (numout,*) 'dsg = ', dsg(ji,jv)
2855	! WRITE (numout,*) 'dsp = ', dsp(ji,jv)
2856	! WRITE (numout,*) 'humrel = ', humrel(ji, jv)
2857	! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv)
2858	! WRITE (numout,*) 'Soil precipitation = ',zpreci(ji,jv)
2859	! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji)
2860	! WRITE (numout,*) 'returnflow = ',returnflow(ji)
2861	! WRITE (numout,*) '============================'
2862	! ENDIF
2863	! ENDDO
2864	!ENDDO
2865	ENDIF
2866
2867	!! 5. Top layers that are very large or very dry can be deadlock situations for the Choisnel scheme
2868
2869	!! Top layers that are very large or that hold a tiny amount of water .
2870	!! They are handled here.
2871	IF (printlev>=3) WRITE(numout,*) 'hydro_soil 2.0 : Resolve deadlocks'
2872	DO jv=1,nvm
2873	DO ji=1,kjpindex
2874
2875	!! 5.1 If the top layer is almost dry, we merge the two layers
2876	IF ( ABS(dsp(ji,jv)-dsg(ji,jv)) .LT. min_sechiba ) THEN ! min_sechiba = 1.e-8 (constantes.f90)
2877	bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv)
2878	dsp(ji,jv) = zmaxh - bqsb(ji,jv) / ruu_ch(ji) ! *** can this make dsp > dpu_cste
2879	gqsb(ji,jv) = zero
2880	dsg(ji,jv) = zero
2881	dss(ji,jv) = dsp(ji,jv)
2882	ENDIF
2883
2884	!! 5.2 Internal drainage from the top to the bottom layer
2885	!! Much faster when the relative moisture of the top layer exceeds 75\% of the water holding capacity.
2886	!! The parameters are exp_drain = 1.5, min_drain = 0.002 mm/h and max_drain = 0.2 mm/h
2887	gqseuil = min_sechiba * deux * ruu_ch(ji) ! min_sechiba = 1.e-8 (constantes.f90)
2888	eausup = dsg(ji,jv) * ruu_ch(ji)
2889	wd1 = .75*eausup
2890
2891	IF (eausup .GT. gqseuil) THEN ! dsg > 2.e-8 m
2892	gdrainage(ji,jv) = min_drain * (gqsb(ji,jv)/eausup)
2893
2894	IF ( gqsb(ji,jv) .GE. wd1 .AND. dsg(ji,jv) .GT. 0.10 ) THEN
2895
2896	gdrainage(ji,jv) = gdrainage(ji,jv) + &
2897	(max_drain-min_drain)((gqsb(ji,jv)-wd1) / (eausup-wd1))*exp_drain
2898
2899	ENDIF
2900
2901	gdrainage(ji,jv)=MIN(gdrainage(ji,jv), MAX(gqsb(ji,jv), zero))
2902
2903	ELSE
2904	gdrainage(ji,jv)=zero ! *** why not remove the top layer in that case ??
2905	ENDIF
2906	!
2907	gqsb(ji,jv) = gqsb(ji,jv) - gdrainage(ji,jv)
2908	bqsb(ji,jv) = bqsb(ji,jv) + gdrainage(ji,jv)
2909	dsg(ji,jv) = dsg(ji,jv) - gdrainage(ji,jv) / ruu_ch(ji)
2910	dsp(ji,jv) = zmaxh - bqsb(ji,jv)/ruu_ch(ji) ! *** this line shouldn't make dsp>dpu_cste
2911	ENDDO
2912	ENDDO
2913
2914	!! 6. We want the vegetation to share a common bottom moisture:
2915	!! Thus, we must account for moisture diffusion between the soil columns.
2916	IF (printlev>=3) WRITE(numout,*) 'hydrolc_soil 3.0 : Vertical diffusion'
2917
2918	!! 6.1 Average of bottom and top moisture over the "veget" fractions
2919	mean_bqsb(:) = zero
2920	mean_gqsb(:) = zero
2921	DO jv = 1, nvm
2922	DO ji = 1, kjpindex
2923	IF ( vegtot(ji) .GT. zero ) THEN
2924	!MM veget(:,1) BUG ??!!!!!!!!!!!
2925	IF (jv .EQ. 1) THEN
2926	mean_bqsb(ji) = mean_bqsb(ji) + tot_bare_soil(ji)/vegtot(ji)*bqsb(ji,jv)
2927	mean_gqsb(ji) = mean_gqsb(ji) + tot_bare_soil(ji)/vegtot(ji)*gqsb(ji,jv)
2928	ELSE
2929	mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv)
2930	mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv)
2931	ENDIF
2932	ENDIF
2933	ENDDO
2934	ENDDO
2935
2936	OnceMore = .TRUE.
2937	niter = 0
2938	lbad_ij(:)=.TRUE.
2939	DO WHILE ( OnceMore .AND. ( niter .LT. nitermax ) ) ! nitermax prevents infinite loops (should actually never occur)
2940
2941	niter = niter + 1
2942
2943	!! 6.2 We assign the mean value in all the soil columns in a grid-cell
2944	DO jv = 1, nvm
2945	DO ji = 1, kjpindex
2946	IF (lbad_ij(ji)) THEN
2947	!MM veget(:,1) BUG ??!!!!!!!!!!!
2948	IF (jv .EQ. 1) THEN
2949	IF (tot_bare_soil(ji).GT.zero) THEN
2950	!
2951	bqsb(ji,jv) = mean_bqsb(ji)
2952	dsp(ji,jv) = zmaxh - bqsb(ji,jv)/ruu_ch(ji)
2953	ENDIF
2954	ELSE IF ( veget(ji,jv) .GT. zero ) THEN
2955	!
2956	bqsb(ji,jv) = mean_bqsb(ji)
2957	dsp(ji,jv) = zmaxh - bqsb(ji,jv)/ruu_ch(ji) ! *** can this line make dsp>dpu_cste ?
2958	ENDIF
2959	ENDIF
2960
2961	ENDDO
2962	ENDDO
2963
2964	!! 6.3 Iterative adjustment if top layer larger than new average dsp
2965	!! After averaging the bottom moistures, dsp becomes the mean depth of soil that is not filled by bottom moisture.
2966	!! In "bad" points where dsg > dsp, we merge the two soil layers.
2967	!! This adjustment needs to be done iteratively (WHILE loop)
2968	!! We diagnose the bad points where dsg>dsp
2969	!MM veget(:,1) BUG ??!!!!!!!!!!!
2970	!!$ lbad(:,:) = ( ( dsp(:,:) .LT. dsg(:,:) ) .AND. &
2971	!!$ ( dsg(:,:) .GT. zero ) .AND. &
2972	!!$ ( veget(:,:) .GT. zero ) )
2973	lbad(:,:) = ( ( dsp(:,:) .LT. dsg(:,:) ) .AND. &
2974	( dsg(:,:) .GT. zero ) )
2975	DO jv = 1, nvm
2976	DO ji = 1, kjpindex
2977	IF (jv .EQ. 1) THEN
2978	lbad(ji,jv) = lbad(ji,jv) .AND. ( tot_bare_soil(ji) .GT. zero )
2979	ELSE
2980	lbad(ji,jv) = lbad(ji,jv) .AND. ( veget(ji,jv) .GT. zero )
2981	ENDIF
2982	ENDDO
2983	ENDDO
2984
2985
2986	!! If there are no such points, we'll do no further iteration
2987	IF ( COUNT( lbad(:,:) ) .EQ. 0 ) OnceMore = .FALSE.
2988	DO ji = 1, kjpindex
2989	IF (COUNT(lbad(ji,:)) == 0 ) lbad_ij(ji)=.FALSE.
2990	ENDDO
2991
2992	!! In the bad PFTs, we merge the two soil layers.
2993	DO jv = 1, nvm
2994	!YM
2995	! !
2996	! DO ji = 1, kjpindex
2997	! IF ( veget(ji,jv) .GT. zero ) THEN
2998	! !
2999	! bqsb(ji,jv) = mean_bqsb(ji)
3000	! dsp(ji,jv) = zmaxh - bqsb(ji,jv)/ruu_ch(ji)
3001	! ENDIF
3002	! !
3003	! ENDDO
3004	!
3005	DO ji = 1, kjpindex
3006	IF ( lbad(ji,jv) ) THEN
3007	!
3008	runoff(ji,jv) = runoff(ji,jv) + &
3009	MAX( bqsb(ji,jv) + gqsb(ji,jv) - mx_eau_var(ji), zero)
3010	!
3011	bqsb(ji,jv) = MIN( bqsb(ji,jv) + gqsb(ji,jv), mx_eau_var(ji))
3012	!
3013	gqsb(ji,jv) = zero
3014	dsp(ji,jv) = zmaxh - bqsb(ji,jv)/ruu_ch(ji) ! *** could this line make dsp>dpu_cste ?
3015	! *** set a max to dpu_cste to be sure !
3016	dss(ji,jv) = dsp(ji,jv)
3017	dsg(ji,jv) = zero
3018
3019	ENDIF
3020	ENDDO
3021	ENDDO
3022
3023	!! New average of bottom and top moisture over the "veget" fractions, for the next iteration
3024	mean_bqsb(:) = zero
3025	mean_gqsb(:) = zero
3026	DO jv = 1, nvm
3027	DO ji = 1, kjpindex
3028	IF ( vegtot(ji) .GT. zero ) THEN
3029	!MM veget(:,1) BUG ??!!!!!!!!!!!
3030	IF (jv .EQ. 1) THEN
3031	mean_bqsb(ji) = mean_bqsb(ji) + tot_bare_soil(ji)/vegtot(ji)*bqsb(ji,jv)
3032	mean_gqsb(ji) = mean_gqsb(ji) + tot_bare_soil(ji)/vegtot(ji)*gqsb(ji,jv)
3033	ELSE
3034	mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv)
3035	mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv)
3036	ENDIF
3037	ENDIF
3038	ENDDO
3039	ENDDO
3040
3041	ENDDO ! end while, but no warning if nitermax has been reached ***
3042
3043	!! 7. Compute total runoff from all soil columns and partition it into drainage and surface runoff
3044
3045	! *** why isn't run_off_tot divided by vegtot ?
3046	IF (printlev>=3) WRITE(numout,*) 'hydrolc_soil 4.0: Computes total runoff'
3047
3048	!! 7.1 Average total runoff
3049	run_off_tot(:) = zero
3050	DO ji = 1, kjpindex
3051	IF ( vegtot(ji) .GT. zero ) THEN
3052	run_off_tot(ji) = runoff(ji,1)tot_bare_soil(ji) + SUM(runoff(ji,2:nvm)veget(ji,2:nvm))
3053	ELSE
3054	run_off_tot(ji) = tot_melt(ji) + irrigation(ji) + reinfiltration(ji)
3055	ENDIF
3056	ENDDO
3057
3058	!! 7.2 Diagnose drainage and surface runoff (95% and 5%)
3059	drainage(:) = 0.95 * run_off_tot(:)
3060	run_off_tot(:) = run_off_tot(:) - drainage(:)
3061	! *** from now on, we lost track of total runoff
3062	! *** the name "run_off_tot" is VERY misleading
3063
3064	!! 8. Diagnostics which are needed to carry information to other modules:
3065
3066	IF (printlev>=3) WRITE(numout,*) 'hydro_soil 5.0: Diagnostics'
3067
3068	! Reset dsg if necessary
3069	WHERE (gqsb(:,:) .LE. zero) dsg(:,:) = zero
3070
3071	! Average moisture profile for shumdiag
3072	DO ji=1,kjpindex
3073	mean_dsg(ji) = mean_gqsb(ji)/ruu_ch(ji) ! mean depth of water in top layer in meters
3074	ENDDO
3075
3076	!! 8.1 Moisture stress factors on vegetation (humrel and vegstress)
3077	!! Moisture stress factors on vegetation:
3078	!! - "humrel" on transpiration (for condveg) exponentially decreases with top dry soil depth;
3079	!! the decay factor is the one of root density with depth (from 1 to 8 depending on PFTs),
3080	!! following de Rosnay and Polcher (1998)
3081	!! - "vegstress" on vegetation growth (for stomate)
3082	IF (printlev>=3) WRITE(numout,*) 'hydro_soil 6.0 : Moisture stress'
3083
3084	a_subgrd(:,:) = zero
3085
3086	DO jv = 1, nvm
3087	DO ji=1,kjpindex
3088	!
3089	! computes relative surface humidity
3090	!
3091	! Only use the standard formulas if total soil moisture is larger than zero.
3092	! Else stress functions are set to zero.
3093	! This will avoid that large negative soil moisture accumulate over time by the
3094	! the creation of small skin reservoirs which evaporate quickly.
3095	!
3096	IF ( gqsb(ji,jv)+bqsb(ji,jv) .GT. zero ) THEN
3097	!
3098	IF (dsg(ji,jv).EQ. zero .OR. gqsb(ji,jv).EQ.zero) THEN
3099	humrel(ji,jv) = EXP( - humcste(jv) * zmaxh * (dsp(ji,jv)/zmaxh) )
3100	dsg(ji,jv) = zero
3101
3102	! humrel = 0 if dsp is larger than its value at the wilting point, or if the bottom layer's soil is negative
3103	! *** the condition based on qwilt (= 5.0) doesn't make much sense: (i) the Choisnel scheme works in available
3104	! *** moisture, i.e. above wilting point (Ducharne & Laval, 2000), so that the moisture at withing point is ZERO;
3105	! *** (ii) with the chosen values in constantes_soil.f90, qwilt/ruu_ch is around 0.033 and dpu_cste-qwilt/ruu_ch
3106	! *** is around dpu
3107	IF (dsp(ji,jv).GT.(zmaxh - min_sechiba) .OR. bqsb(ji,jv).LT.zero) THEN
3108	humrel(ji,jv) = zero
3109	ENDIF
3110
3111	vegstress(ji,jv) = humrel(ji,jv)
3112
3113	ELSE
3114
3115	!! 8.1.2 If there are two soil layers, we need the "transpiring" fraction "a_subgrd"
3116	!! We compute humrel as the average of the values defined by dss and dsp.
3117	!! The weights depend on "a_subgrd", which is defined by redistributing the top layer's
3118	!! moisture in a virtual layer of depth dsg_min (set to 1 mm in hydrolc),
3119	!! what defines a totally dry and a totally saturated fraction (a_subgrd):
3120	!! - if the top layer's moisture > 1mm, then a_subgrd=1, and the humrel is normally deduced from dss only
3121	!! - if the top layer's moisture is very small, a_subgrd decreases from 1 to 0 as the top layer's moisture
3122	!! vanishes, and it serves to describe a smooth transition to the case where the top soil layer has
3123	!! vanished and humrel depends on dsp.
3124	!! \latexonly
3125	!! \includegraphics[scale = 1]{asubgrid.pdf}
3126	!! \endlatexonly
3127	!
3128	zhumrel_lo(ji) = EXP( - humcste(jv) * dsp(ji,jv))
3129	zhumrel_up(ji) = EXP( - humcste(jv) * dss(ji,jv))
3130	a_subgrd(ji,jv)=MIN(MAX(dsg(ji,jv)-dss(ji,jv),zero)/dsg_min,un)
3131	humrel(ji,jv)=a_subgrd(ji,jv)zhumrel_up(ji)+(un-a_subgrd(ji,jv))zhumrel_lo(ji)
3132
3133	!! As we need a slower variable for vegetation growth, vegstress is computed differently from humrel.
3134	! *** la formule ci-dessous est d'un empirisme absolu : on ajoute deux stresses, on en retranche un autre...????
3135	! que veut dire slower variable ??
3136	vegstress(ji,jv) = zhumrel_lo(ji) + zhumrel_up(ji) - EXP( - humcste(jv) * dsg(ji,jv) )
3137
3138	ENDIF
3139
3140	ELSE
3141
3142	!! 8.1.3 If total moisture is negative, the two moisture stress factors are set to zero.
3143	!! This should avoid that large negative soil moisture accumulate over time because of small skin layers
3144	!! which transpire quickly.
3145	! *** Yet, CMIP5 exhibits such behaviors => because of rsol ?? or other reasons ?? The routing shouldn't
3146	! *** create such situations, but couldn't some patches here or there ???
3147	humrel(ji,jv) = zero
3148	vegstress(ji,jv) = zero
3149
3150	ENDIF
3151	!
3152	ENDDO
3153	ENDDO
3154
3155	!! 8.2 Mean relative soil moisture in the diagnostic soil levels: shumdiag (for thermosoil)
3156	!! It is deduced from the mean moisture and depth of the two soil layers in the grid cell,
3157	!! depending on how the soil diagnostic levels intersect the two mean soil depths
3158	DO jd = 1,nbdl
3159	DO ji = 1, kjpindex
3160	IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN ! mean_dsg = mean depth of water in top layer in meters
3161
3162	! If the diagnostic soil level entirely fits into the mean top soil layer depth, they have the same
3163	! relative moisture
3164	shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji)
3165	ELSE
3166	IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN
3167
3168	! If the diagnostic soil level intersects both soil layers, its relative moisture is the weighted
3169	! mean of the ones of two soil layers
3170	gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd))
3171	btr = 1 - gtr
3172	shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + &
3173	& btr*mean_bqsb(ji)/mx_eau_var(ji)
3174	ELSE
3175
3176	! If the diagnostic soil level entirely fits into the mean bottom soil layer depth,
3177	! they have the same relative moisture
3178	shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji)
3179	ENDIF
3180	ENDIF
3181	shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero)
3182	ENDDO
3183	ENDDO
3184
3185	!! 8.3 Fraction of visibly dry soil in the bare soil fraction for soil albedo
3186	!! if we want to account for its dependance on soil moisture in condveg.
3187	!! We redistribute the top layer's moisture in a virtual layer of depth 0.1 m,
3188	!! what defines a totally saturated and a totally dry fraction (drysoil_frac).
3189	!! The latter is thus 1 if dss > 0.1 m.
3190	drysoil_frac(:) = MIN(MAX(dss(:,1),zero)*10._r_std, un)
3191
3192	!! 8.4 Resistance to bare soil evaporation
3193	!! This resistance increases with the height of top dry soil, described by hdry.
3194	!! It depends on both dss and dsp (dry soil heights) in the bare soil fraction,
3195	!! and, as for the soil moisture stress factor on transpiration (humrel),
3196	!! we use "a_subgrd" (see 8.1.2) to describe a smooth transition between cases
3197	!! when the top layer's moisture > 1mm and hdry=dss and cases when the top layer's
3198	!! has vanished and hdy=dsp.
3199	!! The relationship between rsol and dry soil height has been changed since
3200	!! de Rosnay and Polcher (1998), to make rsol becomes VERY large when hdry
3201	!! gets close to dpu_cste
3202	!! Owing to this non-linear dependance, BSE tends to zero when the depth of dry reaches dpu_cste.
3203	! *** if hdry > dpu (several cases might lead to this pathological situation),
3204	! *** we shift to the other limb of the hyperbolic function, and rsol decreases towards hdry*rsol_cste
3205	! *** can this explain the accumulation of very negative soil moisture in arid zones in CMIP5 ???
3206	! *** COULD'NT WE SIMPLY SET THAT hdry CANNOT EXCEED dpu_cste ???
3207	hdry(:) = a_subgrd(:,1)dss(:,1) + (un-a_subgrd(:,1))dsp(:,1)
3208
3209	rsol(:) = -un
3210	DO ji = 1, kjpindex
3211	IF (tot_bare_soil(ji) .GE. min_sechiba) THEN
3212
3213	! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
3214	! on modifie le rsol pour que la resistance croisse subitement si on s'approche
3215	! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70
3216	!rsol(ji) = dss(ji,1) * rsol_cste
3217	rsol(ji) = ( hdry(ji) + un/(10.(zmaxh - hdry(ji))+1.e-10)2 ) rsol_cste
3218	ENDIF
3219	ENDDO
3220
3221	!! 8.5 Litter humidity factor, used in stomate
3222	!! It varies from 1 when mean top dry soil height is 0, and decreases as mean top dry soil height increases
3223	litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter ) ! hcrit_litter=0.08_r_std
3224
3225	!! Special case of it has just been raining a few drops: the top layer exists, but its height
3226	!! is ridiculously small and the mean top dry soil height is almost zero: we assume a dry litter
3227	WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. & ! constantes : min_sechiba = 1.E-8_r_std
3228	( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) )
3229	litterhumdiag(:) = zero
3230	ENDWHERE
3231
3232	IF (printlev>=3) WRITE (numout,*) ' hydrolc_soil done '
3233
3234	END SUBROUTINE hydrolc_soil
3235
3236
3237	SUBROUTINE hydrolc_waterbal_init (kjpindex, qsintveg, snow, snow_nobio)
3238
3239	!! 0. Variable and parameter declaration
3240	!! 0.1 Input variables
3241	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
3242	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Amount of water in the canopy interception
3243	!! reservoir @tex ($kg m^{-2}$) @endtex
3244	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
3245	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !! Snow water equivalent on nobio areas
3246	!! @tex ($kg m^{-2}$) @endtex
3247	!! 0.4 Local variables
3248	INTEGER(i_std) :: ji, jv, jn
3249	REAL(r_std),DIMENSION (kjpindex) :: watveg !! Total amount of intercepted water in a grid-cell
3250	REAL(r_std),DIMENSION (kjpindex) :: sum_snow_nobio !! Total amount of snow from the "nobio" fraction
3251	!! in a grid-cell @tex ($kg m^{-2}$) @endtex
3252	!_ ================================================================================================================================
3253
3254	!! 1. We initialize the total amount of water in terrestrial grid-cells at the beginning of the first time step
3255
3256	tot_water_beg(:) = zero
3257	watveg(:) = zero
3258	sum_snow_nobio(:) = zero
3259
3260	DO jv = 1, nvm
3261	watveg(:) = watveg(:) + qsintveg(:,jv)
3262	ENDDO
3263
3264	DO jn = 1, nnobio ! nnobio=1
3265	sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn)
3266	ENDDO
3267
3268	DO ji = 1, kjpindex
3269	tot_water_beg(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + &
3270	& watveg(ji) + snow(ji) + sum_snow_nobio(ji)
3271	ENDDO
3272	tot_water_end(:) = tot_water_beg(:)
3273
3274	END SUBROUTINE hydrolc_waterbal_init
3275
3276	!! ================================================================================================================================
3277	!! SUBROUTINE : hydrolc_waterbal
3278	!!
3279	!>\BRIEF This subroutine checks the water balance closure in each terrestrial grid-cell
3280	!! over a time step.
3281	!!
3282	!! DESCRIPTION : The change in total water amount over the time step must equal the algebraic sum of the
3283	!! water fluxes causing the change, integrated over the timestep. The equality is checked to some precision,
3284	!! set for double precision calculations (REAL*8). Note that this precision depends on the time step.
3285	!! This verification does not make much sense in REAL*4 as the precision is the same as some of the fluxes.
3286	!! The computation is only done over the soil area, as over glaciers (and lakes?) we do not have
3287	!! water conservation.
3288	!! *** The above sentence is not consistent with what I understand from the code, since snow_nobio is accounted for.
3289	!! *** what is the value of epsilon(1) ?
3290	!!
3291	!! RECENT CHANGE(S) : None
3292	!!
3293	!! MAIN OUTPUT VARIABLE(S) : None
3294	!!
3295	!! REFERENCE(S) : None
3296	!!
3297	!! FLOWCHART : None
3298	!! \n
3299	!_ ================================================================================================================================
3300
3301	SUBROUTINE hydrolc_waterbal (kjpindex, index, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
3302	& precip_rain, precip_snow, returnflow, reinfiltration, irrigation, tot_melt, vevapwet, transpir, vevapnu,&
3303	& vevapsno, vevapflo, floodout, run_off_tot, drainage)
3304
3305	!! 0. Variable and parameter declaration
3306
3307	!! 0.1 Input variables
3308
3309	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
3310	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the points on the map
3311	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by
3312	!! vegetation for each PFT, except for soil
3313	!! (0-1, unitless)
3314	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of terrestrial ice+lakes+...
3315	!! (0-1, unitless)
3316	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Amount of water in the canopy interception
3317	!! reservoir @tex ($kg m^{-2}$) @endtex
3318	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
3319	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout):: snow_nobio !! Snow water equivalent on nobio areas
3320	!! @tex ($kg m^{-2}$) @endtex
3321	! Ice water balance
3322	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rainfall @tex ($kg m^{-2}$) @endtex
3323	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snowfall @tex ($kg m^{-2}$) @endtex
3324	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil
3325	!! @tex ($kg m^{-2}$) @endtex
3326	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: reinfiltration !! Water returning from routing to the top reservoir
3327	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Irrigation water applied to soils
3328	!! @tex ($kg m^{-2}$) @endtex
3329	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt @tex ($kg m^{-2}$) @endtex
3330	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss over each PFT
3331	!! @tex ($kg m^{-2}$) @endtex
3332	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration over each PFT
3333	!! @tex ($kg m^{-2}$) @endtex
3334	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapnu !! Bare soil evaporation @tex ($kg m^{-2}$) @endtex
3335	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapflo !! Floodplains evaporation
3336	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapsno !! Snow evaporation @tex ($kg m^{-2}$) @endtex
3337	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: floodout !! Outflow from floodplains
3338	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: run_off_tot !! Diagnosed surface runoff @tex ($kg m^{-2}$) @endtex
3339	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: drainage !! Diagnosed rainage at the bottom of soil
3340	!! @tex ($kg m^{-2}$) @endtex
3341
3342	!! 0.2 Output variables
3343
3344	!! 0.3 Modified variables
3345
3346	!! 0.4 Local variables
3347
3348	INTEGER(i_std) :: ji, jv, jn !! Grid-cell, PFT and "nobio" fraction indices
3349	!! (unitless)
3350	REAL(r_std) :: allowed_err !! Allowed error in the water budget closure
3351	!! @tex ($kg m^{-2}$) @endtex
3352	REAL(r_std),DIMENSION (kjpindex) :: watveg !! Total amount of intercepted water in a grid-cell
3353	!! @tex ($kg m^{-2}$) @endtex
3354	REAL(r_std),DIMENSION (kjpindex) :: delta_water !! Change in total water amount
3355	!! @tex ($kg m^{-2}$) @endtex
3356	REAL(r_std),DIMENSION (kjpindex) :: tot_flux !! Algebraic sum of the water fluxes integrated over
3357	!! the timestep @tex ($kg m^{-2}$) @endtex
3358	REAL(r_std),DIMENSION (kjpindex) :: sum_snow_nobio !! Total amount of snow from the "nobio" fraction
3359	!! in a grid-cell @tex ($kg m^{-2}$) @endtex
3360	REAL(r_std),DIMENSION (kjpindex) :: sum_vevapwet !! Sum of interception loss in the grid-cell
3361	!! @tex ($kg m^{-2}$) @endtex
3362	REAL(r_std),DIMENSION (kjpindex) :: sum_transpir !! Sum of transpiration in the grid-cell
3363	!! @tex ($kg m^{-2}$) @endtex
3364	!_ ================================================================================================================================
3365
3366	!! 1. We check the water balance : this is done at the end of a time step
3367	tot_water_end(:) = zero
3368
3369	!! 1.1 If the "nobio" fraction does not complement the "bio" fraction, we issue a warning
3370	!! If the "nobio" fraction (ice, lakes, etc.) does not complement the "bio" fraction (vegetation and bare soil),
3371	!! we do not need to go any further, and we output a warning ! *** how are oceans treated ?
3372	DO ji = 1, kjpindex
3373
3374	! Modif Nathalie
3375	! IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. EPSILON(un) ) THEN
3376	IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. (100*EPSILON(un)) ) THEN
3377	WRITE(numout,*) 'HYDROL problem in vegetation or frac_nobio on point ', ji
3378	WRITE(numout,*) 'totfrac_nobio : ', totfrac_nobio(ji)
3379	WRITE(numout,*) 'vegetation fraction : ', vegtot(ji)
3380	!STOP 'in hydrolc_waterbal'
3381	ENDIF
3382	ENDDO
3383
3384	!! 1.2 We calculate the total amount of water in grid-cells at the end the time step
3385	watveg(:) = zero
3386	sum_vevapwet(:) = zero
3387	sum_transpir(:) = zero
3388	sum_snow_nobio(:) = zero
3389
3390	!cdir NODEP
3391	DO jv = 1,nvm
3392	watveg(:) = watveg(:) + qsintveg(:,jv)
3393	sum_vevapwet(:) = sum_vevapwet(:) + vevapwet(:,jv)
3394	sum_transpir(:) = sum_transpir(:) + transpir(:,jv)
3395	ENDDO
3396
3397	!cdir NODEP
3398	DO jn = 1,nnobio
3399	sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn)
3400	ENDDO
3401
3402	!cdir NODEP
3403	DO ji = 1, kjpindex
3404	tot_water_end(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + &
3405	& watveg(ji) + snow(ji) + sum_snow_nobio(ji)
3406	ENDDO
3407
3408	!! 2.3 Calculate the change in total water amount during the time step
3409	!! Calculate the change in total water amount during the time, stepand the algebraic sum
3410	!! of the water fluxes supposed to cause the total water amount change.
3411	!! If the model conserves water, they should be equal, since the fluxes are used in their integrated form,
3412	!! over the time step dt_sechiba, with a unit of kg m^{-2}.
3413	DO ji = 1, kjpindex
3414
3415	delta_water(ji) = tot_water_end(ji) - tot_water_beg(ji)
3416
3417	tot_flux(ji) = precip_rain(ji) + precip_snow(ji) + returnflow(ji) + reinfiltration(ji) + irrigation(ji) - &
3418	& sum_vevapwet(ji) - sum_transpir(ji) - vevapnu(ji) - vevapsno(ji) - vevapflo(ji) + &
3419	& floodout(ji)- run_off_tot(ji) - drainage(ji)
3420
3421	ENDDO
3422
3423	!! 1.4 We do the check, given some minimum required precision
3424	!! This is a wild guess and corresponds to what works on an IEEE machine under double precision (REAL*8).
3425	!! If the water balance is not closed at this precision, we issue a warning with some diagnostics.
3426	allowed_err = 50000EPSILON(un) ! ** how much is epsilon(1) ? where is it defined ?
3427
3428	DO ji = 1, kjpindex
3429	IF ( ABS(delta_water(ji)-tot_flux(ji)) .GT. allowed_err ) THEN
3430	WRITE(numout,*) 'HYDROL does not conserve water. The erroneous point is : ', ji
3431	WRITE(numout,) 'The error in mm/d is :', (delta_water(ji)-tot_flux(ji))/dt_sechibaone_day, &
3432	& ' and in mm/dt : ', delta_water(ji)-tot_flux(ji)
3433	WRITE(numout,*) 'delta_water : ', delta_water(ji), ' tot_flux : ', tot_flux(ji)
3434	WRITE(numout,*) 'Actual and allowed error : ', ABS(delta_water(ji)-tot_flux(ji)), allowed_err
3435	WRITE(numout,*) 'vegtot : ', vegtot(ji)
3436	WRITE(numout,*) 'precip_rain : ', precip_rain(ji)
3437	WRITE(numout,*) 'precip_snow : ', precip_snow(ji)
3438	WRITE(numout,*) 'Water from irrigation, floodplains:', reinfiltration(ji), returnflow(ji), irrigation(ji)
3439	WRITE(numout,*) 'Total water in soil :', mean_bqsb(ji) + mean_gqsb(ji)
3440	WRITE(numout,*) 'Water on vegetation :', watveg(ji)
3441	WRITE(numout,*) 'Snow mass :', snow(ji)
3442	WRITE(numout,*) 'Snow mass on ice :', sum_snow_nobio(ji)
3443	WRITE(numout,*) 'Melt water :', tot_melt(ji)
3444	WRITE(numout,*) 'evapwet : ', vevapwet(ji,:)
3445	WRITE(numout,*) 'transpir : ', transpir(ji,:)
3446	WRITE(numout,*) 'evapnu, evapsno, evapflo : ', vevapnu(ji), vevapsno(ji), vevapflo(ji)
3447	WRITE(numout,*) 'floodout : ', floodout(ji)
3448	WRITE(numout,*) 'drainage : ', drainage(ji)
3449	STOP 'in hydrolc_waterbal'
3450	ENDIF
3451
3452	ENDDO
3453
3454	!! 2. Transfer the total water amount at the end of the current timestep to the begining of the next one
3455
3456	tot_water_beg = tot_water_end
3457
3458	END SUBROUTINE hydrolc_waterbal
3459
3460
3461	!! ================================================================================================================================
3462	!! SUBROUTINE : hydrolc_alma_init
3463	!!
3464	!>\BRIEF Initialize variables needed in hydrolc_alma
3465	!!
3466	!! DESCRIPTION : None
3467	!!
3468	!! RECENT CHANGE(S) : None
3469	!!
3470	!! REFERENCE(S) : None
3471	!!
3472	!! FLOWCHART : None
3473	!! \n
3474	!_ ================================================================================================================================
3475
3476	SUBROUTINE hydrolc_alma_init (kjpindex, index, qsintveg, snow, snow_nobio)
3477
3478	!! 0. Variable and parameter declaration
3479
3480	!! 0.1 Input variables
3481
3482	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
3483	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the points on the map (unitless)
3484	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Amount of water in the canopy interception reservoir
3485	!! @tex ($kg m^{-2}$) @endtex
3486	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
3487	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in):: snow_nobio !! Snow water equivalent on nobio areas
3488	!! @tex ($kg m^{-2}$) @endtex
3489
3490	!! 0.4 Local variabless
3491
3492	INTEGER(i_std) :: ji !! Grid-cell indices (unitless)
3493	REAL(r_std) :: watveg !! Total amount of intercepted water in a grid-cell
3494	!! @tex ($kg m^{-2}$) @endtex
3495	!_ ================================================================================================================================
3496
3497	!! 1. Initialize water in terrestrial grid-cells at the beginning of the first time step
3498
3499	!! Initialize the amounts of water in terrestrial grid-cells at the beginning of the first time step,
3500	!! for the three water reservoirs (interception, soil and snow).
3501	tot_watveg_beg(:) = zero
3502	tot_watsoil_beg(:) = zero
3503	snow_beg(:) = zero
3504
3505	DO ji = 1, kjpindex
3506	watveg = SUM(qsintveg(ji,:))
3507	tot_watveg_beg(ji) = watveg
3508	tot_watsoil_beg(ji) = mean_bqsb(ji) + mean_gqsb(ji)
3509	snow_beg(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
3510	ENDDO
3511
3512	tot_watveg_end(:) = tot_watveg_beg(:)
3513	tot_watsoil_end(:) = tot_watsoil_beg(:)
3514	snow_end(:) = snow_beg(:)
3515
3516	END SUBROUTINE hydrolc_alma_init
3517
3518	!! ================================================================================================================================
3519	!! SUBROUTINE : hydrolc_alma
3520	!!
3521	!>\BRIEF This routine computes diagnostic variables required under the ALMA standards:
3522	!! changes in water amounts over the time steps (in the soil, snow and interception reservoirs); total water amount
3523	!! at the end of the time steps; soil relative humidity.
3524	!!
3525	!! DESCRIPTION : None
3526	!!
3527	!! RECENT CHANGE(S) : None
3528	!!
3529	!! MAIN OUTPUT VARIABLE(S) :
3530	!! - soilwet: mean relative humidity of soil in the grid-cells (0-1, unitless)
3531	!! - delsoilmoist, delswe, delintercept: changes in water amount during the time step @tex ($kg m^{-2}$) @endtex,
3532	!! for the three water reservoirs (soil,snow,interception)
3533	!! - tot_watsoil_end: total water amount at the end of the time steps aincluding the three water resrevoirs
3534	!! @tex ($kg m^{-2}$) @endtex.
3535	!!
3536	!! REFERENCE(S) : None
3537	!!
3538	!! FLOWCHART : None
3539	!! \n
3540	!_ ================================================================================================================================
3541
3542	SUBROUTINE hydrolc_alma (kjpindex, index, qsintveg, snow, snow_nobio, soilwet)
3543
3544	!! 0. Variable and parameter declaration
3545
3546	!! 0.1 Input variables
3547
3548	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
3549	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indices of the points on the map (unitless)
3550	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Amount of water in the canopy interception reservoir
3551	!! @tex ($kg m^{-2}$) @endtex
3552	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent @tex ($kg m^{-2}$) @endtex
3553	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in):: snow_nobio !! Snow water equivalent on nobio areas
3554	!! @tex ($kg m^{-2}$) @endtex
3555
3556	!! 0.2 Output variables
3557
3558	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilwet !! Mean relative humidity of soil in the grid-cells
3559	!! (0-1, unitless)
3560
3561	!! 0.3 Modified variables
3562
3563	!! 0.4 Local variabless
3564
3565	INTEGER(i_std) :: ji !! Grid-cell indices (unitless)
3566	REAL(r_std) :: watveg !! Total amount of intercepted water in a grid-cell
3567	!! @tex ($kg m^{-2}$) @endtex
3568	!_ ================================================================================================================================
3569
3570	!! 1. We calculate the required variables at the end of each time step
3571
3572	tot_watveg_end(:) = zero
3573	tot_watsoil_end(:) = zero
3574	snow_end(:) = zero
3575	delintercept(:) = zero
3576	delsoilmoist(:) = zero
3577	delswe(:) = zero
3578
3579	DO ji = 1, kjpindex
3580	watveg = SUM(qsintveg(ji,:))
3581	tot_watveg_end(ji) = watveg
3582	tot_watsoil_end(ji) = mean_bqsb(ji) + mean_gqsb(ji)
3583	snow_end(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
3584	!
3585	delintercept(ji) = tot_watveg_end(ji) - tot_watveg_beg(ji)
3586	delsoilmoist(ji) = tot_watsoil_end(ji) - tot_watsoil_beg(ji)
3587	delswe(ji) = snow_end(ji) - snow_beg(ji)
3588	!
3589	ENDDO
3590
3591	!! 2. Transfer the total water amount at the end of the current timestep to the begining of the next one.
3592
3593	tot_watveg_beg = tot_watveg_end
3594	tot_watsoil_beg = tot_watsoil_end
3595	snow_beg(:) = snow_end(:)
3596
3597	!! 3. We calculate the mean relative humidity of soil in the grid-cells
3598
3599	DO ji = 1,kjpindex
3600	soilwet(ji) = tot_watsoil_end(ji) / mx_eau_var(ji)
3601	ENDDO
3602
3603	END SUBROUTINE hydrolc_alma
3604
3605
3606	!! ================================================================================================================================
3607	!! SUBROUTINE : hydrolc_hdiff
3608	!!
3609	!>\BRIEF This subroutine performs horizontal diffusion of water between each PFT/soil column,
3610	!! if the flag ok_hdiff is true.
3611	!!
3612	!! DESCRIPTION : The diffusion is realized on the water contents of both the top and bottom
3613	!! soil layers, but the bottom soil layers share the same soil moisture since the end of hydrolc_soil (step 6).
3614	!! This subroutine thus only modifies the moisture of the top soil layer.
3615	!! \latexonly
3616	!! \input{hdiff.tex}
3617	!! \endlatexonly
3618	!!
3619	!! RECENT CHANGE(S) : None
3620	!!
3621	!! MAIN OUTPUT VARIABLE(S) : gqsb, bqsb, dsg, dss, dsp
3622	!!
3623	!! REFERENCE(S) : None
3624	!!
3625	!! FLOWCHART : None
3626	!! \n
3627	!_ ================================================================================================================================
3628
3629	SUBROUTINE hydrolc_hdiff(kjpindex, veget, tot_bare_soil, ruu_ch, gqsb, bqsb, dsg, dss, dsp)
3630
3631	!! 0. Variable and parameter declaration
3632
3633	!! 0.1 Input variables
3634
3635	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size (number of grid cells) (unitless)
3636	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Grid-cell fraction effectively covered by
3637	!! vegetation for each PFT, except for soil
3638	!! (0-1, unitless)
3639	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_bare_soil !! Total evaporating bare soil fraction
3640	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ruu_ch !! Volumetric soil water capacity
3641	!! @tex ($kg m^{-3}$) @endtex
3642
3643	!! 0.2 Output variables
3644
3645	!! 0.3 Modified variables
3646
3647	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Water content in the top layer
3648	!! @tex ($kg m^{-2}$) @endtex
3649	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Water content in the bottom layer
3650	!! @tex ($kg m^{-2}$) @endtex
3651	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Depth of the top layer (m)
3652	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Depth of dry soil at the top, whether in the top
3653	!! or bottom layer (m)
3654	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Depth of dry soil in the bottom layer plus depth
3655	!! of top layer (m)
3656
3657	!! 0.4 Local variables
3658
3659	REAL(r_std), DIMENSION (kjpindex) :: bqsb_mean !! Mean value of bqsb in each grid-cell
3660	!! @tex ($kg m^{-2}$) @endtex
3661	REAL(r_std), DIMENSION (kjpindex) :: gqsb_mean !! Mean value of gqsb in each grid-cell
3662	!! @tex ($kg m^{-2}$) @endtex
3663	REAL(r_std), DIMENSION (kjpindex) :: dss_mean !! Mean value of dss in each grid-cell (m)
3664	REAL(r_std), DIMENSION (kjpindex) :: vegtot !! Total fraction of grid-cell covered by PFTs (bare
3665	!! soil + vegetation) (0-1, unitless)
3666	! *** this variable is declared at the module level
3667	REAL(r_std) :: x !! Coefficient of diffusion by time step
3668	!! (0-1, unitless)
3669	INTEGER(i_std) :: ji,jv !! Indices for grid-cells and PFTs (unitless)
3670	REAL(r_std), SAVE :: tau_hdiff !! Time scale for horizontal water diffusion (s)
3671	!$OMP THREADPRIVATE(tau_hdiff)
3672	LOGICAL, SAVE :: firstcall_hydrolc=.TRUE. !! Flag to set tau_hdiff at the first time step
3673	!! ! *** why not define this parameter in
3674	!! constantes_soil.f90 ?
3675	!$OMP THREADPRIVATE(firstcall_hydrolc)
3676	!_ ================================================================================================================================
3677
3678	!! 1. First time step: we assign a value to the time scale for horizontal water diffusion (in seconds)
3679
3680	IF ( firstcall_hydrolc ) THEN
3681
3682	!Config Key = HYDROL_TAU_HDIFF
3683	!Config Desc = time scale (s) for horizontal diffusion of water
3684	!Config Def = one_day
3685	!Config If = HYDROL_OK_HDIFF
3686	!Config Help = Defines how fast diffusion occurs horizontally between
3687	!Config the individual PFTs' water reservoirs. If infinite, no
3688	!Config diffusion.
3689	!Config Units = [seconds]
3690	tau_hdiff = one_day ! 86400 s
3691	CALL getin_p('HYDROL_TAU_HDIFF',tau_hdiff)
3692
3693	WRITE (numout,*) 'Hydrol: Horizontal diffusion, tau (s)=',tau_hdiff
3694
3695	firstcall_hydrolc = .FALSE.
3696
3697	ENDIF
3698
3699	!! 2. We calculate mean values of bqsb, gqsb and dss over the grid-cell ("bio" fraction).
3700
3701	! Calculate mean values of bqsb, gqsb and dss over the grid-cell ("bio" fraction)
3702	! This could be done with SUM instruction but this kills vectorization
3703	! * We compute here vegtot=SUM(veget), but it is already defined at a higher level (declared in the module) *
3704	bqsb_mean(:) = zero
3705	gqsb_mean(:) = zero
3706	dss_mean(:) = zero
3707	vegtot(:) = zero
3708	!
3709	DO jv = 1, nvm
3710	DO ji = 1, kjpindex
3711	!MM veget(:,1) BUG ??!!!!!!!!!!!
3712	IF (jv .EQ. 1) THEN
3713	bqsb_mean(ji) = bqsb_mean(ji) + tot_bare_soil(ji)*bqsb(ji,jv)
3714	gqsb_mean(ji) = gqsb_mean(ji) + tot_bare_soil(ji)*gqsb(ji,jv)
3715	dss_mean(ji) = dss_mean(ji) + tot_bare_soil(ji)*dss(ji,jv)
3716	vegtot(ji) = vegtot(ji) + tot_bare_soil(ji)
3717	ELSE
3718	bqsb_mean(ji) = bqsb_mean(ji) + veget(ji,jv)*bqsb(ji,jv)
3719	gqsb_mean(ji) = gqsb_mean(ji) + veget(ji,jv)*gqsb(ji,jv)
3720	dss_mean(ji) = dss_mean(ji) + veget(ji,jv)*dss(ji,jv)
3721	vegtot(ji) = vegtot(ji) + veget(ji,jv)
3722	ENDIF
3723	ENDDO
3724	ENDDO
3725
3726	DO ji = 1, kjpindex
3727	IF (vegtot(ji) .GT. zero) THEN
3728	bqsb_mean(ji) = bqsb_mean(ji)/vegtot(ji)
3729	gqsb_mean(ji) = gqsb_mean(ji)/vegtot(ji)
3730	dss_mean(ji) = dss_mean(ji)/vegtot(ji)
3731	ENDIF
3732	ENDDO
3733
3734	!! 3. Relax PFT values towards grid-cell mean
3735
3736	!! Relax values towards mean : the diffusion is proportional to the deviation to the mean,
3737	!! and inversely proportional to the timescale tau_hdiff
3738	!! We change the moisture of two soil layers, but also the depth of dry in the top soil layer.
3739	!! Therefore, we need to accordingly adjust the top soil layer depth
3740	!! and the dry soil depth above the bottom moisture (dsp).
3741	! *** This sequence does not seem to be fully consistent with the variable definitions,
3742	! *** especially for the dry soil depths dss and dsp:
3743	! *** (i) is the "diffusion" of dss correct is some PFTs only have one soil layer (dss=dsp),
3744	! *** given that hydrolc_soil acted to merge the bottom soil moistures bqsb
3745	! *** (ii) if gqsb is very small, we let the top layer vanish, but dss is not updated to dsp
3746	!
3747	! *** Also, it would be make much more sense to perform the diffusion in hydrolc_soil, when the
3748	! *** bottom soil moistures are homogenized. It would benefit from the iterative consistency
3749	! *** check between the dsg and dsp, and it would prevent from doing some calculations twice.
3750	! *** Finally, it would be better to keep the water balance calculation for the real end of the
3751	! *** hydrological processes.
3752	!
3753	! *** FORTUNATELY, the default value of ok_hdiff is false (hydrolc_init.f90)
3754	x = MAX( zero, MIN( dt_sechiba/tau_hdiff, un ) )
3755
3756	DO jv = 1, nvm
3757	DO ji = 1, kjpindex
3758
3759	! *** bqsb does not change as bqsb(ji,jv)=bqsb_mean(ji) since hydrolc_soil, step 6
3760	bqsb(ji,jv) = (un-x) * bqsb(ji,jv) + x * bqsb_mean(ji)
3761	gqsb(ji,jv) = (un-x) * gqsb(ji,jv) + x * gqsb_mean(ji)
3762
3763	! *** is it meaningful to apply the diffusion equation to dss ?
3764	dss(ji,jv) = (un-x) * dss(ji,jv) + x * dss_mean(ji)
3765
3766	IF (gqsb(ji,jv) .LT. min_sechiba) THEN
3767	dsg(ji,jv) = zero ! *** in this case, we should also set dss to dsp (defined below)
3768	ELSE
3769	dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) / ruu_ch(ji)
3770	ENDIF
3771	dsp(ji,jv) = zmaxh - bqsb(ji,jv) / ruu_ch(ji) ! no change here since bqsb does not change
3772
3773	ENDDO
3774	ENDDO
3775
3776	END SUBROUTINE hydrolc_hdiff
3777
3778	END MODULE hydrolc

Note: See TracBrowser for help on using the repository browser.

Download in other formats: