Context Navigation

source: tags/ORCHIDEE_1_9_5_1/ORCHIDEE_OL/readdim2.f90 @ 56

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

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

Download in other formats: