Context Navigation

forcesoil.f90 @ 8730

Last change on this file since 8730 was 1451, checked in by sebastiaan.luyssaert, 11 years ago
DEV: trying to close the mass balance. Runs only for 12 years
File size: 27.4 KB

Line
1	! =================================================================================================================================
2	! PROGRAM : forcesoil
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	!>\BRIEF This subroutine runs the soilcarbon submodel using specific initial conditions
9	!! and driving variables in order to obtain soil carbon stocks closed to the steady-state values
10	!! quicker than when using the ''full'' ORCHIDEE.
11	!!
12	!!\n DESCRIPTION: None
13	!! This subroutine computes the soil carbon stocks by calling the soilcarbon routine at each time step. \n
14	!! The aim is to obtain soil carbon stocks closed to the steady-state values and ultimately to create
15	!! an updated stomate restart file for the stomate component. The state variables of the subsystem are the clay content
16	!! (fixed value) and the soil carbon stocks. Initial conditions for the state variables are read in an
17	!! input stomate restart file. Driving variables are Soil carbon input, Water and Temperature stresses on
18	!! Organic Matter decomposition. Driving variables are read from a specific forcing file produced by a former run of ORCHIDEE
19	!! (SECHIBA+STOMATE). \n
20	!! The FORCESOIL program first consists in reading a set of input files, allocating variables and
21	!! preparing output stomate restart file. \n
22	!! Then, a loop over time is performed in which the soilcarbon routine is called at each time step. \n
23	!! Last, final values of the soil carbon stocks are written into the output stomate restart file. \n
24	!! No flag is associated with the use of the FORCESOIL program. \n
25	!!
26	!! RECENT CHANGE(S): None
27	!!
28	!! REFERENCE(S) : None
29	!!
30	!! FLOWCHART : None
31	!!
32	!! SVN :
33	!! $HeadURL: $
34	!! $Date: $
35	!! $Revision: $
36	!! \n
37	!_ =================================================================================================================================
38
39	PROGRAM forcesoil
40
41	USE netcdf
42	!-
43	USE defprec
44	USE constantes
45	USE constantes_mtc
46	USE pft_parameters
47	USE stomate_data
48	USE ioipsl_para
49	USE mod_orchidee_para
50	USE stomate_soilcarbon
51
52	! USE parallel
53	!-
54	IMPLICIT NONE
55	!-
56	CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out
57	INTEGER(i_std) :: iim,jjm !! Indices (unitless)
58
59	INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless)
60	INTEGER(i_std) :: kjpindex !! Domain size (unitless)
61
62	INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?)
63	!! and number of time steps of the simulation (unitless)
64	CHARACTER(LEN=30) :: time_str !! Length of the simulation (year)
65	REAL(r_std) :: dt_files !! time step between two successive itaus (?)
66	!! (requested by restini routine) (seconds)
67	REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?)
68	INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless)
69	CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file
70	!! (unitless)
71	!-
72	CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless)
73	INTEGER :: Cforcing_id !! ID of the forcing file (unitless)
74	INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file
75	!! (unitless)
76	REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days)
77	INTEGER :: nparan !! Number of values stored per year in the forcing file
78	!! (unitless)
79	INTEGER :: nbyear
80	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless)
81	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless)
82	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless)
83	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined
84	!! in lat/lon (2D) (degrees)
85	REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless)
86
87
88	INTEGER :: i,m,iatt,iv,iyear !! counters (unitless)
89
90	CHARACTER(LEN=80) :: var_name
91	CHARACTER(LEN=800) :: taboo_vars !! string used for storing the name of the variables
92	!! of the stomate restart file that are not automatically
93	!! duplicated from input to output restart file (unitless)
94	REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless)
95	INTEGER(i_std),PARAMETER :: nbvarmax=300 !! maximum # of variables assumed in the stomate restart file
96	!! (unitless)
97	INTEGER(i_std) :: nbvar !! # of variables effectively present
98	!! in the stomate restart file (unitless)
99	CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored
100	!! in the stomate restart file (unitless)
101	INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable
102	!! of the stomate restart file
103	INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable
104	!! of the stomate restart file
105	INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable
106	!! of the stomate restart file
107	LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable
108	!! of the stomate restart file are of length 1 (ie scalar)
109	!! (unitless)
110	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless)
111	REAL(r_std) :: x_tmp !! temporary variable used to store return value
112	!! from nf90_get_att (unitless)
113	CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT
114	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless)
115	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g !! soil carbon input (nbpglob,ncarb,nvm,time)
116	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
117	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g !! Temperature control (nbp_glo,above/below,time) on OM decomposition
118	!! (unitless)
119	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g !! Moisture control (nbp_glo,abo/below,time) on OM decomposition
120	!! ?? Should be defined per PFT as well (unitless)
121	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$)
122
123	REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless)
124	REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) !! soil carbon input (nbp_loc,ncarb,nvm,time)
125	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
126	REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) !! Temperature control (nbp_loc,above/below,time) on OM decomposition
127	!! (unitless)
128	REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) !! Moisture control (nbp_loc,abo/below,time) on OM decomposition
129	!! ?? Should be defined per PFT as well (unitless)
130	REAL(r_std),ALLOCATABLE :: carbon(:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$)
131	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
132	!! (requested by soilcarbon routine but not used here)
133
134	INTEGER(i_std) :: ier,iret !! Used for hangling errors
135
136	CHARACTER(LEN=30) :: temp_name
137	LOGICAL :: debug !! boolean used for printing messages
138	LOGICAL :: l_error !! boolean for memory allocation
139	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: matrixA !! Carbon fluxes matrix*
140	REAL(r_std), DIMENSION(:,:),ALLOCATABLE :: dummy_veget_max !! the subroutin soil carbon needs veget_max to check for the
141	!! closure of the C_balance. This is of interest in this
142	!! routine. A dummy variable is passed
143
144	!_ =================================================================================================================================
145
146	CALL Init_orchidee_para
147	CALL init_timer
148
149	!-
150	! Configure the number of PFTS
151	!-
152
153	! 1. Read the number of PFTs
154	!
155	!Config Key = NVM
156	!Config Desc = number of PFTs
157	!Config If = OK_SECHIBA or OK_STOMATE
158	!Config Def = 13
159	!Config Help = The number of vegetation types define by the user
160	!Config Units = [-]
161	CALL getin_p('NVM',nvm)
162
163	! 2. Allocation
164	l_error = .FALSE.
165	ALLOCATE(pft_to_mtc(nvm),stat=ier)
166	l_error = l_error .OR. (ier .NE. 0)
167	IF (l_error) THEN
168	STOP 'pft_to_mtc (forcesoil only) : error in memory allocation'
169	ENDIF
170
171	! 3. Initialisation of the correspondance table
172	pft_to_mtc(:) = undef_int
173
174	! 4.Reading of the conrrespondance table in the .def file
175	!
176	!Config Key = PFT_TO_MTC
177	!Config Desc = correspondance array linking a PFT to MTC
178	!Config if = OK_SECHIBA or OK_STOMATE
179	!Config Def = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13
180	!Config Help =
181	!Config Units = [-]
182	CALL getin_p('PFT_TO_MTC',pft_to_mtc)
183
184	! 4.1 if nothing is found, we use the standard configuration
185	IF(nvm == nvmc ) THEN
186	IF(pft_to_mtc(1) == undef_int) THEN
187	WRITE(numout,*) 'Note to the user : we will use ORCHIDEE to its standard configuration'
188	pft_to_mtc(:) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 /)
189	ENDIF
190	ELSE
191	IF(pft_to_mtc(1) == undef_int) THEN
192	WRITE(numout,*)' The array PFT_TO_MTC is empty : we stop'
193	ENDIF
194	ENDIF
195
196	! 4.2 What happened if pft_to_mtc(j) > nvmc (if the mtc doesn't exist)?
197	DO i = 1, nvm
198	IF(pft_to_mtc(i) > nvmc) THEN
199	WRITE(numout,*) "the MTC you chose doesn't exist"
200	STOP 'we stop reading pft_to_mtc'
201	ENDIF
202	ENDDO
203
204	! 4.3 Check if pft_to_mtc(1) = 1
205	IF(pft_to_mtc(1) /= 1) THEN
206	WRITE(numout,*) 'the first pft has to be the bare soil'
207	STOP 'we stop reading next values of pft_to_mtc'
208	ELSE
209	DO i = 2,nvm
210	IF(pft_to_mtc(i) == 1) THEN
211	WRITE(numout,*) 'only pft_to_mtc(1) has to be the bare soil'
212	STOP 'we stop reading pft_to_mtc'
213	ENDIF
214	ENDDO
215	ENDIF
216
217	! 5. Allocate and initialize natural ans is_c4
218
219	! 5.1 Memory allocation
220	l_error = .FALSE.
221	ALLOCATE(natural(nvm),stat=ier)
222	l_error = l_error .OR. (ier .NE. 0)
223	ALLOCATE(is_c4(nvm),stat=ier)
224
225	IF (l_error) THEN
226	STOP 'natural or is_c4 (forcesoil only) : error in memory allocation'
227	ENDIF
228
229	! 5.2 Initialisation
230	DO i = 1, nvm
231	natural(i) = natural_mtc(pft_to_mtc(i))
232	is_c4(i) = is_c4_mtc(pft_to_mtc(i))
233	ENDDO
234
235	!-
236	!-
237	! set debug to have more information
238	!-
239	!Config Key = DEBUG_INFO
240	!Config Desc = Flag for debug information
241	!Config If =
242	!Config Def = n
243	!Config Help = This option allows to switch on the output of debug
244	!Config information without recompiling the code.
245	!Config Units = [FLAG]
246	!-
247	debug = .FALSE.
248	CALL getin_p('DEBUG_INFO',debug)
249	!!-
250	!! 1. Initialisation stage
251	!! Reading a set of input files, allocating variables and preparing output restart file.
252	!!-
253	! Define restart file name
254	! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var).
255	! User values are used if present in the .def file.
256	! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used.
257	!-
258	IF (is_root_prc) THEN
259	sto_restname_in = 'stomate_start.nc'
260	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
261	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
262	sto_restname_out = 'stomate_rest_out.nc'
263	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
264	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
265	!-
266	! Open the input file and Get some Dimension and Attributes ID's
267	!-
268	iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto)
269	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g)
270	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g)
271	iret = NF90_INQ_VARID (rest_id_sto, "time", iv)
272	iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar)
273	iret = NF90_CLOSE (rest_id_sto)
274	i=INDEX(thecalendar,ACHAR(0))
275	IF ( i > 0 ) THEN
276	thecalendar(i:20)=' '
277	ENDIF
278	!-
279	! Allocate longitudes and latitudes
280	!-
281	ALLOCATE (lon(iim_g,jjm_g))
282	ALLOCATE (lat(iim_g,jjm_g))
283	lon(:,:) = zero
284	lat(:,:) = zero
285	lev(1) = zero
286	!-
287	CALL restini &
288	& (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
289	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
290	ENDIF
291
292	CALL bcast(date0)
293	CALL bcast(thecalendar)
294	WRITE(numout,*) "calendar = ",thecalendar
295	!-
296	! calendar
297	!-
298	CALL ioconf_calendar (thecalendar)
299	CALL ioget_calendar (one_year,one_day)
300	CALL ioconf_startdate(date0)
301	!
302	!! For master process only
303	!
304	IF (is_root_prc) THEN
305	!-
306	! define forcing file's name (Cforcing_name var)
307	! User value is used if present in the .def file
308	! If not, default (NONE) is used
309	!-
310	Cforcing_name = 'NONE'
311	CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
312	!-
313	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
314	! and allocate variables.
315	!-
316	iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
317	IF (iret /= NF90_NOERR) THEN
318	CALL ipslerr (3,'forcesoil', &
319	& 'Could not open file : ', &
320	& Cforcing_name,'(Do you have forget it ?)')
321	ENDIF
322	!-
323	! Total Domain size is stored in nbp_glo variable
324	!-
325	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
326	nbp_glo = NINT(x_tmp)
327	!-
328	! Number of values stored per year in the forcing file is stored in nparan var.
329	!-
330	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
331	nparan = NINT(x_tmp)
332	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nbyear',x_tmp)
333	nbyear = NINT(x_tmp)
334	!-
335	ALLOCATE (indices_g(nbp_glo))
336	ALLOCATE (clay_g(nbp_glo))
337	!-
338	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
339	ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
340	ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g)
341	indices_g(:) = NINT(x_indices_g(:))
342	WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
343	DEALLOCATE (x_indices_g)
344	!-
345	ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
346	ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g)
347	!-
348	! time step of forcesoil program (in days)
349	!-
350	dt_forcesoil = one_year / FLOAT(nparan)
351	WRITE(numout,*) 'time step (d): ',dt_forcesoil
352	WRITE(numout,*) 'nparan: ',nparan
353	WRITE(numout,*) 'nbyear: ',nbyear
354	!-
355	! read and write the variables in the output restart file we do not modify within the Forcesoil program
356	! ie all variables stored in the input restart file except those stored in taboo_vars
357	!-
358	taboo_vars ='$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
359	& '$day_counter$ $dt_days$ $date$ '
360	!-
361	DO m = 1,nvm
362	WRITE(part_str,'(I2)') m
363	IF (m < 10) part_str(1:1) = '0'
364	temp_name = '$carbon_'//part_str(1:LEN_TRIM(part_str))//'$'
365	taboo_vars = TRIM(taboo_vars)//' '//TRIM(temp_name)
366	ENDDO
367	!-
368	CALL ioget_vname(rest_id_sto, nbvar, varnames)
369	!-
370	! read and write some special variables (1D or variables that we need)
371	!-
372	var_name = 'day_counter'
373	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
374	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
375	!-
376	var_name = 'dt_days'
377	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
378	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
379	!-
380	var_name = 'date'
381	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
382	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
383	!-
384	DO iv=1,nbvar
385	!-- check if the variable is to be written here
386	IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
387	!---- get variable dimensions, especially 3rd dimension
388	CALL ioget_vdim &
389	& (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
390	l1d = ALL(vardims(1:varnbdim) == 1)
391	!---- read it
392	IF (l1d) THEN
393	CALL restget &
394	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
395	& 1, itau_dep, .TRUE., xtmp)
396	ELSE
397	ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
398	IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
399	!----
400	CALL restget &
401	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
402	& 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g)
403	ENDIF
404	!---- write it
405	IF (l1d) THEN
406	CALL restput &
407	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
408	& 1, itau_dep, xtmp)
409	ELSE
410	CALL restput &
411	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
412	& 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g)
413	!----
414	DEALLOCATE(var_3d)
415	ENDIF
416	ENDIF
417	ENDDO
418	!-
419	! read soil carbon stocks values stored in the input restart file
420	!-
421	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm))
422	carbon_g(:,:,:) = val_exp
423	DO m = 1, nvm
424	WRITE (part_str, '(I2)') m
425	IF (m<10) part_str(1:1)='0'
426	var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
427	CALL restget &
428	& (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
429	& .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g)
430	IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero
431	!-- do not write this variable: it will be modified.
432	ENDDO
433	WRITE(numout,*) "date0 : ",date0, itau_dep
434	!-
435	! Analytical spinup is set to false
436	!
437	spinup_analytic = .FALSE.
438
439	! Length of the run (in Years)
440	! User value is used if present in the .def file
441	! If not, default value (10000 Years) is used
442	!-
443	WRITE(time_str,'(a)') '10000Y'
444	CALL getin('TIME_LENGTH', time_str)
445	write(numout,*) 'Number of years for carbon spinup : ',time_str
446	! transform into itau
447	CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
448	write(numout,*) 'Number of time steps to do: ',itau_len
449	!-
450	! read soil carbon inputs, water and temperature stresses on OM decomposition
451	! into the forcing file - We read an average year.
452	!-
453	ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan*nbyear))
454	ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan*nbyear))
455	ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan*nbyear))
456	!-
457	ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
458	ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g)
459	ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id)
460	ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_g)
461	ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id)
462	ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_g)
463	!-
464	ier = NF90_CLOSE (Cforcing_id)
465	!-
466	ENDIF
467	CALL bcast(nparan)
468	CALL bcast(nbyear)
469	CALL bcast(dt_forcesoil)
470	CALL bcast(iim_g)
471	CALL bcast(jjm_g)
472	CALL bcast(nbp_glo)
473	CALL bcast(itau_dep)
474	CALL bcast(itau_len)
475	!
476	! We must initialize data_para :
477	!
478	!
479	CALL init_orchidee_data_para_driver(nbp_glo,indices_g)
480
481	kjpindex=nbp_loc
482	jjm=jj_nb
483	iim=iim_g
484	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
485
486	!---
487	!--- Create the index table
488	!---
489	!--- This job returns a LOCAL kindex.
490	!---
491	ALLOCATE (indices(kjpindex),stat=ier)
492	!
493	!! scattering to all processes in parallel mode
494	!
495	CALL scatter(indices_g,indices)
496	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
497	IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
498	!-
499	! Allocation of the variables for a processor
500	!-
501	ALLOCATE(clay(kjpindex))
502	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan*nbyear))
503	ALLOCATE(control_temp(kjpindex,nlevs,nparan*nbyear))
504	ALLOCATE(control_moist(kjpindex,nlevs,nparan*nbyear))
505	ALLOCATE(carbon(kjpindex,ncarb,nvm))
506	ALLOCATE(resp_hetero_soil(kjpindex,nvm))
507	ALLOCATE(matrixA(kjpindex,nvm,nbpools,nbpools))
508	DO i = 1,nbpools
509	matrixA(:,:,i,i) = un
510	ENDDO
511	iatt = 0
512
513	!-
514	! Initialization of the variables for a processor
515	!-
516	CALL Scatter(clay_g,clay)
517	CALL Scatter(soilcarbon_input_g,soilcarbon_input)
518	CALL Scatter(control_temp_g,control_temp)
519	CALL Scatter(control_moist_g,control_moist)
520	CALL Scatter(carbon_g,carbon)
521
522	!-
523	! Configuration of the parameters
524	!-
525
526	!Config Key = FRAC_CARB_AP
527	!Config Desc = frac carb coefficients from active pool: depends on clay content
528	!Config if = OK_STOMATE
529	!Config Def = 0.004
530	!Config Help = fraction of the active pool going to the passive pool
531	!Config Units = [-]
532	CALL getin_p('FRAC_CARB_AP',frac_carb_ap)
533	!
534	!Config Key = FRAC_CARB_SA
535	!Config Desc = frac_carb_coefficients from slow pool
536	!Config if = OK_STOMATE
537	!Config Def = 0.42
538	!Config Help = fraction of the slow pool going to the active pool
539	!Config Units = [-]
540	CALL getin_p('FRAC_CARB_SA',frac_carb_sa)
541	!
542	!Config Key = FRAC_CARB_SP
543	!Config Desc = frac_carb_coefficients from slow pool
544	!Config if = OK_STOMATE
545	!Config Def = 0.03
546	!Config Help = fraction of the slow pool going to the passive pool
547	!Config Units = [-]
548	CALL getin_p('FRAC_CARB_SP',frac_carb_sp)
549	!
550	!Config Key = FRAC_CARB_PA
551	!Config Desc = frac_carb_coefficients from passive pool
552	!Config if = OK_STOMATE
553	!Config Def = 0.45
554	!Config Help = fraction of the passive pool going to the passive pool
555	!Config Units = [-]
556	CALL getin_p('FRAC_CARB_PA',frac_carb_pa)
557	!
558	!Config Key = FRAC_CARB_PS
559	!Config Desc = frac_carb_coefficients from passive pool
560	!Config if = OK_STOMATE
561	!Config Def = 0.0
562	!Config Help = fraction of the passive pool going to the passive pool
563	!Config Units = [-]
564	CALL getin_p('FRAC_CARB_PS',frac_carb_ps)
565	!
566	!Config Key = ACTIVE_TO_PASS_CLAY_FRAC
567	!Config Desc =
568	!Config if = OK_STOMATE
569	!Config Def = .68
570	!Config Help =
571	!Config Units = [-]
572	CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac)
573	!
574	!Config Key = CARBON_TAU_IACTIVE
575	!Config Desc = residence times in carbon pools
576	!Config if = OK_STOMATE
577	!Config Def = 0.149
578	!Config Help =
579	!Config Units = [days]
580	CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive)
581	!
582	!Config Key = CARBON_TAU_ISLOW
583	!Config Desc = residence times in carbon pools
584	!Config if = OK_STOMATE
585	!Config Def = 5.48
586	!Config Help =
587	!Config Units = [days]
588	CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow)
589	!
590	!Config Key = CARBON_TAU_IPASSIVE
591	!Config Desc = residence times in carbon pools
592	!Config if = OK_STOMATE
593	!Config Def = 241.
594	!Config Help =
595	!Config Units = [days]
596	CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive)
597	!
598	!Config Key = FLUX_TOT_COEFF
599	!Config Desc =
600	!Config if = OK_STOMATE
601	!Config Def = 1.2, 1.4,.75
602	!Config Help =
603	!Config Units = [days]
604	CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff)
605
606	!!-
607	!! 2. Computational step
608	!! Loop over time - Call of soilcarbon routine at each time step
609	!! Updated soil carbon stocks are stored into carbon variable
610	!! We only keep the last value of carbon variable (no time dimension).
611	!!-
612	iyear=1
613
614	! Initialize the dummy variable
615	ALLOCATE(dummy_veget_max(kjpindex,nvm))
616	dummy_veget_max(:,:) = un
617
618	DO i=1,itau_len
619	iatt = iatt+1
620	IF (iatt > nparan*nbyear) THEN
621	IF (debug) WRITE(numout,*) iyear
622	iatt = 1
623	iyear=iyear+1
624	ENDIF
625	CALL soilcarbon &
626	& (kjpindex, dt_forcesoil, clay, &
627	& soilcarbon_input(:,:,:,iatt), &
628	& control_temp(:,:,iatt), control_moist(:,:,iatt), &
629	& carbon, resp_hetero_soil, &
630	& matrixA, dummy_veget_max)
631	ENDDO
632	WRITE(numout,*) "End of soilcarbon LOOP."
633
634	!
635	!! Gathering of variables towards main processor in parallel mode
636	!
637	CALL Gather(carbon,carbon_g)
638	!!-
639	!! 3. write new carbon stocks into the ouput restart file
640	!!-
641	IF (is_root_prc) THEN
642	DO m=1,nvm
643	WRITE (part_str, '(I2)') m
644	IF (m<10) part_str(1:1)='0'
645	var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
646	CALL restput &
647	& (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
648	& carbon_g(:,:,m), 'scatter', nbp_glo, indices_g)
649	ENDDO
650	!-
651	CALL getin_dump
652	CALL restclo
653	ENDIF
654	#ifdef CPP_PARA
655	CALL MPI_FINALIZE(ier)
656	#endif
657	WRITE(numout,*) "End of forcesoil."
658	!--------------------
659	END PROGRAM forcesoil

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_driver/forcesoil.f90 @ 8730

Download in other formats: