Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_driver/dim2_driver.f90

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

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

Download in other formats: