Context Navigation

source: branches/publications/ORCHIDEE_v2_r5968/src_sechiba/hydrolc.f90

Last change on this file was 5091, checked in by josefine.ghattas, 7 years ago
Change unit for output variable snowliq and snowliqtot from m to kg/m2 and take into acount contfrac.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 197.6 KB

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

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

Download in other formats: