Context Navigation

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, 3 years 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	!
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
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_airtemp_air(:)+cte_gravzlev_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	!
893	CONTAINS
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	!
968	END PROGRAM orchidedriver

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

Download in other formats: