Context Navigation

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

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

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

Download in other formats: