Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

forcesoil.f90 @ 855

Last change on this file since 855 was 687, checked in by fabienne.maignan, 12 years ago
Editing done
File size: 21.7 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_veg
46	USE constantes_co2
47	USE stomate_constants
48	USE ioipsl
49	USE stomate_soilcarbon
50	USE parallel
51	!-
52	IMPLICIT NONE
53	!-
54	CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out
55	INTEGER(i_std) :: iim,jjm !! Indices (unitless)
56
57	INTEGER(i_std),PARAMETER :: llm = 1 !! Vertical Layers (requested by restini routine) (unitless)
58	INTEGER(i_std) :: kjpindex !! Domain size (unitless)
59
60	INTEGER(i_std) :: itau_dep,itau_len !! Time step read in the restart file (?)
61	!! and number of time steps of the simulation (unitless)
62	CHARACTER(LEN=30) :: time_str !! Length of the simulation (year)
63	REAL(r_std) :: dt_files !! time step between two successive itaus (?)
64	!! (requested by restini routine) (seconds)
65	REAL(r_std) :: date0 !! Time at which itau = 0 (requested by restini routine) (?)
66	INTEGER(i_std) :: rest_id_sto !! ID of the input restart file (unitless)
67	CHARACTER(LEN=20), SAVE :: thecalendar = 'noleap' !! Type of calendar defined in the input restart file
68	!! (unitless)
69	!-
70	CHARACTER(LEN=100) :: Cforcing_name !! Name of the forcing file (unitless)
71	INTEGER :: Cforcing_id !! ID of the forcing file (unitless)
72	INTEGER :: v_id !! ID of the variable 'Index' stored in the forcing file
73	!! (unitless)
74	REAL(r_std) :: dt_forcesoil !! Time step at which soilcarbon routine is called (days)
75	INTEGER :: nparan !! Number of values stored per year in the forcing file
76	!! (unitless)
77
78	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! Grid Point Index used per processor (unitless)
79	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Grid Point Index for all processor (unitless)
80	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Grid Point Index for all processor (unitless)
81	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat !! Longitude and Latitude of each grid point defined
82	!! in lat/lon (2D) (degrees)
83	REAL(r_std),DIMENSION(llm) :: lev !! Number of level (requested by restini routine) (unitless)
84
85
86	INTEGER :: i,m,iatt,iv,iyear !! counters (unitless)
87
88	CHARACTER(LEN=80) :: var_name !! temporary variable to store string (unitless)
89	CHARACTER(LEN=400) :: taboo_vars !! string used for storing the name of the variables
90	!! of the stomate restart file that are not automatically
91	!! duplicated from input to output restart file (unitless)
92	REAL(r_std),DIMENSION(1) :: xtmp !! scalar read/written in restget/restput routines (unitless)
93	INTEGER(i_std),PARAMETER :: nbvarmax=300 !! maximum # of variables assumed in the stomate restart file
94	!! (unitless)
95	INTEGER(i_std) :: nbvar !! # of variables effectively present
96	!! in the stomate restart file (unitless)
97	CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames !! list of the names of the variables stored
98	!! in the stomate restart file (unitless)
99	INTEGER(i_std) :: varnbdim !! # of dimensions of a given variable
100	!! of the stomate restart file
101	INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! maximal # of dimensions assumed for any variable
102	!! of the stomate restart file
103	INTEGER,DIMENSION(varnbdim_max) :: vardims !! length of each dimension of a given variable
104	!! of the stomate restart file
105	LOGICAL :: l1d !! boolean : TRUE if all dimensions of a given variable
106	!! of the stomate restart file are of length 1 (ie scalar)
107	!! (unitless)
108	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d !! matrix read/written in restget/restput routines (unitless)
109	REAL(r_std) :: x_tmp !! temporary variable used to store return value
110	!! from nf90_get_att (unitless)
111	CHARACTER(LEN=10) :: part_str !! string suffix indicating the index of a PFT
112	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g !! clay fraction (nbpglo) (unitless)
113	REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g !! soil carbon input (nbpglob,ncarb,nvm,time)
114	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
115	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g !! Temperature control (nbp_glo,above/below,time) on OM decomposition
116	!! (unitless)
117	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g !! Moisture control (nbp_glo,abo/below,time) on OM decomposition
118	!! ?? Should be defined per PFT as well (unitless)
119	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g !! Soil carbon stocks (nbp_glo,ncarb,nvm) (\f$gC m^{-2}\f$)
120
121	REAL(r_std),ALLOCATABLE :: clay(:) !! clay fraction (nbp_loc) (unitless)
122	REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) !! soil carbon input (nbp_loc,ncarb,nvm,time)
123	!! (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
124	REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) !! Temperature control (nbp_loc,above/below,time) on OM decomposition
125	!! (unitless)
126	REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) !! Moisture control (nbp_loc,abo/below,time) on OM decomposition
127	!! ?? Should be defined per PFT as well (unitless)
128	REAL(r_std),ALLOCATABLE :: carbon(:,:,:) !! Soil carbon stocks (nbp_loc,ncarb,nvm) (\f$gC m^{-2}\f$)
129	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil !! Heterotrophic respiration (\f$gC m^{-2} dt_forcesoil^{-1}\f$)
130	!! (requested by soilcarbon routine but not used here)
131
132	INTEGER(i_std) :: ier,iret !! Used for hangling errors
133
134	LOGICAL :: debug !! boolean used for printing messages
135
136	CALL Init_para(.FALSE.)
137	CALL init_timer
138
139	!---------------------------------------------------------------------
140	!-
141	! set debug to have more information
142	!-
143	!Config Key = DEBUG_INFO
144	!Config Desc = Flag for debug information
145	!Config Def = n
146	!Config Help = This option allows to switch on the output of debug
147	!Config information without recompiling the code.
148	!-
149	debug = .FALSE.
150	CALL getin_p('DEBUG_INFO',debug)
151	!!-
152	!! 1. Initialisation stage
153	!! Reading a set of input files, allocating variables and preparing output restart file.
154	!!-
155	! Define restart file name
156	! for reading initial conditions (sto_restname_in var) and for writting final conditions (sto_restname_out var).
157	! User values are used if present in the .def file.
158	! If not present, default values (stomate_start.nc and stomate_rest_out.c) are used.
159	!-
160	IF (is_root_prc) THEN
161	sto_restname_in = 'stomate_start.nc'
162	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
163	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
164	sto_restname_out = 'stomate_rest_out.nc'
165	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
166	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
167	!-
168	! Open the input file and Get some Dimension and Attributes ID's
169	!-
170	iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto)
171	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g)
172	iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g)
173	iret = NF90_INQ_VARID (rest_id_sto, "time", iv)
174	iret = NF90_GET_ATT (rest_id_sto, iv, 'calendar',thecalendar)
175	iret = NF90_CLOSE (rest_id_sto)
176	i=INDEX(thecalendar,ACHAR(0))
177	IF ( i > 0 ) THEN
178	thecalendar(i:20)=' '
179	ENDIF
180	!-
181	! Allocate longitudes and latitudes
182	!-
183	ALLOCATE (lon(iim_g,jjm_g))
184	ALLOCATE (lat(iim_g,jjm_g))
185	lon(:,:) = 0.0
186	lat(:,:) = 0.0
187	lev(1) = 0.0
188	!-
189	CALL restini &
190	& (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, &
191	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
192	ENDIF
193
194	CALL bcast(date0)
195	CALL bcast(thecalendar)
196	WRITE(numout,*) "calendar = ",thecalendar
197	!-
198	! calendar
199	!-
200	CALL ioconf_calendar (thecalendar)
201	CALL ioget_calendar (one_year,one_day)
202	CALL ioconf_startdate(date0)
203	!
204	!! For master process only
205	!
206	IF (is_root_prc) THEN
207	!-
208	! define forcing file's name (Cforcing_name var)
209	! User value is used if present in the .def file
210	! If not, default (NONE) is used
211	!-
212	Cforcing_name = 'NONE'
213	CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name)
214	!-
215	! Open FORCESOIL's forcing file to read some basic info (dimensions, variable ID's)
216	! and allocate variables.
217	!-
218	iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id)
219	IF (iret /= NF90_NOERR) THEN
220	CALL ipslerr (3,'forcesoil', &
221	& 'Could not open file : ', &
222	& Cforcing_name,'(Do you have forget it ?)')
223	ENDIF
224	!-
225	! Total Domain size is stored in nbp_glo variable
226	!-
227	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp)
228	nbp_glo = NINT(x_tmp)
229	!-
230	! Number of values stored per year in the forcing file is stored in nparan var.
231	!-
232	ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp)
233	nparan = NINT(x_tmp)
234	!-
235	ALLOCATE (indices_g(nbp_glo))
236	ALLOCATE (clay_g(nbp_glo))
237	!-
238	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
239	ier = NF90_INQ_VARID (Cforcing_id,'index',v_id)
240	ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g)
241	indices_g(:) = NINT(x_indices_g(:))
242	WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g
243	DEALLOCATE (x_indices_g)
244	!-
245	ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id)
246	ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g)
247	!-
248	! time step of forcesoil program (in days)
249	!-
250	dt_forcesoil = one_year / FLOAT(nparan)
251	WRITE(numout,*) 'time step (d): ',dt_forcesoil
252	!-
253	! read and write the variables in the output restart file we do not modify within the Forcesoil program
254	! ie all variables stored in the input restart file except those stored in taboo_vars
255	!-
256	taboo_vars = '$lon$ $lat$ $lev$ $nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$ '// &
257	& '$day_counter$ $dt_days$ $date$ '// &
258	& '$carbon_01$ $carbon_02$ $carbon_03$ $carbon_04$ $carbon_05$'// &
259	& '$carbon_06$ $carbon_07$ $carbon_08$ $carbon_09$ $carbon_10$'// &
260	& '$carbon_11$ $carbon_12$ $carbon_13$'
261	!-
262	CALL ioget_vname(rest_id_sto, nbvar, varnames)
263	!-
264	! read and write some special variables (1D or variables that we need)
265	!-
266	var_name = 'day_counter'
267	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
268	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
269	!-
270	var_name = 'dt_days'
271	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
272	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
273	!-
274	var_name = 'date'
275	CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
276	CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
277	!-
278	DO iv=1,nbvar
279	!-- check if the variable is to be written here
280	IF (INDEX(taboo_vars,'$'//TRIM(varnames(iv))//'$') == 0 ) THEN
281	!---- get variable dimensions, especially 3rd dimension
282	CALL ioget_vdim &
283	& (rest_id_sto, varnames(iv), varnbdim_max, varnbdim, vardims)
284	l1d = ALL(vardims(1:varnbdim) == 1)
285	!---- read it
286	IF (l1d) THEN
287	CALL restget &
288	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
289	& 1, itau_dep, .TRUE., xtmp)
290	ELSE
291	ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier)
292	IF (ier /= 0) STOP 'ALLOCATION PROBLEM'
293	!----
294	CALL restget &
295	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
296	& 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g)
297	ENDIF
298	!---- write it
299	IF (l1d) THEN
300	CALL restput &
301	& (rest_id_sto, TRIM(varnames(iv)), 1, vardims(3), &
302	& 1, itau_dep, xtmp)
303	ELSE
304	CALL restput &
305	& (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), &
306	& 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g)
307	!----
308	DEALLOCATE(var_3d)
309	ENDIF
310	ENDIF
311	ENDDO
312	!-
313	! read soil carbon stocks values stored in the input restart file
314	!-
315	ALLOCATE(carbon_g(nbp_glo,ncarb,nvm))
316	carbon_g(:,:,:) = val_exp
317	DO m = 1, nvm
318	WRITE (part_str, '(I2)') m
319	IF (m<10) part_str(1:1)='0'
320	var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
321	CALL restget &
322	& (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
323	& .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g)
324	IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero
325	!-- do not write this variable: it will be modified.
326	ENDDO
327	WRITE(numout,*) "date0 : ",date0, itau_dep
328	!-
329	! Length of the run (in Years)
330	! User value is used if present in the .def file
331	! If not, default value (10000 Years) is used
332	!-
333	WRITE(time_str,'(a)') '10000Y'
334	CALL getin('TIME_LENGTH', time_str)
335	write(numout,*) 'Number of years for carbon spinup : ',time_str
336	! transform into itau
337	CALL tlen2itau(time_str, dt_forcesoil*one_day, date0, itau_len)
338	write(numout,*) 'Number of time steps to do: ',itau_len
339	!-
340	! read soil carbon inputs, water and temperature stresses on OM decomposition
341	! into the forcing file - We read an average year.
342	!?! It should not be an average year
343	!?! ORCHIDEE is capable to produce multi-annual Cforcing files
344	!?! time dimension for soilcarbon_input_g, control_moist_g, control_temp_g
345	!?! should not be nparan but nparan*nb_year
346	!!?? I completely agree.
347	!-
348	ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan))
349	ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan))
350	ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan))
351	!-
352	ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id)
353	ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g)
354	ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id)
355	ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_g)
356	ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id)
357	ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_g)
358	!-
359	ier = NF90_CLOSE (Cforcing_id)
360	!-
361	ENDIF
362	CALL bcast(nparan)
363	CALL bcast(dt_forcesoil)
364	CALL bcast(iim_g)
365	CALL bcast(jjm_g)
366	call bcast(nbp_glo)
367	CALL bcast(itau_dep)
368	CALL bcast(itau_len)
369	!
370	! We must initialize data_para :
371	!
372	!
373	CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
374
375	kjpindex=nbp_loc
376	jjm=jj_nb
377	iim=iim_g
378	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
379
380	!---
381	!--- Create the index table
382	!---
383	!--- This job returns a LOCAL kindex.
384	!---
385	ALLOCATE (indices(kjpindex),stat=ier)
386	!
387	!! scattering to all processes in parallel mode
388	!
389	CALL scatter(indices_g,indices)
390	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
391	IF (debug) WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex)
392	!-
393	! Allocation of the variables for a processor
394	!-
395	ALLOCATE(clay(kjpindex))
396	!?! It should not be an average year
397	!?! ORCHIDEE is capable to produce multi-annual Cforcing files
398	!?! time dimension for soilcarbon_input, control_moist, control_temp
399	!?! should not be nparan but nparan*nb_year
400	ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan))
401	ALLOCATE(control_temp(kjpindex,nlevs,nparan))
402	ALLOCATE(control_moist(kjpindex,nlevs,nparan))
403	ALLOCATE(carbon(kjpindex,ncarb,nvm))
404	ALLOCATE(resp_hetero_soil(kjpindex,nvm))
405	iatt = 0
406
407	!-
408	! Initialization of the variables for a processor
409	!-
410	CALL Scatter(clay_g,clay)
411	CALL Scatter(soilcarbon_input_g,soilcarbon_input)
412	CALL Scatter(control_temp_g,control_temp)
413	CALL Scatter(control_moist_g,control_moist)
414	CALL Scatter(carbon_g,carbon)
415
416	!!-
417	!! 2. Computational step
418	!! Loop over time - Call of soilcarbon routine at each time step
419	!! Updated soil carbon stocks are stored into carbon variable
420	!! We only keep the last value of carbon variable (no time dimension).
421	!!-
422	iyear=1
423	DO i=1,itau_len
424	iatt = iatt+1
425	IF (iatt > nparan) THEN
426	IF (debug) WRITE(numout,*) iyear
427	iatt = 1
428	iyear=iyear+1
429	ENDIF
430	!?! Time index should not be iatt
431	!?! but iatt+(iyear-1)*nparan or index should be i
432	CALL soilcarbon &
433	& (kjpindex, dt_forcesoil, clay, &
434	& soilcarbon_input(:,:,:,iatt), &
435	& control_temp(:,:,iatt), control_moist(:,:,iatt), &
436	& carbon, resp_hetero_soil)
437	ENDDO
438	WRITE(numout,*) "End of soilcarbon LOOP."
439
440	!
441	!! Gathering of variables towards main processor in parallel mode
442	!
443	CALL Gather(carbon,carbon_g)
444	!!-
445	!! 3. write new carbon stocks into the ouput restart file
446	!!-
447	IF (is_root_prc) THEN
448	DO m=1,nvm
449	WRITE (part_str, '(I2)') m
450	IF (m<10) part_str(1:1)='0'
451	var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str))
452	CALL restput &
453	& (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, &
454	& carbon_g(:,:,m), 'scatter', nbp_glo, indices_g)
455	ENDDO
456	!-
457	CALL getin_dump
458	CALL restclo
459	ENDIF
460	#ifdef CPP_PARA
461	CALL MPI_FINALIZE(ier)
462	#endif
463	WRITE(numout,*) "End of forcesoil."
464	!--------------------
465	END PROGRAM forcesoil

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

Context Navigation

source: perso/abdelouhab.djerrah/ORCHIDEE_OL/forcesoil.f90 @ 855

Download in other formats: