Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_driver/dim2_driver.f90 @ 5520

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

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

Download in other formats: