source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/orchideedriver.f90 @ 7262

Last change on this file since 7262 was 7262, checked in by agnes.ducharne, 6 months ago

As in r7065 and r76158 for the new driver. Runs OK on jean-zay with the old driver, but not yet tested with the new driver. The corresponding configurations need to be added following r7081, r7082, and r7719.

File size: 40.3 KB
Line 
1! =================================================================================================================================
2! PROGRAM       : orchideedriver
3!
4! CONTACT       : jan.polcher@lmd.jussieu.fr
5!
6! LICENCE      : IPSL (2016)
7! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8!
9!>\BRIEF      This is the main program for the new driver. This only organises the data and calls sechiba_main.
10!!            The main work is done in glogrid.f90 and forcing_tools.f90.
11!!
12!!\n DESCRIPTION: Call the various modules to get the forcing data and provide it to SECHIBA. The only complexity
13!!                is setting-up the domain decomposition and distributing the grid information.
14!!                The code is parallel from tip to toe using the domain decomposition inherited from LMDZ.
15!!
16!! RECENT CHANGE(S): None
17!!
18!! REFERENCE(S) :
19!!
20!! SVN          :
21!! $HeadURL:  $
22!! $Date:  $
23!! $Revision: $
24!! \n
25!_
26!================================================================================================================================
27!
28PROGRAM orchidedriver
29  !---------------------------------------------------------------------
30  !-
31  !-
32  !---------------------------------------------------------------------
33  USE defprec
34  USE netcdf
35  !
36  !
37  USE ioipsl_para
38  USE mod_orchidee_para
39  !
40  USE grid
41  USE time
42  USE timer
43  USE constantes
44  USE constantes_soil
45  USE forcing_tools
46  USE globgrd
47  !
48  USE sechiba
49  USE control
50  USE ioipslctrl
51  USE xios_orchidee
52 
53  !
54  !-
55  IMPLICIT NONE
56  !-
57  CHARACTER(LEN=80) :: gridfilename
58  CHARACTER(LEN=80), DIMENSION(100) :: forfilename
59  INTEGER(i_std) :: nb_forcefile
60  CHARACTER(LEN=8)  :: model_guess
61  INTEGER(i_std)    :: iim_glo, jjm_glo, file_id
62  !-
63  INTEGER(i_std)    :: nbseg
64  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon_glo, lat_glo, area_glo
65  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: mask_glo
66  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: corners_glo
67  REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: corners_lon, corners_lat
68  INTEGER(i_std) :: nbindex_g, kjpindex
69  INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: kindex, kindex_g
70  REAL(r_std), DIMENSION(2) :: zoom_lon, zoom_lat
71  !
72  ! Variables for the global grid available on all procs and used
73  ! to fill the ORCHIDEE variable on the root_proc
74  !
75  REAL(r_std), ALLOCATABLE, DIMENSION(:,:)    :: lalo_glo
76  REAL(r_std), ALLOCATABLE, DIMENSION(:)      :: contfrac_glo
77  CHARACTER(LEN=20)                           :: calendar
78  !-
79  !- Variables local to each processors.
80  !-
81  INTEGER(i_std) :: i, j, ik, in, nbdt, first_point
82  INTEGER(i_std) :: itau, itau_offset, itau_sechiba
83  REAL(r_std)    :: date0, date0_shifted, dt, julian
84  REAL(r_std)    :: date0_tmp, dt_tmp
85  INTEGER(i_std) :: nbdt_tmp
86  REAL(r_std)    :: timestep_interval(2), timestep_int_next(2)
87  !
88  INTEGER(i_std) :: rest_id, rest_id_stom
89  INTEGER(i_std) ::  hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC
90  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lalo_loc
91  INTEGER(i_std) :: iim, jjm, ier
92  REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lon, lat
93  !-
94  !- input fields
95  !-
96  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: u             !! Lowest level wind speed
97  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: v             !! Lowest level wind speed
98  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: zlev_uv       !! Height of first layer
99  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: zlev_tq       !! Height of first layer
100  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: qair          !! Lowest level specific humidity
101  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: precip_rain   !! Rain precipitation
102  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: precip_snow   !! Snow precipitation
103  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: lwdown        !! Down-welling long-wave flux
104  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: swdown        !! Downwelling surface short-wave flux
105  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: sinang        !! cosine of solar zenith angle
106  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: temp_air      !! Air temperature in Kelvin
107  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: epot_air      !! Air potential energy
108  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: ccanopy       !! CO2 concentration in the canopy
109  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: petAcoef      !! Coeficients A from the PBL resolution
110  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: peqAcoef      !! One for T and another for q
111  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: petBcoef      !! Coeficients B from the PBL resolution
112  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: peqBcoef      !! One for T and another for q
113  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: cdrag         !! Cdrag
114  REAL(r_std), ALLOCATABLE, DIMENSION (:)             :: pb            !! Lowest level pressure
115  !-
116  !- output fields
117  !-
118  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0m           !! Surface roughness for momentum (m)
119  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: z0h           !! Surface roughness for heat (m)
120  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: coastalflow   !! Diffuse flow of water into the ocean (m^3/dt)
121  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: riverflow     !! Largest rivers flowing into the ocean (m^3/dt)
122  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: tsol_rad      !! Radiative surface temperature
123  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: vevapp        !! Total of evaporation
124  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: temp_sol_new  !! New soil temperature
125  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: qsurf         !! Surface specific humidity
126  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: albedo        !! Albedo
127  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: fluxsens      !! Sensible chaleur flux
128  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: fluxlat       !! Latent chaleur flux
129  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: emis          !! Emissivity
130  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: netco2        !! netco2flux: Sum CO2 flux over PFTs
131  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: carblu        !! fco2_lu: Land Cover Change CO2 flux
132  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: carbwh        !! fco2_wh: Wood harvest CO2 flux
133  REAL(r_std), ALLOCATABLE, DIMENSION (:)            :: carbha        !! fco2_ha: Crop harvest CO2 flux
134  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: veget_diag    !! Fraction of vegetation type (unitless, 0-1)
135  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: lai_diag      !! Leaf area index (m^2 m^{-2}
136  REAL(r_std), ALLOCATABLE, DIMENSION (:,:)          :: height_diag   !! Vegetation Height (m)
137  !-
138  !-
139  !-
140  REAL(r_std) :: atmco2
141  REAL(r_std), ALLOCATABLE, DIMENSION (:)  :: u_tq, v_tq, swnet
142  LOGICAL :: lrestart_read = .TRUE. !! Logical for _restart_ file to read
143  LOGICAL :: lrestart_write = .FALSE. !! Logical for _restart_ file to write'
144  !
145  ! Timer variables
146  !
147  LOGICAL, PARAMETER :: timemeasure=.TRUE.
148  REAL(r_std) :: waitput_cputime=0.0, waitget_cputime=0.0, orchidee_cputime=0.0
149  REAL(r_std) :: waitput_walltime=0.0, waitget_walltime=0.0, orchidee_walltime=0.0
150  !
151  !
152  ! Print point
153  !
154!!  REAL(r_std), DIMENSION(2) :: testpt=(/44.8,-25.3/)
155!!  REAL(r_std), DIMENSION(2) :: testpt=(/44.8,-18.3/)
156!!  REAL(r_std), DIMENSION(2) :: testpt=(/-60.25,-5.25/)
157!!  REAL(r_std), DIMENSION(2) :: testpt=(/46.7,10.3/)
158!!  REAL(r_std), DIMENSION(2) :: testpt=(/0.25,49.25/)
159  ! Case when no ouput is desired.
160  REAL(r_std), DIMENSION(2) :: testpt=(/9999.99,9999.99/)
161  INTEGER(i_std) :: ktest
162  INTEGER :: printlev_loc  !! local write level
163
164  OFF_LINE_MODE = .TRUE. 
165
166
167  !-
168  !---------------------------------------------------------------------------------------
169  !-
170  !- Define MPI communicator
171  !-
172  !---------------------------------------------------------------------------------------
173  !-
174  !
175  ! Set parallel processing in ORCHIDEE
176  !
177  CALL Init_orchidee_para()
178  !
179  !====================================================================================
180  !
181  ! Start timer now that the paralelisation is initialized.
182  !
183  IF ( timemeasure ) THEN
184     CALL init_timer
185     CALL start_timer(timer_global)
186     CALL start_timer(timer_mpi)
187  ENDIF
188  !
189  !
190  !---------------------------------------------------------------------------------------
191  !-
192  !- Start the getconf processes
193  !-
194  !---------------------------------------------------------------------------------------
195  !-
196
197  ! Set specific write level to orchideedriver using PRINTLEV_orchideedriver=[0-4]
198  ! in run.def. The global printlev is used as default value.
199  printlev_loc=get_printlev('orchideedriver')
200
201  !-
202  !Config Key   = GRID_FILE
203  !Config Desc  = Name of file containing the forcing data
204  !Config If    = [-]
205  !Config Def   = grid_file.nc
206  !Config Help  = This is the name of the file from which we will read
207  !Config         or write into it the description of the grid from
208  !Config         the forcing file.
209  !Config         compliant.
210  !Config Units = [FILE]
211  !-
212  gridfilename='NONE'
213  CALL getin_p('GRID_FILE', gridfilename)
214  !-
215  forfilename(:)=" "
216  forfilename(1)='forcing_file.nc'
217  CALL getin_p('FORCING_FILE', forfilename)
218  !-
219  !- Define the zoom
220  !-
221  zoom_lon=(/-180,180/)
222  zoom_lat=(/-90,90/)
223  !
224  !Config Key   = LIMIT_WEST
225  !Config Desc  = Western limit of region
226  !Config If    = [-]
227  !Config Def   = -180.
228  !Config Help  = Western limit of the region we are
229  !Config         interested in. Between -180 and +180 degrees
230  !Config         The model will use the smalest regions from
231  !Config         region specified here and the one of the forcing file.
232  !Config Units = [Degrees]
233  !-
234  CALL getin_p('LIMIT_WEST',zoom_lon(1))
235  !-
236  !Config Key   = LIMIT_EAST
237  !Config Desc  = Eastern limit of region
238  !Config If    = [-]
239  !Config Def   = 180.
240  !Config Help  = Eastern limit of the region we are
241  !Config         interested in. Between -180 and +180 degrees
242  !Config         The model will use the smalest regions from
243  !Config         region specified here and the one of the forcing file.
244  !Config Units = [Degrees]
245  !-
246  CALL getin_p('LIMIT_EAST',zoom_lon(2))
247  !-
248  !Config Key   = LIMIT_NORTH
249  !Config Desc  = Northern limit of region
250  !Config If    = [-]
251  !Config Def   = 90.
252  !Config Help  = Northern limit of the region we are
253  !Config         interested in. Between +90 and -90 degrees
254  !Config         The model will use the smalest regions from
255  !Config         region specified here and the one of the forcing file.
256  !Config Units = [Degrees]
257  !-
258  CALL getin_p('LIMIT_NORTH',zoom_lat(2))
259  !-
260  !Config Key   = LIMIT_SOUTH
261  !Config Desc  = Southern limit of region
262  !Config If    = [-]
263  !Config Def   = -90.
264  !Config Help  = Southern limit of the region we are
265  !Config         interested in. Between 90 and -90 degrees
266  !Config         The model will use the smalest regions from
267  !Config         region specified here and the one of the forcing file.
268  !Config Units = [Degrees]
269  !-
270  CALL getin_p('LIMIT_SOUTH',zoom_lat(1))
271  IF ( (zoom_lon(1)+180 < EPSILON(zoom_lon(1))) .AND. (zoom_lon(2)-180 < EPSILON(zoom_lon(2))) .AND.&
272       &(zoom_lat(1)+90 < EPSILON(zoom_lat(1))) .AND. (zoom_lat(2)-90 < EPSILON(zoom_lat(2))) ) THEN
273
274     ! We are here only if zoom_lon and zoom_lat have there original values which
275     ! means that they have not been modified by the getin LIMIT_ above.
276     ! Read WEST_EAST and SOUTH_NORTH from run.def.
277     
278     !Config Key   = WEST_EAST
279     !Config Desc  = Longitude interval to use from the forcing data
280     !Config If    = [-]
281     !Config Def   = -180, 180
282     !Config Help  = This function allows to zoom into the forcing data
283     !Config Units = [degrees east]
284     !-
285     CALL getin_p('WEST_EAST', zoom_lon)
286     !
287     !Config Key   = SOUTH_NORTH
288     !Config Desc  = Latitude interval to use from the forcing data
289     !Config If    = [-]
290     !Config Def   = -90, 90
291     !Config Help  = This function allows to zoom into the forcing data
292     !Config Units = [degrees north]
293     !-
294     CALL getin_p('SOUTH_NORTH', zoom_lat)
295  ENDIF
296  !-
297  !-
298  !- Get some basic variables from the run.def
299  !-
300  atmco2=350.
301  CALL getin_p('ATM_CO2',atmco2)
302  !
303  !====================================================================================
304  !-
305  !-
306  !- Get the grid on all processors.
307  !-
308  !---------------------------------------------------------------------------------------
309  !-
310  !- Read the grid, only on the root proc. from the forcing file using tools in the globgrd module.
311  !- The grid is then broadcast to all other broadcast.
312  !
313  nb_forcefile = 0
314  DO ik=1,100
315     IF ( INDEX(forfilename(ik), '.nc') > 0 ) nb_forcefile = nb_forcefile+1
316  ENDDO
317  !
318  IF ( is_root_prc) THEN
319     CALL globgrd_getdomsz(gridfilename, iim_glo, jjm_glo, nbindex_g, model_guess, file_id, forfilename, zoom_lon, zoom_lat)
320     write(numout,*) "nbindex_g after calling globgrd_getdomsz in orchideedriver", nbindex_g
321     nbseg = 4
322  ENDIF
323  !
324  CALL bcast(iim_glo)
325  CALL bcast(jjm_glo)
326  CALL bcast(nbindex_g)
327  CALL bcast(nbseg)
328  !-
329  !- Allocation of memory
330  !- variables over the entire grid (thus in x,y)
331  ALLOCATE(lon_glo(iim_glo, jjm_glo))
332  ALLOCATE(lat_glo(iim_glo, jjm_glo))
333  ALLOCATE(mask_glo(iim_glo, jjm_glo))
334  ALLOCATE(area_glo(iim_glo, jjm_glo))
335  ALLOCATE(corners_glo(iim_glo, jjm_glo, nbseg, 2))
336  !
337  ! Gathered variables
338  ALLOCATE(kindex_g(nbindex_g))
339  ALLOCATE(contfrac_glo(nbindex_g))
340  !-
341  IF ( is_root_prc) THEN
342     CALL globgrd_getgrid(file_id, iim_glo, jjm_glo, nbindex_g, model_guess, &
343          &               lon_glo, lat_glo, mask_glo, area_glo, corners_glo,&
344          &               kindex_g, contfrac_glo, calendar)
345  ENDIF
346 
347  !
348  CALL bcast(lon_glo)
349  CALL bcast(lat_glo)
350  CALL bcast(mask_glo)
351  CALL bcast(area_glo)
352  CALL bcast(corners_glo)
353  CALL bcast(kindex_g)
354  CALL bcast(contfrac_glo)
355  CALL bcast(calendar)
356  CALL bcast(model_guess)
357  !
358  ALLOCATE(lalo_glo(nbindex_g,2))
359  DO ik=1,nbindex_g
360     !
361     j = ((kindex_g(ik)-1)/iim_glo)+1
362     i = (kindex_g(ik)-(j-1)*iim_glo)
363     !
364     IF ( i > iim_glo .OR. j > jjm_glo ) THEN
365        WRITE(100+mpi_rank,*) "Error in the indexing (ik, kindex, i, j) : ", ik, kindex(ik), i, j
366        STOP "ERROR in orchideedriver"
367     ENDIF
368     !
369     lalo_glo(ik,1) = lat_glo(i,j)
370     lalo_glo(ik,2) = lon_glo(i,j)
371     !
372  ENDDO
373  !
374 
375  WRITE(*,*) "Rank", mpi_rank, " Before parallel region All land points : ",  nbindex_g
376  WRITE(*,*) "Rank", mpi_rank, " from ", iim_glo, " point in Lon. and ", jjm_glo, "in Lat."
377  !-
378  !---------------------------------------------------------------------------------------
379  !-
380  !- Now that the grid is distributed on all procs we can start
381  !- initialise the ORCHIDEE domain on each proc (longitude, latitude, indices)
382  !-
383  !---------------------------------------------------------------------------------------
384  !-
385  !- init_data_para also transfers kindex_g to index_g (the variable used in ORCHIDEE)
386  !-
387  CALL grid_set_glo(iim_glo, jjm_glo, nbindex_g)
388 
389
390  CALL grid_allocate_glo(nbseg)
391  ! Copy the list of indexes of land points into index_g used by ORCHIDEE and then broacast to all
392  ! processors
393  CALL bcast(nbindex_g)
394  IF ( is_root_prc) index_g = kindex_g
395  CALL bcast(index_g)
396  !
397  WRITE(numout,*) "Rank", mpi_rank, "Into Init_orchidee_data_para_driver with ", nbindex_g
398  WRITE(numout,*) "Rank", mpi_rank, "Into ", index_g(1), index_g(nbindex_g)
399  !
400  CALL Init_orchidee_data_para_driver(nbindex_g,index_g)
401  CALL init_ioipsl_para 
402  !
403  WRITE(numout,*) "Rank", mpi_rank, "After init_data_para global size : ",  nbp_glo, SIZE(index_g), iim_g, iim_glo, jjm_g, jjm_glo
404  WRITE(numout,'("After init_data_para local : ij_nb, jj_nb",2I4)') iim_glo, jj_nb
405  !
406  ! Allocate grid on the local processor
407  !
408  IF ( model_guess == "regular") THEN
409     CALL grid_init (nbp_loc, nbseg, regular_lonlat, "ForcingGrid")
410  ELSE IF ( model_guess == "WRF") THEN
411     CALL grid_init (nbp_loc, nbseg, regular_xy, "WRFGrid")
412  ELSE
413     CALL ipslerr(3, "orchidedriver", "The grid found in the GRID_FILE is not supported by ORCHIDEE", "", "")
414  ENDIF
415  !
416  !
417  ! Transfer the global grid variables to the ORCHIDEE version on the root proc
418  ! *_glo -> *_g
419  ! Variables *_g were allocated with the CALL init_grid
420  !
421  !
422  lalo_g(:,:) = lalo_glo(:,:)
423  contfrac_g(:) = contfrac_glo(:)
424  lon_g(:,:) = lon_glo(:,:)
425  lat_g(:,:) = lat_glo(:,:)
426  !
427  !
428  ! Set the local dimensions of the fields
429  !
430  iim = iim_glo
431  jjm = jj_nb
432  kjpindex = nbp_loc
433  !
434  WRITE(numout,*) mpi_rank, "DIMENSIONS of grid on processor : iim, jjm, kjpindex = ", iim, jjm, kjpindex
435  !
436  !  Allocate the local arrays we need :
437  !
438  ALLOCATE(lon(iim,jjm), lat(iim,jjm))
439  ALLOCATE(corners_lon(nbseg,iim,jjm), corners_lat(nbseg,iim,jjm))
440  ALLOCATE(kindex(kjpindex))
441  !
442  lon=lon_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
443  lat=lat_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
444  DO in=1,nbseg
445     corners_lon(in,:,:)=corners_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank),in,1)
446     corners_lat(in,:,:)=corners_glo(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank),in,2)
447  ENDDO
448  !
449  !
450  ! Redistribute the indeces on all procs (apple distribution of land points)
451  !
452  CALL bcast(lon_g)
453  CALL bcast(lat_g)
454  CALL scatter(index_g, kindex)
455  !
456  !
457  ! Apply the offset needed so that kindex refers to the index of the land point
458  ! on the current region, i.e. the local lon lat domain.
459  !
460  kindex(1:kjpindex)=kindex(1:kjpindex)-(jj_begin-1)*iim_glo
461  !
462  ! This routine transforms the global grid into a series of polygons for all land
463  ! points identified by index_g.
464  !
465  CALL grid_stuff(nbindex_g, iim_g, jjm_g, lon_g, lat_g, index_g, contfrac_glo)
466  !
467  ! Distribute the global lalo to the local processor level lalo
468  !
469  ALLOCATE(lalo_loc(kjpindex,2))
470  CALL scatter(lalo_glo, lalo_loc)
471  lalo(:,:) = lalo_loc(:,:)
472  !
473  !====================================================================================
474  !-
475  !- Prepare the time for the simulation
476  !-
477  !- Set the calendar and get some information
478  !-
479  CALL ioconf_calendar(calendar)
480  CALL ioget_calendar(one_year, one_day)
481  !-
482  !- get the time period for the run
483  !-
484  CALL forcing_integration_time(date0, dt, nbdt)
485  write(numout, *) "orchideedriver date0", date0, dt, nbdt
486  !
487  !-
488  !- Set the start date in IOIPSL for the calendar and initialize the module time
489  !-
490  CALL ioconf_startdate(date0)
491  CALL time_initialize(0, date0, dt, "END")
492  !
493  !
494  !====================================================================================
495  !-
496  !- Initialize the forcing files and prepare the time stepping through the data.
497  !-
498  !
499  CALL forcing_open(forfilename, iim_glo,  jjm_glo, lon_glo, lat_glo, nbindex_g, zoom_lon, zoom_lat, &
500       &            index_g, kjpindex, numout)
501  !
502  !-
503  !====================================================================================
504  !-
505  !- Initialise the ORCHIDEE system in 4 steps :
506  !- 1 The control flags,
507  !- 2 Allocate memory (needs to be done after initializing the control flags because of nvm).
508  !- 2 the restart system of IOIPSL
509  !- 3 The history mechanism
510  !- 4 Finally the first call to SECHIBA will initialise all the internal variables
511  !
512  ! 1 Setting flags and some parameters (nvm)
513  !
514  CALL control_initialize
515  !
516  ! 2 Allocation
517  !
518  ALLOCATE(zlev_tq(kjpindex), zlev_uv(kjpindex))
519  ALLOCATE(u(kjpindex), v(kjpindex), pb(kjpindex))
520  ALLOCATE(temp_air(kjpindex))
521  ALLOCATE(qair(kjpindex))
522  ALLOCATE(petAcoef(kjpindex), peqAcoef(kjpindex), petBcoef(kjpindex), peqBcoef(kjpindex))
523  ALLOCATE(ccanopy(kjpindex))
524  ALLOCATE(cdrag(kjpindex))
525  ALLOCATE(precip_rain(kjpindex))
526  ALLOCATE(precip_snow(kjpindex))
527  ALLOCATE(swdown(kjpindex))
528  ALLOCATE(swnet(kjpindex))
529  ALLOCATE(lwdown(kjpindex))
530  ALLOCATE(sinang(kjpindex))
531  ALLOCATE(vevapp(kjpindex))
532  ALLOCATE(fluxsens(kjpindex))
533  ALLOCATE(fluxlat(kjpindex))
534  ALLOCATE(coastalflow(kjpindex))
535  ALLOCATE(riverflow(kjpindex))
536  ALLOCATE(netco2(kjpindex))
537  ALLOCATE(carblu(kjpindex))
538  ALLOCATE(carbwh(kjpindex))
539  ALLOCATE(carbha(kjpindex))
540  ALLOCATE(tsol_rad(kjpindex))
541  ALLOCATE(temp_sol_new(kjpindex))
542  ALLOCATE(qsurf(kjpindex))
543  ALLOCATE(albedo(kjpindex,2))
544  ALLOCATE(emis(kjpindex))
545  ALLOCATE(epot_air(kjpindex))
546  ALLOCATE(u_tq(kjpindex), v_tq(kjpindex))
547  ALLOCATE(z0m(kjpindex))
548  ALLOCATE(z0h(kjpindex))
549  ALLOCATE(veget_diag(kjpindex,nvm))
550  ALLOCATE(lai_diag(kjpindex,nvm))
551  ALLOCATE(height_diag(kjpindex,nvm))
552  !-
553  !---------------------------------------------------------------------------------------
554  !-
555  !- Get a first set of forcing data
556  !-
557  !---------------------------------------------------------------------------------------
558  !-
559  !- Some default values so that the operations before the ORCHIDEE initialisation do not fail.
560  !-
561  z0m(:) = 0.1
562  albedo(:,:) = 0.13
563  !
564  itau = 0
565  !
566  CALL ioipslctrl_restini(itau, date0, dt, rest_id, rest_id_stom, itau_offset, date0_shifted)
567  !
568  ! To ensure that itau starts with 0 at date0 for the restart, we have to set an off-set to achieve this.
569  ! itau_offset will get used to prduce itau_sechiba.
570  !
571  itau_offset=-itau_offset-1
572  !
573  ! Get the date of the first time step
574  !
575  WRITE(*,*) "itau_offset : date0 : ", year_end, month_end, day_end, sec_end
576 
577  !!- Initialize module for output with XIOS
578  CALL xios_orchidee_init( MPI_COMM_ORCH,                &
579       date0,    year_end, month_end, day_end, julian_diff,  &
580       lon,      lat,      znt)
581  CALL sechiba_xios_initialize 
582
583  CALL xios_orchidee_close_definition
584  IF (printlev_loc >= 2) WRITE(numout,*) 'After xios_orchidee_close_definition'
585
586  !- Initialize IOIPSL sechiba output files
587  itau_sechiba = itau+itau_offset
588  CALL ioipslctrl_history(iim, jjm, lon, lat,  kindex, kjpindex, itau_sechiba, &
589       date0, dt, hist_id, hist2_id, hist_id_stom, hist_id_stom_IPCC)
590  WRITE(*,*) "HISTORY : Defined for ", itau_sechiba, date0, dt
591  !
592  !-
593  !---------------------------------------------------------------------------------------
594  !-
595  !- Go into the time loop
596  !-
597  !---------------------------------------------------------------------------------------
598  !-
599  DO itau = 1,nbdt
600     !
601     CALL time_nextstep( itau )
602     !     
603     timestep_interval(1) = julian_start
604     timestep_interval(2) = julian_end
605     julian = (julian_start + julian_end) /2.0  !julian_end
606     
607     !
608     ! Get the forcing data
609     !
610     CALL forcing_getvalues(timestep_interval, dt, zlev_tq, zlev_uv, temp_air, qair, &
611          &                 precip_rain, precip_snow, swdown, lwdown, sinang, u, v, pb)
612     
613     !-
614     !
615     IF ( itau == nbdt ) lrestart_write = .TRUE.
616     !
617     ! Adaptation of the forcing data to SECHIBA's needs
618     !
619     ! Contrary to what the documentation says, ORCHIDEE expects surface pressure in hPa.
620     pb(:) = pb(:)/100.
621     epot_air(:) = cp_air*temp_air(:)+cte_grav*zlev_tq(:)
622     ccanopy(:) = atmco2
623     cdrag(:) = 0.0
624     !
625     petBcoef(:) = epot_air(:)
626     peqBcoef(:) = qair(:)
627     petAcoef(:) = zero
628     peqAcoef(:) = zero
629     !
630     ! Interpolate the wind (which is at hight zlev_uv) to the same height
631     ! as the temperature and humidity (at zlev_tq).
632     !
633     u_tq(:) = u(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:))
634     v_tq(:) = v(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:))
635     !
636     !
637     swnet(:) =(1.-(albedo(:,1)+albedo(:,2))/2.)*swdown(:)
638     !
639     !
640     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, temp_air, "RECEIVED Air temperature")
641     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, qair, "RECEIVED Air humidity")
642     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, precip_rain*one_day, "RECEIVED Rainfall")
643     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, precip_snow*one_day, "RECEIVED Snowfall")
644     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, swnet, "RECEIVED net solar")
645     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, lwdown, "RECEIVED lwdown")
646     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, u, "RECEIVED East-ward wind")
647     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, v, "RECEIVED North-ward wind")
648     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, pb*100, "RECEIVED surface pressure")
649     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, zlev_uv, "RECEIVED UV height")
650     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, zlev_tq, "RECEIVED TQ height")
651     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, sinang, "RECEIVED sinang")
652     !
653     IF ( itau .NE. 1 ) THEN
654        IF ( timemeasure ) THEN
655           waitget_cputime = waitget_cputime + Get_cpu_Time(timer_global)
656           waitget_walltime = waitget_walltime + Get_real_Time(timer_global)
657           CALL stop_timer(timer_global)
658           CALL start_timer(timer_global)
659        ENDIF
660     ENDIF
661     !
662     !---------------------------------------------------------------------------------------
663     !-
664     !- IF first time step : Call to SECHIBA_initialize to set-up ORCHIDEE before doing an actual call
665     !- which will provide the first fluxes.
666     !-
667     !---------------------------------------------------------------------------------------
668     !
669     itau_sechiba = itau+itau_offset
670     !
671     ! Update the calendar in xios by sending the new time step
672     !CALL xios_orchidee_update_calendar(itau_sechiba)
673     CALL xios_orchidee_update_calendar(itau_sechiba)
674     !
675     IF ( itau == 1 ) THEN
676        !
677        IF ( timemeasure ) THEN
678           WRITE(numout,*) '------> CPU Time for start-up of main : ',Get_cpu_Time(timer_global)
679           WRITE(numout,*) '------> Real Time for start-up of main : ',Get_real_Time(timer_global)
680           CALL stop_timer(timer_global)
681           CALL start_timer(timer_global)
682        ENDIF
683        !
684        CALL sechiba_initialize( &
685             itau_sechiba,  iim*jjm,      kjpindex,      kindex,                   &
686             lalo_loc,     contfrac,     neighbours,    resolution,  zlev_tq,      &
687             u_tq,         v_tq,         qair,          temp_air,                  &
688             petAcoef,     peqAcoef,     petBcoef,      peqBcoef,                  &
689             precip_rain,  precip_snow,  lwdown,        swnet,       swdown,       &
690             pb,           rest_id,      hist_id,       hist2_id,                  &
691             rest_id_stom, hist_id_stom, hist_id_stom_IPCC,                        &
692             coastalflow,  riverflow,    tsol_rad,      vevapp,      qsurf,        &
693             z0m,          z0h,          albedo,        fluxsens,    fluxlat,  emis, &
694             temp_sol_new, cdrag)
695        !
696        CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, temp_sol_new, "Init temp_sol_new")
697        !
698        ! Net solar and the wind at the right hight are recomputed with the correct values.
699        !
700        swnet(:) =(1.-(albedo(:,1)+albedo(:,2))/2.)*swdown(:)
701        u_tq(:) = u(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:))
702        v_tq(:) = v(:)*LOG(zlev_tq(:)/z0m(:))/LOG(zlev_uv(:)/z0m(:))
703        !
704        lrestart_read = .FALSE.
705        !
706        CALL histwrite_p(hist_id, 'LandPoints',  itau+1, (/ REAL(kindex) /), kjpindex, kindex)
707        CALL histwrite_p(hist_id, 'Areas',  itau+1, area, kjpindex, kindex)
708        CALL histwrite_p(hist_id, 'Contfrac',  itau+1, contfrac, kjpindex, kindex)
709        !
710        IF ( timemeasure ) THEN
711           WRITE(numout,*) '------> CPU Time for set-up of ORCHIDEE : ',Get_cpu_Time(timer_global)
712           WRITE(numout,*) '------> Real Time for set-up of ORCHIDEE : ',Get_real_Time(timer_global)
713           CALL stop_timer(timer_global)
714           CALL start_timer(timer_global)
715        ENDIF
716        !
717     ENDIF
718     !
719     !---------------------------------------------------------------------------------------
720     !-
721     !- Main call to SECHIBA 
722     !-
723     !---------------------------------------------------------------------------------------
724     !
725     !
726     !
727     CALL sechiba_main (itau_sechiba, iim*jjm, kjpindex, kindex, &
728          & lrestart_read, lrestart_write, &
729          & lalo_loc, contfrac, neighbours, resolution, &
730          ! First level conditions
731          & zlev_tq, u_tq, v_tq, qair, temp_air, epot_air, ccanopy, &
732          ! Variables for the implicit coupling
733          & cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
734          ! Rain, snow, radiation and surface pressure
735          & precip_rain ,precip_snow, lwdown, swnet, swdown, sinang, pb, &
736          ! Output : Fluxes
737          & vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
738          ! CO2 fluxes
739          & netco2, carblu, carbwh, carbha, &
740          ! Surface temperatures and surface properties
741          & tsol_rad, temp_sol_new, qsurf, albedo, emis, z0m, z0h, &
742          ! Vegetation, lai and vegetation height
743          & veget_diag, lai_diag, height_diag, &
744          ! File ids
745          & rest_id, hist_id, hist2_id, rest_id_stom, hist_id_stom, hist_id_stom_IPCC)
746     !
747     !
748     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, temp_sol_new, "Produced temp_sol_new")
749     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, fluxsens, "Produced fluxsens")
750     CALL forcing_printpoint(julian, testpt(1), testpt(2), kjpindex, lalo_loc, fluxlat, "Produced fluxlat")
751     !
752     IF ( timemeasure ) THEN
753        orchidee_cputime = orchidee_cputime + Get_cpu_Time(timer_global)
754        orchidee_walltime = orchidee_walltime + Get_real_Time(timer_global)
755        CALL stop_timer(timer_global)
756        CALL start_timer(timer_global)
757     ENDIF
758     !
759     !---------------------------------------------------------------------------------------
760     !-
761     !- Write diagnostics
762     !-
763     !---------------------------------------------------------------------------------------
764     !
765     CALL xios_orchidee_send_field("LandPoints" ,(/ ( REAL(ik), ik=1,kjpindex ) /))
766     CALL xios_orchidee_send_field("areas", area)
767     CALL xios_orchidee_send_field("contfrac",contfrac)
768     CALL xios_orchidee_send_field("temp_air",temp_air)
769     CALL xios_orchidee_send_field("qair",qair)
770     CALL xios_orchidee_send_field("swnet",swnet)
771     CALL xios_orchidee_send_field("swdown",swdown)
772     ! zpb in hPa, output in Pa
773     CALL xios_orchidee_send_field("pb",pb)
774     !
775     IF ( .NOT. almaoutput ) THEN
776        !
777        !  ORCHIDEE INPUT variables
778        !
779        CALL histwrite_p (hist_id, 'swdown',   itau_sechiba, swdown,   kjpindex, kindex)
780        CALL histwrite_p (hist_id, 'tair',     itau_sechiba, temp_air, kjpindex, kindex)
781        CALL histwrite_p (hist_id, 'qair',     itau_sechiba, qair, kjpindex, kindex)
782        CALL histwrite_p (hist_id, 'evap',     itau_sechiba, vevapp, kjpindex, kindex)
783        CALL histwrite_p (hist_id, 'coastalflow',itau_sechiba, coastalflow, kjpindex, kindex)
784        CALL histwrite_p (hist_id, 'riverflow',itau_sechiba, riverflow, kjpindex, kindex)
785        !
786        CALL histwrite_p (hist_id, 'temp_sol', itau_sechiba, temp_sol_new, kjpindex, kindex)
787        CALL histwrite_p (hist_id, 'tsol_max', itau_sechiba, temp_sol_new, kjpindex, kindex)
788        CALL histwrite_p (hist_id, 'tsol_min', itau_sechiba, temp_sol_new, kjpindex, kindex)
789        CALL histwrite_p (hist_id, 'fluxsens', itau_sechiba, fluxsens, kjpindex, kindex)
790        CALL histwrite_p (hist_id, 'fluxlat',  itau_sechiba, fluxlat,  kjpindex, kindex)
791        CALL histwrite_p (hist_id, 'swnet',    itau_sechiba, swnet,    kjpindex, kindex)
792        CALL histwrite_p (hist_id, 'alb_vis',  itau_sechiba, albedo(:,1), kjpindex, kindex)
793        CALL histwrite_p (hist_id, 'alb_nir',  itau_sechiba, albedo(:,2), kjpindex, kindex)
794        !
795        IF ( hist2_id > 0 ) THEN
796           CALL histwrite_p (hist2_id, 'swdown',   itau_sechiba, swdown, kjpindex, kindex)
797           CALL histwrite_p (hist2_id, 'tair',     itau_sechiba, temp_air, kjpindex, kindex)
798           CALL histwrite_p (hist2_id, 'qair',     itau_sechiba, qair, kjpindex, kindex)
799           !
800           CALL histwrite_p (hist2_id, 'evap',     itau_sechiba, vevapp, kjpindex, kindex)
801           CALL histwrite_p (hist2_id, 'coastalflow',itau_sechiba, coastalflow, kjpindex, kindex)
802           CALL histwrite_p (hist2_id, 'riverflow',itau_sechiba, riverflow, kjpindex, kindex)
803           !
804           CALL histwrite_p (hist2_id, 'temp_sol', itau_sechiba, temp_sol_new, kjpindex, kindex)
805           CALL histwrite_p (hist2_id, 'tsol_max', itau_sechiba, temp_sol_new, kjpindex, kindex)
806           CALL histwrite_p (hist2_id, 'tsol_min', itau_sechiba, temp_sol_new, kjpindex, kindex)
807           CALL histwrite_p (hist2_id, 'fluxsens', itau_sechiba, fluxsens, kjpindex, kindex)
808           CALL histwrite_p (hist2_id, 'fluxlat',  itau_sechiba, fluxlat,  kjpindex, kindex)
809           CALL histwrite_p (hist2_id, 'swnet',    itau_sechiba, swnet,    kjpindex, kindex)
810           !
811           CALL histwrite_p (hist2_id, 'alb_vis',  itau_sechiba, albedo(:,1), kjpindex, kindex)
812           CALL histwrite_p (hist2_id, 'alb_nir',  itau_sechiba, albedo(:,2), kjpindex, kindex)
813        ENDIF
814     ELSE
815        !
816        ! Input variables
817        !
818        CALL histwrite_p (hist_id, 'SinAng', itau_sechiba, sinang, kjpindex, kindex)
819        CALL histwrite_p (hist_id, 'LWdown', itau_sechiba, lwdown, kjpindex, kindex)
820        CALL histwrite_p (hist_id, 'SWdown', itau_sechiba, swdown, kjpindex, kindex)
821        CALL histwrite_p (hist_id, 'Tair', itau_sechiba, temp_air, kjpindex, kindex)
822        CALL histwrite_p (hist_id, 'Qair', itau_sechiba, qair, kjpindex, kindex)
823        CALL histwrite_p (hist_id, 'SurfP', itau_sechiba, pb, kjpindex, kindex)
824        CALL histwrite_p (hist_id, 'Windu', itau_sechiba, u_tq, kjpindex, kindex)
825        CALL histwrite_p (hist_id, 'Windv', itau_sechiba, v_tq, kjpindex, kindex)
826        !
827        CALL histwrite_p (hist_id, 'Evap', itau_sechiba, vevapp, kjpindex, kindex)
828        CALL histwrite_p (hist_id, 'SWnet',    itau_sechiba, swnet, kjpindex, kindex)
829        CALL histwrite_p (hist_id, 'Qh', itau_sechiba, fluxsens, kjpindex, kindex)
830        CALL histwrite_p (hist_id, 'Qle',  itau_sechiba, fluxlat, kjpindex, kindex)
831        CALL histwrite_p (hist_id, 'AvgSurfT', itau_sechiba, temp_sol_new, kjpindex, kindex)
832        CALL histwrite_p (hist_id, 'RadT', itau_sechiba, temp_sol_new, kjpindex, kindex)
833        !
834        ! There is a mess with the units passed to the coupler. To be checked with Marc
835        !
836        IF ( river_routing ) THEN
837           CALL histwrite_p (hist_id, 'CoastalFlow',itau_sechiba, coastalflow, kjpindex, kindex)
838           CALL histwrite_p (hist_id, 'RiverFlow',itau_sechiba, riverflow, kjpindex, kindex)
839        ENDIF
840        !
841        IF ( hist2_id > 0 ) THEN
842           CALL histwrite_p (hist2_id, 'Evap', itau_sechiba, vevapp, kjpindex, kindex)
843           CALL histwrite_p (hist2_id, 'SWnet',    itau_sechiba, swnet, kjpindex, kindex)
844           CALL histwrite_p (hist2_id, 'Qh', itau_sechiba, fluxsens, kjpindex, kindex)
845           CALL histwrite_p (hist2_id, 'Qle',  itau_sechiba, fluxlat, kjpindex, kindex)
846           CALL histwrite_p (hist2_id, 'AvgSurfT', itau_sechiba, temp_sol_new, kjpindex, kindex)
847           CALL histwrite_p (hist2_id, 'RadT', itau_sechiba, temp_sol_new, kjpindex, kindex)
848        ENDIF
849     ENDIF
850     !
851     !
852  ENDDO
853  !-
854  !-
855  !---------------------------------------------------------------------------------------
856  !-
857  !- Close eveything
858  !-
859  !--
860  ! Close IOIPSL history files
861  CALL histclo
862  IF(is_root_prc) THEN
863     ! Close restart files
864     CALL restclo
865     CALL getin_dump
866  ENDIF
867  !-
868  !- Deallocate all variables and reset initialization variables
869  !-
870  CALL orchideedriver_clear()
871  !
872  WRITE(numout,*) "End at proc ", mpi_rank
873  !
874  !
875  !---------------------------------------------------------------------------------------
876  !-
877  !- Get time and close IOIPSL, OASIS and MPI
878  !-
879  !---------------------------------------------------------------------------------------
880  !-
881  IF ( timemeasure ) THEN
882     WRITE(numout,*) '------> Total CPU Time waiting to get forcing : ',waitget_cputime
883     WRITE(numout,*) '------> Total Real Time waiting to get forcing : ',waitget_walltime
884     WRITE(numout,*) '------> Total CPU Time for ORCHIDEE : ', orchidee_cputime
885     WRITE(numout,*) '------> Total Real Time for ORCHIDEE : ', orchidee_walltime
886     WRITE(numout,*) '------> Total CPU Time waiting to put fluxes : ',waitput_cputime
887     WRITE(numout,*) '------> Total Real Time waiting to put fluxes : ',waitput_walltime
888     WRITE(numout,*) '------> Total CPU Time for closing : ',  Get_cpu_Time(timer_global)
889     WRITE(numout,*) '------> Total Real Time for closing : ', Get_real_Time(timer_global)
890     WRITE(numout,*) '------> Total without MPI : CPU Time  : ', Get_cpu_Time(timer_mpi)
891     WRITE(numout,*) '------> Total without MPI : Real Time : ', Get_real_Time(timer_mpi)
892     CALL stop_timer(timer_global)
893     CALL stop_timer(timer_mpi)
894  ENDIF
895  !
896  ! Finalize MPI and XIOS
897  CALL Finalize_mpi
898  !
899CONTAINS
900!
901!! ================================================================================================================================
902!! SUBROUTINE   : orchideedriver_clear
903!!
904!>\BRIEF         Clear orchideedriver
905!!
906!! DESCRIPTION  :  Deallocate memory and reset initialization variables to there original values
907!!                 This subroutine calls forcing_tools_clear and sechiba_clear.
908!!
909!_ ================================================================================================================================
910!
911  SUBROUTINE orchideedriver_clear
912    !- Deallocate all variables existing on all procs
913    !-
914    !- Deallocate all variables existing on all procs (list still incomplete)
915    !-
916    IF ( ALLOCATED(lon_glo) ) DEALLOCATE(lon_glo)
917    IF ( ALLOCATED(lat_glo) ) DEALLOCATE(lat_glo)
918    IF ( ALLOCATED(mask_glo) ) DEALLOCATE(mask_glo)
919    IF ( ALLOCATED(area_glo) ) DEALLOCATE(area_glo)
920    IF ( ALLOCATED(corners_glo) ) DEALLOCATE(corners_glo)
921    IF ( ALLOCATED(corners_lon) ) DEALLOCATE(corners_lon)
922    IF ( ALLOCATED(corners_lat) ) DEALLOCATE(corners_lat)
923    IF ( ALLOCATED(kindex_g) ) DEALLOCATE(kindex_g)
924    IF ( ALLOCATED(contfrac_glo) ) DEALLOCATE(contfrac_glo)
925    IF ( ALLOCATED(lalo_glo) ) DEALLOCATE(lalo_glo)
926    IF ( ALLOCATED(lon) ) DEALLOCATE(lon)
927    IF ( ALLOCATED(lat) ) DEALLOCATE(lat)
928    IF ( ALLOCATED(kindex) ) DEALLOCATE(kindex)
929    IF ( ALLOCATED(lalo_loc) ) DEALLOCATE(lalo_loc)
930    IF ( ALLOCATED(zlev_tq) ) DEALLOCATE(zlev_tq)
931    IF ( ALLOCATED(zlev_uv) ) DEALLOCATE(zlev_uv)
932    IF ( ALLOCATED(u) ) DEALLOCATE(u)
933    IF ( ALLOCATED(v) ) DEALLOCATE(v)
934    IF ( ALLOCATED(pb) ) DEALLOCATE(pb)
935    IF ( ALLOCATED(temp_air) ) DEALLOCATE(temp_air)
936    IF ( ALLOCATED(qair) ) DEALLOCATE(qair)
937    IF ( ALLOCATED(precip_rain) ) DEALLOCATE(precip_rain)
938    IF ( ALLOCATED(precip_snow) ) DEALLOCATE(precip_snow)
939    IF ( ALLOCATED(swdown) ) DEALLOCATE(swdown)
940    IF ( ALLOCATED(swnet) ) DEALLOCATE(swnet)
941    IF ( ALLOCATED(lwdown) ) DEALLOCATE(lwdown)
942    IF ( ALLOCATED(sinang) ) DEALLOCATE(sinang)
943    IF ( ALLOCATED(epot_air) ) DEALLOCATE(epot_air)
944    IF ( ALLOCATED(ccanopy) ) DEALLOCATE(ccanopy)
945    IF ( ALLOCATED(cdrag) ) DEALLOCATE(cdrag)
946    IF ( ALLOCATED(swnet) ) DEALLOCATE(swnet)
947    IF ( ALLOCATED(petAcoef) ) DEALLOCATE(petAcoef)
948    IF ( ALLOCATED(peqAcoef) ) DEALLOCATE(peqAcoef)
949    IF ( ALLOCATED(petBcoef) ) DEALLOCATE(petBcoef)
950    IF ( ALLOCATED(peqBcoef) ) DEALLOCATE(peqBcoef)
951    IF ( ALLOCATED(u_tq) ) DEALLOCATE(u_tq)
952    IF ( ALLOCATED(v_tq) ) DEALLOCATE(v_tq)
953    IF ( ALLOCATED(vevapp) ) DEALLOCATE(vevapp)
954    IF ( ALLOCATED(fluxsens) ) DEALLOCATE(fluxsens)
955    IF ( ALLOCATED(fluxlat) ) DEALLOCATE(fluxlat)
956    IF ( ALLOCATED(coastalflow) ) DEALLOCATE(coastalflow)
957    IF ( ALLOCATED(riverflow) ) DEALLOCATE(riverflow)
958    IF ( ALLOCATED(netco2) ) DEALLOCATE(netco2)
959    IF ( ALLOCATED(carblu) ) DEALLOCATE(carblu)
960    IF ( ALLOCATED(carbwh) ) DEALLOCATE(carbwh)
961    IF ( ALLOCATED(carbha) ) DEALLOCATE(carbha)
962    IF ( ALLOCATED(tsol_rad) ) DEALLOCATE(tsol_rad)
963    IF ( ALLOCATED(temp_sol_new) ) DEALLOCATE(temp_sol_new)
964    IF ( ALLOCATED(qsurf) ) DEALLOCATE(qsurf)
965    IF ( ALLOCATED(albedo) ) DEALLOCATE(albedo)
966    IF ( ALLOCATED(emis) ) DEALLOCATE(emis)
967    IF ( ALLOCATED(z0m) ) DEALLOCATE(z0m)
968    IF ( ALLOCATED(z0h) ) DEALLOCATE(z0h)
969    !
970    WRITE(numout,*) "Memory deallocated"
971    !
972  END SUBROUTINE orchideedriver_clear
973  !
974END PROGRAM orchidedriver
Note: See TracBrowser for help on using the repository browser.