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

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

As in r6651 and r6652 for new driver.

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