Context Navigation

source: branches/ORCHIDEE_EXT/ORCHIDEE_OL/readdim2.f90 @ 65

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

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

Download in other formats: