Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

hydrolc.f90 @ 7746

Last change on this file since 7746 was 3965, checked in by jan.polcher, 8 years ago

Merge with trunk at revision3959.
This includes all the developments made for CMIP6 and passage to XIOS2.
All conflicts are resolved and the code compiles.

But it still does not link because of an "undefined reference to `_intel_fast_memmove'"

Property svn:keywords set to Revision Date HeadURL Date Author Revision

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_sechiba/hydrolc.f90 @ 7746

Download in other formats: