Context Navigation

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

Last change on this file since 257 was 257, checked in by didier.solyga, 13 years ago
Externalized version merged with the trunk
File size: 104.8 KB

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

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

Download in other formats: