Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_driver/readdim2.f90

Last change on this file was 2088, checked in by yiying.chen, 11 years ago
drivefixed
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 91.1 KB

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

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

Download in other formats: