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

Last change on this file since 3493 was 2577, checked in by emillour, 3 years ago

Mars GCM utilities:
Fixes in the utilities for the picky gfortran 10+ compiler.
JL+EM
File size: 43.9 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=100)  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=100), 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(*,'(a100)') tmpvar
202do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
203   nbvar=nbvar+1
204   var(nbvar)=tmpvar
205   read(*,'(a100)') 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(*,*) '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
745ierr=NF_CLOSE(nout)
746
747write(*,*) "Program completed !"
748
749contains
750
751!******************************************************************************
752Subroutine 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
764implicit none
765
766include "netcdf.inc" ! NetCDF definitions
767
768!==============================================================================
769! Arguments:
770!==============================================================================
771character (len=*), intent(in):: filename
772! filename(): the file's name
773integer,intent(in) ::latlen,lonlen,altlen,ctllen
774real, intent(in):: lat(latlen)
775! lat(): latitude
776real, intent(in):: lon(lonlen)
777! lon(): longitude
778real, intent(in):: alt(altlen)
779! alt(): altitude
780real, intent(in):: ctl(ctllen)
781! ctl(): controle
782integer,intent(in) :: GCM_layers ! number of GCM layers
783character (len=*), intent(in) :: altlong_name
784! altlong_name(): [netcdf] altitude "long_name" attribute
785character (len=*), intent(in) :: altunits
786! altunits(): [netcdf] altitude "units" attribute
787character (len=*), intent(in) :: altpositive
788! altpositive(): [netcdf] altitude "positive" attribute
789integer, intent(out):: nout
790! nout: [netcdf] file ID
791integer, intent(out):: latdimout
792! latdimout: [netcdf] lat() (i.e.: latitude)  ID
793integer, intent(out):: londimout
794! londimout: [netcdf] lon()  ID
795integer, intent(out):: altdimout
796! altdimout: [netcdf] alt()  ID
797integer, intent(out):: timedimout
798! timedimout: [netcdf] "Time"  ID
799integer,intent(out) :: layerdimout
800! layerdimout: [netcdf] "GCM_layers" ID
801integer, intent(out):: timevarout
802! timevarout: [netcdf] "Time" (considered as a variable) ID
803
804!==============================================================================
805! Local variables:
806!==============================================================================
807!integer :: latdim,londim,altdim,timedim
808integer :: nvarid,ierr
809integer :: 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!==============================================================================
816write(*,*) "creating "//trim(adjustl(filename))//'...'
817ierr = 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
819if (ierr.NE.NF_NOERR) then
820   WRITE(*,*)'ERROR: Impossible to create the file ',trim(filename)
821   write(*,*) NF_STRERROR(ierr)
822   stop
823endif
824
825!==============================================================================
826! 2. Define/write "dimensions" and get their IDs
827!==============================================================================
828
829ierr = NF_DEF_DIM(nout, "latitude", latlen, latdimout)
830if (ierr.NE.NF_NOERR) then
831   WRITE(*,*)'initiate: error failed to define dimension <latitude>'
832   write(*,*) NF_STRERROR(ierr)
833   stop
834endif
835ierr = NF_DEF_DIM(nout, "longitude", lonlen, londimout)
836if (ierr.NE.NF_NOERR) then
837   WRITE(*,*)'initiate: error failed to define dimension <longitude>'
838   write(*,*) NF_STRERROR(ierr)
839   stop
840endif
841ierr = NF_DEF_DIM(nout, "altitude", altlen, altdimout)
842if (ierr.NE.NF_NOERR) then
843   WRITE(*,*)'initiate: error failed to define dimension <altitude>'
844   write(*,*) NF_STRERROR(ierr)
845   stop
846endif
847if (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
854endif
855ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
856if (ierr.NE.NF_NOERR) then
857   WRITE(*,*)'initiate: error failed to define dimension <Time>'
858   write(*,*) NF_STRERROR(ierr)
859   stop
860endif
861
862ierr = NF_DEF_DIM(nout, "GCM_layers", GCM_layers, layerdimout)
863
864! End netcdf define mode
865ierr = NF_ENDDEF(nout)
866if (ierr.NE.NF_NOERR) then
867   WRITE(*,*)'initiate: error failed to switch out of define mode'
868   write(*,*) NF_STRERROR(ierr)
869   stop
870endif
871
872!==============================================================================
873! 3. Write "Time" and its attributes
874!==============================================================================
875
876call 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
883call def_var(nout,"latitude","latitude","degrees_north",1,&
884             (/latdimout/),nvarid,ierr)
885
886ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
887if (ierr.NE.NF_NOERR) then
888   WRITE(*,*)'initiate: error failed writing variable <latitude>'
889   write(*,*) NF_STRERROR(ierr)
890   stop
891endif
892
893!==============================================================================
894! 4. Write "longitude" (data and attributes)
895!==============================================================================
896
897call def_var(nout,"longitude","East longitude","degrees_east",1,&
898             (/londimout/),nvarid,ierr)
899
900
901ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
902if (ierr.NE.NF_NOERR) then
903   WRITE(*,*)'initiate: error failed writing variable <longitude>'
904   write(*,*) NF_STRERROR(ierr)
905   stop
906endif
907
908!==============================================================================
909! 5. Write "altitude" (data and attributes)
910!==============================================================================
911
912! Switch to netcdf define mode
913ierr = NF_REDEF (nout)
914
915ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,[altdimout],nvarid)
916
917ierr = NF_PUT_ATT_TEXT (nout,nvarid,'long_name',len_trim(adjustl(altlong_name)),adjustl(altlong_name))
918ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',len_trim(adjustl(altunits)),adjustl(altunits))
919ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',len_trim(adjustl(altpositive)),adjustl(altpositive))
920
921! End netcdf define mode
922ierr = NF_ENDDEF(nout)
923
924ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
925if (ierr.NE.NF_NOERR) then
926   WRITE(*,*)'initiate: error failed writing variable <altitude>'
927   write(*,*) NF_STRERROR(ierr)
928   stop
929endif
930
931!==============================================================================
932! 6. Write "controle" (data and attributes)
933!==============================================================================
934
935if (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
952endif
953
954end Subroutine initiate
955!******************************************************************************
956subroutine 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
967implicit none
968
969include "netcdf.inc" ! NetCDF definitions
970
971!==============================================================================
972! Arguments:
973!==============================================================================
974integer, intent(in) :: infid  ! NetCDF output file ID
975integer, intent(in) :: lonlen ! # of grid points along longitude
976integer, intent(in) :: latlen ! # of grid points along latitude
977integer, intent(in) :: altlen ! # of grid points along altitude
978integer, intent(in) :: GCM_layers ! # of GCM atmospheric layers
979integer, intent(in) :: outfid ! NetCDF output file ID
980integer, intent(in) :: londimout ! longitude dimension ID
981integer, intent(in) :: latdimout ! latitude dimension ID
982integer, intent(in) :: altdimout ! altitude dimension ID
983integer, intent(in) :: layerdimout ! GCM_layers dimension ID
984!==============================================================================
985! Local variables:
986!==============================================================================
987real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
988real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
989integer :: apsid,bpsid,phisinitid
990integer :: ierr
991integer :: tmpdimid ! temporary dimension ID
992integer :: tmpvarid ! temporary variable ID
993logical :: 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
1002allocate(aps(GCM_layers),stat=ierr)
1003if (ierr.ne.0) then
1004  write(*,*) "init2: failed to allocate aps!"
1005  stop
1006endif
1007ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
1008if (ierr.ne.NF_NOERR) then
1009  write(*,*) "Ooops. Failed to get aps ID. OK."
1010  hybrid=.false.
1011else
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
1038endif
1039
1040! ground geopotential phisinit
1041allocate(phisinit(lonlen,latlen),stat=ierr)
1042if (ierr.ne.0) then
1043  write(*,*) "init2: failed to allocate phisinit(lonlen,latlen)"
1044  stop
1045endif
1046ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
1047if (ierr.ne.NF_NOERR) then
1048  write(*,*) "Failed to get phisinit ID. OK"
1049  phisinit = 0.
1050  phis = .false.
1051else
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.
1058endif
1059
1060
1061!==============================================================================
1062! 2. Write
1063!==============================================================================
1064
1065!==============================================================================
1066! 2.2. Hybrid coordinates aps() and bps()
1067!==============================================================================
1068
1069IF (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)
1102ENDIF ! of IF (hybrid)
1103
1104!==============================================================================
1105! 2.2. phisinit()
1106!==============================================================================
1107
1108IF (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)
1126ENDIF ! of IF (phis)
1127
1128
1129end subroutine init2
1130
1131!******************************************************************************
1132subroutine 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
1142implicit none
1143
1144include "netcdf.inc" ! NetCDF definitions
1145
1146!==============================================================================
1147! Arguments:
1148!==============================================================================
1149integer, intent(in) :: nid
1150! nid: [netcdf] file ID #
1151character (len=*), intent(in) :: name
1152! name(): [netcdf] variable's name
1153character (len=*), intent(in) :: long_name
1154! long_name(): [netcdf] variable's "long_name" attribute
1155character (len=*), intent(in) :: units
1156! unit(): [netcdf] variable's "units" attribute
1157integer, intent(in) :: nbdim
1158! nbdim: number of dimensions of the variable
1159integer, dimension(nbdim), intent(in) :: dim
1160! dim(nbdim): [netcdf] dimension(s) ID(s)
1161integer, intent(out) :: nvarid
1162! nvarid: [netcdf] ID # of the variable
1163integer, intent(out) :: ierr
1164! ierr: [netcdf] subroutines returned error code
1165
1166! Switch to netcdf define mode
1167ierr=NF_REDEF(nid)
1168
1169! Insert the definition of the variable
1170ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
1171
1172! Write the attributes
1173ierr=NF_PUT_ATT_TEXT(nid,nvarid,"long_name",len_trim(adjustl(long_name)),adjustl(long_name))
1174ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
1175
1176! End netcdf define mode
1177ierr=NF_ENDDEF(nid)
1178
1179end subroutine def_var
1180!******************************************************************************
1181subroutine  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
1192implicit none
1193
1194include "netcdf.inc"  ! NetCDF definitions
1195
1196!==============================================================================
1197! Arguments:
1198!==============================================================================
1199integer, intent(in) :: nout
1200! nout: [netcdf] file ID #
1201integer, intent(in) :: nvarid
1202! varid: [netcdf] variable ID #
1203integer, intent(in) :: validr
1204! validr : [netcdf] routines return code for "valid_range" attribute
1205integer, intent(in) :: miss
1206! miss : [netcdf] routines return code for "missing_value" attribute
1207real, dimension(2), intent(in) :: valid_range
1208! valid_range(2): [netcdf] "valid_range" attribute (min and max)
1209real, intent(in) :: missing
1210! missing: [netcdf] "missing_value" attribute
1211
1212!==============================================================================
1213! Local variables:
1214!==============================================================================
1215integer :: 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
1222ierr = NF_REDEF (nout)
1223if (ierr.ne.NF_NOERR) then
1224   print*,'missing_value: '
1225   print*, NF_STRERROR(ierr)
1226endif
1227
1228! Write valid_range() attribute
1229if (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
1239endif
1240if (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
1251endif
1252
1253! End netcdf dataset define mode
1254ierr = NF_ENDDEF(nout)
1255
1256end subroutine  missing_value
1257
1258!*****************************************************************************
1259subroutine 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!==============================================================================
1269implicit none
1270!==============================================================================
1271! Arguments:
1272!==============================================================================
1273real,intent(in) :: x ! abscissa at which interpolation is to be done
1274real,intent(in) :: missing ! missing value (if no interpolation is performed)
1275integer :: nd ! size of arrays
1276real,dimension(nd),intent(in) :: xd ! array of known absissas
1277real,dimension(nd),intent(in) :: yd ! array of correponding values
1278
1279real,intent(out) :: y ! interpolated value
1280!==============================================================================
1281! Local variables:
1282!==============================================================================
1283integer :: i
1284
1285! default: set y to 'missing'
1286y=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
1303end subroutine interpolf
1304
1305!******************************************************************************
1306
1307end program localtime
1308
Note: See TracBrowser for help on using the repository browser.