Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE_OL/teststomate.f90 @ 2369

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

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

Download in other formats: