Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

hydrolc.f90 @ 8787

Last change on this file since 8787 was 5080, checked in by chunjing.qiu, 7 years ago
soil freezing, soil moisture, fwet bugs fixed
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 204.4 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-PEAT_r5488/src_sechiba/hydrolc.f90 @ 8787

Download in other formats: