Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/dim2_driver.f90 @ 7262

Last change on this file since 7262 was 7256, checked in by agnes.ducharne, 3 years ago
Reduced xios compression for speed and correction to run leap years with dim2_driver, following r5855 and r5867 of the trunk.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 64.8 KB

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

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

Download in other formats: