Context Navigation

source: tags/ORCHIDEE_OL/readdim2.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago
import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 61.3 KB

Line
1	!
2	MODULE readdim2
3	!---------------------------------------------------------------------
4	!- $Header: /home/ssipsl/CVSREP/ORCHIDEE_OL/readdim2.f90,v 1.23 2010/04/22 13:11:24 ssipsl Exp $
5	!- IPSL (2006)
6	!- This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
7	!-
8	USE ioipsl
9	USE weather
10	USE TIMER
11	USE constantes
12	USE constantes_veg
13	USE solar
14	USE grid
15	!-
16	IMPLICIT NONE
17	!-
18	PRIVATE
19	PUBLIC :: forcing_read, forcing_info, forcing_grid
20	!-
21	INTEGER, SAVE :: iim_full, jjm_full, llm_full, ttm_full
22	INTEGER, SAVE :: iim_zoom, jjm_zoom
23	INTEGER, SAVE :: iim_g_begin,jjm_g_begin,iim_g_end,jjm_g_end
24	REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: data_full, lon_full, lat_full
25	REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: lev_full
26	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: itau, i_index, j_index,j_index_g
27	INTEGER, SAVE :: i_test, j_test
28	LOGICAL, SAVE :: allow_weathergen, interpol
29	LOGICAL, SAVE, PUBLIC :: weathergen, is_watchout
30	REAL, SAVE :: merid_res, zonal_res
31	LOGICAL, SAVE :: have_zaxis=.FALSE.
32	!-
33	CONTAINS
34	!-
35	!=====================================================================
36	!-
37	SUBROUTINE forcing_info(filename, iim, jjm, llm, tm, date0, dt_force,&
38	& force_id)
39
40	!---------------------------------------------------------------------
41	!
42	!- This subroutine will get all the info from the forcing file and
43	!- prepare for the zoom if needed.
44	!- It returns to the caller the sizes of the data it will receive at
45	!- the forcing_read call. This is important so that the caller can
46	!- allocate the right space.
47	!-
48	!- filename : name of the file to be opened
49	!- iim : size in x of the forcing data
50	!- jjm : size in y of the forcing data
51	!- llm : number of levels in the forcing data (not yet used)
52	!- tm : Time dimension of the forcing
53	!- date0 : The date at which the forcing file starts (julian days)
54	!- dt_force : time-step of the forcing file in seconds
55	!- force_id : ID of the forcing file
56	!-
57	!- ARGUMENTS
58	!-
59	USE parallel
60	IMPLICIT NONE
61	!-
62	CHARACTER(LEN=*) :: filename
63	INTEGER :: iim, jjm, llm, tm
64	REAL :: date0, dt_force
65	INTEGER, INTENT(OUT) :: force_id
66	!- LOCAL
67	CHARACTER(LEN=20) :: calendar_str
68	REAL :: juld_1, juld_2
69	LOGICAL :: debug = .FALSE.
70	REAL, ALLOCATABLE, DIMENSION(:,:) :: fcontfrac
71	REAL, ALLOCATABLE, DIMENSION(:,:) :: tair
72	LOGICAL :: contfrac_exists=.FALSE.
73	INTEGER :: NbPoint
74	INTEGER :: i_test,j_test
75	INTEGER :: i,j,ind
76	INTEGER, ALLOCATABLE, DIMENSION(:) :: index_l
77	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
78	REAL, ALLOCATABLE, DIMENSION(:) :: lev, levuv
79
80	!-
81	CALL flininfo(filename, iim_full, jjm_full, llm_full, ttm_full, force_id)
82	CALL flinclo(force_id)
83	!-
84	IF ( debug ) WRITE(numout,*) 'forcing_info : Details from forcing file :', &
85	iim_full, jjm_full, llm_full, ttm_full
86	!-
87	IF ( llm_full < 1 ) THEN
88	have_zaxis = .FALSE.
89	ELSE
90	have_zaxis = .TRUE.
91	ENDIF
92	WRITE(numout,*) 'have_zaxis : ', llm_full, have_zaxis
93	!-
94	ALLOCATE(itau(ttm_full))
95	ALLOCATE(data_full(iim_full, jjm_full),lon_full(iim_full, jjm_full),&
96	& lat_full(iim_full, jjm_full))
97	ALLOCATE(lev_full(llm_full))
98	ALLOCATE(fcontfrac(iim_full,jjm_full))
99	ALLOCATE(i_index(iim_full), j_index(jjm_full),j_index_g(jjm_full))
100	!-
101	lev_full(:) = zero
102	!-
103	dt_force=zero
104	CALL flinopen &
105	& (filename, .FALSE., iim_full, jjm_full, llm_full, lon_full, lat_full, &
106	& lev_full, ttm_full, itau, date0, dt_force, force_id)
107	IF ( dt_force == zero ) THEN
108	dt_force = itau(2) - itau(1)
109	itau(:) = itau(:) / dt_force
110	ENDIF
111	! WRITE(numout,*) "forcing_info : Forcing time step out of flinopen : ",dt_force
112	!-
113	!- What are the alowed options for the temportal interpolation
114	!-
115	!Config Key = ALLOW_WEATHERGEN
116	!Config Desc = Allow weather generator to create data
117	!Config Def = n
118	!Config Help = This flag allows the forcing-reader to generate
119	!Config synthetic data if the data in the file is too sparse
120	!Config and the temporal resolution would not be enough to
121	!Config run the model.
122	!-
123	allow_weathergen = .FALSE.
124	CALL getin_p('ALLOW_WEATHERGEN',allow_weathergen)
125	!-
126	!- The calendar was set by the forcing file. If no "calendar" attribute was
127	!- found then it is assumed to be gregorian,
128	!MM => FALSE !! it is NOT assumed anything !
129	!- else it is what ever is written in this attribute.
130	!-
131	CALL ioget_calendar(calendar_str)
132	! WRITE(numout,*) "forcing_info : Calendar used : ",calendar_str
133	IF ( calendar_str == 'XXXX' ) THEN
134	WRITE(numout,*) "forcing_info : The calendar was not found in the forcing file."
135	IF (allow_weathergen) THEN
136	! Then change the calendar
137	CALL ioconf_calendar("noleap")
138	ELSE
139	WRITE(numout,*) "forcing_info : We will force it to gregorian by default."
140	CALL ioconf_calendar("gregorian") !! = 365.2425 ; "noleap" = 365.0; "360d"; "julian"=365.25
141	ENDIF
142	ENDIF
143	WRITE(numout,*) "readdim2 : Calendar used by the model : ",calendar_str
144	IF (ttm_full .GE. 2) THEN
145	juld_1 = itau2date(itau(1), date0, dt_force)
146	juld_2 = itau2date(itau(2), date0, dt_force)
147	ELSE
148	juld_1 = 0
149	juld_2 = 0
150	CALL ipslerr ( 3, 'forcing_info','What is that only one time step in the forcing file ?', &
151	& ' That can not be right.','verify forcing file.')
152	STOP
153	ENDIF
154	!-
155	!- Initialize one_year / one_day
156	CALL ioget_calendar (one_year, one_day)
157	!-
158	!- What is the distance between the two first states. From this we will deduce what is
159	!- to be done.
160	weathergen = .FALSE.
161	interpol = .FALSE.
162	is_watchout = .FALSE.
163	!-
164	IF ( ABS(ABS(juld_2-juld_1)-30.) .LE. 2.) THEN
165	IF ( allow_weathergen ) THEN
166	weathergen = .TRUE.
167	WRITE(numout,*) 'Using weather generator.'
168	ELSE
169	CALL ipslerr ( 3, 'forcing_info', &
170	& 'This seems to be a monthly file.', &
171	& 'We should use a weather generator with this file.', &
172	& 'This should be allowed in the run.def')
173	ENDIF
174	ELSEIF ( ABS(juld_1-juld_2) .LE. 1./4.) THEN
175	interpol = .TRUE.
176	WRITE(numout,*) 'We will interpolate between the forcing data time steps.'
177	ELSE
178	! Using the weather generator with data other than monthly ones probably
179	! needs some thinking.
180	WRITE(numout,*) 'The time step is not suitable:',ABS(juld_1-juld_2),' days.'
181	CALL ipslerr ( 3, 'forcing_info','The time step is not suitable.', &
182	& '','We cannot do anything with these forcing data.')
183	ENDIF
184	!-
185	!- redefine the forcing time step if the weather generator is activated
186	!-
187	IF ( weathergen ) THEN
188	!Config Key = DT_WEATHGEN
189	!Config Desc = Calling frequency of weather generator (s)
190	!Config If = ALLOW_WEATHERGEN
191	!Config Def = 1800.
192	!Config Help = Determines how often the weather generator
193	!Config is called (time step in s). Should be equal
194	!Config to or larger than Sechiba's time step (say,
195	!Config up to 6 times Sechiba's time step or so).
196	dt_force = 1800.
197	CALL getin_p('DT_WEATHGEN',dt_force)
198	ENDIF
199	!-
200	!- Define the zoom
201	!-
202	!Config Key = LIMIT_WEST
203	!Config Desc = Western limit of region
204	!Config Def = -180.
205	!Config Help = Western limit of the region we are
206	!Config interested in. Between -180 and +180 degrees
207	!Config The model will use the smalest regions from
208	!Config region specified here and the one of the forcing file.
209	!-
210	limit_west = -180.
211	CALL getin_p('LIMIT_WEST',limit_west)
212	!-
213	!Config Key = LIMIT_EAST
214	!Config Desc = Eastern limit of region
215	!Config Def = 180.
216	!Config Help = Eastern limit of the region we are
217	!Config interested in. Between -180 and +180 degrees
218	!Config The model will use the smalest regions from
219	!Config region specified here and the one of the forcing file.
220	!-
221	limit_east = 180.
222	CALL getin_p('LIMIT_EAST',limit_east)
223	!-
224	!Config Key = LIMIT_NORTH
225	!Config Desc = Northern limit of region
226	!Config Def = 90.
227	!Config Help = Northern limit of the region we are
228	!Config interested in. Between +90 and -90 degrees
229	!Config The model will use the smalest regions from
230	!Config region specified here and the one of the forcing file.
231	!-
232	limit_north = 90.
233	CALL getin_p('LIMIT_NORTH',limit_north)
234	!-
235	!Config Key = LIMIT_SOUTH
236	!Config Desc = Southern limit of region
237	!Config Def = -90.
238	!Config Help = Southern limit of the region we are
239	!Config interested in. Between 90 and -90 degrees
240	!Config The model will use the smalest regions from
241	!Config region specified here and the one of the forcing file.
242	!-
243	limit_south = -90.
244	CALL getin_p('LIMIT_SOUTH',limit_south)
245	!-
246	!- Calculate domain size
247	!-
248	IF ( interpol ) THEN
249	!-
250	!- If we use temporal interpolation, then we cannot change the resolution (yet?)
251	!-
252	IF (is_root_prc) THEN
253
254	CALL domain_size (limit_west, limit_east, limit_north, limit_south,&
255	& iim_full, jjm_full, lon_full, lat_full, iim_zoom, jjm_zoom,&
256	& i_index, j_index_g)
257
258	j_index(:)=j_index_g(:)
259
260	ALLOCATE(tair(iim_full,jjm_full))
261	CALL flinget (force_id,'Tair',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
262	CALL forcing_zoom(data_full, tair)
263
264	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
265	IF ( contfrac_exists ) THEN
266	WRITE(numout,*) "contfrac exist in the forcing file."
267	CALL flinget (force_id,'contfrac',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
268	CALL forcing_zoom(data_full, fcontfrac)
269	WRITE(numout,*) "fcontfrac min/max :",MINVAL(fcontfrac(1:iim_zoom,1:jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,1:jjm_zoom))
270	ELSE
271	fcontfrac(:,:)=1.
272	ENDIF
273
274
275	DO i=1,iim_zoom
276	DO j=1,jjm_zoom
277	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
278	tair(i,j) = 999999.
279	ENDIF
280	ENDDO
281	ENDDO
282
283	ALLOCATE(index_l(iim_zoom*jjm_zoom))
284	!- In this point is returning the global NbPoint with the global index
285	CALL forcing_landind(iim_zoom,jjm_zoom,tair,.TRUE.,NbPoint,index_l,i_test,j_test)
286	ELSE
287	ALLOCATE(index_l(1))
288	ENDIF
289
290	CALL init_data_para(iim_zoom,jjm_zoom,NbPoint,index_l)
291
292	!
293	!- global index index_g is the index_l of root proc
294	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
295
296	DEALLOCATE(index_l)
297
298	CALL bcast(jjm_zoom)
299	CALL bcast(i_index)
300	CALL bcast(j_index_g)
301
302	ind=0
303	DO j=1,jjm_zoom
304	IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
305	ind=ind+1
306	j_index(ind)=j_index_g(j)
307	ENDIF
308	ENDDO
309
310	jjm_zoom=jj_nb
311	iim_zoom=iim_g
312
313	!-
314	!- If we use the weather generator, then we read zonal and meridional resolutions
315	!- This should be unified one day...
316	!-
317	ELSEIF ( weathergen ) THEN
318	!-
319	!Config Key = MERID_RES
320	!Config Desc = North-South Resolution
321	!Config Def = 2.
322	!Config If = ALLOW_WEATHERGEN
323	!Config Help = North-South Resolution of the region we are
324	!Config interested in. In degrees
325	!-
326	merid_res = 2.
327	CALL getin_p('MERID_RES',merid_res)
328	!-
329	!Config Key = ZONAL_RES
330	!Config Desc = East-West Resolution
331	!Config Def = 2.
332	!Config If = ALLOW_WEATHERGEN
333	!Config Help = East-West Resolution of the region we are
334	!Config interested in. In degrees
335	!-
336	zonal_res = 2.
337	CALL getin_p('ZONAL_RES',zonal_res)
338	!-
339	!- Number of time steps is meaningless in this case
340	!-
341	! ttm_full = HUGE( ttm_full )
342	!MM Number (realistic) of time steps for half hour dt
343	ttm_full = NINT(one_year * 86400. / dt_force)
344	!-
345	IF (is_root_prc) THEN
346
347	!- In this point is returning the global NbPoint with the global index
348	CALL weathgen_domain_size (limit_west,limit_east,limit_north,limit_south, &
349	zonal_res,merid_res,iim_zoom,jjm_zoom)
350	ALLOCATE(index_l(iim_zoom*jjm_zoom))
351	ENDIF
352	CALL bcast(iim_zoom)
353	CALL bcast(jjm_zoom)
354
355	ALLOCATE(lon(iim_zoom,jjm_zoom))
356	ALLOCATE(lat(iim_zoom,jjm_zoom))
357	ALLOCATE(lev(llm_full),levuv(llm_full))
358
359	! We need lon and lat now for weathgen_init
360	CALL forcing_grid (iim_zoom,jjm_zoom,llm_full,lon,lat,lev,levuv,init_f=.TRUE.)
361
362	IF (is_root_prc) THEN
363	CALL weathgen_init &
364	& (filename,dt_force,force_id,iim_zoom,jjm_zoom, &
365	& zonal_res,merid_res,lon,lat,index_l,NbPoint,&
366	& i_index,j_index_g)
367	ELSE
368	ALLOCATE(index_l(1))
369	index_l(1)=1
370	ENDIF
371
372	CALL init_data_para(iim_zoom,jjm_zoom,NbPoint,index_l)
373
374	!
375	!- global index index_g is the index_l of root proc
376	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
377
378	DEALLOCATE(index_l)
379
380	CALL bcast(i_index)
381	CALL bcast(j_index_g)
382
383	ind=0
384	DO j=1,jjm_zoom
385	IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
386	ind=ind+1
387	j_index(ind)=j_index_g(j)
388	ENDIF
389	ENDDO
390
391	jjm_zoom=jj_nb
392	iim_zoom=iim_g
393	!
394	CALL weathgen_read_file(force_id,iim_zoom,jjm_zoom)
395
396	!-
397	ELSE
398	!-
399	STOP 'ERROR: Neither interpolation nor weather generator is specified.'
400	!-
401	ENDIF
402	!-
403	!- Transfer the right information to the caller
404	!-
405	iim = iim_zoom
406	jjm = jjm_zoom
407	llm = 1
408	tm = ttm_full
409	!-
410	IF ( debug ) WRITE(numout,*) 'forcing_info : end : ', iim,jjm, llm,tm
411	!-
412	END SUBROUTINE forcing_info
413	!-
414	!-
415	!=====================================================================
416	SUBROUTINE forcing_read &
417	& (filename, rest_id, lrstread, lrstwrite, &
418	& itauin, istp, itau_split, split, nb_spread, netrad_cons, date0, &
419	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
420	& swdown, precip, snowf, tair, u, v, qair, pb, lwdown, &
421	& fcontfrac, fneighbours, fresolution, &
422	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
423	& kindex, nbindex, force_id)
424	!---------------------------------------------------------------------
425	!- filename : name of the file to be opened
426	!- rest_id : ID of restart file
427	!- lrstread : read restart file?
428	!- lrstwrite : write restart file?
429	!- itauin : time step for which we need the data
430	!- istp : time step for restart file
431	!- itau_split : Where are we within the splitting
432	!- of the time-steps of the forcing files
433	!- (it decides IF we READ or not)
434	!- split : The number of time steps we do
435	!- between two time-steps of the forcing
436	!- nb_spread : Over how many time steps do we spread the precipitation
437	!- netrad_cons: flag that decides if net radiation should be conserved.
438	!- date0 : The date at which the forcing file starts (julian days)
439	!- dt_force : time-step of the forcing file in seconds
440	!- iim : Size of the grid in x
441	!- jjm : size of the grid in y
442	!- lon : Longitudes
443	!- lat : Latitudes
444	!- zlev : First Levels if it exists (ie if watchout file)
445	!- zlevuv : First Levels of the wind (equal precedent, if it exists)
446	!- ttm : number of time steps in all in the forcing file
447	!- swdown : Downward solar radiation (W/m^2)
448	!- precip : Precipitation (Rainfall) (kg/m^2s)
449	!- snowf : Snowfall (kg/m^2s)
450	!- tair : 1st level (2m ? in off-line) air temperature (K)
451	!- u and v : 1st level (2m/10m ? in off-line) (in theory !) wind speed (m/s)
452	!- qair : 1st level (2m ? in off-line) humidity (kg/kg)
453	!- pb : Surface pressure (Pa)
454	!- lwdown : Downward long wave radiation (W/m^2)
455	!- fcontfrac : Continental fraction (no unit)
456	!- fneighbours: land neighbours
457	!- fresolution: resolution in x and y dimensions for each points
458	!-
459	!- From a WATCHOUT file :
460	!- SWnet : Net surface short-wave flux
461	!- Eair : Air potential energy
462	!- petAcoef : Coeficients A from the PBL resolution for T
463	!- peqAcoef : Coeficients A from the PBL resolution for q
464	!- petBcoef : Coeficients B from the PBL resolution for T
465	!- peqBcoef : Coeficients B from the PBL resolution for q
466	!- cdrag : Surface drag
467	!- ccanopy : CO2 concentration in the canopy
468	!-
469	!- kindex : Index of all land-points in the data
470	!- (used for the gathering)
471	!- nbindex : Number of land points
472	!- force_id : FLINCOM file id.
473	!- It is used to close the file at the end of the run.
474	!-
475	!---------------------------------------------------------------------
476	IMPLICIT NONE
477	!-
478	CHARACTER(LEN=*) :: filename
479	INTEGER, INTENT(IN) :: force_id
480	INTEGER, INTENT(INOUT) :: nbindex
481	INTEGER :: rest_id
482	LOGICAL :: lrstread, lrstwrite
483	INTEGER :: itauin, istp, itau_split, split, nb_spread
484	LOGICAL :: netrad_cons
485	REAL :: date0, dt_force
486	INTEGER :: iim, jjm, ttm
487	REAL,DIMENSION(iim,jjm) :: lon, lat, zlev, zlevuv, &
488	& swdown, precip, snowf, tair, u, v, qair, pb, lwdown, &
489	& fcontfrac
490	REAL,DIMENSION(iim,jjm,2) :: fresolution
491	INTEGER,DIMENSION(iim,jjm,8) :: fneighbours
492	! for watchout files
493	REAL,DIMENSION(iim,jjm) :: &
494	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
495	INTEGER,DIMENSION(iim*jjm), INTENT(OUT) :: kindex
496	!-
497	INTEGER :: ik,i,j
498	!
499	IF ( interpol ) THEN
500	!-
501	CALL forcing_read_interpol &
502	(filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, &
503	dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
504	swdown, 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	ELSEIF ( weathergen ) THEN
510	!-
511	IF (lrstread) THEN
512	fcontfrac(:,:) = 1.0
513	WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
514	ENDIF
515
516	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
517	CALL weathgen_main (istp, istp, filename, force_id, iim, jjm, 1, &
518	rest_id, lrstread, lrstwrite, &
519	limit_west, limit_east, limit_north, limit_south, &
520	zonal_res, merid_res, lon, lat, date0, dt_force, &
521	kindex, nbindex, &
522	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
523	ELSE
524	CALL weathgen_main (itauin, istp, filename, force_id, iim, jjm, 1, &
525	rest_id, lrstread, lrstwrite, &
526	limit_west, limit_east, limit_north, limit_south, &
527	zonal_res, merid_res, lon, lat, date0, dt_force, &
528	kindex, nbindex, &
529	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
530	ENDIF
531	!-
532	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
533	!---
534	!--- Allocate grid stuff
535	!---
536	CALL init_grid ( nbindex )
537	!---
538	!--- Compute
539	!---
540	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, kindex)
541	!CALL grid_stuff (nbindex, iim, jjm, lon, lat, kindex)
542	DO ik=1,nbindex
543
544	j = ((kindex(ik)-1)/iim) + 1
545	i = (kindex(ik) - (j-1)*iim)
546	!-
547	!- Store variable to help describe the grid
548	!- once the points are gathered.
549	!-
550	fneighbours(i,j,:) = neighbours(ik,:)
551	!
552	fresolution(i,j,:) = resolution(ik,:)
553	ENDDO
554	ENDIF
555	ELSE
556	!-
557	STOP 'ERROR: Neither interpolation nor weather generator is specified.'
558	!-
559	ENDIF
560	!-
561	IF (.NOT. is_watchout) THEN
562	! We have to compute Potential air energy
563	WHERE(tair(:,:) < val_exp)
564	eair(:,:) = cp_airtair(:,:)+cte_gravzlev(:,:)
565	ENDWHERE
566	ENDIF
567	!-
568	!
569	!-------------------------
570	END SUBROUTINE forcing_read
571	!=====================================================================
572	!-
573	!-
574	!=====================================================================
575	SUBROUTINE forcing_read_interpol &
576	& (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0, &
577	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, swdown, rainf, snowf, tair, &
578	& u, v, qair, pb, lwdown, &
579	& fcontfrac, fneighbours, fresolution, &
580	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
581	& kindex, nbindex, force_id)
582	!---------------------------------------------------------------------
583	!- filename : name of the file to be opened
584	!- itauin : time step for which we need the data
585	!- itau_split : Where are we within the splitting
586	!- of the time-steps of the forcing files
587	!- (it decides IF we READ or not)
588	!- split : The number of time steps we do
589	!- between two time-steps of the forcing
590	!- nb_spread : Over how many time steps do we spread the precipitation
591	!- netrad_cons: flag that decides if net radiation should be conserved.
592	!- date0 : The date at which the forcing file starts (julian days)
593	!- dt_force : time-step of the forcing file in seconds
594	!- iim : Size of the grid in x
595	!- jjm : size of the grid in y
596	!- lon : Longitudes
597	!- lat : Latitudes
598	!- zlev : First Levels if it exists (ie if watchout file)
599	!- zlevuv : First Levels of the wind (equal precedent, if it exists)
600	!- ttm : number of time steps in all in the forcing file
601	!- swdown : Downward solar radiation (W/m^2)
602	!- rainf : Rainfall (kg/m^2s)
603	!- snowf : Snowfall (kg/m^2s)
604	!- tair : 2m air temperature (K)
605	!- u and v : 2m (in theory !) wind speed (m/s)
606	!- qair : 2m humidity (kg/kg)
607	!- pb : Surface pressure (Pa)
608	!- lwdown : Downward long wave radiation (W/m^2)
609	!- fcontfrac : Continental fraction (no unit)
610	!- fneighbours: land neighbours
611	!- fresolution: resolution in x and y dimensions for each points
612	!-
613	!- From a WATCHOUT file :
614	!- SWnet : Net surface short-wave flux
615	!- Eair : Air potential energy
616	!- petAcoef : Coeficients A from the PBL resolution for T
617	!- peqAcoef : Coeficients A from the PBL resolution for q
618	!- petBcoef : Coeficients B from the PBL resolution for T
619	!- peqBcoef : Coeficients B from the PBL resolution for q
620	!- cdrag : Surface drag
621	!- ccanopy : CO2 concentration in the canopy
622	!-
623	!- kindex : Index of all land-points in the data
624	!- (used for the gathering)
625	!- nbindex : Number of land points
626	!- force_id : FLINCOM file id.
627	!- It is used to close the file at the end of the run.
628	!---------------------------------------------------------------------
629	USE parallel
630	IMPLICIT NONE
631	!-
632	INTEGER,PARAMETER :: lm=1
633	!-
634	!- Input variables
635	!-
636	CHARACTER(LEN=*) :: filename
637	INTEGER :: itauin, itau_split, split, nb_spread
638	LOGICAL :: netrad_cons
639	REAL :: date0, dt_force
640	INTEGER :: iim, jjm, ttm
641	REAL,DIMENSION(:,:),INTENT(IN) :: lon, lat !- LOCAL data array (size=iim,jjm)
642	INTEGER, INTENT(IN) :: force_id
643	!-
644	!- Output variables
645	!-
646	REAL,DIMENSION(:,:),INTENT(OUT) :: zlev, zlevuv, & !- LOCAL data array (size=iim,jjm)
647	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown, &
648	& fcontfrac
649	REAL,DIMENSION(:,:,:),INTENT(OUT) :: fresolution !- LOCAL data array (size=iim,jjm,2)
650	INTEGER,DIMENSION(:,:,:),INTENT(OUT) :: fneighbours !- LOCAL data array (size=iim,jjm,8)
651	! for watchout files
652	REAL,DIMENSION(:,:),INTENT(OUT) :: &
653	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
654	INTEGER,DIMENSION(:),INTENT(INOUT) :: kindex !- LOCAL index of the map
655	INTEGER, INTENT(INOUT) :: nbindex
656	!-
657	!- Local variables
658	!-
659	INTEGER, SAVE :: last_read=0
660	INTEGER, SAVE :: itau_read, itau_read_nm1=0, itau_read_n=0
661	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
662	& zlev_nm1, swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
663	& zlev_n, swdown_n, rainf_n, snowf_n, tair_n, u_n, v_n, qair_n, &
664	& pb_n, lwdown_n, mean_sinang, sinang
665	! just for grid stuff if the forcing file is a watchout file
666	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmpdata
667	! variables to be read in watchout files
668	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
669	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
670	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n
671	REAL, SAVE :: julian_for ! Date of the forcing to be read
672	REAL :: julian, ss, rw
673	!jur,
674	REAL, SAVE :: julian0 ! First day of this year
675	INTEGER :: yy, mm, dd, is, i, j, ik
676	! if Wind_N and Wind_E are in the file (and not just Wind)
677	LOGICAL, SAVE :: wind_N_exists=.FALSE.
678	LOGICAL :: wind_E_exists=.FALSE.
679	LOGICAL, SAVE :: contfrac_exists=.FALSE.
680	LOGICAL, SAVE :: neighbours_exists=.FALSE.
681	LOGICAL, SAVE :: initialized = .FALSE.
682	LOGICAL :: check=.FALSE.
683	! to bypass FPE problem on integer convertion with missing_value heigher than precision
684	INTEGER, PARAMETER :: undef_int = 999999999
685	!---------------------------------------------------------------------
686	IF (check) THEN
687	WRITE(numout,*) &
688	& " FORCING READ : itau_read, itau_split : ",itauin,itau_split
689	ENDIF
690	!-
691	itau_read = itauin
692	!-
693	!- This part initializes the reading of the forcing. As you can see
694	!- we only go through here if both time steps are zero.
695	!-
696	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
697	!-
698	!- Tests on forcing file type
699	CALL flinquery_var(force_id, 'Wind_N', wind_N_exists)
700	IF ( wind_N_exists ) THEN
701	CALL flinquery_var(force_id, 'Wind_E', wind_E_exists)
702	IF ( .NOT. wind_E_exists ) THEN
703	CALL ipslerr(3,'forcing_read_interpol', &
704	& 'Variable Wind_E does not exist in forcing file', &
705	& 'But variable Wind_N exists.','Please, rename Wind_N to Wind;')
706	ENDIF
707	ENDIF
708	CALL flinquery_var(force_id, 'levels', is_watchout)
709	IF ( is_watchout ) THEN
710	WRITE(numout,*) "Read a Watchout File."
711	ENDIF
712	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
713	!-
714	IF (check) WRITE(numout,*) 'ALLOCATE all the memory needed'
715	!-
716	ALLOCATE &
717	& (swdown_nm1(iim,jjm), rainf_nm1(iim,jjm), snowf_nm1(iim,jjm), &
718	& tair_nm1(iim,jjm), u_nm1(iim,jjm), v_nm1(iim,jjm), qair_nm1(iim,jjm), &
719	& pb_nm1(iim,jjm), lwdown_nm1(iim,jjm))
720	ALLOCATE &
721	& (swdown_n(iim,jjm), rainf_n(iim,jjm), snowf_n(iim,jjm), &
722	& tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), &
723	& pb_n(iim,jjm), lwdown_n(iim,jjm))
724	! to read of watchout files
725	IF (is_watchout) THEN
726	ALLOCATE &
727	& (zlev_nm1(iim,jjm), zlev_n(iim,jjm), &
728	& SWnet_nm1(iim,jjm), Eair_nm1(iim,jjm), cdrag_nm1(iim,jjm), ccanopy_nm1(iim,jjm), &
729	& petAcoef_nm1(iim,jjm), peqAcoef_nm1(iim,jjm), petBcoef_nm1(iim,jjm), peqBcoef_nm1(iim,jjm), &
730	& SWnet_n(iim,jjm), Eair_n(iim,jjm), cdrag_n(iim,jjm), ccanopy_n(iim,jjm), &
731	& petAcoef_n(iim,jjm), peqAcoef_n(iim,jjm), petBcoef_n(iim,jjm), peqBcoef_n(iim,jjm))
732	ENDIF
733	ALLOCATE &
734	& (mean_sinang(iim,jjm), sinang(iim,jjm))
735	!-
736	IF (check) WRITE(numout,*) 'Memory ALLOCATED'
737	!-
738	!- We need for the driver surface air temperature and humidity before the
739	!- the loop starts. Thus we read it here.
740	!-
741	CALL forcing_just_read (iim, jjm, zlev, ttm, 1, 1, &
742	& swdown, rainf, snowf, tair, &
743	& u, v, qair, pb, lwdown, &
744	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
745	& force_id, wind_N_exists, check)
746	!----
747	IF (is_watchout) THEN
748	zlevuv(:,:) = zlev(:,:)
749	ENDIF
750	!-- First in case it's not a watchout file
751	IF ( .NOT. is_watchout ) THEN
752	IF ( contfrac_exists ) THEN
753	WRITE(numout,*) "contfrac exist in the forcing file."
754	CALL flinget (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
755	CALL forcing_zoom(data_full, fcontfrac)
756	WRITE(numout,*) "fcontfrac min/max :",MINVAL(fcontfrac(1:iim_zoom,jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,jjm_zoom))
757	!
758	! We need to make sure that when we gather the points we pick all
759	! the points where contfrac is above 0. Thus we prepare tair for
760	! subroutine forcing_landind
761	!
762	DO i=1,iim
763	DO j=1,jjm
764	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0. ! bande de recouvrement du scatter2D
765	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0. ! => on mets les données qu'on veut pas du noeud à missing_value
766	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
767	tair(i,j) = 999999.
768	ENDIF
769	ENDDO
770	ENDDO
771	ELSE
772	fcontfrac(:,:) = 1.0
773	ENDIF
774	!---
775	!--- Create the index table
776	!---
777	!--- This job return a LOCAL kindex
778	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
779	#ifdef CPP_PARA
780	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
781	! Force nbindex points for parallel computation
782	nbindex=nbp_loc
783	CALL scatter(index_g,kindex)
784	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
785	#endif
786	ik=MAX(nbindex/2,1)
787	j_test = (((kindex(ik)-1)/iim) + 1)
788	i_test = (kindex(ik) - (j_test-1)*iim)
789	IF (check) THEN
790	WRITE(numout,*) 'New test point is : ', i_test, j_test
791	ENDIF
792	!---
793	!--- Allocate grid stuff
794	!---
795	CALL init_grid ( nbindex )
796	!---
797	!--- All grid variables
798	!---
799	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, kindex)
800	DO ik=1,nbindex
801	!
802	j = ((kindex(ik)-1)/iim) + 1
803	i = (kindex(ik) - (j-1)*iim)
804	!-
805	!- Store variable to help describe the grid
806	!- once the points are gathered.
807	!-
808	fneighbours(i,j,:) = neighbours(ik,:)
809	!
810	fresolution(i,j,:) = resolution(ik,:)
811	ENDDO
812	ELSE
813	!-- Second, in case it is a watchout file
814	ALLOCATE (tmpdata(iim,jjm))
815	tmpdata(:,:) = 0.0
816	!--
817	IF ( .NOT. contfrac_exists ) THEN
818	CALL ipslerr (3,'forcing_read_interpol', &
819	& 'Could get contfrac variable in a watchout file :',TRIM(filename), &
820	& '(Problem with file ?)')
821	ENDIF
822	CALL flinget (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
823	CALL forcing_zoom(data_full, fcontfrac)
824	!
825	! We need to make sure that when we gather the points we pick all
826	! the points where contfrac is above 0. Thus we prepare tair for
827	! subroutine forcing_landind
828	!
829	DO i=1,iim
830	DO j=1,jjm
831	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0.
832	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0.
833	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
834	tair(i,j) = 999999.
835	ENDIF
836	ENDDO
837	ENDDO
838	!---
839	!--- Create the index table
840	!---
841	!--- This job return a LOCAL kindex
842	CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
843	#ifdef CPP_PARA
844	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
845	! Force nbindex points for parallel computation
846	nbindex=nbp_loc
847	CALL scatter(index_g,kindex)
848	kindex(:)=kindex(:)-(jj_begin-1)*iim_g
849	#endif
850	ik=MAX(nbindex/2,1)
851	j_test = (((kindex(ik)-1)/iim) + 1)
852	i_test = (kindex(ik) - (j_test-1)*iim)
853	IF (check) THEN
854	WRITE(numout,*) 'New test point is : ', i_test, j_test
855	ENDIF
856	!---
857	!--- Allocate grid stuff
858	!---
859	CALL init_grid ( nbindex )
860	neighbours(:,:) = -1
861	resolution(:,:) = 0.
862	min_resol(:) = 1.e6
863	max_resol(:) = -1.
864	!---
865	!--- All grid variables
866	!---
867	!-
868	!- Get variables to help describe the grid
869	CALL flinquery_var(force_id, 'neighboursNN', neighbours_exists)
870	IF ( .NOT. neighbours_exists ) THEN
871	CALL ipslerr (3,'forcing_read_interpol', &
872	& 'Could get neighbours in a watchout file :',TRIM(filename), &
873	& '(Problem with file ?)')
874	ELSE
875	WRITE(numout,*) "Watchout file contains neighbours and resolutions."
876	ENDIF
877	!
878	fneighbours(:,:,:) = undef_int
879	!
880	!- once the points are gathered.
881	CALL flinget (force_id,'neighboursNN',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
882	CALL forcing_zoom(data_full, tmpdata)
883	WHERE(tmpdata(:,:) < undef_int)
884	fneighbours(:,:,1) = NINT(tmpdata(:,:))
885	ENDWHERE
886	!
887	CALL flinget (force_id,'neighboursNE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
888	CALL forcing_zoom(data_full, tmpdata)
889	WHERE(tmpdata(:,:) < undef_int)
890	fneighbours(:,:,2) = NINT(tmpdata(:,:))
891	ENDWHERE
892	!
893	CALL flinget (force_id,'neighboursEE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
894	CALL forcing_zoom(data_full, tmpdata)
895	WHERE(tmpdata(:,:) < undef_int)
896	fneighbours(:,:,3) = NINT(tmpdata(:,:))
897	ENDWHERE
898	!
899	CALL flinget (force_id,'neighboursSE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
900	CALL forcing_zoom(data_full, tmpdata)
901	WHERE(tmpdata(:,:) < undef_int)
902	fneighbours(:,:,4) = NINT(tmpdata(:,:))
903	ENDWHERE
904	!
905	CALL flinget (force_id,'neighboursSS',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
906	CALL forcing_zoom(data_full, tmpdata)
907	WHERE(tmpdata(:,:) < undef_int)
908	fneighbours(:,:,5) = NINT(tmpdata(:,:))
909	ENDWHERE
910	!
911	CALL flinget (force_id,'neighboursSW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
912	CALL forcing_zoom(data_full, tmpdata)
913	WHERE(tmpdata(:,:) < undef_int)
914	fneighbours(:,:,6) = NINT(tmpdata(:,:))
915	ENDWHERE
916	!
917	CALL flinget (force_id,'neighboursWW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
918	CALL forcing_zoom(data_full, tmpdata)
919	WHERE(tmpdata(:,:) < undef_int)
920	fneighbours(:,:,7) = NINT(tmpdata(:,:))
921	ENDWHERE
922	!
923	CALL flinget (force_id,'neighboursNW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
924	CALL forcing_zoom(data_full, tmpdata)
925	WHERE(tmpdata(:,:) < undef_int)
926	fneighbours(:,:,8) = NINT(tmpdata(:,:))
927	ENDWHERE
928	!
929	! now, resolution of the grid
930	CALL flinget (force_id,'resolutionX',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
931	CALL forcing_zoom(data_full, tmpdata)
932	fresolution(:,:,1) = tmpdata(:,:)
933	!
934	CALL flinget (force_id,'resolutionY',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
935	CALL forcing_zoom(data_full, tmpdata)
936	fresolution(:,:,2) = tmpdata(:,:)
937	!
938	DO ik=1,nbindex
939	!
940	j = ((kindex(ik)-1)/iim) + 1
941	i = (kindex(ik) - (j-1)*iim)
942	!-
943	!- Store variable to help describe the grid
944	!- once the points are gathered.
945	!-
946	neighbours(ik,:) = fneighbours(i,j,:)
947	!
948	resolution(ik,:) = fresolution(i,j,:)
949	!
950
951	ENDDO
952	CALL gather(neighbours,neighbours_g)
953	CALL gather(resolution,resolution_g)
954	min_resol(1) = MINVAL(resolution(:,1))
955	min_resol(2) = MAXVAL(resolution(:,2))
956	max_resol(1) = MAXVAL(resolution(:,1))
957	max_resol(2) = MAXVAL(resolution(:,2))
958	!
959	area(:) = resolution(:,1)*resolution(:,2)
960	CALL gather(area,area_g)
961	!--
962	DEALLOCATE (tmpdata)
963	ENDIF
964	WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
965	!---
966	ENDIF
967	!---
968	IF (check) THEN
969	WRITE(numout,*) &
970	& 'The dates : ',itau_read,itau_split,itau_read_nm1,itau_read_n
971	ENDIF
972	!---
973	!--- Here we do the work in case only interpolation is needed.
974	!---
975	IF ( initialized .AND. interpol ) THEN
976	!---
977	IF (itau_read /= last_read) THEN
978	!---
979	!----- Start or Restart
980	IF (itau_read_n == 0) THEN
981	! Case of a restart or a shift in the forcing file.
982	IF (itau_read > 1) THEN
983	itau_read_nm1=itau_read-1
984	CALL forcing_just_read (iim, jjm, zlev_nm1, ttm, itau_read_nm1, itau_read_nm1, &
985	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
986	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
987	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
988	& force_id, wind_N_exists, check)
989	! Case of a simple start.
990	ELSE IF (dt_force*ttm > one_day-1. ) THEN
991	! if the forcing file contains at least 24 hours,
992	! we will use the last forcing step of the first day
993	! as initiale condition to prevent first shift off reading.
994	itau_read_nm1 = NINT (one_day/dt_force)
995	WRITE(numout,) "the forcing file contains 24 hours :",dt_forcettm,one_day-1.
996	WRITE(numout,*) "we will use the last forcing step of the first day : itau_read_nm1 ",itau_read_nm1
997	CALL forcing_just_read (iim, jjm, zlev_nm1, ttm, itau_read_nm1, itau_read_nm1, &
998	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
999	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1000	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1001	& force_id, wind_N_exists, check)
1002	ELSE
1003	! if the forcing file contains less than 24 hours,
1004	! just say error !
1005	CALL ipslerr(3,'forcing_read_interpol', &
1006	& 'The forcing file contains less than 24 hours !', &
1007	& 'We can''t intialize interpolation with such a file.','')
1008	ENDIF
1009	ELSE
1010	!----- Normal mode : copy old step
1011	swdown_nm1(:,:) = swdown_n(:,:)
1012	rainf_nm1(:,:) = rainf_n(:,:)
1013	snowf_nm1(:,:) = snowf_n(:,:)
1014	tair_nm1(:,:) = tair_n(:,:)
1015	u_nm1(:,:) = u_n(:,:)
1016	v_nm1(:,:) = v_n(:,:)
1017	qair_nm1(:,:) = qair_n(:,:)
1018	pb_nm1(:,:) = pb_n(:,:)
1019	lwdown_nm1(:,:) = lwdown_n(:,:)
1020	IF (is_watchout) THEN
1021	zlev_nm1(:,:) = zlev_n(:,:)
1022	! Net surface short-wave flux
1023	SWnet_nm1(:,:) = SWnet_n(:,:)
1024	! Air potential energy
1025	Eair_nm1(:,:) = Eair_n(:,:)
1026	! Coeficients A from the PBL resolution for T
1027	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1028	! Coeficients A from the PBL resolution for q
1029	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1030	! Coeficients B from the PBL resolution for T
1031	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1032	! Coeficients B from the PBL resolution for q
1033	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1034	! Surface drag
1035	cdrag_nm1(:,:) = cdrag_n(:,:)
1036	! CO2 concentration in the canopy
1037	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1038	ENDIF
1039	itau_read_nm1 = itau_read_n
1040	ENDIF
1041	!-----
1042	itau_read_n = itau_read
1043	IF (itau_read_n > ttm) THEN
1044	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1045	WRITE(numout,*) &
1046	& 'WARNING : We are going back to the start of the file'
1047	itau_read_n =1
1048	ENDIF
1049	IF (check) THEN
1050	WRITE(numout,*) &
1051	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1052	ENDIF
1053	!-----
1054	!----- Get a reduced julian day !
1055	!----- This is needed because we lack the precision on 32 bit machines.
1056	!-----
1057	IF ( dt_force .GT. 3600. ) THEN
1058	julian_for = itau2date(itau_read-1, date0, dt_force)
1059	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1060
1061	! first day of this year
1062	CALL ymds2ju (yy,1,1,0.0, julian0)
1063	!-----
1064	IF (check) THEN
1065	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1066	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1067	ENDIF
1068	ENDIF
1069	!-----
1070	CALL forcing_just_read (iim, jjm, zlev_n, ttm, itau_read_n, itau_read_n, &
1071	& swdown_n, rainf_n, snowf_n, tair_n, &
1072	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1073	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1074	& force_id, wind_N_exists, check)
1075	!---
1076	last_read = itau_read_n
1077	!-----
1078	!----- Compute mean solar angle for the comming period
1079	!-----
1080	IF (check) WRITE(numout,*) 'Going into solarang', split, one_day
1081	!-----
1082	IF ( dt_force .GT. 3600. ) THEN
1083	mean_sinang(:,:) = 0.0
1084	DO is=1,split
1085	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1086	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1087	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1088	CALL solarang (julian, julian0, iim, jjm, lon, lat, sinang)
1089	mean_sinang(:,:) = mean_sinang(:,:)+sinang(:,:)
1090	ENDDO
1091	mean_sinang(:,:) = mean_sinang(:,:)/split
1092	! WRITE(,) "mean_sinang =",MAXVAL(mean_sinang)
1093	ENDIF
1094	!-----
1095	ENDIF
1096	!---
1097	!--- Do the interpolation
1098	IF (check) WRITE(numout,*) 'Doing the interpolation between time steps'
1099	!---
1100	IF (split > 1) THEN
1101	rw = REAL(itau_split)/split
1102	ELSE
1103	rw = 1.
1104	ENDIF
1105	IF (check) WRITE(numout,*) 'Coeff of interpollation : ',rw
1106	!---
1107	qair(:,:) = (qair_n(:,:)-qair_nm1(:,:))*rw + qair_nm1(:,:)
1108	tair(:,:) = (tair_n(:,:)-tair_nm1(:,:))*rw + tair_nm1(:,:)
1109	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1110	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1111	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1112	IF (is_watchout) THEN
1113	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1114	zlevuv(:,:) = zlev(:,:)
1115	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1116	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1117	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1118	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1119	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1120	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1121	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1122	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1123	ENDIF
1124	!---
1125	!--- Here we need to allow for an option
1126	!--- where radiative energy is conserved
1127	!---
1128	IF ( netrad_cons ) THEN
1129	lwdown(:,:) = lwdown_n(:,:)
1130	ELSE
1131	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1132	ENDIF
1133	!---
1134	!--- For the solar radiation we decompose the mean value
1135	!--- using the zenith angle of the sun if the time step in the forcing data is
1136	!---- more than an hour. Else we use the standard linera interpolation
1137	!----
1138	IF (check) WRITE(numout,*) 'Ready to deal with the solar radiation'
1139	!----
1140	IF ( dt_force .GT. 3600. ) THEN
1141	!---
1142	IF ( netrad_cons ) THEN
1143	WRITE(numout,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
1144	ENDIF
1145	!---
1146	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1147	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1148	!!$ julian = julian_for + rw*dt_force/one_day
1149	IF (check) THEN
1150	WRITE(numout,'(a,f20.10,2I3)') &
1151	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1152	ENDIF
1153	!---
1154	CALL solarang(julian, julian0, iim, jjm, lon, lat, sinang)
1155	!---
1156	WHERE (mean_sinang(:,:) > 0.)
1157	swdown(:,:) = swdown_n(:,:) *sinang(:,:)/mean_sinang(:,:)
1158	ELSEWHERE
1159	swdown(:,:) = 0.0
1160	END WHERE
1161	!---
1162	WHERE (swdown(:,:) > 2000. )
1163	swdown(:,:) = 2000.
1164	END WHERE
1165	!---
1166	ELSE
1167	!!$ IF ( .NOT. is_watchout ) THEN
1168	!---
1169	IF ( netrad_cons ) THEN
1170	swdown(:,:) = swdown_n(:,:)
1171	ELSE
1172	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1173	ENDIF
1174	!---
1175	ENDIF
1176	!---
1177	IF (check) THEN
1178	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1179	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1180	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1181	IF (is_watchout) THEN
1182	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1183	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1184	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1185	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1186	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1187	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1188	ENDIF
1189	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1190	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1191	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1192	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1193	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1194	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1195	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1196	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1197	ENDIF
1198	!---
1199	!--- For precip we suppose that the rain
1200	!--- is the sum over the next 6 hours
1201	!---
1202	IF (itau_split <= nb_spread) THEN
1203	rainf(:,:) = rainf_n(:,:)*(split/nb_spread)
1204	snowf(:,:) = snowf_n(:,:)*(split/nb_spread)
1205	ELSE
1206	rainf(:,:) = 0.0
1207	snowf(:,:) = 0.0
1208	ENDIF
1209	IF (check) THEN
1210	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1211	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1212	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1213	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1214	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1215	ENDIF
1216	!---
1217	ENDIF
1218	!---
1219	!--- Here we might put the call to the weather generator ... one day.
1220	!--- Pour le moment, le branchement entre interpolation et generateur de temps
1221	!--- est fait au-dessus.
1222	!---
1223	!- IF ( initialized .AND. weathergen ) THEN
1224	!- ....
1225	!- ENDIF
1226	!---
1227	!--- At this point the code should be initialized. If not we have a problem !
1228	!---
1229	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
1230	!---
1231	initialized = .TRUE.
1232	!---
1233	ELSE
1234	IF ( .NOT. initialized ) THEN
1235	WRITE(numout,*) 'Why is the code forcing_read not initialized ?'
1236	WRITE(numout,*) 'Have you called it with both time-steps set to zero ?'
1237	STOP
1238	ENDIF
1239	ENDIF
1240	!
1241	!-------------------------
1242	END SUBROUTINE forcing_read_interpol
1243	!=====================================================================
1244	!-
1245	!=====================================================================
1246	SUBROUTINE forcing_just_read &
1247	& (iim, jjm, zlev, ttm, itb, ite, &
1248	& swdown, rainf, snowf, tair, &
1249	& u, v, qair, pb, lwdown, &
1250	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
1251	& force_id, wind_N_exists, check)
1252	!---------------------------------------------------------------------
1253	!- iim : Size of the grid in x
1254	!- jjm : size of the grid in y
1255	!- zlev : First Levels if it exists (ie if watchout file)
1256	!- ttm : number of time steps in all in the forcing file
1257	!- itb, ite : index of respectively begin and end of read for each variable
1258	!- swdown : Downward solar radiation (W/m^2)
1259	!- rainf : Rainfall (kg/m^2s)
1260	!- snowf : Snowfall (kg/m^2s)
1261	!- tair : 2m air temperature (K)
1262	!- u and v : 2m (in theory !) wind speed (m/s)
1263	!- qair : 2m humidity (kg/kg)
1264	!- pb : Surface pressure (Pa)
1265	!- lwdown : Downward long wave radiation (W/m^2)
1266	!-
1267	!- From a WATCHOUT file :
1268	!- SWnet : Net surface short-wave flux
1269	!- Eair : Air potential energy
1270	!- petAcoef : Coeficients A from the PBL resolution for T
1271	!- peqAcoef : Coeficients A from the PBL resolution for q
1272	!- petBcoef : Coeficients B from the PBL resolution for T
1273	!- peqBcoef : Coeficients B from the PBL resolution for q
1274	!- cdrag : Surface drag
1275	!- ccanopy : CO2 concentration in the canopy
1276	!- force_id : FLINCOM file id.
1277	!- It is used to close the file at the end of the run.
1278	!- wind_N_exists : if Wind_N and Wind_E are in the file (and not just Wind)
1279	!- check : Prompt for reading
1280	!---------------------------------------------------------------------
1281	IMPLICIT NONE
1282	!-
1283	INTEGER, INTENT(in) :: iim, jjm, ttm
1284	INTEGER, INTENT(in) :: itb, ite
1285	REAL, DIMENSION(iim,jjm), INTENT(out) :: zlev, &
1286	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown
1287	! for watchout files
1288	REAL, DIMENSION(iim,jjm), INTENT(out) :: &
1289	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
1290	INTEGER, INTENT(in) :: force_id
1291	! if Wind_N and Wind_E are in the file (and not just Wind)
1292	LOGICAL, INTENT(in) :: wind_N_exists
1293	LOGICAL :: check
1294	!-
1295	!---------------------------------------------------------------------
1296	CALL flinget (force_id,'Tair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1297	CALL forcing_zoom(data_full, tair)
1298	CALL flinget (force_id,'SWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1299	CALL forcing_zoom(data_full, swdown)
1300	CALL flinget (force_id,'LWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1301	CALL forcing_zoom(data_full, lwdown)
1302	CALL flinget (force_id,'Snowf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1303	CALL forcing_zoom(data_full, snowf)
1304	CALL flinget (force_id,'Rainf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1305	CALL forcing_zoom(data_full, rainf)
1306	!SZ FLUXNET input file correction
1307	! rainf=rainf/1800.
1308	!MM Rainf and not Snowf ?
1309	CALL flinget (force_id,'PSurf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1310	CALL forcing_zoom(data_full, pb)
1311	CALL flinget (force_id,'Qair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1312	CALL forcing_zoom(data_full, qair)
1313	!---
1314	IF ( wind_N_exists ) THEN
1315	CALL flinget (force_id,'Wind_N', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1316	CALL forcing_zoom(data_full, u)
1317	CALL flinget (force_id,'Wind_E', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1318	CALL forcing_zoom(data_full, v)
1319	ELSE
1320	CALL flinget (force_id,'Wind', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1321	CALL forcing_zoom(data_full, u)
1322	v=0.0
1323	ENDIF
1324	!----
1325	IF ( is_watchout ) THEN
1326	CALL flinget (force_id,'levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1327	CALL forcing_zoom(data_full, zlev)
1328	! Net surface short-wave flux
1329	CALL flinget (force_id,'SWnet', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1330	CALL forcing_zoom(data_full, SWnet)
1331	! Air potential energy
1332	CALL flinget (force_id,'Eair', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1333	CALL forcing_zoom(data_full, Eair)
1334	! Coeficients A from the PBL resolution for T
1335	CALL flinget (force_id,'petAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1336	CALL forcing_zoom(data_full, petAcoef)
1337	! Coeficients A from the PBL resolution for q
1338	CALL flinget (force_id,'peqAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1339	CALL forcing_zoom(data_full, peqAcoef)
1340	! Coeficients B from the PBL resolution for T
1341	CALL flinget (force_id,'petBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1342	CALL forcing_zoom(data_full, petBcoef)
1343	! Coeficients B from the PBL resolution for q
1344	CALL flinget (force_id,'peqBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1345	CALL forcing_zoom(data_full, peqBcoef)
1346	! Surface drag
1347	CALL flinget (force_id,'cdrag', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1348	CALL forcing_zoom(data_full, cdrag)
1349	! CO2 concentration in the canopy
1350	CALL flinget (force_id,'ccanopy', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1351	CALL forcing_zoom(data_full, ccanopy)
1352	ENDIF
1353	!
1354	!----
1355	IF (check) WRITE(numout,*) 'Variables have been extracted between ',itb,' and ',ite,' iterations of the forcing file.'
1356	!-------------------------
1357	END SUBROUTINE forcing_just_read
1358	!=====================================================================
1359	!-
1360	SUBROUTINE forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
1361	!---
1362	!--- This subroutine finds the indices of the land points over which the land
1363	!--- surface scheme is going to run.
1364	!---
1365	IMPLICIT NONE
1366	!-
1367	!- ARGUMENTS
1368	!-
1369	INTEGER, INTENT(IN) :: iim, jjm
1370	REAL, INTENT(IN) :: tair(iim,jjm)
1371	INTEGER, INTENT(OUT) :: i_test, j_test, nbindex
1372	INTEGER, INTENT(OUT) :: kindex(iim*jjm)
1373	LOGICAL :: check
1374	!-
1375	!- LOCAL
1376	INTEGER :: i, j, ig
1377	!-
1378	!-
1379	ig = 0
1380	i_test = 0
1381	j_test = 0
1382	!---
1383	IF (MINVAL(tair(:,:)) < 100.) THEN
1384	!----- In this case the 2m temperature is in Celsius
1385	DO j=1,jjm
1386	DO i=1,iim
1387	IF (tair(i,j) < 100.) THEN
1388	ig = ig+1
1389	kindex(ig) = (j-1)*iim+i
1390	!
1391	! Here we find at random a land-point on which we can do
1392	! some printouts for test.
1393	!
1394	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1395	i_test = i
1396	j_test = j
1397	IF (check) THEN
1398	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1399	ENDIF
1400	ENDIF
1401	ENDIF
1402	ENDDO
1403	ENDDO
1404	ELSE
1405	!----- 2m temperature is in Kelvin
1406	DO j=1,jjm
1407	DO i=1,iim
1408	IF (tair(i,j) < 500.) THEN
1409	ig = ig+1
1410	kindex(ig) = (j-1)*iim+i
1411	!
1412	! Here we find at random a land-point on which we can do
1413	! some printouts for test.
1414	!
1415	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1416	i_test = i
1417	j_test = j
1418	IF (check) THEN
1419	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
1420	ENDIF
1421	ENDIF
1422	ENDIF
1423	ENDDO
1424	ENDDO
1425	ENDIF
1426	!---
1427	nbindex = ig
1428	!---
1429	END SUBROUTINE forcing_landind
1430	!-
1431	!=====================================================================
1432	!-
1433	SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,lev,levuv,init_f)
1434	!-
1435	!- This subroutine calculates the longitudes and latitudes of the model grid.
1436	!-
1437	USE parallel
1438	IMPLICIT NONE
1439	!-
1440	INTEGER, INTENT(in) :: iim, jjm, llm
1441	LOGICAL, INTENT(in) :: init_f
1442	REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat
1443	REAL, DIMENSION(llm), INTENT(out) :: lev, levuv
1444	!-
1445	INTEGER :: i,j
1446	REAL :: zlev, wlev
1447	!-
1448	LOGICAL :: debug = .FALSE.
1449	!-
1450	!- Should be unified one day
1451	!-
1452	IF ( debug ) WRITE(numout,*) 'forcing_grid : options : ', weathergen, interpol
1453	!-
1454	!Config Key = HEIGHT_LEV1
1455	!Config Desc = Height at which T and Q are given
1456	!Config Def = 2.0
1457	!Config Help = The atmospheric variables (temperature and specific
1458	!Config humidity) are measured at a specific level.
1459	!Config The height of this level is needed to compute
1460	!Config correctly the turbulent transfer coefficients.
1461	!Config Look at the description of the forcing
1462	!Config DATA for the correct value.
1463	!-
1464	zlev = 2.0
1465	CALL getin_p('HEIGHT_LEV1', zlev)
1466	!-
1467	!Config Key = HEIGHT_LEVW
1468	!Config Desc = Height at which the wind is given
1469	!Config Def = 10.0
1470	!Config Help = The height at which wind is needed to compute
1471	!Config correctly the turbulent transfer coefficients.
1472	!-
1473	wlev = 10.0
1474	CALL getin_p('HEIGHT_LEVW', wlev)
1475	!-
1476	IF ( weathergen ) THEN
1477	IF (init_f) THEN
1478	DO i = 1, iim
1479	lon(i,:) = limit_west + merid_res/2. + &
1480	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim)
1481	ENDDO
1482	!-
1483	DO j = 1, jjm
1484	lat(:,j) = limit_north - zonal_res/2. - &
1485	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm)
1486	ENDDO
1487	ELSE
1488	IF (is_root_prc) THEN
1489	DO i = 1, iim_g
1490	lon_g(i,:) = limit_west + merid_res/2. + &
1491	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim_g)
1492	ENDDO
1493	!-
1494	DO j = 1, jjm_g
1495	lat_g(:,j) = limit_north - zonal_res/2. - &
1496	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm_g)
1497	ENDDO
1498	ELSE
1499	ALLOCATE(lon_g(iim_g, jjm_g), lat_g(iim_g, jjm_g))
1500	ENDIF
1501	CALL bcast(lon_g)
1502	CALL bcast(lat_g)
1503	lon=lon_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
1504	lat=lat_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
1505	ENDIF
1506	!-
1507	lev(:) = zlev
1508	levuv(:) = wlev
1509	!-
1510	ELSEIF ( interpol ) THEN
1511	!-
1512	CALL forcing_zoom(lon_full, lon)
1513	IF ( debug ) WRITE(numout,*) 'forcing_grid : out of zoom on lon'
1514	CALL forcing_zoom(lat_full, lat)
1515	IF ( debug ) WRITE(numout,*) 'forcing_grid : out of zoom on lat'
1516	!
1517	IF ( have_zaxis ) THEN
1518	lev(:) = lev_full(:)
1519	levuv(:) = lev_full(:)
1520	ELSE
1521	lev(:) = zlev
1522	levuv(:) = wlev
1523	ENDIF
1524	IF ( debug ) WRITE(numout,*) 'forcing_grid : levels : ', lev(:), levuv(:)
1525	!-
1526	ELSE
1527	!-
1528	STOP 'Neither weather generator nor temporal interpolation is specified.'
1529	!-
1530	ENDIF
1531	!-
1532	IF ( debug ) WRITE(numout,*) 'forcing_grid : ended'
1533	!-
1534	END SUBROUTINE forcing_grid
1535	!-
1536	!=====================================================================
1537	!-
1538	SUBROUTINE forcing_zoom(x_f, x_z)
1539	!-
1540	!- This subroutine takes the slab of data over which we wish to run the model.
1541	!-
1542	IMPLICIT NONE
1543	!-
1544	REAL, INTENT(IN) :: x_f(iim_full, jjm_full)
1545	REAL, INTENT(OUT) :: x_z(iim_zoom, jjm_zoom)
1546	!-
1547	INTEGER :: i,j
1548	!-
1549	DO i=1,iim_zoom
1550	DO j=1,jjm_zoom
1551	x_z(i,j) = x_f(i_index(i),j_index(j))
1552	ENDDO
1553	ENDDO
1554	!-
1555	END SUBROUTINE forcing_zoom
1556	!
1557	! ---------------------------------------------------------------------
1558	!
1559
1560	SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, &
1561	& iim_f, jjm_f, lon, lat, iim, jjm, iind, jind)
1562
1563	IMPLICIT NONE
1564	!
1565	! ARGUMENTS
1566	!
1567	REAL, INTENT(inout) :: limit_west,limit_east,limit_north,limit_south
1568	INTEGER, INTENT(in) :: iim_f, jjm_f
1569	REAL, INTENT(in) :: lon(iim_f, jjm_f), lat(iim_f, jjm_f)
1570	INTEGER, INTENT(out) :: iim,jjm
1571	INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f)
1572	!
1573	! LOCAL
1574	!
1575	INTEGER :: i, j
1576	REAL :: lolo
1577	LOGICAL :: over_dateline = .FALSE.
1578	!
1579	!
1580	IF ( ( ABS(limit_east) .GT. 180. ) .OR. &
1581	( ABS(limit_west) .GT. 180. ) ) THEN
1582	WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
1583	CALL ipslerr (3,'domain_size', &
1584	& 'Longitudes problem.','In run.def file :', &
1585	& 'limit_east > 180. or limit_west > 180.')
1586	ENDIF
1587	!
1588	IF ( limit_west .GT. limit_east ) over_dateline = .TRUE.
1589	!
1590	IF ( ( limit_south .LT. -90. ) .OR. &
1591	( limit_north .GT. 90. ) .OR. &
1592	( limit_south .GE. limit_north ) ) THEN
1593	WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
1594	CALL ipslerr (3,'domain_size', &
1595	& 'Latitudes problem.','In run.def file :', &
1596	& 'limit_south < -90. or limit_north > 90. or limit_south >= limit_north')
1597	ENDIF
1598	!
1599	! Here we assume that the grid of the forcing data is regular. Else we would have
1600	! to do more work to find the index table.
1601	!
1602	iim = 0
1603	DO i=1,iim_f
1604	!
1605	lolo = lon(i,1)
1606	IF ( lon(i,1) .GT. 180. ) lolo = lon(i,1) - 360.
1607	IF ( lon(i,1) .LT. -180. ) lolo = lon(i,1) + 360.
1608	!
1609	IF (lon(i,1) < limit_west) iim_g_begin = i+1
1610	IF (lon(i,1) < limit_east) iim_g_end = i
1611	!
1612	IF ( over_dateline ) THEN
1613	IF ( lolo .LE. limit_west .OR. lolo .GE. limit_east ) THEN
1614	iim = iim + 1
1615	iind(iim) = i
1616	ENDIF
1617	ELSE
1618	IF ( lolo .GE. limit_west .AND. lolo .LE. limit_east ) THEN
1619	iim = iim + 1
1620	iind(iim) = i
1621	ENDIF
1622	ENDIF
1623	!
1624	ENDDO
1625	!
1626	jjm = 0
1627	DO j=1,jjm_f
1628	IF (lat(1,j) > limit_north) jjm_g_begin = j+1
1629	IF (lat(1,j) > limit_south) jjm_g_end = j
1630	!
1631	IF ( lat(1,j) .GE. limit_south .AND. lat(1,j) .LE. limit_north) THEN
1632	jjm = jjm + 1
1633	jind(jjm) = j
1634	ENDIF
1635	ENDDO
1636	!
1637	WRITE(numout,*) 'Domain zoom size: iim, jjm = ', iim, jjm
1638	END SUBROUTINE domain_size
1639
1640	!------------------
1641	END MODULE readdim2

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

Download in other formats: