Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE/src_sechiba/sechiba.f90 @ 112

Last change on this file since 112 was 106, checked in by didier.solyga, 14 years ago
Minor modifications
File size: 69.7 KB

Line
1	!! This module computes continental processes SECHIBA
2	!!
3	!! See also this graph <a href="http:shema.gif"><img SRC=" http:shema.gif" WIDTH="25" LENGTH="25"></A>
4	!!
5	!! @author Marie-Alice Foujols and Jan Polcher
6	!! @Version : $Revision: 1.46 $, $Date: 2010/05/07 08:28:13 $
7	!!
8	!! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_sechiba/sechiba.f90,v 1.46 2010/05/07 08:28:13 ssipsl Exp $
9	!! IPSL (2006)
10	!! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
11	!!
12	MODULE sechiba
13
14	! routines called : diffuco_main, enerbil_main, hydrolc_main, enrbil_fusion, condveg_main, thermosoil_main
15	!
16	USE ioipsl
17	!
18	! modules used :
19	USE constantes
20	USE pft_parameters
21	USE diffuco
22	USE condveg
23	USE enerbil
24	USE hydrol
25	USE hydrolc
26	USE thermosoil
27	USE sechiba_io
28	USE slowproc
29	USE routing
30	! USE write_field_p, only : WriteFieldI_p
31
32
33	IMPLICIT NONE
34
35	! public routines :
36	! sechiba_main
37	! sechiba_clear
38	PRIVATE
39	PUBLIC sechiba_main,sechiba_clear
40
41	! Index arrays we need internaly
42	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexveg !! indexing array for the 3D fields of vegetation
43	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexnobio !! indexing array for the 3D fields of other surfaces (ice, lakes, ...)
44	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexsoil !! indexing array for the 3D fields of soil types
45	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexgrnd !! indexing array for the 3D ground heat profiles
46	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexlayer !! indexing array for the 3D fields of soil layers
47	INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: indexalb !! indexing array for the 2 fields of albedo
48
49	! three dimensions array allocated, computed, saved and got in sechiba module
50
51	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:):: assim_param !! min+max+opt temps, vcmax, vjmax for photosynthesis
52
53	! two dimensions array allocated, computed, saved and got in sechiba module
54
55	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget !! Fraction of vegetation type
56	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: veget_max !! Max. fraction of vegetation type (LAI -> infty)
57	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: height !! Vegetation Height (m)
58	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: lai !! Surface foliere
59	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: gpp !! STOMATE: GPP. gC/m**2 of total area
60	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: humrel !! Relative humidity
61	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vegstress !! Vegetation moisture stress (only for vegetation growth)
62	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: soiltype !! Map of soil types, created in slowproc in the
63	!! order : silt, sand and clay.
64
65	!
66	! one dimension array allocated, computed and used in sechiba module and passed to other
67	! modules called
68
69	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta1 !! Snow resistance
70	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta4 !! Bare soil resistance
71	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilcap !! Soil calorific capacity
72	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: soilflx !! Soil flux
73	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: temp_sol !! Soil temperature
74	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: qsurf !! near soil air moisture
75	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow !! Snow mass [Kg/m^2]
76	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: snow_age !! Snow age
77	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drysoil_frac !! Fraction of visibly (albedo) Dry soil (Between 0 and 1)
78	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rsol !! resistance to bare soil evaporation
79	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evap_bare_lim !! Bare soil stress
80	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot !! Soil Potential Evaporation
81	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: evapot_corr !! Soil Potential Evaporation Correction (Milly 1992)
82	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapsno !! Snow evaporation
83	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vevapnu !! Bare soil evaporation
84	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: t2mdiag !! 2 meter temperature
85	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: tot_melt !! Total melt
86	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: vbeta !! Resistance coefficient
87	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: valpha !! Resistance coefficient
88	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: fusion !!
89	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: rau !! Density
90	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: deadleaf_cover !! Fraction of soil covered by dead leaves
91	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: co2_flux !! CO2 flux (gC/m**2 of average ground/s)
92	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: totfrac_nobio !! Total fraction of continental ice+lakes+cities+...
93	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: runoff !! Surface runoff generated by hydrol
94	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: drainage !! Deep drainage generated by hydrol
95	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: returnflow !! Routed water which returns into the soil
96	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: irrigation !! irrigation going back into the soils
97	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: emis !! Surface emissivity
98	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: z0 !! Surface roughness
99	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: roughheight !! Effective height for roughness
100	!
101	! Arrays which are diagnostics from the physical processes for
102	! for the biological processes. They are not saved in the restart file because at the first time step,
103	! they are recalculated. However, they must be saved in memory as they are in slowproc which is called
104	! before the modules which calculate them.
105	!
106	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: shumdiag !! Relative soil moisture
107	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:) :: litterhumdiag!! litter humidity
108	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: stempdiag !! Temperature which controls canopy evolution
109
110	! two dimensions array allocated, computed and used in sechiba module and passed to other
111	! modules called
112
113	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintveg !! Water on vegetation due to interception
114	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta2 !! Interception resistance
115	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbeta3 !! Vegetation resistance
116	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vbetaco2 !! STOMATE: Vegetation resistance to CO2
117	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: cimean !! STOMATE: Mean ci
118	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: vevapwet !! Interception
119	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: transpir !! Transpiration
120	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: qsintmax !! Maximum water on vegetation for interception
121	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: rveget !! Vegetation resistance
122	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio !! Water balance over other surface types (that is snow !)
123	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: snow_nobio_age !! Snow age on other surface types
124	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: frac_nobio !! Fraction of continental ice, lakes, ...
125	REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:) :: albedo !! Surface albedo for visible and NIR
126	!
127	! variables used inside sechiba module : declaration and initialisation
128	!
129	LOGICAL, SAVE :: l_first_sechiba = .TRUE.!! Initialisation has to be done one time
130	CHARACTER(LEN=80) , SAVE :: var_name !! To store variables names for I/O
131
132	LOGICAL, SAVE :: river_routing !! Flag that decides if we route.
133	LOGICAL, SAVE :: hydrol_cwrr !! Selects the CWRR hydrology.
134
135	LOGICAL, SAVE :: myfalse=.FALSE.
136	LOGICAL, SAVE :: mytrue=.TRUE.
137
138	CONTAINS
139
140	!! Main routine for sechiba module.
141	!!
142	!! Should be called two times:
143	!! - first time to have initial values
144	!! - second time to have complete algorithm
145	!!
146	!! Algorithm:
147	!! 3 series of called SECHIBA processes
148	!! - initialisation (first time only)
149	!! - time step evolution (every time step)
150	!! - prepares output and storage of restart arrays (last time only)
151	!!
152	!! One serie consists of:
153	!! - call sechiba_var_init to do some initialisation
154	!! - call slowproc_main to do some daily initialisation
155	!! - call diffuco_main for diffusion coefficient calculation
156	!! - call enerbil_main for energy bilan calculation
157	!! - call hydrolc_main for hydrologic processes calculation
158	!! - call enerbil_fusion : last part with fusion
159	!! - call condveg_main for surface conditions such as roughness, albedo, and emmisivity
160	!! - call thermosoil_main for soil thermodynamic calculation
161	!! - call sechiba_end to swap new fields to previous
162	!!
163	!! @call sechiba_var_init
164	!! @call slowproc_main
165	!! @call diffuco_main
166	!! @call enerbil_main
167	!! @call hydrolc_main
168	!! @call enerbil_fusion
169	!! @call condveg_main
170	!! @call thermosoil_main
171	!! @call sechiba_end
172	!!
173	SUBROUTINE sechiba_main (kjit, kjpij, kjpindex, index, dtradia, date0, &
174	& ldrestart_read, ldrestart_write, control_in, &
175	& lalo, contfrac, neighbours, resolution,&
176	! First level conditions
177	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
178	! & zlev, u, v, qair, temp_air, epot_air, ccanopy, &
179	& zlev, u, v, qair, q2m, t2m, temp_air, epot_air, ccanopy, &
180	! Variables for the implicit coupling
181	& tq_cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
182	! Rain, snow, radiation and surface pressure
183	& precip_rain, precip_snow, lwdown, swnet, swdown, pb, &
184	! Output : Fluxes
185	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
186	! Surface temperatures and surface properties
187	& tsol_rad, temp_sol_new, qsurf_out, albedo_out, emis_out, z0_out, &
188	! File ids
189	& rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
190
191	! interface description for dummy arguments
192	! input scalar
193	INTEGER(i_std), INTENT(in) :: kjit !! Time step number
194	INTEGER(i_std), INTENT(in) :: kjpij !! Total size of size. This is the un-compressed grid
195	INTEGER(i_std), INTENT(in) :: kjpindex !! Domain size
196	INTEGER(i_std),INTENT (in) :: rest_id !! _Restart_ file identifier
197	INTEGER(i_std),INTENT (in) :: hist_id !! _History_ file identifier
198	INTEGER(i_std),INTENT (in) :: hist2_id !! _History_ file 2 identifier
199	INTEGER(i_std),INTENT (in) :: rest_id_stom !! STOMATE's _Restart_ file identifier
200	INTEGER(i_std),INTENT (in) :: hist_id_stom !! STOMATE's _History_ file identifier
201	INTEGER(i_std),INTENT(in) :: hist_id_stom_IPCC !! STOMATE's IPCC _history_ file file identifier
202	REAL(r_std), INTENT (in) :: dtradia !! Time step in seconds
203	REAL(r_std), INTENT (in) :: date0 !! initial date
204	LOGICAL, INTENT(in) :: ldrestart_read !! Logical for _restart_ file to read
205	LOGICAL, INTENT(in) :: ldrestart_write !! Logical for _restart_ file to write
206	TYPE(control_type), INTENT(in) :: control_in !! Flags that (de)activate parts of the model
207	! input fields
208	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
209	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: contfrac !! Fraction of continent in the grid
210	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
211	!
212	INTEGER(i_std), DIMENSION (kjpindex,8), INTENT(in) :: neighbours !! neighoring grid points if land
213	REAL(r_std), DIMENSION (kjpindex,2), INTENT(in) :: resolution !! size in x an y of the grid (m)
214	!
215	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: u !! Lowest level wind speed
216	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: v !! Lowest level wind speed
217	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: zlev !! Height of first layer
218	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: qair !! Lowest level specific humidity
219	! Ajout Nathalie - Juin 2006 - Q2M/t2m pour calcul Rveget
220	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: q2m !! 2m specific humidity
221	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: t2m !! 2m air temperature
222	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_rain !! Rain precipitation
223	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: precip_snow !! Snow precipitation
224	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: lwdown !! Down-welling long-wave flux
225	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swnet !! Net surface short-wave flux
226	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: swdown !! Down-welling surface short-wave flux
227	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature in Kelvin
228	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: epot_air !! Air potential energy
229	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: ccanopy !! CO2 concentration in the canopy
230	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petAcoef !! Coeficients A from the PBL resolution
231	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqAcoef !! One for T and another for q
232	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: petBcoef !! Coeficients B from the PBL resolution
233	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: peqBcoef !! One for T and another for q
234	REAL(r_std),DIMENSION (kjpindex), INTENT (inout) :: tq_cdrag !! This is the cdrag without the wind multiplied
235	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
236	! output fields
237	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: coastalflow !! Diffuse water flow to the oceans
238	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: riverflow !! River outflow to the oceans
239	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: tsol_rad !! Radiative surface temperature
240	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: vevapp !! Total of evaporation
241	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol_new !! New soil temperature
242	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: qsurf_out !! Surface specific humidity
243	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: z0_out !! Surface roughness (output diagnostic)
244	REAL(r_std),DIMENSION (kjpindex,2), INTENT (out) :: albedo_out !! Albedo (output diagnostic)
245	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxsens !! Sensible chaleur flux
246	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: fluxlat !! Latent chaleur flux
247	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: emis_out !! Emissivity
248
249	REAL(r_std), ALLOCATABLE, DIMENSION (:) :: runoff1,drainage1, soilcap1,soilflx1
250	REAL(r_std), ALLOCATABLE, DIMENSION (:,:) :: shumdiag1
251
252	REAL(r_std), DIMENSION(kjpindex) :: histvar !! computations for history files
253
254	INTEGER(i_std) :: j
255
256
257
258
259	IF (long_print) WRITE(numout,*) ' kjpindex =',kjpindex
260
261	! do SECHIBA'S first call initialisation
262
263	IF (l_first_sechiba) THEN
264
265	CALL sechiba_init (kjit, ldrestart_read, kjpij, kjpindex, index, rest_id, control_in, lalo)
266
267	ALLOCATE(runoff1 (kjpindex),drainage1 (kjpindex), soilcap1 (kjpindex),soilflx1 (kjpindex))
268	ALLOCATE(shumdiag1(kjpindex,nbdl))
269	!
270	! computes initialisation of energy bilan
271	!
272	IF (ldrestart_read) THEN
273
274	IF (long_print) WRITE (numout,*) ' we have to read a restart file for SECHIBA variables'
275	var_name='soilcap' ;
276	CALL ioconf_setatt('UNITS', '-')
277	CALL ioconf_setatt('LONG_NAME','Soil calorific capacity')
278	soilcap1=val_exp
279	IF ( ok_var(var_name) ) THEN
280	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., soilcap1, "gather", nbp_glo, index_g)
281	IF (MINVAL(soilcap1) < MAXVAL(soilcap1) .OR. MAXVAL(soilcap1) < val_exp) THEN
282	soilcap(:) = soilcap1(:)
283	ENDIF
284	ENDIF
285	!
286	var_name='soilflx' ;
287	CALL ioconf_setatt('UNITS', '-')
288	CALL ioconf_setatt('LONG_NAME','Soil flux')
289	soilflx1=val_exp
290	IF ( ok_var(var_name) ) THEN
291	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., soilflx1, "gather", nbp_glo, index_g)
292	IF (MINVAL(soilflx1) < MAXVAL(soilflx1) .OR. MAXVAL(soilflx1) < val_exp) THEN
293	soilflx(:) = soilflx1(:)
294	ENDIF
295	ENDIF
296	!
297	var_name='shumdiag' ;
298	CALL ioconf_setatt('UNITS', '-')
299	CALL ioconf_setatt('LONG_NAME','Relative soil moisture')
300	shumdiag1=val_exp
301	IF ( ok_var(var_name) ) THEN
302	CALL restget_p (rest_id, var_name, nbp_glo, nbdl, 1, kjit, .TRUE., shumdiag1, "gather", nbp_glo, index_g)
303	IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
304	shumdiag(:,:) = shumdiag1(:,:)
305	ENDIF
306	ENDIF
307	ENDIF
308
309	!
310	! computes slow variables
311	!
312	CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
313	ldrestart_read, ldrestart_write, control%ok_co2, control%ok_stomate, &
314	index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
315	t2mdiag, t2mdiag, temp_sol, stempdiag, &
316	vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
317	deadleaf_cover, &
318	assim_param, &
319	lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
320	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
321	co2_flux)
322	!
323	! computes initialisation of diffusion coeff
324	!
325
326	CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
327	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
328	! & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb , &
329	& zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb , &
330	& rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
331	& swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
332	& vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
333	!
334	! computes initialisation of energy bilan
335	!
336	CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
337	& index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef,&
338	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
339	& cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
340	& vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
341	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id)
342	!
343	! computes initialisation of hydrologie
344	!
345	IF (ldrestart_read) THEN
346
347	var_name='runoff' ;
348	CALL ioconf_setatt('UNITS', 'mm/d')
349	CALL ioconf_setatt('LONG_NAME','Complete runoff')
350	runoff1=val_exp
351	IF ( ok_var(var_name) ) THEN
352	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff1, "gather", nbp_glo, index_g)
353	IF (MINVAL(runoff1) < MAXVAL(runoff1) .OR. MAXVAL(runoff1) < val_exp) THEN
354	runoff(:) = runoff1(:)
355	ENDIF
356	ENDIF
357
358	var_name='drainage' ;
359	CALL ioconf_setatt('UNITS', 'mm/d')
360	CALL ioconf_setatt('LONG_NAME','Deep drainage')
361	drainage1=val_exp
362	IF ( ok_var(var_name) ) THEN
363	CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage1, "gather", nbp_glo, index_g)
364	IF (MINVAL(drainage1) < MAXVAL(drainage1) .OR. MAXVAL(drainage1) < val_exp) THEN
365	drainage(:) = drainage1(:)
366	ENDIF
367	ENDIF
368
369	IF ( ok_var("shumdiag") ) THEN
370	IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
371	shumdiag(:,:) = shumdiag1(:,:)
372	ENDIF
373	ENDIF
374	ENDIF
375	!
376	IF ( .NOT. hydrol_cwrr ) THEN
377	CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
378	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
379	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
380	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
381	& vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id)
382	evap_bare_lim(:) = -un
383	ELSE
384	CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, indexsoil, indexlayer, &
385	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
386	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
387	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
388	& vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
389	& shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
390	ENDIF
391	IF (ldrestart_read) THEN
392
393	IF ( ok_var("runoff") ) THEN
394	IF (MINVAL(runoff1) < MAXVAL(runoff1) .OR. MAXVAL(runoff1) < val_exp) THEN
395	runoff(:) = runoff1(:)
396	ENDIF
397	ENDIF
398
399	IF ( ok_var("drainage") ) THEN
400	IF (MINVAL(drainage1) < MAXVAL(drainage1) .OR. MAXVAL(drainage1) < val_exp) THEN
401	drainage(:) = drainage1(:)
402	ENDIF
403	ENDIF
404
405	IF ( ok_var("shumdiag") ) THEN
406	IF (MINVAL(shumdiag1) < MAXVAL(shumdiag1) .OR. MAXVAL(shumdiag1) < val_exp) THEN
407	shumdiag(:,:) = shumdiag1(:,:)
408	ENDIF
409	ENDIF
410	ENDIF
411
412	!
413	! computes initialisation of condveg
414	!
415	CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
416	& lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
417	& zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
418	& drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
419	!
420	! computes initialisation of Soil Thermodynamic
421	!
422	CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
423	& index, indexgrnd, temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
424
425	IF (ldrestart_read) THEN
426	IF ( ok_var("soilcap") ) THEN
427	IF (MINVAL(soilcap1) < MAXVAL(soilcap1) .OR. MAXVAL(soilcap1) < val_exp) THEN
428	soilcap(:) = soilcap1(:)
429	ENDIF
430	ENDIF
431	!
432	IF ( ok_var("soilflx") ) THEN
433	IF (MINVAL(soilflx1) < MAXVAL(soilflx1) .OR. MAXVAL(soilflx1) < val_exp) THEN
434	soilflx(:) = soilflx1(:)
435	ENDIF
436	ENDIF
437	ENDIF
438
439	!
440	! If we chose to route the water then we call the module. Else variables
441	! are set to zero.
442	!
443	!
444	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
445	CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
446	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
447	& drainage, evapot_corr, precip_rain, humrel, &
448	& stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
449	ELSE
450	riverflow(:) = zero
451	coastalflow(:) = zero
452	returnflow(:) = zero
453	irrigation(:) = zero
454	ENDIF
455	!
456	! Write the internal variables into the output fields
457	!
458	z0_out(:) = z0(:)
459	emis_out(:) = emis(:)
460	albedo_out(:,:) = albedo(:,:)
461	qsurf_out(:) = qsurf(:)
462
463	DEALLOCATE(runoff1,drainage1,soilcap1,soilflx1)
464	DEALLOCATE(shumdiag1)
465
466	!
467	! This line should remain last as it ends the initialisation and returns to the caling
468	! routine.
469	!
470	RETURN
471	!
472	ENDIF
473
474	!
475	! computes some initialisation every SECHIBA's call
476	!
477	CALL sechiba_var_init (kjpindex, rau, pb, temp_air)
478
479	!
480	! computes diffusion coeff
481	!
482	CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, u, v, &
483	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
484	! & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb , &
485	& zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb , &
486	& rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
487	& swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
488	& vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
489
490	!
491	! computes energy bilan
492	!
493
494	CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, &
495	& index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef, &
496	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
497	& cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
498	& vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
499	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id)
500	!
501	! computes hydrologie
502	!
503	IF ( .NOT. hydrol_cwrr ) THEN
504	CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, &
505	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
506	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
507	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
508	& vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, rest_id, hist_id, hist2_id)
509	evap_bare_lim(:) = -un
510	ELSE
511	CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexveg, indexsoil, indexlayer, &
512	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
513	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
514	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
515	& vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
516	& shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
517	ENDIF
518	!
519	! computes last part of energy bilan
520	!
521	CALL enerbil_fusion (kjpindex, dtradia, tot_melt, soilcap, temp_sol_new, fusion)
522
523	!
524	! computes condveg
525	!
526	CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index,&
527	& lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
528	& zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
529	& drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
530	!
531	! computes Soil Thermodynamic
532	!
533	CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, indexgrnd, &
534	& temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
535
536	!
537	! If we chose to route the water then we call the module. Else variables
538	! are set to zero.
539	!
540	!
541	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
542	CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, myfalse, index, &
543	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
544	& drainage, evapot_corr, precip_rain, humrel, &
545	& stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
546	! returnflow(:) = returnflow(:) * 100.
547	ELSE
548	riverflow(:) = zero
549	coastalflow(:) = zero
550	returnflow(:) = zero
551	irrigation(:) = zero
552	ENDIF
553
554
555	!
556	! computes slow variables
557	! ok_co2 and ok_stomate are interpreted as flags that determine whether the
558	! forcing files are written.
559	!
560	CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
561	ldrestart_read, myfalse, control%ok_co2, control%ok_stomate, &
562	index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
563	t2mdiag, t2mdiag, temp_sol, stempdiag, &
564	vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
565	deadleaf_cover, &
566	assim_param, &
567	lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
568	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
569	co2_flux)
570
571	!
572	! call swap from new computed variables
573	!
574	CALL sechiba_end (kjpindex, dtradia, temp_sol, temp_sol_new)
575	!
576	! Write the internal variables into the output fields
577	!
578	z0_out(:) = z0(:)
579	emis_out(:) = emis(:)
580	albedo_out(:,:) = albedo(:,:)
581	qsurf_out(:) = qsurf(:)
582	!
583	! Writing the global variables on the history tape
584	!
585	!
586	IF ( .NOT. almaoutput ) THEN
587	CALL histwrite(hist_id, 'beta', kjit, vbeta, kjpindex, index)
588	CALL histwrite(hist_id, 'z0', kjit, z0, kjpindex, index)
589	CALL histwrite(hist_id, 'roughheight', kjit, roughheight, kjpindex, index)
590	CALL histwrite(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
591	CALL histwrite(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
592	CALL histwrite(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
593	CALL histwrite(hist_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
594	CALL histwrite(hist_id, 'subli', kjit, vevapsno, kjpindex, index)
595	CALL histwrite(hist_id, 'evapnu', kjit, vevapnu, kjpindex, index)
596	CALL histwrite(hist_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
597	CALL histwrite(hist_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
598	CALL histwrite(hist_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
599	CALL histwrite(hist_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
600	CALL histwrite(hist_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
601	CALL histwrite(hist_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
602	CALL histwrite(hist_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
603	CALL histwrite(hist_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
604	CALL histwrite(hist_id, 'rsol', kjit, rsol, kjpindex, index)
605	CALL histwrite(hist_id, 'snow', kjit, snow, kjpindex, index)
606	CALL histwrite(hist_id, 'snowage', kjit, snow_age, kjpindex, index)
607	CALL histwrite(hist_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
608	CALL histwrite(hist_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
609	CALL histwrite(hist_id, 'soiltype', kjit, soiltype, kjpindex*nstm, indexsoil)
610	IF ( control%ok_co2 ) THEN
611	CALL histwrite(hist_id, 'vbetaco2', kjit, vbetaco2, kjpindex*nvm, indexveg)
612	CALL histwrite(hist_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
613	CALL histwrite(hist_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
614	ENDIF
615	IF ( control%ok_stomate ) THEN
616	CALL histwrite(hist_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)
617	ENDIF
618
619	histvar(:)=SUM(vevapwet(:,:),dim=2)/86400
620	CALL histwrite(hist_id, 'evspsblveg', kjit, histvar, kjpindex, index)
621
622	histvar(:)=(vevapnu(:)+vevapsno(:))/86400
623	CALL histwrite(hist_id, 'evspsblsoi', kjit, histvar, kjpindex, index)
624
625	histvar(:)=SUM(transpir(:,:),dim=2)/86400
626	CALL histwrite(hist_id, 'tran', kjit, histvar, kjpindex, index)
627
628	!------------------------------------
629
630	! histvar(:)=SUM(veget_max(:,2:9),dim=2)100contfrac(:)
631	! CALL histwrite(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
632
633	! histvar(:)=SUM(veget_max(:,10:11),dim=2)100contfrac(:)
634	! CALL histwrite(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
635
636	! histvar(:)=SUM(veget_max(:,12:13),dim=2)100contfrac(:)
637	! CALL histwrite(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
638
639	!$$ 25/10/10 Modif DS & NViovy
640
641	!!$ DO j = 2, nvm
642	!!$ IF (is_tree(j) .AND. natural(j)) THEN
643	!!$ histvar(:)=SUM(veget_max(:,j))100contfrac(:)
644	!!$ CALL histwrite(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
645	!!$ ELSE IF ((.NOT. is_tree(j)) .AND. natural(j)) THEN
646	!!$ histvar(:)=SUM(veget_max(:,j))100contfrac(:)
647	!!$ CALL histwrite(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
648	!!$ ELSE IF (.NOT. natural(j)) THEN
649	!!$ histvar(:)=SUM(veget_max(:,j))100contfrac(:)
650	!!$ CALL histwrite(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
651	!!$ ELSE
652	!!$ WRITE(numout,*) 'Something wrong happened in histvar of veget_max in sechiba.f90'
653	!!$ ENDIF
654	!!$ ENDDO
655	!!$
656	WHERE (is_tree(:) .AND. natural(:))
657	histvar(:)=SUM(veget_max(:,j))100contfrac(:)
658	ENDWHERE
659	CALL histwrite(hist_id, 'treeFrac', kjit, histvar, kjpindex, index)
660
661	WHERE (.NOT. is_tree(:) .AND. natural(:))
662	histvar(:)=SUM(veget_max(:,j))100contfrac(:)
663	ENDWHERE
664	CALL histwrite(hist_id, 'grassFrac', kjit, histvar, kjpindex, index)
665
666	WHERE (.NOT. natural(:))
667	histvar(:)=SUM(veget_max(:,j))100contfrac(:)
668	ENDWHERE
669	CALL histwrite(hist_id, 'cropFrac', kjit, histvar, kjpindex, index)
670	!!$
671
672	!--------------------------------------------
673	histvar(:)=veget_max(:,1)100contfrac(:)
674	CALL histwrite(hist_id, 'baresoilFrac', kjit, histvar, kjpindex, index)
675
676	histvar(:)=SUM(frac_nobio(:,1:nnobio),dim=2)100contfrac(:)
677	CALL histwrite(hist_id, 'residualFrac', kjit, histvar, kjpindex, index)
678
679	ELSE
680	CALL histwrite(hist_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
681	CALL histwrite(hist_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
682	CALL histwrite(hist_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
683	CALL histwrite(hist_id, 'Qf', kjit, fusion, kjpindex, index)
684	CALL histwrite(hist_id, 'SWE', kjit, snow, kjpindex, index)
685	CALL histwrite(hist_id, 'ESoil', kjit, vevapnu, kjpindex, index)
686	CALL histwrite(hist_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
687	CALL histwrite(hist_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
688	CALL histwrite(hist_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
689	CALL histwrite(hist_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
690	IF ( control%ok_co2 ) THEN
691	CALL histwrite(hist_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
692	ENDIF
693	IF ( control%ok_stomate ) THEN
694	CALL histwrite(hist_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
695	ENDIF
696	ENDIF
697	!
698	IF ( hist2_id > 0 ) THEN
699	IF ( .NOT. almaoutput ) THEN
700	CALL histwrite(hist2_id, 'tsol_rad', kjit, tsol_rad, kjpindex, index)
701	CALL histwrite(hist2_id, 'qsurf', kjit, qsurf, kjpindex, index)
702	CALL histwrite(hist2_id, 'albedo', kjit, albedo, kjpindex*2, indexalb)
703	CALL histwrite(hist2_id, 'emis', kjit, emis, kjpindex, index)
704	!
705	CALL histwrite(hist2_id, 'beta', kjit, vbeta, kjpindex, index)
706	CALL histwrite(hist2_id, 'z0', kjit, z0, kjpindex, index)
707	CALL histwrite(hist2_id, 'roughheight', kjit, roughheight, kjpindex, index)
708	CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
709	CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
710	CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
711	CALL histwrite(hist2_id, 'lai', kjit, lai, kjpindex*nvm, indexveg)
712	CALL histwrite(hist2_id, 'subli', kjit, vevapsno, kjpindex, index)
713	CALL histwrite(hist2_id, 'vevapnu', kjit, vevapnu, kjpindex, index)
714	CALL histwrite(hist2_id, 'transpir', kjit, transpir, kjpindex*nvm, indexveg)
715	CALL histwrite(hist2_id, 'inter', kjit, vevapwet, kjpindex*nvm, indexveg)
716	CALL histwrite(hist2_id, 'vbeta1', kjit, vbeta1, kjpindex, index)
717	CALL histwrite(hist2_id, 'vbeta2', kjit, vbeta2, kjpindex*nvm, indexveg)
718	CALL histwrite(hist2_id, 'vbeta3', kjit, vbeta3, kjpindex*nvm, indexveg)
719	CALL histwrite(hist2_id, 'vbeta4', kjit, vbeta4, kjpindex, index)
720	CALL histwrite(hist2_id, 'drysoil_frac', kjit, drysoil_frac, kjpindex, index)
721	CALL histwrite(hist2_id, 'rveget', kjit, rveget, kjpindex*nvm, indexveg)
722	CALL histwrite(hist2_id, 'rsol', kjit, rsol, kjpindex, index)
723	CALL histwrite(hist2_id, 'snow', kjit, snow, kjpindex, index)
724	CALL histwrite(hist2_id, 'snowage', kjit, snow_age, kjpindex, index)
725	CALL histwrite(hist2_id, 'snownobio', kjit, snow_nobio, kjpindex*nnobio, indexnobio)
726	CALL histwrite(hist2_id, 'snownobioage', kjit, snow_nobio_age, kjpindex*nnobio, indexnobio)
727	CALL histwrite(hist2_id, 'soiltype', kjit, soiltype, kjpindex*nstm, indexsoil)
728	IF ( control%ok_co2 ) THEN
729	CALL histwrite(hist2_id, 'vbetaco2', kjit, vbetaco2, kjpindex*nvm, indexveg)
730	CALL histwrite(hist2_id, 'gpp', kjit, gpp, kjpindex*nvm, indexveg)
731	CALL histwrite(hist2_id, 'cimean', kjit, cimean, kjpindex*nvm, indexveg)
732	ENDIF
733	IF ( control%ok_stomate ) THEN
734	CALL histwrite(hist2_id, 'nee', kjit, co2_flux, kjpindex*nvm, indexveg)
735	ENDIF
736	ELSE
737	CALL histwrite(hist2_id, 'vegetfrac', kjit, veget, kjpindex*nvm, indexveg)
738	CALL histwrite(hist2_id, 'maxvegetfrac', kjit, veget_max, kjpindex*nvm, indexveg)
739	CALL histwrite(hist2_id, 'nobiofrac', kjit, frac_nobio, kjpindex*nnobio, indexnobio)
740	CALL histwrite(hist2_id, 'Qf', kjit, fusion, kjpindex, index)
741	CALL histwrite(hist2_id, 'SWE', kjit, snow, kjpindex, index)
742	CALL histwrite(hist2_id, 'ESoil', kjit, vevapnu, kjpindex, index)
743	CALL histwrite(hist2_id, 'TVeg', kjit, transpir, kjpindex*nvm, indexveg)
744	CALL histwrite(hist2_id, 'ECanop', kjit, vevapwet, kjpindex*nvm, indexveg)
745	CALL histwrite(hist2_id, 'ACond', kjit, tq_cdrag, kjpindex, index)
746	CALL histwrite(hist2_id, 'SnowFrac', kjit, vbeta1, kjpindex, index)
747	IF ( control%ok_co2 ) THEN
748	CALL histwrite(hist2_id, 'GPP', kjit, gpp, kjpindex*nvm, indexveg)
749	ENDIF
750	IF ( control%ok_stomate ) THEN
751	CALL histwrite(hist2_id, 'NEE', kjit, co2_flux, kjpindex*nvm, indexveg)
752	ENDIF
753	ENDIF
754	ENDIF
755	!
756	! prepares restart file for the next simulation from SECHIBA and from other modules
757	!
758	IF (ldrestart_write) THEN
759
760	IF (long_print) WRITE (numout,*) ' we have to write a restart file '
761
762	!
763	! call diffuco_main to write restart files
764	!
765
766	CALL diffuco_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, u, v, &
767	! Ajout Nathalie - Juin 2006 - passage q2m/t2m pour calcul Rveget
768	! & zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, pb , &
769	& zlev, z0, roughheight, temp_sol, temp_air, rau, tq_cdrag, qsurf, qair, q2m, t2m, pb , &
770	& rsol, evap_bare_lim, evapot, snow, frac_nobio, snow_nobio, totfrac_nobio, &
771	& swnet, swdown, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
772	& vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, rveget, cimean, rest_id, hist_id, hist2_id)
773
774	!
775	! call energy bilan to write restart files
776	!
777	CALL enerbil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, &
778	& index, zlev, lwdown, swnet, epot_air, temp_air, u, v, petAcoef, petBcoef,&
779	& qair, peqAcoef, peqBcoef, pb, rau, vbeta, valpha, vbeta1, vbeta2, vbeta3, vbeta4, vbetaco2, &
780	& cimean, ccanopy, emis, soilflx, soilcap, tq_cdrag, humrel, fluxsens, fluxlat, &
781	& vevapp, transpir, gpp, vevapnu, vevapwet, vevapsno, t2mdiag, temp_sol, tsol_rad, &
782	& temp_sol_new, qsurf, evapot, evapot_corr, rest_id, hist_id, hist2_id)
783
784	!
785	! call hydrologie to write restart files
786	!
787	IF ( .NOT. hydrol_cwrr ) THEN
788	CALL hydrolc_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, &
789	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
790	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
791	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, &
792	& humrel, vegstress, rsol, drysoil_frac, evapot, evapot_corr, shumdiag, litterhumdiag, soilcap, &
793	& rest_id, hist_id, hist2_id)
794	evap_bare_lim(:) = -un
795	ELSE
796	CALL hydrol_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexveg, indexsoil, indexlayer, &
797	& temp_sol_new, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max,&
798	& qsintmax, qsintveg, vevapnu, vevapsno, snow, snow_age, snow_nobio, snow_nobio_age,&
799	& tot_melt, transpir, precip_rain, precip_snow, returnflow, irrigation, humrel, &
800	& vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, &
801	& shumdiag, litterhumdiag, soilcap, soiltype, rest_id, hist_id, hist2_id)
802	rsol(:) = -un
803	ENDIF
804	!
805	! call condveg to write restart files
806	!
807	CALL condveg_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
808	& lalo, neighbours, resolution, contfrac, veget, veget_max, frac_nobio, totfrac_nobio, &
809	& zlev, snow, snow_age, snow_nobio, snow_nobio_age, &
810	& drysoil_frac, height, deadleaf_cover, emis, albedo, z0, roughheight, soiltype, rest_id, hist_id, hist2_id)
811
812	!
813	! call Soil Thermodynamic to write restart files
814	!
815	CALL thermosoil_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, indexgrnd, &
816	& temp_sol_new, snow, soilcap, soilflx, shumdiag, stempdiag, rest_id, hist_id, hist2_id)
817
818	!
819	! If we chose to route the water then we call the module. Else variables
820	! are set to zero.
821	!
822	!
823	IF ( river_routing .AND. nbp_glo .GT. 1) THEN
824	CALL routing_main (kjit, kjpindex, dtradia, ldrestart_read, ldrestart_write, index, &
825	& lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, runoff, &
826	& drainage, evapot_corr, precip_rain, humrel, &
827	& stempdiag, returnflow, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
828	ELSE
829	riverflow(:) = zero
830	coastalflow(:) = zero
831	returnflow(:) = zero
832	irrigation(:) = zero
833	ENDIF
834
835	!
836	! call slowproc_main to write restart files
837	!
838
839	CALL slowproc_main (kjit, kjpij, kjpindex, dtradia, date0, &
840	ldrestart_read, ldrestart_write, control%ok_co2, control%ok_stomate, &
841	index, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
842	t2mdiag, t2mdiag, temp_sol, stempdiag, &
843	vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
844	deadleaf_cover, &
845	assim_param, &
846	lai, height, veget, frac_nobio, veget_max, totfrac_nobio, qsintmax, &
847	rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC, &
848	co2_flux)
849
850
851	var_name= 'shumdiag'
852	CALL restput_p(rest_id, var_name, nbp_glo, nbdl, 1, kjit, shumdiag, 'scatter', nbp_glo, index_g)
853
854	var_name= 'runoff'
855	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, runoff, 'scatter', nbp_glo, index_g)
856	var_name= 'drainage'
857	CALL restput_p(rest_id, var_name, nbp_glo, 1, 1, kjit, drainage, 'scatter', nbp_glo, index_g)
858	END IF
859
860	IF (long_print) WRITE (numout,*) ' sechiba_main done '
861
862	END SUBROUTINE sechiba_main
863
864	!! Initialisation for SECHIBA processes
865	!! - does dynamic allocation for local arrays
866	!! - reads _restart_ file or set initial values to a raisonable value
867	!! - reads initial map
868	!!
869	SUBROUTINE sechiba_init (kjit, ldrestart_read, kjpij, kjpindex, index, rest_id, control_in, lalo)
870
871	! interface description
872	! input scalar
873	INTEGER(i_std), INTENT (in) :: kjit !! Time step number
874	LOGICAL,INTENT (in) :: ldrestart_read !! Logical for restart file to read
875	INTEGER(i_std), INTENT (in) :: kjpij !! Size of full domaine
876	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain size
877	INTEGER(i_std),DIMENSION (kjpindex), INTENT (in) :: index !! Indeces of the points on the map
878	INTEGER(i_std), INTENT (in) :: rest_id !! _Restart_ file identifier
879	TYPE(control_type), INTENT(in) :: control_in !! Flags that (de)activate parts of the model
880	! input fields
881	REAL(r_std),DIMENSION (kjpindex,2), INTENT (in) :: lalo !! Geographical coordinates
882	! output scalar
883	! output fields
884
885	! local declaration
886	INTEGER(i_std) :: ier, ji, jv
887	!
888	! initialisation
889	!
890	IF (l_first_sechiba) THEN
891	l_first_sechiba=.FALSE.
892	ELSE
893	WRITE (numout,*) ' l_first_sechiba false . we stop '
894	STOP 'sechiba_init'
895	ENDIF
896
897	! 1. make dynamic allocation with good dimension
898
899	! 1.0 The 3D vegetation indexation table
900
901	ALLOCATE (indexveg(kjpindex*nvm),stat=ier)
902	IF (ier.NE.0) THEN
903	WRITE (numout,) ' error in indexveg allocation. We stop. We need kjpindex words = ',kjpindexnvm
904	STOP 'sechiba_init'
905	END IF
906
907	ALLOCATE (indexsoil(kjpindex*nstm),stat=ier)
908	IF (ier.NE.0) THEN
909	WRITE (numout,) ' error in indexsoil allocation. We stop. We need kjpindex words = ',kjpindexnstm
910	STOP 'sechiba_init'
911	END IF
912
913	ALLOCATE (indexnobio(kjpindex*nnobio),stat=ier)
914	IF (ier.NE.0) THEN
915	WRITE (numout,) ' error in indexnobio allocation. We stop. We need kjpindex words = ',kjpindexnnobio
916	STOP 'sechiba_init'
917	END IF
918
919	ALLOCATE (indexgrnd(kjpindex*ngrnd),stat=ier)
920	IF (ier.NE.0) THEN
921	WRITE (numout,) ' error in indexgrnd allocation. We stop. We need kjpindex words = ',kjpindexngrnd
922	STOP 'sechiba_init'
923	END IF
924
925	ALLOCATE (indexlayer(kjpindex*nslm),stat=ier)
926	IF (ier.NE.0) THEN
927	WRITE (numout,) ' error in indexlayer allocation. We stop. We need kjpindex words = ',kjpindexnslm
928	STOP 'sechiba_init'
929	END IF
930
931	ALLOCATE (indexalb(kjpindex*2),stat=ier)
932	IF (ier.NE.0) THEN
933	WRITE (numout,) ' error in indexalb allocation. We stop. We need kjpindex words = ',kjpindex2
934	STOP 'sechiba_init'
935	END IF
936
937	! 1.1 one dimension array allocation with restartable value
938
939	IF (long_print) WRITE (numout,*) 'Allocation of 1D variables. We need for each kjpindex words = ',kjpindex
940	ALLOCATE (snow(kjpindex),stat=ier)
941	IF (ier.NE.0) THEN
942	WRITE (numout,*) ' error in snow allocation. We stop. We need kjpindex words = ',kjpindex
943	STOP 'sechiba_init'
944	END IF
945	snow(:) = undef_sechiba
946
947	ALLOCATE (snow_age(kjpindex),stat=ier)
948	IF (ier.NE.0) THEN
949	WRITE (numout,*) ' error in snow_age allocation. We stop. We need kjpindex words = ',kjpindex
950	STOP 'sechiba_init'
951	END IF
952	snow_age(:) = undef_sechiba
953
954	ALLOCATE (drysoil_frac(kjpindex),stat=ier)
955	IF (ier.NE.0) THEN
956	WRITE (numout,*) ' error in drysoil_frac allocation. We stop. We need kjpindex words = ',kjpindex
957	STOP 'sechiba_init'
958	END IF
959	drysoil_frac(:) = zero
960
961	ALLOCATE (rsol(kjpindex),stat=ier)
962	IF (ier.NE.0) THEN
963	WRITE (numout,*) ' error in rsol allocation. We stop. We need kjpindex words = ',kjpindex
964	STOP 'sechiba_init'
965	END IF
966
967	ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
968	IF (ier.NE.0) THEN
969	WRITE (numout,*) ' error in evap_bare_lim allocation. We stop. We need kjpindex words = ',kjpindex
970	STOP 'sechiba_init'
971	END IF
972
973	ALLOCATE (evapot(kjpindex),stat=ier)
974	IF (ier.NE.0) THEN
975	WRITE (numout,*) ' error in evapot allocation. We stop. We need kjpindex words = ',kjpindex
976	STOP 'sechiba_init'
977	END IF
978	evapot(:) = undef_sechiba
979
980	ALLOCATE (evapot_corr(kjpindex),stat=ier)
981	IF (ier.NE.0) THEN
982	WRITE (numout,*) ' error in evapot_corr allocation. We stop. We need kjpindex words = ',kjpindex
983	STOP 'sechiba_init'
984	END IF
985
986	ALLOCATE (humrel(kjpindex,nvm),stat=ier)
987	IF (ier.NE.0) THEN
988	WRITE (numout,*) ' error in humrel allocation. We stop. we need kjpindex words = ',kjpindex
989	STOP 'sechiba_init'
990	END IF
991	humrel(:,:) = undef_sechiba
992
993	ALLOCATE (vegstress(kjpindex,nvm),stat=ier)
994	IF (ier.NE.0) THEN
995	WRITE (numout,*) ' error in vegstress allocation. We stop. we need kjpindex words = ',kjpindex
996	STOP 'sechiba_init'
997	END IF
998	vegstress(:,:) = undef_sechiba
999
1000	ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
1001	IF (ier.NE.0) THEN
1002	WRITE (numout,*) ' error in soiltype allocation. We stop. we need kjpindex words = ',kjpindex
1003	STOP 'sechiba_init'
1004	END IF
1005	soiltype(:,:)=undef_sechiba
1006
1007	ALLOCATE (vbeta1(kjpindex),stat=ier)
1008	IF (ier.NE.0) THEN
1009	WRITE (numout,*) ' error in vbeta1 allocation. We stop. We need kjpindex words = ',kjpindex
1010	STOP 'sechiba_init'
1011	END IF
1012
1013	ALLOCATE (vbeta4(kjpindex),stat=ier)
1014	IF (ier.NE.0) THEN
1015	WRITE (numout,*) ' error in vbeta4 allocation. We stop. We need kjpindex words = ',kjpindex
1016	STOP 'sechiba_init'
1017	END IF
1018
1019	ALLOCATE (soilcap(kjpindex),stat=ier)
1020	IF (ier.NE.0) THEN
1021	WRITE (numout,*) ' error in soilcap allocation. We stop. We need kjpindex words = ',kjpindex
1022	STOP 'sechiba_init'
1023	END IF
1024
1025	ALLOCATE (soilflx(kjpindex),stat=ier)
1026	IF (ier.NE.0) THEN
1027	WRITE (numout,*) ' error in soilflx allocation. We stop. We need kjpindex words = ',kjpindex
1028	STOP 'sechiba_init'
1029	END IF
1030
1031	ALLOCATE (temp_sol(kjpindex),stat=ier)
1032	IF (ier.NE.0) THEN
1033	WRITE (numout,*) ' error in temp_sol allocation. We stop. We need kjpindex words = ',kjpindex
1034	STOP 'sechiba_init'
1035	END IF
1036	temp_sol(:) = undef_sechiba
1037
1038	ALLOCATE (qsurf(kjpindex),stat=ier)
1039	IF (ier.NE.0) THEN
1040	WRITE (numout,*) ' error in qsurf allocation. We stop. We need kjpindex words = ',kjpindex
1041	STOP 'sechiba_init'
1042	END IF
1043	qsurf(:) = undef_sechiba
1044
1045	! 1.2 two dimensions array allocation with restartable value
1046
1047	ALLOCATE (qsintveg(kjpindex,nvm),stat=ier)
1048	IF (ier.NE.0) THEN
1049	WRITE (numout,*) ' error in qsintveg allocation. We stop. We need kjpindex x nvm words = ',&
1050	& kjpindex,' x ' ,nvm,' = ',kjpindex*nvm
1051	STOP 'sechiba_init'
1052	END IF
1053	qsintveg(:,:) = undef_sechiba
1054
1055	ALLOCATE (vbeta2(kjpindex,nvm),stat=ier)
1056	IF (ier.NE.0) THEN
1057	WRITE (numout,*) ' error in vbeta2 allocation. We stop. We need kjpindex x nvm words = ',&
1058	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1059	STOP 'sechiba_init'
1060	END IF
1061
1062	ALLOCATE (vbeta3(kjpindex,nvm),stat=ier)
1063	IF (ier.NE.0) THEN
1064	WRITE (numout,*) ' error in vbeta3 allocation. We stop.We need kjpindex x nvm words = ',&
1065	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1066	STOP 'sechiba_init'
1067	END IF
1068
1069	ALLOCATE (vbetaco2(kjpindex,nvm),stat=ier)
1070	IF (ier.NE.0) THEN
1071	WRITE (numout,*) ' error in vbetaco2 allocation. We stop.We need kjpindex x nvm words = ',&
1072	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1073	STOP 'sechiba_init'
1074	END IF
1075
1076	ALLOCATE (cimean(kjpindex,nvm),stat=ier)
1077	IF (ier.NE.0) THEN
1078	WRITE (numout,*) ' error in cimean allocation. We stop.We need kjpindex x nvm words = ',&
1079	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1080	STOP 'sechiba_init'
1081	END IF
1082
1083	ALLOCATE (gpp(kjpindex,nvm),stat=ier)
1084	IF (ier.NE.0) THEN
1085	WRITE (numout,*) ' error in gpp allocation. We stop.We need kjpindex x nvm words = ',&
1086	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1087	STOP 'sechiba_init'
1088	END IF
1089	gpp(:,:) = undef_sechiba
1090
1091	ALLOCATE (veget(kjpindex,nvm),stat=ier)
1092	IF (ier.NE.0) THEN
1093	WRITE (numout,*) ' error in veget allocation. We stop. We need kjpindex x nvm words = ',&
1094	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1095	STOP 'sechiba_init'
1096	END IF
1097	veget(:,:)=undef_sechiba
1098
1099	ALLOCATE (veget_max(kjpindex,nvm),stat=ier)
1100	IF (ier.NE.0) THEN
1101	WRITE (numout,*) ' error in veget_max allocation. We stop. We need kjpindex x nvm words = ',&
1102	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1103	STOP 'sechiba_init'
1104	END IF
1105	veget_max(:,:)=undef_sechiba
1106
1107	ALLOCATE (lai(kjpindex,nvm),stat=ier)
1108	IF (ier.NE.0) THEN
1109	WRITE (numout,*) ' error in lai allocation. We stop. We need kjpindex x nvm words = ',&
1110	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1111	STOP 'sechiba_init'
1112	END IF
1113	lai(:,:)=undef_sechiba
1114
1115	ALLOCATE (height(kjpindex,nvm),stat=ier)
1116	IF (ier.NE.0) THEN
1117	WRITE (numout,*) ' error in height allocation. We stop. We need kjpindex x nvm words = ',&
1118	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1119	STOP 'sechiba_init'
1120	END IF
1121	height(:,:)=undef_sechiba
1122
1123	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
1124	IF (ier.NE.0) THEN
1125	WRITE (numout,*) ' error in frac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1126	STOP 'sechiba_init'
1127	END IF
1128	frac_nobio(:,:) = undef_sechiba
1129
1130	ALLOCATE (albedo(kjpindex,2),stat=ier)
1131	IF (ier.NE.0) THEN
1132	WRITE (numout,*) ' error in albedo allocation. We stop. We need kjpindex words = ',kjpindex
1133	STOP 'sechiba_init'
1134	END IF
1135
1136	ALLOCATE (snow_nobio(kjpindex,nnobio),stat=ier)
1137	IF (ier.NE.0) THEN
1138	WRITE (numout,*) ' error in snow_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1139	STOP 'sechiba_init'
1140	END IF
1141	snow_nobio(:,:) = undef_sechiba
1142
1143	ALLOCATE (snow_nobio_age(kjpindex,nnobio),stat=ier)
1144	IF (ier.NE.0) THEN
1145	WRITE (numout,*) ' error in snow_nobio_age allocation. We stop. We need kjpindex words = ',kjpindex
1146	STOP 'sechiba_init'
1147	END IF
1148	snow_nobio_age(:,:) = undef_sechiba
1149
1150	ALLOCATE (assim_param(kjpindex,nvm,npco2),stat=ier)
1151	IF (ier.NE.0) THEN
1152	WRITE (numout,*) ' error in assim_param allocation. We stop. We need kjpindex x nvm x npco2 words = ',&
1153	& kjpindex,' x ' ,nvm,' x ',npco2, ' = ',kjpindexnvmnpco2
1154	STOP 'sechiba_init'
1155	END IF
1156
1157	! 1.3 one dimension array allocation
1158
1159	ALLOCATE (vevapsno(kjpindex),stat=ier)
1160	IF (ier.NE.0) THEN
1161	WRITE (numout,*) ' error in vevapsno allocation. We stop. We need kjpindex words = ',kjpindex
1162	STOP 'sechiba_init'
1163	END IF
1164
1165	ALLOCATE (vevapnu(kjpindex),stat=ier)
1166	IF (ier.NE.0) THEN
1167	WRITE (numout,*) ' error in vevapnu allocation. We stop. We need kjpindex words = ',kjpindex
1168	STOP 'sechiba_init'
1169	END IF
1170
1171	ALLOCATE (t2mdiag(kjpindex),stat=ier)
1172	IF (ier.NE.0) THEN
1173	WRITE (numout,*) ' error in t2mdiag allocation. We stop. We need kjpindex words = ',kjpindex
1174	STOP 'sechiba_init'
1175	END IF
1176
1177	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
1178	IF (ier.NE.0) THEN
1179	WRITE (numout,*) ' error in totfrac_nobio allocation. We stop. We need kjpindex words = ',kjpindex
1180	STOP 'sechiba_init'
1181	END IF
1182
1183	ALLOCATE (runoff(kjpindex),stat=ier)
1184	IF (ier.NE.0) THEN
1185	WRITE (numout,*) ' error in runoff allocation. We stop. We need kjpindex words = ',kjpindex
1186	STOP 'sechiba_init'
1187	END IF
1188
1189	ALLOCATE (drainage(kjpindex),stat=ier)
1190	IF (ier.NE.0) THEN
1191	WRITE (numout,*) ' error in drainage allocation. We stop. We need kjpindex words = ',kjpindex
1192	STOP 'sechiba_init'
1193	END IF
1194
1195	ALLOCATE (returnflow(kjpindex),stat=ier)
1196	IF (ier.NE.0) THEN
1197	WRITE (numout,*) ' error in returnflow allocation. We stop. We need kjpindex words = ',kjpindex
1198	STOP 'sechiba_init'
1199	END IF
1200	returnflow(:) = zero
1201
1202	ALLOCATE (irrigation(kjpindex),stat=ier)
1203	IF (ier.NE.0) THEN
1204	WRITE (numout,*) ' error in irrigation allocation. We stop. We need kjpindex words = ',kjpindex
1205	STOP 'sechiba_init'
1206	END IF
1207	irrigation(:) = zero
1208
1209	ALLOCATE (z0(kjpindex),stat=ier)
1210	IF (ier.NE.0) THEN
1211	WRITE (numout,*) ' error in z0 allocation. We stop. We need kjpindex words = ',kjpindex
1212	STOP 'sechiba_init'
1213	END IF
1214
1215	ALLOCATE (roughheight(kjpindex),stat=ier)
1216	IF (ier.NE.0) THEN
1217	WRITE (numout,*) ' error in roughheight allocation. We stop. We need kjpindex words = ',kjpindex
1218	STOP 'sechiba_init'
1219	END IF
1220
1221	ALLOCATE (emis(kjpindex),stat=ier)
1222	IF (ier.NE.0) THEN
1223	WRITE (numout,*) ' error in emis allocation. We stop. We need kjpindex words = ',kjpindex
1224	STOP 'sechiba_init'
1225	END IF
1226
1227	ALLOCATE (tot_melt(kjpindex),stat=ier)
1228	IF (ier.NE.0) THEN
1229	WRITE (numout,*) ' error in tot_melt allocation. We stop. We need kjpindex words = ',kjpindex
1230	STOP 'sechiba_init'
1231	END IF
1232
1233	ALLOCATE (valpha(kjpindex),stat=ier)
1234	IF (ier.NE.0) THEN
1235	WRITE (numout,*) ' error in valpha allocation. We stop. We need kjpindex words = ',kjpindex
1236	STOP 'sechiba_init'
1237	END IF
1238
1239	ALLOCATE (vbeta(kjpindex),stat=ier)
1240	IF (ier.NE.0) THEN
1241	WRITE (numout,*) ' error in vbeta allocation. We stop. We need kjpindex words = ',kjpindex
1242	STOP 'sechiba_init'
1243	END IF
1244
1245	ALLOCATE (fusion(kjpindex),stat=ier)
1246	IF (ier.NE.0) THEN
1247	WRITE (numout,*) ' error in fusion allocation. We stop. We need kjpindex words = ',kjpindex
1248	STOP 'sechiba_init'
1249	END IF
1250
1251	ALLOCATE (rau(kjpindex),stat=ier)
1252	IF (ier.NE.0) THEN
1253	WRITE (numout,*) ' error in rau allocation. We stop. We need kjpindex words = ',kjpindex
1254	STOP 'sechiba_init'
1255	END IF
1256
1257	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
1258	IF (ier.NE.0) THEN
1259	WRITE (numout,*) ' error in deadleaf_cover allocation. We stop. We need kjpindex words = ',kjpindex
1260	STOP 'sechiba_init'
1261	END IF
1262
1263	ALLOCATE (stempdiag(kjpindex, nbdl),stat=ier)
1264	IF (ier.NE.0) THEN
1265	WRITE (numout,) ' error in stempdiag allocation. We stop. We need kjpindexnbdl words = ',&
1266	& kjpindex*nbdl
1267	STOP 'sechiba_init'
1268	END IF
1269	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
1270	IF (ier.NE.0) THEN
1271	WRITE (numout,) ' error in co2_flux allocation. We stop. We need kjpindexnvm words = ' ,kjpindex*nvm
1272	STOP 'sechiba_init'
1273	END IF
1274	co2_flux(:,:)=zero
1275
1276	ALLOCATE (shumdiag(kjpindex,nbdl),stat=ier)
1277	IF (ier.NE.0) THEN
1278	WRITE (numout,) ' error in shumdiag allocation. We stop. We need kjpindexnbdl words = ',&
1279	& kjpindex*nbdl
1280	STOP 'sechiba_init'
1281	END IF
1282
1283	ALLOCATE (litterhumdiag(kjpindex),stat=ier)
1284	IF (ier.NE.0) THEN
1285	WRITE (numout,*) ' error in litterhumdiag allocation. We stop. We need kjpindex words = ',kjpindex
1286	STOP 'sechiba_init'
1287	END IF
1288
1289	! 1.4 two dimensions array allocation
1290
1291	ALLOCATE (vevapwet(kjpindex,nvm),stat=ier)
1292	IF (ier.NE.0) THEN
1293	WRITE (numout,*) ' error in vevapwet allocation. We stop. We need kjpindex x nvm words = ',&
1294	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1295	STOP 'sechiba_init'
1296	END IF
1297
1298	ALLOCATE (transpir(kjpindex,nvm),stat=ier)
1299	IF (ier.NE.0) THEN
1300	WRITE (numout,*) ' error in transpir allocation. We stop. We need kjpindex x nvm words = ',&
1301	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1302	STOP 'sechiba_init'
1303	END IF
1304
1305	ALLOCATE (qsintmax(kjpindex,nvm),stat=ier)
1306	IF (ier.NE.0) THEN
1307	WRITE (numout,*) ' error in qsintmax allocation. We stop. We need kjpindex x nvm words = ',&
1308	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1309	STOP 'sechiba_init'
1310	END IF
1311
1312	ALLOCATE (rveget(kjpindex,nvm),stat=ier)
1313	IF (ier.NE.0) THEN
1314	WRITE (numout,*) ' error in rveget allocation. We stop. We need kjpindex x nvm words = ',&
1315	& kjpindex,' x ' ,nvm, ' = ',kjpindex*nvm
1316	STOP 'sechiba_init'
1317	END IF
1318
1319	!
1320	! 1.5 Get the indexing table for the vegetation fields. In SECHIBA we work on reduced grids but to store in the
1321	! full 3D filed vegetation variable we need another index table : indexveg, indexsoil, indexnobio and
1322	! indexgrnd
1323	!
1324	DO ji = 1, kjpindex
1325	!
1326	DO jv = 1, nvm
1327	indexveg((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1328	ENDDO
1329	!
1330	DO jv = 1, nstm
1331	indexsoil((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1332	ENDDO
1333	!
1334	DO jv = 1, nnobio
1335	indexnobio((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1336	ENDDO
1337	!
1338	DO jv = 1, ngrnd
1339	indexgrnd((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1340	ENDDO
1341	!
1342	DO jv = 1, nslm
1343	indexlayer((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1344	ENDDO
1345	!
1346	DO jv = 1, 2
1347	indexalb((jv-1)kjpindex + ji) = INDEX(ji) + (jv-1)kjpij
1348	ENDDO
1349	!
1350	ENDDO
1351
1352
1353	!
1354	! 2. restart value
1355	!
1356	! open restart input file if needed
1357	! and read data from restart input file
1358	!
1359
1360	IF (ldrestart_read) THEN
1361
1362	IF (long_print) WRITE (numout,*) ' we have to read a restart file for SECHIBA variables'
1363	!
1364	! Read the default value that will be put into variable which are not in the restart file
1365	!
1366	CALL ioget_expval(val_exp)
1367
1368	ENDIF
1369
1370	!
1371	river_routing = control_in%river_routing
1372	hydrol_cwrr = control_in%hydrol_cwrr
1373	!
1374	! 4. run control: store flags in a common variable
1375	!
1376
1377	control%ok_co2 = control_in%ok_co2
1378	control%ok_sechiba = control_in%ok_sechiba
1379	control%ok_stomate = control_in%ok_stomate
1380	control%ok_dgvm = control_in%ok_dgvm
1381	control%ok_pheno = control_in%ok_pheno
1382	control%stomate_watchout = control_in%stomate_watchout
1383
1384	IF (long_print) WRITE (numout,*) ' sechiba_init done '
1385
1386	END SUBROUTINE sechiba_init
1387	!
1388	!------------------------------------------------------------------
1389	!
1390	SUBROUTINE sechiba_clear (forcing_name,cforcing_name)
1391
1392	CHARACTER(LEN=100), INTENT(in) :: forcing_name
1393	CHARACTER(LEN=100), INTENT(in) :: cforcing_name
1394
1395	!
1396	! initialisation
1397	!
1398	l_first_sechiba=.TRUE.
1399
1400	! 1. Deallocate all dynamic variables
1401
1402	IF ( ALLOCATED (indexveg)) DEALLOCATE (indexveg)
1403	IF ( ALLOCATED (indexsoil)) DEALLOCATE (indexsoil)
1404	IF ( ALLOCATED (indexnobio)) DEALLOCATE (indexnobio)
1405	IF ( ALLOCATED (indexgrnd)) DEALLOCATE (indexgrnd)
1406	IF ( ALLOCATED (indexlayer)) DEALLOCATE (indexlayer)
1407	IF ( ALLOCATED (indexalb)) DEALLOCATE (indexalb)
1408	IF ( ALLOCATED (snow)) DEALLOCATE (snow)
1409	IF ( ALLOCATED (snow_age)) DEALLOCATE (snow_age)
1410	IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
1411	IF ( ALLOCATED (rsol)) DEALLOCATE (rsol)
1412	IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
1413	IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
1414	IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
1415	IF ( ALLOCATED (humrel)) DEALLOCATE (humrel)
1416	IF ( ALLOCATED (vegstress)) DEALLOCATE (vegstress)
1417	IF ( ALLOCATED (soiltype)) DEALLOCATE (soiltype)
1418	IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
1419	IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
1420	IF ( ALLOCATED (soilcap)) DEALLOCATE (soilcap)
1421	IF ( ALLOCATED (soilflx)) DEALLOCATE (soilflx)
1422	IF ( ALLOCATED (temp_sol)) DEALLOCATE (temp_sol)
1423	IF ( ALLOCATED (qsurf)) DEALLOCATE (qsurf)
1424	IF ( ALLOCATED (qsintveg)) DEALLOCATE (qsintveg)
1425	IF ( ALLOCATED (vbeta2)) DEALLOCATE (vbeta2)
1426	IF ( ALLOCATED (vbeta3)) DEALLOCATE (vbeta3)
1427	IF ( ALLOCATED (vbetaco2)) DEALLOCATE (vbetaco2)
1428	IF ( ALLOCATED (cimean)) DEALLOCATE (cimean)
1429	IF ( ALLOCATED (gpp)) DEALLOCATE (gpp)
1430	IF ( ALLOCATED (veget)) DEALLOCATE (veget)
1431	IF ( ALLOCATED (veget_max)) DEALLOCATE (veget_max)
1432	IF ( ALLOCATED (lai)) DEALLOCATE (lai)
1433	IF ( ALLOCATED (height)) DEALLOCATE (height)
1434	IF ( ALLOCATED (roughheight)) DEALLOCATE (roughheight)
1435	IF ( ALLOCATED (frac_nobio)) DEALLOCATE (frac_nobio)
1436	IF ( ALLOCATED (snow_nobio)) DEALLOCATE (snow_nobio)
1437	IF ( ALLOCATED (snow_nobio_age)) DEALLOCATE (snow_nobio_age)
1438	IF ( ALLOCATED (assim_param)) DEALLOCATE (assim_param)
1439	IF ( ALLOCATED (vevapsno)) DEALLOCATE (vevapsno)
1440	IF ( ALLOCATED (vevapnu)) DEALLOCATE (vevapnu)
1441	IF ( ALLOCATED (t2mdiag)) DEALLOCATE (t2mdiag)
1442	IF ( ALLOCATED (totfrac_nobio)) DEALLOCATE (totfrac_nobio)
1443	IF ( ALLOCATED (runoff)) DEALLOCATE (runoff)
1444	IF ( ALLOCATED (drainage)) DEALLOCATE (drainage)
1445	IF ( ALLOCATED (returnflow)) DEALLOCATE (returnflow)
1446	IF ( ALLOCATED (irrigation)) DEALLOCATE (irrigation)
1447	IF ( ALLOCATED (tot_melt)) DEALLOCATE (tot_melt)
1448	IF ( ALLOCATED (valpha)) DEALLOCATE (valpha)
1449	IF ( ALLOCATED (vbeta)) DEALLOCATE (vbeta)
1450	IF ( ALLOCATED (fusion)) DEALLOCATE (fusion)
1451	IF ( ALLOCATED (rau)) DEALLOCATE (rau)
1452	IF ( ALLOCATED (deadleaf_cover)) DEALLOCATE (deadleaf_cover)
1453	IF ( ALLOCATED (stempdiag)) DEALLOCATE (stempdiag)
1454	IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
1455	IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
1456	IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
1457	IF ( ALLOCATED (vevapwet)) DEALLOCATE (vevapwet)
1458	IF ( ALLOCATED (transpir)) DEALLOCATE (transpir)
1459	IF ( ALLOCATED (qsintmax)) DEALLOCATE (qsintmax)
1460	IF ( ALLOCATED (rveget)) DEALLOCATE (rveget)
1461	! 2. clear all modules
1462	CALL pft_clear
1463	CALL slowproc_clear
1464	CALL diffuco_clear
1465	CALL enerbil_clear
1466	IF ( hydrol_cwrr ) THEN
1467	CALL hydrol_clear
1468	ELSE
1469	CALL hydrolc_clear
1470	ENDIF
1471	CALL condveg_clear
1472	CALL thermosoil_clear
1473	CALL routing_clear
1474	!3. give name to next block
1475	stomate_forcing_name=forcing_name
1476	stomate_Cforcing_name=Cforcing_name
1477
1478	END SUBROUTINE sechiba_clear
1479
1480	!! SECHIBA's variables initialisation
1481	!! called every time step
1482	!!
1483	SUBROUTINE sechiba_var_init (kjpindex, rau, pb, temp_air)
1484
1485	! interface description
1486	! input scalar
1487	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain dimension
1488	! input fields
1489	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: pb !! Lowest level pressure
1490	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_air !! Air temperature
1491	! output fields
1492	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: rau !! Density
1493
1494	! local declaration
1495	INTEGER(i_std) :: ji
1496
1497	!
1498	! initialisation
1499	!
1500
1501	!
1502	! 1. calcul of rau: air density
1503	!
1504
1505	DO ji = 1,kjpindex
1506	rau(ji) = pa_par_hpa * pb(ji) / (cte_molr*temp_air(ji))
1507	END DO
1508
1509	IF (long_print) WRITE (numout,*) ' sechiba_var_init done '
1510
1511	END SUBROUTINE sechiba_var_init
1512
1513	!!
1514	!! Swap new fields to previous fields
1515	!!
1516	SUBROUTINE sechiba_end (kjpindex, dtradia, temp_sol, temp_sol_new)
1517
1518	! interface description
1519	! input scalar
1520	INTEGER(i_std), INTENT (in) :: kjpindex !! Domain dimension
1521	REAL(r_std),INTENT (in) :: dtradia !! Time step in seconds
1522	! input fields
1523	REAL(r_std),DIMENSION (kjpindex), INTENT (in) :: temp_sol_new !! New soil temperature
1524	! output fields
1525	REAL(r_std),DIMENSION (kjpindex), INTENT (out) :: temp_sol !! Soil temperature
1526
1527	!
1528	! swap
1529	!
1530	temp_sol(:) = temp_sol_new(:)
1531
1532	IF (long_print) WRITE (numout,*) ' sechiba_end done '
1533
1534	END SUBROUTINE sechiba_end
1535
1536	END MODULE sechiba

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

Download in other formats: