Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_driver/dim2_driver.f90 @ 5925

Last change on this file since 5925 was 3876, checked in by ronny.lauerwald, 8 years ago
Merge ORCHIDOC into TRUNK
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 65.8 KB

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

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

Download in other formats: