Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/hydrol.f90 @ 311

Last change on this file since 311 was 311, checked in by didier.solyga, 13 years ago
Move and clean the rest of the externalized parameters from sechiba and stomate to src_parameters. Add two subroutines in constantes. Correct Olson type number 79 in vegcorr
File size: 157.1 KB

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

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

Download in other formats: