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

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

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