source: trunk/LMDZ.MARS/util/localtime.F90 @ 2434

Last change on this file since 2434 was 2434, checked in by abierjon, 4 years ago

Mars GCM:
Update utilities :
concat.F90

  • rewrite and simplify the handling of time and offset so that any file can be concatenated, including files from different years or stats file
  • use modulo to shift starting sols in concat.nc to values between 0 and 669
  • add a "adls" option in addition to "sol" or "ls" to add a Ls 1D variable

while keeping "sol" as the Time axis

  • add conservation of altitude attributes (long_name,units,positive)
  • enable absence of both hybrid (aps,bps) and sigma levels

solzenangle.F90

  • improve calculation for 1 sol file (stats) to use all local time data
  • read the starting sol in the file instead of asking it to the user (except for stats file) ; keep the possibility for user to change it
  • ask mean Ls value for stats file instead of sol number
  • fix crash when using "all variable" option (ticket #66)
  • fix bug on aps,bps by using GCM_layers dimension instead of altitude

localtime.F90

  • fix crash when using "all variable" option (ticket #66)
  • fix bug on aps,bps by using GCM_layers dimension instead of altitude

For all 3 utilities :

  • handle all ndim cases for the variables (ticket #52)
  • change "title" attribute into "long_name" by default (ticket #48)
  • extend size of long_name character string (ticket #57)
  • remove the use of #ifdef NC_DOUBLE (ticket #67)

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