Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE_OL/dim2_driver.f90 @ 4639

Last change on this file since 4639 was 847, checked in by didier.solyga, 12 years ago
Formatted labels so a script can automatically generate the orchidee.default file for ORCHIDEE_OL modules.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 64.0 KB

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

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

Download in other formats: