Context Navigation

source: tags/ORCHIDEE_1_9_6/ORCHIDEE_OL/readdim2.f90 @ 3850

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

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

Download in other formats: