Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_driver/teststomate.f90 @ 8398

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

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

Download in other formats: