Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE_OL/dim2_driver.f90 @ 183

Last change on this file since 183 was 65, checked in by didier.solyga, 14 years ago
Import first version of ORCHIDEE_EXT
File size: 62.2 KB

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

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

Download in other formats: