Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_driver/dim2_driver.f90 @ 7852

Last change on this file since 7852 was 7460, checked in by josefine.ghattas, 2 years ago
Uniformized all parameters SPRED_PREC and SPREAD_PREC. Now the A is always included. This was already done in ORCHIDEE_2_2 rev [7264], see ticket #821
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 58.9 KB

Line
1	! ==============================================================================================================================\n
2	! PROGRAM : driver
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Main program for the offline driver dim2_driver
10	!!
11	!!\n DESCRIPTION : This PROGRAM is the driver of the dim2 version of ORCHIDEE.
12	!!
13	!! The forcing data needs to be in netCDF format and should contain the following variables :
14	!! - Incoming SW radiation
15	!! - Incoming LW radiation
16	!! - Precipitation
17	!! - Air temperature at a reference level
18	!! - Air humidity at the same level
19	!! - wind at the same level
20	!! - surface pressure
21	!!
22	!! RECENT CHANGE(S): None
23	!!
24	!! REFERENCE(S) : None
25	!!
26	!! SVN :
27	!! $HeadURL$
28	!! $Date$
29	!! $Revision$
30	!! \n
31	!_ ================================================================================================================================
32
33	PROGRAM driver
34
35	USE netcdf
36	USE ioipsl_para
37	USE grid
38	USE intersurf, ONLY : intersurf_main_2d, intersurf_initialize_2d, intersurf_clear
39	USE constantes
40	USE time
41	USE readdim2
42	USE mod_orchidee_para
43	USE timer
44
45	IMPLICIT NONE
46
47	INTEGER :: iim, jjm, llm
48	INTEGER :: im, jm, lm, tm, is, force_id, itest, jtest
49	REAL :: dt, dt_force, dt_rest, date0, date0_rest
50	REAL :: date_cur ! Starting date of the current execution
51	REAL :: zlflu
52	REAL :: alpha
53
54	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
55	& swdown, coszang, precip_rain, precip_snow, tair_obs, u, v, &
56	& qair_obs, pb, lwdown, &
57	& eair_obs, zlev_vec, zlevuv_vec, relax
58	!- Variables which are forcings for SECHIBA
59	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
60	& petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, &
61	& for_u, for_v, for_swnet, for_swdown, for_coszang, for_lwdown, &
62	& for_psurf, for_qair, for_tair, for_eair, &
63	& for_ccanopy, for_rau
64
65	REAL, ALLOCATABLE, DIMENSION(:,:) :: &
66	& for_contfrac, old_zlev, old_qair, old_eair, tsol_rad, vevapp, &
67	& temp_sol_NEW, temp_sol_old, qsurf, dtdt, coastalflow, riverflow, &
68	& fluxsens, fluxlat, emis, z0, tmp_z0
69
70	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: for_neighbours
71
72	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: for_resolution
73
74	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: albedo
75	REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_vis
76	REAL, ALLOCATABLE, DIMENSION(:,:) :: albedo_nir
77
78	INTEGER, ALLOCATABLE, DIMENSION(:) :: kindex
79	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
80	REAL, ALLOCATABLE, DIMENSION(:) :: tmplev
81
82	REAL :: old_tair
83	REAL :: atmco2,co2inc
84	LOGICAL :: CO2_varying
85	INTEGER :: nbindex
86	REAL :: julian, ss
87	INTEGER :: yy, mm, dd, yy_prev
88
89	LOGICAL :: relaxation
90
91	CHARACTER(LEN=512) :: filename, driv_restname_in, driv_restname_out
92	CHARACTER(LEN=30) :: time_str, var_name
93
94	CHARACTER(LEN=100) :: str, str1 ! temporary variables to print data
95
96	INTEGER :: it, istp, istp_old, rest_id, it_force
97
98	INTEGER :: split, split_start, nb_spread, for_offset
99	INTEGER :: itau_dep, itau_dep_rest, itau_fin, itau_skip, itau_len
100
101	INTEGER,DIMENSION(2) :: ml
102
103	LOGICAL :: lstep_init, lstep_last
104	LOGICAL :: lwdown_cons !! Flag to conserve lwdown radiation from forcing
105	LOGICAL :: swdown_cons !! Flag to conserve swdown radiation from forcing
106
107	! to check variables passed to intersurf
108	INTEGER :: ik
109	INTEGER :: i,j
110	INTEGER :: printlev_loc !! local write level
111
112	LOGICAL :: Flag
113	LOGICAL :: driver_reset_time
114
115	REAL :: fill_init
116
117	fill_init=REAL(nf90_fill_real,r_std)
118	CALL ioconf_expval(val_exp)
119	!-
120	! Init parallelism
121
122	CALL Init_orchidee_para
123	CALL init_timer
124	CALL start_timer(timer_global)
125	CALL start_timer(timer_mpi)
126
127	! Set specific write level to dim2_driver using PRINTLEV_dim2_driver=[0-4]
128	! in run.def. The global printlev is used as default value.
129	printlev_loc=get_printlev('dim2_driver')
130
131
132	!=====================================================================
133	!- 1.0 This section defines the general set-up of the experiment :
134	!- - Forcing data to be used with its time step
135	!- - Restart file to be used
136	!- - The time step that will be used
137	!- - The starting date of the simulation
138	!- - Length of integration
139	!- - Basic flags for SSIPSL
140
141	!Config Key = LWDOWN_CONS
142	!Config Desc = Conserve longwave downwelling radiation in the forcing
143	!Config Def = n
144	!Config If =
145	!Config Help = If LWDOWN_CONS=False a non conservative interpolation is used
146	!Config Units = [FLAG]
147	lwdown_cons = .FALSE.
148	CALL getin_p('LWDOWN_CONS', lwdown_cons)
149
150	!Config Key = SWDOWN_CONS
151	!Config Desc = Conserve shortwave downwelling radiation in the forcing
152	!Config Def = LWDOWN_CONS
153	!Config If =
154	!Config Help = If SWDOWN_CONS=False a non conservative interpolation is used
155	!Config Units = [FLAG]
156	swdown_cons = lwdown_cons
157	CALL getin_p('SWDOWN_CONS', swdown_cons)
158
159	!=====================================================================
160	!- 1.1 Initialize the driving variables. It essentialy gets the mode
161	!- used and the size of the driving variables.
162	!=====================================================================
163	IF (printlev_loc>=4) WRITE(numout,*) 'Reading name of the forcing file'
164	!-
165	!Config Key = FORCING_FILE
166	!Config Desc = Name of file containing the forcing data
167	!Config If = [-]
168	!Config Def = forcing_file.nc
169	!Config Help = This is the name of the file which should be opened
170	!Config for reading the forcing data of the dim0 model.
171	!Config The format of the file has to be netCDF and COADS
172	!Config compliant.
173	!Config Units = [FILE]
174	!-
175	filename='forcing_file.nc'
176	CALL getin_p('FORCING_FILE',filename)
177	!-
178	IF (printlev_loc>=4) WRITE(numout,*) 'Opening forcing file'
179	!-
180	! We call flininfo to obtain the dimensions
181	! of iim, jjm and others.
182	! This information will allow us to allocate all the space needed.
183	!-
184	CALL forcing_info &
185	& (filename, iim, jjm, llm, tm, date0, dt_force, force_id)
186	!-
187	IF (printlev>=1) WRITE(numout,*) 'Information about forcing file : date0 ', date0, &
188	'iim, jjm, llm, tm',iim,jjm,llm,tm,' dt_force ',dt_force
189	!-
190	CALL init_ioipsl_para
191	!-
192	IF (printlev_loc>=4) THEN
193	WRITE(numout,*) 'Allocate memory for the driver :', iim, jjm, llm
194	ENDIF
195	!-
196	ALLOCATE (tmplev(llm))
197	tmplev(:)=0.
198
199	ALLOCATE &
200	& (swdown(iim,jjm), coszang(iim,jjm), precip_rain(iim,jjm), precip_snow(iim,jjm), &
201	& tair_obs(iim,jjm), u(iim,jjm), v(iim,jjm), qair_obs(iim,jjm), &
202	& pb(iim,jjm), lwdown(iim,jjm), &
203	& eair_obs(iim,jjm), zlev_vec(iim,jjm), zlevuv_vec(iim,jjm), relax(iim,jjm))
204	!-
205	ALLOCATE &
206	& (petAcoef(iim,jjm), peqAcoef(iim,jjm), &
207	& petBcoef(iim,jjm), peqBcoef(iim,jjm), &
208	& cdrag(iim,jjm), for_u(iim,jjm), for_v(iim,jjm), &
209	& for_swnet(iim,jjm), for_swdown(iim,jjm), for_coszang(iim,jjm), for_lwdown(iim,jjm), &
210	& for_psurf(iim,jjm), for_qair(iim,jjm), for_tair(iim,jjm), &
211	& for_eair(iim,jjm), for_ccanopy(iim,jjm), for_rau(iim,jjm))
212
213	ALLOCATE &
214	& (for_contfrac(iim,jjm), for_neighbours(iim,jjm,8), for_resolution(iim,jjm,2), &
215	& old_zlev(iim,jjm), old_qair(iim,jjm), old_eair(iim,jjm), &
216	& tsol_rad(iim,jjm), vevapp(iim,jjm), &
217	& temp_sol_NEW(iim,jjm), temp_sol_old(iim,jjm), &
218	& dtdt(iim,jjm), coastalflow(iim,jjm), riverflow(iim,jjm), &
219	& fluxsens(iim,jjm), fluxlat(iim,jjm), emis(iim,jjm), &
220	& z0(iim,jjm), tmp_z0(iim,jjm), qsurf(iim,jjm))
221
222	ALLOCATE(albedo(iim,jjm,2))
223	ALLOCATE(albedo_vis(iim,jjm),albedo_nir(iim,jjm))
224	ALLOCATE(kindex(iim*jjm))
225	ALLOCATE(lon(iim,jjm), lat(iim,jjm))
226
227	swdown(:,:) = fill_init
228	precip_rain(:,:) = 0.0
229	precip_snow(:,:) = 0.0
230	tair_obs(:,:) = 0.0
231	u(:,:) = fill_init
232	v(:,:) = fill_init
233	qair_obs(:,:) = fill_init
234	pb(:,:) = fill_init
235	lwdown(:,:) = fill_init
236	eair_obs(:,:) = fill_init
237	zlev_vec(:,:) = 0.0
238	zlevuv_vec(:,:) = 0.0
239	relax(:,:) = 0.0
240	petAcoef(:,:) = 0.0
241	peqAcoef(:,:) = 0.0
242	petBcoef(:,:) = 0.0
243	peqBcoef(:,:) = 0.0
244	cdrag(:,:) = 0.0
245	for_u(:,:) = fill_init
246	for_v(:,:) = fill_init
247	for_swnet(:,:) = fill_init
248	for_swdown(:,:) = fill_init
249	for_lwdown(:,:) = fill_init
250	for_psurf(:,:) = fill_init
251	for_qair(:,:) = fill_init
252	for_tair(:,:) = fill_init
253	for_eair(:,:) = fill_init
254	for_ccanopy(:,:) = 0.0
255	for_rau(:,:) = fill_init
256	for_contfrac(:,:) = 0.0
257	for_neighbours(:,:,:) = 0
258	for_resolution(:,:,:) = 0.0
259	old_zlev(:,:) = 0.0
260	old_qair(:,:) = 0.0
261	old_eair(:,:) = 0.0
262	tsol_rad(:,:) = 0.0
263	vevapp(:,:) = 0.0
264	temp_sol_NEW(:,:) = fill_init
265	temp_sol_old(:,:) = fill_init
266	dtdt(:,:) = 0.0
267	coastalflow(:,:) = 0.0
268	riverflow(:,:) = 0.0
269	fluxsens(:,:) = fill_init
270	fluxlat(:,:) = fill_init
271	emis(:,:) = 0.0
272	z0(:,:) = fill_init
273	tmp_z0(:,:) = fill_init
274	qsurf(:,:) = 0.0
275	albedo(:,:,:) = fill_init
276	albedo_vis(:,:) = fill_init
277	albedo_nir(:,:) = fill_init
278	kindex(:) = 0
279	lon(:,:) = 0.0
280	lat(:,:) = 0.0
281	!-
282	! We need to know the grid.
283	! Then we can initialize the restart files, and then we
284	! can give read the restart files in the forcing subroutine.
285	!-
286	CALL forcing_grid (iim,jjm,llm,lon,lat,init_f=.FALSE.)
287	!=====================================================================
288	!- 1.2 Time step to be used.
289	!- This is relative to the time step of the forcing data
290	!=====================================================================
291	IF ( .NOT. weathergen ) THEN
292	!Config Key = DT_SECHIBA
293	!Config Desc = Time-step of the SECHIBA component
294	!Config If = NOT(WEATHERGEN)
295	!Config Def = 1800.
296	!Config Help = Determines the time resolution at which
297	!Config the calculations in the SECHIBA component
298	!Config are done
299	!Config Units = [seconds]
300	dt = 1800
301	CALL getin_p('DT_SECHIBA', dt)
302	split = INT(dt_force/dt)
303	IF (printlev_loc >= 1) WRITE(numout,*) 'Time step in forcing file: dt_force=',dt_force
304	IF (printlev_loc >= 1) WRITE(numout,*) 'Time step in sechiba component: dt_sechiba=',dt
305	IF (printlev_loc >= 1) WRITE(numout,*) 'Splitting of each forcing time step: split=',split
306
307	IF ( split .LT. 1. ) THEN
308	CALL ipslerr_p ( 3, 'dim2_driver',&
309	'Time step of the forcing file is higher than the time step in sechiba',&
310	'Please, modify DT_SECHIBA parameter value !','')
311	END IF
312	ELSE
313	! Case weathergen:
314	! The model time step in sechiba is always the same as the forcing time step
315	dt = dt_force
316	split = 1
317	ENDIF
318
319	!=====================================================================
320	!- 1.3 Initialize the restart file for the driver
321	!=====================================================================
322	!Config Key = RESTART_FILEIN
323	!Config Desc = Name of restart to READ for initial conditions
324	!Config If = [-]
325	!Config Def = NONE
326	!Config Help = This is the name of the file which will be opened
327	!Config to extract the initial values of all prognostic
328	!Config values of the model. This has to be a netCDF file.
329	!Config Not truly COADS compliant. NONE will mean that
330	!Config no restart file is to be expected.
331	!Config Units = [FILE]
332	!-
333	driv_restname_in = 'NONE'
334	CALL getin_p('RESTART_FILEIN',driv_restname_in)
335	if (printlev_loc>=4) WRITE(numout,*) 'INPUT RESTART_FILE : ',TRIM(driv_restname_in)
336	!-
337	!Config Key = RESTART_FILEOUT
338	!Config Desc = Name of restart files to be created by the driver
339	!Config If = [-]
340	!Config Def = driver_rest_out.nc
341	!Config Help = This variable give the name for
342	!Config the restart files. The restart software within
343	!Config IOIPSL will add .nc if needed
344	!Config Units = [FILE]
345	!-
346	driv_restname_out = 'driver_rest_out.nc'
347	CALL getin_p('RESTART_FILEOUT', driv_restname_out)
348	if (printlev_loc>=4) WRITE(numout,*) 'OUTPUT RESTART_FILE : ',TRIM(driv_restname_out)
349	!-
350	! Set default values for the start and end of the simulation
351	! in the forcing chronology.
352
353	CALL gather2D_mpi(lon,lon_g)
354	CALL gather2D_mpi(lat,lat_g)
355
356
357	!Config Key = DRIVER_reset_time
358	!Config Desc = Overwrite time values from the driver restart file
359	!Config If = [-]
360	!Config Def = n
361	!Config Units = [FLAG]
362
363	driver_reset_time=.FALSE.
364	CALL getin_p('DRIVER_reset_time', driver_reset_time)
365	IF (printlev_loc>=4) WRITE(numout,*) 'driver_reset_time=',driver_reset_time
366
367	IF (is_root_prc) THEN
368	! Set default values for the time variables
369	itau_dep_rest = 0
370	date0_rest = date0
371	dt_rest = dt
372
373	IF (printlev_loc>=4) WRITE(numout,*) &
374	'Before driver restart file initialization : itau_dep_rest, date0_rest, dt_rest = ', &
375	itau_dep_rest, date0_rest, dt_rest
376
377	CALL restini &
378	(driv_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, tmplev, &
379	driv_restname_out, itau_dep_rest, date0_rest, dt_rest, rest_id, driver_reset_time, &
380	use_compression=nc_restart_compression)
381
382	IF (printlev_loc>=4) WRITE(numout,*) &
383	'After driver restart file initialization : itau_dep_rest, date0_rest, dt_rest = ', &
384	itau_dep_rest, date0_rest, dt_rest
385
386	IF (itau_dep_rest == 0 .OR. driver_reset_time) THEN
387	! Restart file did not exist or we decide to overwrite time values in it.
388	! Set time values to read the begining of the forcing file.
389	itau_dep=0
390	itau_fin=tm
391	date0_rest = date0
392	ELSE
393	! Take time values from restart file
394	itau_dep = itau_dep_rest
395	itau_fin = itau_dep+tm
396	ENDIF
397
398	IF (printlev_loc >= 1) WRITE(numout,*) &
399	'Restart file initialized : itau_dep, itau_fin, date0_rest, dt_rest = ', &
400	itau_dep, itau_fin, date0_rest, dt_rest
401	ENDIF
402
403	! Communicate values from root_prc to all the others
404	CALL bcast (itau_dep_rest)
405	CALL bcast (itau_dep)
406	CALL bcast (itau_fin)
407	CALL bcast (date0_rest)
408	CALL bcast (dt_rest)
409
410	!=====================================================================
411	!- 1.4 Here we do the first real reading of the driving. It only
412	!- gets a few variables.
413	!=====================================================================
414
415	! prepares kindex table from the information obtained
416	! from the forcing data and reads some initial values for
417	! temperature, etc.
418	!-
419	kindex(1) = 1
420	!-
421	CALL forcing_READ &
422	& (filename, rest_id, .TRUE., .FALSE., &
423	& 0, itau_dep, 0, split, nb_spread, lwdown_cons, swdown_cons, &
424	& date0, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, &
425	& swdown, coszang, precip_rain, precip_snow, tair_obs, &
426	& u, v, qair_obs, pb, lwdown, for_contfrac, for_neighbours, for_resolution, &
427	& for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, &
428	& kindex, nbindex, force_id)
429	!-
430	IF (printlev_loc >= 2) WRITE (numout,*) ">> Number of land points =",nbindex
431	IF (nbindex == 0) THEN
432	WRITE(numout,*) "Limits : (W,E / N,S)", limit_west, limit_east, &
433	& limit_north, limit_south
434	CALL ipslerr_p ( 3, 'dim2_driver','number of land points error.', &
435	& ' is zero !','stop driver')
436	ENDIF
437	!-
438	DO ik=1,nbindex
439	jlandindex(ik) = (((kindex(ik)-1)/iim) + 1)
440	ilandindex(ik) = (kindex(ik) - (jlandindex(ik)-1)*iim)
441	ENDDO
442	IF (printlev_loc>=4) THEN
443	WRITE(numout,*) "kindex of land points : ", kindex(1:nbindex)
444	WRITE(numout,*) "index i of land points : ", ilandindex
445	WRITE(numout,*) "index j of land points : ", jlandindex
446	ENDIF
447
448	im = iim; jm = jjm; lm = llm;
449	IF ( (iim > 1).AND.(jjm > 1) ) THEN
450	jtest = INT((kindex(INT(nbindex/2))-1)/iim)+1
451	itest = MAX( 1, kindex(INT(nbindex/2))-(jtest-1)*iim )
452	ELSE
453	jtest = 1
454	itest = 1
455	ENDIF
456	IF (printlev_loc>=3) WRITE(numout,*) "test point in dim2_driver : ",itest,jtest
457	!-
458	IF ((im /= iim) .AND. (jm /= jjm) .AND. (lm /= llm)) THEN
459	WRITE(numout,*) ' dimensions are not good. Verify FILE :'
460	WRITE(numout,*) ' filename = ',filename
461	WRITE(numout,*) ' im, jm, lm lus = ', im, jm, lm
462	WRITE(numout,*) ' iim, jjm, llm demandes = ', iim, jjm, llm
463	CALL ipslerr_p(3,'dim2_driver','Pb in dimensions','','')
464	ENDIF
465	!=====================================================================
466	!- 1.5 Configures the time-steps and other parameters
467	!- of the run.
468	!=====================================================================
469	!-
470	! If the time steping of the restart is different from the one
471	! of the forcing we need to convert the itau_dep into the
472	! chronology of the forcing. This ensures that the forcing will
473	! start at the date of the restart file. Obviously the itau_fin
474	! needs to be shifted as well !
475	!-
476	IF ( (dt_rest /= dt_force).AND.(itau_dep > 1) ) THEN
477	itau_dep = NINT((itau_dep*dt_rest )/dt_force)
478	itau_fin = itau_dep+tm
479	if (printlev_loc>=3) WRITE(numout,*) &
480	& 'The time steping of the restart is different from the one ',&
481	& 'of the forcing we need to convert the itau_dep into the ',&
482	& 'chronology of the forcing. This ensures that the forcing will ',&
483	& 'start at the date of the restart file. Obviously the itau_fin ',&
484	& 'needs to be shifted as well : itau_dep, itau_fin ', &
485	& itau_dep, itau_fin
486	ENDIF
487	!-
488	! Same things if the starting dates are not the same.
489	! Everything should look as if we had only one forcing file !
490	!-
491	IF (date0 /= date0_rest) THEN
492	WRITE(numout,*) 'date0_rest , date0 : ',date0_rest , date0
493	for_offset = NINT((date0_rest-date0)*one_day/dt_force)
494	ELSE IF (date0 .gt. date0_rest) THEN
495	write(str,*) "date0=", date0
496	write(str1,*) "date0_rest=", date0_rest
497	CALL ipslerr_p(3, 'dim2_driver', &
498	'restart date (date0_rest) must be bigger or equal than the forcing file date (date0)',&
499	trim(str),trim(str1))
500	ELSE
501	for_offset = 0
502	ENDIF
503	IF (printlev_loc >= 3) WRITE(numout,*) 'OFFSET FOR THE data read :', for_offset
504
505	CALL ioconf_startdate(date0_rest)
506	!-
507	!Config Key = TIME_SKIP
508	!Config Desc = Time in the forcing file at which the model is started.
509	!Config If = [-]
510	!Config Def = 0
511	!Config Help = This time give the point in time at which the model
512	!Config should be started. If exists, the date of the restart file is use.
513	!Config The FORMAT of this date can be either of the following :
514	!Config n : time step n within the forcing file
515	!Config nS : n seconds after the first time-step in the file
516	!Config nD : n days after the first time-step
517	!Config nM : n month after the first time-step (year of 365 days)
518	!Config nY : n years after the first time-step (year of 365 days)
519	!Config Or combinations :
520	!Config nYmM: n years and m month
521	!Config Units = [seconds, days, months, years]
522	!-
523	itau_skip = 0
524	WRITE(time_str,'(I10)') itau_skip
525	CALL getin_p('TIME_SKIP', time_str)
526	!-
527	! Transform into itau
528	!-
529	CALL tlen2itau (time_str, dt_force, date0, itau_skip)
530	!-
531	itau_dep = itau_dep+itau_skip
532	!-
533	! We need to select the right position of the splited time steps.
534	!-
535	istp = itau_dep*split+1
536	IF (MOD(istp-1,split) /= 0) THEN
537	split_start = MOD(istp-1,split)+1
538	ELSE
539	split_start = 1
540	ENDIF
541	istp_old = itau_dep_rest
542	itau_len = itau_fin-itau_dep
543
544	!Config Key = TIME_LENGTH
545	!Config Desc = Length of the integration in time.
546	!Config If = [-]
547	!Config Def = Full length of the forcing file
548	!Config Help = Length of integration. By default the entire length
549	!Config of the forcing is used. The FORMAT of this date can
550	!Config be either of the following :
551	!Config n : time step n within the forcing file
552	!Config nS : n seconds after the first time-step in the file
553	!Config nD : n days after the first time-step
554	!Config nM : n month after the first time-step (year of 365 days)
555	!Config nY : n years after the first time-step (year of 365 days)
556	!Config Or combinations :
557	!Config nYmM: n years and m month
558	!Config Units = [seconds, days, months, years]
559	!-
560	WRITE(time_str,'(I10)') itau_len
561	CALL getin_p('TIME_LENGTH', time_str)
562	!-
563	! Transform time lenght into number of time step (itau_len) using date of current year
564	!-
565	date_cur = itau2date(itau_dep, date0, dt_force)
566	CALL tlen2itau (time_str, dt_force, date_cur, itau_len)
567	!-
568	itau_fin = itau_dep+itau_len
569	!-
570	IF (printlev_loc >= 1) THEN
571	WRITE(numout,*) '>> Time origine in the forcing file :', date0
572	WRITE(numout,*) '>> Time origine in the restart file :', date0_rest
573	WRITE(numout,*) '>> Simulate starts at forcing time-step : ', itau_dep
574	WRITE(numout,*) '>> The splited time-steps start at (Sets the '
575	WRITE(numout,*) '>> chronology for the history and restart files):',istp
576	WRITE(numout,*) '>> The time spliting starts at :', split_start
577	WRITE(numout,*) '>> Simulation ends at forcing time-step: ', itau_fin
578	WRITE(numout,*) '>> Length of the simulation is thus :', itau_len
579	WRITE(numout,*) '>> Length of the forcing data is in time-steps : ', tm
580	WRITE(numout,*) '>> Time steps : true, forcing and restart : ', dt,dt_force,dt_rest
581	END IF
582
583	IF (tm < itau_len) THEN
584	CALL ipslerr_p ( 2, 'dim2_driver','Length of the simulation is greater than.', &
585	' Length of the forcing data is in time-steps','verify TIME_LENGTH parameter.')
586	ENDIF
587
588
589	!=====================================================================
590	!- 2.0 This section is going to define the details by which
591	!- the input data is going to be used to force the
592	!- land-surface scheme. The tasks are the following :
593	!- - Is the coupling going to be explicit or implicit
594	!- - Type of interpolation to be used.
595	!- - At which height are the atmospheric forcings going to be used ?
596	!- - Is a relaxation method going to be used on the forcing
597	!- - Does net radiation in the interpolated data need to be conserved
598	!- - How do we distribute precipitation.
599	!=====================================================================
600	!Config Key = RELAXATION
601	!Config Desc = method of forcing
602	!Config If = [-]
603	!Config Def = n
604	!Config Help = A method is proposed by which the first atmospheric
605	!Config level is not directly forced by observations but
606	!Config relaxed with a time constant towards observations.
607	!Config For the moment the methods tends to smooth too much
608	!Config the diurnal cycle and introduces a time shift.
609	!Config A more sophisticated method is needed.
610	!Config Units = [FLAG]
611	!-
612	relaxation = .FALSE.
613	CALL getin_p('RELAXATION', relaxation)
614	IF ( relaxation ) THEN
615	WRITE(numout,*) 'dim2_driver : The relaxation option is temporarily disabled as it does not'
616	WRITE(numout,*) ' produce energy conservation in ORCHIDEE. If you intend to use it'
617	WRITE(numout,*) ' you should validate it.'
618	CALL ipslerr_p(3,'dim2_driver','relaxation option is not activated.','This option needs to be validated.','')
619
620	!Config Key = RELAX_A
621	!Config Desc = Time constant of the relaxation layer
622	!Config If = RELAXATION
623	!Config Def = 1.0
624	!Config Help = The time constant associated to the atmospheric
625	!Config conditions which are going to be computed
626	!Config in the relaxed layer. To avoid too much
627	!Config damping the value should be larger than 1000.
628	!Config Units = [days?]
629	!-
630	alpha = 1000.0
631	CALL getin_p('RELAX_A', alpha)
632	ENDIF
633
634	!Config Key = SPREAD_PREC
635	!Config Desc = On how long we spread the precipitation, value in nb of dt_sechiba
636	!Config If = orchidee_ol
637	!Config Def = Half of the forcing time step or uniform, depending on dt_force and dt_sechiba
638	!Config Help = Spread the precipitation over SPREAD_PREC steps of the splited forcing
639	!Config time step. This ONLY applied if the forcing time step has been splited.
640	!Config If the value indicated is greater than SPLIT_DT, SPLIT_DT is used for it.
641	!Config Units = [-]
642	!-
643	IF ( dt_force >= 3*one_hour) THEN
644	! Distribut the precipitations over the half of the forcing time step
645	nb_spread = INT(0.5 * (dt_force/dt))
646	ELSE
647	! Distribut the precipitations uniformly over the forcing time step
648	nb_spread = dt_force/dt
649	END IF
650
651	CALL getin_p('SPREAD_PREC', nb_spread)
652	IF (nb_spread > split) THEN
653	WRITE(numout,*) 'WARNING : nb_spread is too large it will be '
654	WRITE(numout,*) ' set to the value of split'
655	nb_spread = split
656	ELSE IF (split == 1) THEN
657	nb_spread = 1
658	ENDIF
659
660
661	!=====================================================================
662	!- 3.0 Finaly we can prepare all the variables before launching
663	!- the simulation !
664	!=====================================================================
665	! Initialize LOGICAL and the length of the integration
666	!-
667	lstep_init = .TRUE.
668	lstep_last = .FALSE.
669
670	temp_sol_NEW(:,:) = tp_00
671	!-
672	!Config Key = ATM_CO2
673	!Config Desc = Value to precribe atmosoheric CO2
674	!Config If = [FORCE_CO2_VEG=y or Offline mode]
675	!Config Def = 350.
676	!Config Help = Used in offline mode or in coupled mode if FORCE_CO2_VEG=y
677	!Config Units = [ppm]
678	atmco2=350.
679	CALL getin_p('ATM_CO2',atmco2)
680	for_ccanopy(:,:)=atmco2
681
682	!Config Key = CO2_varying
683	!Config Desc = A flag to specify if CO2 level will vary within the simulation
684	!Config If = [FORCE_CO2_VEG=y or Offline mode]
685	!Config Def = .FALSE.
686	!Config Help =
687	!Config Units = [y/n]
688	CO2_varying=.FALSE.
689	CALL getin_p('CO2_varying',CO2_varying)
690
691
692	!Config Key = CO2_inc
693	!Config Desc = Relative yearly increase of the CO2 level
694	!Config If = [FORCE_CO2_VEG=y or Offline mode]
695	!Config Def = 1.
696	!Config Help =
697	!Config Units = [-]
698	co2inc=1.
699	IF (CO2_varying) THEN
700	CALL getin_p('CO2_inc',co2inc)
701	ENDIF
702
703	!-
704	! Preparing for the implicit scheme.
705	! This means loading the prognostic variables from the restart file.
706	!-
707	fluxsens = val_exp
708	CALL restget_p (rest_id, 'fluxsens', iim_g, jjm_g, 1, istp_old, .TRUE., fluxsens)
709	IF (ALL(fluxsens(:,:) == val_exp)) fluxsens(:,:) = zero
710	!-
711	vevapp = val_exp
712	CALL restget_p(rest_id, 'vevapp', iim_g, jjm_g, 1, istp_old, .TRUE., vevapp)
713	IF (ALL(vevapp(:,:) == val_exp)) vevapp(:,:) = zero
714	!-
715	old_zlev = val_exp
716	CALL restget_p (rest_id, 'zlev_old', iim_g, jjm_g, 1, istp_old, .TRUE., old_zlev)
717	IF (ALL(old_zlev(:,:) == val_exp)) old_zlev(:,:)=zlev_vec(:,:)
718	!-
719	old_qair = val_exp
720	CALL restget_p (rest_id, 'qair_old', iim_g, jjm_g, 1, istp_old, .TRUE., old_qair)
721	IF (ALL(old_qair(:,:) == val_exp)) old_qair(:,:) = qair_obs(:,:)
722	!-
723	old_eair = val_exp
724	CALL restget_p (rest_id, 'eair_old', iim_g, jjm_g, 1, istp_old, .TRUE., old_eair)
725	IF (ALL(old_eair(:,:) == val_exp)) THEN
726	old_eair = 0.0 ! Init value
727	DO ik=1,nbindex
728	i=ilandindex(ik)
729	j=jlandindex(ik)
730	old_eair(i,j) = cp_air * tair_obs(i,j) + cte_grav*zlev_vec(i,j)
731	ENDDO
732	ENDIF
733	!-
734	! old density is also needed because we do not yet have the right pb
735	!-
736	!=> obsolete ??!! (tjrs calcul apres forcing_read)
737	for_rau = val_exp
738	CALL restget_p (rest_id, 'rau_old', iim_g, jjm_g, 1, istp_old, .TRUE., for_rau)
739	IF (ALL(for_rau(:,:) == val_exp)) THEN
740	for_rau = fill_init
741	DO ik=1,nbindex
742	i=ilandindex(ik)
743	j=jlandindex(ik)
744	for_rau(i,j) = pb(i,j)/(cte_molr*(tair_obs(i,j)))
745	ENDDO
746	ENDIF
747	!-
748	! For this variable the restart is extracted by SECHIBA
749	!-
750	temp_sol_NEW(:,:) = tair_obs(:,:)
751	!-
752	if (.NOT. is_watchout) THEN
753	!-
754	! This does not yield a correct restart in the case of relaxation
755	!-
756	petAcoef = val_exp
757	CALL restget_p (rest_id, 'petAcoef', iim_g, jjm_g, 1, istp_old, .TRUE., petAcoef)
758	IF (ALL(petAcoef(:,:) == val_exp)) petAcoef(:,:) = zero
759	!--
760	petBcoef = val_exp
761	CALL restget_p (rest_id, 'petBcoef', iim_g, jjm_g, 1, istp_old, .TRUE., petBcoef)
762	IF (ALL(petBcoef(:,:) == val_exp)) petBcoef(:,:) = old_eair(:,:)
763	!--
764	peqAcoef = val_exp
765	CALL restget_p (rest_id, 'peqAcoef', iim_g, jjm_g, 1, istp_old, .TRUE., peqAcoef)
766	IF (ALL(peqAcoef(:,:) == val_exp)) peqAcoef(:,:) = zero
767	!--
768	peqBcoef = val_exp
769	CALL restget_p (rest_id, 'peqBcoef', iim_g, jjm_g, 1, istp_old, .TRUE., peqBcoef)
770	IF (ALL(peqBcoef(:,:) == val_exp)) peqBcoef(:,:) = old_qair(:,:)
771	ENDIF
772	!-
773	! And other variables which need initial variables. These variables
774	! will get properly initialized by ORCHIDEE when it is called for
775	! the first time.
776	!-
777	albedo(:,:,:) = 0.13
778	emis(:,:) = 1.0
779	z0(:,:) = 0.1
780	!--
781	!=====================================================================
782	!- 4.0 START THE TIME LOOP
783	!=====================================================================
784
785	julian = itau2date(istp, date0_rest, dt)
786	CALL ju2ymds(julian, yy, mm, dd, ss)
787
788
789	it = itau_dep+1
790	DO WHILE ( it <= itau_fin )
791	!----
792	it_force = it+for_offset
793	IF (it_force < 0) THEN
794	WRITE(numout,*) 'The day is not in the forcing file :', &
795	& it_force, it, for_offset
796	CALL ipslerr_p(3,'dim2_driver','Pb in forcing file','','')
797	ENDIF
798	is=split_start
799	DO WHILE ( is <= split )
800	!-----
801	yy_prev=yy
802	julian = itau2date(istp, date0_rest, dt)
803	CALL ju2ymds(julian, yy, mm, dd, ss)
804	IF (CO2_varying .AND. (yy /= yy_prev)) THEN
805	for_ccanopy(:,:)=for_ccanopy(:,:)*co2inc
806	ENDIF
807
808
809	IF (printlev_loc>=3) THEN
810	WRITE(numout,*) "=============================================================="
811	WRITE(numout,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") &
812	& yy,mm,dd,ss/3600.
813	#ifdef CPP_PARA
814	IF (is_root_prc) THEN
815	WRITE(,) "=============================================================="
816	WRITE(*,"('New iteration at date : ',I4,'-',I2.2,'-',I2.2,':',F8.4)") &
817	& yy,mm,dd,ss/3600.
818	ENDIF
819	#endif
820	ENDIF
821	!-----
822	IF ( (it == itau_fin).AND.(is == split) ) THEN
823	lstep_last = .TRUE.
824	ENDIF
825	!-----
826	IF (printlev_loc>=3) WRITE(numout,*) 'Into forcing_read'
827	!-----
828	CALL forcing_READ &
829	& (filename, rest_id, .FALSE., lstep_last, &
830	& it_force, istp, is, split, nb_spread, lwdown_cons, swdown_cons, &
831	& date0_rest, dt_force, iim, jjm, lon, lat, zlev_vec, zlevuv_vec, tm, &
832	& swdown, coszang, precip_rain, precip_snow, tair_obs, &
833	& u, v, qair_obs, pb, for_lwdown, for_contfrac, for_neighbours, for_resolution, &
834	& for_swnet, eair_obs, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, for_ccanopy, &
835	& kindex, nbindex, force_id)
836
837	!-----
838	!---- SECHIBA expects surface pressure in hPa
839	!-----
840	for_psurf(:,:) = pb(:,:)/100.
841
842	IF (printlev_loc>=4) THEN
843	WRITE(numout,*) "dim2_driver 0 ",it_force
844	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
845	WRITE(numout,*) "Lowest level wind speed North = ", &
846	& (/ ( u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
847	WRITE(numout,*) "Lowest level wind speed East = ", &
848	& (/ ( v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
849	WRITE(numout,*) "z0 ; Surface roughness = ", &
850	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
851	WRITE(numout,*) "Height of first layer = ", &
852	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
853	WRITE(numout,*) "Lowest level specific humidity = ", &
854	& (/ ( qair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
855	WRITE(numout,*) "Rain precipitation = ", &
856	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
857	WRITE(numout,*) "Snow precipitation = ", &
858	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
859	WRITE(numout,*) "Down-welling long-wave flux = ", &
860	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
861	WRITE(numout,*) "Net surface short-wave flux = ", &
862	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
863	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
864	& (/ ( swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
865	WRITE(numout,*) "Air temperature in Kelvin = ", &
866	& (/ ( tair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
867	WRITE(numout,*) "Air potential energy = ", &
868	& (/ ( eair_obs(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
869	WRITE(numout,*) "CO2 concentration in the canopy = ", &
870	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
871	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
872	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
873	WRITE(numout,*) "One for T and another for q = ", &
874	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
875	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
876	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
877	WRITE(numout,*) "One for T and another for q = ", &
878	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
879	WRITE(numout,*) "Cdrag = ", &
880	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
881	WRITE(numout,*) "CO2 concentration in the canopy = ", &
882	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
883	WRITE(numout,*) "Lowest level pressure = ", &
884	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
885	WRITE(numout,*) "Geographical coordinates lon = ", &
886	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
887	WRITE(numout,*) "Geographical coordinates lat = ", &
888	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
889	WRITE(numout,*) "Fraction of continent in the grid = ", &
890	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
891	ENDIF
892	!-----
893	!---- Prepare : tmp_qair, tmp_eair, tmp_tair, tmp_pb
894	!---- and : for_u, for_v, for_lwdown, for_swnet, for_swdown
895	!---- All the work is done in forcing_read
896	IF (printlev_loc>=3) WRITE(numout,*) 'Prepare the atmospheric forcing'
897	!-----
898	IF (.NOT. is_watchout) THEN
899	DO ik=1,nbindex
900	i=ilandindex(ik)
901	j=jlandindex(ik)
902	eair_obs(i,j) = cp_airtair_obs(i,j)+cte_gravzlev_vec(i,j)
903	for_swnet(i,j) = (1.-(albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j)
904	ENDDO
905	ENDIF
906	DO ik=1,nbindex
907	i=ilandindex(ik)
908	j=jlandindex(ik)
909	for_swdown(i,j) = swdown(i,j)
910	for_coszang(i,j) = coszang(i,j)
911	ENDDO
912	!-----
913	!---- Computing the buffer zone !
914	!-----
915	IF (relaxation) THEN
916	DO ik=1,nbindex
917	i=ilandindex(ik)
918	j=jlandindex(ik)
919	for_qair(i,j) = peqAcoef(i,j)(-1.) vevapp(i,j)*dt+peqBcoef(i,j)
920	!-------
921	for_eair(i,j) = petAcoef(i,j)(-1.) fluxsens(i,j)+petBcoef(i,j)
922	!-------
923	ENDDO
924	DO ik=1,nbindex
925	i=ilandindex(ik)
926	j=jlandindex(ik)
927	for_tair(i,j) = (for_eair(i,j) - cte_grav*zlev_vec(i,j))/cp_air
928	!-------
929	!!$ if (.NOT. is_watchout) &
930	!!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:)
931	!-------
932	ENDDO
933	DO ik=1,nbindex
934	i=ilandindex(ik)
935	j=jlandindex(ik)
936	for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j))
937	!-------
938	relax(i,j) = for_rau(i,j)*alpha
939	ENDDO
940
941	DO ik=1,nbindex
942	i=ilandindex(ik)
943	j=jlandindex(ik)
944	zlflu = zlev_vec(i,j)/2.0*dt
945	peqAcoef(i,j) = 1.0/(zlflu+relax(i,j))
946	peqBcoef(i,j) = (relax(i,j) * qair_obs(i,j)/(zlflu+relax(i,j))) + &
947	& for_qair(i,j)/(1.0+relax(i,j)/zlflu)
948	ENDDO
949	!-------
950	! relax(:,:) = for_rau(:,:)*alpha
951	DO ik=1,nbindex
952	i=ilandindex(ik)
953	j=jlandindex(ik)
954	petAcoef(i,j) = 1.0/(zlflu+relax(i,j))
955	petBcoef(i,j) = ( relax(i,j) * eair_obs(i,j) / (zlflu+relax(i,j)) ) &
956	& + for_eair(i,j)/(1.0+relax(i,j)/zlflu)
957	ENDDO
958	ELSE
959	for_qair(:,:) = fill_init
960	for_eair(:,:) = fill_init
961	for_tair(:,:) = fill_init
962	DO ik=1,nbindex
963	i=ilandindex(ik)
964	j=jlandindex(ik)
965	for_qair(i,j) = qair_obs(i,j)
966	for_eair(i,j) = eair_obs(i,j)
967	for_tair(i,j) = tair_obs(i,j)
968	ENDDO
969	!-------
970	!!$ if (.NOT. is_watchout) &
971	!!$ epot_sol(:,:) = cp_air*temp_sol_NEW(:,:)
972	!-------
973	DO ik=1,nbindex
974	i=ilandindex(ik)
975	j=jlandindex(ik)
976	for_rau(i,j) = pb(i,j) / (cte_molr*for_tair(i,j))
977	ENDDO
978	!-------
979	IF (.NOT. is_watchout) THEN
980	petAcoef(:,:) = 0.0
981	peqAcoef(:,:) = 0.0
982	DO ik=1,nbindex
983	i=ilandindex(ik)
984	j=jlandindex(ik)
985	petBcoef(i,j) = eair_obs(i,j)
986	peqBcoef(i,j) = qair_obs(i,j)
987	ENDDO
988	ENDIF
989	ENDIF
990	!-----
991	IF (.NOT. is_watchout) &
992	cdrag(:,:) = 0.0
993	! for_ccanopy(:,:)=atmco2
994	!-----
995	!---- SECHIBA expects wind, temperature and humidity at the same height.
996	!---- If this is not the case then we need to correct for that.
997	!-----
998	DO ik=1,nbindex
999	i=ilandindex(ik)
1000	j=jlandindex(ik)
1001	for_u(i,j) = u(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / &
1002	LOG(zlevuv_vec(i,j)/z0(i,j))
1003	for_v(i,j) = v(i,j)*LOG(zlev_vec(i,j)/z0(i,j)) / &
1004	LOG(zlevuv_vec(i,j)/z0(i,j))
1005	END DO
1006
1007	!-----
1008	!---- Prepare the other variables WITH the special CASE
1009	!---- of splited time steps
1010	!----
1011	!---- PRINT input value for printlev_loc>=3
1012	!-----
1013	IF (printlev_loc>=3) THEN
1014	WRITE(numout,*) ' >>>>>> time step it_force = ',it_force
1015	WRITE(numout,*) &
1016	& ' tair, qair, eair = ', &
1017	& for_tair(itest,jtest),for_qair(itest,jtest), &
1018	& for_eair(itest,jtest)
1019	WRITE(numout,*) &
1020	& ' OBS : tair, qair, eair = ', &
1021	& tair_obs(itest,jtest),qair_obs(itest,jtest), &
1022	& eair_obs(itest,jtest)
1023	WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest)
1024	WRITE(numout,*) ' precip rain et snow = ', &
1025	& precip_rain(itest,jtest),precip_snow(itest,jtest)
1026	WRITE(numout,*) ' lwdown et swnet = ', &
1027	& for_lwdown(itest,jtest),for_swnet(itest,jtest)
1028	WRITE(numout,*) ' petAcoef et peqAcoef = ', &
1029	& petAcoef(itest,jtest), peqAcoef(itest,jtest)
1030	WRITE(numout,*) ' petBcoef et peqAcoef = ', &
1031	& petBcoef(itest,jtest),peqBcoef(itest,jtest)
1032	WRITE(numout,*) ' zlev = ',zlev_vec(itest,jtest)
1033	ENDIF
1034	!-----
1035
1036	IF (lstep_init) THEN
1037
1038	DO ik=1,nbindex
1039	i=ilandindex(ik)
1040	j=jlandindex(ik)
1041	for_swdown(i,j) = swdown(i,j)
1042	for_coszang(i,j) = coszang(i,j)
1043	ENDDO
1044	IF (printlev_loc>=4) THEN
1045	WRITE(numout,*) "dim2_driver lstep_init ",it_force
1046	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
1047	WRITE(numout,*) "Lowest level wind speed North = ", &
1048	& (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1049	WRITE(numout,*) "Lowest level wind speed East = ", &
1050	& (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1051	WRITE(numout,*) "z0 ; Surface roughness = ", &
1052	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1053	WRITE(numout,*) "Height of first layer = ", &
1054	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1055	WRITE(numout,*) "Lowest level specific humidity = ", &
1056	& (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1057	WRITE(numout,*) "Rain precipitation = ", &
1058	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
1059	WRITE(numout,*) "Snow precipitation = ", &
1060	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
1061	WRITE(numout,*) "Down-welling long-wave flux = ", &
1062	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1063	WRITE(numout,*) "Net surface short-wave flux = ", &
1064	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1065	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
1066	& (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1067	WRITE(numout,*) "Air temperature in Kelvin = ", &
1068	& (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1069	WRITE(numout,*) "Air potential energy = ", &
1070	& (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1071	WRITE(numout,*) "CO2 concentration in the canopy = ", &
1072	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1073	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
1074	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1075	WRITE(numout,*) "One for T and another for q = ", &
1076	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1077	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
1078	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1079	WRITE(numout,*) "One for T and another for q = ", &
1080	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1081	WRITE(numout,*) "Cdrag = ", &
1082	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1083	WRITE(numout,*) "Lowest level pressure = ", &
1084	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1085	WRITE(numout,*) "Geographical coordinates lon = ", &
1086	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1087	WRITE(numout,*) "Geographical coordinates lat = ", &
1088	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1089	WRITE(numout,*) "Fraction of continent in the grid = ", &
1090	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1091	ENDIF
1092	!-------
1093	!------ CALL sechiba to initialize fields
1094	!------ and have some initial results: emis, albedo, z0
1095	!-------
1096	CALL intersurf_initialize_2d &
1097	& (istp_old, iim, jjm, nbindex, kindex, dt, &
1098	& lstep_init, .FALSE., lon, lat, for_contfrac, for_resolution, date0_rest, &
1099	! first level conditions
1100	& zlev_vec, for_u, for_v, &
1101	& for_qair, for_tair, for_eair, for_ccanopy, &
1102	! Variables for the implicit coupling
1103	& cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1104	! Rain, snow, radiation and surface pressure
1105	& precip_rain, precip_snow, &
1106	& for_lwdown, for_swnet, for_swdown, for_psurf, &
1107	! Output : Fluxes
1108	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
1109	! Surface temperatures and surface properties
1110	& tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0, coszang )
1111
1112	CALL Stop_timer(timer_global)
1113	CALL Stop_timer(timer_mpi)
1114	CALL Start_timer(timer_global)
1115	CALL Start_timer(timer_mpi)
1116	!
1117	lstep_init = .FALSE.
1118	!
1119	! Get Restart values for albedo and z0,
1120	! as they modify forcing variables swnet and wind.
1121	!-------
1122	! albedo
1123	albedo_vis = val_exp
1124	CALL restget_p (rest_id, 'albedo_vis', iim_g, jjm_g, 1, istp_old, .TRUE., albedo_vis)
1125	IF (ALL(albedo_vis(:,:) == val_exp)) THEN
1126	albedo_vis(:,:)=albedo(:,:,1)
1127	ELSE
1128	albedo(:,:,1)=albedo_vis(:,:)
1129	ENDIF
1130	!
1131	albedo_nir = val_exp
1132	CALL restget_p (rest_id, 'albedo_nir', iim_g, jjm_g, 1, istp_old, .TRUE., albedo_nir)
1133	IF (ALL(albedo_nir(:,:) == val_exp)) THEN
1134	albedo_nir(:,:)=albedo(:,:,2)
1135	ELSE
1136	albedo(:,:,2)=albedo_nir(:,:)
1137	ENDIF
1138
1139	CALL restget_p (rest_id, 'z0', iim_g, jjm_g, 1, istp_old, .TRUE., tmp_z0)
1140	IF (.NOT. ALL(tmp_z0 == val_exp)) THEN
1141	! 'z0' was found in restart file. Redefine z0 with this value.
1142	z0 = tmp_z0
1143	END IF
1144	!-------
1145	DO ik=1,nbindex
1146	i=ilandindex(ik)
1147	j=jlandindex(ik)
1148	temp_sol_old(i,j) = temp_sol_NEW(i,j)
1149	for_swnet(i,j) = (1.- (albedo(i,j,1)+albedo(i,j,2))/2.)*swdown(i,j)
1150	for_swdown(i,j) = swdown(i,j)
1151	for_coszang(i,j) = coszang(i,j)
1152	ENDDO
1153	!
1154	! MM : z0 have been modified then we must lower the wind again
1155	!-----
1156	!---- SECHIBA expects wind, temperature and humidity at the same height.
1157	!---- If this is not the case then we need to correct for that.
1158	!-----
1159	DO ik=1,nbindex
1160	i=ilandindex(ik)
1161	j=jlandindex(ik)
1162	for_u(i,j) = u(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / &
1163	LOG(zlevuv_vec(i,j)/z0(i,j))
1164	for_v(i,j) = v(i,j) * LOG(zlev_vec(i,j)/z0(i,j)) / &
1165	LOG(zlevuv_vec(i,j)/z0(i,j))
1166	END DO
1167
1168	!-----
1169	!---- PRINT input value after lstep_init for printlev_loc>=3
1170	!-----
1171	IF (printlev_loc>=3) THEN
1172	WRITE(numout,*) ' >>>>>> after lstep_init = ',lstep_init
1173	WRITE(numout,*) ' u et v = ',for_u(itest,jtest),for_v(itest,jtest)
1174	WRITE(numout,*) ' swnet = ', for_swnet(itest,jtest)
1175	ENDIF
1176	!-------
1177	IF (printlev_loc>=4) THEN
1178	WRITE(numout,*) "dim2_driver lstep_init outputs"
1179	! Output : Fluxes
1180	WRITE(numout,*) "vevapp ; Total of evaporation = ", &
1181	& (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1182	WRITE(numout,*) "Sensible heat flux = ", &
1183	& (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1184	WRITE(numout,*) "Latent heat flux = ", &
1185	& (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1186	WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", &
1187	& (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1188	WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", &
1189	& (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1190	! Surface temperatures and surface properties
1191	WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", &
1192	& (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1193	WRITE(numout,*) "temp_sol_new ; New soil temperature = ", &
1194	& (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1195	WRITE(numout,*) "qsurf ; Surface specific humidity = ", &
1196	& (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1197	WRITE(numout,*) "albedoVIS = ", &
1198	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /)
1199	WRITE(numout,*) "albedoNIR = ", &
1200	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /)
1201	WRITE(numout,*) "emis ; Emissivity = ", &
1202	& (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1203	WRITE(numout,*) "z0 ; Surface roughness = ", &
1204	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1205	ENDIF
1206	!-------
1207	IF (printlev_loc>=3) THEN
1208	WRITE(numout,*) &
1209	& ' OUT rest : z0, albedoVIS, albedoNIR, emis = ', &
1210	& z0(itest,jtest),albedo(itest,jtest,1), &
1211	& albedo(itest,jtest,2),emis(itest,jtest)
1212	WRITE(numout,*) ' OUT rest : coastal and river flow = ', &
1213	& coastalflow(itest,jtest), riverflow(itest,jtest)
1214	WRITE(numout,*) ' OUT rest : tsol_rad, vevapp = ', &
1215	& tsol_rad(itest,jtest), vevapp(itest,jtest)
1216	WRITE(numout,*) ' OUT rest : temp_sol_new =', &
1217	& temp_sol_NEW(itest,jtest)
1218	ENDIF
1219
1220	ENDIF ! lstep_init
1221	!-----
1222	!---- Calling SECHIBA and doing the number crunching.
1223	!---- Note that for the first time step SECHIBA is called twice.
1224	!----
1225	!---- All H_2O fluxes are now in Kg/m^2s
1226	!-----
1227	IF (printlev_loc>=4) THEN
1228	WRITE(numout,*) "dim2_driver ",it_force
1229	WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex)
1230	WRITE(numout,*) "Lowest level wind speed North = ", &
1231	& (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1232	WRITE(numout,*) "Lowest level wind speed East = ", &
1233	& (/ ( for_v(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1234	WRITE(numout,*) "z0 ; Surface roughness = ", &
1235	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1236	WRITE(numout,*) "Height of first layer = ", &
1237	& (/ ( zlev_vec(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1238	WRITE(numout,*) "Lowest level specific humidity = ", &
1239	& (/ ( for_qair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1240	WRITE(numout,*) "Rain precipitation = ", &
1241	& (/ ( precip_rain(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
1242	WRITE(numout,*) "Snow precipitation = ", &
1243	& (/ ( precip_snow(ilandindex(ik), jlandindex(ik))*dt,ik=1,nbindex ) /)
1244	WRITE(numout,*) "Down-welling long-wave flux = ", &
1245	& (/ ( for_lwdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1246	WRITE(numout,*) "Net surface short-wave flux = ", &
1247	& (/ ( for_swnet(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1248	WRITE(numout,*) "Downwelling surface short-wave flux = ", &
1249	& (/ ( for_swdown(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1250	WRITE(numout,*) "Air temperature in Kelvin = ", &
1251	& (/ ( for_tair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1252	WRITE(numout,*) "Air potential energy = ", &
1253	& (/ ( for_eair(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1254	WRITE(numout,*) "CO2 concentration in the canopy = ", &
1255	& (/ ( for_ccanopy(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1256	WRITE(numout,*) "Coeficients A from the PBL resolution = ", &
1257	& (/ ( petAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1258	WRITE(numout,*) "One for T and another for q = ", &
1259	& (/ ( peqAcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1260	WRITE(numout,*) "Coeficients B from the PBL resolution = ", &
1261	& (/ ( petBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1262	WRITE(numout,*) "One for T and another for q = ", &
1263	& (/ ( peqBcoef(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1264	WRITE(numout,*) "Cdrag = ", &
1265	& (/ ( cdrag(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1266	WRITE(numout,*) "Lowest level pressure = ", &
1267	& (/ ( for_psurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1268	WRITE(numout,*) "Geographical coordinates lon = ", &
1269	& (/ ( lon(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1270	WRITE(numout,*) "Geographical coordinates lat = ", &
1271	& (/ ( lat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1272	WRITE(numout,*) "Fraction of continent in the grid = ", &
1273	& (/ ( for_contfrac(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1274	ENDIF
1275
1276	CALL intersurf_main_2d &
1277	& (istp, iim, jjm, nbindex, kindex, dt, &
1278	& lstep_init, lstep_last, lon, lat, for_contfrac, for_resolution, date0_rest, &
1279	! first level conditions
1280	& zlev_vec, for_u, for_v, &
1281	& for_qair, for_tair, for_eair, for_ccanopy, &
1282	! Variables for the implicit coupling
1283	& cdrag, petAcoef, peqAcoef, petBcoef, peqBcoef, &
1284	! Rain, snow, radiation and surface pressure
1285	& precip_rain, precip_snow, &
1286	& for_lwdown, for_swnet, for_swdown, for_psurf, &
1287	! Output : Fluxes
1288	& vevapp, fluxsens, fluxlat, coastalflow, riverflow, &
1289	! Surface temperatures and surface properties
1290	& tsol_rad, temp_sol_NEW, qsurf, albedo, emis, z0, &
1291	! VOC : radiation
1292	& for_coszang)
1293
1294	!-------
1295	IF (printlev_loc>=4) THEN
1296	WRITE(numout,*) "dim2_driver outputs"
1297	! Output : Fluxes
1298	WRITE(numout,*) "vevapp ; Total of evaporation = ", &
1299	& (/ ( vevapp(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1300	WRITE(numout,*) "Sensible heat flux = ", &
1301	& (/ ( fluxsens(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1302	WRITE(numout,*) "Latent heat flux = ", &
1303	& (/ ( fluxlat(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1304	WRITE(numout,*) "coastalflow ; Diffuse flow of water into the ocean (m^3/dt) = ", &
1305	& (/ ( coastalflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1306	WRITE(numout,*) "riverflow ; Largest rivers flowing into the ocean (m^3/dt) = ", &
1307	& (/ ( riverflow(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1308	! Surface temperatures and surface properties
1309	WRITE(numout,*) "tsol_rad ; Radiative surface temperature = ", &
1310	& (/ ( tsol_rad(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1311	WRITE(numout,*) "temp_sol_new ; New soil temperature = ", &
1312	& (/ ( temp_sol_NEW(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1313	WRITE(numout,*) "qsurf ; Surface specific humidity = ", &
1314	& (/ ( qsurf(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1315	WRITE(numout,*) "albedoVIS = ", &
1316	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 1),ik=1,nbindex ) /)
1317	WRITE(numout,*) "albedoNIR = ", &
1318	& (/ ( albedo(ilandindex(ik), jlandindex(ik), 2),ik=1,nbindex ) /)
1319	WRITE(numout,*) "emis ; Emissivity = ", &
1320	& (/ ( emis(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1321	WRITE(numout,*) "z0 ; Surface roughness = ", &
1322	& (/ ( z0(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /)
1323	ENDIF
1324	!-----
1325	dtdt(:,:) = zero
1326	DO ik=1,nbindex
1327	i=ilandindex(ik)
1328	j=jlandindex(ik)
1329	dtdt(i,j) = ABS(temp_sol_NEW(i,j)-temp_sol_old(i,j))/dt
1330	ENDDO
1331	!-----
1332	!---- Test if the point with the largest change has more than 5K per dt
1333	!-----
1334	IF (printlev_loc >=3) THEN
1335	IF (MAXVAL(dtdt(:,:)) > 5./dt) THEN
1336	ml = MAXLOC(dtdt)
1337	CALL ju2ymds(julian, yy, mm, dd, ss)
1338	WRITE(numout,"('ATT :',I5,' big temperature jumps on ', &
1339	I4,'-',I2.2,'-',I2.2,':',F8.4)") &
1340	COUNT(dtdt(:,:) > 5./dt),yy,mm,dd,ss/3600.
1341	WRITE(numout,*) &
1342	'Maximum change of surface temperature located at :', &
1343	lon(ml(1),ml(2)),lat(ml(1),ml(2))
1344	WRITE(numout,*) 'Coordinates in grid space: ',ml(1),ml(2)
1345	WRITE(numout,*) 'Change from ',temp_sol_old(ml(1),ml(2)), &
1346	' to ',temp_sol_new(ml(1),ml(2)),&
1347	'with sw_in = ',for_swnet(ml(1),ml(2))
1348	old_tair = &
1349	(old_eair(ml(1),ml(2))-cte_grav*old_zlev(ml(1),ml(2)))/cp_air
1350	WRITE(numout,*) 'Air temperature change from ',old_tair, &
1351	' to ',for_tair(ml(1),ml(2))
1352	WRITE(numout,*) 'Max of dtdt : ',dtdt(ml(1),ml(2)),' with dt = ',dt
1353	ENDIF
1354	END IF
1355
1356	temp_sol_old(:,:) = temp_sol_NEW(:,:)
1357	!-----
1358	!---- PRINT output value for printlev_loc>=3
1359	!-----
1360	IF (printlev_loc>=3) THEN
1361	WRITE(numout,*) ' OUT : z0, albedoVIS, albedoNIR, emis = ', &
1362	& z0(itest,jtest),albedo(itest,jtest,1), &
1363	& albedo(itest,jtest,2),emis(itest,jtest)
1364	WRITE(numout,*) ' OUT : coastal and river flow = ',&
1365	& coastalflow(itest,jtest), riverflow(itest,jtest)
1366	WRITE(numout,*) ' OUT : tsol_rad, vevapp = ', &
1367	& tsol_rad(itest,jtest), vevapp(itest,jtest)
1368	WRITE(numout,*) ' OUT : temp_sol_new =', temp_sol_NEW(itest,jtest)
1369	ENDIF
1370	!-----
1371	!---- Give some variables to the output package
1372	!---- for saving on the history tape
1373	!-----
1374	IF (printlev_loc>=3) WRITE(numout,*) 'history written for ', istp
1375	!-----
1376	istp_old = istp
1377	istp = istp+1
1378	!-----
1379	old_zlev(:,:) = zlev_vec(:,:)
1380	old_qair(:,:) = for_qair(:,:)
1381	old_eair(:,:) = for_eair(:,:)
1382	!-----
1383	is = is + 1
1384	ENDDO ! DO WHILE (is <= split)
1385	split_start = 1
1386
1387	IF (it==itau_fin-1) THEN
1388
1389	CALL Write_Load_Balance(REAL(Get_cpu_time(timer_mpi),r_std))
1390
1391	ENDIF
1392	it = it + 1
1393	ENDDO ! DO WHILE (it <= itau_fin)
1394
1395	!-
1396	! Archive in restart file the prognostic variables
1397	!-
1398	IF (printlev_loc>=3) WRITE(numout,*) 'Write the restart for the driver', istp_old
1399	!-
1400	CALL restput_p (rest_id, 'fluxsens', iim_g, jjm_g, 1, istp_old, fluxsens)
1401	CALL restput_p (rest_id, 'vevapp', iim_g, jjm_g, 1, istp_old, vevapp)
1402	CALL restput_p (rest_id, 'zlev_old', iim_g, jjm_g, 1, istp_old, old_zlev)
1403	CALL restput_p (rest_id, 'qair_old', iim_g, jjm_g, 1, istp_old, old_qair)
1404	CALL restput_p (rest_id, 'eair_old', iim_g, jjm_g, 1, istp_old, old_eair)
1405	CALL restput_p (rest_id, 'rau_old', iim_g, jjm_g, 1, istp_old, for_rau)
1406	CALL restput_p (rest_id, 'albedo_vis', iim_g, jjm_g, 1, istp_old, albedo(:,:,1))
1407	CALL restput_p (rest_id, 'albedo_nir', iim_g, jjm_g, 1, istp_old, albedo(:,:,2))
1408	CALL restput_p (rest_id, 'z0', iim_g, jjm_g, 1, istp_old, z0)
1409
1410	if (.NOT. is_watchout) THEN
1411	CALL restput_p (rest_id, 'petAcoef', iim_g, jjm_g, 1, istp_old, petAcoef)
1412	CALL restput_p (rest_id, 'petBcoef', iim_g, jjm_g, 1, istp_old, petBcoef)
1413	CALL restput_p (rest_id, 'peqAcoef', iim_g, jjm_g, 1, istp_old, peqAcoef)
1414	CALL restput_p (rest_id, 'peqBcoef', iim_g, jjm_g, 1, istp_old, peqBcoef)
1415	ENDIF
1416	!-
1417	IF (printlev_loc>=3) WRITE(numout,*) 'Restart for the driver written'
1418	!=====================================================================
1419	!- 5.0 Closing all files
1420	!=====================================================================
1421	CALL flinclo(force_id)
1422	IF ( printlev_loc>=3 ) WRITE(numout,*) 'FLIN CLOSED'
1423	CALL histclo
1424	IF ( printlev_loc>=3 ) WRITE(numout,*) 'HIST CLOSED'
1425
1426	IF(is_root_prc) THEN
1427	CALL restclo
1428	IF ( printlev_loc>=3 ) WRITE(numout,*) 'REST CLOSED'
1429	CALL getin_dump
1430	IF ( printlev_loc>=3 ) WRITE(numout,*) 'GETIN CLOSED'
1431	ENDIF
1432
1433
1434
1435	WRITE(numout,*) '-------------------------------------------'
1436	WRITE(numout,*) '------> CPU Time Global ',Get_cpu_Time(timer_global)
1437	WRITE(numout,*) '------> CPU Time without mpi ',Get_cpu_Time(timer_mpi)
1438	WRITE(numout,*) '------> Real Time Global ',Get_real_Time(timer_global)
1439	WRITE(numout,*) '------> real Time without mpi ',Get_real_Time(timer_mpi)
1440	WRITE(numout,*) '-------------------------------------------'
1441
1442	! Call driver_clear for deallocation and reset of initialization variables
1443	CALL driver_clear
1444
1445	CALL Finalize_mpi
1446
1447
1448	WRITE(numout,*) 'END of dim2_driver'
1449
1450
1451	CONTAINS
1452
1453	!! ================================================================================================================================
1454	!! SUBROUTINE : driver_clear
1455	!!
1456	!>\BRIEF Clear driver main program and call clear funcions for underlaying module intersurf
1457	!!
1458	!! DESCRIPTION : Deallocate memory and reset initialization variables to there original values
1459	!! This subroutine call intersurf_clear which will call sechiba_clear.
1460	!!
1461	!_ ================================================================================================================================
1462	SUBROUTINE driver_clear
1463
1464	CALL intersurf_clear
1465
1466	END SUBROUTINE driver_clear
1467
1468	END PROGRAM driver

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

Download in other formats: