Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/teststomate.f90 @ 7262

Last change on this file since 7262 was 5559, checked in by josefine.ghattas, 6 years ago
Enhancement on grid module: merge variables grid_type and GridType? meaning the same. Only grid_type is used now, this one is choosen because it is the same as used in LMDZ.
File size: 49.2 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	ok_sechiba = .FALSE.
296	ok_stomate = .TRUE.
297
298	! Deactivate writing of stomate_forcing.nc file
299	allow_forcing_write=.FALSE.
300
301	!-
302	! is DGVM activated?
303	!-
304	ok_dgvm = .FALSE.
305	CALL getin_p('STOMATE_OK_DGVM',ok_dgvm)
306	WRITE(numout,*) 'LPJ is activated: ',ok_dgvm
307
308	! Initialize parameter for off-line use : no coupling to atmospheric model
309	OFF_LINE_MODE=.TRUE.
310	!-
311	! Configuration
312	!-
313	! 1. Number of PFTs defined by the user
314
315	!Config Key = NVM
316	!Config Desc = number of PFTs
317	!Config If = OK_SECHIBA or OK_STOMATE
318	!Config Def = 13
319	!Config Help = The number of vegetation types define by the user
320	!Config Units = [-]
321	!
322	CALL getin_p('NVM',nvm)
323	WRITE(numout,*)'the number of pfts is : ', nvm
324
325	! 2. Should we read the parameters in the run.def file ?
326
327	!Config Key = IMPOSE_PARAM
328	!Config Desc = Do you impose the values of the parameters?
329	!Config if = OK_SECHIBA or OK_STOMATE
330	!Config Def = y
331	!Config Help = This flag can deactivate the reading of some parameters.
332	!Config Useful if you want to use the standard values without commenting the run.def
333	!Config Units = [FLAG]
334	!
335	CALL getin_p('IMPOSE_PARAM',impose_param)
336
337	! 3. Allocate and intialize the pft parameters
338
339	CALL pft_parameters_main()
340
341	! 4. Activation sub-models of ORCHIDEE
342
343	CALL activate_sub_models()
344
345	! 5. Vegetation configuration (impose_veg, land_use, lcchange...previously in slowproc)
346
347	CALL veget_config
348
349	! 6. Read the parameters in the run.def file
350
351	IF (impose_param) THEN
352	CALL config_pft_parameters
353	CALL config_stomate_pft_parameters
354	CALL config_co2_parameters
355	CALL config_stomate_parameters
356	ENDIF
357	!-
358	IF (ok_dgvm) THEN
359	IF ( impose_param ) THEN
360	CALL config_dgvm_parameters
361	ENDIF
362	ENDIF
363	!-
364	!-
365	! restart files
366	!-
367	IF (is_root_prc) THEN
368	! Sechiba's restart files
369	sec_restname_in = 'sechiba_start.nc'
370	CALL getin('SECHIBA_restart_in',sec_restname_in)
371	WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
372	IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
373	STOP 'Need a restart file for Sechiba'
374	ENDIF
375	sec_restname_out = 'sechiba_rest_out.nc'
376	CALL getin('SECHIBA_rest_out',sec_restname_out)
377	WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
378	! Stomate's restart files
379	sto_restname_in = 'stomate_start.nc'
380	CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
381	WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
382	sto_restname_out = 'stomate_rest_out.nc'
383	CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
384	WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
385
386	!-
387	! We need to know iim, jjm.
388	! Get them from the restart files themselves.
389	!-
390	iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
391	IF (iret /= NF90_NOERR) THEN
392	CALL ipslerr (3,'teststomate', &
393	& 'Could not open file : ', &
394	& sec_restname_in,'(Do you have forget it ?)')
395	ENDIF
396	iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
397	iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
398	iret = NF90_INQ_VARID (ncfid, "time", iv)
399	iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
400	iret = NF90_CLOSE (ncfid)
401	i=INDEX(thecalendar,ACHAR(0))
402	IF ( i > 0 ) THEN
403	thecalendar(i:20)=' '
404	ENDIF
405	ENDIF
406	CALL bcast(iim_g)
407	CALL bcast(jjm_g)
408	CALL bcast(thecalendar)
409	!-
410	! calendar
411	!-
412	CALL ioconf_calendar (thecalendar)
413	CALL ioget_calendar (one_year,one_day)
414	!
415	! Parallelization :
416	!
417	! CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
418	CALL grid_set_glo(iim_g,jjm_g,nbp_glo)
419	CALL grid_allocate_glo(4)
420
421	! Initialize index_g needed for module grid
422	! index_g is declared in mod_orchidee_para and allocated by grid_allocate_glo
423	! index_g and indices_g are the same but indeces_g only declared on root proc.
424	IF (is_root_prc) THEN
425	index_g(:)=indices_g(:)
426	! Only use index_g from now and on => Deallocate indices_g.
427	DEALLOCATE(indices_g)
428	END IF
429	CALL bcast(index_g)
430
431	CALL init_orchidee_data_para_driver(nbp_glo, index_g)
432	kjpindex=nbp_loc
433	jjm=jj_nb
434	iim=iim_g
435	kjpij=iim*jjm
436	IF (printlev_loc>=3) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
437	!-
438	!-
439	! read info about grids
440	!-
441	!-
442	llm=1
443	ALLOCATE(lev(llm))
444	IF (is_root_prc) THEN
445	!-
446	ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
447	ier = NF90_GET_VAR (forcing_id,vid,lalo_g)
448	!-
449	ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
450	ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
451	ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g)
452	neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
453	DEALLOCATE (x_neighbours_g)
454	!-
455	ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
456	ier = NF90_GET_VAR (forcing_id,vid,resolution_g)
457	!-
458	ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
459	ier = NF90_GET_VAR (forcing_id,vid,contfrac_g)
460
461	lon_g(:,:) = zero
462	lat_g(:,:) = zero
463	lev(1) = zero
464	!-
465	CALL restini &
466	& (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
467	& sec_restname_out, itau_dep, date0, dt, rest_id_sec)
468	!-
469	IF ( dt .NE. dtradia ) THEN
470	WRITE(numout,*) 'dt',dt
471	WRITE(numout,*) 'dtradia',dtradia
472	CALL ipslerr (3,'teststomate', &
473	& 'PROBLEM with time steps.', &
474	& sec_restname_in,'(dt .NE. dtradia)')
475	ENDIF
476	!-
477	CALL restini &
478	& (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
479	& sto_restname_out, itau_dep, date0, dt, rest_id_sto)
480	!-
481	IF ( dt .NE. dtradia ) THEN
482	WRITE(numout,*) 'dt',dt
483	WRITE(numout,*) 'dtradia',dtradia
484	CALL ipslerr (3,'teststomate', &
485	& 'PROBLEM with time steps.', &
486	& sto_restname_in,'(dt .NE. dtradia)')
487	ENDIF
488	ENDIF
489	CALL bcast(rest_id_sec)
490	CALL bcast(rest_id_sto)
491	CALL bcast(itau_dep)
492	CALL bcast(date0)
493	CALL bcast(dt)
494	CALL bcast(lev)
495	!---
496	!--- Create the index table
497	!---
498	!--- This job return a LOCAL kindex
499	!---
500	ALLOCATE (indices(kjpindex),stat=ier)
501	IF (printlev_loc>=3 .AND. is_root_prc) WRITE(numout,*) "index_g = ",index_g(1:nbp_glo)
502	CALL scatter(index_g,indices)
503	indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
504	IF (printlev_loc>=3) WRITE(numout,*) "indices = ",indices(1:kjpindex)
505
506	!---
507	!--- initialize global grid
508	!---
509	CALL grid_init( kjpindex, 4, regular_lonlat, "ForcingGrid" )
510	CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, index_g)
511
512	!---
513	!--- initialize local grid
514	!---
515	jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
516	if (printlev_loc>=3) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
517	ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
518	IF (printlev_loc>=3) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
519	ALLOCATE(lon(iim,jjm))
520	ALLOCATE(lat(iim,jjm))
521	lon=zero
522	lat=zero
523	CALL scatter2D_mpi(lon_g,lon)
524	CALL scatter2D_mpi(lat_g,lat)
525
526	DO ji=1,kjpindex
527
528	j = jlandindex(ji)
529	i = ilandindex(ji)
530
531	!- Create the internal coordinate table
532	!-
533	lalo(ji,1) = lat(i,j)
534	lalo(ji,2) = lon(i,j)
535	ENDDO
536	CALL scatter(neighbours_g,neighbours)
537	CALL scatter(resolution_g,resolution)
538	CALL scatter(contfrac_g,contfrac)
539	CALL scatter(area_g,area)
540	!-
541	!- Check if we have by any chance a rectilinear grid. This would allow us to
542	!- simplify the output files.
543	!
544	rectilinear = .FALSE.
545	IF (is_root_prc) THEN
546	IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
547	& ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
548	rectilinear = .TRUE.
549	ENDIF
550	ENDIF
551	CALL bcast(rectilinear)
552	IF (rectilinear) THEN
553	ALLOCATE(lon_rect(iim),stat=ier)
554	IF (ier .NE. 0) THEN
555	WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
556	STOP 'intersurf_history'
557	ENDIF
558	ALLOCATE(lat_rect(jjm),stat=ier)
559	IF (ier .NE. 0) THEN
560	WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
561	STOP 'intersurf_history'
562	ENDIF
563	lon_rect(:) = lon(:,1)
564	lat_rect(:) = lat(1,:)
565	ENDIF
566	!-
567	! allocate arrays
568	!-
569	!
570	a_er = .FALSE.
571	ALLOCATE (hist_PFTaxis(nvm),stat=ier)
572	a_er = a_er .OR. (ier.NE.0)
573	ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
574	a_er = a_er .OR. (ier.NE.0)
575	ALLOCATE (soiltile(kjpindex,nstm),stat=ier)
576	a_er = a_er .OR. (ier.NE.0)
577	ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
578	a_er = a_er .OR. (ier.NE.0)
579	ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
580	a_er = a_er .OR. (ier.NE.0)
581	ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
582	a_er = a_er .OR. (ier.NE.0)
583	ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
584	a_er = a_er .OR. (ier.NE.0)
585	ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
586	a_er = a_er .OR. (ier.NE.0)
587	ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
588	a_er = a_er .OR. (ier.NE.0)
589	ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
590	a_er = a_er .OR. (ier.NE.0)
591	ALLOCATE (njsc(kjpindex),stat=ier)
592	a_er = a_er .OR. (ier.NE.0)
593	ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
594	a_er = a_er .OR. (ier.NE.0)
595	ALLOCATE (reinf_slope(kjpindex),stat=ier)
596	a_er = a_er .OR. (ier.NE.0)
597	ALLOCATE (t2m(kjpindex),stat=ier)
598	a_er = a_er .OR. (ier.NE.0)
599	ALLOCATE (t2m_min(kjpindex),stat=ier)
600	a_er = a_er .OR. (ier.NE.0)
601	ALLOCATE (temp_sol(kjpindex),stat=ier)
602	a_er = a_er .OR. (ier.NE.0)
603	ALLOCATE (soiltemp(kjpindex,nslm),stat=ier)
604	a_er = a_er .OR. (ier.NE.0)
605	ALLOCATE (soilhum(kjpindex,nslm),stat=ier)
606	a_er = a_er .OR. (ier.NE.0)
607	ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
608	a_er = a_er .OR. (ier.NE.0)
609	ALLOCATE (litterhum(kjpindex),stat=ier)
610	a_er = a_er .OR. (ier.NE.0)
611	ALLOCATE (precip_rain(kjpindex),stat=ier)
612	a_er = a_er .OR. (ier.NE.0)
613	ALLOCATE (precip_snow(kjpindex),stat=ier)
614	a_er = a_er .OR. (ier.NE.0)
615	ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
616	a_er = a_er .OR. (ier.NE.0)
617	ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
618	a_er = a_er .OR. (ier.NE.0)
619	ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
620	a_er = a_er .OR. (ier.NE.0)
621	ALLOCATE (height_x(kjpindex,nvm),stat=ier)
622	a_er = a_er .OR. (ier.NE.0)
623	ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
624	a_er = a_er .OR. (ier.NE.0)
625	ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
626	a_er = a_er .OR. (ier.NE.0)
627	ALLOCATE (fco2_lu(kjpindex),stat=ier)
628	a_er = a_er .OR. (ier.NE.0)
629	ALLOCATE (temp_growth(kjpindex),stat=ier)
630	a_er = a_er .OR. (ier.NE.0)
631	ALLOCATE (tot_bare_soil(kjpindex),stat=ier)
632	a_er = a_er .OR. (ier.NE.0)
633	ALLOCATE (frac_age_x(kjpindex,nvm,nleafages),stat=ier)
634	a_er = a_er .OR. (ier.NE.0)
635	IF (a_er) THEN
636	CALL ipslerr (3,'teststomate', &
637	& 'PROBLEM WITH ALLOCATION', &
638	& 'for local variables 1','')
639	ENDIF
640	!
641	! prepare forcing
642	! forcing file of STOMATE is read by block
643	! in order to minimize the reading frequency
644	! Here is done the calculation of the number
645	! of time steps we load in memory
646	!
647	max_totsize = 50
648	CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
649	max_totsize = max_totsize * 1000000
650
651	totsize_1step = SIZE(soiltile(:,3))*KIND(soiltile(:,3)) + &
652	SIZE(humrel_x)*KIND(humrel_x) + &
653	SIZE(litterhum)*KIND(litterhum) + &
654	SIZE(t2m)*KIND(t2m) + &
655	SIZE(t2m_min)*KIND(t2m_min) + &
656	SIZE(temp_sol)*KIND(temp_sol) + &
657	SIZE(soiltemp)*KIND(soiltemp) + &
658	SIZE(soilhum)*KIND(soilhum) + &
659	SIZE(precip_rain)*KIND(precip_rain) + &
660	SIZE(precip_snow)*KIND(precip_snow) + &
661	SIZE(gpp_x)*KIND(gpp_x) + &
662	SIZE(veget_force_x)*KIND(veget_force_x) + &
663	SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
664	SIZE(lai_force_x)*KIND(lai_force_x)
665	CALL reduce_sum(totsize_1step,totsize_tmp)
666	CALL bcast(totsize_tmp)
667	totsize_1step=totsize_tmp
668
669	! check for consistency of the total number of forcing steps
670	IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
671	CALL ipslerr_p (3,'teststomate', &
672	& 'stomate: error with total number of forcing steps', &
673	& 'nsft','teststomate computation different with forcing file value.')
674	ENDIF
675	! number of forcing steps in memory
676	nsfm = MIN(nsft, &
677	& MAX(1,NINT( REAL(max_totsize,r_std) &
678	& /REAL(totsize_1step,r_std))))
679	!-
680	WRITE(numout,*) 'Offline forcing of Stomate:'
681	WRITE(numout,*) ' Total number of forcing steps:',nsft
682	WRITE(numout,*) ' Number of forcing steps in memory:',nsfm
683	!-
684	! init_forcing defined into stomate.f90
685	! allocate and set to zero driving variables of STOMATE
686	! ie variables that are read in the forcing file
687	!-
688	CALL init_forcing(kjpindex,nsfm,nsft)
689	!-
690	! ensure that we read all new forcing states
691	iisf = nsfm
692	! initialize the table that contains the indices
693	! of the forcing states that will be in memory
694	isf(:) = (/ (i,i=1,nsfm) /)
695
696	nf_written(:) = .FALSE.
697	nf_written(isf(:)) = .TRUE.
698
699	!-
700	! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
701	!-
702	itau_step = NINT(dt_force/dtradia)
703	IF (printlev_loc>=3) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
704	!
705	CALL ioconf_startdate(date0)
706	!
707	CALL time_initialize(itau_dep, date0, dt_sechiba, "START")
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	CALL getin_p ("DEPTH_MAX_H", zmaxh)
920	!
921	!MM Problem here with dpu which depends on soil type
922	DO l = 1, nslm-1
923	! first 2.0 is dpu
924	! second 2.0 is average
925	diaglev(l) = zmaxh/(2*(nslm-1) -1) ( ( 2(l-1) -1) + ( 2(l) -1) ) / 2.0
926	ENDDO
927	diaglev(nslm) = zmaxh
928	!
929	!-
930	! Read the parameters in the "*.def" files
931	!-
932	!
933	!Config Key = CLAYFRACTION_DEFAULT
934	!Config Desc =
935	!Config If = OK_SECHIBA
936	!Config Def = 0.2
937	!Config Help =
938	!Config Units = [-]
939	CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
940	!
941	!Config Key = MIN_VEGFRAC
942	!Config Desc = Minimal fraction of mesh a vegetation type can occupy
943	!Config If = OK_SECHIBA
944	!Config Def = 0.001
945	!Config Help =
946	!Config Units = [-]
947	CALL getin_p('MIN_VEGFRAC',min_vegfrac)
948	!
949	!Config Key = STEMPDIAG_BID
950	!Config Desc = only needed for an initial LAI if there is no restart file
951	!Config If = OK_SECHIBA
952	!Config Def = 280.
953	!Config Help =
954	!Config Units = [K]
955	CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
956	!
957	!-
958	CALL ioget_expval(val_exp)
959	ldrestart_read = .TRUE.
960	ldrestart_write = .FALSE.
961	!-
962	! read some variables we need from SECHIBA's restart file
963	!-
964	CALL slowproc_initialize (itau, kjpij, kjpindex, &
965	rest_id_sec, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, &
966	indices, indexveg, lalo, neighbours, &
967	resolution, contfrac, t2m, &
968	soiltile, reinf_slope, deadleaf_cover, assim_param_x, &
969	lai_x, frac_age_x, height_x, veget_x, &
970	frac_nobio, njsc, veget_max_x, tot_bare_soil, &
971	totfrac_nobio, qsintmax_x, co2_flux, fco2_lu, temp_growth)
972
973	! Old interface to slowproc_main, before revision 2581
974	! CALL slowproc_main &
975	! & (itau, kjpij, kjpindex, dt_force, date0, &
976	! & ldrestart_read, ldrestart_write, &
977	! & indices, indexveg, lalo, neighbours, resolution, contfrac, soiltile, reinf_slope, &
978	! & t2m, t2m_min, temp_sol, soiltemp, &
979	! & humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
980	! & deadleaf_cover, assim_param_x, lai_x, frac_age_x, height_x, veget_x, &
981	! & frac_nobio, njsc, veget_max_x, totfrac_nobio, qsintmax_x, &
982	! & 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)
983	! correct date
984
985
986	days_since_beg=1
987
988	CALL time_nextstep( itau_dep+itau_step )
989
990	lstep_init_intersurf=.FALSE.
991	!-
992	! 4. Time loop
993	!â»
994	DO itau = itau_dep+itau_step,itau_fin,itau_step
995	!
996	!
997	CALL time_nextstep( itau )
998	!-
999	! next forcing state
1000	iisf = iisf+1
1001	IF (printlev_loc>=3) WRITE(numout,*) "itau,iisf : ",itau,iisf
1002	!---
1003	IF (iisf .GT. nsfm) THEN
1004	!-----
1005	!---- we have to read new forcing states from the file
1006	!-----
1007	!---- determine blocks of forcing states that are contiguous in memory
1008	!-----
1009	CALL stomate_forcing_read(forcing_id,nsfm)
1010
1011	!--------------------------
1012
1013	!-----
1014	!---- determine which forcing states must be read next time
1015	!-----
1016	isf(1) = isf(nsfm)+1
1017	IF ( isf(1) .GT. nsft ) isf(1) = 1
1018	DO iiisf = 2, nsfm
1019	isf(iiisf) = isf(iiisf-1)+1
1020	IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
1021	ENDDO
1022	nf_written(isf(:)) = .TRUE.
1023	!---- start again at first forcing state
1024	iisf = 1
1025	ENDIF
1026	humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
1027	litterhum(:) = litterhum_daily_fm(:,iisf)
1028	t2m(:) = t2m_daily_fm(:,iisf)
1029	t2m_min(:) = t2m_min_daily_fm(:,iisf)
1030	temp_sol(:) = tsurf_daily_fm(:,iisf)
1031	soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
1032	soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
1033	precip_rain(:) = precip_fm(:,iisf)
1034	gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
1035	veget_force_x(:,:) = veget_fm(:,:,iisf)
1036	veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
1037	lai_force_x(:,:) = lai_fm(:,:,iisf)
1038	WHERE ( t2m(:) .LT. ZeroCelsius )
1039	precip_snow(:) = precip_rain(:)
1040	precip_rain(:) = zero
1041	ELSEWHERE
1042	precip_snow(:) = zero
1043	ENDWHERE
1044	!---
1045	!-- scale GPP to new lai and veget_max
1046	!---
1047	WHERE ( lai_x(:,:) .EQ. zero ) gpp_x(:,:) = zero
1048	!-- scale GPP to new LAI
1049	WHERE (lai_force_x(:,:) .GT. zero )
1050	gpp_x(:,:) = gpp_x(:,:)ATAN(2.lai_x(:,:)) &
1051	& /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
1052	ENDWHERE
1053	!- scale GPP to new veget_max
1054	WHERE (veget_max_force_x(:,:) .GT. zero )
1055	gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
1056	ENDWHERE
1057	!-
1058	!- number crunching
1059	!-
1060	ldrestart_read = .FALSE.
1061	ldrestart_write = .FALSE.
1062	CALL slowproc_main &
1063	& (itau, kjpij, kjpindex, &
1064	& indices, indexveg, lalo, neighbours, resolution, contfrac, soiltile, &
1065	& t2m, temp_sol, soiltemp, &
1066	& humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1067	& deadleaf_cover, assim_param_x, lai_x, frac_age_x, height_x, veget_x, &
1068	& frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
1069	& rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, &
1070	& co2_flux, fco2_lu, temp_growth, tot_bare_soil)
1071
1072	ENDDO ! end of the time loop
1073
1074
1075	!-
1076	! 5. write restart files
1077	!-
1078	IF (is_root_prc) THEN
1079	! first, read and write variables that are not managed otherwise
1080	taboo_vars = &
1081	& '$lat$ $lon$ $lev$ $veget_year$ '// &
1082	& '$height$ $veget$ $veget_max$ $frac_nobio$ '// &
1083	& '$lai$ $soiltile_frac$ $clay_frac$ '// &
1084	& '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
1085	!-
1086	CALL ioget_vname(rest_id_sec, nbvar, varnames)
1087	!-
1088	DO iv = 1, nbvar
1089	!-- check if the variable is to be written here
1090	IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
1091	IF (printlev_loc>=3) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
1092
1093	!---- get variable dimensions, especially 3rd dimension
1094	CALL ioget_vdim &
1095	& (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
1096	l1d = ALL(vardims(1:varnbdim) .EQ. 1)
1097	!---- read it
1098	IF (l1d) THEN
1099	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
1100	1,1,1,itau_dep,.TRUE.,var_1d)
1101	ELSE
1102	ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
1103	IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
1104	CALL restget (rest_id_sec,TRIM(varnames(iv)), &
1105	nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
1106	"gather",nbp_glo,index_g)
1107	ENDIF
1108	!---- write it
1109	IF (l1d) THEN
1110	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
1111	1,1,1,itau_fin,var_1d)
1112	ELSE
1113	CALL restput (rest_id_sec,TRIM(varnames(iv)), &
1114	nbp_glo,vardims(3),1,itau_fin,var_3d, &
1115	'scatter',nbp_glo,index_g)
1116	DEALLOCATE (var_3d)
1117	ENDIF
1118	ENDIF
1119	ENDDO
1120	ENDIF
1121
1122	#ifdef CPP_PARA
1123	CALL MPI_BARRIER(MPI_COMM_ORCH,ierr)
1124	#endif
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.

Download in other formats: