Context Navigation

teststomate.f90 @ 6715

Last change on this file since 6715 was 4646, checked in by josefine.ghattas, 7 years ago

Ticket #185 :
Created new module containing time variables and subroutines to calculate them. Variables in the time module are public and can be used elswhere in the model. They should not be modified outside the time module.

Note that now each time step have a time interval. xx_start or xx_end values for the date can be used. Previously the date values correspondend to the end of the interval. See description in module time for more information.

These changes do not make any difference in results for couled mode with LMDZ or dim2_driver offline driver.

These modifications also corrects the problem with do_slow related to orchideedriver, ticket #327

File size: 49.7 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_driver/teststomate.f90 @ 6715

Download in other formats: