Context Navigation

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

Last change on this file since 367 was 367, checked in by didier.solyga, 13 years ago
Synchronize the externalized version with revisions 353, 355 and 356 of the trunk
File size: 34.2 KB

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

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

Download in other formats: