Context Navigation

source: tags/ORCHIDEE_1_9_5_2/ORCHIDEE_OL/readdim2.f90 @ 6268

Last change on this file since 6268 was 379, checked in by martial.mancip, 13 years ago
Correct calendars in all drivers.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 61.4 KB

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

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

Download in other formats: