Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

teststomate.f90 @ 7325

Last change on this file since 7325 was 6339, checked in by elodie.salmon, 5 years ago
New: include methane emissions for peatlands and other pft
File size: 58.2 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_mict_peat_ch4/src_driver/teststomate.f90 @ 7325

Download in other formats: