Context Navigation

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

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

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

Download in other formats: