Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_driver/teststomate.f90 @ 7911

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

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

Download in other formats: