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

Last change on this file since 626 was 410, checked in by acolaitis, 13 years ago

Modified NCDF norm of our files from classic to 64bit offset to support variables indices of more than integer*4 maximum length. This is a priori retrocompatible with classic format. I have tested it, it works for all file outputs, newstart (on classic or 64-bit offset files). One can check the format of his .nc with ncdump -k file.nc
File size: 33.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=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
29integer :: nid,ierr,miss
30! nid: [netcdf] file ID #
31! ierr: [netcdf] subroutine returned error code
32! miss: [netcdf] subroutine returned error code
33integer :: i,j,k,inter
34! for various loops
35integer :: varid
36! varid: [netcdf] variable ID #
37real, dimension(:), allocatable:: lat,lon,alt,time
38! lat(): array, stores latitude coordinates
39! lon(): array, stores longitude coordinates
40! alt(): array, stores altitude coordinates
41! time(): array, stores time coordinates
42integer :: nbvar,nbvarfile,ndim
43! nbvar: # of variables to concatenate
44! nbfile: # number of input file(s)
45! nbvarfile: total # of variables in an input file
46! ndim: [netcdf] # (3 or 4) of dimensions (for variables)
47integer :: latdim,londim,altdim,timedim
48! latdim: [netcdf] "latitude" dim ID
49! londim: [netcdf] "longitude" dim ID
50! altdim: [netcdf] "altdim" dim ID
51! timedim: [netcdf] "timedim" dim ID
52integer :: latvar,lonvar,altvar,timevar
53! latvar: [netcdf] ID of "latitude" variable
54! lonvar: [netcdf] ID of "longitude" variable
55! altvar: [netcdf] ID of "altitude" variable
56! timevar: [netcdf] ID of "Time" variable
57integer :: latlen,lonlen,altlen,timelen,timelen_lt,timelen_tot
58integer :: ilat,ilon,ialt,it
59! latlen: # of elements of lat() array
60! lonlen: # of elements of lon() array
61! altvar: # of elements of alt() array
62! timelen: # of elemnets of time() array
63! timelen_tot:# =timelen or timelen+1 (if 1 more time to interpolate needed)
64! timelen_lt: # of elemnets of time() array in output
65integer :: nout,latdimout,londimout,altdimout,timedimout,timevarout
66integer :: nhour,nsol
67! nout: [netcdf] output file ID
68! latdimout: [netcdf] output latitude (dimension) ID
69! londimout: [netcdf] output longitude (dimension) ID
70! altdimout: [netcdf] output altitude (dimension) ID
71! timedimout: [netcdf] output time (dimension) ID
72! timevarout: [netcdf] ID of output "Time" variable
73integer :: varidout
74! varidout: [netcdf] variable ID # (of a variable to write to the output file)
75integer :: Nnotconcat,var_ok
76! Nnotconcat: # of (leading)variables that won't be concatenated
77! var_ok: flag (0 or 1)
78integer, dimension(4) :: corner,edges,dim
79! corner: [netcdf]
80! edges: [netcdf]
81! dim: [netcdf]
82real, dimension(:,:,:,:), allocatable :: var3d
83! var3D(,,,): 4D array to store a field
84
85real, dimension(:,:,:,:), allocatable :: var3d_lt
86! var3D_lt(,,,): 4D array to store a field in local time coordinate
87real, dimension(:), allocatable :: lt_gcm,lt_hour
88real, dimension(:), allocatable :: lt_out,lt_outc
89real, dimension(:), allocatable :: var_gcm
90
91real :: missing
92!PARAMETER(missing=1E+20)
93! missing: [netcdf] to handle "missing" values when reading/writing files
94real, dimension(2) :: valid_range
95! valid_range(2): [netcdf] interval in which a value is considered valid
96logical :: stats   ! stats=T when reading a "stats" kind of ile
97
98
99!==============================================================================
100! 1.1. Get input file name(s)
101!==============================================================================
102write(*,*)
103write(*,*) "which file do you want to use?  (diagfi... stats...  concat...)"
104
105read(*,'(a50)') file
106
107if (len_trim(file).eq.0) then
108   write(*,*) "no file... game over"
109   stop ""
110endif
111
112!==============================================================================
113! 1.3. Open the first input file
114!==============================================================================
115
116ierr = NF_OPEN(file,NF_NOWRITE,nid)
117if (ierr.NE.NF_NOERR) then
118   write(*,*) 'ERROR: Pb opening file '//file
119   stop ""
120endif
121
122ierr=NF_INQ_NVARS(nid,nbvarfile)
123! nbvarfile now set to be the (total) number of variables in file
124
125!==============================================================================
126! 1.4. Ask for (output) "Time" axis type
127!==============================================================================
128
129notconcat(1)='Time'
130notconcat(2)='controle'
131notconcat(3)='rlonu'
132notconcat(4)='latitude'
133notconcat(5)='longitude'
134notconcat(6)='altitude'
135notconcat(7)='rlatv'
136notconcat(8)='aps'
137notconcat(9)='bps'
138notconcat(10)='ap'
139notconcat(11)='bp'
140notconcat(12)='cu'
141notconcat(13)='cv'
142notconcat(14)='aire'
143notconcat(15)='phisinit'
144
145!==============================================================================
146! 1.5. Get (and check) list of variables to concatenate
147!==============================================================================
148write(*,*)
149   Nnotconcat=0
150do i=1,nbvarfile
151   ierr=NF_INQ_VARNAME(nid,i,vartmp)
152   ! vartmp now contains the "name" of variable of ID # i
153   var_ok=0
154   do inter=1,15
155      if (vartmp.eq.notconcat(inter)) then
156         var_ok=1
157         Nnotconcat=Nnotconcat+1
158      endif
159   enddo       
160   if (var_ok.eq.0)  write(*,*) trim(vartmp)
161enddo
162
163! Nnotconcat: # of variables that won't be concatenated
164! nbvarfile: total # of variables in file
165allocate(var(nbvarfile-Nnotconcat))
166
167
168write(*,*)
169write(*,*) "which variables do you want to redistribute ?"
170write(*,*) "all / list of <variables> (separated by <Enter>s)"
171write(*,*) "(an empty line , i.e: just <Enter>, implies end of list)"
172nbvar=0
173read(*,'(a50)') tmpvar
174do while ((tmpvar/=' ').AND.(trim(tmpvar)/='all'))
175   nbvar=nbvar+1
176   var(nbvar)=tmpvar
177   read(*,'(a50)') tmpvar
178enddo
179
180if (tmpvar=="all") then
181         nbvar=nbvarfile-Nnotconcat
182         do j=Nnotconcat+1,nbvarfile
183            ierr=nf_inq_varname(nid,j,var(j-Nnotconcat))
184         enddo
185! Variables names from the file are catched
186   nbvar=nbvarfile-Nnotconcat
187   do i=1,nbvar
188      ierr=nf_inq_varname(nid,i+Nnotconcat,var(i))
189      write(*,'(a9,1x,i2,1x,a1,1x,a50)') "variable ",i,":",var(i)
190   enddo
191else if(nbvar==0) then
192   write(*,*) "no variable... game over"
193   stop ""
194endif ! of if (tmpvar=="all")
195
196!==============================================================================
197! 1.6. Get output file name
198!==============================================================================
199 filename=file(1:len_trim(file)-3)//"_LT.nc"
200 write(*,*) filename
201
202!==============================================================================
203! 2.1. Open input file
204!==============================================================================
205
206   if (i/=1) then
207      write(*,*)
208      write(*,*) "opening "//trim(file)//"..."
209      ierr = NF_OPEN(file,NF_NOWRITE,nid)
210      if (ierr.NE.NF_NOERR) then
211         write(*,*) 'ERROR: Pb opening file '//file
212         write(*,*) NF_STRERROR(ierr)
213         stop ""
214      endif
215   endif
216
217!==============================================================================
218! 2.2. Read (and check) dimensions of variables from input file
219!==============================================================================
220
221   ierr=NF_INQ_DIMID(nid,"latitude",latdim)
222   ierr=NF_INQ_VARID(nid,"latitude",latvar)
223   if (ierr.NE.NF_NOERR) then
224      write(*,*) 'ERROR: Field <latitude> is missing in file'//file
225      stop "" 
226   endif
227   ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
228!  write(*,*) "latlen: ",latlen
229
230   ierr=NF_INQ_DIMID(nid,"longitude",londim)
231   ierr=NF_INQ_VARID(nid,"longitude",lonvar)
232   if (ierr.NE.NF_NOERR) then
233      write(*,*) 'ERROR: Field <longitude> is missing in file'//file
234      stop ""
235   endif
236   ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
237!  write(*,*) "lonlen: ",lonlen
238
239   ierr=NF_INQ_DIMID(nid,"altitude",altdim)
240   ierr=NF_INQ_VARID(nid,"altitude",altvar)
241   if (ierr.NE.NF_NOERR) then
242      write(*,*) 'ERROR: Field <altitude> is missing in file'//file
243      stop ""
244   endif
245   ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
246!  write(*,*) "altlen: ",altlen
247
248!==============================================================================
249! 2.3. Read (and check compatibility of) dimensions of
250!       variables from input file
251!==============================================================================
252
253! First call; initialize/allocate
254      allocate(lat(latlen))
255      allocate(lon(lonlen))
256      allocate(alt(altlen))
257#ifdef NC_DOUBLE
258      ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
259      ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
260      ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
261#else
262      ierr = NF_GET_VAR_REAL(nid,latvar,lat)
263      ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
264      ierr = NF_GET_VAR_REAL(nid,altvar,alt)
265#endif
266!==============================================================================
267! 2.4. Handle "Time" dimension from input file
268!==============================================================================
269
270!==============================================================================
271! 2.4.0 Read "Time" dimension from input file
272!==============================================================================
273   ierr=NF_INQ_DIMID(nid,"Time",timedim)
274   if (ierr.NE.NF_NOERR) then
275      write(*,*) 'ERROR: Dimension <Time> is missing in file'//file
276      stop ""
277   endif
278   ierr=NF_INQ_VARID(nid,"Time",timevar)
279   if (ierr.NE.NF_NOERR) then
280      write(*,*) 'ERROR: Field <Time> is missing in file'//file
281      stop ""
282   endif
283   ierr=NF_INQ_DIMLEN(nid,timedim,timelen)
284!  write(*,*) "timelen: ",timelen
285
286   ! allocate time() array and fill it with values from input file
287   allocate(time(timelen))
288
289#ifdef NC_DOUBLE
290   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
291#else
292   ierr = NF_GET_VAR_REAL(nid,timevar,time)
293#endif
294
295
296!==============================================================================
297! 2.4.1 Select local times to be saved "Time" in output file
298!==============================================================================
299      write(*,*) 'number of local time to be used ?'
300      read(*,*) nhour
301      allocate(lt_hour(nhour))
302      write(*,*) 'list of selected local time (0<t<24)'
303      do it=1,nhour
304         read(*,*) lt_hour(it)
305      end do
306
307      if ((timelen.eq.12).and.(time(1).eq.2).and.(time(timelen).eq.24))then
308         write(*,*) 'We detect a ""stats"" file'
309         stats=.true.
310         timelen_lt=nhour
311         allocate(lt_out(timelen_lt))
312         do it=1,nhour
313           lt_out(it)=lt_hour(it)/24.
314         end do
315!        We rewrite the time from "stats" from 0 to 1 sol...
316         do it=1,timelen  ! loop temps in
317               time(it) = time(it)/24.
318         end do
319         nsol =1
320      else   ! case of a diagfi or concatnc file
321        stats=.false.
322!       Number of sol in input file
323        nsol = int(time(timelen)) - int(time(1))
324        timelen_lt=nhour*nsol
325        allocate(lt_out(timelen_lt))
326        i=0
327        do k=1,nsol
328          do it=1,nhour
329            i=i+1
330            lt_out(i)=int(time(1)) + k-1  + lt_hour(it)/24.
331          end do
332        end do
333        end if
334
335      if (nsol.gt.1) then
336         timelen_tot=timelen
337      else
338!        if only 1 sol available, we must add 1 timestep for the interpolation
339         timelen_tot=timelen+1
340      endif     
341      allocate(lt_gcm(timelen_tot))
342      allocate(var_gcm(timelen_tot))
343
344      allocate(lt_outc(timelen_lt))
345         
346
347!==============================================================================
348! 2.4.1.5 Initiate dimension in output file
349!==============================================================================
350
351
352   ! Initialize output file's lat,lon,alt and time dimensions
353      call initiate (filename,lat,lon,alt,nout,&
354           latdimout,londimout,altdimout,timedimout,timevarout)
355   ! Initialize output file's aps,bps and phisinit variables
356     call init2(nid,lonlen,latlen,altlen,&
357                nout,londimout,latdimout,altdimout)
358
359
360!==============================================================================
361! 2.4.2 Write/extend "Time" dimension/values in output file
362!==============================================================================
363
364   if (stats) then
365
366     do it=1,timelen_lt
367#ifdef NC_DOUBLE
368        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_hour(it))
369#else
370        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_hour(it))
371#endif
372     enddo
373   else
374     do it=1,timelen_lt
375#ifdef NC_DOUBLE
376        ierr= NF_PUT_VARA_DOUBLE(nout,timevarout,it,1,lt_out(it))
377#else
378        ierr= NF_PUT_VARA_REAL(nout,timevarout,it,1,lt_out(it))
379#endif
380     enddo
381  end if
382
383!==============================================================================
384! 2.5 Read/write variables
385!==============================================================================
386
387   do j=1,nbvar ! loop on variables to read/write
388
389!==============================================================================
390! 2.5.1 Check that variable to be read existe in input file
391!==============================================================================
392
393      write(*,*) "variable ",var(j)
394      ierr=nf_inq_varid(nid,var(j),varid)
395      if (ierr.NE.NF_NOERR) then
396         write(*,*) 'ERROR: Field <',var(j),'> not found in file'//file
397         stop ""
398      endif
399      ierr=nf_inq_varndims(nid,varid,ndim)
400
401!==============================================================================
402! 2.5.2 Prepare things in order to read/write the variable
403!==============================================================================
404
405      ! build dim(),corner() and edges() arrays
406      ! and allocate var3d() array
407      if (ndim==3) then
408         allocate(var3d(lonlen,latlen,timelen,1))
409         allocate(var3d_lt(lonlen,latlen,timelen_lt,1))
410         dim(1)=londimout
411         dim(2)=latdimout
412         dim(3)=timedimout
413
414         ! start indexes (where data values will be written)
415         corner(1)=1
416         corner(2)=1
417         corner(3)=1
418         corner(4)=1
419
420         ! length (along dimensions) of block of data to be written
421         edges(1)=lonlen
422         edges(2)=latlen
423         edges(3)=timelen_lt
424         edges(4)=1
425   
426      else if (ndim==4) then
427         allocate(var3d(lonlen,latlen,altlen,timelen))
428         allocate(var3d_lt(lonlen,latlen,altlen,timelen_lt))
429         dim(1)=londimout
430         dim(2)=latdimout
431         dim(3)=altdimout
432         dim(4)=timedimout
433
434         ! start indexes (where data values will be written)
435         corner(1)=1
436         corner(2)=1
437         corner(3)=1
438         corner(4)=1
439
440         ! length (along dimensions) of block of data to be written
441         edges(1)=lonlen
442         edges(2)=latlen
443         edges(3)=altlen
444         edges(4)=timelen_lt
445      endif
446
447         units="                                                    "
448         title="                                                    "
449         ierr=nf_get_att_text(nid,varid,"title",title)
450         ierr=nf_get_att_text(nid,varid,"units",units)
451         call def_var(nout,var(j),title,units,ndim,dim,varidout,ierr)
452
453     
454
455!==============================================================================
456! 2.5.3 Read from input file
457!==============================================================================
458
459#ifdef NC_DOUBLE
460      ierr = NF_GET_VAR_DOUBLE(nid,varid,var3d)
461      miss=NF_GET_ATT_DOUBLE(nid,varid,"missing_value",missing)
462      miss=NF_GET_ATT_DOUBLE(nid,varid,"valid_range",valid_range)
463#else
464      ierr = NF_GET_VAR_REAL(nid,varid,var3d)
465      miss=NF_GET_ATT_REAL(nid,varid,"missing_value",missing)
466      miss=NF_GET_ATT_REAL(nid,varid,"valid_range",valid_range)
467#endif
468
469!==============================================================================
470! 2.5.4 Read from input file
471!==============================================================================
472! Interpolation in local time :
473
474        do ilon=1,lonlen
475!          write(*,*) 'processing longitude =', lon(ilon)
476!          Local time at each longitude :
477           do it=1,timelen  ! loop temps in
478               lt_gcm(it) = time(it) +0.5*lon(ilon)/180.
479           end do
480           if (nsol.eq.1) lt_gcm(timelen+1) = lt_gcm(1)+1
481
482           do it=1,timelen_lt ! loop time local time
483              lt_outc(it)=lt_out(it)
484!             Correction to use same local time previous or following
485!             day where data are missing :
486
487              if(lt_outc(it).gt.lt_gcm(timelen_tot))lt_outc(it)=lt_out(it)-1
488              if(lt_outc(it).lt.lt_gcm(1))lt_outc(it)=lt_out(it)+1
489           end do
490
491           if (ndim==3) then
492             do ilat=1,latlen
493               do it=1,timelen  ! loop temps in
494                  var_gcm(it) = var3d(ilon,ilat,it,1)
495               end do
496               if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
497               do it=1,timelen_lt ! loop time local time
498                  call interpolf(lt_outc(it),var3d_lt(ilon,ilat,it,1),&
499                                 missing,lt_gcm,var_gcm,timelen_tot)
500               end do   
501             enddo
502
503           else if  (ndim==4) then
504             do ilat=1,latlen
505               do ialt=1,altlen
506                 do it=1,timelen ! loop temps in
507                    var_gcm(it) = var3d(ilon,ilat,ialt,it)
508                 end do
509                 if (nsol.eq.1) var_gcm(timelen+1) = var_gcm(1)
510                 do it=1,timelen_lt ! loop time local time
511                    call interpolf(lt_outc(it),var3d_lt(ilon,ilat,ialt,it),&
512                                   missing,lt_gcm,var_gcm,timelen_tot)
513                 end do   
514               enddo
515             enddo
516           end if
517
518        end do
519
520
521!==============================================================================
522! 2.5.5 write (append) to the output file
523!==============================================================================
524
525#ifdef NC_DOUBLE
526      ierr= NF_PUT_VARA_DOUBLE(nout,varidout,corner,edges,var3d_lt)
527#else
528      ierr= NF_PUT_VARA_REAL(nout,varidout,corner,edges,var3d_lt)
529#endif
530
531      if (ierr.ne.NF_NOERR) then
532         write(*,*) 'PUT_VAR ERROR: ',NF_STRERROR(ierr)
533         stop ""
534      endif
535
536! In case there is a "valid_range" attribute
537      ! Write "valid_range" and "missing_value" attributes in output file
538      if (miss.eq.NF_NOERR) then
539         call missing_value(nout,varidout,valid_range,missing)
540      endif
541
542      ! free var3d() array
543      deallocate(var3d)
544      deallocate(var3d_lt)
545
546   enddo ! of do j=1,nbvar
547
548   deallocate(time)
549   deallocate(lt_gcm)
550   deallocate(lt_out)
551   deallocate(lt_outc)
552   deallocate(var_gcm)
553
554   ! Close input file
555   ierr=nf_close(nid)
556
557! Close output file
558ierr=NF_CLOSE(nout)
559
560contains
561
562!******************************************************************************
563Subroutine initiate (filename,lat,lon,alt,&
564                     nout,latdimout,londimout,altdimout,timedimout,timevarout)
565!==============================================================================
566! Purpose:
567! Create and initialize a data file (NetCDF format)
568!==============================================================================
569! Remarks:
570! The NetCDF file (created in this subroutine) remains open
571!==============================================================================
572
573implicit none
574
575include "netcdf.inc" ! NetCDF definitions
576
577!==============================================================================
578! Arguments:
579!==============================================================================
580character (len=*), intent(in):: filename
581! filename(): the file's name
582real, dimension(:), intent(in):: lat
583! lat(): latitude
584real, dimension(:), intent(in):: lon
585! lon(): longitude
586real, dimension(:), intent(in):: alt
587! alt(): altitude
588integer, intent(out):: nout
589! nout: [netcdf] file ID
590integer, intent(out):: latdimout
591! latdimout: [netcdf] lat() (i.e.: latitude)  ID
592integer, intent(out):: londimout
593! londimout: [netcdf] lon()  ID
594integer, intent(out):: altdimout
595! altdimout: [netcdf] alt()  ID
596integer, intent(out):: timedimout
597! timedimout: [netcdf] "Time"  ID
598integer, intent(out):: timevarout
599! timevarout: [netcdf] "Time" (considered as a variable) ID
600
601!==============================================================================
602! Local variables:
603!==============================================================================
604!integer :: latdim,londim,altdim,timedim
605integer :: nvarid,ierr
606! nvarid: [netcdf] ID of a variable
607! ierr: [netcdf]  return error code (from called subroutines)
608
609!==============================================================================
610! 1. Create (and open) output file
611!==============================================================================
612write(*,*) "creating "//trim(adjustl(filename))//'...'
613ierr = NF_CREATE(filename,IOR(NF_CLOBBER,NF_64BIT_OFFSET),nout)
614! NB: setting NF_CLOBBER mode means that it's OK to overwrite an existing file
615if (ierr.NE.NF_NOERR) then
616   WRITE(*,*)'ERROR: Impossible to create the file.'
617   stop ""
618endif
619
620!==============================================================================
621! 2. Define/write "dimensions" and get their IDs
622!==============================================================================
623
624ierr = NF_DEF_DIM(nout, "latitude", size(lat), latdimout)
625ierr = NF_DEF_DIM(nout, "longitude", size(lon), londimout)
626ierr = NF_DEF_DIM(nout, "altitude", size(alt), altdimout)
627ierr = NF_DEF_DIM(nout, "Time", NF_UNLIMITED, timedimout)
628
629! End netcdf define mode
630ierr = NF_ENDDEF(nout)
631
632!==============================================================================
633! 3. Write "Time" and its attributes
634!==============================================================================
635
636call def_var(nout,"Time","Time","years since 0000-00-0 00:00:00",1,&
637             (/timedimout/),timevarout,ierr)
638
639!==============================================================================
640! 4. Write "latitude" (data and attributes)
641!==============================================================================
642
643call def_var(nout,"latitude","latitude","degrees_north",1,&
644             (/latdimout/),nvarid,ierr)
645
646#ifdef NC_DOUBLE
647ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
648#else
649ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
650#endif
651
652!==============================================================================
653! 4. Write "longitude" (data and attributes)
654!==============================================================================
655
656call def_var(nout,"longitude","East longitude","degrees_east",1,&
657             (/londimout/),nvarid,ierr)
658
659#ifdef NC_DOUBLE
660ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
661#else
662ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
663#endif
664
665!==============================================================================
666! 4. Write "altitude" (data and attributes)
667!==============================================================================
668
669! Switch to netcdf define mode
670ierr = NF_REDEF (nout)
671
672#ifdef NC_DOUBLE
673ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,altdimout,nvarid)
674#else
675ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,altdimout,nvarid)
676#endif
677
678ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",8,"altitude")
679ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"km")
680ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
681
682! End netcdf define mode
683ierr = NF_ENDDEF(nout)
684
685#ifdef NC_DOUBLE
686ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
687#else
688ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
689#endif
690
691end Subroutine initiate
692!******************************************************************************
693subroutine init2(infid,lonlen,latlen,altlen, &
694                 outfid,londimout,latdimout,altdimout)
695!==============================================================================
696! Purpose:
697! Copy aps(), bps() and phisinit() from input file to outpout file
698!==============================================================================
699! Remarks:
700! The NetCDF files must be open
701!==============================================================================
702
703implicit none
704
705include "netcdf.inc" ! NetCDF definitions
706
707!==============================================================================
708! Arguments:
709!==============================================================================
710integer, intent(in) :: infid  ! NetCDF output file ID
711integer, intent(in) :: lonlen ! # of grid points along longitude
712integer, intent(in) :: latlen ! # of grid points along latitude
713integer, intent(in) :: altlen ! # of grid points along latitude
714integer, intent(in) :: outfid ! NetCDF output file ID
715integer, intent(in) :: londimout ! longitude dimension ID
716integer, intent(in) :: latdimout ! latitude dimension ID
717integer, intent(in) :: altdimout ! altitude dimension ID
718!==============================================================================
719! Local variables:
720!==============================================================================
721real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
722real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
723integer :: apsid,bpsid,phisinitid
724integer :: ierr
725integer :: tmpvarid ! temporary variable ID
726logical :: phis, aps_ok, bps_ok ! are "phisinit" "aps" "bps" available ?
727
728
729!==============================================================================
730! 1. Read data from input file
731!==============================================================================
732
733! hybrid coordinate aps
734  allocate(aps(altlen))
735ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
736if (ierr.ne.NF_NOERR) then
737  write(*,*) "oops Failed to get aps ID. OK"
738  aps_ok=.false.
739else
740  ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
741  if (ierr.ne.NF_NOERR) then
742   stop "error: Failed reading aps"
743  endif
744  aps_ok=.true.
745endif
746
747! hybrid coordinate bps
748  allocate(bps(altlen))
749ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
750if (ierr.ne.NF_NOERR) then
751  write(*,*) "oops: Failed to get bps ID. OK"
752  bps_ok=.false.
753else
754  ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
755  if (ierr.ne.NF_NOERR) then
756    stop "Error: Failed reading bps"
757  endif
758  bps_ok=.true.
759endif
760
761! ground geopotential phisinit
762ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
763allocate(phisinit(lonlen,latlen))
764if (ierr.ne.NF_NOERR) then
765  write(*,*) "Failed to get phisinit ID. OK"
766  phisinit = 0.
767  phis = .false.
768else
769  ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
770  if (ierr.ne.NF_NOERR) then
771    stop "Error: Failed reading phisinit"
772  endif
773  phis = .true.
774endif
775
776
777!==============================================================================
778! 2. Write
779!==============================================================================
780
781!==============================================================================
782! 2.2. Hybrid coordinates aps() and bps()
783!==============================================================================
784
785IF (aps_ok) then
786! define aps
787call def_var(nout,"aps","hybrid pressure at midlayers"," ",1,&
788             (/altdimout/),apsid,ierr)
789if (ierr.ne.NF_NOERR) then
790  stop "Error: Failed to def_var aps"
791endif
792
793! write aps
794#ifdef NC_DOUBLE
795ierr=NF_PUT_VAR_DOUBLE(outfid,apsid,aps)
796#else
797ierr=NF_PUT_VAR_REAL(outfid,apsid,aps)
798#endif
799if (ierr.ne.NF_NOERR) then
800  stop "Error: Failed to write aps"
801endif
802ENDIF
803
804IF (bps_ok) then
805! define bps
806call def_var(nout,"bps","hybrid sigma at midlayers"," ",1,&
807             (/altdimout/),bpsid,ierr)
808if (ierr.ne.NF_NOERR) then
809  stop "Error: Failed to def_var bps"
810endif
811
812! write bps
813#ifdef NC_DOUBLE
814ierr=NF_PUT_VAR_DOUBLE(outfid,bpsid,bps)
815#else
816ierr=NF_PUT_VAR_REAL(outfid,bpsid,bps)
817#endif
818if (ierr.ne.NF_NOERR) then
819  stop "Error: Failed to write bps"
820endif
821ENDIF
822
823!==============================================================================
824! 2.2. phisinit()
825!==============================================================================
826
827
828IF (phis) THEN
829
830!define phisinit
831 call def_var(nout,"phisinit","Ground level geopotential"," ",2,&
832            (/londimout,latdimout/),phisinitid,ierr)
833 if (ierr.ne.NF_NOERR) then
834     stop "Error: Failed to def_var phisinit"
835  endif
836
837! write phisinit
838#ifdef NC_DOUBLE
839ierr=NF_PUT_VAR_DOUBLE(outfid,phisinitid,phisinit)
840#else
841ierr=NF_PUT_VAR_REAL(outfid,phisinitid,phisinit)
842#endif
843if (ierr.ne.NF_NOERR) then
844  stop "Error: Failed to write phisinit"
845endif
846
847END IF
848
849
850! Cleanup
851deallocate(aps)
852deallocate(bps)
853deallocate(phisinit)
854
855end subroutine init2
856!******************************************************************************
857subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
858!==============================================================================
859! Purpose: Write a variable (i.e: add a variable to a dataset)
860! called "name"; along with its attributes "title", "units"...
861! to a file (following the NetCDF format)
862!==============================================================================
863! Remarks:
864! The NetCDF file must be open
865!==============================================================================
866
867implicit none
868
869include "netcdf.inc" ! NetCDF definitions
870
871!==============================================================================
872! Arguments:
873!==============================================================================
874integer, intent(in) :: nid
875! nid: [netcdf] file ID #
876character (len=*), intent(in) :: name
877! name(): [netcdf] variable's name
878character (len=*), intent(in) :: title
879! title(): [netcdf] variable's "title" attribute
880character (len=*), intent(in) :: units
881! unit(): [netcdf] variable's "units" attribute
882integer, intent(in) :: nbdim
883! nbdim: number of dimensions of the variable
884integer, dimension(nbdim), intent(in) :: dim
885! dim(nbdim): [netcdf] dimension(s) ID(s)
886integer, intent(out) :: nvarid
887! nvarid: [netcdf] ID # of the variable
888integer, intent(out) :: ierr
889! ierr: [netcdf] subroutines returned error code
890
891! Switch to netcdf define mode
892ierr=NF_REDEF(nid)
893
894! Insert the definition of the variable
895#ifdef NC_DOUBLE
896ierr=NF_DEF_VAR(nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
897#else
898ierr=NF_DEF_VAR(nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
899#endif
900
901! Write the attributes
902ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",len_trim(adjustl(title)),adjustl(title))
903ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",len_trim(adjustl(units)),adjustl(units))
904
905! End netcdf define mode
906ierr=NF_ENDDEF(nid)
907
908end subroutine def_var
909!******************************************************************************
910subroutine  missing_value(nout,nvarid,valid_range,missing)
911!==============================================================================
912! Purpose:
913! Write "valid_range" and "missing_value" attributes (of a given
914! variable) to a netcdf file
915!==============================================================================
916! Remarks:
917! NetCDF file must be open
918! Variable (nvarid) ID must be set
919!==============================================================================
920
921implicit none
922
923include "netcdf.inc"  ! NetCDF definitions
924
925!==============================================================================
926! Arguments:
927!==============================================================================
928integer, intent(in) :: nout
929! nout: [netcdf] file ID #
930integer, intent(in) :: nvarid
931! varid: [netcdf] variable ID #
932real, dimension(2), intent(in) :: valid_range
933! valid_range(2): [netcdf] "valid_range" attribute (min and max)
934real, intent(in) :: missing
935! missing: [netcdf] "missing_value" attribute
936
937!==============================================================================
938! Local variables:
939!==============================================================================
940integer :: ierr
941! ierr: [netcdf] subroutine returned error code
942!      INTEGER nout,nvarid,ierr
943!      REAL missing
944!      REAL valid_range(2)
945
946! Switch to netcdf dataset define mode
947ierr = NF_REDEF (nout)
948if (ierr.ne.NF_NOERR) then
949   print*,'missing_value: '
950   print*, NF_STRERROR(ierr)
951endif
952
953! Write valid_range() attribute
954#ifdef NC_DOUBLE
955ierr=NF_PUT_ATT_DOUBLE(nout,nvarid,'valid_range',NF_DOUBLE,2,valid_range)
956#else
957ierr=NF_PUT_ATT_REAL(nout,nvarid,'valid_range',NF_FLOAT,2,valid_range)
958#endif
959
960if (ierr.ne.NF_NOERR) then
961   print*,'missing_value: valid_range attribution failed'
962   print*, NF_STRERROR(ierr)
963   !write(*,*) 'NF_NOERR', NF_NOERR
964   !CALL abort
965   stop ""
966endif
967
968! Write "missing_value" attribute
969#ifdef NC_DOUBLE
970ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
971#else
972ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
973#endif
974
975if (ierr.NE.NF_NOERR) then
976   print*, 'missing_value: missing value attribution failed'
977   print*, NF_STRERROR(ierr)
978!    WRITE(*,*) 'NF_NOERR', NF_NOERR
979!    CALL abort
980   stop ""
981endif
982
983! End netcdf dataset define mode
984ierr = NF_ENDDEF(nout)
985
986end subroutine  missing_value
987!******************************************************************************
988
989end program localtime
990
991!*****************************************************************************
992subroutine interpolf(x,y,missing,xd,yd,nd)
993!==============================================================================
994! Purpose:
995! Yield y=f(x), where xd() end yd() are arrays of known values,
996! using linear interpolation
997! If x is not included in the interval spaned by xd(), then y is set
998! to a default value 'missing'
999! Note:
1000! Array xd() should contain ordered (either increasing or decreasing) abscissas
1001!==============================================================================
1002implicit none
1003!==============================================================================
1004! Arguments:
1005!==============================================================================
1006real,intent(in) :: x ! abscissa at which interpolation is to be done
1007real,intent(in) :: missing ! missing value (if no interpolation is performed)
1008integer :: nd ! size of arrays
1009real,dimension(nd),intent(in) :: xd ! array of known absissas
1010real,dimension(nd),intent(in) :: yd ! array of correponding values
1011
1012real,intent(out) :: y ! interpolated value
1013!==============================================================================
1014! Local variables:
1015!==============================================================================
1016integer :: i
1017
1018! default: set y to 'missing'
1019y=missing
1020
1021   do i=1,nd-1
1022     if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
1023          ((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
1024        y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1025        exit
1026     endif
1027   enddo
1028
1029
1030end subroutine interpolf
1031
Note: See TracBrowser for help on using the repository browser.