Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE/src_sechiba/hydrolc.f90 @ 55

Last change on this file since 55 was 45, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 104.6 KB

Line
1	!!
2	!! This module computes hydrologic processes on continental points.
3	!!
4	!! @author Marie-Alice Foujols and Jan Polcher
5	!! @Version : $Revision$, $Date$
6	!!
7	!< $HeadURL$
8	!< $Date$
9	!< $Author$
10	!< $Revision$
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 constantes_soil
24	USE constantes_veg
25	USE sechiba_io
26	USE grid
27	USE parallel
28	! USE Write_Field_p
29
30	IMPLICIT NONE
31
32	! public routines :
33	! hydrol
34	PRIVATE
35	PUBLIC :: hydrolc_main,hydrolc_clear
36
37	!
38	! variables used inside hydrol module : declaration and initialisation
39	!
40	LOGICAL, SAVE :: l_first_hydrol=.TRUE. !! Initialisation has to be done one time
41	!
42	LOGICAL, SAVE :: check_waterbal=.FALSE. !! The check the water balance
43	LOGICAL, SAVE :: ok_hdiff !! do horizontal diffusion?
44
45	CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O
46
47	! one dimension array allocated, computed, saved and got in hydrol module
48
49	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: bqsb !! Hauteur d'eau dans le reservoir profond
50	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gqsb !! Hauteur d'eau dans le reservoir de surface
51	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsg !! Hauteur du reservoir de surface
52	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dsp !! Hauteur au dessus du reservoir profond
53	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_bqsb !! diagnostique du reservoir profond
54	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mean_gqsb !! diagnostique du reservoir de surface
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_beg !! Total amount of water at start of time step
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_end !! Total amount of water at end of time step
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_beg !! Total amount of water on vegetation at start of time step
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_end !! Total amount of water on vegetation at end of time step
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_beg !! Total amount of water in the soil at start of time step
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_end !! Total amount of water in the soil at end of time step
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_beg !! Total amount of snow at start of time step
62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_end !! Total amount of snow at end of time step
63	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delsoilmoist !! Change in soil moisture
64	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delintercept !! Change in interception storage
65	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delswe !! Change in SWE
66
67	! one dimension array allocated, computed and used in hydrol module exclusively
68
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dss !! Hauteur au dessus du reservoir de surface
70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: hdry !! Dry soil heigth in meters
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Eau tombee sur le sol
72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsnowveg !! Sublimation of snow on vegetation
73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
74	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snowmelt !! Quantite de neige fondue
75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: icemelt !! Quantite de glace fondue
76	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gdrainage !! Drainage between reservoirs
77
78
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vegtot !! Total vegetation
80	! The last vegetation map which was used to distribute the reservoirs
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: resdist !! Distribution of reservoirs
82
83	!! profondeur du reservoir contenant le maximum d'eau
84	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var
85	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: ruu_ch !! Quantite d'eau maximum
86	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: runoff !! Ruissellement
87
88	! Ajout Nathalie - le 28 Mars 2006 - sur conseils Fred Hourdin
89	! Modifs stabilite
90	REAL(r_std), PARAMETER :: dsg_min = 0.001
91
92	CONTAINS
93
94	!!
95	!! Main routine for hydrol module
96	!! - called only one time for initialisation
97	!! - called every time step
98	!! - called one more time at last time step for writing _restart_ file
99	!!
100	!! Algorithm:
101	!! - call hydrolc_snow for snow process (including age of snow)
102	!! - call hydrolc_canop for canopy process
103	!! - call hydrolc_soil for bare soil process
104	!!
105	!! @call hydrolc_snow
106	!! @call hydrolc_canop
107	!! @call hydrolc_soil
108	!!
109	SUBROUTINE hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
110	& temp_sol_new, run_off_tot, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
111	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, tot_melt, transpir, &
112	& precip_rain, precip_snow, returnflow, irrigation, humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr,&
113	& shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id)
114
115	! interface description
116	! input scalar
117	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
118	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
119	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
120	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
121	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
122	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
123	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
124	! input fields
125	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
126	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in) :: indexveg !! Indeces of the points on the 3D map
127	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
128	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
129	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil
130	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation returning to soil moisture
131	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
132	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, ...
133	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+...
134	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
135	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss
136	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
137	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty)
138	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
139	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
140	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation
141	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_corr !! Soil Potential Evaporation Correction
142	! modified fields
143	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapnu !! Bare soil evaporation
144	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapsno !! Snow evaporation
145	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow !! Snow mass [Kg/m^2]
146	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow_age !! Snow age
147	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2]
148	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
149	!! We consider that any water on the ice is snow and we only peforme a water balance to have consistency.
150	!! The water balance is limite to + or - 10^6 so that accumulation is not endless
151	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Relative humidity
152	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vegstress !! Veg. moisture stress (only for vegetation growth)
153	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation due to interception
154	! output fields
155	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: run_off_tot !! Complete runoff
156	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: drainage !! Drainage
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,nbdl), INTENT (inout):: shumdiag !! relative soil moisture
161	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
162	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt
163
164	!
165	! local declaration
166	!
167	REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness
168	REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Depth of snow layer
169
170	INTEGER(i_std) :: ji,jv
171	REAL(r_std), DIMENSION(kjpindex) :: histvar !! computations for history files
172
173	!
174	! do initialisation
175	!
176	IF (l_first_hydrol) THEN
177
178	IF (long_print) WRITE (numout,*) ' l_first_hydrol : call hydrolc_init '
179
180	CALL hydrolc_init (kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, &
181	& snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
182
183	CALL hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag)
184	!
185	! If we check the water balance we first save the total amount of water
186	!
187	IF (check_waterbal) THEN
188	CALL hydrolc_waterbal(kjpindex, index, .TRUE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
189	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,&
190	& run_off_tot, drainage)
191	ENDIF
192	!
193	IF (almaoutput) THEN
194	CALL hydrolc_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwet)
195	ENDIF
196
197	!
198	! shared time step
199	!
200	IF (long_print) WRITE (numout,*) 'hydrolc pas de temps = ',dtradia
201
202	RETURN
203
204	ENDIF
205
206	!
207	! prepares restart file for the next simulation
208	!
209	IF (ldrestart_write) THEN
210
211	IF (long_print) WRITE (numout,*) ' we have to complete restart file with HYDROLOGIC variables '
212
213	var_name= 'humrel'
214	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, humrel, 'scatter', nbp_glo, index_g)
215	!
216	var_name= 'vegstress'
217	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, vegstress, 'scatter', nbp_glo, index_g)
218	!
219	var_name= 'snow'
220	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow, 'scatter', nbp_glo, index_g)
221	!
222	var_name= 'snow_age'
223	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow_age, 'scatter', nbp_glo, index_g)
224	!
225	var_name= 'snow_nobio'
226	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio, 'scatter', nbp_glo, index_g)
227	!
228	var_name= 'snow_nobio_age'
229	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio_age, 'scatter', nbp_glo, index_g)
230	!
231	var_name= 'bqsb'
232	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, bqsb, 'scatter', nbp_glo, index_g)
233	!
234	var_name= 'gqsb'
235	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, gqsb, 'scatter', nbp_glo, index_g)
236	!
237	var_name= 'dsg'
238	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsg, 'scatter', nbp_glo, index_g)
239	!
240	var_name= 'dsp'
241	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dsp, 'scatter', nbp_glo, index_g)
242	!
243	var_name= 'dss'
244	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, dss, 'scatter', nbp_glo, index_g)
245	!
246	var_name= 'hdry'
247	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, hdry, 'scatter', nbp_glo, index_g)
248	!
249	var_name= 'qsintveg'
250	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, qsintveg, 'scatter', nbp_glo, index_g)
251	!
252	var_name= 'resdist'
253	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, resdist, 'scatter', nbp_glo, index_g)
254	RETURN
255	!
256	END IF
257
258	!
259	! computes snow
260	!
261
262	CALL hydrolc_snow(kjpindex, dtradia, precip_rain, precip_snow, temp_sol_new, soilcap, &
263	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
264	& tot_melt, snowdepth)
265	!
266	! computes canopy
267	!
268	!
269	CALL hydrolc_vegupd(kjpindex, veget, ruu_ch, qsintveg, gqsb, bqsb, dsg, dss,dsp, resdist)
270	!
271	CALL hydrolc_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol)
272	!
273	! computes hydro_soil
274	!
275	CALL hydrolc_soil(kjpindex, vevapnu, precisol, returnflow, irrigation, tot_melt, mx_eau_var, veget, ruu_ch, transpir,&
276	& gqsb, bqsb, dsg, dss, rsol, drysoil_frac, hdry, dsp, runoff, run_off_tot, drainage, humrel, vegstress, &
277	& shumdiag, litterhumdiag)
278	!
279	! computes horizontal diffusion between the water reservoirs
280	!
281	IF ( ok_hdiff ) THEN
282	CALL hydrolc_hdiff(kjpindex, dtradia, veget, ruu_ch, gqsb, bqsb, dsg, dss, dsp)
283	ENDIF
284	!
285	! If we check the water balance we end with the comparison of total water change and fluxes
286	!
287	IF (check_waterbal) THEN
288	CALL hydrolc_waterbal(kjpindex, index, .FALSE., dtradia, veget, totfrac_nobio, qsintveg, snow, snow_nobio,&
289	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno,&
290	& run_off_tot, drainage )
291	ENDIF
292	!
293	! If we use the ALMA standards
294	!
295	IF (almaoutput) THEN
296	CALL hydrolc_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet)
297	ENDIF
298
299	!
300	IF ( .NOT. almaoutput ) THEN
301	CALL histwrite(hist_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
302	CALL histwrite(hist_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
303	CALL histwrite(hist_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
304	CALL histwrite(hist_id, 'runoff', kjit, run_off_tot, kjpindex, index)
305	CALL histwrite(hist_id, 'drainage', kjit, drainage, kjpindex, index)
306	CALL histwrite(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
307	CALL histwrite(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
308	CALL histwrite(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
309	CALL histwrite(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
310	CALL histwrite(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
311	CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
312
313	histvar(:)=zero
314	DO jv = 1, nvm
315	DO ji = 1, kjpindex
316	IF ( vegtot(ji) .GT. zero ) THEN
317	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))mx_eau_eau, zero)
318	ENDIF
319	ENDDO
320	ENDDO
321	CALL histwrite(hist_id, 'mrsos', kjit, histvar, kjpindex, index)
322
323	histvar(:)=mean_bqsb(:)+mean_gqsb(:)
324	CALL histwrite(hist_id, 'mrso', kjit, histvar, kjpindex, index)
325
326	histvar(:)=run_off_tot(:)/one_day
327	CALL histwrite(hist_id, 'mrros', kjit, histvar, kjpindex, index)
328
329	histvar(:)=(run_off_tot(:)+drainage(:))/one_day
330	CALL histwrite(hist_id, 'mrro', kjit, histvar, kjpindex, index)
331
332	histvar(:)=(precip_rain(:)-SUM(precisol(:,:),dim=2))/one_day
333	CALL histwrite(hist_id, 'prveg', kjit, histvar, kjpindex, index)
334
335	ELSE
336	CALL histwrite(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
337	CALL histwrite(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
338	CALL histwrite(hist_id, 'Qs', kjit, run_off_tot, kjpindex, index)
339	CALL histwrite(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
340	CALL histwrite(hist_id, 'Qsm', kjit, tot_melt, kjpindex, index)
341	CALL histwrite(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
342	CALL histwrite(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
343	CALL histwrite(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
344	!
345	CALL histwrite(hist_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
346	CALL histwrite(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
347	CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
348	!
349	CALL histwrite(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
350	CALL histwrite(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
351	!
352	CALL histwrite(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
353	!
354	ENDIF
355	IF ( hist2_id > 0 ) THEN
356	IF ( .NOT. almaoutput ) THEN
357	CALL histwrite(hist2_id, 'dss', kjit, dss, kjpindex*nvm, indexveg)
358	CALL histwrite(hist2_id, 'bqsb', kjit, mean_bqsb, kjpindex, index)
359	CALL histwrite(hist2_id, 'gqsb', kjit, mean_gqsb, kjpindex, index)
360	CALL histwrite(hist2_id, 'runoff', kjit, run_off_tot, kjpindex, index)
361	CALL histwrite(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
362	CALL histwrite(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
363	CALL histwrite(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
364	CALL histwrite(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
365	CALL histwrite(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
366	CALL histwrite(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
367	CALL histwrite(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
368
369	histvar(:)=zero
370	DO jv = 1, nvm
371	DO ji = 1, kjpindex
372	IF ( vegtot(ji) .GT. zero ) THEN
373	histvar(ji)=histvar(ji)+veget(ji,jv)/vegtot(ji)MAX((0.1-dss(ji,jv))mx_eau_eau, zero)
374	ENDIF
375	ENDDO
376	ENDDO
377	CALL histwrite(hist2_id, 'mrsos', kjit, histvar, kjpindex, index)
378
379	histvar(:)=(run_off_tot(:)+drainage(:))/one_day
380	CALL histwrite(hist2_id, 'mrro', kjit, histvar, kjpindex, index)
381
382	ELSE
383	CALL histwrite(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
384	CALL histwrite(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
385	CALL histwrite(hist2_id, 'Qs', kjit, run_off_tot, kjpindex, index)
386	CALL histwrite(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
387	CALL histwrite(hist2_id, 'Qsm', kjit, tot_melt, kjpindex, index)
388	CALL histwrite(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
389	CALL histwrite(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
390	CALL histwrite(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
391	!
392	CALL histwrite(hist2_id, 'SoilMoist', kjit, tot_watsoil_end, kjpindex, index)
393	CALL histwrite(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
394	!
395	CALL histwrite(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
396	CALL histwrite(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
397	!
398	CALL histwrite(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
399	!
400	ENDIF
401	ENDIF
402
403	!
404	IF (long_print) WRITE (numout,*) ' hydrolc_main Done '
405
406	END SUBROUTINE hydrolc_main
407
408	!! Algorithm:
409	!! - dynamic allocation for local array
410	!! - _restart_ file reading for HYDROLOGIC variables
411	!!
412	SUBROUTINE hydrolc_init(kjit, ldrestart_read, kjpindex, index, rest_id, veget, humrel, vegstress, &
413	& snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
414
415	! interface description
416	! input scalar
417	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
418	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read
419	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
420	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
421	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
422	! input fields
423	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Carte de vegetation
424	! output fields
425	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Stress hydrique, relative humidity
426	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
427	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow mass [Kg/m^2]
428	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age
429	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio !! Snow on ice, lakes, ...
430	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age !! Snow age on ice, lakes, ...
431	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception
432
433	! local declaration
434	INTEGER(i_std) :: ier
435	INTEGER(i_std) :: ji,jv,ik
436	REAL(r_std), DIMENSION (kjpindex,nvm) :: zdsp, tmpdss
437
438	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: dsp_g
439	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: zdsp_g
440
441	REAL(r_std), DIMENSION(kjpindex) :: a_subgrd
442
443	! initialisation
444	IF (l_first_hydrol) THEN
445	l_first_hydrol=.FALSE.
446	ELSE
447	WRITE (numout,*) ' l_first_hydrol false . we stop '
448	STOP 'hydrolc_init'
449	ENDIF
450
451	!Config Key = HYDROL_OK_HDIFF
452	!Config Desc = do horizontal diffusion?
453	!Config Def = n
454	!Config Help = If TRUE, then water can diffuse horizontally between
455	!Config the PFTs' water reservoirs.
456
457	ok_hdiff = .FALSE.
458	CALL getin_p('HYDROL_OK_HDIFF',ok_hdiff)
459
460	! make dynamic allocation with good dimension
461
462	! one dimension array allocation with possible restart value
463
464	ALLOCATE (bqsb(kjpindex,nvm),stat=ier)
465	IF (ier.NE.0) THEN
466	WRITE (numout,*) ' error in bqsb allocation. We stop. We need kjpindex words = ',kjpindex
467	STOP 'hydrolc_init'
468	END IF
469	bqsb(:,:) = zero
470
471	ALLOCATE (gqsb(kjpindex,nvm),stat=ier)
472	IF (ier.NE.0) THEN
473	WRITE (numout,*) ' error in gqsb allocation. We stop. We need kjpindex words = ',kjpindex
474	STOP 'hydrolc_init'
475	END IF
476	gqsb(:,:) = zero
477
478	ALLOCATE (dsg(kjpindex,nvm),stat=ier)
479	IF (ier.NE.0) THEN
480	WRITE (numout,*) ' error in dsg allocation. We stop. We need kjpindex words = ',kjpindex
481	STOP 'hydrolc_init'
482	END IF
483	dsg(:,:) = zero
484
485	ALLOCATE (dsp(kjpindex,nvm),stat=ier)
486	IF (ier.NE.0) THEN
487	WRITE (numout,*) ' error in dsp allocation. We stop. We need kjpindex words = ',kjpindex
488	STOP 'hydrolc_init'
489	END IF
490	dsp(:,:) = zero
491
492	! one dimension array allocation
493
494	ALLOCATE (mean_bqsb(kjpindex),stat=ier)
495	IF (ier.NE.0) THEN
496	WRITE (numout,*) ' error in mean_bqsb allocation. We stop. We need kjpindex words = ',kjpindex
497	STOP 'hydrolc_init'
498	END IF
499	mean_bqsb(:) = zero
500
501	ALLOCATE (mean_gqsb(kjpindex),stat=ier)
502	IF (ier.NE.0) THEN
503	WRITE (numout,*) ' error in mean_gqsb allocation. We stop. We need kjpindex words = ',kjpindex
504	STOP 'hydrolc_init'
505	END IF
506	mean_gqsb(:) = zero
507
508	ALLOCATE (dss(kjpindex,nvm),stat=ier)
509	IF (ier.NE.0) THEN
510	WRITE (numout,*) ' error in dss allocation. We stop. We need kjpindex words = ',kjpindex
511	STOP 'hydrolc_init'
512	END IF
513	dss(:,:) = zero
514
515	ALLOCATE (hdry(kjpindex),stat=ier)
516	IF (ier.NE.0) THEN
517	WRITE (numout,*) ' error in hdry allocation. We stop. We need kjpindex words = ',kjpindex
518	STOP 'hydrolc_init'
519	END IF
520	hdry(:) = zero
521
522	ALLOCATE (precisol(kjpindex,nvm),stat=ier)
523	IF (ier.NE.0) THEN
524	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
525	STOP 'hydrolc_init'
526	END IF
527	precisol(:,:) = zero
528
529	ALLOCATE (gdrainage(kjpindex,nvm),stat=ier)
530	IF (ier.NE.0) THEN
531	WRITE (numout,*) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindex
532	STOP 'hydrolc_init'
533	END IF
534	gdrainage(:,:) = zero
535
536	ALLOCATE (subsnowveg(kjpindex),stat=ier)
537	IF (ier.NE.0) THEN
538	WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex
539	STOP 'hydrolc_init'
540	END IF
541	subsnowveg(:) = zero
542
543	ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier)
544	IF (ier.NE.0) THEN
545	WRITE (numout,) ' error in subsnownobio allocation. We stop. We need kjpindexnnobio words = ', &
546	kjpindex*nnobio
547	STOP 'hydrolc_init'
548	END IF
549	subsnownobio(:,:) = zero
550
551	ALLOCATE (snowmelt(kjpindex),stat=ier)
552	IF (ier.NE.0) THEN
553	WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex
554	STOP 'hydrolc_init'
555	END IF
556	snowmelt(:) = zero
557
558	ALLOCATE (icemelt(kjpindex),stat=ier)
559	IF (ier.NE.0) THEN
560	WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex
561	STOP 'hydrolc_init'
562	END IF
563	icemelt(:) = zero
564
565	ALLOCATE (mx_eau_var(kjpindex),stat=ier)
566	IF (ier.NE.0) THEN
567	WRITE (numout,*) ' error in mx_eau_var allocation. We stop. We need kjpindex words = ',kjpindex
568	STOP 'hydrolc_init'
569	END IF
570	mx_eau_var(:) = zero
571
572	ALLOCATE (ruu_ch(kjpindex),stat=ier)
573	IF (ier.NE.0) THEN
574	WRITE (numout,*) ' error in ruu_ch allocation. We stop. We need kjpindex words = ',kjpindex
575	STOP 'hydrolc_init'
576	END IF
577	ruu_ch(:) = zero
578
579	ALLOCATE (vegtot(kjpindex),stat=ier)
580	IF (ier.NE.0) THEN
581	WRITE (numout,) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindexnvm
582	STOP 'hydrolc_init'
583	END IF
584	vegtot(:) = zero
585
586	ALLOCATE (resdist(kjpindex,nvm),stat=ier)
587	IF (ier.NE.0) THEN
588	WRITE (numout,) ' error in resdist allocation. We stop. We need kjpindex words = ',kjpindexnvm
589	STOP 'hydrolc_init'
590	END IF
591	resdist(:,:) = zero
592
593	ALLOCATE (runoff(kjpindex,nvm),stat=ier)
594	IF (ier.NE.0) THEN
595	WRITE (numout,) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindexnvm
596	STOP 'hydrolc_init'
597	END IF
598	runoff(:,:) = zero
599	!
600	! If we check the water balance we need two more variables
601	!
602	IF ( check_waterbal ) THEN
603
604	ALLOCATE (tot_water_beg(kjpindex),stat=ier)
605	IF (ier.NE.0) THEN
606	WRITE (numout,*) ' error in tot_water_beg allocation. We stop. We need kjpindex words = ',kjpindex
607	STOP 'hydrolc_init'
608	END IF
609
610	ALLOCATE (tot_water_end(kjpindex),stat=ier)
611	IF (ier.NE.0) THEN
612	WRITE (numout,*) ' error in tot_water_end allocation. We stop. We need kjpindex words = ',kjpindex
613	STOP 'hydrolc_init'
614	END IF
615
616	ENDIF
617	!
618	! If we use the almaoutputs we need four more variables
619	!
620	IF ( almaoutput ) THEN
621
622	ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
623	IF (ier.NE.0) THEN
624	WRITE (numout,*) ' error in tot_watveg_beg allocation. We stop. We need kjpindex words = ',kjpindex
625	STOP 'hydrolc_init'
626	END IF
627
628	ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
629	IF (ier.NE.0) THEN
630	WRITE (numout,*) ' error in tot_watveg_end allocation. We stop. We need kjpindex words = ',kjpindex
631	STOP 'hydrolc_init'
632	END IF
633
634	ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
635	IF (ier.NE.0) THEN
636	WRITE (numout,*) ' error in tot_watsoil_beg allocation. We stop. We need kjpindex words = ',kjpindex
637	STOP 'hydrolc_init'
638	END IF
639
640	ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
641	IF (ier.NE.0) THEN
642	WRITE (numout,*) ' error in tot_watsoil_end allocation. We stop. We need kjpindex words = ',kjpindex
643	STOP 'hydrolc_init'
644	END IF
645
646	ALLOCATE (delsoilmoist(kjpindex),stat=ier)
647	IF (ier.NE.0) THEN
648	WRITE (numout,*) ' error in delsoilmoist allocation. We stop. We need kjpindex words = ',kjpindex
649	STOP 'hydrolc_init'
650	END IF
651
652	ALLOCATE (delintercept(kjpindex),stat=ier)
653	IF (ier.NE.0) THEN
654	WRITE (numout,*) ' error in delintercept. We stop. We need kjpindex words = ',kjpindex
655	STOP 'hydrolc_init'
656	END IF
657
658	ALLOCATE (delswe(kjpindex),stat=ier)
659	IF (ier.NE.0) THEN
660	WRITE (numout,*) ' error in delswe. We stop. We need kjpindex words = ',kjpindex
661	STOP 'hydrolc_init'
662	ENDIF
663
664	ALLOCATE (snow_beg(kjpindex),stat=ier)
665	IF (ier.NE.0) THEN
666	WRITE (numout,*) ' error in snow_beg allocation. We stop. We need kjpindex words =',kjpindex
667	STOP 'hydrolc_init'
668	END IF
669
670	ALLOCATE (snow_end(kjpindex),stat=ier)
671	IF (ier.NE.0) THEN
672	WRITE (numout,*) ' error in snow_end allocation. We stop. We need kjpindex words =',kjpindex
673	STOP 'hydrolc_init'
674	END IF
675
676	ENDIF
677
678	! open restart input file done by sechiba_init
679	! and read data from restart input file for HYDROLOGIC process
680
681	IF (ldrestart_read) THEN
682
683	IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
684
685	var_name= 'snow'
686	CALL ioconf_setatt('UNITS', 'kg/m^2')
687	CALL ioconf_setatt('LONG_NAME','Snow mass')
688	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
689	!!$ ! correction for old restart
690	!!$ DO ik=1, kjpindex
691	!!$ if(snow(ik).gt.maxmass_glacier) snow(ik)=maxmass_glacier
692	!!$ ENDDO
693	!
694	var_name= 'snow_age'
695	CALL ioconf_setatt('UNITS', 'd')
696	CALL ioconf_setatt('LONG_NAME','Snow age')
697	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
698	!
699	var_name= 'snow_nobio'
700	CALL ioconf_setatt('UNITS', 'kg/m^2')
701	CALL ioconf_setatt('LONG_NAME','Snow on other surface types')
702	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
703	!!$ ! correction for old restart
704	!!$ DO ik=1, kjpindex
705	!!$ if(snow_nobio(ik,iice).gt.maxmass_glacier) snow_nobio(ik,iice)=maxmass_glacier
706	!!$ ENDDO
707	!
708	var_name= 'snow_nobio_age'
709	CALL ioconf_setatt('UNITS', 'd')
710	CALL ioconf_setatt('LONG_NAME','Snow age on other surface types')
711	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
712	!
713	var_name= 'humrel'
714	CALL ioconf_setatt('UNITS', '-')
715	CALL ioconf_setatt('LONG_NAME','Soil moisture stress')
716	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., humrel, "gather", nbp_glo, index_g)
717	!
718	var_name= 'vegstress'
719	CALL ioconf_setatt('UNITS', '-')
720	CALL ioconf_setatt('LONG_NAME','Vegetation growth moisture stress')
721	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
722	!
723	var_name= 'bqsb'
724	CALL ioconf_setatt('UNITS', 'kg/m^2')
725	CALL ioconf_setatt('LONG_NAME','Deep soil moisture')
726	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., bqsb, "gather", nbp_glo, index_g)
727	!
728	var_name= 'gqsb'
729	CALL ioconf_setatt('UNITS', 'kg/m^2')
730	CALL ioconf_setatt('LONG_NAME','Surface soil moisture')
731	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., gqsb, "gather", nbp_glo, index_g)
732	!
733	var_name= 'dsg'
734	CALL ioconf_setatt('UNITS', 'm')
735	CALL ioconf_setatt('LONG_NAME','Depth of upper reservoir')
736	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsg, "gather", nbp_glo, index_g)
737	!
738	var_name= 'dsp'
739	CALL ioconf_setatt('UNITS', 'm')
740	CALL ioconf_setatt('LONG_NAME','Depth to lower reservoir')
741	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dsp, "gather", nbp_glo, index_g)
742	!
743	var_name= 'dss'
744	CALL ioconf_setatt('UNITS', 'm')
745	CALL ioconf_setatt('LONG_NAME','Hauteur au dessus du reservoir de surface')
746	IF ( ok_var(var_name) ) THEN
747	CALL restget_p (rest_id, var_name, nbp_glo, nvm , 1, kjit, .TRUE., dss, "gather", nbp_glo, index_g)
748	ENDIF
749	!
750	var_name= 'hdry'
751	CALL ioconf_setatt('UNITS', 'm')
752	CALL ioconf_setatt('LONG_NAME','Dry soil heigth in meters')
753	IF ( ok_var(var_name) ) THEN
754	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., hdry, "gather", nbp_glo, index_g)
755	ENDIF
756	!
757	var_name= 'qsintveg'
758	CALL ioconf_setatt('UNITS', 'kg/m^2')
759	CALL ioconf_setatt('LONG_NAME','Intercepted moisture')
760	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
761	!
762	var_name= 'resdist'
763	CALL ioconf_setatt('UNITS', '-')
764	CALL ioconf_setatt('LONG_NAME','Distribution of reservoirs')
765	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
766	!
767	! get restart values if non were found in the restart file
768	!
769	!Config Key = HYDROL_SNOW
770	!Config Desc = Initial snow mass if not found in restart
771	!Config Def = 0.0
772	!Config Help = The initial value of snow mass if its value is not found
773	!Config in the restart file. This should only be used if the model is
774	!Config started without a restart file.
775	!
776	CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', zero)
777	!
778	!Config Key = HYDROL_SNOWAGE
779	!Config Desc = Initial snow age if not found in restart
780	!Config Def = 0.0
781	!Config Help = The initial value of snow age if its value is not found
782	!Config in the restart file. This should only be used if the model is
783	!Config started without a restart file.
784	!
785	CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', zero)
786	!
787	!Config Key = HYDROL_SNOW_NOBIO
788	!Config Desc = Initial snow amount on ice, lakes, etc. if not found in restart
789	!Config Def = 0.0
790	!Config Help = The initial value of snow if its value is not found
791	!Config in the restart file. This should only be used if the model is
792	!Config started without a restart file.
793	!
794	CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero)
795	!
796	!Config Key = HYDROL_SNOW_NOBIO_AGE
797	!Config Desc = Initial snow age on ice, lakes, etc. if not found in restart
798	!Config Def = 0.0
799	!Config Help = The initial value of snow age if its value is not found
800	!Config in the restart file. This should only be used if the model is
801	!Config started without a restart file.
802	!
803	CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero)
804	!
805	!Config Key = HYDROL_HUMR
806	!Config Desc = Initial soil moisture stress if not found in restart
807	!Config Def = 1.0
808	!Config Help = The initial value of soil moisture stress if its value is not found
809	!Config in the restart file. This should only be used if the model is
810	!Config started without a restart file.
811	!
812	CALL setvar_p (humrel, val_exp,'HYDROL_HUMR', un)
813	CALL setvar_p (vegstress, val_exp,'HYDROL_HUMR', un)
814	!
815	!Config Key = HYDROL_BQSB
816	!Config Desc = Initial restart deep soil moisture if not found in restart
817	!Config Def = DEF
818	!Config Help = The initial value of deep soil moisture if its value is not found
819	!Config in the restart file. This should only be used if the model is
820	!Config started without a restart file. Default behaviour is a saturated soil.
821	!
822	CALL setvar_p (bqsb, val_exp, 'HYDROL_BQSB', mx_eau_eau*dpu_cste)
823	!
824	!Config Key = HYDROL_GQSB
825	!Config Desc = Initial upper soil moisture if not found in restart
826	!Config Def = 0.0
827	!Config Help = The initial value of upper soil moisture if its value is not found
828	!Config in the restart file. This should only be used if the model is
829	!Config started without a restart file.
830	!
831	CALL setvar_p (gqsb, val_exp, 'HYDROL_GQSB', zero)
832	!
833	!Config Key = HYDROL_DSG
834	!Config Desc = Initial upper reservoir depth if not found in restart
835	!Config Def = 0.0
836	!Config Help = The initial value of upper reservoir depth if its value is not found
837	!Config in the restart file. This should only be used if the model is
838	!Config started without a restart file.
839	!
840	CALL setvar_p (dsg, val_exp, 'HYDROL_DSG', zero)
841
842	! set inital value for dsp if needed
843	!
844	!Config Key = HYDROL_DSP
845	!Config Desc = Initial dry soil above upper reservoir if not found in restart
846	!Config Def = DEF
847	!Config Help = The initial value of dry soil above upper reservoir if its value
848	!Config is not found in the restart file. This should only be used if
849	!Config the model is started without a restart file. The default behaviour
850	!Config is to compute it from the variables above. Should be OK most of
851	!Config the time.
852	!
853	zdsp(:,:) = dpu_cste - bqsb(:,:) / mx_eau_eau
854	dsp(1,1) = val_exp
855	call getin_p('HYDROL_DSP', dsp(1,1))
856	IF (dsp(1,1) == val_exp) THEN
857	dsp(:,:) = zdsp(:,:)
858	ELSE
859	IF (is_root_prc) &
860	ALLOCATE(zdsp_g(nbp_glo,nvm),dsp_g(nbp_glo,nvm))
861	CALL gather(zdsp,zdsp_g)
862	IF (is_root_prc) &
863	CALL setvar (dsp_g, val_exp, 'HYDROL_DSP', zdsp_g)
864	CALL scatter(dsp_g,dsp)
865	IF (is_root_prc) &
866	DEALLOCATE(zdsp_g, dsp_g)
867	ENDIF
868	!
869	!Config Key = HYDROL_QSV
870	!Config Desc = Initial water on canopy if not found in restart
871	!Config Def = 0.0
872	!Config Help = The initial value of moisture on canopy if its value
873	!Config is not found in the restart file. This should only be used if
874	!Config the model is started without a restart file.
875	!
876	CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', zero)
877	!
878	tmpdss = dsg - gqsb / mx_eau_eau
879	IF ( ok_var("dss") ) THEN
880	CALL setvar_p (dss, val_exp, 'NO_KEYWORD', tmpdss)
881	ELSE
882	dss(:,:)=tmpdss(:,:)
883	ENDIF
884
885	IF ( ok_var("hdry") ) THEN
886	IF (MINVAL(hdry) .EQ. MAXVAL(hdry) .AND. MAXVAL(hdry) .EQ. val_exp) THEN
887	a_subgrd(:) = zero
888	DO ji=1,kjpindex
889	IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN
890	!
891	IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN
892	! Ajouts Nathalie - Fred - le 28 Mars 2006
893	a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),zero)/dsg_min,un)
894	!
895	ENDIF
896	ENDIF
897	!
898	ENDDO
899	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
900	! revu 22 novembre 2007
901	hdry(:) = a_subgrd(:)dss(:,1) + (1.-a_subgrd(:))dsp(:,1)
902	ENDIF
903	ELSE
904	a_subgrd(:) = zero
905	DO ji=1,kjpindex
906	IF ( gqsb(ji,1)+bqsb(ji,1) .GT. zero ) THEN
907	!
908	IF (.NOT. (dsg(ji,1).EQ. zero .OR. gqsb(ji,1).EQ.zero)) THEN
909	! Ajouts Nathalie - Fred - le 28 Mars 2006
910	a_subgrd(ji)=MIN(MAX(dsg(ji,1)-dss(ji,1),zero)/dsg_min,un)
911	!
912	ENDIF
913	ENDIF
914	!
915	ENDDO
916
917	! Correction Nathalie - le 28 Mars 2006 - re-ecriture drysoil_frac/hdry - Fred Hourdin
918	! revu 22 novembre 2007
919	hdry(:) = a_subgrd(:)dss(:,1) + (un-a_subgrd(:))dsp(:,1)
920	ENDIF
921	!
922	! There is no need to configure the initialisation of resdist. If not available it is the vegetation map
923	!
924	IF ( MINVAL(resdist) .EQ. MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
925	resdist = veget
926	ENDIF
927	!
928	! Remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1. Thus we need vegtot
929	!
930	DO ji = 1, kjpindex
931	vegtot(ji) = SUM(veget(ji,:))
932	ENDDO
933	!
934	ENDIF
935
936	!
937	! Where vegetation fraction is zero, set water to that of bare soil.
938	! This does not create any additional water.
939	!
940	DO jv = 2, nvm
941	DO ji = 1, kjpindex
942	IF ( veget(ji,jv) .LT. EPSILON(un) ) THEN
943	gqsb(ji,jv) = gqsb(ji,1)
944	bqsb(ji,jv) = bqsb(ji,1)
945	dsg(ji,jv) = dsg(ji,1)
946	dss(ji,jv) = dss(ji,1)
947	dsp(ji,jv) = dsp(ji,1)
948	ENDIF
949	ENDDO
950	ENDDO
951	!
952	DO ik=1, kjpindex
953	if(snow(ik).gt.maxmass_glacier) then
954	WRITE(numout,*)' il faut diminuer le stock de neige car snow > maxmass_glacier dans restart'
955	snow(ik)=maxmass_glacier
956	endif
957	if(snow_nobio(ik,iice).gt.maxmass_glacier) then
958	WRITE(numout,*)' il faut diminuer le stock de neige car snow_nobio > maxmass_glacier dans restart'
959	snow_nobio(ik,iice)=maxmass_glacier
960	endif
961	ENDDO
962	!
963	IF (long_print) WRITE (numout,*) ' hydrolc_init done '
964	!
965	END SUBROUTINE hydrolc_init
966	!
967	!-------------------------------------
968	!
969	SUBROUTINE hydrolc_clear()
970
971	l_first_hydrol=.TRUE.
972
973	IF (ALLOCATED (bqsb)) DEALLOCATE (bqsb)
974	IF (ALLOCATED (gqsb)) DEALLOCATE (gqsb)
975	IF (ALLOCATED (dsg)) DEALLOCATE (dsg)
976	IF (ALLOCATED (dsp)) DEALLOCATE (dsp)
977	IF (ALLOCATED (mean_bqsb)) DEALLOCATE (mean_bqsb)
978	IF (ALLOCATED (mean_gqsb)) DEALLOCATE (mean_gqsb)
979	IF (ALLOCATED (dss)) DEALLOCATE (dss)
980	IF (ALLOCATED (hdry)) DEALLOCATE (hdry)
981	IF (ALLOCATED (precisol)) DEALLOCATE (precisol)
982	IF (ALLOCATED (gdrainage)) DEALLOCATE (gdrainage)
983	IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg)
984	IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio)
985	IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt)
986	IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt)
987	IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var)
988	IF (ALLOCATED (ruu_ch)) DEALLOCATE (ruu_ch)
989	IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot)
990	IF (ALLOCATED (resdist)) DEALLOCATE (resdist)
991	IF (ALLOCATED (runoff)) DEALLOCATE (runoff)
992	IF (ALLOCATED (tot_water_beg)) DEALLOCATE (tot_water_beg)
993	IF (ALLOCATED (tot_water_end)) DEALLOCATE (tot_water_end)
994	IF (ALLOCATED (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
995	IF (ALLOCATED (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
996	IF (ALLOCATED (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
997	IF (ALLOCATED (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
998	IF (ALLOCATED (delsoilmoist)) DEALLOCATE (delsoilmoist)
999	IF (ALLOCATED (delintercept)) DEALLOCATE (delintercept)
1000	IF (ALLOCATED (snow_beg)) DEALLOCATE (snow_beg)
1001	IF (ALLOCATED (snow_end)) DEALLOCATE (snow_end)
1002	IF (ALLOCATED (delswe)) DEALLOCATE (delswe)
1003	!
1004	END SUBROUTINE hydrolc_clear
1005
1006	!! This routine initializes HYDROLOGIC variables
1007	!! - mx_eau_var
1008	!! - ruu_ch
1009	!!
1010	SUBROUTINE hydrolc_var_init (kjpindex, veget, veget_max, rsol, drysoil_frac, mx_eau_var, ruu_ch, shumdiag, litterhumdiag)
1011
1012	! interface description
1013	! input scalar
1014	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1015	! input fields
1016	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
1017	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type
1018	! output fields
1019	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rsol !! Resistance to bare soil evaporation
1020	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Fraction of visible dry soil
1021	!! Profondeur du reservoir contenant le maximum d'eau
1022	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: mx_eau_var !!
1023	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: ruu_ch !! Quantite d'eau maximum
1024	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out):: shumdiag !! relative soil moisture
1025	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
1026	! local declaration
1027	INTEGER(i_std) :: ji,jv, jd
1028	REAL(r_std), DIMENSION(kjpindex) :: mean_dsg
1029	REAL(r_std) :: gtr, btr
1030	REAL(r_std), DIMENSION(nbdl+1) :: tmp_dl
1031	!
1032	! initialisation
1033	tmp_dl(1) = 0
1034	tmp_dl(2:nbdl+1) = diaglev(1:nbdl)
1035	!
1036	mx_eau_var(:) = zero
1037	!
1038	DO ji = 1,kjpindex
1039	DO jv = 1,nvm
1040	mx_eau_var(ji) = mx_eau_var(ji) + veget(ji,jv)wmax_veg(jv)dpu_cste
1041	END DO
1042	IF (vegtot(ji) .GT. zero) THEN
1043	mx_eau_var(ji) = mx_eau_var(ji)/vegtot(ji)
1044	ELSE
1045	mx_eau_var(ji) = mx_eau_eau*dpu_cste
1046	ENDIF
1047	ruu_ch(ji) = mx_eau_var(ji) / dpu_cste
1048	END DO
1049	!
1050	!
1051	! could be done with SUM instruction but this kills vectorization
1052	mean_bqsb(:) = zero
1053	mean_gqsb(:) = zero
1054	mean_dsg(:) = zero
1055	DO jv = 1, nvm
1056	DO ji = 1, kjpindex
1057	mean_bqsb(ji) = mean_bqsb(ji) + resdist(ji,jv)*bqsb(ji,jv)
1058	mean_gqsb(ji) = mean_gqsb(ji) + resdist(ji,jv)*gqsb(ji,jv)
1059	mean_dsg(ji) = mean_dsg(ji) + resdist(ji,jv)*dsg(ji,jv)
1060	ENDDO
1061	ENDDO
1062	mean_dsg(:) = MAX( mean_dsg(:), mean_gqsb(:)/ruu_ch(:) )
1063
1064	DO ji = 1, kjpindex
1065	IF (vegtot(ji) .GT. zero) THEN
1066	mean_bqsb(ji) = mean_bqsb(ji)/vegtot(ji)
1067	mean_gqsb(ji) = mean_gqsb(ji)/vegtot(ji)
1068	mean_dsg(ji) = mean_dsg(ji)/vegtot(ji)
1069	ENDIF
1070	ENDDO
1071
1072	DO jd = 1,nbdl
1073	!
1074	DO ji = 1,kjpindex
1075	!!$ !
1076	!!$ DO jd = 1,nbdl
1077	IF ( tmp_dl(jd+1) .LT. mean_dsg(ji)) THEN
1078	shumdiag(ji,jd) = mean_gqsb(ji)/mx_eau_var(ji)
1079	ELSE
1080	IF ( tmp_dl(jd) .LT. mean_dsg(ji)) THEN
1081	gtr = (mean_dsg(ji)-tmp_dl(jd))/(tmp_dl(jd+1)-tmp_dl(jd))
1082	btr = 1 - gtr
1083	shumdiag(ji,jd) = gtr*mean_gqsb(ji)/mx_eau_var(ji) + &
1084	& btr*mean_bqsb(ji)/mx_eau_var(ji)
1085	ELSE
1086	shumdiag(ji,jd) = mean_bqsb(ji)/mx_eau_var(ji)
1087	ENDIF
1088	ENDIF
1089	shumdiag(ji,jd) = MAX(MIN(shumdiag(ji,jd), un), zero)
1090	ENDDO
1091	ENDDO
1092
1093	! The fraction of soil which is visibly dry (dry when dss = 0.1 m)
1094	drysoil_frac(:) = MIN(MAX(dss(:,1),zero)*10._r_std, un)
1095	!
1096	! Compute the resistance to bare soil evaporation
1097	!
1098	rsol(:) = -un
1099	DO ji = 1, kjpindex
1100	IF (veget(ji,1) .GE. min_sechiba) THEN
1101	!
1102	! Correction Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
1103	! on modifie le rsol pour que la resistance croisse subitement si on s'approche
1104	! du fond. En gros, rsol=hdry*rsol_cste pour hdry < 1m70
1105	!rsol(ji) = dss(ji,1) * rsol_cste
1106	!rsol(ji) = ( drysoil_frac(ji) + un/(10.(dpu_cste - drysoil_frac(ji))+1.e-10)2 ) rsol_cste
1107	rsol(ji) = ( hdry(ji) + un/(10.(dpu_cste - hdry(ji))+1.e-10)2 ) rsol_cste
1108	ENDIF
1109	ENDDO
1110
1111	!
1112	! litter humidity.
1113	!
1114
1115	!!$ DO ji = 1, kjpindex
1116	!!$ litterhumdiag(ji) = EXP( - dss(ji,1) / hcrit_litter )
1117	!!$ ENDDO
1118	litterhumdiag(:) = EXP( - hdry(:) / hcrit_litter )
1119
1120	! special case: it has just been raining a few drops. The upper soil
1121	! reservoir is ridiculously small, but the dry soil height is zero.
1122	! Don't take it into account.
1123
1124	!!$ DO ji = 1, kjpindex
1125	!!$ IF ( ( dss(ji,1) .LT. min_sechiba ) .AND. &
1126	!!$ ( mean_dsg(ji) .GT. min_sechiba ) .AND. &
1127	!!$ ( mean_dsg(ji) .LT. 5.E-4 ) ) THEN
1128	!!$ litterhumdiag(ji) = zero
1129	!!$ ENDIF
1130	!!$ ENDDO
1131	WHERE ( ( hdry(:) .LT. min_sechiba ) .AND. &
1132	( mean_dsg(:) .GT. min_sechiba ) .AND. ( mean_dsg(:) .LT. 5.E-4 ) )
1133	litterhumdiag(:) = zero
1134	ENDWHERE
1135
1136	!
1137	IF (long_print) WRITE (numout,*) ' hydrolc_var_init done '
1138
1139	END SUBROUTINE hydrolc_var_init
1140
1141	!! This routine computes snow processes
1142	!!
1143	SUBROUTINE hydrolc_snow (kjpindex, dtradia, precip_rain, precip_snow , temp_sol_new, soilcap,&
1144	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
1145	& tot_melt, snowdepth)
1146
1147	!
1148	! interface description
1149	! input scalar
1150	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1151	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
1152	! input fields
1153	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1154	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation
1155	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! New soil temperature
1156	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
1157	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
1158	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
1159	! modified fields
1160	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation
1161	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation
1162	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
1163	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
1164	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
1165	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age!! Snow age on ice, lakes, ...
1166	! output fields
1167	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt
1168	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowdepth !! Snow depth
1169	!
1170	! local declaration
1171	!
1172	INTEGER(i_std) :: ji, jv
1173	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
1174	REAL(r_std), DIMENSION (kjpindex) :: xx !! temporary
1175	REAL(r_std) :: snowmelt_tmp !! The name says it all !
1176	LOGICAL, DIMENSION (kjpindex) :: warnings
1177	LOGICAL :: any_warning
1178	!
1179	! for continental points
1180	!
1181
1182	!
1183	! 0. initialisation
1184	!
1185	DO jv = 1, nnobio
1186	DO ji=1,kjpindex
1187	subsnownobio(ji,jv) = zero
1188	ENDDO
1189	ENDDO
1190	DO ji=1,kjpindex
1191	subsnowveg(ji) = zero
1192	snowmelt(ji) = zero
1193	icemelt(ji) = zero
1194	tot_melt(ji) = zero
1195	ENDDO
1196	!
1197	! 1. On vegetation
1198	!
1199	!cdir NODEP
1200	DO ji=1,kjpindex
1201	!
1202	! 1.1. It is snowing
1203	!
1204	snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji)
1205	ENDDO
1206	!
1207	DO ji=1,kjpindex
1208	!
1209	! 1.2. Sublimation - separate between vegetated and no-veget fractions
1210	! Care has to be taken as we might have sublimation from the
1211	! the frac_nobio while there is no snow on the rest of the grid.
1212	!
1213	IF ( snow(ji) > snowcri ) THEN
1214	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
1215	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
1216	ELSE
1217	! Correction Nathalie - Juillet 2006.
1218	! On doit d'abord tester s'il existe un frac_nobio!
1219	! Pour le moment je ne regarde que le iice
1220	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1221	subsnownobio(ji,iice) = vevapsno(ji)
1222	subsnowveg(ji) = zero
1223	ELSE
1224	subsnownobio(ji,iice) = zero
1225	subsnowveg(ji) = vevapsno(ji)
1226	ENDIF
1227	ENDIF
1228	ENDDO
1229	!
1230	warnings(:) = .FALSE.
1231	any_warning = .FALSE.
1232	!cdir NODEP
1233	DO ji=1,kjpindex
1234	!
1235	! 1.2.1 Check that sublimation on the vegetated fraction is possible.
1236	!
1237	IF (subsnowveg(ji) .GT. snow(ji)) THEN
1238	! What could not be sublimated goes into bare soil evaporation
1239	! Nathalie - Juillet 2006 - il faut avant tout tester s'il existe du
1240	! frac_nobio sur ce pixel pour eviter de puiser dans le sol!
1241	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1242	subsnownobio(ji,iice) = subsnownobio(ji,iice) + (subsnowveg(ji) - snow(ji))
1243	ELSE
1244	vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji))
1245	warnings(ji) = .TRUE.
1246	any_warning = .TRUE.
1247	ENDIF
1248	! Sublimation is thus limited to what is available
1249	subsnowveg(ji) = snow(ji)
1250	snow(ji) = zero
1251	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
1252	ELSE
1253	snow(ji) = snow(ji) - subsnowveg(ji)
1254	ENDIF
1255	ENDDO
1256	IF ( any_warning ) THEN
1257	WRITE(numout,*)' ATTENTION on prend de l eau au sol nu sur au moins un point car evapsno est trop fort!'
1258	!!$ DO ji=1,kjpindex
1259	!!$ IF ( warnings(ji) ) THEN
1260	!!$ WRITE(numout,*)' ATTENTION on prend de l eau au sol nu car evapsno est trop fort!'
1261	!!$ WRITE(numout,)' ',ji,' vevapnu (en mm/jour) = ',vevapnu(ji)one_day/dtradia
1262	!!$ ENDIF
1263	!!$ ENDDO
1264	ENDIF
1265	!
1266	warnings(:) = .FALSE.
1267	any_warning = .FALSE.
1268	!cdir NODEP
1269	DO ji=1,kjpindex
1270	!
1271	! 1.3. snow melt only if temperature positive
1272	!
1273	IF (temp_sol_new(ji).GT.tp_00) THEN
1274	!
1275	IF (snow(ji).GT.sneige) THEN
1276	!
1277	snowmelt(ji) = (un - frac_nobio(ji,iice))(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1278	!
1279	! 1.3.1.1 enough snow for melting or not
1280	!
1281	IF (snowmelt(ji).LT.snow(ji)) THEN
1282	snow(ji) = snow(ji) - snowmelt(ji)
1283	ELSE
1284	snowmelt(ji) = snow(ji)
1285	snow(ji) = zero
1286	END IF
1287	!
1288	ELSEIF (snow(ji).GE.zero) THEN
1289	!
1290	! 1.3.2 not enough snow
1291	!
1292	snowmelt(ji) = snow(ji)
1293	snow(ji) = zero
1294	ELSE
1295	!
1296	! 1.3.3 negative snow - now snow melt
1297	!
1298	snow(ji) = zero
1299	snowmelt(ji) = zero
1300	warnings(ji) = .TRUE.
1301	any_warning = .TRUE.
1302	!
1303	END IF
1304
1305	ENDIF
1306	ENDDO
1307	IF ( any_warning ) THEN
1308	DO ji=1,kjpindex
1309	IF ( warnings(ji) ) THEN
1310	WRITE(numout,*) 'hydrolc_snow: WARNING! snow was negative and was reset to zero for point ',ji,'. '
1311	ENDIF
1312	ENDDO
1313	ENDIF
1314	!
1315	DO ji=1,kjpindex
1316	!
1317	! 1.4. Ice melt only if there is more than a given mass : maxmass_glacier,
1318	! i.e. only weight melts glaciers !
1319	! Ajouts Edouard Davin / Nathalie de Noblet add extra to melting
1320	!
1321	IF ( snow(ji) .GT. maxmass_glacier ) THEN
1322	snowmelt(ji) = snowmelt(ji) + (snow(ji) - maxmass_glacier)
1323	snow(ji) = maxmass_glacier
1324	ENDIF
1325	!
1326	END DO
1327	!
1328	! 2. On Land ice
1329	!
1330	DO ji=1,kjpindex
1331	!
1332	! 2.1. It is snowing
1333	!
1334	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
1335	& frac_nobio(ji,iice)*precip_rain(ji)
1336	!
1337	! 2.2. Sublimation - was calculated before it can give us negative snow_nobio but that is OK
1338	! Once it goes below a certain values (-maxmass_glacier for instance) we should kill
1339	! the frac_nobio(ji,iice) !
1340	!
1341	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
1342	!
1343	! 2.3. snow melt only for continental ice fraction
1344	!
1345	snowmelt_tmp = zero
1346	IF (temp_sol_new(ji) .GT. tp_00) THEN
1347	!
1348	! 2.3.1 If there is snow on the ice-fraction it can melt
1349	!
1350	snowmelt_tmp = frac_nobio(ji,iice)(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1351	!
1352	IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN
1353	snowmelt_tmp = MAX( zero, snow_nobio(ji,iice))
1354	ENDIF
1355	snowmelt(ji) = snowmelt(ji) + snowmelt_tmp
1356	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp
1357	!
1358	ENDIF
1359	!
1360	! 2.4 Ice melt only if there is more than a given mass : maxmass_glacier,
1361	! i.e. only weight melts glaciers !
1362	!
1363	IF ( snow_nobio(ji,iice) .GT. maxmass_glacier ) THEN
1364	icemelt(ji) = snow_nobio(ji,iice) - maxmass_glacier
1365	snow_nobio(ji,iice) = maxmass_glacier
1366	ENDIF
1367	!
1368	END DO
1369
1370	!
1371	! 3. On other surface types - not done yet
1372	!
1373	IF ( nnobio .GT. 1 ) THEN
1374	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
1375	CALL ipslerr (3,'hydrolc_snow', '', &
1376	& 'CANNOT TREAT SNOW ON THESE SURFACE TYPES', '')
1377	ENDIF
1378
1379	!
1380	! 4. computes total melt (snow and ice)
1381	!
1382
1383	DO ji = 1, kjpindex
1384	tot_melt(ji) = icemelt(ji) + snowmelt(ji)
1385	ENDDO
1386
1387	!
1388	! 5. computes snow age on veg and ice (for albedo)
1389	!
1390	DO ji = 1, kjpindex
1391	!
1392	! 5.1 Snow age on vegetation
1393	!
1394	IF (snow(ji) .LE. zero) THEN
1395	snow_age(ji) = zero
1396	ELSE
1397	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dtradia/one_day) &
1398	& * EXP(-precip_snow(ji) / snow_trans)
1399	ENDIF
1400	!
1401	! 5.2 Snow age on ice
1402	!
1403	! age of snow on ice: a little bit different because in cold regions, we really
1404	! cannot negect the effect of cold temperatures on snow metamorphism any more.
1405	!
1406	IF (snow_nobio(ji,iice) .LE. zero) THEN
1407	snow_nobio_age(ji,iice) = zero
1408	ELSE
1409	!
1410	d_age(ji) = ( snow_nobio_age(ji,iice) + &
1411	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dtradia/one_day ) * &
1412	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
1413	IF (d_age(ji) .GT. zero ) THEN
1414	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
1415	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
1416	d_age(ji) = d_age(ji) / (un+xx(ji))
1417	ENDIF
1418	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
1419	!
1420	ENDIF
1421
1422	ENDDO
1423
1424	!
1425	! 6.0 Diagnose the depth of the snow layer
1426	!
1427
1428	DO ji = 1, kjpindex
1429	snowdepth(ji) = snow(ji) /sn_dens
1430	ENDDO
1431
1432	IF (long_print) WRITE (numout,*) ' hydrolc_snow done '
1433
1434	END SUBROUTINE hydrolc_snow
1435
1436	!! This routine computes canopy processes
1437	!!
1438	SUBROUTINE hydrolc_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg, precisol)
1439
1440	!
1441	! interface description
1442	!
1443	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1444	! input fields
1445	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation
1446	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: vevapwet !! Interception loss
1447	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation type
1448	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum water on vegetation for interception
1449	! modified fields
1450	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Water on vegetation due to interception
1451	! output fields
1452	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Eau tombee sur le sol
1453
1454	!
1455	! local declaration
1456	!
1457	INTEGER(i_std) :: ji, jv
1458	REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew
1459	LOGICAL, SAVE :: firstcall=.TRUE.
1460	REAL(r_std), SAVE, DIMENSION(nvm) :: throughfall_by_pft
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: