Context Navigation

source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_driver/dim2_driver.f90 @ 8398

Last change on this file since 8398 was 4287, checked in by josefine.ghattas, 7 years ago
Enhencement on clear functions. Added call to clear functions from all offline drivers. See ticket #232
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 65.4 KB

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

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

Download in other formats: