Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE_OL/teststomate.f90 @ 6268

Last change on this file since 6268 was 396, checked in by martial.mancip, 13 years ago
Correct teststomate time steps again. Correct a bug in restart writes. Now teststomate can run with routing activated in SPINUP.
File size: 32.5 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, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time
26	USE parallel
27	!-
28	IMPLICIT NONE
29	!-
30	! Declarations
31	!-
32	INTEGER(i_std) :: kjpij,kjpindex
33	REAL(r_std) :: dtradia,dt_force
34
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 :: t2m,t2m_min,temp_sol
46	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum
47	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: humrel_x
48	REAL(r_std),DIMENSION(:),ALLOCATABLE :: litterhum
49	REAL(r_std),DIMENSION(:),ALLOCATABLE :: precip_rain,precip_snow
50	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: gpp_x
51	REAL(r_std),DIMENSION(:),ALLOCATABLE :: deadleaf_cover
52	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: assim_param_x
53	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: height_x
54	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x
55	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux
56	REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu
57
58	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g
59	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g
60	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g
61
62	INTEGER :: ier,iret
63	INTEGER :: ncfid
64	CHARACTER(LEN=20),SAVE :: thecalendar='noleap'
65
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, &
71	& stom_histname, stom_ipcc_histname
72	INTEGER(i_std) :: iim,jjm,llm
73	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
74	REAL, ALLOCATABLE, DIMENSION(:) :: lev
75	LOGICAL :: rectilinear
76	REAL, ALLOCATABLE, DIMENSION(:) :: lon_rect, lat_rect
77	REAL(r_std) :: dt
78	INTEGER(i_std) :: itau_dep,itau_fin,itau,itau_len,itau_step
79	REAL(r_std) :: date0
80	INTEGER(i_std) :: rest_id_sec,rest_id_sto
81	INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC
82	CHARACTER(LEN=30) :: time_str
83	REAL(r_std) :: dt_slow_
84	REAL :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc
85	REAL,DIMENSION(nvm) :: hist_PFTaxis
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	! activate CO2, STOMATE, but not sechiba
218	!-
219	control%river_routing = .FALSE.
220	control%hydrol_cwrr = .FALSE.
221	control%ok_sechiba = .FALSE.
222	!
223	control%stomate_watchout = .TRUE.
224	control%ok_co2 = .TRUE.
225	control%ok_stomate = .TRUE.
226	!-
227	! is DGVM activated?
228	!-
229	control%ok_dgvm = .FALSE.
230	CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm)
231	WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm
232
233	!-
234	! restart files
235	!-
236	IF (is_root_prc) THEN
237	! Sechiba's restart files
238	sec_restname_in = 'sechiba_start.nc'
239	CALL getin('SECHIBA_restart_in',sec_restname_in)
240	WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
241	IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
242	STOP 'Need a restart file for Sechiba'
243	ENDIF
244	sec_restname_out = 'sechiba_rest_out.nc'
245	CALL getin('SECHIBA_rest_out',sec_restname_out)
246	WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
247	! Stomate's restart files
248	sto_restname_in = 'stomate_start.nc'
249	CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
250	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
251	sto_restname_out = 'stomate_rest_out.nc'
252	CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
253	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
254
255	!-
256	! We need to know iim, jjm.
257	! Get them from the restart files themselves.
258	!-
259	iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
260	IF (iret /= NF90_NOERR) THEN
261	CALL ipslerr (3,'teststomate', &
262	& 'Could not open file : ', &
263	& sec_restname_in,'(Do you have forget it ?)')
264	ENDIF
265	iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
266	iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
267	iret = NF90_INQ_VARID (ncfid, "time", iv)
268	iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
269	iret = NF90_CLOSE (ncfid)
270	i=INDEX(thecalendar,ACHAR(0))
271	IF ( i > 0 ) THEN
272	thecalendar(i:20)=' '
273	ENDIF
274	ENDIF
275	CALL bcast(iim_g)
276	CALL bcast(jjm_g)
277	CALL bcast(thecalendar)
278	!-
279	! calendar
280	!-
281	CALL ioconf_calendar (thecalendar)
282	CALL ioget_calendar (one_year,one_day)
283	!
284	! Parallelization :
285	!
286	CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
287	kjpindex=nbp_loc
288	jjm=jj_nb
289	iim=iim_g
290	kjpij=iim*jjm
291	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
292	!-
293	!-
294	! read info about grids
295	!-
296	!-
297	llm=1
298	ALLOCATE(lev(llm))
299	IF (is_root_prc) THEN
300	!-
301	ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
302	ier = NF90_GET_VAR (forcing_id,vid,lalo_g)
303	!-
304	ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
305	ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
306	ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g)
307	neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
308	DEALLOCATE (x_neighbours_g)
309	!-
310	ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
311	ier = NF90_GET_VAR (forcing_id,vid,resolution_g)
312	!-
313	ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
314	ier = NF90_GET_VAR (forcing_id,vid,contfrac_g)
315
316	lon_g(:,:) = 0.0
317	lat_g(:,:) = 0.0
318	lev(1) = 0.0
319	!-
320	CALL restini &
321	& (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
322	& sec_restname_out, itau_dep, date0, dt, rest_id_sec)
323	!-
324	IF ( dt .NE. dtradia ) THEN
325	WRITE(numout,*) 'dt',dt
326	WRITE(numout,*) 'dtradia',dtradia
327	CALL ipslerr (3,'teststomate', &
328	& 'PROBLEM with time steps.', &
329	& sec_restname_in,'(dt .NE. dtradia)')
330	ENDIF
331	!-
332	CALL restini &
333	& (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
334	& sto_restname_out, itau_dep, date0, dt, rest_id_sto)
335	!-
336	IF ( dt .NE. dtradia ) THEN
337	WRITE(numout,*) 'dt',dt
338	WRITE(numout,*) 'dtradia',dtradia
339	CALL ipslerr (3,'teststomate', &
340	& 'PROBLEM with time steps.', &
341	& sto_restname_in,'(dt .NE. dtradia)')
342	ENDIF
343	ENDIF
344	CALL bcast(rest_id_sec)
345	CALL bcast(rest_id_sto)
346	CALL bcast(itau_dep)
347	CALL bcast(date0)
348	CALL bcast(dt)
349	CALL bcast(lev)
350	!---
351	!--- Create the index table
352	!---
353	!--- This job return a LOCAL kindex
354	!---
355	ALLOCATE (indices(kjpindex),stat=ier)
356	IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo)
357	CALL scatter(indices_g,indices)
358	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
359	IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex)
360
361	!---
362	!--- initialize global grid
363	!---
364	CALL init_grid( kjpindex )
365	CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g)
366
367	!---
368	!--- initialize local grid
369	!---
370	jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
371	if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
372	ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
373	IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
374	ALLOCATE(lon(iim,jjm))
375	ALLOCATE(lat(iim,jjm))
376	lon=zero
377	lat=zero
378	CALL scatter2D(lon_g,lon)
379	CALL scatter2D(lat_g,lat)
380
381	DO ji=1,kjpindex
382
383	j = jlandindex(ji)
384	i = ilandindex(ji)
385
386	!- Create the internal coordinate table
387	!-
388	lalo(ji,1) = lat(i,j)
389	lalo(ji,2) = lon(i,j)
390	ENDDO
391	CALL scatter(neighbours_g,neighbours)
392	CALL scatter(resolution_g,resolution)
393	CALL scatter(contfrac_g,contfrac)
394	CALL scatter(area_g,area)
395	!-
396	!- Check if we have by any change a rectilinear grid. This would allow us to
397	!- simplify the output files.
398	!
399	rectilinear = .FALSE.
400	IF (is_root_prc) THEN
401	IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
402	& ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
403	rectilinear = .TRUE.
404	ENDIF
405	ENDIF
406	CALL bcast(rectilinear)
407	IF (rectilinear) THEN
408	ALLOCATE(lon_rect(iim),stat=ier)
409	IF (ier .NE. 0) THEN
410	WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
411	STOP 'intersurf_history'
412	ENDIF
413	ALLOCATE(lat_rect(jjm),stat=ier)
414	IF (ier .NE. 0) THEN
415	WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
416	STOP 'intersurf_history'
417	ENDIF
418	lon_rect(:) = lon(:,1)
419	lat_rect(:) = lat(1,:)
420	ENDIF
421	!-
422	! allocate arrays
423	!-
424	!
425	a_er = .FALSE.
426	ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
427	a_er = a_er .OR. (ier.NE.0)
428	ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
429	a_er = a_er .OR. (ier.NE.0)
430	ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
431	a_er = a_er .OR. (ier.NE.0)
432	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
433	a_er = a_er .OR. (ier.NE.0)
434	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
435	a_er = a_er .OR. (ier.NE.0)
436	ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
437	a_er = a_er .OR. (ier.NE.0)
438	ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
439	a_er = a_er .OR. (ier.NE.0)
440	ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
441	a_er = a_er .OR. (ier.NE.0)
442	ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
443	a_er = a_er .OR. (ier.NE.0)
444	ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
445	a_er = a_er .OR. (ier.NE.0)
446	ALLOCATE (t2m(kjpindex),stat=ier)
447	a_er = a_er .OR. (ier.NE.0)
448	ALLOCATE (t2m_min(kjpindex),stat=ier)
449	a_er = a_er .OR. (ier.NE.0)
450	ALLOCATE (temp_sol(kjpindex),stat=ier)
451	a_er = a_er .OR. (ier.NE.0)
452	ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
453	a_er = a_er .OR. (ier.NE.0)
454	ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
455	a_er = a_er .OR. (ier.NE.0)
456	ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
457	a_er = a_er .OR. (ier.NE.0)
458	ALLOCATE (litterhum(kjpindex),stat=ier)
459	a_er = a_er .OR. (ier.NE.0)
460	ALLOCATE (precip_rain(kjpindex),stat=ier)
461	a_er = a_er .OR. (ier.NE.0)
462	ALLOCATE (precip_snow(kjpindex),stat=ier)
463	a_er = a_er .OR. (ier.NE.0)
464	ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
465	a_er = a_er .OR. (ier.NE.0)
466	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
467	a_er = a_er .OR. (ier.NE.0)
468	ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
469	a_er = a_er .OR. (ier.NE.0)
470	ALLOCATE (height_x(kjpindex,nvm),stat=ier)
471	a_er = a_er .OR. (ier.NE.0)
472	ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
473	a_er = a_er .OR. (ier.NE.0)
474	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
475	a_er = a_er .OR. (ier.NE.0)
476	ALLOCATE (fco2_lu(kjpindex),stat=ier)
477	a_er = a_er .OR. (ier.NE.0)
478	IF (a_er) THEN
479	CALL ipslerr (3,'teststomate', &
480	& 'PROBLEM WITH ALLOCATION', &
481	& 'for local variables 1','')
482	ENDIF
483	!
484	! prepare forcing
485	!
486	max_totsize = 50
487	CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
488	max_totsize = max_totsize * 1000000
489
490	totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
491	SIZE(humrel_x)*KIND(humrel_x) + &
492	SIZE(litterhum)*KIND(litterhum) + &
493	SIZE(t2m)*KIND(t2m) + &
494	SIZE(t2m_min)*KIND(t2m_min) + &
495	SIZE(temp_sol)*KIND(temp_sol) + &
496	SIZE(soiltemp)*KIND(soiltemp) + &
497	SIZE(soilhum)*KIND(soilhum) + &
498	SIZE(precip_rain)*KIND(precip_rain) + &
499	SIZE(precip_snow)*KIND(precip_snow) + &
500	SIZE(gpp_x)*KIND(gpp_x) + &
501	SIZE(veget_force_x)*KIND(veget_force_x) + &
502	SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
503	SIZE(lai_force_x)*KIND(lai_force_x)
504	CALL reduce_sum(totsize_1step,totsize_tmp)
505	CALL bcast(totsize_tmp)
506	totsize_1step=totsize_tmp
507
508	! total number of forcing steps
509	IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
510	CALL ipslerr (3,'teststomate', &
511	& 'stomate: error with total number of forcing steps', &
512	& 'nsft','teststomate computation different with forcing file value.')
513	ENDIF
514	! number of forcing steps in memory
515	nsfm = MIN(nsft, &
516	& MAX(1,NINT( REAL(max_totsize,r_std) &
517	& /REAL(totsize_1step,r_std))))
518	!-
519	WRITE(numout,*) 'Offline forcing of Stomate:'
520	WRITE(numout,*) ' Total number of forcing steps:',nsft
521	WRITE(numout,*) ' Number of forcing steps in memory:',nsfm
522	!-
523	CALL init_forcing(kjpindex,nsfm,nsft)
524	!-
525	! ensure that we read all new forcing states
526	iisf = nsfm
527	! initialize that table that contains the indices
528	! of the forcing states that will be in memory
529	isf(:) = (/ (i,i=1,nsfm) /)
530
531	nf_written(:) = .FALSE.
532	nf_written(isf(:)) = .TRUE.
533
534	!-
535	! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
536	!-
537	itau_step = NINT(dt_force/dtradia)
538	IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
539	!
540	CALL ioconf_startdate(date0)
541	CALL intsurf_time( itau_dep, date0, dtradia )
542	!-
543	! Length of integration
544	!-
545	WRITE(time_str,'(a)') '1Y'
546	CALL getin_p ('TIME_LENGTH', time_str)
547	! transform into itau
548	CALL tlen2itau(time_str, dt, date0, itau_len)
549	! itau_len*dtradia must be a multiple of dt_force
550	itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
551	& *dt_force/dtradia)
552	!-
553	itau_fin = itau_dep+itau_len
554	!-
555	! set up STOMATE history file
556	!-
557	!Config Key = STOMATE_OUTPUT_FILE
558	!Config Desc = Name of file in which STOMATE's output is going
559	!Config to be written
560	!Config Def = stomate_history.nc
561	!Config Help = This file is going to be created by the model
562	!Config and will contain the output from the model.
563	!Config This file is a truly COADS compliant netCDF file.
564	!Config It will be generated by the hist software from
565	!Config the IOIPSL package.
566	!-
567	stom_histname='stomate_history.nc'
568	CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
569	WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
570	!-
571	!Config Key = STOMATE_HIST_DT
572	!Config Desc = STOMATE history time step (d)
573	!Config Def = 10.
574	!Config Help = Time step of the STOMATE history file
575	!-
576	hist_days_stom = 10.
577	CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
578	IF ( hist_days_stom == -1. ) THEN
579	hist_dt_stom = -1.
580	WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
581	ELSE
582	hist_dt_stom = NINT( hist_days_stom ) * one_day
583	WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
584	hist_dt_stom/one_day
585	ENDIF
586	!-
587	!-
588	!- initialize
589	WRITE(numout,*) "before histbeg : ",date0,dt
590	IF ( rectilinear ) THEN
591	#ifdef CPP_PARA
592	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
593	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
594	#else
595	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
596	& itau_dep, date0, dt, hori_id, hist_id_stom)
597	#endif
598	ELSE
599	#ifdef CPP_PARA
600	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
601	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
602	#else
603	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
604	& itau_dep, date0, dt, hori_id, hist_id_stom)
605	#endif
606	ENDIF
607	!- define PFT axis
608	hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
609	!- declare this axis
610	CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
611	& '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
612	!- define Pool_10 axis
613	hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
614	!- declare this axis
615	CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
616	& '1', 10, hist_pool_10axis, hist_pool_10axis_id)
617
618	!- define Pool_100 axis
619	hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
620	!- declare this axis
621	CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
622	& '1', 100, hist_pool_100axis, hist_pool_100axis_id)
623
624	!- define Pool_11 axis
625	hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
626	!- declare this axis
627	CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
628	& '1', 11, hist_pool_11axis, hist_pool_11axis_id)
629	!- define Pool_101 axis
630	hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
631	!- declare this axis
632	CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
633	& '1', 101, hist_pool_101axis, hist_pool_101axis_id)
634
635	!- define STOMATE history file
636	CALL stom_define_history (hist_id_stom, nvm, iim, jjm, &
637	& dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
638	& hist_pool_10axis_id, hist_pool_100axis_id, &
639	& hist_pool_11axis_id, hist_pool_101axis_id)
640
641	!- end definition
642	CALL histend(hist_id_stom)
643	!-
644	!-
645	! STOMATE IPCC OUTPUTS IS ACTIVATED
646	!-
647	!Config Key = STOMATE_IPCC_OUTPUT_FILE
648	!Config Desc = Name of file in which STOMATE's output is going
649	!Config to be written
650	!Config Def = stomate_ipcc_history.nc
651	!Config Help = This file is going to be created by the model
652	!Config and will contain the output from the model.
653	!Config This file is a truly COADS compliant netCDF file.
654	!Config It will be generated by the hist software from
655	!Config the IOIPSL package.
656	!-
657	stom_ipcc_histname='stomate_ipcc_history.nc'
658	CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)
659	WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
660	!-
661	!Config Key = STOMATE_IPCC_HIST_DT
662	!Config Desc = STOMATE IPCC history time step (d)
663	!Config Def = 0.
664	!Config Help = Time step of the STOMATE IPCC history file
665	!-
666	hist_days_stom_ipcc = zero
667	CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)
668	IF ( hist_days_stom_ipcc == moins_un ) THEN
669	hist_dt_stom_ipcc = moins_un
670	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
671	ELSE
672	hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
673	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
674	hist_dt_stom_ipcc/one_day
675	ENDIF
676
677	! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
678	dt_slow_ = one_day
679	CALL getin_p('DT_SLOW', dt_slow_)
680	IF ( hist_days_stom_ipcc > zero ) THEN
681	IF (dt_slow_ > hist_dt_stom_ipcc) THEN
682	WRITE(numout,*) "DT_SLOW = ",dt_slow_," , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
683	CALL ipslerr (3,'intsurf_history', &
684	& 'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
685	& '(must be less or equal)')
686	ENDIF
687	ENDIF
688
689	IF ( hist_dt_stom_ipcc == 0 ) THEN
690	hist_id_stom_ipcc = -1
691	ELSE
692	!-
693	!- initialize
694	IF ( rectilinear ) THEN
695	#ifdef CPP_PARA
696	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
697	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
698	#else
699	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
700	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
701	#endif
702	ELSE
703	#ifdef CPP_PARA
704	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
705	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
706	#else
707	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
708	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
709	#endif
710	ENDIF
711	!- declare this axis
712	CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
713	& '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)
714
715	!- define STOMATE history file
716	CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
717	& dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)
718
719	!- end definition
720	CALL histend(hist_id_stom_IPCC)
721
722	ENDIF
723	!
724	CALL histwrite(hist_id_stom, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
725	IF ( hist_id_stom_IPCC > 0 ) THEN
726	CALL histwrite(hist_id_stom_IPCC, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
727	ENDIF
728	!
729	hist_id_sec = -1
730	hist_id_sec2 = -1
731	!-
732	! first call of slowproc to initialize variables
733	!-
734	itau = itau_dep
735	!
736	DO ji=1,kjpindex
737	DO jv=1,nvm
738	indexveg((jv-1)kjpindex + ji) = indices(ji) + (jv-1)kjpij
739	ENDDO
740	ENDDO
741	!-
742	!MM Problem here with dpu which depends on soil type
743	DO l = 1, nbdl-1
744	! first 2.0 is dpu
745	! second 2.0 is average
746	diaglev(l) = dpu_cste/(2*(nbdl-1) -1) ( ( 2(l-1) -1) + ( 2(l) -1) ) / 2.0
747	ENDDO
748	diaglev(nbdl) = dpu_cste
749	!
750	CALL ioget_expval(val_exp)
751	ldrestart_read = .TRUE.
752	ldrestart_write = .FALSE.
753	!-
754	! read some variables we need from SECHIBA's restart file
755	!-
756	CALL slowproc_main &
757	& (itau, kjpij, kjpindex, dt_force, date0, &
758	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
759	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
760	& t2m, t2m_min, temp_sol, soiltemp, &
761	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
762	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
763	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
764	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
765	! correct date
766	day_counter = one_day - dt_force
767	date=1
768	!-
769	! time loop
770	!-
771	IF (debug) check_time=.TRUE.
772	CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
773	l_first_intersurf=.FALSE.
774	!
775	DO itau = itau_dep+itau_step,itau_fin,itau_step
776	!
777	CALL intsurf_time( itau, date0, dtradia )
778	!
779	!-- next forcing state
780	iisf = iisf+1
781	IF (debug) WRITE(numout,*) "itau,iisf : ",itau,iisf
782	!---
783	IF (iisf .GT. nsfm) THEN
784	!-----
785	!---- we have to read new forcing states from the file
786	!-----
787	!---- determine blocks of forcing states that are contiguous in memory
788	!-----
789	CALL forcing_read(forcing_id,nsfm)
790
791	!--------------------------
792
793	!-----
794	!---- determine which forcing states must be read next time
795	!-----
796	isf(1) = isf(nsfm)+1
797	IF ( isf(1) .GT. nsft ) isf(1) = 1
798	DO iiisf = 2, nsfm
799	isf(iiisf) = isf(iiisf-1)+1
800	IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
801	ENDDO
802	nf_written(isf(:)) = .TRUE.
803	!---- start again at first forcing state
804	iisf = 1
805	ENDIF
806	! Bug here ! soiltype(:,3) != clay
807	! soiltype(:,3) = clay_fm(:,iisf)
808	humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
809	litterhum(:) = litterhum_daily_fm(:,iisf)
810	t2m(:) = t2m_daily_fm(:,iisf)
811	t2m_min(:) = t2m_min_daily_fm(:,iisf)
812	temp_sol(:) = tsurf_daily_fm(:,iisf)
813	soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
814	soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
815	precip_rain(:) = precip_fm(:,iisf)
816	gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
817	veget_force_x(:,:) = veget_fm(:,:,iisf)
818	veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
819	lai_force_x(:,:) = lai_fm(:,:,iisf)
820	WHERE ( t2m(:) .LT. ZeroCelsius )
821	precip_snow(:) = precip_rain(:)
822	precip_rain(:) = 0.
823	ELSEWHERE
824	precip_snow(:) = 0.
825	ENDWHERE
826	!---
827	!-- scale GPP to new lai and veget_max
828	!---
829	WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
830	!-- scale GPP to new LAI
831	WHERE (lai_force_x(:,:) .GT. 0.0 )
832	gpp_x(:,:) = gpp_x(:,:)ATAN(2.lai_x(:,:)) &
833	& /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
834	ENDWHERE
835	!-- scale GPP to new veget_max
836	WHERE (veget_max_force_x(:,:) .GT. 0.0 )
837	gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
838	ENDWHERE
839	!---
840	!-- number crunching
841	!---
842	ldrestart_read = .FALSE.
843	ldrestart_write = .FALSE.
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_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
853	day_counter = one_day - dt_force
854	ENDDO ! end of the time loop
855	!-
856	! write restart files
857	!-
858	IF (is_root_prc) THEN
859	! first, read and write variables that are not managed otherwise
860	taboo_vars = &
861	& '$lat$ $lon$ $lev$ $veget_year$ '// &
862	& '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
863	& '$lai$ $soiltype_frac$ $clay_frac$ '// &
864	& '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
865	!-
866	CALL ioget_vname(rest_id_sec, nbvar, varnames)
867	!-
868	DO iv = 1, nbvar
869	!-- check if the variable is to be written here
870	IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
871	IF (debug) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
872
873	!---- get variable dimensions, especially 3rd dimension
874	CALL ioget_vdim &
875	& (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
876	l1d = ALL(vardims(1:varnbdim) .EQ. 1)
877	!---- read it
878	IF (l1d) THEN
879	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
880	1,1,1,itau_dep,.TRUE.,var_1d)
881	ELSE
882	ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
883	IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
884	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
885	nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
886	"gather",nbp_glo,indices_g)
887	ENDIF
888	!---- write it
889	IF (l1d) THEN
890	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
891	1,1,1,itau_fin,var_1d)
892	ELSE
893	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
894	nbp_glo,vardims(3),1,itau_fin,var_3d, &
895	'scatter',nbp_glo,indices_g)
896	DEALLOCATE (var_3d)
897	ENDIF
898	ENDIF
899	ENDDO
900	ENDIF
901	CALL barrier_para()
902
903	! call slowproc to write restart files
904	ldrestart_read = .FALSE.
905	ldrestart_write = .TRUE.
906	!-
907	IF (debug) WRITE(numout,*) "Call slowproc for restart."
908	CALL slowproc_main &
909	& (itau_fin, kjpij, kjpindex, dt_force, date0, &
910	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
911	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
912	& t2m, t2m_min, temp_sol, soiltemp, &
913	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
914	& deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
915	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
916	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
917	!-
918	! close files
919	!-
920	IF (is_root_prc) THEN
921	CALL restclo
922	IF ( debug ) WRITE(numout,*) 'REST CLOSED'
923	ENDIF
924	CALL histclo
925
926	IF (is_root_prc) &
927	ier = NF90_CLOSE (forcing_id)
928
929	IF (is_root_prc) THEN
930	CALL getin_dump()
931	ENDIF
932	#ifdef CPP_PARA
933	CALL MPI_FINALIZE(ier)
934	#endif
935	WRITE(numout,*) "End of teststomate."
936
937	!---------------
938	END PROGRAM teststomate
939	!
940	!===
941	!

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

Download in other formats: