Context Navigation

hydrolc.f90 @ 7326

Last change on this file since 7326 was 2608, checked in by sebastiaan.luyssaert, 9 years ago
DEV: code compiles but needs to be tested. Trunk changes up to and including r2223 thus MICT with the bug. Note that the explicitsnow code was copied but has some conflicts with the DOFOCO code
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 193.3 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_sechiba/hydrolc.f90 @ 7326

Download in other formats: