Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

teststomate.f90 @ 7442

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

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

Context Navigation

source: branches/publications/ORCHIDEE-MICT-BIOENERGY_r7298/src_driver/teststomate.f90 @ 7442

Download in other formats: