Context Navigation

source: branches/publications/ORCHIDEE-LEAK-r5919/src_driver/readdim2.f90 @ 7346

Last change on this file since 7346 was 3283, checked in by josefine.ghattas, 8 years ago
Ticket #233 : Stop if more processors than land points Albert Jornet-Puig
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 97.4 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(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 set_grid_glo(iim_zoom,jjm_zoom,NbPoint)
328	CALL allocate_grid_glo
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(-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 set_grid_glo(iim_zoom,jjm_zoom,NbPoint)
451	CALL allocate_grid_glo
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 init_grid ( nbindex )
622	!---
623	!--- Compute
624	!---
625	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, kindex)
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	! if Wind_N and Wind_E are in the file (and not just Wind)
770	LOGICAL, SAVE :: wind_N_exists=.FALSE.
771	LOGICAL :: wind_E_exists=.FALSE.
772	LOGICAL, SAVE :: contfrac_exists=.FALSE.
773	LOGICAL, SAVE :: neighbours_exists=.FALSE.
774	LOGICAL, SAVE :: initialized = .FALSE.
775	LOGICAL :: check=.FALSE.
776	! to bypass FPE problem on integer convertion with missing_value heigher than precision
777	INTEGER, PARAMETER :: undef_int = 999999999
778	!---------------------------------------------------------------------
779
780	itau_read = MOD((itauin-1),ttm)+1
781
782	IF (check) THEN
783	WRITE(numout,*) &
784	" FORCING READ : itauin, itau_read, itau_split : ",&
785	itauin, itau_read, itau_split
786	ENDIF
787
788	!-
789	!- This part initializes the reading of the forcing. As you can see
790	!- we only go through here if both time steps are zero.
791	!-
792	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
793	!-
794	!- Tests on forcing file type
795	CALL flinquery_var(force_id, 'Wind_N', wind_N_exists)
796	IF ( wind_N_exists ) THEN
797	CALL flinquery_var(force_id, 'Wind_E', wind_E_exists)
798	IF ( .NOT. wind_E_exists ) THEN
799	CALL ipslerr_p(3,'forcing_read_interpol', &
800	& 'Variable Wind_E does not exist in forcing file', &
801	& 'But variable Wind_N exists.','Please, rename Wind_N to Wind;')
802	ENDIF
803	ENDIF
804	CALL flinquery_var(force_id, 'levels', is_watchout)
805	IF ( is_watchout ) THEN
806	WRITE(numout,*) "Read a Watchout File."
807	ENDIF
808	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
809	!-
810	IF (check) WRITE(numout,*) 'ALLOCATE all the memory needed'
811	!-
812	ALLOCATE &
813	& (swdown_nm1(iim,jjm), rainf_nm1(iim,jjm), snowf_nm1(iim,jjm), &
814	& tair_nm1(iim,jjm), u_nm1(iim,jjm), v_nm1(iim,jjm), qair_nm1(iim,jjm), &
815	& pb_nm1(iim,jjm), lwdown_nm1(iim,jjm))
816	ALLOCATE &
817	& (swdown_n(iim,jjm), rainf_n(iim,jjm), snowf_n(iim,jjm), &
818	& tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), &
819	& pb_n(iim,jjm), lwdown_n(iim,jjm))
820
821	IF(daily_interpol) THEN
822	ALLOCATE &
823	& (startday_n(iim,jjm), startday_nm1(iim,jjm), daylength_n(iim,jjm), &
824	& daylength_nm1(iim,jjm), tmax_n(iim,jjm), tmax_nm1(iim,jjm), tmin_n(iim,jjm), &
825	& tmin_nm1(iim,jjm), tmin_nm2(iim,jjm), qsatta(iim,jjm), qsattmin_n(iim,jjm), qsattmin_nm1(iim,jjm), &
826	& qmin_n(iim,jjm), qmin_nm1(iim,jjm), qmax_n(iim,jjm), qmax_nm1(iim,jjm), qsa(iim,jjm) )
827	ENDIF
828
829
830	ALLOCATE &
831	& (zlev_nm1(iim,jjm), zlev_n(iim,jjm), zlevuv_nm1(iim,jjm), zlevuv_n(iim,jjm), &
832	& SWnet_nm1(iim,jjm), Eair_nm1(iim,jjm), cdrag_nm1(iim,jjm), ccanopy_nm1(iim,jjm), &
833	& petAcoef_nm1(iim,jjm), peqAcoef_nm1(iim,jjm), petBcoef_nm1(iim,jjm), peqBcoef_nm1(iim,jjm), &
834	& SWnet_n(iim,jjm), Eair_n(iim,jjm), cdrag_n(iim,jjm), ccanopy_n(iim,jjm), &
835	& petAcoef_n(iim,jjm), peqAcoef_n(iim,jjm), petBcoef_n(iim,jjm), peqBcoef_n(iim,jjm))
836	ALLOCATE &
837	& (mean_coszang(iim,jjm))
838	!-
839	IF (check) WRITE(numout,*) 'Memory ALLOCATED'
840	!-
841	!- We need for the driver surface air temperature and humidity before the
842	!- the loop starts. Thus we read it here.
843	!-
844	CALL forcing_just_read (iim, jjm, zlev, zlevuv, ttm, 1, 1, &
845	& swdown, rainf, snowf, tair, &
846	& u, v, qair, pb, lwdown, &
847	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
848	& force_id, wind_N_exists, check)
849	!----
850
851	!-- First in case it's not a watchout file
852	IF ( .NOT. is_watchout ) THEN
853	IF ( contfrac_exists ) THEN
854	WRITE(numout,*) "contfrac exist in the forcing file."
855	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
856	CALL forcing_zoom(data_full, fcontfrac)
857	WRITE(numout,*) "fcontfrac min/max :",MINVAL(fcontfrac(1:iim_zoom,jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,jjm_zoom))
858	!
859	! We need to make sure that when we gather the points we pick all
860	! the points where contfrac is above 0. Thus we prepare tair for
861	! subroutine forcing_landind
862	!
863	DO i=1,iim
864	DO j=1,jjm
865	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0. ! bande de recouvrement du scatter2D
866	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0. ! => on mets les données qu'on veut pas du noeud à missing_value
867	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
868	tair(i,j) = 999999.
869	ENDIF
870	ENDDO
871	ENDDO
872	ELSE
873	fcontfrac(:,:) = 1.0
874	ENDIF
875	!---
876	!--- Create the index table
877	!---
878	!--- This job return a LOCAL kindex
879	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
880	#ifdef CPP_PARA
881	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
882	! Force nbindex points for parallel computation
883	nbindex=nbp_loc
884	CALL scatter(index_g,kindex(1:nbindex))
885	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
886	#endif
887	ik=MAX(nbindex/2,1)
888	j_test = (((kindex(ik)-1)/iim) + 1)
889	i_test = (kindex(ik) - (j_test-1)*iim)
890	IF (check) THEN
891	WRITE(numout,*) 'New test point is : ', i_test, j_test
892	ENDIF
893	!---
894	!--- Allocate grid stuff
895	!---
896	CALL init_grid ( nbindex )
897	!---
898	!--- All grid variables
899	!---
900	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, kindex)
901	DO ik=1,nbindex
902	!
903	j = ((kindex(ik)-1)/iim) + 1
904	i = (kindex(ik) - (j-1)*iim)
905	!-
906	!- Store variable to help describe the grid
907	!- once the points are gathered.
908	!-
909	fneighbours(i,j,:) = neighbours(ik,:)
910	!
911	fresolution(i,j,:) = resolution(ik,:)
912	ENDDO
913	ELSE
914	!-- Second, in case it is a watchout file
915	ALLOCATE (tmpdata(iim,jjm))
916	tmpdata(:,:) = 0.0
917	!--
918	IF ( .NOT. contfrac_exists ) THEN
919	CALL ipslerr_p (3,'forcing_read_interpol', &
920	& 'Could get contfrac variable in a watchout file :',TRIM(filename), &
921	& '(Problem with file ?)')
922	ENDIF
923	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
924	CALL forcing_zoom(data_full, fcontfrac)
925	!
926	! We need to make sure that when we gather the points we pick all
927	! the points where contfrac is above 0. Thus we prepare tair for
928	! subroutine forcing_landind
929	!
930	DO i=1,iim
931	DO j=1,jjm
932	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0.
933	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0.
934	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
935	tair(i,j) = 999999.
936	ENDIF
937	ENDDO
938	ENDDO
939	!---
940	!--- Create the index table
941	!---
942	!--- This job return a LOCAL kindex
943	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
944	#ifdef CPP_PARA
945	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
946	! Force nbindex points for parallel computation
947	nbindex=nbp_loc
948	CALL scatter(index_g,kindex)
949	kindex(:)=kindex(:)-offset
950	! kindex(:)=kindex(:)-(jj_begin-1)*iim_g
951	#endif
952	ik=MAX(nbindex/2,1)
953	j_test = (((kindex(ik)-1)/iim) + 1)
954	i_test = (kindex(ik) - (j_test-1)*iim)
955	IF (check) THEN
956	WRITE(numout,*) 'New test point is : ', i_test, j_test
957	ENDIF
958	!---
959	!--- Allocate grid stuff
960	!---
961	CALL init_grid ( nbindex )
962	neighbours(:,:) = -1
963	resolution(:,:) = 0.
964	min_resol(:) = 1.e6
965	max_resol(:) = -1.
966	!---
967	!--- All grid variables
968	!---
969	!-
970	!- Get variables to help describe the grid
971	CALL flinquery_var(force_id, 'neighboursNN', neighbours_exists)
972	IF ( .NOT. neighbours_exists ) THEN
973	CALL ipslerr_p (3,'forcing_read_interpol', &
974	& 'Could get neighbours in a watchout file :',TRIM(filename), &
975	& '(Problem with file ?)')
976	ELSE
977	WRITE(numout,*) "Watchout file contains neighbours and resolutions."
978	ENDIF
979	!
980	fneighbours(:,:,:) = undef_int
981	!
982	!- once the points are gathered.
983	CALL flinget_buffer (force_id,'neighboursNN',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
984	CALL forcing_zoom(data_full, tmpdata)
985	WHERE(tmpdata(:,:) < undef_int)
986	fneighbours(:,:,1) = NINT(tmpdata(:,:))
987	ENDWHERE
988	!
989	CALL flinget_buffer (force_id,'neighboursNE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
990	CALL forcing_zoom(data_full, tmpdata)
991	WHERE(tmpdata(:,:) < undef_int)
992	fneighbours(:,:,2) = NINT(tmpdata(:,:))
993	ENDWHERE
994	!
995	CALL flinget_buffer (force_id,'neighboursEE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
996	CALL forcing_zoom(data_full, tmpdata)
997	WHERE(tmpdata(:,:) < undef_int)
998	fneighbours(:,:,3) = NINT(tmpdata(:,:))
999	ENDWHERE
1000	!
1001	CALL flinget_buffer (force_id,'neighboursSE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1002	CALL forcing_zoom(data_full, tmpdata)
1003	WHERE(tmpdata(:,:) < undef_int)
1004	fneighbours(:,:,4) = NINT(tmpdata(:,:))
1005	ENDWHERE
1006	!
1007	CALL flinget_buffer (force_id,'neighboursSS',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1008	CALL forcing_zoom(data_full, tmpdata)
1009	WHERE(tmpdata(:,:) < undef_int)
1010	fneighbours(:,:,5) = NINT(tmpdata(:,:))
1011	ENDWHERE
1012	!
1013	CALL flinget_buffer (force_id,'neighboursSW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1014	CALL forcing_zoom(data_full, tmpdata)
1015	WHERE(tmpdata(:,:) < undef_int)
1016	fneighbours(:,:,6) = NINT(tmpdata(:,:))
1017	ENDWHERE
1018	!
1019	CALL flinget_buffer (force_id,'neighboursWW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1020	CALL forcing_zoom(data_full, tmpdata)
1021	WHERE(tmpdata(:,:) < undef_int)
1022	fneighbours(:,:,7) = NINT(tmpdata(:,:))
1023	ENDWHERE
1024	!
1025	CALL flinget_buffer (force_id,'neighboursNW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1026	CALL forcing_zoom(data_full, tmpdata)
1027	WHERE(tmpdata(:,:) < undef_int)
1028	fneighbours(:,:,8) = NINT(tmpdata(:,:))
1029	ENDWHERE
1030	!
1031	! now, resolution of the grid
1032	CALL flinget_buffer (force_id,'resolutionX',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1033	CALL forcing_zoom(data_full, tmpdata)
1034	fresolution(:,:,1) = tmpdata(:,:)
1035	!
1036	CALL flinget_buffer (force_id,'resolutionY',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1037	CALL forcing_zoom(data_full, tmpdata)
1038	fresolution(:,:,2) = tmpdata(:,:)
1039	!
1040	DO ik=1,nbindex
1041	!
1042	j = ((kindex(ik)-1)/iim) + 1
1043	i = (kindex(ik) - (j-1)*iim)
1044	!-
1045	!- Store variable to help describe the grid
1046	!- once the points are gathered.
1047	!-
1048	neighbours(ik,:) = fneighbours(i,j,:)
1049	!
1050	resolution(ik,:) = fresolution(i,j,:)
1051	!
1052
1053	ENDDO
1054	CALL gather(neighbours,neighbours_g)
1055	CALL gather(resolution,resolution_g)
1056	min_resol(1) = MINVAL(resolution(:,1))
1057	min_resol(2) = MAXVAL(resolution(:,2))
1058	max_resol(1) = MAXVAL(resolution(:,1))
1059	max_resol(2) = MAXVAL(resolution(:,2))
1060	!
1061	area(:) = resolution(:,1)*resolution(:,2)
1062	CALL gather(area,area_g)
1063	!--
1064	DEALLOCATE (tmpdata)
1065	ENDIF
1066	WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
1067	!---
1068	ENDIF
1069	!---
1070	IF (check) THEN
1071	WRITE(numout,*) &
1072	& 'The dates : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1073	ENDIF
1074	!---
1075	!--- Here we do the work in case only interpolation is needed.
1076	!---
1077	IF ( initialized .AND. interpol ) THEN
1078	!---
1079	IF ( daily_interpol ) THEN
1080
1081	IF (split > 1) THEN
1082	IF ( itau_split <= (split/2.) ) THEN
1083	rw = REAL(itau_split+split/2.)/split
1084	ELSE
1085	rw = REAL(itau_split-split/2.)/split
1086	ENDIF
1087	ELSE
1088	rw = 1.
1089	ENDIF
1090
1091	IF ((last_read == 0) .OR. ( rw==(1./split)) ) THEN
1092	!---
1093	!----- Start or Restart
1094	IF (last_read == 0) THEN
1095	! Case of a restart or a shift in the forcing file.
1096	IF (itau_read > 1) THEN
1097	itau_read_nm1=itau_read-1
1098	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1099	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1100	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1101	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1102	& force_id, wind_N_exists, check)
1103	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1104	! Case of a simple start.
1105	ELSE
1106	itau_read_nm1 = un
1107	WRITE(numout,*) "we will use the forcing of the first day to initialize "
1108	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1109	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1110	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1111	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1112	& force_id, wind_N_exists, check)
1113	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1114	ENDIF
1115	tmin_nm2(:,:)=tmin_nm1(:,:)
1116	IF ( dt_force .GT. 3600. ) THEN
1117	mean_coszang(:,:) = 0.0
1118	daylength_n(:,:) = 0.
1119	DO is=1,split
1120	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1121	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1122	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1123	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1124	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1125	WHERE( coszang(:,:) > 0. )
1126	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1127	ENDWHERE
1128	ENDDO
1129	mean_coszang(:,:) = mean_coszang(:,:)/split
1130	daylength_nm1(:,:)=daylength_n(:,:)
1131	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1132	ENDIF
1133	ELSE
1134	!----- Normal mode : copy old step
1135	swdown_nm1(:,:) = swdown_n(:,:)
1136	rainf_nm1(:,:) = rainf_n(:,:)
1137	snowf_nm1(:,:) = snowf_n(:,:)
1138	tair_nm1(:,:) = tair_n(:,:)
1139	u_nm1(:,:) = u_n(:,:)
1140	v_nm1(:,:) = v_n(:,:)
1141	qair_nm1(:,:) = qair_n(:,:)
1142	pb_nm1(:,:) = pb_n(:,:)
1143	lwdown_nm1(:,:) = lwdown_n(:,:)
1144	tmin_nm2(:,:) = tmin_nm1(:,:)
1145	tmin_nm1(:,:) = tmin_n(:,:)
1146	tmax_nm1(:,:) = tmax_n(:,:)
1147
1148	IF (is_watchout) THEN
1149	zlev_nm1(:,:) = zlev_n(:,:)
1150	zlevuv_nm1(:,:) = zlevuv_n(:,:)
1151	! Net surface short-wave flux
1152	SWnet_nm1(:,:) = SWnet_n(:,:)
1153	! Air potential energy
1154	Eair_nm1(:,:) = Eair_n(:,:)
1155	! Coeficients A from the PBL resolution for T
1156	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1157	! Coeficients A from the PBL resolution for q
1158	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1159	! Coeficients B from the PBL resolution for T
1160	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1161	! Coeficients B from the PBL resolution for q
1162	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1163	! Surface drag
1164	cdrag_nm1(:,:) = cdrag_n(:,:)
1165	! CO2 concentration in the canopy
1166	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1167	ENDIF
1168	itau_read_nm1 = itau_read_n
1169	ENDIF
1170	!-----
1171	!-----
1172	IF(last_read==0)THEN
1173	itau_read_n = itau_read
1174	ELSE
1175	itau_read_n = itau_read+1
1176	ENDIF
1177
1178	IF (itau_read_n > ttm) THEN
1179	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1180	WRITE(numout,*) &
1181	& 'WARNING : We are going back to the start of the file'
1182	itau_read_n =1
1183	ENDIF
1184	IF (check) THEN
1185	WRITE(numout,*) &
1186	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1187	ENDIF
1188	!-----
1189	!----- Get a reduced julian day !
1190	!----- This is needed because we lack the precision on 32 bit machines.
1191	!-----
1192	IF ( dt_force .GT. 3600. ) THEN
1193	julian_for = itau2date(itau_read-1, date0, dt_force)
1194	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1195
1196	! first day of this year
1197	CALL ymds2ju (yy,1,1,0.0, julian0)
1198	!-----
1199	IF (check) THEN
1200	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1201	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1202	ENDIF
1203	ENDIF
1204	!-----
1205	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1206	& swdown_n, rainf_n, snowf_n, tmin_n, &
1207	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1208	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1209	& force_id, wind_N_exists, check)
1210	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_n, itau_read_n, tmax_n, force_id )
1211
1212	!---
1213	last_read = itau_read_n
1214	!-----
1215	!----- Compute mean solar angle for the comming period
1216	!-----
1217	IF (check) WRITE(numout,*) 'Going into solarang', split, one_day
1218	!-----
1219
1220	!-----
1221	ENDIF
1222	!---
1223	IF ( itau_split == 1. ) THEN
1224	IF ( dt_force .GT. 3600. ) THEN
1225	mean_coszang(:,:) = 0.0
1226	daylength_nm1(:,:)=daylength_n(:,:)
1227	daylength_n(:,:) = 0.
1228	DO is=1,split
1229	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1230	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1231	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1232	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1233	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1234	WHERE( coszang(:,:) > 0. )
1235	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1236	ENDWHERE
1237	ENDDO
1238	mean_coszang(:,:) = mean_coszang(:,:)/split
1239	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1240	ENDIF
1241	ENDIF
1242
1243	!--- Do the interpolation
1244	IF (check) WRITE(numout,*) 'Doing the interpolation between time steps'
1245	!---
1246
1247	IF (check) WRITE(numout,*) 'Coeff of interpollation : ',rw
1248	!---
1249
1250	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1251	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1252	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1253
1254	!--- Take care of the height of the vertical levels
1255	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1256	zlevuv(:,:) = (zlevuv_n(:,:)-zlevuv_nm1(:,:))*rw + zlevuv_nm1(:,:)
1257
1258	hour=REAL(itau_split)/split*24
1259	startday_n(:,:)=12.-daylength_n(:,:)/2.
1260	startday_nm1(:,:)=12.-daylength_nm1(:,:)/2.
1261
1262	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1263	tair(:,:)=(tmax_nm1(:,:)-tmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1264	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_nm1(:,:)+tmin_nm1(:,:))/2.
1265	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1266	tair(:,:)=(tmax_n(:,:)-tmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1267	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_n(:,:)+tmin_n(:,:))/2.
1268	ELSEWHERE ( hour < startday_n(:,:) )
1269	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1270	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1271	ELSEWHERE
1272	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1273	& (24.-14.+startday_n(:,:))-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1274	ENDWHERE
1275
1276	CALL weathgen_qsat_2d (iim,jjm,tmin_n,pb,qsattmin_n)
1277	CALL weathgen_qsat_2d (iim,jjm,tmin_nm1,pb,qsattmin_nm1)
1278	CALL weathgen_qsat_2d (iim,jjm,tair,pb,qsatta)
1279
1280	!---
1281	qmin_nm1(:,:) = MIN(qair_nm1(:,:),0.99*qsattmin_nm1(:,:))
1282	qmin_n(:,:) = MIN(qair_n(:,:),0.99*qsattmin_n(:,:))
1283	qmax_nm1(:,:) = (qair_nm1(:,:)-qmin_nm1(:,:)) + qair_nm1(:,:)
1284	qmax_n(:,:) = (qair_n(:,:)-qmin_n(:,:)) + qair_n(:,:)
1285
1286	qsa(:,:) = 0.99*qsatta(:,:)
1287
1288
1289	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1290	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1291	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_nm1(:,:)+qmin_nm1(:,:))/2.)
1292	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1293	qair(:,:)=MIN(qsa(:,:),(qmax_n(:,:)-qmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1294	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_n(:,:)+qmin_n(:,:))/2.)
1295	ELSEWHERE ( hour < startday_n(:,:) )
1296	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1297	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1298	ELSEWHERE
1299	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1300	& (24.-14.+startday_n(:,:))-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1301	ENDWHERE
1302
1303	IF (is_watchout) THEN
1304	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1305	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1306	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1307	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1308	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1309	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1310	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1311	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1312	ENDIF
1313	!---
1314	!--- Here we need to allow for an option
1315	!--- where radiative energy is conserved
1316	!---
1317	IF ( netrad_cons ) THEN
1318	lwdown(:,:) = lwdown_n(:,:)
1319	ELSE
1320	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1321	ENDIF
1322	!---
1323	!--- For the solar radiation we decompose the mean value
1324	!--- using the zenith angle of the sun if the time step in the forcing data is
1325	!---- more than an hour. Else we use the standard linera interpolation
1326	!----
1327	IF (check) WRITE(numout,*) 'Ready to deal with the solar radiation'
1328	!----
1329	IF ( dt_force .GT. 3600. ) THEN
1330	!---
1331	IF ( netrad_cons ) THEN
1332	WRITE(numout,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
1333	ENDIF
1334	!---
1335	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1336	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1337	!!$ julian = julian_for + rw*dt_force/one_day
1338	IF (check) THEN
1339	WRITE(numout,'(a,f20.10,2I3)') &
1340	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1341	ENDIF
1342	!---
1343	CALL solarang(julian, julian0, iim, jjm, lon*0, lat, coszang)
1344	!---
1345
1346	WHERE ((mean_coszang(:,:) > 0.) .AND. (hour <= 12 ))
1347	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1348	ELSEWHERE ((mean_coszang(:,:) > 0.) .AND. (hour > 12 ))
1349	swdown(:,:) = swdown_nm1(:,:) *coszang(:,:)/mean_coszang(:,:)
1350	ELSEWHERE
1351	swdown(:,:) = 0.0
1352	END WHERE
1353	!---
1354	WHERE (swdown(:,:) > 2000. )
1355	swdown(:,:) = 2000.
1356	END WHERE
1357	!---
1358	ELSE
1359	!---
1360	IF ( netrad_cons ) THEN
1361	swdown(:,:) = swdown_n(:,:)
1362	ELSE
1363	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1364	ENDIF
1365	!---
1366	ENDIF
1367	!---
1368	IF (check) THEN
1369	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1370	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1371	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1372	IF (is_watchout) THEN
1373	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1374	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1375	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1376	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1377	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1378	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1379	ENDIF
1380	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1381	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1382	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1383	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1384	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1385	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1386	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1387	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1388	ENDIF
1389	!---
1390	!--- For precip we suppose that the rain
1391	!--- is the sum over the next 6 hours
1392	!---
1393	WHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)>=273.15))
1394	rainf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1395	snowf(:,:) = 0.0
1396	ELSEWHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)<273.15))
1397	snowf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1398	rainf(:,:) = 0.0
1399	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)>=273.15))
1400	rainf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1401	snowf(:,:) = 0.0
1402	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)<273.15))
1403	snowf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1404	rainf(:,:) = 0.0
1405	ELSEWHERE
1406	snowf(:,:) = 0.0
1407	rainf(:,:) = 0.0
1408	ENDWHERE
1409
1410	IF (check) THEN
1411	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1412	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1413	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1414	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1415	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1416	ENDIF
1417	!---
1418
1419
1420	ELSE
1421
1422	IF (itau_read /= last_read) THEN
1423	!---
1424	!----- Start or Restart
1425	IF (itau_read_n == 0) THEN
1426	! Case of a restart or a shift in the forcing file.
1427	IF (itau_read > 1) THEN
1428	itau_read_nm1=itau_read-1
1429	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1430	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1431	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1432	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1433	& force_id, wind_N_exists, check)
1434	! Case of a simple start.
1435	ELSE IF (dt_force*ttm > one_day-1. ) THEN
1436	! if the forcing file contains at least 24 hours,
1437	! we will use the last forcing step of the first day
1438	! as initiale condition to prevent first shift off reading.
1439	itau_read_nm1 = NINT (one_day/dt_force)
1440	WRITE(numout,) "the forcing file contains 24 hours :",dt_forcettm,one_day-1.
1441	WRITE(numout,*) "we will use the last forcing step of the first day : itau_read_nm1 ",itau_read_nm1
1442	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1443	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1444	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1445	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1446	& force_id, wind_N_exists, check)
1447	ELSE
1448	! if the forcing file contains less than 24 hours,
1449	! just say error !
1450	CALL ipslerr_p(3,'forcing_read_interpol', &
1451	& 'The forcing file contains less than 24 hours !', &
1452	& 'We can''t intialize interpolation with such a file.','')
1453	ENDIF
1454	ELSE
1455	!----- Normal mode : copy old step
1456	swdown_nm1(:,:) = swdown_n(:,:)
1457	rainf_nm1(:,:) = rainf_n(:,:)
1458	snowf_nm1(:,:) = snowf_n(:,:)
1459	tair_nm1(:,:) = tair_n(:,:)
1460	u_nm1(:,:) = u_n(:,:)
1461	v_nm1(:,:) = v_n(:,:)
1462	qair_nm1(:,:) = qair_n(:,:)
1463	pb_nm1(:,:) = pb_n(:,:)
1464	lwdown_nm1(:,:) = lwdown_n(:,:)
1465	IF (is_watchout) THEN
1466	zlev_nm1(:,:) = zlev_n(:,:)
1467	! Net surface short-wave flux
1468	SWnet_nm1(:,:) = SWnet_n(:,:)
1469	! Air potential energy
1470	Eair_nm1(:,:) = Eair_n(:,:)
1471	! Coeficients A from the PBL resolution for T
1472	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1473	! Coeficients A from the PBL resolution for q
1474	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1475	! Coeficients B from the PBL resolution for T
1476	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1477	! Coeficients B from the PBL resolution for q
1478	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1479	! Surface drag
1480	cdrag_nm1(:,:) = cdrag_n(:,:)
1481	! CO2 concentration in the canopy
1482	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1483	ENDIF
1484	itau_read_nm1 = itau_read_n
1485	ENDIF
1486	!-----
1487	itau_read_n = itau_read
1488	IF (itau_read_n > ttm) THEN
1489	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1490	WRITE(numout,*) &
1491	& 'WARNING : We are going back to the start of the file'
1492	itau_read_n =1
1493	ENDIF
1494	IF (check) THEN
1495	WRITE(numout,*) &
1496	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1497	ENDIF
1498	!-----
1499	!----- Get a reduced julian day !
1500	!----- This is needed because we lack the precision on 32 bit machines.
1501	!-----
1502	IF ( dt_force .GT. 3600. ) THEN
1503	julian_for = itau2date(itau_read-1, date0, dt_force)
1504	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1505
1506	! first day of this year
1507	CALL ymds2ju (yy,1,1,0.0, julian0)
1508	!-----
1509	IF (check) THEN
1510	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1511	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1512	ENDIF
1513	ENDIF
1514	!-----
1515	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1516	& swdown_n, rainf_n, snowf_n, tair_n, &
1517	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1518	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1519	& force_id, wind_N_exists, check)
1520	!---
1521	last_read = itau_read_n
1522	!-----
1523	!----- Compute mean solar angle for the comming period
1524	!-----
1525	IF (check) WRITE(numout,*) 'Going into solarang', split, one_day
1526	!-----
1527	IF ( dt_force .GT. 3600. ) THEN
1528	mean_coszang(:,:) = 0.0
1529	DO is=1,split
1530	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1531	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1532	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1533	CALL solarang (julian, julian0, iim, jjm, lon, lat, coszang)
1534	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1535	ENDDO
1536	mean_coszang(:,:) = mean_coszang(:,:)/split
1537	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1538	ENDIF
1539	!-----
1540	ENDIF
1541	!---
1542	!--- Do the interpolation
1543	IF (check) WRITE(numout,*) 'Doing the interpolation between time steps'
1544	!---
1545	IF (split > 1) THEN
1546	rw = REAL(itau_split)/split
1547	ELSE
1548	rw = 1.
1549	ENDIF
1550	IF (check) WRITE(numout,*) 'Coeff of interpollation : ',rw
1551	!---
1552	qair(:,:) = (qair_n(:,:)-qair_nm1(:,:))*rw + qair_nm1(:,:)
1553	tair(:,:) = (tair_n(:,:)-tair_nm1(:,:))*rw + tair_nm1(:,:)
1554	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1555	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1556	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1557	IF (is_watchout) THEN
1558	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1559	zlevuv(:,:) = zlev(:,:)
1560	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1561	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1562	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1563	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1564	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1565	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1566	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1567	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1568	ENDIF
1569	!---
1570	!--- Here we need to allow for an option
1571	!--- where radiative energy is conserved
1572	!---
1573	IF ( netrad_cons ) THEN
1574	lwdown(:,:) = lwdown_n(:,:)
1575	ELSE
1576	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1577	ENDIF
1578	!---
1579	!--- For the solar radiation we decompose the mean value
1580	!--- using the zenith angle of the sun if the time step in the forcing data is
1581	!---- more than an hour. Else we use the standard linera interpolation
1582	!----
1583	IF (check) WRITE(numout,*) 'Ready to deal with the solar radiation'
1584	!----
1585	IF ( dt_force .GT. 3600. ) THEN
1586	!---
1587	IF ( netrad_cons ) THEN
1588	WRITE(numout,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
1589	ENDIF
1590	!---
1591	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1592	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1593	!!$ julian = julian_for + rw*dt_force/one_day
1594	IF (check) THEN
1595	WRITE(numout,'(a,f20.10,2I3)') &
1596	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1597	ENDIF
1598	!---
1599	CALL solarang(julian, julian0, iim, jjm, lon, lat, coszang)
1600	!---
1601	WHERE (mean_coszang(:,:) > 0.)
1602	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1603	ELSEWHERE
1604	swdown(:,:) = 0.0
1605	END WHERE
1606	!---
1607	WHERE (swdown(:,:) > 2000. )
1608	swdown(:,:) = 2000.
1609	END WHERE
1610	!---
1611	ELSE
1612	!---
1613	IF ( netrad_cons ) THEN
1614	swdown(:,:) = swdown_n(:,:)
1615	ELSE
1616	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1617	ENDIF
1618	!---
1619	ENDIF
1620	!---
1621	IF (check) THEN
1622	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1623	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1624	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1625	IF (is_watchout) THEN
1626	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1627	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1628	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1629	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1630	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1631	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1632	ENDIF
1633	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1634	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1635	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1636	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1637	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1638	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1639	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1640	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1641	ENDIF
1642	!---
1643	!--- For precip we suppose that the rain
1644	!--- is the sum over the next 6 hours
1645	!---
1646	IF (itau_split <= nb_spread) THEN
1647	rainf(:,:) = rainf_n(:,:)*(split/REAL(nb_spread))
1648	snowf(:,:) = snowf_n(:,:)*(split/REAL(nb_spread))
1649	ELSE
1650	rainf(:,:) = 0.0
1651	snowf(:,:) = 0.0
1652	ENDIF
1653	IF (check) THEN
1654	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1655	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1656	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1657	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1658	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1659	ENDIF
1660	!---
1661	ENDIF ! (daily_interpol)
1662	ENDIF
1663	!---
1664	!--- Here we might put the call to the weather generator ... one day.
1665	!--- Pour le moment, le branchement entre interpolation et generateur de temps
1666	!--- est fait au-dessus.
1667	!---
1668	!- IF ( initialized .AND. weathergen ) THEN
1669	!- ....
1670	!- ENDIF
1671	!---
1672	!--- At this point the code should be initialized. If not we have a problem !
1673	!---
1674	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
1675	!---
1676	initialized = .TRUE.
1677	!---
1678	ELSE
1679	IF ( .NOT. initialized ) THEN
1680	WRITE(numout,*) 'Why is the code forcing_read not initialized ?'
1681	WRITE(numout,*) 'Have you called it with both time-steps set to zero ?'
1682	CALL ipslerr_p(3,'forcing_read_interpol','Pb in initialization','','')
1683	ENDIF
1684	ENDIF
1685	!
1686	!-------------------------
1687	END SUBROUTINE forcing_read_interpol
1688	!=====================================================================
1689	!-
1690	!=====================================================================
1691	SUBROUTINE forcing_just_read &
1692	& (iim, jjm, zlev, zlev_uv, ttm, itb, ite, &
1693	& swdown, rainf, snowf, tair, &
1694	& u, v, qair, pb, lwdown, &
1695	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
1696	& force_id, wind_N_exists, check)
1697	!---------------------------------------------------------------------
1698	!- iim : Size of the grid in x
1699	!- jjm : size of the grid in y
1700	!- zlev : height of the varibales T and Q
1701	!- zlev_uv : height of the varibales U and V
1702	!- ttm : number of time steps in all in the forcing file
1703	!- itb, ite : index of respectively begin and end of read for each variable
1704	!- swdown : Downward solar radiation (W/m^2)
1705	!- rainf : Rainfall (kg/m^2s)
1706	!- snowf : Snowfall (kg/m^2s)
1707	!- tair : 2m air temperature (K)
1708	!- u and v : 2m (in theory !) wind speed (m/s)
1709	!- qair : 2m humidity (kg/kg)
1710	!- pb : Surface pressure (Pa)
1711	!- lwdown : Downward long wave radiation (W/m^2)
1712	!-
1713	!- From a WATCHOUT file :
1714	!- SWnet : Net surface short-wave flux
1715	!- Eair : Air potential energy
1716	!- petAcoef : Coeficients A from the PBL resolution for T
1717	!- peqAcoef : Coeficients A from the PBL resolution for q
1718	!- petBcoef : Coeficients B from the PBL resolution for T
1719	!- peqBcoef : Coeficients B from the PBL resolution for q
1720	!- cdrag : Surface drag
1721	!- ccanopy : CO2 concentration in the canopy
1722	!- force_id : FLINCOM file id.
1723	!- It is used to close the file at the end of the run.
1724	!- wind_N_exists : if Wind_N and Wind_E are in the file (and not just Wind)
1725	!- check : Prompt for reading
1726	!---------------------------------------------------------------------
1727	IMPLICIT NONE
1728	!-
1729	INTEGER, INTENT(in) :: iim, jjm, ttm
1730	INTEGER, INTENT(in) :: itb, ite
1731	REAL, DIMENSION(iim,jjm), INTENT(out) :: zlev, zlev_uv, &
1732	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown
1733	! for watchout files
1734	REAL, DIMENSION(iim,jjm), INTENT(out) :: &
1735	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
1736	INTEGER, INTENT(in) :: force_id
1737	! if Wind_N and Wind_E are in the file (and not just Wind)
1738	LOGICAL, INTENT(in) :: wind_N_exists
1739	LOGICAL :: check
1740	INTEGER :: i, j
1741	REAL :: rau
1742
1743	!-
1744	!---------------------------------------------------------------------
1745	IF ( daily_interpol ) THEN
1746	CALL flinget_buffer (force_id,'Tmin' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1747	CALL forcing_zoom(data_full, tair)
1748	CALL flinget_buffer (force_id,'precip' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1749	CALL forcing_zoom(data_full, rainf)
1750	ELSE
1751	CALL flinget_buffer (force_id,'Tair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1752	CALL forcing_zoom(data_full, tair)
1753	CALL flinget_buffer (force_id,'Snowf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1754	CALL forcing_zoom(data_full, snowf)
1755	CALL flinget_buffer (force_id,'Rainf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1756	CALL forcing_zoom(data_full, rainf)
1757	ENDIF
1758
1759
1760	CALL flinget_buffer (force_id,'SWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1761	CALL forcing_zoom(data_full, swdown)
1762	CALL flinget_buffer (force_id,'LWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1763	CALL forcing_zoom(data_full, lwdown)
1764
1765	CALL flinget_buffer (force_id,'PSurf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1766	CALL forcing_zoom(data_full, pb)
1767	CALL flinget_buffer (force_id,'Qair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1768	CALL forcing_zoom(data_full, qair)
1769	!---
1770	IF ( wind_N_exists ) THEN
1771	CALL flinget_buffer (force_id,'Wind_N', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1772	CALL forcing_zoom(data_full, u)
1773	CALL flinget_buffer (force_id,'Wind_E', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1774	CALL forcing_zoom(data_full, v)
1775	ELSE
1776	CALL flinget_buffer (force_id,'Wind', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1777	CALL forcing_zoom(data_full, u)
1778	v=0.0
1779	ENDIF
1780
1781	!-
1782	!- Deal with the height of the atmospheric forcing varibles
1783	!-
1784	!----
1785	IF ( zheight ) THEN
1786	zlev(:,:) = zlev_fixed
1787	ELSE IF ( zsigma .OR. zhybrid ) THEN
1788	DO i=1,iim
1789	DO j=1,jjm
1790	IF ( tair(i,j) < val_exp ) THEN
1791	rau = pb(i,j)/(cte_molr*tair(i,j))
1792
1793	zlev(i,j) = (pb(i,j) - (zhybrid_a + zhybrid_bpb(i,j)))/(rau cte_grav)
1794	ELSE
1795	zlev(i,j) = 0.0
1796	ENDIF
1797	ENDDO
1798	ENDDO
1799	ELSE IF ( zlevels ) THEN
1800	CALL flinget_buffer (force_id,'Levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1801	CALL forcing_zoom(data_full, zlev)
1802	ELSE
1803	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1804	& 'We cannot determine the height for T and Q.','stop readdim2')
1805	ENDIF
1806
1807	IF ( zsamelev_uv ) THEN
1808	zlev_uv(:,:) = zlev(:,:)
1809	ELSE
1810	IF ( zheight ) THEN
1811	zlev_uv(:,:) = zlevuv_fixed
1812	ELSE IF ( zsigma .OR. zhybrid ) THEN
1813	DO i=1,iim
1814	DO j=1,jjm
1815	IF ( tair(i,j) < val_exp ) THEN
1816	rau = pb(i,j)/(cte_molr*tair(i,j))
1817
1818	zlev_uv(i,j) = (pb(i,j) - (zhybriduv_a + zhybriduv_bpb(i,j)))/(rau cte_grav)
1819	ELSE
1820	zlev_uv(i,j) = 0.0
1821	ENDIF
1822	ENDDO
1823	ENDDO
1824	ELSE IF ( zlevels ) THEN
1825	CALL flinget_buffer (force_id,'Levels_uv', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1826	CALL forcing_zoom(data_full, zlev_uv)
1827	ELSE
1828	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1829	& 'We cannot determine the height for U and V.','stop readdim2')
1830	ENDIF
1831	ENDIF
1832	!----
1833	IF ( is_watchout ) THEN
1834	CALL flinget_buffer (force_id,'levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1835	CALL forcing_zoom(data_full, zlev)
1836	!
1837	! If we are in WATHCOUT it means T,Q are at the same height as U,V
1838	!
1839	zlev_uv(:,:) = zlev(:,:)
1840	! Net surface short-wave flux
1841	CALL flinget_buffer (force_id,'SWnet', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1842	CALL forcing_zoom(data_full, SWnet)
1843	! Air potential energy
1844	CALL flinget_buffer (force_id,'Eair', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1845	CALL forcing_zoom(data_full, Eair)
1846	! Coeficients A from the PBL resolution for T
1847	CALL flinget_buffer (force_id,'petAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1848	CALL forcing_zoom(data_full, petAcoef)
1849	! Coeficients A from the PBL resolution for q
1850	CALL flinget_buffer (force_id,'peqAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1851	CALL forcing_zoom(data_full, peqAcoef)
1852	! Coeficients B from the PBL resolution for T
1853	CALL flinget_buffer (force_id,'petBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1854	CALL forcing_zoom(data_full, petBcoef)
1855	! Coeficients B from the PBL resolution for q
1856	CALL flinget_buffer (force_id,'peqBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1857	CALL forcing_zoom(data_full, peqBcoef)
1858	! Surface drag
1859	CALL flinget_buffer (force_id,'cdrag', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1860	CALL forcing_zoom(data_full, cdrag)
1861	! CO2 concentration in the canopy
1862	CALL flinget_buffer (force_id,'ccanopy', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1863	CALL forcing_zoom(data_full, ccanopy)
1864	ENDIF
1865	!
1866	!----
1867	IF (check) WRITE(numout,*) 'Variables have been extracted between ',itb,' and ',ite,' iterations of the forcing file.'
1868	!-------------------------
1869	END SUBROUTINE forcing_just_read
1870	!=====================================================================
1871
1872	!-
1873	SUBROUTINE forcing_just_read_tmax &
1874	& (iim, jjm, ttm, itb, ite, tmax, force_id )
1875	!---------------------------------------------------------------------
1876	!- iim : Size of the grid in x
1877	!- jjm : size of the grid in y
1878	!- ttm : number of time steps in all in the forcing file
1879	!- itb, ite : index of respectively begin and end of read for each variable
1880	!- tmax : 2m air temperature (K)
1881	!- force_id : FLINCOM file id.
1882	!- It is used to close the file at the end of the run.
1883	!---------------------------------------------------------------------
1884	IMPLICIT NONE
1885	!-
1886	INTEGER, INTENT(in) :: iim, jjm, ttm
1887	INTEGER, INTENT(in) :: itb, ite
1888	REAL, DIMENSION(iim,jjm), INTENT(out) :: tmax
1889	INTEGER, INTENT(in) :: force_id
1890	!-
1891	!---------------------------------------------------------------------
1892	CALL flinget_buffer (force_id,'Tmax' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1893	CALL forcing_zoom(data_full, tmax)
1894	!-------------------------
1895	END SUBROUTINE forcing_just_read_tmax
1896	!=====================================================================
1897
1898	!-
1899	SUBROUTINE forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
1900	!---
1901	!--- This subroutine finds the indices of the land points over which the land
1902	!--- surface scheme is going to run.
1903	!---
1904	IMPLICIT NONE
1905	!-
1906	!- ARGUMENTS
1907	!-
1908	INTEGER, INTENT(IN) :: iim, jjm
1909	REAL, INTENT(IN) :: tair(iim,jjm)
1910	INTEGER, INTENT(OUT) :: i_test, j_test, nbindex
1911	INTEGER, INTENT(OUT) :: kindex(iim*jjm)
1912	LOGICAL :: check
1913	!-
1914	!- LOCAL
1915	INTEGER :: i, j, ig
1916	!-
1917	!-
1918	ig = 0
1919	i_test = 0
1920	j_test = 0
1921	!---
1922	IF (MINVAL(tair(:,:)) < 100.) THEN
1923	!----- In this case the 2m temperature is in Celsius
1924	DO j=1,jjm
1925	DO i=1,iim
1926	IF (tair(i,j) < 100.) THEN
1927	ig = ig+1
1928	kindex(ig) = (j-1)*iim+i
1929	!
1930	! Here we find at random a land-point on which we can do
1931	! some printouts for test.
1932	!
1933	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1934	i_test = i
1935	j_test = j
1936	IF (check) THEN
1937	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1938	ENDIF
1939	ENDIF
1940	ENDIF
1941	ENDDO
1942	ENDDO
1943	ELSE
1944	!----- 2m temperature is in Kelvin
1945	DO j=1,jjm
1946	DO i=1,iim
1947	IF (tair(i,j) < 500.) THEN
1948	ig = ig+1
1949	kindex(ig) = (j-1)*iim+i
1950	!
1951	! Here we find at random a land-point on which we can do
1952	! some printouts for test.
1953	!
1954	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1955	i_test = i
1956	j_test = j
1957	IF (check) THEN
1958	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1959	ENDIF
1960	ENDIF
1961	ENDIF
1962	ENDDO
1963	ENDDO
1964	ENDIF
1965	!---
1966	nbindex = ig
1967	!---
1968	END SUBROUTINE forcing_landind
1969	!-
1970	!=====================================================================
1971	!-
1972	SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,init_f)
1973	!-
1974	!- This subroutine calculates the longitudes and latitudes of the model grid.
1975	!-
1976	IMPLICIT NONE
1977	!-
1978	INTEGER, INTENT(in) :: iim, jjm, llm
1979	LOGICAL, INTENT(in) :: init_f
1980	REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat
1981	!-
1982	INTEGER :: i,j
1983	!-
1984	!- Should be unified one day
1985	!-
1986	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : options : ', weathergen, interpol
1987	!-
1988	IF ( weathergen ) THEN
1989	IF (init_f) THEN
1990	DO i = 1, iim
1991	lon(i,:) = limit_west + merid_res/2. + &
1992	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim)
1993	ENDDO
1994	!-
1995	DO j = 1, jjm
1996	lat(:,j) = limit_north - zonal_res/2. - &
1997	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm)
1998	ENDDO
1999	ELSE
2000	IF (is_root_prc) THEN
2001	DO i = 1, iim_g
2002	lon_g(i,:) = limit_west + merid_res/2. + &
2003	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim_g)
2004	ENDDO
2005	!-
2006	DO j = 1, jjm_g
2007	lat_g(:,j) = limit_north - zonal_res/2. - &
2008	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm_g)
2009	ENDDO
2010	ENDIF
2011	CALL bcast(lon_g)
2012	CALL bcast(lat_g)
2013	lon=lon_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2014	lat=lat_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2015	ENDIF
2016	!-
2017	ELSEIF ( interpol ) THEN
2018	!-
2019	CALL forcing_zoom(lon_full, lon)
2020	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lon'
2021	CALL forcing_zoom(lat_full, lat)
2022	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lat'
2023
2024	ELSE
2025	CALL ipslerr_p(3,'forcing_grid','Neither interpolation nor weather generator is specified.','','')
2026	ENDIF
2027
2028	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : ended'
2029
2030	END SUBROUTINE forcing_grid
2031	!-
2032	!=====================================================================
2033	!-
2034	SUBROUTINE forcing_zoom(x_f, x_z)
2035	!-
2036	!- This subroutine takes the slab of data over which we wish to run the model.
2037	!-
2038	IMPLICIT NONE
2039	!-
2040	REAL, INTENT(IN) :: x_f(iim_full, jjm_full)
2041	REAL, INTENT(OUT) :: x_z(iim_zoom, jjm_zoom)
2042	!-
2043	INTEGER :: i,j
2044	!-
2045	DO i=1,iim_zoom
2046	DO j=1,jjm_zoom
2047	x_z(i,j) = x_f(i_index(i),j_index(j))
2048	ENDDO
2049	ENDDO
2050	!-
2051	END SUBROUTINE forcing_zoom
2052
2053	!-
2054	!=====================================================================
2055	!-
2056
2057	SUBROUTINE forcing_vertical(force_id)
2058	!
2059	!- This subroutine explores the forcing file to decide what information is available
2060	!- on the vertical coordinate.
2061	!
2062	INTEGER, INTENT(IN) :: force_id
2063	!
2064	LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
2065	LOGICAL :: foundvar
2066	LOGICAL :: levlegacy
2067
2068	!
2069	!- Set all the defaults
2070	!
2071	zfixed=.FALSE.
2072	zsigma=.FALSE.
2073	zhybrid=.FALSE.
2074	zlevels=.FALSE.
2075	zheight=.FALSE.
2076	zsamelev_uv = .TRUE.
2077	levlegacy = .FALSE.
2078	!
2079	foundvar = .FALSE.
2080	!
2081	!- We have a forcing file to explore so let us see if we find any of the conventions
2082	!- which allow us to find the height of T,Q,U and V.
2083	!
2084	IF ( force_id > 0 ) THEN
2085	!
2086	! Case for sigma levels
2087	!
2088	IF ( .NOT. foundvar ) THEN
2089	CALL flinquery_var(force_id, 'Sigma', var_exists)
2090	CALL flinquery_var(force_id, 'Sigma_uv', varuv_exists)
2091	IF ( var_exists ) THEN
2092	foundvar = .TRUE.
2093	zsigma = .TRUE.
2094	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2095	ENDIF
2096	ENDIF
2097	!
2098	! Case for Hybrid levels
2099	!
2100	IF ( .NOT. foundvar ) THEN
2101	CALL flinquery_var(force_id, 'HybSigA', vara_exists)
2102	IF ( vara_exists ) THEN
2103	CALL flinquery_var(force_id, 'HybSigB', varb_exists)
2104	IF ( varb_exists ) THEN
2105	var_exists=.TRUE.
2106	ELSE
2107	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
2108	& 'Hybrid vertical levels for T and Q','stop readdim2')
2109	ENDIF
2110	ENDIF
2111	CALL flinquery_var(force_id, 'HybSigA_uv', vara_exists)
2112	IF ( vara_exists ) THEN
2113	CALL flinquery_var(force_id, 'HybSigB_uv', varb_exists)
2114	IF ( varb_exists ) THEN
2115	varuv_exists=.TRUE.
2116	ELSE
2117	CALL ipslerr ( 3, 'forcing_vertical','Missing the B coefficient for', &
2118	& 'Hybrid vertical levels for U and V','stop readdim2')
2119	ENDIF
2120	ENDIF
2121	IF ( var_exists ) THEN
2122	foundvar = .TRUE.
2123	zhybrid = .TRUE.
2124	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2125	ENDIF
2126	ENDIF
2127	!
2128	! Case for levels (i.e. a 2d time varying field with the height in meters)
2129	!
2130	IF ( .NOT. foundvar ) THEN
2131	CALL flinquery_var(force_id, 'Levels', var_exists)
2132	CALL flinquery_var(force_id, 'Levels_uv', varuv_exists)
2133	IF ( var_exists ) THEN
2134	foundvar = .TRUE.
2135	zlevels = .TRUE.
2136	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2137	ENDIF
2138	ENDIF
2139	!
2140	! Case where a fixed height is provided in meters
2141	!
2142	IF ( .NOT. foundvar ) THEN
2143	CALL flinquery_var(force_id, 'Height_Lev1', var_exists)
2144	CALL flinquery_var(force_id, 'Height_Levuv', varuv_exists)
2145	IF ( var_exists ) THEN
2146	foundvar = .TRUE.
2147	zheight = .TRUE.
2148	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2149	ENDIF
2150	ENDIF
2151	!
2152	! Case where a fixed height is provided in meters in the lev variable
2153	!
2154	IF ( .NOT. foundvar ) THEN
2155	CALL flinquery_var(force_id, 'lev', var_exists)
2156	IF ( var_exists ) THEN
2157	foundvar = .TRUE.
2158	zheight = .TRUE.
2159	levlegacy = .TRUE.
2160	ENDIF
2161	ENDIF
2162	!
2163	ENDIF
2164	!
2165	! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
2166	! except the case zlevels
2167	!
2168	IF ( foundvar .AND. .NOT. zlevels ) THEN
2169	!
2170	IF ( zheight ) THEN
2171	!
2172	! Constant height
2173	!
2174	IF ( levlegacy ) THEN
2175	CALL flinget (force_id,'lev',1, 1, 1, 1, 1, 1, zlev_fixed)
2176	ELSE
2177	CALL flinget (force_id,'Height_Lev1',1, 1, 1, 1, 1, 1, zlev_fixed)
2178	IF ( .NOT. zsamelev_uv ) THEN
2179	CALL flinget (force_id,'Height_Levuv',1, 1, 1, 1, 1, 1, zlevuv_fixed)
2180	ENDIF
2181	ENDIF
2182	WRITE(numout,*) "forcing_vertical : case ZLEV : Read from forcing file :", zlev_fixed, zlevuv_fixed
2183	!
2184	ELSE IF ( zsigma .OR. zhybrid ) THEN
2185	!
2186	! Sigma or hybrid levels
2187	!
2188	IF ( zsigma ) THEN
2189	CALL flinget (force_id,'Sigma',1, 1, 1, 1, 1, 1, zhybrid_b)
2190	zhybrid_a = zero
2191	IF ( .NOT. zsamelev_uv ) THEN
2192	CALL flinget (force_id,'Sigma_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2193	zhybriduv_a = zero
2194	ENDIF
2195	ELSE
2196	CALL flinget (force_id,'HybSigB',1, 1, 1, 1, 1, 1, zhybrid_b)
2197	CALL flinget (force_id,'HybSigA',1, 1, 1, 1, 1, 1, zhybrid_a)
2198	IF ( .NOT. zsamelev_uv ) THEN
2199	CALL flinget (force_id,'HybSigB_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2200	CALL flinget (force_id,'HybSigA_uv',1, 1, 1, 1, 1, 1, zhybriduv_a)
2201	ENDIF
2202	ENDIF
2203	WRITE(numout,*) "forcing_vertical : case Pressure coordinates : "
2204	WRITE(numout,*) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
2205	ELSE
2206	!
2207	! Why are we here ???
2208	!
2209	CALL ipslerr ( 3, 'forcing_vertical','What is the option used to describe the height of', &
2210	& 'the atmospheric forcing ?','Please check your forcing file.')
2211	ENDIF
2212	ENDIF
2213	!
2214	!- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
2215	!- read what has been specified by the user.
2216	!
2217	IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
2218	!
2219	!-
2220	!Config Key = HEIGHT_LEV1
2221	!Config Desc = Height at which T and Q are given
2222	!Config Def = 2.0
2223	!Config If = offline mode
2224	!Config Help = The atmospheric variables (temperature and specific
2225	!Config humidity) are measured at a specific level.
2226	!Config The height of this level is needed to compute
2227	!Config correctly the turbulent transfer coefficients.
2228	!Config Look at the description of the forcing
2229	!Config DATA for the correct value.
2230	!Config Units = [m]
2231	!-
2232	zlev_fixed = 2.0
2233	CALL getin('HEIGHT_LEV1', zlev_fixed)
2234
2235	!-
2236	!Config Key = HEIGHT_LEVW
2237	!Config Desc = Height at which the wind is given
2238	!Config Def = 10.0
2239	!Config If = offline mode
2240	!Config Help = The height at which wind is needed to compute
2241	!Config correctly the turbulent transfer coefficients.
2242	!Config Units= [m]
2243	!-
2244	zlevuv_fixed = 10.0
2245	CALL getin('HEIGHT_LEVW', zlevuv_fixed)
2246
2247	zheight = .TRUE.
2248
2249	IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
2250	zsamelev_uv = .FALSE.
2251	ENDIF
2252
2253	CALL ipslerr ( 2, 'forcing_vertical','The height of the atmospheric forcing variables', &
2254	& 'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
2255	ENDIF
2256
2257	END SUBROUTINE forcing_vertical
2258
2259	!-
2260	!=====================================================================
2261	!-
2262
2263	SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, &
2264	& iim_f, jjm_f, lon, lat, iim, jjm, iind, jind)
2265
2266	IMPLICIT NONE
2267	!
2268	! ARGUMENTS
2269	!
2270	REAL, INTENT(inout) :: limit_west,limit_east,limit_north,limit_south
2271	INTEGER, INTENT(in) :: iim_f, jjm_f
2272	REAL, INTENT(in) :: lon(iim_f, jjm_f), lat(iim_f, jjm_f)
2273	INTEGER, INTENT(out) :: iim,jjm
2274	INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f)
2275	!
2276	! LOCAL
2277	!
2278	INTEGER :: i, j
2279	REAL :: lolo
2280	LOGICAL :: over_dateline = .FALSE.
2281	!
2282	!
2283	IF ( ( ABS(limit_east) .GT. 180. ) .OR. &
2284	( ABS(limit_west) .GT. 180. ) ) THEN
2285	WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
2286	CALL ipslerr_p (3,'domain_size', &
2287	& 'Longitudes problem.','In run.def file :', &
2288	& 'limit_east > 180. or limit_west > 180.')
2289	ENDIF
2290	!
2291	IF ( limit_west .GT. limit_east ) over_dateline = .TRUE.
2292	!
2293	IF ( ( limit_south .LT. -90. ) .OR. &
2294	( limit_north .GT. 90. ) .OR. &
2295	( limit_south .GE. limit_north ) ) THEN
2296	WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
2297	CALL ipslerr_p (3,'domain_size', &
2298	& 'Latitudes problem.','In run.def file :', &
2299	& 'limit_south < -90. or limit_north > 90. or limit_south >= limit_north')
2300	ENDIF
2301	!
2302	! Here we assume that the grid of the forcing data is regular. Else we would have
2303	! to do more work to find the index table.
2304	!
2305	iim = 0
2306	DO i=1,iim_f
2307	!
2308	lolo = lon(i,1)
2309	IF ( lon(i,1) .GT. 180. ) lolo = lon(i,1) - 360.
2310	IF ( lon(i,1) .LT. -180. ) lolo = lon(i,1) + 360.
2311	!
2312	IF (lon(i,1) < limit_west) iim_g_begin = i+1
2313	IF (lon(i,1) < limit_east) iim_g_end = i
2314	!
2315	IF ( over_dateline ) THEN
2316	IF ( lolo .LE. limit_west .OR. lolo .GE. limit_east ) THEN
2317	iim = iim + 1
2318	iind(iim) = i
2319	ENDIF
2320	ELSE
2321	IF ( lolo .GE. limit_west .AND. lolo .LE. limit_east ) THEN
2322	iim = iim + 1
2323	iind(iim) = i
2324	ENDIF
2325	ENDIF
2326	!
2327	ENDDO
2328	!
2329	jjm = 0
2330	DO j=1,jjm_f
2331	IF (lat(1,j) > limit_north) jjm_g_begin = j+1
2332	IF (lat(1,j) > limit_south) jjm_g_end = j
2333	!
2334	IF ( lat(1,j) .GE. limit_south .AND. lat(1,j) .LE. limit_north) THEN
2335	jjm = jjm + 1
2336	jind(jjm) = j
2337	ENDIF
2338	ENDDO
2339	!
2340	WRITE(numout,*) 'Domain zoom size: iim, jjm = ', iim, jjm
2341	END SUBROUTINE domain_size
2342	!-
2343	!=====================================================================
2344	!-
2345	SUBROUTINE flinget_buffer(force_id, varname, iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
2346
2347	!! This subroutine is a wrap of flinget/IOIPSL. The arguments are the same.
2348	!! flinget_buffer will call flinget and buffer the forcing data localy in this subroutine.
2349	!! According to the variable NBUFF set in run.def, several time steps can be read at the same time
2350	!! from the forcing file. If NBUFF=0, the full forcing file is read.
2351	!! The output, data_full, from this subroutine is always only one time step of corresponding to itb.
2352	!! itb must be equal to ite.
2353
2354	!! Input arguments
2355	INTEGER, INTENT(in) :: force_id !! Id for forcing file
2356	CHARACTER(len=*), INTENT(in) :: varname !! Name of current variable to be read
2357	INTEGER, INTENT(in) :: iim_full, jjm_full, llm_full !! Horizontal and vertical domaine
2358	INTEGER, INTENT(in) :: ttm !! Full lenght of forcing file
2359	INTEGER, INTENT(in) :: itb, ite !! Time step to be read from forcing file. itb must be equal to ite
2360
2361	!! Output argument
2362	REAL, DIMENSION(iim_full, jjm_full), INTENT(out) :: data_full !! Data for time step itb.
2363
2364	!! Define specific type to buffer data together with name and index
2365	TYPE buffer_type
2366	CHARACTER(len=20) :: name !! Name of variable in forcing file
2367	INTEGER :: istart !! Start index of current buffered data
2368	INTEGER :: iend !! End index of current buffered data
2369	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: data !! Data read from forcing file for intervall [istart,iend]
2370	END TYPE buffer_type
2371
2372	!! Local variables
2373	INTEGER, PARAMETER :: maxvar=20 !! Max number of variables to be buffered
2374	TYPE(buffer_type), DIMENSION(maxvar),SAVE :: data_buffer !! Containing all variables and the current buffered data
2375	INTEGER, SAVE :: nbuff !! Number of time steps to be buffered
2376	INTEGER, SAVE :: lastindex=0 !! Current number of variables stored in data_buffer
2377	INTEGER, SAVE :: ttm0 !! Time lenght of forcing file, stored for test purpose
2378	LOGICAL, SAVE :: first=.TRUE. !! First call to this subroutine
2379	INTEGER :: index !! Index in data_buffer for current variable
2380	INTEGER :: i, ierr !! Loop and error variables
2381
2382
2383	!! 1. Initialization
2384	IF (first) THEN
2385	data_buffer(:)%name='undef'
2386	! Read NBUFF from run.def
2387	! Note that getin_p is not used because this subroutine might be called only by master process
2388
2389	!Config Key = NBUFF
2390	!Config Desc = Number of time steps of data to buffer between each reading of the forcing file
2391	!Config If = OFF_LINE
2392	!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.
2393	!Config Def = 1
2394	!Config Units= -
2395
2396	nbuff=1
2397	CALL getin('NBUFF', nbuff)
2398
2399	IF (nbuff == 0 .OR. nbuff >ttm) THEN
2400	! Set nbuff as the full forcing file lenght
2401	nbuff=ttm
2402	ELSE IF (nbuff < 0) THEN
2403	! Negativ nbuff not possible
2404	CALL ipslerr_p(3,'flinget_buffer','NBUFF must be a positiv number','Set NBUFF=0 for full simulation lenght','')
2405	END IF
2406	WRITE(numout,*)'flinget_buffer: NBUFF=',nbuff,' number of time step will be buffered'
2407	WRITE(numout,*)'flinget_buffer: Choose a lower value for NBUFF if problem with memory'
2408
2409	! Save dimensions to check following timesteps
2410	! ttm is the full lenght of forcing file
2411	ttm0=ttm
2412
2413	first=.FALSE.
2414	END IF
2415
2416	!! 2. Coeherence tests on input arguments
2417	IF (ttm /= ttm0) THEN
2418	WRITE(numout,*)'Problem with ttm=',ttm,' ttm0=',ttm0
2419	CALL ipslerr_p(3,'flinget_buffer','Error with ttm and ttm0','','')
2420	END IF
2421	IF (itb /= ite) THEN
2422	WRITE(numout,*) 'There is a problem. Why is itb not equal ite ?'
2423	WRITE(numout,*) 'itb=',itb,' ite=',ite,' varname=',varname
2424	CALL ipslerr_p(3,'flinget_buffer','ite not equal itb','','')
2425	END IF
2426
2427
2428	!! 3. Find index for current variable
2429	index=0
2430	DO i=1, maxvar
2431	IF ( trim(varname) == data_buffer(i)%name ) THEN
2432	index=i
2433	CYCLE
2434	END IF
2435	END DO
2436
2437	!! 4. Initialize and allocate if necesary the current variable
2438	IF ( index == 0 ) THEN
2439	! The variable was not found
2440	! This must be the first time for current variable
2441	index=lastindex+1
2442	lastindex=index
2443	IF (index > maxvar) CALL ipslerr_p(3,'flinget_buffer','to many variables','maxvar is too small','')
2444
2445	! Initialize the data_buffer for this index
2446	data_buffer(index)%name=trim(varname)
2447	ALLOCATE(data_buffer(index)%data(iim_full,jjm_full,nbuff),stat=ierr)
2448	IF (ierr /= 0) CALL ipslerr_p(3,'flinget_buffer','pb alloc data_buffer%data','for variable=',varname)
2449	data_buffer(index)%istart=0
2450	data_buffer(index)%iend=0
2451	END IF
2452
2453
2454	!! 5. Call flinget if current time step (itb) is outside the buffered intervall
2455	IF (( itb > data_buffer(index)%iend ) .OR. ( itb < data_buffer(index)%istart )) THEN
2456	! itb is not in the time slice previously read or it is the first time to read
2457	! Reading of forcing file will now be done
2458	! First recalculate index to be read
2459	data_buffer(index)%istart = itb
2460	data_buffer(index)%iend = itb + nbuff - 1
2461
2462	! WRITE(numout,*) 'Now do flinget for ',varname,', itb=',itb,', istart=',&
2463	! data_buffer(index)%istart,', iend=',data_buffer(index)%iend
2464	CALL flinget (force_id,varname, iim_full, jjm_full, llm_full, ttm, data_buffer(index)%istart, &
2465	data_buffer(index)%iend, data_buffer(index)%data(:,:,:))
2466	END IF
2467
2468	!! 6. Initialize the output variable with data from buffered variable
2469	! Find index for the time step corrsponding to itb in the time slice previously read from forcing file
2470	i=itb-data_buffer(index)%istart+1
2471	! Initialize output variable
2472	data_full(:,:) = data_buffer(index)%data(:,:,i)
2473
2474
2475	END SUBROUTINE flinget_buffer
2476	!------------------
2477	END MODULE readdim2

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

Download in other formats: