Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

teststomate.f90 @ 7325

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_driver/teststomate.f90 @ 7325

Download in other formats: