Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

forcesoil.f90

Last change on this file was 6849, checked in by yidi.xu, 4 years ago
ORCHIDEE-MICT-OP for oil palm growth modelling
File size: 52.9 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	INTEGER(i_std) :: nc_opts
201	!_ =================================================================================================================================
202
203	CALL Init_orchidee_para
204	CALL init_timer
205
206	! Set specific write level to forcesoil using WRITELEVEL_forcesoil=[0-4] in run.def.
207	! The global printlev is used as default value.
208	printlev_loc=get_printlev('forcesoil')
209
210	!-
211	! Configure the number of PFTS
212	!-
213	ok_pc=.FALSE.
214	CALL getin_p('OK_PC',ok_pc)
215	! 1. Read the number of PFTs
216	!
217	!Config Key = NVM
218	!Config Desc = number of PFTs
219	!Config If = OK_SECHIBA or OK_STOMATE
220	!Config Def = 13
221	!Config Help = The number of vegetation types define by the user
222	!Config Units = [-]
223	CALL getin_p('NVM',nvm)
224
225	! 2. Allocation
226	ALLOCATE(pft_to_mtc(nvm),stat=ier)
227	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'pft_to_mtc : error in memory allocation', '', '')
228
229	! 3. Initialisation of the correspondance table
230	pft_to_mtc(:) = undef_int
231
232	! 4.Reading of the conrrespondance table in the .def file
233	!
234	!Config Key = PFT_TO_MTC
235	!Config Desc = correspondance array linking a PFT to MTC
236	!Config if = OK_SECHIBA or OK_STOMATE
237	!Config Def = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13
238	!Config Help =
239	!Config Units = [-]
240	CALL getin_p('PFT_TO_MTC',pft_to_mtc)
241
242	! 4.1 if nothing is found, we use the standard configuration
243	IF(nvm <= nvmc ) THEN
244	IF(pft_to_mtc(1) == undef_int) THEN
245	WRITE(numout,*) 'Note to the user : we will use ORCHIDEE to its standard configuration'
246	pft_to_mtc(:) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 /)
247	ENDIF
248	ELSE
249	IF(pft_to_mtc(1) == undef_int) THEN
250	WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop'
251	ENDIF
252	ENDIF
253
254	! 4.2 What happened if pft_to_mtc(j) > nvmc (if the mtc doesn't exist)?
255	DO i = 1, nvm
256	IF(pft_to_mtc(i) > nvmc) THEN
257	CALL ipslerr_p(3, 'forcesoil', 'the MTC you chose doesnt exist', 'we stop reading pft_to_mtc', '')
258	ENDIF
259	ENDDO
260
261	! 4.3 Check if pft_to_mtc(1) = 1
262	IF(pft_to_mtc(1) /= 1) THEN
263	CALL ipslerr_p(3, 'forcesoil', 'the first pft has to be the bare soil', 'we stop reading next values of pft_to_mtc', '')
264	ENDIF
265
266	DO i = 2,nvm
267	IF(pft_to_mtc(i) == 1) THEN
268	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', '')
269	ENDIF
270	ENDDO
271
272	! 5. Allocate and initialize natural and is_c4
273
274	! 5.1 Memory allocation
275	ALLOCATE(natural(nvm),stat=ier)
276	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'natural : error in memory allocation', '', '')
277
278	ALLOCATE(is_c4(nvm),stat=ier)
279	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'is_c4 : error in memory allocation', '', '')
280
281	ALLOCATE(permafrost_veg_exists(nvm),stat=ier)
282	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'permafrost_veg_exists : error in memory allocation', '', '')
283
284	! 5.2 Initialisation
285	DO i = 1, nvm
286	natural(i) = natural_mtc(pft_to_mtc(i))
287	is_c4(i) = is_c4_mtc(pft_to_mtc(i))
288	ENDDO
289
290	DO i = 1, nvm
291	permafrost_veg_exists(i) = permafrost_veg_exists_mtc(pft_to_mtc(i))
292	ENDDO
293
294	!!!! yidi Allocate and initialize is_oilpalm/is_oilpalm_ffbharvest
295	ALLOCATE(is_oilpalm(nvm),stat=ier)
296	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'is_oilpalm : error in memory allocation', '', '')
297
298	DO i = 1, nvm
299	is_oilpalm(i) = is_oilpalm_mtc(pft_to_mtc(i))
300	ENDDO
301
302	ALLOCATE(is_oilpalm_ffbharvest(nvm),stat=ier)
303	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'is_oilpalm_ffbharvest : error in memory allocation', '', '')
304
305	DO i = 1, nvm
306	is_oilpalm_ffbharvest(i) = is_oilpalm_ffbharvest_mtc(pft_to_mtc(i))
307	ENDDO
308	!!!! yidi
309
310	!!-
311	!! 1. Initialisation stage
312	!! Reading a set of input files, allocating variables and preparing output restart file.
313	!!-
314	! Define restart file name
315	! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var).
316	! User values are used if present in the .def file.
317	! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used.
318	!-
319	IF (is_root_prc) THEN
320	sto_restname_in = 'stomate_start.nc'
321	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
322	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
323	sto_restname_out = 'stomate_rest_out.nc'
324	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
325	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
326	IF (ok_pc) CALL getin ('satsoil', satsoil)
327	!-
328	! Open the input file and Get some Dimension and Attributes ID's
329	!-
330	CALL nccheck( NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto))
331	CALL nccheck( NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g))
332	CALL nccheck( NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g))
333	CALL nccheck( NF90_INQ_VARID (rest_id_sto, "time", iv))
334	CALL nccheck( NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar))
335	CALL nccheck( NF90_CLOSE (rest_id_sto))
336	i=INDEX(thecalendar,ACHAR(0))
337	IF ( i > 0 ) THEN
338	thecalendar(i:20)=' '
339	ENDIF
340	!-
341	! Allocate longitudes and latitudes
342	!-
343	ALLOCATE (lon(iim_g,jjm_g), stat=ier)
344	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lon : error in memory allocation', '', '')
345	ALLOCATE (lat(iim_g,jjm_g), stat=ier)
346	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lat : error in memory allocation', '', '')
347	lon(:,:) = zero
348	lat(:,:) = zero
349	lev(1) = zero
350	!-
351	CALL restini &
352	& (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
353	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto, &
354	& use_compression=NC_COMPRESSION_ENABLE )
355	ENDIF
356
357	CALL bcast(date0)
358	CALL bcast(thecalendar)
359	WRITE(numout,*) "calendar = ",thecalendar
360	!-
361	! calendar
362	!-
363	CALL ioconf_calendar (thecalendar)
364	CALL ioget_calendar (one_year,one_day)
365	CALL ioconf_startdate(date0)
366	!
367	IF (ok_pc) THEN
368	!- Permafrost variables (zz_deep and zz_coef_deep are constants)
369	ALLOCATE (zz_deep(ndeep), stat=ier)
370	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'zz_deep : error in memory allocation', '', '')
371	ALLOCATE (zz_coef_deep(ndeep), stat=ier)
372	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'zz_coef_deep : error in memory allocation', '', '')
373	ENDIF
374	!-
375	! define forcing file's name (Cforcing_name var)
376	! User value is used if present in the .def file
377	! If not, default (NONE) is used
378	!-
379	IF (.NOT.ok_pc) THEN
380	Cforcing_name = 'NONE'
381	CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
382	ELSE
383	Cforcing_name = 'stomate_Cforcing_permafrost.nc'
384	CALL getin ('STOMATE_CFORCING_PF_NM',Cforcing_name)
385	ENDIF
386	!
387	!! For master process only
388	!
389	IF (is_root_prc) THEN
390	!-
391	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
392	! and allocate variables.
393	!-
394	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name),IOR(NF90_NOWRITE,NF90_NETCDF4),Cforcing_id))
395	!-
396	! Total Domain size is stored in nbp_glo variable
397	!-
398	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp))
399	nbp_glo = NINT(x_tmp)
400	!-
401	! Number of values stored per year in the forcing file is stored in nparan var.
402	!-
403	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp))
404	nparan = NINT(x_tmp)
405	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp))
406	nbyear = NINT(x_tmp)
407	!-
408	ALLOCATE (indices_g(nbp_glo), stat=ier)
409	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'indices_g : error in memory allocation', '', '')
410	ALLOCATE (clay_g(nbp_glo), stat=ier)
411	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'clay_g : error in memory allocation', '', '')
412	!-
413	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
414	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'x_indices_g : error in memory allocation', '', '')
415	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'index',v_id))
416	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,x_indices_g))
417	indices_g(:) = NINT(x_indices_g(:))
418	WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
419	DEALLOCATE (x_indices_g)
420	!-
421	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'clay',v_id))
422	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,clay_g))
423	!-
424	IF (ok_pc) THEN
425	!- Permafrost variables (zz_deep and zz_coef_deep are constants)
426	! ALLOCATE (zz_deep(ndeep))
427	! ALLOCATE (zz_coef_deep(ndeep))
428	ALLOCATE (z_organic_g(nbp_glo), stat=ier)
429	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'z_organic_g : error in memory allocation', '', '')
430	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'zz_deep',v_id))
431	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,zz_deep))
432	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'zz_coef_deep',v_id))
433	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,zz_coef_deep))
434	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'z_organic',v_id))
435	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,z_organic_g))
436	ENDIF
437	CALL nccheck( NF90_CLOSE(Cforcing_id) )
438	! time step of forcesoil program (in days)
439	!-
440	dt_forcesoil = one_year / FLOAT(nparan)
441	WRITE(numout,*) 'time step (d): ',dt_forcesoil
442	WRITE(numout,*) 'nparan: ',nparan
443	WRITE(numout,*) 'nbyear: ',nbyear
444	!-
445	! read and write the variables in the output restart file we do not modify within the Forcesoil program
446	! ie all variables stored in the input restart file except those stored in taboo_vars
447	!-
448	IF (.NOT. ok_pc) THEN
449	!-
450	taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
451	& '$day_counter$ $dt_days$ $date$ $carbon$ '
452	!-
453	ELSE
454	!-
455	taboo_vars = '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
456	& '$day_counter$ $dt_days$ $date$ $deepC_a$ $deepC_s$ '// &
457	& '$deepC_p$ $O2_soil$ $CH4_soil$ $O2_snow$ $CH4_snow$ '// &
458	& '$altmax$ '
459	!-
460	ENDIF ! ok_pc
461	!-
462	CALL ioget_vname(rest_id_sto, nbvar, varnames)
463	!-
464	! read and write some special variables (1D or variables that we need)
465	!-
466	var_name = 'day_counter'
467	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
468	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
469	!-
470	var_name = 'dt_days'
471	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
472	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
473	!-
474	var_name = 'date'
475	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
476	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
477	!-
478	DO iv=1,nbvar
479	!-- check if the variable is to be written here
480	IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
481	!---- get variable dimensions, especially 3rd dimension
482	CALL ioget_vdim &
483	& (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
484	l1d = ALL(vardims(1:varnbdim) == 1)
485	WRITE(,) TRIM(varnames(iv)),": ", varnbdim, "-", vardims(1:varnbdim)," l1d=" ,l1d
486	!---- read it
487	IF (l1d) THEN
488	CALL restget &
489	& (rest_id_sto, TRIM(varnames(iv)), 1, 1, &
490	& 1, itau_dep, .TRUE., xtmp)
491	CALL restput &
492	& (rest_id_sto, TRIM(varnames(iv)), 1, 1, &
493	& 1, itau_dep, xtmp)
494	ELSE
495	orch_vardims = varnbdim - 1 ! exclude time dimension introduced by IOIPSL
496	! Deal for different number of dimension
497	IF (orch_vardims == 2) THEN
498	!----
499	! vardims X(1), Y(2), remaining variables, time(last position)
500	ALLOCATE( var_2d(nbp_glo), stat=ier)
501	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_4d : error in memory allocation', '', '')
502	!----
503	CALL restget &
504	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, 1, &
505	& 1, itau_dep, .TRUE., var_2d, "gather", nbp_glo, indices_g)
506	CALL restput &
507	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, 1, &
508	& 1, itau_dep, var_2d, 'scatter', nbp_glo, indices_g)
509	!----
510	DEALLOCATE(var_2d)
511	ELSE IF (orch_vardims == 3) THEN
512	!----
513	! vardims X(1), Y(2), remaining variables, time(last position)
514	ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
515	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_3d : error in memory allocation', '', '')
516	!----
517	CALL restget &
518	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
519	& 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g)
520	CALL restput &
521	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
522	& 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g)
523	!----
524	DEALLOCATE(var_3d)
525	ELSE IF (orch_vardims == 4) THEN
526	!----
527	! vardims X(1), Y(2), remaining variables, time(last position)
528	ALLOCATE( var_4d(nbp_glo,vardims(3),vardims(4)), stat=ier)
529	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_4d : error in memory allocation', '', '')
530	!----
531	CALL restget &
532	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
533	& vardims(4), itau_dep, .TRUE., var_4d, "gather", nbp_glo, indices_g)
534	CALL restput &
535	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
536	& vardims(4), itau_dep, var_4d, 'scatter', nbp_glo, indices_g)
537	!----
538	DEALLOCATE(var_4d)
539	ELSE IF (orch_vardims == 5) THEN
540	!----
541	! vardims X(1), Y(2), remaining variables, time(last position)
542	ALLOCATE( var_5d(nbp_glo,vardims(3),vardims(4),vardims(5)), stat=ier)
543	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_5d : error in memory allocation', '', '')
544	!----
545	CALL restget &
546	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
547	& vardims(4), vardims(5), itau_dep, .TRUE., var_5d, "gather", nbp_glo, indices_g)
548	CALL restput &
549	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
550	& vardims(4), vardims(5), itau_dep, var_5d, 'scatter', nbp_glo, indices_g)
551	!----
552	DEALLOCATE(var_5d)
553	ELSE IF (orch_vardims == 6) THEN
554	!----
555	! vardims X(1), Y(2), remaining variables, time(last position)
556	ALLOCATE( var_6d(nbp_glo,vardims(3),vardims(4),vardims(5),vardims(6)), stat=ier)
557	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'var_6d : error in memory allocation', '', '')
558	!----
559	CALL restget &
560	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
561	& vardims(4), vardims(5), vardims(6), itau_dep, .TRUE., var_6d, "gather", nbp_glo, indices_g)
562	CALL restput &
563	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
564	& vardims(4), vardims(5), vardims(6), itau_dep, var_6d, 'scatter', nbp_glo, indices_g)
565	!----
566	DEALLOCATE(var_6d)
567	ELSE
568	WRITE( msg3, '(i5)' ) orch_vardims
569	CALL ipslerr(3, 'forcesoil', 'Varialbes(1) Restart Read/Write not implement for N dimensions(2)', &
570	& TRIM(varnames(iv)), TRIM(msg3))
571	ENDIF
572	ENDIF
573	ENDIF
574	ENDDO
575	! Length of the run (in Years)
576	! User value is used if present in the .def file
577	! If not, default value (10000 Years) is used
578	!-
579	WRITE(time_str,'(a)') '10000Y'
580	CALL getin('TIME_LENGTH', time_str)
581	write(numout,*) 'Number of years for carbon spinup : ',time_str
582	! transform into itau
583	CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
584	write(numout,*) 'Number of time steps to do: ',itau_len
585
586	! read soil carbon stocks values stored in the input restart file
587	!-
588	IF (.NOT. ok_pc) THEN
589	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm), stat=ier)
590	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_g : error in memory allocation', '', '')
591	carbon_g(:,:,:) = val_exp
592	CALL restget &
593	& (rest_id_sto, 'carbon', nbp_glo, ncarb , nvm, itau_dep, &
594	& .TRUE., carbon_g, 'gather', nbp_glo, indices_g)
595	IF (ALL(carbon_g == val_exp)) carbon_g = zero
596	WRITE(numout,*) "date0 : ",date0, itau_dep
597	ELSE
598	!-
599	! Permafrost carbon
600	!-
601	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm), stat=ier)
602	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_g : error in memory allocation', '', '')
603	carbon_g(:,:,:) = 0.
604	ALLOCATE(carbon_surf_g(nbp_glo,ncarb,nvm), stat=ier)
605	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'carbon_surf_g : error in memory allocation', '', '')
606	carbon_surf_g(:,:,:) = 0.
607	ALLOCATE(deepC_a_g(nbp_glo,ndeep,nvm), stat=ier)
608	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_a_g : error in memory allocation', '', '')
609	ALLOCATE(deepC_s_g(nbp_glo,ndeep,nvm), stat=ier)
610	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_s_g : error in memory allocation', '', '')
611	ALLOCATE(deepC_p_g(nbp_glo,ndeep,nvm), stat=ier)
612	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'deepC_p_g : error in memory allocation', '', '')
613	ALLOCATE(O2_soil_g(nbp_glo,ndeep,nvm), stat=ier)
614	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'O2_soil_g : error in memory allocation', '', '')
615	ALLOCATE(CH4_soil_g(nbp_glo,ndeep,nvm), stat=ier)
616	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'CH4_soil_g : error in memory allocation', '', '')
617	ALLOCATE(O2_snow_g(nbp_glo,nsnow,nvm), stat=ier)
618	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'O2_snow_g : error in memory allocation', '', '')
619	ALLOCATE(CH4_snow_g(nbp_glo,nsnow,nvm), stat=ier)
620	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'CH4_snow_g : error in memory allocation', '', '')
621	ALLOCATE(altmax_g(nbp_glo,nvm), stat=ier)
622	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'altmax_g : error in memory allocation', '', '')
623
624	deepC_a_g(:,:,:) = val_exp
625	CALL restget (rest_id_sto, 'deepC_a', nbp_glo, ndeep, nvm, itau_dep, &
626	& .TRUE., deepC_a_g, 'gather', nbp_glo, indices_g)
627	IF (ALL(deepC_a_g == val_exp)) deepC_a_g = zero
628
629	deepC_s_g(:,:,:) = val_exp
630	CALL restget (rest_id_sto, 'deepC_s', nbp_glo, ndeep, nvm, itau_dep, &
631	& .TRUE., deepC_s_g, 'gather', nbp_glo, indices_g)
632	IF (ALL(deepC_s_g == val_exp)) deepC_s_g = zero
633
634	deepC_p_g(:,:,:) = val_exp
635	CALL restget (rest_id_sto, 'deepC_p', nbp_glo, ndeep, nvm, itau_dep, &
636	& .TRUE., deepC_p_g, 'gather', nbp_glo, indices_g)
637	IF (ALL(deepC_p_g == val_exp)) deepC_p_g = zero
638
639	var_name= 'altmax'
640	altmax_g(:,:) = val_exp
641	CALL restget (rest_id_sto, var_name, nbp_glo, nvm, 1, itau_dep, .TRUE., altmax_g, "gather", nbp_glo, indices_g)
642	IF ( ALL( altmax_g(:,:) .EQ. val_exp ) ) THEN
643	CALL ipslerr(3, 'forcesoil', 'altmax is not found in stomate restart file', '', '')
644	END IF
645
646	CALL getin('reset_soilc', reset_soilc)
647	IF (reset_soilc) THEN
648	CALL ipslerr(1, 'forcesoil', 'deepC_a, deepC_s and deeC_p', &
649	'are ignored and set to zero value due to', 'reset_soilc option')
650	deepC_a_g(:,:,:) = zero
651	deepC_s_g(:,:,:) = zero
652	deepC_p_g(:,:,:) = zero
653	ENDIF
654
655	O2_soil_g(:,:,:) = val_exp
656	CALL restget (rest_id_sto, 'O2_soil', nbp_glo, ndeep, nvm, itau_dep, &
657	& .TRUE., O2_soil_g, 'gather', nbp_glo, indices_g)
658	IF (ALL(O2_soil_g == val_exp)) O2_soil_g = O2_init_conc
659
660	CH4_soil_g(:,:,:) = val_exp
661	CALL restget (rest_id_sto, 'CH4_soil', nbp_glo, ndeep, nvm, itau_dep, &
662	& .TRUE., CH4_soil_g, 'gather', nbp_glo, indices_g)
663	IF (ALL(CH4_soil_g == val_exp)) CH4_soil_g = CH4_init_conc
664
665	O2_snow_g(:,:,:) = val_exp
666	CALL restget (rest_id_sto, 'O2_snow', nbp_glo, nsnow, nvm, itau_dep, &
667	& .TRUE., O2_snow_g, 'gather', nbp_glo, indices_g)
668	IF (ALL(O2_snow_g == val_exp)) O2_snow_g = O2_init_conc
669
670	CH4_snow_g(:,:,:) = val_exp
671	CALL restget (rest_id_sto, 'CH4_snow', nbp_glo, nsnow, nvm, itau_dep, &
672	& .TRUE., CH4_snow_g, 'gather', nbp_glo, indices_g)
673	IF (ALL(CH4_snow_g == val_exp)) CH4_snow_g = CH4_init_conc
674	ENDIF ! ok_pc
675	ENDIF ! is_root_prc
676	!
677	CALL bcast(nbp_glo)
678	CALL bcast(iim_g)
679	CALL bcast(jjm_g)
680	IF (.NOT. ALLOCATED(indices_g)) ALLOCATE (indices_g(nbp_glo))
681	CALL bcast(indices_g)
682	CALL bcast(nparan)
683	CALL bcast(nbyear)
684	CALL bcast(dt_forcesoil)
685	CALL bcast(itau_dep)
686	CALL bcast(itau_len)
687	IF (ok_pc) THEN
688	CALL bcast(zz_deep)
689	CALL bcast(zz_coef_deep)
690	ENDIF
691	!
692	! we must initialize data_para :
693	CALL init_orchidee_data_para_driver(nbp_glo,indices_g)
694
695	kjpindex=nbp_loc
696	jjm=jj_nb
697	iim=iim_g
698	IF (printlev_loc>=3) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
699	!-
700	! Analytical spinup is set to false
701	!
702	spinup_analytic = .FALSE.
703	!-
704	! read soil carbon inputs, water and temperature stresses on OM
705	! decomposition
706	! into the forcing file - We read an average year.
707	!-
708	IF (.NOT. ok_pc) THEN
709	!-
710	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
711	! and allocate variables.
712	!-
713	#ifdef CPP_PARA
714	nc_opts = IOR(NF90_NOWRITE, NF90_MPIIO)
715	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name),IOR(nc_opts, NF90_NETCDF4),Cforcing_id, &
716	& comm = MPI_COMM_ORCH, info = MPI_INFO_NULL ))
717	#else
718	CALL nccheck( NF90_OPEN (TRIM(Cforcing_name), IOR(NF90_NOWRITE, NF90_NETCDF4),Cforcing_id) )
719	#endif
720	ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear), stat=ier)
721	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'control_temp : error in memory allocation', '', '')
722	ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear), stat=ier)
723	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'control_moist : error in memory allocation', '', '')
724	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear), stat=ier)
725	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'soilcarbon_input : error in memory allocation', '', '')
726	ALLOCATE(veget_max(kjpindex,nvm,nparan*nbyear), stat=ier)
727	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'veget_max : error in memory allocation', 'No ok_pc', '')
728	!
729	CALL nccheck( NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'dt_sechiba',dt_sechiba))
730	!-
731	start_4d = (/ nbp_mpi_para_begin(mpi_rank), 1, 1, 1 /)
732	count_4d = (/ nbp_mpi_para(mpi_rank), ncarb, nvm, nparan*nbyear /)
733	CALL nccheck( NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id))
734	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input, &
735	& start = start_4d, count = count_4d ))
736	!
737	start_3d = (/ nbp_mpi_para_begin(mpi_rank), 1, 1 /)
738	count_3d = (/ nbp_mpi_para(mpi_rank), nlevs, nparan*nbyear /)
739	CALL nccheck( NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id))
740	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,control_moist, &
741	& start = start_3d, count = count_3d))
742	CALL nccheck( NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id))
743	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,control_temp, &
744	& start = start_3d, count = count_3d))
745	count_3d = (/ nbp_mpi_para(mpi_rank), nvm, nparan*nbyear /)
746	CALL nccheck( NF90_INQ_VARID (Cforcing_id, 'veget_max',v_id))
747	CALL nccheck( NF90_GET_VAR (Cforcing_id,v_id,veget_max, &
748	& start = start_3d, count = count_3d))
749	!- Close Netcdf carbon permafrost file reference
750	CALL nccheck( NF90_CLOSE (Cforcing_id))
751	ELSE
752	!-
753	! Read permafrost-related soil carbon from Cforcing_id
754	!-
755	ALLOCATE(pb(kjpindex,nparan*nbyear), stat=ier)
756	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'pb : error in memory allocation', '', '')
757	ALLOCATE(snow(kjpindex,nparan*nbyear), stat=ier)
758	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snow : error in memory allocation', '', '')
759	ALLOCATE(tprof(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
760	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'tprof : error in memory allocation', '', '')
761	ALLOCATE(fbact(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
762	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'fbact : error in memory allocation', '', '')
763	ALLOCATE(hslong(kjpindex,ndeep,nvm,nparan*nbyear), stat=ier)
764	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'hslong : error in memory allocation', '', '')
765	ALLOCATE(veget_max(kjpindex,nvm,nparan*nbyear), stat=ier)
766	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'veget_max : error in memory allocation', '', '')
767	ALLOCATE(rprof(kjpindex,nvm,nparan*nbyear), stat=ier)
768	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'rprof : error in memory allocation', '', '')
769	ALLOCATE(tsurf(kjpindex,nparan*nbyear), stat=ier)
770	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'tsurf : error in memory allocation', '', '')
771	ALLOCATE(lalo(kjpindex,2), stat=ier)
772	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'lalo : error in memory allocation', '', '')
773	ALLOCATE(snowdz(kjpindex,nsnow,nparan*nbyear), stat=ier)
774	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snowdz_ : error in memory allocation', '', '')
775	ALLOCATE(snowrho(kjpindex,nsnow,nparan*nbyear), stat=ier)
776	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'snowrho : error in memory allocation', '', '')
777	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear), stat=ier)
778	IF (ier /= 0) CALL ipslerr(3, 'forcesoil', 'soilcarbon_input : error in memory allocation', '', '')
779	!-
780	CALL stomate_io_carbon_permafrost_read(Cforcing_name, nparan, nbyear,&
781	nbp_mpi_para_begin(mpi_rank), nbp_mpi_para(mpi_rank), &
782	soilcarbon_input, pb, snow, tsurf, &
783	tprof, fbact, hslong, rprof, &
784	lalo, snowdz, snowrho, veget_max )
785	ENDIF ! ok_pc
786	!---
787	!--- Create the index table
788	!---
789	!--- This job returns a LOCAL kindex.
790	!---
791	ALLOCATE (indices(kjpindex),stat=ier)
792	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'indices : error in memory allocation', '', '')
793	!
794	!! scattering to all processes in parallel mode
795	!
796	CALL scatter(indices_g,indices)
797	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
798	IF (printlev_loc>=3) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
799	!
800	! Initialize _index variables
801	CALL stomate_init_index(nbp_glo, kjpindex, indices)
802	!-
803	! Allocation of the variables for a processor
804	!-
805	ALLOCATE(clay(kjpindex), stat=ier)
806	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'clay : error in memory allocation', '', '')
807	ALLOCATE(carbon(kjpindex,ncarb,nvm), stat=ier)
808	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'indices : error in memory allocation', '', '')
809	!-
810	IF (.NOT. ok_pc) THEN
811	ALLOCATE(resp_hetero_soil(kjpindex,nvm), stat=ier)
812	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'resp_hetero_soil : error in memory allocation', '', '')
813	ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools), stat=ier)
814	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'matrixA : error in memory allocation', '', '')
815	DO i = 1,nbpools
816	matrixA(:,:,i,i) = un
817	ENDDO
818	ELSE
819	ALLOCATE(carbon_surf(kjpindex,ncarb,nvm), stat=ier)
820	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'carbon_surf : error in memory allocation', '', '')
821	ALLOCATE(deepC_a(kjpindex,ndeep,nvm), stat=ier)
822	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_a : error in memory allocation', '', '')
823	ALLOCATE(deepC_s(kjpindex,ndeep,nvm), stat=ier)
824	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_s : error in memory allocation', '', '')
825	ALLOCATE(deepC_p(kjpindex,ndeep,nvm), stat=ier)
826	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'deepC_p : error in memory allocation', '', '')
827	ALLOCATE(O2_soil(kjpindex,ndeep,nvm), stat=ier)
828	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'O2_soil : error in memory allocation', '', '')
829	ALLOCATE(CH4_soil(kjpindex,ndeep,nvm), stat=ier)
830	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'CH4_soil : error in memory allocation', '', '')
831	ALLOCATE(O2_snow(kjpindex,nsnow,nvm), stat=ier)
832	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'O2_snow : error in memory allocation', '', '')
833	ALLOCATE(CH4_snow(kjpindex,nsnow,nvm), stat=ier)
834	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'CH4_snow : error in memory allocation', '', '')
835	ALLOCATE(altmax(kjpindex,nvm), stat=ier)
836	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'altmax : error in memory allocation', '', '')
837	ALLOCATE(z_organic(kjpindex), stat=ier)
838	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'z_organic : error in memory allocation', '', '')
839	ENDIF
840	!-
841	! Initialization of the variables for a processor
842	!-
843	CALL Scatter(clay_g,clay)
844	CALL Scatter(carbon_g,carbon)
845	IF (ok_pc) THEN
846	CALL Scatter(z_organic_g,z_organic)
847	CALL Scatter(carbon_surf_g,carbon_surf)
848	CALL Scatter(deepC_a_g,deepC_a)
849	CALL Scatter(deepC_s_g,deepC_s)
850	CALL Scatter(deepC_p_g,deepC_p)
851	CALL Scatter(O2_soil_g,O2_soil)
852	CALL Scatter(CH4_soil_g,CH4_soil)
853	CALL Scatter(O2_snow_g,O2_snow)
854	CALL Scatter(CH4_snow_g,CH4_snow)
855	CALL Scatter(altmax_g,altmax)
856	ENDIF
857	!-
858	! Configuration of the parameters
859	!-
860	!Config Key = FRAC_CARB_AP
861	!Config Desc = frac carb coefficients from active pool: depends on clay
862	!content
863	!Config if = OK_STOMATE
864	!Config Def = 0.004
865	!Config Help = fraction of the active pool going to the passive pool
866	!Config Units = [-]
867	CALL getin_p('FRAC_CARB_AP',frac_carb_ap)
868	!
869	!Config Key = FRAC_CARB_SA
870	!Config Desc = frac_carb_coefficients from slow pool
871	!Config if = OK_STOMATE
872	!Config Def = 0.42
873	!Config Help = fraction of the slow pool going to the active pool
874	!Config Units = [-]
875	CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
876	!
877	!Config Key = FRAC_CARB_SP
878	!Config Desc = frac_carb_coefficients from slow pool
879	!Config if = OK_STOMATE
880	!Config Def = 0.03
881	!Config Help = fraction of the slow pool going to the passive pool
882	!Config Units = [-]
883	CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
884	!
885	!Config Key = FRAC_CARB_PA
886	!Config Desc = frac_carb_coefficients from passive pool
887	!Config if = OK_STOMATE
888	!Config Def = 0.45
889	!Config Help = fraction of the passive pool going to the passive pool
890	!Config Units = [-]
891	CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
892	!
893	!Config Key = FRAC_CARB_PS
894	!Config Desc = frac_carb_coefficients from passive pool
895	!Config if = OK_STOMATE
896	!Config Def = 0.0
897	!Config Help = fraction of the passive pool going to the passive pool
898	!Config Units = [-]
899	CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
900	!
901	!Config Key = ACTIVE_TO_PASS_CLAY_FRAC
902	!Config Desc =
903	!Config if = OK_STOMATE
904	!Config Def = .68
905	!Config Help =
906	!Config Units = [-]
907	CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
908	!
909	!Config Key = CARBON_TAU_IACTIVE
910	!Config Desc = residence times in carbon pools
911	!Config if = OK_STOMATE
912	!Config Def = 0.149
913	!Config Help =
914	!Config Units = [days]
915	CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
916	!
917	!Config Key = CARBON_TAU_ISLOW
918	!Config Desc = residence times in carbon pools
919	!Config if = OK_STOMATE
920	!Config Def = 5.48
921	!Config Help =
922	!Config Units = [days]
923	CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
924	!
925	!Config Key = CARBON_TAU_IPASSIVE
926	!Config Desc = residence times in carbon pools
927	!Config if = OK_STOMATE
928	!Config Def = 241.
929	!Config Help =
930	!Config Units = [days]
931	CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
932	!
933	!Config Key = FLUX_TOT_COEFF
934	!Config Desc =
935	!Config if = OK_STOMATE
936	!Config Def = 1.2, 1.4,.75
937	!Config Help =
938	!Config Units = [days]
939	CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)
940
941	!
942	!! 2. Computational step
943	!! Loop over time - Call of soilcarbon routine at each time step
944	!! Updated soil carbon stocks are stored into carbon variable
945	!! We only keep the last value of carbon variable (no time dimension).
946	!!-
947	IF (.NOT.ok_pc) THEN
948	iyear=1
949	iatt = 0
950	DO i=1,itau_len
951	iatt = iatt+1
952	IF (iatt > nparan*nbyear) THEN
953	IF (printlev>=3) WRITE(numout,*) iyear
954	iatt = 1
955	iyear=iyear+1
956	ENDIF
957	CALL soilcarbon &
958	& (kjpindex, dt_forcesoil, clay, &
959	& soilcarbon_input(:,:,:,iatt), &
960	& control_temp(:,:,iatt), control_moist(:,:,iatt), veget_max(:,:,iatt), &
961	& carbon, resp_hetero_soil, &
962	& matrixA)
963	ENDDO
964	WRITE(numout,*) "End of soilcarbon LOOP."
965	ELSE
966	IF ( satsoil ) hslong(:,:,:,:) = 1.
967	!these variables are only ouputs from deep_carbcycle (thus not necessary for
968	!Gather and Scatter)
969	ALLOCATE(heat_Zimov(kjpindex,ndeep,nvm), stat=ier)
970	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'heat_Zimov : error in memory allocation', '', '')
971	ALLOCATE(sfluxCH4_deep(kjpindex), stat=ier)
972	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'sfluxCH4_deep : error in memory allocation', '', '')
973	ALLOCATE(sfluxCO2_deep(kjpindex), stat=ier)
974	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'sfluxCO2_deep : error in memory allocation', '', '')
975	ALLOCATE(fixed_cryoturbation_depth(kjpindex,nvm), stat=ier)
976	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'fixed_cryoturbation_depth : error in memory allocation', '', '')
977	ALLOCATE(resp_hetero_soil(kjpindex, nvm), stat=ier)
978	IF (ier /= 0) CALL ipslerr_p(3, 'forcesoil', 'resp_hetero_soil : error in memory allocation', '', '')
979
980	iatt = 0
981	iyear=1
982	DO i=1,itau_len
983	iatt = iatt+1
984	IF (iatt > nparan*nbyear) THEN
985	IF (printlev>=3) WRITE(numout,*) iyear
986	iatt = 1
987	iyear=iyear+1
988	ENDIF
989	WRITE(numout, *) "Forcesoil:: deep_carbcycle, iyear=", iyear
990	CALL deep_carbcycle(kjpindex, indices, iatt, dt_forcesoil*one_day, lalo, clay, &
991	tsurf(:,iatt), tprof(:,:,:,iatt), hslong(:,:,:,iatt), snow(:,iatt), heat_Zimov, pb(:,iatt), &
992	sfluxCH4_deep, sfluxCO2_deep, &
993	deepC_a, deepC_s, deepC_p, O2_soil, CH4_soil, O2_snow, CH4_snow, &
994	zz_deep, zz_coef_deep, z_organic, soilcarbon_input(:,:,:,iatt), &
995	veget_max(:,:,iatt), rprof(:,:,iatt), altmax, carbon, carbon_surf, resp_hetero_soil, &
996	fbact(:,:,:,iatt), fixed_cryoturbation_depth, snowdz(:,:,iatt), snowrho(:,:,iatt))
997	ENDDO
998
999	ENDIF
1000	!!-
1001	!! 3. write new carbon stocks into the ouput restart file
1002	!!-
1003	CALL restput_p (rest_id_sto, 'carbon', nbp_glo, ncarb , nvm, itau_dep, &
1004	& carbon, 'scatter', nbp_glo, indices_g)
1005
1006	IF (ok_pc) THEN
1007	CALL restput_p (rest_id_sto, 'deepC_a', nbp_glo, ndeep, nvm, itau_dep, &
1008	& deepC_a, 'scatter', nbp_glo, indices_g)
1009	CALL restput_p (rest_id_sto, 'deepC_s', nbp_glo, ndeep, nvm, itau_dep, &
1010	& deepC_s, 'scatter', nbp_glo, indices_g)
1011	CALL restput_p (rest_id_sto, 'deepC_p', nbp_glo, ndeep, nvm, itau_dep, &
1012	& deepC_p, 'scatter', nbp_glo, indices_g)
1013	CALL restput_p (rest_id_sto, 'O2_soil', nbp_glo, ndeep, nvm, itau_dep, &
1014	& O2_soil, 'scatter', nbp_glo, indices_g)
1015	CALL restput_p (rest_id_sto, 'CH4_soil', nbp_glo, ndeep, nvm, itau_dep, &
1016	& CH4_soil, 'scatter', nbp_glo, indices_g)
1017	CALL restput_p (rest_id_sto, 'O2_snow', nbp_glo, nsnow, nvm, itau_dep, &
1018	& O2_snow, 'scatter', nbp_glo, indices_g)
1019	CALL restput_p (rest_id_sto, 'CH4_snow', nbp_glo, nsnow, nvm, itau_dep, &
1020	& CH4_snow, 'scatter', nbp_glo, indices_g)
1021	CALL restput_p (rest_id_sto, 'altmax', nbp_glo, nvm, 1, itau_dep, &
1022	& altmax, 'scatter', nbp_glo, indices_g)
1023	ENDIF
1024	!-
1025	IF (is_root_prc) THEN
1026	!- Close restart files
1027	CALL getin_dump
1028	CALL restclo
1029	ENDIF
1030	!-
1031	#ifdef CPP_PARA
1032	CALL MPI_FINALIZE(ier)
1033	#endif
1034	WRITE(numout,*) "End of forcesoil."
1035	!--------------------
1036	END PROGRAM forcesoil

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-OP-r6850/src_driver/forcesoil.f90

Download in other formats: