Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_sechiba/hydrolc.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

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

File size: 204.5 KB

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

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

Download in other formats: