Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

hydrol.f90 @ 8

Last change on this file since 8 was 8, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 157.9 KB

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5/ORCHIDEE/src_sechiba/hydrol.f90 @ 8

Download in other formats: