Context Navigation

dim2_driver.f90

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

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

Context Navigation

source: tags/ORCHIDEE_1_9_5/ORCHIDEE_OL/dim2_driver.f90

Download in other formats: