Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE_OL/teststomate.f90 @ 313

Last change on this file since 313 was 313, checked in by didier.solyga, 13 years ago
Update teststomate.f90 and forcesoil.f90 in order to take into account the modifications of the revisions 311 and 312
File size: 35.7 KB

Line
1	PROGRAM teststomate
2	!---------------------------------------------------------------------
3	! This program reads a forcing file for STOMATE
4	! which was created with STOMATE coupled to SECHIBA.
5	! It then integrates STOMATE for a
6	! given number of years using this forcing.
7	! The GPP is scaled to the new calculated LAI...
8	! Nothing is done for soil moisture which should
9	! actually evolve with the vegetation.
10	!---------------------------------------------------------------------
11	!- IPSL (2006)
12	!- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13	USE netcdf
14	!-
15	USE defprec
16	USE constantes
17	USE pft_parameters
18	USE stomate_data
19	USE ioipsl
20	USE grid
21	USE slowproc
22	USE stomate
23	USE intersurf, ONLY: stom_define_history , intsurf_time
24	USE parallel
25	!-
26	IMPLICIT NONE
27	!-
28	! Declarations
29	!-
30	INTEGER(i_std) :: vegax_id
31	INTEGER(i_std) :: kjpij,kjpindex
32	REAL(r_std) :: dtradia,dt_force
33	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices
34	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indexveg
35	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype, veget_x
36	REAL(r_std),DIMENSION(:),ALLOCATABLE :: totfrac_nobio
37	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: frac_nobio
38	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_x
39	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_x
40	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_force_x
41	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_force_x
42	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_force_x
43	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon,lat
44	REAL(r_std),DIMENSION(:),ALLOCATABLE :: t2m,t2m_min,temp_sol
45	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum
46	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: humrel_x
47	REAL(r_std),DIMENSION(:),ALLOCATABLE :: litterhum
48	REAL(r_std),DIMENSION(:),ALLOCATABLE :: precip_rain,precip_snow
49	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: gpp_x
50	REAL(r_std),DIMENSION(:),ALLOCATABLE :: deadleaf_cover
51	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: assim_param_x
52	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: height_x
53	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x
54	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux
55	!!$ >> DS add for externalised version
56	REAL(r_std),DIMENSION(:),ALLOCATABLE :: hist_PFTaxis
57	!!$ >> DS
58	!!$ >> Add for correct 1.9.5.1 version 01/07/2011
59	TYPE(control_type), SAVE :: control_flags
60	REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu
61	!!$>> DS
62
63
64	INTEGER :: ier,iret
65	INTEGER :: ncfid
66	LOGICAL :: a_er
67	CHARACTER(LEN=80) :: &
68	& dri_restname_in,dri_restname_out, &
69	& sec_restname_in,sec_restname_out, &
70	& sto_restname_in,sto_restname_out, stom_histname
71	INTEGER(i_std) :: iim,jjm
72	INTEGER(i_std),PARAMETER :: llm = 1
73	REAL(r_std),DIMENSION(llm) :: lev
74	REAL(r_std) :: dt_files
75	INTEGER(i_std) :: itau_dep,itau,itau_len,itau_step
76	REAL(r_std) :: date0
77	INTEGER(i_std) :: rest_id_sec,rest_id_sto
78	INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_sto,hist_id_stom_IPCC
79	CHARACTER(LEN=30) :: time_str
80	REAL :: hist_days_stom,hist_dt_stom
81	REAL(r_std),DIMENSION(10) :: hist_pool_10axis
82	REAL(r_std),DIMENSION(100) :: hist_pool_100axis
83	REAL(r_std),DIMENSION(11) :: hist_pool_11axis
84	REAL(r_std),DIMENSION(101) :: hist_pool_101axis
85	INTEGER :: hist_PFTaxis_id,hori_id
86	INTEGER :: hist_pool_10axis_id
87	INTEGER :: hist_pool_100axis_id
88	INTEGER :: hist_pool_11axis_id
89	INTEGER :: hist_pool_101axis_id
90	INTEGER :: hist_level
91	INTEGER,PARAMETER :: max_hist_level = 10
92	INTEGER(i_std) :: i,j,iv,id
93	CHARACTER*80 :: var_name
94	CHARACTER(LEN=40),DIMENSION(10) :: fluxop
95	LOGICAL :: ldrestart_read,ldrestart_write
96	LOGICAL :: l1d
97	INTEGER(i_std),PARAMETER :: nbvarmax=200
98	INTEGER(i_std) :: nbvar
99	CHARACTER*50,DIMENSION(nbvarmax) :: varnames
100	INTEGER(i_std) :: varnbdim
101	INTEGER(i_std),PARAMETER :: varnbdim_max=20
102	INTEGER,DIMENSION(varnbdim_max) :: vardims
103	CHARACTER*200 :: taboo_vars
104	REAL(r_std),DIMENSION(1) :: xtmp
105	INTEGER :: nsfm,nsft
106	INTEGER :: iisf
107	INTEGER(i_std) :: max_totsize,totsize_1step
108	INTEGER(i_std),DIMENSION(0:2) :: ifirst,ilast
109	INTEGER(i_std) :: iblocks,nblocks
110	INTEGER,PARAMETER :: ndm = 10
111	INTEGER,DIMENSION(ndm) :: start,count
112	INTEGER :: ndim,v_id
113	INTEGER :: force_id
114	CHARACTER(LEN=100) :: forcing_name
115	REAL :: x
116	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices
117	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours
118	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d
119	!-
120	REAL(r_std) :: time_sec,time_step_sec
121	REAL(r_std) :: time_dri,time_step_dri
122	REAL(r_std),DIMENSION(1) :: r1d
123	REAL(r_std) :: julian,djulian
124	! REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype
125	!-
126	! the following variables contain the forcing data
127	!-
128	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: clay_fm
129	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: humrel_x_fm
130	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: litterhum_fm
131	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_fm
132	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_min_fm
133	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: temp_sol_fm
134	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soiltemp_fm
135	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soilhum_fm
136	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: precip_fm
137	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: gpp_x_fm
138	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_force_x_fm
139	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_max_force_x_fm
140	REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lai_force_x_fm
141	INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: isf
142	LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:) :: nf_written
143	INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nf_cumul
144	INTEGER(i_std) :: ji,jv
145	!---------------------------------------------------------------------
146	!- No parallelisation yet in teststomate !
147	#ifdef CPP_PARA
148	CALL ipslerr (3,'teststomate', &
149	& 'You try to run testsomate compiled with parallelisation. (CPP_PARA key)', &
150	& 'But it wasn''t programmed yet and teststomate will stop.','You must compiled it without CPP_PARA key.')
151	#endif
152
153	!-
154	! calendar
155	!-
156	CALL ioconf_calendar ('noleap')
157	CALL ioget_calendar (one_year,one_day)
158	!-
159	! open STOMATE's forcing file to read some basic info
160	!-
161	forcing_name = 'stomate_forcing.nc'
162	CALL getin ('STOMATE_FORCING_NAME',forcing_name)
163	iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,force_id)
164	IF (iret /= NF90_NOERR) THEN
165	CALL ipslerr (3,'teststomate', &
166	& 'Could not open file : ', &
167	& forcing_name,'(Do you have forget it ?)')
168	ENDIF
169	ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dtradia',dtradia)
170	ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dt_slow',dt_force)
171	ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'nsft',x)
172	nsft = NINT(x)
173	ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpij',x)
174	kjpij = NINT(x)
175	ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpindex',x)
176	kjpindex = NINT(x)
177	CALL init_grid( kjpindex )
178	!-
179	write(,) 'ATTENTION',dtradia,dt_force
180
181	!!$ >> DS added for the externalised version
182	!!$ we need to know the pft parameters values used for stomate
183
184	! 1. Read the number of PFTs
185	CALL getin('NVM',nvm)
186	! 2. Allocate and read the pft parameters stomate specific
187	CALL pft_parameters_main
188	CALL getin_stomate_pft_parameters
189	!!$ DS
190
191	!-
192	! allocate arrays
193	!-
194	a_er = .FALSE.
195	!!$ >> DS added for the externalised version
196	ALLOCATE (hist_PFTaxis(nvm),stat=ier)
197	a_er = a_er .OR. (ier.NE.0)
198	! 01/07/2011 allocate fco2_lu
199	ALLOCATE (fco2_lu(kjpindex),stat=ier)
200	a_er = a_er .OR. (ier.NE.0)
201	!!$ end add DS
202	ALLOCATE (indices(kjpindex),stat=ier)
203	a_er = a_er .OR. (ier.NE.0)
204	ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
205	a_er = a_er .OR. (ier.NE.0)
206	ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
207	a_er = a_er .OR. (ier.NE.0)
208	ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
209	a_er = a_er .OR. (ier.NE.0)
210	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
211	a_er = a_er .OR. (ier.NE.0)
212	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
213	a_er = a_er .OR. (ier.NE.0)
214	ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
215	a_er = a_er .OR. (ier.NE.0)
216	ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
217	a_er = a_er .OR. (ier.NE.0)
218	ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
219	a_er = a_er .OR. (ier.NE.0)
220	ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
221	a_er = a_er .OR. (ier.NE.0)
222	ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
223	a_er = a_er .OR. (ier.NE.0)
224	ALLOCATE (t2m(kjpindex),stat=ier)
225	a_er = a_er .OR. (ier.NE.0)
226	ALLOCATE (t2m_min(kjpindex),stat=ier)
227	a_er = a_er .OR. (ier.NE.0)
228	ALLOCATE (temp_sol(kjpindex),stat=ier)
229	a_er = a_er .OR. (ier.NE.0)
230	ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
231	a_er = a_er .OR. (ier.NE.0)
232	ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
233	a_er = a_er .OR. (ier.NE.0)
234	ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
235	a_er = a_er .OR. (ier.NE.0)
236	ALLOCATE (litterhum(kjpindex),stat=ier)
237	a_er = a_er .OR. (ier.NE.0)
238	ALLOCATE (precip_rain(kjpindex),stat=ier)
239	a_er = a_er .OR. (ier.NE.0)
240	ALLOCATE (precip_snow(kjpindex),stat=ier)
241	a_er = a_er .OR. (ier.NE.0)
242	ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
243	a_er = a_er .OR. (ier.NE.0)
244	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
245	a_er = a_er .OR. (ier.NE.0)
246	ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
247	a_er = a_er .OR. (ier.NE.0)
248	ALLOCATE (height_x(kjpindex,nvm),stat=ier)
249	a_er = a_er .OR. (ier.NE.0)
250	ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
251	a_er = a_er .OR. (ier.NE.0)
252	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
253	a_er = a_er .OR. (ier.NE.0)
254	IF ( a_er ) STOP 'PROBLEM WITH ALLOCATION'
255	!
256	! prepare forcing
257	!
258	max_totsize = 50
259	CALL getin ('STOMATE_FORCING_MEMSIZE',max_totsize)
260	max_totsize = max_totsize * 1000000
261	totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
262	SIZE(humrel_x)*KIND(humrel_x) + &
263	SIZE(litterhum)*KIND(litterhum) + &
264	SIZE(t2m)*KIND(t2m) + &
265	SIZE(t2m_min)*KIND(t2m_min) + &
266	SIZE(temp_sol)*KIND(temp_sol) + &
267	SIZE(soiltemp)*KIND(soiltemp) + &
268	SIZE(soilhum)*KIND(soilhum) + &
269	SIZE(precip_rain)*KIND(precip_rain) + &
270	SIZE(precip_snow)*KIND(precip_snow) + &
271	SIZE(gpp_x)*KIND(gpp_x) + &
272	SIZE(veget_force_x)*KIND(veget_force_x) + &
273	SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
274	SIZE(lai_force_x)*KIND(lai_force_x)
275	! total number of forcing steps
276	nsft = INT(one_year/(dt_force/one_day))
277	! number of forcing steps in memory
278	nsfm = MIN(nsft,MAX(1,NINT(FLOAT(max_totsize)/FLOAT(totsize_1step))))
279	!-
280	WRITE(numout,*) 'Offline forcing of Stomate:'
281	WRITE(numout,*) ' Total number of forcing steps:',nsft
282	WRITE(numout,*) ' Number of forcing steps in memory:',nsfm
283	!-
284	a_er = .FALSE.
285	ALLOCATE (clay_fm(kjpindex,nsfm),stat=ier)
286	a_er = a_er.OR.(ier.NE.0)
287	ALLOCATE (humrel_x_fm(kjpindex,nvm,nsfm),stat=ier)
288	a_er = a_er.OR.(ier.NE.0)
289	ALLOCATE (litterhum_fm(kjpindex,nsfm),stat=ier)
290	a_er = a_er.OR.(ier.NE.0)
291	ALLOCATE (t2m_fm(kjpindex,nsfm),stat=ier)
292	a_er = a_er.OR.(ier.NE.0)
293	ALLOCATE (t2m_min_fm(kjpindex,nsfm),stat=ier)
294	a_er = a_er.OR.(ier.NE.0)
295	ALLOCATE (temp_sol_fm(kjpindex,nsfm),stat=ier)
296	a_er = a_er.OR.(ier.NE.0)
297	ALLOCATE (soiltemp_fm(kjpindex,nbdl,nsfm),stat=ier)
298	a_er = a_er.OR.(ier.NE.0)
299	ALLOCATE (soilhum_fm(kjpindex,nbdl,nsfm),stat=ier)
300	a_er = a_er.OR.(ier.NE.0)
301	ALLOCATE (precip_fm(kjpindex,nsfm),stat=ier)
302	a_er = a_er.OR.(ier.NE.0)
303	ALLOCATE (gpp_x_fm(kjpindex,nvm,nsfm),stat=ier)
304	a_er = a_er.OR.(ier.NE.0)
305	ALLOCATE (veget_force_x_fm(kjpindex,nvm,nsfm),stat=ier)
306	a_er = a_er.OR.(ier.NE.0)
307	ALLOCATE (veget_max_force_x_fm(kjpindex,nvm,nsfm),stat=ier)
308	a_er = a_er.OR.(ier.NE.0)
309	ALLOCATE (lai_force_x_fm(kjpindex,nvm,nsfm),stat=ier)
310	a_er = a_er.OR.(ier.NE.0)
311	ALLOCATE (isf(nsfm),stat=ier)
312	a_er = a_er.OR.(ier.NE.0)
313	ALLOCATE (nf_written(nsft),stat=ier)
314	a_er = a_er.OR.(ier.NE.0)
315	ALLOCATE (nf_cumul(nsft),stat=ier)
316	a_er = a_er.OR.(ier.NE.0)
317	IF (a_er) THEN
318	STOP 'stomate: error in memory allocation for forcing data'
319	ENDIF
320	! ensure that we read all new forcing states
321	iisf = nsfm
322	! initialize that table that contains the indices
323	! of the forcing states that will be in memory
324	isf(:) = (/ (i,i=1,nsfm) /)
325	!-
326	! read info about grids
327	!-
328	contfrac(:) = 1.0
329	!-
330	ALLOCATE (x_indices(kjpindex),stat=ier)
331	ier = NF90_INQ_VARID (force_id,'index',v_id)
332	ier = NF90_GET_VAR (force_id,v_id,x_indices)
333	indices(:) = NINT(x_indices(:))
334	DEALLOCATE (x_indices)
335	!-
336	DO ji=1,kjpindex
337	DO jv=1,nvm
338	indexveg((jv-1)kjpindex+ji) = indices(ji)+(jv-1)kjpij
339	ENDDO
340	ENDDO
341	!-
342	ier = NF90_INQ_VARID (force_id,'lalo',v_id)
343	ier = NF90_GET_VAR (force_id,v_id,lalo)
344	!-
345	ALLOCATE (x_neighbours(kjpindex,8),stat=ier)
346	ier = NF90_INQ_VARID (force_id,'neighbours',v_id)
347	ier = NF90_GET_VAR (force_id,v_id,x_neighbours)
348	neighbours(:,:) = NINT(x_neighbours(:,:))
349	DEALLOCATE (x_neighbours)
350	!-
351	ier = NF90_INQ_VARID (force_id,'resolution',v_id)
352	ier = NF90_GET_VAR (force_id,v_id,resolution)
353	!-
354	ier = NF90_INQ_VARID (force_id,'contfrac',v_id)
355	ier = NF90_GET_VAR (force_id,v_id,contfrac)
356	!-
357	! activate CO2, STOMATE, but not sechiba
358	!-
359	control%river_routing = .FALSE.
360	control%hydrol_cwrr = .FALSE.
361	control%ok_sechiba = .FALSE.
362	!
363	control%stomate_watchout = .TRUE.
364	control%ok_co2 = .TRUE.
365	control%ok_stomate = .TRUE.
366	!!$ >> DS now we search the values for the scalar parameters in some .def files
367	! 07/2011
368	CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate)
369	CALL veget_config
370	!
371	CALL getin_co2_parameters
372	CALL getin_stomate_parameters
373	!!$ >> DS
374	!-
375	! is DGVM activated?
376	!-
377	control%ok_dgvm = .FALSE.
378	CALL getin('STOMATE_OK_DGVM',control%ok_dgvm)
379	WRITE(,) 'LPJ is activated: ',control%ok_dgvm
380	!!$ >> DS for the externalisation reading the scalar parameters when ok_dgvm is set to true
381	IF (control%ok_dgvm) THEN
382	CALL getin_dgvm_parameters
383	ENDIF
384	!!$ >> DS
385	!-
386	! restart files
387	!-
388	! Sechiba's restart files
389	sec_restname_in = 'sechiba_start.nc'
390	CALL getin('SECHIBA_restart_in',sec_restname_in)
391	WRITE(,) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
392	IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
393	STOP 'Need a restart file for Sechiba'
394	ENDIF
395	sec_restname_out = 'sechiba_restart.nc'
396	CALL getin('SECHIBA_rest_out',sec_restname_out)
397	WRITE(,) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
398	! Stomate's restart files
399	sto_restname_in = 'stomate_start.nc'
400	CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in)
401	WRITE(,) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
402	sto_restname_out = 'stomate_restart.nc'
403	CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out)
404	WRITE(,) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
405	!-
406	! We need to know iim, jjm.
407	! Get them from the restart files themselves.
408	!-
409	iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
410	IF (iret /= NF90_NOERR) THEN
411	CALL ipslerr (3,'teststomate', &
412	& 'Could not open file : ', &
413	& sec_restname_in,'(Do you have forget it ?)')
414	ENDIF
415	iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim)
416	iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm)
417	iret = NF90_CLOSE (ncfid)
418	! Allocate longitudes and latitudes
419	ALLOCATE (lon(iim,jjm),stat=ier)
420	a_er = a_er.OR.(ier.NE.0)
421	ALLOCATE (lat(iim,jjm),stat=ier)
422	a_er = a_er.OR.(ier.NE.0)
423	lon(:,:) = 0.0
424	lat(:,:) = 0.0
425	lev(1) = 0.0
426	!-
427	CALL restini &
428	& (sec_restname_in, iim, jjm, lon, lat, llm, lev, &
429	& sec_restname_out, itau_dep, date0, dt_files, rest_id_sec)
430	!-
431	CALL restini &
432	& (sto_restname_in, iim, jjm, lon, lat, llm, lev, &
433	& sto_restname_out, itau_dep, date0, dt_files, rest_id_sto)
434	!-
435	IF ( dt_files .NE. dtradia ) THEN
436	WRITE(,) 'dt_files',dt_files
437	WRITE(,) 'dtradia',dtradia
438	STOP 'PROBLEM with time steps.'
439	ENDIF
440	!-
441	! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
442	!-
443	itau_step = NINT(dt_force/dtradia)
444	!
445	CALL ioconf_startdate(date0)
446	CALL intsurf_time( itau_dep, date0, dtradia )
447	!-
448	! Length of integration
449	!-
450	WRITE(time_str,'(a)') '1Y'
451	CALL getin ('TIME_LENGTH', time_str)
452	! transform into itau
453	CALL tlen2itau(time_str, dt_files, date0, itau_len)
454	! itau_len*dtradia must be a multiple of dt_force
455	itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
456	& *dt_force/dtradia)
457	!-
458	! set up STOMATE history file
459	!-
460	!Config Key = STOMATE_OUTPUT_FILE
461	!Config Desc = Name of file in which STOMATE's output is going
462	!Config to be written
463	!Config Def = stomate_history.nc
464	!Config Help = This file is going to be created by the model
465	!Config and will contain the output from the model.
466	!Config This file is a truly COADS compliant netCDF file.
467	!Config It will be generated by the hist software from
468	!Config the IOIPSL package.
469	!-
470	stom_histname='stomate_history.nc'
471	CALL getin ('STOMATE_OUTPUT_FILE', stom_histname)
472	WRITE(,) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
473	!-
474	!Config Key = STOMATE_HIST_DT
475	!Config Desc = STOMATE history time step (d)
476	!Config Def = 10.
477	!Config Help = Time step of the STOMATE history file
478	!-
479	hist_days_stom = 10.
480	CALL getin ('STOMATE_HIST_DT', hist_days_stom)
481	IF ( hist_days_stom == -1. ) THEN
482	hist_dt_stom = -1.
483	WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
484	ELSE
485	hist_dt_stom = NINT( hist_days_stom ) * one_day
486	WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
487	hist_dt_stom/one_day
488	ENDIF
489	!-
490	! initialize
491	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
492	& itau_dep, date0, dt_files, hori_id, hist_id_sto)
493	! define PFT axis
494	hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
495	! declare this axis
496	CALL histvert (hist_id_sto, 'PFT', 'Plant functional type', &
497	& '-', nvm, hist_PFTaxis, hist_PFTaxis_id)
498	!- define Pool_10 axis
499	hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
500	!- declare this axis
501	CALL histvert (hist_id_sto, 'P10', 'Pool 10 years', &
502	& '-', 10, hist_pool_10axis, hist_pool_10axis_id)
503	!- define Pool_100 axis
504	hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
505	!- declare this axis
506	CALL histvert (hist_id_sto, 'P100', 'Pool 100 years', &
507	& '-', 100, hist_pool_100axis, hist_pool_100axis_id)
508	!- define Pool_11 axis
509	hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
510	!- declare this axis
511	CALL histvert (hist_id_sto, 'P11', 'Pool 10 years + 1', &
512	& '-', 11, hist_pool_11axis, hist_pool_11axis_id)
513	!- define Pool_101 axis
514	hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
515	!- declare this axis
516	CALL histvert (hist_id_sto, 'P101', 'Pool 100 years + 1', &
517	& '-', 101, hist_pool_101axis, hist_pool_101axis_id)
518	! define STOMATE history file
519	CALL stom_define_history (hist_id_sto, nvm, iim, jjm, &
520	& dt_files, hist_dt_stom, hori_id, hist_PFTaxis_id, &
521	& hist_pool_10axis_id, hist_pool_100axis_id, &
522	& hist_pool_11axis_id, hist_pool_101axis_id,control_flags)
523	! end definition
524	CALL histend(hist_id_sto)
525	hist_id_sec = -1
526	hist_id_sec2 = -1
527	hist_id_stom_IPCC = -1
528	!-
529	! read some variables we need from SECHIBA's restart file
530	!-
531	CALL ioget_expval(val_exp)
532	!-
533	! first call of slowproc to initialize variables
534	!-
535	itau = itau_dep
536	ldrestart_read = .TRUE.
537	ldrestart_write = .FALSE.
538	!-
539	!MM Problem here with dpu which depends on soil type
540	DO jv = 1, nbdl-1
541	! first 2.0 is dpu
542	! second 2.0 is average
543	diaglev(jv) = 2.0/(2*(nbdl-1) -1) ( ( 2(jv-1) -1) + ( 2(jv) -1) ) / 2.0
544	ENDDO
545	diaglev(nbdl) = 2.0
546	!-
547	! For sequential use only, we must initialize data_para :
548	!
549	CALL init_para(.FALSE.)
550	!
551	CALL init_data_para(iim,jjm,kjpindex,indices)
552	!
553	!- global index index_g is the index_l of root proc
554	IF (is_root_prc) index_g(:)=indices(1:kjpindex)
555	!-
556
557	! >> DS : getin the value of slowproc parameters
558	CALL getin('CLAYFRACTION_DEFAULT',clayfraction_default)
559	CALL getin('SOILTYPE_DEFAULT',soiltype_default)
560	! >> DS
561
562	CALL slowproc_main &
563	& (itau, kjpij, kjpindex, dt_force, date0, &
564	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
565	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
566	& t2m, t2m_min, temp_sol, soiltemp, &
567	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
568	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
569	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
570	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, &
571	!!>> DS add 01/07/2011 correction due to some modifications in the trunk
572	& fco2_lu)
573	! correct date
574	day_counter = one_day - dt_force
575	date=1
576	!-
577	! time loop
578	!-
579	DO itau = itau_dep+itau_step,itau_dep+itau_len,itau_step
580	!-- next forcing state
581	iisf = iisf+1
582	!---
583	IF (iisf .GT. nsfm) THEN
584	!-----
585	!---- we have to read new forcing states from the file
586	!-----
587	!---- determine blocks of forcing states that are contiguous in memory
588	!-----
589	nblocks = 0
590	ifirst(:) = 1
591	ilast(:) = 1
592	!-----
593	DO iisf=1,nsfm
594	IF ( (nblocks .NE. 0) ) THEN
595	IF ( (isf(iisf) .EQ. isf(ilast(nblocks))+1) ) THEN
596	!-------- element is contiguous with last element found
597	ilast(nblocks) = iisf
598	ELSE
599	!-------- found first element of new block
600	nblocks = nblocks+1
601	IF (nblocks .GT. nsfm) THEN
602	! IF (nblocks .GT. 2) THEN
603	STOP 'Problem in teststomate'
604	ENDIF
605	ifirst(nblocks) = iisf
606	ilast(nblocks) = iisf
607	ENDIF
608	ELSE
609	!-------- found first element of new block
610	nblocks = nblocks+1
611	IF (nblocks .GT. nsfm) THEN
612	! IF (nblocks .GT. 2) THEN
613	STOP 'Problem in teststomate'
614	ENDIF
615	ifirst(nblocks) = iisf
616	ilast(nblocks) = iisf
617	ENDIF
618	ENDDO
619	!-----
620	DO iblocks=1,nblocks
621	IF (ifirst(iblocks) .NE. ilast(iblocks)) THEN
622	ndim = 2;
623	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
624	count(1:ndim) = SHAPE(clay_fm)
625	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
626	ier = NF90_INQ_VARID (force_id,'clay',v_id)
627	a_er = a_er.OR.(ier.NE.0)
628	ier = NF90_GET_VAR (force_id,v_id, &
629	& clay_fm(:,ifirst(iblocks):ilast(iblocks)), &
630	& start=start(1:ndim), count=count(1:ndim))
631	a_er = a_er.OR.(ier.NE.0)
632	!---------
633	ndim = 3;
634	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
635	count(1:ndim) = SHAPE(humrel_x_fm)
636	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
637	ier = NF90_INQ_VARID (force_id,'humrel',v_id)
638	a_er = a_er.OR.(ier.NE.0)
639	ier = NF90_GET_VAR (force_id,v_id, &
640	& humrel_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
641	& start=start(1:ndim), count=count(1:ndim))
642	a_er = a_er.OR.(ier.NE.0)
643	!---------
644	ndim = 2;
645	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
646	count(1:ndim) = SHAPE(litterhum_fm)
647	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
648	ier = NF90_INQ_VARID (force_id,'litterhum',v_id)
649	a_er = a_er.OR.(ier.NE.0)
650	ier = NF90_GET_VAR (force_id,v_id, &
651	& litterhum_fm(:,ifirst(iblocks):ilast(iblocks)), &
652	& start=start(1:ndim), count=count(1:ndim))
653	a_er = a_er.OR.(ier.NE.0)
654	!---------
655	ndim = 2;
656	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
657	count(1:ndim) = SHAPE(t2m_fm)
658	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
659	ier = NF90_INQ_VARID (force_id,'t2m',v_id)
660	a_er = a_er.OR.(ier.NE.0)
661	ier = NF90_GET_VAR (force_id,v_id, &
662	& t2m_fm(:,ifirst(iblocks):ilast(iblocks)), &
663	& start=start(1:ndim), count=count(1:ndim))
664	a_er = a_er.OR.(ier.NE.0)
665	!---------
666	ndim = 2;
667	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
668	count(1:ndim) = SHAPE(t2m_min_fm)
669	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
670	ier = NF90_INQ_VARID (force_id,'t2m_min',v_id)
671	a_er = a_er.OR.(ier.NE.0)
672	ier = NF90_GET_VAR (force_id,v_id, &
673	& t2m_min_fm(:,ifirst(iblocks):ilast(iblocks)), &
674	& start=start(1:ndim), count=count(1:ndim))
675	a_er = a_er.OR.(ier.NE.0)
676	!---------
677	ndim = 2;
678	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
679	count(1:ndim) = SHAPE(temp_sol_fm)
680	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
681	ier = NF90_INQ_VARID (force_id,'tsurf',v_id)
682	a_er = a_er.OR.(ier.NE.0)
683	ier = NF90_GET_VAR (force_id,v_id, &
684	& temp_sol_fm(:,ifirst(iblocks):ilast(iblocks)), &
685	& start=start(1:ndim), count=count(1:ndim))
686	a_er = a_er.OR.(ier.NE.0)
687	!---------
688	ndim = 3;
689	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
690	count(1:ndim) = SHAPE(soiltemp_fm)
691	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
692	ier = NF90_INQ_VARID (force_id,'tsoil',v_id)
693	a_er = a_er.OR.(ier.NE.0)
694	ier = NF90_GET_VAR (force_id,v_id, &
695	& soiltemp_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
696	& start=start(1:ndim), count=count(1:ndim))
697	a_er = a_er.OR.(ier.NE.0)
698	!---------
699	ndim = 3;
700	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
701	count(1:ndim) = SHAPE(soilhum_fm)
702	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
703	ier = NF90_INQ_VARID (force_id,'soilhum',v_id)
704	a_er = a_er.OR.(ier.NE.0)
705	ier = NF90_GET_VAR (force_id,v_id, &
706	& soilhum_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
707	& start=start(1:ndim), count=count(1:ndim))
708	a_er = a_er.OR.(ier.NE.0)
709	!---------
710	ndim = 2;
711	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
712	count(1:ndim) = SHAPE(precip_fm)
713	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
714	ier = NF90_INQ_VARID (force_id,'precip',v_id)
715	a_er = a_er.OR.(ier.NE.0)
716	ier = NF90_GET_VAR (force_id,v_id, &
717	& precip_fm(:,ifirst(iblocks):ilast(iblocks)), &
718	& start=start(1:ndim), count=count(1:ndim))
719	a_er = a_er.OR.(ier.NE.0)
720	!---------
721	ndim = 3;
722	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
723	count(1:ndim) = SHAPE(gpp_x_fm)
724	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
725	ier = NF90_INQ_VARID (force_id,'gpp',v_id)
726	a_er = a_er.OR.(ier.NE.0)
727	ier = NF90_GET_VAR (force_id,v_id, &
728	& gpp_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
729	& start=start(1:ndim), count=count(1:ndim))
730	a_er = a_er.OR.(ier.NE.0)
731	!---------
732	ndim = 3;
733	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
734	count(1:ndim) = SHAPE(veget_force_x_fm)
735	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
736	ier = NF90_INQ_VARID (force_id,'veget',v_id)
737	a_er = a_er.OR.(ier.NE.0)
738	ier = NF90_GET_VAR (force_id,v_id, &
739	& veget_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
740	& start=start(1:ndim), count=count(1:ndim))
741	a_er = a_er.OR.(ier.NE.0)
742	!---------
743	ndim = 3;
744	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
745	count(1:ndim) = SHAPE(veget_max_force_x_fm)
746	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
747	ier = NF90_INQ_VARID (force_id,'veget_max',v_id)
748	a_er = a_er.OR.(ier.NE.0)
749	ier = NF90_GET_VAR (force_id,v_id, &
750	& veget_max_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
751	& start=start(1:ndim), count=count(1:ndim))
752	a_er = a_er.OR.(ier.NE.0)
753	!---------
754	ndim = 3;
755	start(:) = 1; start(ndim) = isf(ifirst(iblocks));
756	count(1:ndim) = SHAPE(lai_force_x_fm)
757	count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1
758	ier = NF90_INQ_VARID (force_id,'lai',v_id)
759	a_er = a_er.OR.(ier.NE.0)
760	ier = NF90_GET_VAR (force_id,v_id, &
761	& lai_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), &
762	& start=start(1:ndim), count=count(1:ndim))
763	a_er = a_er.OR.(ier.NE.0)
764	ENDIF
765	ENDDO
766	!-----
767	!---- determine which forcing states must be read next time
768	!-----
769	isf(1) = isf(nsfm)+1
770	IF ( isf(1) .GT. nsft ) isf(1) = 1
771	DO iisf = 2, nsfm
772	isf(iisf) = isf(iisf-1)+1
773	IF ( isf(iisf) .GT. nsft ) isf(iisf) = 1
774	ENDDO
775	!---- start again at first forcing state
776	iisf = 1
777	ENDIF
778	soiltype(:,3) = clay_fm(:,iisf)
779	humrel_x(:,:) = humrel_x_fm(:,:,iisf)
780	litterhum(:) = litterhum_fm(:,iisf)
781	t2m(:) = t2m_fm(:,iisf)
782	t2m_min(:) = t2m_min_fm(:,iisf)
783	temp_sol(:) = temp_sol_fm(:,iisf)
784	soiltemp(:,:) = soiltemp_fm(:,:,iisf)
785	soilhum(:,:) = soilhum_fm(:,:,iisf)
786	precip_rain(:) = precip_fm(:,iisf)
787	gpp_x(:,:) = gpp_x_fm(:,:,iisf)
788	veget_force_x(:,:) = veget_force_x_fm(:,:,iisf)
789	veget_max_force_x(:,:) = veget_max_force_x_fm(:,:,iisf)
790	lai_force_x(:,:) = lai_force_x_fm(:,:,iisf)
791	WHERE ( t2m(:) .LT. ZeroCelsius )
792	precip_snow(:) = precip_rain(:)
793	precip_rain(:) = 0.
794	ELSEWHERE
795	precip_snow(:) = 0.
796	ENDWHERE
797	!---
798	!-- scale GPP to new lai and veget_max
799	!---
800	WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
801	!-- scale GPP to new LAI
802	WHERE (lai_force_x(:,:) .GT. 0.0 )
803	gpp_x(:,:) = gpp_x(:,:)atan(2.lai_x(:,:)) &
804	& /atan( 2.*MAX(lai_force_x(:,:),0.01))
805	ENDWHERE
806	!-- scale GPP to new veget_max
807	WHERE (veget_max_force_x(:,:) .GT. 0.0 )
808	gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
809	ENDWHERE
810	!---
811	!-- number crunching
812	!---
813	CALL intsurf_time( itau, date0, dtradia )
814	ldrestart_read = .FALSE.
815	ldrestart_write = .FALSE.
816	CALL slowproc_main &
817	& (itau, kjpij, kjpindex, dt_force, date0, &
818	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
819	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
820	& t2m, t2m_min, temp_sol, soiltemp, &
821	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
822	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
823	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
824	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, &
825	!!>> DS add 01/07/2011 correction due to some modifications in the trunk
826	& fco2_lu)
827	day_counter = one_day - dt_force
828	ENDDO ! end of the time loop
829	!-
830	itau = itau -itau_step
831	!-
832	! write restart files
833	!-
834	! first, read and write variables that are not managed otherwise
835	taboo_vars = &
836	& '$lat$ $lon$ $lev$ $veget_year$ '// &
837	& '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
838	& '$lai$ $soiltype_frac$ $clay_frac$ '// &
839	& '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
840	!-
841	CALL ioget_vname(rest_id_sec, nbvar, varnames)
842	!!$!-
843	!!$! read and write some special variables (1D or variables that we need)
844	!!$!-
845	!!$ var_name = 'day_counter'
846	!!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
847	!!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
848	!!$!-
849	!!$ var_name = 'dt_days'
850	!!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
851	!!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
852	!!$!-
853	!!$ var_name = 'date'
854	!!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp)
855	!!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp)
856	!-
857	DO iv = 1, nbvar
858	!-- check if the variable is to be written here
859	IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
860	!---- get variable dimensions, especially 3rd dimension
861	CALL ioget_vdim &
862	& (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
863	l1d = ALL(vardims(1:varnbdim) .EQ. 1)
864	ALLOCATE(var_3d(kjpindex,vardims(3)),stat=ier)
865	IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
866	!---- read it
867	IF (l1d) THEN
868	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
869	1,vardims(3),1,itau_dep,.TRUE.,var_3d)
870	ELSE
871	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
872	kjpindex,vardims(3),1,itau_dep,.TRUE.,var_3d, &
873	"gather",kjpindex,indices)
874	ENDIF
875	!---- write it
876	IF (l1d) THEN
877	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
878	1,vardims(3),1,itau,var_3d)
879	ELSE
880	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
881	kjpindex,vardims(3),1,itau,var_3d, &
882	'scatter',kjpindex,indices)
883	ENDIF
884	DEALLOCATE (var_3d)
885	ENDIF
886	ENDDO
887	! call slowproc to write restart files
888	ldrestart_read = .FALSE.
889	ldrestart_write = .TRUE.
890	!-
891	CALL slowproc_main &
892	& (itau, kjpij, kjpindex, dt_force, date0, &
893	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
894	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
895	& t2m, t2m_min, temp_sol, soiltemp, &
896	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
897	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
898	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
899	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, &
900	!!>> DS add 01/07/2011 correction due to some modifications in the trunk
901	& fco2_lu)
902
903	!-
904	! close files
905	!-
906	CALL restclo
907	CALL histclo
908	ier = NF90_CLOSE (force_id)
909	!-
910	! write a new driver restart file with correct time step
911	!-
912	write(,) 'teststomate: writing driver restart file with correct time step.'
913	dri_restname_in = 'driver_start.nc'
914	CALL getin ('RESTART_FILEIN',dri_restname_in)
915	dri_restname_out = 'driver_restart.nc'
916	CALL getin ('RESTART_FILEOUT',dri_restname_out)
917	CALL SYSTEM &
918	& ('cp '//TRIM(dri_restname_in)//' '//TRIM(dri_restname_out))
919	!-
920	iret = NF90_OPEN (TRIM(sec_restname_out),NF90_NOWRITE,ncfid)
921	iret = NF90_INQ_VARID (ncfid,'time',v_id)
922	iret = NF90_GET_VAR (ncfid,v_id,r1d)
923	time_sec = r1d(1)
924	iret = NF90_INQ_VARID (ncfid,'time_steps',v_id)
925	iret = NF90_GET_VAR (ncfid,v_id,time_step_sec)
926	iret = NF90_CLOSE (ncfid)
927	!-
928	iret = NF90_OPEN (TRIM(dri_restname_out),NF90_WRITE,ncfid)
929	iret = NF90_INQ_VARID (ncfid,'time',v_id)
930	iret = NF90_GET_VAR (ncfid,v_id,r1d)
931	time_dri = r1d(1)
932	r1d(1) = time_sec
933	iret = NF90_PUT_VAR (ncfid,v_id,r1d)
934	iret = NF90_INQ_VARID (ncfid,'time_steps',v_id)
935	iret = NF90_GET_VAR (ncfid,v_id,time_step_dri)
936	iret = NF90_PUT_VAR (ncfid,v_id,time_step_sec)
937	iret = NF90_INQ_VARID (ncfid,'julian',v_id)
938	iret = NF90_GET_VAR (ncfid,v_id,r1d)
939	julian = r1d(1)
940	djulian = (time_step_sec-time_step_dri)*dtradia/one_day
941	julian = julian &
942	& +djulian-FLOAT(INT((julian+djulian)/one_year))*one_year
943	r1d(1) = julian
944	iret = NF90_PUT_VAR (ncfid,v_id,r1d)
945	iret = NF90_CLOSE (ncfid)
946
947	CALL getin_dump
948	!---------------
949	END PROGRAM teststomate
950	!
951	!===
952	!

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

Download in other formats: