Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_driver/teststomate.f90 @ 5242

Last change on this file since 5242 was 2091, checked in by matthew.mcgrath, 11 years ago
DEV: Fixed a bug for teststomate
File size: 55.9 KB

Line
1	! =================================================================================================================================
2	! PROGRAM : teststomate
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006). This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!
8	!>\BRIEF This program runs the STOMATE submodel using specific initial conditions
9	!! and driving variables in order to obtain the state variables of STOMATE closed to the steady-state values
10	!! quicker than when using the 'full' ORCHIDEE.
11	!!
12	!! \n DESCRIPTION (functional, design, flags):
13	!! It integrates STOMATE for a given number of years using a specific forcing file.
14	!! The initial GPP from this forcing is scaled in agreement with the updated calculated LAI in STOMATE
15	!! Nothing is done for soil moisture which should actually evolve with the vegetation.
16	!! - 1. It first reads a set of parameters, allocates variables and initializes them.
17	!! - 2. Then, a section on the set up of the STOMATE history file
18	!! - 3. A first call to slowproc subroutine in order to initialize SECHIBA variables
19	!! - 4. A loop over time in which slowproc is called at each time step
20	!! - 5. A restart file is written
21	!!
22	!! RECENT CHANGE(S) : None
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!! SVN :
27	!! $HeadURL$
28	!! $Date$
29	!! $Revision$
30	!! \n
31	!_ ================================================================================================================================
32
33	PROGRAM teststomate
34
35	! modules used;
36
37	USE netcdf
38	USE defprec
39	USE constantes
40	USE constantes_soil
41	USE pft_parameters
42	USE structures
43	USE stomate_data
44	USE ioipsl_para
45	USE mod_orchidee_para
46	use tools_para
47	USE grid
48	USE slowproc
49	USE stomate
50	USE intersurf, ONLY: stom_define_history, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time, impose_param
51	! USE parallel
52
53	IMPLICIT NONE
54
55	!! 0. Variable and parameter declaration
56
57	INTEGER(i_std) :: kjpij,kjpindex !! # of total/land points
58	REAL(r_std) :: dtradia,dt_force !! time step of sechiba and stomate components
59	!! read into the forcing file [second]
60
61	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices !! index over land points
62	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indexveg
63	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltile, veget_x !! texture and fraction of vegetation type
64	!! (accounts for LAI dynamic)
65	REAL(r_std),DIMENSION(:),ALLOCATABLE :: totfrac_nobio !! Total fraction of ice+lakes+cities etc. in
66	!! the mesh
67	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: frac_nobio !! Fraction of ice, lakes, cities etc. in the
68	!! mesh
69	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_x !! Fraction of vegetation type
70	INTEGER(i_std), DIMENSION(:),ALLOCATABLE :: njsc
71	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_x !! Leaf Area Index as calculated in STOMATE
72	!! [m2/m2]
73	REAL(r_std),DIMENSION(:),ALLOCATABLE :: sinang_x !! Sine of the sun angle
74	REAL(r_std), DIMENSION(:,:),ALLOCATABLE :: ind_x !! density of individuals
75	!! @tex $(m^{-2})$ @endtex
76	REAL(r_std), DIMENSION(:,:,:,:,:), ALLOCATABLE :: circ_class_biomass_x !! Stem diameter
77	!! @tex $(m)$ @endtex
78	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: circ_class_n_x !! Number of trees within each circumference
79	!! class @tex $(m^{-2})$ @endtex
80
81	REAL(r_std), DIMENSION(:,:,:,:), ALLOCATABLE :: biomass_x !! Total biomass
82
83
84	REAL(r_std), DIMENSION(:,:,:), ALLOCATABLE :: lai_per_level_x !!
85	TYPE(laieff_type), ALLOCATABLE, DIMENSION (:,:,:) :: laieff_fit_x !! The parameters for fitting the effective
86	!! LAI function
87	!$OMP THREADPRIVATE(laieff_fit)
88
89	REAL(r_std),DIMENSION (:,:,:), ALLOCATABLE:: frac_age_x !! Age efficacity from STOMATE for isoprene
90	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_force_x !! Fraction of vegetation of PFTs in each grid
91	!! cell (read in the forcinng file) [-]
92	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_force_x !! fraction of maximum vegetation of PFTs in
93	!! each grid cell (read in the forcinng file)
94	!! [-]
95	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_force_x !! Leaf Area Index as read in the forcing file
96	!! [m2/m2]
97	REAL(r_std),DIMENSION(:),ALLOCATABLE :: reinf_slope
98	REAL(r_std),DIMENSION(:),ALLOCATABLE :: t2m,t2m_min,temp_sol !! 2 m air temperature, min. 2 m air temp.
99	!! during forcing time step and Surface
100	!! temperature [K]
101	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum !! Soil temperature and Relative soil moisture
102	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: humrel_x
103	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: transpir_x !! transpiration
104	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: transpir_mod_x !! transpiration divided by veget_max (per m2)
105	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: pot_transpir_mod_x !! potential transpiration divided by veget_max (per m2)
106	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: transpir_supply_x !! Water supply to the leaves for transpiration
107	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: vir_transpir_supply_x !! Virtual supply of water for transpiration to deal
108	!! with water stress when PFT1 becomes vegetated in LCC
109	!! @tex $(mm dt^{-1})$ @endtex
110	REAL(r_std),DIMENSION(:),ALLOCATABLE :: litterhum !! Litter humidity
111	REAL(r_std),DIMENSION(:),ALLOCATABLE :: precip_rain,precip_snow !! Rainfall, Snowfall
112	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: gpp_x !! GPP (gC/(m**2 of total ground)/time step)
113	REAL(r_std),DIMENSION(:),ALLOCATABLE :: deadleaf_cover !! fraction of soil covered by dead leaves
114	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: assim_param_x !! min/max/opt temperature for photosynthesis
115	!! and vcmax/vjmax parameters
116	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: height_x !! height (m)
117	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x !! Maximum water on vegetation for interception
118	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux !! CO2 flux in gC/m**2 of ground/second
119	REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu !! Land Cover Change CO2 flux (gC/m**2 of
120	REAL(r_std),DIMENSION(:),ALLOCATABLE :: temp_growth !! Growth temperature - Is equal to t2m_month (Â°C)
121	!! average ground/s)
122	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: Isotrop_Tran_Tot_p_x !! Transmitted radiation per level (unitless)
123
124	INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g !! Vector of indeces of land points
125	REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g !! Vector of indeces of land points
126	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g !! Indeces of neighbours for each land point
127
128	INTEGER :: ier,iret !! Return value of functions used to catch
129	!! errors
130	INTEGER :: ncfid !! Id of the restart file of SECHIBA used in
131	!! input
132	CHARACTER(LEN=20),SAVE :: thecalendar='noleap' !! Type of calendar used
133	LOGICAL :: a_er !! Flag used to check errors when allocating
134	CHARACTER(LEN=80) :: & !! Name of the restart files used for
135	& dri_restname_in,dri_restname_out, & !! the driver in/out
136	& sec_restname_in,sec_restname_out, & !! the sechiba component in/out
137	& sto_restname_in,sto_restname_out, & !! the stomate component in/out
138	& stom_histname, stom_ipcc_histname !! Name of the history files for stomate
139	!! (standard and ippc format)
140	INTEGER(i_std) :: iim,jjm,llm !! # of lon,lat and time-step in the restart
141	!! files
142	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat !! Arrays of Longitude and latitude (iim,jjm)
143	!! dimension for each processor
144	REAL, ALLOCATABLE, DIMENSION(:) :: lev !! height level (required by the restini
145	!! interface)
146	LOGICAL :: rectilinear !! Is the grid rectilinear
147	REAL, ALLOCATABLE, DIMENSION(:) :: lon_rect, lat_rect !! Vectors of Longitude and latitude (iim or
148	!! jjm) dimension used in case of a
149	!! rectilinear grid
150	REAL(r_std) :: dt !! Time step of sechiba component read in the
151	!! restart file [seconds]
152	INTEGER(i_std) :: itau_dep,itau_fin,itau !! First, last and current time step of SECHIBA
153	!! component
154	INTEGER(i_std) :: itau_len,itau_step !! Total length of the simulation and length
155	!! between 2 calls of slowproc (expreseed in
156	!! time steps of SECHIBA component)
157	REAL(r_std) :: date0 !! Initial date
158	INTEGER(i_std) :: rest_id_sec,rest_id_sto !! ID's of the restart files for the SECHIBA
159	!! and STOMATE components
160	INTEGER(i_std) :: hist_id_sec,hist_id_sec2 !! ID's of the history files of SECHIBA
161	!! component (required by the interface of
162	!! slowproc but not used)
163	INTEGER(i_std) :: hist_id_stom,hist_id_stom_IPCC !! ID's of the history files of STOMATE
164	!! component
165	CHARACTER(LEN=30) :: time_str !! String used for reading the length of
166	!! simulation in the .def file
167	REAL(r_std) :: dt_slow_ !!
168	REAL :: hist_days_stom,hist_days_stom_ipcc !! Time frequency at which variables are
169	!! written in the history file for the STOMATE
170	!! component (standard and ipcc format) [day]
171	REAL :: hist_dt_stom,hist_dt_stom_ipcc !! Time frequency at which variables are
172	!! written in the history file for the STOMATE
173	!! component (standard and ipcc format)
174	!! [second]
175	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: hist_PFTaxis !! Vector with PFT indeces used as axis in the
176	!! history file
177	REAL(r_std), DIMENSION(nshort) :: hist_pool_s_axis !! Deforestation axis
178	REAL(r_std), DIMENSION(nmedium) :: hist_pool_m_axis !! Deforestation axis
179	REAL(r_std), DIMENSION(nlong) :: hist_pool_l_axis !! Deforestation axis
180	REAL(r_std), DIMENSION(nshort+1) :: hist_pool_ss_axis !! Deforestation axis
181	REAL(r_std), DIMENSION(nmedium+1) :: hist_pool_mm_axis !! Deforestation axis
182	REAL(r_std), DIMENSION(nlong+1) :: hist_pool_ll_axis !! Deforestation axis
183	! REAL(r_std),DIMENSION(10) :: hist_pool_10axis !! Vector with 10-year indeces used as axis in
184	!! the history file (for land-use change)
185	! REAL(r_std),DIMENSION(100) :: hist_pool_100axis !! Vector with 100-year indeces used as axis in
186	!! the history file (for land-use change)
187	! REAL(r_std),DIMENSION(11) :: hist_pool_11axis !! Vector with 11-year indeces used as axis in
188	!! the history file (for land-use change)
189	! REAL(r_std),DIMENSION(101) :: hist_pool_101axis !! Vector with 101-year indeces used as axis in
190	!! the history file (for land-use change)
191	INTEGER :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id !! Id of the axis for PFT in the standard/IPCC
192	!! format
193	INTEGER :: hori_id
194	INTEGER(i_std) :: hist_pool_s_axis_id
195	INTEGER(i_std) :: hist_pool_m_axis_id
196	INTEGER(i_std) :: hist_pool_l_axis_id
197	INTEGER(i_std) :: hist_pool_ss_axis_id
198	INTEGER(i_std) :: hist_pool_mm_axis_id
199	INTEGER(i_std) :: hist_pool_ll_axis_id
200	INTEGER(i_std) :: canx_id
201	INTEGER(i_std) :: canx_tot_stom_id
202	! INTEGER :: hist_pool_10axis_id !! Id of the axis for 10-year vector
203	! INTEGER :: hist_pool_100axis_id !! Id of the axis for the 100-year vector
204	! INTEGER :: hist_pool_11axis_id !! Id of the axis for the 11-year vector
205	! INTEGER :: hist_pool_101axis_id !! Id of the axis for the 101-year vector
206	INTEGER(i_std) :: i,j,iv !! used as counters
207	LOGICAL :: ldrestart_read,ldrestart_write !! Flags to activate reading/writing of the
208	!! restart file
209	LOGICAL :: l1d !! Are vector elements 1-dimension size ?
210	INTEGER(i_std),PARAMETER :: nbvarmax=200 !! Maximum number of variables assumed in the
211	!! restart file of SECHIBA component used as
212	!! input
213	INTEGER(i_std) :: nbvar !! Number of variables present in the restart
214	!! file of SECHIBA component used as input
215	CHARACTER(LEN=50),DIMENSION(nbvarmax) :: varnames !! Name of the variables present in the restart
216	!! file of SECHIBA component used as input
217	INTEGER(i_std) :: varnbdim !! Number of dimensions of a variable
218	!! varnames(i)
219	INTEGER(i_std),PARAMETER :: varnbdim_max=20 !! Maximum number of dimensions of a variable
220	!! varnames(i)
221	INTEGER,DIMENSION(varnbdim_max) :: vardims
222	CHARACTER(LEN=200) :: taboo_vars
223	REAL(r_std),DIMENSION(1) :: xtmp
224	INTEGER :: nsfm,nsft
225	INTEGER :: iisf,iiisf
226	INTEGER(i_std) :: max_totsize,totsize_1step
227	INTEGER(i_std) :: totsize_tmp
228	INTEGER :: vid
229	CHARACTER(LEN=100) :: forcing_name
230	REAL :: x
231
232	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d
233	REAL(r_std) :: var_1d(1)
234
235	!-
236	REAL(r_std) :: time_sec,time_step_sec
237	REAL(r_std) :: time_dri,time_step_dri
238	REAL(r_std),DIMENSION(1) :: r1d
239	REAL(r_std) :: julian,djulian
240
241	INTEGER(i_std) :: ji,jv,l
242
243	LOGICAL :: debug
244
245	REAL(r_std), DIMENSION(:),ALLOCATABLE :: z_levels_tot !! A dummy array to label canopy levels
246	INTEGER(i_std) :: ilev,ivm !! Index
247
248	!---------------------------------------------------------------------
249
250	!-
251	! 1. Reading parameters, Allocating variables and Initializations
252	!-
253
254	CALL Init_orchidee_para
255	!ANNE CALL init_para(.FALSE.)
256	CALL init_timer
257
258	IF (is_root_prc) THEN
259	!-
260	! open STOMATE's forcing file to read some basic info
261	!-
262	forcing_name = 'stomate_forcing.nc'
263	CALL getin ('STOMATE_FORCING_NAME',forcing_name)
264	iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id)
265	IF (iret /= NF90_NOERR) THEN
266	CALL ipslerr (3,'teststomate', &
267	& 'Could not open file : ', &
268	& forcing_name,'(Do you have forget it ?)')
269	ENDIF
270	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia)
271	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force)
272	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x)
273	nsft = NINT(x)
274	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x)
275	kjpij = NINT(x)
276	ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x)
277	nbp_glo = NINT(x)
278	ENDIF
279
280	CALL bcast(dtradia)
281	CALL bcast(dt_force)
282	CALL bcast(nsft)
283	CALL bcast(nbp_glo)
284	!-
285	write(numout,*) 'ATTENTION',dtradia,dt_force
286	!-
287	! read info about land points
288	!-
289	IF (is_root_prc) THEN
290	a_er=.FALSE.
291	ALLOCATE (indices_g(nbp_glo),stat=ier)
292	a_er = a_er .OR. (ier.NE.0)
293	IF (a_er) THEN
294	CALL ipslerr (3,'teststomate', &
295	& 'PROBLEM WITH ALLOCATION', &
296	& 'for local variables 1','')
297	ENDIF
298	!
299	ALLOCATE (x_indices_g(nbp_glo),stat=ier)
300	a_er = a_er .OR. (ier.NE.0)
301	IF (a_er) THEN
302	CALL ipslerr (3,'teststomate', &
303	& 'PROBLEM WITH ALLOCATION', &
304	& 'for global variables 1','')
305	ENDIF
306	ier = NF90_INQ_VARID (forcing_id,'index',vid)
307	IF (ier .NE. 0) THEN
308	CALL ipslerr (3,'teststomate', &
309	& 'PROBLEM WITH ALLOCATION', &
310	& 'for global variables 1','')
311	ENDIF
312	ier = NF90_GET_VAR (forcing_id,vid,x_indices_g)
313	IF (iret /= NF90_NOERR) THEN
314	CALL ipslerr (3,'teststomate', &
315	& 'PROBLEM WITH variable "index" in file ', &
316	& forcing_name,'(check this file)')
317	ENDIF
318	indices_g(:) = NINT(x_indices_g(:))
319	DEALLOCATE (x_indices_g)
320	ELSE
321	ALLOCATE (indices_g(0))
322	ENDIF
323
324	!---------------------------------------------------------------------
325	!-
326	! set debug to have more information
327	!-
328	!Config Key = DEBUG_INFO
329	!Config Desc = Flag for debug information
330	!Config If = [-]
331	!Config Def = n
332	!Config Help = This option allows to switch on the output of debug
333	!Config information without recompiling the code.
334	!Config Units = [FLAG]
335	!-
336	debug = .FALSE.
337	CALL getin_p('DEBUG_INFO',debug)
338	!
339	!Config Key = LONGPRINT
340	!Config Desc = ORCHIDEE will print more messages
341	!Config If = [-]
342	!Config Def = n
343	!Config Help = This flag permits to print more debug messages in the run.
344	!Config Units = [FLAG]
345	!
346	long_print = .FALSE.
347	CALL getin_p('LONGPRINT',long_print)
348	!-
349	! activate CO2, STOMATE, but not sechiba
350	!-
351	control%river_routing = .FALSE.
352	control%hydrol_cwrr = .FALSE.
353	control%ok_sechiba = .FALSE.
354	!
355	control%stomate_watchout = .TRUE.
356	control%ok_co2 = .TRUE.
357	control%ok_stomate = .TRUE.
358	control%forest_management = .FALSE.
359	!-
360	! is DGVM activated?
361	!-
362	control%ok_dgvm = .FALSE.
363	CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm)
364	WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm
365
366	!-
367	! Configuration
368	!-
369	! 1. Number of PFTs defined by the user
370
371	!Config Key = NVM
372	!Config Desc = number of PFTs
373	!Config If = OK_SECHIBA or OK_STOMATE
374	!Config Def = 13
375	!Config Help = The number of vegetation types define by the user
376	!Config Units = [-]
377	!
378	CALL getin_p('NVM',nvm)
379	WRITE(numout,*)'the number of pfts is : ', nvm
380
381	! 2. Should we read the parameters in the run.def file ?
382
383	!Config Key = IMPOSE_PARAM
384	!Config Desc = Do you impose the values of the parameters?
385	!Config if = OK_SECHIBA or OK_STOMATE
386	!Config Def = y
387	!Config Help = This flag can deactivate the reading of some parameters.
388	!Config Useful if you want to use the standard values without commenting the run.def
389	!Config Units = [FLAG]
390	!
391	CALL getin_p('IMPOSE_PARAM',impose_param)
392
393	! 3. Allocate and intialize the pft parameters
394
395	CALL pft_parameters_main(control)
396
397	! 4. Activation sub-models of ORCHIDEE
398
399	CALL activate_sub_models(control)
400
401	! 5. Vegetation configuration (impose_veg, land_use, lcchange...previously in slowproc)
402
403	CALL veget_config
404
405	! 6. Read the parameters in the run.def file
406
407	IF (impose_param) THEN
408	CALL config_pft_parameters
409	CALL config_stomate_pft_parameters
410	CALL config_co2_parameters
411	CALL config_stomate_parameters
412	ENDIF
413	!-
414	IF (control%ok_dgvm) THEN
415	IF ( impose_param ) THEN
416	CALL config_dgvm_parameters
417	ENDIF
418	ENDIF
419	!-
420	!-
421	! restart files
422	!-
423	IF (is_root_prc) THEN
424	! Sechiba's restart files
425	sec_restname_in = 'sechiba_start.nc'
426	CALL getin('SECHIBA_restart_in',sec_restname_in)
427	WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
428	IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
429	STOP 'Need a restart file for Sechiba'
430	ENDIF
431	sec_restname_out = 'sechiba_rest_out.nc'
432	CALL getin('SECHIBA_rest_out',sec_restname_out)
433	WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
434	! Stomate's restart files
435	sto_restname_in = 'stomate_start.nc'
436	CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
437	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
438	sto_restname_out = 'stomate_rest_out.nc'
439	CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
440	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
441
442	!-
443	! We need to know iim, jjm.
444	! Get them from the restart files themselves.
445	!-
446	iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
447	IF (iret /= NF90_NOERR) THEN
448	CALL ipslerr (3,'teststomate', &
449	& 'Could not open file : ', &
450	& sec_restname_in,'(Do you have forget it ?)')
451	ENDIF
452	iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
453	iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
454	iret = NF90_INQ_VARID (ncfid, "time", iv)
455	iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
456	iret = NF90_CLOSE (ncfid)
457	i=INDEX(thecalendar,ACHAR(0))
458	IF ( i > 0 ) THEN
459	thecalendar(i:20)=' '
460	ENDIF
461	ENDIF
462	CALL bcast(iim_g)
463	CALL bcast(jjm_g)
464	CALL bcast(thecalendar)
465	!-
466	! calendar
467	!-
468	CALL ioconf_calendar (thecalendar)
469	CALL ioget_calendar (one_year,one_day)
470	!
471	! Parallelization :
472	!
473	! CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
474	call init_orchidee_data_para_driver(nbp_glo, indices_g)
475	kjpindex=nbp_loc
476	jjm=jj_nb
477	iim=iim_g
478	kjpij=iim*jjm
479	IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
480	!-
481	!-
482	! read info about grids
483	!-
484	!-
485	llm=1
486	ALLOCATE(lev(llm))
487	IF (is_root_prc) THEN
488	!-
489	ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
490	ier = NF90_GET_VAR (forcing_id,vid,lalo_g)
491	!-
492	ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
493	ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
494	ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g)
495	neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
496	DEALLOCATE (x_neighbours_g)
497	!-
498	ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
499	ier = NF90_GET_VAR (forcing_id,vid,resolution_g)
500	!-
501	ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
502	ier = NF90_GET_VAR (forcing_id,vid,contfrac_g)
503
504	lon_g(:,:) = zero
505	lat_g(:,:) = zero
506	lev(1) = zero
507	!-
508	CALL restini &
509	& (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
510	& sec_restname_out, itau_dep, date0, dt, rest_id_sec)
511	!-
512	IF ( dt .NE. dtradia ) THEN
513	WRITE(numout,*) 'dt',dt
514	WRITE(numout,*) 'dtradia',dtradia
515	CALL ipslerr (3,'teststomate', &
516	& 'PROBLEM with time steps.', &
517	& sec_restname_in,'(dt .NE. dtradia)')
518	ENDIF
519	!-
520	CALL restini &
521	& (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
522	& sto_restname_out, itau_dep, date0, dt, rest_id_sto)
523	!-
524	IF ( dt .NE. dtradia ) THEN
525	WRITE(numout,*) 'dt',dt
526	WRITE(numout,*) 'dtradia',dtradia
527	CALL ipslerr (3,'teststomate', &
528	& 'PROBLEM with time steps.', &
529	& sto_restname_in,'(dt .NE. dtradia)')
530	ENDIF
531	ENDIF
532	CALL bcast(rest_id_sec)
533	CALL bcast(rest_id_sto)
534	CALL bcast(itau_dep)
535	CALL bcast(date0)
536	CALL bcast(dt)
537	CALL bcast(lev)
538	!---
539	!--- Create the index table
540	!---
541	!--- This job return a LOCAL kindex
542	!---
543	ALLOCATE (indices(kjpindex),stat=ier)
544	IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo)
545	CALL scatter(indices_g,indices)
546	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
547	IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex)
548
549	!---
550	!--- initialize global grid
551	!---
552	CALL init_grid( kjpindex )
553	CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g)
554
555	!---
556	!--- initialize local grid
557	!---
558	jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
559	if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
560	ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
561	IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
562	ALLOCATE(lon(iim,jjm))
563	ALLOCATE(lat(iim,jjm))
564	lon=zero
565	lat=zero
566	CALL scatter2D_mpi(lon_g,lon)
567	CALL scatter2D_mpi(lat_g,lat)
568
569	DO ji=1,kjpindex
570
571	j = jlandindex(ji)
572	i = ilandindex(ji)
573
574	!- Create the internal coordinate table
575	!-
576	lalo(ji,1) = lat(i,j)
577	lalo(ji,2) = lon(i,j)
578	ENDDO
579	CALL scatter(neighbours_g,neighbours)
580	CALL scatter(resolution_g,resolution)
581	CALL scatter(contfrac_g,contfrac)
582	CALL scatter(area_g,area)
583	!-
584	!- Check if we have by any chance a rectilinear grid. This would allow us to
585	!- simplify the output files.
586	!
587	rectilinear = .FALSE.
588	IF (is_root_prc) THEN
589	IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
590	& ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
591	rectilinear = .TRUE.
592	ENDIF
593	ENDIF
594	CALL bcast(rectilinear)
595	IF (rectilinear) THEN
596	ALLOCATE(lon_rect(iim),stat=ier)
597	IF (ier .NE. 0) THEN
598	WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
599	STOP 'intersurf_history'
600	ENDIF
601	ALLOCATE(lat_rect(jjm),stat=ier)
602	IF (ier .NE. 0) THEN
603	WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
604	STOP 'intersurf_history'
605	ENDIF
606	lon_rect(:) = lon(:,1)
607	lat_rect(:) = lat(1,:)
608	ENDIF
609	!-
610	! allocate arrays
611	!-
612	!
613	a_er = .FALSE.
614	ALLOCATE (hist_PFTaxis(nvm),stat=ier)
615	a_er = a_er .OR. (ier.NE.0)
616	ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
617	a_er = a_er .OR. (ier.NE.0)
618	ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
619	a_er = a_er .OR. (ier.NE.0)
620	ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
621	a_er = a_er .OR. (ier.NE.0)
622	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
623	a_er = a_er .OR. (ier.NE.0)
624	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
625	a_er = a_er .OR. (ier.NE.0)
626	ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
627	a_er = a_er .OR. (ier.NE.0)
628	ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
629	a_er = a_er .OR. (ier.NE.0)
630
631	! the following are variables needed for the effective LAI calculation in stomate, but they do not
632	! change any of the pools, so their values don't matter in teststomate
633	ALLOCATE (sinang_x(kjpindex),stat=ier)
634	a_er = a_er .OR. (ier.NE.0)
635	sinang_x=0.5_r_std
636	ALLOCATE (circ_class_n_x(kjpindex,nvm,ncirc),stat=ier)
637	a_er = a_er .OR. (ier.NE.0)
638	ALLOCATE (circ_class_biomass_x(kjpindex,nvm,ncirc,nparts,nelements),stat=ier)
639	a_er = a_er .OR. (ier.NE.0)
640	! end effective LAI variables
641
642	ALLOCATE (ind_x(kjpindex,nvm),stat=ier)
643	a_er = a_er .OR. (ier.NE.0)
644	ALLOCATE (biomass_x(kjpindex,nvm,nparts,nelements),stat=ier)
645	a_er = a_er .OR. (ier.NE.0)
646	ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
647	a_er = a_er .OR. (ier.NE.0)
648	ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
649	a_er = a_er .OR. (ier.NE.0)
650	ALLOCATE (njsc(kjpindex),stat=ier)
651	a_er = a_er .OR. (ier.NE.0)
652	ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
653	a_er = a_er .OR. (ier.NE.0)
654	ALLOCATE (reinf_slope(kjpindex),stat=ier)
655	a_er = a_er .OR. (ier.NE.0)
656	ALLOCATE (t2m(kjpindex),stat=ier)
657	a_er = a_er .OR. (ier.NE.0)
658	ALLOCATE (t2m_min(kjpindex),stat=ier)
659	a_er = a_er .OR. (ier.NE.0)
660	ALLOCATE (temp_sol(kjpindex),stat=ier)
661	a_er = a_er .OR. (ier.NE.0)
662	ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
663	a_er = a_er .OR. (ier.NE.0)
664	ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
665	a_er = a_er .OR. (ier.NE.0)
666	ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
667	a_er = a_er .OR. (ier.NE.0)
668	ALLOCATE (transpir_x(kjpindex,nvm),stat=ier)
669	a_er = a_er .OR. (ier.NE.0)
670	ALLOCATE (transpir_mod_x(kjpindex,nvm),stat=ier)
671	a_er = a_er .OR. (ier.NE.0)
672	ALLOCATE (pot_transpir_mod_x(kjpindex,nvm),stat=ier)
673	a_er = a_er .OR. (ier.NE.0)
674	ALLOCATE (transpir_supply_x(kjpindex,nvm),stat=ier)
675	a_er = a_er .OR. (ier.NE.0)
676	ALLOCATE (vir_transpir_supply_x(kjpindex,nvm),stat=ier)
677	a_er = a_er .OR. (ier.NE.0)
678	ALLOCATE (litterhum(kjpindex),stat=ier)
679	a_er = a_er .OR. (ier.NE.0)
680	ALLOCATE (precip_rain(kjpindex),stat=ier)
681	a_er = a_er .OR. (ier.NE.0)
682	ALLOCATE (precip_snow(kjpindex),stat=ier)
683	a_er = a_er .OR. (ier.NE.0)
684	ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
685	a_er = a_er .OR. (ier.NE.0)
686	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
687	a_er = a_er .OR. (ier.NE.0)
688	ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
689	a_er = a_er .OR. (ier.NE.0)
690	ALLOCATE (height_x(kjpindex,nvm),stat=ier)
691	a_er = a_er .OR. (ier.NE.0)
692	ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
693	a_er = a_er .OR. (ier.NE.0)
694	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
695	a_er = a_er .OR. (ier.NE.0)
696	ALLOCATE (fco2_lu(kjpindex),stat=ier)
697	a_er = a_er .OR. (ier.NE.0)
698	ALLOCATE (temp_growth(kjpindex),stat=ier)
699	a_er = a_er .OR. (ier.NE.0)
700	ALLOCATE (frac_age_x(kjpindex,nvm,nleafages),stat=ier)
701	a_er = a_er .OR. (ier.NE.0)
702	IF (a_er) THEN
703	CALL ipslerr (3,'teststomate', &
704	& 'PROBLEM WITH ALLOCATION', &
705	& 'for local variables 1','')
706	ENDIF
707	!
708	! prepare forcing
709	! forcing file of STOMATE is read by block
710	! in order to minimize the reading frequency
711	! Here is done the calculation of the number
712	! of time steps we load in memory
713	!
714	max_totsize = 50
715	CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
716	max_totsize = max_totsize * 1000000
717
718	totsize_1step = SIZE(soiltile(:,3))*KIND(soiltile(:,3)) + &
719	SIZE(humrel_x)*KIND(humrel_x) + &
720	SIZE(litterhum)*KIND(litterhum) + &
721	SIZE(t2m)*KIND(t2m) + &
722	SIZE(t2m_min)*KIND(t2m_min) + &
723	SIZE(temp_sol)*KIND(temp_sol) + &
724	SIZE(soiltemp)*KIND(soiltemp) + &
725	SIZE(soilhum)*KIND(soilhum) + &
726	SIZE(precip_rain)*KIND(precip_rain) + &
727	SIZE(precip_snow)*KIND(precip_snow) + &
728	SIZE(gpp_x)*KIND(gpp_x) + &
729	SIZE(veget_force_x)*KIND(veget_force_x) + &
730	SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
731	SIZE(lai_force_x)*KIND(lai_force_x)
732	CALL reduce_sum(totsize_1step,totsize_tmp)
733	CALL bcast(totsize_tmp)
734	totsize_1step=totsize_tmp
735
736	! check for consistency of the total number of forcing steps
737	IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
738	CALL ipslerr_p (3,'teststomate', &
739	& 'stomate: error with total number of forcing steps', &
740	& 'nsft','teststomate computation different with forcing file value.')
741	ENDIF
742	! number of forcing steps in memory
743	nsfm = MIN(nsft, &
744	& MAX(1,NINT( REAL(max_totsize,r_std) &
745	& /REAL(totsize_1step,r_std))))
746	!-
747	WRITE(numout,*) 'Offline forcing of Stomate:'
748	WRITE(numout,*) ' Total number of forcing steps:',nsft
749	WRITE(numout,*) ' Number of forcing steps in memory:',nsfm
750	!-
751	! init_forcing defined into stomate.f90
752	! allocate and set to zero driving variables of STOMATE
753	! ie variables that are read in the forcing file
754	!-
755	CALL init_forcing(kjpindex,nsfm,nsft)
756	!-
757	! ensure that we read all new forcing states
758	iisf = nsfm
759	! initialize the table that contains the indices
760	! of the forcing states that will be in memory
761	isf(:) = (/ (i,i=1,nsfm) /)
762
763	nf_written(:) = .FALSE.
764	nf_written(isf(:)) = .TRUE.
765
766	!-
767	! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
768	!-
769	itau_step = NINT(dt_force/dtradia)
770	IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
771	!
772	CALL ioconf_startdate(date0)
773	CALL intsurf_time( itau_dep, date0, dtradia )
774	!-
775	! Length of integration
776	!-
777	WRITE(time_str,'(a)') '1Y'
778	CALL getin_p ('TIME_LENGTH', time_str)
779	! transform into itau
780	CALL tlen2itau(time_str, dt, date0, itau_len)
781	! itau_len*dtradia must be a multiple of dt_force
782	itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
783	& *dt_force/dtradia)
784	!-
785	itau_fin = itau_dep+itau_len
786	!-
787	! 2. set up STOMATE history file
788	!-
789	!Config Key = STOMATE_OUTPUT_FILE
790	!Config Desc = Name of file in which STOMATE's output is going to be written
791	!Config If =
792	!Config Def = stomate_history.nc
793	!Config Help = This file is going to be created by the model
794	!Config and will contain the output from the model.
795	!Config This file is a truly COADS compliant netCDF file.
796	!Config It will be generated by the hist software from
797	!Config the IOIPSL package.
798	!Config Units = FILE
799	!-
800	stom_histname='stomate_history.nc'
801	CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
802	WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
803	!-
804	!Config Key = STOMATE_HIST_DT
805	!Config Desc = STOMATE history time step
806	!Config If =
807	!Config Def = 10.
808	!Config Help = Time step of the STOMATE history file
809	!Config Units = days [d]
810	!-
811	hist_days_stom = 10.
812	CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
813	IF ( hist_days_stom == -1. ) THEN
814	hist_dt_stom = -1.
815	WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
816	ELSE
817	hist_dt_stom = NINT( hist_days_stom ) * one_day
818	WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
819	hist_dt_stom/one_day
820	ENDIF
821	!-
822	!- initialize
823	WRITE(numout,*) "before histbeg : ",date0,dt
824	IF ( rectilinear ) THEN
825	#ifdef CPP_PARA
826	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
827	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
828	#else
829	CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
830	& itau_dep, date0, dt, hori_id, hist_id_stom)
831	#endif
832	ELSE
833	#ifdef CPP_PARA
834	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
835	& itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
836	#else
837	CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
838	& itau_dep, date0, dt, hori_id, hist_id_stom)
839	#endif
840	ENDIF
841	!- define PFT axis
842	hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
843	!- declare this axis
844	CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
845	& '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
846
847	!!$!- define Pool_10 axis
848	!!$ hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
849	!!$!- declare this axis
850	!!$ CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
851	!!$ & '1', 10, hist_pool_10axis, hist_pool_10axis_id)
852	!!$
853	!!$!- define Pool_100 axis
854	!!$ hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
855	!!$!- declare this axis
856	!!$ CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
857	!!$ & '1', 100, hist_pool_100axis, hist_pool_100axis_id)
858	!!$
859	!!$!- define Pool_11 axis
860	!!$ hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
861	!!$!- declare this axis
862	!!$ CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
863	!!$ & '1', 11, hist_pool_11axis, hist_pool_11axis_id)
864	!!$
865	!!$!- define Pool_101 axis
866	!!$ hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
867	!!$!- declare this axis
868	!!$ CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
869	!!$ & '1', 101, hist_pool_101axis, hist_pool_101axis_id)
870
871
872	! deforestation
873	!- define Pool_short axis
874	hist_pool_s_axis = (/ ( REAL(i,r_std), i=1,nshort ) /)
875	!- declare this axis
876	CALL histvert (hist_id_stom, 'P_S', 'Short lived pool', &
877	& '1', nshort, hist_pool_s_axis, hist_pool_s_axis_id)
878
879	!- define Pool_medium axis
880	hist_pool_m_axis = (/ ( REAL(i,r_std), i=1,nmedium ) /)
881	!- declare this axis
882	CALL histvert (hist_id_stom, 'P_M', 'Medium lived pool', &
883	& '1', nmedium, hist_pool_m_axis, hist_pool_m_axis_id)
884
885	!- define Pool_long axis
886	hist_pool_l_axis = (/ ( REAL(i,r_std), i=1,nlong ) /)
887	!- declare this axis
888	CALL histvert (hist_id_stom, 'P_L', 'Long lived pool', &
889	& '1', nlong, hist_pool_l_axis, hist_pool_l_axis_id)
890
891	!- define Pool_ss axis
892	hist_pool_ss_axis = (/ ( REAL(i,r_std), i=1,nshort+1 ) /)
893	!- declare this axis
894	CALL histvert (hist_id_stom, 'P_SS', 'short lived pool years + 1', &
895	& '1', nshort+1, hist_pool_ss_axis, hist_pool_ss_axis_id)
896
897	!- define Pool_mm axis
898	hist_pool_mm_axis = (/ ( REAL(i,r_std), i=1,nmedium+1 ) /)
899	!- declare this axis
900	CALL histvert (hist_id_stom, 'P_MM', 'medium lived pool years + 1', &
901	& '1', nmedium+1, hist_pool_mm_axis, hist_pool_mm_axis_id)
902
903	!- define Pool_ll axis
904	hist_pool_ll_axis = (/ ( REAL(i,r_std), i=1,nlong+1 ) /)
905	!- declare this axis
906	CALL histvert (hist_id_stom, 'P_LL', 'long lived pool years + 1', &
907	& '1', nlong+1, hist_pool_ll_axis, hist_pool_ll_axis_id)
908
909	! For the number of canopy levels
910	CALL histvert(hist_id_stom, 'canlevel', 'Canopy levels', 'm', &
911	& nlevels, z_level, canx_id)
912
913	! For the number of canopy levels including all those used in photosynthesis.
914	! Notice we cannot use the real heights here, since they can change every day
915	! as the canopy develops.
916	ALLOCATE(z_levels_tot(nlevels_tot))
917	DO ilev=1,nlevels_tot
918	z_levels_tot(ilev)=REAL(ilev,r_std)
919	ENDDO
920	CALL histvert(hist_id_stom, 'canleveltot', 'Total canopy levels', '-', &
921	& nlevels_tot, z_levels_tot, canx_tot_stom_id)
922	DEALLOCATE(z_levels_tot)
923
924	!- define STOMATE history file
925	CALL stom_define_history (hist_id_stom, nvm, control, iim, jjm, &
926	& dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
927	& hist_pool_s_axis_id, hist_pool_m_axis_id, &
928	& hist_pool_l_axis_id, hist_pool_ss_axis_id, &
929	& hist_pool_mm_axis_id, hist_pool_ll_axis_id, canx_id,&
930	canx_tot_stom_id)
931
932	!- end definition
933	CALL histend(hist_id_stom)
934	!-
935	!-
936	! STOMATE IPCC OUTPUTS IS ACTIVATED
937	!-
938	!Config Key = STOMATE_IPCC_OUTPUT_FILE
939	!Config Desc = Name of file in which STOMATE's output is going to be written
940	!Config If =
941	!Config Def = stomate_ipcc_history.nc
942	!Config Help = This file is going to be created by the model
943	!Config and will contain the output from the model.
944	!Config This file is a truly COADS compliant netCDF file.
945	!Config It will be generated by the hist software from
946	!Config the IOIPSL package.
947	!Config Units = FILE
948	!-
949	stom_ipcc_histname='stomate_ipcc_history.nc'
950	CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)
951	WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
952	!-
953	!Config Key = STOMATE_IPCC_HIST_DT
954	!Config Desc = STOMATE IPCC history time step
955	!Config If =
956	!Config Def = 0.
957	!Config Help = Time step of the STOMATE IPCC history file
958	!Config Units = days [d]
959	!-
960	hist_days_stom_ipcc = zero
961	CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)
962	IF ( hist_days_stom_ipcc == moins_un ) THEN
963	hist_dt_stom_ipcc = moins_un
964	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
965	ELSE
966	hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
967	WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
968	hist_dt_stom_ipcc/one_day
969	ENDIF
970
971	! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
972	dt_slow_ = one_day
973	CALL getin_p('DT_SLOW', dt_slow_)
974	IF ( hist_days_stom_ipcc > zero ) THEN
975	IF (dt_slow_ > hist_dt_stom_ipcc) THEN
976	WRITE(numout,*) "DT_SLOW = ",dt_slow_," , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
977	CALL ipslerr_p (3,'intsurf_history', &
978	& 'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
979	& '(must be less or equal)')
980	ENDIF
981	ENDIF
982
983	IF ( hist_dt_stom_ipcc == 0 ) THEN
984	hist_id_stom_ipcc = -1
985	ELSE
986	!-
987	!- initialize
988	IF ( rectilinear ) THEN
989	#ifdef CPP_PARA
990	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
991	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
992	#else
993	CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, &
994	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
995	#endif
996	ELSE
997	#ifdef CPP_PARA
998	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
999	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
1000	#else
1001	CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, &
1002	& itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
1003	#endif
1004	ENDIF
1005	!- declare this axis
1006	CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
1007	& '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)
1008
1009	!- define STOMATE history file
1010	CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
1011	& dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)
1012
1013	!- end definition
1014	CALL histend(hist_id_stom_IPCC)
1015
1016	ENDIF
1017	!
1018	CALL histwrite_p(hist_id_stom, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
1019	IF ( hist_id_stom_IPCC > 0 ) THEN
1020	CALL histwrite_p(hist_id_stom_IPCC, 'Areas', itau_dep+itau_step, area, kjpindex, indices)
1021	ENDIF
1022	!
1023	hist_id_sec = -1
1024	hist_id_sec2 = -1
1025	!-
1026	! 3. first call of slowproc to initialize variables
1027	!-
1028	itau = itau_dep
1029	!
1030	DO ji=1,kjpindex
1031	DO jv=1,nvm
1032	indexveg((jv-1)kjpindex + ji) = indices(ji) + (jv-1)kjpij
1033	ENDDO
1034	ENDDO
1035	!-
1036	!
1037	!Config key = HYDROL_CWRR
1038	!Config Desc = Allows to switch on the multilayer hydrology of CWRR
1039	!Config If = OK_SECHIBA
1040	!Config Def = n
1041	!Config Help = This flag allows the user to decide if the vertical
1042	!Config hydrology should be treated using the multi-layer
1043	!Config diffusion scheme adapted from CWRR by Patricia de Rosnay.
1044	!Config by default the Choisnel hydrology is used.
1045	!Config Units = [FLAG]
1046	CALL getin_p ("HYDROL_SOIL_DEPTH", dpu_max)
1047	!
1048	!MM Problem here with dpu which depends on soil type
1049	DO l = 1, nbdl-1
1050	! first 2.0 is dpu
1051	! second 2.0 is average
1052	diaglev(l) = dpu_max/(2*(nbdl-1) -1) ( ( 2(l-1) -1) + ( 2(l) -1) ) / 2.0
1053	ENDDO
1054	diaglev(nbdl) = dpu_max
1055	!
1056	!-
1057	! Read the parameters in the "*.def" files
1058	!-
1059	!
1060	!Config Key = CLAYFRACTION_DEFAULT
1061	!Config Desc =
1062	!Config If = OK_SECHIBA
1063	!Config Def = 0.2
1064	!Config Help =
1065	!Config Units = [-]
1066	CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
1067	!
1068	!Config Key = MIN_VEGFRAC
1069	!Config Desc = Minimal fraction of mesh a vegetation type can occupy
1070	!Config If = OK_SECHIBA
1071	!Config Def = 0.001
1072	!Config Help =
1073	!Config Units = [-]
1074	CALL getin_p('MIN_VEGFRAC',min_vegfrac)
1075	!
1076	!Config Key = STEMPDIAG_BID
1077	!Config Desc = only needed for an initial LAI if there is no restart file
1078	!Config If = OK_SECHIBA
1079	!Config Def = 280.
1080	!Config Help =
1081	!Config Units = [K]
1082	CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
1083	!
1084	!-
1085	CALL ioget_expval(val_exp)
1086	ldrestart_read = .TRUE.
1087	ldrestart_write = .FALSE.
1088	!-
1089	! read some variables we need from SECHIBA's restart file
1090	!-
1091	CALL slowproc_main &
1092	& (itau, kjpij, kjpindex, dt_force, date0, &
1093	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1094	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
1095	& t2m, t2m_min, temp_sol, soiltemp, &
1096	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1097	& deadleaf_cover, assim_param_x, lai_x, frac_age_x, height_x, veget_x, &
1098	& frac_nobio, njsc, veget_max_x, totfrac_nobio, qsintmax_x, &
1099	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu, &
1100	sinang_x, circ_class_biomass_x, circ_class_n_x, transpir_x, transpir_mod_x, &
1101	pot_transpir_mod_x, transpir_supply_x, biomass_x, ind_x, lai_per_level_x, laieff_fit_x, &
1102	vir_transpir_supply_x, temp_growth, Isotrop_Tran_Tot_p_x)
1103
1104	! correct date
1105
1106
1107	day_counter = one_day - dt_force
1108	date=1
1109	IF (debug) check_time=.TRUE.
1110	CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
1111	l_first_intersurf=.FALSE.
1112	!-
1113	! 4. Time loop
1114	!â»
1115	DO itau = itau_dep+itau_step,itau_fin,itau_step
1116	!
1117	CALL intsurf_time( itau, date0, dtradia )
1118	!-
1119	! next forcing state
1120	iisf = iisf+1
1121	IF (debug) WRITE(numout,*) "itau,iisf : ",itau,iisf
1122	!---
1123	IF (iisf .GT. nsfm) THEN
1124	!-----
1125	!---- we have to read new forcing states from the file
1126	!-----
1127	!---- determine blocks of forcing states that are contiguous in memory
1128	!-----
1129	CALL forcing_read(forcing_id,nsfm)
1130
1131	!--------------------------
1132
1133	!-----
1134	!---- determine which forcing states must be read next time
1135	!-----
1136	isf(1) = isf(nsfm)+1
1137	IF ( isf(1) .GT. nsft ) isf(1) = 1
1138	DO iiisf = 2, nsfm
1139	isf(iiisf) = isf(iiisf-1)+1
1140	IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
1141	ENDDO
1142	nf_written(isf(:)) = .TRUE.
1143	!---- start again at first forcing state
1144	iisf = 1
1145	ENDIF
1146	humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
1147	litterhum(:) = litterhum_daily_fm(:,iisf)
1148	t2m(:) = t2m_daily_fm(:,iisf)
1149	t2m_min(:) = t2m_min_daily_fm(:,iisf)
1150	temp_sol(:) = tsurf_daily_fm(:,iisf)
1151	soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
1152	soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
1153	precip_rain(:) = precip_fm(:,iisf)
1154	gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
1155	veget_force_x(:,:) = veget_fm(:,:,iisf)
1156	veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
1157	lai_force_x(:,:) = lai_fm(:,:,iisf)
1158	WHERE ( t2m(:) .LT. ZeroCelsius )
1159	precip_snow(:) = precip_rain(:)
1160	precip_rain(:) = zero
1161	ELSEWHERE
1162	precip_snow(:) = zero
1163	ENDWHERE
1164	!---
1165	!-- scale GPP to new lai and veget_max
1166	!---
1167	WHERE ( lai_x(:,:) .EQ. zero ) gpp_x(:,:) = zero
1168	!-- scale GPP to new LAI
1169	WHERE (lai_force_x(:,:) .GT. zero )
1170	gpp_x(:,:) = gpp_x(:,:)ATAN(2.lai_x(:,:)) &
1171	& /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
1172	ENDWHERE
1173	!- scale GPP to new veget_max
1174	WHERE (veget_max_force_x(:,:) .GT. zero )
1175	gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
1176	ENDWHERE
1177	!-
1178	!- number crunching
1179	!-
1180	ldrestart_read = .FALSE.
1181	ldrestart_write = .FALSE.
1182	CALL slowproc_main &
1183	& (itau, kjpij, kjpindex, dt_force, date0, &
1184	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1185	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
1186	& t2m, t2m_min, temp_sol, soiltemp, &
1187	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1188	& deadleaf_cover, assim_param_x, lai_x, frac_age_x, height_x, veget_x, &
1189	& frac_nobio, njsc, veget_max_x, totfrac_nobio, qsintmax_x, &
1190	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu, &
1191	sinang_x, circ_class_biomass_x, circ_class_n_x, transpir_x, transpir_mod_x, &
1192	pot_transpir_mod_x, transpir_supply_x, biomass_x,ind_x, lai_per_level_x, laieff_fit_x, &
1193	vir_transpir_supply_x, temp_growth, Isotrop_Tran_Tot_p_x)
1194
1195	day_counter = one_day - dt_force
1196	ENDDO ! end of the time loop
1197
1198
1199	!-
1200	! 5. write restart files
1201	!-
1202	IF (is_root_prc) THEN
1203	! first, read and write variables that are not managed otherwise
1204	taboo_vars = &
1205	& '$lat$ $lon$ $lev$ $veget_year$ '// &
1206	& '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
1207	& '$lai$ $soiltile_frac$ $clay_frac$ '// &
1208	& '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
1209	!-
1210	CALL ioget_vname(rest_id_sec, nbvar, varnames)
1211	!-
1212	DO iv = 1, nbvar
1213	!-- check if the variable is to be written here
1214	IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
1215	IF (debug) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
1216
1217	!---- get variable dimensions, especially 3rd dimension
1218	CALL ioget_vdim &
1219	& (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
1220	l1d = ALL(vardims(1:varnbdim) .EQ. 1)
1221	!---- read it
1222	IF (l1d) THEN
1223	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
1224	1,1,1,itau_dep,.TRUE.,var_1d)
1225	ELSE
1226	ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
1227	IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
1228	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
1229	nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
1230	"gather",nbp_glo,indices_g)
1231	ENDIF
1232	!---- write it
1233	IF (l1d) THEN
1234	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
1235	1,1,1,itau_fin,var_1d)
1236	ELSE
1237	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
1238	nbp_glo,vardims(3),1,itau_fin,var_3d, &
1239	'scatter',nbp_glo,indices_g)
1240	DEALLOCATE (var_3d)
1241	ENDIF
1242	ENDIF
1243	ENDDO
1244	ENDIF
1245	CALL barrier_para()
1246
1247	! call slowproc to write restart files
1248	ldrestart_read = .FALSE.
1249	ldrestart_write = .TRUE.
1250	!-
1251	IF (debug) WRITE(numout,*) "Call slowproc for restart."
1252	CALL slowproc_main &
1253	& (itau_fin, kjpij, kjpindex, dt_force, date0, &
1254	& ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1255	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
1256	& t2m, t2m_min, temp_sol, soiltemp, &
1257	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1258	& deadleaf_cover, assim_param_x, lai_x, frac_age_x, height_x, veget_x, &
1259	& frac_nobio, njsc, veget_max_x, totfrac_nobio, qsintmax_x, &
1260	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu, &
1261	sinang_x, circ_class_biomass_x, circ_class_n_x, transpir_x, transpir_mod_x, &
1262	pot_transpir_mod_x, transpir_supply_x, biomass_x,ind_x, lai_per_level_x, laieff_fit_x, &
1263	vir_transpir_supply_x, temp_growth, Isotrop_Tran_Tot_p_x)
1264
1265
1266	!-
1267	! close files
1268	!-
1269	IF (is_root_prc) THEN
1270	CALL restclo
1271	IF ( debug ) WRITE(numout,*) 'REST CLOSED'
1272	ENDIF
1273	CALL histclo
1274
1275	IF (is_root_prc) &
1276	ier = NF90_CLOSE (forcing_id)
1277
1278	IF (is_root_prc) THEN
1279	CALL getin_dump()
1280	ENDIF
1281	#ifdef CPP_PARA
1282	CALL MPI_FINALIZE(ier)
1283	#endif
1284	WRITE(numout,*) "End of teststomate."
1285
1286	!---------------
1287	END PROGRAM teststomate
1288	!
1289	!===
1290	!

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

Download in other formats: