Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

forcesoil.f90 @ 7325

Last change on this file since 7325 was 4977, checked in by simon.bowring, 6 years ago
Currently running (13/02/2018) version includes all necessarily changes to include DOC in MICT code... further parametrisation necessary to equate soil pools with those of normal forcesoil restarts
File size: 51.7 KB

Line
1	! =================================================================================================================================
2	! PROGRAM : forcesoil
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.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 utils
44	USE defprec
45	USE constantes
46	USE constantes_var
47	USE constantes_mtc
48	USE constantes_soil
49	USE pft_parameters
50	USE stomate_data
51	USE ioipsl_para
52	USE mod_orchidee_para
53	USE stomate_soilcarbon
54	USE stomate_permafrost_soilcarbon
55	USE constantes_soil_var
56	USE stomate
57	USE stomate_io_carbon_permafrost
58	#ifdef CPP_PARA
59	USE mpi
60	#endif
61	!-
62	IMPLICIT NONE
63	!-
64	!-
65	CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out
66	INTEGER(i_std) :: iim,jjm !! Indices (unitless)
67
68	INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless)
69	INTEGER(i_std) :: kjpindex !! Domain size (unitless)
70
71	INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?)
72	!! and number of time steps of the simulation (unitless)
73	CHARACTER(LEN=30) :: time_str !! Length of the simulation (year)
74	REAL(r_std) :: dt_files !! time step between two successive itaus (?)
75	!! (requested by restini routine) (seconds)
76	REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?)
77	INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless)
78	CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file
79	!! (unitless)
80	!-
81	CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless)
82	INTEGER :: Cforcing_id !! ID of the forcing file (unitless)
83	INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file
84	!! (unitless)
85	REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days)
86	INTEGER :: nparan !! Number of values stored per year in the forcing file
87	!! (unitless)
88	INTEGER :: nbyear
89	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless)
90	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless)
91	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless)
92	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined
93	!! in lat/lon (2D) (degrees)
94	REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless)
95
96
97	INTEGER :: i,m,iatt,iv,iyear !! counters (unitless)
98
99	CHARACTER(LEN=100) :: var_name
100	CHARACTER(LEN=8000) :: taboo_vars !! string used for storing the name of the variables
101	!! of the stomate restart file that are not automatically
102	!! duplicated from input to output restart file (unitless)
103	REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless)
104	INTEGER(i_std),PARAMETER :: nbvarmax=1000 !! maximum # of variables assumed in the stomate restart file
105	!! (unitless)
106	INTEGER(i_std) :: nbvar !! # of variables effectively present
107	!! in the stomate restart file (unitless)
108	CHARACTER(LEN=1000),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored
109	!! in the stomate restart file (unitless)
110	INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable
111	!! of the stomate restart file
112	INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable
113	!! of the stomate restart file
114	INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable
115	!! of the stomate restart file
116	LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable
117	!! of the stomate restart file are of length 1 (ie scalar)
118	!! (unitless)
119	REAL(r_std),DIMENSION(:),ALLOCATABLE :: var_2d !! matrix read/written in restget/restput routines (unitless)
120	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless)
121	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: var_4d !! matrix read/written in restget/restput routines (unitless)
122	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: var_5d !! matrix read/written in restget/restput routines (unitless)
123	REAL(r_std),DIMENSION(:,:,:,:,:),ALLOCATABLE :: var_6d !! matrix read/written in restget/restput routines (unitless)
124	REAL(r_std) :: x_tmp !! temporary variable used to store return value
125	INTEGER(i_std) :: orch_vardims !! Orchidee dimensions (different to IOIPSL -exclude time dim-)
126	!! from nf90_get_att (unitless)
127	CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT
128	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless)
129	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
130	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g !! Temperature control (nbp_glo,above/below,time) on OM decomposition
131	!! (unitless)
132	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g !! Moisture control (nbp_glo,abo/below,time) on OM decomposition
133	!! ?? Should be defined per PFT as well (unitless)
134	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$)
135
136	REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless)
137	REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) !! soil carbon input (nbp_loc,ncarb,nvm,time)
138	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
139	REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) !! Temperature control (nbp_loc,above/below,time) on OM decomposition
140	!! (unitless)
141	REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) !! Moisture control (nbp_loc,abo/below,time) on OM decomposition
142	!! ?? Should be defined per PFT as well (unitless)
143	REAL(r_std),ALLOCATABLE :: carbon(:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$)
144	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
145	!! (requested by soilcarbon routine but not used here)
146
147	INTEGER(i_std) :: printlev_loc !! Local write level
148	INTEGER(i_std) :: ier,iret !! Used for hangling errors
149
150	CHARACTER(LEN=50) :: temp_name
151	CHARACTER(LEN=100) :: msg3
152	LOGICAL :: debug !! boolean used for printing messages
153	LOGICAL :: l_error !! boolean for memory allocation
154	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: matrixA !! Carbon fluxes matrix
155	! allocateable arrays needed for permafrost carbon
156	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_a_g !! active carbon concentration
157	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_s_g !! slow carbon concentration
158	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_p_g !! passive carbon concentration
159	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: O2_soil_g !! oxygen in the soil
160	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: CH4_soil_g !! methane in the soil
161	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: O2_snow_g !! oxygen in the snow
162	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: CH4_snow_g !! methane in the snow
163	REAL(r_std),DIMENSION(:),ALLOCATABLE :: z_organic_g !! organic carbon depth
164	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: snowdz_g !! snow depth at each layer
165	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: snowrho_g !! snow density at each layer
166	REAL(r_std),DIMENSION(:),ALLOCATABLE :: zz_deep !! deep vertical profile
167	REAL(r_std),DIMENSION(:),ALLOCATABLE :: zz_coef_deep !! deep vertical profile
168	REAL(r_std),DIMENSION(:),ALLOCATABLE :: z_organic !! organic carbon depth
169	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_a !! active carbon concentration
170	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_s !! slow carbon concentration
171	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: deepC_p !! passive carbon concentration
172	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: O2_soil !! oxygen in the soil
173	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: CH4_soil !! methane in the soil
174	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: O2_snow !! oxygen in the snow
175	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: CH4_snow !! methane in the snow
176	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: pb !! surface pressure
177	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: snow !! snow mass
178	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: tprof !! deep soil temperature profile
179	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: fbact !! factor for soil carbon decomposition
180	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: hslong !! deep soil humidity
181	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: veget_max !! maximum vegetation fraction
182	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: rprof !! PFT rooting depth
183	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: tsurf !! surface temperature
184	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: snowrho !! snow density
185	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: snowdz !! snow depth
186	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lalo !! Geogr. coordinates (latitude,longitude) (degrees)
187	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: heat_Zimov !! heating associated with decomposition [W/m**3 soil]
188	REAL(R_STD), DIMENSION(:), ALLOCATABLE :: sfluxCH4_deep, sfluxCO2_deep !! [g / m**2]
189	REAL(R_STD), DIMENSION(:,:), ALLOCATABLE :: altmax !! active layer thickness (m)
190	REAL(R_STD), DIMENSION(:,:), ALLOCATABLE :: altmax_g !! global active layer thickness (m)
191	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: carbon_surf_g
192	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: carbon_surf !! vertically-integrated (diagnostic) soil carbon pool: active, slow, or passive, (gC/(m**2 of ground))
193	REAL(R_STD), ALLOCATABLE, DIMENSION(:,:) :: fixed_cryoturbation_depth !! depth to hold cryoturbation to for fixed runs
194	LOGICAL, SAVE :: satsoil = .FALSE.
195	LOGICAL :: reset_soilc = .false.
196
197	INTEGER(i_std) :: start_2d(2), count_2d(2)
198	INTEGER(i_std) :: start_4d(4), count_4d(4), start_3d(3), count_3d(3)
199	!_ =================================================================================================================================
200
201	CALL Init_orchidee_para
202	CALL init_timer
203
204	! Set specific write level to forcesoil using WRITELEVEL_forcesoil=[0-4] in run.def.
205	! The global printlev is used as default value.
206	printlev_loc=get_printlev('forcesoil')
207
208	!-
209	! Configure the number of PFTS
210	!-
211	ok_pc=.FALSE.
212	CALL getin_p('OK_PC',ok_pc)
213	! 1. Read the number of PFTs
214	!
215	!Config Key = NVM
216	!Config Desc = number of PFTs
217	!Config If = OK_SECHIBA or OK_STOMATE
218	!Config Def = 13
219	!Config Help = The number of vegetation types define by the user
220	!Config Units = [-]
221	CALL getin_p('NVM',nvm)
222
223	! 2. Allocation
224	ALLOCATE(pft_to_mtc(nvm),stat=ier)
225	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'pft_to_mtc : error in memory allocation', '', '')
226
227	! 3. Initialisation of the correspondance table
228	pft_to_mtc(:) = undef_int
229
230	! 4.Reading of the conrrespondance table in the .def file
231	!
232	!Config Key = PFT_TO_MTC
233	!Config Desc = correspondance array linking a PFT to MTC
234	!Config if = OK_SECHIBA or OK_STOMATE
235	!Config Def = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13
236	!Config Help =
237	!Config Units = [-]
238	CALL getin_p('PFT_TO_MTC',pft_to_mtc)
239
240	! 4.1 if nothing is found, we use the standard configuration
241	IF(nvm <= nvmc ) THEN
242	IF(pft_to_mtc(1) == undef_int) THEN
243	WRITE(numout,*) 'Note to the user : we will use ORCHIDEE to its standard configuration'
244	pft_to_mtc(:) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 /)
245	ENDIF
246	ELSE
247	IF(pft_to_mtc(1) == undef_int) THEN
248	WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop'
249	ENDIF
250	ENDIF
251
252	! 4.2 What happened if pft_to_mtc(j) > nvmc (if the mtc doesn't exist)?
253	DO i = 1, nvm
254	IF(pft_to_mtc(i) > nvmc) THEN
255	CALL ipslerr_p(3, 'forcesoil', 'the MTC you chose doesnt exist', 'we stop reading pft_to_mtc', '')
256	ENDIF
257	ENDDO
258
259	! 4.3 Check if pft_to_mtc(1) = 1
260	IF(pft_to_mtc(1) /= 1) THEN
261	CALL ipslerr_p(3, 'forcesoil', 'the first pft has to be the bare soil', 'we stop reading next values of pft_to_mtc', '')
262	ENDIF
263
264	DO i = 2,nvm
265	IF(pft_to_mtc(i) == 1) THEN
266	CALL ipslerr_p(3, 'forcesoil', 'only pft_to_mtc(1) has to be the bare soil', 'we stop reading next values of pft_to_mtc', '')
267	ENDIF
268	ENDDO
269
270	! 5. Allocate and initialize natural and is_c4
271
272	! 5.1 Memory allocation
273	ALLOCATE(natural(nvm),stat=ier)
274	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'natural : error in memory allocation', '', '')
275
276	ALLOCATE(is_c4(nvm),stat=ier)
277	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'is_c4 : error in memory allocation', '', '')
278
279	ALLOCATE(permafrost_veg_exists(nvm),stat=ier)
280	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'permafrost_veg_exists : error in memory allocation', '', '')
281
282	! 5.2 Initialisation
283	DO i = 1, nvm
284	natural(i) = natural_mtc(pft_to_mtc(i))
285	is_c4(i) = is_c4_mtc(pft_to_mtc(i))
286	ENDDO
287
288	DO i = 1, nvm
289	permafrost_veg_exists(i) = permafrost_veg_exists_mtc(pft_to_mtc(i))
290	ENDDO
291	!!-
292	!! 1. Initialisation stage
293	!! Reading a set of input files, allocating variables and preparing output restart file.
294	!!-
295	! Define restart file name
296	! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var).
297	! User values are used if present in the .def file.
298	! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used.
299	!-
300	IF (is_root_prc) THEN
301	sto_restname_in = 'stomate_start.nc'
302	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
303	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
304	sto_restname_out = 'stomate_rest_out.nc'
305	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
306	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
307	IF (ok_pc) CALL getin ('satsoil', satsoil)
308	!-
309	! Open the input file and Get some Dimension and Attributes ID's
310	!-
311	CALL nccheck( NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto))
312	CALL nccheck( NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g))
313	CALL nccheck( NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g))
314	CALL nccheck( NF90_INQ_VARID (rest_id_sto, "time", iv))
315	CALL nccheck( NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar))
316	CALL nccheck( NF90_CLOSE (rest_id_sto))
317	i=INDEX(thecalendar,ACHAR(0))
318	IF ( i > 0 ) THEN
319	thecalendar(i:20)=' '
320	ENDIF
321	!-
322	! Allocate longitudes and latitudes
323	!-
324	ALLOCATE (lon(iim_g,jjm_g), stat=ier)
325	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lon : error in memory allocation', '', '')
326	ALLOCATE (lat(iim_g,jjm_g), stat=ier)
327	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lat : error in memory allocation', '', '')
328	lon(:,:) = zero
329	lat(:,:) = zero
330	lev(1) = zero
331	!-
332	CALL restini &
333	& (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
334	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto, &
335	& NC_COMPRESSION_ENABLE )
336	ENDIF
337
338	CALL bcast(date0)
339	CALL bcast(thecalendar)
340	WRITE(numout,*) "calendar = ",thecalendar
341	!-
342	! calendar
343	!-
344	CALL ioconf_calendar (thecalendar)
345	CALL ioget_calendar (one_year,one_day)
346	CALL ioconf_startdate(date0)
347	!
348	IF (ok_pc) THEN
349	!- Permafrost variables (zz_deep and zz_coef_deep are constants)
350	ALLOCATE (zz_deep(ndeep), stat=ier)
351	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'zz_deep : error in memory allocation', '', '')
352	ALLOCATE (zz_coef_deep(ndeep), stat=ier)
353	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'zz_coef_deep : error in memory allocation', '', '')
354	ENDIF
355	!-
356	! define forcing file's name (Cforcing_name var)
357	! User value is used if present in the .def file
358	! If not, default (NONE) is used
359	!-
360	IF (.NOT.ok_pc) THEN
361	Cforcing_name = 'NONE'
362	CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
363	ELSE
364	Cforcing_name = 'stomate_Cforcing_permafrost.nc'
365	CALL getin ('STOMATE_CFORCING_PF_NM',Cforcing_name)
366	ENDIF
367	!
368	!! For master process only
369	!
370	IF (is_root_prc) THEN
371	!-
372	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
373	! and allocate variables.
374	!-
375	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id))
376	!-
377	! Total Domain size is stored in nbp_glo variable
378	!-
379	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp))
380	nbp_glo = NINT(x_tmp)
381	!-
382	! Number of values stored per year in the forcing file is stored in nparan var.
383	!-
384	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp))
385	nparan = NINT(x_tmp)
386	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp))
387	nbyear = NINT(x_tmp)
388	!-
389	ALLOCATE (indices_g(nbp_glo), stat=ier)
390	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'indices_g : error in memory allocation', '', '')
391	ALLOCATE (clay_g(nbp_glo), stat=ier)
392	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'clay_g : error in memory allocation', '', '')
393	!-
394	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
395	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'x_indices_g : error in memory allocation', '', '')
396	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'index',v_id))
397	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,x_indices_g))
398	indices_g(:) = NINT(x_indices_g(:))
399	WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
400	DEALLOCATE (x_indices_g)
401	!-
402	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'clay',v_id))
403	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,clay_g))
404	!-
405	IF (ok_pc) THEN
406	!- Permafrost variables (zz_deep and zz_coef_deep are constants)
407	! ALLOCATE (zz_deep(ndeep))
408	! ALLOCATE (zz_coef_deep(ndeep))
409	ALLOCATE (z_organic_g(nbp_glo), stat=ier)
410	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'z_organic_g : error in memory allocation', '', '')
411	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'zz_deep',v_id))
412	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,zz_deep))
413	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'zz_coef_deep',v_id))
414	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,zz_coef_deep))
415	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'z_organic',v_id))
416	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,z_organic_g))
417	ENDIF
418	CALL nccheck( NF90_CLOSE(Cforcing_id) )
419	! time step of forcesoil program (in days)
420	!-
421	dt_forcesoil = one_year / FLOAT(nparan)
422	WRITE(numout,*) 'time step (d): ',dt_forcesoil
423	WRITE(numout,*) 'nparan: ',nparan
424	WRITE(numout,*) 'nbyear: ',nbyear
425	!-
426	! read and write the variables in the output restart file we do not modify within the Forcesoil program
427	! ie all variables stored in the input restart file except those stored in taboo_vars
428	!-
429	IF (.NOT. ok_pc) THEN
430	!-
431	taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
432	& '$day_counter$ $dt_days$ $date$ $carbon$ '
433	!-
434	ELSE
435	!-
436	taboo_vars = '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
437	& '$day_counter$ $dt_days$ $date$ $deepC_a$ $deepC_s$ '// &
438	& '$deepC_p$ $O2_soil$ $CH4_soil$ $O2_snow$ $CH4_snow$ '// &
439	& '$altmax$ '
440	!-
441	ENDIF ! ok_pc
442	!-
443	CALL ioget_vname(rest_id_sto, nbvar, varnames)
444	!-
445	! read and write some special variables (1D or variables that we need)
446	!-
447	var_name = 'day_counter'
448	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
449	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
450	!-
451	var_name = 'dt_days'
452	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
453	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
454	!-
455	var_name = 'date'
456	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
457	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
458	!-
459	DO iv=1,nbvar
460	!-- check if the variable is to be written here
461	IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
462	!---- get variable dimensions, especially 3rd dimension
463	CALL ioget_vdim &
464	& (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
465	l1d = ALL(vardims(1:varnbdim) == 1)
466	WRITE(,) TRIM(varnames(iv)),": ", varnbdim, "-", vardims(1:varnbdim)," l1d=" ,l1d
467	!---- read it
468	IF (l1d) THEN
469	CALL restget &
470	& (rest_id_sto, TRIM(varnames(iv)), 1, 1, &
471	& 1, itau_dep, .TRUE., xtmp)
472	CALL restput &
473	& (rest_id_sto, TRIM(varnames(iv)), 1, 1, &
474	& 1, itau_dep, xtmp)
475	ELSE
476	orch_vardims = varnbdim - 1 ! exclude time dimension introduced by IOIPSL
477	! Deal for different number of dimension
478	IF (orch_vardims == 2) THEN
479	!----
480	! vardims X(1), Y(2), remaining variables, time(last position)
481	ALLOCATE( var_2d(nbp_glo), stat=ier)
482	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_4d : error in memory allocation', '', '')
483	!----
484	CALL restget &
485	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, 1, &
486	& 1, itau_dep, .TRUE., var_2d, "gather", nbp_glo, indices_g)
487	CALL restput &
488	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, 1, &
489	& 1, itau_dep, var_2d, 'scatter', nbp_glo, indices_g)
490	!----
491	DEALLOCATE(var_2d)
492	ELSE IF (orch_vardims == 3) THEN
493	!----
494	! vardims X(1), Y(2), remaining variables, time(last position)
495	ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
496	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_3d : error in memory allocation', '', '')
497	!----
498	CALL restget &
499	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
500	& 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g)
501	CALL restput &
502	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
503	& 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g)
504	!----
505	DEALLOCATE(var_3d)
506	ELSE IF (orch_vardims == 4) THEN
507	!----
508	! vardims X(1), Y(2), remaining variables, time(last position)
509	ALLOCATE( var_4d(nbp_glo,vardims(3),vardims(4)), stat=ier)
510	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_4d : error in memory allocation', '', '')
511	!----
512	CALL restget &
513	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
514	& vardims(4), itau_dep, .TRUE., var_4d, "gather", nbp_glo, indices_g)
515	CALL restput &
516	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
517	& vardims(4), itau_dep, var_4d, 'scatter', nbp_glo, indices_g)
518	!----
519	DEALLOCATE(var_4d)
520	ELSE IF (orch_vardims == 5) THEN
521	!----
522	! vardims X(1), Y(2), remaining variables, time(last position)
523	ALLOCATE( var_5d(nbp_glo,vardims(3),vardims(4),vardims(5)), stat=ier)
524	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_5d : error in memory allocation', '', '')
525	!----
526	CALL restget &
527	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
528	& vardims(4), vardims(5), itau_dep, .TRUE., var_5d, "gather", nbp_glo, indices_g)
529	CALL restput &
530	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
531	& vardims(4), vardims(5), itau_dep, var_5d, 'scatter', nbp_glo, indices_g)
532	!----
533	DEALLOCATE(var_5d)
534	ELSE IF (orch_vardims == 6) THEN
535	!----
536	! vardims X(1), Y(2), remaining variables, time(last position)
537	ALLOCATE( var_6d(nbp_glo,vardims(3),vardims(4),vardims(5),vardims(6)), stat=ier)
538	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_6d : error in memory allocation', '', '')
539	!----
540	CALL restget &
541	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
542	& vardims(4), vardims(5), vardims(6), itau_dep, .TRUE., var_6d, "gather", nbp_glo, indices_g)
543	CALL restput &
544	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
545	& vardims(4), vardims(5), vardims(6), itau_dep, var_6d, 'scatter', nbp_glo, indices_g)
546	!----
547	DEALLOCATE(var_6d)
548	ELSE
549	WRITE( msg3, '(i5)' ) orch_vardims
550	CALL ipslerr(3, 'forcesoil', 'Varialbes(1) Restart Read/Write not implement for N dimensions(2)', &
551	& TRIM(varnames(iv)), TRIM(msg3))
552	ENDIF
553	ENDIF
554	ENDIF
555	ENDDO
556	! Length of the run (in Years)
557	! User value is used if present in the .def file
558	! If not, default value (10000 Years) is used
559	!-
560	WRITE(time_str,'(a)') '10000Y'
561	CALL getin('TIME_LENGTH', time_str)
562	write(numout,*) 'Number of years for carbon spinup : ',time_str
563	! transform into itau
564	CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
565	write(numout,*) 'Number of time steps to do: ',itau_len
566
567	! read soil carbon stocks values stored in the input restart file
568	!-
569	IF (.NOT. ok_pc) THEN
570	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm), stat=ier)
571	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_g : error in memory allocation', '', '')
572	carbon_g(:,:,:) = val_exp
573	CALL restget &
574	& (rest_id_sto, 'carbon', nbp_glo, ncarb , nvm, itau_dep, &
575	& .TRUE., carbon_g, 'gather', nbp_glo, indices_g)
576	IF (ALL(carbon_g == val_exp)) carbon_g = zero
577	WRITE(numout,*) "date0 : ",date0, itau_dep
578	ELSE
579	!-
580	! Permafrost carbon
581	!-
582	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm), stat=ier)
583	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_g : error in memory allocation', '', '')
584	carbon_g(:,:,:) = 0.
585	ALLOCATE(carbon_surf_g(nbp_glo,ncarb,nvm), stat=ier)
586	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_surf_g : error in memory allocation', '', '')
587	carbon_surf_g(:,:,:) = 0.
588	ALLOCATE(deepC_a_g(nbp_glo,ndeep,nvm), stat=ier)
589	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_a_g : error in memory allocation', '', '')
590	ALLOCATE(deepC_s_g(nbp_glo,ndeep,nvm), stat=ier)
591	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_s_g : error in memory allocation', '', '')
592	ALLOCATE(deepC_p_g(nbp_glo,ndeep,nvm), stat=ier)
593	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_p_g : error in memory allocation', '', '')
594	ALLOCATE(O2_soil_g(nbp_glo,ndeep,nvm), stat=ier)
595	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'O2_soil_g : error in memory allocation', '', '')
596	ALLOCATE(CH4_soil_g(nbp_glo,ndeep,nvm), stat=ier)
597	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'CH4_soil_g : error in memory allocation', '', '')
598	ALLOCATE(O2_snow_g(nbp_glo,nsnow,nvm), stat=ier)
599	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'O2_snow_g : error in memory allocation', '', '')
600	ALLOCATE(CH4_snow_g(nbp_glo,nsnow,nvm), stat=ier)
601	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'CH4_snow_g : error in memory allocation', '', '')
602	ALLOCATE(altmax_g(nbp_glo,nvm), stat=ier)
603	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'altmax_g : error in memory allocation', '', '')
604
605	deepC_a_g(:,:,:) = val_exp
606	CALL restget (rest_id_sto, 'deepC_a', nbp_glo, ndeep, nvm, itau_dep, &
607	& .TRUE., deepC_a_g, 'gather', nbp_glo, indices_g)
608	IF (ALL(deepC_a_g == val_exp)) deepC_a_g = zero
609
610	deepC_s_g(:,:,:) = val_exp
611	CALL restget (rest_id_sto, 'deepC_s', nbp_glo, ndeep, nvm, itau_dep, &
612	& .TRUE., deepC_s_g, 'gather', nbp_glo, indices_g)
613	IF (ALL(deepC_s_g == val_exp)) deepC_s_g = zero
614
615	deepC_p_g(:,:,:) = val_exp
616	CALL restget (rest_id_sto, 'deepC_p', nbp_glo, ndeep, nvm, itau_dep, &
617	& .TRUE., deepC_p_g, 'gather', nbp_glo, indices_g)
618	IF (ALL(deepC_p_g == val_exp)) deepC_p_g = zero
619
620	var_name= 'altmax'
621	altmax_g(:,:) = val_exp
622	CALL restget (rest_id_sto, var_name, nbp_glo, nvm, 1, itau_dep, .TRUE., altmax_g, "gather", nbp_glo, indices_g)
623	IF ( ALL( altmax_g(:,:) .EQ. val_exp ) ) THEN
624	CALL ipslerr(3, 'forcesoil', 'altmax is not found in stomate restart file', '', '')
625	END IF
626
627	CALL getin('reset_soilc', reset_soilc)
628	IF (reset_soilc) THEN
629	CALL ipslerr(1, 'forcesoil', 'deepC_a, deepC_s and deeC_p', &
630	'are ignored and set to zero value due to', 'reset_soilc option')
631	deepC_a_g(:,:,:) = zero
632	deepC_s_g(:,:,:) = zero
633	deepC_p_g(:,:,:) = zero
634	ENDIF
635
636	O2_soil_g(:,:,:) = val_exp
637	CALL restget (rest_id_sto, 'O2_soil', nbp_glo, ndeep, nvm, itau_dep, &
638	& .TRUE., O2_soil_g, 'gather', nbp_glo, indices_g)
639	IF (ALL(O2_soil_g == val_exp)) O2_soil_g = O2_init_conc
640
641	CH4_soil_g(:,:,:) = val_exp
642	CALL restget (rest_id_sto, 'CH4_soil', nbp_glo, ndeep, nvm, itau_dep, &
643	& .TRUE., CH4_soil_g, 'gather', nbp_glo, indices_g)
644	IF (ALL(CH4_soil_g == val_exp)) CH4_soil_g = CH4_init_conc
645
646	O2_snow_g(:,:,:) = val_exp
647	CALL restget (rest_id_sto, 'O2_snow', nbp_glo, nsnow, nvm, itau_dep, &
648	& .TRUE., O2_snow_g, 'gather', nbp_glo, indices_g)
649	IF (ALL(O2_snow_g == val_exp)) O2_snow_g = O2_init_conc
650
651	CH4_snow_g(:,:,:) = val_exp
652	CALL restget (rest_id_sto, 'CH4_snow', nbp_glo, nsnow, nvm, itau_dep, &
653	& .TRUE., CH4_snow_g, 'gather', nbp_glo, indices_g)
654	IF (ALL(CH4_snow_g == val_exp)) CH4_snow_g = CH4_init_conc
655	ENDIF ! ok_pc
656	ENDIF ! is_root_prc
657	!
658	CALL bcast(nbp_glo)
659	CALL bcast(iim_g)
660	CALL bcast(jjm_g)
661	IF (.NOT. ALLOCATED(indices_g)) ALLOCATE (indices_g(nbp_glo))
662	CALL bcast(indices_g)
663	CALL bcast(nparan)
664	CALL bcast(nbyear)
665	CALL bcast(dt_forcesoil)
666	CALL bcast(itau_dep)
667	CALL bcast(itau_len)
668	CALL bcast(zz_deep)
669	CALL bcast(zz_coef_deep)
670	!
671	! we must initialize data_para :
672	CALL init_orchidee_data_para_driver(nbp_glo,indices_g)
673
674	kjpindex=nbp_loc
675	jjm=jj_nb
676	iim=iim_g
677	IF (printlev_loc>=3) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
678	!-
679	! Analytical spinup is set to false
680	!
681	spinup_analytic = .FALSE.
682	!-
683	! read soil carbon inputs, water and temperature stresses on OM
684	! decomposition
685	! into the forcing file - We read an average year.
686	!-
687	IF (.NOT. ok_pc) THEN
688	!-
689	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
690	! and allocate variables.
691	!-
692	#ifdef CPP_PARA
693	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name),IOR(NF90_NOWRITE, NF90_MPIIO),Cforcing_id, &
694	& comm = MPI_COMM_WORLD, info = MPI_INFO_NULL ))
695	#else
696	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name), NF90_NOWRITE,Cforcing_id) )
697	#endif
698	ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear), stat=ier)
699	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'control_temp : error in memory allocation', '', '')
700	ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear), stat=ier)
701	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'control_moist : error in memory allocation', '', '')
702	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear), stat=ier)
703	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'soilcarbon_input : error in memory allocation', '', '')
704	!-
705	start_4d = (/ nbp_mpi_para_begin(mpi_rank), 1, 1, 1 /)
706	count_4d = (/ nbp_mpi_para(mpi_rank), ncarb, nvm, nparan*nbyear /)
707	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id))
708	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input, &
709	& start = start_4d, count = count_4d ))
710	!
711	start_3d = (/ nbp_mpi_para_begin(mpi_rank), 1, 1 /)
712	count_3d = (/ nbp_mpi_para(mpi_rank), nlevs, nparan*nbyear /)
713	CALL nccheck( NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id))
714	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,control_moist, &
715	& start = start_3d, count = count_3d))
716	CALL nccheck( NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id))
717	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,control_temp, &
718	& start = start_3d, count = count_3d))
719	!- Close Netcdf carbon permafrost file reference
720	CALL nccheck( NF90_CLOSE (Cforcing_id))
721	ELSE
722	!-
723	! Read permafrost-related soil carbon from Cforcing_id
724	!-
725	ALLOCATE(pb(kjpindex,nparan*nbyear), stat=ier)
726	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'pb : error in memory allocation', '', '')
727	ALLOCATE(snow(kjpindex,nparan*nbyear), stat=ier)
728	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snow : error in memory allocation', '', '')
729	ALLOCATE(tprof(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
730	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'tprof : error in memory allocation', '', '')
731	ALLOCATE(fbact(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
732	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'fbact : error in memory allocation', '', '')
733	ALLOCATE(hslong(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
734	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'hslong : error in memory allocation', '', '')
735	ALLOCATE(veget_max(kjpindex,nvm,nparan*nbyear), stat=ier)
736	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'veget_max : error in memory allocation', '', '')
737	ALLOCATE(rprof(kjpindex,nvm,nparan*nbyear), stat=ier)
738	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'rprof : error in memory allocation', '', '')
739	ALLOCATE(tsurf(kjpindex,nparan*nbyear), stat=ier)
740	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'tsurf : error in memory allocation', '', '')
741	ALLOCATE(lalo(kjpindex,2), stat=ier)
742	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lalo : error in memory allocation', '', '')
743	ALLOCATE(snowdz(kjpindex,nsnow,nparan*nbyear), stat=ier)
744	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snowdz_ : error in memory allocation', '', '')
745	ALLOCATE(snowrho(kjpindex,nsnow,nparan*nbyear), stat=ier)
746	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snowrho : error in memory allocation', '', '')
747	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear), stat=ier)
748	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'soilcarbon_input : error in memory allocation', '', '')
749	!-
750	CALL stomate_io_carbon_permafrost_read(Cforcing_name, nparan, nbyear,&
751	nbp_mpi_para_begin(mpi_rank), nbp_mpi_para(mpi_rank), &
752	soilcarbon_input, pb, snow, tsurf, &
753	tprof, fbact, hslong, rprof, &
754	lalo, snowdz, snowrho, veget_max )
755	ENDIF ! ok_pc
756	!---
757	!--- Create the index table
758	!---
759	!--- This job returns a LOCAL kindex.
760	!---
761	ALLOCATE (indices(kjpindex),stat=ier)
762	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'indices : error in memory allocation', '', '')
763	!
764	!! scattering to all processes in parallel mode
765	!
766	CALL scatter(indices_g,indices)
767	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
768	IF (printlev_loc>=3) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
769	!
770	! Initialize _index variables
771	CALL stomate_init_index(nbp_glo, kjpindex, indices)
772	!-
773	! Allocation of the variables for a processor
774	!-
775	ALLOCATE(clay(kjpindex), stat=ier)
776	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'clay : error in memory allocation', '', '')
777	ALLOCATE(carbon(kjpindex,ncarb,nvm), stat=ier)
778	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'indices : error in memory allocation', '', '')
779	!-
780	IF (.NOT. ok_pc) THEN
781	ALLOCATE(resp_hetero_soil(kjpindex,nvm), stat=ier)
782	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'resp_hetero_soil : error in memory allocation', '', '')
783	ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools), stat=ier)
784	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'matrixA : error in memory allocation', '', '')
785	DO i = 1,nbpools
786	matrixA(:,:,i,i) = un
787	ENDDO
788	ELSE
789	ALLOCATE(carbon_surf(kjpindex,ncarb,nvm), stat=ier)
790	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'carbon_surf : error in memory allocation', '', '')
791	ALLOCATE(deepC_a(kjpindex,ndeep,nvm), stat=ier)
792	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_a : error in memory allocation', '', '')
793	ALLOCATE(deepC_s(kjpindex,ndeep,nvm), stat=ier)
794	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_s : error in memory allocation', '', '')
795	ALLOCATE(deepC_p(kjpindex,ndeep,nvm), stat=ier)
796	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_p : error in memory allocation', '', '')
797	ALLOCATE(O2_soil(kjpindex,ndeep,nvm), stat=ier)
798	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'O2_soil : error in memory allocation', '', '')
799	ALLOCATE(CH4_soil(kjpindex,ndeep,nvm), stat=ier)
800	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'CH4_soil : error in memory allocation', '', '')
801	ALLOCATE(O2_snow(kjpindex,nsnow,nvm), stat=ier)
802	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'O2_snow : error in memory allocation', '', '')
803	ALLOCATE(CH4_snow(kjpindex,nsnow,nvm), stat=ier)
804	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'CH4_snow : error in memory allocation', '', '')
805	ALLOCATE(altmax(kjpindex,nvm), stat=ier)
806	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'altmax : error in memory allocation', '', '')
807	ALLOCATE(z_organic(kjpindex), stat=ier)
808	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'z_organic : error in memory allocation', '', '')
809	ENDIF
810	!-
811	! Initialization of the variables for a processor
812	!-
813	CALL Scatter(clay_g,clay)
814	CALL Scatter(carbon_g,carbon)
815	IF (ok_pc) THEN
816	CALL Scatter(z_organic_g,z_organic)
817	CALL Scatter(carbon_surf_g,carbon_surf)
818	CALL Scatter(deepC_a_g,deepC_a)
819	CALL Scatter(deepC_s_g,deepC_s)
820	CALL Scatter(deepC_p_g,deepC_p)
821	CALL Scatter(O2_soil_g,O2_soil)
822	CALL Scatter(CH4_soil_g,CH4_soil)
823	CALL Scatter(O2_snow_g,O2_snow)
824	CALL Scatter(CH4_snow_g,CH4_snow)
825	CALL Scatter(altmax_g,altmax)
826	ENDIF
827	!-
828	! Configuration of the parameters
829	!-
830	!Config Key = FRAC_CARB_AP
831	!Config Desc = frac carb coefficients from active pool: depends on clay
832	!content
833	!Config if = OK_STOMATE
834	!Config Def = 0.004
835	!Config Help = fraction of the active pool going to the passive pool
836	!Config Units = [-]
837	CALL getin_p('FRAC_CARB_AP',frac_carb_ap)
838	!
839	!Config Key = FRAC_CARB_SA
840	!Config Desc = frac_carb_coefficients from slow pool
841	!Config if = OK_STOMATE
842	!Config Def = 0.42
843	!Config Help = fraction of the slow pool going to the active pool
844	!Config Units = [-]
845	CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
846	!
847	!Config Key = FRAC_CARB_SP
848	!Config Desc = frac_carb_coefficients from slow pool
849	!Config if = OK_STOMATE
850	!Config Def = 0.03
851	!Config Help = fraction of the slow pool going to the passive pool
852	!Config Units = [-]
853	CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
854	!
855	!Config Key = FRAC_CARB_PA
856	!Config Desc = frac_carb_coefficients from passive pool
857	!Config if = OK_STOMATE
858	!Config Def = 0.45
859	!Config Help = fraction of the passive pool going to the passive pool
860	!Config Units = [-]
861	CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
862	!
863	!Config Key = FRAC_CARB_PS
864	!Config Desc = frac_carb_coefficients from passive pool
865	!Config if = OK_STOMATE
866	!Config Def = 0.0
867	!Config Help = fraction of the passive pool going to the passive pool
868	!Config Units = [-]
869	CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
870	!
871	!Config Key = ACTIVE_TO_PASS_CLAY_FRAC
872	!Config Desc =
873	!Config if = OK_STOMATE
874	!Config Def = .68
875	!Config Help =
876	!Config Units = [-]
877	CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
878	!
879	!Config Key = CARBON_TAU_IACTIVE
880	!Config Desc = residence times in carbon pools
881	!Config if = OK_STOMATE
882	!Config Def = 0.149
883	!Config Help =
884	!Config Units = [days]
885	CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
886	!
887	!Config Key = CARBON_TAU_ISLOW
888	!Config Desc = residence times in carbon pools
889	!Config if = OK_STOMATE
890	!Config Def = 5.48
891	!Config Help =
892	!Config Units = [days]
893	CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
894	!
895	!Config Key = CARBON_TAU_IPASSIVE
896	!Config Desc = residence times in carbon pools
897	!Config if = OK_STOMATE
898	!Config Def = 241.
899	!Config Help =
900	!Config Units = [days]
901	CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
902	!
903	!Config Key = FLUX_TOT_COEFF
904	!Config Desc =
905	!Config if = OK_STOMATE
906	!Config Def = 1.2, 1.4,.75
907	!Config Help =
908	!Config Units = [days]
909	CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)
910
911	!
912	!! 2. Computational step
913	!! Loop over time - Call of soilcarbon routine at each time step
914	!! Updated soil carbon stocks are stored into carbon variable
915	!! We only keep the last value of carbon variable (no time dimension).
916	!!-
917	IF (.NOT.ok_pc) THEN
918	iyear=1
919	iatt = 0
920	DO i=1,itau_len
921	iatt = iatt+1
922	IF (iatt > nparan*nbyear) THEN
923	IF (printlev>=3) WRITE(numout,*) iyear
924	iatt = 1
925	iyear=iyear+1
926	ENDIF
927	CALL soilcarbon &
928	& (kjpindex, dt_forcesoil, clay, &
929	& soilcarbon_input(:,:,:,iatt), &
930	& control_temp(:,:,iatt), control_moist(:,:,iatt), &
931	& carbon, resp_hetero_soil, &
932	& matrixA)
933	ENDDO
934	WRITE(numout,*) "End of soilcarbon LOOP."
935	ELSE
936	IF ( satsoil ) hslong(:,:,:,:) = 1.
937	!these variables are only ouputs from deep_carbcycle (thus not necessary for
938	!Gather and Scatter)
939	ALLOCATE(heat_Zimov(kjpindex,ndeep,nvm), stat=ier)
940	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'heat_Zimov : error in memory allocation', '', '')
941	ALLOCATE(sfluxCH4_deep(kjpindex), stat=ier)
942	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'sfluxCH4_deep : error in memory allocation', '', '')
943	ALLOCATE(sfluxCO2_deep(kjpindex), stat=ier)
944	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'sfluxCO2_deep : error in memory allocation', '', '')
945	ALLOCATE(fixed_cryoturbation_depth(kjpindex,nvm), stat=ier)
946	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'fixed_cryoturbation_depth : error in memory allocation', '', '')
947	ALLOCATE(resp_hetero_soil(kjpindex, nvm), stat=ier)
948	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'resp_hetero_soil : error in memory allocation', '', '')
949
950	iatt = 0
951	iyear=1
952	DO i=1,itau_len
953	iatt = iatt+1
954	IF (iatt > nparan*nbyear) THEN
955	IF (printlev>=3) WRITE(numout,*) iyear
956	iatt = 1
957	iyear=iyear+1
958	ENDIF
959	WRITE(numout, *) "Forcesoil:: deep_carbcycle, iyear=", iyear
960	CALL deep_carbcycle(kjpindex, indices, iatt, dt_forcesoil*one_day, lalo, clay, &
961	tsurf(:,iatt), tprof(:,:,:,iatt), hslong(:,:,:,iatt), snow(:,iatt), heat_Zimov, pb(:,iatt), &
962	sfluxCH4_deep, sfluxCO2_deep, &
963	deepC_a, deepC_s, deepC_p, O2_soil, CH4_soil, O2_snow, CH4_snow, &
964	zz_deep, zz_coef_deep, z_organic, soilcarbon_input(:,:,:,iatt), &
965	veget_max(:,:,iatt), rprof(:,:,iatt), altmax, carbon, carbon_surf, resp_hetero_soil, &
966	fbact(:,:,:,iatt), fixed_cryoturbation_depth, snowdz(:,:,iatt), snowrho(:,:,iatt))
967	ENDDO
968
969	ENDIF
970	!!-
971	!! 3. write new carbon stocks into the ouput restart file
972	!!-
973	CALL restput_p (rest_id_sto, 'carbon', nbp_glo, ncarb , nvm, itau_dep, &
974	& carbon, 'scatter', nbp_glo, indices_g)
975
976	IF (ok_pc) THEN
977	CALL restput_p (rest_id_sto, 'deepC_a', nbp_glo, ndeep, nvm, itau_dep, &
978	& deepC_a, 'scatter', nbp_glo, indices_g)
979	CALL restput_p (rest_id_sto, 'deepC_s', nbp_glo, ndeep, nvm, itau_dep, &
980	& deepC_s, 'scatter', nbp_glo, indices_g)
981	CALL restput_p (rest_id_sto, 'deepC_p', nbp_glo, ndeep, nvm, itau_dep, &
982	& deepC_p, 'scatter', nbp_glo, indices_g)
983	CALL restput_p (rest_id_sto, 'O2_soil', nbp_glo, ndeep, nvm, itau_dep, &
984	& O2_soil, 'scatter', nbp_glo, indices_g)
985	CALL restput_p (rest_id_sto, 'CH4_soil', nbp_glo, ndeep, nvm, itau_dep, &
986	& CH4_soil, 'scatter', nbp_glo, indices_g)
987	CALL restput_p (rest_id_sto, 'O2_snow', nbp_glo, nsnow, nvm, itau_dep, &
988	& O2_snow, 'scatter', nbp_glo, indices_g)
989	CALL restput_p (rest_id_sto, 'CH4_snow', nbp_glo, nsnow, nvm, itau_dep, &
990	& CH4_snow, 'scatter', nbp_glo, indices_g)
991	CALL restput_p (rest_id_sto, 'altmax', nbp_glo, nvm, 1, itau_dep, &
992	& altmax, 'scatter', nbp_glo, indices_g)
993	ENDIF
994	!-
995	IF (is_root_prc) THEN
996	!- Close restart files
997	CALL getin_dump
998	CALL restclo
999	ENDIF
1000	!-
1001	#ifdef CPP_PARA
1002	CALL MPI_FINALIZE(ier)
1003	#endif
1004	WRITE(numout,*) "End of forcesoil."
1005	!--------------------
1006	END PROGRAM forcesoil

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_driver/forcesoil.f90 @ 7325

Download in other formats: