Context Navigation

source: tags/ORCHIDEE/src_sechiba/hydrolc.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 104.7 KB

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

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

Download in other formats: