Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/hydrolc.f90 @ 116

Last change on this file since 116 was 112, checked in by didier.solyga, 13 years ago
uncomment the division of throughfall_by_pft, could be wrong if division made after getin in the previous commit
File size: 105.0 KB

Line
1	!!
2	!! This module computes hydrologic processes on continental points.
3	!!
4	!! @author Marie-Alice Foujols and Jan Polcher
5	!! @Version : $Revision: 1.21 $, $Date: 2010/05/07 08:28:13 $
6	!!
7	!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/hydrolc.f90,v 1.21 2010/05/07 08:28:13 ssipsl Exp $
8	!! IPSL (2006)
9	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!!
11	MODULE hydrolc
12	!
13	!
14	! routines called : restput, restget
15	!
16	USE ioipsl
17	!
18	! modules used :
19	USE constantes
20	USE pft_parameters
21	USE sechiba_io
22	USE grid
23	USE parallel
24	! USE Write_Field_p
25
26	IMPLICIT NONE
27
28	! public routines :
29	! hydrol
30	PRIVATE
31	PUBLIC :: hydrolc_main,hydrolc_clear
32
33	!
34	! variables used inside hydrol module : declaration and initialisation
35	!
36	LOGICAL, SAVE :: l_first_hydrol=.TRUE. !! Initialisation has to be done one time
37	!
38	LOGICAL, SAVE :: check_waterbal=.FALSE. !! The check the water balance
39	LOGICAL, SAVE :: ok_hdiff !! do horizontal diffusion?
40
41	CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O
42
43	! one dimension array allocated, computed, saved and got in hydrol module
44
45	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: bqsb !! Hauteur d'eau dans le reservoir profond
46	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gqsb !! Hauteur d'eau dans le reservoir de surface
47	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsg !! Hauteur du reservoir de surface
48	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsp !! Hauteur au dessus du reservoir profond
49	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_bqsb !! diagnostique du reservoir profond
50	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_gqsb !! diagnostique du reservoir de surface
51	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_beg !! Total amount of water at start of time step
52	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_end !! Total amount of water at end of time step
53	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_beg !! Total amount of water on vegetation at start of time step
54	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_end !! Total amount of water on vegetation at end of time step
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_beg !! Total amount of water in the soil at start of time step
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_end !! Total amount of water in the soil at end of time step
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_beg !! Total amount of snow at start of time step
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_end !! Total amount of snow at end of time step
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delsoilmoist !! Change in soil moisture
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delintercept !! Change in interception storage
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delswe !! Change in SWE
62
63	! one dimension array allocated, computed and used in hydrol module exclusively
64
65	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dss !! Hauteur au dessus du reservoir de surface
66	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: hdry !! Dry soil heigth in meters
67	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Eau tombee sur le sol
68	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsnowveg !! Sublimation of snow on vegetation
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snowmelt !! Quantite de neige fondue
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: icemelt !! Quantite de glace fondue
72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gdrainage !! Drainage between reservoirs
73
74
75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vegtot !! Total vegetation
76	! The last vegetation map which was used to distribute the reservoirs
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: resdist !! Distribution of reservoirs
78
79	!! profondeur du reservoir contenant le maximum d'eau
80	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ruu_ch !! Quantite d'eau maximum
82	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: runoff !! Ruissellement
83
84	! Ajout Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
85	! Modifs stabilite
86	REAL(r_std), PARAMETER :: dsg_min = 0.001
87
88	CONTAINS
89
90	!!
91	!! Main routine for hydrol module
92	!! - called only one time for initialisation
93	!! - called every time step
94	!! - called one more time at last time step for writing _restart_ file
95	!!
96	!! Algorithm:
97	!! - call hydrolc_snow for snow process (including age of snow)
98	!! - call hydrolc_canop for canopy process
99	!! - call hydrolc_soil for bare soil process
100	!!
101	!! @call hydrolc_snow
102	!! @call hydrolc_canop
103	!! @call hydrolc_soil
104	!!
105	SUBROUTINE hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
106	& temp_sol_new, run_off_tot, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
107	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, &
108	& precip_rain, precip_snow, returnflow, irrigation, humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr,&
109	& shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id)
110
111	! interface description
112	! input scalar
113	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
114	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
115	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
116	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
117	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
118	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
119	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
120	! input fields
121	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
122	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map
123	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
124	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
125	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil
126	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation returning to soil moisture
127	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
128	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, ...
129	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+...
130	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
131	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss
132	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
133	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty)
134	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
135	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
136	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation
137	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction
138	! modified fields
139	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapnu !! Bare soil evaporation
140	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapsno !! Snow evaporation
141	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow !! Snow mass [Kg/m^2]
142	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow_age !! Snow age
143	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2]
144	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
145	!! We consider that any water on the ice is snow and we only peforme a water balance to have consistency.
146	!! The water balance is limite to + or - 10^6 so that accumulation is not endless
147	! output fields
148	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: run_off_tot !! Complete runoff
149	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drainage !! Drainage
150	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Relative humidity
151	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
152	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistence to bare soil evaporation
153	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visibly dry soil (between 0 and 1)
154	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out):: shumdiag !! relative soil moisture
155	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
156	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt
157	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception
158
159	!
160	! local declaration
161	!
162	REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness
163	REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Depth of snow layer
164
165	INTEGER(i_std) :: ji,jv
166	REAL(r_std), DIMENSION(kjpindex) :: histvar !! computations for history files
167
168	!
169	! do initialisation
170	!
171	IF (l_first_hydrol) THEN
172
173	sneige = snowcri/mille
174
175	IF (long_print) WRITE (numout,*) ' l_first_hydrol : call hydrolc_init '
176
177	CALL hydrolc_init (kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, &
178	& snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
179
180	CALL hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag)
181	!
182	! If we check the water balance we first save the total amount of water
183	!
184	IF (check_waterbal) THEN
185	CALL hydrolc_waterbal(kjpindex, index, .TRUE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
186	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,&
187	& run_off_tot, drainage)
188	ENDIF
189	!
190	IF (almaoutput) THEN
191	CALL hydrolc_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwet)
192	ENDIF
193
194	!
195	! shared time step
196	!
197	IF (long_print) WRITE (numout,*) 'hydrolc pas de temps = ',dtradia
198
199	RETURN
200
201	ENDIF
202
203	!
204	! prepares restart file for the next simulation
205	!
206	IF (ldrestart_write) THEN
207
208	IF (long_print) WRITE (numout,*) ' we have to complete restart file with HYDROLOGIC variables '
209
210	var_name= 'humrel'
211	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, humrel, 'scatter', nbp_glo, index_g)
212	!
213	var_name= 'vegstress'
214	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, vegstress, 'scatter', nbp_glo, index_g)
215	!
216	var_name= 'snow'
217	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow, 'scatter', nbp_glo, index_g)
218	!
219	var_name= 'snow_age'
220	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow_age, 'scatter', nbp_glo, index_g)
221	!
222	var_name= 'snow_nobio'
223	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio, 'scatter', nbp_glo, index_g)
224	!
225	var_name= 'snow_nobio_age'
226	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio_age, 'scatter', nbp_glo, index_g)
227	!
228	var_name= 'bqsb'
229	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, bqsb, 'scatter', nbp_glo, index_g)
230	!
231	var_name= 'gqsb'
232	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, gqsb, 'scatter', nbp_glo, index_g)
233	!
234	var_name= 'dsg'
235	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsg, 'scatter', nbp_glo, index_g)
236	!
237	var_name= 'dsp'
238	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsp, 'scatter', nbp_glo, index_g)
239	!
240	var_name= 'dss'
241	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dss, 'scatter', nbp_glo, index_g)
242	!
243	var_name= 'hdry'
244	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, hdry, 'scatter', nbp_glo, index_g)
245	!
246	var_name= 'qsintveg'
247	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, qsintveg, 'scatter', nbp_glo, index_g)
248	!
249	var_name= 'resdist'
250	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, resdist, 'scatter', nbp_glo, index_g)
251	RETURN
252	!
253	END IF
254
255	!
256	! computes snow
257	!
258
259	CALL hydrolc_snow(kjpindex, dtradia, precip_rain, precip_snow, temp_sol_new, soilcap, &
260	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
261	& tot_melt, snowdepth)
262	!
263	! computes canopy
264	!
265	!
266	CALL hydrolc_vegupd(kjpindex, veget, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss,dsp, resdist)
267	!
268	CALL hydrolc_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol)
269	!
270	! computes hydro_soil
271	!
272	CALL hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, irrigation, tot_melt, mx_eau_var, veget, ruu_ch, transpir,&
273	& gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, runoff, run_off_tot, drainage, humrel, vegstress, &
274	& shumdiag, litterhumdiag)
275	!
276	! computes horizontal diffusion between the water reservoirs
277	!
278	IF ( ok_hdiff ) THEN
279	CALL hydrolc_hdiff(kjpindex, dtradia, veget, ruu_ch, gqsb, bqsb, dsg, dss, dsp)
280	ENDIF
281	!
282	! If we check the water balance we end with the comparison of total water change and fluxes
283	!
284	IF (check_waterbal) THEN
285	CALL hydrolc_waterbal(kjpindex, index, .FALSE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
286	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,&
287	& run_off_tot, drainage )
288	ENDIF
289	!
290	! If we use the ALMA standards
291	!
292	IF (almaoutput) THEN
293	CALL hydrolc_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet)
294	ENDIF
295
296
297	!
298	IF ( .NOT. almaoutput ) THEN
299	CALL histwrite(hist_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
300	CALL histwrite(hist_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
301	CALL histwrite(hist_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
302	CALL histwrite(hist_id, 'runoff', kjit, run_off_tot, kjpindex, index)
303	CALL histwrite(hist_id, 'drainage', kjit, drainage, kjpindex, index)
304	CALL histwrite(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
305	CALL histwrite(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
306	CALL histwrite(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
307	CALL histwrite(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
308	CALL histwrite(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
309	CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
310
311	histvar(:)=zero
312	DO jv = 1, nvm
313	DO ji = 1, kjpindex
314	IF ( vegtot(ji) .GT. zero ) THEN
315	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))mx_eau_eau, 0.0)
316	ENDIF
317	ENDDO
318	ENDDO
319	CALL histwrite(hist_id, 'mrsos', kjit, histvar, kjpindex, index)
320
321	histvar(:)=mean_bqsb(:)+mean_gqsb(:)
322	CALL histwrite(hist_id, 'mrso', kjit, histvar, kjpindex, index)
323
324	histvar(:)=run_off_tot(:)/86400.
325	CALL histwrite(hist_id, 'mrros', kjit, histvar, kjpindex, index)
326
327	histvar(:)=(run_off_tot(:)+drainage(:))/86400.
328	CALL histwrite(hist_id, 'mrro', kjit, histvar, kjpindex, index)
329
330	histvar(:)=(precip_rain(:)-SUM(precisol(:,:),dim=2))/86400.
331	CALL histwrite(hist_id, 'prveg', kjit, histvar, kjpindex, index)
332
333	ELSE
334	CALL histwrite(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
335	CALL histwrite(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
336	CALL histwrite(hist_id, 'Qs', kjit, run_off_tot, kjpindex, index)
337	CALL histwrite(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
338	CALL histwrite(hist_id, 'Qsm', kjit, tot_melt, kjpindex, index)
339	CALL histwrite(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
340	CALL histwrite(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
341	CALL histwrite(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
342	!
343	CALL histwrite(hist_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
344	CALL histwrite(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
345	CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
346	!
347	CALL histwrite(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
348	CALL histwrite(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
349	!
350	CALL histwrite(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
351	!
352	ENDIF
353	IF ( hist2_id > 0 ) THEN
354	IF ( .NOT. almaoutput ) THEN
355	CALL histwrite(hist2_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
356	CALL histwrite(hist2_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
357	CALL histwrite(hist2_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
358	CALL histwrite(hist2_id, 'runoff', kjit, run_off_tot, kjpindex, index)
359	CALL histwrite(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
360	CALL histwrite(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
361	CALL histwrite(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
362	CALL histwrite(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
363	CALL histwrite(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
364	CALL histwrite(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
365	CALL histwrite(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
366
367	histvar(:)=zero
368	DO jv = 1, nvm
369	DO ji = 1, kjpindex
370	IF ( vegtot(ji) .GT. zero ) THEN
371	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))mx_eau_eau, 0.0)
372	ENDIF
373	ENDDO
374	ENDDO
375	CALL histwrite(hist2_id, 'mrsos', kjit, histvar, kjpindex, index)
376
377	histvar(:)=(run_off_tot(:)+drainage(:))/86400.
378	CALL histwrite(hist2_id, 'mrro', kjit, histvar, kjpindex, index)
379
380	ELSE
381	CALL histwrite(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
382	CALL histwrite(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
383	CALL histwrite(hist2_id, 'Qs', kjit, run_off_tot, kjpindex, index)
384	CALL histwrite(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
385	CALL histwrite(hist2_id, 'Qsm', kjit, tot_melt, kjpindex, index)
386	CALL histwrite(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
387	CALL histwrite(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
388	CALL histwrite(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
389	!
390	CALL histwrite(hist2_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
391	CALL histwrite(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
392	!
393	CALL histwrite(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
394	CALL histwrite(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
395	!
396	CALL histwrite(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
397	!
398	ENDIF
399	ENDIF
400
401	!
402	IF (long_print) WRITE (numout,*) ' hydrolc_main Done '
403
404	END SUBROUTINE hydrolc_main
405
406	!! Algorithm:
407	!! - dynamic allocation for local array
408	!! - _restart_ file reading for HYDROLOGIC variables
409	!!
410	SUBROUTINE hydrolc_init(kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, &
411	& snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
412
413	! interface description
414	! input scalar
415	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
416	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read
417	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
418	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
419	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
420	! input fields
421	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Carte de vegetation
422	! output fields
423	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Stress hydrique, relative humidity
424	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
425	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow mass [Kg/m^2]
426	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age
427	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio !! Snow on ice, lakes, ...
428	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age !! Snow age on ice, lakes, ...
429	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception
430
431	! local declaration
432	INTEGER(i_std) :: ier
433	INTEGER(i_std) :: ji,jv,ik
434	REAL(r_std), DIMENSION (kjpindex,nvm) :: zdsp, tmpdss
435
436	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: dsp_g
437	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: zdsp_g
438
439	REAL(r_std), DIMENSION(kjpindex) :: a_subgrd
440
441	! initialisation
442	IF (l_first_hydrol) THEN
443	l_first_hydrol=.FALSE.
444	ELSE
445	WRITE (numout,*) ' l_first_hydrol false . we stop '
446	STOP 'hydrolc_init'
447	ENDIF
448
449	!Config Key = HYDROL_OK_HDIFF
450	!Config Desc = do horizontal diffusion?
451	!Config Def = n
452	!Config Help = If TRUE, then water can diffuse horizontally between
453	!Config the PFTs' water reservoirs.
454
455	ok_hdiff = .FALSE.
456	CALL getin_p('HYDROL_OK_HDIFF',ok_hdiff)
457
458	! make dynamic allocation with good dimension
459
460	! one dimension array allocation with possible restart value
461
462	ALLOCATE (bqsb(kjpindex,nvm),stat=ier)
463	IF (ier.NE.0) THEN
464	WRITE (numout,*) ' error in bqsb allocation. We stop. We need kjpindex words = ',kjpindex
465	STOP 'hydrolc_init'
466	END IF
467	bqsb(:,:) = zero
468
469	ALLOCATE (gqsb(kjpindex,nvm),stat=ier)
470	IF (ier.NE.0) THEN
471	WRITE (numout,*) ' error in gqsb allocation. We stop. We need kjpindex words = ',kjpindex
472	STOP 'hydrolc_init'
473	END IF
474	gqsb(:,:) = zero
475
476	ALLOCATE (dsg(kjpindex,nvm),stat=ier)
477	IF (ier.NE.0) THEN
478	WRITE (numout,*) ' error in dsg allocation. We stop. We need kjpindex words = ',kjpindex
479	STOP 'hydrolc_init'
480	END IF
481	dsg(:,:) = zero
482
483	ALLOCATE (dsp(kjpindex,nvm),stat=ier)
484	IF (ier.NE.0) THEN
485	WRITE (numout,*) ' error in dsp allocation. We stop. We need kjpindex words = ',kjpindex
486	STOP 'hydrolc_init'
487	END IF
488	dsp(:,:) = zero
489
490	! one dimension array allocation
491
492	ALLOCATE (mean_bqsb(kjpindex),stat=ier)
493	IF (ier.NE.0) THEN
494	WRITE (numout,*) ' error in mean_bqsb allocation. We stop. We need kjpindex words = ',kjpindex
495	STOP 'hydrolc_init'
496	END IF
497	mean_bqsb(:) = zero
498
499	ALLOCATE (mean_gqsb(kjpindex),stat=ier)
500	IF (ier.NE.0) THEN
501	WRITE (numout,*) ' error in mean_gqsb allocation. We stop. We need kjpindex words = ',kjpindex
502	STOP 'hydrolc_init'
503	END IF
504	mean_gqsb(:) = zero
505
506	ALLOCATE (dss(kjpindex,nvm),stat=ier)
507	IF (ier.NE.0) THEN
508	WRITE (numout,*) ' error in dss allocation. We stop. We need kjpindex words = ',kjpindex
509	STOP 'hydrolc_init'
510	END IF
511	dss(:,:) = zero
512
513	ALLOCATE (hdry(kjpindex),stat=ier)
514	IF (ier.NE.0) THEN
515	WRITE (numout,*) ' error in hdry allocation. We stop. We need kjpindex words = ',kjpindex
516	STOP 'hydrolc_init'
517	END IF
518	hdry(:) = zero
519
520	ALLOCATE (precisol(kjpindex,nvm),stat=ier)
521	IF (ier.NE.0) THEN
522	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
523	STOP 'hydrolc_init'
524	END IF
525	precisol(:,:) = zero
526
527	ALLOCATE (gdrainage(kjpindex,nvm),stat=ier)
528	IF (ier.NE.0) THEN
529	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
530	STOP 'hydrolc_init'
531	END IF
532	gdrainage(:,:) = zero
533
534	ALLOCATE (subsnowveg(kjpindex),stat=ier)
535	IF (ier.NE.0) THEN
536	WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex
537	STOP 'hydrolc_init'
538	END IF
539	subsnowveg(:) = zero
540
541	ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier)
542	IF (ier.NE.0) THEN
543	WRITE (numout,) ' error in subsnownobio allocation. We stop. We need kjpindexnnobio words = ', &
544	kjpindex*nnobio
545	STOP 'hydrolc_init'
546	END IF
547	subsnownobio(:,:) = zero
548
549	ALLOCATE (snowmelt(kjpindex),stat=ier)
550	IF (ier.NE.0) THEN
551	WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex
552	STOP 'hydrolc_init'
553	END IF
554	snowmelt(:) = zero
555
556	ALLOCATE (icemelt(kjpindex),stat=ier)
557	IF (ier.NE.0) THEN
558	WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex
559	STOP 'hydrolc_init'
560	END IF
561	icemelt(:) = zero
562
563	ALLOCATE (mx_eau_var(kjpindex),stat=ier)
564	IF (ier.NE.0) THEN
565	WRITE (numout,*) ' error in mx_eau_var allocation. We stop. We need kjpindex words = ',kjpindex
566	STOP 'hydrolc_init'
567	END IF
568	mx_eau_var(:) = zero
569
570	ALLOCATE (ruu_ch(kjpindex),stat=ier)
571	IF (ier.NE.0) THEN
572	WRITE (numout,*) ' error in ruu_ch allocation. We stop. We need kjpindex words = ',kjpindex
573	STOP 'hydrolc_init'
574	END IF
575	ruu_ch(:) = zero
576
577	ALLOCATE (vegtot(kjpindex),stat=ier)
578	IF (ier.NE.0) THEN
579	WRITE (numout,) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindexnvm
580	STOP 'hydrolc_init'
581	END IF
582	vegtot(:) = zero
583
584	ALLOCATE (resdist(kjpindex,nvm),stat=ier)
585	IF (ier.NE.0) THEN
586	WRITE (numout,) ' error in resdist allocation. We stop. We need kjpindex words = ',kjpindexnvm
587	STOP 'hydrolc_init'
588	END IF
589	resdist(:,:) = zero
590
591	ALLOCATE (runoff(kjpindex,nvm),stat=ier)
592	IF (ier.NE.0) THEN
593	WRITE (numout,) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindexnvm
594	STOP 'hydrolc_init'
595	END IF
596	runoff(:,:) = zero
597	!
598	! If we check the water balance we need two more variables
599	!
600	IF ( check_waterbal ) THEN
601
602	ALLOCATE (tot_water_beg(kjpindex),stat=ier)
603	IF (ier.NE.0) THEN
604	WRITE (numout,*) ' error in tot_water_beg allocation. We stop. We need kjpindex words = ',kjpindex
605	STOP 'hydrolc_init'
606	END IF
607
608	ALLOCATE (tot_water_end(kjpindex),stat=ier)
609	IF (ier.NE.0) THEN
610	WRITE (numout,*) ' error in tot_water_end allocation. We stop. We need kjpindex words = ',kjpindex
611	STOP 'hydrolc_init'
612	END IF
613
614	ENDIF
615	!
616	! If we use the almaoutputs we need four more variables
617	!
618	IF ( almaoutput ) THEN
619
620	ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
621	IF (ier.NE.0) THEN
622	WRITE (numout,*) ' error in tot_watveg_beg allocation. We stop. We need kjpindex words = ',kjpindex
623	STOP 'hydrolc_init'
624	END IF
625
626	ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
627	IF (ier.NE.0) THEN
628	WRITE (numout,*) ' error in tot_watveg_end allocation. We stop. We need kjpindex words = ',kjpindex
629	STOP 'hydrolc_init'
630	END IF
631
632	ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
633	IF (ier.NE.0) THEN
634	WRITE (numout,*) ' error in tot_watsoil_beg allocation. We stop. We need kjpindex words = ',kjpindex
635	STOP 'hydrolc_init'
636	END IF
637
638	ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
639	IF (ier.NE.0) THEN
640	WRITE (numout,*) ' error in tot_watsoil_end allocation. We stop. We need kjpindex words = ',kjpindex
641	STOP 'hydrolc_init'
642	END IF
643
644	ALLOCATE (delsoilmoist(kjpindex),stat=ier)
645	IF (ier.NE.0) THEN
646	WRITE (numout,*) ' error in delsoilmoist allocation. We stop. We need kjpindex words = ',kjpindex
647	STOP 'hydrolc_init'
648	END IF
649
650	ALLOCATE (delintercept(kjpindex),stat=ier)
651	IF (ier.NE.0) THEN
652	WRITE (numout,*) ' error in delintercept. We stop. We need kjpindex words = ',kjpindex
653	STOP 'hydrolc_init'
654	END IF
655
656	ALLOCATE (delswe(kjpindex),stat=ier)
657	IF (ier.NE.0) THEN
658	WRITE (numout,*) ' error in delswe. We stop. We need kjpindex words = ',kjpindex
659	STOP 'hydrolc_init'
660	ENDIF
661
662	ALLOCATE (snow_beg(kjpindex),stat=ier)
663	IF (ier.NE.0) THEN
664	WRITE (numout,*) ' error in snow_beg allocation. We stop. We need kjpindex words =',kjpindex
665	STOP 'hydrolc_init'
666	END IF
667
668	ALLOCATE (snow_end(kjpindex),stat=ier)
669	IF (ier.NE.0) THEN
670	WRITE (numout,*) ' error in snow_end allocation. We stop. We need kjpindex words =',kjpindex
671	STOP 'hydrolc_init'
672	END IF
673
674	ENDIF
675
676	! open restart input file done by sechiba_init
677	! and read data from restart input file for HYDROLOGIC process
678
679	IF (ldrestart_read) THEN
680
681	IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
682
683	var_name= 'snow'
684	CALL ioconf_setatt('UNITS', 'kg/m^2')
685	CALL ioconf_setatt('LONG_NAME','Snow mass')
686	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
687	!!$ ! correction for old restart
688	!!$ DO ik=1, kjpindex
689	!!$ if(snow(ik).gt.maxmass_glacier) snow(ik)=maxmass_glacier
690	!!$ ENDDO
691	!
692	var_name= 'snow_age'
693	CALL ioconf_setatt('UNITS', 'd')
694	CALL ioconf_setatt('LONG_NAME','Snow age')
695	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
696	!
697	var_name= 'snow_nobio'
698	CALL ioconf_setatt('UNITS', 'kg/m^2')
699	CALL ioconf_setatt('LONG_NAME','Snow on other surface types')
700	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
701	!!$ ! correction for old restart
702	!!$ DO ik=1, kjpindex
703	!!$ if(snow_nobio(ik,iice).gt.maxmass_glacier) snow_nobio(ik,iice)=maxmass_glacier
704	!!$ ENDDO
705	!
706	var_name= 'snow_nobio_age'
707	CALL ioconf_setatt('UNITS', 'd')
708	CALL ioconf_setatt('LONG_NAME','Snow age on other surface types')
709	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
710	!
711	var_name= 'humrel'
712	CALL ioconf_setatt('UNITS', '-')
713	CALL ioconf_setatt('LONG_NAME','Soil moisture stress')
714	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g)
715	!
716	var_name= 'vegstress'
717	CALL ioconf_setatt('UNITS', '-')
718	CALL ioconf_setatt('LONG_NAME','Vegetation growth moisture stress')
719	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
720	!
721	var_name= 'bqsb'
722	CALL ioconf_setatt('UNITS', 'kg/m^2')
723	CALL ioconf_setatt('LONG_NAME','Deep soil moisture')
724	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., bqsb, "gather", nbp_glo, index_g)
725	!
726	var_name= 'gqsb'
727	CALL ioconf_setatt('UNITS', 'kg/m^2')
728	CALL ioconf_setatt('LONG_NAME','Surface soil moisture')
729	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., gqsb, "gather", nbp_glo, index_g)
730	!
731	var_name= 'dsg'
732	CALL ioconf_setatt('UNITS', 'm')
733	CALL ioconf_setatt('LONG_NAME','Depth of upper reservoir')
734	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsg, "gather", nbp_glo, index_g)
735	!
736	var_name= 'dsp'
737	CALL ioconf_setatt('UNITS', 'm')
738	CALL ioconf_setatt('LONG_NAME','Depth to lower reservoir')
739	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsp, "gather", nbp_glo, index_g)
740	!
741	var_name= 'dss'
742	CALL ioconf_setatt('UNITS', 'm')
743	CALL ioconf_setatt('LONG_NAME','Hauteur au dessus du reservoir de surface')
744	IF ( ok_var(var_name) ) THEN
745	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dss, "gather", nbp_glo, index_g)
746	ENDIF
747	!
748	var_name= 'hdry'
749	CALL ioconf_setatt('UNITS', 'm')
750	CALL ioconf_setatt('LONG_NAME','Dry soil heigth in meters')
751	IF ( ok_var(var_name) ) THEN
752	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., hdry, "gather", nbp_glo, index_g)
753	ENDIF
754	!
755	var_name= 'qsintveg'
756	CALL ioconf_setatt('UNITS', 'kg/m^2')
757	CALL ioconf_setatt('LONG_NAME','Intercepted moisture')
758	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
759	!
760	var_name= 'resdist'
761	CALL ioconf_setatt('UNITS', '-')
762	CALL ioconf_setatt('LONG_NAME','Distribution of reservoirs')
763	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
764	!
765	! get restart values if non were found in the restart file
766	!
767	!Config Key = HYDROL_SNOW
768	!Config Desc = Initial snow mass if not found in restart
769	!Config Def = 0.0
770	!Config Help = The initial value of snow mass if its value is not found
771	!Config in the restart file. This should only be used if the model is
772	!Config started without a restart file.
773	!
774	CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', 0.0_r_std)
775	!
776	!Config Key = HYDROL_SNOWAGE
777	!Config Desc = Initial snow age if not found in restart
778	!Config Def = 0.0
779	!Config Help = The initial value of snow age if its value is not found
780	!Config in the restart file. This should only be used if the model is
781	!Config started without a restart file.
782	!
783	CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', 0.0_r_std)
784	!
785	!Config Key = HYDROL_SNOW_NOBIO
786	!Config Desc = Initial snow amount on ice, lakes, etc. if not found in restart
787	!Config Def = 0.0
788	!Config Help = The initial value of snow if its value is not found
789	!Config in the restart file. This should only be used if the model is
790	!Config started without a restart file.
791	!
792	CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', 0.0_r_std)
793	!
794	!Config Key = HYDROL_SNOW_NOBIO_AGE
795	!Config Desc = Initial snow age on ice, lakes, etc. if not found in restart
796	!Config Def = 0.0
797	!Config Help = The initial value of snow age if its value is not found
798	!Config in the restart file. This should only be used if the model is
799	!Config started without a restart file.
800	!
801	CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', 0.0_r_std)
802	!
803	!Config Key = HYDROL_HUMR
804	!Config Desc = Initial soil moisture stress if not found in restart
805	!Config Def = 1.0
806	!Config Help = The initial value of soil moisture stress if its value is not found
807	!Config in the restart file. This should only be used if the model is
808	!Config started without a restart file.
809	!
810	CALL setvar_p (humrel, val_exp,'HYDROL_HUMR', 1.0_r_std)
811	CALL setvar_p (vegstress, val_exp,'HYDROL_HUMR', 1.0_r_std)
812	!
813	!Config Key = HYDROL_BQSB
814	!Config Desc = Initial restart deep soil moisture if not found in restart
815	!Config Def = DEF
816	!Config Help = The initial value of deep soil moisture if its value is not found
817	!Config in the restart file. This should only be used if the model is
818	!Config started without a restart file. Default behaviour is a saturated soil.
819	!
820	CALL setvar_p (bqsb, val_exp, 'HYDROL_BQSB', mx_eau_eau*dpu_cste)
821	!
822	!Config Key = HYDROL_GQSB
823	!Config Desc = Initial upper soil moisture if not found in restart
824	!Config Def = 0.0
825	!Config Help = The initial value of upper soil moisture if its value is not found
826	!Config in the restart file. This should only be used if the model is
827	!Config started without a restart file.
828	!
829	CALL setvar_p (gqsb, val_exp, 'HYDROL_GQSB', 0.0_r_std)
830	!
831	!Config Key = HYDROL_DSG
832	!Config Desc = Initial upper reservoir depth if not found in restart
833	!Config Def = 0.0
834	!Config Help = The initial value of upper reservoir depth if its value is not found
835	!Config in the restart file. This should only be used if the model is
836	!Config started without a restart file.
837	!
838	CALL setvar_p (dsg, val_exp, 'HYDROL_DSG', 0.0_r_std)
839
840	! set inital value for dsp if needed
841	!
842	!Config Key = HYDROL_DSP
843	!Config Desc = Initial dry soil above upper reservoir if not found in restart
844	!Config Def = DEF
845	!Config Help = The initial value of dry soil above upper reservoir if its value
846	!Config is not found in the restart file. This should only be used if
847	!Config the model is started without a restart file. The default behaviour
848	!Config is to compute it from the variables above. Should be OK most of
849	!Config the time.
850	!
851	zdsp(:,:) = dpu_cste - bqsb(:,:) / mx_eau_eau
852	dsp(1,1) = val_exp
853	call getin_p('HYDROL_DSP', dsp(1,1))
854	IF (dsp(1,1) == val_exp) THEN
855	dsp(:,:) = zdsp(:,:)
856	ELSE
857	IF (is_root_prc) &
858	ALLOCATE(zdsp_g(nbp_glo,nvm),dsp_g(nbp_glo,nvm))
859	CALL gather(zdsp,zdsp_g)
860	IF (is_root_prc) &
861	CALL setvar (dsp_g, val_exp, 'HYDROL_DSP', zdsp_g)
862	CALL scatter(dsp_g,dsp)
863	IF (is_root_prc) &
864	DEALLOCATE(zdsp_g, dsp_g)
865	ENDIF
866	!
867	!Config Key = HYDROL_QSV
868	!Config Desc = Initial water on canopy if not found in restart
869	!Config Def = 0.0
870	!Config Help = The initial value of moisture on canopy if its value
871	!Config is not found in the restart file. This should only be used if
872	!Config the model is started without a restart file.
873	!
874	CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', 0.0_r_std)
875	!
876	tmpdss = dsg - gqsb / mx_eau_eau
877	IF ( ok_var("dss") ) THEN
878	CALL setvar_p (dss, val_exp, 'NO_KEYWORD', tmpdss)
879	ELSE
880	dss(:,:)=tmpdss(:,:)
881	ENDIF
882
883	IF ( ok_var("hdry") ) THEN
884	IF (MINVAL(hdry) .EQ. MAXVAL(hdry) .AND. MAXVAL(hdry) .EQ. val_exp) THEN
885	a_subgrd(:) = zero
886	DO ji=1,kjpindex
887	IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN
888	!
889	IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN
890	! Ajouts Nathalie - Fred - le 28 Mars 2006
891	a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),0.)/dsg_min,1.)
892	!
893	ENDIF
894	ENDIF
895	!
896	ENDDO
897	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
898	! revu 22 novembre 2007
899	hdry(:) = a_subgrd(:)dss(:,1) + (1.-a_subgrd(:))dsp(:,1)
900	ENDIF
901	ELSE
902	a_subgrd(:) = zero
903	DO ji=1,kjpindex
904	IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN
905	!
906	IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN
907	! Ajouts Nathalie - Fred - le 28 Mars 2006
908	a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),0.)/dsg_min,1.)
909	!
910	ENDIF
911	ENDIF
912	!
913	ENDDO
914
915	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
916	! revu 22 novembre 2007
917	hdry(:) = a_subgrd(:)dss(:,1) + (1.-a_subgrd(:))dsp(:,1)
918	ENDIF
919	!
920	! There is no need to configure the initialisation of resdist. If not available it is the vegetation map
921	!
922	IF ( MINVAL(resdist) .EQ. MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
923	resdist = veget
924	ENDIF
925	!
926	! Remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1. Thus we need vegtot
927	!
928	DO ji = 1, kjpindex
929	vegtot(ji) = SUM(veget(ji,:))
930	ENDDO
931	!
932	ENDIF
933
934	!
935	! Where vegetation fraction is zero, set water to that of bare soil.
936	! This does not create any additional water.
937	!
938	DO jv = 2, nvm
939	DO ji = 1, kjpindex
940	IF ( veget(ji,jv) .LT. EPSILON(un) ) THEN
941	gqsb(ji,jv) = gqsb(ji,1)
942	bqsb(ji,jv) = bqsb(ji,1)
943	dsg(ji,jv) = dsg(ji,1)
944	dss(ji,jv) = dss(ji,1)
945	dsp(ji,jv) = dsp(ji,1)
946	ENDIF
947	ENDDO
948	ENDDO
949	!
950	DO ik=1, kjpindex
951	if(snow(ik).gt.maxmass_glacier) then
952	WRITE(numout,*)' il faut diminuer le stock de neige car snow > maxmass_glacier dans restart'
953	snow(ik)=maxmass_glacier
954	endif
955	if(snow_nobio(ik,iice).gt.maxmass_glacier) then
956	WRITE(numout,*)' il faut diminuer le stock de neige car snow_nobio > maxmass_glacier dans restart'
957	snow_nobio(ik,iice)=maxmass_glacier
958	endif
959	ENDDO
960	!
961	IF (long_print) WRITE (numout,*) ' hydrolc_init done '
962	!
963	END SUBROUTINE hydrolc_init
964	!
965	!-------------------------------------
966	!
967	SUBROUTINE hydrolc_clear()
968
969	l_first_hydrol=.TRUE.
970
971	IF (ALLOCATED (bqsb)) DEALLOCATE (bqsb)
972	IF (ALLOCATED (gqsb)) DEALLOCATE (gqsb)
973	IF (ALLOCATED (dsg)) DEALLOCATE (dsg)
974	IF (ALLOCATED (dsp)) DEALLOCATE (dsp)
975	IF (ALLOCATED (mean_bqsb)) DEALLOCATE (mean_bqsb)
976	IF (ALLOCATED (mean_gqsb)) DEALLOCATE (mean_gqsb)
977	IF (ALLOCATED (dss)) DEALLOCATE (dss)
978	IF (ALLOCATED (hdry)) DEALLOCATE (hdry)
979	IF (ALLOCATED (precisol)) DEALLOCATE (precisol)
980	IF (ALLOCATED (gdrainage)) DEALLOCATE (gdrainage)
981	IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg)
982	IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio)
983	IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt)
984	IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt)
985	IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var)
986	IF (ALLOCATED (ruu_ch)) DEALLOCATE (ruu_ch)
987	IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot)
988	IF (ALLOCATED (resdist)) DEALLOCATE (resdist)
989	IF (ALLOCATED (runoff)) DEALLOCATE (runoff)
990	IF (ALLOCATED (tot_water_beg)) DEALLOCATE (tot_water_beg)
991	IF (ALLOCATED (tot_water_end)) DEALLOCATE (tot_water_end)
992	IF (ALLOCATED (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
993	IF (ALLOCATED (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
994	IF (ALLOCATED (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
995	IF (ALLOCATED (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
996	IF (ALLOCATED (delsoilmoist)) DEALLOCATE (delsoilmoist)
997	IF (ALLOCATED (delintercept)) DEALLOCATE (delintercept)
998	IF (ALLOCATED (snow_beg)) DEALLOCATE (snow_beg)
999	IF (ALLOCATED (snow_end)) DEALLOCATE (snow_end)
1000	IF (ALLOCATED (delswe)) DEALLOCATE (delswe)
1001	!
1002	END SUBROUTINE hydrolc_clear
1003
1004	!! This routine initializes HYDROLOGIC variables
1005	!! - mx_eau_var
1006	!! - ruu_ch
1007	!!
1008	SUBROUTINE hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag)
1009
1010	! interface description
1011	! input scalar
1012	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1013	! input fields
1014	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
1015	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type
1016	! output fields
1017	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistance to bare soil evaporation
1018	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visible dry soil
1019	!! Profondeur du reservoir contenant le maximum d'eau
1020	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: mx_eau_var !!
1021	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ruu_ch !! Quantite d'eau maximum
1022	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out):: shumdiag !! relative soil moisture
1023	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
1024	! local declaration
1025	INTEGER(i_std) :: ji,jv, jd
1026	REAL(r_std), DIMENSION(kjpindex) :: mean_dsg
1027	REAL(r_std) :: gtr, btr
1028	REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl
1029	!
1030	! initialisation
1031	tmp_dl(1) = 0
1032	tmp_dl(2:nbdl+1) = diaglev(1:nbdl)
1033	!
1034	mx_eau_var(:) = zero
1035	!
1036	DO ji = 1,kjpindex
1037	DO jv = 1,nvm
1038	mx_eau_var(ji) = mx_eau_var(ji) + veget(ji,jv)wmax_veg(jv)dpu_cste
1039	END DO
1040	IF (vegtot(ji) .GT. zero) THEN
1041	mx_eau_var(ji) = mx_eau_var(ji)/vegtot(ji)
1042	ELSE
1043	mx_eau_var(ji) = mx_eau_eau*dpu_cste
1044	ENDIF
1045	ruu_ch(ji) = mx_eau_var(ji) / dpu_cste
1046	END DO
1047	!
1048	!
1049	! could be done with SUM instruction but this kills vectorization
1050	mean_bqsb(:) = zero
1051	mean_gqsb(:) = zero
1052	mean_dsg(:) = zero
1053	DO jv = 1, nvm
1054	DO ji = 1, kjpindex
1055	mean_bqsb(ji) = mean_bqsb(ji) + resdist(ji,jv)*bqsb(ji,jv)
1056	mean_gqsb(ji) = mean_gqsb(ji) + resdist(ji,jv)*gqsb(ji,jv)
1057	mean_dsg(ji) = mean_dsg(ji) + resdist(ji,jv)*dsg(ji,jv)
1058	ENDDO
1059	ENDDO
1060	mean_dsg(:) = MAX( mean_dsg(:), mean_gqsb(:)/ruu_ch(:) )
1061
1062	DO ji = 1, kjpindex
1063	IF (vegtot(ji) .GT. zero) THEN
1064	mean_bqsb(ji) = mean_bqsb(ji)/vegtot(ji)
1065	mean_gqsb(ji) = mean_gqsb(ji)/vegtot(ji)
1066	mean_dsg(ji) = mean_dsg(ji)/vegtot(ji)
1067	ENDIF
1068	ENDDO
1069
1070	DO jd = 1,nbdl
1071	!
1072	DO ji = 1,kjpindex
1073	!!$ !
1074	!!$ DO jd = 1,nbdl
1075	IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN
1076	shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji)
1077	ELSE
1078	IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN
1079	gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd))
1080	btr = 1 - gtr
1081	shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + &
1082	& btr*mean_bqsb(ji)/mx_eau_var(ji)
1083	ELSE
1084	shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji)
1085	ENDIF
1086	ENDIF
1087	shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero)
1088	ENDDO
1089	ENDDO
1090
1091	! The fraction of soil which is visibly dry (dry when dss = 0.1 m)
1092	drysoil_frac(:) = MIN(MAX(dss(:,1),0.)*10._r_std, un)
1093	!
1094	! Compute the resistance to bare soil evaporation
1095	!
1096	rsol(:) = -un
1097	DO ji = 1, kjpindex
1098	IF (veget(ji,1) .GE. min_sechiba) THEN
1099	!
1100	! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
1101	! on modifie le rsol pour que la resistance croisse subitement si on s'approche
1102	! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70
1103	!rsol(ji) = dss(ji,1) * rsol_cste
1104	!rsol(ji) = ( drysoil_frac(ji) + 1./(10.(dpu_cste - drysoil_frac(ji))+1.e-10)2 ) rsol_cste
1105	rsol(ji) = ( hdry(ji) + 1./(10.(dpu_cste - hdry(ji))+1.e-10)2 ) rsol_cste
1106	ENDIF
1107	ENDDO
1108
1109	!
1110	! litter humidity.
1111	!
1112
1113	!!$ DO ji = 1, kjpindex
1114	!!$ litterhumdiag(ji) = EXP( - dss(ji,1) / hcrit_litter )
1115	!!$ ENDDO
1116	litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter )
1117
1118	! special case: it has just been raining a few drops. The upper soil
1119	! reservoir is ridiculously small, but the dry soil height is zero.
1120	! Don't take it into account.
1121
1122	!!$ DO ji = 1, kjpindex
1123	!!$ IF ( ( dss(ji,1) .LT. min_sechiba ) .AND. &
1124	!!$ ( mean_dsg(ji) .GT. min_sechiba ) .AND. &
1125	!!$ ( mean_dsg(ji) .LT. 5.E-4 ) ) THEN
1126	!!$ litterhumdiag(ji) = 0.0
1127	!!$ ENDIF
1128	!!$ ENDDO
1129	WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. &
1130	( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) )
1131	litterhumdiag(:) = zero
1132	ENDWHERE
1133
1134	!
1135	IF (long_print) WRITE (numout,*) ' hydrolc_var_init done '
1136
1137	END SUBROUTINE hydrolc_var_init
1138
1139	!! This routine computes snow processes
1140	!!
1141	SUBROUTINE hydrolc_snow (kjpindex, dtradia, precip_rain, precip_snow , temp_sol_new, soilcap,&
1142	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
1143	& tot_melt, snowdepth)
1144
1145	!
1146	! interface description
1147	! input scalar
1148	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1149	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
1150	! input fields
1151	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1152	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation
1153	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! New soil temperature
1154	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
1155	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
1156	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
1157	! modified fields
1158	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation
1159	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation
1160	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
1161	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
1162	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
1163	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age!! Snow age on ice, lakes, ...
1164	! output fields
1165	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt
1166	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowdepth !! Snow depth
1167	!
1168	! local declaration
1169	!
1170	INTEGER(i_std) :: ji, jv
1171	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
1172	REAL(r_std), DIMENSION (kjpindex) :: xx !! temporary
1173	REAL(r_std) :: snowmelt_tmp !! The name says it all !
1174	LOGICAL, DIMENSION (kjpindex) :: warnings
1175	LOGICAL :: any_warning
1176	!
1177	! for continental points
1178	!
1179
1180	!
1181	! 0. initialisation
1182	!
1183	DO jv = 1, nnobio
1184	DO ji=1,kjpindex
1185	subsnownobio(ji,jv) = zero
1186	ENDDO
1187	ENDDO
1188	DO ji=1,kjpindex
1189	subsnowveg(ji) = zero
1190	snowmelt(ji) = zero
1191	icemelt(ji) = zero
1192	tot_melt(ji) = zero
1193	ENDDO
1194	!
1195	! 1. On vegetation
1196	!
1197	!cdir NODEP
1198	DO ji=1,kjpindex
1199	!
1200	! 1.1. It is snowing
1201	!
1202	snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji)
1203	ENDDO
1204	!
1205	DO ji=1,kjpindex
1206	!
1207	! 1.2. Sublimation - separate between vegetated and no-veget fractions
1208	! Care has to be taken as we might have sublimation from the
1209	! the frac_nobio while there is no snow on the rest of the grid.
1210	!
1211	IF ( snow(ji) > snowcri ) THEN
1212	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
1213	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
1214	ELSE
1215	! Correction Nathalie - Juillet 2006.
1216	! On doit d'abord tester s'il existe un frac_nobio!
1217	! Pour le moment je ne regarde que le iice
1218	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1219	subsnownobio(ji,iice) = vevapsno(ji)
1220	subsnowveg(ji) = zero
1221	ELSE
1222	subsnownobio(ji,iice) = zero
1223	subsnowveg(ji) = vevapsno(ji)
1224	ENDIF
1225	ENDIF
1226	ENDDO
1227	!
1228	warnings(:) = .FALSE.
1229	any_warning = .FALSE.
1230	!cdir NODEP
1231	DO ji=1,kjpindex
1232	!
1233	! 1.2.1 Check that sublimation on the vegetated fraction is possible.
1234	!
1235	IF (subsnowveg(ji) .GT. snow(ji)) THEN
1236	! What could not be sublimated goes into bare soil evaporation
1237	! Nathalie - Juillet 2006 - il faut avant tout tester s'il existe du
1238	! frac_nobio sur ce pixel pour eviter de puiser dans le sol!
1239	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1240	subsnownobio(ji,iice) = subsnownobio(ji,iice) + (subsnowveg(ji) - snow(ji))
1241	ELSE
1242	vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji))
1243	warnings(ji) = .TRUE.
1244	any_warning = .TRUE.
1245	ENDIF
1246	! Sublimation is thus limited to what is available
1247	subsnowveg(ji) = snow(ji)
1248	snow(ji) = zero
1249	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
1250	ELSE
1251	snow(ji) = snow(ji) - subsnowveg(ji)
1252	ENDIF
1253	ENDDO
1254	IF ( any_warning ) THEN
1255	WRITE(numout,*)' ATTENTION on prend de l eau au sol nu sur au moins un point car evapsno est trop fort!'
1256	!!$ DO ji=1,kjpindex
1257	!!$ IF ( warnings(ji) ) THEN
1258	!!$ WRITE(numout,*)' ATTENTION on prend de l eau au sol nu car evapsno est trop fort!'
1259	!!$ WRITE(numout,)' ',ji,' vevapnu (en mm/jour) = ',vevapnu(ji)one_day/dtradia
1260	!!$ ENDIF
1261	!!$ ENDDO
1262	ENDIF
1263	!
1264	warnings(:) = .FALSE.
1265	any_warning = .FALSE.
1266	!cdir NODEP
1267	DO ji=1,kjpindex
1268	!
1269	! 1.3. snow melt only if temperature positive
1270	!
1271	IF (temp_sol_new(ji).GT.tp_00) THEN
1272	!
1273	IF (snow(ji).GT.sneige) THEN
1274	!
1275	snowmelt(ji) = (1. - frac_nobio(ji,iice))(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1276	!
1277	! 1.3.1.1 enough snow for melting or not
1278	!
1279	IF (snowmelt(ji).LT.snow(ji)) THEN
1280	snow(ji) = snow(ji) - snowmelt(ji)
1281	ELSE
1282	snowmelt(ji) = snow(ji)
1283	snow(ji) = zero
1284	END IF
1285	!
1286	ELSEIF (snow(ji).GE.zero) THEN
1287	!
1288	! 1.3.2 not enough snow
1289	!
1290	snowmelt(ji) = snow(ji)
1291	snow(ji) = zero
1292	ELSE
1293	!
1294	! 1.3.3 negative snow - now snow melt
1295	!
1296	snow(ji) = zero
1297	snowmelt(ji) = zero
1298	warnings(ji) = .TRUE.
1299	any_warning = .TRUE.
1300	!
1301	END IF
1302
1303	ENDIF
1304	ENDDO
1305	IF ( any_warning ) THEN
1306	DO ji=1,kjpindex
1307	IF ( warnings(ji) ) THEN
1308	WRITE(numout,*) 'hydrolc_snow: WARNING! snow was negative and was reset to zero for point ',ji,'. '
1309	ENDIF
1310	ENDDO
1311	ENDIF
1312	!
1313	DO ji=1,kjpindex
1314	!
1315	! 1.4. Ice melt only if there is more than a given mass : maxmass_glacier,
1316	! i.e. only weight melts glaciers !
1317	! Ajouts Edouard Davin / Nathalie de Noblet add extra to melting
1318	!
1319	IF ( snow(ji) .GT. maxmass_glacier ) THEN
1320	snowmelt(ji) = snowmelt(ji) + (snow(ji) - maxmass_glacier)
1321	snow(ji) = maxmass_glacier
1322	ENDIF
1323	!
1324	END DO
1325	!
1326	! 2. On Land ice
1327	!
1328	DO ji=1,kjpindex
1329	!
1330	! 2.1. It is snowing
1331	!
1332	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
1333	& frac_nobio(ji,iice)*precip_rain(ji)
1334	!
1335	! 2.2. Sublimation - was calculated before it can give us negative snow_nobio but that is OK
1336	! Once it goes below a certain values (-maxmass_glacier for instance) we should kill
1337	! the frac_nobio(ji,iice) !
1338	!
1339	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
1340	!
1341	! 2.3. snow melt only for continental ice fraction
1342	!
1343	snowmelt_tmp = zero
1344	IF (temp_sol_new(ji) .GT. tp_00) THEN
1345	!
1346	! 2.3.1 If there is snow on the ice-fraction it can melt
1347	!
1348	snowmelt_tmp = frac_nobio(ji,iice)(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1349	!
1350	IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN
1351	snowmelt_tmp = MAX( zero, snow_nobio(ji,iice))
1352	ENDIF
1353	snowmelt(ji) = snowmelt(ji) + snowmelt_tmp
1354	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp
1355	!
1356	ENDIF
1357	!
1358	! 2.4 Ice melt only if there is more than a given mass : maxmass_glacier,
1359	! i.e. only weight melts glaciers !
1360	!
1361	IF ( snow_nobio(ji,iice) .GT. maxmass_glacier ) THEN
1362	icemelt(ji) = snow_nobio(ji,iice) - maxmass_glacier
1363	snow_nobio(ji,iice) = maxmass_glacier
1364	ENDIF
1365	!
1366	END DO
1367
1368	!
1369	! 3. On other surface types - not done yet
1370	!
1371	IF ( nnobio .GT. 1 ) THEN
1372	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
1373	CALL ipslerr (3,'hydrolc_snow', '', &
1374	& 'CANNOT TREAT SNOW ON THESE SURFACE TYPES', '')
1375	ENDIF
1376
1377	!
1378	! 4. computes total melt (snow and ice)
1379	!
1380
1381	DO ji = 1, kjpindex
1382	tot_melt(ji) = icemelt(ji) + snowmelt(ji)
1383	ENDDO
1384
1385	!
1386	! 5. computes snow age on veg and ice (for albedo)
1387	!
1388	DO ji = 1, kjpindex
1389	!
1390	! 5.1 Snow age on vegetation
1391	!
1392	IF (snow(ji) .LE. zero) THEN
1393	snow_age(ji) = zero
1394	ELSE
1395	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dtradia/one_day) &
1396	& * EXP(-precip_snow(ji) / snow_trans)
1397	ENDIF
1398	!
1399	! 5.2 Snow age on ice
1400	!
1401	! age of snow on ice: a little bit different because in cold regions, we really
1402	! cannot negect the effect of cold temperatures on snow metamorphism any more.
1403	!
1404	IF (snow_nobio(ji,iice) .LE. zero) THEN
1405	snow_nobio_age(ji,iice) = zero
1406	ELSE
1407	!
1408	d_age(ji) = ( snow_nobio_age(ji,iice) + &
1409	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dtradia/one_day ) * &
1410	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
1411	IF (d_age(ji) .GT. 0. ) THEN
1412	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
1413	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
1414	d_age(ji) = d_age(ji) / (un+xx(ji))
1415	ENDIF
1416	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
1417	!
1418	ENDIF
1419
1420	ENDDO
1421
1422	!
1423	! 6.0 Diagnose the depth of the snow layer
1424	!
1425
1426	DO ji = 1, kjpindex
1427	snowdepth(ji) = snow(ji) /sn_dens
1428	ENDDO
1429
1430	IF (long_print) WRITE (numout,*) ' hydrolc_snow done '
1431
1432	END SUBROUTINE hydrolc_snow
1433
1434	!! This routine computes canopy processes
1435	!!
1436	SUBROUTINE hydrolc_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol)
1437
1438	!
1439	! interface description
1440	!
1441	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1442	! input fields
1443	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation
1444	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: vevapwet !! Interception loss
1445	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation type
1446	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum water on vegetation for interception
1447	! modified fields
1448	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Water on vegetation due to interception
1449	! output fields
1450	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Eau tombee sur le sol
1451
1452	!
1453	! local declaration
1454	!
1455	INTEGER(i_std) :: ji, jv
1456	REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew
1457	LOGICAL, SAVE :: firstcall=.TRUE.
1458	! REAL(r_std), SAVE, DIMENSION(nvm) :: throughfall_by_pft
1459
1460	IF ( firstcall ) THEN
1461	!!$ !Config Key = PERCENT_THROUGHFALL_PFT
1462	!!$ !Config Desc = Percent by PFT of precip that is not intercepted by the canopy
1463	!!$ !Config Def = 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30. 30.
1464	!!$ !Config Help = During one rainfall event, PERCENT_THROUGHFALL_PFT% of the incident rainfall
1465	!!$ !Config will get directly to the ground without being intercepted, for each PFT.
1466	!!$
1467	!!$! throughfall_by_pft = (/ 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30., 30. /)
1468	!!$ CALL getin_p('PERCENT_THROUGHFALL_PFT',throughfall_by_pft)
1469	throughfall_by_pft = throughfall_by_pft / 100.
1470
1471	firstcall=.FALSE.
1472	ENDIF
1473
1474	! calcul de qsintmax a prevoir a chaque pas de temps
1475	! dans ini_sechiba
1476	! boucle sur les points continentaux
1477	! calcul de qsintveg au pas de temps suivant
1478	! par ajout du flux interception loss
1479	! calcule par enerbil en fonction
1480	! des calculs faits dans diffuco
1481	! calcul de ce qui tombe sur le sol
1482	! avec accumulation dans precisol
1483	! essayer d'harmoniser le traitement du sol nu
1484	! avec celui des differents types de vegetation
1485	! fait si on impose qsintmax ( ,1) = 0.0
1486	!
1487	! loop for continental subdomain
1488	!
1489
1490	!
1491	! 1. evaporation off the continents
1492	!
1493	! 1.1 The interception loss is take off the canopy.
1494
1495	DO jv=1,nvm
1496	qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv)
1497	END DO
1498
1499	! 1.2 It is raining : precip_rain is shared for each vegetation
1500	! type
1501	! sum (veget (1,nvm)) must be egal to 1-totfrac_nobio.
1502	! iniveget computes veget each day
1503	!
1504	DO jv=1,nvm
1505	! Correction Nathalie - Juin 2006 - une partie de la pluie arrivera toujours sur le sol
1506	! sorte de throughfall supplementaire
1507	!qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * precip_rain(:)
1508	qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
1509	END DO
1510
1511	!
1512	! 1.3 Limits the effect and sum what receives soil
1513	!
1514	precisol(:,:) = zero
1515	DO jv=1,nvm
1516	DO ji = 1, kjpindex
1517	zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv))
1518	! correction throughfall Nathalie - Juin 2006
1519	!precisol(ji,jv) = qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
1520	precisol(ji,jv) = (veget(ji,jv)throughfall_by_pft(jv)precip_rain(ji)) + qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
1521	ENDDO
1522	ENDDO
1523	!
1524	! 1.4 swap qsintveg to the new value
1525	!
1526
1527	DO jv=1,nvm
1528	qsintveg(:,jv) = zqsintvegnew (:,jv)
1529	END DO
1530
1531	IF (long_print) WRITE (numout,*) ' hydrolc_canop done '
1532
1533	END SUBROUTINE hydrolc_canop
1534	!!
1535	!!
1536	!!
1537	SUBROUTINE hydrolc_vegupd(kjpindex, veget, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss, dsp, resdist)
1538	!
1539	!
1540	! The vegetation cover has changed and we need to adapt the reservoir distribution.
1541	! You may note that this occurs after evaporation and so on have been computed. It is
1542	! not a problem as a new vegetation fraction will start with humrel=0 and thus will have no
1543	! evaporation. If this is not the case it should have been caught above.
1544	!
1545	! input scalar
1546	INTEGER(i_std), INTENT(in) :: kjpindex
1547	! input fields
1548	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! New vegetation map
1549	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Quantite d'eau maximum
1550	! modified fields
1551	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation due to interception
1552	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface
1553	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond
1554	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface
1555	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface
1556	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond
1557	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(inout) :: resdist !! Old vegetation map
1558	!
1559	!
1560	! local declaration
1561	!
1562	INTEGER(i_std) :: ji,jv
1563	!
1564	REAL(r_std),DIMENSION (kjpindex,nvm) :: qsintveg2 !! Water on vegetation due to interception over old veget
1565	REAL(r_std), DIMENSION (kjpindex,nvm) :: bdq, gdq, qsdq
1566	REAL(r_std), DIMENSION (kjpindex,nvm) :: vmr !! variation of veget
1567	REAL(r_std), DIMENSION(kjpindex) :: gtr, btr, vtr, qstr, fra
1568	REAL(r_std), DIMENSION(kjpindex) :: vegchtot
1569	REAL(r_std), PARAMETER :: EPS1 = EPSILON(un)
1570	!
1571	!
1572	DO jv = 1, nvm
1573	DO ji = 1, kjpindex
1574	IF ( ABS(veget(ji,jv)-resdist(ji,jv)) .GT. EPS1 ) THEN
1575	vmr(ji,jv) = veget(ji,jv)-resdist(ji,jv)
1576	ELSE
1577	vmr(ji,jv) = zero
1578	ENDIF
1579	!
1580	IF (resdist(ji,jv) .GT. 0.) THEN
1581	qsintveg2(ji,jv) = qsintveg(ji,jv)/resdist(ji,jv)
1582	ELSE
1583	qsintveg2(ji,jv) = zero
1584	ENDIF
1585	ENDDO
1586	ENDDO
1587	!
1588	vegchtot(:) = 0.
1589	DO jv = 1, nvm
1590	DO ji = 1, kjpindex
1591	vegchtot(ji) = vegchtot(ji) + ABS( vmr(ji,jv) )
1592	ENDDO
1593	ENDDO
1594	!
1595	DO jv = 1, nvm
1596	DO ji = 1, kjpindex
1597	IF ( vegchtot(ji) .GT. 0. ) THEN
1598	gdq(ji,jv) = ABS(vmr(ji,jv)) * gqsb(ji,jv)
1599	bdq(ji,jv) = ABS(vmr(ji,jv)) * bqsb(ji,jv)
1600	qsdq(ji,jv) = ABS(vmr(ji,jv)) * qsintveg2(ji,jv)
1601	ENDIF
1602	ENDDO
1603	ENDDO
1604	!
1605	! calculate water mass that we have to redistribute
1606	!
1607	gtr(:) = zero
1608	btr(:) = zero
1609	qstr(:) = zero
1610	vtr(:) = zero
1611	!
1612	!
1613	DO jv = 1, nvm
1614	DO ji = 1, kjpindex
1615	IF ( ( vegchtot(ji) .GT. 0. ) .AND. ( vmr(ji,jv) .LT. 0. ) ) THEN
1616	gtr(ji) = gtr(ji) + gdq(ji,jv)
1617	btr(ji) = btr(ji) + bdq(ji,jv)
1618	qstr(ji) = qstr(ji) + qsdq(ji,jv)
1619	vtr(ji) = vtr(ji) - vmr(ji,jv)
1620	ENDIF
1621	ENDDO
1622	ENDDO
1623	!
1624	! put it into reservoir of plant whose surface area has grown
1625	DO jv = 1, nvm
1626	DO ji = 1, kjpindex
1627	IF ( vegchtot(ji) .GT. 0. .AND. ABS(vtr(ji)) .GT. EPS1) THEN
1628	fra(ji) = vmr(ji,jv) / vtr(ji)
1629	IF ( vmr(ji,jv) .GT. 0.) THEN
1630	IF (veget(ji,jv) .GT. 0.) THEN
1631	gqsb(ji,jv) = (resdist(ji,jv)gqsb(ji,jv) + fra(ji)gtr(ji))/veget(ji,jv)
1632	bqsb(ji,jv) = (resdist(ji,jv)bqsb(ji,jv) + fra(ji)btr(ji))/veget(ji,jv)
1633	ENDIF
1634	qsintveg(ji,jv) = qsintveg(ji,jv) + fra(ji)* qstr(ji)
1635	ELSE
1636	qsintveg(ji,jv) = qsintveg(ji,jv) - qsdq(ji,jv)
1637	ENDIF
1638	!
1639	! dss is not altered so that this transfer of moisture does not directly
1640	! affect transpiration
1641	IF (gqsb(ji,jv) .LT. min_sechiba) THEN
1642	dsg(ji,jv) = zero
1643	ELSE
1644	dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) &
1645	/ ruu_ch(ji)
1646	ENDIF
1647	dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji)
1648	ENDIF
1649	ENDDO
1650	ENDDO
1651
1652	! Now that the work is done resdist needs an update !
1653	DO jv = 1, nvm
1654	resdist(:,jv) = veget(:,jv)
1655	ENDDO
1656
1657	!
1658	! Where vegetation fraction is zero, set water to that of bare soil.
1659	! This does not create any additional water.
1660	!
1661	DO jv = 2, nvm
1662	DO ji = 1, kjpindex
1663	IF ( veget(ji,jv) .LT. EPS1 ) THEN
1664	gqsb(ji,jv) = gqsb(ji,1)
1665	bqsb(ji,jv) = bqsb(ji,1)
1666	dsg(ji,jv) = dsg(ji,1)
1667	dss(ji,jv) = dss(ji,1)
1668	dsp(ji,jv) = dsp(ji,1)
1669	ENDIF
1670	ENDDO
1671	ENDDO
1672
1673	RETURN
1674	!
1675	END SUBROUTINE hydrolc_vegupd
1676	!!
1677	!! This routines computes soil processes
1678	!!
1679	SUBROUTINE hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, irrigation, tot_melt, mx_eau_var, veget, ruu_ch, transpir,&
1680	& gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, runoff, run_off_tot, drainage, humrel, vegstress, &
1681	& shumdiag, litterhumdiag)
1682	!
1683	! interface description
1684	! input scalar
1685	INTEGER(i_std), INTENT(in) :: kjpindex
1686	! input fields
1687	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: vevapnu !! Bare soil evaporation
1688	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: precisol !! Eau tombee sur le sol
1689	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Water returning to the deep reservoir
1690	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: irrigation !! Irrigation
1691	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt !! Total melt
1692	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: mx_eau_var !! Profondeur du reservoir contenant le
1693	!! maximum d'eau
1694	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! Vegetation map
1695	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: ruu_ch !! Quantite d'eau maximum
1696	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! Transpiration
1697	! modified fields
1698	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface
1699	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond
1700	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface
1701	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface
1702	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond
1703	! output fields
1704	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: runoff !! Ruissellement
1705	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: run_off_tot !! Complete runoff
1706	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! Drainage
1707	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: humrel !! Relative humidity
1708	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
1709	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! relative soil moisture
1710	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
1711	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: rsol !! resistance to bare soil evaporation
1712	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drysoil_frac !! Fraction of visible fry soil
1713	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: hdry !! Dry soil heigth in meters
1714
1715	!
1716	! local declaration
1717	!
1718	INTEGER(i_std) :: ji,jv, jd
1719	REAL(r_std), DIMENSION(kjpindex) :: zhumrel_lo, zhumrel_up
1720	REAL(r_std), DIMENSION(kjpindex,nvm) :: zeflux !! Soil evaporation
1721	REAL(r_std), DIMENSION(kjpindex,nvm) :: zpreci !! Soil precipitation
1722	LOGICAL, DIMENSION(kjpindex,nvm) :: warning
1723	REAL(r_std) :: gtr, btr
1724	REAL(r_std), DIMENSION(kjpindex) :: mean_dsg
1725	LOGICAL :: OnceMore
1726	INTEGER(i_std), PARAMETER :: nitermax = 100
1727	INTEGER(i_std) :: niter
1728	INTEGER(i_std) :: nbad
1729	LOGICAL, DIMENSION(kjpindex,nvm) :: lbad
1730	LOGICAL, DIMENSION(kjpindex) :: lbad_ij
1731	REAL(r_std) :: gqseuil , eausup, wd1
1732	REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl
1733	! Ajout Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
1734	! Modifs stabilite
1735	REAL(r_std), DIMENSION(kjpindex,nvm) :: a_subgrd
1736	!
1737	! 0. we have only one flux field corresponding to water evaporated from the surface
1738	!
1739	DO jv=1,nvm
1740	DO ji = 1, kjpindex
1741	IF ( veget(ji,jv) .GT. zero ) THEN
1742	zeflux(ji,jv) = transpir(ji,jv)/veget(ji,jv)
1743	zpreci(ji,jv) = precisol(ji,jv)/veget(ji,jv)
1744	ELSE
1745	zeflux(ji,jv) = zero
1746	zpreci(ji,jv) = zero
1747	ENDIF
1748	ENDDO
1749	ENDDO
1750	!
1751	! We need a test on the bare soil fraction because we can have bare soil evaporation even when
1752	! there is no bare soil because of transfers (snow for instance). This should only apply if there
1753	! is vegetation but we do not test this case.
1754	!
1755
1756	! case 1 - there is vegetation and bare soil
1757	DO ji = 1, kjpindex
1758	IF ( (vegtot(ji) .GT. zero) .AND. (veget(ji,1) .GT. min_sechiba) ) THEN
1759	zeflux(ji,1) = vevapnu(ji)/veget(ji,1)
1760	ENDIF
1761	ENDDO
1762
1763	! case 2 - there is vegetation but no bare soil
1764	DO jv = 1, nvm
1765	DO ji = 1, kjpindex
1766	IF ( (vegtot(ji) .GT. zero) .AND. (veget(ji,1) .LE. min_sechiba)&
1767	& .AND. (veget(ji,jv) .GT. min_sechiba)) THEN
1768	zeflux(ji,jv) = zeflux(ji,jv) + vevapnu(ji)/vegtot(ji)
1769	ENDIF
1770	ENDDO
1771	ENDDO
1772
1773	!
1774	! 0.1 Other temporary variables
1775	!
1776	tmp_dl(1) = 0
1777	tmp_dl(2:nbdl+1) = diaglev(1:nbdl)
1778	!
1779
1780	!
1781	! 1.1 Transpiration for each vegetation type
1782	! Evaporated water is taken out of the ground.
1783	!
1784	!
1785	DO jv=1,nvm
1786	DO ji=1,kjpindex
1787	!
1788	gqsb(ji,jv) = gqsb(ji,jv) - zeflux(ji,jv)
1789	!
1790	! 1.2 Add snow and ice melt, troughfall from vegetation and irrigation.
1791	!
1792	IF(vegtot(ji) .NE. zero) THEN
1793	! snow and ice melt and troughfall from vegetation
1794	gqsb(ji,jv) = gqsb(ji,jv) + zpreci(ji,jv) + tot_melt(ji)/vegtot(ji)
1795	!
1796	! We take care to add the irrigation only to the vegetated part if possible
1797	!
1798	IF (ABS(vegtot(ji)-veget(ji,1)) .LE. min_sechiba) THEN
1799	gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/vegtot(ji)
1800	ELSE
1801	IF ( jv > 1 ) THEN
1802	! Only add the irrigation to the upper soil if there is a reservoir.
1803	! Without this the water evaporates right away.
1804	IF ( gqsb(ji,jv) > zero ) THEN
1805	gqsb(ji,jv) = gqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-veget(ji,1))
1806	ELSE
1807	bqsb(ji,jv) = bqsb(ji,jv) + irrigation(ji)/(vegtot(ji)-veget(ji,1))
1808	ENDIF
1809	ENDIF
1810	ENDIF
1811	!
1812	! 1.3 We add the water returned from rivers to the lower reservoir.
1813	!
1814	bqsb(ji,jv) = bqsb(ji,jv) + returnflow(ji)/vegtot(ji)
1815	!
1816	ENDIF
1817	!
1818	END DO
1819	ENDDO
1820	!
1821	! 1.3 Computes run-off
1822	!
1823	runoff(:,:) = zero
1824	!
1825	! 1.4 Soil moisture is updated
1826	!
1827	warning(:,:) = .FALSE.
1828	DO jv=1,nvm
1829	DO ji = 1, kjpindex
1830	!
1831	runoff(ji,jv) = MAX(gqsb(ji,jv) + bqsb(ji,jv) - mx_eau_var(ji), zero)
1832	!
1833	IF (mx_eau_var(ji) .LE. (gqsb(ji,jv) + bqsb(ji,jv))) THEN
1834	!
1835	! 1.4.1 Plus de reservoir de surface: le sol est sature
1836	! d'eau. Le reservoir de surface est inexistant
1837	! Tout est dans le reservoir de fond.
1838	! Le ruissellement correspond a ce qui deborde.
1839	!
1840	gqsb(ji,jv) = zero
1841	dsg(ji,jv) = zero
1842	bqsb(ji,jv) = mx_eau_var(ji)
1843	dsp(ji,jv) = zero
1844	dss(ji,jv) = dsp (ji,jv)
1845
1846	ELSEIF ((gqsb(ji,jv) + bqsb(ji,jv)).GE.zero) THEN
1847	!
1848	IF (gqsb(ji,jv) .GT. dsg(ji,jv) * ruu_ch(ji)) THEN
1849	!
1850	! 1.4.2 On agrandit le reservoir de surface
1851	! car il n y a pas eu ruissellement
1852	! et toute l'eau du reservoir de surface
1853	! tient dans un reservoir dont la taille
1854	! est plus importante que l actuel.
1855	! La hauteur de sol sec dans le reservoir
1856	! de surface est alors nulle.
1857	! Le reste ne bouge pas.
1858	!
1859	dsg(ji,jv) = gqsb(ji,jv) / ruu_ch(ji)
1860	dss(ji,jv) = zero
1861	ELSEIF (gqsb(ji,jv) .GT. zero ) THEN
1862	!
1863	! 1.4.3 L eau tient dans le reservoir de surface
1864	! tel qu il existe.
1865	! Calcul de la nouvelle hauteur de sol sec
1866	! dans la couche de surface.
1867	! Le reste ne bouge pas.
1868	!
1869	dss(ji,jv) = ((dsg(ji,jv) * ruu_ch(ji)) - gqsb(ji,jv)) / ruu_ch(ji)
1870	ELSE
1871	!
1872	! 1.4.4 La quantite d eau du reservoir de surface
1873	! est negative. Cela revient a enlever
1874	! cette quantite au reservoir profond.
1875	! Le reservoir de surface est alors vide.
1876	! (voir aussi 1.4.1)
1877	!
1878	bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv)
1879	dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji)
1880	gqsb(ji,jv) = zero
1881	dsg(ji,jv) = zero
1882	dss(ji,jv) = dsp(ji,jv)
1883	END IF
1884
1885	ELSE
1886	!
1887	! 1.4.5 Le reservoir profond est aussi asseche.
1888	! La quantite d eau a enlever depasse la quantite
1889	! disponible dans le reservoir profond.
1890	!
1891	!
1892	! Ceci ne devrait jamais arriver plus d'une fois par point. C-a-d une fois la valeur negative
1893	! atteinte les flux doivent etre nuls. On ne signal que ce cas donc.
1894	!
1895	IF ( ( zeflux(ji,jv) .GT. zero ) .AND. &
1896	( gqsb(ji,jv) + bqsb(ji,jv) .LT. -1.e-10 ) ) THEN
1897	warning(ji,jv) = .TRUE.
1898	! WRITE (numout,*) 'WARNING! Soil Moisture will be negative'
1899	! WRITE (numout,*) 'ji, jv = ', ji,jv
1900	! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji)
1901	! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv)
1902	! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv)
1903	! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv)
1904	! WRITE (numout,*) 'dss = ', dss(ji,jv)
1905	! WRITE (numout,*) 'dsg = ', dsg(ji,jv)
1906	! WRITE (numout,*) 'dsp = ', dsp(ji,jv)
1907	! WRITE (numout,*) 'humrel = ', humrel(ji, jv)
1908	! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv)
1909	! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji)
1910	! WRITE (numout,*) '============================'
1911	ENDIF
1912	!
1913	bqsb(ji,jv) = gqsb(ji,jv) + bqsb(ji,jv)
1914	dsp(ji,jv) = dpu_cste
1915	gqsb(ji,jv) = zero
1916	dsg(ji,jv) = zero
1917	dss(ji,jv) = dsp(ji,jv)
1918	!
1919	ENDIF
1920	!
1921	ENDDO
1922	ENDDO
1923	!
1924	nbad = COUNT( warning(:,:) .EQV. .TRUE. )
1925	IF ( nbad .GT. 0 ) THEN
1926	WRITE(numout,*) 'hydrolc_soil: WARNING! Soil moisture was negative at', &
1927	nbad, ' points:'
1928	!DO jv = 1, nvm
1929	! DO ji = 1, kjpindex
1930	! IF ( warning(ji,jv) ) THEN
1931	! WRITE(numout,*) ' ji,jv = ', ji,jv
1932	! WRITE (numout,*) 'mx_eau_var = ', mx_eau_var(ji)
1933	! WRITE (numout,*) 'veget, resdist =', veget(ji,jv), resdist(ji,jv)
1934	! WRITE (numout,*) 'bqsb = ', bqsb(ji,jv)
1935	! WRITE (numout,*) 'gqsb = ', gqsb(ji,jv)
1936	! WRITE (numout,*) 'runoff = ',runoff(ji,jv)
1937	! WRITE (numout,*) 'dss = ', dss(ji,jv)
1938	! WRITE (numout,*) 'dsg = ', dsg(ji,jv)
1939	! WRITE (numout,*) 'dsp = ', dsp(ji,jv)
1940	! WRITE (numout,*) 'humrel = ', humrel(ji, jv)
1941	! WRITE (numout,*) 'Soil evaporation = ', zeflux(ji,jv)
1942	! WRITE (numout,*) 'Soil precipitation = ',zpreci(ji,jv)
1943	! WRITE (numout,*) 'input = ',precisol(ji, jv), tot_melt(ji)
1944	! WRITE (numout,*) 'returnflow = ',returnflow(ji)
1945	! WRITE (numout,*) '============================'
1946	! ENDIF
1947	! ENDDO
1948	!ENDDO
1949	ENDIF
1950	!
1951	! 2.0 Very large upper reservoirs or very close upper and lower reservoirs
1952	! can be deadlock situations for Choisnel. They are handled here
1953	!
1954	IF (long_print) WRITE(numout,*) 'hydro_soil 2.0 : Resolve deadlocks'
1955	DO jv=1,nvm
1956	DO ji=1,kjpindex
1957	!
1958	! 2.1 the two reservoirs are very close to each other
1959	!
1960	IF ( ABS(dsp(ji,jv)-dsg(ji,jv)) .LT. min_resdis ) THEN
1961	bqsb(ji,jv) = bqsb(ji,jv) + gqsb(ji,jv)
1962	dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji)
1963	gqsb(ji,jv) = zero
1964	dsg(ji,jv) = zero
1965	dss(ji,jv) = dsp(ji,jv)
1966	ENDIF
1967	!
1968	! 2.2 Draine some water from the upper to the lower reservoir
1969	!
1970	gqseuil = min_resdis * deux * ruu_ch(ji)
1971	eausup = dsg(ji,jv) * ruu_ch(ji)
1972	wd1 = .75*eausup
1973	!
1974	IF (eausup .GT. gqseuil) THEN
1975	gdrainage(ji,jv) = min_drain * (gqsb(ji,jv)/eausup)
1976	!
1977	IF ( gqsb(ji,jv) .GE. wd1 .AND. dsg(ji,jv) .GT. 0.10 ) THEN
1978	!
1979	gdrainage(ji,jv) = gdrainage(ji,jv) + &
1980	(max_drain-min_drain)((gqsb(ji,jv)-wd1) / (eausup-wd1))*exp_drain
1981	!
1982	ENDIF
1983	!
1984	gdrainage(ji,jv)=MIN(gdrainage(ji,jv), MAX(gqsb(ji,jv), zero))
1985	!
1986	ELSE
1987	gdrainage(ji,jv)=zero
1988	ENDIF
1989	!
1990	gqsb(ji,jv) = gqsb(ji,jv) - gdrainage(ji,jv)
1991	bqsb(ji,jv) = bqsb(ji,jv) + gdrainage(ji,jv)
1992	dsg(ji,jv) = dsg(ji,jv) - gdrainage(ji,jv) / ruu_ch(ji)
1993	dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji)
1994	!
1995	!
1996	ENDDO
1997	!
1998	ENDDO
1999
2000	!
2001	! 3.0 Diffusion of water between the reservoirs of the different plants
2002	!
2003	IF (long_print) WRITE(numout,*) 'hydrolc_soil 3.0 : Vertical diffusion'
2004
2005	mean_bqsb(:) = 0.
2006	mean_gqsb(:) = 0.
2007	DO jv = 1, nvm
2008	DO ji = 1, kjpindex
2009	IF ( vegtot(ji) .GT. zero ) THEN
2010	mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv)
2011	mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv)
2012	ENDIF
2013	ENDDO
2014	ENDDO
2015
2016	OnceMore = .TRUE.
2017	niter = 0
2018
2019	!YM
2020	lbad_ij(:)=.TRUE.
2021
2022	! nitermax prevents infinite loops (should actually never occur)
2023	DO WHILE ( OnceMore .AND. ( niter .LT. nitermax ) )
2024	!
2025	niter = niter + 1
2026
2027	!YM
2028	DO jv = 1, nvm
2029	!
2030	DO ji = 1, kjpindex
2031	IF (lbad_ij(ji)) THEN
2032	IF ( veget(ji,jv) .GT. 0. ) THEN
2033	!
2034	bqsb(ji,jv) = mean_bqsb(ji)
2035	dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji)
2036	ENDIF
2037	ENDIF
2038	!
2039	ENDDO
2040	ENDDO
2041	!
2042	! where do we have to do something?
2043	lbad(:,:) = ( ( dsp(:,:) .LT. dsg(:,:) ) .AND. &
2044	( dsg(:,:) .GT. zero ) .AND. &
2045	( veget(:,:) .GT. zero ) )
2046	!
2047	! if there are no such points any more, we'll do no more iteration
2048	IF ( COUNT( lbad(:,:) ) .EQ. 0 ) OnceMore = .FALSE.
2049
2050	DO ji = 1, kjpindex
2051	IF (COUNT(lbad(ji,:)) == 0 ) lbad_ij(ji)=.FALSE.
2052	ENDDO
2053	!
2054	DO jv = 1, nvm
2055	!YM
2056	! !
2057	! DO ji = 1, kjpindex
2058	! IF ( veget(ji,jv) .GT. 0. ) THEN
2059	! !
2060	! bqsb(ji,jv) = mean_bqsb(ji)
2061	! dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji)
2062	! ENDIF
2063	! !
2064	! ENDDO
2065	!
2066	DO ji = 1, kjpindex
2067	IF ( lbad(ji,jv) ) THEN
2068	!
2069	runoff(ji,jv) = runoff(ji,jv) + &
2070	MAX( bqsb(ji,jv) + gqsb(ji,jv) - mx_eau_var(ji), zero)
2071	!
2072	bqsb(ji,jv) = MIN( bqsb(ji,jv) + gqsb(ji,jv), mx_eau_var(ji))
2073	!
2074	gqsb(ji,jv) = zero
2075	dsp(ji,jv) = dpu_cste - bqsb(ji,jv)/ruu_ch(ji)
2076	dss(ji,jv) = dsp(ji,jv)
2077	dsg(ji,jv) = zero
2078	!
2079	ENDIF
2080	ENDDO
2081	!
2082	ENDDO
2083	!
2084	mean_bqsb(:) = 0.
2085	mean_gqsb(:) = 0.
2086	DO jv = 1, nvm
2087	DO ji = 1, kjpindex
2088	IF ( vegtot(ji) .GT. zero ) THEN
2089	mean_bqsb(ji) = mean_bqsb(ji) + veget(ji,jv)/vegtot(ji)*bqsb(ji,jv)
2090	mean_gqsb(ji) = mean_gqsb(ji) + veget(ji,jv)/vegtot(ji)*gqsb(ji,jv)
2091	ENDIF
2092	ENDDO
2093	ENDDO
2094	!
2095	ENDDO
2096
2097	!
2098	! 4. computes total runoff
2099	!
2100	IF (long_print) WRITE(numout,*) 'hydrolc_soil 4.0: Computes total runoff'
2101
2102	run_off_tot(:) = zero
2103	DO ji = 1, kjpindex
2104	IF ( vegtot(ji) .GT. zero ) THEN
2105	DO jv = 1, nvm
2106	run_off_tot(ji) = run_off_tot(ji) + (runoff(ji,jv)*veget(ji,jv))
2107	ENDDO
2108	ELSE
2109	run_off_tot(ji) = tot_melt(ji) + irrigation(ji)
2110	ENDIF
2111	ENDDO
2112	!
2113	! 4.1 We estimate some drainage !
2114	!
2115	drainage(:) = 0.95 * run_off_tot(:)
2116	run_off_tot(:) = run_off_tot(:) - drainage(:)
2117	!
2118	! 5. Some averaged diagnostics
2119	!
2120	IF (long_print) WRITE(numout,*) 'hydro_soil 5.0: Diagnostics'
2121
2122	!
2123	! 5.1 reset dsg if necessary
2124	!
2125	WHERE (gqsb(:,:) .LE. zero) dsg(:,:) = zero
2126	!
2127	DO ji=1,kjpindex
2128	!
2129	! 5.2 Compute an average moisture profile
2130	!
2131	mean_dsg(ji) = mean_gqsb(ji)/ruu_ch(ji)
2132	!
2133	ENDDO
2134
2135	!
2136	! 6. Compute the moisture stress on vegetation
2137	!
2138	IF (long_print) WRITE(numout,*) 'hydro_soil 6.0 : Moisture stress'
2139
2140	a_subgrd(:,:) = zero
2141
2142	DO jv = 1, nvm
2143	DO ji=1,kjpindex
2144	!
2145	! computes relative surface humidity
2146	!
2147	! Only use the standard formulas if total soil moisture is larger than zero.
2148	! Else stress functions are set to zero.
2149	! This will avoid that large negative soil moisture accumulate over time by the
2150	! the creation of small skin reservoirs which evaporate quickly.
2151	!
2152	IF ( gqsb(ji,jv)+bqsb(ji,jv) .GT. zero ) THEN
2153	!
2154	IF (dsg(ji,jv).EQ. zero .OR. gqsb(ji,jv).EQ.zero) THEN
2155	humrel(ji,jv) = EXP( - humcste(jv) * dpu_cste * (dsp(ji,jv)/dpu_cste) )
2156	dsg(ji,jv) = zero
2157	!
2158	! if the dry soil height is larger than the one corresponding
2159	! to the wilting point, or negative lower soil moisture : humrel is 0.0
2160	!
2161	IF (dsp(ji,jv).GT.(dpu_cste - (qwilt / ruu_ch(ji))) .OR. bqsb(ji,jv).LT.zero) THEN
2162	humrel(ji,jv) = zero
2163	ENDIF
2164	!
2165	! In this case we can take for vegetation growth the same values for humrel and vegstress
2166	!
2167	vegstress(ji,jv) = humrel(ji,jv)
2168	!
2169	ELSE
2170
2171	! Corrections Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
2172	!zhumrel_lo(ji) = EXP( - humcste(jv) * dpu_cste * (dsp(ji,jv)/dpu_cste) )
2173	!zhumrel_up(ji) = EXP( - humcste(jv) * dpu_cste * (dss(ji,jv)/dsg(ji,jv)) )
2174	!humrel(ji,jv) = MAX(zhumrel_lo(ji),zhumrel_up(ji))
2175	!
2176	! As we need a slower variable for vegetation growth the stress is computed
2177	! differently than in humrel.
2178	!
2179	zhumrel_lo(ji) = EXP( - humcste(jv) * dsp(ji,jv))
2180	zhumrel_up(ji) = EXP( - humcste(jv) * dss(ji,jv))
2181	! Ajouts Nathalie - Fred - le 28 Mars 2006
2182	a_subgrd(ji,jv)=MIN(MAX(dsg(ji,jv)-dss(ji,jv),0.)/dsg_min,1.)
2183	humrel(ji,jv)=a_subgrd(ji,jv)zhumrel_up(ji)+(1.-a_subgrd(ji,jv))zhumrel_lo(ji)
2184	!
2185	vegstress(ji,jv) = zhumrel_lo(ji) + zhumrel_up(ji) - EXP( - humcste(jv) * dsg(ji,jv) )
2186	!
2187	ENDIF
2188	!
2189	ELSE
2190	!
2191	humrel(ji,jv) = zero
2192	vegstress(ji,jv) = zero
2193	!
2194	ENDIF
2195	!
2196	ENDDO
2197	ENDDO
2198
2199
2200	!!$DS Debug humrel 14/02/2011
2201	!!$ DO jv = 1, nvm
2202	!!$ DO ji=1,kjpindex
2203	!!$ WRITE(numout,*)'the value of humrel for the PFT 10 is : ',humrel(ji,10)
2204	!!$ WRITE(numout,*)'the value of humrel for the PFT 12 is : ',humrel(ji,12)
2205	!!$ ENDDO
2206	!!$ ENDDO
2207
2208	!
2209	! 7. Diagnostics which are needed to carry information to other modules
2210	!
2211	! 7.2 Relative soil moisture
2212	!
2213	DO jd = 1,nbdl
2214	DO ji = 1, kjpindex
2215	IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN
2216	shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji)
2217	ELSE
2218	IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN
2219	gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd))
2220	btr = 1 - gtr
2221	shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + &
2222	& btr*mean_bqsb(ji)/mx_eau_var(ji)
2223	ELSE
2224	shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji)
2225	ENDIF
2226	ENDIF
2227	shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero)
2228	ENDDO
2229	ENDDO
2230
2231	! The fraction of visibly dry soil (dry when dss(:,1) = 0.1 m)
2232	drysoil_frac(:) = MIN(MAX(dss(:,1),0.)*10._r_std, un)
2233
2234	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
2235	! revu 22 novembre 2007
2236	hdry(:) = a_subgrd(:,1)dss(:,1) + (1.-a_subgrd(:,1))dsp(:,1)
2237	!
2238	! Compute the resistance to bare soil evaporation.
2239	!
2240	rsol(:) = -un
2241	DO ji = 1, kjpindex
2242	IF (veget(ji,1) .GE. min_sechiba) THEN
2243	!
2244	! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
2245	! on modifie le rsol pour que la resistance croisse subitement si on s'approche
2246	! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70
2247	!rsol(ji) = dss(ji,1) * rsol_cste
2248	rsol(ji) = ( hdry(ji) + 1./(10.(dpu_cste - hdry(ji))+1.e-10)2 ) rsol_cste
2249	ENDIF
2250	ENDDO
2251
2252	!
2253	! 8. litter humidity.
2254	!
2255
2256	litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter )
2257
2258	! special case: it has just been raining a few drops. The upper soil
2259	! reservoir is ridiculously small, but the dry soil height is zero.
2260	! Don't take it into account.
2261
2262	WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. &
2263	( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) )
2264	litterhumdiag(:) = zero
2265	ENDWHERE
2266
2267	!
2268	IF (long_print) WRITE (numout,*) ' hydrolc_soil done '
2269
2270	END SUBROUTINE hydrolc_soil
2271	!!
2272	!! This routines checks the water balance. First it gets the total
2273	!! amount of water and then it compares the increments with the fluxes.
2274	!! The computation is only done over the soil area as over glaciers (and lakes?)
2275	!! we do not have water conservation.
2276	!!
2277	!! This verification does not make much sense in REAL*4 as the precision is the same as some
2278	!! of the fluxes
2279	!!
2280	SUBROUTINE hydrolc_waterbal (kjpindex, index, first_call, dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
2281	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno, run_off_tot, &
2282	& drainage)
2283	!
2284	!
2285	!
2286	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
2287	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
2288	LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ?
2289	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
2290	!
2291	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
2292	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio!! Total fraction of continental ice+lakes+...
2293	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
2294	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
2295	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !!Ice water balance
2296	!
2297	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
2298	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
2299	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Water returning from routing to the deep reservoir
2300	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation
2301	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt
2302	!
2303	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss
2304	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
2305	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapnu !! Bare soil evaporation
2306	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapsno !! Snow evaporation
2307	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: run_off_tot !! Total runoff
2308	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: drainage !! Drainage
2309	!
2310	! LOCAL
2311	!
2312	INTEGER(i_std) :: ji, jv, jn
2313	REAL(r_std) :: allowed_err
2314	REAL(r_std),DIMENSION (kjpindex) :: watveg, delta_water, tot_flux, sum_snow_nobio, sum_vevapwet, sum_transpir
2315	!
2316	!
2317	!
2318	IF ( first_call ) THEN
2319
2320	tot_water_beg(:) = zero
2321	watveg(:) = zero
2322	sum_snow_nobio(:) = zero
2323	!cdir NODEP
2324	DO jv = 1, nvm
2325	watveg(:) = watveg(:) + qsintveg(:,jv)
2326	ENDDO
2327	!cdir NODEP
2328	DO jn = 1, nnobio
2329	sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn)
2330	ENDDO
2331	!cdir NODEP
2332	DO ji = 1, kjpindex
2333	tot_water_beg(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + &
2334	& watveg(ji) + snow(ji) + sum_snow_nobio(ji)
2335	ENDDO
2336	tot_water_end(:) = tot_water_beg(:)
2337
2338	RETURN
2339
2340	ENDIF
2341	!
2342	! Check the water balance
2343	!
2344	tot_water_end(:) = zero
2345	!
2346	DO ji = 1, kjpindex
2347	!
2348	! If the fraction of ice, lakes, etc. does not complement the vegetation fraction then we do not
2349	! need to go any further
2350	!
2351	! Modif Nathalie
2352	! IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. EPSILON(un) ) THEN
2353	IF ( (un - (totfrac_nobio(ji) + vegtot(ji))) .GT. (100*EPSILON(un)) ) THEN
2354	WRITE(numout,*) 'HYDROL problem in vegetation or frac_nobio on point ', ji
2355	WRITE(numout,*) 'totfrac_nobio : ', totfrac_nobio(ji)
2356	WRITE(numout,*) 'vegetation fraction : ', vegtot(ji)
2357	!STOP 'in hydrolc_waterbal'
2358	ENDIF
2359	ENDDO
2360	!
2361	watveg(:) = zero
2362	sum_vevapwet(:) = zero
2363	sum_transpir(:) = zero
2364	sum_snow_nobio(:) = zero
2365	!cdir NODEP
2366	DO jv = 1,nvm
2367	watveg(:) = watveg(:) + qsintveg(:,jv)
2368	sum_vevapwet(:) = sum_vevapwet(:) + vevapwet(:,jv)
2369	sum_transpir(:) = sum_transpir(:) + transpir(:,jv)
2370	ENDDO
2371	!cdir NODEP
2372	DO jn = 1,nnobio
2373	sum_snow_nobio(:) = sum_snow_nobio(:) + snow_nobio(:,jn)
2374	ENDDO
2375	!
2376	!cdir NODEP
2377	DO ji = 1, kjpindex
2378	tot_water_end(ji) = (mean_bqsb(ji) + mean_gqsb(ji))*vegtot(ji) + &
2379	& watveg(ji) + snow(ji) + sum_snow_nobio(ji)
2380	ENDDO
2381	!
2382	DO ji = 1, kjpindex
2383	!
2384	delta_water(ji) = tot_water_end(ji) - tot_water_beg(ji)
2385	!
2386	tot_flux(ji) = precip_rain(ji) + precip_snow(ji) + returnflow(ji) + irrigation(ji) - &
2387	& sum_vevapwet(ji) - sum_transpir(ji) - vevapnu(ji) - vevapsno(ji) - &
2388	& run_off_tot(ji) - drainage(ji)
2389	!
2390	ENDDO
2391	!
2392	! Set some precision ! This is a wild guess and corresponds to what works on an IEEE machine
2393	! under double precision (REAL*8).
2394	!
2395	allowed_err = 50000*EPSILON(un)
2396	!
2397	DO ji = 1, kjpindex
2398	IF ( ABS(delta_water(ji)-tot_flux(ji)) .GT. allowed_err ) THEN
2399	WRITE(numout,*) 'HYDROL does not conserve water. The erroneous point is : ', ji
2400	WRITE(numout,*) 'The error in mm/d is :', (delta_water(ji)-tot_flux(ji))/dtradia, &
2401	& ' and in mm/dt : ', delta_water(ji)-tot_flux(ji)
2402	WRITE(numout,*) 'delta_water : ', delta_water(ji), ' tot_flux : ', tot_flux(ji)
2403	WRITE(numout,*) 'Actual and allowed error : ', ABS(delta_water(ji)-tot_flux(ji)), allowed_err
2404	WRITE(numout,*) 'vegtot : ', vegtot(ji)
2405	WRITE(numout,*) 'precip_rain : ', precip_rain(ji)
2406	WRITE(numout,*) 'precip_snow : ', precip_snow(ji)
2407	WRITE(numout,*) 'Water from irrigation : ', returnflow(ji), irrigation(ji)
2408	WRITE(numout,*) 'Total water in soil :', mean_bqsb(ji) + mean_gqsb(ji)
2409	WRITE(numout,*) 'Water on vegetation :', watveg(ji)
2410	WRITE(numout,*) 'Snow mass :', snow(ji)
2411	WRITE(numout,*) 'Snow mass on ice :', sum_snow_nobio(ji)
2412	WRITE(numout,*) 'Melt water :', tot_melt(ji)
2413	WRITE(numout,*) 'evapwet : ', vevapwet(ji,:)
2414	WRITE(numout,*) 'transpir : ', transpir(ji,:)
2415	WRITE(numout,*) 'evapnu, evapsno : ', vevapnu(ji), vevapsno(ji)
2416
2417	WRITE(numout,*) 'drainage : ', drainage(ji)
2418	STOP 'in hydrolc_waterbal'
2419	ENDIF
2420	!
2421	ENDDO
2422	!
2423	! Transfer the total water amount at the end of the current timestep top the begining of the next one.
2424	!
2425	tot_water_beg = tot_water_end
2426	!
2427	END SUBROUTINE hydrolc_waterbal
2428	!
2429	! This routine computes the changes in soil moisture and interception storage for the ALMA outputs
2430	!
2431	SUBROUTINE hydrolc_alma (kjpindex, index, first_call, qsintveg, snow, snow_nobio, soilwet)
2432	!
2433	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
2434	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
2435	LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ?
2436	!
2437	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
2438	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent
2439	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilwet !! Soil wetness
2440	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2]
2441	!
2442	! LOCAL
2443	!
2444	INTEGER(i_std) :: ji
2445	REAL(r_std) :: watveg
2446	!
2447	!
2448	!
2449	IF ( first_call ) THEN
2450
2451	tot_watveg_beg(:) = zero
2452	tot_watsoil_beg(:) = zero
2453	snow_beg(:) = zero
2454	!
2455	DO ji = 1, kjpindex
2456	watveg = SUM(qsintveg(ji,:))
2457	tot_watveg_beg(ji) = watveg
2458	tot_watsoil_beg(ji) = mean_bqsb(ji) + mean_gqsb(ji)
2459	snow_beg(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
2460	ENDDO
2461	!
2462	tot_watveg_end(:) = tot_watveg_beg(:)
2463	tot_watsoil_end(:) = tot_watsoil_beg(:)
2464	snow_end(:) = snow_beg(:)
2465
2466	RETURN
2467
2468	ENDIF
2469	!
2470	! Calculate the values for the end of the time step
2471	!
2472	tot_watveg_end(:) = zero
2473	tot_watsoil_end(:) = zero
2474	snow_end(:) = zero
2475	delintercept(:) = zero
2476	delsoilmoist(:) = zero
2477	delswe(:) = zero
2478	!
2479	DO ji = 1, kjpindex
2480	watveg = SUM(qsintveg(ji,:))
2481	tot_watveg_end(ji) = watveg
2482	tot_watsoil_end(ji) = mean_bqsb(ji) + mean_gqsb(ji)
2483	snow_end(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
2484	!
2485	delintercept(ji) = tot_watveg_end(ji) - tot_watveg_beg(ji)
2486	delsoilmoist(ji) = tot_watsoil_end(ji) - tot_watsoil_beg(ji)
2487	delswe(ji) = snow_end(ji) - snow_beg(ji)
2488	!
2489	!
2490	ENDDO
2491	!
2492	!
2493	! Transfer the total water amount at the end of the current timestep top the begining of the next one.
2494	!
2495	tot_watveg_beg = tot_watveg_end
2496	tot_watsoil_beg = tot_watsoil_end
2497	snow_beg(:) = snow_end(:)
2498	!
2499	DO ji = 1,kjpindex
2500	soilwet(ji) = tot_watsoil_end(ji) / mx_eau_var(ji)
2501	ENDDO
2502	!
2503	END SUBROUTINE hydrolc_alma
2504	!-
2505	SUBROUTINE hydrolc_hdiff(kjpindex, dtradia, veget, ruu_ch, gqsb, bqsb, dsg, dss, dsp)
2506
2507	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
2508	REAL(r_std), INTENT (in) :: dtradia !! time step (s)
2509	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
2510	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ruu_ch !! Quantite d'eau maximum
2511
2512	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: gqsb !! Hauteur d'eau dans le reservoir de surface
2513	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: bqsb !! Hauteur d'eau dans le reservoir profond
2514	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsg !! Hauteur du reservoir de surface
2515	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dss !! Hauteur au dessus du reservoir de surface
2516	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: dsp !! Hauteur au dessus du reservoir profond
2517
2518	REAL(r_std), DIMENSION (kjpindex) :: bqsb_mean
2519	REAL(r_std), DIMENSION (kjpindex) :: gqsb_mean
2520	REAL(r_std), DIMENSION (kjpindex) :: dss_mean
2521	REAL(r_std), DIMENSION (kjpindex) :: vegtot
2522	REAL(r_std) :: x
2523	INTEGER(i_std) :: ji,jv
2524	REAL(r_std), SAVE :: tau_hdiff
2525	LOGICAL, SAVE :: firstcall=.TRUE.
2526
2527	IF ( firstcall ) THEN
2528
2529	!Config Key = HYDROL_TAU_HDIFF
2530	!Config Desc = time scale (s) for horizontal diffusion of water
2531	!Config Def = 86400.
2532	!Config If = HYDROL_OK_HDIFF
2533	!Config Help = Defines how fast diffusion occurs horizontally between
2534	!Config the individual PFTs' water reservoirs. If infinite, no
2535	!Config diffusion.
2536
2537	tau_hdiff = 86400.
2538	CALL getin_p('HYDROL_TAU_HDIFF',tau_hdiff)
2539
2540	WRITE (numout,*) 'Hydrol: Horizontal diffusion, tau (s)=',tau_hdiff
2541
2542	firstcall = .FALSE.
2543
2544	ENDIF
2545
2546	! Calculate mean values
2547	! could be done with SUM instruction but this kills vectorization
2548	!
2549	bqsb_mean(:) = zero
2550	gqsb_mean(:) = zero
2551	dss_mean(:) = zero
2552	vegtot(:) = zero
2553	!
2554	DO jv = 1, nvm
2555	DO ji = 1, kjpindex
2556	bqsb_mean(ji) = bqsb_mean(ji) + veget(ji,jv)*bqsb(ji,jv)
2557	gqsb_mean(ji) = gqsb_mean(ji) + veget(ji,jv)*gqsb(ji,jv)
2558	dss_mean(ji) = dss_mean(ji) + veget(ji,jv)*dss(ji,jv)
2559	vegtot(ji) = vegtot(ji) + veget(ji,jv)
2560	ENDDO
2561	ENDDO
2562
2563	DO ji = 1, kjpindex
2564	IF (vegtot(ji) .GT. zero) THEN
2565	bqsb_mean(ji) = bqsb_mean(ji)/vegtot(ji)
2566	gqsb_mean(ji) = gqsb_mean(ji)/vegtot(ji)
2567	dss_mean(ji) = dss_mean(ji)/vegtot(ji)
2568	ENDIF
2569	ENDDO
2570
2571	! relax values towards mean.
2572	!
2573	x = MAX( zero, MIN( dtradia/tau_hdiff, un ) )
2574	!
2575	DO jv = 1, nvm
2576	DO ji = 1, kjpindex
2577	!
2578	bqsb(ji,jv) = (un-x) * bqsb(ji,jv) + x * bqsb_mean(ji)
2579	gqsb(ji,jv) = (un-x) * gqsb(ji,jv) + x * gqsb_mean(ji)
2580	dss(ji,jv) = (un-x) * dss(ji,jv) + x * dss_mean(ji)
2581	!
2582	IF (gqsb(ji,jv) .LT. min_sechiba) THEN
2583	dsg(ji,jv) = zero
2584	ELSE
2585	dsg(ji,jv) = (dss(ji,jv) * ruu_ch(ji) + gqsb(ji,jv)) / ruu_ch(ji)
2586	ENDIF
2587	dsp(ji,jv) = dpu_cste - bqsb(ji,jv) / ruu_ch(ji)
2588	!
2589	ENDDO
2590	ENDDO
2591
2592	END SUBROUTINE hydrolc_hdiff
2593	!
2594	END MODULE hydrolc

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

Download in other formats: