Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE_OL/dim2_driver.f90 @ 118

Last change on this file since 118 was 48, checked in by mmaipsl, 14 years ago

MM: Tests with lf95 compiler : correct f95 strict norm problems.

There is no change in numerical result after these commits.

Property svn:keywords set to HeadURL Date Author Revision

File size: 62.2 KB

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

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

Download in other formats: