Context Navigation

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

Last change on this file since 382 was 382, checked in by didier.solyga, 13 years ago
Merge the revisions 373 + 376 to 389 from the trunk in the externalized version
File size: 61.5 KB

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

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

Download in other formats: