Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_driver/readdim2.f90 @ 8692

Last change on this file since 8692 was 4333, checked in by josefine.ghattas, 8 years ago
Merge with trunk for revision r3623 - r3977 Some corrections for gfortran
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 97.6 KB

Line
1
2	MODULE readdim2
3	!---------------------------------------------------------------------
4	!< $HeadURL$
5	!< $Date$
6	!< $Author$
7	!< $Revision$
8	!- IPSL (2006)
9	!- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10	!-
11	USE ioipsl_para
12	USE weather
13	USE TIMER
14	USE constantes
15	USE solar
16	USE grid
17	USE mod_orchidee_para
18	!-
19	IMPLICIT NONE
20	!-
21	PRIVATE
22	PUBLIC :: forcing_read, forcing_info, forcing_grid
23	!-
24	INTEGER, SAVE :: iim_full, jjm_full, llm_full, ttm_full
25	INTEGER, SAVE :: iim_zoom, jjm_zoom
26	INTEGER, SAVE :: iim_g_begin,jjm_g_begin,iim_g_end,jjm_g_end
27	REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: data_full, lon_full, lat_full
28	REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: lev_full
29	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: itau, i_index, j_index,j_index_g
30	INTEGER, SAVE :: i_test, j_test
31	INTEGER, SAVE :: printlev_loc !! Local printlev
32	LOGICAL, SAVE :: allow_weathergen, interpol, daily_interpol
33	LOGICAL, SAVE, PUBLIC :: weathergen, is_watchout
34	REAL, SAVE :: merid_res, zonal_res
35	LOGICAL, SAVE :: have_zaxis=.FALSE.
36	!-
37	!- Heigh controls and data
38	!-
39	LOGICAL, SAVE :: zfixed, zsigma, zhybrid, zlevels, zheight
40	LOGICAL, SAVE :: zsamelev_uv
41	REAL, SAVE :: zlev_fixed, zlevuv_fixed
42	REAL, SAVE :: zhybrid_a, zhybrid_b
43	REAL, SAVE :: zhybriduv_a, zhybriduv_b
44
45	CONTAINS
46	!-
47	!=====================================================================
48	!-
49	SUBROUTINE forcing_info(filename, iim, jjm, llm, tm, date0, dt_force,&
50	& force_id)
51
52	!---------------------------------------------------------------------
53	!
54	!- This subroutine will get all the info from the forcing file and
55	!- prepare for the zoom if needed.
56	!- It returns to the caller the sizes of the data it will receive at
57	!- the forcing_read call. This is important so that the caller can
58	!- allocate the right space.
59	!-
60	!- filename : name of the file to be opened
61	!- iim : size in x of the forcing data
62	!- jjm : size in y of the forcing data
63	!- llm : number of levels in the forcing data (not yet used)
64	!- tm : Time dimension of the forcing
65	!- date0 : The date at which the forcing file starts (julian days)
66	!- dt_force : time-step of the forcing file in seconds
67	!- force_id : ID of the forcing file
68	!-
69	!- ARGUMENTS
70	!-
71	IMPLICIT NONE
72	!-
73	CHARACTER(LEN=*) :: filename
74	INTEGER :: iim, jjm, llm, tm
75	REAL :: date0, dt_force
76	INTEGER, INTENT(OUT) :: force_id
77	!- LOCAL
78	CHARACTER(LEN=20) :: calendar_str
79	CHARACTER(LEN=200):: printstr !! temporary character string to contain error message
80	REAL :: juld_1, juld_2
81	REAL, ALLOCATABLE, DIMENSION(:,:) :: fcontfrac
82	REAL, ALLOCATABLE, DIMENSION(:,:) :: qair
83	LOGICAL :: contfrac_exists=.FALSE.
84	INTEGER :: NbPoint
85	INTEGER :: i_test,j_test
86	INTEGER :: i,j,ind,ttm_part
87	INTEGER, ALLOCATABLE, DIMENSION(:) :: index_l
88	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
89	REAL, ALLOCATABLE, DIMENSION(:) :: lev, levuv
90
91	!-
92	CALL flininfo(filename, iim_full, jjm_full, llm_full, ttm_full, force_id)
93	!-
94	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_info : Details from forcing file :', &
95	iim_full, jjm_full, llm_full, ttm_full
96	!-
97	IF ( llm_full < 1 ) THEN
98	have_zaxis = .FALSE.
99	ELSE
100	have_zaxis = .TRUE.
101	ENDIF
102	WRITE(numout,*) 'have_zaxis : ', llm_full, have_zaxis
103	!-
104	ttm_part = 2
105	ALLOCATE(itau(ttm_part))
106	ALLOCATE(data_full(iim_full, jjm_full),lon_full(iim_full, jjm_full),&
107	& lat_full(iim_full, jjm_full))
108	ALLOCATE(lev_full(llm_full))
109	ALLOCATE(fcontfrac(iim_full,jjm_full))
110	!-
111	lev_full(:) = zero
112	!-
113	dt_force=zero
114	CALL flinopen &
115	& (filename, .FALSE., iim_full, jjm_full, llm_full, lon_full, lat_full, &
116	& lev_full, ttm_part, itau, date0, dt_force, force_id)
117	IF ( dt_force == zero ) THEN
118	dt_force = itau(2) - itau(1)
119	itau(:) = itau(:) / dt_force
120	ENDIF
121	! WRITE(numout,*) "forcing_info : Forcing time step out of flinopen : ",dt_force
122	!-
123	!- What are the alowed options for the temportal interpolation
124	!-
125	!Config Key = ALLOW_WEATHERGEN
126	!Config Desc = Allow weather generator to create data
127	!Config If = [-]
128	!Config Def = n
129	!Config Help = This flag allows the forcing-reader to generate
130	!Config synthetic data if the data in the file is too sparse
131	!Config and the temporal resolution would not be enough to
132	!Config run the model.
133	!Config Units = [FLAG]
134	!-
135	allow_weathergen = .FALSE.
136	CALL getin_p('ALLOW_WEATHERGEN',allow_weathergen)
137	!-
138	!- The calendar was set by the forcing file. If no "calendar" attribute was
139	!- found then it is assumed to be gregorian,
140	!MM => FALSE !! it is NOT assumed anything !
141	!- else it is what ever is written in this attribute.
142	!-
143	CALL ioget_calendar(calendar_str)
144	i=INDEX(calendar_str,ACHAR(0))
145	IF ( i > 0 ) THEN
146	calendar_str(i:20)=' '
147	ENDIF
148	! WRITE(numout,*) "forcing_info : Calendar used : ",calendar_str
149	IF ( calendar_str == 'XXXX' ) THEN
150	WRITE(numout,*) "forcing_info : The calendar was not found in the forcing file."
151	IF (allow_weathergen) THEN
152	! Then change the calendar
153	CALL ioconf_calendar("noleap")
154	ELSE
155	WRITE(numout,*) "forcing_info : We will force it to gregorian by default."
156	CALL ioconf_calendar("gregorian") !! = 365.2425 ; "noleap" = 365.0; "360d"; "julian"=365.25
157	ENDIF
158	ENDIF
159	WRITE(numout,*) "readdim2 : Calendar used by the model : ",calendar_str
160	IF (ttm_full .GE. 2) THEN
161	juld_1 = itau2date(itau(1), date0, dt_force)
162	juld_2 = itau2date(itau(2), date0, dt_force)
163	ELSE
164	juld_1 = 0
165	juld_2 = 0
166	CALL ipslerr_p ( 3, 'forcing_info','What is that only one time step in the forcing file ?', &
167	& ' That can not be right.','verify forcing file.')
168	ENDIF
169	!-
170	!- Initialize one_year / one_day
171	CALL ioget_calendar (one_year, one_day)
172	!-
173	!- What is the distance between the two first states. From this we will deduce what is
174	!- to be done.
175	weathergen = .FALSE.
176	interpol = .FALSE.
177	daily_interpol = .FALSE.
178	is_watchout = .FALSE.
179	!-
180	IF ( ABS(ABS(juld_2-juld_1)-30.) .LE. 2.) THEN
181	IF ( allow_weathergen ) THEN
182	weathergen = .TRUE.
183	WRITE(numout,*) 'Using weather generator.'
184	ELSE
185	CALL ipslerr_p ( 3, 'forcing_info', &
186	& 'This seems to be a monthly file.', &
187	& 'We should use a weather generator with this file.', &
188	& 'This should be allowed in the run.def')
189	ENDIF
190	ELSEIF (( ABS(juld_1-juld_2) .LE. 1./4.) .OR. ( ABS(juld_1-juld_2) .EQ. 1.)) THEN
191	interpol = .TRUE.
192	WRITE(numout,*) 'We will interpolate between the forcing data time steps.'
193	IF ( ABS(juld_1-juld_2) .EQ. 1.) THEN
194	daily_interpol = .TRUE.
195	ENDIF
196	ELSE
197	! Using the weather generator with data other than monthly ones probably
198	! needs some thinking.
199	WRITE(numout,*) 'The time step is not suitable:',ABS(juld_1-juld_2),' days.'
200	CALL ipslerr_p ( 3, 'forcing_info','The time step is not suitable.', &
201	& '','We cannot do anything with these forcing data.')
202	ENDIF
203	!-
204	!- redefine the forcing time step if the weather generator is activated
205	!-
206	IF ( weathergen ) THEN
207	!Config Key = DT_WEATHGEN
208	!Config Desc = Calling frequency of weather generator
209	!Config If = ALLOW_WEATHERGEN
210	!Config Def = 1800.
211	!Config Help = Determines how often the weather generator
212	!Config is called (time step in s). Should be equal
213	!Config to or larger than Sechiba's time step (say,
214	!Config up to 6 times Sechiba's time step or so).
215	!Config Units = [seconds]
216	dt_force = 1800.
217	CALL getin_p('DT_WEATHGEN',dt_force)
218	ENDIF
219	!-
220	!- Define the zoom
221	!-
222	!Config Key = LIMIT_WEST
223	!Config Desc = Western limit of region
224	!Config If = [-]
225	!Config Def = -180.
226	!Config Help = Western limit of the region we are
227	!Config interested in. Between -180 and +180 degrees
228	!Config The model will use the smalest regions from
229	!Config region specified here and the one of the forcing file.
230	!Config Units = [Degrees]
231	!-
232	limit_west = -180.
233	CALL getin_p('LIMIT_WEST',limit_west)
234	!-
235	!Config Key = LIMIT_EAST
236	!Config Desc = Eastern limit of region
237	!Config If = [-]
238	!Config Def = 180.
239	!Config Help = Eastern limit of the region we are
240	!Config interested in. Between -180 and +180 degrees
241	!Config The model will use the smalest regions from
242	!Config region specified here and the one of the forcing file.
243	!Config Units = [Degrees]
244	!-
245	limit_east = 180.
246	CALL getin_p('LIMIT_EAST',limit_east)
247	!-
248	!Config Key = LIMIT_NORTH
249	!Config Desc = Northern limit of region
250	!Config If = [-]
251	!Config Def = 90.
252	!Config Help = Northern limit of the region we are
253	!Config interested in. Between +90 and -90 degrees
254	!Config The model will use the smalest regions from
255	!Config region specified here and the one of the forcing file.
256	!Config Units = [Degrees]
257	!-
258	limit_north = 90.
259	CALL getin_p('LIMIT_NORTH',limit_north)
260	!-
261	!Config Key = LIMIT_SOUTH
262	!Config Desc = Southern limit of region
263	!Config If = [-]
264	!Config Def = -90.
265	!Config Help = Southern limit of the region we are
266	!Config interested in. Between 90 and -90 degrees
267	!Config The model will use the smalest regions from
268	!Config region specified here and the one of the forcing file.
269	!Config Units = [Degrees]
270	!-
271	limit_south = -90.
272	CALL getin_p('LIMIT_SOUTH',limit_south)
273	!-
274	!- Calculate domain size
275	!-
276	IF ( interpol ) THEN
277	!-
278	!- If we use temporal interpolation, then we cannot change the resolution (yet?)
279	!-
280	ALLOCATE(i_index(iim_full), j_index(jjm_full),j_index_g(jjm_full))
281	IF (is_root_prc) THEN
282
283	CALL domain_size (limit_west, limit_east, limit_north, limit_south,&
284	& iim_full, jjm_full, lon_full, lat_full, iim_zoom, jjm_zoom,&
285	& i_index, j_index_g)
286
287	j_index(:)=j_index_g(:)
288
289	ALLOCATE(qair(iim_full,jjm_full))
290	CALL flinget_buffer (force_id,'Qair',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
291	CALL forcing_zoom(data_full, qair)
292
293	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
294	IF ( contfrac_exists ) THEN
295	WRITE(numout,*) "contfrac exist in the forcing file."
296	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
297	CALL forcing_zoom(data_full, fcontfrac)
298	WRITE(numout,*) "fcontfrac min/max :",MINVAL(fcontfrac(1:iim_zoom,1:jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,1:jjm_zoom))
299	ELSE
300	fcontfrac(:,:)=1.
301	ENDIF
302
303
304	DO i=1,iim_zoom
305	DO j=1,jjm_zoom
306	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
307	qair(i,j) = 999999.
308	ENDIF
309	ENDDO
310	ENDDO
311
312	ALLOCATE(index_l(iim_zoom*jjm_zoom))
313	!- In this point is returning the global NbPoint with the global index
314	CALL forcing_landind(iim_zoom,jjm_zoom,qair,.TRUE.,NbPoint,index_l,i_test,j_test)
315	!
316	! Work out the vertical layers to be used
317	!
318	CALL forcing_vertical_ioipsl(force_id)
319	ELSE
320	ALLOCATE(index_l(1))
321	ENDIF
322
323	! Initiate global grid and parallelism
324	CALL bcast(iim_zoom)
325	CALL bcast(jjm_zoom)
326	CALL bcast(NbPoint)
327	CALL grid_set_glo(iim_zoom,jjm_zoom,NbPoint)
328	CALL grid_allocate_glo(4)
329
330	! Check consistency in the number of mpi processors and the number of land points
331	! in order to prevent an exception
332	IF (NbPoint < mpi_size) THEN
333	WRITE(printstr,*) 'The number of landpoints found (', NbPoint, &
334	') is less than the number of processors selected (', mpi_size,')'
335	CALL ipslerr_p(3, 'forcing_info', 'Wrong parallelization options', &
336	TRIM(printstr), &
337	'Decrease the number of processors for the current grid')
338	ENDIF
339
340	!
341	!- global index index_g is the index_l of root proc
342	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
343
344	DEALLOCATE(index_l)
345
346	!
347	! Distribute to all processors the information on the forcing
348	!
349	CALL bcast(index_g)
350	CALL Init_orchidee_data_para_driver(nbp_glo,index_g)
351	CALL init_ioipsl_para
352
353	! Initialize printlev_loc
354	printlev_loc=get_printlev('readdim2')
355	WRITE(numout,*) 'In readdim2, : standard PRINTLEV= ', printlev
356	WRITE(numout,*) 'In readdim2, : local PRINTLEV_readdim2= ', printlev_loc
357
358	! CALL Init_writeField_p
359
360	CALL bcast(jjm_zoom)
361	CALL bcast(i_index)
362	CALL bcast(j_index_g)
363	CALL bcast(zfixed)
364	CALL bcast(zsigma)
365	CALL bcast(zhybrid)
366	CALL bcast(zlevels)
367	CALL bcast(zheight)
368	CALL bcast(zsamelev_uv)
369	CALL bcast(zlev_fixed)
370	CALL bcast(zlevuv_fixed)
371	CALL bcast(zhybrid_a)
372	CALL bcast(zhybrid_b)
373	CALL bcast(zhybriduv_a)
374	CALL bcast(zhybriduv_b)
375	ind=0
376	DO j=1,jjm_zoom
377	IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
378	ind=ind+1
379	j_index(ind)=j_index_g(j)
380	ENDIF
381	ENDDO
382
383	jjm_zoom=jj_nb
384	iim_zoom=iim_g
385
386	!-
387	!- If we use the weather generator, then we read zonal and meridional resolutions
388	!- This should be unified one day...
389	!-
390	ELSEIF ( weathergen ) THEN
391	!-
392	!Config Key = MERID_RES
393	!Config Desc = North-South Resolution
394	!Config Def = 2.
395	!Config If = ALLOW_WEATHERGEN
396	!Config Help = North-South Resolution of the region we are
397	!Config interested in.
398	!Config Units = [Degrees]
399	!-
400	merid_res = 2.
401	CALL getin_p('MERID_RES',merid_res)
402	!-
403	!Config Key = ZONAL_RES
404	!Config Desc = East-West Resolution
405	!Config Def = 2.
406	!Config If = ALLOW_WEATHERGEN
407	!Config Help = East-West Resolution of the region we are
408	!Config interested in. In degrees
409	!Config Units = [Degrees]
410	!-
411	zonal_res = 2.
412	CALL getin_p('ZONAL_RES',zonal_res)
413	!-
414	!- Number of time steps is meaningless in this case
415	!-
416	! ttm_full = HUGE( ttm_full )
417	!MM Number (realistic) of time steps for half hour dt
418	ttm_full = NINT(one_year * 86400. / dt_force)
419	!-
420	IF (is_root_prc) THEN
421
422	!- In this point is returning the global NbPoint with the global index
423	CALL weathgen_domain_size (limit_west,limit_east,limit_north,limit_south, &
424	zonal_res,merid_res,iim_zoom,jjm_zoom)
425	ALLOCATE(index_l(iim_zoom*jjm_zoom))
426	ENDIF
427	CALL bcast(iim_zoom)
428	CALL bcast(jjm_zoom)
429
430	ALLOCATE(lon(iim_zoom,jjm_zoom))
431	ALLOCATE(lat(iim_zoom,jjm_zoom))
432	ALLOCATE(lev(llm_full),levuv(llm_full))
433
434	! We need lon and lat now for weathgen_init
435	CALL forcing_grid (iim_zoom,jjm_zoom,llm_full,lon,lat,init_f=.TRUE.)
436	CALL forcing_vertical_ioipsl(-1)
437
438	IF (is_root_prc) THEN
439	CALL weathgen_init &
440	& (filename,dt_force,force_id,iim_zoom,jjm_zoom, &
441	& zonal_res,merid_res,lon,lat,index_l,NbPoint)
442	!!$,&
443	!!$ & i_index,j_index_g)
444	ELSE
445	ALLOCATE(index_l(1))
446	index_l(1)=1
447	ENDIF
448
449	CALL bcast(NbPoint)
450	CALL grid_set_glo(iim_zoom,jjm_zoom,NbPoint)
451	CALL grid_allocate_glo(4)
452
453	!
454	!- global index index_g is the index_l of root proc
455	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
456
457	DEALLOCATE(index_l)
458
459	CALL bcast(index_g)
460	CALL Init_orchidee_data_para_driver(nbp_glo,index_g)
461	CALL init_ioipsl_para
462	! CALL Init_writeField_p
463	CALL bcast(jjm_zoom)
464	!!$ CALL bcast(i_index)
465	!!$ CALL bcast(j_index_g)
466
467	!!$ ind=0
468	!!$ DO j=1,jjm_zoom
469	!!$ IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
470	!!$ ind=ind+1
471	!!$ j_index(ind)=j_index_g(j)
472	!!$ ENDIF
473	!!$ ENDDO
474
475	jjm_zoom=jj_nb
476	iim_zoom=iim_g
477	!
478	CALL weathgen_read_file(force_id,iim_zoom,jjm_zoom)
479
480	!-
481	ELSE
482	!-
483	CALL ipslerr_p(3,'forcing_info','Neither interpolation nor weather generator is specified.','','')
484	!-
485	ENDIF
486	!-
487	!- Transfer the right information to the caller
488	!-
489	iim = iim_zoom
490	jjm = jjm_zoom
491	llm = 1
492	tm = ttm_full
493	!-
494	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_info : end : ', iim,jjm, llm,tm
495	!-
496	END SUBROUTINE forcing_info
497	!-
498	!-
499	!=====================================================================
500	SUBROUTINE forcing_read &
501	& (filename, rest_id, lrstread, lrstwrite, &
502	& itauin, istp, itau_split, split, nb_spread, netrad_cons, date0, &
503	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
504	& swdown, coszang, precip, snowf, tair, u, v, qair, pb, lwdown, &
505	& fcontfrac, fneighbours, fresolution, &
506	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
507	& kindex, nbindex, force_id)
508	!---------------------------------------------------------------------
509	!- filename : name of the file to be opened
510	!- rest_id : ID of restart file
511	!- lrstread : read restart file?
512	!- lrstwrite : write restart file?
513	!- itauin : time step for which we need the data
514	!- istp : time step for restart file
515	!- itau_split : Where are we within the splitting
516	!- of the time-steps of the forcing files
517	!- (it decides IF we READ or not)
518	!- split : The number of time steps we do
519	!- between two time-steps of the forcing
520	!- nb_spread : Over how many time steps do we spread the precipitation
521	!- netrad_cons: flag that decides if net radiation should be conserved.
522	!- date0 : The date at which the forcing file starts (julian days)
523	!- dt_force : time-step of the forcing file in seconds
524	!- iim : Size of the grid in x
525	!- jjm : size of the grid in y
526	!- lon : Longitudes
527	!- lat : Latitudes
528	!- zlev : First Levels if it exists (ie if watchout file)
529	!- zlevuv : First Levels of the wind (equal precedent, if it exists)
530	!- ttm : number of time steps in all in the forcing file
531	!- swdown : Downward solar radiation (W/m^2)
532	!- coszang : Cosine of the solar zenith angle (unitless)
533	!- precip : Precipitation (Rainfall) (kg/m^2s)
534	!- snowf : Snowfall (kg/m^2s)
535	!- tair : 1st level (2m ? in off-line) air temperature (K)
536	!- u and v : 1st level (2m/10m ? in off-line) (in theory !) wind speed (m/s)
537	!- qair : 1st level (2m ? in off-line) humidity (kg/kg)
538	!- pb : Surface pressure (Pa)
539	!- lwdown : Downward long wave radiation (W/m^2)
540	!- fcontfrac : Continental fraction (no unit)
541	!- fneighbours: land neighbours
542	!- fresolution: resolution in x and y dimensions for each points
543	!-
544	!- From a WATCHOUT file :
545	!- SWnet : Net surface short-wave flux
546	!- Eair : Air potential energy
547	!- petAcoef : Coeficients A from the PBL resolution for T
548	!- peqAcoef : Coeficients A from the PBL resolution for q
549	!- petBcoef : Coeficients B from the PBL resolution for T
550	!- peqBcoef : Coeficients B from the PBL resolution for q
551	!- cdrag : Surface drag
552	!- ccanopy : CO2 concentration in the canopy
553	!-
554	!- kindex : Index of all land-points in the data
555	!- (used for the gathering)
556	!- nbindex : Number of land points
557	!- force_id : FLINCOM file id.
558	!- It is used to close the file at the end of the run.
559	!-
560	!---------------------------------------------------------------------
561	IMPLICIT NONE
562	!-
563	CHARACTER(LEN=*) :: filename
564	INTEGER, INTENT(IN) :: force_id
565	INTEGER, INTENT(INOUT) :: nbindex
566	INTEGER :: rest_id
567	LOGICAL :: lrstread, lrstwrite
568	INTEGER :: itauin, istp, itau_split, split, nb_spread
569	LOGICAL :: netrad_cons
570	REAL :: date0, dt_force
571	INTEGER :: iim, jjm, ttm
572	REAL,DIMENSION(iim,jjm) :: lon, lat, zlev, zlevuv, &
573	& swdown, coszang, precip, snowf, tair, u, v, qair, pb, lwdown, &
574	& fcontfrac
575	REAL,DIMENSION(iim,jjm,2) :: fresolution
576	INTEGER,DIMENSION(iim,jjm,8) :: fneighbours
577	! for watchout files
578	REAL,DIMENSION(iim,jjm) :: &
579	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
580	INTEGER,DIMENSION(iim*jjm), INTENT(INOUT) :: kindex
581	!-
582	INTEGER :: ik,i,j
583	!
584	IF ( interpol ) THEN
585	!-
586	CALL forcing_read_interpol &
587	(filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, &
588	dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
589	swdown, coszang, precip, snowf, tair, u, v, qair, pb, lwdown, &
590	fcontfrac, fneighbours, fresolution, &
591	SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
592	kindex, nbindex, force_id)
593	!-
594	ELSEIF ( weathergen ) THEN
595	!-
596	IF (lrstread) THEN
597	fcontfrac(:,:) = 1.0
598	WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
599	ENDIF
600
601	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
602	CALL weathgen_main (istp, istp, filename, force_id, iim, jjm, 1, &
603	rest_id, lrstread, lrstwrite, &
604	limit_west, limit_east, limit_north, limit_south, &
605	zonal_res, merid_res, lon, lat, date0, dt_force, &
606	kindex, nbindex, &
607	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
608	ELSE
609	CALL weathgen_main (itauin, istp, filename, force_id, iim, jjm, 1, &
610	rest_id, lrstread, lrstwrite, &
611	limit_west, limit_east, limit_north, limit_south, &
612	zonal_res, merid_res, lon, lat, date0, dt_force, &
613	kindex, nbindex, &
614	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
615	ENDIF
616	!-
617	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
618	!---
619	!--- Allocate grid stuff
620	!---
621	CALL grid_init ( nbindex, 4, "RegLonLat", "ForcingGrid" )
622	!---
623	!--- Compute
624	!---
625	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, index_g)
626	!CALL grid_stuff (nbindex, iim, jjm, lon, lat, kindex)
627	DO ik=1,nbindex
628
629	j = ((kindex(ik)-1)/iim) + 1
630	i = (kindex(ik) - (j-1)*iim)
631	!-
632	!- Store variable to help describe the grid
633	!- once the points are gathered.
634	!-
635	fneighbours(i,j,:) = neighbours(ik,:)
636	!
637	fresolution(i,j,:) = resolution(ik,:)
638	ENDDO
639	ENDIF
640	ELSE
641	!-
642
643	CALL ipslerr_p(3,'forcing_read','Neither interpolation nor weather generator is specified.','','')
644
645	ENDIF
646	!-
647	IF (.NOT. is_watchout) THEN
648	! We have to compute Potential air energy
649	WHERE(tair(:,:) < val_exp)
650	eair(:,:) = cp_airtair(:,:)+cte_gravzlev(:,:)
651	ENDWHERE
652	ENDIF
653	!-
654	!
655	!-------------------------
656	END SUBROUTINE forcing_read
657	!=====================================================================
658	!-
659	!-
660	!=====================================================================
661	SUBROUTINE forcing_read_interpol &
662	& (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, &
663	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, swdown, coszang, rainf, snowf, tair, &
664	& u, v, qair, pb, lwdown, &
665	& fcontfrac, fneighbours, fresolution, &
666	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
667	& kindex, nbindex, force_id)
668	!---------------------------------------------------------------------
669	!- filename : name of the file to be opened
670	!- itauin : time step for which we need the data
671	!- itau_split : Where are we within the splitting
672	!- of the time-steps of the forcing files
673	!- (it decides IF we READ or not)
674	!- split : The number of time steps we do
675	!- between two time-steps of the forcing
676	!- nb_spread : Over how many time steps do we spread the precipitation
677	!- netrad_cons: flag that decides if net radiation should be conserved.
678	!- date0 : The date at which the forcing file starts (julian days)
679	!- dt_force : time-step of the forcing file in seconds
680	!- iim : Size of the grid in x
681	!- jjm : size of the grid in y
682	!- lon : Longitudes
683	!- lat : Latitudes
684	!- zlev : First Levels if it exists (ie if watchout file)
685	!- zlevuv : First Levels of the wind (equal precedent, if it exists)
686	!- ttm : number of time steps in all in the forcing file
687	!- swdown : Downward solar radiation (W/m^2)
688	!- coszang : Cosine of the solar zenith angle (unitless)
689	!- rainf : Rainfall (kg/m^2s)
690	!- snowf : Snowfall (kg/m^2s)
691	!- tair : 2m air temperature (K)
692	!- u and v : 2m (in theory !) wind speed (m/s)
693	!- qair : 2m humidity (kg/kg)
694	!- pb : Surface pressure (Pa)
695	!- lwdown : Downward long wave radiation (W/m^2)
696	!- fcontfrac : Continental fraction (no unit)
697	!- fneighbours: land neighbours
698	!- fresolution: resolution in x and y dimensions for each points
699	!-
700	!- From a WATCHOUT file :
701	!- SWnet : Net surface short-wave flux
702	!- Eair : Air potential energy
703	!- petAcoef : Coeficients A from the PBL resolution for T
704	!- peqAcoef : Coeficients A from the PBL resolution for q
705	!- petBcoef : Coeficients B from the PBL resolution for T
706	!- peqBcoef : Coeficients B from the PBL resolution for q
707	!- cdrag : Surface drag
708	!- ccanopy : CO2 concentration in the canopy
709	!-
710	!- kindex : Index of all land-points in the data
711	!- (used for the gathering)
712	!- nbindex : Number of land points
713	!- force_id : FLINCOM file id.
714	!- It is used to close the file at the end of the run.
715	!---------------------------------------------------------------------
716	IMPLICIT NONE
717	!-
718	INTEGER,PARAMETER :: lm=1
719	!-
720	!- Input variables
721	!-
722	CHARACTER(LEN=*) :: filename
723	INTEGER :: itauin, itau_split, split, nb_spread
724	LOGICAL :: netrad_cons
725	REAL :: date0, dt_force
726	INTEGER :: iim, jjm, ttm
727	REAL,DIMENSION(:,:),INTENT(IN) :: lon, lat !- LOCAL data array (size=iim,jjm)
728	INTEGER, INTENT(IN) :: force_id
729	!-
730	!- Output variables
731	!-
732	REAL,DIMENSION(:,:),INTENT(OUT) :: zlev, zlevuv, & !- LOCAL data array (size=iim,jjm)
733	& swdown, coszang, rainf, snowf, tair, u, v, qair, pb, lwdown, &
734	& fcontfrac
735	REAL,DIMENSION(:,:,:),INTENT(OUT) :: fresolution !- LOCAL data array (size=iim,jjm,2)
736	INTEGER,DIMENSION(:,:,:),INTENT(OUT) :: fneighbours !- LOCAL data array (size=iim,jjm,8)
737	! for watchout files
738	REAL,DIMENSION(:,:),INTENT(OUT) :: &
739	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
740	INTEGER,DIMENSION(:),INTENT(INOUT) :: kindex !- LOCAL index of the map
741	INTEGER, INTENT(INOUT) :: nbindex
742	!-
743	!- Local variables
744	!-
745	INTEGER, SAVE :: last_read=0
746	INTEGER, SAVE :: itau_read, itau_read_nm1=0, itau_read_n=0
747	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
748	& zlev_nm1, zlevuv_nm1, swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, u_nm1, v_nm1, qair_nm1, &
749	& pb_nm1, lwdown_nm1, &
750	& zlev_n, zlevuv_n, swdown_n, rainf_n, snowf_n, tair_n, u_n, v_n, qair_n, &
751	& pb_n, lwdown_n, mean_coszang
752
753	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
754	& startday_n, startday_nm1, daylength_n, daylength_nm1, tmax_n, tmax_nm1, tmin_nm1, tmin_nm2, tmin_n, &
755	& qsatta, qsattmin_n, qsattmin_nm1, qmin_n, qmin_nm1, qmax_n, qmax_nm1, qsa
756	REAL,SAVE :: hour
757
758	! just for grid stuff if the forcing file is a watchout file
759	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmpdata
760	! variables to be read in watchout files
761	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
762	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
763	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n
764	REAL, SAVE :: julian_for ! Date of the forcing to be read
765	REAL :: julian, ss, rw
766	!jur,
767	REAL, SAVE :: julian0 ! First day of this year
768	INTEGER :: yy, mm, dd, is, i, j, ik
769	REAL(r_std), DIMENSION(2) :: min_resol, max_resol
770	! if Wind_N and Wind_E are in the file (and not just Wind)
771	LOGICAL, SAVE :: wind_N_exists=.FALSE.
772	LOGICAL :: wind_E_exists=.FALSE.
773	LOGICAL, SAVE :: contfrac_exists=.FALSE.
774	LOGICAL, SAVE :: neighbours_exists=.FALSE.
775	LOGICAL, SAVE :: initialized = .FALSE.
776	LOGICAL :: check=.FALSE.
777	! to bypass FPE problem on integer convertion with missing_value heigher than precision
778	INTEGER, PARAMETER :: undef_int = 999999999
779	!---------------------------------------------------------------------
780
781	itau_read = MOD((itauin-1),ttm)+1
782
783	IF (check) THEN
784	WRITE(numout,*) &
785	" FORCING READ : itauin, itau_read, itau_split : ",&
786	itauin, itau_read, itau_split
787	ENDIF
788
789	!-
790	!- This part initializes the reading of the forcing. As you can see
791	!- we only go through here if both time steps are zero.
792	!-
793	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
794	!-
795	!- Tests on forcing file type
796	CALL flinquery_var(force_id, 'Wind_N', wind_N_exists)
797	IF ( wind_N_exists ) THEN
798	CALL flinquery_var(force_id, 'Wind_E', wind_E_exists)
799	IF ( .NOT. wind_E_exists ) THEN
800	CALL ipslerr_p(3,'forcing_read_interpol', &
801	& 'Variable Wind_E does not exist in forcing file', &
802	& 'But variable Wind_N exists.','Please, rename Wind_N to Wind;')
803	ENDIF
804	ENDIF
805	CALL flinquery_var(force_id, 'levels', is_watchout)
806	IF ( is_watchout ) THEN
807	WRITE(numout,*) "Read a Watchout File."
808	ENDIF
809	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
810	!-
811	IF (check) WRITE(numout,*) 'ALLOCATE all the memory needed'
812	!-
813	ALLOCATE &
814	& (swdown_nm1(iim,jjm), rainf_nm1(iim,jjm), snowf_nm1(iim,jjm), &
815	& tair_nm1(iim,jjm), u_nm1(iim,jjm), v_nm1(iim,jjm), qair_nm1(iim,jjm), &
816	& pb_nm1(iim,jjm), lwdown_nm1(iim,jjm))
817	ALLOCATE &
818	& (swdown_n(iim,jjm), rainf_n(iim,jjm), snowf_n(iim,jjm), &
819	& tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), &
820	& pb_n(iim,jjm), lwdown_n(iim,jjm))
821
822	IF(daily_interpol) THEN
823	ALLOCATE &
824	& (startday_n(iim,jjm), startday_nm1(iim,jjm), daylength_n(iim,jjm), &
825	& daylength_nm1(iim,jjm), tmax_n(iim,jjm), tmax_nm1(iim,jjm), tmin_n(iim,jjm), &
826	& tmin_nm1(iim,jjm), tmin_nm2(iim,jjm), qsatta(iim,jjm), qsattmin_n(iim,jjm), qsattmin_nm1(iim,jjm), &
827	& qmin_n(iim,jjm), qmin_nm1(iim,jjm), qmax_n(iim,jjm), qmax_nm1(iim,jjm), qsa(iim,jjm) )
828	ENDIF
829
830
831	ALLOCATE &
832	& (zlev_nm1(iim,jjm), zlev_n(iim,jjm), zlevuv_nm1(iim,jjm), zlevuv_n(iim,jjm), &
833	& SWnet_nm1(iim,jjm), Eair_nm1(iim,jjm), cdrag_nm1(iim,jjm), ccanopy_nm1(iim,jjm), &
834	& petAcoef_nm1(iim,jjm), peqAcoef_nm1(iim,jjm), petBcoef_nm1(iim,jjm), peqBcoef_nm1(iim,jjm), &
835	& SWnet_n(iim,jjm), Eair_n(iim,jjm), cdrag_n(iim,jjm), ccanopy_n(iim,jjm), &
836	& petAcoef_n(iim,jjm), peqAcoef_n(iim,jjm), petBcoef_n(iim,jjm), peqBcoef_n(iim,jjm))
837	ALLOCATE &
838	& (mean_coszang(iim,jjm))
839	!-
840	IF (check) WRITE(numout,*) 'Memory ALLOCATED'
841	!-
842	!- We need for the driver surface air temperature and humidity before the
843	!- the loop starts. Thus we read it here.
844	!-
845	CALL forcing_just_read (iim, jjm, zlev, zlevuv, ttm, 1, 1, &
846	& swdown, rainf, snowf, tair, &
847	& u, v, qair, pb, lwdown, &
848	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
849	& force_id, wind_N_exists, check)
850	!----
851
852	!-- First in case it's not a watchout file
853	IF ( .NOT. is_watchout ) THEN
854	IF ( contfrac_exists ) THEN
855	WRITE(numout,*) "contfrac exist in the forcing file."
856	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
857	CALL forcing_zoom(data_full, fcontfrac)
858	WRITE(numout,*) "fcontfrac min/max :",MINVAL(fcontfrac(1:iim_zoom,jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,jjm_zoom))
859	!
860	! We need to make sure that when we gather the points we pick all
861	! the points where contfrac is above 0. Thus we prepare tair for
862	! subroutine forcing_landind
863	!
864	DO i=1,iim
865	DO j=1,jjm
866	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0. ! bande de recouvrement du scatter2D
867	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0. ! => on mets les données qu'on veut pas du noeud à missing_value
868	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
869	tair(i,j) = 999999.
870	ENDIF
871	ENDDO
872	ENDDO
873	ELSE
874	fcontfrac(:,:) = 1.0
875	ENDIF
876	!---
877	!--- Create the index table
878	!---
879	!--- This job return a LOCAL kindex
880	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
881	#ifdef CPP_PARA
882	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
883	! Force nbindex points for parallel computation
884	nbindex=nbp_loc
885	CALL scatter(index_g,kindex(1:nbindex))
886	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
887	#endif
888	ik=MAX(nbindex/2,1)
889	j_test = (((kindex(ik)-1)/iim) + 1)
890	i_test = (kindex(ik) - (j_test-1)*iim)
891	IF (check) THEN
892	WRITE(numout,*) 'New test point is : ', i_test, j_test
893	ENDIF
894	!---
895	!--- Allocate grid stuff
896	!---
897	CALL grid_init ( nbindex, 4, "RegLonLat", "ForcingGrid" )
898	!---
899	!--- All grid variables
900	!---
901	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, index_g)
902	DO ik=1,nbindex
903	!
904	j = ((kindex(ik)-1)/iim) + 1
905	i = (kindex(ik) - (j-1)*iim)
906	!-
907	!- Store variable to help describe the grid
908	!- once the points are gathered.
909	!-
910	fneighbours(i,j,:) = neighbours(ik,:)
911	!
912	fresolution(i,j,:) = resolution(ik,:)
913	ENDDO
914	ELSE
915	!-- Second, in case it is a watchout file
916	ALLOCATE (tmpdata(iim,jjm))
917	tmpdata(:,:) = 0.0
918	!--
919	IF ( .NOT. contfrac_exists ) THEN
920	CALL ipslerr_p (3,'forcing_read_interpol', &
921	& 'Could get contfrac variable in a watchout file :',TRIM(filename), &
922	& '(Problem with file ?)')
923	ENDIF
924	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
925	CALL forcing_zoom(data_full, fcontfrac)
926	!
927	! We need to make sure that when we gather the points we pick all
928	! the points where contfrac is above 0. Thus we prepare tair for
929	! subroutine forcing_landind
930	!
931	DO i=1,iim
932	DO j=1,jjm
933	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0.
934	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0.
935	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
936	tair(i,j) = 999999.
937	ENDIF
938	ENDDO
939	ENDDO
940	!---
941	!--- Create the index table
942	!---
943	!--- This job return a LOCAL kindex
944	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
945	#ifdef CPP_PARA
946	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
947	! Force nbindex points for parallel computation
948	nbindex=nbp_loc
949	CALL scatter(index_g,kindex)
950	kindex(:)=kindex(:)-offset
951	! kindex(:)=kindex(:)-(jj_begin-1)*iim_g
952	#endif
953	ik=MAX(nbindex/2,1)
954	j_test = (((kindex(ik)-1)/iim) + 1)
955	i_test = (kindex(ik) - (j_test-1)*iim)
956	IF (check) THEN
957	WRITE(numout,*) 'New test point is : ', i_test, j_test
958	ENDIF
959	!---
960	!--- Allocate grid stuff
961	!---
962	CALL grid_init ( nbindex, 4, "RegLonLat", "ForcingGrid" )
963	neighbours(:,:) = -1
964	resolution(:,:) = 0.
965	min_resol(:) = 1.e6
966	max_resol(:) = -1.
967	!---
968	!--- All grid variables
969	!---
970	!-
971	!- Get variables to help describe the grid
972	CALL flinquery_var(force_id, 'neighboursNN', neighbours_exists)
973	IF ( .NOT. neighbours_exists ) THEN
974	CALL ipslerr_p (3,'forcing_read_interpol', &
975	& 'Could get neighbours in a watchout file :',TRIM(filename), &
976	& '(Problem with file ?)')
977	ELSE
978	WRITE(numout,*) "Watchout file contains neighbours and resolutions."
979	ENDIF
980	!
981	fneighbours(:,:,:) = undef_int
982	!
983	!- once the points are gathered.
984	CALL flinget_buffer (force_id,'neighboursNN',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
985	CALL forcing_zoom(data_full, tmpdata)
986	WHERE(tmpdata(:,:) < undef_int)
987	fneighbours(:,:,1) = NINT(tmpdata(:,:))
988	ENDWHERE
989	!
990	CALL flinget_buffer (force_id,'neighboursNE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
991	CALL forcing_zoom(data_full, tmpdata)
992	WHERE(tmpdata(:,:) < undef_int)
993	fneighbours(:,:,2) = NINT(tmpdata(:,:))
994	ENDWHERE
995	!
996	CALL flinget_buffer (force_id,'neighboursEE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
997	CALL forcing_zoom(data_full, tmpdata)
998	WHERE(tmpdata(:,:) < undef_int)
999	fneighbours(:,:,3) = NINT(tmpdata(:,:))
1000	ENDWHERE
1001	!
1002	CALL flinget_buffer (force_id,'neighboursSE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1003	CALL forcing_zoom(data_full, tmpdata)
1004	WHERE(tmpdata(:,:) < undef_int)
1005	fneighbours(:,:,4) = NINT(tmpdata(:,:))
1006	ENDWHERE
1007	!
1008	CALL flinget_buffer (force_id,'neighboursSS',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1009	CALL forcing_zoom(data_full, tmpdata)
1010	WHERE(tmpdata(:,:) < undef_int)
1011	fneighbours(:,:,5) = NINT(tmpdata(:,:))
1012	ENDWHERE
1013	!
1014	CALL flinget_buffer (force_id,'neighboursSW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1015	CALL forcing_zoom(data_full, tmpdata)
1016	WHERE(tmpdata(:,:) < undef_int)
1017	fneighbours(:,:,6) = NINT(tmpdata(:,:))
1018	ENDWHERE
1019	!
1020	CALL flinget_buffer (force_id,'neighboursWW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1021	CALL forcing_zoom(data_full, tmpdata)
1022	WHERE(tmpdata(:,:) < undef_int)
1023	fneighbours(:,:,7) = NINT(tmpdata(:,:))
1024	ENDWHERE
1025	!
1026	CALL flinget_buffer (force_id,'neighboursNW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1027	CALL forcing_zoom(data_full, tmpdata)
1028	WHERE(tmpdata(:,:) < undef_int)
1029	fneighbours(:,:,8) = NINT(tmpdata(:,:))
1030	ENDWHERE
1031	!
1032	! now, resolution of the grid
1033	CALL flinget_buffer (force_id,'resolutionX',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1034	CALL forcing_zoom(data_full, tmpdata)
1035	fresolution(:,:,1) = tmpdata(:,:)
1036	!
1037	CALL flinget_buffer (force_id,'resolutionY',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1038	CALL forcing_zoom(data_full, tmpdata)
1039	fresolution(:,:,2) = tmpdata(:,:)
1040	!
1041	DO ik=1,nbindex
1042	!
1043	j = ((kindex(ik)-1)/iim) + 1
1044	i = (kindex(ik) - (j-1)*iim)
1045	!-
1046	!- Store variable to help describe the grid
1047	!- once the points are gathered.
1048	!-
1049	neighbours(ik,:) = fneighbours(i,j,:)
1050	!
1051	resolution(ik,:) = fresolution(i,j,:)
1052	!
1053
1054	ENDDO
1055	CALL gather(neighbours,neighbours_g)
1056	CALL gather(resolution,resolution_g)
1057	min_resol(1) = MINVAL(resolution(:,1))
1058	min_resol(2) = MAXVAL(resolution(:,2))
1059	max_resol(1) = MAXVAL(resolution(:,1))
1060	max_resol(2) = MAXVAL(resolution(:,2))
1061	!
1062	area(:) = resolution(:,1)*resolution(:,2)
1063	CALL gather(area,area_g)
1064	!--
1065	DEALLOCATE (tmpdata)
1066	ENDIF
1067	WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
1068	!---
1069	ENDIF
1070	!---
1071	IF (check) THEN
1072	WRITE(numout,*) &
1073	& 'The dates : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1074	ENDIF
1075	!---
1076	!--- Here we do the work in case only interpolation is needed.
1077	!---
1078	IF ( initialized .AND. interpol ) THEN
1079	!---
1080	IF ( daily_interpol ) THEN
1081
1082	IF (split > 1) THEN
1083	IF ( itau_split <= (split/2.) ) THEN
1084	rw = REAL(itau_split+split/2.)/split
1085	ELSE
1086	rw = REAL(itau_split-split/2.)/split
1087	ENDIF
1088	ELSE
1089	rw = 1.
1090	ENDIF
1091
1092	IF ((last_read == 0) .OR. ( rw==(1./split)) ) THEN
1093	!---
1094	!----- Start or Restart
1095	IF (last_read == 0) THEN
1096	! Case of a restart or a shift in the forcing file.
1097	IF (itau_read > 1) THEN
1098	itau_read_nm1=itau_read-1
1099	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1100	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1101	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1102	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1103	& force_id, wind_N_exists, check)
1104	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1105	! Case of a simple start.
1106	ELSE
1107	itau_read_nm1 = un
1108	WRITE(numout,*) "we will use the forcing of the first day to initialize "
1109	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1110	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1111	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1112	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1113	& force_id, wind_N_exists, check)
1114	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1115	ENDIF
1116	tmin_nm2(:,:)=tmin_nm1(:,:)
1117	IF ( dt_force .GT. 3600. ) THEN
1118	mean_coszang(:,:) = 0.0
1119	daylength_n(:,:) = 0.
1120	DO is=1,split
1121	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1122	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1123	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1124	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1125	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1126	WHERE( coszang(:,:) > 0. )
1127	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1128	ENDWHERE
1129	ENDDO
1130	mean_coszang(:,:) = mean_coszang(:,:)/split
1131	daylength_nm1(:,:)=daylength_n(:,:)
1132	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1133	ENDIF
1134	ELSE
1135	!----- Normal mode : copy old step
1136	swdown_nm1(:,:) = swdown_n(:,:)
1137	rainf_nm1(:,:) = rainf_n(:,:)
1138	snowf_nm1(:,:) = snowf_n(:,:)
1139	tair_nm1(:,:) = tair_n(:,:)
1140	u_nm1(:,:) = u_n(:,:)
1141	v_nm1(:,:) = v_n(:,:)
1142	qair_nm1(:,:) = qair_n(:,:)
1143	pb_nm1(:,:) = pb_n(:,:)
1144	lwdown_nm1(:,:) = lwdown_n(:,:)
1145	tmin_nm2(:,:) = tmin_nm1(:,:)
1146	tmin_nm1(:,:) = tmin_n(:,:)
1147	tmax_nm1(:,:) = tmax_n(:,:)
1148
1149	IF (is_watchout) THEN
1150	zlev_nm1(:,:) = zlev_n(:,:)
1151	zlevuv_nm1(:,:) = zlevuv_n(:,:)
1152	! Net surface short-wave flux
1153	SWnet_nm1(:,:) = SWnet_n(:,:)
1154	! Air potential energy
1155	Eair_nm1(:,:) = Eair_n(:,:)
1156	! Coeficients A from the PBL resolution for T
1157	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1158	! Coeficients A from the PBL resolution for q
1159	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1160	! Coeficients B from the PBL resolution for T
1161	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1162	! Coeficients B from the PBL resolution for q
1163	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1164	! Surface drag
1165	cdrag_nm1(:,:) = cdrag_n(:,:)
1166	! CO2 concentration in the canopy
1167	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1168	ENDIF
1169	itau_read_nm1 = itau_read_n
1170	ENDIF
1171	!-----
1172	!-----
1173	IF(last_read==0)THEN
1174	itau_read_n = itau_read
1175	ELSE
1176	itau_read_n = itau_read+1
1177	ENDIF
1178
1179	IF (itau_read_n > ttm) THEN
1180	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1181	WRITE(numout,*) &
1182	& 'WARNING : We are going back to the start of the file'
1183	itau_read_n =1
1184	ENDIF
1185	IF (check) THEN
1186	WRITE(numout,*) &
1187	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1188	ENDIF
1189	!-----
1190	!----- Get a reduced julian day !
1191	!----- This is needed because we lack the precision on 32 bit machines.
1192	!-----
1193	IF ( dt_force .GT. 3600. ) THEN
1194	julian_for = itau2date(itau_read-1, date0, dt_force)
1195	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1196
1197	! first day of this year
1198	CALL ymds2ju (yy,1,1,0.0, julian0)
1199	!-----
1200	IF (check) THEN
1201	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1202	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1203	ENDIF
1204	ENDIF
1205	!-----
1206	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1207	& swdown_n, rainf_n, snowf_n, tmin_n, &
1208	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1209	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1210	& force_id, wind_N_exists, check)
1211	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_n, itau_read_n, tmax_n, force_id )
1212
1213	!---
1214	last_read = itau_read_n
1215	!-----
1216	!----- Compute mean solar angle for the comming period
1217	!-----
1218	IF (check) WRITE(numout,*) 'Going into solarang', split, one_day
1219	!-----
1220
1221	!-----
1222	ENDIF
1223	!---
1224	IF ( itau_split == 1. ) THEN
1225	IF ( dt_force .GT. 3600. ) THEN
1226	mean_coszang(:,:) = 0.0
1227	daylength_nm1(:,:)=daylength_n(:,:)
1228	daylength_n(:,:) = 0.
1229	DO is=1,split
1230	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1231	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1232	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1233	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1234	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1235	WHERE( coszang(:,:) > 0. )
1236	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1237	ENDWHERE
1238	ENDDO
1239	mean_coszang(:,:) = mean_coszang(:,:)/split
1240	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1241	ENDIF
1242	ENDIF
1243
1244	!--- Do the interpolation
1245	IF (check) WRITE(numout,*) 'Doing the interpolation between time steps'
1246	!---
1247
1248	IF (check) WRITE(numout,*) 'Coeff of interpollation : ',rw
1249	!---
1250
1251	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1252	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1253	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1254
1255	!--- Take care of the height of the vertical levels
1256	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1257	zlevuv(:,:) = (zlevuv_n(:,:)-zlevuv_nm1(:,:))*rw + zlevuv_nm1(:,:)
1258
1259	hour=REAL(itau_split)/split*24
1260	startday_n(:,:)=12.-daylength_n(:,:)/2.
1261	startday_nm1(:,:)=12.-daylength_nm1(:,:)/2.
1262
1263	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1264	tair(:,:)=(tmax_nm1(:,:)-tmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1265	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_nm1(:,:)+tmin_nm1(:,:))/2.
1266	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1267	tair(:,:)=(tmax_n(:,:)-tmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1268	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_n(:,:)+tmin_n(:,:))/2.
1269	ELSEWHERE ( hour < startday_n(:,:) )
1270	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1271	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1272	ELSEWHERE
1273	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1274	& (24.-14.+startday_n(:,:))-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1275	ENDWHERE
1276
1277	CALL weathgen_qsat_2d (iim,jjm,tmin_n,pb,qsattmin_n)
1278	CALL weathgen_qsat_2d (iim,jjm,tmin_nm1,pb,qsattmin_nm1)
1279	CALL weathgen_qsat_2d (iim,jjm,tair,pb,qsatta)
1280
1281	!---
1282	qmin_nm1(:,:) = MIN(qair_nm1(:,:),0.99*qsattmin_nm1(:,:))
1283	qmin_n(:,:) = MIN(qair_n(:,:),0.99*qsattmin_n(:,:))
1284	qmax_nm1(:,:) = (qair_nm1(:,:)-qmin_nm1(:,:)) + qair_nm1(:,:)
1285	qmax_n(:,:) = (qair_n(:,:)-qmin_n(:,:)) + qair_n(:,:)
1286
1287	qsa(:,:) = 0.99*qsatta(:,:)
1288
1289
1290	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1291	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1292	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_nm1(:,:)+qmin_nm1(:,:))/2.)
1293	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1294	qair(:,:)=MIN(qsa(:,:),(qmax_n(:,:)-qmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1295	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_n(:,:)+qmin_n(:,:))/2.)
1296	ELSEWHERE ( hour < startday_n(:,:) )
1297	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1298	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1299	ELSEWHERE
1300	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1301	& (24.-14.+startday_n(:,:))-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1302	ENDWHERE
1303
1304	IF (is_watchout) THEN
1305	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1306	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1307	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1308	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1309	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1310	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1311	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1312	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1313	ENDIF
1314	!---
1315	!--- Here we need to allow for an option
1316	!--- where radiative energy is conserved
1317	!---
1318	IF ( netrad_cons ) THEN
1319	lwdown(:,:) = lwdown_n(:,:)
1320	ELSE
1321	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1322	ENDIF
1323	!---
1324	!--- For the solar radiation we decompose the mean value
1325	!--- using the zenith angle of the sun if the time step in the forcing data is
1326	!---- more than an hour. Else we use the standard linera interpolation
1327	!----
1328	IF (check) WRITE(numout,*) 'Ready to deal with the solar radiation'
1329	!----
1330	IF ( dt_force .GT. 3600. ) THEN
1331	!---
1332	IF ( netrad_cons ) THEN
1333	WRITE(numout,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
1334	ENDIF
1335	!---
1336	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1337	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1338	!!$ julian = julian_for + rw*dt_force/one_day
1339	IF (check) THEN
1340	WRITE(numout,'(a,f20.10,2I3)') &
1341	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1342	ENDIF
1343	!---
1344	CALL solarang(julian, julian0, iim, jjm, lon*0, lat, coszang)
1345	!---
1346
1347	WHERE ((mean_coszang(:,:) > 0.) .AND. (hour <= 12 ))
1348	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1349	ELSEWHERE ((mean_coszang(:,:) > 0.) .AND. (hour > 12 ))
1350	swdown(:,:) = swdown_nm1(:,:) *coszang(:,:)/mean_coszang(:,:)
1351	ELSEWHERE
1352	swdown(:,:) = 0.0
1353	END WHERE
1354	!---
1355	WHERE (swdown(:,:) > 2000. )
1356	swdown(:,:) = 2000.
1357	END WHERE
1358	!---
1359	ELSE
1360	!---
1361	IF ( netrad_cons ) THEN
1362	swdown(:,:) = swdown_n(:,:)
1363	ELSE
1364	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1365	ENDIF
1366	!---
1367	ENDIF
1368	!---
1369	IF (check) THEN
1370	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1371	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1372	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1373	IF (is_watchout) THEN
1374	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1375	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1376	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1377	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1378	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1379	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1380	ENDIF
1381	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1382	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1383	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1384	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1385	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1386	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1387	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1388	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1389	ENDIF
1390	!---
1391	!--- For precip we suppose that the rain
1392	!--- is the sum over the next 6 hours
1393	!---
1394	WHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)>=273.15))
1395	rainf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1396	snowf(:,:) = 0.0
1397	ELSEWHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)<273.15))
1398	snowf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1399	rainf(:,:) = 0.0
1400	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)>=273.15))
1401	rainf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1402	snowf(:,:) = 0.0
1403	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)<273.15))
1404	snowf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1405	rainf(:,:) = 0.0
1406	ELSEWHERE
1407	snowf(:,:) = 0.0
1408	rainf(:,:) = 0.0
1409	ENDWHERE
1410
1411	IF (check) THEN
1412	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1413	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1414	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1415	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1416	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1417	ENDIF
1418	!---
1419
1420
1421	ELSE
1422
1423	IF (itau_read /= last_read) THEN
1424	!---
1425	!----- Start or Restart
1426	IF (itau_read_n == 0) THEN
1427	! Case of a restart or a shift in the forcing file.
1428	IF (itau_read > 1) THEN
1429	itau_read_nm1=itau_read-1
1430	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1431	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1432	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1433	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1434	& force_id, wind_N_exists, check)
1435	! Case of a simple start.
1436	ELSE IF (dt_force*ttm > one_day-1. ) THEN
1437	! if the forcing file contains at least 24 hours,
1438	! we will use the last forcing step of the first day
1439	! as initiale condition to prevent first shift off reading.
1440	itau_read_nm1 = NINT (one_day/dt_force)
1441	WRITE(numout,) "the forcing file contains 24 hours :",dt_forcettm,one_day-1.
1442	WRITE(numout,*) "we will use the last forcing step of the first day : itau_read_nm1 ",itau_read_nm1
1443	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1444	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1445	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1446	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1447	& force_id, wind_N_exists, check)
1448	ELSE
1449	! if the forcing file contains less than 24 hours,
1450	! just say error !
1451	CALL ipslerr_p(3,'forcing_read_interpol', &
1452	& 'The forcing file contains less than 24 hours !', &
1453	& 'We can''t intialize interpolation with such a file.','')
1454	ENDIF
1455	ELSE
1456	!----- Normal mode : copy old step
1457	swdown_nm1(:,:) = swdown_n(:,:)
1458	rainf_nm1(:,:) = rainf_n(:,:)
1459	snowf_nm1(:,:) = snowf_n(:,:)
1460	tair_nm1(:,:) = tair_n(:,:)
1461	u_nm1(:,:) = u_n(:,:)
1462	v_nm1(:,:) = v_n(:,:)
1463	qair_nm1(:,:) = qair_n(:,:)
1464	pb_nm1(:,:) = pb_n(:,:)
1465	lwdown_nm1(:,:) = lwdown_n(:,:)
1466	IF (is_watchout) THEN
1467	zlev_nm1(:,:) = zlev_n(:,:)
1468	! Net surface short-wave flux
1469	SWnet_nm1(:,:) = SWnet_n(:,:)
1470	! Air potential energy
1471	Eair_nm1(:,:) = Eair_n(:,:)
1472	! Coeficients A from the PBL resolution for T
1473	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1474	! Coeficients A from the PBL resolution for q
1475	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1476	! Coeficients B from the PBL resolution for T
1477	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1478	! Coeficients B from the PBL resolution for q
1479	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1480	! Surface drag
1481	cdrag_nm1(:,:) = cdrag_n(:,:)
1482	! CO2 concentration in the canopy
1483	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1484	ENDIF
1485	itau_read_nm1 = itau_read_n
1486	ENDIF
1487	!-----
1488	itau_read_n = itau_read
1489	IF (itau_read_n > ttm) THEN
1490	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1491	WRITE(numout,*) &
1492	& 'WARNING : We are going back to the start of the file'
1493	itau_read_n =1
1494	ENDIF
1495	IF (check) THEN
1496	WRITE(numout,*) &
1497	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1498	ENDIF
1499	!-----
1500	!----- Get a reduced julian day !
1501	!----- This is needed because we lack the precision on 32 bit machines.
1502	!-----
1503	IF ( dt_force .GT. 3600. ) THEN
1504	julian_for = itau2date(itau_read-1, date0, dt_force)
1505	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1506
1507	! first day of this year
1508	CALL ymds2ju (yy,1,1,0.0, julian0)
1509	!-----
1510	IF (check) THEN
1511	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1512	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1513	ENDIF
1514	ENDIF
1515	!-----
1516	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1517	& swdown_n, rainf_n, snowf_n, tair_n, &
1518	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1519	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1520	& force_id, wind_N_exists, check)
1521	!---
1522	last_read = itau_read_n
1523	!-----
1524	!----- Compute mean solar angle for the comming period
1525	!-----
1526	IF (check) WRITE(numout,*) 'Going into solarang', split, one_day
1527	!-----
1528	IF ( dt_force .GT. 3600. ) THEN
1529	mean_coszang(:,:) = 0.0
1530	DO is=1,split
1531	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1532	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1533	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1534	CALL solarang (julian, julian0, iim, jjm, lon, lat, coszang)
1535	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1536	ENDDO
1537	mean_coszang(:,:) = mean_coszang(:,:)/split
1538	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1539	ENDIF
1540	!-----
1541	ENDIF
1542	!---
1543	!--- Do the interpolation
1544	IF (check) WRITE(numout,*) 'Doing the interpolation between time steps'
1545	!---
1546	IF (split > 1) THEN
1547	rw = REAL(itau_split)/split
1548	ELSE
1549	rw = 1.
1550	ENDIF
1551	IF (check) WRITE(numout,*) 'Coeff of interpollation : ',rw
1552	!---
1553	qair(:,:) = (qair_n(:,:)-qair_nm1(:,:))*rw + qair_nm1(:,:)
1554	tair(:,:) = (tair_n(:,:)-tair_nm1(:,:))*rw + tair_nm1(:,:)
1555	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1556	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1557	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1558	IF (is_watchout) THEN
1559	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1560	zlevuv(:,:) = zlev(:,:)
1561	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1562	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1563	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1564	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1565	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1566	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1567	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1568	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1569	ENDIF
1570	!---
1571	!--- Here we need to allow for an option
1572	!--- where radiative energy is conserved
1573	!---
1574	IF ( netrad_cons ) THEN
1575	lwdown(:,:) = lwdown_n(:,:)
1576	ELSE
1577	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1578	ENDIF
1579	!---
1580	!--- For the solar radiation we decompose the mean value
1581	!--- using the zenith angle of the sun if the time step in the forcing data is
1582	!---- more than an hour. Else we use the standard linera interpolation
1583	!----
1584	IF (check) WRITE(numout,*) 'Ready to deal with the solar radiation'
1585	!----
1586	IF ( dt_force .GT. 3600. ) THEN
1587	!---
1588	IF ( netrad_cons ) THEN
1589	WRITE(numout,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
1590	ENDIF
1591	!---
1592	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1593	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1594	!!$ julian = julian_for + rw*dt_force/one_day
1595	IF (check) THEN
1596	WRITE(numout,'(a,f20.10,2I3)') &
1597	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1598	ENDIF
1599	!---
1600	CALL solarang(julian, julian0, iim, jjm, lon, lat, coszang)
1601	!---
1602	WHERE (mean_coszang(:,:) > 0.)
1603	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1604	ELSEWHERE
1605	swdown(:,:) = 0.0
1606	END WHERE
1607	!---
1608	WHERE (swdown(:,:) > 2000. )
1609	swdown(:,:) = 2000.
1610	END WHERE
1611	!---
1612	ELSE
1613	!---
1614	IF ( netrad_cons ) THEN
1615	swdown(:,:) = swdown_n(:,:)
1616	ELSE
1617	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1618	ENDIF
1619	!---
1620	ENDIF
1621	!---
1622	IF (check) THEN
1623	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1624	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1625	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1626	IF (is_watchout) THEN
1627	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1628	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1629	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1630	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1631	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1632	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1633	ENDIF
1634	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1635	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1636	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1637	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1638	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1639	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1640	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1641	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1642	ENDIF
1643	!---
1644	!--- For precip we suppose that the rain
1645	!--- is the sum over the next 6 hours
1646	!---
1647	IF (itau_split <= nb_spread) THEN
1648	rainf(:,:) = rainf_n(:,:)*(split/REAL(nb_spread))
1649	snowf(:,:) = snowf_n(:,:)*(split/REAL(nb_spread))
1650	ELSE
1651	rainf(:,:) = 0.0
1652	snowf(:,:) = 0.0
1653	ENDIF
1654	IF (check) THEN
1655	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1656	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1657	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1658	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1659	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1660	ENDIF
1661	!---
1662	ENDIF ! (daily_interpol)
1663	ENDIF
1664	!---
1665	!--- Here we might put the call to the weather generator ... one day.
1666	!--- Pour le moment, le branchement entre interpolation et generateur de temps
1667	!--- est fait au-dessus.
1668	!---
1669	!- IF ( initialized .AND. weathergen ) THEN
1670	!- ....
1671	!- ENDIF
1672	!---
1673	!--- At this point the code should be initialized. If not we have a problem !
1674	!---
1675	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
1676	!---
1677	initialized = .TRUE.
1678	!---
1679	ELSE
1680	IF ( .NOT. initialized ) THEN
1681	WRITE(numout,*) 'Why is the code forcing_read not initialized ?'
1682	WRITE(numout,*) 'Have you called it with both time-steps set to zero ?'
1683	CALL ipslerr_p(3,'forcing_read_interpol','Pb in initialization','','')
1684	ENDIF
1685	ENDIF
1686	!
1687	!-------------------------
1688	END SUBROUTINE forcing_read_interpol
1689	!=====================================================================
1690	!-
1691	!=====================================================================
1692	SUBROUTINE forcing_just_read &
1693	& (iim, jjm, zlev, zlev_uv, ttm, itb, ite, &
1694	& swdown, rainf, snowf, tair, &
1695	& u, v, qair, pb, lwdown, &
1696	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
1697	& force_id, wind_N_exists, check)
1698	!---------------------------------------------------------------------
1699	!- iim : Size of the grid in x
1700	!- jjm : size of the grid in y
1701	!- zlev : height of the varibales T and Q
1702	!- zlev_uv : height of the varibales U and V
1703	!- ttm : number of time steps in all in the forcing file
1704	!- itb, ite : index of respectively begin and end of read for each variable
1705	!- swdown : Downward solar radiation (W/m^2)
1706	!- rainf : Rainfall (kg/m^2s)
1707	!- snowf : Snowfall (kg/m^2s)
1708	!- tair : 2m air temperature (K)
1709	!- u and v : 2m (in theory !) wind speed (m/s)
1710	!- qair : 2m humidity (kg/kg)
1711	!- pb : Surface pressure (Pa)
1712	!- lwdown : Downward long wave radiation (W/m^2)
1713	!-
1714	!- From a WATCHOUT file :
1715	!- SWnet : Net surface short-wave flux
1716	!- Eair : Air potential energy
1717	!- petAcoef : Coeficients A from the PBL resolution for T
1718	!- peqAcoef : Coeficients A from the PBL resolution for q
1719	!- petBcoef : Coeficients B from the PBL resolution for T
1720	!- peqBcoef : Coeficients B from the PBL resolution for q
1721	!- cdrag : Surface drag
1722	!- ccanopy : CO2 concentration in the canopy
1723	!- force_id : FLINCOM file id.
1724	!- It is used to close the file at the end of the run.
1725	!- wind_N_exists : if Wind_N and Wind_E are in the file (and not just Wind)
1726	!- check : Prompt for reading
1727	!---------------------------------------------------------------------
1728	IMPLICIT NONE
1729	!-
1730	INTEGER, INTENT(in) :: iim, jjm, ttm
1731	INTEGER, INTENT(in) :: itb, ite
1732	REAL, DIMENSION(iim,jjm), INTENT(out) :: zlev, zlev_uv, &
1733	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown
1734	! for watchout files
1735	REAL, DIMENSION(iim,jjm), INTENT(out) :: &
1736	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
1737	INTEGER, INTENT(in) :: force_id
1738	! if Wind_N and Wind_E are in the file (and not just Wind)
1739	LOGICAL, INTENT(in) :: wind_N_exists
1740	LOGICAL :: check
1741	INTEGER :: i, j
1742	REAL :: rau
1743
1744	!-
1745	!---------------------------------------------------------------------
1746	IF ( daily_interpol ) THEN
1747	CALL flinget_buffer (force_id,'Tmin' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1748	CALL forcing_zoom(data_full, tair)
1749	CALL flinget_buffer (force_id,'precip' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1750	CALL forcing_zoom(data_full, rainf)
1751	ELSE
1752	CALL flinget_buffer (force_id,'Tair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1753	CALL forcing_zoom(data_full, tair)
1754	CALL flinget_buffer (force_id,'Snowf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1755	CALL forcing_zoom(data_full, snowf)
1756	CALL flinget_buffer (force_id,'Rainf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1757	CALL forcing_zoom(data_full, rainf)
1758	ENDIF
1759
1760
1761	CALL flinget_buffer (force_id,'SWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1762	CALL forcing_zoom(data_full, swdown)
1763	CALL flinget_buffer (force_id,'LWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1764	CALL forcing_zoom(data_full, lwdown)
1765
1766	CALL flinget_buffer (force_id,'PSurf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1767	CALL forcing_zoom(data_full, pb)
1768	CALL flinget_buffer (force_id,'Qair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1769	CALL forcing_zoom(data_full, qair)
1770	!---
1771	IF ( wind_N_exists ) THEN
1772	CALL flinget_buffer (force_id,'Wind_N', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1773	CALL forcing_zoom(data_full, u)
1774	CALL flinget_buffer (force_id,'Wind_E', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1775	CALL forcing_zoom(data_full, v)
1776	ELSE
1777	CALL flinget_buffer (force_id,'Wind', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1778	CALL forcing_zoom(data_full, u)
1779	v=0.0
1780	ENDIF
1781
1782	!-
1783	!- Deal with the height of the atmospheric forcing varibles
1784	!-
1785	!----
1786	IF ( zheight ) THEN
1787	zlev(:,:) = zlev_fixed
1788	ELSE IF ( zsigma .OR. zhybrid ) THEN
1789	DO i=1,iim
1790	DO j=1,jjm
1791	IF ( tair(i,j) < val_exp ) THEN
1792	rau = pb(i,j)/(cte_molr*tair(i,j))
1793
1794	zlev(i,j) = (pb(i,j) - (zhybrid_a + zhybrid_bpb(i,j)))/(rau cte_grav)
1795	ELSE
1796	zlev(i,j) = 0.0
1797	ENDIF
1798	ENDDO
1799	ENDDO
1800	ELSE IF ( zlevels ) THEN
1801	CALL flinget_buffer (force_id,'Levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1802	CALL forcing_zoom(data_full, zlev)
1803	ELSE
1804	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1805	& 'We cannot determine the height for T and Q.','stop readdim2')
1806	ENDIF
1807
1808	IF ( zsamelev_uv ) THEN
1809	zlev_uv(:,:) = zlev(:,:)
1810	ELSE
1811	IF ( zheight ) THEN
1812	zlev_uv(:,:) = zlevuv_fixed
1813	ELSE IF ( zsigma .OR. zhybrid ) THEN
1814	DO i=1,iim
1815	DO j=1,jjm
1816	IF ( tair(i,j) < val_exp ) THEN
1817	rau = pb(i,j)/(cte_molr*tair(i,j))
1818
1819	zlev_uv(i,j) = (pb(i,j) - (zhybriduv_a + zhybriduv_bpb(i,j)))/(rau cte_grav)
1820	ELSE
1821	zlev_uv(i,j) = 0.0
1822	ENDIF
1823	ENDDO
1824	ENDDO
1825	ELSE IF ( zlevels ) THEN
1826	CALL flinget_buffer (force_id,'Levels_uv', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1827	CALL forcing_zoom(data_full, zlev_uv)
1828	ELSE
1829	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1830	& 'We cannot determine the height for U and V.','stop readdim2')
1831	ENDIF
1832	ENDIF
1833	!----
1834	IF ( is_watchout ) THEN
1835	CALL flinget_buffer (force_id,'levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1836	CALL forcing_zoom(data_full, zlev)
1837	!
1838	! If we are in WATHCOUT it means T,Q are at the same height as U,V
1839	!
1840	zlev_uv(:,:) = zlev(:,:)
1841	! Net surface short-wave flux
1842	CALL flinget_buffer (force_id,'SWnet', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1843	CALL forcing_zoom(data_full, SWnet)
1844	! Air potential energy
1845	CALL flinget_buffer (force_id,'Eair', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1846	CALL forcing_zoom(data_full, Eair)
1847	! Coeficients A from the PBL resolution for T
1848	CALL flinget_buffer (force_id,'petAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1849	CALL forcing_zoom(data_full, petAcoef)
1850	! Coeficients A from the PBL resolution for q
1851	CALL flinget_buffer (force_id,'peqAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1852	CALL forcing_zoom(data_full, peqAcoef)
1853	! Coeficients B from the PBL resolution for T
1854	CALL flinget_buffer (force_id,'petBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1855	CALL forcing_zoom(data_full, petBcoef)
1856	! Coeficients B from the PBL resolution for q
1857	CALL flinget_buffer (force_id,'peqBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1858	CALL forcing_zoom(data_full, peqBcoef)
1859	! Surface drag
1860	CALL flinget_buffer (force_id,'cdrag', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1861	CALL forcing_zoom(data_full, cdrag)
1862	! CO2 concentration in the canopy
1863	CALL flinget_buffer (force_id,'ccanopy', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1864	CALL forcing_zoom(data_full, ccanopy)
1865	ENDIF
1866	!
1867	!----
1868	IF (check) WRITE(numout,*) 'Variables have been extracted between ',itb,' and ',ite,' iterations of the forcing file.'
1869	!-------------------------
1870	END SUBROUTINE forcing_just_read
1871	!=====================================================================
1872
1873	!-
1874	SUBROUTINE forcing_just_read_tmax &
1875	& (iim, jjm, ttm, itb, ite, tmax, force_id )
1876	!---------------------------------------------------------------------
1877	!- iim : Size of the grid in x
1878	!- jjm : size of the grid in y
1879	!- ttm : number of time steps in all in the forcing file
1880	!- itb, ite : index of respectively begin and end of read for each variable
1881	!- tmax : 2m air temperature (K)
1882	!- force_id : FLINCOM file id.
1883	!- It is used to close the file at the end of the run.
1884	!---------------------------------------------------------------------
1885	IMPLICIT NONE
1886	!-
1887	INTEGER, INTENT(in) :: iim, jjm, ttm
1888	INTEGER, INTENT(in) :: itb, ite
1889	REAL, DIMENSION(iim,jjm), INTENT(out) :: tmax
1890	INTEGER, INTENT(in) :: force_id
1891	!-
1892	!---------------------------------------------------------------------
1893	CALL flinget_buffer (force_id,'Tmax' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1894	CALL forcing_zoom(data_full, tmax)
1895	!-------------------------
1896	END SUBROUTINE forcing_just_read_tmax
1897	!=====================================================================
1898
1899	!-
1900	SUBROUTINE forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
1901	!---
1902	!--- This subroutine finds the indices of the land points over which the land
1903	!--- surface scheme is going to run.
1904	!---
1905	IMPLICIT NONE
1906	!-
1907	!- ARGUMENTS
1908	!-
1909	INTEGER, INTENT(IN) :: iim, jjm
1910	REAL, INTENT(IN) :: tair(iim,jjm)
1911	INTEGER, INTENT(OUT) :: i_test, j_test, nbindex
1912	INTEGER, INTENT(OUT) :: kindex(iim*jjm)
1913	LOGICAL :: check
1914	!-
1915	!- LOCAL
1916	INTEGER :: i, j, ig
1917	!-
1918	!-
1919	ig = 0
1920	i_test = 0
1921	j_test = 0
1922	!---
1923	IF (MINVAL(tair(:,:)) < 100.) THEN
1924	!----- In this case the 2m temperature is in Celsius
1925	DO j=1,jjm
1926	DO i=1,iim
1927	IF (tair(i,j) < 100.) THEN
1928	ig = ig+1
1929	kindex(ig) = (j-1)*iim+i
1930	!
1931	! Here we find at random a land-point on which we can do
1932	! some printouts for test.
1933	!
1934	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1935	i_test = i
1936	j_test = j
1937	IF (check) THEN
1938	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1939	ENDIF
1940	ENDIF
1941	ENDIF
1942	ENDDO
1943	ENDDO
1944	ELSE
1945	!----- 2m temperature is in Kelvin
1946	DO j=1,jjm
1947	DO i=1,iim
1948	IF (tair(i,j) < 500.) THEN
1949	ig = ig+1
1950	kindex(ig) = (j-1)*iim+i
1951	!
1952	! Here we find at random a land-point on which we can do
1953	! some printouts for test.
1954	!
1955	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1956	i_test = i
1957	j_test = j
1958	IF (check) THEN
1959	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1960	ENDIF
1961	ENDIF
1962	ENDIF
1963	ENDDO
1964	ENDDO
1965	ENDIF
1966	!---
1967	nbindex = ig
1968	!---
1969	END SUBROUTINE forcing_landind
1970	!-
1971	!=====================================================================
1972	!-
1973	SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,init_f)
1974	!-
1975	!- This subroutine calculates the longitudes and latitudes of the model grid.
1976	!-
1977	IMPLICIT NONE
1978	!-
1979	INTEGER, INTENT(in) :: iim, jjm, llm
1980	LOGICAL, INTENT(in) :: init_f
1981	REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat
1982	!-
1983	INTEGER :: i,j
1984	!-
1985	!- Should be unified one day
1986	!-
1987	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : options : ', weathergen, interpol
1988	!-
1989	IF ( weathergen ) THEN
1990	IF (init_f) THEN
1991	DO i = 1, iim
1992	lon(i,:) = limit_west + merid_res/2. + &
1993	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim)
1994	ENDDO
1995	!-
1996	DO j = 1, jjm
1997	lat(:,j) = limit_north - zonal_res/2. - &
1998	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm)
1999	ENDDO
2000	ELSE
2001	IF (is_root_prc) THEN
2002	DO i = 1, iim_g
2003	lon_g(i,:) = limit_west + merid_res/2. + &
2004	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim_g)
2005	ENDDO
2006	!-
2007	DO j = 1, jjm_g
2008	lat_g(:,j) = limit_north - zonal_res/2. - &
2009	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm_g)
2010	ENDDO
2011	ENDIF
2012	CALL bcast(lon_g)
2013	CALL bcast(lat_g)
2014	lon=lon_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2015	lat=lat_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2016	ENDIF
2017	!-
2018	ELSEIF ( interpol ) THEN
2019	!-
2020	CALL forcing_zoom(lon_full, lon)
2021	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lon'
2022	CALL forcing_zoom(lat_full, lat)
2023	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lat'
2024
2025	ELSE
2026	CALL ipslerr_p(3,'forcing_grid','Neither interpolation nor weather generator is specified.','','')
2027	ENDIF
2028
2029	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : ended'
2030
2031	END SUBROUTINE forcing_grid
2032	!-
2033	!=====================================================================
2034	!-
2035	SUBROUTINE forcing_zoom(x_f, x_z)
2036	!-
2037	!- This subroutine takes the slab of data over which we wish to run the model.
2038	!-
2039	IMPLICIT NONE
2040	!-
2041	REAL, INTENT(IN) :: x_f(iim_full, jjm_full)
2042	REAL, INTENT(OUT) :: x_z(iim_zoom, jjm_zoom)
2043	!-
2044	INTEGER :: i,j
2045	!-
2046	DO i=1,iim_zoom
2047	DO j=1,jjm_zoom
2048	x_z(i,j) = x_f(i_index(i),j_index(j))
2049	ENDDO
2050	ENDDO
2051	!-
2052	END SUBROUTINE forcing_zoom
2053
2054	!-
2055	!=====================================================================
2056	!-
2057
2058	SUBROUTINE forcing_vertical_ioipsl(force_id)
2059	!
2060	!- This subroutine explores the forcing file to decide what information is available
2061	!- on the vertical coordinate.
2062	!
2063	INTEGER, INTENT(IN) :: force_id
2064	!
2065	LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
2066	LOGICAL :: foundvar
2067	LOGICAL :: levlegacy
2068
2069	!
2070	!- Set all the defaults
2071	!
2072	zfixed=.FALSE.
2073	zsigma=.FALSE.
2074	zhybrid=.FALSE.
2075	zlevels=.FALSE.
2076	zheight=.FALSE.
2077	zsamelev_uv = .TRUE.
2078	levlegacy = .FALSE.
2079	!
2080	foundvar = .FALSE.
2081	!
2082	!- We have a forcing file to explore so let us see if we find any of the conventions
2083	!- which allow us to find the height of T,Q,U and V.
2084	!
2085	IF ( force_id > 0 ) THEN
2086	!
2087	! Case for sigma levels
2088	!
2089	IF ( .NOT. foundvar ) THEN
2090	CALL flinquery_var(force_id, 'Sigma', var_exists)
2091	CALL flinquery_var(force_id, 'Sigma_uv', varuv_exists)
2092	IF ( var_exists ) THEN
2093	foundvar = .TRUE.
2094	zsigma = .TRUE.
2095	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2096	ENDIF
2097	ENDIF
2098	!
2099	! Case for Hybrid levels
2100	!
2101	IF ( .NOT. foundvar ) THEN
2102	CALL flinquery_var(force_id, 'HybSigA', vara_exists)
2103	IF ( vara_exists ) THEN
2104	CALL flinquery_var(force_id, 'HybSigB', varb_exists)
2105	IF ( varb_exists ) THEN
2106	var_exists=.TRUE.
2107	ELSE
2108	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','Missing the B coefficient for', &
2109	& 'Hybrid vertical levels for T and Q','stop readdim2')
2110	ENDIF
2111	ENDIF
2112	CALL flinquery_var(force_id, 'HybSigA_uv', vara_exists)
2113	IF ( vara_exists ) THEN
2114	CALL flinquery_var(force_id, 'HybSigB_uv', varb_exists)
2115	IF ( varb_exists ) THEN
2116	varuv_exists=.TRUE.
2117	ELSE
2118	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','Missing the B coefficient for', &
2119	& 'Hybrid vertical levels for U and V','stop readdim2')
2120	ENDIF
2121	ENDIF
2122	IF ( var_exists ) THEN
2123	foundvar = .TRUE.
2124	zhybrid = .TRUE.
2125	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2126	ENDIF
2127	ENDIF
2128	!
2129	! Case for levels (i.e. a 2d time varying field with the height in meters)
2130	!
2131	IF ( .NOT. foundvar ) THEN
2132	CALL flinquery_var(force_id, 'Levels', var_exists)
2133	CALL flinquery_var(force_id, 'Levels_uv', varuv_exists)
2134	IF ( var_exists ) THEN
2135	foundvar = .TRUE.
2136	zlevels = .TRUE.
2137	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2138	ENDIF
2139	ENDIF
2140	!
2141	! Case where a fixed height is provided in meters
2142	!
2143	IF ( .NOT. foundvar ) THEN
2144	CALL flinquery_var(force_id, 'Height_Lev1', var_exists)
2145	CALL flinquery_var(force_id, 'Height_Levuv', varuv_exists)
2146	IF ( var_exists ) THEN
2147	foundvar = .TRUE.
2148	zheight = .TRUE.
2149	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2150	ENDIF
2151	ENDIF
2152	!
2153	! Case where a fixed height is provided in meters in the lev variable
2154	!
2155	IF ( .NOT. foundvar ) THEN
2156	CALL flinquery_var(force_id, 'lev', var_exists)
2157	IF ( var_exists ) THEN
2158	foundvar = .TRUE.
2159	zheight = .TRUE.
2160	levlegacy = .TRUE.
2161	ENDIF
2162	ENDIF
2163	!
2164	ENDIF
2165	!
2166	! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
2167	! except the case zlevels
2168	!
2169	IF ( foundvar .AND. .NOT. zlevels ) THEN
2170	!
2171	IF ( zheight ) THEN
2172	!
2173	! Constant height
2174	!
2175	IF ( levlegacy ) THEN
2176	CALL flinget (force_id,'lev',1, 1, 1, 1, 1, 1, zlev_fixed)
2177	ELSE
2178	CALL flinget (force_id,'Height_Lev1',1, 1, 1, 1, 1, 1, zlev_fixed)
2179	IF ( .NOT. zsamelev_uv ) THEN
2180	CALL flinget (force_id,'Height_Levuv',1, 1, 1, 1, 1, 1, zlevuv_fixed)
2181	ENDIF
2182	ENDIF
2183	WRITE(numout,*) "forcing_vertical_ioipsl : case ZLEV : Read from forcing file :", zlev_fixed, zlevuv_fixed
2184	!
2185	ELSE IF ( zsigma .OR. zhybrid ) THEN
2186	!
2187	! Sigma or hybrid levels
2188	!
2189	IF ( zsigma ) THEN
2190	CALL flinget (force_id,'Sigma',1, 1, 1, 1, 1, 1, zhybrid_b)
2191	zhybrid_a = zero
2192	IF ( .NOT. zsamelev_uv ) THEN
2193	CALL flinget (force_id,'Sigma_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2194	zhybriduv_a = zero
2195	ENDIF
2196	ELSE
2197	CALL flinget (force_id,'HybSigB',1, 1, 1, 1, 1, 1, zhybrid_b)
2198	CALL flinget (force_id,'HybSigA',1, 1, 1, 1, 1, 1, zhybrid_a)
2199	IF ( .NOT. zsamelev_uv ) THEN
2200	CALL flinget (force_id,'HybSigB_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2201	CALL flinget (force_id,'HybSigA_uv',1, 1, 1, 1, 1, 1, zhybriduv_a)
2202	ENDIF
2203	ENDIF
2204	WRITE(numout,*) "forcing_vertical_ioipsl : case Pressure coordinates : "
2205	WRITE(numout,*) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
2206	ELSE
2207	!
2208	! Why are we here ???
2209	!
2210	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','What is the option used to describe the height of', &
2211	& 'the atmospheric forcing ?','Please check your forcing file.')
2212	ENDIF
2213	ENDIF
2214	!
2215	!- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
2216	!- read what has been specified by the user.
2217	!
2218	IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
2219	!
2220	!-
2221	!Config Key = HEIGHT_LEV1
2222	!Config Desc = Height at which T and Q are given
2223	!Config Def = 2.0
2224	!Config If = offline mode
2225	!Config Help = The atmospheric variables (temperature and specific
2226	!Config humidity) are measured at a specific level.
2227	!Config The height of this level is needed to compute
2228	!Config correctly the turbulent transfer coefficients.
2229	!Config Look at the description of the forcing
2230	!Config DATA for the correct value.
2231	!Config Units = [m]
2232	!-
2233	zlev_fixed = 2.0
2234	CALL getin('HEIGHT_LEV1', zlev_fixed)
2235
2236	!-
2237	!Config Key = HEIGHT_LEVW
2238	!Config Desc = Height at which the wind is given
2239	!Config Def = 10.0
2240	!Config If = offline mode
2241	!Config Help = The height at which wind is needed to compute
2242	!Config correctly the turbulent transfer coefficients.
2243	!Config Units= [m]
2244	!-
2245	zlevuv_fixed = 10.0
2246	CALL getin('HEIGHT_LEVW', zlevuv_fixed)
2247
2248	zheight = .TRUE.
2249
2250	IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
2251	zsamelev_uv = .FALSE.
2252	ENDIF
2253
2254	CALL ipslerr ( 2, 'forcing_vertical_ioipsl','The height of the atmospheric forcing variables', &
2255	& 'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
2256	ENDIF
2257
2258	END SUBROUTINE forcing_vertical_ioipsl
2259
2260	!-
2261	!=====================================================================
2262	!-
2263
2264	SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, &
2265	& iim_f, jjm_f, lon, lat, iim, jjm, iind, jind)
2266
2267	IMPLICIT NONE
2268	!
2269	! ARGUMENTS
2270	!
2271	REAL, INTENT(inout) :: limit_west,limit_east,limit_north,limit_south
2272	INTEGER, INTENT(in) :: iim_f, jjm_f
2273	REAL, INTENT(in) :: lon(iim_f, jjm_f), lat(iim_f, jjm_f)
2274	INTEGER, INTENT(out) :: iim,jjm
2275	INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f)
2276	!
2277	! LOCAL
2278	!
2279	INTEGER :: i, j
2280	REAL :: lolo
2281	LOGICAL :: over_dateline = .FALSE.
2282	!
2283	!
2284	IF ( ( ABS(limit_east) .GT. 180. ) .OR. &
2285	( ABS(limit_west) .GT. 180. ) ) THEN
2286	WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
2287	CALL ipslerr_p (3,'domain_size', &
2288	& 'Longitudes problem.','In run.def file :', &
2289	& 'limit_east > 180. or limit_west > 180.')
2290	ENDIF
2291	!
2292	IF ( limit_west .GT. limit_east ) over_dateline = .TRUE.
2293	!
2294	IF ( ( limit_south .LT. -90. ) .OR. &
2295	( limit_north .GT. 90. ) .OR. &
2296	( limit_south .GE. limit_north ) ) THEN
2297	WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
2298	CALL ipslerr_p (3,'domain_size', &
2299	& 'Latitudes problem.','In run.def file :', &
2300	& 'limit_south < -90. or limit_north > 90. or limit_south >= limit_north')
2301	ENDIF
2302	!
2303	! Here we assume that the grid of the forcing data is regular. Else we would have
2304	! to do more work to find the index table.
2305	!
2306	iim = 0
2307	DO i=1,iim_f
2308	!
2309	lolo = lon(i,1)
2310	IF ( lon(i,1) .GT. 180. ) lolo = lon(i,1) - 360.
2311	IF ( lon(i,1) .LT. -180. ) lolo = lon(i,1) + 360.
2312	!
2313	IF (lon(i,1) < limit_west) iim_g_begin = i+1
2314	IF (lon(i,1) < limit_east) iim_g_end = i
2315	!
2316	IF ( over_dateline ) THEN
2317	IF ( lolo .LE. limit_west .OR. lolo .GE. limit_east ) THEN
2318	iim = iim + 1
2319	iind(iim) = i
2320	ENDIF
2321	ELSE
2322	IF ( lolo .GE. limit_west .AND. lolo .LE. limit_east ) THEN
2323	iim = iim + 1
2324	iind(iim) = i
2325	ENDIF
2326	ENDIF
2327	!
2328	ENDDO
2329	!
2330	jjm = 0
2331	DO j=1,jjm_f
2332	IF (lat(1,j) > limit_north) jjm_g_begin = j+1
2333	IF (lat(1,j) > limit_south) jjm_g_end = j
2334	!
2335	IF ( lat(1,j) .GE. limit_south .AND. lat(1,j) .LE. limit_north) THEN
2336	jjm = jjm + 1
2337	jind(jjm) = j
2338	ENDIF
2339	ENDDO
2340	!
2341	WRITE(numout,*) 'Domain zoom size: iim, jjm = ', iim, jjm
2342	END SUBROUTINE domain_size
2343	!-
2344	!=====================================================================
2345	!-
2346	SUBROUTINE flinget_buffer(force_id, varname, iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
2347
2348	!! This subroutine is a wrap of flinget/IOIPSL. The arguments are the same.
2349	!! flinget_buffer will call flinget and buffer the forcing data localy in this subroutine.
2350	!! According to the variable NBUFF set in run.def, several time steps can be read at the same time
2351	!! from the forcing file. If NBUFF=0, the full forcing file is read.
2352	!! The output, data_full, from this subroutine is always only one time step of corresponding to itb.
2353	!! itb must be equal to ite.
2354
2355	!! Input arguments
2356	INTEGER, INTENT(in) :: force_id !! Id for forcing file
2357	CHARACTER(len=*), INTENT(in) :: varname !! Name of current variable to be read
2358	INTEGER, INTENT(in) :: iim_full, jjm_full, llm_full !! Horizontal and vertical domaine
2359	INTEGER, INTENT(in) :: ttm !! Full lenght of forcing file
2360	INTEGER, INTENT(in) :: itb, ite !! Time step to be read from forcing file. itb must be equal to ite
2361
2362	!! Output argument
2363	REAL, DIMENSION(iim_full, jjm_full), INTENT(out) :: data_full !! Data for time step itb.
2364
2365	!! Define specific type to buffer data together with name and index
2366	TYPE buffer_type
2367	CHARACTER(len=20) :: name !! Name of variable in forcing file
2368	INTEGER :: istart !! Start index of current buffered data
2369	INTEGER :: iend !! End index of current buffered data
2370	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: data !! Data read from forcing file for intervall [istart,iend]
2371	END TYPE buffer_type
2372
2373	!! Local variables
2374	INTEGER, PARAMETER :: maxvar=20 !! Max number of variables to be buffered
2375	TYPE(buffer_type), DIMENSION(maxvar),SAVE :: data_buffer !! Containing all variables and the current buffered data
2376	INTEGER, SAVE :: nbuff !! Number of time steps to be buffered
2377	INTEGER, SAVE :: lastindex=0 !! Current number of variables stored in data_buffer
2378	INTEGER, SAVE :: ttm0 !! Time lenght of forcing file, stored for test purpose
2379	LOGICAL, SAVE :: first=.TRUE. !! First call to this subroutine
2380	INTEGER :: index !! Index in data_buffer for current variable
2381	INTEGER :: i, ierr !! Loop and error variables
2382
2383
2384	!! 1. Initialization
2385	IF (first) THEN
2386	data_buffer(:)%name='undef'
2387	! Read NBUFF from run.def
2388	! Note that getin_p is not used because this subroutine might be called only by master process
2389
2390	!Config Key = NBUFF
2391	!Config Desc = Number of time steps of data to buffer between each reading of the forcing file
2392	!Config If = OFF_LINE
2393	!Config Help = The full simulation time length will be read if NBUFF equal 0. NBUFF > 1 can be used for smaller regions or site simulations only.
2394	!Config Def = 1
2395	!Config Units= -
2396
2397	nbuff=1
2398	CALL getin('NBUFF', nbuff)
2399
2400	IF (nbuff == 0 .OR. nbuff >ttm) THEN
2401	! Set nbuff as the full forcing file lenght
2402	nbuff=ttm
2403	ELSE IF (nbuff < 0) THEN
2404	! Negativ nbuff not possible
2405	CALL ipslerr_p(3,'flinget_buffer','NBUFF must be a positiv number','Set NBUFF=0 for full simulation lenght','')
2406	END IF
2407	WRITE(numout,*)'flinget_buffer: NBUFF=',nbuff,' number of time step will be buffered'
2408	WRITE(numout,*)'flinget_buffer: Choose a lower value for NBUFF if problem with memory'
2409
2410	! Save dimensions to check following timesteps
2411	! ttm is the full lenght of forcing file
2412	ttm0=ttm
2413
2414	first=.FALSE.
2415	END IF
2416
2417	!! 2. Coeherence tests on input arguments
2418	IF (ttm /= ttm0) THEN
2419	WRITE(numout,*)'Problem with ttm=',ttm,' ttm0=',ttm0
2420	CALL ipslerr_p(3,'flinget_buffer','Error with ttm and ttm0','','')
2421	END IF
2422	IF (itb /= ite) THEN
2423	WRITE(numout,*) 'There is a problem. Why is itb not equal ite ?'
2424	WRITE(numout,*) 'itb=',itb,' ite=',ite,' varname=',varname
2425	CALL ipslerr_p(3,'flinget_buffer','ite not equal itb','','')
2426	END IF
2427
2428
2429	!! 3. Find index for current variable
2430	index=0
2431	DO i=1, maxvar
2432	IF ( trim(varname) == data_buffer(i)%name ) THEN
2433	index=i
2434	CYCLE
2435	END IF
2436	END DO
2437
2438	!! 4. Initialize and allocate if necesary the current variable
2439	IF ( index == 0 ) THEN
2440	! The variable was not found
2441	! This must be the first time for current variable
2442	index=lastindex+1
2443	lastindex=index
2444	IF (index > maxvar) CALL ipslerr_p(3,'flinget_buffer','to many variables','maxvar is too small','')
2445
2446	! Initialize the data_buffer for this index
2447	data_buffer(index)%name=trim(varname)
2448	ALLOCATE(data_buffer(index)%data(iim_full,jjm_full,nbuff),stat=ierr)
2449	IF (ierr /= 0) CALL ipslerr_p(3,'flinget_buffer','pb alloc data_buffer%data','for variable=',varname)
2450	data_buffer(index)%istart=0
2451	data_buffer(index)%iend=0
2452	END IF
2453
2454
2455	!! 5. Call flinget if current time step (itb) is outside the buffered intervall
2456	IF (( itb > data_buffer(index)%iend ) .OR. ( itb < data_buffer(index)%istart )) THEN
2457	! itb is not in the time slice previously read or it is the first time to read
2458	! Reading of forcing file will now be done
2459	! First recalculate index to be read
2460	data_buffer(index)%istart = itb
2461	data_buffer(index)%iend = itb + nbuff - 1
2462
2463	! WRITE(numout,*) 'Now do flinget for ',varname,', itb=',itb,', istart=',&
2464	! data_buffer(index)%istart,', iend=',data_buffer(index)%iend
2465	CALL flinget (force_id,varname, iim_full, jjm_full, llm_full, ttm, data_buffer(index)%istart, &
2466	data_buffer(index)%iend, data_buffer(index)%data(:,:,:))
2467	END IF
2468
2469	!! 6. Initialize the output variable with data from buffered variable
2470	! Find index for the time step corrsponding to itb in the time slice previously read from forcing file
2471	i=itb-data_buffer(index)%istart+1
2472	! Initialize output variable
2473	data_full(:,:) = data_buffer(index)%data(:,:,i)
2474
2475
2476	END SUBROUTINE flinget_buffer
2477	!------------------
2478	END MODULE readdim2

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

Download in other formats: