Context Navigation

readdim2.f90

Last change on this file was 7751, checked in by josefine.ghattas, 22 months ago
Bug correction for use of daily forcing, see ticket #861
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 105.5 KB

Line
1	! ==============================================================================================================================\n
2	! MODULE : readdim2
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module contains subroutines for reading the forcing file for the dim2_driver.
10	!!
11	!!
12	!!\n DESCRIPTION : This module contains subroutines for reading the forcing file for the dim2_driver.
13	!! Following subroutines are public and called from dim2_driver :
14	!! - forcing_info : Open the forcing file and return information about the grid in the forcing file.
15	!! Prepare for a zoom if needed.
16	!! Initialization of parallelization related to the grid.
17	!! - forcing_read : Return the forcing data for the current time step of the model. The forcing file will
18	!! be read if it has not already been done for the current time-step in the forcing file.
19	!! - forcing_grid : Calculate the longitudes and latitudes of the model grid.
20	!!
21	!! RECENT CHANGE(S): None
22	!!
23	!! REFERENCE(S) : None
24	!!
25	!! SVN :
26	!! $HeadURL$
27	!! $Date$
28	!! $Revision$
29	!! \n
30	!_ ================================================================================================================================
31
32	MODULE readdim2
33
34	USE ioipsl_para
35	USE weather
36	USE TIMER
37	USE constantes
38	USE time
39	USE solar
40	USE grid
41	USE mod_orchidee_para
42
43	IMPLICIT NONE
44
45	PRIVATE
46	PUBLIC :: forcing_read, forcing_info, forcing_grid
47
48	INTEGER, SAVE :: iim_full, jjm_full, llm_full, ttm_full
49	INTEGER, SAVE :: iim_zoom, jjm_zoom
50	INTEGER, SAVE :: iim_g_begin,jjm_g_begin,iim_g_end,jjm_g_end
51	REAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: data_full, lon_full, lat_full
52	REAL, SAVE, ALLOCATABLE, DIMENSION(:) :: lev_full
53	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: itau, i_index, j_index,j_index_g
54	INTEGER, SAVE :: i_test, j_test
55	INTEGER, SAVE :: printlev_loc !! Local printlev
56	LOGICAL, SAVE :: allow_weathergen, interpol, daily_interpol
57	LOGICAL, SAVE, PUBLIC :: weathergen, is_watchout
58	REAL, SAVE :: merid_res, zonal_res
59	LOGICAL, SAVE :: have_zaxis=.FALSE.
60	!-
61	!- Heigh controls and data
62	!-
63	LOGICAL, SAVE :: zfixed, zsigma, zhybrid, zlevels, zheight
64	LOGICAL, SAVE :: zsamelev_uv
65	REAL, SAVE :: zlev_fixed, zlevuv_fixed
66	REAL, SAVE :: zhybrid_a, zhybrid_b
67	REAL, SAVE :: zhybriduv_a, zhybriduv_b
68
69	CONTAINS
70
71	!! ==============================================================================================================================\n
72	!! SUBROUTINE : forcing_info
73	!!
74	!>\BRIEF Open the forcing file and return information about the grid in the forcing file.
75	!!
76	!!\n DESCRIPTION : This subroutine will get all the information from the forcing file and prepare for the zoom if needed.
77	!! It returns to the caller the sizes of the data it will receive at the forcing_read call.
78	!! This is important so that the caller can allocate the right space.
79	!!
80	!!
81	!! RECENT CHANGE(S): None
82	!!
83	!! MAIN OUTPUT VARIABLE(S):
84	!!
85	!! REFERENCE(S) :
86	!!
87	!_ ================================================================================================================================
88	SUBROUTINE forcing_info(filename, iim, jjm, llm, tm, date0, dt_force, force_id)
89
90	IMPLICIT NONE
91
92	!! 0.1 Input variables
93	CHARACTER(LEN=*), INTENT(in) :: filename !! Name of the file to be opened
94
95	!! 0.2 Output variables
96	INTEGER, INTENT(out) :: iim !! Size in x of the forcing data
97	INTEGER, INTENT(out) :: jjm !! Size in y of the forcing data
98	INTEGER, INTENT(out) :: llm !! Number of levels in the forcing data (not yet used)
99	INTEGER, INTENT(out) :: tm !! Time dimension of the forcing
100	REAL, INTENT(out) :: date0 !! The date at which the forcing file starts (julian days)
101	REAL, INTENT(out) :: dt_force !! Time-step of the forcing file in seconds
102	INTEGER, INTENT(out) :: force_id !! Id of the forcing file
103
104	!! 0.3 Local variables
105	CHARACTER(LEN=20) :: calendar_str
106	CHARACTER(LEN=200) :: printstr !! temporary character string to contain error message
107	REAL :: juld_1, juld_2
108	REAL, ALLOCATABLE, DIMENSION(:,:) :: fcontfrac
109	REAL, ALLOCATABLE, DIMENSION(:,:) :: qair
110	LOGICAL :: contfrac_exists=.FALSE.
111	INTEGER :: NbPoint
112	INTEGER :: i_test,j_test
113	INTEGER :: i,j,ind,ttm_part
114	INTEGER, ALLOCATABLE, DIMENSION(:) :: index_l
115	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
116	REAL, ALLOCATABLE, DIMENSION(:) :: lev, levuv
117
118	!-
119	CALL flininfo(filename, iim_full, jjm_full, llm_full, ttm_full, force_id)
120	!-
121	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_info : Details from forcing file :', &
122	iim_full, jjm_full, llm_full, ttm_full
123	!-
124	IF ( llm_full < 1 ) THEN
125	have_zaxis = .FALSE.
126	ELSE
127	have_zaxis = .TRUE.
128	ENDIF
129	IF ( printlev_loc>=3 ) WRITE(numout,*) 'have_zaxis : ', llm_full, have_zaxis
130	!-
131	ttm_part = 2
132	ALLOCATE(itau(ttm_part))
133	ALLOCATE(data_full(iim_full, jjm_full),lon_full(iim_full, jjm_full),&
134	& lat_full(iim_full, jjm_full))
135	ALLOCATE(lev_full(llm_full))
136	!-
137	lev_full(:) = zero
138	!-
139	dt_force=zero
140	CALL flinopen &
141	& (filename, .FALSE., iim_full, jjm_full, llm_full, lon_full, lat_full, &
142	& lev_full, ttm_part, itau, date0, dt_force, force_id)
143	IF ( dt_force == zero ) THEN
144	dt_force = itau(2) - itau(1)
145	itau(:) = itau(:) / dt_force
146	ENDIF
147	! WRITE(numout,*) "forcing_info : Forcing time step out of flinopen : ",dt_force
148	!-
149	!- What are the alowed options for the temportal interpolation
150	!-
151	!Config Key = ALLOW_WEATHERGEN
152	!Config Desc = Allow weather generator to create data
153	!Config If = [-]
154	!Config Def = n
155	!Config Help = This flag allows the forcing-reader to generate
156	!Config synthetic data if the data in the file is too sparse
157	!Config and the temporal resolution would not be enough to
158	!Config run the model.
159	!Config Units = [FLAG]
160	!-
161	allow_weathergen = .FALSE.
162	CALL getin_p('ALLOW_WEATHERGEN',allow_weathergen)
163	!-
164	!- The calendar was set by the forcing file. If no "calendar" attribute was
165	!- found then it is assumed to be gregorian,
166	!MM => FALSE !! it is NOT assumed anything !
167	!- else it is what ever is written in this attribute.
168	!-
169	CALL ioget_calendar(calendar_str)
170	i=INDEX(calendar_str,ACHAR(0))
171	IF ( i > 0 ) THEN
172	calendar_str(i:20)=' '
173	ENDIF
174	! WRITE(numout,*) "forcing_info : Calendar used : ",calendar_str
175	IF ( calendar_str == 'XXXX' ) THEN
176	IF (printlev_loc >= 1) WRITE(numout,*) "forcing_info : The calendar was not found in the forcing file."
177	IF (allow_weathergen) THEN
178	! Then change the calendar
179	CALL ioconf_calendar("noleap")
180	ELSE
181	IF ( printlev_loc>=1 ) WRITE(numout,*) "forcing_info : We will force it to gregorian by default."
182	CALL ioconf_calendar("gregorian") !! = 365.2425 ; "noleap" = 365.0; "360d"; "julian"=365.25
183	ENDIF
184	ENDIF
185	IF (printlev_loc >= 1) WRITE(numout,*) "readdim2 : Calendar used by the model : ",calendar_str
186	IF (ttm_full .GE. 2) THEN
187	juld_1 = itau2date(itau(1), date0, dt_force)
188	juld_2 = itau2date(itau(2), date0, dt_force)
189	ELSE
190	juld_1 = 0
191	juld_2 = 0
192	CALL ipslerr_p ( 3, 'forcing_info','What is that only one time step in the forcing file ?', &
193	& ' That can not be right.','verify forcing file.')
194	ENDIF
195	!-
196	!- Initialize one_year / one_day
197	CALL ioget_calendar (one_year, one_day)
198	!-
199	!- What is the distance between the two first states. From this we will deduce what is
200	!- to be done.
201	weathergen = .FALSE.
202	interpol = .FALSE.
203	daily_interpol = .FALSE.
204	is_watchout = .FALSE.
205	!-
206	IF ( ABS(ABS(juld_2-juld_1)-30.) .LE. 2.) THEN
207	IF ( allow_weathergen ) THEN
208	weathergen = .TRUE.
209	IF (printlev_loc >= 1) WRITE(numout,*) 'Using weather generator.'
210	ELSE
211	CALL ipslerr_p ( 3, 'forcing_info', &
212	& 'This seems to be a monthly file.', &
213	& 'We should use a weather generator with this file.', &
214	& 'This should be allowed in the run.def')
215	ENDIF
216	ELSEIF (( ABS(juld_1-juld_2) .LE. 1./4.) .OR. ( ABS(juld_1-juld_2) .EQ. 1.)) THEN
217	interpol = .TRUE.
218	IF (printlev_loc >= 1) WRITE(numout,*) 'We will interpolate between the forcing data time steps.'
219	IF ( ABS(juld_1-juld_2) .EQ. 1.) THEN
220	daily_interpol = .TRUE.
221	ENDIF
222	ELSE
223	! Using the weather generator with data other than monthly ones probably
224	! needs some thinking.
225	WRITE(numout,*) 'The time step is not suitable:',ABS(juld_1-juld_2),' days.'
226	CALL ipslerr_p ( 3, 'forcing_info','The time step is not suitable.', &
227	& '','We cannot do anything with these forcing data.')
228	ENDIF
229	!-
230	!- redefine the forcing time step if the weather generator is activated
231	!-
232	IF ( weathergen ) THEN
233	!Config Key = DT_WEATHGEN
234	!Config Desc = Calling frequency of weather generator
235	!Config If = ALLOW_WEATHERGEN
236	!Config Def = 1800.
237	!Config Help = Determines how often the weather generator
238	!Config is called (time step in s). Should be equal
239	!Config to or larger than Sechiba's time step (say,
240	!Config up to 6 times Sechiba's time step or so).
241	!Config Units = [seconds]
242	dt_force = 1800.
243	CALL getin_p('DT_WEATHGEN',dt_force)
244	ENDIF
245	!-
246	!- Define the zoom
247	!-
248	!Config Key = LIMIT_WEST
249	!Config Desc = Western limit of region
250	!Config If = [-]
251	!Config Def = -180.
252	!Config Help = Western limit of the region we are
253	!Config interested in. Between -180 and +180 degrees
254	!Config The model will use the smalest regions from
255	!Config region specified here and the one of the forcing file.
256	!Config Units = [Degrees]
257	!-
258	limit_west = -180.
259	CALL getin_p('LIMIT_WEST',limit_west)
260	!-
261	!Config Key = LIMIT_EAST
262	!Config Desc = Eastern limit of region
263	!Config If = [-]
264	!Config Def = 180.
265	!Config Help = Eastern limit of the region we are
266	!Config interested in. Between -180 and +180 degrees
267	!Config The model will use the smalest regions from
268	!Config region specified here and the one of the forcing file.
269	!Config Units = [Degrees]
270	!-
271	limit_east = 180.
272	CALL getin_p('LIMIT_EAST',limit_east)
273	!-
274	!Config Key = LIMIT_NORTH
275	!Config Desc = Northern limit of region
276	!Config If = [-]
277	!Config Def = 90.
278	!Config Help = Northern limit of the region we are
279	!Config interested in. Between +90 and -90 degrees
280	!Config The model will use the smalest regions from
281	!Config region specified here and the one of the forcing file.
282	!Config Units = [Degrees]
283	!-
284	limit_north = 90.
285	CALL getin_p('LIMIT_NORTH',limit_north)
286	!-
287	!Config Key = LIMIT_SOUTH
288	!Config Desc = Southern limit of region
289	!Config If = [-]
290	!Config Def = -90.
291	!Config Help = Southern limit of the region we are
292	!Config interested in. Between 90 and -90 degrees
293	!Config The model will use the smalest regions from
294	!Config region specified here and the one of the forcing file.
295	!Config Units = [Degrees]
296	!-
297	limit_south = -90.
298	CALL getin_p('LIMIT_SOUTH',limit_south)
299	!-
300	!- Calculate domain size
301	!-
302
303	IF ( interpol ) THEN
304	!-
305	!- If we use temporal interpolation, then we cannot change the resolution (yet?)
306	!-
307	ALLOCATE(i_index(iim_full), j_index(jjm_full),j_index_g(jjm_full))
308	IF (is_root_prc) THEN
309
310	CALL domain_size (limit_west, limit_east, limit_north, limit_south,&
311	& iim_full, jjm_full, lon_full, lat_full, iim_zoom, jjm_zoom,&
312	& i_index, j_index_g)
313
314	j_index(:)=j_index_g(:)
315
316	ALLOCATE(qair(iim_full,jjm_full))
317	CALL flinget_buffer (force_id,'Qair',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
318	CALL forcing_zoom(data_full, qair)
319
320	ALLOCATE(fcontfrac(iim_zoom,jjm_zoom))
321	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
322	IF ( contfrac_exists ) THEN
323	IF (printlev_loc >= 1) WRITE(numout,*) "contfrac exist in the forcing file."
324	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, 1, ttm_full, 1, 1, data_full)
325	CALL forcing_zoom(data_full, fcontfrac)
326	IF (printlev_loc >= 2) WRITE(numout,*) "fcontfrac min/max :", &
327	MINVAL(fcontfrac(1:iim_zoom,1:jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,1:jjm_zoom))
328	ELSE
329	fcontfrac(:,:)=1.
330	ENDIF
331
332	DO i=1,iim_zoom
333	DO j=1,jjm_zoom
334	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
335	qair(i,j) = 999999.
336	ENDIF
337	ENDDO
338	ENDDO
339
340	DEALLOCATE(fcontfrac)
341
342	ALLOCATE(index_l(iim_zoom*jjm_zoom))
343	!- In this point is returning the global NbPoint with the global index
344	CALL forcing_landind(iim_zoom,jjm_zoom,qair,NbPoint,index_l,i_test,j_test)
345	!
346	! Work out the vertical layers to be used
347	!
348	CALL forcing_vertical_ioipsl(force_id)
349	ELSE
350	ALLOCATE(index_l(1))
351	ENDIF
352
353	! Initiate global grid and parallelism
354	CALL bcast(iim_zoom)
355	CALL bcast(jjm_zoom)
356	CALL bcast(NbPoint)
357	CALL grid_set_glo(iim_zoom,jjm_zoom,NbPoint)
358	CALL grid_allocate_glo(4)
359
360	! Check consistency in the number of mpi processors and the number of land points
361	! in order to prevent an exception
362	IF (NbPoint < mpi_size) THEN
363	WRITE(printstr,*) 'The number of landpoints found (', NbPoint, &
364	') is less than the number of processors selected (', mpi_size,')'
365	CALL ipslerr_p(3, 'forcing_info', 'Wrong parallelization options', &
366	TRIM(printstr), &
367	'Decrease the number of processors for the current grid')
368	ENDIF
369
370	!
371	!- global index index_g is the index_l of root proc
372	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
373
374	DEALLOCATE(index_l)
375
376	!
377	! Distribute to all processors the information on the forcing
378	!
379	CALL bcast(index_g)
380	CALL Init_orchidee_data_para_driver(nbp_glo,index_g)
381	CALL init_ioipsl_para
382
383	! Initialize printlev_loc
384	printlev_loc=get_printlev('readdim2')
385	IF (printlev_loc >= 2) WRITE(numout,*) 'Standard PRINTLEV= ', printlev
386	IF (printlev_loc >= 2) WRITE(numout,*) 'Local PRINTLEV_readdim2= ', printlev_loc
387
388	! CALL Init_writeField_p
389
390	CALL bcast(jjm_zoom)
391	CALL bcast(i_index)
392	CALL bcast(j_index_g)
393	CALL bcast(zfixed)
394	CALL bcast(zsigma)
395	CALL bcast(zhybrid)
396	CALL bcast(zlevels)
397	CALL bcast(zheight)
398	CALL bcast(zsamelev_uv)
399	CALL bcast(zlev_fixed)
400	CALL bcast(zlevuv_fixed)
401	CALL bcast(zhybrid_a)
402	CALL bcast(zhybrid_b)
403	CALL bcast(zhybriduv_a)
404	CALL bcast(zhybriduv_b)
405	ind=0
406	DO j=1,jjm_zoom
407	IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
408	ind=ind+1
409	j_index(ind)=j_index_g(j)
410	ENDIF
411	ENDDO
412
413	jjm_zoom=jj_nb
414	iim_zoom=iim_g
415
416	!-
417	!- If we use the weather generator, then we read zonal and meridional resolutions
418	!- This should be unified one day...
419	!-
420	ELSEIF ( weathergen ) THEN
421	!-
422	!Config Key = MERID_RES
423	!Config Desc = North-South Resolution
424	!Config Def = 2.
425	!Config If = ALLOW_WEATHERGEN
426	!Config Help = North-South Resolution of the region we are
427	!Config interested in.
428	!Config Units = [Degrees]
429	!-
430	merid_res = 2.
431	CALL getin_p('MERID_RES',merid_res)
432	!-
433	!Config Key = ZONAL_RES
434	!Config Desc = East-West Resolution
435	!Config Def = 2.
436	!Config If = ALLOW_WEATHERGEN
437	!Config Help = East-West Resolution of the region we are
438	!Config interested in. In degrees
439	!Config Units = [Degrees]
440	!-
441	zonal_res = 2.
442	CALL getin_p('ZONAL_RES',zonal_res)
443	!-
444	!- Number of time steps is meaningless in this case
445	!-
446	! ttm_full = HUGE( ttm_full )
447	!MM Number (realistic) of time steps for half hour dt
448	ttm_full = NINT(one_year * 86400. / dt_force)
449	!-
450	IF (is_root_prc) THEN
451
452	!- In this point is returning the global NbPoint with the global index
453	CALL weathgen_domain_size (limit_west,limit_east,limit_north,limit_south, &
454	zonal_res,merid_res,iim_zoom,jjm_zoom)
455	ALLOCATE(index_l(iim_zoom*jjm_zoom))
456	ENDIF
457	CALL bcast(iim_zoom)
458	CALL bcast(jjm_zoom)
459
460	ALLOCATE(lon(iim_zoom,jjm_zoom))
461	ALLOCATE(lat(iim_zoom,jjm_zoom))
462	ALLOCATE(lev(llm_full),levuv(llm_full))
463
464	! We need lon and lat now for weathgen_init
465	CALL forcing_grid (iim_zoom,jjm_zoom,llm_full,lon,lat,init_f=.TRUE.)
466	CALL forcing_vertical_ioipsl(-1)
467
468	IF (is_root_prc) THEN
469	CALL weathgen_init &
470	& (filename,dt_force,force_id,iim_zoom,jjm_zoom, &
471	& zonal_res,merid_res,lon,lat,index_l,NbPoint)
472	!!$,&
473	!!$ & i_index,j_index_g)
474	ELSE
475	ALLOCATE(index_l(1))
476	index_l(1)=1
477	ENDIF
478
479	CALL bcast(NbPoint)
480	CALL grid_set_glo(iim_zoom,jjm_zoom,NbPoint)
481	CALL grid_allocate_glo(4)
482
483	!
484	!- global index index_g is the index_l of root proc
485	IF (is_root_prc) index_g(:)=index_l(1:NbPoint)
486
487	DEALLOCATE(index_l)
488
489	CALL bcast(index_g)
490	CALL Init_orchidee_data_para_driver(nbp_glo,index_g)
491	CALL init_ioipsl_para
492	! CALL Init_writeField_p
493	CALL bcast(jjm_zoom)
494	!!$ CALL bcast(i_index)
495	!!$ CALL bcast(j_index_g)
496
497	!!$ ind=0
498	!!$ DO j=1,jjm_zoom
499	!!$ IF ( (j >= jj_begin) .AND. (j <= jj_end) ) THEN
500	!!$ ind=ind+1
501	!!$ j_index(ind)=j_index_g(j)
502	!!$ ENDIF
503	!!$ ENDDO
504
505	jjm_zoom=jj_nb
506	iim_zoom=iim_g
507	!
508	CALL weathgen_read_file(force_id,iim_zoom,jjm_zoom)
509
510	!-
511	ELSE
512	!-
513	CALL ipslerr_p(3,'forcing_info','Neither interpolation nor weather generator is specified.','','')
514	!-
515	ENDIF
516	!-
517	!- Transfer the right information to the caller
518	!-
519	iim = iim_zoom
520	jjm = jjm_zoom
521	llm = 1
522	tm = ttm_full
523	!-
524	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_info : end : ', iim,jjm, llm,tm
525	!-
526	END SUBROUTINE forcing_info
527
528
529	!! ==============================================================================================================================\n
530	!! SUBROUTINE : forcing_read
531	!!
532	!>\BRIEF Return forcing data for the current time step
533	!!
534	!!\n DESCRIPTION : Return the forcing data for the current time step of the model. The forcing file will
535	!! be read if it has not already been done for the current time-step in the forcing file.
536	!!
537	!! RECENT CHANGE(S): None
538	!!
539	!! MAIN OUTPUT VARIABLE(S):
540	!!
541	!! REFERENCE(S) :
542	!!
543	!_ ================================================================================================================================
544	SUBROUTINE forcing_read &
545	& (filename, rest_id, lrstread, lrstwrite, &
546	& itauin, istp, itau_split, split, nb_spread, lwdown_cons, swdown_cons, date0, &
547	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
548	& swdown, coszang, precip, snowf, tair, u, v, qair, pb, lwdown, &
549	& fcontfrac, fneighbours, fresolution, &
550	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
551	& kindex, nbindex, force_id)
552
553	IMPLICIT NONE
554
555	!! 0. Variable and parameter declaration
556	!! 0.1 Input variables
557	CHARACTER(LEN=*), INTENT(IN) :: filename !! name of the file to be opened
558	INTEGER, INTENT(IN) :: force_id !! FLINCOM file id. It is used to close the file at the end of the run.
559	INTEGER, INTENT(IN) :: rest_id !! ID of restart file
560	LOGICAL, INTENT(IN) :: lrstread !! read restart file?
561	LOGICAL, INTENT(IN) :: lrstwrite !! write restart file?
562	INTEGER, INTENT(IN) :: itauin !! time step for which we need the data
563	INTEGER, INTENT(IN) :: istp !! time step for restart file
564	INTEGER, INTENT(IN) :: itau_split !! Current step between 2 forcing times-step (it decides if it is time to read)
565	INTEGER, INTENT(IN) :: split !! The number of time steps between two time-steps of the forcing
566	INTEGER, INTENT(IN) :: nb_spread !! Over how many time steps do we spread the precipitation
567	LOGICAL, INTENT(IN) :: lwdown_cons !! Flag to conserve lwdown radiation from forcing
568	LOGICAL, INTENT(IN) :: swdown_cons !! Flag to conserve swdown radiation from forcing
569	REAL, INTENT(IN) :: date0 !! The date at which the forcing file starts (julian days)
570	REAL, INTENT(IN) :: dt_force !! time-step of the forcing file in seconds
571	INTEGER, INTENT(IN) :: iim !! Size of the grid in x
572	INTEGER, INTENT(IN) :: jjm !! Size of the grid in y
573	INTEGER, INTENT(IN) :: ttm !! number of time steps in all in the forcing file
574	REAL,DIMENSION(iim,jjm), INTENT(IN) :: lon !! Longitudes
575	REAL,DIMENSION(iim,jjm), INTENT(IN) :: lat !! Latitudes
576
577	!! 0.2 Output variables
578	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: zlev !! First Levels if it exists (ie if watchout file)
579	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: zlevuv !! First Levels of the wind (equal precedent, if it exists)
580	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: swdown !! Downward solar radiation (W/m^2)
581	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: coszang !! Cosine of the solar zenith angle (unitless)
582	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: precip !! Precipitation (Rainfall) (kg/m^2s)
583	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: snowf !! Snowfall (kg/m^2s)
584	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: tair !! 1st level (2m ? in off-line) air temperature (K)
585	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: u !! 1st level (2m/10m ? in off-line) (in theory !) wind speed (m/s)
586	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: v !! 1st level (2m/10m ? in off-line) (in theory !) wind speed (m/s)
587	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: qair !! 1st level (2m ? in off-line) humidity (kg/kg)
588	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: pb !! Surface pressure (Pa)
589	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: lwdown !! Downward long wave radiation (W/m^2)
590	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: fcontfrac !! Continental fraction (no unit)
591	REAL,DIMENSION(iim,jjm,2), INTENT(OUT) :: fresolution !! resolution in x and y dimensions for each points
592	INTEGER,DIMENSION(iim,jjm,8), INTENT(OUT) :: fneighbours !! land neighbours
593
594	!! From a WATCHOUT file :
595	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: SWnet !! Net surface short-wave flux
596	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: Eair !! Air potential energy
597	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: petAcoef !! Coeficients A from the PBL resolution for T
598	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: peqAcoef !! Coeficients A from the PBL resolution for q
599	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: petBcoef !! Coeficients B from the PBL resolution for T
600	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: peqBcoef !! Coeficients B from the PBL resolution for q
601	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: cdrag !! Surface drag
602	REAL,DIMENSION(iim,jjm), INTENT(OUT) :: ccanopy !! CO2 concentration in the canopy
603
604	!! 0.3 Modified variable
605	INTEGER, INTENT(INOUT) :: nbindex !! Number of land points
606	INTEGER,DIMENSION(iim*jjm), INTENT(INOUT) :: kindex !! Index of all land-points in the data (used for the gathering)
607
608	!! 0.4 Local variables
609	INTEGER :: ik,i,j
610
611	IF ( interpol ) THEN
612
613	CALL forcing_read_interpol &
614	(filename, itauin, itau_split, split, nb_spread, lwdown_cons, swdown_cons, date0, &
615	dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, &
616	swdown, coszang, precip, snowf, tair, u, v, qair, pb, lwdown, &
617	fcontfrac, fneighbours, fresolution, &
618	SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
619	kindex, nbindex, force_id)
620
621	ELSEIF ( weathergen ) THEN
622
623	IF (lrstread) THEN
624	fcontfrac(:,:) = 1.0
625	IF (printlev_loc >= 2) WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
626	ENDIF
627
628	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
629	CALL weathgen_main (istp, istp, filename, force_id, iim, jjm, 1, &
630	rest_id, lrstread, lrstwrite, &
631	limit_west, limit_east, limit_north, limit_south, &
632	zonal_res, merid_res, lon, lat, date0, dt_force, &
633	kindex, nbindex, &
634	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
635	ELSE
636	CALL weathgen_main (itauin, istp, filename, force_id, iim, jjm, 1, &
637	rest_id, lrstread, lrstwrite, &
638	limit_west, limit_east, limit_north, limit_south, &
639	zonal_res, merid_res, lon, lat, date0, dt_force, &
640	kindex, nbindex, &
641	swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
642	ENDIF
643
644	IF ( (itauin == 0).AND.(itau_split == 0) ) THEN
645	!---
646	!--- Allocate grid stuff
647	!---
648	CALL grid_init ( nbindex, 4, regular_lonlat, "ForcingGrid" )
649	!---
650	!--- Compute
651	!---
652	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, index_g)
653	!CALL grid_stuff (nbindex, iim, jjm, lon, lat, kindex)
654	DO ik=1,nbindex
655
656	j = ((kindex(ik)-1)/iim) + 1
657	i = (kindex(ik) - (j-1)*iim)
658	!-
659	!- Store variable to help describe the grid
660	!- once the points are gathered.
661	!-
662	fneighbours(i,j,:) = neighbours(ik,:)
663	!
664	fresolution(i,j,:) = resolution(ik,:)
665	ENDDO
666	ENDIF
667	ELSE
668
669	CALL ipslerr_p(3,'forcing_read','Neither interpolation nor weather generator is specified.','','')
670
671	ENDIF
672
673	IF (.NOT. is_watchout) THEN
674	! We have to compute Potential air energy
675	WHERE(tair(:,:) < val_exp)
676	eair(:,:) = cp_airtair(:,:)+cte_gravzlev(:,:)
677	ENDWHERE
678	ENDIF
679
680	END SUBROUTINE forcing_read
681
682	!! ==============================================================================================================================\n
683	!! SUBROUTINE : forcing_read_interpol
684	!!
685	!>\BRIEF
686	!!
687	!!\n DESCRIPTION :
688	!!
689	!! RECENT CHANGE(S): None
690	!!
691	!! MAIN OUTPUT VARIABLE(S):
692	!!
693	!! REFERENCE(S) :
694	!!
695	!_ ================================================================================================================================
696	SUBROUTINE forcing_read_interpol &
697	& (filename, itauin, itau_split, split, nb_spread, lwdown_cons, swdown_cons, date0, &
698	& dt_force, iim, jjm, lon, lat, zlev, zlevuv, ttm, swdown, coszang, rainf, snowf, tair, &
699	& u, v, qair, pb, lwdown, &
700	& fcontfrac, fneighbours, fresolution, &
701	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
702	& kindex, nbindex, force_id)
703	!---------------------------------------------------------------------
704	!- filename : name of the file to be opened
705	!- itauin : time step for which we need the data
706	!- itau_split : Where are we within the splitting
707	!- of the time-steps of the forcing files
708	!- (it decides IF we READ or not)
709	!- split : The number of time steps we do
710	!- between two time-steps of the forcing
711	!- nb_spread : Over how many time steps do we spread the precipitation
712	!- lwdown_cons: flag that decides if lwdown radiation should be conserved.
713	!- swdown_cons: flag that decides if swdown radiation should be conserved.
714	!- date0 : The date at which the forcing file starts (julian days)
715	!- dt_force : time-step of the forcing file in seconds
716	!- iim : Size of the grid in x
717	!- jjm : size of the grid in y
718	!- lon : Longitudes
719	!- lat : Latitudes
720	!- zlev : First Levels if it exists (ie if watchout file)
721	!- zlevuv : First Levels of the wind (equal precedent, if it exists)
722	!- ttm : number of time steps in all in the forcing file
723	!- swdown : Downward solar radiation (W/m^2)
724	!- coszang : Cosine of the solar zenith angle (unitless)
725	!- rainf : Rainfall (kg/m^2s)
726	!- snowf : Snowfall (kg/m^2s)
727	!- tair : 2m air temperature (K)
728	!- u and v : 2m (in theory !) wind speed (m/s)
729	!- qair : 2m humidity (kg/kg)
730	!- pb : Surface pressure (Pa)
731	!- lwdown : Downward long wave radiation (W/m^2)
732	!- fcontfrac : Continental fraction (no unit)
733	!- fneighbours: land neighbours
734	!- fresolution: resolution in x and y dimensions for each points
735	!-
736	!- From a WATCHOUT file :
737	!- SWnet : Net surface short-wave flux
738	!- Eair : Air potential energy
739	!- petAcoef : Coeficients A from the PBL resolution for T
740	!- peqAcoef : Coeficients A from the PBL resolution for q
741	!- petBcoef : Coeficients B from the PBL resolution for T
742	!- peqBcoef : Coeficients B from the PBL resolution for q
743	!- cdrag : Surface drag
744	!- ccanopy : CO2 concentration in the canopy
745	!-
746	!- kindex : Index of all land-points in the data
747	!- (used for the gathering)
748	!- nbindex : Number of land points
749	!- force_id : FLINCOM file id.
750	!- It is used to close the file at the end of the run.
751	!---------------------------------------------------------------------
752	IMPLICIT NONE
753	!-
754	INTEGER,PARAMETER :: lm=1
755	!-
756	!- Input variables
757	!-
758	CHARACTER(LEN=*) :: filename
759	INTEGER :: itauin, itau_split, split, nb_spread
760	LOGICAL, INTENT(IN) :: lwdown_cons, swdown_cons
761	REAL :: date0, dt_force
762	INTEGER :: iim, jjm, ttm
763	REAL,DIMENSION(:,:),INTENT(IN) :: lon, lat !- LOCAL data array (size=iim,jjm)
764	INTEGER, INTENT(IN) :: force_id
765	!-
766	!- Output variables
767	!-
768	REAL,DIMENSION(:,:),INTENT(OUT) :: zlev, zlevuv, & !- LOCAL data array (size=iim,jjm)
769	& swdown, coszang, rainf, snowf, tair, u, v, qair, pb, lwdown, &
770	& fcontfrac
771	REAL,DIMENSION(:,:,:),INTENT(OUT) :: fresolution !- LOCAL data array (size=iim,jjm,2)
772	INTEGER,DIMENSION(:,:,:),INTENT(OUT) :: fneighbours !- LOCAL data array (size=iim,jjm,8)
773	! for watchout files
774	REAL,DIMENSION(:,:),INTENT(OUT) :: &
775	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
776	INTEGER,DIMENSION(:),INTENT(INOUT) :: kindex !- LOCAL index of the map
777	INTEGER, INTENT(INOUT) :: nbindex
778	!-
779	!- Local variables
780	!-
781	INTEGER, SAVE :: last_read=0
782	INTEGER, SAVE :: itau_read, itau_read_nm1=0, itau_read_n=0
783	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
784	& zlev_nm1, zlevuv_nm1, swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, u_nm1, v_nm1, qair_nm1, &
785	& pb_nm1, lwdown_nm1, &
786	& zlev_n, zlevuv_n, swdown_n, rainf_n, snowf_n, tair_n, u_n, v_n, qair_n, &
787	& pb_n, lwdown_n, mean_coszang
788
789	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
790	& startday_n, startday_nm1, daylength_n, daylength_nm1, tmax_n, tmax_nm1, tmin_nm1, tmin_nm2, tmin_n, &
791	& qsatta, qsattmin_n, qsattmin_nm1, qmin_n, qmin_nm1, qmax_n, qmax_nm1, qsa
792	REAL,SAVE :: hour
793
794	! just for grid stuff if the forcing file is a watchout file
795	REAL, ALLOCATABLE, DIMENSION(:,:) :: tmpdata
796	! variables to be read in watchout files
797	REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
798	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
799	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n
800	REAL, SAVE :: julian_for ! Date of the forcing to be read
801	REAL :: julian, ss, rw
802	!jur,
803	REAL, SAVE :: julian0 ! First day of this year
804	INTEGER :: yy, mm, dd, is, i, j, ik
805	REAL(r_std), DIMENSION(2) :: min_resol, max_resol
806	! if Wind_N and Wind_E are in the file (and not just Wind)
807	LOGICAL, SAVE :: wind_N_exists=.FALSE.
808	LOGICAL :: wind_E_exists=.FALSE.
809	LOGICAL, SAVE :: contfrac_exists=.FALSE.
810	LOGICAL, SAVE :: neighbours_exists=.FALSE.
811	LOGICAL, SAVE :: initialized = .FALSE.
812	! to bypass FPE problem on integer convertion with missing_value heigher than precision
813	INTEGER, PARAMETER :: undef_int = 999999999
814	!---------------------------------------------------------------------
815
816	itau_read = MOD((itauin-1),ttm)+1
817
818	IF (printlev_loc >= 5) THEN
819	WRITE(numout,*) &
820	" FORCING READ : itauin, itau_read, itau_split : ",&
821	itauin, itau_read, itau_split
822	ENDIF
823
824	!-
825	!- This part initializes the reading of the forcing. As you can see
826	!- we only go through here if both time steps are zero.
827	!-
828	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
829	!-
830	!- Tests on forcing file type
831	CALL flinquery_var(force_id, 'Wind_N', wind_N_exists)
832	IF ( wind_N_exists ) THEN
833	CALL flinquery_var(force_id, 'Wind_E', wind_E_exists)
834	IF ( .NOT. wind_E_exists ) THEN
835	CALL ipslerr_p(3,'forcing_read_interpol', &
836	& 'Variable Wind_E does not exist in forcing file', &
837	& 'But variable Wind_N exists.','Please, rename Wind_N to Wind;')
838	ENDIF
839	ENDIF
840	CALL flinquery_var(force_id, 'levels', is_watchout)
841	IF ( is_watchout ) THEN
842	WRITE(numout,*) "Read a Watchout File."
843	ENDIF
844	CALL flinquery_var(force_id, 'contfrac', contfrac_exists)
845	!-
846	IF (printlev_loc >= 5) WRITE(numout,*) 'ALLOCATE all the memory needed'
847	!-
848	ALLOCATE &
849	& (swdown_nm1(iim,jjm), rainf_nm1(iim,jjm), snowf_nm1(iim,jjm), &
850	& tair_nm1(iim,jjm), u_nm1(iim,jjm), v_nm1(iim,jjm), qair_nm1(iim,jjm), &
851	& pb_nm1(iim,jjm), lwdown_nm1(iim,jjm))
852	ALLOCATE &
853	& (swdown_n(iim,jjm), rainf_n(iim,jjm), snowf_n(iim,jjm), &
854	& tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), &
855	& pb_n(iim,jjm), lwdown_n(iim,jjm))
856
857	IF(daily_interpol) THEN
858	ALLOCATE &
859	& (startday_n(iim,jjm), startday_nm1(iim,jjm), daylength_n(iim,jjm), &
860	& daylength_nm1(iim,jjm), tmax_n(iim,jjm), tmax_nm1(iim,jjm), tmin_n(iim,jjm), &
861	& tmin_nm1(iim,jjm), tmin_nm2(iim,jjm), qsatta(iim,jjm), qsattmin_n(iim,jjm), qsattmin_nm1(iim,jjm), &
862	& qmin_n(iim,jjm), qmin_nm1(iim,jjm), qmax_n(iim,jjm), qmax_nm1(iim,jjm), qsa(iim,jjm) )
863	ENDIF
864
865
866	ALLOCATE &
867	& (zlev_nm1(iim,jjm), zlev_n(iim,jjm), zlevuv_nm1(iim,jjm), zlevuv_n(iim,jjm), &
868	& SWnet_nm1(iim,jjm), Eair_nm1(iim,jjm), cdrag_nm1(iim,jjm), ccanopy_nm1(iim,jjm), &
869	& petAcoef_nm1(iim,jjm), peqAcoef_nm1(iim,jjm), petBcoef_nm1(iim,jjm), peqBcoef_nm1(iim,jjm), &
870	& SWnet_n(iim,jjm), Eair_n(iim,jjm), cdrag_n(iim,jjm), ccanopy_n(iim,jjm), &
871	& petAcoef_n(iim,jjm), peqAcoef_n(iim,jjm), petBcoef_n(iim,jjm), peqBcoef_n(iim,jjm))
872	ALLOCATE &
873	& (mean_coszang(iim,jjm))
874	!-
875	IF (printlev_loc >= 5) WRITE(numout,*) 'Memory ALLOCATED'
876	!-
877	!- We need for the driver surface air temperature and humidity before the
878	!- the loop starts. Thus we read it here.
879	!-
880	CALL forcing_just_read (iim, jjm, zlev, zlevuv, ttm, 1, 1, &
881	& swdown, rainf, snowf, tair, &
882	& u, v, qair, pb, lwdown, &
883	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
884	& force_id, wind_N_exists)
885	!----
886
887	!-- First in case it's not a watchout file
888	IF ( .NOT. is_watchout ) THEN
889	IF ( contfrac_exists ) THEN
890	IF (printlev_loc >= 1) WRITE(numout,*) "contfrac exist in the forcing file."
891	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
892	CALL forcing_zoom(data_full, fcontfrac)
893	IF (printlev_loc >= 2) WRITE(numout,*) "fcontfrac min/max :", &
894	MINVAL(fcontfrac(1:iim_zoom,jjm_zoom)),MAXVAL(fcontfrac(1:iim_zoom,jjm_zoom))
895	!
896	! We need to make sure that when we gather the points we pick all
897	! the points where contfrac is above 0. Thus we prepare tair for
898	! subroutine forcing_landind
899	!
900	DO i=1,iim
901	DO j=1,jjm
902	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0. ! bande de recouvrement du scatter2D
903	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0. ! => on mets les données qu'on veut pas du noeud à missing_value
904	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
905	tair(i,j) = 999999.
906	ENDIF
907	ENDDO
908	ENDDO
909	ELSE
910	fcontfrac(:,:) = 1.0
911	ENDIF
912	!---
913	!--- Create the index table
914	!---
915	!--- This job return a LOCAL kindex
916	CALL forcing_landind(iim, jjm, tair, nbindex, kindex, i_test, j_test)
917	#ifdef CPP_PARA
918	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
919	! Force nbindex points for parallel computation
920	nbindex=nbp_loc
921	CALL scatter(index_g,kindex(1:nbindex))
922	kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
923	#endif
924	ik=MAX(nbindex/2,1)
925	j_test = (((kindex(ik)-1)/iim) + 1)
926	i_test = (kindex(ik) - (j_test-1)*iim)
927	IF (printlev_loc >= 5) THEN
928	WRITE(numout,*) 'New test point is : ', i_test, j_test
929	ENDIF
930	!---
931	!--- Allocate grid stuff
932	!---
933	CALL grid_init ( nbindex, 4, regular_lonlat, "ForcingGrid" )
934	!---
935	!--- All grid variables
936	!---
937	CALL grid_stuff(nbp_glo, iim_g, jjm_g, lon_g, lat_g, index_g)
938	DO ik=1,nbindex
939	!
940	j = ((kindex(ik)-1)/iim) + 1
941	i = (kindex(ik) - (j-1)*iim)
942	!-
943	!- Store variable to help describe the grid
944	!- once the points are gathered.
945	!-
946	fneighbours(i,j,:) = neighbours(ik,:)
947	!
948	fresolution(i,j,:) = resolution(ik,:)
949	ENDDO
950	ELSE
951	!-- Second, in case it is a watchout file
952	ALLOCATE (tmpdata(iim,jjm))
953	tmpdata(:,:) = 0.0
954	!--
955	IF ( .NOT. contfrac_exists ) THEN
956	CALL ipslerr_p (3,'forcing_read_interpol', &
957	& 'Could get contfrac variable in a watchout file :',TRIM(filename), &
958	& '(Problem with file ?)')
959	ENDIF
960	CALL flinget_buffer (force_id,'contfrac',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
961	CALL forcing_zoom(data_full, fcontfrac)
962	!
963	! We need to make sure that when we gather the points we pick all
964	! the points where contfrac is above 0. Thus we prepare tair for
965	! subroutine forcing_landind
966	!
967	DO i=1,iim
968	DO j=1,jjm
969	IF ( j==1 .AND. i<ii_begin) fcontfrac(i,j)=0.
970	IF ( j==jjm .AND. i>ii_end) fcontfrac(i,j)=0.
971	IF ( fcontfrac(i,j) <= EPSILON(1.) ) THEN
972	tair(i,j) = 999999.
973	ENDIF
974	ENDDO
975	ENDDO
976
977	!---
978	!--- Create the index table
979	!---
980	!--- This job return a LOCAL kindex
981	CALL forcing_landind(iim, jjm, tair, nbindex, kindex, i_test, j_test)
982	#ifdef CPP_PARA
983	! We keep previous function forcing_landind, only to get a valid (i_test,j_test)
984	! Force nbindex points for parallel computation
985	nbindex=nbp_loc
986	CALL scatter(index_g,kindex)
987	kindex(:)=kindex(:)-offset
988	! kindex(:)=kindex(:)-(jj_begin-1)*iim_g
989	#endif
990	ik=MAX(nbindex/2,1)
991	j_test = (((kindex(ik)-1)/iim) + 1)
992	i_test = (kindex(ik) - (j_test-1)*iim)
993	IF (printlev_loc >= 5) THEN
994	WRITE(numout,*) 'New test point is : ', i_test, j_test
995	ENDIF
996	!---
997	!--- Allocate grid stuff
998	!---
999	CALL grid_init ( nbindex, 4, regular_lonlat, "ForcingGrid" )
1000	neighbours(:,:) = -1
1001	resolution(:,:) = 0.
1002	min_resol(:) = 1.e6
1003	max_resol(:) = -1.
1004	!---
1005	!--- All grid variables
1006	!---
1007	!-
1008	!- Get variables to help describe the grid
1009	CALL flinquery_var(force_id, 'neighboursNN', neighbours_exists)
1010	IF ( .NOT. neighbours_exists ) THEN
1011	CALL ipslerr_p (3,'forcing_read_interpol', &
1012	& 'Could get neighbours in a watchout file :',TRIM(filename), &
1013	& '(Problem with file ?)')
1014	ELSE
1015	IF (printlev_loc >= 2) WRITE(numout,*) "Watchout file contains neighbours and resolutions."
1016	ENDIF
1017	!
1018	fneighbours(:,:,:) = undef_int
1019	!
1020	!- once the points are gathered.
1021	CALL flinget_buffer (force_id,'neighboursNN',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1022	CALL forcing_zoom(data_full, tmpdata)
1023	WHERE(tmpdata(:,:) < undef_int)
1024	fneighbours(:,:,1) = NINT(tmpdata(:,:))
1025	ENDWHERE
1026	!
1027	CALL flinget_buffer (force_id,'neighboursNE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1028	CALL forcing_zoom(data_full, tmpdata)
1029	WHERE(tmpdata(:,:) < undef_int)
1030	fneighbours(:,:,2) = NINT(tmpdata(:,:))
1031	ENDWHERE
1032	!
1033	CALL flinget_buffer (force_id,'neighboursEE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1034	CALL forcing_zoom(data_full, tmpdata)
1035	WHERE(tmpdata(:,:) < undef_int)
1036	fneighbours(:,:,3) = NINT(tmpdata(:,:))
1037	ENDWHERE
1038	!
1039	CALL flinget_buffer (force_id,'neighboursSE',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1040	CALL forcing_zoom(data_full, tmpdata)
1041	WHERE(tmpdata(:,:) < undef_int)
1042	fneighbours(:,:,4) = NINT(tmpdata(:,:))
1043	ENDWHERE
1044	!
1045	CALL flinget_buffer (force_id,'neighboursSS',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1046	CALL forcing_zoom(data_full, tmpdata)
1047	WHERE(tmpdata(:,:) < undef_int)
1048	fneighbours(:,:,5) = NINT(tmpdata(:,:))
1049	ENDWHERE
1050	!
1051	CALL flinget_buffer (force_id,'neighboursSW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1052	CALL forcing_zoom(data_full, tmpdata)
1053	WHERE(tmpdata(:,:) < undef_int)
1054	fneighbours(:,:,6) = NINT(tmpdata(:,:))
1055	ENDWHERE
1056	!
1057	CALL flinget_buffer (force_id,'neighboursWW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1058	CALL forcing_zoom(data_full, tmpdata)
1059	WHERE(tmpdata(:,:) < undef_int)
1060	fneighbours(:,:,7) = NINT(tmpdata(:,:))
1061	ENDWHERE
1062	!
1063	CALL flinget_buffer (force_id,'neighboursNW',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1064	CALL forcing_zoom(data_full, tmpdata)
1065	WHERE(tmpdata(:,:) < undef_int)
1066	fneighbours(:,:,8) = NINT(tmpdata(:,:))
1067	ENDWHERE
1068	!
1069	! now, resolution of the grid
1070	CALL flinget_buffer (force_id,'resolutionX',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1071	CALL forcing_zoom(data_full, tmpdata)
1072	fresolution(:,:,1) = tmpdata(:,:)
1073	!
1074	CALL flinget_buffer (force_id,'resolutionY',iim_full, jjm_full, llm_full, ttm, 1, 1, data_full)
1075	CALL forcing_zoom(data_full, tmpdata)
1076	fresolution(:,:,2) = tmpdata(:,:)
1077	!
1078	DO ik=1,nbindex
1079	!
1080	j = ((kindex(ik)-1)/iim) + 1
1081	i = (kindex(ik) - (j-1)*iim)
1082	!-
1083	!- Store variable to help describe the grid
1084	!- once the points are gathered.
1085	!-
1086	neighbours(ik,:) = fneighbours(i,j,:)
1087	!
1088	resolution(ik,:) = fresolution(i,j,:)
1089	!
1090
1091	ENDDO
1092	CALL gather(neighbours,neighbours_g)
1093	CALL gather(resolution,resolution_g)
1094	min_resol(1) = MINVAL(resolution(:,1))
1095	min_resol(2) = MAXVAL(resolution(:,2))
1096	max_resol(1) = MAXVAL(resolution(:,1))
1097	max_resol(2) = MAXVAL(resolution(:,2))
1098	!
1099	area(:) = resolution(:,1)*resolution(:,2)
1100	CALL gather(area,area_g)
1101	!--
1102	DEALLOCATE (tmpdata)
1103	ENDIF
1104	IF (printlev_loc >= 2) WRITE(numout,*) 'contfrac : ', MINVAL(fcontfrac), MAXVAL(fcontfrac)
1105	!---
1106	ENDIF
1107	!---
1108	IF (printlev_loc >= 5) THEN
1109	WRITE(numout,*) &
1110	& 'The dates : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1111	ENDIF
1112	!---
1113	!--- Here we do the work in case only interpolation is needed.
1114	!---
1115	IF ( initialized .AND. interpol ) THEN
1116	!---
1117	IF ( daily_interpol ) THEN
1118
1119	IF (split > 1) THEN
1120	IF ( itau_split <= (split/2.) ) THEN
1121	rw = REAL(itau_split+split/2.)/split
1122	ELSE
1123	rw = REAL(itau_split-split/2.)/split
1124	ENDIF
1125	ELSE
1126	rw = 1.
1127	ENDIF
1128
1129	IF ((last_read == 0) .OR. ( rw==(1./split)) ) THEN
1130	!---
1131	!----- Start or Restart
1132	IF (last_read == 0) THEN
1133	! Case of a restart or a shift in the forcing file.
1134	IF (itau_read > 1) THEN
1135	itau_read_nm1=itau_read-1
1136	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1137	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1138	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1139	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1140	& force_id, wind_N_exists)
1141	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1142	! Case of a simple start.
1143	ELSE
1144	itau_read_nm1 = un
1145	IF (printlev_loc >= 2) WRITE(numout,*) "we will use the forcing of the first day to initialize "
1146	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1147	& swdown_nm1, rainf_nm1, snowf_nm1, tmin_nm1, &
1148	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1149	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1150	& force_id, wind_N_exists)
1151	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_nm1, itau_read_nm1, tmax_nm1, force_id )
1152	ENDIF
1153	tmin_nm2(:,:)=tmin_nm1(:,:)
1154	IF ( dt_force .GT. 3600. ) THEN
1155	mean_coszang(:,:) = 0.0
1156	daylength_n(:,:) = 0.
1157	DO is=1,split
1158	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1159	julian_for = itau2date(itau_read-1, date0, dt_force)
1160	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1161	! first day of this year
1162	CALL ymds2ju (yy,1,1,0.0, julian0)
1163	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1164	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1165	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1166	WHERE( coszang(:,:) > 0. )
1167	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1168	ENDWHERE
1169	ENDDO
1170	mean_coszang(:,:) = mean_coszang(:,:)/split
1171	daylength_nm1(:,:)=daylength_n(:,:)
1172	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1173	ENDIF
1174	ELSE
1175	!----- Normal mode : copy old step
1176	swdown_nm1(:,:) = swdown_n(:,:)
1177	rainf_nm1(:,:) = rainf_n(:,:)
1178	snowf_nm1(:,:) = snowf_n(:,:)
1179	tair_nm1(:,:) = tair_n(:,:)
1180	u_nm1(:,:) = u_n(:,:)
1181	v_nm1(:,:) = v_n(:,:)
1182	qair_nm1(:,:) = qair_n(:,:)
1183	pb_nm1(:,:) = pb_n(:,:)
1184	lwdown_nm1(:,:) = lwdown_n(:,:)
1185	tmin_nm2(:,:) = tmin_nm1(:,:)
1186	tmin_nm1(:,:) = tmin_n(:,:)
1187	tmax_nm1(:,:) = tmax_n(:,:)
1188
1189	IF (is_watchout) THEN
1190	zlev_nm1(:,:) = zlev_n(:,:)
1191	zlevuv_nm1(:,:) = zlevuv_n(:,:)
1192	! Net surface short-wave flux
1193	SWnet_nm1(:,:) = SWnet_n(:,:)
1194	! Air potential energy
1195	Eair_nm1(:,:) = Eair_n(:,:)
1196	! Coeficients A from the PBL resolution for T
1197	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1198	! Coeficients A from the PBL resolution for q
1199	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1200	! Coeficients B from the PBL resolution for T
1201	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1202	! Coeficients B from the PBL resolution for q
1203	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1204	! Surface drag
1205	cdrag_nm1(:,:) = cdrag_n(:,:)
1206	! CO2 concentration in the canopy
1207	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1208	ENDIF
1209	itau_read_nm1 = itau_read_n
1210	ENDIF
1211	!-----
1212	!-----
1213	IF(last_read==0)THEN
1214	itau_read_n = itau_read
1215	ELSE
1216	itau_read_n = itau_read+1
1217	ENDIF
1218
1219	IF (itau_read_n > ttm) THEN
1220	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1221	WRITE(numout,*) &
1222	& 'WARNING : We are going back to the start of the file'
1223	itau_read_n =1
1224	ENDIF
1225	IF (printlev_loc >= 5) THEN
1226	WRITE(numout,*) &
1227	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1228	ENDIF
1229	!-----
1230	!----- Get a reduced julian day !
1231	!----- This is needed because we lack the precision on 32 bit machines.
1232	!-----
1233	IF ( dt_force .GT. 3600. ) THEN
1234	julian_for = itau2date(itau_read-1, date0, dt_force)
1235	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1236
1237	! first day of this year
1238	CALL ymds2ju (yy,1,1,0.0, julian0)
1239	!-----
1240	IF (printlev_loc >= 5) THEN
1241	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1242	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1243	ENDIF
1244	ENDIF
1245	!-----
1246	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1247	& swdown_n, rainf_n, snowf_n, tmin_n, &
1248	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1249	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1250	& force_id, wind_N_exists)
1251	CALL forcing_just_read_tmax (iim, jjm, ttm, itau_read_n, itau_read_n, tmax_n, force_id )
1252
1253	!---
1254	last_read = itau_read_n
1255	!-----
1256	!----- Compute mean solar angle for the comming period
1257	!-----
1258	IF (printlev_loc >= 5) WRITE(numout,*) 'Going into solarang', split, one_day
1259	!-----
1260
1261	!-----
1262	ENDIF
1263	!---
1264	IF ( itau_split == 1. ) THEN
1265	IF ( dt_force .GT. 3600. ) THEN
1266	mean_coszang(:,:) = 0.0
1267	daylength_nm1(:,:)=daylength_n(:,:)
1268	daylength_n(:,:) = 0.
1269	DO is=1,split
1270	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1271	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1272	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1273	! first day of this year
1274	CALL ymds2ju (yy,1,1,0.0, julian0)
1275	CALL solarang (julian, julian0, iim, jjm, lon*0, lat, coszang)
1276	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1277	WHERE( coszang(:,:) > 0. )
1278	daylength_n(:,:)=daylength_n(:,:)+1./split*24
1279	ENDWHERE
1280	ENDDO
1281	mean_coszang(:,:) = mean_coszang(:,:)/split
1282	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1283	ENDIF
1284	ENDIF
1285
1286	!--- Do the interpolation
1287	IF (printlev_loc >= 5) WRITE(numout,*) 'Doing the interpolation between time steps'
1288	!---
1289
1290	IF (printlev_loc >= 5) WRITE(numout,*) 'Coeff of interpollation : ',rw
1291	!---
1292
1293	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1294	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1295	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1296
1297	!--- Take care of the height of the vertical levels
1298	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1299	zlevuv(:,:) = (zlevuv_n(:,:)-zlevuv_nm1(:,:))*rw + zlevuv_nm1(:,:)
1300
1301	hour=REAL(itau_split)/split*24
1302	startday_n(:,:)=12.-daylength_n(:,:)/2.
1303	startday_nm1(:,:)=12.-daylength_nm1(:,:)/2.
1304
1305	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1306	tair(:,:)=(tmax_nm1(:,:)-tmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1307	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_nm1(:,:)+tmin_nm1(:,:))/2.
1308	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1309	tair(:,:)=(tmax_n(:,:)-tmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1310	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (tmax_n(:,:)+tmin_n(:,:))/2.
1311	ELSEWHERE ( hour < startday_n(:,:) )
1312	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1313	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1314	ELSEWHERE
1315	tair(:,:)=(tmax_nm1(:,:)-tmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1316	& (24.-14.+startday_n(:,:))-14.))+(tmax_nm1(:,:)+tmin_n(:,:))/2.
1317	ENDWHERE
1318
1319	CALL weathgen_qsat_2d (iim,jjm,tmin_n,pb,qsattmin_n)
1320	CALL weathgen_qsat_2d (iim,jjm,tmin_nm1,pb,qsattmin_nm1)
1321	CALL weathgen_qsat_2d (iim,jjm,tair,pb,qsatta)
1322
1323	!---
1324	qmin_nm1(:,:) = MIN(qair_nm1(:,:),0.99*qsattmin_nm1(:,:))
1325	qmin_n(:,:) = MIN(qair_n(:,:),0.99*qsattmin_n(:,:))
1326	qmax_nm1(:,:) = (qair_nm1(:,:)-qmin_nm1(:,:)) + qair_nm1(:,:)
1327	qmax_n(:,:) = (qair_n(:,:)-qmin_n(:,:)) + qair_n(:,:)
1328
1329	qsa(:,:) = 0.99*qsatta(:,:)
1330
1331
1332	WHERE ( ( hour >= startday_n(:,:) ) .AND. ( hour > 12) .AND. ( hour <= 14) )
1333	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_nm1(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1334	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_nm1(:,:)+qmin_nm1(:,:))/2.)
1335	ELSEWHERE( ( hour >= startday_n(:,:) ) .AND. ( hour <= 12) )
1336	qair(:,:)=MIN(qsa(:,:),(qmax_n(:,:)-qmin_n(:,:))/2 * ( sin(pi/(14-startday_n(:,:))(hour-0.5 &
1337	& (14.-startday_n(:,:))-startday_n(:,:))) )+ (qmax_n(:,:)+qmin_n(:,:))/2.)
1338	ELSEWHERE ( hour < startday_n(:,:) )
1339	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_nm1(:,:) ) &
1340	& (hour + 24.+0.5*(24.-14.+startday_nm1(:,:) )-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1341	ELSEWHERE
1342	qair(:,:)=MIN(qsa(:,:),(qmax_nm1(:,:)-qmin_n(:,:))/2.sin(pi/(24.-14.+startday_n(:,:))(hour+0.5* &
1343	& (24.-14.+startday_n(:,:))-14.))+(qmax_nm1(:,:)+qmin_n(:,:))/2.)
1344	ENDWHERE
1345
1346	IF (is_watchout) THEN
1347	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1348	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1349	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1350	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1351	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1352	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1353	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1354	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1355	ENDIF
1356	!---
1357	!--- Here we need to allow for an option
1358	!--- where radiative energy is conserved
1359	!---
1360	IF ( lwdown_cons ) THEN
1361	lwdown(:,:) = lwdown_n(:,:)
1362	ELSE
1363	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1364	ENDIF
1365	!---
1366	!--- For the solar radiation we decompose the mean value
1367	!--- using the zenith angle of the sun, conservative approach under 2000W/m2
1368	!----
1369	IF (printlev_loc >= 5) WRITE(numout,*) 'Ready to deal with the solar radiation'
1370	!----
1371	! We compute mean SWdown between t and t+Dt then we take t+Dt/2.
1372	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1373	!!$ julian = julian_for + rw*dt_force/one_day
1374	IF (printlev_loc >= 5) THEN
1375	WRITE(numout,'(a,f20.10,2I3)') &
1376	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1377	ENDIF
1378
1379	CALL solarang(julian, julian0, iim, jjm, lon*0, lat, coszang)
1380
1381	WHERE ((mean_coszang(:,:) > 0.) .AND. (hour <= 12 ))
1382	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1383	ELSEWHERE ((mean_coszang(:,:) > 0.) .AND. (hour > 12 ))
1384	swdown(:,:) = swdown_nm1(:,:) *coszang(:,:)/mean_coszang(:,:)
1385	ELSEWHERE
1386	swdown(:,:) = 0.0
1387	END WHERE
1388
1389	WHERE (swdown(:,:) > 2000. )
1390	swdown(:,:) = 2000.
1391	END WHERE
1392
1393	IF (printlev_loc >= 5) THEN
1394	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1395	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1396	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1397	IF (is_watchout) THEN
1398	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1399	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1400	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1401	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1402	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1403	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1404	ENDIF
1405	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1406	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1407	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1408	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1409	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1410	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1411	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1412	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1413	ENDIF
1414	!---
1415	!--- For precip we suppose that the rain
1416	!--- is the sum over the next 6 hours
1417	!---
1418	WHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)>=273.15))
1419	rainf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1420	snowf(:,:) = 0.0
1421	ELSEWHERE ((itau_split <= nb_spread).AND.(hour<=12).AND.(tair(:,:)<273.15))
1422	snowf(:,:) = rainf_n(:,:) *(split/REAL(nb_spread))
1423	rainf(:,:) = 0.0
1424	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)>=273.15))
1425	rainf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1426	snowf(:,:) = 0.0
1427	ELSEWHERE ((itau_split <= nb_spread).AND.(hour>12).AND.(tair(:,:)<273.15))
1428	snowf(:,:) = rainf_nm1(:,:) *(split/REAL(nb_spread))
1429	rainf(:,:) = 0.0
1430	ELSEWHERE
1431	snowf(:,:) = 0.0
1432	rainf(:,:) = 0.0
1433	ENDWHERE
1434
1435	IF (printlev_loc >= 5) THEN
1436	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1437	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1438	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1439	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1440	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1441	ENDIF
1442	!---
1443
1444	ELSE ! If not daily_interpol
1445
1446	IF (itau_read /= last_read) THEN
1447	!---
1448	!----- Start or Restart
1449	IF (itau_read_n == 0) THEN
1450	! Case of a restart or a shift in the forcing file.
1451	IF (itau_read > 1) THEN
1452	itau_read_nm1=itau_read-1
1453	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1454	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1455	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1456	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1457	& force_id, wind_N_exists)
1458	! Case of a simple start.
1459	ELSE IF (dt_force*ttm > one_day-1. ) THEN
1460	! if the forcing file contains at least 24 hours,
1461	! we will use the last forcing step of the first day
1462	! as initiale condition to prevent first shift off reading.
1463	itau_read_nm1 = NINT (one_day/dt_force)
1464	IF (printlev_loc >= 1) WRITE(numout,) "The forcing file contains 24 hours :",dt_forcettm,one_day-1.
1465	IF (printlev_loc >= 1) WRITE(numout,*) "We will use the last forcing step of the first day : itau_read_nm1 ",&
1466	itau_read_nm1
1467	CALL forcing_just_read (iim, jjm, zlev_nm1, zlevuv_nm1, ttm, itau_read_nm1, itau_read_nm1, &
1468	& swdown_nm1, rainf_nm1, snowf_nm1, tair_nm1, &
1469	& u_nm1, v_nm1, qair_nm1, pb_nm1, lwdown_nm1, &
1470	& SWnet_nm1, Eair_nm1, petAcoef_nm1, peqAcoef_nm1, petBcoef_nm1, peqBcoef_nm1, cdrag_nm1, ccanopy_nm1, &
1471	& force_id, wind_N_exists)
1472	ELSE
1473	! if the forcing file contains less than 24 hours,
1474	! just say error !
1475	CALL ipslerr_p(3,'forcing_read_interpol', &
1476	& 'The forcing file contains less than 24 hours !', &
1477	& 'We can''t intialize interpolation with such a file.','')
1478	ENDIF
1479	ELSE
1480	!----- Normal mode : copy old step
1481	swdown_nm1(:,:) = swdown_n(:,:)
1482	rainf_nm1(:,:) = rainf_n(:,:)
1483	snowf_nm1(:,:) = snowf_n(:,:)
1484	tair_nm1(:,:) = tair_n(:,:)
1485	u_nm1(:,:) = u_n(:,:)
1486	v_nm1(:,:) = v_n(:,:)
1487	qair_nm1(:,:) = qair_n(:,:)
1488	pb_nm1(:,:) = pb_n(:,:)
1489	lwdown_nm1(:,:) = lwdown_n(:,:)
1490	IF (is_watchout) THEN
1491	zlev_nm1(:,:) = zlev_n(:,:)
1492	! Net surface short-wave flux
1493	SWnet_nm1(:,:) = SWnet_n(:,:)
1494	! Air potential energy
1495	Eair_nm1(:,:) = Eair_n(:,:)
1496	! Coeficients A from the PBL resolution for T
1497	petAcoef_nm1(:,:) = petAcoef_n(:,:)
1498	! Coeficients A from the PBL resolution for q
1499	peqAcoef_nm1(:,:) = peqAcoef_n(:,:)
1500	! Coeficients B from the PBL resolution for T
1501	petBcoef_nm1(:,:) = petBcoef_n(:,:)
1502	! Coeficients B from the PBL resolution for q
1503	peqBcoef_nm1(:,:) = peqBcoef_n(:,:)
1504	! Surface drag
1505	cdrag_nm1(:,:) = cdrag_n(:,:)
1506	! CO2 concentration in the canopy
1507	ccanopy_nm1(:,:) = ccanopy_n(:,:)
1508	ENDIF
1509	itau_read_nm1 = itau_read_n
1510	ENDIF
1511	!-----
1512	itau_read_n = itau_read
1513	IF (itau_read_n > ttm) THEN
1514	WRITE(numout,*) 'WARNING --WARNING --WARNING --WARNING '
1515	WRITE(numout,*) &
1516	& 'WARNING : We are going back to the start of the file'
1517	itau_read_n =1
1518	ENDIF
1519	IF (printlev_loc >= 5) THEN
1520	WRITE(numout,*) &
1521	& 'The dates 2 : ',itau_read,itau_split,itau_read_nm1,itau_read_n
1522	ENDIF
1523	!-----
1524	!----- Get a reduced julian day !
1525	!----- This is needed because we lack the precision on 32 bit machines.
1526	!-----
1527	IF ( dt_force .GT. 3600. ) THEN
1528	julian_for = itau2date(itau_read-1, date0, dt_force)
1529	CALL ju2ymds (julian_for, yy, mm, dd, ss)
1530
1531	! first day of this year
1532	CALL ymds2ju (yy,1,1,0.0, julian0)
1533	!-----
1534	IF (printlev_loc >= 5) THEN
1535	WRITE(numout,*) 'Forcing for Julian day ',julian_for,'is read'
1536	WRITE(numout,*) 'Date for this day ',yy,' / ',mm,' / ',dd," ",ss
1537	ENDIF
1538	ENDIF
1539	!-----
1540	CALL forcing_just_read (iim, jjm, zlev_n, zlevuv_n, ttm, itau_read_n, itau_read_n, &
1541	& swdown_n, rainf_n, snowf_n, tair_n, &
1542	& u_n, v_n, qair_n, pb_n, lwdown_n, &
1543	& SWnet_n, Eair_n, petAcoef_n, peqAcoef_n, petBcoef_n, peqBcoef_n, cdrag_n, ccanopy_n, &
1544	& force_id, wind_N_exists)
1545	!---
1546	last_read = itau_read_n
1547	!-----
1548	!----- Compute mean solar angle for the comming period
1549	!-----
1550	IF (printlev_loc >= 5) WRITE(numout,*) 'Going into solarang', split, one_day
1551	!-----
1552	IF ( dt_force .GT. 3600. ) THEN
1553	mean_coszang(:,:) = 0.0
1554	DO is=1,split
1555	!MM we compute mean SWdown between t and t+Dt then I take t+Dt/2.
1556	julian = julian_for+((is-0.5)/split)*dt_force/one_day
1557	! first day of this year
1558	CALL ymds2ju (yy,1,1,0.0, julian0)
1559	!!$ julian = julian_for+(FLOAT(is)/split)*dt_force/one_day
1560	CALL solarang (julian, julian0, iim, jjm, lon, lat, coszang)
1561	mean_coszang(:,:) = mean_coszang(:,:)+coszang(:,:)
1562	ENDDO
1563	mean_coszang(:,:) = mean_coszang(:,:)/split
1564	! WRITE(,) "mean_coszang =",MAXVAL(mean_coszang)
1565	ENDIF
1566	!-----
1567	ENDIF
1568	!---
1569	!--- Do the interpolation
1570	IF (printlev_loc >= 5) WRITE(numout,*) 'Doing the interpolation between time steps'
1571	!---
1572	IF (split > 1) THEN
1573	rw = REAL(itau_split)/split
1574	ELSE
1575	rw = 1.
1576	ENDIF
1577	IF (printlev_loc >= 5) WRITE(numout,*) 'Coeff of interpollation : ',rw
1578	!---
1579	qair(:,:) = (qair_n(:,:)-qair_nm1(:,:))*rw + qair_nm1(:,:)
1580	tair(:,:) = (tair_n(:,:)-tair_nm1(:,:))*rw + tair_nm1(:,:)
1581	pb(:,:) = (pb_n(:,:)-pb_nm1(:,:))*rw + pb_nm1(:,:)
1582	u(:,:) = (u_n(:,:)-u_nm1(:,:))*rw + u_nm1(:,:)
1583	v(:,:) = (v_n(:,:)-v_nm1(:,:))*rw + v_nm1(:,:)
1584	IF (is_watchout) THEN
1585	zlev(:,:) = (zlev_n(:,:)-zlev_nm1(:,:))*rw + zlev_nm1(:,:)
1586	zlevuv(:,:) = zlev(:,:)
1587	SWnet(:,:) = (SWnet_n(:,:)-SWnet_nm1(:,:))*rw + SWnet_nm1(:,:)
1588	Eair(:,:) = (Eair_n(:,:)-Eair_nm1(:,:))*rw + Eair_nm1(:,:)
1589	petAcoef(:,:) = (petAcoef_n(:,:)-petAcoef_nm1(:,:))*rw + petAcoef_nm1(:,:)
1590	peqAcoef(:,:) = (peqAcoef_n(:,:)-peqAcoef_nm1(:,:))*rw + peqAcoef_nm1(:,:)
1591	petBcoef(:,:) = (petBcoef_n(:,:)-petBcoef_nm1(:,:))*rw + petBcoef_nm1(:,:)
1592	peqBcoef(:,:) = (peqBcoef_n(:,:)-peqBcoef_nm1(:,:))*rw + peqBcoef_nm1(:,:)
1593	cdrag(:,:) = (cdrag_n(:,:)-cdrag_nm1(:,:))*rw + cdrag_nm1(:,:)
1594	ccanopy(:,:) = (ccanopy_n(:,:)-ccanopy_nm1(:,:))*rw + ccanopy_nm1(:,:)
1595	ENDIF
1596	!---
1597	!--- Here we need to allow for an option
1598	!--- where radiative energy is conserved
1599	!---
1600	IF ( lwdown_cons ) THEN
1601	lwdown(:,:) = lwdown_n(:,:)
1602	ELSE
1603	lwdown(:,:) = (lwdown_n(:,:)-lwdown_nm1(:,:))*rw + lwdown_nm1(:,:)
1604	ENDIF
1605	!---
1606	!--- For the solar radiation we decompose the mean value
1607	!--- using the zenith angle of the sun if the time step in the forcing data is
1608	!---- more than an hour. Else we use the standard linera interpolation
1609	!----
1610	IF (printlev_loc >= 5) WRITE(numout,*) 'Ready to deal with the solar radiation'
1611	!----
1612	IF ( dt_force .GT. 3600. ) THEN
1613
1614	! In this case solar radiation will be conserved
1615
1616	! We compute mean SWdown between t and t+Dt then we take t+Dt/2.
1617	julian = julian_for + (itau_split-0.5)/split*dt_force/one_day
1618	!!$ julian = julian_for + rw*dt_force/one_day
1619	IF (printlev_loc >= 5) THEN
1620	WRITE(numout,'(a,f20.10,2I3)') &
1621	& 'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
1622	ENDIF
1623	!---
1624	CALL solarang(julian, julian0, iim, jjm, lon, lat, coszang)
1625	!---
1626	WHERE (mean_coszang(:,:) > 0.)
1627	swdown(:,:) = swdown_n(:,:) *coszang(:,:)/mean_coszang(:,:)
1628	ELSEWHERE
1629	swdown(:,:) = 0.0
1630	END WHERE
1631	!---
1632	WHERE (swdown(:,:) > 2000. )
1633	swdown(:,:) = 2000.
1634	END WHERE
1635	!---
1636	ELSE ! If dt_force < 3600 (1h)
1637
1638	IF ( swdown_cons ) THEN
1639	! Conserve swdown radiation
1640	swdown(:,:) = swdown_n(:,:)
1641	ELSE
1642	swdown(:,:) = (swdown_n(:,:)-swdown_nm1(:,:))*rw + swdown_nm1(:,:)
1643	ENDIF
1644	!---
1645	ENDIF
1646	!---
1647	IF (printlev_loc >= 5) THEN
1648	WRITE(numout,*) '__ Forcing read at ',itau_split,' :',i_test, j_test
1649	WRITE(numout,*) 'SWdown : ',swdown_nm1(i_test, j_test), &
1650	& ' < ', swdown(i_test, j_test), ' < ', swdown_n(i_test, j_test)
1651	IF (is_watchout) THEN
1652	WRITE(numout,*) 'SWnet : ',swnet_nm1(i_test, j_test), &
1653	& ' < ', swnet(i_test, j_test), ' < ', swnet_n(i_test, j_test)
1654	WRITE(numout,*) 'levels :',zlev_nm1(i_test, j_test), &
1655	& ' < ', zlev(i_test, j_test), ' < ', zlev_n(i_test, j_test)
1656	WRITE(numout,*) 'EAIR :',Eair_nm1(i_test, j_test), &
1657	& ' < ', eair(i_test, j_test), ' < ', Eair_n(i_test, j_test)
1658	ENDIF
1659	WRITE(numout,*) 'TAIR :',tair_nm1(i_test, j_test), &
1660	& ' < ', tair(i_test, j_test), ' < ', tair_n(i_test, j_test)
1661	WRITE(numout,*) 'QAIR :',qair_nm1(i_test, j_test), &
1662	& ' < ', qair(i_test, j_test), ' < ', qair_n(i_test, j_test)
1663	WRITE(numout,*) 'U :',u_nm1(i_test, j_test), &
1664	& ' < ', u(i_test, j_test), ' < ', u_n(i_test, j_test)
1665	WRITE(numout,*) 'V :',v_nm1(i_test, j_test), &
1666	& ' < ', v(i_test, j_test), ' < ', v_n(i_test, j_test)
1667	ENDIF
1668	!---
1669	!--- For precip we suppose that the rain
1670	!--- is the sum over the next 6 hours
1671	!---
1672	IF (itau_split <= nb_spread) THEN
1673	rainf(:,:) = rainf_n(:,:)*(split/REAL(nb_spread))
1674	snowf(:,:) = snowf_n(:,:)*(split/REAL(nb_spread))
1675	ELSE
1676	rainf(:,:) = 0.0
1677	snowf(:,:) = 0.0
1678	ENDIF
1679	IF (printlev_loc >= 5) THEN
1680	WRITE(numout,*) '__ Forcing read at ',itau_split,' :'
1681	WRITE(numout,*) 'Rainf : ',rainf_nm1(i_test, j_test), &
1682	& ' < ', rainf(i_test, j_test), ' < ', rainf_n(i_test, j_test)
1683	WRITE(numout,*) 'Snowf : ',snowf_nm1(i_test, j_test), &
1684	& ' < ', snowf(i_test, j_test), ' < ', snowf_n(i_test, j_test)
1685	ENDIF
1686	!---
1687	ENDIF ! (daily_interpol)
1688	ENDIF
1689	!---
1690	!--- Here we might put the call to the weather generator ... one day.
1691	!--- Pour le moment, le branchement entre interpolation et generateur de temps
1692	!--- est fait au-dessus.
1693	!---
1694	!- IF ( initialized .AND. weathergen ) THEN
1695	!- ....
1696	!- ENDIF
1697	!---
1698	!--- At this point the code should be initialized. If not we have a problem !
1699	!---
1700	IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
1701	!---
1702	initialized = .TRUE.
1703	!---
1704	ELSE
1705	IF ( .NOT. initialized ) THEN
1706	WRITE(numout,*) 'Why is the code forcing_read not initialized ?'
1707	WRITE(numout,*) 'Have you called it with both time-steps set to zero ?'
1708	CALL ipslerr_p(3,'forcing_read_interpol','Pb in initialization','','')
1709	ENDIF
1710	ENDIF
1711
1712	END SUBROUTINE forcing_read_interpol
1713
1714
1715	!! ==============================================================================================================================\n
1716	!! SUBROUTINE : forcing_just_read
1717	!!
1718	!>\BRIEF
1719	!!
1720	!!\n DESCRIPTION :
1721	!!
1722	!! RECENT CHANGE(S): None
1723	!!
1724	!! MAIN OUTPUT VARIABLE(S):
1725	!!
1726	!! REFERENCE(S) :
1727	!!
1728	!_ ================================================================================================================================
1729	SUBROUTINE forcing_just_read &
1730	& (iim, jjm, zlev, zlev_uv, ttm, itb, ite, &
1731	& swdown, rainf, snowf, tair, &
1732	& u, v, qair, pb, lwdown, &
1733	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy, &
1734	& force_id, wind_N_exists)
1735	!---------------------------------------------------------------------
1736	!- iim : Size of the grid in x
1737	!- jjm : size of the grid in y
1738	!- zlev : height of the varibales T and Q
1739	!- zlev_uv : height of the varibales U and V
1740	!- ttm : number of time steps in all in the forcing file
1741	!- itb, ite : index of respectively begin and end of read for each variable
1742	!- swdown : Downward solar radiation (W/m^2)
1743	!- rainf : Rainfall (kg/m^2s)
1744	!- snowf : Snowfall (kg/m^2s)
1745	!- tair : 2m air temperature (K)
1746	!- u and v : 2m (in theory !) wind speed (m/s)
1747	!- qair : 2m humidity (kg/kg)
1748	!- pb : Surface pressure (Pa)
1749	!- lwdown : Downward long wave radiation (W/m^2)
1750	!-
1751	!- From a WATCHOUT file :
1752	!- SWnet : Net surface short-wave flux
1753	!- Eair : Air potential energy
1754	!- petAcoef : Coeficients A from the PBL resolution for T
1755	!- peqAcoef : Coeficients A from the PBL resolution for q
1756	!- petBcoef : Coeficients B from the PBL resolution for T
1757	!- peqBcoef : Coeficients B from the PBL resolution for q
1758	!- cdrag : Surface drag
1759	!- ccanopy : CO2 concentration in the canopy
1760	!- force_id : FLINCOM file id.
1761	!- It is used to close the file at the end of the run.
1762	!- wind_N_exists : if Wind_N and Wind_E are in the file (and not just Wind)
1763	!---------------------------------------------------------------------
1764	IMPLICIT NONE
1765	!-
1766	INTEGER, INTENT(in) :: iim, jjm, ttm
1767	INTEGER, INTENT(in) :: itb, ite
1768	REAL, DIMENSION(iim,jjm), INTENT(out) :: zlev, zlev_uv, &
1769	& swdown, rainf, snowf, tair, u, v, qair, pb, lwdown
1770	! for watchout files
1771	REAL, DIMENSION(iim,jjm), INTENT(out) :: &
1772	& SWnet, Eair, petAcoef, peqAcoef, petBcoef, peqBcoef, cdrag, ccanopy
1773	INTEGER, INTENT(in) :: force_id
1774	! if Wind_N and Wind_E are in the file (and not just Wind)
1775	LOGICAL, INTENT(in) :: wind_N_exists
1776	INTEGER :: i, j
1777	REAL :: rau
1778
1779	!-
1780	!---------------------------------------------------------------------
1781	IF ( daily_interpol ) THEN
1782	CALL flinget_buffer (force_id,'Tmin' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1783	CALL forcing_zoom(data_full, tair)
1784	CALL flinget_buffer (force_id,'precip' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1785	CALL forcing_zoom(data_full, rainf)
1786	ELSE
1787	CALL flinget_buffer (force_id,'Tair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1788	CALL forcing_zoom(data_full, tair)
1789	CALL flinget_buffer (force_id,'Snowf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1790	CALL forcing_zoom(data_full, snowf)
1791	CALL flinget_buffer (force_id,'Rainf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1792	CALL forcing_zoom(data_full, rainf)
1793	ENDIF
1794
1795
1796	CALL flinget_buffer (force_id,'SWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1797	CALL forcing_zoom(data_full, swdown)
1798	CALL flinget_buffer (force_id,'LWdown', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1799	CALL forcing_zoom(data_full, lwdown)
1800
1801	CALL flinget_buffer (force_id,'PSurf' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1802	CALL forcing_zoom(data_full, pb)
1803	CALL flinget_buffer (force_id,'Qair' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1804	CALL forcing_zoom(data_full, qair)
1805	!---
1806	IF ( wind_N_exists ) THEN
1807	CALL flinget_buffer (force_id,'Wind_N', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1808	CALL forcing_zoom(data_full, u)
1809	CALL flinget_buffer (force_id,'Wind_E', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1810	CALL forcing_zoom(data_full, v)
1811	ELSE
1812	CALL flinget_buffer (force_id,'Wind', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1813	CALL forcing_zoom(data_full, u)
1814	v=0.0
1815	ENDIF
1816
1817	!-
1818	!- Deal with the height of the atmospheric forcing varibles
1819	!-
1820	!----
1821	IF ( zheight ) THEN
1822	zlev(:,:) = zlev_fixed
1823	ELSE IF ( zsigma .OR. zhybrid ) THEN
1824	DO i=1,iim
1825	DO j=1,jjm
1826	IF ( tair(i,j) < val_exp ) THEN
1827	rau = pb(i,j)/(cte_molr*tair(i,j))
1828
1829	zlev(i,j) = (pb(i,j) - (zhybrid_a + zhybrid_bpb(i,j)))/(rau cte_grav)
1830	ELSE
1831	zlev(i,j) = 0.0
1832	ENDIF
1833	ENDDO
1834	ENDDO
1835	ELSE IF ( zlevels ) THEN
1836	CALL flinget_buffer (force_id,'Levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1837	CALL forcing_zoom(data_full, zlev)
1838	ELSE
1839	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1840	& 'We cannot determine the height for T and Q.','stop readdim2')
1841	ENDIF
1842
1843	IF ( zsamelev_uv ) THEN
1844	zlev_uv(:,:) = zlev(:,:)
1845	ELSE
1846	IF ( zheight ) THEN
1847	zlev_uv(:,:) = zlevuv_fixed
1848	ELSE IF ( zsigma .OR. zhybrid ) THEN
1849	DO i=1,iim
1850	DO j=1,jjm
1851	IF ( tair(i,j) < val_exp ) THEN
1852	rau = pb(i,j)/(cte_molr*tair(i,j))
1853
1854	zlev_uv(i,j) = (pb(i,j) - (zhybriduv_a + zhybriduv_bpb(i,j)))/(rau cte_grav)
1855	ELSE
1856	zlev_uv(i,j) = 0.0
1857	ENDIF
1858	ENDDO
1859	ENDDO
1860	ELSE IF ( zlevels ) THEN
1861	CALL flinget_buffer (force_id,'Levels_uv', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1862	CALL forcing_zoom(data_full, zlev_uv)
1863	ELSE
1864	CALL ipslerr(3, 'forcing_just_read','No case for the vertical levels was specified.', &
1865	& 'We cannot determine the height for U and V.','stop readdim2')
1866	ENDIF
1867	ENDIF
1868	!----
1869	IF ( is_watchout ) THEN
1870	CALL flinget_buffer (force_id,'levels', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1871	CALL forcing_zoom(data_full, zlev)
1872	!
1873	! If we are in WATHCOUT it means T,Q are at the same height as U,V
1874	!
1875	zlev_uv(:,:) = zlev(:,:)
1876	! Net surface short-wave flux
1877	CALL flinget_buffer (force_id,'SWnet', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1878	CALL forcing_zoom(data_full, SWnet)
1879	! Air potential energy
1880	CALL flinget_buffer (force_id,'Eair', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1881	CALL forcing_zoom(data_full, Eair)
1882	! Coeficients A from the PBL resolution for T
1883	CALL flinget_buffer (force_id,'petAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1884	CALL forcing_zoom(data_full, petAcoef)
1885	! Coeficients A from the PBL resolution for q
1886	CALL flinget_buffer (force_id,'peqAcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1887	CALL forcing_zoom(data_full, peqAcoef)
1888	! Coeficients B from the PBL resolution for T
1889	CALL flinget_buffer (force_id,'petBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1890	CALL forcing_zoom(data_full, petBcoef)
1891	! Coeficients B from the PBL resolution for q
1892	CALL flinget_buffer (force_id,'peqBcoef', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1893	CALL forcing_zoom(data_full, peqBcoef)
1894	! Surface drag
1895	CALL flinget_buffer (force_id,'cdrag', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1896	CALL forcing_zoom(data_full, cdrag)
1897	! CO2 concentration in the canopy
1898	CALL flinget_buffer (force_id,'ccanopy', iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1899	CALL forcing_zoom(data_full, ccanopy)
1900	ENDIF
1901	!
1902	!----
1903	IF (printlev_loc >= 5) WRITE(numout,*) 'Variables have been extracted between ',itb,&
1904	' and ',ite,' iterations of the forcing file.'
1905	!-------------------------
1906	END SUBROUTINE forcing_just_read
1907
1908
1909	!! ==============================================================================================================================\n
1910	!! SUBROUTINE : forcing_just_read_tmax
1911	!!
1912	!>\BRIEF
1913	!!
1914	!!\n DESCRIPTION :
1915	!!
1916	!! RECENT CHANGE(S): None
1917	!!
1918	!! MAIN OUTPUT VARIABLE(S):
1919	!!
1920	!! REFERENCE(S) :
1921	!!
1922	!_ ================================================================================================================================
1923	SUBROUTINE forcing_just_read_tmax &
1924	& (iim, jjm, ttm, itb, ite, tmax, force_id )
1925	!---------------------------------------------------------------------
1926	!- iim : Size of the grid in x
1927	!- jjm : size of the grid in y
1928	!- ttm : number of time steps in all in the forcing file
1929	!- itb, ite : index of respectively begin and end of read for each variable
1930	!- tmax : 2m air temperature (K)
1931	!- force_id : FLINCOM file id.
1932	!- It is used to close the file at the end of the run.
1933	!---------------------------------------------------------------------
1934	IMPLICIT NONE
1935	!-
1936	INTEGER, INTENT(in) :: iim, jjm, ttm
1937	INTEGER, INTENT(in) :: itb, ite
1938	REAL, DIMENSION(iim,jjm), INTENT(out) :: tmax
1939	INTEGER, INTENT(in) :: force_id
1940	!-
1941	!---------------------------------------------------------------------
1942	CALL flinget_buffer (force_id,'Tmax' , iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
1943	CALL forcing_zoom(data_full, tmax)
1944	!-------------------------
1945	END SUBROUTINE forcing_just_read_tmax
1946
1947
1948	!! ==============================================================================================================================\n
1949	!! SUBROUTINE : forcing_landind
1950	!!
1951	!>\BRIEF
1952	!!
1953	!!\n DESCRIPTION : This subroutine finds the indices of the land points over which the land
1954	!! surface scheme is going to run.
1955	!!
1956	!! RECENT CHANGE(S): None
1957	!!
1958	!! MAIN OUTPUT VARIABLE(S):
1959	!!
1960	!! REFERENCE(S) :
1961	!!
1962	!_ ================================================================================================================================
1963	SUBROUTINE forcing_landind(iim, jjm, tair, nbindex, kindex, i_test, j_test)
1964	!---
1965	!---
1966	!---
1967	IMPLICIT NONE
1968	!-
1969	!- ARGUMENTS
1970	!-
1971	INTEGER, INTENT(IN) :: iim, jjm
1972	REAL, INTENT(IN) :: tair(iim,jjm)
1973	INTEGER, INTENT(OUT) :: i_test, j_test, nbindex
1974	INTEGER, INTENT(OUT) :: kindex(iim*jjm)
1975	!-
1976	!- LOCAL
1977	INTEGER :: i, j, ig
1978	!-
1979	!-
1980	ig = 0
1981	i_test = 0
1982	j_test = 0
1983	!---
1984	IF (MINVAL(tair(:,:)) < 100.) THEN
1985	!----- In this case the 2m temperature is in Celsius
1986	DO j=1,jjm
1987	DO i=1,iim
1988	IF (tair(i,j) < 100.) THEN
1989	ig = ig+1
1990	kindex(ig) = (j-1)*iim+i
1991	!
1992	! Here we find at random a land-point on which we can do
1993	! some printouts for test.
1994	!
1995	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
1996	i_test = i
1997	j_test = j
1998	IF (printlev_loc >= 5) THEN
1999	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
2000	ENDIF
2001	ENDIF
2002	ENDIF
2003	ENDDO
2004	ENDDO
2005	ELSE
2006	!----- 2m temperature is in Kelvin
2007	DO j=1,jjm
2008	DO i=1,iim
2009	IF (tair(i,j) < 500.) THEN
2010	ig = ig+1
2011	kindex(ig) = (j-1)*iim+i
2012	!
2013	! Here we find at random a land-point on which we can do
2014	! some printouts for test.
2015	!
2016	IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
2017	i_test = i
2018	j_test = j
2019	IF (printlev_loc >= 5) THEN
2020	WRITE(numout,*) 'The test point chosen for output is : ', i_test, j_test
2021	ENDIF
2022	ENDIF
2023	ENDIF
2024	ENDDO
2025	ENDDO
2026	ENDIF
2027
2028	nbindex = ig
2029
2030	END SUBROUTINE forcing_landind
2031
2032
2033	!! ==============================================================================================================================\n
2034	!! SUBROUTINE : forcing_grid
2035	!!
2036	!>\BRIEF
2037	!!
2038	!!\n DESCRIPTION : This subroutine calculates the longitudes and latitudes of the model grid.
2039	!!
2040	!! RECENT CHANGE(S): None
2041	!!
2042	!! MAIN OUTPUT VARIABLE(S):
2043	!!
2044	!! REFERENCE(S) :
2045	!!
2046	!_ ================================================================================================================================
2047	SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,init_f)
2048
2049	IMPLICIT NONE
2050
2051	INTEGER, INTENT(in) :: iim, jjm, llm
2052	LOGICAL, INTENT(in) :: init_f
2053	REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat
2054	!-
2055	INTEGER :: i,j
2056	!-
2057	!- Should be unified one day
2058	!-
2059	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : options : ', weathergen, interpol
2060	!-
2061	IF ( weathergen ) THEN
2062	IF (init_f) THEN
2063	DO i = 1, iim
2064	lon(i,:) = limit_west + merid_res/2. + &
2065	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim)
2066	ENDDO
2067	!-
2068	DO j = 1, jjm
2069	lat(:,j) = limit_north - zonal_res/2. - &
2070	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm)
2071	ENDDO
2072	ELSE
2073	IF (is_root_prc) THEN
2074	DO i = 1, iim_g
2075	lon_g(i,:) = limit_west + merid_res/2. + &
2076	FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim_g)
2077	ENDDO
2078	!-
2079	DO j = 1, jjm_g
2080	lat_g(:,j) = limit_north - zonal_res/2. - &
2081	FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm_g)
2082	ENDDO
2083	ENDIF
2084	CALL bcast(lon_g)
2085	CALL bcast(lat_g)
2086	lon=lon_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2087	lat=lat_g(:,jj_para_begin(mpi_rank):jj_para_end(mpi_rank))
2088	ENDIF
2089	!-
2090	ELSEIF ( interpol ) THEN
2091	!-
2092	CALL forcing_zoom(lon_full, lon)
2093	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lon'
2094	CALL forcing_zoom(lat_full, lat)
2095	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : out of zoom on lat'
2096
2097	ELSE
2098	CALL ipslerr_p(3,'forcing_grid','Neither interpolation nor weather generator is specified.','','')
2099	ENDIF
2100
2101	IF ( printlev_loc>=3 ) WRITE(numout,*) 'forcing_grid : ended'
2102
2103	END SUBROUTINE forcing_grid
2104
2105
2106	!! ==============================================================================================================================\n
2107	!! SUBROUTINE : forcing_zoom
2108	!!
2109	!>\BRIEF This subroutine extracts the region of data over which we wish to run the model.
2110	!!
2111	!!\n DESCRIPTION : This subroutine extracts the region of data over which we wish to run the model.
2112	!!
2113	!! RECENT CHANGE(S): None
2114	!!
2115	!! MAIN OUTPUT VARIABLE(S):
2116	!!
2117	!! REFERENCE(S) :
2118	!!
2119	!_ ================================================================================================================================
2120	SUBROUTINE forcing_zoom(x_f, x_z)
2121
2122	IMPLICIT NONE
2123	!-
2124	REAL, DIMENSION(iim_full, jjm_full), INTENT(IN) :: x_f
2125	REAL, DIMENSION(iim_zoom, jjm_zoom), INTENT(OUT) :: x_z
2126
2127	INTEGER :: i,j
2128	!-
2129	DO i=1,iim_zoom
2130	DO j=1,jjm_zoom
2131	x_z(i,j) = x_f(i_index(i),j_index(j))
2132	ENDDO
2133	ENDDO
2134	!-
2135	END SUBROUTINE forcing_zoom
2136
2137
2138	!! ==============================================================================================================================\n
2139	!! SUBROUTINE : forcing_vertical_ioipsl
2140	!!
2141	!>\BRIEF
2142	!!
2143	!!\n DESCRIPTION : This subroutine explores the forcing file to decide what information is available
2144	!! on the vertical coordinate.
2145	!!
2146	!! RECENT CHANGE(S): None
2147	!!
2148	!! MAIN OUTPUT VARIABLE(S):
2149	!!
2150	!! REFERENCE(S) :
2151	!!
2152	!_ ================================================================================================================================
2153	SUBROUTINE forcing_vertical_ioipsl(force_id)
2154
2155	INTEGER, INTENT(IN) :: force_id
2156
2157	LOGICAL :: var_exists, vara_exists, varb_exists, varuv_exists
2158	LOGICAL :: foundvar
2159	LOGICAL :: levlegacy
2160
2161	!
2162	!- Set all the defaults
2163	!
2164	zfixed=.FALSE.
2165	zsigma=.FALSE.
2166	zhybrid=.FALSE.
2167	zlevels=.FALSE.
2168	zheight=.FALSE.
2169	zsamelev_uv = .TRUE.
2170	levlegacy = .FALSE.
2171	!
2172	foundvar = .FALSE.
2173	!
2174	!- We have a forcing file to explore so let us see if we find any of the conventions
2175	!- which allow us to find the height of T,Q,U and V.
2176	!
2177	IF ( force_id > 0 ) THEN
2178	!
2179	! Case for sigma levels
2180	!
2181	IF ( .NOT. foundvar ) THEN
2182	CALL flinquery_var(force_id, 'Sigma', var_exists)
2183	CALL flinquery_var(force_id, 'Sigma_uv', varuv_exists)
2184	IF ( var_exists ) THEN
2185	foundvar = .TRUE.
2186	zsigma = .TRUE.
2187	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2188	ENDIF
2189	ENDIF
2190	!
2191	! Case for Hybrid levels
2192	!
2193	IF ( .NOT. foundvar ) THEN
2194	CALL flinquery_var(force_id, 'HybSigA', vara_exists)
2195	IF ( vara_exists ) THEN
2196	CALL flinquery_var(force_id, 'HybSigB', varb_exists)
2197	IF ( varb_exists ) THEN
2198	var_exists=.TRUE.
2199	ELSE
2200	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','Missing the B coefficient for', &
2201	& 'Hybrid vertical levels for T and Q','stop readdim2')
2202	ENDIF
2203	ENDIF
2204	CALL flinquery_var(force_id, 'HybSigA_uv', vara_exists)
2205	IF ( vara_exists ) THEN
2206	CALL flinquery_var(force_id, 'HybSigB_uv', varb_exists)
2207	IF ( varb_exists ) THEN
2208	varuv_exists=.TRUE.
2209	ELSE
2210	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','Missing the B coefficient for', &
2211	& 'Hybrid vertical levels for U and V','stop readdim2')
2212	ENDIF
2213	ENDIF
2214	IF ( var_exists ) THEN
2215	foundvar = .TRUE.
2216	zhybrid = .TRUE.
2217	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2218	ENDIF
2219	ENDIF
2220	!
2221	! Case for levels (i.e. a 2d time varying field with the height in meters)
2222	!
2223	IF ( .NOT. foundvar ) THEN
2224	CALL flinquery_var(force_id, 'Levels', var_exists)
2225	CALL flinquery_var(force_id, 'Levels_uv', varuv_exists)
2226	IF ( var_exists ) THEN
2227	foundvar = .TRUE.
2228	zlevels = .TRUE.
2229	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2230	ENDIF
2231	ENDIF
2232	!
2233	! Case where a fixed height is provided in meters
2234	!
2235	IF ( .NOT. foundvar ) THEN
2236	CALL flinquery_var(force_id, 'Height_Lev1', var_exists)
2237	CALL flinquery_var(force_id, 'Height_Levuv', varuv_exists)
2238	IF ( var_exists ) THEN
2239	foundvar = .TRUE.
2240	zheight = .TRUE.
2241	IF ( varuv_exists ) zsamelev_uv = .FALSE.
2242	ENDIF
2243	ENDIF
2244	!
2245	! Case where a fixed height is provided in meters in the lev variable
2246	!
2247	IF ( .NOT. foundvar ) THEN
2248	CALL flinquery_var(force_id, 'lev', var_exists)
2249	IF ( var_exists ) THEN
2250	foundvar = .TRUE.
2251	zheight = .TRUE.
2252	levlegacy = .TRUE.
2253	ENDIF
2254	ENDIF
2255	!
2256	ENDIF
2257	!
2258	! We found forcing variables so we need to extract the values if we are dealing with constant values (i.e. all
2259	! except the case zlevels
2260	!
2261	IF ( foundvar .AND. .NOT. zlevels ) THEN
2262	!
2263	IF ( zheight ) THEN
2264	!
2265	! Constant height
2266	!
2267	IF ( levlegacy ) THEN
2268	CALL flinget (force_id,'lev',1, 1, 1, 1, 1, 1, zlev_fixed)
2269	ELSE
2270	CALL flinget (force_id,'Height_Lev1',1, 1, 1, 1, 1, 1, zlev_fixed)
2271	IF ( .NOT. zsamelev_uv ) THEN
2272	CALL flinget (force_id,'Height_Levuv',1, 1, 1, 1, 1, 1, zlevuv_fixed)
2273	ENDIF
2274	ENDIF
2275	IF (printlev_loc >= 1) WRITE(numout,*) "forcing_vertical_ioipsl : case ZLEV : Read from forcing file :", &
2276	zlev_fixed, zlevuv_fixed
2277	!
2278	ELSE IF ( zsigma .OR. zhybrid ) THEN
2279	!
2280	! Sigma or hybrid levels
2281	!
2282	IF ( zsigma ) THEN
2283	CALL flinget (force_id,'Sigma',1, 1, 1, 1, 1, 1, zhybrid_b)
2284	zhybrid_a = zero
2285	IF ( .NOT. zsamelev_uv ) THEN
2286	CALL flinget (force_id,'Sigma_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2287	zhybriduv_a = zero
2288	ENDIF
2289	ELSE
2290	CALL flinget (force_id,'HybSigB',1, 1, 1, 1, 1, 1, zhybrid_b)
2291	CALL flinget (force_id,'HybSigA',1, 1, 1, 1, 1, 1, zhybrid_a)
2292	IF ( .NOT. zsamelev_uv ) THEN
2293	CALL flinget (force_id,'HybSigB_uv',1, 1, 1, 1, 1, 1, zhybriduv_b)
2294	CALL flinget (force_id,'HybSigA_uv',1, 1, 1, 1, 1, 1, zhybriduv_a)
2295	ENDIF
2296	ENDIF
2297	IF (printlev_loc >= 1) WRITE(numout,*) "forcing_vertical_ioipsl : case Pressure coordinates : "
2298	IF (printlev_loc >= 1) WRITE(numout,*) "Read from forcing file :", zhybrid_b, zhybrid_a, zhybriduv_b, zhybriduv_a
2299	ELSE
2300	!
2301	! Why are we here ???
2302	!
2303	CALL ipslerr ( 3, 'forcing_vertical_ioipsl','What is the option used to describe the height of', &
2304	& 'the atmospheric forcing ?','Please check your forcing file.')
2305	ENDIF
2306	ENDIF
2307	!
2308	!- We have no forcing file to explore or we did not find anything. So revert back to the run.def and
2309	!- read what has been specified by the user.
2310	!
2311	IF ( force_id < 0 .OR. .NOT. foundvar ) THEN
2312	!
2313	!-
2314	!Config Key = HEIGHT_LEV1
2315	!Config Desc = Height at which T and Q are given
2316	!Config Def = 2.0
2317	!Config If = offline mode
2318	!Config Help = The atmospheric variables (temperature and specific
2319	!Config humidity) are measured at a specific level.
2320	!Config The height of this level is needed to compute
2321	!Config correctly the turbulent transfer coefficients.
2322	!Config Look at the description of the forcing
2323	!Config DATA for the correct value.
2324	!Config Units = [m]
2325	!-
2326	zlev_fixed = 2.0
2327	CALL getin('HEIGHT_LEV1', zlev_fixed)
2328
2329	!-
2330	!Config Key = HEIGHT_LEVW
2331	!Config Desc = Height at which the wind is given
2332	!Config Def = 10.0
2333	!Config If = offline mode
2334	!Config Help = The height at which wind is needed to compute
2335	!Config correctly the turbulent transfer coefficients.
2336	!Config Units= [m]
2337	!-
2338	zlevuv_fixed = 10.0
2339	CALL getin('HEIGHT_LEVW', zlevuv_fixed)
2340
2341	zheight = .TRUE.
2342
2343	IF ( ABS(zlevuv_fixed-zlev_fixed) > EPSILON(zlev_fixed)) THEN
2344	zsamelev_uv = .FALSE.
2345	ELSE
2346	zsamelev_uv = .TRUE.
2347	ENDIF
2348
2349	CALL ipslerr ( 2, 'forcing_vertical_ioipsl','The height of the atmospheric forcing variables', &
2350	& 'was not found in the netCDF file.','Thus the values in run.def were used ... or their defaults.')
2351	ENDIF
2352
2353	END SUBROUTINE forcing_vertical_ioipsl
2354
2355
2356	!! ==============================================================================================================================\n
2357	!! SUBROUTINE : domain_size
2358	!!
2359	!>\BRIEF
2360	!!
2361	!!\n DESCRIPTION :
2362	!!
2363	!! RECENT CHANGE(S): None
2364	!!
2365	!! MAIN OUTPUT VARIABLE(S):
2366	!!
2367	!! REFERENCE(S) :
2368	!!
2369	!_ ================================================================================================================================
2370	SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, &
2371	& iim_f, jjm_f, lon, lat, iim, jjm, iind, jind)
2372
2373	IMPLICIT NONE
2374	!
2375	! ARGUMENTS
2376	!
2377	REAL, INTENT(inout) :: limit_west,limit_east,limit_north,limit_south
2378	INTEGER, INTENT(in) :: iim_f, jjm_f
2379	REAL, INTENT(in) :: lon(iim_f, jjm_f), lat(iim_f, jjm_f)
2380	INTEGER, INTENT(out) :: iim,jjm
2381	INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f)
2382	!
2383	! LOCAL
2384	!
2385	INTEGER :: i, j
2386	REAL :: lolo
2387	LOGICAL :: over_dateline = .FALSE.
2388	!
2389	!
2390	IF ( ( ABS(limit_east) .GT. 180. ) .OR. &
2391	( ABS(limit_west) .GT. 180. ) ) THEN
2392	WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
2393	CALL ipslerr_p (3,'domain_size', &
2394	& 'Longitudes problem.','In run.def file :', &
2395	& 'limit_east > 180. or limit_west > 180.')
2396	ENDIF
2397	!
2398	IF ( limit_west .GT. limit_east ) over_dateline = .TRUE.
2399	!
2400	IF ( ( limit_south .LT. -90. ) .OR. &
2401	( limit_north .GT. 90. ) .OR. &
2402	( limit_south .GE. limit_north ) ) THEN
2403	WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
2404	CALL ipslerr_p (3,'domain_size', &
2405	& 'Latitudes problem.','In run.def file :', &
2406	& 'limit_south < -90. or limit_north > 90. or limit_south >= limit_north')
2407	ENDIF
2408	!
2409	! Here we assume that the grid of the forcing data is regular. Else we would have
2410	! to do more work to find the index table.
2411	!
2412	iim = 0
2413	DO i=1,iim_f
2414	!
2415	lolo = lon(i,1)
2416	IF ( lon(i,1) .GT. 180. ) lolo = lon(i,1) - 360.
2417	IF ( lon(i,1) .LT. -180. ) lolo = lon(i,1) + 360.
2418	!
2419	IF (lon(i,1) < limit_west) iim_g_begin = i+1
2420	IF (lon(i,1) < limit_east) iim_g_end = i
2421	!
2422	IF ( over_dateline ) THEN
2423	IF ( lolo .LE. limit_west .OR. lolo .GE. limit_east ) THEN
2424	iim = iim + 1
2425	iind(iim) = i
2426	ENDIF
2427	ELSE
2428	IF ( lolo .GE. limit_west .AND. lolo .LE. limit_east ) THEN
2429	iim = iim + 1
2430	iind(iim) = i
2431	ENDIF
2432	ENDIF
2433	!
2434	ENDDO
2435	!
2436	jjm = 0
2437	DO j=1,jjm_f
2438	IF (lat(1,j) > limit_north) jjm_g_begin = j+1
2439	IF (lat(1,j) > limit_south) jjm_g_end = j
2440	!
2441	IF ( lat(1,j) .GE. limit_south .AND. lat(1,j) .LE. limit_north) THEN
2442	jjm = jjm + 1
2443	jind(jjm) = j
2444	ENDIF
2445	ENDDO
2446
2447	IF (printlev_loc >= 1) WRITE(numout,*) 'Domain zoom size: iim, jjm = ', iim, jjm
2448
2449	END SUBROUTINE domain_size
2450
2451
2452	!! ==============================================================================================================================\n
2453	!! SUBROUTINE : flinget_buffer
2454	!!
2455	!>\BRIEF
2456	!!
2457	!!\n DESCRIPTION : This subroutine is a wrap of flinget/IOIPSL. The arguments are the same.
2458	!! flinget_buffer will call flinget and buffer the forcing data localy in this subroutine.
2459	!! According to the variable NBUFF set in run.def, several time steps can be read at the same time
2460	!! from the forcing file. If NBUFF=0, the full forcing file is read.
2461	!! The output, data_full, from this subroutine is always only one time step of corresponding to itb.
2462	!! itb must be equal to ite.
2463	!!
2464	!! RECENT CHANGE(S): None
2465	!!
2466	!! MAIN OUTPUT VARIABLE(S):
2467	!!
2468	!! REFERENCE(S) :
2469	!!
2470	!_ ================================================================================================================================
2471	SUBROUTINE flinget_buffer(force_id, varname, iim_full, jjm_full, llm_full, ttm, itb, ite, data_full)
2472
2473	!! Input arguments
2474	INTEGER, INTENT(in) :: force_id !! Id for forcing file
2475	CHARACTER(len=*), INTENT(in) :: varname !! Name of current variable to be read
2476	INTEGER, INTENT(in) :: iim_full, jjm_full, llm_full !! Horizontal and vertical domaine
2477	INTEGER, INTENT(in) :: ttm !! Full lenght of forcing file
2478	INTEGER, INTENT(in) :: itb, ite !! Time step to be read from forcing file. itb must be equal to ite
2479
2480	!! Output argument
2481	REAL, DIMENSION(iim_full, jjm_full), INTENT(out) :: data_full !! Data for time step itb.
2482
2483	!! Define specific type to buffer data together with name and index
2484	TYPE buffer_type
2485	CHARACTER(len=20) :: name !! Name of variable in forcing file
2486	INTEGER :: istart !! Start index of current buffered data
2487	INTEGER :: iend !! End index of current buffered data
2488	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: data !! Data read from forcing file for intervall [istart,iend]
2489	END TYPE buffer_type
2490
2491	!! Local variables
2492	INTEGER, PARAMETER :: maxvar=20 !! Max number of variables to be buffered
2493	TYPE(buffer_type), DIMENSION(maxvar),SAVE :: data_buffer !! Containing all variables and the current buffered data
2494	INTEGER, SAVE :: nbuff !! Number of time steps to be buffered
2495	INTEGER, SAVE :: lastindex=0 !! Current number of variables stored in data_buffer
2496	INTEGER, SAVE :: ttm0 !! Time lenght of forcing file, stored for test purpose
2497	LOGICAL, SAVE :: first=.TRUE. !! First call to this subroutine
2498	INTEGER :: index !! Index in data_buffer for current variable
2499	INTEGER :: i, ierr !! Loop and error variables
2500	INTEGER :: nbuff_new !! Temporary variable for nbuff used in the end of the forcing file
2501
2502	!! 1. Initialization
2503	IF (first) THEN
2504	data_buffer(:)%name='undef'
2505	! Read NBUFF from run.def
2506	! Note that getin_p is not used because this subroutine might be called only by master process
2507
2508	!Config Key = NBUFF
2509	!Config Desc = Number of time steps of data to buffer between each reading of the forcing file
2510	!Config If = OFF_LINE
2511	!Config Help = The full simulation time length will be read if NBUFF equal 0.
2512	!Config NBUFF > 1 can be used for smaller regions or site simulations only.
2513	!Config Def = 1
2514	!Config Units= -
2515
2516	nbuff=1
2517	CALL getin('NBUFF', nbuff)
2518
2519	IF (nbuff == 0 .OR. nbuff >ttm) THEN
2520	! Set nbuff as the full forcing file lenght
2521	nbuff=ttm
2522	ELSE IF (nbuff < 0) THEN
2523	! Negativ nbuff not possible
2524	CALL ipslerr_p(3,'flinget_buffer','NBUFF must be a positiv number','Set NBUFF=0 for full simulation lenght','')
2525	END IF
2526	IF (printlev_loc >= 1) WRITE(numout,*)'flinget_buffer: NBUFF=',nbuff,' number of time step will be buffered'
2527	IF (printlev_loc >= 1) WRITE(numout,*)'flinget_buffer: Choose a lower value for NBUFF if problem with memory'
2528
2529	! Save dimensions to check following timesteps
2530	! ttm is the full lenght of forcing file
2531	ttm0=ttm
2532
2533	first=.FALSE.
2534	END IF
2535
2536	!! 2. Coeherence tests on input arguments
2537	IF (ttm /= ttm0) THEN
2538	WRITE(numout,*)'Problem with ttm=',ttm,' ttm0=',ttm0
2539	CALL ipslerr_p(3,'flinget_buffer','Error with ttm and ttm0','','')
2540	END IF
2541	IF (itb /= ite) THEN
2542	WRITE(numout,*) 'There is a problem. Why is itb not equal ite ?'
2543	WRITE(numout,*) 'itb=',itb,' ite=',ite,' varname=',varname
2544	CALL ipslerr_p(3,'flinget_buffer','ite not equal itb','','')
2545	END IF
2546
2547
2548	!! 3. Find index for current variable
2549	index=0
2550	DO i=1, maxvar
2551	IF ( trim(varname) == data_buffer(i)%name ) THEN
2552	index=i
2553	CYCLE
2554	END IF
2555	END DO
2556
2557	!! 4. Initialize and allocate if necesary the current variable
2558	IF ( index == 0 ) THEN
2559	! The variable was not found
2560	! This must be the first time for current variable
2561	index=lastindex+1
2562	lastindex=index
2563	IF (index > maxvar) CALL ipslerr_p(3,'flinget_buffer','to many variables','maxvar is too small','')
2564
2565	! Initialize the data_buffer for this index
2566	data_buffer(index)%name=trim(varname)
2567	ALLOCATE(data_buffer(index)%data(iim_full,jjm_full,nbuff),stat=ierr)
2568	IF (ierr /= 0) CALL ipslerr_p(3,'flinget_buffer','pb alloc data_buffer%data','for variable=',varname)
2569	data_buffer(index)%istart=0
2570	data_buffer(index)%iend=0
2571	END IF
2572
2573
2574	!! 5. Call flinget if current time step (itb) is outside the buffered intervall
2575	IF (( itb > data_buffer(index)%iend ) .OR. ( itb < data_buffer(index)%istart )) THEN
2576	! itb is not in the time slice previously read or it is the first time to read
2577	! Reading of forcing file will now be done
2578	! First recalculate index to be read
2579	data_buffer(index)%istart = itb
2580	data_buffer(index)%iend = itb + nbuff - 1
2581
2582	! Check and correct if data_buffer(index)%iend is exceeding file size
2583	IF (data_buffer(index)%iend > ttm) THEN
2584	! iend is exceeding the limit of the file. Change iend to the last time step in the file.
2585	data_buffer(index)%iend = ttm
2586
2587	! Calculate a new smaller nbuff
2588	nbuff_new = ttm - itb + 1
2589
2590	! Resize data buffer
2591	DEALLOCATE(data_buffer(index)%data)
2592	ALLOCATE(data_buffer(index)%data(iim_full,jjm_full, nbuff_new), stat=ierr )
2593	IF (ierr /= 0) CALL ipslerr_p(3,'flinget_buffer','pb realloc data_buffer%data with new nbuff','','')
2594	END IF
2595
2596	! WRITE(numout,*) 'Now do flinget for ',varname,', itb=',itb,', istart=',&
2597	! data_buffer(index)%istart,', iend=',data_buffer(index)%iend
2598	CALL flinget (force_id,varname, iim_full, jjm_full, llm_full, ttm, data_buffer(index)%istart, &
2599	data_buffer(index)%iend, data_buffer(index)%data(:,:,:))
2600	END IF
2601
2602	!! 6. Initialize the output variable with data from buffered variable
2603	! Find index for the time step corrsponding to itb in the time slice previously read from forcing file
2604	i=itb-data_buffer(index)%istart+1
2605	! Initialize output variable
2606	data_full(:,:) = data_buffer(index)%data(:,:,i)
2607
2608
2609	END SUBROUTINE flinget_buffer
2610
2611	END MODULE readdim2

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

Context Navigation

source: branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/readdim2.f90

Download in other formats: