Context Navigation

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

Last change on this file since 116 was 112, checked in by didier.solyga, 13 years ago
uncomment the division of throughfall_by_pft, could be wrong if division made after getin in the previous commit
File size: 157.5 KB

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

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

Download in other formats: