Context Navigation

localtime.F90 @ 2227

Last change on this file since 2227 was 2227, checked in by abierjon, 5 years ago

Mars GCM:
Bug fixes in localtime.F90 : 1) Following of the correction by EM (r2215), separation of the "valid_range" and "missing_value" indicators, which prevented localtime to write the missing_value attribute in the output file and was making r2215 useless ; 2) Altitude attributes written in the output file are no more specified by default but are read in input file (important for zrecasted files).
AB+EM

File size: 43.0 KB

Line
1	program localtime
2
3	! ----------------------------------------------------------------------------
4	! Program to redistribute and interpolate the variable a the same
5	! local times everywhere
6	! input : diagfi.nc / concat.nc / stats.nc kind of files
7	! author: F. Forget
8	! ----------------------------------------------------------------------------
9
10	implicit none
11
12	include "netcdf.inc" ! NetCDF definitions
13
14	character (len=50) file
15	! file(): input file(s) names(s)
16	character (len=30), dimension(15) :: notconcat
17	! notconcat(): names of the (15) variables that won't be concatenated
18	character (len=50), dimension(:), allocatable :: var
19	! var(): name(s) of variable(s) that will be concatenated
20	character (len=50) :: tmpvar,title,units
21	! tmpvar(): used to temporarily store a variable name
22	! title(): [netcdf] title attribute
23	! units(): [netcdf] units attribute
24	character (len=100) :: filename,vartmp
25	! filename(): output file name
26	! vartmp(): temporary variable name (read from netcdf input file)
27	!character (len=1) :: ccopy
28	! ccpy: 'y' or 'n' answer
29
30	character (len=50) :: altlong_name,altunits,altpositive
31	! altlong_name(): [netcdf] altitude "long_name" attribute
32	! altunits(): [netcdf] altitude "units" attribute
33	! altpositive(): [netcdf] altitude "positive" attribute
34
35	integer :: nid,ierr,miss,validr
36	! nid: [netcdf] file ID #
37	! ierr: [netcdf] subroutine returned error code
38	! miss: [netcdf] subroutine returned error code
39	integer :: i,j,k,inter
40	! for various loops
41	integer :: varid
42	! varid: [netcdf] variable ID #
43	real, dimension(:), allocatable:: lat,lon,alt,ctl,time
44	! lat(): array, stores latitude coordinates
45	! lon(): array, stores longitude coordinates
46	! alt(): array, stores altitude coordinates
47	! ctl(): array, stores controle array
48	! time(): array, stores time coordinates
49	integer :: nbvar,nbvarfile,ndim
50	! nbvar: # of variables to concatenate
51	! nbfile: # number of input file(s)
52	! nbvarfile: total # of variables in an input file
53	! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
54	integer :: latdim,londim,altdim,ctldim,timedim
55	! latdim: [netcdf] "latitude" dim ID
56	! londim: [netcdf] "longitude" dim ID
57	! altdim: [netcdf] "altdim" dim ID
58	! ctldim: [netcdf] "controle" dim ID
59	! timedim: [netcdf] "timedim" dim ID
60	integer :: latvar,lonvar,altvar,ctlvar,timevar
61	! latvar: [netcdf] ID of "latitude" variable
62	! lonvar: [netcdf] ID of "longitude" variable
63	! altvar: [netcdf] ID of "altitude" variable
64	! ctlvar: [netcdf] ID of "controle" variable
65	! timevar: [netcdf] ID of "Time" variable
66	integer :: latlen,lonlen,altlen,ctllen,timelen,timelen_lt,timelen_tot
67	integer :: ilat,ilon,ialt,it
68	! latlen: # of elements of lat() array
69	! lonlen: # of elements of lon() array
70	! altvar: # of elements of alt() array
71	! ctlvar: # of elements of ctl() array
72	! timelen: # of elemnets of time() array
73	! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
74	! timelen_lt: # of elemnets of time() array in output
75	integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
76	integer :: nhour,nsol
77	! nout: [netcdf] output file ID
78	! latdimout: [netcdf] output latitude (dimension) ID
79	! londimout: [netcdf] output longitude (dimension) ID
80	! altdimout: [netcdf] output altitude (dimension) ID
81	! timedimout: [netcdf] output time (dimension) ID
82	! timevarout: [netcdf] ID of output "Time" variable
83	integer :: varidout
84	! varidout: [netcdf] variable ID # (of a variable to write to the output file)
85	integer :: Nnotconcat,var_ok
86	! Nnotconcat: # of (leading)variables that won't be concatenated
87	! var_ok: flag (0 or 1)
88	integer, dimension(4) :: corner,edges,dim
89	! corner: [netcdf]
90	! edges: [netcdf]
91	! dim: [netcdf]
92	real, dimension(:,:,:,:), allocatable :: var3d
93	! var3D(,,,): 4D array to store a field
94
95	real, dimension(:,:,:,:), allocatable :: var3d_lt
96	! var3D_lt(,,,): 4D array to store a field in local time coordinate
97	real, dimension(:), allocatable :: lt_gcm,lt_hour
98	real, dimension(:), allocatable :: lt_out,lt_outc
99	real, dimension(:), allocatable :: var_gcm
100
101	real :: missing
102	!PARAMETER(missing=1E+20)
103	! missing: [netcdf] to handle "missing" values when reading/writing files
104	real, dimension(2) :: valid_range
105	! valid_range(2): [netcdf] interval in which a value is considered valid
106	logical :: stats ! stats=T when reading a "stats" kind of ile
107
108
109	!==============================================================================
110	! 1.1. Get input file name(s)
111	!==============================================================================
112	write(,)
113	write(,) "which file do you want to use? (diagfi... stats... concat...)"
114
115	read(*,'(a50)') file
116
117	if (len_trim(file).eq.0) then
118	write(,) "no file... game over"
119	stop ""
120	endif
121
122	!==============================================================================
123	! 1.3. Open the input file
124	!==============================================================================
125
126	ierr = NF_OPEN(file,NF_NOWRITE,nid)
127	if (ierr.NE.NF_NOERR) then
128	write(,) 'ERROR: Pb opening file '//trim(file)
129	write(,) NF_STRERROR(ierr)
130	stop ""
131	endif
132
133	ierr=NF_INQ_NVARS(nid,nbvarfile)
134	! nbvarfile now set to be the (total) number of variables in file
135	if (ierr.NE.NF_NOERR) then
136	write(,) 'ERROR: Pb with NF_INQ_NVARS'
137	write(,) NF_STRERROR(ierr)
138	stop ""
139	endif
140
141	!==============================================================================
142	! 1.4. List of variables that should not be processed
143	!==============================================================================
144
145	notconcat(1)='Time'
146	notconcat(2)='controle'
147	notconcat(3)='rlonu'
148	notconcat(4)='latitude'
149	notconcat(5)='longitude'
150	notconcat(6)='altitude'
151	notconcat(7)='rlatv'
152	notconcat(8)='aps'
153	notconcat(9)='bps'
154	notconcat(10)='ap'
155	notconcat(11)='bp'
156	notconcat(12)='cu'
157	notconcat(13)='cv'
158	notconcat(14)='aire'
159	notconcat(15)='phisinit'
160
161	!==============================================================================
162	! 1.5. Get (and check) list of variables to process
163	!==============================================================================
164	write(,)
165	Nnotconcat=0
166	do i=1,nbvarfile
167	ierr=NF_INQ_VARNAME(nid,i,vartmp)
168	! vartmp now contains the "name" of variable of ID # i
169	var_ok=0
170	do inter=1,15
171	if (vartmp.eq.notconcat(inter)) then
172	var_ok=1
173	Nnotconcat=Nnotconcat+1
174	endif
175	enddo
176	if (var_ok.eq.0) write(,) trim(vartmp)
177	enddo
178
179	! Nnotconcat: # of variables that won't be processed
180	! nbvarfile: total # of variables in file
181	allocate(var(nbvarfile-Nnotconcat),stat=ierr)
182	if (ierr.ne.0) then
183	write(,) "Error: failed allocation of var(nbvarfile-Nnotconcat)"
184	write(,) " nbvarfile=",nbvarfile
185	write(,) " Nnotconcat=",Nnotconcat
186	stop
187	endif
188
189
190	write(,)
191	write(,) "which variables do you want to redistribute ?"
192	write(,) "all / list of <variables> (separated by <Enter>s)"
193	write(,) "(an empty line , i.e: just <Enter>, implies end of list)"
194	nbvar=0
195	read(*,'(a50)') tmpvar
196	do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
197	nbvar=nbvar+1
198	var(nbvar)=tmpvar
199	read(*,'(a50)') tmpvar
200	enddo
201
202	if (tmpvar=="all") then
203	nbvar=nbvarfile-Nnotconcat
204	do j=Nnotconcat+1,nbvarfile
205	ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
206	enddo
207	! Variables names from the file are stored in var()
208	nbvar=nbvarfile-Nnotconcat
209	do i=1,nbvar
210	ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
211	write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
212	enddo
213	else if(nbvar==0) then
214	write(,) "no variable... game over"
215	stop ""
216	endif ! of if (tmpvar=="all")
217
218	!==============================================================================
219	! 1.6. Get output file name
220	!==============================================================================
221	filename=file(1:len_trim(file)-3)//"_LT.nc"
222	write(,) trim(filename)
223
224	!==============================================================================
225	! 2.1. Open input file
226	!==============================================================================
227
228	write(,)
229	write(,) "opening "//trim(file)//"..."
230	ierr = NF_OPEN(file,NF_NOWRITE,nid)
231	if (ierr.NE.NF_NOERR) then
232	write(,) 'ERROR: Pb opening file '//trim(file)
233	write(,) NF_STRERROR(ierr)
234	stop ""
235	endif
236
237	!==============================================================================
238	! 2.2. Read (and check) dimensions of variables from input file
239	!==============================================================================
240
241	ierr=NF_INQ_DIMID(nid,"latitude",latdim)
242	if (ierr.NE.NF_NOERR) then
243	write(,) 'ERROR: Dimension <latitude> is missing in file '//trim(file)
244	stop ""
245	endif
246	ierr=NF_INQ_VARID(nid,"latitude",latvar)
247	if (ierr.NE.NF_NOERR) then
248	write(,) 'ERROR: Field <latitude> is missing in file '//trim(file)
249	stop ""
250	endif
251	ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
252	! write(,) "latlen: ",latlen
253
254	ierr=NF_INQ_DIMID(nid,"longitude",londim)
255	if (ierr.NE.NF_NOERR) then
256	write(,) 'ERROR: Dimension <longitude> is missing in file '//trim(file)
257	stop ""
258	endif
259	ierr=NF_INQ_VARID(nid,"longitude",lonvar)
260	if (ierr.NE.NF_NOERR) then
261	write(,) 'ERROR: Field <longitude> is missing in file'//trim(file)
262	stop ""
263	endif
264	ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
265	! write(,) "lonlen: ",lonlen
266
267	ierr=NF_INQ_DIMID(nid,"altitude",altdim)
268	if (ierr.NE.NF_NOERR) then
269	write(,) 'ERROR: Dimension <altitude> is missing in file '//trim(file)
270	stop ""
271	endif
272	ierr=NF_INQ_VARID(nid,"altitude",altvar)
273	if (ierr.NE.NF_NOERR) then
274	write(,) 'ERROR: Field <altitude> is missing in file'//trim(file)
275	stop ""
276	endif
277	ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
278	! write(,) "altlen: ",altlen
279
280	ierr=NF_INQ_DIMID(nid,"index",ctldim)
281	if (ierr.NE.NF_NOERR) then
282	write(,) 'ERROR: Dimension <index> is missing in file '//trim(file)
283	stop ""
284	endif
285	ierr=NF_INQ_VARID(nid,"controle",ctlvar)
286	if (ierr.NE.NF_NOERR) then
287	write(,) 'Field <controle> is missing in file'//trim(file)
288	ctllen=0
289	!stop ""
290	else
291	ierr=NF_INQ_DIMLEN(nid,ctldim,ctllen)
292	endif
293	! write(,) "controle: ",controle
294
295	!==============================================================================
296	! 2.3. Read (and check compatibility of) dimensions of
297	! variables from input file
298	!==============================================================================
299
300	! First call; initialize/allocate
301	allocate(lat(latlen),stat=ierr)
302	if (ierr.ne.0) then
303	write(,) "Error: failed to allocate lat(latlen)"
304	stop
305	endif
306	allocate(lon(lonlen),stat=ierr)
307	if (ierr.ne.0) then
308	write(,) "Error: failed to allocate lon(lonlen)"
309	stop
310	endif
311	allocate(alt(altlen),stat=ierr)
312	if (ierr.ne.0) then
313	write(,) "Error: failed to allocate alt(altlen)"
314	stop
315	endif
316	allocate(ctl(ctllen),stat=ierr)
317	if (ierr.ne.0) then
318	write(,) "Error: failed to allocate ctl(ctllen)"
319	stop
320	endif
321
322	#ifdef NC_DOUBLE
323	ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
324	#else
325	ierr = NF_GET_VAR_REAL(nid,latvar,lat)
326	#endif
327	if (ierr.ne.0) then
328	write(,) "Error: failed to load latitude"
329	write(,) NF_STRERROR(ierr)
330	stop
331	endif
332
333	#ifdef NC_DOUBLE
334	ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
335	#else
336	ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
337	#endif
338	if (ierr.ne.0) then
339	write(,) "Error: failed to load longitude"
340	write(,) NF_STRERROR(ierr)
341	stop
342	endif
343
344	#ifdef NC_DOUBLE
345	ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
346	#else
347	ierr = NF_GET_VAR_REAL(nid,altvar,alt)
348	#endif
349	if (ierr.ne.0) then
350	write(,) "Error: failed to load altitude"
351	write(,) NF_STRERROR(ierr)
352	stop
353	endif
354	! Get altitude attributes to handle files with any altitude type
355	ierr=nf_get_att_text(nid,altvar,'long_name',altlong_name)
356	ierr=nf_get_att_text(nid,altvar,'units',altunits)
357	ierr=nf_get_att_text(nid,altvar,'positive',altpositive)
358
359	if (ctllen .ne. 0) then
360	#ifdef NC_DOUBLE
361	ierr = NF_GET_VAR_DOUBLE(nid,ctlvar,ctl)
362	#else
363	ierr = NF_GET_VAR_REAL(nid,ctlvar,ctl)
364	#endif
365	if (ierr.ne.0) then
366	write(,) "Error: failed to load controle"
367	write(,) NF_STRERROR(ierr)
368	stop
369	endif
370	endif ! of if (ctllen .ne. 0)
371
372	!==============================================================================
373	! 2.4. Handle "Time" dimension from input file
374	!==============================================================================
375
376	!==============================================================================
377	! 2.4.0 Read "Time" dimension from input file
378	!==============================================================================
379	ierr=NF_INQ_DIMID(nid,"Time",timedim)
380	if (ierr.NE.NF_NOERR) then
381	write(,) 'ERROR: Dimension <Time> is missing in file'//trim(file)
382	stop ""
383	endif
384	ierr=NF_INQ_VARID(nid,"Time",timevar)
385	if (ierr.NE.NF_NOERR) then
386	write(,) 'ERROR: Field <Time> is missing in file'//trim(file)
387	stop ""
388	endif
389	ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
390	! write(,) "timelen: ",timelen
391
392	! Sanity Check: Does "Time" has a "title" attribute
393	! and is it "Solar longitude" ?
394	ierr=nf_get_att_text(nid,timevar,"title",title)
395	if ((ierr.EQ.NF_NOERR).and.(title.eq."Solar longitude")) then
396	write(,) "ERROR: Time axis in input file is in Solar Longitude!"
397	write(,) " localtime requires sols as time axis!"
398	write(,) " Might as well stop here."
399	stop
400	endif
401
402	! allocate time() array and fill it with values from input file
403	allocate(time(timelen),stat=ierr)
404	if (ierr.ne.0) then
405	write(,) "Error: failed to allocate time(timelen)"
406	stop
407	endif
408
409	#ifdef NC_DOUBLE
410	ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
411	#else
412	ierr = NF_GET_VAR_REAL(nid,timevar,time)
413	#endif
414	if (ierr.NE.NF_NOERR) then
415	write(,) "Error , failed to load Time"
416	write(,) NF_STRERROR(ierr)
417	stop
418	endif
419
420
421	!==============================================================================
422	! 2.4.1 Select local times to be saved "Time" in output file
423	!==============================================================================
424	write(,) 'number of local time to be used ?'
425	read(,) nhour
426	allocate(lt_hour(nhour),stat=ierr)
427	if (ierr.ne.0) then
428	write(,) "Error: failed to allocate lt_hour(nhour)"
429	stop
430	endif
431	write(,) 'list of selected local time (0<t<24)'
432	do it=1,nhour
433	read(,) lt_hour(it)
434	end do
435
436	if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
437	write(,) 'We detect a ""stats"" file'
438	stats=.true.
439	timelen_lt=nhour
440	allocate(lt_out(timelen_lt),stat=ierr)
441	if (ierr.ne.0) then
442	write(,) "Error: failed to allocate lt_hour(nhour)"
443	stop
444	endif
445	do it=1,nhour
446	lt_out(it)=lt_hour(it)/24.
447	end do
448	! We rewrite the time from "stats" from 0 to 1 sol...
449	do it=1,timelen ! loop temps in
450	time(it) = time(it)/24.
451	end do
452	nsol =1
453	else ! case of a diagfi or concatnc file
454	stats=.false.
455	! Number of sol in input file
456	nsol = int(time(timelen)) - int(time(1))
457	timelen_lt=nhour*nsol
458	allocate(lt_out(timelen_lt),stat=ierr)
459	if (ierr.ne.0) then
460	write(,) "Error: failed to allocate lt_hour(nhour)"
461	stop
462	endif
463	i=0
464	do k=1,nsol
465	do it=1,nhour
466	i=i+1
467	lt_out(i)=int(time(1)) + k-1 + lt_hour(it)/24.
468	end do
469	end do
470	end if
471
472	if (nsol.gt.1) then
473	timelen_tot=timelen
474	else
475	! if only 1 sol available, we must add 1 timestep for the interpolation
476	timelen_tot=timelen+1
477	endif
478	allocate(lt_gcm(timelen_tot),stat=ierr)
479	if (ierr.ne.0) then
480	write(,) "Error: failed to allocate lt_gcm(timelen_tot)"
481	stop
482	endif
483	allocate(var_gcm(timelen_tot),stat=ierr)
484	if (ierr.ne.0) then
485	write(,) "Error: failed to allocate var_gcm(timelen_tot)"
486	stop
487	endif
488	allocate(lt_outc(timelen_lt),stat=ierr)
489	if (ierr.ne.0) then
490	write(,) "Error: failed to allocate lt_outc(timelen_lt)"
491	stop
492	endif
493
494	!==============================================================================
495	! 2.4.1.5 Initiate dimension in output file
496	!==============================================================================
497
498
499	! Initialize output file's lat,lon,alt and time dimensions
500	call initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
501	altlong_name,altunits,altpositive,&
502	nout,latdimout,londimout,altdimout,timedimout,timevarout)
503
504	! Initialize output file's aps,bps and phisinit variables
505	call init2(nid,lonlen,latlen,altlen,altdim,&
506	nout,londimout,latdimout,altdimout)
507
508
509	!==============================================================================
510	! 2.4.2 Write/extend "Time" dimension/values in output file
511	!==============================================================================
512
513	if (stats) then
514
515	do it=1,timelen_lt
516	#ifdef NC_DOUBLE
517	ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_hour(it))
518	#else
519	ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
520	#endif
521	enddo
522	else
523	do it=1,timelen_lt
524	#ifdef NC_DOUBLE
525	ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_out(it))
526	#else
527	ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
528	#endif
529	if (ierr.NE.NF_NOERR) then
530	write(,) "Error , failed to write Time"
531	stop
532	endif
533	enddo
534	end if
535
536	!==============================================================================
537	! 2.5 Read/write variables
538	!==============================================================================
539
540	do j=1,nbvar ! loop on variables to read/write
541
542	!==============================================================================
543	! 2.5.1 Check that variable to be read existe in input file
544	!==============================================================================
545
546	write(,) "variable ",var(j)
547	ierr=nf_inq_varid(nid,var(j),varid)
548	if (ierr.NE.NF_NOERR) then
549	write(,) 'ERROR: Field <',var(j),'> not found in file'//file
550	stop ""
551	endif
552	ierr=nf_inq_varndims(nid,varid,ndim)
553
554	!==============================================================================
555	! 2.5.2 Prepare things in order to read/write the variable
556	!==============================================================================
557
558	! build dim(),corner() and edges() arrays
559	! and allocate var3d() array
560	if (ndim==3) then
561	allocate(var3d(lonlen,latlen,timelen,1))
562	allocate(var3d_lt(lonlen,latlen,timelen_lt,1))
563	dim(1)=londimout
564	dim(2)=latdimout
565	dim(3)=timedimout
566
567	! start indexes (where data values will be written)
568	corner(1)=1
569	corner(2)=1
570	corner(3)=1
571	corner(4)=1
572
573	! length (along dimensions) of block of data to be written
574	edges(1)=lonlen
575	edges(2)=latlen
576	edges(3)=timelen_lt
577	edges(4)=1
578
579	else if (ndim==4) then
580	allocate(var3d(lonlen,latlen,altlen,timelen))
581	allocate(var3d_lt(lonlen,latlen,altlen,timelen_lt))
582	dim(1)=londimout
583	dim(2)=latdimout
584	dim(3)=altdimout
585	dim(4)=timedimout
586
587	! start indexes (where data values will be written)
588	corner(1)=1
589	corner(2)=1
590	corner(3)=1
591	corner(4)=1
592
593	! length (along dimensions) of block of data to be written
594	edges(1)=lonlen
595	edges(2)=latlen
596	edges(3)=altlen
597	edges(4)=timelen_lt
598	endif
599
600	units=" "
601	title=" "
602	ierr=nf_get_att_text(nid,varid,"title",title)
603	ierr=nf_get_att_text(nid,varid,"units",units)
604	call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
605
606
607
608	!==============================================================================
609	! 2.5.3 Read from input file
610	!==============================================================================
611
612	#ifdef NC_DOUBLE
613	ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
614	miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
615	validr=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
616	#else
617	ierr = NF_GET_VAR_REAL(nid,varid,var3d)
618	miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
619	validr=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
620	#endif
621
622	!==============================================================================
623	! 2.5.4 Read from input file
624	!==============================================================================
625	! Interpolation in local time :
626
627	do ilon=1,lonlen
628	! write(,) 'processing longitude =', lon(ilon)
629	! Local time at each longitude :
630	do it=1,timelen ! loop temps in
631	lt_gcm(it) = time(it) +0.5*lon(ilon)/180.
632	end do
633	if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
634
635	do it=1,timelen_lt ! loop time local time
636	lt_outc(it)=lt_out(it)
637	! Correction to use same local time previous or following
638	! day where data are missing :
639
640	if(lt_outc(it).gt.lt_gcm(timelen_tot))lt_outc(it)=lt_out(it)-1
641	if(lt_outc(it).lt.lt_gcm(1))lt_outc(it)=lt_out(it)+1
642	end do
643
644	if (ndim==3) then
645	do ilat=1,latlen
646	do it=1,timelen ! loop temps in
647	var_gcm(it) = var3d(ilon,ilat,it,1)
648	end do
649	if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
650	do it=1,timelen_lt ! loop time local time
651	call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
652	missing,lt_gcm,var_gcm,timelen_tot)
653	end do
654	enddo
655
656	else if (ndim==4) then
657	do ilat=1,latlen
658	do ialt=1,altlen
659	do it=1,timelen ! loop temps in
660	var_gcm(it) = var3d(ilon,ilat,ialt,it)
661	end do
662	if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
663	do it=1,timelen_lt ! loop time local time
664	call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
665	missing,lt_gcm,var_gcm,timelen_tot)
666	end do
667	enddo
668	enddo
669	end if
670
671	end do
672
673
674	!==============================================================================
675	! 2.5.5 write (append) to the output file
676	!==============================================================================
677
678	#ifdef NC_DOUBLE
679	ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
680	#else
681	ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
682	#endif
683
684	if (ierr.ne.NF_NOERR) then
685	write(,) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
686	stop ""
687	endif
688
689	! In case there is a "valid_range" attribute
690	! Write "valid_range" and "missing_value" attributes in output file
691	if (miss.eq.NF_NOERR) then
692	call missing_value(nout,varidout,validr,miss,valid_range,missing)
693	endif
694
695	! free var3d() array
696	deallocate(var3d)
697	deallocate(var3d_lt)
698
699	enddo ! of do j=1,nbvar
700
701	deallocate(time)
702	deallocate(lt_gcm)
703	deallocate(lt_out)
704	deallocate(lt_outc)
705	deallocate(var_gcm)
706
707	! Close input file
708	ierr=nf_close(nid)
709
710	! Close output file
711	ierr=NF_CLOSE(nout)
712
713	contains
714
715	!******************************************************************************
716	Subroutine initiate(filename,lat,lon,alt,ctl,latlen,lonlen,altlen,ctllen,&
717	altlong_name,altunits,altpositive,&
718	nout,latdimout,londimout,altdimout,timedimout,timevarout)
719	!==============================================================================
720	! Purpose:
721	! Create and initialize a data file (NetCDF format)
722	!==============================================================================
723	! Remarks:
724	! The NetCDF file (created in this subroutine) remains open
725	!==============================================================================
726
727	implicit none
728
729	include "netcdf.inc" ! NetCDF definitions
730
731	!==============================================================================
732	! Arguments:
733	!==============================================================================
734	character (len=*), intent(in):: filename
735	! filename(): the file's name
736	integer,intent(in) ::latlen,lonlen,altlen,ctllen
737	real, intent(in):: lat(latlen)
738	! lat(): latitude
739	real, intent(in):: lon(lonlen)
740	! lon(): longitude
741	real, intent(in):: alt(altlen)
742	! alt(): altitude
743	real, intent(in):: ctl(ctllen)
744	! ctl(): controle
745	character (len=*), intent(in) :: altlong_name
746	! altlong_name(): [netcdf] altitude "long_name" attribute
747	character (len=*), intent(in) :: altunits
748	! altunits(): [netcdf] altitude "units" attribute
749	character (len=*), intent(in) :: altpositive
750	! altpositive(): [netcdf] altitude "positive" attribute
751	integer, intent(out):: nout
752	! nout: [netcdf] file ID
753	integer, intent(out):: latdimout
754	! latdimout: [netcdf] lat() (i.e.: latitude) ID
755	integer, intent(out):: londimout
756	! londimout: [netcdf] lon() ID
757	integer, intent(out):: altdimout
758	! altdimout: [netcdf] alt() ID
759	integer, intent(out):: timedimout
760	! timedimout: [netcdf] "Time" ID
761	integer, intent(out):: timevarout
762	! timevarout: [netcdf] "Time" (considered as a variable) ID
763
764	!==============================================================================
765	! Local variables:
766	!==============================================================================
767	!integer :: latdim,londim,altdim,timedim
768	integer :: nvarid,ierr
769	integer :: ctldimout
770	! nvarid: [netcdf] ID of a variable
771	! ierr: [netcdf] return error code (from called subroutines)
772
773	!==============================================================================
774	! 1. Create (and open) output file
775	!==============================================================================
776	write(,) "creating "//trim(adjustl(filename))//'...'
777	ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
778	! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
779	if (ierr.NE.NF_NOERR) then
780	WRITE(,)'ERROR: Impossible to create the file ',trim(filename)
781	write(,) NF_STRERROR(ierr)
782	stop ""
783	endif
784
785	!==============================================================================
786	! 2. Define/write "dimensions" and get their IDs
787	!==============================================================================
788
789	ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
790	if (ierr.NE.NF_NOERR) then
791	WRITE(,)'initiate: error failed to define dimension <latitude>'
792	write(,) NF_STRERROR(ierr)
793	stop ""
794	endif
795	ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
796	if (ierr.NE.NF_NOERR) then
797	WRITE(,)'initiate: error failed to define dimension <longitude>'
798	write(,) NF_STRERROR(ierr)
799	stop ""
800	endif
801	ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
802	if (ierr.NE.NF_NOERR) then
803	WRITE(,)'initiate: error failed to define dimension <altitude>'
804	write(,) NF_STRERROR(ierr)
805	stop ""
806	endif
807	if (size(ctl).ne.0) then
808	ierr = NF_DEF_DIM(nout, "index", ctllen, ctldimout)
809	if (ierr.NE.NF_NOERR) then
810	WRITE(,)'initiate: error failed to define dimension <index>'
811	write(,) NF_STRERROR(ierr)
812	stop ""
813	endif
814	endif
815	ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
816	if (ierr.NE.NF_NOERR) then
817	WRITE(,)'initiate: error failed to define dimension <Time>'
818	write(,) NF_STRERROR(ierr)
819	stop ""
820	endif
821
822	! End netcdf define mode
823	ierr = NF_ENDDEF(nout)
824	if (ierr.NE.NF_NOERR) then
825	WRITE(,)'initiate: error failed to switch out of define mode'
826	write(,) NF_STRERROR(ierr)
827	stop ""
828	endif
829
830	!==============================================================================
831	! 3. Write "Time" and its attributes
832	!==============================================================================
833
834	call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
835	(/timedimout/),timevarout,ierr)
836
837	!==============================================================================
838	! 4. Write "latitude" (data and attributes)
839	!==============================================================================
840
841	call def_var(nout,"latitude","latitude","degrees_north",1,&
842	(/latdimout/),nvarid,ierr)
843
844	#ifdef NC_DOUBLE
845	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
846	#else
847	ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
848	#endif
849	if (ierr.NE.NF_NOERR) then
850	WRITE(,)'initiate: error failed writing variable <latitude>'
851	write(,) NF_STRERROR(ierr)
852	stop ""
853	endif
854
855	!==============================================================================
856	! 4. Write "longitude" (data and attributes)
857	!==============================================================================
858
859	call def_var(nout,"longitude","East longitude","degrees_east",1,&
860	(/londimout/),nvarid,ierr)
861
862	#ifdef NC_DOUBLE
863	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
864	#else
865	ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
866	#endif
867	if (ierr.NE.NF_NOERR) then
868	WRITE(,)'initiate: error failed writing variable <longitude>'
869	write(,) NF_STRERROR(ierr)
870	stop ""
871	endif
872
873	!==============================================================================
874	! 5. Write "altitude" (data and attributes)
875	!==============================================================================
876
877	! Switch to netcdf define mode
878	ierr = NF_REDEF (nout)
879
880	#ifdef NC_DOUBLE
881	ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
882	#else
883	ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
884	#endif
885
886	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
887	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
888	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
889
890	! End netcdf define mode
891	ierr = NF_ENDDEF(nout)
892
893	#ifdef NC_DOUBLE
894	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
895	#else
896	ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
897	#endif
898	if (ierr.NE.NF_NOERR) then
899	WRITE(,)'initiate: error failed writing variable <altitude>'
900	write(,) NF_STRERROR(ierr)
901	stop ""
902	endif
903
904	!==============================================================================
905	! 6. Write "controle" (data and attributes)
906	!==============================================================================
907
908	if (size(ctl).ne.0) then
909	! Switch to netcdf define mode
910	ierr = NF_REDEF (nout)
911
912	#ifdef NC_DOUBLE
913	ierr = NF_DEF_VAR (nout,"controle",NF_DOUBLE,1,ctldimout,nvarid)
914	#else
915	ierr = NF_DEF_VAR (nout,"controle",NF_FLOAT,1,ctldimout,nvarid)
916	#endif
917
918	ierr = NF_PUT_ATT_TEXT (nout,nvarid,"title",18,"Control parameters")
919
920	! End netcdf define mode
921	ierr = NF_ENDDEF(nout)
922
923	#ifdef NC_DOUBLE
924	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,ctl)
925	#else
926	ierr = NF_PUT_VAR_REAL (nout,nvarid,ctl)
927	#endif
928	if (ierr.NE.NF_NOERR) then
929	WRITE(,)'initiate: error failed writing variable <controle>'
930	write(,) NF_STRERROR(ierr)
931	stop ""
932	endif
933	endif
934
935	end Subroutine initiate
936	!******************************************************************************
937	subroutine init2(infid,lonlen,latlen,altlen,altdimid, &
938	outfid,londimout,latdimout,altdimout)
939	!==============================================================================
940	! Purpose:
941	! Copy aps(), bps() and phisinit() from input file to output file
942	!==============================================================================
943	! Remarks:
944	! The NetCDF files must be open
945	!==============================================================================
946
947	implicit none
948
949	include "netcdf.inc" ! NetCDF definitions
950
951	!==============================================================================
952	! Arguments:
953	!==============================================================================
954	integer, intent(in) :: infid ! NetCDF output file ID
955	integer, intent(in) :: lonlen ! # of grid points along longitude
956	integer, intent(in) :: latlen ! # of grid points along latitude
957	integer, intent(in) :: altlen ! # of grid points along altitude
958	integer, intent(in) :: altdimid ! ID of altitude dimension
959	integer, intent(in) :: outfid ! NetCDF output file ID
960	integer, intent(in) :: londimout ! longitude dimension ID
961	integer, intent(in) :: latdimout ! latitude dimension ID
962	integer, intent(in) :: altdimout ! altitude dimension ID
963	!==============================================================================
964	! Local variables:
965	!==============================================================================
966	real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
967	real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
968	integer :: apsid,bpsid,phisinitid
969	integer :: ierr
970	integer :: tmpdimid ! temporary dimension ID
971	integer :: tmpvarid ! temporary variable ID
972	logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
973
974
975	!==============================================================================
976	! 1. Read data from input file
977	!==============================================================================
978
979	! hybrid coordinate aps
980	allocate(aps(altlen),stat=ierr)
981	if (ierr.ne.0) then
982	write(,) "init2: failed to allocate aps(altlen)"
983	stop
984	endif
985
986	ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
987	if (ierr.ne.NF_NOERR) then
988	write(,) "oops Failed to get aps ID. OK"
989	aps_ok=.false.
990	else
991	! Check that aps() is the right size (it most likely won't be if
992	! zrecast has be used to generate the input file)
993	ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
994	if (tmpdimid.ne.altdimid) then
995	write(,) "init2: wrong dimension size for aps(), skipping it ..."
996	aps_ok=.false.
997	else
998	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
999	if (ierr.ne.NF_NOERR) then
1000	stop "init2 error: Failed reading aps"
1001	endif
1002	aps_ok=.true.
1003	endif
1004	endif
1005
1006	! hybrid coordinate bps
1007	allocate(bps(altlen),stat=ierr)
1008	if (ierr.ne.0) then
1009	write(,) "init2: failed to allocate bps(altlen)"
1010	stop
1011	endif
1012
1013	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
1014	if (ierr.ne.NF_NOERR) then
1015	write(,) "oops: Failed to get bps ID. OK"
1016	bps_ok=.false.
1017	else
1018	! Check that bps() is the right size
1019	ierr=NF_INQ_VARDIMID(infid,tmpvarid,tmpdimid)
1020	if (tmpdimid.ne.altdimid) then
1021	write(,) "init2: wrong dimension size for bps(), skipping it ..."
1022	bps_ok=.false.
1023	else
1024	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
1025	if (ierr.ne.NF_NOERR) then
1026	stop "init2 Error: Failed reading bps"
1027	endif
1028	bps_ok=.true.
1029	endif
1030	endif
1031
1032	! ground geopotential phisinit
1033	allocate(phisinit(lonlen,latlen),stat=ierr)
1034	if (ierr.ne.0) then
1035	write(,) "init2: failed to allocate phisinit(lonlen,latlen)"
1036	stop
1037	endif
1038	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1039	if (ierr.ne.NF_NOERR) then
1040	write(,) "Failed to get phisinit ID. OK"
1041	phisinit = 0.
1042	phis = .false.
1043	else
1044	ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
1045	if (ierr.ne.NF_NOERR) then
1046	stop "init2 Error: Failed reading phisinit"
1047	endif
1048	phis = .true.
1049	endif
1050
1051
1052	!==============================================================================
1053	! 2. Write
1054	!==============================================================================
1055
1056	!==============================================================================
1057	! 2.2. Hybrid coordinates aps() and bps()
1058	!==============================================================================
1059
1060	IF (aps_ok) then
1061	! define aps
1062	call def_var(outfid,"aps","hybrid pressure at midlayers"," ",1,&
1063	(/altdimout/),apsid,ierr)
1064	if (ierr.ne.NF_NOERR) then
1065	stop "Error: Failed to def_var aps"
1066	endif
1067
1068	! write aps
1069	#ifdef NC_DOUBLE
1070	ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
1071	#else
1072	ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
1073	#endif
1074	if (ierr.ne.NF_NOERR) then
1075	stop "Error: Failed to write aps"
1076	endif
1077	ENDIF ! of IF (aps_ok)
1078
1079	IF (bps_ok) then
1080	! define bps
1081	call def_var(outfid,"bps","hybrid sigma at midlayers"," ",1,&
1082	(/altdimout/),bpsid,ierr)
1083	if (ierr.ne.NF_NOERR) then
1084	stop "Error: Failed to def_var bps"
1085	endif
1086
1087	! write bps
1088	#ifdef NC_DOUBLE
1089	ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
1090	#else
1091	ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
1092	#endif
1093	if (ierr.ne.NF_NOERR) then
1094	stop "Error: Failed to write bps"
1095	endif
1096	ENDIF ! of IF (bps_ok)
1097
1098	!==============================================================================
1099	! 2.2. phisinit()
1100	!==============================================================================
1101
1102	IF (phis) THEN
1103
1104	!define phisinit
1105	call def_var(outfid,"phisinit","Ground level geopotential"," ",2,&
1106	(/londimout,latdimout/),phisinitid,ierr)
1107	if (ierr.ne.NF_NOERR) then
1108	stop "Error: Failed to def_var phisinit"
1109	endif
1110
1111	! write phisinit
1112	#ifdef NC_DOUBLE
1113	ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
1114	#else
1115	ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
1116	#endif
1117	if (ierr.ne.NF_NOERR) then
1118	stop "Error: Failed to write phisinit"
1119	endif
1120
1121	ENDIF ! of IF (phis)
1122
1123
1124	! Cleanup
1125	deallocate(aps)
1126	deallocate(bps)
1127	deallocate(phisinit)
1128
1129	end subroutine init2
1130	!******************************************************************************
1131	subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
1132	!==============================================================================
1133	! Purpose: Write a variable (i.e: add a variable to a dataset)
1134	! called "name"; along with its attributes "title", "units"...
1135	! to a file (following the NetCDF format)
1136	!==============================================================================
1137	! Remarks:
1138	! The NetCDF file must be open
1139	!==============================================================================
1140
1141	implicit none
1142
1143	include "netcdf.inc" ! NetCDF definitions
1144
1145	!==============================================================================
1146	! Arguments:
1147	!==============================================================================
1148	integer, intent(in) :: nid
1149	! nid: [netcdf] file ID #
1150	character (len=*), intent(in) :: name
1151	! name(): [netcdf] variable's name
1152	character (len=*), intent(in) :: title
1153	! title(): [netcdf] variable's "title" attribute
1154	character (len=*), intent(in) :: units
1155	! unit(): [netcdf] variable's "units" attribute
1156	integer, intent(in) :: nbdim
1157	! nbdim: number of dimensions of the variable
1158	integer, dimension(nbdim), intent(in) :: dim
1159	! dim(nbdim): [netcdf] dimension(s) ID(s)
1160	integer, intent(out) :: nvarid
1161	! nvarid: [netcdf] ID # of the variable
1162	integer, intent(out) :: ierr
1163	! ierr: [netcdf] subroutines returned error code
1164
1165	! Switch to netcdf define mode
1166	ierr=NF_REDEF(nid)
1167
1168	! Insert the definition of the variable
1169	#ifdef NC_DOUBLE
1170	ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
1171	#else
1172	ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1173	#endif
1174
1175	! Write the attributes
1176	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
1177	ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1178
1179	! End netcdf define mode
1180	ierr=NF_ENDDEF(nid)
1181
1182	end subroutine def_var
1183	!******************************************************************************
1184	subroutine missing_value(nout,nvarid,validr,miss,valid_range,missing)
1185	!==============================================================================
1186	! Purpose:
1187	! Write "valid_range" and "missing_value" attributes (of a given
1188	! variable) to a netcdf file
1189	!==============================================================================
1190	! Remarks:
1191	! NetCDF file must be open
1192	! Variable (nvarid) ID must be set
1193	!==============================================================================
1194
1195	implicit none
1196
1197	include "netcdf.inc" ! NetCDF definitions
1198
1199	!==============================================================================
1200	! Arguments:
1201	!==============================================================================
1202	integer, intent(in) :: nout
1203	! nout: [netcdf] file ID #
1204	integer, intent(in) :: nvarid
1205	! varid: [netcdf] variable ID #
1206	integer, intent(in) :: validr
1207	! validr : [netcdf] routines return code for "valid_range" attribute
1208	integer, intent(in) :: miss
1209	! miss : [netcdf] routines return code for "missing_value" attribute
1210	real, dimension(2), intent(in) :: valid_range
1211	! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1212	real, intent(in) :: missing
1213	! missing: [netcdf] "missing_value" attribute
1214
1215	!==============================================================================
1216	! Local variables:
1217	!==============================================================================
1218	integer :: ierr
1219	! ierr: [netcdf] subroutine returned error code
1220	! INTEGER nout,nvarid,ierr
1221	! REAL missing
1222	! REAL valid_range(2)
1223
1224	! Switch to netcdf dataset define mode
1225	ierr = NF_REDEF (nout)
1226	if (ierr.ne.NF_NOERR) then
1227	print*,'missing_value: '
1228	print*, NF_STRERROR(ierr)
1229	endif
1230
1231	! Write valid_range() attribute
1232	if (validr.eq.NF_NOERR) then
1233	#ifdef NC_DOUBLE
1234	ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
1235	#else
1236	ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
1237	#endif
1238
1239	if (ierr.ne.NF_NOERR) then
1240	print*,'missing_value: valid_range attribution failed'
1241	print*, NF_STRERROR(ierr)
1242	!write(,) 'NF_NOERR', NF_NOERR
1243	!CALL abort
1244	stop ""
1245	endif
1246	endif
1247
1248	! Write "missing_value" attribute
1249	if (miss.eq.NF_NOERR) then
1250	#ifdef NC_DOUBLE
1251	ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
1252	#else
1253	ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
1254	#endif
1255
1256	if (ierr.NE.NF_NOERR) then
1257	print*, 'missing_value: missing value attribution failed'
1258	print*, NF_STRERROR(ierr)
1259	! WRITE(,) 'NF_NOERR', NF_NOERR
1260	! CALL abort
1261	stop ""
1262	endif
1263	endif
1264
1265	! End netcdf dataset define mode
1266	ierr = NF_ENDDEF(nout)
1267
1268	end subroutine missing_value
1269
1270	!*****************************************************************************
1271	subroutine interpolf(x,y,missing,xd,yd,nd)
1272	!==============================================================================
1273	! Purpose:
1274	! Yield y=f(x), where xd() end yd() are arrays of known values,
1275	! using linear interpolation
1276	! If x is not included in the interval spaned by xd(), then y is set
1277	! to a default value 'missing'
1278	! Note:
1279	! Array xd() should contain ordered (either increasing or decreasing) abscissas
1280	!==============================================================================
1281	implicit none
1282	!==============================================================================
1283	! Arguments:
1284	!==============================================================================
1285	real,intent(in) :: x ! abscissa at which interpolation is to be done
1286	real,intent(in) :: missing ! missing value (if no interpolation is performed)
1287	integer :: nd ! size of arrays
1288	real,dimension(nd),intent(in) :: xd ! array of known absissas
1289	real,dimension(nd),intent(in) :: yd ! array of correponding values
1290
1291	real,intent(out) :: y ! interpolated value
1292	!==============================================================================
1293	! Local variables:
1294	!==============================================================================
1295	integer :: i
1296
1297	! default: set y to 'missing'
1298	y=missing
1299
1300	do i=1,nd-1
1301	if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1302	((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1303	if ((yd(i).eq.missing).or.(yd(i+1).eq.missing)) then
1304	! cannot perform the interpolation if an encompasing value
1305	! is already set to 'missing'
1306	else
1307	!linear interpolation based on encompassing values
1308	y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1309	endif
1310	exit
1311	endif
1312	enddo
1313
1314
1315	end subroutine interpolf
1316
1317	!******************************************************************************
1318
1319	end program localtime
1320

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

Download in other formats:

Original Format