Context Navigation

source: tags/ORCHIDEE_2_1/ORCHIDEE/src_driver/readdim2.f90 @ 6593

Last change on this file since 6593 was 5609, checked in by josefine.ghattas, 6 years ago
Integrated correction done in changeset [5329] in branches/ORCHIDEE-MICT/ORCHIDEE. See ticket #446
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 105.5 KB

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

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

Download in other formats: