Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE/src_sechiba/hydrol.f90 @ 880

Last change on this file since 880 was 566, checked in by didier.solyga, 13 years ago
Update Parameters labels : add Config Units and Config If for all parameters
Property svn:keywords set to `Revision Date HeadURL Date Author Revision`
File size: 157.2 KB

Line
1	!!
2	!! This module computes hydrologic processes on continental points.
3	!!
4	!! @author Marie-Alice Foujols and Jan Polcher
5	!! @Version : $Revision$, $Date$
6	!!
7	!< $HeadURL$
8	!< $Date$
9	!< $Author$
10	!< $Revision$
11	!! IPSL (2006)
12	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13	!!
14	MODULE hydrol
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
26	! for debug :
27	USE grid
28
29	IMPLICIT NONE
30
31	! public routines :
32	! hydrol
33	PRIVATE
34	PUBLIC :: hydrol_main,hydrol_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=.TRUE. !! The check the water balance
42	LOGICAL, SAVE :: check_cwrr=.TRUE. !! The check the water balance
43	REAL(r_std), PARAMETER :: allowed_err = 1.0E-8_r_std
44
45	! one dimension array allocated, computed, saved and got in hydrol module
46
47	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_beg !! Total amount of water at start of time step
48	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_water_end !! Total amount of water at end of time step
49	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_beg !! Total amount of water on vegetation at start of time step
50	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watveg_end !! Total amount of water on vegetation at end of time step
51	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_beg !! Total amount of water in the soil at start of time step
52	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_watsoil_end !! Total amount of water in the soil at end of time step
53	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_beg !! Total amount of snow at start of time step
54	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_end !! Total amount of snow at end of time step
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delsoilmoist !! Change in soil moisture
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delintercept !! Change in interception storage
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: delswe !! Change in SWE^Q
58
59	! array allocated, computed, saved and got in hydrol module
60	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_veget !! zero/one when veget fraction is zero/higher
61	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: mask_soiltype !! zero/one where soil fraction is zero/higher
62	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mask_corr_veg_soil !! zero/one where veg frac on a soil type is zero/higher
63	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mask_return !! zero/one where there is no/is returnflow
64	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: index_nsat !! Indices of the non-saturated layers
65	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: index_sat !! Indices of the saturated layers
66	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: n_nsat !! Number of non-saturated layers
67	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: n_sat !! Number of saturated layers
68	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: nslme !! last efficient layer
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: humrelv !! humrel for each soil type
70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: vegstressv !! vegstress for each soil type
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:,:):: us !! relative humidity
72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol !! Eau tombee sur le sol
73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: precisol_ns !! Eau tombee sur le sol par type de sol
74	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ae_ns !! Evaporation du sol nu par type de sol
75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: evap_bare_lim_ns !! limitation of bare soil evaporation on each soil column (used to deconvoluate vevapnu)
76	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: free_drain_coef !! Coefficient for free drainage at bottom
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: rootsink !! stress racinaire par niveau et type de sol
78	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsnowveg !! Sublimation of snow on vegetation
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: subsnownobio !! Sublimation of snow on other surface types (ice, lakes, ...)
80	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snowmelt !! Quantite de neige fondue
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: icemelt !! Quantite de glace fondue
82	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: subsinksoil !! Excess of sublimation as a sink for the soil
83	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vegtot !! Total vegetation
84	! The last vegetation map which was used to distribute the reservoirs
85	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: resdist !! Distribution of reservoirs
86	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: mx_eau_var
87
88	! arrays used by cwrr scheme
89	REAL(r_std), SAVE, DIMENSION (nslm+1,nstm) :: zz !!
90	REAL(r_std), SAVE, DIMENSION (nslm+1,nstm) :: dz !!
91	REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: mc_lin !!
92	REAL(r_std), SAVE, DIMENSION (nstm) :: v1r !! Residual soil water content of the first layer
93	REAL(r_std), SAVE, DIMENSION (nstm) :: vBs !! Saturated soil water content of the bottom layer
94
95	REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: k_lin !!
96	REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: d_lin !!
97	REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: a_lin !!
98	REAL(r_std), SAVE, DIMENSION (imin:imax,nstm) :: b_lin !!
99
100	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: humtot !! (:)
101	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: flux !! (:)
102	LOGICAL, ALLOCATABLE, SAVE, DIMENSION (:) :: resolv !! (:)
103
104
105	!! linarization coefficients of hydraulic conductivity K
106	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: a !! (:,nslm)
107	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: b !!
108	!! linarization coefficients of hydraulic diffusivity D
109	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: d !!
110
111	!! matrix coefficients
112	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: e !!
113	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: f !!
114	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: g1 !!
115
116
117	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ep !!
118	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: fp !!
119	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gp !!
120	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rhs !!
121	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: srhs !!
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gam !!
123
124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc !! (:,nstm) Total moisture content (mm)
125	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmcs !! (nstm) Total moisture constent at saturation (mm)
126	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter !! (:,nstm) Total moisture in the litter by soil type
127	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_mea !! Total moisture in the litter over the grid
128
129	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_wilt !! (:,nstm) Moisture of litter at wilt pt
130	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_field !! (:,nstm) Moisture of litter at field cap.
131	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_res !! (:,nstm) Moisture of litter at residual moisture.
132	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_sat !! (:,nstm) Moisture of litter at saturatiion
133	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_awet !! (:,nstm) Moisture of litter at mc_awet
134	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tmc_litter_adry !! (:,nstm) Moisture of litter at mc_dry
135	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_wet_mea !! Total moisture in the litter over the grid below which albedo is fixed
136	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tmc_litt_dry_mea !! Total moisture in the litter over the grid above which albedo is fixed
137
138	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: v1 !! (:)
139	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vB !! (:)
140	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qflux00 !! flux at the top of the soil column
141
142	!! par type de sol :
143	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: ru_ns !! ruissellement
144	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: dr_ns !! drainage
145	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: tr_ns !! transpiration
146
147	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: corr_veg_soil !! (:,nvm,nstm)
148	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: corr_veg_soil_max !! (:,nvm,nstm)
149	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc !! (:,nslm,nstm) m³ x m³
150	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soilmoist !! (:,nslm)
151	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: soil_wet !! soil wetness
152	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soil_wet_litter !! soil wetness of the litter
153	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: qflux !! fluxes between the soil layers
154
155	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: tmat !!
156	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: stmat !!
157
158	LOGICAL, SAVE :: interpol_diag=.FALSE.
159
160	CONTAINS
161
162	!!
163	!! Main routine for hydrol module
164	!! - called only one time for initialisation
165	!! - called every time step
166	!! - called one more time at last time step for writing _restart_ file
167	!!
168	!! Algorithm:
169	!! - call hydrol_snow for snow process (including age of snow)
170	!! - call hydrol_canop for canopy process
171	!! - call hydrol_soil for bare soil process
172	!!
173	!! @call hydrol_snow
174	!! @call hydrol_canop
175	!! @call hydrol_soil
176	!!
177	SUBROUTINE hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
178	& index, indexveg, indexsoil, indexlayer,&
179	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, &
180	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
181	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, &
182	& humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, &
183	& shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
184
185	! interface description
186	! input scalar
187	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
188	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
189	INTEGER(i_std),INTENT (in) :: rest_id,hist_id !! _Restart_ file and _history_ file identifier
190	INTEGER(i_std),INTENT (in) :: hist2_id !! _history_ file 2 identifier
191	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
192	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
193	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
194	! input fields
195	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
196	INTEGER(i_std),DIMENSION (kjpindex*nvm), INTENT (in):: indexveg !! Indeces of the points on the 3D map for veg
197	INTEGER(i_std),DIMENSION (kjpindex*nstm), INTENT (in):: indexsoil !! Indeces of the points on the 3D map for soil
198	INTEGER(i_std),DIMENSION (kjpindex*nslm), INTENT (in):: indexlayer !! Indeces of the points on the 3D map for soil layers
199	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
200	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
201	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Routed water which comes back into the soil
202	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation returning to soil moisture
203	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
204
205	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: frac_nobio !! Fraction of ice, lakes, ...
206	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio !! Total fraction of ice+lakes+...
207	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: soilcap !! Soil capacity
208	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
209	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss
210	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
211	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type (LAI -> infty)
212	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintmax !! Maximum water on vegetation for interception
213	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
214	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot !! Soil Potential Evaporation
215	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: evapot_penm !! Soil Potential Evaporation Correction
216	! modified fields
217	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !!
218	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapnu !! Bare soil evaporation
219	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: vevapsno !! Snow evaporation
220	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow !! Snow mass [Kg/m^2]
221	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: snow_age !! Snow age
222	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2]
223	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (inout) :: snow_nobio_age !! Snow age on ice, lakes, ...
224	!! We consider that any water on the ice is snow and we only peforme a water balance to have consistency.
225	!! The water balance is limite to + or - 10^6 so that accumulation is not endless
226	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: runoff !! Complete runoff
227	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: drainage !! Drainage
228	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (inout):: shumdiag !! relative soil moisture
229	! output fields
230	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Relative humidity
231	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
232	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! function of litter wetness
233	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
234	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tot_melt !! Total melt
235	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception
236
237	!
238	! local declaration
239	!
240	INTEGER(i_std) :: jst, jsl
241	REAL(r_std),DIMENSION (kjpindex) :: soilwet !! A temporary diagnostic of soil wetness
242	REAL(r_std),DIMENSION (kjpindex) :: snowdepth !! Depth of snow layer
243	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
244
245	!
246	! do initialisation
247	!
248	IF (l_first_hydrol) THEN
249
250	IF (long_print) WRITE (numout,*) ' l_first_hydrol : call hydrol_init '
251
252	CALL hydrol_init (kjit, ldrestart_read, kjpindex, index, rest_id, veget, soiltype, humrel,&
253	& vegstress, snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
254
255	CALL hydrol_var_init (kjpindex, veget, soiltype, mx_eau_var, shumdiag, litterhumdiag, &
256	& drysoil_frac, evap_bare_lim)
257	!
258	! If we check the water balance we first save the total amount of water
259	!
260	IF (check_waterbal) THEN
261	CALL hydrol_waterbal(kjpindex, index, .TRUE., dtradia, veget, &
262	& totfrac_nobio, qsintveg, snow, snow_nobio,&
263	& precip_rain, precip_snow, returnflow, irrigation, tot_melt, &
264	& vevapwet, transpir, vevapnu, vevapsno, runoff,drainage)
265	ENDIF
266	!
267	IF (almaoutput) THEN
268	CALL hydrol_alma(kjpindex, index, .TRUE., qsintveg, snow, snow_nobio, soilwet)
269	ENDIF
270
271	RETURN
272
273	ENDIF
274
275	!
276	! prepares restart file for the next simulation
277	!
278	IF (ldrestart_write) THEN
279
280	IF (long_print) WRITE (numout,*) ' we have to complete restart file with HYDROLOGIC variables '
281
282	DO jst=1,nstm
283	! var_name= "mc_1" ... "mc_3"
284	WRITE (var_name,"('moistc_',i1)") jst
285	CALL restput_p(rest_id, var_name, nbp_glo, nslm, 1, kjit, mc(:,:,jst), 'scatter', nbp_glo, index_g)
286	END DO
287	!
288	DO jst=1,nstm
289	DO jsl=1,nslm
290	! var_name= "us_1_01" ... "us_3_11"
291	WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
292	CALL restput_p(rest_id, var_name, nbp_glo,nvm, 1,kjit,us(:,:,jst,jsl),'scatter',nbp_glo,index_g)
293	END DO
294	END DO
295	!
296	var_name= 'free_drain_coef'
297	CALL restput_p(rest_id, var_name, nbp_glo, nstm, 1, kjit, free_drain_coef, 'scatter', nbp_glo, index_g)
298	!
299	var_name= 'ae_ns'
300	CALL restput_p(rest_id, var_name, nbp_glo, nstm, 1, kjit, ae_ns, 'scatter', nbp_glo, index_g)
301	!
302	var_name= 'vegstress'
303	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, vegstress, 'scatter', nbp_glo, index_g)
304	!
305	var_name= 'snow'
306	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow, 'scatter', nbp_glo, index_g)
307	!
308	var_name= 'snow_age'
309	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, snow_age, 'scatter', nbp_glo, index_g)
310	!
311	var_name= 'snow_nobio'
312	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio, 'scatter', nbp_glo, index_g)
313	!
314	var_name= 'snow_nobio_age'
315	CALL restput_p(rest_id, var_name, nbp_glo, nnobio, 1, kjit, snow_nobio_age, 'scatter', nbp_glo, index_g)
316	!
317	var_name= 'qsintveg'
318	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, qsintveg, 'scatter', nbp_glo, index_g)
319	!
320	var_name= 'resdist'
321	CALL restput_p(rest_id, var_name, nbp_glo, nvm, 1, kjit, resdist, 'scatter', nbp_glo, index_g)
322	RETURN
323	!
324	END IF
325
326	!
327	! shared time step
328	!
329	IF (long_print) WRITE (numout,*) 'hydrol pas de temps = ',dtradia
330
331	!
332	! computes snow
333	!
334
335	CALL hydrol_snow(kjpindex, dtradia, precip_rain, precip_snow, temp_sol_new, soilcap, &
336	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
337	& tot_melt, snowdepth)
338
339	!
340	! computes canopy
341	!
342	!
343	CALL hydrol_vegupd(kjpindex, veget, veget_max, soiltype, qsintveg,resdist)
344	!
345
346	CALL hydrol_canop(kjpindex, precip_rain, vevapwet, veget, qsintmax, qsintveg,precisol,tot_melt)
347
348	! computes hydro_soil
349	!
350
351	CALL hydrol_soil(kjpindex, dtradia, veget, veget_max, soiltype, transpir, vevapnu, evapot, &
352	& evapot_penm, runoff, drainage, returnflow, irrigation, &
353	& tot_melt,evap_bare_lim, shumdiag, litterhumdiag, humrel, vegstress, drysoil_frac)
354	!
355	! If we check the water balance we end with the comparison of total water change and fluxes
356	!
357	IF (check_waterbal) THEN
358	CALL hydrol_waterbal(kjpindex, index, .FALSE., dtradia, veget, totfrac_nobio, &
359	& qsintveg, snow,snow_nobio, precip_rain, precip_snow, returnflow, &
360	& irrigation, tot_melt, vevapwet, transpir, vevapnu, vevapsno, runoff, drainage)
361	ENDIF
362	!
363	! If we use the ALMA standards
364	!
365	IF (almaoutput) THEN
366	CALL hydrol_alma(kjpindex, index, .FALSE., qsintveg, snow, snow_nobio, soilwet)
367	ENDIF
368
369	!
370	IF ( .NOT. almaoutput ) THEN
371	DO jst=1,nstm
372	! var_name= "mc_1" ... "mc_3"
373	WRITE (var_name,"('moistc_',i1)") jst
374	CALL histwrite(hist_id, trim(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
375
376	! var_name= "vegetsoil_1" ... "vegetsoil_3"
377	WRITE (var_name,"('vegetsoil_',i1)") jst
378	CALL histwrite(hist_id, trim(var_name), kjit,corr_veg_soil(:,:,jst), kjpindex*nvm, indexveg)
379	ENDDO
380	CALL histwrite(hist_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
381	CALL histwrite(hist_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
382	CALL histwrite(hist_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
383	CALL histwrite(hist_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
384	CALL histwrite(hist_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
385	CALL histwrite(hist_id, 'humtot', kjit, humtot, kjpindex, index)
386	CALL histwrite(hist_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
387	CALL histwrite(hist_id, 'drainage', kjit, drainage, kjpindex, index)
388	CALL histwrite(hist_id, 'runoff', kjit, runoff, kjpindex, index)
389	CALL histwrite(hist_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
390	CALL histwrite(hist_id, 'rain', kjit, precip_rain, kjpindex, index)
391	CALL histwrite(hist_id, 'snowf', kjit, precip_snow, kjpindex, index)
392	CALL histwrite(hist_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
393	CALL histwrite(hist_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
394	IF ( hist2_id > 0 ) THEN
395	DO jst=1,nstm
396	! var_name= "mc_1" ... "mc_3"
397	WRITE (var_name,"('moistc_',i1)") jst
398	CALL histwrite(hist2_id, trim(var_name), kjit,mc(:,:,jst), kjpindex*nslm, indexlayer)
399
400	! var_name= "vegetsoil_1" ... "vegetsoil_3"
401	WRITE (var_name,"('vegetsoil_',i1)") jst
402	CALL histwrite(hist2_id, trim(var_name), kjit,corr_veg_soil(:,:,jst), kjpindex*nvm, indexveg)
403	ENDDO
404	CALL histwrite(hist2_id, 'evapnu_soil', kjit, ae_ns, kjpindex*nstm, indexsoil)
405	CALL histwrite(hist2_id, 'drainage_soil', kjit, dr_ns, kjpindex*nstm, indexsoil)
406	CALL histwrite(hist2_id, 'transpir_soil', kjit, tr_ns, kjpindex*nstm, indexsoil)
407	CALL histwrite(hist2_id, 'runoff_soil', kjit, ru_ns, kjpindex*nstm, indexsoil)
408	CALL histwrite(hist2_id, 'humtot_soil', kjit, tmc, kjpindex*nstm, indexsoil)
409	CALL histwrite(hist2_id, 'humtot', kjit, humtot, kjpindex, index)
410	CALL histwrite(hist2_id, 'humrel', kjit, humrel, kjpindex*nvm, indexveg)
411	CALL histwrite(hist2_id, 'drainage', kjit, drainage, kjpindex, index)
412	CALL histwrite(hist2_id, 'runoff', kjit, runoff, kjpindex, index)
413	CALL histwrite(hist2_id, 'precisol', kjit, precisol, kjpindex*nvm, indexveg)
414	CALL histwrite(hist2_id, 'rain', kjit, precip_rain, kjpindex, index)
415	CALL histwrite(hist2_id, 'snowf', kjit, precip_snow, kjpindex, index)
416	CALL histwrite(hist2_id, 'qsintmax', kjit, qsintmax, kjpindex*nvm, indexveg)
417	CALL histwrite(hist2_id, 'qsintveg', kjit, qsintveg, kjpindex*nvm, indexveg)
418	ENDIF
419	ELSE
420	CALL histwrite(hist_id, 'Snowf', kjit, precip_snow, kjpindex, index)
421	CALL histwrite(hist_id, 'Rainf', kjit, precip_rain, kjpindex, index)
422	CALL histwrite(hist_id, 'Qs', kjit, runoff, kjpindex, index)
423	CALL histwrite(hist_id, 'Qsb', kjit, drainage, kjpindex, index)
424	CALL histwrite(hist_id, 'Qsm', kjit, tot_melt, kjpindex, index)
425	CALL histwrite(hist_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
426	CALL histwrite(hist_id, 'DelSWE', kjit, delswe, kjpindex, index)
427	CALL histwrite(hist_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
428	!
429	CALL histwrite(hist_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
430	CALL histwrite(hist_id, 'SoilWet', kjit, soilwet, kjpindex, index)
431	!
432	CALL histwrite(hist_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
433	CALL histwrite(hist_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
434	!
435	CALL histwrite(hist_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
436	!
437	IF ( hist2_id > 0 ) THEN
438	CALL histwrite(hist2_id, 'Snowf', kjit, precip_snow, kjpindex, index)
439	CALL histwrite(hist2_id, 'Rainf', kjit, precip_rain, kjpindex, index)
440	CALL histwrite(hist2_id, 'Qs', kjit, runoff, kjpindex, index)
441	CALL histwrite(hist2_id, 'Qsb', kjit, drainage, kjpindex, index)
442	CALL histwrite(hist2_id, 'Qsm', kjit, tot_melt, kjpindex, index)
443	CALL histwrite(hist2_id, 'DelSoilMoist', kjit, delsoilmoist, kjpindex, index)
444	CALL histwrite(hist2_id, 'DelSWE', kjit, delswe, kjpindex, index)
445	CALL histwrite(hist2_id, 'DelIntercept', kjit, delintercept, kjpindex, index)
446	!
447	CALL histwrite(hist2_id, 'SoilMoist', kjit, soilmoist, kjpindex*nslm, indexlayer)
448	CALL histwrite(hist2_id, 'SoilWet', kjit, soilwet, kjpindex, index)
449	!
450	CALL histwrite(hist2_id, 'RootMoist', kjit, tot_watsoil_end, kjpindex, index)
451	CALL histwrite(hist2_id, 'SubSnow', kjit, vevapsno, kjpindex, index)
452	!
453	CALL histwrite(hist2_id, 'SnowDepth', kjit, snowdepth, kjpindex, index)
454	ENDIF
455	ENDIF
456
457	IF (long_print) WRITE (numout,*) ' hydrol_main Done '
458
459	END SUBROUTINE hydrol_main
460
461	!! Algorithm:
462	!! - dynamic allocation for local array
463	!! - _restart_ file reading for HYDROLOGIC variables
464	!!
465	SUBROUTINE hydrol_init(kjit, ldrestart_read, kjpindex, index, rest_id, veget, soiltype, humrel,&
466	& vegstress, snow, snow_age, snow_nobio, snow_nobio_age, qsintveg)
467
468	! interface description
469	! input scalar
470	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
471	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for _restart_ file to read
472	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
473	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
474	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
475	! input fields
476	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Carte de vegetation
477	REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
478	! output fields
479	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Stress hydrique, relative humidity
480	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
481	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow !! Snow mass [Kg/m^2]
482	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: snow_age !! Snow age
483	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio !! Snow on ice, lakes, ...
484	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (out) :: snow_nobio_age !! Snow age on ice, lakes, ...
485	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: qsintveg !! Water on vegetation due to interception
486
487	! local declaration
488	INTEGER(i_std) :: ier, ierror, ipdt
489	INTEGER(i_std) :: ji, jv, jst, jsl, ik
490	CHARACTER(LEN=80) :: var_name !! To store variables names for I/O
491
492	! initialisation
493	IF (l_first_hydrol) THEN
494	l_first_hydrol=.FALSE.
495	ELSE
496	WRITE (numout,*) ' l_first_hydrol false . we stop '
497	STOP 'hydrol_init'
498	ENDIF
499
500	sneige = snowcri/mille
501	throughfall_by_pft = throughfall_by_pft / 100.
502
503	! make dynamic allocation with good dimension
504
505	! one dimension array allocation with possible restart value
506
507	ALLOCATE (mask_corr_veg_soil(kjpindex,nvm,nstm),stat=ier)
508	IF (ier.NE.0) THEN
509	WRITE (numout,) ' error in mask_corr_veg_soil allocation. We stop. We need kjpindex words = ',kjpindexnvm*nstm
510	STOP 'hydrol_init'
511	END IF
512
513	ALLOCATE (mask_veget(kjpindex,nvm),stat=ier)
514	IF (ier.NE.0) THEN
515	WRITE (numout,) ' error in mask_veget allocation. We stop. We need kjpindex words = ',kjpindexnvm
516	STOP 'hydrol_init'
517	END IF
518
519	ALLOCATE (mask_soiltype(kjpindex,nstm),stat=ier)
520	IF (ier.NE.0) THEN
521	WRITE (numout,) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindexnstm
522	STOP 'hydrol_init'
523	END IF
524
525	ALLOCATE (mask_return(kjpindex),stat=ier)
526	IF (ier.NE.0) THEN
527	WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindex
528	STOP 'hydrol_init'
529	END IF
530
531	ALLOCATE (index_nsat(kjpindex,nstm),stat=ier)
532	IF (ier.NE.0) THEN
533	WRITE (numout,) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindexnstm
534	STOP 'hydrol_init'
535	END IF
536
537	ALLOCATE (index_sat(kjpindex,nstm),stat=ier)
538	IF (ier.NE.0) THEN
539	WRITE (numout,) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindexnstm
540	STOP 'hydrol_init'
541	END IF
542
543	ALLOCATE (n_nsat(nstm),stat=ier)
544	IF (ier.NE.0) THEN
545	WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',nstm
546	STOP 'hydrol_init'
547	END IF
548
549	ALLOCATE (n_sat(nstm),stat=ier)
550	IF (ier.NE.0) THEN
551	WRITE (numout,*) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',nstm
552	STOP 'hydrol_init'
553	END IF
554
555	ALLOCATE (nslme(kjpindex,nstm),stat=ier)
556	IF (ier.NE.0) THEN
557	WRITE (numout,) ' error in mask_soiltype allocation. We stop. We need kjpindex words = ',kjpindexnstm
558	STOP 'hydrol_init'
559	END IF
560
561	ALLOCATE (humrelv(kjpindex,nvm,nstm),stat=ier)
562	IF (ier.NE.0) THEN
563	WRITE (numout,) ' error in humrelv allocation. We stop. We need kjpindex words = ',kjpindexnvm*nstm
564	STOP 'hydrol_init'
565	END IF
566
567	ALLOCATE (vegstressv(kjpindex,nvm,nstm),stat=ier)
568	IF (ier.NE.0) THEN
569	WRITE (numout,) ' error in vegstressv allocation. We stop. We need kjpindex words = ',kjpindexnvm*nstm
570	STOP 'hydrol_init'
571	END IF
572
573	ALLOCATE (us(kjpindex,nvm,nstm,nslm),stat=ier)
574	IF (ier.NE.0) THEN
575	WRITE (numout,) ' error in us allocation. We stop. We need kjpindex words = ',kjpindexnvmnstmnslm
576	STOP 'hydrol_init'
577	END IF
578
579	ALLOCATE (precisol(kjpindex,nvm),stat=ier)
580	IF (ier.NE.0) THEN
581	WRITE (numout,) ' error in precisol allocation. We stop. We need kjpindex words = ',kjpindexnvm
582	STOP 'hydrol_init'
583	END IF
584
585	ALLOCATE (precisol_ns(kjpindex,nstm),stat=ier)
586	IF (ier.NE.0) THEN
587	WRITE (numout,) ' error in precisol_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
588	STOP 'hydrol_init'
589	END IF
590
591	ALLOCATE (free_drain_coef(kjpindex,nstm),stat=ier)
592	IF (ier.NE.0) THEN
593	WRITE (numout,) ' error in free_drain_coef allocation. We stop. We need kjpindex words = ',kjpindexnstm
594	STOP 'hydrol_init'
595	END IF
596
597	ALLOCATE (ae_ns(kjpindex,nstm),stat=ier)
598	IF (ier.NE.0) THEN
599	WRITE (numout,) ' error in ae_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
600	STOP 'hydrol_init'
601	END IF
602
603	ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier)
604	IF (ier.NE.0) THEN
605	WRITE (numout,) ' error in evap_bare_lim_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
606	STOP 'hydrol_init'
607	END IF
608
609	ALLOCATE (rootsink(kjpindex,nslm,nstm),stat=ier)
610	IF (ier.NE.0) THEN
611	WRITE (numout,) ' error in rootsink allocation. We stop. We need kjpindex words = ',kjpindexnslm*nstm
612	STOP 'hydrol_init'
613	END IF
614
615	ALLOCATE (subsnowveg(kjpindex),stat=ier)
616	IF (ier.NE.0) THEN
617	WRITE (numout,*) ' error in subsnowveg allocation. We stop. We need kjpindex words = ',kjpindex
618	STOP 'hydrol_init'
619	END IF
620
621	ALLOCATE (subsnownobio(kjpindex,nnobio),stat=ier)
622	IF (ier.NE.0) THEN
623	WRITE (numout,) ' error in subsnownobio allocation. We stop. We need kjpindex words = ',kjpindexnnobio
624	STOP 'hydrol_init'
625	END IF
626
627	ALLOCATE (snowmelt(kjpindex),stat=ier)
628	IF (ier.NE.0) THEN
629	WRITE (numout,*) ' error in snowmelt allocation. We stop. We need kjpindex words = ',kjpindex
630	STOP 'hydrol_init'
631	END IF
632
633	ALLOCATE (icemelt(kjpindex),stat=ier)
634	IF (ier.NE.0) THEN
635	WRITE (numout,*) ' error in icemelt allocation. We stop. We need kjpindex words = ',kjpindex
636	STOP 'hydrol_init'
637	END IF
638
639	ALLOCATE (subsinksoil(kjpindex),stat=ier)
640	IF (ier.NE.0) THEN
641	WRITE (numout,*) ' error in subsinksoil allocation. We stop. We need kjpindex words = ',kjpindex
642	STOP 'hydrol_init'
643	END IF
644
645	ALLOCATE (mx_eau_var(kjpindex),stat=ier)
646	IF (ier.NE.0) THEN
647	WRITE (numout,*) ' error in mx_eau_var allocation. We stop. We need kjpindex words = ',kjpindex
648	STOP 'hydrol_init'
649	END IF
650
651	ALLOCATE (vegtot(kjpindex),stat=ier)
652	IF (ier.NE.0) THEN
653	WRITE (numout,*) ' error in vegtot allocation. We stop. We need kjpindex words = ',kjpindex
654	STOP 'hydrol_init'
655	END IF
656
657	ALLOCATE (resdist(kjpindex,nvm),stat=ier)
658	IF (ier.NE.0) THEN
659	WRITE (numout,) ' error in resdist allocation. We stop. We need kjpindex words = ',kjpindexnvm
660	STOP 'hydrol_init'
661	END IF
662
663	ALLOCATE (humtot(kjpindex),stat=ier)
664	IF (ier.NE.0) THEN
665	WRITE (numout,*) ' error in humtot allocation. We stop. We need kjpindex words = ',kjpindex
666	STOP 'hydrol_init'
667	END IF
668
669	ALLOCATE (flux(kjpindex),stat=ier)
670	IF (ier.NE.0) THEN
671	WRITE (numout,*) ' error in flux allocation. We stop. We need kjpindex words = ',kjpindex
672	STOP 'hydrol_init'
673	END IF
674
675	ALLOCATE (resolv(kjpindex),stat=ier)
676	IF (ier.NE.0) THEN
677	WRITE (numout,*) ' error in resolv allocation. We stop. We need kjpindex words = ',kjpindex
678	STOP 'hydrol_init'
679	END IF
680
681	ALLOCATE (a(kjpindex,nslm),stat=ier)
682	IF (ier.NE.0) THEN
683	WRITE (numout,) ' error in a allocation. We stop. We need kjpindex words = ',kjpindexnslm
684	STOP 'hydrol_init'
685	END IF
686
687	ALLOCATE (b(kjpindex,nslm),stat=ier)
688	IF (ier.NE.0) THEN
689	WRITE (numout,) ' error in b allocation. We stop. We need kjpindex words = ',kjpindexnslm
690	STOP 'hydrol_init'
691	END IF
692
693	ALLOCATE (d(kjpindex,nslm),stat=ier)
694	IF (ier.NE.0) THEN
695	WRITE (numout,) ' error in d allocation. We stop. We need kjpindex words = ',kjpindexnslm
696	STOP 'hydrol_init'
697	END IF
698
699	ALLOCATE (e(kjpindex,nslm),stat=ier)
700	IF (ier.NE.0) THEN
701	WRITE (numout,) ' error in e allocation. We stop. We need kjpindex words = ',kjpindexnslm
702	STOP 'hydrol_init'
703	END IF
704
705	ALLOCATE (f(kjpindex,nslm),stat=ier)
706	IF (ier.NE.0) THEN
707	WRITE (numout,) ' error in f allocation. We stop. We need kjpindex words = ',kjpindexnslm
708	STOP 'hydrol_init'
709	END IF
710
711	ALLOCATE (g1(kjpindex,nslm),stat=ier)
712	IF (ier.NE.0) THEN
713	WRITE (numout,) ' error in g1 allocation. We stop. We need kjpindex words = ',kjpindexnslm
714	STOP 'hydrol_init'
715	END IF
716
717	ALLOCATE (ep(kjpindex,nslm),stat=ier)
718	IF (ier.NE.0) THEN
719	WRITE (numout,) ' error in ep allocation. We stop. We need kjpindex words = ',kjpindexnslm
720	STOP 'hydrol_init'
721	END IF
722
723	ALLOCATE (fp(kjpindex,nslm),stat=ier)
724	IF (ier.NE.0) THEN
725	WRITE (numout,) ' error in fp allocation. We stop. We need kjpindex words = ',kjpindexnslm
726	STOP 'hydrol_init'
727	END IF
728
729	ALLOCATE (gp(kjpindex,nslm),stat=ier)
730	IF (ier.NE.0) THEN
731	WRITE (numout,) ' error in gp allocation. We stop. We need kjpindex words = ',kjpindexnslm
732	STOP 'hydrol_init'
733	END IF
734
735	ALLOCATE (rhs(kjpindex,nslm),stat=ier)
736	IF (ier.NE.0) THEN
737	WRITE (numout,) ' error in rhs allocation. We stop. We need kjpindex words = ',kjpindexnslm
738	STOP 'hydrol_init'
739	END IF
740
741	ALLOCATE (srhs(kjpindex,nslm),stat=ier)
742	IF (ier.NE.0) THEN
743	WRITE (numout,) ' error in srhs allocation. We stop. We need kjpindex words = ',kjpindexnslm
744	STOP 'hydrol_init'
745	END IF
746
747	ALLOCATE (gam(kjpindex,nslm),stat=ier)
748	IF (ier.NE.0) THEN
749	WRITE (numout,) ' error in gam allocation. We stop. We need kjpindex words = ',kjpindexnslm
750	STOP 'hydrol_init'
751	END IF
752
753	ALLOCATE (tmc(kjpindex,nstm),stat=ier)
754	IF (ier.NE.0) THEN
755	WRITE (numout,) ' error in tmc allocation. We stop. We need kjpindex words = ',kjpindexnstm
756	STOP 'hydrol_init'
757	END IF
758
759	ALLOCATE (tmcs(nstm),stat=ier)
760	IF (ier.NE.0) THEN
761	WRITE (numout,*) ' error in tmcs allocation. We stop. We need kjpindex words = ',nstm
762	STOP 'hydrol_init'
763	END IF
764
765	ALLOCATE (tmc_litter(kjpindex,nstm),stat=ier)
766	IF (ier.NE.0) THEN
767	WRITE (numout,) ' error in tmc_litter allocation. We stop. We need kjpindex words = ',kjpindexnstm
768	STOP 'hydrol_init'
769	END IF
770
771	ALLOCATE (tmc_litt_mea(kjpindex),stat=ier)
772	IF (ier.NE.0) THEN
773	WRITE (numout,*) ' error in tmc_litt_mea allocation. We stop. We need kjpindex words = ',kjpindex
774	STOP 'hydrol_init'
775	END IF
776
777	ALLOCATE (tmc_litter_res(kjpindex,nstm),stat=ier)
778	IF (ier.NE.0) THEN
779	WRITE (numout,) ' error in tmc_litter_res allocation. We stop. We need kjpindex words = ',kjpindexnstm
780	STOP 'hydrol_init'
781	END IF
782
783	ALLOCATE (tmc_litter_wilt(kjpindex,nstm),stat=ier)
784	IF (ier.NE.0) THEN
785	WRITE (numout,) ' error in tmc_litter_wilt allocation. We stop. We need kjpindex words = ',kjpindexnstm
786	STOP 'hydrol_init'
787	END IF
788
789	ALLOCATE (tmc_litter_field(kjpindex,nstm),stat=ier)
790	IF (ier.NE.0) THEN
791	WRITE (numout,) ' error in tmc_litter_field allocation. We stop. We need kjpindex words = ',kjpindexnstm
792	STOP 'hydrol_init'
793	END IF
794
795	ALLOCATE (tmc_litter_sat(kjpindex,nstm),stat=ier)
796	IF (ier.NE.0) THEN
797	WRITE (numout,) ' error in tmc_litter_sat allocation. We stop. We need kjpindex words = ',kjpindexnstm
798	STOP 'hydrol_init'
799	END IF
800
801	ALLOCATE (tmc_litter_awet(kjpindex,nstm),stat=ier)
802	IF (ier.NE.0) THEN
803	WRITE (numout,) ' error in tmc_litter_awet allocation. We stop. We need kjpindex words = ',kjpindexnstm
804	STOP 'hydrol_init'
805	END IF
806
807	ALLOCATE (tmc_litter_adry(kjpindex,nstm),stat=ier)
808	IF (ier.NE.0) THEN
809	WRITE (numout,) ' error in tmc_litter_adry allocation. We stop. We need kjpindex words = ',kjpindexnstm
810	STOP 'hydrol_init'
811	END IF
812
813	ALLOCATE (tmc_litt_wet_mea(kjpindex),stat=ier)
814	IF (ier.NE.0) THEN
815	WRITE (numout,*) ' error in tmc_litt_wet_mea allocation. We stop. We need kjpindex words = ',kjpindex
816	STOP 'hydrol_init'
817	END IF
818
819	ALLOCATE (tmc_litt_dry_mea(kjpindex),stat=ier)
820	IF (ier.NE.0) THEN
821	WRITE (numout,*) ' error in tmc_litt_dry_mea allocation. We stop. We need kjpindex words = ',kjpindex
822	STOP 'hydrol_init'
823	END IF
824
825	ALLOCATE (v1(kjpindex,nstm),stat=ier)
826	IF (ier.NE.0) THEN
827	WRITE (numout,) ' error in v1 allocation. We stop. We need kjpindex words = ',kjpindexnstm
828	STOP 'hydrol_init'
829	END IF
830
831	ALLOCATE (vB(kjpindex,nstm),stat=ier)
832	IF (ier.NE.0) THEN
833	WRITE (numout,) ' error in vB allocation. We stop. We need kjpindex words = ',kjpindexnstm
834	STOP 'hydrol_init'
835	END IF
836
837	ALLOCATE (qflux00(kjpindex,nstm),stat=ier)
838	IF (ier.NE.0) THEN
839	WRITE (numout,) ' error in qflux00 allocation. We stop. We need kjpindex words = ',kjpindexnstm
840	STOP 'hydrol_init'
841	END IF
842
843	ALLOCATE (ru_ns(kjpindex,nstm),stat=ier)
844	IF (ier.NE.0) THEN
845	WRITE (numout,) ' error in ru_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
846	STOP 'hydrol_init'
847	END IF
848	ru_ns(:,:) = zero
849
850	ALLOCATE (dr_ns(kjpindex,nstm),stat=ier)
851	IF (ier.NE.0) THEN
852	WRITE (numout,) ' error in dr_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
853	STOP 'hydrol_init'
854	END IF
855	dr_ns(:,:) = zero
856
857	ALLOCATE (tr_ns(kjpindex,nstm),stat=ier)
858	IF (ier.NE.0) THEN
859	WRITE (numout,) ' error in tr_ns allocation. We stop. We need kjpindex words = ',kjpindexnstm
860	STOP 'hydrol_init'
861	END IF
862
863	ALLOCATE (corr_veg_soil(kjpindex,nvm,nstm),stat=ier)
864	IF (ier.NE.0) THEN
865	WRITE (numout,) ' error in corr_veg_soil allocation. We stop. We need kjpindex words = ',kjpindexnvm*nstm
866	STOP 'hydrol_init'
867	END IF
868
869	ALLOCATE (corr_veg_soil_max(kjpindex,nvm,nstm),stat=ier)
870	IF (ier.NE.0) THEN
871	WRITE (numout,) ' error in corr_veg_soil_max allocation. We stop. We need kjpindex words = ',kjpindexnvm*nstm
872	STOP 'hydrol_init'
873	END IF
874
875
876	ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier)
877	IF (ier.NE.0) THEN
878	WRITE (numout,) ' error in mc allocation. We stop. We need kjpindex words = ',kjpindexnslm*nstm
879	STOP 'hydrol_init'
880	END IF
881
882	ALLOCATE (soilmoist(kjpindex,nslm),stat=ier)
883	IF (ier.NE.0) THEN
884	WRITE (numout,) ' error in soilmoist allocation. We stop. We need kjpindex words = ',kjpindexnslm
885	STOP 'hydrol_init'
886	END IF
887
888	ALLOCATE (soil_wet(kjpindex,nslm,nstm),stat=ier)
889	IF (ier.NE.0) THEN
890	WRITE (numout,) ' error in soil_wet allocation. We stop. We need kjpindex words = ',kjpindexnslm*nstm
891	STOP 'hydrol_init'
892	END IF
893
894	ALLOCATE (soil_wet_litter(kjpindex,nstm),stat=ier)
895	IF (ier.NE.0) THEN
896	WRITE (numout,) ' error in soil_wet allocation. We stop. We need kjpindex words = ',kjpindexnstm
897	STOP 'hydrol_init'
898	END IF
899
900	ALLOCATE (qflux(kjpindex,nslm,nstm),stat=ier)
901	IF (ier.NE.0) THEN
902	WRITE (numout,) ' error in qflux allocation. We stop. We need kjpindex words = ',kjpindexnslm*nstm
903	STOP 'hydrol_init'
904	END IF
905
906	ALLOCATE (tmat(kjpindex,nslm,3),stat=ier)
907	IF (ier.NE.0) THEN
908	WRITE (numout,) ' error in tmat allocation. We stop. We need kjpindex words = ',kjpindexnslm*trois
909	STOP 'hydrol_init'
910	END IF
911
912	ALLOCATE (stmat(kjpindex,nslm,3),stat=ier)
913	IF (ier.NE.0) THEN
914	WRITE (numout,) ' error in stmat allocation. We stop. We need kjpindex words = ',kjpindexnslm*trois
915	STOP 'hydrol_init'
916	END IF
917
918	!
919	! If we check the water balance we need two more variables
920	!
921	IF ( check_waterbal ) THEN
922
923	ALLOCATE (tot_water_beg(kjpindex),stat=ier)
924	IF (ier.NE.0) THEN
925	WRITE (numout,*) ' error in tot_water_beg allocation. We stop. We need kjpindex words = ',kjpindex
926	STOP 'hydrol_init'
927	END IF
928
929	ALLOCATE (tot_water_end(kjpindex),stat=ier)
930	IF (ier.NE.0) THEN
931	WRITE (numout,*) ' error in tot_water_end allocation. We stop. We need kjpindex words = ',kjpindex
932	STOP 'hydrol_init'
933	END IF
934
935	ENDIF
936	!
937	! If we use the almaoutputs we need four more variables
938	!
939	IF ( almaoutput ) THEN
940
941	ALLOCATE (tot_watveg_beg(kjpindex),stat=ier)
942	IF (ier.NE.0) THEN
943	WRITE (numout,*) ' error in tot_watveg_beg allocation. We stop. We need kjpindex words = ',kjpindex
944	STOP 'hydrol_init'
945	END IF
946
947	ALLOCATE (tot_watveg_end(kjpindex),stat=ier)
948	IF (ier.NE.0) THEN
949	WRITE (numout,*) ' error in tot_watveg_end allocation. We stop. We need kjpindex words = ',kjpindex
950	STOP 'hydrol_init'
951	END IF
952
953	ALLOCATE (tot_watsoil_beg(kjpindex),stat=ier)
954	IF (ier.NE.0) THEN
955	WRITE (numout,*) ' error in tot_watsoil_beg allocation. We stop. We need kjpindex words = ',kjpindex
956	STOP 'hydrol_init'
957	END IF
958
959	ALLOCATE (tot_watsoil_end(kjpindex),stat=ier)
960	IF (ier.NE.0) THEN
961	WRITE (numout,*) ' error in tot_watsoil_end allocation. We stop. We need kjpindex words = ',kjpindex
962	STOP 'hydrol_init'
963	END IF
964
965	ALLOCATE (delsoilmoist(kjpindex),stat=ier)
966	IF (ier.NE.0) THEN
967	WRITE (numout,*) ' error in delsoilmoist allocation. We stop. We need kjpindex words = ',kjpindex
968	STOP 'hydrol_init'
969	END IF
970
971	ALLOCATE (delintercept(kjpindex),stat=ier)
972	IF (ier.NE.0) THEN
973	WRITE (numout,*) ' error in delintercept. We stop. We need kjpindex words = ',kjpindex
974	STOP 'hydrol_init'
975	END IF
976
977	ALLOCATE (delswe(kjpindex),stat=ier)
978	IF (ier.NE.0) THEN
979	WRITE (numout,*) ' error in delswe. We stop. We need kjpindex words = ',kjpindex
980	STOP 'hydrol_init'
981	ENDIF
982
983	ALLOCATE (snow_beg(kjpindex),stat=ier)
984	IF (ier.NE.0) THEN
985	WRITE (numout,*) ' error in snow_beg allocation. We stop. We need kjpindex words =',kjpindex
986	STOP 'hydrol_init'
987	END IF
988
989	ALLOCATE (snow_end(kjpindex),stat=ier)
990	IF (ier.NE.0) THEN
991	WRITE (numout,*) ' error in snow_end allocation. We stop. We need kjpindex words =',kjpindex
992	STOP 'hydrol_init'
993	END IF
994
995	ENDIF
996
997	! open restart input file done by sechiba_init
998	! and read data from restart input file for HYDROLOGIC process
999
1000	IF (ldrestart_read) THEN
1001
1002	IF (long_print) WRITE (numout,*) ' we have to read a restart file for HYDROLOGIC variables'
1003
1004	IF (is_root_prc) CALL ioconf_setatt('UNITS', '-')
1005	!
1006	DO jst=1,nstm
1007	! var_name= "mc_1" ... "mc_3"
1008	WRITE (var_name,"('moistc_',I1)") jst
1009	CALL ioconf_setatt('LONG_NAME',var_name)
1010	CALL restget_p (rest_id, var_name, nbp_glo, nslm , 1, kjit, .TRUE., mc(:,:,jst), "gather", nbp_glo, index_g)
1011	END DO
1012	!
1013	CALL ioconf_setatt('UNITS', '-')
1014	DO jst=1,nstm
1015	DO jsl=1,nslm
1016	! var_name= "us_1_01" ... "us_3_11"
1017	WRITE (var_name,"('us_',i1,'_',i2.2)") jst,jsl
1018	CALL ioconf_setatt('LONG_NAME',var_name)
1019	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., us(:,:,jst,jsl), "gather", nbp_glo, index_g)
1020	END DO
1021	END DO
1022	!
1023	var_name= 'free_drain_coef'
1024	CALL ioconf_setatt('UNITS', '-')
1025	CALL ioconf_setatt('LONG_NAME','Coefficient for free drainage at bottom of soil')
1026	CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., free_drain_coef, "gather", nbp_glo, index_g)
1027	!
1028	var_name= 'ae_ns'
1029	CALL ioconf_setatt('UNITS', 'kg/m^2')
1030	CALL ioconf_setatt('LONG_NAME','Bare soil evap on each soil type')
1031	CALL restget_p (rest_id, var_name, nbp_glo, nstm, 1, kjit, .TRUE., ae_ns, "gather", nbp_glo, index_g)
1032	!
1033	var_name= 'snow'
1034	CALL ioconf_setatt('UNITS', 'kg/m^2')
1035	CALL ioconf_setatt('LONG_NAME','Snow mass')
1036	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow, "gather", nbp_glo, index_g)
1037	!
1038	var_name= 'snow_age'
1039	CALL ioconf_setatt('UNITS', 'd')
1040	CALL ioconf_setatt('LONG_NAME','Snow age')
1041	CALL restget_p (rest_id, var_name, nbp_glo, 1 , 1, kjit, .TRUE., snow_age, "gather", nbp_glo, index_g)
1042	!
1043	var_name= 'snow_nobio'
1044	CALL ioconf_setatt('UNITS', 'kg/m^2')
1045	CALL ioconf_setatt('LONG_NAME','Snow on other surface types')
1046	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio, "gather", nbp_glo, index_g)
1047	!
1048	var_name= 'snow_nobio_age'
1049	CALL ioconf_setatt('UNITS', 'd')
1050	CALL ioconf_setatt('LONG_NAME','Snow age on other surface types')
1051	CALL restget_p (rest_id, var_name, nbp_glo, nnobio , 1, kjit, .TRUE., snow_nobio_age, "gather", nbp_glo, index_g)
1052	!
1053	var_name= 'vegstress'
1054	CALL ioconf_setatt('UNITS', '-')
1055	CALL ioconf_setatt('LONG_NAME','Vegetation growth moisture stress')
1056	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., vegstress, "gather", nbp_glo, index_g)
1057	!
1058	var_name= 'qsintveg'
1059	CALL ioconf_setatt('UNITS', 'kg/m^2')
1060	CALL ioconf_setatt('LONG_NAME','Intercepted moisture')
1061	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., qsintveg, "gather", nbp_glo, index_g)
1062	!
1063	var_name= 'resdist'
1064	CALL ioconf_setatt('UNITS', '-')
1065	CALL ioconf_setatt('LONG_NAME','Distribution of reservoirs')
1066	CALL restget_p (rest_id, var_name, nbp_glo, nvm, 1, kjit, .TRUE., resdist, "gather", nbp_glo, index_g)
1067	!
1068	!
1069	!
1070	! get restart values if non were found in the restart file
1071	!
1072	!Config Key = HYDROL_MOISTURE_CONTENT
1073	!Config Desc = Soil moisture on each soil tile and levels
1074	!Config If = OK_CWRR
1075	!Config Def = 0.3
1076	!Config Help = The initial value of mc if its value is not found
1077	!Config in the restart file. This should only be used if the model is
1078	!Config started without a restart file.
1079	!Config Units =
1080	!
1081	CALL setvar_p (mc, val_exp, 'HYDROL_MOISTURE_CONTENT', 0.3_r_std)
1082	!
1083	!Config Key = US_INIT
1084	!Config Desc = US_NVM_NSTM_NSLM
1085	!Config If = OK_CWRR
1086	!Config Def = 0.0
1087	!Config Help = The initial value of us if its value is not found
1088	!Config in the restart file. This should only be used if the model is
1089	!Config started without a restart file.
1090	!Config Units =
1091	!
1092	DO jsl=1,nslm
1093	CALL setvar_p (us(:,:,:,jsl), val_exp, 'US_INIT', zero)
1094	ENDDO
1095	!
1096	!Config Key = FREE_DRAIN_COEF
1097	!Config Desc = Coefficient for free drainage at bottom
1098	!Config If = OK_CWRR
1099	!Config Def = 1.0, 1.0, 1.0
1100	!Config Help = The initial value of free drainage if its value is not found
1101	!Config in the restart file. This should only be used if the model is
1102	!Config started without a restart file.
1103	!Config Units =
1104	!
1105	CALL setvar_p (free_drain_coef, val_exp, 'FREE_DRAIN_COEF', free_drain_max)
1106	!
1107	!Config Key = EVAPNU_SOIL
1108	!Config Desc = Bare soil evap on each soil if not found in restart
1109	!Config If = OK_CWRR
1110	!Config Def = 0.0
1111	!Config Help = The initial value of bare soils evap if its value is not found
1112	!Config in the restart file. This should only be used if the model is
1113	!Config started without a restart file.
1114	!Config Units =
1115	!
1116	CALL setvar_p (ae_ns, val_exp, 'EVAPNU_SOIL', zero)
1117	!
1118	!Config Key = HYDROL_SNOW
1119	!Config Desc = Initial snow mass if not found in restart
1120	!Config If = OK_SECHIBA
1121	!Config Def = 0.0
1122	!Config Help = The initial value of snow mass if its value is not found
1123	!Config in the restart file. This should only be used if the model is
1124	!Config started without a restart file.
1125	!Config Units =
1126	!
1127	CALL setvar_p (snow, val_exp, 'HYDROL_SNOW', zero)
1128	!
1129	!Config Key = HYDROL_SNOWAGE
1130	!Config Desc = Initial snow age if not found in restart
1131	!Config If = OK_SECHIBA
1132	!Config Def = 0.0
1133	!Config Help = The initial value of snow age if its value is not found
1134	!Config in the restart file. This should only be used if the model is
1135	!Config started without a restart file.
1136	!Config Units =
1137	!
1138	CALL setvar_p (snow_age, val_exp, 'HYDROL_SNOWAGE', zero)
1139	!
1140	!Config Key = HYDROL_SNOW_NOBIO
1141	!Config Desc = Initial snow amount on ice, lakes, etc. if not found in restart
1142	!Config If = OK_SECHIBA
1143	!Config Def = 0.0
1144	!Config Help = The initial value of snow if its value is not found
1145	!Config in the restart file. This should only be used if the model is
1146	!Config started without a restart file.
1147	!Config Units =
1148	!
1149	CALL setvar_p (snow_nobio, val_exp, 'HYDROL_SNOW_NOBIO', zero)
1150	!
1151	!Config Key = HYDROL_SNOW_NOBIO_AGE
1152	!Config Desc = Initial snow age on ice, lakes, etc. if not found in restart
1153	!Config If = OK_SECHIBA
1154	!Config Def = 0.0
1155	!Config Help = The initial value of snow age if its value is not found
1156	!Config in the restart file. This should only be used if the model is
1157	!Config started without a restart file.
1158	!Config Units =
1159	!
1160	CALL setvar_p (snow_nobio_age, val_exp, 'HYDROL_SNOW_NOBIO_AGE', zero)
1161	!
1162	!
1163	!
1164	!Config Key = HYDROL_QSV
1165	!Config Desc = Initial water on canopy if not found in restart
1166	!Config If = OK_SECHIBA
1167	!Config Def = 0.0
1168	!Config Help = The initial value of moisture on canopy if its value
1169	!Config is not found in the restart file. This should only be used if
1170	!Config the model is started without a restart file.
1171	!Config Units =
1172	!
1173	CALL setvar_p (qsintveg, val_exp, 'HYDROL_QSV', zero)
1174	!
1175	! There is no need to configure the initialisation of resdist. If not available it is the vegetation map
1176	!
1177	IF ( MINVAL(resdist) .EQ. MAXVAL(resdist) .AND. MINVAL(resdist) .EQ. val_exp) THEN
1178	resdist = veget
1179	ENDIF
1180	!
1181	! Remember that it is only frac_nobio + SUM(veget(,:)) that is equal to 1. Thus we need vegtot
1182	!
1183	DO ji = 1, kjpindex
1184	vegtot(ji) = SUM(veget(ji,:))
1185	ENDDO
1186	!
1187	!
1188	! compute the masks for veget
1189
1190	! mask_veget(:,:) = MIN( un, MAX(zero,veget(:,:)))
1191	! mask_soiltype(:,:) = MIN( un, MAX(zero,soiltype(:,:)))
1192	! mask_corr_veg_soil(:,:,:) = MIN( un, MAX(zero,corr_veg_soil(:,:,:)))
1193
1194	mask_veget(:,:) = 0
1195	mask_soiltype(:,:) = 0
1196	mask_corr_veg_soil(:,:,:) = 0
1197
1198	mask_return(:) = 0
1199	index_nsat(:,:) = 0
1200	index_sat(:,:) = 0
1201	n_nsat(:) = 1
1202	n_sat(:) = 0
1203	nslme(:,:) = nslm
1204
1205	DO ji = 1, kjpindex
1206
1207	DO jst = 1, nstm
1208	IF(soiltype(ji,jst) .GT. min_sechiba) THEN
1209	mask_soiltype(ji,jst) = 1
1210	ENDIF
1211	END DO
1212
1213	DO jv = 1, nvm
1214	IF(veget(ji,jv) .GT. min_sechiba) THEN
1215	mask_veget(ji,jv) = 1
1216	ENDIF
1217
1218	DO jst = 1, nstm
1219	IF(corr_veg_soil(ji,jv,jst) .GT. min_sechiba) THEN
1220	mask_corr_veg_soil(ji,jv,jst) = 1
1221	ENDIF
1222	END DO
1223	END DO
1224
1225	! WRITE(numout,*) 'mask: soiltype,mask_soiltype',soiltype(ji,:),mask_soiltype(ji,:)
1226
1227	END DO
1228	! set humrelv from us
1229
1230	humrelv(:,:,:) = SUM(us,dim=4)
1231	vegstressv(:,:,:) = SUM(us,dim=4)
1232	! set humrel from humrelv
1233
1234	humrel(:,:) = zero
1235
1236	DO jst=1,nstm
1237	DO jv=1,nvm
1238	DO ji=1,kjpindex
1239
1240	vegstress(ji,jv)=vegstress(ji,jv) + vegstressv(ji,jv,jst) * soiltype(ji,jst)
1241
1242	! IF(veget(ji,jv).NE.zero) THEN
1243	humrel(ji,jv)=humrel(ji,jv) + humrelv(ji,jv,jst) * &
1244	& soiltype(ji,jst)
1245	humrel(ji,jv)=MAX(humrel(ji,jv), zero)* mask_veget(ji,jv)
1246	! ELSE
1247	! humrel(ji,jv)= zero
1248	! ENDIF
1249	END DO
1250	END DO
1251	END DO
1252	! vegstress(:,:)=humrel(:,:)
1253	ENDIF
1254	!
1255	!
1256	IF (long_print) WRITE (numout,*) ' hydrol_init done '
1257	!
1258	END SUBROUTINE hydrol_init
1259	!
1260	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1261	!
1262	SUBROUTINE hydrol_clear()
1263
1264	l_first_hydrol=.TRUE.
1265	IF (ALLOCATED (mask_veget)) DEALLOCATE (mask_veget)
1266	IF (ALLOCATED (mask_soiltype)) DEALLOCATE (mask_soiltype)
1267	IF (ALLOCATED (mask_corr_veg_soil)) DEALLOCATE (mask_corr_veg_soil)
1268	IF (ALLOCATED (mask_return)) DEALLOCATE (mask_return)
1269	IF (ALLOCATED (index_nsat)) DEALLOCATE (index_nsat)
1270	IF (ALLOCATED (index_sat)) DEALLOCATE (index_sat)
1271	IF (ALLOCATED (n_nsat)) DEALLOCATE (n_nsat)
1272	IF (ALLOCATED (n_sat)) DEALLOCATE (n_sat)
1273	IF (ALLOCATED (nslme)) DEALLOCATE (nslme)
1274	IF (ALLOCATED (humrelv)) DEALLOCATE (humrelv)
1275	IF (ALLOCATED (vegstressv)) DEALLOCATE (vegstressv)
1276	IF (ALLOCATED (us)) DEALLOCATE (us)
1277	IF (ALLOCATED (precisol)) DEALLOCATE (precisol)
1278	IF (ALLOCATED (precisol_ns)) DEALLOCATE (precisol_ns)
1279	IF (ALLOCATED (free_drain_coef)) DEALLOCATE (free_drain_coef)
1280	IF (ALLOCATED (ae_ns)) DEALLOCATE (ae_ns)
1281	IF (ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
1282	IF (ALLOCATED (rootsink)) DEALLOCATE (rootsink)
1283	IF (ALLOCATED (subsnowveg)) DEALLOCATE (subsnowveg)
1284	IF (ALLOCATED (subsnownobio)) DEALLOCATE (subsnownobio)
1285	IF (ALLOCATED (snowmelt)) DEALLOCATE (snowmelt)
1286	IF (ALLOCATED (icemelt)) DEALLOCATE (icemelt)
1287	IF (ALLOCATED (subsinksoil)) DEALLOCATE (subsinksoil)
1288	IF (ALLOCATED (mx_eau_var)) DEALLOCATE (mx_eau_var)
1289	IF (ALLOCATED (vegtot)) DEALLOCATE (vegtot)
1290	IF (ALLOCATED (resdist)) DEALLOCATE (resdist)
1291	IF (ALLOCATED (tot_water_beg)) DEALLOCATE (tot_water_beg)
1292	IF (ALLOCATED (tot_water_end)) DEALLOCATE (tot_water_end)
1293	IF (ALLOCATED (tot_watveg_beg)) DEALLOCATE (tot_watveg_beg)
1294	IF (ALLOCATED (tot_watveg_end)) DEALLOCATE (tot_watveg_end)
1295	IF (ALLOCATED (tot_watsoil_beg)) DEALLOCATE (tot_watsoil_beg)
1296	IF (ALLOCATED (tot_watsoil_end)) DEALLOCATE (tot_watsoil_end)
1297	IF (ALLOCATED (delsoilmoist)) DEALLOCATE (delsoilmoist)
1298	IF (ALLOCATED (delintercept)) DEALLOCATE (delintercept)
1299	IF (ALLOCATED (snow_beg)) DEALLOCATE (snow_beg)
1300	IF (ALLOCATED (snow_end)) DEALLOCATE (snow_end)
1301	IF (ALLOCATED (delswe)) DEALLOCATE (delswe)
1302	! more allocation for cwrr scheme
1303	IF (ALLOCATED (v1)) DEALLOCATE (v1)
1304	IF (ALLOCATED (vB)) DEALLOCATE (vB)
1305	IF (ALLOCATED (humtot)) DEALLOCATE (humtot)
1306	IF (ALLOCATED (flux)) DEALLOCATE (flux)
1307	IF (ALLOCATED (resolv)) DEALLOCATE (resolv)
1308	IF (ALLOCATED (a)) DEALLOCATE (a)
1309	IF (ALLOCATED (b)) DEALLOCATE (b)
1310	IF (ALLOCATED (d)) DEALLOCATE (d)
1311	IF (ALLOCATED (e)) DEALLOCATE (e)
1312	IF (ALLOCATED (f)) DEALLOCATE (f)
1313	IF (ALLOCATED (g1)) DEALLOCATE (g1)
1314	IF (ALLOCATED (ep)) DEALLOCATE (ep)
1315	IF (ALLOCATED (fp)) DEALLOCATE (fp)
1316	IF (ALLOCATED (gp)) DEALLOCATE (gp)
1317	IF (ALLOCATED (rhs)) DEALLOCATE (rhs)
1318	IF (ALLOCATED (srhs)) DEALLOCATE (srhs)
1319	IF (ALLOCATED (gam)) DEALLOCATE (gam)
1320	IF (ALLOCATED (tmc)) DEALLOCATE (tmc)
1321	IF (ALLOCATED (tmcs)) DEALLOCATE (tmcs)
1322	IF (ALLOCATED (tmc_litter)) DEALLOCATE (tmc_litter)
1323	IF (ALLOCATED (tmc_litt_mea)) DEALLOCATE (tmc_litt_mea)
1324	IF (ALLOCATED (tmc_litter_res)) DEALLOCATE (tmc_litter_res)
1325	IF (ALLOCATED (tmc_litter_wilt)) DEALLOCATE (tmc_litter_wilt)
1326	IF (ALLOCATED (tmc_litter_field)) DEALLOCATE (tmc_litter_field)
1327	IF (ALLOCATED (tmc_litter_sat)) DEALLOCATE (tmc_litter_sat)
1328	IF (ALLOCATED (tmc_litter_awet)) DEALLOCATE (tmc_litter_awet)
1329	IF (ALLOCATED (tmc_litter_adry)) DEALLOCATE (tmc_litter_adry)
1330	IF (ALLOCATED (tmc_litt_wet_mea)) DEALLOCATE (tmc_litt_wet_mea)
1331	IF (ALLOCATED (tmc_litt_dry_mea)) DEALLOCATE (tmc_litt_dry_mea)
1332	IF (ALLOCATED (qflux00)) DEALLOCATE (qflux00)
1333	IF (ALLOCATED (ru_ns)) DEALLOCATE (ru_ns)
1334	IF (ALLOCATED (dr_ns)) DEALLOCATE (dr_ns)
1335	IF (ALLOCATED (tr_ns)) DEALLOCATE (tr_ns)
1336	IF (ALLOCATED (corr_veg_soil)) DEALLOCATE (corr_veg_soil)
1337	IF (ALLOCATED (corr_veg_soil_max)) DEALLOCATE (corr_veg_soil_max)
1338	IF (ALLOCATED (mc)) DEALLOCATE (mc)
1339	IF (ALLOCATED (soilmoist)) DEALLOCATE (soilmoist)
1340	IF (ALLOCATED (soil_wet)) DEALLOCATE (soil_wet)
1341	IF (ALLOCATED (soil_wet_litter)) DEALLOCATE (soil_wet_litter)
1342	IF (ALLOCATED (qflux)) DEALLOCATE (qflux)
1343	IF (ALLOCATED (tmat)) DEALLOCATE (tmat)
1344	IF (ALLOCATED (stmat)) DEALLOCATE (stmat)
1345
1346	RETURN
1347
1348	END SUBROUTINE hydrol_clear
1349
1350	!! This routine initializes HYDROLOGIC variables
1351	!! - mx_eau_var
1352
1353	SUBROUTINE hydrol_var_init (kjpindex, veget, soiltype, mx_eau_var, shumdiag, litterhumdiag, drysoil_frac, evap_bare_lim)
1354
1355	! interface description
1356	! input scalar
1357	INTEGER(i_std), INTENT(in) :: kjpindex !! domain size
1358	! input fields
1359	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! fraction of vegetation type
1360	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
1361	! output fields
1362	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: mx_eau_var !!
1363	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! relative soil moisture
1364	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
1365	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! function of litter humidity
1366	REAL(r_std),DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !!
1367	! local declaration
1368	REAL(r_std), DIMENSION (kjpindex) :: dpu_mean !! mean soil depth
1369
1370	INTEGER(i_std) :: ji, jv, jd, jst, jsl
1371	REAL(r_std) :: m, frac
1372	!
1373	! initialisation
1374	mx_eau_var(:) = zero
1375	dpu_mean(:)= zero
1376	!
1377	DO ji = 1,kjpindex
1378	DO jst = 1,nstm
1379	dpu_mean(ji)=dpu_mean(ji)+dpu(jst)*soiltype(ji,jst)
1380	END DO
1381	END DO
1382	!
1383	DO ji = 1,kjpindex
1384	DO jst = 1,nstm
1385	mx_eau_var(ji) = mx_eau_var(ji) + soiltype(ji,jst)*&
1386	& dpu(jst)millemcs(jst)
1387	END DO
1388	END DO
1389
1390	DO ji = 1,kjpindex
1391	IF (vegtot(ji) .LE. zero) THEN
1392	mx_eau_var(ji) = mx_eau_eau*deux
1393	ENDIF
1394
1395	END DO
1396
1397
1398	!
1399	! Calcul the matrix coef for dublin model:
1400	! pice-wise linearised hydraulic conductivity k_lin=alin * mc_lin + b_lin
1401	! and diffusivity d_lin in each interval of mc, called mc_lin,
1402	! between imin, for residual mcr,
1403	! and imax for saturation mcs.
1404	!
1405	DO jst=1,nstm
1406	m = un - un / nvan(jst)
1407	! WRITE(numout,*) 'jst',jst,imin,imax
1408	mc_lin(imin,jst)=mcr(jst)
1409	mc_lin(imax,jst)=mcs(jst)
1410	tmcs(jst)=dpu(jst)* mille*mcs(jst)
1411	zz(1,jst) = zero
1412	dz(1,jst) = zero
1413	DO jsl=2,nslm
1414	zz(jsl,jst) = dpu(jst)* mille((2(jsl-1))-1)/ ((2*(nslm-1))-1)
1415	dz(jsl,jst) = zz(jsl,jst)-zz(jsl-1,jst)
1416	! WRITE(numout,*) 'jsl, zz,dz',jsl, dpu(jst),zz(jsl,jst),dz(jsl,jst)
1417	ENDDO
1418	zz(nslm+1,jst) = zz(nslm,jst)
1419	dz(nslm+1,jst) = zero
1420	DO ji= imin+1, imax-1
1421	mc_lin(ji,jst) = mcr(jst) + (ji-imin)*(mcs(jst)-mcr(jst))/(imax-imin)
1422	ENDDO
1423	DO ji = imin,imax-1
1424	frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst)))
1425	k_lin(ji,jst) = ks(jst) * (frac*0.5) ( un - ( un - frac (un/m)) m )**2
1426	frac=MIN(un,(mc_lin(ji+1,jst)-mcr(jst))/(mcs(jst)-mcr(jst)))
1427	k_lin(ji+1,jst) = ks(jst) * (frac*0.5) ( un - ( un - frac (un/m)) m )**2
1428	a_lin(ji,jst) = (k_lin(ji+1,jst)-k_lin(ji,jst)) / (mc_lin(ji+1,jst)-mc_lin(ji,jst))
1429	b_lin(ji,jst) = k_lin(ji,jst) - a_lin(ji,jst)*mc_lin(ji,jst)
1430	!- Il faudrait ici definir a et b pour mc > mcs, et mc < mcr car c'est un cas auquel on peut etre confronte... a reflechir
1431
1432	IF(ji.NE.imin.AND.ji.NE.imax-1) THEN
1433	frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst)))
1434	d_lin(ji,jst) =(k_lin(ji,jst) / (avan(jst)mnvan(jst))) * &
1435	& ( (frac*(-un/m))/(mc_lin(ji,jst)-mcr(jst)) ) &
1436	& ( frac(-un/m) -un ) (-m)
1437	frac=MIN(un,(mc_lin(ji+1,jst)-mcr(jst))/(mcs(jst)-mcr(jst)))
1438	d_lin(ji+1,jst) =(k_lin(ji+1,jst) / (avan(jst)mnvan(jst)))*&
1439	& ( (frac*(-un/m))/(mc_lin(ji+1,jst)-mcr(jst)) ) &
1440	& ( frac(-un/m) -un ) (-m)
1441	d_lin(ji,jst) = undemi * (d_lin(ji,jst)+d_lin(ji+1,jst))
1442	ELSEIF(ji.EQ.imin) THEN
1443	d_lin(ji,jst) = zero
1444	ELSEIF(ji.EQ.imax-1) THEN
1445	frac=MIN(un,(mc_lin(ji,jst)-mcr(jst))/(mcs(jst)-mcr(jst)))
1446	d_lin(ji,jst) =(k_lin(ji,jst) / (avan(jst)mnvan(jst))) * &
1447	& ( (frac*(-un/m))/(mc_lin(ji,jst)-mcr(jst)) ) &
1448	& ( frac(-un/m) -un ) (-m)
1449	ENDIF
1450	ENDDO
1451	ENDDO
1452
1453
1454
1455	! Compute the litter humidity, shumdiag and fry
1456
1457	litterhumdiag(:) = zero
1458	tmc_litt_mea(:) = zero
1459	tmc_litt_wet_mea(:) = zero
1460	tmc_litt_dry_mea(:) = zero
1461	shumdiag(:,:) = zero
1462	soilmoist(:,:) = zero
1463	humtot(:) = zero
1464	tmc(:,:) = zero
1465
1466	! Loop on soil types to compute the variables (ji,jst)
1467
1468	DO jst=1,nstm
1469
1470	! the residual 1st layer soil moisture:
1471	v1r(jst) = dz(2,jst)*mcr(jst)/deux
1472
1473	! the saturated Bottom layer soil moisture:
1474	! v A CALCULER SUR TOUT LE PROFIL (vs et vr egalement)
1475	vBs(jst) = dz(nslm,jst)*mcs(jst)/deux
1476
1477	! The total soil moisture for each soil type:
1478
1479	DO ji=1,kjpindex
1480	tmc(ji,jst)= dz(2,jst) * ( trois*mc(ji,1,jst)+ mc(ji,2,jst))/huit
1481	END DO
1482
1483	DO jsl=2,nslm-1
1484	DO ji=1,kjpindex
1485	tmc(ji,jst) = tmc(ji,jst) + dz(jsl,jst) * ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
1486	& + dz(jsl+1,jst)(troismc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
1487	END DO
1488	END DO
1489
1490	DO ji=1,kjpindex
1491	tmc(ji,jst) = tmc(ji,jst) + dz(nslm,jst) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
1492	END DO
1493
1494
1495	! The litter variables:
1496
1497	DO ji=1,kjpindex
1498	tmc_litter(ji,jst) = dz(2,jst) * (trois*mc(ji,1,jst)+mc(ji,2,jst))/huit
1499	tmc_litter_wilt(ji,jst) = dz(2,jst) * mcw(jst) / deux
1500	tmc_litter_res(ji,jst) = dz(2,jst) * mcr(jst) / deux
1501	tmc_litter_field(ji,jst) = dz(2,jst) * mcf(jst) / deux
1502	tmc_litter_sat(ji,jst) = dz(2,jst) * mcs(jst) / deux
1503	tmc_litter_awet(ji,jst) = dz(2,jst) * mc_awet(jst) / deux
1504	tmc_litter_adry(ji,jst) = dz(2,jst) * mc_adry(jst) / deux
1505	END DO
1506
1507
1508	! sum from level 1 to 4
1509
1510	DO jsl=2,4
1511
1512	DO ji=1,kjpindex
1513	tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl,jst) * &
1514	& ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
1515	& + dz(jsl+1,jst)(troismc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
1516
1517	tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + &
1518	&(dz(jsl,jst)+ dz(jsl+1,jst))*&
1519	& mcw(jst)/deux
1520	tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + &
1521	&(dz(jsl,jst)+ dz(jsl+1,jst))*&
1522	& mcr(jst)/deux
1523	tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + &
1524	&(dz(jsl,jst)+ dz(jsl+1,jst))* &
1525	& mcs(jst)/deux
1526	tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + &
1527	& (dz(jsl,jst)+ dz(jsl+1,jst))* &
1528	& mcf(jst)/deux
1529	tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + &
1530	&(dz(jsl,jst)+ dz(jsl+1,jst))* &
1531	& mc_awet(jst)/deux
1532	tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + &
1533	& (dz(jsl,jst)+ dz(jsl+1,jst))* &
1534	& mc_adry(jst)/deux
1535	END DO
1536
1537	END DO
1538
1539
1540	! subsequent calcul of soil_wet_litter (tmc-tmcw)/(tmcf-tmcw)
1541
1542	DO ji=1,kjpindex
1543	soil_wet_litter(ji,jst)=MIN(un, MAX(zero,&
1544	&(tmc_litter(ji,jst)-tmc_litter_wilt(ji,jst))/&
1545	& (tmc_litter_field(ji,jst)-tmc_litter_wilt(ji,jst)) ))
1546	END DO
1547
1548	! Soil wetness profiles (mc-mcw)/(mcs-mcw)
1549
1550	DO ji=1,kjpindex
1551	soil_wet(ji,1,jst) = MIN(un, MAX(zero,&
1552	&(troismc(ji,1,jst) + mc(ji,2,jst) - quatremcw(jst))&
1553	& /(quatre*(mcs(jst)-mcw(jst))) ))
1554	humrelv(ji,1,jst) = zero
1555
1556	END DO
1557
1558	DO jsl=2,nslm-1
1559	DO ji=1,kjpindex
1560	soil_wet(ji,jsl,jst) = MIN(un, MAX(zero,&
1561	& (trois*mc(ji,jsl,jst) + &
1562	& mc(ji,jsl-1,jst) *(dz(jsl,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) &
1563	& + mc(ji,jsl+1,jst)*(dz(jsl+1,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) &
1564	& - quatremcw(jst)) / (quatre(mcs(jst)-mcw(jst))) ))
1565	END DO
1566	END DO
1567
1568	DO ji=1,kjpindex
1569	soil_wet(ji,nslm,jst) = MIN(un, MAX(zero,&
1570	& (trois*mc(ji,nslm,jst) &
1571	& + mc(ji,nslm-1,jst)-quatremcw(jst))/(quatre(mcs(jst)-mcw(jst))) ))
1572	END DO
1573
1574	END DO ! loop on soil type
1575
1576
1577	!Now we compute the grid averaged values:
1578
1579	DO jst=1,nstm
1580	DO ji=1,kjpindex
1581
1582	humtot(ji) = humtot(ji) + soiltype(ji,jst) * tmc(ji,jst)
1583
1584	litterhumdiag(ji) = litterhumdiag(ji) + &
1585	& soil_wet_litter(ji,jst) * soiltype(ji,jst)
1586
1587	tmc_litt_wet_mea(ji) = tmc_litt_wet_mea(ji) + &
1588	& tmc_litter_awet(ji,jst)* soiltype(ji,jst)
1589
1590	tmc_litt_dry_mea(ji) = tmc_litt_dry_mea(ji) + &
1591	& tmc_litter_adry(ji,jst) * soiltype(ji,jst)
1592
1593	tmc_litt_mea(ji) = tmc_litt_mea(ji) + &
1594	& tmc_litter(ji,jst) * soiltype(ji,jst)
1595	END DO
1596
1597
1598
1599	DO jsl=1,nbdl
1600	DO ji=1,kjpindex
1601	shumdiag(ji,jsl)= shumdiag(ji,jsl) + soil_wet(ji,jsl,jst) * &
1602	& ((mcs(jst)-mcw(jst))/(mcf(jst)-mcw(jst))) * &
1603	& soiltype(ji,jst)
1604	soilmoist(ji,jsl) = soilmoist(ji,jsl) + mc(ji,jsl,jst)*soiltype(ji,jst)
1605	shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero)
1606	END DO
1607	END DO
1608
1609	END DO ! loop on soiltype
1610	!
1611	!
1612	!
1613	DO ji=1,kjpindex
1614	drysoil_frac(ji) = un + MAX( MIN( (tmc_litt_dry_mea(ji) - tmc_litt_mea(ji)) / &
1615	& (tmc_litt_wet_mea(ji) - tmc_litt_dry_mea(ji)), zero), - un)
1616	END DO
1617
1618	evap_bare_lim = zero
1619
1620	!!$ IF ( COUNT(diaglev .EQ. undef_sechiba) > 0 ) THEN
1621	!!$
1622	!!$ DO jsl=1,nbdl-1
1623	!!$ diaglev(jsl) = zz(jsl,1) + dz(jsl+1,1)/deux
1624	!!$ END DO
1625	!!$ diaglev(nbdl) = zz(nbdl,1)
1626	!!$ interpol_diag = .FALSE.
1627	!!$
1628	!!$ ENDIF
1629
1630	IF (long_print) WRITE (numout,*) ' hydrol_var_init done '
1631
1632	END SUBROUTINE hydrol_var_init
1633
1634	!! This routine computes snow processes
1635	!!
1636	SUBROUTINE hydrol_snow (kjpindex, dtradia, precip_rain, precip_snow , temp_sol_new, soilcap,&
1637	& frac_nobio, totfrac_nobio, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age, &
1638	& tot_melt, snowdepth)
1639
1640	!
1641	! interface description
1642	! input scalar
1643	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1644	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
1645	! input fields
1646	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rainfall
1647	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_snow !! Snow precipitation
1648	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: temp_sol_new !! New soil temperature
1649	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: soilcap !! Soil capacity
1650	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: frac_nobio !! Fraction of continental ice, lakes, ...
1651	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: totfrac_nobio !! Total fraction of continental ice+lakes+ ...
1652	! modified fields
1653	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !! Bare soil evaporation
1654	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapsno !! Snow evaporation
1655	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow !! Snow mass [Kg/m^2]
1656	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: snow_age !! Snow age
1657	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio !! Ice water balance
1658	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(inout) :: snow_nobio_age!! Snow age on ice, lakes, ...
1659	! output fields
1660	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: tot_melt !! Total melt
1661	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: snowdepth !! Snow depth
1662	!
1663	! local declaration
1664	!
1665	INTEGER(i_std) :: ji, jv
1666	REAL(r_std), DIMENSION (kjpindex) :: d_age !! Snow age change
1667	REAL(r_std), DIMENSION (kjpindex) :: xx !! temporary
1668	REAL(r_std) :: snowmelt_tmp !! The name says it all !
1669
1670	!
1671	! for continental points
1672	!
1673
1674	!
1675	! 0. initialisation
1676	!
1677	DO jv = 1, nnobio
1678	DO ji=1,kjpindex
1679	subsnownobio(ji,jv) = zero
1680	ENDDO
1681	ENDDO
1682	DO ji=1,kjpindex
1683	subsnowveg(ji) = zero
1684	snowmelt(ji) = zero
1685	icemelt(ji) = zero
1686	subsinksoil(ji) = zero
1687	tot_melt(ji) = zero
1688	ENDDO
1689	!
1690	! 1. On vegetation
1691	!
1692	DO ji=1,kjpindex
1693	!
1694	! 1.1. It is snowing
1695	!
1696	snow(ji) = snow(ji) + (un - totfrac_nobio(ji))*precip_snow(ji)
1697	!
1698	!
1699	! 1.2. Sublimation - separate between vegetated and no-veget fractions
1700	! Care has to be taken as we might have sublimation from the
1701	! the frac_nobio while there is no snow on the rest of the grid.
1702	!
1703	IF ( snow(ji) > snowcri ) THEN
1704	subsnownobio(ji,iice) = frac_nobio(ji,iice)*vevapsno(ji)
1705	subsnowveg(ji) = vevapsno(ji) - subsnownobio(ji,iice)
1706	ELSE
1707	! Correction Nathalie - Juillet 2006.
1708	! On doit d'abord tester s'il existe un frac_nobio!
1709	! Pour le moment je ne regarde que le iice
1710	IF ( frac_nobio(ji,iice) .GT. min_sechiba) THEN
1711	subsnownobio(ji,iice) = vevapsno(ji)
1712	subsnowveg(ji) = zero
1713	ELSE
1714	subsnownobio(ji,iice) = zero
1715	subsnowveg(ji) = vevapsno(ji)
1716	ENDIF
1717	ENDIF
1718	!
1719	!
1720	! 1.2.1 Check that sublimation on the vegetated fraction is possible.
1721	!
1722	IF (subsnowveg(ji) .GT. snow(ji)) THEN
1723	! What could not be sublimated goes into soil evaporation
1724	! vevapnu(ji) = vevapnu(ji) + (subsnowveg(ji) - snow(ji))
1725	IF( (un - totfrac_nobio(ji)).GT.min_sechiba) THEN
1726	subsinksoil (ji) = (subsnowveg(ji) - snow(ji))/ (un - totfrac_nobio(ji))
1727	END IF
1728	! Sublimation is thus limited to what is available
1729	subsnowveg(ji) = snow(ji)
1730	snow(ji) = zero
1731	vevapsno(ji) = subsnowveg(ji) + subsnownobio(ji,iice)
1732	ELSE
1733	snow(ji) = snow(ji) - subsnowveg(ji)
1734	ENDIF
1735	!
1736	! 1.3. snow melt only if temperature positive
1737	!
1738	IF (temp_sol_new(ji).GT.tp_00) THEN
1739	!
1740	IF (snow(ji).GT.sneige) THEN
1741	!
1742	snowmelt(ji) = (un - frac_nobio(ji,iice))(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1743	!
1744	! 1.3.1.1 enough snow for melting or not
1745	!
1746	IF (snowmelt(ji).LT.snow(ji)) THEN
1747	snow(ji) = snow(ji) - snowmelt(ji)
1748	ELSE
1749	snowmelt(ji) = snow(ji)
1750	snow(ji) = zero
1751	END IF
1752	!
1753	ELSEIF (snow(ji).GE.zero) THEN
1754	!
1755	! 1.3.2 not enough snow
1756	!
1757	snowmelt(ji) = snow(ji)
1758	snow(ji) = zero
1759	ELSE
1760	!
1761	! 1.3.3 negative snow - now snow melt
1762	!
1763	snow(ji) = zero
1764	snowmelt(ji) = zero
1765	WRITE(numout,*) 'hydrol_snow: WARNING! snow was negative and was reset to zero. '
1766	!
1767	END IF
1768
1769	ENDIF
1770	!
1771	! 1.4. Ice melt only if there is more than a given mass : maxmass_glacier,
1772	! i.e. only weight melts glaciers !
1773	! Ajouts Edouard Davin / Nathalie de Noblet add extra to melting
1774	!
1775	IF ( snow(ji) .GT. maxmass_glacier ) THEN
1776	snowmelt(ji) = snowmelt(ji) + (snow(ji) - maxmass_glacier)
1777	snow(ji) = maxmass_glacier
1778	ENDIF
1779	!
1780	END DO
1781	!
1782	! 2. On Land ice
1783	!
1784	DO ji=1,kjpindex
1785	!
1786	! 2.1. It is snowing
1787	!
1788	snow_nobio(ji,iice) = snow_nobio(ji,iice) + frac_nobio(ji,iice)*precip_snow(ji) + &
1789	& frac_nobio(ji,iice)*precip_rain(ji)
1790	!
1791	! 2.2. Sublimation - was calculated before it can give us negative snow_nobio but that is OK
1792	! Once it goes below a certain values (-maxmass_glacier for instance) we should kill
1793	! the frac_nobio(ji,iice) !
1794	!
1795	snow_nobio(ji,iice) = snow_nobio(ji,iice) - subsnownobio(ji,iice)
1796	!
1797	! 2.3. Snow melt only for continental ice fraction
1798	!
1799	snowmelt_tmp = zero
1800	IF (temp_sol_new(ji) .GT. tp_00) THEN
1801	!
1802	! 2.3.1 If there is snow on the ice-fraction it can melt
1803	!
1804	snowmelt_tmp = frac_nobio(ji,iice)(temp_sol_new(ji) - tp_00) soilcap(ji) / chalfu0
1805	!
1806	IF ( snowmelt_tmp .GT. snow_nobio(ji,iice) ) THEN
1807	snowmelt_tmp = MAX( zero, snow_nobio(ji,iice))
1808	ENDIF
1809	snowmelt(ji) = snowmelt(ji) + snowmelt_tmp
1810	snow_nobio(ji,iice) = snow_nobio(ji,iice) - snowmelt_tmp
1811	!
1812	ENDIF
1813	!
1814	! 2.4 Ice melt only if there is more than a given mass : maxmass_glacier,
1815	! i.e. only weight melts glaciers !
1816	!
1817	IF ( snow_nobio(ji,iice) .GT. maxmass_glacier ) THEN
1818	icemelt(ji) = snow_nobio(ji,iice) - maxmass_glacier
1819	snow_nobio(ji,iice) = maxmass_glacier
1820	ENDIF
1821	!
1822	END DO
1823
1824	!
1825	! 3. On other surface types - not done yet
1826	!
1827	IF ( nnobio .GT. 1 ) THEN
1828	WRITE(numout,*) 'WE HAVE',nnobio-1,' SURFACE TYPES I DO NOT KNOW'
1829	WRITE(numout,*) 'CANNOT TREAT SNOW ON THESE SURFACE TYPES'
1830	STOP 'in hydrol_snow'
1831	ENDIF
1832
1833	!
1834	! 4. computes total melt (snow and ice)
1835	!
1836	DO ji = 1, kjpindex
1837	tot_melt(ji) = icemelt(ji) + snowmelt(ji)
1838	ENDDO
1839
1840	!
1841	! 5. computes snow age on veg and ice (for albedo)
1842	!
1843	DO ji = 1, kjpindex
1844	!
1845	! 5.1 Snow age on vegetation
1846	!
1847	IF (snow(ji) .LE. zero) THEN
1848	snow_age(ji) = zero
1849	ELSE
1850	snow_age(ji) =(snow_age(ji) + (un - snow_age(ji)/max_snow_age) * dtradia/one_day) &
1851	& * EXP(-precip_snow(ji) / snow_trans)
1852	ENDIF
1853	!
1854	! 5.2 Snow age on ice
1855	!
1856	! age of snow on ice: a little bit different because in cold regions, we really
1857	! cannot negect the effect of cold temperatures on snow metamorphism any more.
1858	!
1859	IF (snow_nobio(ji,iice) .LE. zero) THEN
1860	snow_nobio_age(ji,iice) = zero
1861	ELSE
1862	!
1863	d_age(ji) = ( snow_nobio_age(ji,iice) + &
1864	& (un - snow_nobio_age(ji,iice)/max_snow_age) * dtradia/one_day ) * &
1865	& EXP(-precip_snow(ji) / snow_trans) - snow_nobio_age(ji,iice)
1866	IF (d_age(ji) .GT. min_sechiba ) THEN
1867	xx(ji) = MAX( tp_00 - temp_sol_new(ji), zero )
1868	xx(ji) = ( xx(ji) / 7._r_std ) ** 4._r_std
1869	d_age(ji) = d_age(ji) / (un+xx(ji))
1870	ENDIF
1871	snow_nobio_age(ji,iice) = MAX( snow_nobio_age(ji,iice) + d_age(ji), zero )
1872	!
1873	ENDIF
1874
1875	ENDDO
1876
1877	!
1878	! 6.0 Diagnose the depth of the snow layer
1879	!
1880
1881	DO ji = 1, kjpindex
1882	snowdepth(ji) = snow(ji) /sn_dens
1883	ENDDO
1884
1885	IF (long_print) WRITE (numout,*) ' hydrol_snow done '
1886
1887	END SUBROUTINE hydrol_snow
1888
1889	!! This routine computes canopy processes
1890	!!
1891	SUBROUTINE hydrol_canop (kjpindex, precip_rain, vevapwet, veget, qsintmax, &
1892	& qsintveg,precisol,tot_melt)
1893
1894	!
1895	! interface description
1896	!
1897	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
1898	! input fields
1899	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: precip_rain !! Rain precipitation
1900	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: vevapwet !! Interception loss
1901	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: veget !! Fraction of vegetation type
1902	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: qsintmax !! Maximum water on vegetation for interception
1903	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt
1904	! modified fields
1905	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(inout) :: qsintveg !! Water on vegetation due to interception
1906	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(out) :: precisol !! Eau tombee sur le sol
1907	! output fields
1908
1909	!
1910	! local declaration
1911	!
1912	INTEGER(i_std) :: ji, jv
1913	REAL(r_std), DIMENSION (kjpindex,nvm) :: zqsintvegnew
1914
1915	! calcul de qsintmax a prevoir a chaque pas de temps
1916	! dans ini_sechiba
1917	! boucle sur les points continentaux
1918	! calcul de qsintveg au pas de temps suivant
1919	! par ajout du flux interception loss
1920	! calcule par enerbil en fonction
1921	! des calculs faits dans diffuco
1922	! calcul de ce qui tombe sur le sol
1923	! avec accumulation dans precisol
1924	! essayer d'harmoniser le traitement du sol nu
1925	! avec celui des differents types de vegetation
1926	! fait si on impose qsintmax ( ,1) = 0.0
1927	!
1928	! loop for continental subdomain
1929	!
1930	!
1931	! 1. evaporation off the continents
1932	!
1933	! 1.1 The interception loss is take off the canopy.
1934	DO jv=1,nvm
1935	qsintveg(:,jv) = qsintveg(:,jv) - vevapwet(:,jv)
1936	END DO
1937
1938	! 1.2 It is raining : precip_rain is shared for each vegetation
1939	! type
1940	! sum (veget (1,nvm)) must be egal to 1-totfrac_nobio.
1941	! iniveget computes veget each day
1942	!
1943	DO jv=1,nvm
1944	! Correction Nathalie - Juin 2006 - une partie de la pluie arrivera toujours sur le sol
1945	! sorte de throughfall supplementaire
1946	!qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * precip_rain(:)
1947	qsintveg(:,jv) = qsintveg(:,jv) + veget(:,jv) * ((1-throughfall_by_pft(jv))*precip_rain(:))
1948	END DO
1949
1950	!
1951	! 1.3 Limits the effect and sum what receives soil
1952	!
1953	precisol(:,:) = zero
1954	DO jv=1,nvm
1955	DO ji = 1, kjpindex
1956	zqsintvegnew(ji,jv) = MIN (qsintveg(ji,jv),qsintmax(ji,jv))
1957	! correction throughfall Nathalie - Juin 2006
1958	!precisol(ji,jv) = qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
1959	precisol(ji,jv) = (veget(ji,jv)throughfall_by_pft(jv)precip_rain(ji)) + qsintveg(ji,jv ) - zqsintvegnew (ji,jv)
1960	ENDDO
1961	END DO
1962	!
1963	DO jv=1,nvm
1964	DO ji = 1, kjpindex
1965	IF (vegtot(ji).GT.min_sechiba) THEN
1966	precisol(ji,jv) = precisol(ji,jv)+tot_melt(ji)*veget(ji,jv)/vegtot(ji)
1967	ENDIF
1968	ENDDO
1969	END DO
1970	!
1971	!
1972	! 1.4 swap qsintveg to the new value
1973	!
1974
1975	DO jv=1,nvm
1976	qsintveg(:,jv) = zqsintvegnew (:,jv)
1977	END DO
1978
1979	IF (long_print) WRITE (numout,*) ' hydrol_canop done '
1980
1981	END SUBROUTINE hydrol_canop
1982	!!
1983	!!
1984	!!
1985	SUBROUTINE hydrol_vegupd(kjpindex, veget, veget_max, soiltype,qsintveg,resdist)
1986	!
1987	! The vegetation cover has changed and we need to adapt the reservoir distribution
1988	! and the distribution of plants on different soil types.
1989	! You may note that this occurs after evaporation and so on have been computed. It is
1990	! not a problem as a new vegetation fraction will start with humrel=0 and thus will have no
1991	! evaporation. If this is not the case it should have been caught above.
1992	!
1993	! input scalar
1994	INTEGER(i_std), INTENT(in) :: kjpindex
1995	! input fields
1996	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! New vegetation map
1997	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. fraction of vegetation type
1998	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types : proportion of each soil type
1999	! modified fields
2000	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: qsintveg !! Water on vegetation
2001	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(inout) :: resdist !! Old vegetation map
2002	!
2003	! local declaration
2004	!
2005	INTEGER(i_std) :: ji,jv,jst,jst_pref
2006	REAL(r_std), DIMENSION (kjpindex,nstm) :: soil_exist,soil_exist_max
2007	REAL(r_std), DIMENSION (kjpindex,nvm) :: veget_exist,veget_exist_max
2008	REAL(r_std), DIMENSION (kjpindex,nvm) :: qsintveg2 !! Water on vegetation due to interception over old veget
2009	REAL(r_std), DIMENSION (kjpindex,nvm) :: vmr !! variation of veget
2010	REAL(r_std), DIMENSION (kjpindex,nvm) :: qsdq
2011	REAL(r_std), DIMENSION(kjpindex) :: vegchtot,vtr, qstr, fra
2012	REAL(r_std), PARAMETER :: EPS1 = EPSILON(un)
2013	!
2014	DO jv = 1, nvm
2015	DO ji = 1, kjpindex
2016	!mask
2017	! vmr(ji,jv) = MAX ( EPSILON(un), MIN ( veget(ji,jv)-resdist(ji,jv) , MAX( EPSILON(un), veget(ji,jv)-resdist(ji,jv)) ) )
2018
2019	! vmr(ji,jv) = MAX ( EPSILON(un), MAX( EPSILON(un), veget(ji,jv)-resdist(ji,jv)) )
2020	! IF(ABS(veget(ji,jv)-resdist(ji,jv)).gt.epsilon(un)) then
2021	! WRITE(numout,*) '-----------------------------------------------'
2022	! WRITE(numout,*) 'vmr,epsilon(un),veget,resdist',vmr(ji,jv),epsilon(un)
2023	! WRITE(numout,*),veget(ji,jv),resdist(ji,jv)
2024	! WRITE(numout,*) 'ABS(veget -resdist',ABS(veget(ji,jv)-resdist(ji,jv))
2025	! endif
2026	IF ( ABS(veget(ji,jv)-resdist(ji,jv)) .GT. EPS1 ) THEN
2027	vmr(ji,jv) = veget(ji,jv)-resdist(ji,jv)
2028	ELSE
2029	vmr(ji,jv) = zero
2030	ENDIF
2031	!
2032	IF (resdist(ji,jv) .GT. min_sechiba) THEN
2033	qsintveg2(ji,jv) = qsintveg(ji,jv)/resdist(ji,jv)
2034	ELSE
2035	qsintveg2(ji,jv) = zero
2036	ENDIF
2037	ENDDO
2038	ENDDO
2039	!
2040	vegchtot(:) = zero
2041	DO jv = 1, nvm
2042	DO ji = 1, kjpindex
2043	vegchtot(ji) = vegchtot(ji) + ABS( vmr(ji,jv) )
2044	ENDDO
2045	ENDDO
2046	!
2047	DO jv = 1, nvm
2048	DO ji = 1, kjpindex
2049	IF ( vegchtot(ji) .GT. min_sechiba ) THEN
2050	qsdq(ji,jv) = ABS(vmr(ji,jv)) * qsintveg2(ji,jv)
2051	ENDIF
2052	ENDDO
2053	ENDDO
2054	!
2055	! calculate water mass that we have to redistribute
2056	!
2057	qstr(:) = zero
2058	vtr(:) = zero
2059	!
2060	!
2061	DO jv = 1, nvm
2062	DO ji = 1, kjpindex
2063	IF ( ( vegchtot(ji) .GT. min_sechiba ) .AND. ( vmr(ji,jv) .LT. -min_sechiba ) ) THEN
2064	qstr(ji) = qstr(ji) + qsdq(ji,jv)
2065	vtr(ji) = vtr(ji) - vmr(ji,jv)
2066	ENDIF
2067	ENDDO
2068	ENDDO
2069	!
2070	! put it into reservoir of plant whose surface area has grown
2071	DO jv = 1, nvm
2072	DO ji = 1, kjpindex
2073	IF ( vegchtot(ji) .GT. min_sechiba .AND. ABS(vtr(ji)) .GT. EPSILON(un)) THEN
2074	fra(ji) = vmr(ji,jv) / vtr(ji)
2075	IF ( vmr(ji,jv) .GT. min_sechiba) THEN
2076	qsintveg(ji,jv) = qsintveg(ji,jv) + fra(ji)* qstr(ji)
2077	ELSE
2078	qsintveg(ji,jv) = qsintveg(ji,jv) - qsdq(ji,jv)
2079	ENDIF
2080	ENDIF
2081	ENDDO
2082	ENDDO
2083	!MM underflow :
2084	DO jv = 1, nvm
2085	DO ji = 1, kjpindex
2086	IF ( ABS(qsintveg(ji,jv)) > zero .AND. ABS(qsintveg(ji,jv)) < EPS1 ) THEN
2087	qsintveg(ji,jv) = EPS1
2088	ENDIF
2089	ENDDO
2090	ENDDO
2091
2092	! Now that the work is done resdist needs an update !
2093	DO jv = 1, nvm
2094	DO ji = 1, kjpindex
2095	resdist(ji,jv) = veget(ji,jv)
2096	ENDDO
2097	ENDDO
2098
2099
2100	! Distribution of the vegetation depending on the soil type
2101
2102	! DO jst = 1, nstm
2103	!
2104	! DO ji = 1, kjpindex
2105	!
2106	!
2107	! soil_exist(ji,jst)=zero
2108	! IF (soiltype(ji,jst) .NE. zero) THEN
2109	! soil_exist(ji,jst)=un
2110	! soil_exist_max(ji,jst)=un
2111	! ENDIF
2112	!
2113	! ENDDO
2114	!
2115	! ENDDO
2116
2117	soil_exist(:,:) = mask_soiltype(:,:)
2118	soil_exist_max(:,:) = mask_soiltype(:,:)
2119	veget_exist(:,:) = zero
2120	veget_exist_max(:,:) = zero
2121
2122	DO jv = 1, nvm
2123
2124	DO ji = 1, kjpindex
2125	IF(vegtot(ji).GT.min_sechiba) THEN
2126	veget_exist(ji,jv)= veget(ji,jv)/vegtot(ji)
2127	veget_exist_max(ji,jv)= veget_max(ji,jv)/vegtot(ji)
2128	ENDIF
2129	ENDDO
2130	ENDDO
2131
2132	! Compute corr_veg_soil
2133
2134	corr_veg_soil(:,:,:) = zero
2135	corr_veg_soil_max(:,:,:) = zero
2136
2137	IF ( COUNT(pref_soil_veg .EQ. 0) > 0 ) THEN
2138
2139	DO jst = 1, nstm
2140
2141	DO jv = nvm, 1, -1
2142
2143	DO ji=1,kjpindex
2144
2145	IF(vegtot(ji).GT.min_sechiba.AND.soiltype(ji,jst).GT.min_sechiba) THEN
2146	corr_veg_soil(ji,jv,jst) = veget(ji,jv)/vegtot(ji)
2147	corr_veg_soil_max(ji,jv,jst) = veget_max(ji,jv)/vegtot(ji)
2148	END IF
2149
2150	END DO
2151	END DO
2152	END DO
2153
2154	ELSE
2155
2156
2157	DO jst = 1, nstm
2158
2159	DO jv = nvm, 1, -1
2160
2161	jst_pref = pref_soil_veg(jv,jst)
2162
2163	DO ji=1,kjpindex
2164	corr_veg_soil(ji,jv,jst_pref) = zero
2165	corr_veg_soil_max(ji,jv,jst_pref) = zero
2166	!for veget distribution used in sechiba via humrel
2167	IF (soil_exist(ji,jst_pref).GT.min_sechiba) THEN
2168	corr_veg_soil(ji,jv,jst_pref)=MIN(veget_exist(ji,jv)/soiltype(ji,jst_pref),soil_exist(ji,jst_pref))
2169	veget_exist(ji,jv)=MAX(veget_exist(ji,jv)-soil_exist(ji,jst_pref)*soiltype(ji,jst_pref),zero)
2170	soil_exist(ji,jst_pref)=MAX(soil_exist(ji,jst_pref)-corr_veg_soil(ji,jv,jst_pref),zero)
2171	ENDIF
2172	!same for max veget_max used in stomate via vegstress for slowproc
2173	IF (soil_exist_max(ji,jst_pref).GT.min_sechiba) THEN
2174	corr_veg_soil_max(ji,jv,jst_pref)= &
2175	& MIN(veget_exist_max(ji,jv)/soiltype(ji,jst_pref),soil_exist_max(ji,jst_pref))
2176	veget_exist_max(ji,jv)=MAX(veget_exist_max(ji,jv)-soil_exist_max(ji,jst_pref)*soiltype(ji,jst_pref),zero)
2177	soil_exist_max(ji,jst_pref)=MAX(soil_exist_max(ji,jst_pref)-corr_veg_soil_max(ji,jv,jst_pref),zero)
2178	ENDIF
2179	ENDDO
2180
2181	ENDDO
2182
2183	ENDDO
2184	ENDIF
2185	!
2186	! update the corresponding masks
2187	!
2188	! mask_veget(:,:) = MIN( un, MAX(zero,veget(:,:)))
2189	! mask_corr_veg_soil(:,:,:) = MIN( un, MAX(zero,corr_veg_soil(:,:,:)))
2190
2191	mask_veget(:,:) = 0
2192	mask_corr_veg_soil(:,:,:) = 0
2193
2194	DO ji = 1, kjpindex
2195
2196	DO jv = 1, nvm
2197	IF(veget(ji,jv) .GT. min_sechiba) THEN
2198	mask_veget(ji,jv) = 1
2199	ENDIF
2200
2201	DO jst = 1, nstm
2202	IF(corr_veg_soil(ji,jv,jst) .GT. min_sechiba) THEN
2203	mask_corr_veg_soil(ji,jv,jst) = 1
2204	ENDIF
2205	END DO
2206	END DO
2207
2208	! WRITE(numout,*) 'mask: soiltype,mask_soiltype',soiltype(ji,:),mask_soiltype(ji,:)
2209
2210	END DO
2211	!
2212
2213	RETURN
2214	!
2215	END SUBROUTINE hydrol_vegupd
2216	!!
2217	!! this routine computes soil processes with CWRR scheme
2218	!!
2219	SUBROUTINE hydrol_soil (kjpindex, dtradia, veget, veget_max, soiltype, transpir, vevapnu, evapot, &
2220	& evapot_penm, runoff, drainage, returnflow, irrigation, &
2221	& tot_melt, evap_bare_lim, shumdiag, litterhumdiag, humrel,vegstress, drysoil_frac)
2222	!
2223	! interface description
2224	! input scalar
2225	INTEGER(i_std), INTENT(in) :: kjpindex
2226	! input fields
2227	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
2228	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Map of vegetation types
2229	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Map of max vegetation types
2230	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
2231	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in) :: transpir !! transpiration
2232	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu !!
2233	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Water returning to the deep reservoir
2234	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: irrigation !! Irrigation
2235	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: evapot !!
2236	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: evapot_penm !!
2237	! modified fields
2238	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: runoff !! complete runoff
2239	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! complete drainage
2240	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt
2241	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !!
2242	REAL(r_std), DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! relative soil moisture
2243	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
2244	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (inout) :: humrel !! Relative humidity
2245	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
2246	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Function of the litter humidity, that will be used to compute albedo
2247	!
2248	! local declaration
2249	!
2250	INTEGER(i_std) :: ji, jv, jsl, jsl1, jst, ji_nsat !! indices
2251	INTEGER(i_std) :: m_sl0, m_sl1, isat !! mask values
2252	REAL(r_std) :: dztmp !! temporary depth
2253	REAL(r_std) :: temp !! temporary value for fluxes
2254	REAL(r_std) :: dpue !! temporary depth
2255	REAL(r_std), DIMENSION(kjpindex) :: tmcold, tmcint
2256	REAL(r_std), DIMENSION(kjpindex,nslm,nstm) :: moderwilt
2257	REAL(r_std), DIMENSION(kjpindex,nslm) :: mcint !! To save mc values for future use
2258	REAL(r_std), DIMENSION(kjpindex) :: correct_excess !! Corrects the flux at nslme layer in case of under residual moisture
2259	REAL(r_std), DIMENSION(kjpindex) :: mce !! Same use as mcint but jut for last efficient layer nslme
2260	REAL(r_std), DIMENSION(kjpindex) :: under_mcr !! Allows under residual soil moisture due to evap
2261	REAL(r_std), DIMENSION(kjpindex,nstm) :: v2
2262	REAL(r_std), DIMENSION(kjpindex,nstm) :: evap_bare_lim_ns !! limitation of bare soi evaporation on each soil column (used to deconvoluate vevapnu)
2263	REAL(r_std) :: deltahum,diff
2264	REAL(r_std), DIMENSION(kjpindex) :: tsink
2265	REAL(r_std), DIMENSION(kjpindex) :: returnflow_soil
2266	REAL(r_std), DIMENSION(kjpindex) :: irrigation_soil
2267	REAL(r_std), DIMENSION(kjpindex,nstm) :: runoff_excess !! Runoff generated after soil saturation
2268	REAL(r_std) :: excess
2269	LOGICAL :: propagate !! if we propagate an excess
2270	!
2271	!
2272	returnflow_soil(:) = zero
2273	irrigation_soil(:) = zero
2274	qflux(:,:,:) = zero
2275	mask_return(:) = 0
2276	index_nsat(:,:) = 0
2277	index_sat(:,:) = 0
2278	nslme(:,:) = nslm
2279	runoff_excess(:,:) = zero
2280	mce(:) = zero
2281	under_mcr(:) = zero
2282	correct_excess(:) = zero
2283	free_drain_coef(:,:) = zero
2284	n_sat(:) = 0
2285	n_nsat(:) = 1
2286	!
2287	! split 2d variables to 3d variables, per soil type
2288	!
2289	CALL hydrol_split_soil (kjpindex, veget, soiltype, vevapnu, transpir, humrel, evap_bare_lim)
2290	!
2291	! for each soil type
2292	!
2293	DO ji=1,kjpindex
2294	IF(vegtot(ji).GT.min_sechiba) THEN
2295	returnflow_soil(ji) = returnflow(ji)/vegtot(ji)
2296	irrigation_soil(ji) = irrigation(ji)/vegtot(ji)
2297	ENDIF
2298
2299	DO jst=1, nstm
2300	!- A priori on considere qu'on est non sature si soiltype > 0
2301	index_nsat(n_nsat(jst),jst)= ji*mask_soiltype(ji,jst)
2302	n_nsat(jst)=n_nsat(jst)+mask_soiltype(ji,jst)
2303	ENDDO
2304
2305	IF(returnflow_soil(ji).GT.min_sechiba) THEN
2306	mask_return(ji) = 1
2307	DO jst= 1,nstm
2308	isat=1
2309	nslme(ji,jst)=nslm-isat
2310	!
2311	DO jsl= nslm,3,-1
2312	IF(mcs(jst)-mc(ji,jsl,jst).LT.min_sechiba) THEN
2313	nslme(ji,jst) = nslme(ji,jst) - isat
2314	ELSE
2315	isat = 0
2316	ENDIF
2317	ENDDO
2318	! We compute the indeces of the non-saturated points
2319	IF (nslme(ji,jst).LT.2) THEN
2320	nslme(ji,jst) = 0
2321	! En fait on est sature!
2322	n_nsat(jst)=n_nsat(jst)-1
2323	n_sat(jst) = n_sat(jst)+1
2324	index_sat(n_sat(jst),jst)=ji
2325	ENDIF
2326	ENDDO
2327	ENDIF
2328	ENDDO
2329
2330	DO jst = 1,nstm
2331	n_nsat(jst) = n_nsat(jst)-1
2332	ENDDO
2333
2334	DO jst = 1,nstm
2335	!
2336	!- We compute the sum of the sinks for future check-up
2337	DO ji=1,kjpindex
2338	tsink(ji) = SUM(rootsink(ji,:,jst))+MAX(ae_ns(ji,jst),zero)+subsinksoil(ji)
2339	ENDDO
2340
2341	! We save the Total moisture content
2342	tmcold(:) = tmc(:,jst)
2343
2344	DO ji_nsat=1,n_nsat(jst)
2345	ji = index_nsat(ji_nsat,jst)
2346
2347	!- the bare soil evaporation is substracted to the soil moisture profile:
2348
2349	dpue = zz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst) / deux
2350	DO jsl = 1, nslme(ji,jst)
2351	mc(ji,jsl,jst) = mc(ji,jsl,jst) &
2352	& - (MAX(ae_ns(ji,jst),zero) + subsinksoil(ji)) / dpue
2353	ENDDO
2354
2355	!- we add the returnflow to the last efficient layer in the soil
2356	mc(ji,nslme(ji,jst),jst) = mc(ji,nslme(ji,jst),jst) + deux * returnflow_soil(ji) &
2357	& / (dz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst))
2358
2359	ENDDO
2360
2361	!!! SUBROUTINE hydrol_avoid_underres
2362	!-when mc(ji,1,jst)<mcr, we put it to mcr:
2363	! Smooth the profile to avoid negative values of ponctual soil moisture:
2364
2365	DO ji_nsat=1,n_nsat(jst)
2366	ji = index_nsat(ji_nsat,jst)
2367	!
2368	! Shifts water lack from top to bottom for under-residual moisture cases
2369	DO jsl = 1,nslm-1
2370	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2371	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2372	mc(ji,jsl+1,jst) = mc(ji,jsl+1,jst) - excess * &
2373	& (dz(jsl,jst)+dz(jsl+1,jst))/(dz(jsl+1,jst)+dz(jsl+2,jst))
2374	ENDDO
2375
2376	excess = MAX(mcr(jst)-mc(ji,nslm,jst),zero)
2377	mc(ji,nslm,jst) = mc(ji,nslm,jst) + excess
2378
2379	!- Then if the soil moisture at bottom is not sufficient, we try to refill the column from the top
2380	DO jsl = nslm-1,1,-1
2381	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess * &
2382	& (dz(jsl+1,jst)+dz(jsl+2,jst))/(dz(jsl+1,jst)+dz(jsl,jst))
2383	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2384	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2385	ENDDO
2386
2387	excess = excess * mask_soiltype(ji,jst)
2388	mc(ji,:,jst) = mc(ji,:,jst) * mask_soiltype(ji,jst)
2389
2390	! Keep the value in case excess is still positive (due to big change in evapot)
2391	under_mcr(ji) = excess * dz(2,jst)/2
2392	END DO
2393	!!! END SUBROUTINE hydrol_avoid_underres
2394
2395	!-we keep the value of mc in mcint:
2396
2397	DO jsl = 1, nslm
2398	DO ji = 1, kjpindex
2399	mcint(ji,jsl) = mc(ji,jsl,jst) - under_mcr(ji) / (dpu(jst) * mille)
2400	ENDDO
2401	ENDDO
2402
2403	DO ji = 1, kjpindex
2404	tmcint(ji) = dz(2,jst) * ( trois*mcint(ji,1) + mcint(ji,2) )/huit
2405	ENDDO
2406
2407	DO jsl = 2,nslm-1
2408	DO ji = 1, kjpindex
2409	tmcint(ji) = tmcint(ji) + dz(jsl,jst) &
2410	& * (trois*mcint(ji,jsl)+mcint(ji,jsl-1))/huit &
2411	& + dz(jsl+1,jst) * (trois*mcint(ji,jsl)+mcint(ji,jsl+1))/huit
2412	ENDDO
2413	END DO
2414	!
2415	DO ji = 1, kjpindex
2416	tmcint(ji) = tmcint(ji) + dz(nslm,jst) &
2417	& * (trois * mcint(ji,nslm) + mcint(ji,nslm-1))/huit
2418	ENDDO
2419
2420	!- On retire les termes puits de la transpiration des couches inactives a la premiere couche inactive.
2421	DO ji_nsat=1,n_nsat(jst)
2422	ji = index_nsat(ji_nsat,jst)
2423
2424	DO jsl = nslme(ji,jst)+1,nslm ! Allows to use mc(ji,nslme(ji,jst)+1,jst)
2425	mc(ji,nslme(ji,jst)+1,jst) = mc(ji,nslme(ji,jst)+1,jst) &
2426	& - deux * rootsink(ji,jsl,jst) / &
2427	& (dz(nslme(ji,jst)+1,jst) + dz(nslme(ji,jst)+2,jst))
2428	!- ATTENTION: n'est pas traite le cas ou la premiere couche inactive ne peut satisfaire la demande en eau des racines en dessous:
2429	!- Le cas ne devrait pas se poser a priori
2430
2431	ENDDO
2432	ENDDO
2433
2434	!- Some initialisation necessary for the diffusion scheme to work
2435	DO ji_nsat=1,n_nsat(jst)
2436	ji = index_nsat(ji_nsat,jst)
2437	!
2438	!
2439	!
2440	!- We correct rootsink for first two layers so that it is not too low in the first layer
2441	v1(ji,jst) = dz(2,jst)/huit * (trois * mc(ji,1,jst)+ mc(ji,2,jst))
2442	rootsink(ji,2,jst) = rootsink(ji,2,jst) + MAX(rootsink(ji,1,jst)-v1(ji,jst), zero)
2443	rootsink(ji,1,jst) = MIN(rootsink(ji,1,jst),v1(ji,jst))
2444	!- estimate maximum surface flux in mm/step, assuming
2445	!- all available water
2446	flux(ji) = zero
2447
2448	IF(vegtot(ji).GT.min_sechiba) THEN
2449	flux(ji) = precisol_ns(ji,jst) - evapot_penm(ji)*&
2450	& AINT(corr_veg_soil(ji,1,jst)+un-min_sechiba) &
2451	& + irrigation_soil(ji)
2452	ENDIF
2453	!- The incoming flux is first dedicated to fill the soil up to mcr (in case needed)
2454	temp = MAX(MIN(flux(ji),under_mcr(ji)), zero)
2455	flux(ji) = flux(ji) - temp
2456	under_mcr(ji) = under_mcr(ji) - temp
2457	END DO
2458
2459	!- module to implement dublin model for one time-step
2460	!- gravity drainage as lower boundary layer
2461	!- m.bruen, CWRR, ucd.
2462	!
2463	!-step3: matrix resolution
2464	!-calcul of the matrix coefficients
2465
2466	!- coefficient are computed depending on the profile mcint:
2467
2468	CALL hydrol_soil_setup(kjpindex,jst,dtradia)
2469
2470	!- Set the values for diffusion scheme
2471
2472
2473	DO ji_nsat=1,n_nsat(jst)
2474	ji = index_nsat(ji_nsat,jst)
2475
2476	! We only run the scheme in case we are not under mcr with the incoming flux
2477	IF (under_mcr(ji).LT.min_sechiba) THEN
2478	resolv(ji)=.TRUE.
2479	! In under residual case, we equally spread the transpiration over the layers
2480	ELSE
2481	under_mcr(ji) = under_mcr(ji) + SUM(rootsink(ji,:,jst))
2482	ENDIF
2483	!- First layer
2484
2485	tmat(ji,1,1) = zero
2486	tmat(ji,1,2) = f(ji,1)
2487	tmat(ji,1,3) = g1(ji,1)
2488	rhs(ji,1) = fp(ji,1) * mc(ji,1,jst) + gp(ji,1)*mc(ji,2,jst) &
2489	& + flux(ji) - (b(ji,1) + b(ji,2))*(dtradia/one_day)/deux - rootsink(ji,1,jst)
2490
2491	!- soil body
2492	DO jsl=2, nslme(ji,jst)-1
2493	tmat(ji,jsl,1) = e(ji,jsl)
2494	tmat(ji,jsl,2) = f(ji,jsl)
2495	tmat(ji,jsl,3) = g1(ji,jsl)
2496	rhs(ji,jsl) = ep(ji,jsl)mc(ji,jsl-1,jst) + fp(ji,jsl)mc(ji,jsl,jst) &
2497	& + gp(ji,jsl) * mc(ji,jsl+1,jst) &
2498	& + (b(ji,jsl-1) - b(ji,jsl+1)) * (dtradia/one_day) / deux &
2499	& - rootsink(ji,jsl,jst)
2500	ENDDO
2501
2502	!- Last layer
2503	jsl=nslme(ji,jst)
2504	tmat(ji,jsl,1) = e(ji,jsl)
2505	tmat(ji,jsl,2) = f(ji,jsl)
2506	tmat(ji,jsl,3) = zero
2507	rhs(ji,jsl) = ep(ji,jsl)mc(ji,jsl-1,jst) + fp(ji,jsl)mc(ji,jsl,jst) &
2508	& + (b(ji,jsl-1) - b(ji,jsl)) * (dtradia/one_day) / deux &
2509	& - rootsink(ji,jsl,jst)
2510
2511	!- store the equations in case needed again
2512	DO jsl=1,nslm
2513	srhs(ji,jsl) = rhs(ji,jsl)
2514	stmat(ji,jsl,1) = tmat(ji,jsl,1)
2515	stmat(ji,jsl,2) = tmat(ji,jsl,2)
2516	stmat(ji,jsl,3) = tmat(ji,jsl,3)
2517	ENDDO
2518	ENDDO
2519
2520	!
2521	!- step 4 : solve equations assuming atmosphere limiting
2522	!-
2523
2524	CALL hydrol_soil_tridiag(kjpindex,jst)
2525
2526
2527	!
2528	!- step 5 : check if really atmosphere limiting
2529	!-
2530	DO ji_nsat=1,n_nsat(jst)
2531	ji = index_nsat(ji_nsat,jst)
2532
2533	resolv(ji) = .FALSE.
2534	!
2535	!- Prepare to rerun in case of under residual with evaporation or over saturation
2536	!-
2537	IF(mc(ji,1,jst).LT.(mcr(jst)-min_sechiba).AND.evapot_penm(ji).GT.min_sechiba) THEN
2538
2539	!- upper layer dry
2540	!- We only rerun the scheme in case it is possible to reduce the evaporation
2541	resolv(ji) = .TRUE.
2542	DO jsl=1,nslm
2543	rhs(ji,jsl) = srhs(ji,jsl)
2544	tmat(ji,jsl,1) = stmat(ji,jsl,1)
2545	tmat(ji,jsl,2) = stmat(ji,jsl,2)
2546	tmat(ji,jsl,3) = stmat(ji,jsl,3)
2547	END DO
2548	tmat(ji,1,2) = un
2549	tmat(ji,1,3) = zero
2550	rhs(ji,1) = mcr(jst)-dmcr
2551
2552	ELSE
2553	IF(mc(ji,1,jst).GT.(mcs(jst)+0.02)) THEN
2554
2555	!- upper layer saturated
2556	resolv(ji) = .TRUE.
2557	DO jsl=1,nslm
2558	rhs(ji,jsl) = srhs(ji,jsl)
2559	tmat(ji,jsl,1) = stmat(ji,jsl,1)
2560	tmat(ji,jsl,2) = stmat(ji,jsl,2)
2561	tmat(ji,jsl,3) = stmat(ji,jsl,3)
2562	END DO
2563	tmat(ji,1,2) = un
2564	tmat(ji,1,3) = zero
2565	rhs(ji,1) = mcs(jst)+dmcs
2566
2567	ENDIF
2568	ENDIF
2569	ENDDO
2570
2571	!
2572	!- step 6 : resolve the equations with new boundary conditions if necessary
2573	!-
2574
2575	CALL hydrol_soil_tridiag(kjpindex,jst)
2576
2577	!-
2578	!- step 6.5 : initialize qflux at bottom of diffusion and avoid over saturated or under residual soil moisture
2579	!-
2580
2581	DO ji_nsat=1,n_nsat(jst)
2582	ji = index_nsat(ji_nsat,jst)
2583
2584	m_sl0 = mask_return(ji) ! the last efficient layer is above the last layer (nslm)
2585
2586	!- We add the flux from the last active layer to the first inactive one (not taken into account in the diffusion)
2587	!- At the same time, we initialize qflux(ji,jsl,jst) which is useful in case of no returnflow.
2588	!- When the first inactive point is to be saturated by the flux, the last efficient begins to fill up and the flux is changed
2589	jsl=nslme(ji,jst)
2590	qflux(ji,jsl,jst) = (a(ji,jsl)(w_timemc(ji,jsl,jst) &
2591	& + (un-w_time)mcint(ji,jsl)) + b(ji,jsl)) (dtradia/one_day)
2592
2593	mc(ji,jsl+m_sl0,jst) = mc(ji,jsl+m_sl0,jst) + &
2594	& m_sl0 * deux * qflux(ji,jsl,jst) / (dz(jsl+m_sl0,jst) + dz(jsl+m_sl0+1,jst))
2595
2596	excess = m_sl0 * MAX(mc(ji,jsl+m_sl0,jst)-mcs(jst),zero)
2597	mc(ji,jsl+m_sl0,jst) = mc(ji,jsl+m_sl0,jst) - excess
2598
2599	DO jsl1 = jsl*m_sl0,1,-1
2600	mc(ji,jsl1,jst) = mc(ji,jsl1,jst) + excess * &
2601	& (dz(jsl1+1,jst) + dz(jsl1+2,jst))/(dz(jsl1,jst) + dz(jsl1+1,jst))
2602	excess = MAX(mc(ji,jsl1,jst) - mcs(jst),zero)
2603	mc(ji,jsl1,jst) = mc(ji,jsl1,jst) - excess
2604	ENDDO
2605
2606	! Smooth the profile to avoid negative values of ponctual soil moisture
2607	DO jsl = 1,nslm-1
2608	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2609	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2610	mc(ji,jsl+1,jst) = mc(ji,jsl+1,jst) - excess * &
2611	& (dz(jsl,jst)+dz(jsl+1,jst))/(dz(jsl+1,jst)+dz(jsl+2,jst))
2612	ENDDO
2613
2614	excess = MAX(mcr(jst)-mc(ji,nslm,jst),zero)
2615	mc(ji,nslm,jst) = mc(ji,nslm,jst) + excess
2616
2617	!- Then if the soil moisture at bottom is not sufficient, we try to refill the column from the top
2618	DO jsl = nslm-1,1,-1
2619	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess * &
2620	& (dz(jsl+1,jst)+dz(jsl+2,jst))/(dz(jsl+1,jst)+dz(jsl,jst))
2621	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2622	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2623	ENDDO
2624
2625	excess = excess * mask_soiltype(ji,jst)
2626	mc(ji,:,jst) = mc(ji,:,jst) * mask_soiltype(ji,jst)
2627
2628	! Keep the value in case excess is still positive (due to big change in evapot)
2629	under_mcr(ji) = under_mcr(ji) + excess * dz(2,jst)/2
2630
2631	!- We do the opposite thing: in case of over-saturation we put the water where it is possible
2632	DO jsl = 1, nslm-1
2633	excess = MAX(mc(ji,jsl,jst)-mcs(jst),zero)
2634	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess
2635	mc(ji,jsl+1,jst) = mc(ji,jsl+1,jst) + excess * &
2636	& (dz(jsl,jst)+dz(jsl+1,jst))/(dz(jsl+1,jst)+dz(jsl+2,jst))
2637	ENDDO
2638
2639	DO jsl = nslm,2,-1
2640	excess = MAX(mc(ji,jsl,jst)-mcs(jst),zero)
2641	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess
2642	mc(ji,jsl-1,jst) = mc(ji,jsl-1,jst) + excess * &
2643	& (dz(jsl,jst)+dz(jsl+1,jst))/(dz(jsl-1,jst)+dz(jsl,jst))
2644	ENDDO
2645	excess = MAX(mc(ji,1,jst)-mcs(jst),zero)
2646	mc(ji,1,jst) = mc(ji,1,jst) - excess
2647
2648	ENDDO ! loop on grid
2649
2650	! Finaly, for soil that are under-residual, we just equally distribute the lack of water
2651	DO ji_nsat=1,n_nsat(jst)
2652	ji = index_nsat(ji_nsat,jst)
2653	DO jsl = 1, nslm
2654	mc(ji,jsl,jst) = mc(ji,jsl,jst) - under_mcr(ji) / (dpu(jst) * mille)
2655	ENDDO
2656	ENDDO
2657
2658	IF (check_cwrr) THEN
2659	DO ji_nsat=1,n_nsat(jst)
2660	ji = index_nsat(ji_nsat,jst)
2661	IF(qflux(ji,nslm,jst)+returnflow_soil(ji).LT.-min_sechiba.AND.soiltype(ji,jst).GT.min_sechiba) THEN
2662	WRITE(numout,*)'NEGATIVE FLUX AT LAST EFFICIENT LAYER IN SOIL'
2663	WRITE(numout,*)'mc[nlsm]_(t), mc[nslm]_(t-1),jst,soil,ji',&
2664	& mc(ji,nslm,jst),mcint(ji,nslm),jst,soiltype(ji,jst),ji
2665	WRITE(numout,*)'irrigation,returnflow,fdc',irrigation_soil(ji),returnflow_soil(ji)
2666	ENDIF
2667	END DO
2668	ENDIF
2669	!
2670	!- step 7 : close the water balance
2671	!
2672	!- drainage through the lower boundary
2673	!- and fluxes for each soil layer
2674	!- with mass balance computed from the bottom to the top
2675	!- of the soil column
2676
2677	!- Compute the flux at every level from bottom to top (using mc and sink values)
2678
2679	DO ji_nsat=1,n_nsat(jst)
2680	ji = index_nsat(ji_nsat,jst)
2681	m_sl0 = mask_return(ji) ! the last efficient layer is above the last layer
2682
2683	DO jsl=nslm-1,nslme(ji,jst),-1
2684	qflux(ji,jsl,jst) = qflux(ji,jsl+1,jst) &
2685	& + (mc(ji,jsl+1,jst) - mcint(ji,jsl+1)) &
2686	& * (dz(jsl+1,jst)+dz(jsl+2,jst))/deux &
2687	& + rootsink(ji,jsl+1,jst)
2688	ENDDO
2689
2690	jsl = nslme(ji,jst)-1
2691	qflux(ji,jsl,jst) = qflux(ji,jsl+1,jst) &
2692	& + (mc(ji,jsl,jst)-mcint(ji,jsl) &
2693	& + troismc(ji,jsl+1,jst) - troismcint(ji,jsl+1)) &
2694	& * (dz(jsl+1,jst)/huit) &
2695	& + rootsink(ji,jsl+1,jst) &
2696	& + (dz(jsl+2,jst)/deux) & ! zero if nslme=nslm
2697	& * (mc(ji,jsl+1,jst) - mcint(ji,jsl+1))
2698
2699	DO jsl = nslme(ji,jst)-2,1,-1
2700	qflux(ji,jsl,jst) = qflux(ji,jsl+1,jst) &
2701	& + (mc(ji,jsl,jst)-mcint(ji,jsl) &
2702	& + troismc(ji,jsl+1,jst) - troismcint(ji,jsl+1)) &
2703	& * (dz(jsl+1,jst)/huit) &
2704	& + rootsink(ji,jsl+1,jst) &
2705	& + (dz(jsl+2,jst)/huit) &
2706	& * (troismc(ji,jsl+1,jst) - troismcint(ji,jsl+1) &
2707	& + mc(ji,jsl+2,jst)-mcint(ji,jsl+2))
2708	END DO
2709
2710	qflux00(ji,jst) = qflux(ji,1,jst) + (dz(2,jst)/huit) &
2711	& * (trois* (mc(ji,1,jst)-mcint(ji,1)) + (mc(ji,2,jst)-mcint(ji,2))) &
2712	& + rootsink(ji,1,jst)
2713	ENDDO
2714
2715	!- Then computes the water balance (evap-runoff-drainage)
2716
2717	DO ji_nsat=1,n_nsat(jst)
2718	ji = index_nsat(ji_nsat,jst)
2719	!
2720	!
2721	! deduction of ae_ns and ru_ns:
2722	! ae_ns+ru_ns=precisol_ns+irrigation-q0
2723	!
2724	ae_ns(ji,jst) = MAX(MIN((precisol_ns(ji,jst)+irrigation_soil(ji)-qflux00(ji,jst)),evapot_penm(ji)),zero)
2725	ru_ns(ji,jst) = precisol_ns(ji,jst)+irrigation_soil(ji)-qflux00(ji,jst)-ae_ns(ji,jst) !+runoff_excess(ji,jst)
2726	!
2727	! In case of negative runoff, we correct it by taking water from the soil
2728	! Le probleme est que desormais, les qflux ne sont plus justes...
2729	! A corriger plus tard...
2730	IF (ru_ns(ji,jst).LT.-min_sechiba) THEN
2731	WRITE(numout,*) 'Negative runoff corrected', ru_ns(ji,jst), mc(ji,1,jst), SUM(rootsink(ji,:,jst))
2732	WRITE(numout,*) 'At...', ji, jst, mask_soiltype(ji,jst), nslme(ji,jst)
2733	excess = -ru_ns(ji,jst)
2734	ru_ns(ji,jst) = zero
2735	! We correct this by taking water from the whole soil
2736	qflux00(ji,jst) = qflux00(ji,jst) - excess
2737	dpue = zz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst) / deux
2738	DO jsl = 1, nslme(ji,jst)
2739	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess / dpue
2740	ENDDO
2741	! Then we have to check if there is no negative value
2742
2743	DO jsl = 1,nslm-1
2744	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2745	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2746	mc(ji,jsl+1,jst) = mc(ji,jsl+1,jst) - excess * &
2747	& (dz(jsl,jst)+dz(jsl+1,jst))/(dz(jsl+1,jst)+dz(jsl+2,jst))
2748	ENDDO
2749
2750	excess = MAX(mcr(jst)-mc(ji,nslm,jst),zero)
2751	mc(ji,nslm,jst) = mc(ji,nslm,jst) + excess
2752
2753	!- Then if the soil moisture at bottom is not sufficient, we try to refill the column from the top
2754	DO jsl = nslm-1,1,-1
2755	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess * &
2756	& (dz(jsl+1,jst)+dz(jsl+2,jst))/(dz(jsl+1,jst)+dz(jsl,jst))
2757	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2758	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2759	ENDDO
2760
2761	excess = excess * mask_soiltype(ji,jst)
2762	mc(ji,:,jst) = mc(ji,:,jst) * mask_soiltype(ji,jst)
2763
2764	! And if excess is still positive, we put the soil under the residual value:
2765	DO jsl = 1, nslm
2766	mc(ji,jsl,jst) = mc(ji,jsl,jst) - excess / (dpu(jst) * mille)
2767	ENDDO
2768	ENDIF
2769
2770	dr_ns(ji,jst) = qflux(ji,nslm,jst)
2771
2772	IF (ABS(ae_ns(ji,jst)).LT.min_sechiba) THEN
2773	ae_ns(ji,jst) = zero
2774	ENDIF
2775
2776	IF(ABS(ru_ns(ji,jst)).LT.min_sechiba) THEN
2777	ru_ns(ji,jst) = zero
2778	ENDIF
2779
2780	IF(ABS(dr_ns(ji,jst)).LT.min_sechiba) THEN
2781	dr_ns(ji,jst) = zero
2782	ENDIF
2783
2784	! We add the evaporation to the soil profile
2785
2786	dpue = zz(nslme(ji,jst),jst) + dz(nslme(ji,jst)+1,jst) / deux
2787	DO jsl = 1, nslme(ji,jst)
2788	mc(ji,jsl,jst) = mc(ji,jsl,jst) + ae_ns(ji,jst) / dpue
2789	ENDDO
2790	END DO
2791
2792	!
2793	!- Special case for saturated soil
2794	!- We did not use the diffusion and just use a sort of bucket system
2795
2796	DO ji_nsat=1,n_sat(jst)
2797	ji = index_sat(ji_nsat,jst)
2798	dr_ns(ji,jst) = zero ! -returnflow_soil(ji)
2799	m_sl0 = mask_soiltype(ji,jst)
2800	!
2801	!- mc1 and mc2 calculation
2802	! Calculation of mc1 after last timestep evap and after precip and irrigation
2803	mc(ji,1,jst) = mc(ji,1,jst) + (precisol_ns(ji,jst) + irrigation_soil(ji) - ae_ns(ji,jst)) &
2804	& * 2 / dz(2,jst)
2805	! Preparing to take water from bottom in case of under-residual mc1
2806	excess = MAX(mcr(jst)-mc(ji,1,jst),zero)
2807	mc(ji,1,jst) = mc(ji,1,jst) + excess
2808	! Calculation of mc2 with returnflow, transpiration and probable lack of water in first layer
2809	mc(ji,2,jst) = mc(ji,2,jst) + (returnflow_soil(ji) - tsink(ji) + ae_ns(ji,jst) - &
2810	& excess * dz(2,jst) / 2) * 2/(dz(2,jst) + dz(3,jst))
2811	! Shifting water to top in case of saturation (returnflow very strong)
2812	excess = MAX(mc(ji,2,jst)-mcs(jst),zero)
2813	mc(ji,2,jst) = mc(ji,2,jst) - excess
2814	mc(ji,1,jst) = mc(ji,1,jst) + excess * (dz(2,jst) + dz(3,jst)) / dz(2,jst)
2815	! Avoiding under residual soil moisture for mc2 if transpiration was very high
2816	DO jsl = 2, nslm-1
2817	excess = MAX(mcr(jst)-mc(ji,jsl,jst),zero)
2818	mc(ji,jsl,jst) = mc(ji,jsl,jst) + excess
2819	mc(ji,jsl+1,jst) = mc(ji,jsl+1,jst) - excess * &
2820	& (dz(jsl,jst) + dz(jsl+1,jst))/(dz(jsl+1,jst) + dz(jsl+2,jst))
2821	ENDDO
2822	excess = m_sl0 * excess
2823	mc(ji,:,jst) = m_sl0 * mc(ji,:,jst)
2824	IF (excess .GT. min_sechiba) THEN
2825	STOP 'Saturated soil evaporating everything... Oups...'
2826	ENDIF
2827	!
2828	!- Deduction of ae_ns (used for next step) allowing first two layers drying out
2829	! We first calculate the water that can be evaporated from the first layer and deduce the runoff
2830	ae_ns(ji,jst) = m_sl0 * MIN((mc(ji,1,jst)-mcr(jst))*dz(2,jst)/2,evapot_penm(ji))
2831	! Generating runoff in case of remaining over saturation in the first layer
2832	excess = m_sl0 * MIN(MAX(mc(ji,1,jst)-mcs(jst),zero),mc(ji,1,jst)-mcr(jst)-ae_ns(ji,jst)*2/dz(2,jst))
2833	mc(ji,1,jst) = mc(ji,1,jst) - excess
2834	ru_ns(ji,jst) = m_sl0 * excess * dz(2,jst)/2
2835	! If it was not sufficient for ae_ns to reach evapot, we evaporate water from the second layer.
2836	ae_ns(ji,jst) = ae_ns(ji,jst) + m_sl0 * &
2837	& MIN((mc(ji,2,jst)-mcr(jst))*(dz(2,jst)+dz(3,jst))/2, evapot_penm(ji) - ae_ns(ji,jst))
2838
2839	! calculation of qflux from what precedes
2840	qflux00(ji,jst) = m_sl0 * (precisol_ns(ji,jst) + irrigation_soil(ji) &
2841	& - ae_ns(ji,jst) - ru_ns(ji,jst))
2842
2843	IF (ABS(ae_ns(ji,jst)).LT.min_sechiba) THEN
2844	ae_ns(ji,jst) = zero
2845	ENDIF
2846
2847	IF(ABS(ru_ns(ji,jst)).LT.min_sechiba) THEN
2848	ru_ns(ji,jst) = zero
2849	ENDIF
2850	ENDDO
2851
2852
2853	!
2854	!- step8: we make some useful output
2855	!- Total soil moisture, soil moisture at litter levels, soil wetness...
2856	!
2857
2858	!-total soil moisture:
2859
2860	DO ji=1,kjpindex
2861	tmc(ji,jst)= dz(2,jst) * (trois*mc(ji,1,jst) + mc(ji,2,jst))/huit
2862	END DO
2863
2864	DO jsl=2,nslm-1
2865	DO ji=1,kjpindex
2866	tmc(ji,jst) = tmc(ji,jst) + dz(jsl,jst) * ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
2867	& + dz(jsl+1,jst)(troismc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
2868	END DO
2869	END DO
2870
2871	DO ji=1,kjpindex
2872	tmc(ji,jst) = tmc(ji,jst) + dz(nslm,jst) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
2873	END DO
2874
2875	! the litter is the 4 top levels of the soil
2876	! we compute various field of soil moisture for the litter (used for stomate and for albedo)
2877
2878	DO ji=1,kjpindex
2879	tmc_litter(ji,jst) = dz(2,jst) * ( trois*mc(ji,1,jst)+ mc(ji,2,jst))/huit
2880	END DO
2881
2882
2883	! sum from level 1 to 4
2884
2885	DO jsl=2,4
2886
2887	DO ji=1,kjpindex
2888	tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl,jst) * &
2889	& ( trois*mc(ji,jsl,jst) + mc(ji,jsl-1,jst))/huit &
2890	& + dz(jsl+1,jst)(troismc(ji,jsl,jst) + mc(ji,jsl+1,jst))/huit
2891	END DO
2892
2893	END DO
2894
2895
2896	! subsequent calcul of soil_wet_litter (tmc-tmcw)/(tmcf-tmcw)
2897
2898	DO ji=1,kjpindex
2899	soil_wet_litter(ji,jst) = MIN(un, MAX(zero,&
2900	& (tmc_litter(ji,jst)-tmc_litter_wilt(ji,jst)) / &
2901	& (tmc_litter_field(ji,jst)-tmc_litter_wilt(ji,jst)) ))
2902	END DO
2903
2904	! Soil wetness profiles (mc-mcw)/(mcs-mcw)
2905	! soil_wet is the ratio of soil moisture to available soil moisture for plant
2906	! (ie soil moisture at saturation minus soil moisture at wilting point).
2907
2908	DO ji=1,kjpindex
2909	soil_wet(ji,1,jst) = MIN(un, MAX(zero,&
2910	& (troismc(ji,1,jst) + mc(ji,2,jst) - quatremcw(jst))&
2911	& /(quatre*(mcs(jst)-mcw(jst))) ))
2912	humrelv(ji,1,jst) = zero
2913	END DO
2914
2915	DO jsl=2,nslm-1
2916	DO ji=1,kjpindex
2917	soil_wet(ji,jsl,jst) = MIN(un, MAX(zero,&
2918	& (trois*mc(ji,jsl,jst) + &
2919	& mc(ji,jsl-1,jst) *(dz(jsl,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) &
2920	& + mc(ji,jsl+1,jst)*(dz(jsl+1,jst)/(dz(jsl,jst)+dz(jsl+1,jst))) &
2921	& - quatremcw(jst)) / (quatre(mcs(jst)-mcw(jst))) ))
2922	END DO
2923	END DO
2924
2925	DO ji=1,kjpindex
2926	soil_wet(ji,nslm,jst) = MIN(un, MAX(zero,&
2927	& (trois*mc(ji,nslm,jst) &
2928	& + mc(ji,nslm-1,jst)-quatremcw(jst))/(quatre(mcs(jst)-mcw(jst))) ))
2929	END DO
2930
2931	!
2932	!- step8: we make the outputs for sechiba:
2933	!-we compute the moderation of transpiration due to wilting point:
2934	! moderwilt is a factor which is zero if soil moisture is below the wilting point
2935	! and is un if soil moisture is above the wilting point.
2936
2937
2938	DO jsl=1,nslm
2939	DO ji=1,kjpindex
2940	moderwilt(ji,jsl,jst) = INT( MAX(soil_wet(ji,jsl,jst), zero) + un - min_sechiba )
2941	END DO
2942	END DO
2943
2944	!- we compute the new humrelv to use in sechiba:
2945	!- loop on each vegetation type
2946
2947	humrelv(:,1,jst) = zero
2948
2949	DO jv = 2,nvm
2950
2951	!- calcul of us for each layer and vegetation type.
2952
2953	DO ji=1,kjpindex
2954	us(ji,jv,jst,1) = moderwilt(ji,1,jst)MIN(un,((troismc(ji,1,jst) + mc(ji,2,jst)) &
2955	& /(quatremcs(jst)pcent(jst))) )* (un-EXP(-humcste(jv)*dz(2,jst)/mille/deux)) &
2956	& /(un-EXP(-humcste(jv)*zz(nslm,jst)/mille))
2957	us(ji,jv,jst,1) = zero
2958	humrelv(ji,jv,jst) = MAX(us(ji,jv,jst,1),zero)
2959	END DO
2960
2961	DO jsl = 2,nslm-1
2962	DO ji=1,kjpindex
2963	us(ji,jv,jst,jsl) =moderwilt(ji,jsl,jst)* &
2964	& MIN( un, &
2965	& ((trois*mc(ji,jsl,jst)+ &
2966	& mc(ji,jsl-1,jst)*(dz(jsl,jst)/(dz(jsl,jst)+dz(jsl+1,jst)))+ &
2967	& mc(ji,jsl+1,jst)*(dz(jsl+1,jst)/(dz(jsl,jst)+dz(jsl+1,jst)))) &
2968	& /(quatremcs(jst)pcent(jst))) )* &
2969	& (EXP(-humcste(jv)zz(jsl,jst)/mille)) &
2970	& (EXP(humcste(jv)*dz(jsl,jst)/mille/deux) - &
2971	& EXP(-humcste(jv)*dz(jsl+1,jst)/mille/deux))/ &
2972	& (EXP(-humcste(jv)*dz(2,jst)/mille/deux) &
2973	& -EXP(-humcste(jv)*zz(nslm,jst)/mille))
2974
2975	us(ji,jv,jst,jsl) = MAX(us (ji,jv,jst,jsl), zero)
2976	humrelv(ji,jv,jst) = MAX((humrelv(ji,jv,jst) + us(ji,jv,jst,jsl)),zero)
2977	END DO
2978	END DO
2979
2980	DO ji=1,kjpindex
2981	us(ji,jv,jst,nslm) =moderwilt(ji,nslm,jst)* &
2982	& MIN(un, &
2983	& ((trois*mc(ji,nslm,jst) + mc(ji,nslm-1,jst)) &
2984	& / (quatremcs(jst)pcent(jst))) ) * &
2985	& (EXP(humcste(jv)dz(nslm,jst)/mille/deux) -un) &
2986	& EXP(-humcste(jv)*zz(nslm,jst)/mille) / &
2987	& (EXP(-humcste(jv)*dz(2,jst)/mille/deux) &
2988	& -EXP(-humcste(jv)*zz(nslm,jst)/mille))
2989	us(ji,jv,jst,nslm) = MAX(us(ji,jv,jst,nslm), zero)
2990	humrelv(ji,jv,jst) = MAX(zero,MIN(un, humrelv(ji,jv,jst) + us(ji,jv,jst,nslm)))
2991	vegstressv(ji,jv,jst) = humrelv(ji,jv,jst)
2992	humrelv(ji,jv,jst) = humrelv(ji,jv,jst) * mask_corr_veg_soil(ji,jv,jst)
2993	! IF(corr_veg_soil(ji,jv,jst).EQ.zero) THEN
2994	! humrelv(ji,jv,jst) = zero
2995	! ENDIF
2996	END DO
2997	END DO
2998
2999	!before closing the soil water, we check the water balance of soil
3000
3001	IF(check_cwrr) THEN
3002	DO ji = 1,kjpindex
3003
3004	deltahum = (tmc(ji,jst) - tmcold(ji))
3005	diff = precisol_ns(ji,jst)-ru_ns(ji,jst)-dr_ns(ji,jst)-tsink(ji) + irrigation_soil(ji) + returnflow_soil(ji)
3006
3007	IF(abs(deltahum-diff)mask_soiltype(ji,jst).gt.deuxallowed_err) THEN
3008
3009	WRITE(numout,*) 'CWRR pat: bilan non nul',ji,jst,deltahum-diff
3010	WRITE(numout,*) 'tmc,tmcold,diff',tmc(ji,jst),tmcold(ji),deltahum
3011	WRITE(numout,*) 'evapot,evapot_penm,ae_ns',evapot(ji),evapot_penm(ji),ae_ns(ji,jst)
3012	WRITE(numout,*) 'flux,ru_ns,qdrain,tsink,q0,precisol,excess',flux(ji),ru_ns(ji,jst), &
3013	& dr_ns(ji,jst),tsink(ji),qflux00(ji,jst),precisol_ns(ji,jst),runoff_excess(ji,jst)
3014	WRITE(numout,*) 'soiltype',soiltype(ji,jst)
3015	WRITE(numout,*) 'irrigation,returnflow',irrigation_soil(ji),returnflow_soil(ji)
3016	WRITE(numout,*) 'mc',mc(ji,:,jst)
3017	WRITE(numout,*) 'nslme',nslme(ji,jst)
3018	WRITE(numout,*) 'qflux',qflux(ji,:,jst)
3019	WRITE(numout,*) 'correct,mce',correct_excess(ji),mce(ji)
3020	STOP 'in hydrol_soil CWRR water balance check'
3021
3022	ENDIF
3023
3024	IF(MINVAL(mc(ji,:,jst)).LT.-min_sechiba) THEN
3025	WRITE(numout,*) 'CWRR MC NEGATIVE', &
3026	& ji,lalo(ji,:),MINLOC(mc(ji,:,jst)),jst,mc(ji,:,jst)
3027	WRITE(numout,*) 'evapot,ae_ns',evapot(ji),ae_ns(ji,jst)
3028	WRITE(numout,*) 'returnflow,irrigation,nslme',returnflow_soil(ji),&
3029	& irrigation_soil(ji),nslme(ji,jst)
3030	WRITE(numout,*) 'flux,ru_ns,qdrain,tsink,q0',flux(ji),ru_ns(ji,jst), &
3031	& dr_ns(ji,jst),tsink(ji),qflux00(ji,jst)
3032	WRITE(numout,*) 'soiltype',soiltype(ji,jst)
3033	STOP 'in hydrol_soil CWRR MC NEGATIVE'
3034	ENDIF
3035	END DO
3036
3037	DO ji_nsat=1,n_nsat(jst)
3038	ji = index_nsat(ji_nsat,jst)
3039	IF (ru_ns(ji,jst).LT.-min_sechiba) THEN
3040	WRITE(numout,*) 'Negative runoff in non-saturated case', ji,jst, mask_soiltype(ji,jst)
3041	WRITE(numout,*) 'mc1, mc2, nslme', mc(ji,1,jst), mc(ji,2,jst), nslme(ji,jst)
3042	WRITE(numout,*) 'mcint1, mcint2, mce', mcint(ji,1), mcint(ji,2), mce(ji)
3043	WRITE(numout,*) 'qflux1, correct, flux', qflux(ji,nslm,jst), correct_excess(ji), flux(ji)
3044	WRITE(numout,*) 'under_mcr, test', under_mcr(ji), tmc(ji,jst)-tmcint(ji)+qflux(ji,nslm,jst)+SUM(rootsink(ji,:,jst))
3045	WRITE(numout,*) 'mc', mc(ji,:,jst)
3046	WRITE(numout,*) 'mcint', mcint(ji,:)
3047	WRITE(numout,*) 'qflux', qflux(ji,:,jst)
3048	WRITE(numout,*) 'rootsink1,evapot_penm,vegtot', rootsink(ji,1,jst), evapot_penm(ji), vegtot(ji)
3049	WRITE(numout,*) 'ae_ns, tsink, returnflow, precisol_ns, irrigation, qflux0, ru_ns', &
3050	& ae_ns(ji,jst), tsink(ji), returnflow_soil(ji), &
3051	& precisol_ns(ji,jst), irrigation_soil(ji), qflux00(ji,jst), ru_ns(ji,jst)
3052	STOP 'STOP in hydrol_soil: Negative runoff, non-saturated soil'
3053	ENDIF
3054	ENDDO
3055
3056	DO ji_nsat=1,n_sat(jst)
3057	ji = index_sat(ji_nsat,jst)
3058	m_sl0 = mask_soiltype(ji,jst)
3059	IF (ru_ns(ji,jst).LT.-min_sechiba) THEN
3060	WRITE(numout,*) 'Negative runoff in saturated case', ji,jst, mask_soiltype(ji,jst), &
3061	& mc(ji,1,jst), mc(ji,2,jst), ae_ns(ji,jst), tsink(ji), returnflow_soil(ji), &
3062	& precisol_ns(ji,jst), irrigation_soil(ji), qflux00(ji,jst), ru_ns(ji,jst)
3063	STOP 'STOP in hydrol_soil: Negative runoff, saturated soil'
3064	ENDIF
3065	ENDDO
3066	ENDIF
3067
3068	END DO ! end of loop on soiltype
3069
3070	!
3071	! sum 3d variables into 2d variables
3072	!
3073	CALL hydrol_diag_soil (kjpindex, veget, veget_max, soiltype, runoff, drainage, &
3074	& evap_bare_lim, evapot, vevapnu, returnflow, irrigation, &
3075	& shumdiag, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt)
3076	RETURN
3077
3078	END SUBROUTINE hydrol_soil
3079
3080	SUBROUTINE hydrol_soil_tridiag(kjpindex,ins)
3081
3082	!- solves a set of linear equations which has a tridiagonal
3083	!- coefficient matrix.
3084
3085	!- arguments
3086	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
3087	INTEGER(i_std), INTENT(in) :: ins !! number of soil type
3088
3089	! -- variables locales
3090
3091	INTEGER(i_std) :: ji,jt,jsl,ji_nsat
3092	REAL(r_std), DIMENSION(kjpindex) :: bet
3093
3094	DO ji_nsat = 1,n_nsat(ins)
3095	ji = index_nsat(ji_nsat,ins)
3096
3097	IF (resolv(ji)) THEN
3098	bet(ji) = tmat(ji,1,2)
3099	mc(ji,1,ins) = rhs(ji,1)/bet(ji)
3100
3101	DO jsl = 2,nslme(ji,ins)
3102	gam(ji,jsl) = tmat(ji,jsl-1,3)/bet(ji)
3103	bet(ji) = tmat(ji,jsl,2) - tmat(ji,jsl,1)*gam(ji,jsl)
3104	mc(ji,jsl,ins) = (rhs(ji,jsl)-tmat(ji,jsl,1)*mc(ji,jsl-1,ins))/bet(ji)
3105	ENDDO
3106
3107	DO jsl = nslme(ji,ins)-1,1,-1
3108	mc(ji,jsl,ins) = mc(ji,jsl,ins) - gam(ji,jsl+1)*mc(ji,jsl+1,ins)
3109	ENDDO
3110	ENDIF
3111	ENDDO
3112	RETURN
3113	END SUBROUTINE hydrol_soil_tridiag
3114
3115
3116	SUBROUTINE hydrol_soil_setup(kjpindex,ins,dtradia)
3117
3118	!
3119	!**** hydrol_soil_setup -
3120	!**** *routine that computes the matrix coef for dublin model.
3121	!**** *uses the linearised hydraulic conductivity k_lin=a_lin mc_lin+b_lin
3122	!**** *and the linearised diffusivity d_lin
3123	!
3124	IMPLICIT NONE
3125	!
3126	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
3127	! parameters
3128	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
3129	INTEGER(i_std), INTENT(in) :: ins ! index of soil type
3130	! local
3131	INTEGER(i_std) :: jsl,ji,i,m_sl0,m_sl1, ji_nsat
3132	REAL(r_std), DIMENSION (nslm) :: temp0, temp5, temp6, temp7
3133	REAL(r_std) :: temp3, temp4
3134
3135	!-first, we identify the interval i in which the current value of mc is located
3136	!-then, we give the values of the linearized parameters to compute
3137	! conductivity and diffusivity as K=a*mc+b and d
3138
3139	DO jsl=1,nslm
3140	DO ji=1,kjpindex
3141	i= MIN(INT((imax-imin)*(MAX(mc(ji,jsl,ins),mcr(ins))-mcr(ins))&
3142	& / (mcs(ins)-mcr(ins)))+imin , imax-1)
3143	a(ji,jsl) = a_lin(i,ins)
3144	b(ji,jsl) = b_lin(i,ins)
3145	d(ji,jsl) = d_lin(i,ins)
3146	ENDDO ! loop on grid
3147	ENDDO
3148
3149
3150	!-second, we compute tridiag matrix coefficients (LEFT and RIGHT)
3151	! of the system to solve [LEFT]mc_{t+1}=[RIGHT]mc{t}+[add terms]:
3152	! e(nslm),f(nslm),g1(nslm) for the [left] vector
3153	! and ep(nslm),fp(nslm),gp(nslm) for the [right] vector
3154
3155	temp3 = w_time*(dtradia/one_day)/deux
3156	temp4 = (un-w_time)*(dtradia/one_day)/deux
3157
3158	DO ji_nsat=1,n_nsat(ins)
3159	ji = index_nsat(ji_nsat,ins)
3160
3161	!- First layer temporary calc
3162	!- Be careful! The order (first layer before last) is very important in case nslme(ji,jst)=1
3163	temp0(1) = trois * dz(2,ins)/huit
3164	temp5(1) = zero
3165	temp6(1) = (d(ji,1)+d(ji,2))/(dz(2,ins)) + a(ji,1)
3166	temp7(1) = (d(ji,1)+d(ji,2))/(dz(2,ins)) - a(ji,2)
3167
3168	!- Main body
3169	DO jsl = 2, nslme(ji,ins)-1
3170	temp0(jsl) = trois * (dz(jsl,ins) + dz(jsl+1,ins))/huit
3171	temp5(jsl) =(d(ji,jsl)+d(ji,jsl-1))/(dz(jsl,ins))+a(ji,jsl-1)
3172	temp6(jsl) = (d(ji,jsl)+d(ji,jsl-1))/(dz(jsl,ins)) + &
3173	& (d(ji,jsl)+d(ji,jsl+1))/(dz(jsl+1,ins))
3174	temp7(jsl) = (d(ji,jsl)+d(ji,jsl+1))/(dz(jsl+1,ins)) &
3175	& - a(ji,jsl+1)
3176	ENDDO
3177
3178	!- Last layer
3179	jsl = nslme(ji,ins)
3180	temp0(jsl) = trois * (dz(jsl,ins) + dz(jsl+1,ins))/huit &
3181	& + dz(jsl+1,ins)/huit
3182	temp5(jsl) = (d(ji,jsl)+d(ji,jsl-1)) / (dz(jsl,ins)) &
3183	& + a(ji,jsl-1)
3184	temp6(jsl) = (d(ji,jsl)+d(ji,jsl-1))/dz(jsl,ins) &
3185	& + a(ji,jsl)
3186	temp7(jsl) = zero
3187
3188	!- coefficient for every layer
3189	DO jsl = 1, nslme(ji,ins)
3190	e(ji,jsl) = dz(jsl,ins)/(huit) - temp3*temp5(jsl)
3191	f(ji,jsl) = temp0(jsl) + temp3*temp6(jsl)
3192	g1(ji,jsl) = dz(jsl+1,ins)/(huit) - temp3*temp7(jsl)
3193	ep(ji,jsl) = dz(jsl,ins)/(huit) + temp4*temp5(jsl)
3194	fp(ji,jsl) = temp0(jsl) - temp4*temp6(jsl)
3195	gp(ji,jsl) = dz(jsl+1,ins)/(huit) + temp4*temp7(jsl)
3196	ENDDO
3197	ENDDO
3198
3199	RETURN
3200	END SUBROUTINE hydrol_soil_setup
3201
3202	!!! fait la connexion entre l'hydrologie et sechiba :
3203	!!! cherche les variables sechiba pour l'hydrologie
3204	!!! "transforme" ces variables
3205	SUBROUTINE hydrol_split_soil (kjpindex, veget, soiltype, vevapnu, transpir, humrel,evap_bare_lim)
3206	!
3207	! interface description
3208	! input scalar
3209	INTEGER(i_std), INTENT(in) :: kjpindex
3210	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(in) :: veget !! Vegetation map
3211	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
3212	REAL(r_std), DIMENSION (kjpindex), INTENT (in) :: vevapnu !! Bare soil evaporation
3213	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
3214	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: humrel !! Relative humidity
3215	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: evap_bare_lim !!
3216	!
3217	! local declaration
3218	!
3219	INTEGER(i_std) :: ji, jv, jsl, jst
3220	REAL(r_std), Dimension (kjpindex) :: vevapnu_old
3221	REAL(r_std), Dimension (kjpindex) :: tmp_check1
3222	REAL(r_std), Dimension (kjpindex) :: tmp_check2
3223	REAL(r_std), DIMENSION (kjpindex,nstm) :: tmp_check3
3224
3225	!
3226	!
3227	! split 2d variables into 3d variables, per soil type
3228	!
3229	precisol_ns(:,:)=zero
3230	DO jv=1,nvm
3231	DO jst=1,nstm
3232	DO ji=1,kjpindex
3233	IF(veget(ji,jv).GT.min_sechiba) THEN
3234	precisol_ns(ji,jst)=precisol_ns(ji,jst)+precisol(ji,jv)* &
3235	& corr_veg_soil(ji,jv,jst) / veget(ji,jv)
3236	ENDIF
3237	END DO
3238	END DO
3239	END DO
3240	!
3241	!
3242	!
3243	vevapnu_old(:)=zero
3244	DO jst=1,nstm
3245	DO ji=1,kjpindex
3246	vevapnu_old(ji)=vevapnu_old(ji)+ &
3247	& ae_ns(ji,jst)soiltype(ji,jst)vegtot(ji)
3248	END DO
3249	END DO
3250	!
3251	!
3252	!
3253	DO jst=1,nstm
3254	DO ji=1,kjpindex
3255	IF (vevapnu_old(ji).GT.min_sechiba) THEN
3256	IF(evap_bare_lim(ji).GT.min_sechiba) THEN
3257	ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji)
3258	ELSE
3259	ae_ns(ji,jst)=ae_ns(ji,jst) * vevapnu(ji)/vevapnu_old(ji)
3260	ENDIF
3261	ELSEIF(veget(ji,1).GT.min_sechiba.AND.soiltype(ji,jst).GT.min_sechiba) THEN
3262	IF(evap_bare_lim(ji).GT.min_sechiba) THEN
3263	ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji)
3264	ELSE
3265	ae_ns(ji,jst)=vevapnu(ji)*corr_veg_soil(ji,1,jst)/veget(ji,1)
3266	ENDIF
3267	ENDIF
3268	precisol_ns(ji,jst)=precisol_ns(ji,jst)+MAX(-ae_ns(ji,jst),zero)
3269	END DO
3270	END DO
3271
3272	tr_ns(:,:)=zero
3273	DO jv=1,nvm
3274	DO jst=1,nstm
3275	DO ji=1,kjpindex
3276	IF (corr_veg_soil(ji,jv,jst).GT.min_sechiba.AND.humrel(ji,jv).GT.min_sechiba) THEN
3277	tr_ns(ji,jst)=tr_ns(ji,jst)+ corr_veg_soil(ji,jv,jst)humrelv(ji,jv,jst) &
3278	& transpir(ji,jv)/(humrel(ji,jv)*veget(ji,jv))
3279	ENDIF
3280	END DO
3281	END DO
3282	END DO
3283
3284	rootsink(:,:,:)=zero
3285	DO jv=1,nvm
3286	DO jsl=1,nslm
3287	DO jst=1,nstm
3288	DO ji=1,kjpindex
3289	IF ((humrel(ji,jv).GT.min_sechiba).AND.(corr_veg_soil(ji,jv,jst).GT.min_sechiba)) THEN
3290	rootsink(ji,jsl,jst) = rootsink(ji,jsl,jst) &
3291	& + corr_veg_soil(ji,jv,jst)* (transpir(ji,jv)*us(ji,jv,jst,jsl))/ &
3292	& (humrel(ji,jv)*veget(ji,jv))
3293	END IF
3294	END DO
3295	END DO
3296	END DO
3297	END DO
3298
3299	IF(check_cwrr) THEN
3300	DO jsl=1,nslm
3301	DO jst=1,nstm
3302	DO ji=1,kjpindex
3303	IF(mc(ji,jsl,jst).LT.-0.05) THEN
3304	WRITE(numout,*) 'CWRR split-----------------------------------------------'
3305	WRITE(numout,*) 'ji,jst,jsl',ji,jst,jsl
3306	WRITE(numout,*) 'mc',mc(ji,jsl,jst)
3307	WRITE(numout,*) 'rootsink,us',rootsink(ji,:,jst),us(ji,:,jst,jsl)
3308	WRITE(numout,*) 'corr_veg_soil',corr_veg_soil(ji,:,jst)
3309	WRITE(numout,*) 'transpir',transpir(ji,:)
3310	WRITE(numout,*) 'veget',veget(ji,:)
3311	WRITE(numout,*) 'humrel',humrel(ji,:)
3312	WRITE(numout,*) 'humrelv (pour ce jst)',humrelv(ji,:,jst)
3313	WRITE(numout,*) 'ae_ns',ae_ns(ji,jst)
3314	WRITE(numout,*) 'ae_ns',ae_ns(ji,jst)
3315	WRITE(numout,*) 'tr_ns',tr_ns(ji,jst)
3316	WRITE(numout,*) 'vevapnuold',vevapnu_old(ji)
3317	ENDIF
3318	END DO
3319	END DO
3320	END DO
3321	ENDIF
3322
3323
3324	! Now we check if the deconvolution is correct and conserves the fluxes:
3325
3326	IF (check_cwrr) THEN
3327
3328
3329	tmp_check1(:)=zero
3330	tmp_check2(:)=zero
3331
3332	! First we check the precisol and evapnu
3333
3334	DO jst=1,nstm
3335	DO ji=1,kjpindex
3336	tmp_check1(ji)=tmp_check1(ji) + &
3337	& (precisol_ns(ji,jst)-MAX(-ae_ns(ji,jst),zero))* &
3338	& soiltype(ji,jst)*vegtot(ji)
3339	END DO
3340	END DO
3341
3342	DO jv=1,nvm
3343	DO ji=1,kjpindex
3344	tmp_check2(ji)=tmp_check2(ji) + precisol(ji,jv)
3345	END DO
3346	END DO
3347
3348
3349	DO ji=1,kjpindex
3350
3351	IF(ABS(tmp_check1(ji)- tmp_check2(ji)).GT.allowed_err) THEN
3352	WRITE(numout,*) 'PRECISOL SPLIT FALSE:ji=',ji,tmp_check1(ji),tmp_check2(ji)
3353	WRITE(numout,*) 'vegtot',vegtot(ji)
3354
3355	DO jv=1,nvm
3356	WRITE(numout,*) 'jv,veget, precisol',jv,veget(ji,jv),precisol(ji,jv)
3357	DO jst=1,nstm
3358	WRITE(numout,*) 'corr_veg_soil:jst',jst,corr_veg_soil(ji,jv,jst)
3359	END DO
3360	END DO
3361
3362	DO jst=1,nstm
3363	WRITE(numout,*) 'jst,precisol_ns',jst,precisol_ns(ji,jst)
3364	WRITE(numout,*) 'soiltype', soiltype(ji,jst)
3365	END DO
3366	STOP 'in hydrol_split_soil check_cwrr'
3367	ENDIF
3368
3369	END DO
3370
3371
3372	tmp_check1(:)=zero
3373	tmp_check2(:)=zero
3374
3375	DO jst=1,nstm
3376	DO ji=1,kjpindex
3377	tmp_check1(ji)=tmp_check1(ji) + ae_ns(ji,jst)* &
3378	& soiltype(ji,jst)*vegtot(ji)
3379	END DO
3380	END DO
3381
3382	DO ji=1,kjpindex
3383
3384	IF(ABS(tmp_check1(ji)- vevapnu(ji)).GT.allowed_err) THEN
3385	WRITE(numout,*) 'VEVAPNU SPLIT FALSE:ji, Sum(ae_ns), vevapnu =',ji,tmp_check1(ji),vevapnu(ji)
3386	WRITE(numout,*) 'vegtot',vegtot(ji)
3387	WRITE(numout,*) 'evap_bare_lim, evap_bare_lim_ns',evap_bare_lim(ji), evap_bare_lim_ns(ji,:)
3388	WRITE(numout,*) 'vevapnu_old',vevapnu_old(ji)
3389	DO jst=1,nstm
3390	WRITE(numout,*) 'jst,ae_ns',jst,ae_ns(ji,jst)
3391	WRITE(numout,*) 'soiltype', soiltype(ji,jst)
3392	END DO
3393	STOP 'in hydrol_split_soil check_cwrr'
3394	ENDIF
3395	ENDDO
3396
3397	! Second we check the transpiration and root sink
3398
3399	tmp_check1(:)=zero
3400	tmp_check2(:)=zero
3401
3402
3403	DO jst=1,nstm
3404	DO ji=1,kjpindex
3405	tmp_check1(ji)=tmp_check1(ji) + tr_ns(ji,jst)* &
3406	& soiltype(ji,jst)*vegtot(ji)
3407	END DO
3408	END DO
3409
3410	DO jv=1,nvm
3411	DO ji=1,kjpindex
3412	tmp_check2(ji)=tmp_check2(ji) + transpir(ji,jv)
3413	END DO
3414	END DO
3415
3416	DO ji=1,kjpindex
3417
3418	IF(ABS(tmp_check1(ji)- tmp_check2(ji)).GT.allowed_err) THEN
3419	WRITE(numout,*) 'TRANSPIR SPLIT FALSE:ji=',ji,tmp_check1(ji),tmp_check2(ji)
3420	WRITE(numout,*) 'vegtot',vegtot(ji)
3421
3422	DO jv=1,nvm
3423	WRITE(numout,*) 'jv,veget, transpir',jv,veget(ji,jv),transpir(ji,jv)
3424	DO jst=1,nstm
3425	WRITE(numout,*) 'corr_veg_soil:ji,jv,jst',ji,jv,jst,corr_veg_soil(ji,jv,jst)
3426	END DO
3427	END DO
3428
3429	DO jst=1,nstm
3430	WRITE(numout,*) 'jst,tr_ns',jst,tr_ns(ji,jst)
3431	WRITE(numout,*) 'soiltype', soiltype(ji,jst)
3432	END DO
3433
3434	STOP 'in hydrol_split_soil check_cwrr'
3435	ENDIF
3436
3437	END DO
3438
3439
3440	tmp_check3(:,:)=zero
3441
3442	DO jst=1,nstm
3443	DO jsl=1,nslm
3444	DO ji=1,kjpindex
3445	tmp_check3(ji,jst)=tmp_check3(ji,jst) + rootsink(ji,jsl,jst)
3446	END DO
3447	END DO
3448	ENDDO
3449
3450	DO jst=1,nstm
3451	DO ji=1,kjpindex
3452	IF(ABS(tmp_check3(ji,jst)- tr_ns(ji,jst)).GT.allowed_err) THEN
3453	WRITE(numout,*) 'ROOTSINK SPLIT FALSE:ji,jst=', ji,jst,&
3454	& tmp_check3(ji,jst),tr_ns(ji,jst)
3455	WRITE(numout,*) 'HUMREL(jv=1:nvm)',humrel(ji,:)
3456	WRITE(numout,*) 'TRANSPIR',transpir(ji,:)
3457	DO jv=1,nvm
3458	WRITE(numout,*) 'jv=',jv,'us=',us(ji,jv,jst,:)
3459	ENDDO
3460	STOP 'in hydrol_split_soil check_cwrr'
3461	ENDIF
3462	END DO
3463	END DO
3464
3465	ENDIF
3466
3467	RETURN
3468
3469	END SUBROUTINE hydrol_split_soil
3470
3471	SUBROUTINE hydrol_diag_soil (kjpindex, veget, veget_max,soiltype, runoff, drainage, &
3472	& evap_bare_lim, evapot, vevapnu, returnflow, irrigation, &
3473	& shumdiag, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt)
3474	!
3475	! interface description
3476	! input scalar
3477	INTEGER(i_std), INTENT(in) :: kjpindex
3478	REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Map of vegetation types
3479	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget_max !! Max. vegetation type
3480	REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in) :: soiltype !! Map of soil types
3481	REAL(r_std), DIMENSION (kjpindex), INTENT (out) :: drysoil_frac !! Function of litter wetness
3482	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: runoff !! complete runoff
3483	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: drainage !! Drainage
3484	REAL(r_std), DIMENSION (kjpindex), INTENT(out) :: evap_bare_lim !!
3485	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: evapot !!
3486	REAL(r_std), DIMENSION (kjpindex), INTENT(inout) :: vevapnu
3487	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: returnflow !! Water returning to the deep reservoir
3488	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: irrigation !! Water from irrigation
3489	REAL(r_std), DIMENSION (kjpindex), INTENT(in) :: tot_melt
3490	REAL(r_std),DIMENSION (kjpindex,nbdl), INTENT (out) :: shumdiag !! relative soil moisture
3491	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: litterhumdiag !! litter humidity
3492	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out) :: humrel !! Relative humidity
3493	REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out) :: vegstress !! Veg. moisture stress (only for vegetation growth)
3494	!
3495	! local declaration
3496	!
3497	INTEGER(i_std) :: ji, jv, jsl, jst
3498	REAL(r_std), DIMENSION (kjpindex) :: mask_vegtot
3499	!
3500	! Put the prognostics variables of soil to zero if soiltype is zero
3501
3502	DO jst=1,nstm
3503
3504	DO ji=1,kjpindex
3505
3506	! IF(soiltype(ji,jst).EQ.zero) THEN
3507
3508	ae_ns(ji,jst) = ae_ns(ji,jst) * mask_soiltype(ji,jst)
3509	dr_ns(ji,jst) = dr_ns(ji,jst) * mask_soiltype(ji,jst)
3510	ru_ns(ji,jst) = ru_ns(ji,jst) * mask_soiltype(ji,jst)
3511	tmc(ji,jst) = tmc(ji,jst) * mask_soiltype(ji,jst)
3512
3513	DO jv=1,nvm
3514	humrelv(ji,jv,jst) = humrelv(ji,jv,jst) * mask_soiltype(ji,jst)
3515	DO jsl=1,nslm
3516	us(ji,jv,jst,jsl) = us(ji,jv,jst,jsl) * mask_soiltype(ji,jst)
3517	END DO
3518	END DO
3519
3520	DO jsl=1,nslm
3521	mc(ji,jsl,jst) = mc(ji,jsl,jst) * mask_soiltype(ji,jst)
3522	END DO
3523
3524	! ENDIF
3525
3526	END DO
3527	END DO
3528
3529
3530	runoff(:) = zero
3531	drainage(:) = zero
3532	humtot(:) = zero
3533	evap_bare_lim(:) = zero
3534	evap_bare_lim_ns(:,:) = zero
3535	shumdiag(:,:)= zero
3536	litterhumdiag(:) = zero
3537	tmc_litt_mea(:) = zero
3538	soilmoist(:,:) = zero
3539	humrel(:,:) = zero
3540	vegstress(:,:) = zero
3541	!
3542	! sum 3d variables in 2d variables with fraction of vegetation per soil type
3543	!
3544
3545
3546	DO ji = 1, kjpindex
3547	! WRITE(numout,*) ' kjpindex',kjpindex,ji,vegtot(ji)
3548	! mask_vegtot = MIN( un, MAX(zero,vegtot(ji) ) )
3549	mask_vegtot(ji) = 0
3550	IF(vegtot(ji) .GT. min_sechiba) THEN
3551	mask_vegtot(ji) = 1
3552	ENDIF
3553	END DO
3554
3555	DO ji = 1, kjpindex
3556	! WRITE(numout,*) 'vegtot,mask_vegtot',ji,vegtot(ji),mask_vegtot(ji)
3557	ae_ns(ji,:) = mask_vegtot(ji) * ae_ns(ji,:) * corr_veg_soil(ji,1,:)
3558	DO jst = 1, nstm
3559	drainage(ji) = mask_vegtot(ji) * (drainage(ji) + vegtot(ji)soiltype(ji,jst) dr_ns(ji,jst))
3560	runoff(ji) = mask_vegtot(ji) * (runoff(ji) + vegtot(ji)soiltype(ji,jst) ru_ns(ji,jst)) &
3561	& + (1 - mask_vegtot(ji)) * (tot_melt(ji) + irrigation(ji) + returnflow(ji))
3562	humtot(ji) = mask_vegtot(ji) * (humtot(ji) + soiltype(ji,jst) * tmc(ji,jst))
3563	END DO
3564	END DO
3565
3566	DO jst=1,nstm
3567	DO ji=1,kjpindex
3568	IF ((evapot(ji).GT.min_sechiba) .AND. &
3569	& (tmc_litter(ji,jst).GT.(tmc_litter_wilt(ji,jst)))) THEN
3570	evap_bare_lim_ns(ji,jst) = ae_ns(ji,jst) / evapot(ji)
3571	ELSEIF((evapot(ji).GT.min_sechiba).AND. &
3572	& (tmc_litter(ji,jst).GT.(tmc_litter_res(ji,jst)))) THEN
3573	evap_bare_lim_ns(ji,jst) = (un/deux) * ae_ns(ji,jst) / evapot(ji)
3574	END IF
3575
3576	END DO
3577	END DO
3578
3579	DO ji = 1, kjpindex
3580	evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)vegtot(ji)soiltype(ji,:))
3581	IF(evap_bare_lim(ji).GT.un + min_sechiba) THEN
3582	WRITE(numout,*) 'CWRR DIAG EVAP_BARE_LIM TOO LARGE', ji, &
3583	& evap_bare_lim(ji),evap_bare_lim_ns(ji,:)
3584	ENDIF
3585	!print *,'HYDROL_DIAG: ji,evap_bare_lim,evap_bare_lim_ns',ji,evap_bare_lim(ji),evap_bare_lim_ns(ji,:)
3586	ENDDO
3587	! we add the excess of snow sublimation to vevapnu
3588
3589	DO ji = 1,kjpindex
3590	vevapnu(ji) = vevapnu (ji) + subsinksoil(ji)*vegtot(ji)
3591	END DO
3592
3593	DO jst=1,nstm
3594	DO jv=1,nvm
3595	DO ji=1,kjpindex
3596	IF(veget_max(ji,jv).GT.min_sechiba) THEN
3597	vegstress(ji,jv)=vegstress(ji,jv)+vegstressv(ji,jv,jst)*soiltype(ji,jst) &
3598	& * corr_veg_soil_max(ji,jv,jst) *vegtot(ji)/veget_max(ji,jv)
3599	vegstress(ji,jv)= MAX(vegstress(ji,jv),zero)
3600	ENDIF
3601
3602	IF(veget(ji,jv).GT.min_sechiba) THEN
3603	humrel(ji,jv)=humrel(ji,jv)+humrelv(ji,jv,jst)*soiltype(ji,jst) &
3604	& * corr_veg_soil(ji,jv,jst)*vegtot(ji)/veget(ji,jv)
3605	humrel(ji,jv)=MAX(humrel(ji,jv),zero)
3606	ENDIF
3607	END DO
3608	END DO
3609	END DO
3610
3611	! vegstress(:,:) = humrel(:,:)
3612
3613
3614	DO jst=1,nstm
3615
3616	DO ji=1,kjpindex
3617	litterhumdiag(ji) = litterhumdiag(ji) + &
3618	& soil_wet_litter(ji,jst) * soiltype(ji,jst)
3619
3620	tmc_litt_mea(ji) = tmc_litt_mea(ji) + &
3621	& tmc_litter(ji,jst) * soiltype(ji,jst)
3622
3623	END DO
3624
3625
3626	DO jsl=1,nbdl
3627	DO ji=1,kjpindex
3628	shumdiag(ji,jsl)= shumdiag(ji,jsl) + soil_wet(ji,jsl,jst) * &
3629	& ((mcs(jst)-mcw(jst))/(mcf(jst)-mcw(jst))) * &
3630	& soiltype(ji,jst)
3631	soilmoist(ji,jsl)=soilmoist(ji,jsl)+mc(ji,jsl,jst)*soiltype(ji,jst)
3632	shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero)
3633	END DO
3634	END DO
3635
3636
3637	END DO
3638
3639
3640	DO ji=1,kjpindex
3641	drysoil_frac(ji) = un + MAX( MIN( (tmc_litt_dry_mea(ji) - tmc_litt_mea(ji)) / &
3642	& (tmc_litt_wet_mea(ji) - tmc_litt_dry_mea(ji)), zero), - un)
3643	END DO
3644
3645
3646
3647
3648
3649	END SUBROUTINE hydrol_diag_soil
3650
3651	!!
3652	!! This routines checks the water balance. First it gets the total
3653	!! amount of water and then it compares the increments with the fluxes.
3654	!! The computation is only done over the soil area as over glaciers (and lakes?)
3655	!! we do not have water conservation.
3656	!!
3657	!! This verification does not make much sense in REAL*4 as the precision is the same as some
3658	!! of the fluxes
3659	!!
3660	SUBROUTINE hydrol_waterbal (kjpindex, index, first_call, dtradia, veget, totfrac_nobio, &
3661	& qsintveg, snow,snow_nobio, precip_rain, precip_snow, returnflow, irrigation, tot_melt, &
3662	& vevapwet, transpir, vevapnu, vevapsno, runoff, drainage)
3663	!
3664	!
3665	!
3666	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
3667	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
3668	LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ?
3669	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
3670	!
3671	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: veget !! Fraction of vegetation type
3672	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: totfrac_nobio!! Total fraction of continental ice+lakes+...
3673	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
3674	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow mass [Kg/m^2]
3675	REAL(r_std), DIMENSION (kjpindex,nnobio), INTENT(in) :: snow_nobio !!Ice water balance
3676	!
3677	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
3678	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
3679	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: returnflow !! Water from irrigation
3680	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: irrigation !! Water from irrigation
3681	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: tot_melt !! Total melt
3682	!
3683	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: vevapwet !! Interception loss
3684	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: transpir !! Transpiration
3685	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapnu !! Bare soil evaporation
3686	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: vevapsno !! Snow evaporation
3687	REAL(r_std),DIMENSION (kjpindex), INTENT(in) :: runoff !! complete runoff
3688	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: drainage !! Drainage
3689	!
3690	! LOCAL
3691	!
3692	INTEGER(i_std) :: ji
3693	REAL(r_std) :: watveg, delta_water, tot_flux
3694	!
3695	!
3696	!
3697	IF ( first_call ) THEN
3698
3699	tot_water_beg(:) = zero
3700
3701	DO ji = 1, kjpindex
3702	watveg = SUM(qsintveg(ji,:))
3703	tot_water_beg(ji) = humtot(ji)*vegtot(ji) + watveg + snow(ji)&
3704	& + SUM(snow_nobio(ji,:))
3705	ENDDO
3706
3707	tot_water_end(:) = tot_water_beg(:)
3708
3709
3710	RETURN
3711
3712	ENDIF
3713	!
3714	! Check the water balance
3715	!
3716	tot_water_end(:) = zero
3717	!
3718	DO ji = 1, kjpindex
3719	!
3720	! If the fraction of ice, lakes, etc. does not complement the vegetation fraction then we do not
3721	! need to go any further
3722	!
3723	IF ( ABS(un - (totfrac_nobio(ji) + vegtot(ji))) .GT. allowed_err ) THEN
3724	WRITE(numout,*) 'HYDROL problem in vegetation or frac_nobio on point ', ji
3725	WRITE(numout,*) 'totfrac_nobio : ', totfrac_nobio(ji)
3726	WRITE(numout,*) 'vegetation fraction : ', vegtot(ji)
3727	STOP 'in hydrol_waterbal'
3728	ENDIF
3729	!
3730	watveg = SUM(qsintveg(ji,:))
3731	tot_water_end(ji) = humtot(ji)*vegtot(ji) + watveg + &
3732	& snow(ji) + SUM(snow_nobio(ji,:))
3733	!
3734	delta_water = tot_water_end(ji) - tot_water_beg(ji)
3735	!
3736	tot_flux = precip_rain(ji) + precip_snow(ji) + irrigation (ji) - &
3737	& SUM(vevapwet(ji,:)) - SUM(transpir(ji,:)) - vevapnu(ji) - vevapsno(ji) - &
3738	& runoff(ji) - drainage(ji) + returnflow(ji)
3739	!
3740	! Set some precision ! This is a wild guess and corresponds to what works on an IEEE machine
3741	! under double precision (REAL*8).
3742	!
3743	!
3744	IF ( ABS(delta_water-tot_flux) .GT. deux*allowed_err ) THEN
3745	WRITE(numout,*) '------------------------------------------------------------------------- '
3746	WRITE(numout,*) 'HYDROL does not conserve water. The erroneous point is : ', ji
3747	WRITE(numout,*) 'The error in mm/s is :', (delta_water-tot_flux)/dtradia, ' and in mm/dt : ', delta_water-tot_flux
3748	WRITE(numout,*) 'delta_water : ', delta_water, ' tot_flux : ', tot_flux
3749	WRITE(numout,*) 'Actual and allowed error : ', ABS(delta_water-tot_flux), allowed_err
3750	WRITE(numout,*) 'vegtot : ', vegtot(ji)
3751	WRITE(numout,*) 'precip_rain : ', precip_rain(ji)
3752	WRITE(numout,*) 'precip_snow : ', precip_snow(ji)
3753	WRITE(numout,*) 'Water from irrigation : ', returnflow(ji),irrigation(ji)
3754	WRITE(numout,*) 'Total water in soil :', humtot(ji)
3755	WRITE(numout,*) 'Water on vegetation :', watveg
3756	WRITE(numout,*) 'Snow mass :', snow(ji)
3757	WRITE(numout,*) 'Snow mass on ice :', SUM(snow_nobio(ji,:))
3758	WRITE(numout,*) 'Melt water :', tot_melt(ji)
3759	WRITE(numout,*) 'evapwet : ', vevapwet(ji,:)
3760	WRITE(numout,*) 'transpir : ', transpir(ji,:)
3761	WRITE(numout,*) 'evapnu, evapsno : ', vevapnu(ji), vevapsno(ji)
3762	WRITE(numout,*) 'drainage,runoff : ', drainage(ji),runoff(ji)
3763	STOP 'in hydrol_waterbal'
3764	ENDIF
3765	!
3766	ENDDO
3767	!
3768	! Transfer the total water amount at the end of the current timestep top the begining of the next one.
3769	!
3770	tot_water_beg = tot_water_end
3771	!
3772	END SUBROUTINE hydrol_waterbal
3773	!
3774	! This routine computes the changes in soil moisture and interception storage for the ALMA outputs
3775	!
3776	SUBROUTINE hydrol_alma (kjpindex, index, first_call, qsintveg, snow, snow_nobio, soilwet)
3777	!
3778	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
3779	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
3780	LOGICAL, INTENT (in) :: first_call !! At which time is this routine called ?
3781	!
3782	REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in) :: qsintveg !! Water on vegetation due to interception
3783	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: snow !! Snow water equivalent
3784	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: soilwet !! Soil wetness
3785	REAL(r_std),DIMENSION (kjpindex,nnobio), INTENT (in) :: snow_nobio !! Water balance on ice, lakes, .. [Kg/m^2]
3786	!
3787	! LOCAL
3788	!
3789	INTEGER(i_std) :: ji
3790	REAL(r_std) :: watveg
3791	!
3792	!
3793	!
3794	IF ( first_call ) THEN
3795
3796	tot_watveg_beg(:) = zero
3797	tot_watsoil_beg(:) = zero
3798	snow_beg(:) = zero
3799	!
3800	DO ji = 1, kjpindex
3801	watveg = SUM(qsintveg(ji,:))
3802	tot_watveg_beg(ji) = watveg
3803	tot_watsoil_beg(ji) = humtot(ji)
3804	snow_beg(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
3805	ENDDO
3806	!
3807	tot_watveg_end(:) = tot_watveg_beg(:)
3808	tot_watsoil_end(:) = tot_watsoil_beg(:)
3809	snow_end(:) = snow_beg(:)
3810
3811	RETURN
3812
3813	ENDIF
3814	!
3815	! Calculate the values for the end of the time step
3816	!
3817	tot_watveg_end(:) = zero
3818	tot_watsoil_end(:) = zero
3819	snow_end(:) = zero
3820	delintercept(:) = zero
3821	delsoilmoist(:) = zero
3822	delswe(:) = zero
3823	!
3824	DO ji = 1, kjpindex
3825	watveg = SUM(qsintveg(ji,:))
3826	tot_watveg_end(ji) = watveg
3827	tot_watsoil_end(ji) = humtot(ji)
3828	snow_end(ji) = snow(ji)+ SUM(snow_nobio(ji,:))
3829	!
3830	delintercept(ji) = tot_watveg_end(ji) - tot_watveg_beg(ji)
3831	delsoilmoist(ji) = tot_watsoil_end(ji) - tot_watsoil_beg(ji)
3832	delswe(ji) = snow_end(ji) - snow_beg(ji)
3833	!
3834	!
3835	ENDDO
3836	!
3837	!
3838	! Transfer the total water amount at the end of the current timestep top the begining of the next one.
3839	!
3840	tot_watveg_beg = tot_watveg_end
3841	tot_watsoil_beg = tot_watsoil_end
3842	snow_beg(:) = snow_end(:)
3843	!
3844	DO ji = 1,kjpindex
3845	soilwet(ji) = tot_watsoil_end(ji) / mx_eau_var(ji)
3846	ENDDO
3847	!
3848	END SUBROUTINE hydrol_alma
3849	!
3850	!
3851	END MODULE hydrol

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

Download in other formats: