Context Navigation

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

Last change on this file since 472 was 472, checked in by didier.solyga, 13 years ago
Replace ok_sechiba, ok_stomate by the following structure control as a argument for the routine pft_parameters_main
File size: 35.2 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, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time,impose_param
24	USE parallel
25	!-
26	IMPLICIT NONE
27	!-
28	! Declarations
29	!-
30	INTEGER(i_std) :: kjpij,kjpindex
31	REAL(r_std) :: dtradia,dt_force
32
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 :: t2m,t2m_min,temp_sol
44	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum
45	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: humrel_x
46	REAL(r_std),DIMENSION(:),ALLOCATABLE :: litterhum
47	REAL(r_std),DIMENSION(:),ALLOCATABLE :: precip_rain,precip_snow
48	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: gpp_x
49	REAL(r_std),DIMENSION(:),ALLOCATABLE :: deadleaf_cover
50	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: assim_param_x
51	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: height_x
52	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x
53	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux
54	REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu
55
56	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g
57	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g
58	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g
59
60	INTEGER :: ier,iret
61	INTEGER :: ncfid
62	CHARACTER(LEN=20),SAVE :: thecalendar='noleap'
63
64	LOGICAL :: a_er
65	CHARACTER(LEN=80) :: &
66	& dri_restname_in,dri_restname_out, &
67	& sec_restname_in,sec_restname_out, &
68	& sto_restname_in,sto_restname_out, &
69	& stom_histname, stom_ipcc_histname
70	INTEGER(i_std) :: iim,jjm,llm
71	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
72	REAL, ALLOCATABLE, DIMENSION(:) :: lev
73	LOGICAL :: rectilinear
74	REAL, ALLOCATABLE, DIMENSION(:) :: lon_rect, lat_rect
75	REAL(r_std) :: dt
76	INTEGER(i_std) :: itau_dep,itau_fin,itau,itau_len,itau_step
77	REAL(r_std) :: date0
78	INTEGER(i_std) :: rest_id_sec,rest_id_sto
79	INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC
80	CHARACTER(LEN=30) :: time_str
81	REAL(r_std) :: dt_slow_
82	REAL :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc
83	!!$ >> DS for the externalized version
84	REAL(r_std),ALLOCATABLE,DIMENSION(:) :: hist_PFTaxis
85	!!$ >> DS
86	REAL(r_std),DIMENSION(10) :: hist_pool_10axis
87	REAL(r_std),DIMENSION(100) :: hist_pool_100axis
88	REAL(r_std),DIMENSION(11) :: hist_pool_11axis
89	REAL(r_std),DIMENSION(101) :: hist_pool_101axis
90	INTEGER :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id,hori_id
91	INTEGER :: hist_pool_10axis_id
92	INTEGER :: hist_pool_100axis_id
93	INTEGER :: hist_pool_11axis_id
94	INTEGER :: hist_pool_101axis_id
95	INTEGER(i_std) :: i,j,iv
96	LOGICAL :: ldrestart_read,ldrestart_write
97	LOGICAL :: l1d
98	INTEGER(i_std),PARAMETER :: nbvarmax=200
99	INTEGER(i_std) :: nbvar
100	CHARACTER*50,DIMENSION(nbvarmax) :: varnames
101	INTEGER(i_std) :: varnbdim
102	INTEGER(i_std),PARAMETER :: varnbdim_max=20
103	INTEGER,DIMENSION(varnbdim_max) :: vardims
104	CHARACTER*200 :: taboo_vars
105	REAL(r_std),DIMENSION(1) :: xtmp
106	INTEGER :: nsfm,nsft
107	INTEGER :: iisf,iiisf
108	INTEGER(i_std) :: max_totsize,totsize_1step,totsize_tmp
109
110	INTEGER :: vid
111	CHARACTER(LEN=100) :: forcing_name
112	REAL :: x
113
114	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d
115	REAL(r_std) :: var_1d(1)
116
117	!-
118	REAL(r_std) :: time_sec,time_step_sec
119	REAL(r_std) :: time_dri,time_step_dri
120	REAL(r_std),DIMENSION(1) :: r1d
121	REAL(r_std) :: julian,djulian
122
123	INTEGER(i_std) :: ji,jv,l
124
125	LOGICAL :: debug
126
127	!---------------------------------------------------------------------
128
129	CALL init_para(.FALSE.)
130	CALL init_timer
131
132	IF (is_root_prc) THEN
133	!-
134	! open STOMATE's forcing file to read some basic info
135	!-
136	forcing_name = 'stomate_forcing.nc'
137	CALL getin ('STOMATE_FORCING_NAME',forcing_name)
138	iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id)
139	IF (iret /= NF90_NOERR) THEN
140	CALL ipslerr (3,'teststomate', &
141	& 'Could not open file : ', &
142	& forcing_name,'(Do you have forget it ?)')
143	ENDIF
144	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia)
145	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force)
146	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x)
147	nsft = NINT(x)
148	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x)
149	kjpij = NINT(x)
150	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x)
151	nbp_glo = NINT(x)
152	ENDIF
153	CALL bcast(dtradia)
154	CALL bcast(dt_force)
155	CALL bcast(nsft)
156	CALL bcast(nbp_glo)
157	!-
158	write(numout,*) 'ATTENTION',dtradia,dt_force
159	!-
160	! read info about land points
161	!-
162	IF (is_root_prc) THEN
163	a_er=.FALSE.
164	ALLOCATE (indices_g(nbp_glo),stat=ier)
165	a_er = a_er .OR. (ier.NE.0)
166	IF (a_er) THEN
167	CALL ipslerr (3,'teststomate', &
168	& 'PROBLEM WITH ALLOCATION', &
169	& 'for local variables 1','')
170	ENDIF
171	!
172	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
173	a_er = a_er .OR. (ier.NE.0)
174	IF (a_er) THEN
175	CALL ipslerr (3,'teststomate', &
176	& 'PROBLEM WITH ALLOCATION', &
177	& 'for global variables 1','')
178	ENDIF
179	ier = NF90_INQ_VARID (forcing_id,'index',vid)
180	IF (ier .NE. 0) THEN
181	CALL ipslerr (3,'teststomate', &
182	& 'PROBLEM WITH ALLOCATION', &
183	& 'for global variables 1','')
184	ENDIF
185	ier = NF90_GET_VAR (forcing_id,vid,x_indices_g)
186	IF (iret /= NF90_NOERR) THEN
187	CALL ipslerr (3,'teststomate', &
188	& 'PROBLEM WITH variable "index" in file ', &
189	& forcing_name,'(check this file)')
190	ENDIF
191	indices_g(:) = NINT(x_indices_g(:))
192	DEALLOCATE (x_indices_g)
193	ELSE
194	ALLOCATE (indices_g(0))
195	ENDIF
196	!---------------------------------------------------------------------
197	!-
198	! set debug to have more information
199	!-
200	!Config Key = DEBUG_INFO
201	!Config Desc = Flag for debug information
202	!Config Def = n
203	!Config Help = This option allows to switch on the output of debug
204	!Config information without recompiling the code.
205	!-
206	debug = .FALSE.
207	CALL getin_p('DEBUG_INFO',debug)
208	!
209	!Config Key = LONGPRINT
210	!Config Desc = ORCHIDEE will print more messages
211	!Config Def = n
212	!Config Help = This flag permits to print more debug messages in the run.
213	!
214	long_print = .FALSE.
215	CALL getin_p('LONGPRINT',long_print)
216
217	!-
218	! activate CO2, STOMATE, but not sechiba
219	!-
220	control%river_routing = .FALSE.
221	control%hydrol_cwrr = .FALSE.
222	control%ok_sechiba = .FALSE.
223	!
224	control%stomate_watchout = .TRUE.
225	control%ok_co2 = .TRUE.
226	control%ok_stomate = .TRUE.
227	!-
228	! is DGVM activated?
229	!-
230	control%ok_dgvm = .FALSE.
231	CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm)
232	WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm
233
234	!!$ >> DS added for the externalised version
235
236	! 1. Number of PFTs defined by the user
237
238	!Config Key = NVM
239	!Config Desc = number of PFTs
240	!Config if = ANYTIME
241	!Config Def = 13
242	!Config Help = The number of vegetation types define by the user
243	!Config Units = NONE
244	!
245	CALL getin_p('NVM',nvm)
246	WRITE(numout,*)'the number of pfts is : ', nvm
247
248	! 2. Should we read the parameters in the run.def file ?
249
250	!Config Key = IMPOSE_PARAM
251	!Config Desc = Do you impose the values of the parameters?
252	!Config if = ANYTIME
253	!Config Def = y
254	!Config Help = This flag can deactivate the reading of some parameters.
255	! Useful if you want to use the standard values without commenting the run.def
256	!Config Units = NONE
257	!
258	CALL getin_p('IMPOSE_PARAM',impose_param)
259
260	! 3. Allocate and intialize the pft parameters
261
262	CALL pft_parameters_main(control)
263
264	! 4. Activation sub-models of ORCHIDEE
265
266	CALL activate_sub_models(control)
267
268	! 5. Vegetation configuration (impose_veg, land_use, lcchange...previously in slowproc)
269
270	CALL veget_config
271
272	! 6. Read the parameters in the run.def file
273
274	IF ( impose_param) THEN
275	CALL getin_pft_parameters
276	CALL getin_stomate_pft_parameters
277	CALL getin_co2_parameters
278	CALL getin_stomate_parameters
279	ENDIF
280	!-
281	IF (control%ok_dgvm) THEN
282	IF ( impose_param ) THEN
283	CALL getin_dgvm_parameters
284	ENDIF
285	ENDIF
286
287	!!$ >> DS
288
289	!-
290	! restart files
291	!-
292	IF (is_root_prc) THEN
293	! Sechiba's restart files
294	sec_restname_in = 'sechiba_start.nc'
295	CALL getin('SECHIBA_restart_in',sec_restname_in)
296	WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
297	IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
298	STOP 'Need a restart file for Sechiba'
299	ENDIF
300	sec_restname_out = 'sechiba_rest_out.nc'
301	CALL getin('SECHIBA_rest_out',sec_restname_out)
302	WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
303	! Stomate's restart files
304	sto_restname_in = 'stomate_start.nc'
305	CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
306	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
307	sto_restname_out = 'stomate_rest_out.nc'
308	CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
309	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
310
311	!-
312	! We need to know iim, jjm.
313	! Get them from the restart files themselves.
314	!-
315	iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
316	IF (iret /= NF90_NOERR) THEN
317	CALL ipslerr (3,'teststomate', &
318	& 'Could not open file : ', &
319	& sec_restname_in,'(Do you have forget it ?)')
320	ENDIF
321	iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
322	iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
323	iret = NF90_INQ_VARID (ncfid, "time", iv)
324	iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
325	iret = NF90_CLOSE (ncfid)
326	i=INDEX(thecalendar,ACHAR(0))
327	IF ( i > 0 ) THEN
328	thecalendar(i:20)=' '
329	ENDIF
330	ENDIF
331	CALL bcast(iim_g)
332	CALL bcast(jjm_g)
333	CALL bcast(thecalendar)
334	!-
335	! calendar
336	!-
337	CALL ioconf_calendar (thecalendar)
338	CALL ioget_calendar (one_year,one_day)
339	!
340	! Parallelization :
341	!
342	CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
343	kjpindex=nbp_loc
344	jjm=jj_nb
345	iim=iim_g
346	kjpij=iim*jjm
347	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
348	!-
349	!-
350	! read info about grids
351	!-
352	!-
353	llm=1
354	ALLOCATE(lev(llm))
355	IF (is_root_prc) THEN
356	!-
357	ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
358	ier = NF90_GET_VAR (forcing_id,vid,lalo_g)
359	!-
360	ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
361	ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
362	ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g)
363	neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
364	DEALLOCATE (x_neighbours_g)
365	!-
366	ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
367	ier = NF90_GET_VAR (forcing_id,vid,resolution_g)
368	!-
369	ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
370	ier = NF90_GET_VAR (forcing_id,vid,contfrac_g)
371
372	lon_g(:,:) = 0.0
373	lat_g(:,:) = 0.0
374	lev(1) = 0.0
375	!-
376	CALL restini &
377	& (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
378	& sec_restname_out, itau_dep, date0, dt, rest_id_sec)
379	!-
380	IF ( dt .NE. dtradia ) THEN
381	WRITE(numout,*) 'dt',dt
382	WRITE(numout,*) 'dtradia',dtradia
383	CALL ipslerr (3,'teststomate', &
384	& 'PROBLEM with time steps.', &
385	& sec_restname_in,'(dt .NE. dtradia)')
386	ENDIF
387	!-
388	CALL restini &
389	& (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
390	& sto_restname_out, itau_dep, date0, dt, rest_id_sto)
391	!-
392	IF ( dt .NE. dtradia ) THEN
393	WRITE(numout,*) 'dt',dt
394	WRITE(numout,*) 'dtradia',dtradia
395	CALL ipslerr (3,'teststomate', &
396	& 'PROBLEM with time steps.', &
397	& sto_restname_in,'(dt .NE. dtradia)')
398	ENDIF
399	ENDIF
400	CALL bcast(rest_id_sec)
401	CALL bcast(rest_id_sto)
402	CALL bcast(itau_dep)
403	CALL bcast(date0)
404	CALL bcast(dt)
405	CALL bcast(lev)
406	!---
407	!--- Create the index table
408	!---
409	!--- This job return a LOCAL kindex
410	!---
411	ALLOCATE (indices(kjpindex),stat=ier)
412	IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo)
413	CALL scatter(indices_g,indices)
414	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
415	IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex)
416
417	!---
418	!--- initialize global grid
419	!---
420	CALL init_grid( kjpindex )
421	CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g)
422
423	!---
424	!--- initialize local grid
425	!---
426	jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
427	if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
428	ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
429	IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
430	ALLOCATE(lon(iim,jjm))
431	ALLOCATE(lat(iim,jjm))
432	lon=zero
433	lat=zero
434	CALL scatter2D(lon_g,lon)
435	CALL scatter2D(lat_g,lat)
436
437	DO ji=1,kjpindex
438
439	j = jlandindex(ji)
440	i = ilandindex(ji)
441
442	!- Create the internal coordinate table
443	!-
444	lalo(ji,1) = lat(i,j)
445	lalo(ji,2) = lon(i,j)
446	ENDDO
447	CALL scatter(neighbours_g,neighbours)
448	CALL scatter(resolution_g,resolution)
449	CALL scatter(contfrac_g,contfrac)
450	CALL scatter(area_g,area)
451	!-
452	!- Check if we have by any change a rectilinear grid. This would allow us to
453	!- simplify the output files.
454	!
455	rectilinear = .FALSE.
456	IF (is_root_prc) THEN
457	IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
458	& ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
459	rectilinear = .TRUE.
460	ENDIF
461	ENDIF
462	CALL bcast(rectilinear)
463	IF (rectilinear) THEN
464	ALLOCATE(lon_rect(iim),stat=ier)
465	IF (ier .NE. 0) THEN
466	WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
467	STOP 'intersurf_history'
468	ENDIF
469	ALLOCATE(lat_rect(jjm),stat=ier)
470	IF (ier .NE. 0) THEN
471	WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
472	STOP 'intersurf_history'
473	ENDIF
474	lon_rect(:) = lon(:,1)
475	lat_rect(:) = lat(1,:)
476	ENDIF
477	!-
478	! allocate arrays
479	!-
480	!
481	a_er = .FALSE.
482	!!$ >> DS added for the externalised version
483	ALLOCATE (hist_PFTaxis(nvm),stat=ier)
484	a_er = a_er .OR. (ier.NE.0)
485	! 01/07/2011 allocate fco2_lu
486	ALLOCATE (fco2_lu(kjpindex),stat=ier)
487	a_er = a_er .OR. (ier.NE.0)
488	!!$ end add DS
489	ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
490	a_er = a_er .OR. (ier.NE.0)
491	ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
492	a_er = a_er .OR. (ier.NE.0)
493	ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
494	a_er = a_er .OR. (ier.NE.0)
495	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
496	a_er = a_er .OR. (ier.NE.0)
497	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
498	a_er = a_er .OR. (ier.NE.0)
499	ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
500	a_er = a_er .OR. (ier.NE.0)
501	ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
502	a_er = a_er .OR. (ier.NE.0)
503	ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
504	a_er = a_er .OR. (ier.NE.0)
505	ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
506	a_er = a_er .OR. (ier.NE.0)
507	ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
508	a_er = a_er .OR. (ier.NE.0)
509	ALLOCATE (t2m(kjpindex),stat=ier)
510	a_er = a_er .OR. (ier.NE.0)
511	ALLOCATE (t2m_min(kjpindex),stat=ier)
512	a_er = a_er .OR. (ier.NE.0)
513	ALLOCATE (temp_sol(kjpindex),stat=ier)
514	a_er = a_er .OR. (ier.NE.0)
515	ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
516	a_er = a_er .OR. (ier.NE.0)
517	ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
518	a_er = a_er .OR. (ier.NE.0)
519	ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
520	a_er = a_er .OR. (ier.NE.0)
521	ALLOCATE (litterhum(kjpindex),stat=ier)
522	a_er = a_er .OR. (ier.NE.0)
523	ALLOCATE (precip_rain(kjpindex),stat=ier)
524	a_er = a_er .OR. (ier.NE.0)
525	ALLOCATE (precip_snow(kjpindex),stat=ier)
526	a_er = a_er .OR. (ier.NE.0)
527	ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
528	a_er = a_er .OR. (ier.NE.0)
529	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
530	a_er = a_er .OR. (ier.NE.0)
531	ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
532	a_er = a_er .OR. (ier.NE.0)
533	ALLOCATE (height_x(kjpindex,nvm),stat=ier)
534	a_er = a_er .OR. (ier.NE.0)
535	ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
536	a_er = a_er .OR. (ier.NE.0)
537	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
538	a_er = a_er .OR. (ier.NE.0)
539	ALLOCATE (fco2_lu(kjpindex),stat=ier)
540	a_er = a_er .OR. (ier.NE.0)
541	IF (a_er) THEN
542	CALL ipslerr (3,'teststomate', &
543	& 'PROBLEM WITH ALLOCATION', &
544	& 'for local variables 1','')
545	ENDIF
546	!
547	! prepare forcing
548	!
549	max_totsize = 50
550	CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
551	max_totsize = max_totsize * 1000000
552
553	totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
554	SIZE(humrel_x)*KIND(humrel_x) + &
555	SIZE(litterhum)*KIND(litterhum) + &
556	SIZE(t2m)*KIND(t2m) + &
557	SIZE(t2m_min)*KIND(t2m_min) + &
558	SIZE(temp_sol)*KIND(temp_sol) + &
559	SIZE(soiltemp)*KIND(soiltemp) + &
560	SIZE(soilhum)*KIND(soilhum) + &
561	SIZE(precip_rain)*KIND(precip_rain) + &
562	SIZE(precip_snow)*KIND(precip_snow) + &
563	SIZE(gpp_x)*KIND(gpp_x) + &
564	SIZE(veget_force_x)*KIND(veget_force_x) + &
565	SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
566	SIZE(lai_force_x)*KIND(lai_force_x)
567	CALL reduce_sum(totsize_1step,totsize_tmp)
568	CALL bcast(totsize_tmp)
569	totsize_1step=totsize_tmp
570
571	! total number of forcing steps
572	IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
573	CALL ipslerr (3,'teststomate', &
574	& 'stomate: error with total number of forcing steps', &
575	& 'nsft','teststomate computation different with forcing file value.')
576	ENDIF
577	! number of forcing steps in memory
578	nsfm = MIN(nsft, &
579	& MAX(1,NINT( REAL(max_totsize,r_std) &
580	& /REAL(totsize_1step,r_std))))
581	!-
582	WRITE(numout,*) 'Offline forcing of Stomate:'
583	WRITE(numout,*) ' Total number of forcing steps:',nsft
584	WRITE(numout,*) ' Number of forcing steps in memory:',nsfm
585	!-
586	CALL init_forcing(kjpindex,nsfm,nsft)
587	!-
588	! ensure that we read all new forcing states
589	iisf = nsfm
590	! initialize that table that contains the indices
591	! of the forcing states that will be in memory
592	isf(:) = (/ (i,i=1,nsfm) /)
593
594	nf_written(:) = .FALSE.
595	nf_written(isf(:)) = .TRUE.
596
597	!-
598	! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
599	!-
600	itau_step = NINT(dt_force/dtradia)
601	IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
602	!
603	CALL ioconf_startdate(date0)
604	CALL intsurf_time( itau_dep, date0, dtradia )
605	!-
606	! Length of integration
607	!-
608	WRITE(time_str,'(a)') '1Y'
609	CALL getin_p ('TIME_LENGTH', time_str)
610	! transform into itau
611	CALL tlen2itau(time_str, dt, date0, itau_len)
612	! itau_len*dtradia must be a multiple of dt_force
613	itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
614	& *dt_force/dtradia)
615	!-
616	itau_fin = itau_dep+itau_len
617	!-
618	! set up STOMATE history file
619	!-
620	!Config Key = STOMATE_OUTPUT_FILE
621	!Config Desc = Name of file in which STOMATE's output is going
622	!Config to be written
623	!Config Def = stomate_history.nc
624	!Config Help = This file is going to be created by the model
625	!Config and will contain the output from the model.
626	!Config This file is a truly COADS compliant netCDF file.
627	!Config It will be generated by the hist software from
628	!Config the IOIPSL package.
629	!-
630	stom_histname='stomate_history.nc'
631	CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
632	WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
633	!-
634	!Config Key = STOMATE_HIST_DT
635	!Config Desc = STOMATE history time step (d)
636	!Config Def = 10.
637	!Config Help = Time step of the STOMATE history file
638	!-
639	hist_days_stom = 10.
640	CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
641	IF ( hist_days_stom == -1. ) THEN
642	hist_dt_stom = -1.
643	WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
644	ELSE
645	hist_dt_stom = NINT( hist_days_stom ) * one_day
646	WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
647	hist_dt_stom/one_day
648	ENDIF
649	!-
650	!-
651	!- initialize
652	WRITE(numout,*) "before histbeg : ",date0,dt
653	IF ( rectilinear ) THEN
654	#ifdef CPP_PARA
655	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
656	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
657	#else
658	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
659	& itau_dep, date0, dt, hori_id, hist_id_stom)
660	#endif
661	ELSE
662	#ifdef CPP_PARA
663	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
664	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
665	#else
666	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
667	& itau_dep, date0, dt, hori_id, hist_id_stom)
668	#endif
669	ENDIF
670	!- define PFT axis
671	hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
672	!- declare this axis
673	CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
674	& '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
675	!- define Pool_10 axis
676	hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
677	!- declare this axis
678	CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
679	& '1', 10, hist_pool_10axis, hist_pool_10axis_id)
680
681	!- define Pool_100 axis
682	hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
683	!- declare this axis
684	CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
685	& '1', 100, hist_pool_100axis, hist_pool_100axis_id)
686
687	!- define Pool_11 axis
688	hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
689	!- declare this axis
690	CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
691	& '1', 11, hist_pool_11axis, hist_pool_11axis_id)
692
693	!- define Pool_101 axis
694	hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
695	!- declare this axis
696	CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
697	& '1', 101, hist_pool_101axis, hist_pool_101axis_id)
698
699	!- define STOMATE history file
700	CALL stom_define_history (hist_id_stom, nvm, iim, jjm, &
701	& dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
702	& hist_pool_10axis_id, hist_pool_100axis_id, &
703	& hist_pool_11axis_id, hist_pool_101axis_id)
704
705	!- end definition
706	CALL histend(hist_id_stom)
707	!-
708	!-
709	! STOMATE IPCC OUTPUTS IS ACTIVATED
710	!-
711	!Config Key = STOMATE_IPCC_OUTPUT_FILE
712	!Config Desc = Name of file in which STOMATE's output is going
713	!Config to be written
714	!Config Def = stomate_ipcc_history.nc
715	!Config Help = This file is going to be created by the model
716	!Config and will contain the output from the model.
717	!Config This file is a truly COADS compliant netCDF file.
718	!Config It will be generated by the hist software from
719	!Config the IOIPSL package.
720	!-
721	stom_ipcc_histname='stomate_ipcc_history.nc'
722	CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)
723	WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
724	!-
725	!Config Key = STOMATE_IPCC_HIST_DT
726	!Config Desc = STOMATE IPCC history time step (d)
727	!Config Def = 0.
728	!Config Help = Time step of the STOMATE IPCC history file
729	!-
730	hist_days_stom_ipcc = zero
731	CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)
732	IF ( hist_days_stom_ipcc == moins_un ) THEN
733	hist_dt_stom_ipcc = moins_un
734	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
735	ELSE
736	hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
737	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
738	hist_dt_stom_ipcc/one_day
739	ENDIF
740
741	! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
742	dt_slow_ = one_day
743	CALL getin_p('DT_SLOW', dt_slow_)
744	IF ( hist_days_stom_ipcc > zero ) THEN
745	IF (dt_slow_ > hist_dt_stom_ipcc) THEN
746	WRITE(numout,*) "DT_SLOW = ",dt_slow_," , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
747	CALL ipslerr (3,'intsurf_history', &
748	& 'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
749	& '(must be less or equal)')
750	ENDIF
751	ENDIF
752
753	IF ( hist_dt_stom_ipcc == 0 ) THEN
754	hist_id_stom_ipcc = -1
755	ELSE
756	!-
757	!- initialize
758	IF ( rectilinear ) THEN
759	#ifdef CPP_PARA
760	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
761	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
762	#else
763	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
764	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
765	#endif
766	ELSE
767	#ifdef CPP_PARA
768	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
769	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
770	#else
771	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
772	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
773	#endif
774	ENDIF
775	!- declare this axis
776	CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
777	& '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)
778
779	!- define STOMATE history file
780	CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
781	& dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)
782
783	!- end definition
784	CALL histend(hist_id_stom_IPCC)
785
786	ENDIF
787	!
788	CALL histwrite(hist_id_stom, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
789	IF ( hist_id_stom_IPCC > 0 ) THEN
790	CALL histwrite(hist_id_stom_IPCC, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
791	ENDIF
792	!
793	hist_id_sec = -1
794	hist_id_sec2 = -1
795	!-
796	! first call of slowproc to initialize variables
797	!-
798	itau = itau_dep
799	!
800	DO ji=1,kjpindex
801	DO jv=1,nvm
802	indexveg((jv-1)kjpindex + ji) = indices(ji) + (jv-1)kjpij
803	ENDDO
804	ENDDO
805	!-
806	!MM Problem here with dpu which depends on soil type
807	DO l = 1, nbdl-1
808	! first 2.0 is dpu
809	! second 2.0 is average
810	diaglev(l) = dpu_cste/(2*(nbdl-1) -1) ( ( 2(l-1) -1) + ( 2(l) -1) ) / 2.0
811	ENDDO
812	diaglev(nbdl) = dpu_cste
813	!
814
815	! >> DS : getin the value of slowproc parameters
816	!
817	!Config Key = CLAYFRACTION_DEFAULT
818	!Config Desc =
819	!Config If = OK_SECHIBA
820	!Config Def = 0.2
821	!Config Help =
822	!Config Units = NONE
823	CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
824	!
825	!Config Key = MIN_VEGFRAC
826	!Config Desc = Minimal fraction of mesh a vegetation type can occupy
827	!Config If = OK_SECHIBA
828	!Config Def = 0.001
829	!Config Help =
830	!Config Units = NONE
831	CALL getin_p('MIN_VEGFRAC',min_vegfrac)
832	!
833	!Config Key = STEMPDIAG_BID
834	!Config Desc = only needed for an initial LAI if there is no restart file
835	!Config If = OK_SECHIBA
836	!Config Def = 280.
837	!Config Help =
838	!Config Units =
839	CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
840	!
841	!Config Key = SOILTYPE_DEFAULT
842	!Config Desc = Default soil texture distribution in the following order : sand, loam and clay
843	!Config If = OK_SECHIBA
844	!Config Def = 0.0, 1.0, 0.0
845	!Config Help =
846	!Config Units = NONE
847	CALL getin_p('SOILTYPE_DEFAULT',soiltype_default)
848	! >> DS
849
850	CALL ioget_expval(val_exp)
851	ldrestart_read = .TRUE.
852	ldrestart_write = .FALSE.
853	!-
854	! read some variables we need from SECHIBA's restart file
855	!-
856	CALL slowproc_main &
857	& (itau, kjpij, kjpindex, dt_force, date0, &
858	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
859	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
860	& t2m, t2m_min, temp_sol, soiltemp, &
861	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
862	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
863	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
864	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
865	! correct date
866	day_counter = one_day - dt_force
867	date=1
868	!-
869	! time loop
870	!-
871	IF (debug) check_time=.TRUE.
872	CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
873	l_first_intersurf=.FALSE.
874	!
875	DO itau = itau_dep+itau_step,itau_fin,itau_step
876	!
877	CALL intsurf_time( itau, date0, dtradia )
878	!
879	!-- next forcing state
880	iisf = iisf+1
881	IF (debug) WRITE(numout,*) "itau,iisf : ",itau,iisf
882	!---
883	IF (iisf .GT. nsfm) THEN
884	!-----
885	!---- we have to read new forcing states from the file
886	!-----
887	!---- determine blocks of forcing states that are contiguous in memory
888	!-----
889	CALL forcing_read(forcing_id,nsfm)
890
891	!--------------------------
892
893	!-----
894	!---- determine which forcing states must be read next time
895	!-----
896	isf(1) = isf(nsfm)+1
897	IF ( isf(1) .GT. nsft ) isf(1) = 1
898	DO iiisf = 2, nsfm
899	isf(iiisf) = isf(iiisf-1)+1
900	IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
901	ENDDO
902	nf_written(isf(:)) = .TRUE.
903	!---- start again at first forcing state
904	iisf = 1
905	ENDIF
906	! Bug here ! soiltype(:,3) != clay
907	! soiltype(:,3) = clay_fm(:,iisf)
908	humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
909	litterhum(:) = litterhum_daily_fm(:,iisf)
910	t2m(:) = t2m_daily_fm(:,iisf)
911	t2m_min(:) = t2m_min_daily_fm(:,iisf)
912	temp_sol(:) = tsurf_daily_fm(:,iisf)
913	soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
914	soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
915	precip_rain(:) = precip_fm(:,iisf)
916	gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
917	veget_force_x(:,:) = veget_fm(:,:,iisf)
918	veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
919	lai_force_x(:,:) = lai_fm(:,:,iisf)
920	WHERE ( t2m(:) .LT. ZeroCelsius )
921	precip_snow(:) = precip_rain(:)
922	precip_rain(:) = 0.
923	ELSEWHERE
924	precip_snow(:) = 0.
925	ENDWHERE
926	!---
927	!-- scale GPP to new lai and veget_max
928	!---
929	WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
930	!-- scale GPP to new LAI
931	WHERE (lai_force_x(:,:) .GT. 0.0 )
932	gpp_x(:,:) = gpp_x(:,:)ATAN(2.lai_x(:,:)) &
933	& /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
934	ENDWHERE
935	!-- scale GPP to new veget_max
936	WHERE (veget_max_force_x(:,:) .GT. 0.0 )
937	gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
938	ENDWHERE
939	!---
940	!-- number crunching
941	!---
942	ldrestart_read = .FALSE.
943	ldrestart_write = .FALSE.
944	CALL slowproc_main &
945	& (itau, kjpij, kjpindex, dt_force, date0, &
946	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
947	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
948	& t2m, t2m_min, temp_sol, soiltemp, &
949	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
950	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
951	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
952	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
953	day_counter = one_day - dt_force
954	ENDDO ! end of the time loop
955	!-
956	! write restart files
957	!-
958	IF (is_root_prc) THEN
959	! first, read and write variables that are not managed otherwise
960	taboo_vars = &
961	& '$lat$ $lon$ $lev$ $veget_year$ '// &
962	& '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
963	& '$lai$ $soiltype_frac$ $clay_frac$ '// &
964	& '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
965	!-
966	CALL ioget_vname(rest_id_sec, nbvar, varnames)
967	!-
968	DO iv = 1, nbvar
969	!-- check if the variable is to be written here
970	IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
971	IF (debug) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
972
973	!---- get variable dimensions, especially 3rd dimension
974	CALL ioget_vdim &
975	& (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
976	l1d = ALL(vardims(1:varnbdim) .EQ. 1)
977	!---- read it
978	IF (l1d) THEN
979	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
980	1,1,1,itau_dep,.TRUE.,var_1d)
981	ELSE
982	ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
983	IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
984	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
985	nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
986	"gather",nbp_glo,indices_g)
987	ENDIF
988	!---- write it
989	IF (l1d) THEN
990	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
991	1,1,1,itau_fin,var_1d)
992	ELSE
993	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
994	nbp_glo,vardims(3),1,itau_fin,var_3d, &
995	'scatter',nbp_glo,indices_g)
996	DEALLOCATE (var_3d)
997	ENDIF
998	ENDIF
999	ENDDO
1000	ENDIF
1001	CALL barrier_para()
1002
1003	! call slowproc to write restart files
1004	ldrestart_read = .FALSE.
1005	ldrestart_write = .TRUE.
1006	!-
1007	IF (debug) WRITE(numout,*) "Call slowproc for restart."
1008	CALL slowproc_main &
1009	& (itau_fin, kjpij, kjpindex, dt_force, date0, &
1010	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1011	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
1012	& t2m, t2m_min, temp_sol, soiltemp, &
1013	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1014	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
1015	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
1016	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
1017	!-
1018	! close files
1019	!-
1020	IF (is_root_prc) THEN
1021	CALL restclo
1022	IF ( debug ) WRITE(numout,*) 'REST CLOSED'
1023	ENDIF
1024	CALL histclo
1025
1026	IF (is_root_prc) &
1027	ier = NF90_CLOSE (forcing_id)
1028
1029	IF (is_root_prc) THEN
1030	CALL getin_dump()
1031	ENDIF
1032	#ifdef CPP_PARA
1033	CALL MPI_FINALIZE(ier)
1034	#endif
1035	WRITE(numout,*) "End of teststomate."
1036
1037	!---------------
1038	END PROGRAM teststomate
1039	!
1040	!===
1041	!

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

Download in other formats: