Context Navigation

teststomate.f90 @ 5143

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_driver/teststomate.f90 @ 5143

Download in other formats: