Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readdim2.f90 @ 6145

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_driver/readdim2.f90 @ 6145

Download in other formats: