Context Navigation

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

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

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

Download in other formats: