source: branches/ORCHIDEE_EXT/ORCHIDEE_OL/teststomate.f90 @ 472

Last change on this file since 472 was 472, checked in by didier.solyga, 10 years ago

Replace ok_sechiba, ok_stomate by the following structure control as a argument for the routine pft_parameters_main

File size: 35.2 KB
Line 
1PROGRAM teststomate
2!---------------------------------------------------------------------
3! This program reads a forcing file for STOMATE
4! which was created with STOMATE coupled to SECHIBA.
5! It then integrates STOMATE for a
6! given number of years using this forcing.
7! The GPP is scaled to the new calculated LAI...
8! Nothing is done for soil moisture which should
9! actually evolve with the vegetation.
10!---------------------------------------------------------------------
11!- IPSL (2006)
12!-  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13  USE netcdf
14!-
15  USE defprec
16  USE constantes
17  USE pft_parameters
18  USE stomate_data
19  USE ioipsl
20  USE grid
21  USE slowproc
22  USE stomate
23  USE intersurf, ONLY: stom_define_history, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time,impose_param
24  USE parallel
25!-
26  IMPLICIT NONE
27!-
28! Declarations
29!-
30  INTEGER(i_std)                            :: kjpij,kjpindex
31  REAL(r_std)                               :: dtradia,dt_force
32
33  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indices
34  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indexveg
35  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: soiltype, veget_x
36  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: totfrac_nobio
37  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: frac_nobio
38  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_max_x
39  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: lai_x
40  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_force_x
41  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: veget_max_force_x
42  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: lai_force_x
43  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: t2m,t2m_min,temp_sol
44  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: soiltemp,soilhum
45  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: humrel_x
46  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: litterhum
47  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: precip_rain,precip_snow
48  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: gpp_x
49  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: deadleaf_cover
50  REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE  :: assim_param_x
51  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: height_x
52  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: qsintmax_x
53  REAL(r_std),DIMENSION(:,:),ALLOCATABLE    :: co2_flux
54  REAL(r_std),DIMENSION(:),ALLOCATABLE      :: fco2_lu
55
56  INTEGER(i_std),DIMENSION(:),ALLOCATABLE   :: indices_g
57  REAL(r_std),DIMENSION(:),ALLOCATABLE   :: x_indices_g
58  REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g
59
60  INTEGER    :: ier,iret
61  INTEGER    :: ncfid
62  CHARACTER(LEN=20),SAVE                      :: thecalendar='noleap'
63
64  LOGICAL    :: a_er
65  CHARACTER(LEN=80) :: &
66 &  dri_restname_in,dri_restname_out, &
67 &  sec_restname_in,sec_restname_out, &
68 &  sto_restname_in,sto_restname_out, &
69 &  stom_histname, stom_ipcc_histname
70  INTEGER(i_std)                    :: iim,jjm,llm
71  REAL, ALLOCATABLE, DIMENSION(:,:)  :: lon, lat
72  REAL, ALLOCATABLE, DIMENSION(:)    :: lev
73  LOGICAL                            :: rectilinear
74  REAL, ALLOCATABLE, DIMENSION(:)    :: lon_rect, lat_rect
75  REAL(r_std)                       :: dt
76  INTEGER(i_std)                    :: itau_dep,itau_fin,itau,itau_len,itau_step
77  REAL(r_std)                       :: date0
78  INTEGER(i_std)                    :: rest_id_sec,rest_id_sto
79  INTEGER(i_std)                    :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC
80  CHARACTER(LEN=30)                :: time_str
81  REAL(r_std)                       :: dt_slow_
82  REAL                             :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc
83!!$ >>  DS for the externalized version
84  REAL(r_std),ALLOCATABLE,DIMENSION(:) :: hist_PFTaxis
85!!$ >>  DS
86  REAL(r_std),DIMENSION(10)         :: hist_pool_10axis     
87  REAL(r_std),DIMENSION(100)        :: hist_pool_100axis     
88  REAL(r_std),DIMENSION(11)         :: hist_pool_11axis     
89  REAL(r_std),DIMENSION(101)        :: hist_pool_101axis     
90  INTEGER                          :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id,hori_id
91  INTEGER                          :: hist_pool_10axis_id
92  INTEGER                          :: hist_pool_100axis_id
93  INTEGER                          :: hist_pool_11axis_id
94  INTEGER                          :: hist_pool_101axis_id
95  INTEGER(i_std)                    :: i,j,iv
96  LOGICAL                          :: ldrestart_read,ldrestart_write
97  LOGICAL                          :: l1d
98  INTEGER(i_std),PARAMETER          :: nbvarmax=200
99  INTEGER(i_std)                    :: nbvar
100  CHARACTER*50,DIMENSION(nbvarmax) :: varnames
101  INTEGER(i_std)                    :: varnbdim
102  INTEGER(i_std),PARAMETER          :: varnbdim_max=20
103  INTEGER,DIMENSION(varnbdim_max)  :: vardims
104  CHARACTER*200                    :: taboo_vars
105  REAL(r_std),DIMENSION(1)         :: xtmp
106  INTEGER                          :: nsfm,nsft
107  INTEGER                          :: iisf,iiisf
108  INTEGER(i_std)                    :: max_totsize,totsize_1step,totsize_tmp
109
110  INTEGER                          :: vid
111  CHARACTER(LEN=100)               :: forcing_name
112  REAL                             :: x
113
114  REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d
115  REAL(r_std) :: var_1d(1)
116
117!-
118  REAL(r_std)                             :: time_sec,time_step_sec
119  REAL(r_std)                             :: time_dri,time_step_dri
120  REAL(r_std),DIMENSION(1)                :: r1d
121  REAL(r_std)                             :: julian,djulian
122
123  INTEGER(i_std)                                :: ji,jv,l
124
125  LOGICAL :: debug
126
127!---------------------------------------------------------------------
128
129  CALL init_para(.FALSE.)
130  CALL init_timer
131
132  IF (is_root_prc) THEN
133     !-
134     ! open STOMATE's forcing file to read some basic info
135     !-
136     forcing_name = 'stomate_forcing.nc'
137     CALL getin ('STOMATE_FORCING_NAME',forcing_name)
138     iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id)
139     IF (iret /= NF90_NOERR) THEN
140        CALL ipslerr (3,'teststomate', &
141             &        'Could not open file : ', &
142             &          forcing_name,'(Do you have forget it ?)')
143     ENDIF
144     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia)
145     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force)
146     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x)
147     nsft = NINT(x)
148     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x)
149     kjpij = NINT(x)
150     ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x)
151     nbp_glo = NINT(x)
152  ENDIF
153  CALL bcast(dtradia)
154  CALL bcast(dt_force)
155  CALL bcast(nsft)
156  CALL bcast(nbp_glo)
157  !-
158  write(numout,*) 'ATTENTION',dtradia,dt_force
159  !-
160  ! read info about land points
161  !-
162  IF (is_root_prc) THEN
163     a_er=.FALSE.
164     ALLOCATE (indices_g(nbp_glo),stat=ier)
165     a_er = a_er .OR. (ier.NE.0)
166     IF (a_er) THEN
167        CALL ipslerr (3,'teststomate', &
168             &        'PROBLEM WITH ALLOCATION', &
169             &        'for local variables 1','')
170     ENDIF
171     !
172     ALLOCATE (x_indices_g(nbp_glo),stat=ier)
173     a_er = a_er .OR. (ier.NE.0)
174     IF (a_er) THEN
175        CALL ipslerr (3,'teststomate', &
176             &        'PROBLEM WITH ALLOCATION', &
177             &        'for global variables 1','')
178     ENDIF
179     ier = NF90_INQ_VARID (forcing_id,'index',vid)
180     IF (ier .NE. 0) THEN
181        CALL ipslerr (3,'teststomate', &
182             &        'PROBLEM WITH ALLOCATION', &
183             &        'for global variables 1','')
184     ENDIF
185     ier = NF90_GET_VAR   (forcing_id,vid,x_indices_g)
186     IF (iret /= NF90_NOERR) THEN
187        CALL ipslerr (3,'teststomate', &
188             &        'PROBLEM WITH variable "index" in file ', &
189             &        forcing_name,'(check this file)')
190     ENDIF
191     indices_g(:) = NINT(x_indices_g(:))
192     DEALLOCATE (x_indices_g)
193  ELSE
194     ALLOCATE (indices_g(0))
195  ENDIF
196!---------------------------------------------------------------------
197!-
198! set debug to have more information
199!-
200  !Config  Key  = DEBUG_INFO
201  !Config  Desc = Flag for debug information
202  !Config  Def  = n
203  !Config  Help = This option allows to switch on the output of debug
204  !Config         information without recompiling the code.
205!-
206  debug = .FALSE.
207  CALL getin_p('DEBUG_INFO',debug)
208  !
209  !Config Key  = LONGPRINT
210  !Config Desc = ORCHIDEE will print more messages
211  !Config Def  = n
212  !Config Help = This flag permits to print more debug messages in the run.
213  !
214  long_print = .FALSE.
215  CALL getin_p('LONGPRINT',long_print)
216
217  !-
218  ! activate CO2, STOMATE, but not sechiba
219  !-
220  control%river_routing = .FALSE.
221  control%hydrol_cwrr = .FALSE.
222  control%ok_sechiba = .FALSE.
223  !
224  control%stomate_watchout = .TRUE.
225  control%ok_co2 = .TRUE.
226  control%ok_stomate = .TRUE.
227  !-
228  ! is DGVM activated?
229  !-
230  control%ok_dgvm = .FALSE.
231  CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm)
232  WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm
233
234!!$ >> DS added for the externalised version
235
236  ! 1. Number of PFTs defined by the user
237
238  !Config Key  = NVM
239  !Config Desc = number of PFTs 
240  !Config if  = ANYTIME
241  !Config Def  = 13
242  !Config Help = The number of vegetation types define by the user
243  !Config Units = NONE
244  !
245  CALL getin_p('NVM',nvm)
246  WRITE(numout,*)'the number of pfts is : ', nvm
247
248  ! 2. Should we read the parameters in the run.def file ?
249 
250  !Config Key  = IMPOSE_PARAM
251  !Config Desc = Do you impose the values of the parameters?
252  !Config if  = ANYTIME
253  !Config Def  = y
254  !Config Help = This flag can deactivate the reading of some parameters.
255  !               Useful if you want to use the standard values without commenting the run.def
256  !Config Units = NONE
257  !
258  CALL getin_p('IMPOSE_PARAM',impose_param)
259
260  ! 3. Allocate and intialize the pft parameters
261
262  CALL pft_parameters_main(control)
263
264  ! 4. Activation sub-models of ORCHIDEE
265
266  CALL activate_sub_models(control)
267
268  ! 5. Vegetation configuration (impose_veg, land_use, lcchange...previously in slowproc)
269
270  CALL veget_config
271
272  ! 6. Read the parameters in the run.def file
273
274  IF ( impose_param) THEN
275     CALL getin_pft_parameters
276     CALL getin_stomate_pft_parameters
277     CALL getin_co2_parameters
278     CALL getin_stomate_parameters
279  ENDIF
280  !-
281  IF (control%ok_dgvm) THEN
282     IF ( impose_param ) THEN
283        CALL getin_dgvm_parameters
284     ENDIF
285  ENDIF
286
287!!$ >> DS
288
289  !-
290  ! restart files
291  !-
292  IF (is_root_prc) THEN
293     ! Sechiba's restart files
294     sec_restname_in = 'sechiba_start.nc'
295     CALL getin('SECHIBA_restart_in',sec_restname_in)
296     WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in)
297     IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN
298        STOP 'Need a restart file for Sechiba'
299     ENDIF
300     sec_restname_out = 'sechiba_rest_out.nc'
301     CALL getin('SECHIBA_rest_out',sec_restname_out)
302     WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out)
303     ! Stomate's restart files
304     sto_restname_in = 'stomate_start.nc'
305     CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in)
306     WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)
307     sto_restname_out = 'stomate_rest_out.nc'
308     CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out)
309     WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)
310
311     !-
312     ! We need to know iim, jjm.
313     ! Get them from the restart files themselves.
314     !-
315     iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid)
316     IF (iret /= NF90_NOERR) THEN
317        CALL ipslerr (3,'teststomate', &
318             &        'Could not open file : ', &
319             &          sec_restname_in,'(Do you have forget it ?)')
320     ENDIF
321     iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g)
322     iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g)
323     iret = NF90_INQ_VARID (ncfid, "time", iv)
324     iret = NF90_GET_ATT (ncfid, iv, 'calendar',thecalendar)
325     iret = NF90_CLOSE (ncfid)
326     i=INDEX(thecalendar,ACHAR(0))
327     IF ( i > 0 ) THEN
328        thecalendar(i:20)=' '
329     ENDIF
330  ENDIF
331  CALL bcast(iim_g)
332  CALL bcast(jjm_g)
333  CALL bcast(thecalendar)
334  !-
335  ! calendar
336  !-
337  CALL ioconf_calendar (thecalendar)
338  CALL ioget_calendar  (one_year,one_day)
339  !
340  ! Parallelization :
341  !
342  CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g)
343  kjpindex=nbp_loc
344  jjm=jj_nb
345  iim=iim_g
346  kjpij=iim*jjm
347  IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm
348  !-
349  !-
350  ! read info about grids
351  !-
352  !-
353  llm=1
354  ALLOCATE(lev(llm))
355  IF (is_root_prc) THEN
356     !-
357     ier = NF90_INQ_VARID (forcing_id,'lalo',vid)
358     ier = NF90_GET_VAR   (forcing_id,vid,lalo_g)
359     !-
360     ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier)
361     ier = NF90_INQ_VARID (forcing_id,'neighbours',vid)
362     ier = NF90_GET_VAR   (forcing_id,vid,x_neighbours_g)
363     neighbours_g(:,:) = NINT(x_neighbours_g(:,:))
364     DEALLOCATE (x_neighbours_g)
365     !-
366     ier = NF90_INQ_VARID (forcing_id,'resolution',vid)
367     ier = NF90_GET_VAR   (forcing_id,vid,resolution_g)
368     !-
369     ier = NF90_INQ_VARID (forcing_id,'contfrac',vid)
370     ier = NF90_GET_VAR   (forcing_id,vid,contfrac_g)
371
372     lon_g(:,:) = 0.0
373     lat_g(:,:) = 0.0
374     lev(1)   = 0.0
375     !-
376     CALL restini &
377          & (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
378          &  sec_restname_out, itau_dep, date0, dt, rest_id_sec)
379     !-
380     IF ( dt .NE. dtradia ) THEN
381        WRITE(numout,*) 'dt',dt
382        WRITE(numout,*) 'dtradia',dtradia
383        CALL ipslerr (3,'teststomate', &
384             &        'PROBLEM with time steps.', &
385             &          sec_restname_in,'(dt .NE. dtradia)')
386     ENDIF
387     !-
388     CALL restini &
389          & (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, &
390          &  sto_restname_out, itau_dep, date0, dt, rest_id_sto)
391     !-
392     IF ( dt .NE. dtradia ) THEN
393        WRITE(numout,*) 'dt',dt
394        WRITE(numout,*) 'dtradia',dtradia
395        CALL ipslerr (3,'teststomate', &
396             &        'PROBLEM with time steps.', &
397             &          sto_restname_in,'(dt .NE. dtradia)')
398     ENDIF
399  ENDIF
400  CALL bcast(rest_id_sec)
401  CALL bcast(rest_id_sto)
402  CALL bcast(itau_dep)
403  CALL bcast(date0)
404  CALL bcast(dt)
405  CALL bcast(lev)
406  !---
407  !--- Create the index table
408  !---
409  !--- This job return a LOCAL kindex
410  !---
411  ALLOCATE (indices(kjpindex),stat=ier)
412  IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo)
413  CALL scatter(indices_g,indices)
414  indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g
415  IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex)
416
417  !---
418  !--- initialize global grid
419  !---
420  CALL init_grid( kjpindex ) 
421  CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g)
422
423  !---
424  !--- initialize local grid
425  !---
426  jlandindex = (((indices(1:kjpindex)-1)/iim) + 1)
427  if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex)
428  ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim)
429  IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex)
430  ALLOCATE(lon(iim,jjm))
431  ALLOCATE(lat(iim,jjm))
432  lon=zero
433  lat=zero
434  CALL scatter2D(lon_g,lon)
435  CALL scatter2D(lat_g,lat)
436
437  DO ji=1,kjpindex
438
439     j = jlandindex(ji)
440     i = ilandindex(ji)
441
442     !- Create the internal coordinate table
443!-
444     lalo(ji,1) = lat(i,j)
445     lalo(ji,2) = lon(i,j)
446  ENDDO
447  CALL scatter(neighbours_g,neighbours)
448  CALL scatter(resolution_g,resolution)
449  CALL scatter(contfrac_g,contfrac)
450  CALL scatter(area_g,area)
451  !-
452  !- Check if we have by any change a rectilinear grid. This would allow us to
453  !- simplify the output files.
454  !
455  rectilinear = .FALSE.
456  IF (is_root_prc) THEN
457     IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. &
458       & ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN
459        rectilinear = .TRUE.
460     ENDIF
461  ENDIF
462  CALL bcast(rectilinear)
463  IF (rectilinear) THEN
464     ALLOCATE(lon_rect(iim),stat=ier)
465     IF (ier .NE. 0) THEN
466        WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim
467        STOP 'intersurf_history'
468     ENDIF
469     ALLOCATE(lat_rect(jjm),stat=ier)
470     IF (ier .NE. 0) THEN
471        WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm
472        STOP 'intersurf_history'
473     ENDIF
474     lon_rect(:) = lon(:,1)
475     lat_rect(:) = lat(1,:)
476  ENDIF
477  !-
478  ! allocate arrays
479  !-
480  !
481  a_er = .FALSE.
482!!$ >> DS added for the externalised version
483  ALLOCATE (hist_PFTaxis(nvm),stat=ier)
484  a_er = a_er .OR. (ier.NE.0)
485! 01/07/2011 allocate fco2_lu
486  ALLOCATE (fco2_lu(kjpindex),stat=ier)
487  a_er = a_er .OR. (ier.NE.0) 
488!!$ end add DS
489  ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
490  a_er = a_er .OR. (ier.NE.0)
491  ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
492  a_er = a_er .OR. (ier.NE.0)
493  ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
494  a_er = a_er .OR. (ier.NE.0)
495  ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
496  a_er = a_er .OR. (ier.NE.0)
497  ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
498  a_er = a_er .OR. (ier.NE.0)
499  ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
500  a_er = a_er .OR. (ier.NE.0)
501  ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
502  a_er = a_er .OR. (ier.NE.0)
503  ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
504  a_er = a_er .OR. (ier.NE.0)
505  ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
506  a_er = a_er .OR. (ier.NE.0)
507  ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
508  a_er = a_er .OR. (ier.NE.0)
509  ALLOCATE (t2m(kjpindex),stat=ier)
510  a_er = a_er .OR. (ier.NE.0)
511  ALLOCATE (t2m_min(kjpindex),stat=ier)
512  a_er = a_er .OR. (ier.NE.0)
513  ALLOCATE (temp_sol(kjpindex),stat=ier)
514  a_er = a_er .OR. (ier.NE.0)
515  ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
516  a_er = a_er .OR. (ier.NE.0)
517  ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
518  a_er = a_er .OR. (ier.NE.0)
519  ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
520  a_er = a_er .OR. (ier.NE.0)
521  ALLOCATE (litterhum(kjpindex),stat=ier)
522  a_er = a_er .OR. (ier.NE.0)
523  ALLOCATE (precip_rain(kjpindex),stat=ier)
524  a_er = a_er .OR. (ier.NE.0)
525  ALLOCATE (precip_snow(kjpindex),stat=ier)
526  a_er = a_er .OR. (ier.NE.0)
527  ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
528  a_er = a_er .OR. (ier.NE.0)
529  ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
530  a_er = a_er .OR. (ier.NE.0)
531  ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
532  a_er = a_er .OR. (ier.NE.0)
533  ALLOCATE (height_x(kjpindex,nvm),stat=ier)
534  a_er = a_er .OR. (ier.NE.0)
535  ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
536  a_er = a_er .OR. (ier.NE.0)
537  ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
538  a_er = a_er .OR. (ier.NE.0)
539  ALLOCATE (fco2_lu(kjpindex),stat=ier)
540  a_er = a_er .OR. (ier.NE.0)
541  IF (a_er) THEN
542     CALL ipslerr (3,'teststomate', &
543          &        'PROBLEM WITH ALLOCATION', &
544          &        'for local variables 1','')
545  ENDIF
546  !
547  ! prepare forcing
548  !
549  max_totsize = 50
550  CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
551  max_totsize = max_totsize * 1000000
552
553  totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
554       SIZE(humrel_x)*KIND(humrel_x) + &
555       SIZE(litterhum)*KIND(litterhum) + &
556       SIZE(t2m)*KIND(t2m) + &
557       SIZE(t2m_min)*KIND(t2m_min) + &
558       SIZE(temp_sol)*KIND(temp_sol) + &
559       SIZE(soiltemp)*KIND(soiltemp) + &
560       SIZE(soilhum)*KIND(soilhum) + &
561       SIZE(precip_rain)*KIND(precip_rain) + &
562       SIZE(precip_snow)*KIND(precip_snow) + &
563       SIZE(gpp_x)*KIND(gpp_x) + &
564       SIZE(veget_force_x)*KIND(veget_force_x) + &
565       SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
566       SIZE(lai_force_x)*KIND(lai_force_x)
567  CALL reduce_sum(totsize_1step,totsize_tmp)
568  CALL bcast(totsize_tmp)
569  totsize_1step=totsize_tmp 
570
571! total number of forcing steps
572  IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
573     CALL ipslerr (3,'teststomate', &
574          &        'stomate: error with total number of forcing steps', &
575          &        'nsft','teststomate computation different with forcing file value.')
576  ENDIF
577! number of forcing steps in memory
578  nsfm = MIN(nsft, &
579       &       MAX(1,NINT( REAL(max_totsize,r_std) &
580       &                  /REAL(totsize_1step,r_std))))
581!-
582  WRITE(numout,*) 'Offline forcing of Stomate:'
583  WRITE(numout,*) '  Total number of forcing steps:',nsft
584  WRITE(numout,*) '  Number of forcing steps in memory:',nsfm
585!-
586  CALL init_forcing(kjpindex,nsfm,nsft)
587  !-
588! ensure that we read all new forcing states
589  iisf = nsfm
590! initialize that table that contains the indices
591! of the forcing states that will be in memory
592  isf(:) = (/ (i,i=1,nsfm) /)
593
594  nf_written(:) = .FALSE.
595  nf_written(isf(:)) = .TRUE.
596
597!-
598! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
599!-
600  itau_step = NINT(dt_force/dtradia)
601  IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
602  !
603  CALL ioconf_startdate(date0)
604  CALL intsurf_time( itau_dep, date0, dtradia )
605!-
606! Length of integration
607!-
608  WRITE(time_str,'(a)') '1Y'
609  CALL getin_p ('TIME_LENGTH', time_str)
610! transform into itau
611  CALL tlen2itau(time_str, dt, date0, itau_len)
612! itau_len*dtradia must be a multiple of dt_force
613  itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
614       &                *dt_force/dtradia)
615  !-
616  itau_fin = itau_dep+itau_len
617  !-
618  ! set up STOMATE history file
619  !-
620  !Config  Key  = STOMATE_OUTPUT_FILE
621  !Config  Desc = Name of file in which STOMATE's output is going
622  !Config         to be written
623  !Config  Def  = stomate_history.nc
624  !Config  Help = This file is going to be created by the model
625  !Config         and will contain the output from the model.
626  !Config         This file is a truly COADS compliant netCDF file.
627  !Config         It will be generated by the hist software from
628  !Config         the IOIPSL package.
629  !-
630  stom_histname='stomate_history.nc'
631  CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
632  WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
633  !-
634  !Config  Key  = STOMATE_HIST_DT
635  !Config  Desc = STOMATE history time step (d)
636  !Config  Def  = 10.
637  !Config  Help = Time step of the STOMATE history file
638  !-
639  hist_days_stom = 10.
640  CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
641  IF ( hist_days_stom == -1. ) THEN
642     hist_dt_stom = -1.
643     WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
644  ELSE
645     hist_dt_stom = NINT( hist_days_stom ) * one_day
646     WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
647             hist_dt_stom/one_day
648  ENDIF
649!-
650  !-
651  !- initialize
652  WRITE(numout,*) "before histbeg : ",date0,dt
653  IF ( rectilinear ) THEN
654#ifdef CPP_PARA
655     CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
656          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
657#else
658     CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
659          &     itau_dep, date0, dt, hori_id, hist_id_stom)
660#endif
661  ELSE
662#ifdef CPP_PARA
663     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
664          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
665#else
666     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
667          &     itau_dep, date0, dt, hori_id, hist_id_stom)
668#endif
669  ENDIF
670  !- define PFT axis
671  hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
672  !- declare this axis
673  CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
674       & '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
675!- define Pool_10 axis
676   hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
677!- declare this axis
678  CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
679       & '1', 10, hist_pool_10axis, hist_pool_10axis_id)
680
681!- define Pool_100 axis
682   hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
683!- declare this axis
684  CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
685       & '1', 100, hist_pool_100axis, hist_pool_100axis_id)
686
687!- define Pool_11 axis
688   hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
689!- declare this axis
690  CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
691       & '1', 11, hist_pool_11axis, hist_pool_11axis_id)
692
693!- define Pool_101 axis
694   hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
695!- declare this axis
696  CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
697       & '1', 101, hist_pool_101axis, hist_pool_101axis_id)
698
699  !- define STOMATE history file
700  CALL stom_define_history (hist_id_stom, nvm, iim, jjm, &
701       & dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
702 & hist_pool_10axis_id, hist_pool_100axis_id, &
703 & hist_pool_11axis_id, hist_pool_101axis_id)
704
705  !- end definition
706  CALL histend(hist_id_stom)
707  !-
708  !-
709  ! STOMATE IPCC OUTPUTS IS ACTIVATED
710  !-
711  !Config  Key  = STOMATE_IPCC_OUTPUT_FILE
712  !Config  Desc = Name of file in which STOMATE's output is going
713  !Config         to be written
714  !Config  Def  = stomate_ipcc_history.nc
715  !Config  Help = This file is going to be created by the model
716  !Config         and will contain the output from the model.
717  !Config         This file is a truly COADS compliant netCDF file.
718  !Config         It will be generated by the hist software from
719  !Config         the IOIPSL package.
720  !-
721  stom_ipcc_histname='stomate_ipcc_history.nc'
722  CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)       
723  WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
724  !-
725  !Config  Key  = STOMATE_IPCC_HIST_DT
726  !Config  Desc = STOMATE IPCC history time step (d)
727  !Config  Def  = 0.
728  !Config  Help = Time step of the STOMATE IPCC history file
729  !-
730  hist_days_stom_ipcc = zero
731  CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)       
732  IF ( hist_days_stom_ipcc == moins_un ) THEN
733     hist_dt_stom_ipcc = moins_un
734     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
735  ELSE
736     hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
737     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
738          hist_dt_stom_ipcc/one_day
739  ENDIF
740
741  ! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
742  dt_slow_ = one_day
743  CALL getin_p('DT_SLOW', dt_slow_)
744  IF ( hist_days_stom_ipcc > zero ) THEN
745     IF (dt_slow_ > hist_dt_stom_ipcc) THEN
746        WRITE(numout,*) "DT_SLOW = ",dt_slow_,"  , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
747        CALL ipslerr (3,'intsurf_history', &
748             &          'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
749             &          '(must be less or equal)')
750     ENDIF
751  ENDIF
752
753  IF ( hist_dt_stom_ipcc == 0 ) THEN
754     hist_id_stom_ipcc = -1
755  ELSE
756     !-
757     !- initialize
758     IF ( rectilinear ) THEN
759#ifdef CPP_PARA
760        CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
761             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
762#else
763        CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
764             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
765#endif
766     ELSE
767#ifdef CPP_PARA
768        CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
769             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
770#else
771        CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
772             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc)
773#endif
774     ENDIF
775     !- declare this axis
776     CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
777          & '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)
778
779     !- define STOMATE history file
780     CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
781          & dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)
782
783     !- end definition
784     CALL histend(hist_id_stom_IPCC)
785
786  ENDIF
787  !
788  CALL histwrite(hist_id_stom, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
789  IF ( hist_id_stom_IPCC > 0 ) THEN
790     CALL histwrite(hist_id_stom_IPCC, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
791  ENDIF
792  !
793  hist_id_sec = -1
794  hist_id_sec2 = -1
795!-
796! first call of slowproc to initialize variables
797!-
798  itau = itau_dep
799 
800  DO ji=1,kjpindex
801     DO jv=1,nvm
802        indexveg((jv-1)*kjpindex + ji) = indices(ji) + (jv-1)*kjpij
803     ENDDO
804  ENDDO
805  !-
806  !MM Problem here with dpu which depends on soil type           
807  DO l = 1, nbdl-1
808     ! first 2.0 is dpu
809     ! second 2.0 is average
810     diaglev(l) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(l-1) -1) + ( 2**(l) -1) ) / 2.0
811  ENDDO
812  diaglev(nbdl) = dpu_cste
813  !
814
815  ! >> DS : getin the value of slowproc parameters
816  !
817  !Config Key  = CLAYFRACTION_DEFAULT
818  !Config Desc =
819  !Config If   = OK_SECHIBA
820  !Config Def  = 0.2
821  !Config Help =
822  !Config Units = NONE   
823  CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
824  !
825  !Config Key  = MIN_VEGFRAC
826  !Config Desc = Minimal fraction of mesh a vegetation type can occupy
827  !Config If   = OK_SECHIBA
828  !Config Def  = 0.001
829  !Config Help =
830  !Config Units = NONE 
831  CALL getin_p('MIN_VEGFRAC',min_vegfrac)
832  !
833  !Config Key  = STEMPDIAG_BID
834  !Config Desc = only needed for an initial LAI if there is no restart file
835  !Config If   = OK_SECHIBA
836  !Config Def  = 280.
837  !Config Help =
838  !Config Units =
839  CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
840  !
841  !Config Key  = SOILTYPE_DEFAULT
842  !Config Desc = Default soil texture distribution in the following order : sand, loam and clay
843  !Config If   = OK_SECHIBA
844  !Config Def  = 0.0, 1.0, 0.0
845  !Config Help =
846  !Config Units = NONE   
847  CALL getin_p('SOILTYPE_DEFAULT',soiltype_default)
848  ! >> DS
849
850  CALL ioget_expval(val_exp)
851  ldrestart_read = .TRUE.
852  ldrestart_write = .FALSE.
853  !-
854  ! read some variables we need from SECHIBA's restart file
855  !-
856  CALL slowproc_main &
857 &  (itau, kjpij, kjpindex, dt_force, date0, &
858 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
859 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
860 &   t2m, t2m_min, temp_sol, soiltemp, &
861 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
862 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
863 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
864 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
865  ! correct date
866  day_counter = one_day - dt_force
867  date=1
868!-
869! time loop
870!-
871  IF (debug) check_time=.TRUE.
872  CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
873  l_first_intersurf=.FALSE.
874  !
875  DO itau = itau_dep+itau_step,itau_fin,itau_step
876    !
877    CALL intsurf_time( itau, date0, dtradia )
878     !
879!-- next forcing state
880    iisf = iisf+1
881    IF (debug) WRITE(numout,*) "itau,iisf : ",itau,iisf
882!---
883    IF (iisf .GT. nsfm) THEN
884!-----
885!---- we have to read new forcing states from the file
886!-----
887!---- determine blocks of forcing states that are contiguous in memory
888!-----
889        CALL forcing_read(forcing_id,nsfm)
890
891!--------------------------
892
893!-----
894!---- determine which forcing states must be read next time
895!-----
896      isf(1) = isf(nsfm)+1
897      IF ( isf(1) .GT. nsft ) isf(1) = 1
898        DO iiisf = 2, nsfm
899           isf(iiisf) = isf(iiisf-1)+1
900           IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
901        ENDDO
902        nf_written(isf(:)) = .TRUE.
903!---- start again at first forcing state
904        iisf = 1
905     ENDIF
906     ! Bug here ! soiltype(:,3) != clay
907!     soiltype(:,3) = clay_fm(:,iisf)
908     humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
909     litterhum(:) = litterhum_daily_fm(:,iisf)
910     t2m(:) = t2m_daily_fm(:,iisf)
911     t2m_min(:) = t2m_min_daily_fm(:,iisf)
912     temp_sol(:) = tsurf_daily_fm(:,iisf)
913     soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
914     soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
915     precip_rain(:) = precip_fm(:,iisf)
916     gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
917     veget_force_x(:,:) = veget_fm(:,:,iisf)
918     veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
919     lai_force_x(:,:) = lai_fm(:,:,iisf)
920     WHERE ( t2m(:) .LT. ZeroCelsius )
921        precip_snow(:) = precip_rain(:)
922        precip_rain(:) = 0.
923     ELSEWHERE
924        precip_snow(:) = 0.
925     ENDWHERE
926!---
927!-- scale GPP to new lai and veget_max
928!---
929     WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
930!-- scale GPP to new LAI
931     WHERE (lai_force_x(:,:) .GT. 0.0 )
932        gpp_x(:,:) = gpp_x(:,:)*ATAN(2.*lai_x(:,:)) &
933 &                           /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
934     ENDWHERE
935!-- scale GPP to new veget_max
936     WHERE (veget_max_force_x(:,:) .GT. 0.0 )
937        gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
938     ENDWHERE
939!---
940!-- number crunching
941!---
942     ldrestart_read = .FALSE.
943     ldrestart_write = .FALSE.
944     CALL slowproc_main &
945 &    (itau, kjpij, kjpindex, dt_force, date0, &
946 &     ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
947 &     indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
948 &     t2m, t2m_min, temp_sol, soiltemp, &
949 &     humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
950 &     deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
951 &     frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
952 &     rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
953     day_counter = one_day - dt_force
954  ENDDO ! end of the time loop
955!-
956! write restart files
957!-
958  IF (is_root_prc) THEN
959! first, read and write variables that are not managed otherwise
960     taboo_vars = &
961          &  '$lat$ $lon$ $lev$ $veget_year$ '// &
962          &  '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
963          &  '$lai$ $soiltype_frac$ $clay_frac$ '// &
964          &  '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
965!-
966     CALL ioget_vname(rest_id_sec, nbvar, varnames)
967!-
968     DO iv = 1, nbvar
969!-- check if the variable is to be written here
970        IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
971           IF (debug) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
972
973!---- get variable dimensions, especially 3rd dimension
974           CALL ioget_vdim &
975                &      (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
976           l1d = ALL(vardims(1:varnbdim) .EQ. 1)
977!---- read it
978           IF (l1d) THEN
979              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
980                   1,1,1,itau_dep,.TRUE.,var_1d)
981           ELSE
982              ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
983              IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
984              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
985                   nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
986                   "gather",nbp_glo,indices_g)
987           ENDIF
988!---- write it
989           IF (l1d) THEN
990              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
991                   1,1,1,itau_fin,var_1d)
992           ELSE
993              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
994                   nbp_glo,vardims(3),1,itau_fin,var_3d, &
995                   'scatter',nbp_glo,indices_g)
996              DEALLOCATE (var_3d)
997           ENDIF
998        ENDIF
999     ENDDO
1000  ENDIF
1001  CALL barrier_para()
1002
1003! call slowproc to write restart files
1004  ldrestart_read = .FALSE.
1005  ldrestart_write = .TRUE.
1006!-
1007  IF (debug) WRITE(numout,*) "Call slowproc for restart."
1008  CALL slowproc_main &
1009 &  (itau_fin, kjpij, kjpindex, dt_force, date0, &
1010 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1011 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
1012 &   t2m, t2m_min, temp_sol, soiltemp, &
1013 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1014 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
1015 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
1016 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
1017!-
1018! close files
1019!-
1020  IF (is_root_prc) THEN
1021     CALL restclo
1022     IF ( debug )  WRITE(numout,*) 'REST CLOSED'
1023  ENDIF
1024  CALL histclo
1025
1026  IF (is_root_prc) &
1027       ier = NF90_CLOSE (forcing_id)
1028
1029  IF (is_root_prc) THEN
1030     CALL getin_dump()
1031  ENDIF
1032#ifdef CPP_PARA
1033  CALL MPI_FINALIZE(ier)
1034#endif
1035  WRITE(numout,*) "End of teststomate."
1036
1037!---------------
1038END PROGRAM teststomate
1039!
1040!===
1041!
Note: See TracBrowser for help on using the repository browser.