Context Navigation

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

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

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

Download in other formats: