Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE/src_sechiba/hydrol.f90 @ 2061

Last change on this file since 2061 was 45, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

Property svn:keywords set to Revision Date HeadURL Date Author Revision

File size: 157.7 KB

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

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

Download in other formats: