Context Navigation

source: branches/publications/ORCHIDEE_SOM_13C_r4859/src_driver/forcesoil.f90 @ 8076

Last change on this file since 8076 was 4023, checked in by marta.camino, 8 years ago
13C and delta 13C included without 13C discretization in soil
File size: 73.2 KB

Line
1	! =================================================================================================================================
2	! PROGRAM : forcesoil
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	!>\BRIEF This subroutine runs the soilcarbon submodel using specific initial conditions
9	!! and driving variables in order to obtain soil carbon stocks closed to the steady-state values
10	!! quicker than when using the ''full'' ORCHIDEE.
11	!!
12	!!\n DESCRIPTION: None
13	!! This subroutine computes the soil carbon stocks by calling the soilcarbon routine at each time step. \n
14	!! The aim is to obtain soil carbon stocks closed to the steady-state values and ultimately to create
15	!! an updated stomate restart file for the stomate component. The state variables of the subsystem are the clay content
16	!! (fixed value) and the soil carbon stocks. Initial conditions for the state variables are read in an
17	!! input stomate restart file. Driving variables are Soil carbon input, Water and Temperature stresses on
18	!! Organic Matter decomposition. Driving variables are read from a specific forcing file produced by a former run of ORCHIDEE
19	!! (SECHIBA+STOMATE). \n
20	!! The FORCESOIL program first consists in reading a set of input files, allocating variables and
21	!! preparing output stomate restart file. \n
22	!! Then, a loop over time is performed in which the soilcarbon routine is called at each time step. \n
23	!! Last, final values of the soil carbon stocks are written into the output stomate restart file. \n
24	!! No flag is associated with the use of the FORCESOIL program. \n
25	!!
26	!! RECENT CHANGE(S): None
27	!!
28	!! REFERENCE(S) : None
29	!!
30	!! FLOWCHART : None
31	!!
32	!! SVN :
33	!! $HeadURL: $
34	!! $Date: $
35	!! $Revision: $
36	!! \n
37	!_ =================================================================================================================================
38
39	PROGRAM forcesoil
40
41	USE netcdf
42	!-
43	USE defprec
44	USE constantes
45	USE constantes_soil
46	USE constantes_mtc
47	! USE pft_parameters
48	USE stomate_data
49	USE ioipsl_para
50	USE mod_orchidee_para
51	USE stomate_soilcarbon
52	USE vertical_soil
53	USE control
54	! USE constantes_var
55
56	!-
57	IMPLICIT NONE
58	!-
59	CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out
60	INTEGER(i_std) :: iim,jjm !! Indices (unitless)
61
62	INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless)
63	INTEGER(i_std) :: kjpindex !! Domain size (unitless)
64
65	INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?)
66	!! and number of time steps of the simulation (unitless)
67	CHARACTER(LEN=30) :: time_str !! Length of the simulation (year)
68	REAL(r_std) :: dt_files !! time step between two successive itaus (?)
69	!! (requested by restini routine) (seconds)
70	REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?)
71	INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless)
72	CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file
73	!! (unitless)
74	!-
75	CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless)
76	INTEGER :: Cforcing_id !! ID of the forcing file (unitless)
77	INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file
78	!! (unitless)
79	REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days)
80	INTEGER :: nparan !! Number of values stored per year in the forcing file
81	!! (unitless)
82	INTEGER :: nbyear
83	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless)
84	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless)
85	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless)
86	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined
87	!! in lat/lon (2D) (degrees)
88	REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless)
89
90
91	INTEGER :: i,k,j,l,m,iatt,iv,iyear !! counters (unitless)
92	REAL(r_std) :: soil_depth !! Soil depth (m)
93	CHARACTER(LEN=80) :: var_name
94	CHARACTER(LEN=20000) :: taboo_vars !! string used for storing the name of the variables
95	!! of the stomate restart file that are not automatically
96	!! duplicated from input to output restart file (unitless)
97	REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless)
98	INTEGER(i_std),PARAMETER :: nbvarmax=8000 !! maximum # of variables assumed in the stomate restart file
99	!! (unitless)
100	INTEGER(i_std) :: nbvar !! # of variables effectively present
101	!! in the stomate restart file (unitless)
102	CHARACTER(LEN=100),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored
103	!! in the stomate restart file (unitless)
104	INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable
105	!! of the stomate restart file
106	INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable
107	!! of the stomate restart file
108	INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable
109	!! of the stomate restart file
110	LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable
111	!! of the stomate restart file are of length 1 (ie scalar)
112	!! (unitless)
113	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless)
114	REAL(r_std) :: x_tmp !! temporary variable used to store return value
115	!! from nf90_get_att (unitless)
116	CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT
117	CHARACTER(LEN=1),DIMENSION(nelements) :: element_str !! string suffix indicating element
118
119	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: soil_ph_g !! Soil pH (0-14, pH unit)
120	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless)
121	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: bulk_dens_g !! soil bulk density (nbpglo) (g cm-3)
122	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:):: soiltile_g !! Fraction of each soil tile (0-1, unitless)
123	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:):: veget_max_g !! PFT "Maximal" coverage fraction of a PFT defined in
124	!! the input vegetation map
125	!! @tex $(m^2 m^{-2})$ @endtex, parallel computing
126	REAL(r_std),DIMENSION(:,:,:,:,:,:),ALLOCATABLE :: soilcarbon_input_g !! soil carbon input (nbpglob,npool,nvm,time)
127	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
128
129	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: control_temp_soil_g !! Temperature control (nbp_glo,nbdl,time) on OM decomposition
130	!! (unitless)
131	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: control_moist_soil_g !! Moisture control (nbp_glo,nbdl,time) on OM decomposition
132	!! ?? Should be defined per PFT as well (unitless)
133	REAL(r_std),DIMENSION (:,:,:), ALLOCATABLE :: moist_soil_g !! soil moisture content \f($m^3 \times m^3$)\f
134	REAL(r_std),DIMENSION (:,:,:,:), ALLOCATABLE :: soil_mc_g !! soil moisture content per soil type \f($m^3 \times m^3$)\f
135	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: floodout_g !! flux out of floodplains
136	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: wat_flux0_g !! Water flux in the first soil layers exported for soil C calculations
137	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: wat_flux_g !! Water flux in the soil layers exported for soil C calculations
138	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::runoff_per_soil_g !! Runoff per soil type [mm]
139	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::drainage_per_soil_g !! Drainage per soil type [mm]
140	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: DOC_to_topsoil_g !! DOC inputs to top of the soil column, from reinfiltration on
141	!! floodplains and from irrigation
142	!! @tex $(gC m^{-2} day{-1})$ @endtex
143	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: DOC_to_subsoil_g !! DOC inputs to bottom of the soil column, from returnflow
144	!! in swamps and lakes
145	!! @tex $(gC m^{-2} day{-1})$ @endtex
146	REAL(r_std),DIMENSION(:,:,:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$)
147	REAL(r_std), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: DOC_g !! Dissolved Organic Carbon in soil
148	!! The unit is given by m^3 of
149	!! ground @tex $(gC m{-2}) $ @endtex
150	REAL(r_std), ALLOCATABLE ::litter_above_g(:,:,:,:,:) !! Metabolic and structural litter, below ground
151	!! The unit is given by m^2 of
152	!! ground @tex $(gCi m{-2})$ @endtex
153	REAL(r_std), ALLOCATABLE ::litter_below_g(:,:,:,:,:,:) !! Metabolic and structural litter, below ground
154	!! The unit is given by m^2 of
155	!! ground @tex $(gC m^{-2})$ @endtex
156	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lignin_struc_above_g !! Ratio Lignine/Carbon in structural litter for above
157	!! ground compartments (unitless)
158	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: lignin_struc_below_g !! Ratio Lignine/Carbon in structural litter for below
159	!! ground compartments (unitless)
160	REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless)
161	REAL(r_std),ALLOCATABLE :: soil_ph(:) !! Soil pH (0-14 pH units)
162	REAL(r_std),ALLOCATABLE :: bulk_dens(:) !! soil bulk density (nbpglo) (g cm-3)
163	REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:,:,:) !! soil carbon input (nbp_loc,npool,nvm,time)
164	!! water @tex $(gC m{-2} of water)$ @endte
165	REAL(r_std),ALLOCATABLE :: control_temp_soil(:,:,:,:) !! Temperature control (nbp_loc,nbdl,time) on OM decomposition
166	!! (unitless)
167	REAL(r_std),ALLOCATABLE :: control_moist_soil(:,:,:,:) !! Moisture control (nbp_loc,nbdl,time) on OM decomposition
168	!! ?? Should be defined per PFT as well (unitless)
169
170	REAL(r_std),ALLOCATABLE :: carbon(:,:,:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$)
171	REAL(r_std), ALLOCATABLE ::litter_above(:,:,:,:,:) !! Metabolic and structural litter, below ground
172	!! The unit is given by m^2 of
173	!! ground @tex $(gCi m{-2})$ @endtex
174	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lignin_struc_above !! Ratio Lignine/Carbon in structural litter for above
175	!! ground compartments (unitless)
176	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: lignin_struc_below !! Ratio Lignine/Carbon in structural litter for below
177	!! ground compartments (unitless)
178	REAL(r_std),ALLOCATABLE :: veget_max(:,:) !! PFT "Maximal" coverage fraction of a PFT defined in
179	!! the input vegetation map
180	!! @tex $(m^2 m^{-2})$ @endtex
181	REAL(r_std), ALLOCATABLE ::litter_below(:,:,:,:,:,:) !! Metabolic and structural litter, below ground
182	!! The unit is given by m^2 of
183	!! ground @tex $(gC m^{-2})$ @endtex
184	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
185	!! (requested by soilcarbon routine but not used here)
186	REAL(r_std), DIMENSION(:,:), ALLOCATABLE :: soilhum !! Daily soil humidity of each soil layer
187	!! (unitless)
188	REAL(r_std),DIMENSION (:,:,:), ALLOCATABLE :: moist_soil !! soil moisture content \f($m^3 \times m^3$)\f
189	REAL(r_std),DIMENSION (:,:,:,:), ALLOCATABLE :: soil_mc !! soil moisture content \f($m^3 \times m^3$)\f
190	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: floodout !! flux out of floodplains
191	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: wat_flux0 !! Water flux in the first soil layers exported for soil C calculations
192	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:,:):: wat_flux !! Water flux in the soil layers exported for soil C calculations
193	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::runoff_per_soil !! Runoff per soil type [mm]
194	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::drainage_per_soil !! Drainage per soil type [mm]
195	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: DOC_to_topsoil !! DOC inputs to top of the soil column, from reinfiltration on
196	!! floodplains and from irrigation
197	!! @tex $(gC m^{-2} day{-1})$ @endtex
198	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: DOC_to_subsoil !! DOC inputs to bottom of the soil column, from returnflow
199	!! in swamps and lakes
200	!! @tex $(gC m^{-2} day{-1})$ @endtex
201	REAL(r_std),DIMENSION (:,:), ALLOCATABLE :: soiltile !! Fraction of each soil tile (0-1, unitless)
202	REAL(r_std), DIMENSION(:,:,:,:,:,:), ALLOCATABLE :: DOC !! Dissolved Organic Carbon in soil
203	!! The unit is given by m^3 of
204	!! ground @tex $(gC m{-2}) $ @endtex
205	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE :: DOC_EXP !! Dissolved Organic Carbon in soil
206	!! The unit is given by m^3 of
207	!! water @tex $(gCm {-3}) $ @endtex
208
209	INTEGER(i_std) :: ier,iret !! Used for hangling errors
210
211	CHARACTER(LEN=30) :: temp_name
212	LOGICAL :: debug !! boolean used for printing messages
213	LOGICAL :: l_error !! boolean for memory allocation
214
215	!_ =================================================================================================================================
216
217	CALL Init_orchidee_para
218	CALL init_timer
219
220	!-
221	debug = .FALSE.
222	CALL getin_p('DEBUG_INFO',debug)
223	!!-
224	!! 1. Initialisation stage
225	!! Reading a set of input files, allocating variables and preparing output restart file.
226	!!-
227	! Define restart file name
228	! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var).
229	! User values are used if present in the .def file.
230	! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used.
231	!-
232	IF (is_root_prc) THEN
233	sto_restname_in = 'stomate_start.nc'
234	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
235	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
236	sto_restname_out = 'stomate_rest_out.nc'
237	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
238	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
239	!-
240	! Open the input file and Get some Dimension and Attributes ID's
241	!-
242	iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto)
243	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g)
244	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g)
245	iret = NF90_INQ_VARID (rest_id_sto, "time", iv)
246	iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar)
247	iret = NF90_CLOSE (rest_id_sto)
248	i=INDEX(thecalendar,ACHAR(0))
249	IF ( i > 0 ) THEN
250	thecalendar(i:20)=' '
251	ENDIF
252	!-
253	! Allocate longitudes and latitudes
254	!-
255	ALLOCATE (lon(iim_g,jjm_g))
256	ALLOCATE (lat(iim_g,jjm_g))
257	lon(:,:) = zero
258	lat(:,:) = zero
259	lev(1) = zero
260	!-
261	CALL restini &
262	& (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
263	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
264	ENDIF
265	CALL control_initialize(dt_files)
266	CALL bcast(date0)
267	CALL bcast(thecalendar)
268	WRITE(numout,*) "calendar = ",thecalendar
269	!-
270	! calendar
271	!-
272	CALL ioconf_calendar (thecalendar)
273	CALL ioget_calendar (one_year,one_day)
274	CALL ioconf_startdate(date0)
275	!
276	!! For master process only
277	!
278	IF (is_root_prc) THEN
279	!-
280	! define forcing file's name (Cforcing_name var)
281	! User value is used if present in the .def file
282	! If not, default (NONE) is used
283	!-
284	Cforcing_name = 'NONE'
285	CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
286	!-
287	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
288	! and allocate variables.
289	!-
290	iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
291	IF (iret /= NF90_NOERR) THEN
292	CALL ipslerr (3,'forcesoil', &
293	& 'Could not open file : ', &
294	& Cforcing_name,'(Do you have forget it ?)')
295	ENDIF
296	!-
297	! Total Domain size is stored in nbp_glo variable
298	!-
299	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
300	nbp_glo = NINT(x_tmp)
301	!-
302	! Number of values stored per year in the forcing file is stored in nparan var.
303	!-
304	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
305	nparan = NINT(x_tmp)
306	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp)
307	nbyear = NINT(x_tmp)
308	!-
309	ALLOCATE (indices_g(nbp_glo))
310	ALLOCATE (clay_g(nbp_glo))
311	ALLOCATE (soil_ph_g(nbp_glo))
312	ALLOCATE (bulk_dens_g(nbp_glo))
313	ALLOCATE (soiltile_g(nbp_glo,nstm))
314	ALLOCATE (veget_max_g(nbp_glo,nvm))
315	!-
316	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
317	ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
318	ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g)
319	indices_g(:) = NINT(x_indices_g(:))
320	WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
321	DEALLOCATE (x_indices_g)
322	!-
323	ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
324	ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g)
325	ier = NF90_INQ_VARID (Cforcing_id,'soil_ph',v_id)
326	ier = NF90_GET_VAR (Cforcing_id,v_id,soil_ph_g)
327	ier = NF90_INQ_VARID (Cforcing_id,'bulk_dens',v_id)
328	ier = NF90_GET_VAR (Cforcing_id,v_id,bulk_dens_g)
329	ier = NF90_INQ_VARID (Cforcing_id,'soiltile',v_id)
330	ier = NF90_GET_VAR (Cforcing_id,v_id,soiltile_g)
331	ier = NF90_INQ_VARID (Cforcing_id,'veget_max',v_id)
332	ier = NF90_GET_VAR (Cforcing_id,v_id,veget_max_g)
333	!-
334	! time step of forcesoil program (in days)
335	!-
336	dt_forcesoil = one_year / FLOAT(nparan)
337	WRITE(numout,*) 'time step (d): ',dt_forcesoil
338	WRITE(numout,*) 'nparan: ',nparan
339	WRITE(numout,*) 'nbyear: ',nbyear
340	!-
341	! read and write the variables in the output restart file we do not modify within the Forcesoil program
342	! ie all variables stored in the input restart file except those stored in taboo_vars
343	!-
344	taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
345	& '$dt_days$ $date$'
346	! & '$day_counter$ $dt_days$ $date$ '
347	!-
348
349	DO l=1,nelements
350	IF (l == icarbon) THEN
351	element_str(l) = 'c'
352	ELSEIF (l == icarbon13) THEN
353	element_str(l) = '_c13'
354	ELSE
355	CALL ipslerr_p(3,'forcesoil','Define element_str','','')
356	ENDIF
357	ENDDO
358
359	DO k=1,nelements
360	DO m = 1,nvm
361	WRITE(part_str,'(I2)') m
362	IF (m < 10) part_str(1:1) = '0'
363	temp_name = '$carbon_z01_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
364	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
365	ENDDO
366
367	DO m = 1,nvm
368	WRITE(part_str,'(I2)') m
369	IF (m < 10) part_str(1:1) = '0'
370	temp_name = '$carbon_z02_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
371	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
372	ENDDO
373
374	DO m = 1,nvm
375	WRITE(part_str,'(I2)') m
376	IF (m < 10) part_str(1:1) = '0'
377	temp_name = '$carbon_z03_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
378	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
379	ENDDO
380
381	DO m = 1,nvm
382	WRITE(part_str,'(I2)') m
383	IF (m < 10) part_str(1:1) = '0'
384	temp_name = '$carbon_z04_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
385	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
386	ENDDO
387
388	DO m = 1,nvm
389	WRITE(part_str,'(I2)') m
390	IF (m < 10) part_str(1:1) = '0'
391	temp_name = '$carbon_z05_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
392	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
393	ENDDO
394
395	DO m = 1,nvm
396	WRITE(part_str,'(I2)') m
397	IF (m < 10) part_str(1:1) = '0'
398	temp_name = '$carbon_z06_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
399	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
400	ENDDO
401
402	DO m = 1,nvm
403	WRITE(part_str,'(I2)') m
404	IF (m < 10) part_str(1:1) = '0'
405	temp_name = '$carbon_z07_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
406	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
407	ENDDO
408
409	DO m = 1,nvm
410	WRITE(part_str,'(I2)') m
411	IF (m < 10) part_str(1:1) = '0'
412	temp_name = '$carbon_z08_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
413	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
414	ENDDO
415
416	DO m = 1,nvm
417	WRITE(part_str,'(I2)') m
418	IF (m < 10) part_str(1:1) = '0'
419	temp_name = '$carbon_z09_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
420	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
421	ENDDO
422
423	DO m = 1,nvm
424	WRITE(part_str,'(I2)') m
425	IF (m < 10) part_str(1:1) = '0'
426	temp_name = '$carbon_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
427	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
428	ENDDO
429
430	DO m = 1,nvm
431	WRITE(part_str,'(I2)') m
432	IF (m < 10) part_str(1:1) = '0'
433	temp_name = '$carbon_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
434	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
435	ENDDO
436
437	DO m = 1,npool
438	WRITE(part_str,'(I1)') m
439	temp_name = '$freedoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
440	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
441	ENDDO
442
443	DO m = 1,npool
444	WRITE(part_str,'(I1)') m
445	temp_name = '$freedoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
446	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
447	ENDDO
448
449	DO m = 1,npool
450	WRITE(part_str,'(I1)') m
451	temp_name = '$freedoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
452	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
453	ENDDO
454
455	DO m = 1,npool
456	WRITE(part_str,'(I1)') m
457	temp_name = '$freedoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
458	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
459	ENDDO
460
461	DO m = 1,npool
462	WRITE(part_str,'(I1)') m
463	temp_name = '$freedoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
464	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
465	ENDDO
466
467	DO m = 1,npool
468	WRITE(part_str,'(I1)') m
469	temp_name = '$freedoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
470	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
471	ENDDO
472
473	DO m = 1,npool
474	WRITE(part_str,'(I1)') m
475	temp_name = '$freedoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
476	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
477	ENDDO
478
479	DO m = 1,npool
480	WRITE(part_str,'(I1)') m
481	temp_name = '$freedoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
482	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
483	ENDDO
484
485	DO m = 1,npool
486	WRITE(part_str,'(I1)') m
487	temp_name = '$freedoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
488	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
489	ENDDO
490
491	DO m = 1,npool
492	WRITE(part_str,'(I1)') m
493	temp_name = '$freedoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
494	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
495	ENDDO
496
497	DO m = 1,npool
498	WRITE(part_str,'(I1)') m
499	temp_name = '$freedoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
500	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
501	ENDDO
502
503	DO m = 1,npool
504	WRITE(part_str,'(I1)') m
505	temp_name = '$adsdoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
506	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
507	ENDDO
508
509	DO m = 1,npool
510	WRITE(part_str,'(I1)') m
511	temp_name = '$adsdoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
512	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
513	ENDDO
514
515	DO m = 1,npool
516	WRITE(part_str,'(I1)') m
517	temp_name = '$adsdoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
518	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
519	ENDDO
520
521	DO m = 1,npool
522	WRITE(part_str,'(I1)') m
523	temp_name = '$adsdoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
524	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
525	ENDDO
526
527	DO m = 1,npool
528	WRITE(part_str,'(I1)') m
529	temp_name = '$adsdoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
530	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
531	ENDDO
532
533	DO m = 1,npool
534	WRITE(part_str,'(I1)') m
535	temp_name = '$adsdoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
536	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
537	ENDDO
538
539	DO m = 1,npool
540	WRITE(part_str,'(I1)') m
541	temp_name = '$adsdoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
542	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
543	ENDDO
544
545	DO m = 1,npool
546	WRITE(part_str,'(I1)') m
547	temp_name = '$adsdoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
548	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
549	ENDDO
550
551	DO m = 1,npool
552	WRITE(part_str,'(I1)') m
553	temp_name = '$adsdoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
554	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
555	ENDDO
556
557	DO m = 1,npool
558	WRITE(part_str,'(I1)') m
559	temp_name = '$adsdoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
560	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
561	ENDDO
562
563	DO m = 1,npool
564	WRITE(part_str,'(I1)') m
565	temp_name = '$adsdoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)//'$'
566	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
567	ENDDO
568	ENDDO
569
570	!-
571	CALL ioget_vname(rest_id_sto, nbvar, varnames)
572	!-
573	! read and write some special variables (1D or variables that we need)
574	!-
575	!var_name = 'day_counter'
576	!CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
577	!CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
578	!-
579	var_name = 'dt_days'
580	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
581	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
582	!-
583	var_name = 'date'
584	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
585	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
586	!-
587
588	DO iv=1,nbvar
589	!-- check if the variable is to be written here
590	IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
591	!---- get variable dimensions, especially 3rd dimension
592	CALL ioget_vdim &
593	& (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
594	l1d = ALL(vardims(1:varnbdim) == 1)
595	!---- read it
596	IF (l1d) THEN
597	CALL restget &
598	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
599	& 1, itau_dep, .TRUE., xtmp)
600	ELSE
601	ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
602	IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
603	!----
604
605	CALL restget &
606	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
607	& 1, itau_dep, .TRUE., var_3d, 'gather', nbp_glo, indices_g)
608	ENDIF
609	!---- write it
610	IF (l1d) THEN
611	CALL restput &
612	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
613	& 1, itau_dep, xtmp)
614	ELSE
615	CALL restput &
616	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
617	& 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g)
618	!----
619	DEALLOCATE(var_3d)
620	ENDIF
621	ENDIF
622	ENDDO
623	!-
624	! read soil carbon stocks values stored in the input restart file
625	!-
626	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm,nbdl,nelements),stat=ier)
627	IF (ier /= 0) THEN
628	WRITE(numout,*) 'Allocatoin error carbon_g ier = ',ier
629	STOP
630	END IF
631	carbon_g(:,:,:,:,:) = val_exp
632	ALLOCATE(DOC_g(nbp_glo,nvm,nbdl,ndoc,npool,nelements),stat=ier)
633	IF (ier /= 0) THEN
634	WRITE(numout,*) 'Allocatoin error DOC_g ier = ',ier
635	STOP
636	END IF
637
638	DOC_g(:,:,:,:,:,:) = val_exp
639	DO k=1,nelements
640	DO m=1,nvm
641	WRITE (part_str, '(I2)') m
642	IF (m<10) part_str(1:1)='0'
643	var_name = 'carbon_z01_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
644	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
645	& .TRUE., carbon_g(:,:,m,1,k), 'gather', nbp_glo, indices_g)
646	IF (ALL(carbon_g(:,:,m,1,k) == val_exp)) carbon_g(:,:,m,1,k) = zero
647	ENDDO
648
649	DO m=1,nvm
650	WRITE (part_str, '(I2)') m
651	IF (m<10) part_str(1:1)='0'
652	var_name = 'carbon_z02_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
653	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
654	& .TRUE., carbon_g(:,:,m,2,k), 'gather', nbp_glo, indices_g)
655	IF (ALL(carbon_g(:,:,m,2,k) == val_exp)) carbon_g(:,:,m,2,k) = zero
656	ENDDO
657
658	DO m=1,nvm
659	WRITE (part_str, '(I2)') m
660	IF (m<10) part_str(1:1)='0'
661	var_name = 'carbon_z03_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
662	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
663	& .TRUE., carbon_g(:,:,m,3,k), 'gather', nbp_glo, indices_g)
664	IF (ALL(carbon_g(:,:,m,3,k) == val_exp)) carbon_g(:,:,m,3,k) = zero
665	ENDDO
666
667	DO m=1,nvm
668	WRITE (part_str, '(I2)') m
669	IF (m<10) part_str(1:1)='0'
670	var_name = 'carbon_z04_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
671	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
672	& .TRUE., carbon_g(:,:,m,4,k), 'gather', nbp_glo, indices_g)
673	IF (ALL(carbon_g(:,:,m,4,k) == val_exp)) carbon_g(:,:,m,4,k) = zero
674	ENDDO
675
676	DO m=1,nvm
677	WRITE (part_str, '(I2)') m
678	IF (m<10) part_str(1:1)='0'
679	var_name = 'carbon_z05_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
680	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
681	& .TRUE., carbon_g(:,:,m,5,k), 'gather', nbp_glo, indices_g)
682	IF (ALL(carbon_g(:,:,m,5,k) == val_exp)) carbon_g(:,:,m,5,k) = zero
683	ENDDO
684
685	DO m=1,nvm
686	WRITE (part_str, '(I2)') m
687	IF (m<10) part_str(1:1)='0'
688	var_name = 'carbon_z06_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
689	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
690	& .TRUE., carbon_g(:,:,m,6,k), 'gather', nbp_glo, indices_g)
691	IF (ALL(carbon_g(:,:,m,6,k) == val_exp)) carbon_g(:,:,m,6,k) = zero
692	ENDDO
693
694	DO m=1,nvm
695	WRITE (part_str, '(I2)') m
696	IF (m<10) part_str(1:1)='0'
697	var_name = 'carbon_z07_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
698	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
699	& .TRUE., carbon_g(:,:,m,7,k), 'gather', nbp_glo, indices_g)
700	IF (ALL(carbon_g(:,:,m,7,k) == val_exp)) carbon_g(:,:,m,7,k) = zero
701	ENDDO
702
703	DO m=1,nvm
704	WRITE (part_str, '(I2)') m
705	IF (m<10) part_str(1:1)='0'
706	var_name = 'carbon_z08_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
707	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
708	& .TRUE., carbon_g(:,:,m,8,k), 'gather', nbp_glo, indices_g)
709	IF (ALL(carbon_g(:,:,m,8,k) == val_exp)) carbon_g(:,:,m,8,k) = zero
710	ENDDO
711
712	DO m=1,nvm
713	WRITE (part_str, '(I2)') m
714	IF (m<10) part_str(1:1)='0'
715	var_name = 'carbon_z09_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
716	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
717	& .TRUE., carbon_g(:,:,m,9,k), 'gather', nbp_glo, indices_g)
718	IF (ALL(carbon_g(:,:,m,9,k) == val_exp)) carbon_g(:,:,m,9,k) = zero
719	ENDDO
720
721	DO m=1,nvm
722	WRITE (part_str, '(I2)') m
723	IF (m<10) part_str(1:1)='0'
724	var_name = 'carbon_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
725	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
726	& .TRUE., carbon_g(:,:,m,10,k), 'gather', nbp_glo, indices_g)
727	IF (ALL(carbon_g(:,:,m,10,k) == val_exp)) carbon_g(:,:,m,10,k) = zero
728	ENDDO
729
730	DO m=1,nvm
731	WRITE (part_str, '(I2)') m
732	IF (m<10) part_str(1:1)='0'
733	var_name = 'carbon_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
734	CALL restget (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
735	& .TRUE., carbon_g(:,:,m,11,k), 'gather', nbp_glo, indices_g)
736	IF (ALL(carbon_g(:,:,m,11,k) == val_exp)) carbon_g(:,:,m,11,k) = zero
737	ENDDO
738
739	DO m=1,npool
740	WRITE (part_str, '(I1)') m
741	var_name = 'freedoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
742	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
743	& .TRUE., DOC_g(:,:,1,ifree,m,k), 'gather', nbp_glo, indices_g)
744	IF (ALL(DOC_g(:,:,1,ifree,m,k) == val_exp)) DOC_g(:,:,1,ifree,m,k) = zero
745	ENDDO
746	DO m=1,npool
747	WRITE (part_str, '(I1)') m
748	var_name = 'freedoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
749	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
750	& .TRUE., DOC_g(:,:,2,ifree,m,k), 'gather', nbp_glo, indices_g)
751	IF (ALL(DOC_g(:,:,2,ifree,m,k) == val_exp)) DOC_g(:,:,2,ifree,m,k) = zero
752	ENDDO
753
754	DO m=1,npool
755	WRITE (part_str, '(I1)') m
756	var_name = 'freedoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
757	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
758	& .TRUE., DOC_g(:,:,3,ifree,m,k), 'gather', nbp_glo, indices_g)
759	IF (ALL(DOC_g(:,:,3,ifree,m,k) == val_exp)) DOC_g(:,:,3,ifree,m,k) = zero
760	ENDDO
761
762	DO m=1,npool
763	WRITE (part_str, '(I1)') m
764	var_name = 'freedoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
765	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
766	& .TRUE., DOC_g(:,:,4,ifree,m,k), 'gather', nbp_glo, indices_g)
767	IF (ALL(DOC_g(:,:,4,ifree,m,k) == val_exp)) DOC_g(:,:,4,ifree,m,k) = zero
768	ENDDO
769
770	DO m=1,npool
771	WRITE (part_str, '(I1)') m
772	var_name = 'freedoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
773	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
774	& .TRUE., DOC_g(:,:,5,ifree,m,k), 'gather', nbp_glo, indices_g)
775	IF (ALL(DOC_g(:,:,5,ifree,m,k) == val_exp)) DOC_g(:,:,5,ifree,m,k) = zero
776	ENDDO
777
778	DO m=1,npool
779	WRITE (part_str, '(I1)') m
780	var_name = 'freedoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
781	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
782	& .TRUE., DOC_g(:,:,6,ifree,m,k), 'gather', nbp_glo, indices_g)
783	IF (ALL(DOC_g(:,:,6,ifree,m,k) == val_exp)) DOC_g(:,:,6,ifree,m,k) = zero
784	ENDDO
785
786	DO m=1,npool
787	WRITE (part_str, '(I1)') m
788	var_name = 'freedoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
789	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
790	& .TRUE., DOC_g(:,:,7,ifree,m,k), 'gather', nbp_glo, indices_g)
791	IF (ALL(DOC_g(:,:,7,ifree,m,k) == val_exp)) DOC_g(:,:,7,ifree,m,k) = zero
792	ENDDO
793
794	DO m=1,npool
795	WRITE (part_str, '(I1)') m
796	var_name = 'freedoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
797	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
798	& .TRUE., DOC_g(:,:,8,ifree,m,k), 'gather', nbp_glo, indices_g)
799	IF (ALL(DOC_g(:,:,8,ifree,m,k) == val_exp)) DOC_g(:,:,8,ifree,m,k) = zero
800	ENDDO
801
802	DO m=1,npool
803	WRITE (part_str, '(I1)') m
804	var_name = 'freedoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
805	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
806	& .TRUE., DOC_g(:,:,9,ifree,m,k), 'gather', nbp_glo, indices_g)
807	IF (ALL(DOC_g(:,:,9,ifree,m,k) == val_exp)) DOC_g(:,:,9,ifree,m,k) = zero
808	ENDDO
809
810	DO m=1,npool
811	WRITE (part_str, '(I1)') m
812	var_name = 'freedoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
813	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
814	& .TRUE., DOC_g(:,:,10,ifree,m,k), 'gather', nbp_glo, indices_g)
815	IF (ALL(DOC_g(:,:,10,ifree,m,k) == val_exp)) DOC_g(:,:,10,ifree,m,k) = zero
816	ENDDO
817
818	DO m=1,npool
819	WRITE (part_str, '(I1)') m
820	var_name = 'freedoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
821	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
822	& .TRUE., DOC_g(:,:,11,ifree,m,k), 'gather', nbp_glo, indices_g)
823	IF (ALL(DOC_g(:,:,11,ifree,m,k) == val_exp)) DOC_g(:,:,11,ifree,m,k) = zero
824	ENDDO
825	DO m=1,npool
826	WRITE (part_str, '(I1)') m
827	var_name = 'adsdoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
828	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
829	& .TRUE., DOC_g(:,:,1,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
830	IF (ALL(DOC_g(:,:,1,iadsorbed,m,k) == val_exp)) DOC_g(:,:,1,iadsorbed,m,k) = zero
831	ENDDO
832	DO m=1,npool
833	WRITE (part_str, '(I1)') m
834	var_name = 'adsdoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
835	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
836	& .TRUE., DOC_g(:,:,2,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
837	IF (ALL(DOC_g(:,:,2,iadsorbed,m,k) == val_exp)) DOC_g(:,:,2,iadsorbed,m,k) = zero
838	ENDDO
839
840	DO m=1,npool
841	WRITE (part_str, '(I1)') m
842	var_name = 'adsdoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
843	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
844	& .TRUE., DOC_g(:,:,3,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
845	IF (ALL(DOC_g(:,:,3,iadsorbed,m,k) == val_exp)) DOC_g(:,:,3,iadsorbed,m,k) = zero
846	ENDDO
847
848	DO m=1,npool
849	WRITE (part_str, '(I1)') m
850	var_name = 'adsdoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
851	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
852	& .TRUE., DOC_g(:,:,4,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
853	IF (ALL(DOC_g(:,:,4,iadsorbed,m,k) == val_exp)) DOC_g(:,:,4,iadsorbed,m,k) = zero
854	ENDDO
855
856	DO m=1,npool
857	WRITE (part_str, '(I1)') m
858	var_name = 'adsdoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
859	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
860	& .TRUE., DOC_g(:,:,5,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
861	IF (ALL(DOC_g(:,:,5,iadsorbed,m,k) == val_exp)) DOC_g(:,:,5,iadsorbed,m,k) = zero
862	ENDDO
863
864	DO m=1,npool
865	WRITE (part_str, '(I1)') m
866	var_name = 'adsdoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
867	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
868	& .TRUE., DOC_g(:,:,6,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
869	IF (ALL(DOC_g(:,:,6,iadsorbed,m,k) == val_exp)) DOC_g(:,:,6,iadsorbed,m,k) = zero
870	ENDDO
871
872	DO m=1,npool
873	WRITE (part_str, '(I1)') m
874	var_name = 'adsdoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
875	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
876	& .TRUE., DOC_g(:,:,7,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
877	IF (ALL(DOC_g(:,:,7,iadsorbed,m,k) == val_exp)) DOC_g(:,:,7,iadsorbed,m,k) = zero
878	ENDDO
879
880	DO m=1,npool
881	WRITE (part_str, '(I1)') m
882	var_name = 'adsdoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
883	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
884	& .TRUE., DOC_g(:,:,8,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
885	IF (ALL(DOC_g(:,:,8,iadsorbed,m,k) == val_exp)) DOC_g(:,:,8,iadsorbed,m,k) = zero
886	ENDDO
887
888	DO m=1,npool
889	WRITE (part_str, '(I1)') m
890	var_name = 'adsdoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
891	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
892	& .TRUE., DOC_g(:,:,9,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
893	IF (ALL(DOC_g(:,:,9,iadsorbed,m,k) == val_exp)) DOC_g(:,:,9,iadsorbed,m,k) = zero
894	ENDDO
895
896	DO m=1,npool
897	WRITE (part_str, '(I1)') m
898	var_name = 'adsdoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
899	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
900	& .TRUE., DOC_g(:,:,10,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
901	IF (ALL(DOC_g(:,:,10,iadsorbed,m,k) == val_exp)) DOC_g(:,:,10,iadsorbed,m,k) = zero
902	ENDDO
903
904	DO m=1,npool
905	WRITE (part_str, '(I1)') m
906	var_name = 'adsdoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
907	CALL restget (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
908	& .TRUE., DOC_g(:,:,11,iadsorbed,m,k), 'gather', nbp_glo, indices_g)
909	IF (ALL(DOC_g(:,:,11,iadsorbed,m,k) == val_exp)) DOC_g(:,:,11,iadsorbed,m,k) = zero
910	ENDDO
911	ENDDO
912
913	WRITE(numout,*) "date0 : ",date0, itau_dep
914	!-
915	! Analytical spinup is set to false
916	!
917	spinup_analytic = .FALSE.
918
919	! Length of the run (in Years)
920	! User value is used if present in the .def file
921	! If not, default value (10000 Years) is used
922	!-
923	WRITE(time_str,'(a)') '10000Y'
924	CALL getin('TIME_LENGTH', time_str)
925	write(numout,*) 'Number of years for carbon spinup : ',time_str
926	! transform into itau
927	CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
928	write(numout,*) 'Number of time steps to do: ',itau_len
929	!-
930	! read soil carbon inputs, water and temperature stresses on OM decomposition
931	! into the forcing file - We read an average year.
932	!-
933	ALLOCATE(soilcarbon_input_g(nbp_glo,nvm,nbdl,npool,nelements,nparan*nbyear))
934	ALLOCATE(litter_above_g(nbp_glo,nlitt,nvm,nelements,nparan*nbyear))
935	ALLOCATE(litter_below_g(nbp_glo,nlitt,nvm,nbdl,nelements,nparan*nbyear))
936	ALLOCATE(lignin_struc_above_g(nbp_glo,nvm,nparan*nbyear))
937	ALLOCATE(lignin_struc_below_g(nbp_glo,nvm,nbdl,nparan*nbyear))
938	!-
939	ALLOCATE(control_temp_soil_g(nbp_glo,nbdl,npool2,nparannbyear))
940	ALLOCATE(control_moist_soil_g(nbp_glo,nbdl,nvm,nparan*nbyear))
941	ALLOCATE(moist_soil_g(nbp_glo,nbdl,nparan*nbyear))
942	ALLOCATE(soil_mc_g(nbp_glo,nbdl,nstm,nparan*nbyear))
943	ALLOCATE(floodout_g(nbp_glo,nparan*nbyear))
944	ALLOCATE(wat_flux0_g(nbp_glo,nstm,nparan*nbyear))
945	ALLOCATE(wat_flux_g(nbp_glo,nbdl,nstm,nparan*nbyear))
946	ALLOCATE(runoff_per_soil_g(nbp_glo,nstm,nparan*nbyear))
947	ALLOCATE(drainage_per_soil_g(nbp_glo,nstm,nparan*nbyear))
948	ALLOCATE(DOC_to_topsoil_g(nbp_glo,nelements,nparan*nbyear))
949	ALLOCATE(DOC_to_subsoil_g(nbp_glo,nelements,nparan*nbyear))
950	!-
951	ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
952	ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g)
953	ier = NF90_INQ_VARID (Cforcing_id, 'control_moist_soil',v_id)
954	ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_soil_g)
955	ier = NF90_INQ_VARID (Cforcing_id, 'control_temp_soil',v_id)
956	ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_soil_g)
957	ier = NF90_INQ_VARID (Cforcing_id, 'litter_above',v_id)
958	ier = NF90_GET_VAR (Cforcing_id,v_id,litter_above_g)
959	ier = NF90_INQ_VARID (Cforcing_id, 'litter_below',v_id)
960	ier = NF90_GET_VAR (Cforcing_id,v_id,litter_below_g)
961	ier = NF90_INQ_VARID (Cforcing_id, 'moist_soil',v_id)
962	ier = NF90_GET_VAR (Cforcing_id,v_id,moist_soil_g)
963	ier = NF90_INQ_VARID (Cforcing_id, 'soil_mc',v_id)
964	ier = NF90_GET_VAR (Cforcing_id,v_id,soil_mc_g)
965	ier = NF90_INQ_VARID (Cforcing_id, 'floodout',v_id)
966	ier = NF90_GET_VAR (Cforcing_id,v_id,floodout_g)
967	ier = NF90_INQ_VARID (Cforcing_id, 'wat_flux0',v_id)
968	ier = NF90_GET_VAR (Cforcing_id,v_id,wat_flux0_g)
969	ier = NF90_INQ_VARID (Cforcing_id, 'wat_flux',v_id)
970	ier = NF90_GET_VAR (Cforcing_id,v_id,wat_flux_g)
971	ier = NF90_INQ_VARID (Cforcing_id, 'runoff_per_soil',v_id)
972	ier = NF90_GET_VAR (Cforcing_id,v_id,runoff_per_soil_g)
973	ier = NF90_INQ_VARID (Cforcing_id, 'drainage_per_soil',v_id)
974	ier = NF90_GET_VAR (Cforcing_id,v_id,drainage_per_soil_g)
975	ier = NF90_INQ_VARID (Cforcing_id, 'DOC_to_topsoil',v_id)
976	ier = NF90_GET_VAR (Cforcing_id,v_id,DOC_to_topsoil_g)
977	ier = NF90_INQ_VARID (Cforcing_id, 'DOC_to_subsoil',v_id)
978	ier = NF90_GET_VAR (Cforcing_id,v_id,DOC_to_subsoil_g)
979	ier = NF90_INQ_VARID (Cforcing_id, 'lignin_struc_above',v_id)
980	ier = NF90_GET_VAR (Cforcing_id,v_id,lignin_struc_above_g)
981	ier = NF90_INQ_VARID (Cforcing_id, 'lignin_struc_below',v_id)
982	ier = NF90_GET_VAR (Cforcing_id,v_id,lignin_struc_below_g)
983	!-
984	ier = NF90_CLOSE (Cforcing_id)
985	!-
986	ENDIF
987	CALL bcast(nparan)
988	CALL bcast(nbyear)
989	CALL bcast(dt_forcesoil)
990	CALL bcast(iim_g)
991	CALL bcast(jjm_g)
992	CALL bcast(nbp_glo)
993	CALL bcast(itau_dep)
994	CALL bcast(itau_len)
995	IF (.NOT. ALLOCATED(indices_g)) ALLOCATE (indices_g(nbp_glo))
996	CALL bcast(indices_g)
997
998	!
999	! We must initialize data_para :
1000	CALL init_orchidee_data_para_driver(nbp_glo,indices_g)
1001
1002	kjpindex=nbp_loc
1003	jjm=jj_nb
1004	iim=iim_g
1005	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
1006
1007	!---
1008	!--- Create the index table
1009	!---
1010	!--- This job returns a LOCAL kindex.
1011	!---
1012	ALLOCATE (indices(kjpindex),stat=ier)
1013	!
1014	!! scattering to all processes in parallel mode
1015	!
1016	CALL scatter(indices_g,indices)
1017	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
1018	IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
1019	!-
1020	! Allocation of the variables for a processor
1021	!-
1022	ALLOCATE(clay(kjpindex))
1023	ALLOCATE(soil_ph(kjpindex))
1024	ALLOCATE(bulk_dens(kjpindex))
1025	ALLOCATE(soiltile(kjpindex,nstm))
1026	ALLOCATE(veget_max(kjpindex,nvm))
1027	ALLOCATE(soilcarbon_input(kjpindex,nvm,nbdl,npool,nelements,nparan*nbyear))
1028	ALLOCATE(control_temp_soil(kjpindex,nbdl,npool2,nparannbyear))
1029	ALLOCATE(control_moist_soil(kjpindex,nbdl,nvm,nparan*nbyear))
1030	ALLOCATE(carbon(kjpindex,ncarb,nvm,nbdl,nelements))
1031	ALLOCATE(litter_above(kjpindex,nlitt,nvm,nelements,nparan*nbyear))
1032	ALLOCATE(litter_below(kjpindex,nlitt,nvm,nbdl,nelements,nparan*nbyear))
1033	ALLOCATE(resp_hetero_soil(kjpindex,nvm,nelements))
1034	ALLOCATE(soilhum(kjpindex,nbdl))
1035	ALLOCATE(moist_soil(kjpindex,nbdl,nparan*nbyear))
1036	ALLOCATE(soil_mc(kjpindex,nbdl,nstm,nparan*nbyear))
1037	ALLOCATE(floodout(kjpindex,nparan*nbyear))
1038	ALLOCATE(wat_flux0(kjpindex,nstm,nparan*nbyear))
1039	ALLOCATE(wat_flux(kjpindex,nbdl,nstm,nparan*nbyear))
1040	ALLOCATE(runoff_per_soil(kjpindex,nstm,nparan*nbyear))
1041	ALLOCATE(drainage_per_soil(kjpindex,nstm,nparan*nbyear))
1042	ALLOCATE(DOC_to_topsoil(kjpindex,nelements,nparan*nbyear))
1043	ALLOCATE(DOC_to_subsoil(kjpindex,nelements,nparan*nbyear))
1044	ALLOCATE(DOC(kjpindex,nvm,nbdl,ndoc,npool,nelements))
1045	ALLOCATE(DOC_EXP(kjpindex,nvm,nexp,nelements))
1046	ALLOCATE(lignin_struc_above(kjpindex,nvm,nparan*nbyear))
1047	ALLOCATE(lignin_struc_below(kjpindex,nvm,nbdl,nparan*nbyear))
1048	iatt = 0
1049	!-
1050	! Initialization of the variables for a processor
1051	!-
1052	CALL Scatter(clay_g,clay)
1053	CALL Scatter(soil_ph_g,soil_ph)
1054	CALL Scatter(bulk_dens_g,bulk_dens)
1055	CALL Scatter(soiltile_g,soiltile)
1056	CALL Scatter(veget_max_g,veget_max)
1057	DO m =1,nvm
1058	DO l=1,nbdl
1059	DO i=1,npool
1060	CALL Scatter(soilcarbon_input_g(:,m,l,i,:,:),soilcarbon_input(:,m,l,i,:,:))
1061	ENDDO
1062	ENDDO
1063	ENDDO
1064	DO i =1,nlitt
1065	DO j=1,nvm
1066	CALL Scatter(litter_above_g(:,i,j,:,:),litter_above(:,i,j,:,:))
1067	ENDDO
1068	ENDDO
1069
1070	DO i =1,nlitt
1071	DO j=1,nvm
1072	DO k = 1,nbdl
1073	CALL Scatter(litter_below_g(:,i,j,k,:,:),litter_below(:,i,j,k,:,:))
1074	ENDDO
1075	ENDDO
1076	ENDDO
1077	CALL Scatter(control_temp_soil_g,control_temp_soil)
1078	CALL Scatter(control_moist_soil_g,control_moist_soil)
1079	CALL Scatter(moist_soil_g,moist_soil)
1080	CALL Scatter(soil_mc_g,soil_mc)
1081	CALL Scatter(floodout_g,floodout)
1082	CALL Scatter(wat_flux0_g,wat_flux0)
1083	CALL Scatter(wat_flux_g,wat_flux)
1084	CALL Scatter(runoff_per_soil_g,runoff_per_soil)
1085	CALL Scatter(drainage_per_soil_g,drainage_per_soil)
1086	CALL Scatter(DOC_to_topsoil_g,DOC_to_topsoil)
1087	CALL Scatter(DOC_to_subsoil_g,DOC_to_subsoil)
1088	DO m =1,nvm
1089	DO l=1,nbdl
1090	DO k=1,nelements
1091	CALL Scatter(carbon_g(:,:,m,l,k),carbon(:,:,m,l,k))
1092	ENDDO
1093	ENDDO
1094	ENDDO
1095	DO m =1,nbdl
1096	DO l=1,ndoc
1097	DO i= 1,npool
1098	DO k=1,nelements
1099	CALL Scatter(DOC_g(:,:,m,l,i,k),DOC(:,:,m,l,i,k))
1100	ENDDO
1101	ENDDO
1102	ENDDO
1103	ENDDO
1104	DO j=1,nvm
1105	CALL Scatter(lignin_struc_above_g(:,j,:),lignin_struc_above(:,j,:))
1106	ENDDO
1107	DO j=1,nvm
1108	DO k = 1,nbdl
1109	CALL Scatter(lignin_struc_below_g(:,j,k,:),lignin_struc_below(:,j,k,:))
1110	ENDDO
1111	ENDDO
1112	!-
1113	! Configu
1114	!-
1115
1116	!Config Key = FRAC_CARB_AP
1117	!Config Desc = frac carb coefficients from active pool: depends on clay content
1118	!Config if = OK_STOMATE
1119	!Config Def = 0.004
1120	!Config Help = fraction of the active pool going to the passive pool
1121	!Config Units = [-]
1122	CALL getin_p('FRAC_CARB_AP',frac_carb_ap)
1123	!
1124	!Config Key = FRAC_CARB_SA
1125	!Config Desc = frac_carb_coefficients from slow pool
1126	!Config if = OK_STOMATE
1127	!Config Def = 0.42
1128	!Config Help = fraction of the slow pool going to the active pool
1129	!Config Units = [-]
1130	CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
1131	!
1132	!Config Key = FRAC_CARB_SP
1133	!Config Desc = frac_carb_coefficients from slow pool
1134	!Config if = OK_STOMATE
1135	!Config Def = 0.03
1136	!Config Help = fraction of the slow pool going to the passive pool
1137	!Config Units = [-]
1138	!CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
1139	!
1140	!Config Key = FRAC_CARB_PA
1141	!Config Desc = frac_carb_coefficients from passive pool
1142	!Config if = OK_STOMATE
1143	!Config Def = 0.45
1144	!Config Help = fraction of the passive pool going to the passive pool
1145	!Config Units = [-]
1146	CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
1147	!
1148	!Config Key = FRAC_CARB_PS
1149	!Config Desc = frac_carb_coefficients from passive pool
1150	!Config if = OK_STOMATE
1151	!Config Def = 0.0
1152	!Config Help = fraction of the passive pool going to the passive pool
1153	!Config Units = [-]
1154	!CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
1155	!
1156	!Config Key = ACTIVE_TO_PASS_CLAY_FRAC
1157	!Config Desc =
1158	!Config if = OK_STOMATE
1159	!Config Def = .68
1160	!Config Help =
1161	!Config Units = [-]
1162	CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
1163	!
1164	!Config Key = CARBON_TAU_IACTIVE
1165	!Config Desc = residence times in carbon pools
1166	!Config if = OK_STOMATE
1167	!Config Def = 0.3
1168	!Config Help =
1169	!Config Units = [days]
1170	CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
1171	!
1172	!Config Key = CARBON_TAU_ISLOW
1173	!Config Desc = residence times in carbon pools
1174	!Config if = OK_STOMATE
1175	!Config Def = 1.12
1176	!Config Help =
1177	!Config Units = [days]
1178	CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
1179	!
1180	!Config Key = CARBON_TAU_IPASSIVE
1181	!Config Desc = residence times in carbon pools
1182	!Config if = OK_STOMATE
1183	!Config Def = 461.98
1184	!Config Help = residence time in the passive pool
1185	!Config Units = [days]
1186	CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
1187	!
1188	!Config Key = PRIMING_PARAM_IACTIVE
1189	!Config Desc = priming parameter in carbon pools
1190	!Config if = OK_STOMATE
1191	!Config Def = 493.66
1192	!Config Help =
1193	!Config Units = [-]
1194	CALL getin_p('PRIMING_PARAM_IACTIVE',priming_param_iactive)
1195	!
1196	!Config Key = PRIMING_PARAM_ISLOW
1197	!Config Desc = priming parameter in carbon pools
1198	!Config if = OK_STOMATE
1199	!Config Def = 194.03
1200	!Config Help =
1201	!Config Units = [-]
1202	CALL getin_p('PRIMING_PARAM_ISLOW',priming_param_islow)
1203	!
1204	!Config Key = PRIMING_PARAM_IPASSIVE
1205	!Config Desc = priming parameter in carbon pools
1206	!Config if = OK_STOMATE
1207	!Config Def = 136.54
1208	!Config Help =
1209	!Config Units = [-]
1210	CALL getin_p('PRIMING_PARAM_IPASSIVE',priming_param_ipassive)
1211	!
1212	! !Config Key = BULK_DENSITY
1213	! !Config Desc = value of the soil bulk density
1214	! !Config if = READ_BULK_DENSITY
1215	! !Config Def = 1.65
1216	! !Config Help =
1217	! !Config Units = [kg m-3]
1218	! CALL getin_p('BULK_DENSITY',bulk_dens)
1219	!
1220	!Config Key = FLUX_TOT_COEFF
1221	!Config Desc =
1222	!Config if = OK_STOMATE
1223	!Config Def = 1.2, 1.4,.75
1224	!Config Help =
1225	!Config Units = [days]
1226	CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)
1227	!
1228	!Config Key = PREF_SOIL_VEG
1229	!Config Desc = The soil tile number for each vegetation
1230	!Config if = OK_SECHIBA or OK_STOMATE
1231	!Config Def = 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3
1232	!Config Help = Gives the number of the soil tile on which we will
1233	!Config put each vegetation. This allows to divide the
1234	!hydrological column
1235	!Config Units = [-]
1236	CALL getin_p('PREF_SOIL_VEG',pref_soil_veg)
1237	!
1238	!Config Key = DOC_tau_labile
1239	!Config Desc = residence time of labile DOC
1240	!Config if = OK_STOMATE
1241	!Config Def = 3.76
1242	!Config Help = mean from Kalbitz et al., 2003 Geoderma
1243	!Config Units = [days]
1244	CALL getin_p('DOC_TAU_LABILE',DOC_tau_labile)
1245	!
1246	!Config Key = DOC_tau_stable
1247	!Config Desc = residence time of stable DOC
1248	!Config if = OK_STOMATE
1249	!Config Def = 602.41
1250	!Config Help = mean from Kalbitz et al., 2003 Geoderma
1251	!Config Units = [days]
1252	CALL getin_p('DOC_TAU_STABLE',DOC_tau_stable)
1253	!
1254	!!Config Key = D_DOC
1255	!Config Desc = diffusion coefficient of DOC
1256	!Config if = OK_STOMATE
1257	!Config Def = 0.0000004428
1258	!Config Help = from Burdige et al., 1999 in Ota et al., 2013
1259	!Config Units = [m^2 hr-1]
1260	CALL getin_p('D_DOC',D_DOC)
1261	!
1262	!!Config Key = Dif
1263	!Config Desc = diffusion coefficient of POC
1264	!Config if = OK_STOMATE
1265	!Config Def = 1E-4
1266	!Config Help = from Koven et al., 2013 BG
1267	!Config Units = [m^2 day-1]
1268	CALL getin_p('Dif',Dif)
1269	!
1270	!!Config Key = CUE
1271	!Config Desc = Microbial carbon use efficiency
1272	!Config if = OK_STOMATE
1273	!Config Def = 0.5
1274	!Config Help = Assumed value of CUE=0.5
1275	!Config Units = [-]
1276	CALL getin_p('CUE',CUE)
1277
1278	!!-
1279	!! 2. Computational step
1280	!! Loop over time - Call of soilcarbon routine at each time step
1281	!! Updated soil carbon stocks are stored into carbon variable
1282	!! We only keep the last value of carbon variable (no time dimension).
1283	!!-
1284	iyear=1
1285	DO i=1,itau_len
1286	iatt = iatt+1
1287	IF (iatt > nparan*nbyear) THEN
1288	IF (debug) WRITE(numout,*) iyear
1289	iatt = 1
1290	iyear=iyear+1
1291	ENDIF
1292
1293	! Get diaglev from module vertical for CWRR
1294	diaglev=znt(1:nbdl)
1295
1296	CALL soilcarbon &
1297	& (kjpindex, dt_forcesoil, clay, &
1298	& soilcarbon_input(:,:,:,:,:,iatt), control_temp_soil(:,:,:,iatt), control_moist_soil(:,:,:,iatt), &
1299	& carbon, resp_hetero_soil, &
1300	& litter_above(:,:,:,:,iatt),litter_below(:,:,:,:,:,iatt),&
1301	& soilhum, DOC, moist_soil(:,:,iatt), DOC_EXP,&
1302	& lignin_struc_above(:,:,iatt), lignin_struc_below(:,:,:,iatt),&
1303	& floodout(:,iatt), runoff_per_soil(:,:,iatt), drainage_per_soil(:,:,iatt), wat_flux0(:,:,iatt), &
1304	& wat_flux(:,:,:,iatt),bulk_dens, soil_ph,veget_max, soil_mc(:,:,:,iatt), soiltile, &
1305	& DOC_to_topsoil(:,:,iatt), DOC_to_subsoil(:,:,iatt))
1306	ENDDO
1307	WRITE(numout,*) "End of soilcarbon LOOP."
1308
1309	WRITE(numout,*) " carbon pft 4 icarbon layer 3 forcesoil=", carbon(:,:,4,3,icarbon)
1310	WRITE(numout,*) " DOC pft 4 icarbon layer 3 forcesoil=", DOC(:,4,3,ifree,:,icarbon)
1311	WRITE(numout,*) " soilcarbon_input icarbon pft 4 layer 3 forcesoil=",soilcarbon_input(:,4,3,:,icarbon,iatt)
1312	WRITE(numout,*) " carbon pft 4 icarbon13 layer 3 forcesoil=", carbon(:,:,4,3,icarbon13)
1313	WRITE(numout,*) " DOC pft 4 icarbon13 layer 3 forcesoil=", DOC(:,4,3,ifree,:,icarbon13)
1314	WRITE(numout,*) " soilcarbon_input icarbon13 pft 4 layer 3 forcesoil=",soilcarbon_input(:,4,3,:,icarbon13,iatt)
1315
1316
1317	!
1318	!! Gathering of variables towards main processor in parallel mode
1319	!
1320	DO m =1,nvm
1321	DO l=1,nbdl
1322	DO k=1,nelements
1323	CALL Gather(carbon(:,:,m,l,k),carbon_g(:,:,m,l,k))
1324	ENDDO
1325	ENDDO
1326	ENDDO
1327	DO m =1,nbdl
1328	DO l=1,ndoc
1329	DO i= 1,npool
1330	DO k=1,nelements
1331	CALL Gather(DOC(:,:,m,l,i,k),DOC_g(:,:,m,l,i,k))
1332	ENDDO
1333	ENDDO
1334	ENDDO
1335	ENDDO
1336	!!-
1337	!! 3. write new carbon stocks into the ouput restart file
1338	!!-
1339	IF (is_root_prc) THEN
1340	DO k=1,nelements
1341	DO m=1,nvm
1342	WRITE (part_str, '(I2)') m
1343	IF (m<10) part_str(1:1)='0'
1344	var_name = 'carbon_z01_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1345	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1346	& carbon_g(:,:,m,1,k), 'scatter', nbp_glo, indices_g)
1347	ENDDO
1348	DO m=1,nvm
1349	WRITE (part_str, '(I2)') m
1350	IF (m<10) part_str(1:1)='0'
1351	var_name = 'carbon_z02_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1352	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1353	& carbon_g(:,:,m,2,k), 'scatter', nbp_glo, indices_g)
1354	ENDDO
1355
1356	DO m=1,nvm
1357	WRITE (part_str, '(I2)') m
1358	IF (m<10) part_str(1:1)='0'
1359	var_name = 'carbon_z03_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1360	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1361	& carbon_g(:,:,m,3,k), 'scatter', nbp_glo, indices_g)
1362	ENDDO
1363
1364	DO m=1,nvm
1365	WRITE (part_str, '(I2)') m
1366	IF (m<10) part_str(1:1)='0'
1367	var_name = 'carbon_z04_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1368	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1369	& carbon_g(:,:,m,4,k), 'scatter', nbp_glo, indices_g)
1370	ENDDO
1371
1372	DO m=1,nvm
1373	WRITE (part_str, '(I2)') m
1374	IF (m<10) part_str(1:1)='0'
1375	var_name = 'carbon_z05_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1376	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1377	& carbon_g(:,:,m,5,k), 'scatter', nbp_glo, indices_g)
1378	ENDDO
1379
1380	DO m=1,nvm
1381	WRITE (part_str, '(I2)') m
1382	IF (m<10) part_str(1:1)='0'
1383	var_name = 'carbon_z06_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1384	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1385	& carbon_g(:,:,m,6,k), 'scatter', nbp_glo, indices_g)
1386	ENDDO
1387
1388	DO m=1,nvm
1389	WRITE (part_str, '(I2)') m
1390	IF (m<10) part_str(1:1)='0'
1391	var_name = 'carbon_z07_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1392	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1393	& carbon_g(:,:,m,7,k), 'scatter', nbp_glo, indices_g)
1394	ENDDO
1395
1396
1397	DO m=1,nvm
1398	WRITE (part_str, '(I2)') m
1399	IF (m<10) part_str(1:1)='0'
1400	var_name = 'carbon_z08_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1401	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1402	& carbon_g(:,:,m,8,k), 'scatter', nbp_glo, indices_g)
1403	ENDDO
1404
1405	DO m=1,nvm
1406	WRITE (part_str, '(I2)') m
1407	IF (m<10) part_str(1:1)='0'
1408	var_name = 'carbon_z09_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1409	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1410	& carbon_g(:,:,m,9,k), 'scatter', nbp_glo, indices_g)
1411	ENDDO
1412
1413	DO m=1,nvm
1414	WRITE (part_str, '(I2)') m
1415	IF (m<10) part_str(1:1)='0'
1416	var_name = 'carbon_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1417	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1418	& carbon_g(:,:,m,10,k), 'scatter', nbp_glo, indices_g)
1419	ENDDO
1420
1421	DO m=1,nvm
1422	WRITE (part_str, '(I2)') m
1423	IF (m<10) part_str(1:1)='0'
1424	var_name = 'carbon_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1425	CALL restput (rest_id_sto, var_name, nbp_glo, ncarb, 1, itau_dep, &
1426	& carbon_g(:,:,m,11,k), 'scatter', nbp_glo, indices_g)
1427	ENDDO
1428
1429	DO m=1,npool
1430	WRITE (part_str, '(I1)') m
1431	var_name = 'freedoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1432	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1433	& DOC_g(:,:,1,ifree,m,k), 'scatter', nbp_glo, indices_g)
1434	ENDDO
1435
1436	DO m=1,npool
1437	WRITE (part_str, '(I1)') m
1438	var_name = 'freedoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1439	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1440	& DOC_g(:,:,2,ifree,m,k), 'scatter', nbp_glo, indices_g)
1441	ENDDO
1442
1443	DO m=1,npool
1444	WRITE (part_str, '(I1)') m
1445	var_name = 'freedoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1446	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1447	& DOC_g(:,:,3,ifree,m,k), 'scatter', nbp_glo, indices_g)
1448	ENDDO
1449
1450	DO m=1,npool
1451	WRITE (part_str, '(I1)') m
1452	var_name = 'freedoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1453	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1454	& DOC_g(:,:,4,ifree,m,k), 'scatter', nbp_glo, indices_g)
1455	ENDDO
1456
1457	DO m=1,npool
1458	WRITE (part_str, '(I1)') m
1459	var_name = 'freedoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1460	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1461	& DOC_g(:,:,5,ifree,m,k), 'scatter', nbp_glo, indices_g)
1462	ENDDO
1463
1464	DO m=1,npool
1465	WRITE (part_str, '(I1)') m
1466	var_name = 'freedoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1467	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1468	& DOC_g(:,:,6,ifree,m,k), 'scatter', nbp_glo, indices_g)
1469	ENDDO
1470
1471	DO m=1,npool
1472	WRITE (part_str, '(I1)') m
1473	var_name = 'freedoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1474	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1475	& DOC_g(:,:,7,ifree,m,k), 'scatter', nbp_glo, indices_g)
1476	ENDDO
1477
1478	DO m=1,npool
1479	WRITE (part_str, '(I1)') m
1480	var_name = 'freedoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1481	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1482	& DOC_g(:,:,8,ifree,m,k), 'scatter', nbp_glo, indices_g)
1483	ENDDO
1484
1485	DO m=1,npool
1486	WRITE (part_str, '(I1)') m
1487	var_name = 'freedoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1488	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1489	& DOC_g(:,:,9,ifree,m,k), 'scatter', nbp_glo, indices_g)
1490	ENDDO
1491
1492	DO m=1,npool
1493	WRITE (part_str, '(I1)') m
1494	var_name = 'freedoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1495	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1496	& DOC_g(:,:,10,ifree,m,k), 'scatter', nbp_glo, indices_g)
1497	ENDDO
1498
1499	DO m=1,npool
1500	WRITE (part_str, '(I1)') m
1501	var_name = 'freedoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1502	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1503	& DOC_g(:,:,11,ifree,m,k), 'scatter', nbp_glo, indices_g)
1504	ENDDO
1505
1506	DO m=1,npool
1507	WRITE (part_str, '(I1)') m
1508	var_name = 'adsdoc_z1_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1509	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1510	& DOC_g(:,:,1,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1511	ENDDO
1512
1513	DO m=1,npool
1514	WRITE (part_str, '(I1)') m
1515	var_name = 'adsdoc_z2_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1516	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1517	& DOC_g(:,:,2,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1518	ENDDO
1519
1520	DO m=1,npool
1521	WRITE (part_str, '(I1)') m
1522	var_name = 'adsdoc_z3_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1523	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1524	& DOC_g(:,:,3,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1525	ENDDO
1526
1527	DO m=1,npool
1528	WRITE (part_str, '(I1)') m
1529	var_name = 'adsdoc_z4_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1530	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1531	& DOC_g(:,:,4,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1532	ENDDO
1533
1534	DO m=1,npool
1535	WRITE (part_str, '(I1)') m
1536	var_name = 'adsdoc_z5_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1537	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1538	& DOC_g(:,:,5,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1539	ENDDO
1540
1541	DO m=1,npool
1542	WRITE (part_str, '(I1)') m
1543	var_name = 'adsdoc_z6_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1544	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1545	& DOC_g(:,:,6,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1546	ENDDO
1547
1548	DO m=1,npool
1549	WRITE (part_str, '(I1)') m
1550	var_name = 'adsdoc_z7_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1551	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1552	& DOC_g(:,:,7,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1553	ENDDO
1554
1555	DO m=1,npool
1556	WRITE (part_str, '(I1)') m
1557	var_name = 'adsdoc_z8_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1558	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1559	& DOC_g(:,:,8,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1560	ENDDO
1561
1562	DO m=1,npool
1563	WRITE (part_str, '(I1)') m
1564	var_name = 'adsdoc_z9_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1565	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1566	& DOC_g(:,:,9,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1567	ENDDO
1568
1569	DO m=1,npool
1570	WRITE (part_str, '(I1)') m
1571	var_name = 'adsdoc_z10_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1572	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1573	& DOC_g(:,:,10,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1574	ENDDO
1575
1576	DO m=1,npool
1577	WRITE (part_str, '(I1)') m
1578	var_name = 'adsdoc_z11_'//part_str(1:LEN_TRIM(part_str))//element_str(k)
1579	CALL restput (rest_id_sto, var_name, nbp_glo, nvm , 1, itau_dep, &
1580	& DOC_g(:,:,11,iadsorbed,m,k), 'scatter', nbp_glo, indices_g)
1581	ENDDO
1582	ENDDO
1583	!-
1584	CALL getin_dump
1585	CALL restclo
1586	! ier = NF90_CLOSE(rest_id_sto)
1587	! IF (ier /= 0) THEN
1588	! WRITE(numout,*) 'Error closing netcdf file = ',ier
1589	! STOP
1590	! END IF
1591
1592
1593	ENDIF
1594
1595
1596	#ifdef CPP_PARA
1597	CALL MPI_FINALIZE(ier)
1598	#endif
1599	WRITE(numout,*) "End of forcesoil."
1600	!--------------------
1601	END PROGRAM forcesoil

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

Download in other formats: