Context Navigation

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

Last change on this file since 282 was 282, checked in by didier.solyga, 13 years ago
Move some labels associated to externalized parameters in sechiba to pft_parameters.f90
File size: 104.2 KB

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

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

Download in other formats: