Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

teststomate.f90 @ 7746

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

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

Context Navigation

source: branches/publications/ORCHIDEE_Biochar/src_driver/teststomate.f90 @ 7746

Download in other formats: