Context Navigation

source: branches/publications/ORCHIDEE_GLUC_r6545/src_driver/readdim2.f90 @ 6737

Last change on this file since 6737 was 4719, checked in by albert.jornet, 7 years ago

Merge: from revisions [4491:4695/trunk/ORCHIDEE]

Merge done in [4671:4718/perso/albert.jornet/MICT_MERGE]

Property svn:keywords set to HeadURL Date Author Revision

File size: 105.9 KB

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

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

Download in other formats: