Context Navigation

localtime.F90 @ 2567

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

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

Download in other formats:

Original Format