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

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

Correct teststomate bug : double allocation of the variable fco2_lu

File size: 35.1 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!!$ end add DS
486  ALLOCATE (indexveg(kjpindex*nvm), stat=ier)
487  a_er = a_er .OR. (ier.NE.0)
488  ALLOCATE (soiltype(kjpindex,nstm),stat=ier)
489  a_er = a_er .OR. (ier.NE.0)
490  ALLOCATE (veget_x(kjpindex,nvm),stat=ier)
491  a_er = a_er .OR. (ier.NE.0)
492  ALLOCATE (totfrac_nobio(kjpindex),stat=ier)
493  a_er = a_er .OR. (ier.NE.0)
494  ALLOCATE (frac_nobio(kjpindex,nnobio),stat=ier)
495  a_er = a_er .OR. (ier.NE.0)
496  ALLOCATE (veget_max_x(kjpindex,nvm),stat=ier)
497  a_er = a_er .OR. (ier.NE.0)
498  ALLOCATE (lai_x(kjpindex,nvm),stat=ier)
499  a_er = a_er .OR. (ier.NE.0)
500  ALLOCATE (veget_force_x(kjpindex,nvm),stat=ier)
501  a_er = a_er .OR. (ier.NE.0)
502  ALLOCATE (veget_max_force_x(kjpindex,nvm),stat=ier)
503  a_er = a_er .OR. (ier.NE.0)
504  ALLOCATE (lai_force_x(kjpindex,nvm),stat=ier)
505  a_er = a_er .OR. (ier.NE.0)
506  ALLOCATE (t2m(kjpindex),stat=ier)
507  a_er = a_er .OR. (ier.NE.0)
508  ALLOCATE (t2m_min(kjpindex),stat=ier)
509  a_er = a_er .OR. (ier.NE.0)
510  ALLOCATE (temp_sol(kjpindex),stat=ier)
511  a_er = a_er .OR. (ier.NE.0)
512  ALLOCATE (soiltemp(kjpindex,nbdl),stat=ier)
513  a_er = a_er .OR. (ier.NE.0)
514  ALLOCATE (soilhum(kjpindex,nbdl),stat=ier)
515  a_er = a_er .OR. (ier.NE.0)
516  ALLOCATE (humrel_x(kjpindex,nvm),stat=ier)
517  a_er = a_er .OR. (ier.NE.0)
518  ALLOCATE (litterhum(kjpindex),stat=ier)
519  a_er = a_er .OR. (ier.NE.0)
520  ALLOCATE (precip_rain(kjpindex),stat=ier)
521  a_er = a_er .OR. (ier.NE.0)
522  ALLOCATE (precip_snow(kjpindex),stat=ier)
523  a_er = a_er .OR. (ier.NE.0)
524  ALLOCATE (gpp_x(kjpindex,nvm),stat=ier)
525  a_er = a_er .OR. (ier.NE.0)
526  ALLOCATE (deadleaf_cover(kjpindex),stat=ier)
527  a_er = a_er .OR. (ier.NE.0)
528  ALLOCATE (assim_param_x(kjpindex,nvm,npco2),stat=ier)
529  a_er = a_er .OR. (ier.NE.0)
530  ALLOCATE (height_x(kjpindex,nvm),stat=ier)
531  a_er = a_er .OR. (ier.NE.0)
532  ALLOCATE (qsintmax_x(kjpindex,nvm),stat=ier)
533  a_er = a_er .OR. (ier.NE.0)
534  ALLOCATE (co2_flux(kjpindex,nvm),stat=ier)
535  a_er = a_er .OR. (ier.NE.0)
536  ALLOCATE (fco2_lu(kjpindex),stat=ier)
537  a_er = a_er .OR. (ier.NE.0)
538  IF (a_er) THEN
539     CALL ipslerr (3,'teststomate', &
540          &        'PROBLEM WITH ALLOCATION', &
541          &        'for local variables 1','')
542  ENDIF
543  !
544  ! prepare forcing
545  !
546  max_totsize = 50
547  CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize)
548  max_totsize = max_totsize * 1000000
549
550  totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + &
551       SIZE(humrel_x)*KIND(humrel_x) + &
552       SIZE(litterhum)*KIND(litterhum) + &
553       SIZE(t2m)*KIND(t2m) + &
554       SIZE(t2m_min)*KIND(t2m_min) + &
555       SIZE(temp_sol)*KIND(temp_sol) + &
556       SIZE(soiltemp)*KIND(soiltemp) + &
557       SIZE(soilhum)*KIND(soilhum) + &
558       SIZE(precip_rain)*KIND(precip_rain) + &
559       SIZE(precip_snow)*KIND(precip_snow) + &
560       SIZE(gpp_x)*KIND(gpp_x) + &
561       SIZE(veget_force_x)*KIND(veget_force_x) + &
562       SIZE(veget_max_force_x)*KIND(veget_max_force_x) + &
563       SIZE(lai_force_x)*KIND(lai_force_x)
564  CALL reduce_sum(totsize_1step,totsize_tmp)
565  CALL bcast(totsize_tmp)
566  totsize_1step=totsize_tmp 
567
568! total number of forcing steps
569  IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN
570     CALL ipslerr (3,'teststomate', &
571          &        'stomate: error with total number of forcing steps', &
572          &        'nsft','teststomate computation different with forcing file value.')
573  ENDIF
574! number of forcing steps in memory
575  nsfm = MIN(nsft, &
576       &       MAX(1,NINT( REAL(max_totsize,r_std) &
577       &                  /REAL(totsize_1step,r_std))))
578!-
579  WRITE(numout,*) 'Offline forcing of Stomate:'
580  WRITE(numout,*) '  Total number of forcing steps:',nsft
581  WRITE(numout,*) '  Number of forcing steps in memory:',nsfm
582!-
583  CALL init_forcing(kjpindex,nsfm,nsft)
584  !-
585! ensure that we read all new forcing states
586  iisf = nsfm
587! initialize that table that contains the indices
588! of the forcing states that will be in memory
589  isf(:) = (/ (i,i=1,nsfm) /)
590
591  nf_written(:) = .FALSE.
592  nf_written(isf(:)) = .TRUE.
593
594!-
595! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA
596!-
597  itau_step = NINT(dt_force/dtradia)
598  IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step
599  !
600  CALL ioconf_startdate(date0)
601  CALL intsurf_time( itau_dep, date0, dtradia )
602!-
603! Length of integration
604!-
605  WRITE(time_str,'(a)') '1Y'
606  CALL getin_p ('TIME_LENGTH', time_str)
607! transform into itau
608  CALL tlen2itau(time_str, dt, date0, itau_len)
609! itau_len*dtradia must be a multiple of dt_force
610  itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) &
611       &                *dt_force/dtradia)
612  !-
613  itau_fin = itau_dep+itau_len
614  !-
615  ! set up STOMATE history file
616  !-
617  !Config  Key  = STOMATE_OUTPUT_FILE
618  !Config  Desc = Name of file in which STOMATE's output is going
619  !Config         to be written
620  !Config  Def  = stomate_history.nc
621  !Config  Help = This file is going to be created by the model
622  !Config         and will contain the output from the model.
623  !Config         This file is a truly COADS compliant netCDF file.
624  !Config         It will be generated by the hist software from
625  !Config         the IOIPSL package.
626  !-
627  stom_histname='stomate_history.nc'
628  CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname)
629  WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)
630  !-
631  !Config  Key  = STOMATE_HIST_DT
632  !Config  Desc = STOMATE history time step (d)
633  !Config  Def  = 10.
634  !Config  Help = Time step of the STOMATE history file
635  !-
636  hist_days_stom = 10.
637  CALL getin_p ('STOMATE_HIST_DT', hist_days_stom)
638  IF ( hist_days_stom == -1. ) THEN
639     hist_dt_stom = -1.
640     WRITE(numout,*) 'output frequency for STOMATE history file (d): one month.'
641  ELSE
642     hist_dt_stom = NINT( hist_days_stom ) * one_day
643     WRITE(numout,*) 'output frequency for STOMATE history file (d): ', &
644             hist_dt_stom/one_day
645  ENDIF
646!-
647  !-
648  !- initialize
649  WRITE(numout,*) "before histbeg : ",date0,dt
650  IF ( rectilinear ) THEN
651#ifdef CPP_PARA
652     CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
653          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
654#else
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)
657#endif
658  ELSE
659#ifdef CPP_PARA
660     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
661          &     itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id)
662#else
663     CALL histbeg(stom_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
664          &     itau_dep, date0, dt, hori_id, hist_id_stom)
665#endif
666  ENDIF
667  !- define PFT axis
668  hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /)
669  !- declare this axis
670  CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', &
671       & '1', nvm, hist_PFTaxis, hist_PFTaxis_id)
672!- define Pool_10 axis
673   hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /)
674!- declare this axis
675  CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', &
676       & '1', 10, hist_pool_10axis, hist_pool_10axis_id)
677
678!- define Pool_100 axis
679   hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /)
680!- declare this axis
681  CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', &
682       & '1', 100, hist_pool_100axis, hist_pool_100axis_id)
683
684!- define Pool_11 axis
685   hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /)
686!- declare this axis
687  CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', &
688       & '1', 11, hist_pool_11axis, hist_pool_11axis_id)
689
690!- define Pool_101 axis
691   hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /)
692!- declare this axis
693  CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', &
694       & '1', 101, hist_pool_101axis, hist_pool_101axis_id)
695
696  !- define STOMATE history file
697  CALL stom_define_history (hist_id_stom, nvm, iim, jjm, &
698       & dt, hist_dt_stom, hori_id, hist_PFTaxis_id, &
699 & hist_pool_10axis_id, hist_pool_100axis_id, &
700 & hist_pool_11axis_id, hist_pool_101axis_id)
701
702  !- end definition
703  CALL histend(hist_id_stom)
704  !-
705  !-
706  ! STOMATE IPCC OUTPUTS IS ACTIVATED
707  !-
708  !Config  Key  = STOMATE_IPCC_OUTPUT_FILE
709  !Config  Desc = Name of file in which STOMATE's output is going
710  !Config         to be written
711  !Config  Def  = stomate_ipcc_history.nc
712  !Config  Help = This file is going to be created by the model
713  !Config         and will contain the output from the model.
714  !Config         This file is a truly COADS compliant netCDF file.
715  !Config         It will be generated by the hist software from
716  !Config         the IOIPSL package.
717  !-
718  stom_ipcc_histname='stomate_ipcc_history.nc'
719  CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname)       
720  WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname)
721  !-
722  !Config  Key  = STOMATE_IPCC_HIST_DT
723  !Config  Desc = STOMATE IPCC history time step (d)
724  !Config  Def  = 0.
725  !Config  Help = Time step of the STOMATE IPCC history file
726  !-
727  hist_days_stom_ipcc = zero
728  CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc)       
729  IF ( hist_days_stom_ipcc == moins_un ) THEN
730     hist_dt_stom_ipcc = moins_un
731     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.'
732  ELSE
733     hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day
734     WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', &
735          hist_dt_stom_ipcc/one_day
736  ENDIF
737
738  ! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters
739  dt_slow_ = one_day
740  CALL getin_p('DT_SLOW', dt_slow_)
741  IF ( hist_days_stom_ipcc > zero ) THEN
742     IF (dt_slow_ > hist_dt_stom_ipcc) THEN
743        WRITE(numout,*) "DT_SLOW = ",dt_slow_,"  , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc
744        CALL ipslerr (3,'intsurf_history', &
745             &          'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', &
746             &          '(must be less or equal)')
747     ENDIF
748  ENDIF
749
750  IF ( hist_dt_stom_ipcc == 0 ) THEN
751     hist_id_stom_ipcc = -1
752  ELSE
753     !-
754     !- initialize
755     IF ( rectilinear ) THEN
756#ifdef CPP_PARA
757        CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect,  1, iim, 1, jjm, &
758             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
759#else
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)
762#endif
763     ELSE
764#ifdef CPP_PARA
765        CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat,  1, iim, 1, jjm, &
766             &     itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id)
767#else
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)
770#endif
771     ENDIF
772     !- declare this axis
773     CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', &
774          & '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id)
775
776     !- define STOMATE history file
777     CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, &
778          & dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id)
779
780     !- end definition
781     CALL histend(hist_id_stom_IPCC)
782
783  ENDIF
784  !
785  CALL histwrite(hist_id_stom, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
786  IF ( hist_id_stom_IPCC > 0 ) THEN
787     CALL histwrite(hist_id_stom_IPCC, 'Areas',  itau_dep+itau_step, area, kjpindex, indices)
788  ENDIF
789  !
790  hist_id_sec = -1
791  hist_id_sec2 = -1
792!-
793! first call of slowproc to initialize variables
794!-
795  itau = itau_dep
796 
797  DO ji=1,kjpindex
798     DO jv=1,nvm
799        indexveg((jv-1)*kjpindex + ji) = indices(ji) + (jv-1)*kjpij
800     ENDDO
801  ENDDO
802  !-
803  !MM Problem here with dpu which depends on soil type           
804  DO l = 1, nbdl-1
805     ! first 2.0 is dpu
806     ! second 2.0 is average
807     diaglev(l) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(l-1) -1) + ( 2**(l) -1) ) / 2.0
808  ENDDO
809  diaglev(nbdl) = dpu_cste
810  !
811
812  ! >> DS : getin the value of slowproc parameters
813  !
814  !Config Key  = CLAYFRACTION_DEFAULT
815  !Config Desc =
816  !Config If   = OK_SECHIBA
817  !Config Def  = 0.2
818  !Config Help =
819  !Config Units = NONE   
820  CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default)
821  !
822  !Config Key  = MIN_VEGFRAC
823  !Config Desc = Minimal fraction of mesh a vegetation type can occupy
824  !Config If   = OK_SECHIBA
825  !Config Def  = 0.001
826  !Config Help =
827  !Config Units = NONE 
828  CALL getin_p('MIN_VEGFRAC',min_vegfrac)
829  !
830  !Config Key  = STEMPDIAG_BID
831  !Config Desc = only needed for an initial LAI if there is no restart file
832  !Config If   = OK_SECHIBA
833  !Config Def  = 280.
834  !Config Help =
835  !Config Units =
836  CALL getin_p('STEMPDIAG_BID',stempdiag_bid)
837  !
838  !Config Key  = SOILTYPE_DEFAULT
839  !Config Desc = Default soil texture distribution in the following order : sand, loam and clay
840  !Config If   = OK_SECHIBA
841  !Config Def  = 0.0, 1.0, 0.0
842  !Config Help =
843  !Config Units = NONE   
844  CALL getin_p('SOILTYPE_DEFAULT',soiltype_default)
845  ! >> DS
846
847  CALL ioget_expval(val_exp)
848  ldrestart_read = .TRUE.
849  ldrestart_write = .FALSE.
850  !-
851  ! read some variables we need from SECHIBA's restart file
852  !-
853  CALL slowproc_main &
854 &  (itau, kjpij, kjpindex, dt_force, date0, &
855 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
856 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
857 &   t2m, t2m_min, temp_sol, soiltemp, &
858 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
859 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
860 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
861 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
862  ! correct date
863  day_counter = one_day - dt_force
864  date=1
865!-
866! time loop
867!-
868  IF (debug) check_time=.TRUE.
869  CALL intsurf_time( itau_dep+itau_step, date0, dtradia )
870  l_first_intersurf=.FALSE.
871  !
872  DO itau = itau_dep+itau_step,itau_fin,itau_step
873    !
874    CALL intsurf_time( itau, date0, dtradia )
875     !
876!-- next forcing state
877    iisf = iisf+1
878    IF (debug) WRITE(numout,*) "itau,iisf : ",itau,iisf
879!---
880    IF (iisf .GT. nsfm) THEN
881!-----
882!---- we have to read new forcing states from the file
883!-----
884!---- determine blocks of forcing states that are contiguous in memory
885!-----
886        CALL forcing_read(forcing_id,nsfm)
887
888!--------------------------
889
890!-----
891!---- determine which forcing states must be read next time
892!-----
893      isf(1) = isf(nsfm)+1
894      IF ( isf(1) .GT. nsft ) isf(1) = 1
895        DO iiisf = 2, nsfm
896           isf(iiisf) = isf(iiisf-1)+1
897           IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1
898        ENDDO
899        nf_written(isf(:)) = .TRUE.
900!---- start again at first forcing state
901        iisf = 1
902     ENDIF
903     ! Bug here ! soiltype(:,3) != clay
904!     soiltype(:,3) = clay_fm(:,iisf)
905     humrel_x(:,:) = humrel_daily_fm(:,:,iisf)
906     litterhum(:) = litterhum_daily_fm(:,iisf)
907     t2m(:) = t2m_daily_fm(:,iisf)
908     t2m_min(:) = t2m_min_daily_fm(:,iisf)
909     temp_sol(:) = tsurf_daily_fm(:,iisf)
910     soiltemp(:,:) = tsoil_daily_fm(:,:,iisf)
911     soilhum(:,:) = soilhum_daily_fm(:,:,iisf)
912     precip_rain(:) = precip_fm(:,iisf)
913     gpp_x(:,:) = gpp_daily_fm(:,:,iisf)
914     veget_force_x(:,:) = veget_fm(:,:,iisf)
915     veget_max_force_x(:,:) = veget_max_fm(:,:,iisf)
916     lai_force_x(:,:) = lai_fm(:,:,iisf)
917     WHERE ( t2m(:) .LT. ZeroCelsius )
918        precip_snow(:) = precip_rain(:)
919        precip_rain(:) = 0.
920     ELSEWHERE
921        precip_snow(:) = 0.
922     ENDWHERE
923!---
924!-- scale GPP to new lai and veget_max
925!---
926     WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0
927!-- scale GPP to new LAI
928     WHERE (lai_force_x(:,:) .GT. 0.0 )
929        gpp_x(:,:) = gpp_x(:,:)*ATAN(2.*lai_x(:,:)) &
930 &                           /ATAN( 2.*MAX(lai_force_x(:,:),0.01))
931     ENDWHERE
932!-- scale GPP to new veget_max
933     WHERE (veget_max_force_x(:,:) .GT. 0.0 )
934        gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)
935     ENDWHERE
936!---
937!-- number crunching
938!---
939     ldrestart_read = .FALSE.
940     ldrestart_write = .FALSE.
941     CALL slowproc_main &
942 &    (itau, kjpij, kjpindex, dt_force, date0, &
943 &     ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
944 &     indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
945 &     t2m, t2m_min, temp_sol, soiltemp, &
946 &     humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
947 &     deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
948 &     frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
949 &     rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
950     day_counter = one_day - dt_force
951  ENDDO ! end of the time loop
952!-
953! write restart files
954!-
955  IF (is_root_prc) THEN
956! first, read and write variables that are not managed otherwise
957     taboo_vars = &
958          &  '$lat$ $lon$ $lev$ $veget_year$ '// &
959          &  '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// &
960          &  '$lai$ $soiltype_frac$ $clay_frac$ '// &
961          &  '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$'
962!-
963     CALL ioget_vname(rest_id_sec, nbvar, varnames)
964!-
965     DO iv = 1, nbvar
966!-- check if the variable is to be written here
967        IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN
968           IF (debug) WRITE(numout,*) "restart var : ",TRIM(varnames(iv)),itau_dep,itau_fin
969
970!---- get variable dimensions, especially 3rd dimension
971           CALL ioget_vdim &
972                &      (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)
973           l1d = ALL(vardims(1:varnbdim) .EQ. 1)
974!---- read it
975           IF (l1d) THEN
976              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
977                   1,1,1,itau_dep,.TRUE.,var_1d)
978           ELSE
979              ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier)
980              IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'
981              CALL restget (rest_id_sec,TRIM(varnames(iv)), &
982                   nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, &
983                   "gather",nbp_glo,indices_g)
984           ENDIF
985!---- write it
986           IF (l1d) THEN
987              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
988                   1,1,1,itau_fin,var_1d)
989           ELSE
990              CALL restput (rest_id_sec,TRIM(varnames(iv)), &
991                   nbp_glo,vardims(3),1,itau_fin,var_3d, &
992                   'scatter',nbp_glo,indices_g)
993              DEALLOCATE (var_3d)
994           ENDIF
995        ENDIF
996     ENDDO
997  ENDIF
998  CALL barrier_para()
999
1000! call slowproc to write restart files
1001  ldrestart_read = .FALSE.
1002  ldrestart_write = .TRUE.
1003!-
1004  IF (debug) WRITE(numout,*) "Call slowproc for restart."
1005  CALL slowproc_main &
1006 &  (itau_fin, kjpij, kjpindex, dt_force, date0, &
1007 &   ldrestart_read, ldrestart_write, .FALSE., .TRUE., &
1008 &   indices, indexveg, lalo, neighbours, resolution, contfrac, soiltype, &
1009 &   t2m, t2m_min, temp_sol, soiltemp, &
1010 &   humrel_x, soilhum, litterhum, precip_rain, precip_snow, gpp_x, &
1011 &   deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, &
1012 &   frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, &
1013 &   rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu)
1014!-
1015! close files
1016!-
1017  IF (is_root_prc) THEN
1018     CALL restclo
1019     IF ( debug )  WRITE(numout,*) 'REST CLOSED'
1020  ENDIF
1021  CALL histclo
1022
1023  IF (is_root_prc) &
1024       ier = NF90_CLOSE (forcing_id)
1025
1026  IF (is_root_prc) THEN
1027     CALL getin_dump()
1028  ENDIF
1029#ifdef CPP_PARA
1030  CALL MPI_FINALIZE(ier)
1031#endif
1032  WRITE(numout,*) "End of teststomate."
1033
1034!---------------
1035END PROGRAM teststomate
1036!
1037!===
1038!
Note: See TracBrowser for help on using the repository browser.