Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dim2_driver.f90 @ 6715

Last change on this file since 6715 was 6714, checked in by yuan.zhang, 4 years ago
ORCHIDEE_DF with light partitioning at site level
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 64.8 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE_DFv1.0_site/src_driver/dim2_driver.f90 @ 6715

Download in other formats: