Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dim2_driver.f90 @ 8367

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

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

Context Navigation

source: branches/ORCHIDEE_3_CMIP6/ORCHIDEE/src_driver/dim2_driver.f90 @ 8367

Download in other formats: