Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_driver/orchideedriver.f90 @ 6715

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

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

Download in other formats: