Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE_OL/dim2_driver.f90 @ 6268

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

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

Download in other formats: