Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

hydrolc.f90 @ 7325

Last change on this file since 7325 was 7297, checked in by wei.li, 3 years ago
updated code for publication, 2021,9,25
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 204.7 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-BIOENERGY_r7298/src_sechiba/hydrolc.f90 @ 7325

Download in other formats: