Context Navigation

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

Last change on this file since 5804 was 4981, checked in by josefine.ghattas, 6 years ago
Changes in comments to avoid ATM_CO2 be written twice in orchidee.default
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 64.7 KB

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

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

Download in other formats: