Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/hydrolc.f90 @ 6268

Last change on this file since 6268 was 843, checked in by didier.solyga, 12 years ago
Improved parameters documentation and presentation. The presentation is the same for all paramters so a script can automatically build a orchidee.default file.
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 104.2 KB

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

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

Download in other formats: