Context Navigation

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

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

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

Download in other formats: