WebSVN – lagranto.ecmwf – Diff – /trunk/lib/iotra.f


c     ****************************************************************
c     * This package provides IO routines for trajectories. A file   *
c     * is characterised by the filename <filename> and the file     *
c     * identifier <fid>. Different modes <mode> are supported:      *
c     *     mode=1: ascii, sorted by trajectory;                     *
c     *     mode=2: ascii, sorted by time;                           *
c     *     mode=3: fortran (unformatted)                            * 
c     *     mode=4: IVE netcdf (for compatibiltzy reasons)           *
c     *     mode=5: KML                                              *
c     *     mode=6: as mode 1, but with pressure as realk number     *
c     * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)>  *
c     * where <ntra> is the number of trajectories, <ntim> the       *
c     * number of times of each trajectory and <ncol> the number of  *
c     * columns of the trajectory. The first 4 columns are: time,    *
c     * longitude, latitude, pressure. The other columns are traced  *
c     * fields. The complete list of all columns is given in the     *
c     * array <vars(ncol)>. Finally, the reference date is given in  *
c     * the array <time(6)=year,month,day,hour,time length of the    *
c     * trajectory (hour,min)>.                                      *
c     *                                                              *
c     * Author: Michael Sprenger, September 2008                     *
c     ****************************************************************
 
c     ----------------------------------------------------------------
c     Open a trajectory file for reading
c     ----------------------------------------------------------------
      
      subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer      fid
      character*80 filename
      integer      mode
      integer      ntra,ntim,ncol
      integer      time(6)
      character*80 vars(ncol)
 
c     Auxiliary variables
      integer      vardim(4)
      real         varmin(4),varmax(4),stag(4)
      real         mdv
      character*80 cfn
      integer      ierr
      integer      i
      integer      nvars
 
c     Open file
      if (mode.eq.1) then
         fid = 10
         open(fid,file=filename)
      elseif (mode.eq.2) then
         fid = 10
         open(fid,file=filename)
      elseif (mode.eq.3) then
         open(fid,file=filename,form='unformatted')
      elseif (mode.eq.4) then
         call cdfopn(filename,fid,ierr)
      elseif (mode.eq.5) then
         print*,' ERROR: Reading KML not supported'
         stop
      elseif (mode.eq.6) then
         fid = 10
         open(fid,file=filename)   
      endif
 
c     Read header information
      call read_hea(fid,time,vars,ntra,ntim,ncol,mode)
 
      end
      
 
c     ----------------------------------------------------------------
c     Open a trajectory file for wrinting
c     ----------------------------------------------------------------
      
      subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)
      
      implicit none
      
c     Declaration of subroutine parameters
      integer      fid
      character*80 filename
      integer      mode
      integer      ntra,ntim,ncol
      integer      time(6) 
      character*80 vars(ncol)
 
c     Auxiliary variables
      integer      vardim(4)
      real         varmin(4),varmax(4),stag(4)
      real         mdv
      character*80 cfn
      integer      ierr
      integer      i
      character*80 varname
      real         rtime(6)
 
c     Open file
      if (mode.eq.1) then
         fid = 10
         open(fid,file=filename)
      elseif (mode.eq.2) then
         fid = 10
         open(fid,file=filename)
      elseif (mode.eq.3) then
         open(fid,file=filename,form='unformatted')
      elseif (mode.eq.4) then
         vardim(1)=ntra
         vardim(2)=1
         vardim(3)=1
         vardim(4)=1
         cfn      =trim(filename)//'_cst'
         mdv      =-999.98999
         call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)
      elseif (mode.eq.5) then
         fid = 10
         open(fid,file=filename)
      elseif (mode.eq.6) then
         fid = 10
         open(fid,file=filename)
      endif
 
c     Write header information
      call write_hea(fid,time,vars,ntra,ntim,ncol,mode)
 
      end
 
 
c     ----------------------------------------------------------------
c     Read a trajectory
c     ----------------------------------------------------------------
 
      subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer   fid
      integer   ntim
      integer   ncol
      integer   ntra
      real      tra(ntra,ntim,ncol)
      integer   mode
 
c     Auxiliary variables
      integer      i,j,n
      real         arr(ntra)
      integer      ntimes
      real         times(10000)
      integer      ierr
      character*80 vars(ncol+2)
      integer      nvars
 
c     Read ascii mode, sorted by trajectory (mode=1)
      if (mode.eq.1) then
         read(fid,*,end=100)
         do n=1,ntra
            do i=1,ntim
               read(fid,*,end=110) (tra(n,i,j),j=1,ncol)
            enddo
         enddo
 
c     Read ascii mode, sorted by time (mode=2)
      elseif (mode.eq.2) then
         read(fid,*,end=100)
         do i=1,ntim
            do n=1,ntra
               read(fid,*,end=100) (tra(n,i,j),j=1,ncol)
            enddo
         enddo
 
c     Read fortran mode (mode=3)
      elseif (mode.eq.3) then
         read(fid) tra
 
c     Read IVE netcdf mode (mode=4)
      elseif (mode.eq.4) then
          call gettimes(fid,times,ntimes,ierr)
          call getvars(fid,nvars,vars,ierr)
          do i=1,ntim
             do j=1,ncol
                if (j.eq.1) then
                   do n=1,ntra
                      tra(n,i,1)=times(i)
                   enddo
                else
                   call getdat(fid,vars(j),times(i),0,arr,ierr)
                   do n=1,ntra
                      tra(n,i,j)=arr(n)
                   enddo
                endif
            enddo
         enddo
         
c     Read ascii trajectory with pressure as real         
      elseif (mode.eq.6) then
         read(fid,*,end=100)
         do n=1,ntra
            do i=1,ntim
               read(fid,*,end=110) (tra(n,i,j),j=1,ncol)
            enddo
         enddo
 
      endif
 
      return
 
c     End of file has been reached: set negative <fid>
 100  fid=-fid
      return
 
c     Error: incomplete trajectory
 110  print*,'<read_tra>: Incomplete trajectory... Stop'
      stop
      
      end
 
 
c     ----------------------------------------------------------------
c     Write a trajectory
c     ----------------------------------------------------------------
 
      subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer   fid
      integer   ntim
      integer   ncol
      integer   ntra
      real      tra(ntra,ntim,ncol)
      integer   mode
 
c     Auxiliary variables
      integer      i,j,n
      real         arr(ntra)
      integer      ierr
      real         time
      character*80 vars(ncol+2)
      integer      nvars
      character*20 lonstr,latstr,levstr
      character*80 outstr
      real         ref_z(3000),ref_p(3000),ref_t(3000)
      real         lev
      character*80 path
 
c     Write ascii mode, sorted by trajectory (mode=1)
      if (mode.eq.1) then
         do n=1,ntra
            write(fid,*)
            do i=1,ntim
 
c              Avoid ugly *s or missing space in output
               do j=5,ncol
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
                     tra(n,i,j) = -999.99
                  endif
               enddo
 
               write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)') 
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
     >               nint(tra(n,i,4)),               ! p
     >               (tra(n,i,j),j=5,ncol)           ! fields
            enddo
         enddo
 
c     Write ascii mode, sorted by time (mode=2)
      elseif (mode.eq.2) then
         do i=1,ntim
            write(fid,*)
            do n=1,ntra
 
c              Avoid ugly *s or missing space in output
               do j=5,ncol
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
                     tra(n,i,j) = -999.99
                  endif
               enddo
 
               write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)') 
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
     >               nint(tra(n,i,4)),               ! p
     >               (tra(n,i,j),j=5,ncol)           ! fields
            enddo
         enddo
 
c     Write fortran mode (mode=3)
      elseif (mode.eq.3) then
         write(fid) tra                              
 
c     Write netcdf mode (mode=4)
      elseif (mode.eq.4) then
         call getvars(fid,nvars,vars,ierr)
         do i=1,ntim
            time=tra(1,i,1)
            do j=2,ncol
               do n=1,ntra
                  arr(n)=tra(n,i,j)
               enddo
               call putdat(fid,vars(j),time,0,arr,ierr)
            enddo
         enddo
 
c     Write KML mode (mode=5)
      elseif (mode.eq.5) then
 
         call getenv('DYN_TOOLS',path)
         path = trim(path)//'/lagranto.ecmwf/goodies/'
 
         open(fid+1,file=trim(path)//'reformat.refprof')
           
           do n=1,6
              read(fid+1,*)
           enddo
           do n=1,3000
              read(fid+1,*)  ref_z(n),ref_t(n),ref_p(n)
              ref_p(n) = 0.01 * ref_p(n)
           enddo
 
         close(fid+1)
 
         do n=1,ntra
           write(fid,"(A)") '<Placemark>'
           write(fid,"(A)") '<name>Absolute Extruded</name>'
           write(fid,"(A)") '<styleUrl>#yellowkLineGreenPoly</styleUrl>'
           write(fid,"(A)") '<LineString>'
           write(fid,"(A)") '<extrude>1</extrude>'
           write(fid,"(A)") '<tessellate>1</tessellate>'
           write(fid,"(A)") '<altitudeMode>absolute</altitudeMode>'
           write(fid,"(A)") '<coordinates>'
 
           do i=1,ntim
             write(lonstr,*) tra(n,i,2)
             write(latstr,*) tra(n,i,3)
             
             call binary(lev,tra(n,i,4),ref_z,ref_p)
             write(levstr,*) lev
 
             outstr = trim(adjustl(lonstr))//','//
     >                trim(adjustl(latstr))//','//
     >                trim(adjustl(levstr))
 
             write(fid,"(A)") outstr
 
           enddo
 
           write(fid,*) '</coordinates>'
           write(fid,*) '</LineString>'
           write(fid,*) '</Placemark>'
         enddo
         
c     Write ASCII trajectory with pressure as real         
      elseif (mode.eq.6) then
         do n=1,ntra
            write(fid,*)
            do i=1,ntim
 
c              Avoid ugly *s or missing space in output
               do j=5,ncol
                  if ( abs(tra(n,i,j)).gt.9999.) then                     
                     print*,'Format problem : ',tra(n,i,j),' -> -999.99'
                     tra(n,i,j) = -999.99
                  endif
               enddo
 
               write(fid,'(1f7.2,f9.2,f8.2,f8.2,100f10.3)') 
     >               (tra(n,i,j),j=1,3),             ! time, lon, lat
     >               tra(n,i,4),                     ! p
     >               (tra(n,i,j),j=5,ncol)           ! fields
            enddo
         enddo
      endif
 
      end
 
 
c     ----------------------------------------------------------------
c     Read header from trajectory file
c     ----------------------------------------------------------------
 
      subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer       fid
      integer       time(6)
      integer       ntra,ntim,ncol
      character*80  vars(ncol)
      integer       mode
 
c     Auxiliary variables
      integer       i
      character     ch(500)
      character*500 str
      integer       ich(500)
      integer       isstr,ileft,iright
      character*80  varname
      real          rtime(6)
      integer       ierr
      integer       nvars
      character*15  str1
      character     str2
      character*13  str3
      character*4   str4
      character*80  linestr
      integer       itmp1,itmp2
 
c     Read ascii format (mode=1,2,6)
      if ( (mode.eq.1).or.(mode.eq.2).or.(mode.eq.6) ) then
 
c        Read the time specification (old and new format)
         read(fid,'(a80)') linestr
         
         if ( linestr(1:14).eq.'Reference date' ) then
            read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)') 
     >           str1,
     >           time(1),time(2),time(3),str2,time(4),time(5),
     >           str3,time(6),str4
                       
         elseif ( linestr(1:11).eq.'time period' ) then
            read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)') 
     >           str1,
     >           time(1),time(2),time(3),str2,time(4),
     >           str3,itmp1,str3,itmp2,str4
            time(5) = 0
            time(6) = itmp1 * 60 + itmp2
 
         endif
 
c        Skip the empty line and read field names
         read(fid,*)
         read(fid,'(a500)',end=100) str
         do i=1,500
            ch(i)=str(i:i)
         enddo
 
c        Split the input string
         isstr=0
         nvars=0
         do i=1,500
            if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then
               isstr=1
               ileft=i
            elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then
               isstr=0
               iright=i-1
               nvars=nvars+1
               vars(nvars)=str(ileft:iright)
            endif
         enddo
 
c        Skip the empty line
         read(fid,*,end=100)
 
c     Read fortran mode (mode=3)
      elseif (mode.eq.3) then
         read(fid) ntra,ntim,ncol
         read(fid) time
         read(fid) vars
 
c     Read IVE netcdf mode (mode=4)
      elseif (mode.eq.4) then
         call getvars(fid,nvars,vars,ierr)  
         varname='BASEDATE'
         call getdat(fid,varname,0.,0,rtime,ierr)
         do i=1,6
            time(i)=nint(rtime(i))
         enddo
      endif
 
      return
 
c     End of file has been reached
 100  fid=-fid
      return
 
c     Excetion handling
 110  print*,'<read_hea>: Unexspected time format.... Stop'
      stop
   
      end
      
 
c     ----------------------------------------------------------------
c     Write header to trajectory file (in ascii mode)
c     ----------------------------------------------------------------
 
      subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer       fid
      integer       time(6)
      integer       ntra,ntim,ncol
      character*80  vars(ncol)
      integer       mode
 
c     Auxiliary variables
      integer       i
      character*500 str
      character*4   str1
      character*2   str2,str3,str4,str5,str6
      integer       vardim(4)
      real          varmin(4),varmax(4),stag(4)
      real          mdv
      integer       ierr
      character*80  varname
      real          rtime(6)
      integer       nvars
 
c     Write ascii format (mode=1,2)
      if ( (mode.eq.1).or.(mode.eq.2) ) then
 
c        Get the strings for output
         write(str1,'(i4)') time(1)
         write(str2,'(i2)') time(2)
         write(str3,'(i2)') time(3)
         write(str4,'(i2)') time(4)
         write(str5,'(i2)') time(5)
         if (time(2).eq. 0) str2(1:1)='0'
         if (time(3).eq. 0) str3(1:1)='0'
         if (time(4).eq. 0) str4(1:1)='0'
         if (time(5).eq. 0) str5(1:1)='0'
         if (time(2).lt.10) str2(1:1)='0'
         if (time(3).lt.10) str3(1:1)='0'
         if (time(4).lt.10) str4(1:1)='0'
         if (time(5).lt.10) str5(1:1)='0'
 
c        Write the time specification
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
     >          'Reference date ',
     >           str1,str2,str3,'_',str4,str5,
     >          ' / Time range',time(6), ' min'
         write(fid,*)
 
c        Write variable names
         str=''
         do i=1,ncol
            varname = vars(i)
            vars(i) = varname(1:9)
            str     = trim(str)//trim(vars(i))
         enddo
         write(fid,'(a6,a9,a8,a6,100a10)') (trim(vars(i)),i=1,ncol)
         write(fid,'(a6,a9,a8,a6,100a10)') 
     >              '------','---------','--------','------',
     >              ('----------',i=5,ncol)
 
c     Write fortran mode (mode=3)
      elseif (mode.eq.3) then
         write(fid) ntra,ntim,ncol
         write(fid) time
         write(fid) vars
 
c     Write IVE netcdf format (mode=4)
      elseif (mode.eq.4) then
         vardim(1)=ntra
         vardim(2)=1
         vardim(3)=1
         vardim(4)=1
         mdv      =-999.98999
         do i=2,ncol
            call putdef(fid,vars(i),4,mdv,vardim,
     >                  varmin,varmax,stag,ierr)
         enddo
         varname='BASEDATE'
         vardim(1)=6
         call putdef(fid,varname,4,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         do i=1,6
            rtime(i)=real(time(i))
         enddo
         call putdat(fid,varname,0.,0,rtime,ierr)
 
c     Write KML format (mode=5)
      elseif (mode.eq.5) then
 
      write(fid,"(A)") '<?xml version="1.0" encoding="UTF-8"?>'
      write(fid,"(A)") '<kml xmlns="http://www.opengis.net/kml/2.2">'
      write(fid,"(A)") '<Document>'
      write(fid,"(A)") '<name>Paths</name>'
      write(fid,"(A)") '<Style id="yellowLineGreenPoly">'
      write(fid,"(A)") '<LineStyle>'
c      write(fid,*) '<color>7f00ffff</color>'    ! Yellow
      write(fid,"(A)") '<color>500A0A0A</color>'     ! Black
      write(fid,"(A)") '<width>4</width>'
      write(fid,"(A)") '</LineStyle>'
      write(fid,"(A)") '<PolyStyle>'
      write(fid,"(A)") '<color>7f00ff00</color>'
      write(fid,"(A)") '</PolyStyle>'
      write(fid,"(A)") '</Style>'
      
c     Write header for ASCII mode iwth real pressure      
      elseif ( mode.eq.6 ) then
 
c        Get the strings for output
         write(str1,'(i4)') time(1)
         write(str2,'(i2)') time(2)
         write(str3,'(i2)') time(3)
         write(str4,'(i2)') time(4)
         write(str5,'(i2)') time(5)
         if (time(2).eq. 0) str2(1:1)='0'
         if (time(3).eq. 0) str3(1:1)='0'
         if (time(4).eq. 0) str4(1:1)='0'
         if (time(5).eq. 0) str5(1:1)='0'
         if (time(2).lt.10) str2(1:1)='0'
         if (time(3).lt.10) str3(1:1)='0'
         if (time(4).lt.10) str4(1:1)='0'
         if (time(5).lt.10) str5(1:1)='0'
 
c        Write the time specification
         write(fid,'(a15,a4,a2,a2,a1,a2,a2,a13,i8,a4)') 
     >          'Reference date ',
     >           str1,str2,str3,'_',str4,str5,
     >          ' / Time range',time(6), ' min'
         write(fid,*)
 
c        Write variable names
         str=''
         do i=1,ncol
            varname = vars(i)
            vars(i) = varname(1:9)
            str     = trim(str)//trim(vars(i))
         enddo
         write(fid,'(a6,a9,a8,a8,100a10)') (trim(vars(i)),i=1,ncol)
         write(fid,'(a6,a9,a8,a8,100a10)') 
     >              '------','---------','--------','--------',
     >              ('----------',i=5,ncol)
 
      endif
 
      end
      
      
c     ----------------------------------------------------------------
c     Close a trajectory file
c     ----------------------------------------------------------------
      
      subroutine close_tra(fid,mode)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer      fid
      integer      mode
      
c     Auxiliary variables
      integer      ierr
 
c     Close file
      if (mode.eq.1) then
         close(abs(fid))
 
      elseif (mode.eq.2) then
         close(abs(fid))
 
      elseif (mode.eq.3) then
         close(fid)
 
      elseif (mode.eq.4) then
         call clscdf(fid,ierr)
 
      elseif (mode.eq.5) then
          write(fid,"(A)") '</Document>'
          write(fid,"(A)") '</kml>'
         close(abs(fid))
         
      elseif (mode.eq.6) then
         close(abs(fid))
         
      endif
 
      end
      
 
c     ----------------------------------------------------------------
c     Determine the mode of a trajectory file
c     ----------------------------------------------------------------
 
      subroutine mode_tra(mode,filename)
      
      implicit none
 
c     Declaration of subroutine parameters
      integer        mode
      character*80   filename
 
c     Auxiliary variables
      integer        len
      character      char0,char1,char2,char3,char4
 
c     Get mode
      mode=-1
      
      len  = len_trim(filename)
 
c     Mode specified by number
      char0 = filename((len-1):(len-1))
      char1 = filename(len:len)
 
      if ( (char0.eq.'.').and.(char1.eq.'1') ) mode=1
      if ( (char0.eq.'.').and.(char1.eq.'2') ) mode=2
      if ( (char0.eq.'.').and.(char1.eq.'3') ) mode=3
      if ( (char0.eq.'.').and.(char1.eq.'4') ) mode=4
      if ( (char0.eq.'.').and.(char1.eq.'5') ) mode=5
      if ( (char0.eq.'.').and.(char1.eq.'6') ) mode=6
 
      if ( mode.gt.0 ) return
 
c     Mode specified by appendix
      char0 = filename((len-3):(len-3))
      char1 = filename((len-2):(len-2))
      char2 = filename((len-1):(len-1))
      char3 = filename(len:len)
      if ( (char1.eq.'.').and.(char2.eq.'l').and.(char3.eq.'s') ) mode=1
      if ( (char1.eq.'.').and.(char2.eq.'t').and.(char3.eq.'i') ) mode=2
      if ( (char1.eq.'.').and.(char2.eq.'d').and.(char3.eq.'u') ) mode=3
      if ( (char1.eq.'.').and.(char2.eq.'n').and.(char3.eq.'c') ) mode=4
 
      if ( (char0.eq.'.').and.(char1.eq.'k').and.
     >                        (char2.eq.'m').and.
     >                        (char3.eq.'l') ) mode = 5
 
      if ( (char1.eq.'.').and.(char2.eq.'r').and.(char3.eq.'p') ) mode=6
 
      end
 
 
c     ----------------------------------------------------------------
c     Get dimension of a trajectory file
c     ----------------------------------------------------------------
    
      subroutine info_tra(filename,ntra,ntim,ncol,mode)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer      fid
      character*80 filename
      integer      mode
      integer      ntra,ntim,ncol
 
c     Auxiliary variables
      integer       vardim(4)
      real          varmin(4),varmax(4),stag(4)
      real          mdv
      character*80  cfn
      integer       ierr
      integer       i,ndim
      integer       nvars
      integer       ntimes
      real          times(100)
      character*500 str
      integer       nline0,nline1,nline2
      integer       isstr,isok
      character     ch
      character*80  vars(200)
      integer       ios
 
c     Open file
      if (mode.eq.1) then
         fid=10
         open(fid,file=filename)
      elseif (mode.eq.2) then
         fid=10
         open(fid,file=filename)
      elseif (mode.eq.3) then
         fid=10
         open(fid,file=filename,form='unformatted')
      elseif (mode.eq.4) then
         call cdfopn(filename,fid,ierr)
      elseif (mode.eq.6) then
         fid=10
         open(fid,file=filename)
      endif
 
c     Get dimension information
      if ( (mode.eq.1).or.(mode.eq.2).or.(mode.eq.6) ) then
         read(fid,*)
         read(fid,*)
         read(fid,'(a500)') str
         read(fid,*)
 
c        Get the number of columns
         isstr=0
         ncol =0
         do i=1,500
            ch = str(i:i)
            if ( (isstr.eq.0).and.(ch.ne.' ') ) then
               isstr=1
            elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
               isstr=0
               ncol=ncol+1
            endif
         enddo
 
c        Init the dimensions for empty trajectory files
         ntra = 0
         ntim = 0
 
c        Get the first data block
         nline0  = 5
         nline1  = 5
         read(fid,*,end=140,iostat=ios)
 100     read(fid,'(a500)',end=110,iostat=ios) str
         if (str.ne.'') then
            nline1 = nline1 + 1
            goto 100
         endif
 110     continue
         
c        Get the total numbers of lines in the data block
         nline2 = nline1
         if ( ios.eq.0 ) then
 120        read(fid,*,end=130)
            nline2 = nline2 + 1
            goto 120
 130        nline2 = nline2 + 1
         endif
 
c        Set the dimensions
         if (mode.eq.1) then
            ntim = nline1 - nline0
            ntra = (nline2-nline0+1)/(ntim+1)
         else
            ntra = nline1 - nline0
            ntim = (nline2-nline0+1)/(ntra+1)
         endif
 
      elseif (mode.eq.3) then
         read(fid) ntra,ntim,ncol
 
      elseif (mode.eq.4) then
         call gettimes(fid,times,ntimes,ierr)
         call getvars(fid,nvars,vars,ierr)
         call getdef(fid,trim(vars(2)),ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         ntra = vardim(1)
         ntim = ntimes
         ncol = nvars-1
 
      endif
 
c     Close file
 140  continue
      if (mode.eq.1) then
         close(fid)
      elseif (mode.eq.2) then
         close(fid)
      elseif (mode.eq.3) then
         close(fid)
      elseif (mode.eq.4) then
         call clscdf(fid,ierr)
      endif
 
      end
     
 
c     ----------------------------------------------------------------
c     Binary search algorithm
c     ----------------------------------------------------------------
 
      subroutine binary (z,p,ref_z,ref_p)
 
      implicit none
 
c     Declaration of subroutine parameters
      real   z
      real   p
      real   ref_z(3000)
      real   ref_p(3000)
 
c     Auxiliary variables
      integer i0,i1,im
 
c     Binary search
      i0 = 1
      i1 = 3000
 100  continue
        im = (i0 + i1) / 2
        if ( p.lt.ref_p(im) ) then
           i0 = im
        else 
           i1 = im
        endif
      if ( (i1-i0).gt.1 ) goto 100
 
c     Linear interpolation in between
      z = ref_z(i0) + ( p - ref_p(i0) ) / ( ref_p(i1) - ref_p(i0) ) *
     >                ( ref_z(i1) - ref_z(i0) ) 
 
      end
 
 
 
 
 

Rev 15	Rev 39
1	`c ****************************************************************`	1	`c ****************************************************************`
2	`c * This package provides IO routines for trajectories. A file *`	2	`c * This package provides IO routines for trajectories. A file *`
3	`c * is characterised by the filename <filename> and the file *`	3	`c * is characterised by the filename <filename> and the file *`
4	`c * identifier <fid>. Different modes <mode> are supported: *`	4	`c * identifier <fid>. Different modes <mode> are supported: *`
5	`c * mode=1: ascii, sorted by trajectory; *`	5	`c * mode=1: ascii, sorted by trajectory; *`
6	`c * mode=2: ascii, sorted by time; *`	6	`c * mode=2: ascii, sorted by time; *`
7	`c * mode=3: fortran (unformatted) *`	7	`c * mode=3: fortran (unformatted) *`
8	`c * mode=4: IVE netcdf (for compatibiltzy reasons) *`	8	`c * mode=4: IVE netcdf (for compatibiltzy reasons) *`
-		9	`c * mode=5: KML *`
-		10	`c * mode=6: as mode 1, but with pressure as realk number *`
9	`c * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)> *`	11	`c * A trajectory set is given as 3d array <tra(ntra,ntim,ncol)> *`
10	`c * where <ntra> is the number of trajectories, <ntim> the *`	12	`c * where <ntra> is the number of trajectories, <ntim> the *`
11	`c * number of times of each trajectory and <ncol> the number of *`	13	`c * number of times of each trajectory and <ncol> the number of *`
12	`c * columns of the trajectory. The first 4 columns are: time, *`	14	`c * columns of the trajectory. The first 4 columns are: time, *`
13	`c * longitude, latitude, pressure. The other columns are traced *`	15	`c * longitude, latitude, pressure. The other columns are traced *`
14	`c * fields. The complete list of all columns is given in the *`	16	`c * fields. The complete list of all columns is given in the *`
15	`c * array <vars(ncol)>. Finally, the reference date is given in *`	17	`c * array <vars(ncol)>. Finally, the reference date is given in *`
16	`c * the array <time(6)=year,month,day,hour,time length of the *`	18	`c * the array <time(6)=year,month,day,hour,time length of the *`
17	`c * trajectory (hour,min)>. *`	19	`c * trajectory (hour,min)>. *`
18	`c * *`	20	`c * *`
19	`c * Author: Michael Sprenger, September 2008 *`	21	`c * Author: Michael Sprenger, September 2008 *`
20	`c ****************************************************************`	22	`c ****************************************************************`
21		23
22	`c ----------------------------------------------------------------`	24	`c ----------------------------------------------------------------`
23	`c Open a trajectory file for reading`	25	`c Open a trajectory file for reading`
24	`c ----------------------------------------------------------------`	26	`c ----------------------------------------------------------------`
25		27
26	`subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`	28	`subroutine ropen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`
27		29
28	`implicit none`	30	`implicit none`
29		31
30	`c Declaration of subroutine parameters`	32	`c Declaration of subroutine parameters`
31	`integer fid`	33	`integer fid`
32	`character*80 filename`	34	`character*80 filename`
33	`integer mode`	35	`integer mode`
34	`integer ntra,ntim,ncol`	36	`integer ntra,ntim,ncol`
35	`integer time(6)`	37	`integer time(6)`
36	`character*80 vars(ncol)`	38	`character*80 vars(ncol)`
37		39
38	`c Auxiliary variables`	40	`c Auxiliary variables`
39	`integer vardim(4)`	41	`integer vardim(4)`
40	`real varmin(4),varmax(4),stag(4)`	42	`real varmin(4),varmax(4),stag(4)`
41	`real mdv`	43	`real mdv`
42	`character*80 cfn`	44	`character*80 cfn`
43	`integer ierr`	45	`integer ierr`
44	`integer i`	46	`integer i`
45	`integer nvars`	47	`integer nvars`
46		48
47	`c Open file`	49	`c Open file`
48	`if (mode.eq.1) then`	50	`if (mode.eq.1) then`
49	`fid = 10`	51	`fid = 10`
50	`open(fid,file=filename)`	52	`open(fid,file=filename)`
51	`elseif (mode.eq.2) then`	53	`elseif (mode.eq.2) then`
52	`fid = 10`	54	`fid = 10`
53	`open(fid,file=filename)`	55	`open(fid,file=filename)`
54	`elseif (mode.eq.3) then`	56	`elseif (mode.eq.3) then`
55	`open(fid,file=filename,form='unformatted')`	57	`open(fid,file=filename,form='unformatted')`
56	`elseif (mode.eq.4) then`	58	`elseif (mode.eq.4) then`
57	`call cdfopn(filename,fid,ierr)`	59	`call cdfopn(filename,fid,ierr)`
58	`elseif (mode.eq.5) then`	60	`elseif (mode.eq.5) then`
59	`print*,' ERROR: Reading KML not supported'`	61	`print*,' ERROR: Reading KML not supported'`
60	`stop`	62	`stop`
-		63	`elseif (mode.eq.6) then`
-		64	`fid = 10`
-		65	`open(fid,file=filename)`
61	`endif`	66	`endif`
62		67
63	`c Read header information`	68	`c Read header information`
64	`call read_hea(fid,time,vars,ntra,ntim,ncol,mode)`	69	`call read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
65		70
66	`end`	71	`end`
67		72
68		73
69	`c ----------------------------------------------------------------`	74	`c ----------------------------------------------------------------`
70	`c Open a trajectory file for wrinting`	75	`c Open a trajectory file for wrinting`
71	`c ----------------------------------------------------------------`	76	`c ----------------------------------------------------------------`
72		77
73	`subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`	78	`subroutine wopen_tra(fid,filename,ntra,ntim,ncol,time,vars,mode)`
74		79
75	`implicit none`	80	`implicit none`
76		81
77	`c Declaration of subroutine parameters`	82	`c Declaration of subroutine parameters`
78	`integer fid`	83	`integer fid`
79	`character*80 filename`	84	`character*80 filename`
80	`integer mode`	85	`integer mode`
81	`integer ntra,ntim,ncol`	86	`integer ntra,ntim,ncol`
82	`integer time(6)`	87	`integer time(6)`
83	`character*80 vars(ncol)`	88	`character*80 vars(ncol)`
84		89
85	`c Auxiliary variables`	90	`c Auxiliary variables`
86	`integer vardim(4)`	91	`integer vardim(4)`
87	`real varmin(4),varmax(4),stag(4)`	92	`real varmin(4),varmax(4),stag(4)`
88	`real mdv`	93	`real mdv`
89	`character*80 cfn`	94	`character*80 cfn`
90	`integer ierr`	95	`integer ierr`
91	`integer i`	96	`integer i`
92	`character*80 varname`	97	`character*80 varname`
93	`real rtime(6)`	98	`real rtime(6)`
94		99
95	`c Open file`	100	`c Open file`
96	`if (mode.eq.1) then`	101	`if (mode.eq.1) then`
97	`fid = 10`	102	`fid = 10`
98	`open(fid,file=filename)`	103	`open(fid,file=filename)`
99	`elseif (mode.eq.2) then`	104	`elseif (mode.eq.2) then`
100	`fid = 10`	105	`fid = 10`
101	`open(fid,file=filename)`	106	`open(fid,file=filename)`
102	`elseif (mode.eq.3) then`	107	`elseif (mode.eq.3) then`
103	`open(fid,file=filename,form='unformatted')`	108	`open(fid,file=filename,form='unformatted')`
104	`elseif (mode.eq.4) then`	109	`elseif (mode.eq.4) then`
105	`vardim(1)=ntra`	110	`vardim(1)=ntra`
106	`vardim(2)=1`	111	`vardim(2)=1`
107	`vardim(3)=1`	112	`vardim(3)=1`
108	`vardim(4)=1`	113	`vardim(4)=1`
109	`cfn =trim(filename)//'_cst'`	114	`cfn =trim(filename)//'_cst'`
110	`mdv =-999.98999`	115	`mdv =-999.98999`
111	`call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)`	116	`call crecdf(filename,fid,varmin,varmax,3,cfn,ierr)`
112	`elseif (mode.eq.5) then`	117	`elseif (mode.eq.5) then`
113	`fid = 10`	118	`fid = 10`
114	`open(fid,file=filename)`	119	`open(fid,file=filename)`
-		120	`elseif (mode.eq.6) then`
-		121	`fid = 10`
-		122	`open(fid,file=filename)`
115	`endif`	123	`endif`
116		124
117	`c Write header information`	125	`c Write header information`
118	`call write_hea(fid,time,vars,ntra,ntim,ncol,mode)`	126	`call write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
119		127
120	`end`	128	`end`
121		129
122		130
123	`c ----------------------------------------------------------------`	131	`c ----------------------------------------------------------------`
124	`c Read a trajectory`	132	`c Read a trajectory`
125	`c ----------------------------------------------------------------`	133	`c ----------------------------------------------------------------`
126		134
127	`subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)`	135	`subroutine read_tra(fid,tra,ntra,ntim,ncol,mode)`
128		136
129	`implicit none`	137	`implicit none`
130		138
131	`c Declaration of subroutine parameters`	139	`c Declaration of subroutine parameters`
132	`integer fid`	140	`integer fid`
133	`integer ntim`	141	`integer ntim`
134	`integer ncol`	142	`integer ncol`
135	`integer ntra`	143	`integer ntra`
136	`real tra(ntra,ntim,ncol)`	144	`real tra(ntra,ntim,ncol)`
137	`integer mode`	145	`integer mode`
138		146
139	`c Auxiliary variables`	147	`c Auxiliary variables`
140	`integer i,j,n`	148	`integer i,j,n`
141	`real arr(ntra)`	149	`real arr(ntra)`
142	`integer ntimes`	150	`integer ntimes`
143	`real times(10000)`	151	`real times(10000)`
144	`integer ierr`	152	`integer ierr`
145	`character*80 vars(ncol+2)`	153	`character*80 vars(ncol+2)`
146	`integer nvars`	154	`integer nvars`
147		155
148	`c Read ascii mode, sorted by trajectory (mode=1)`	156	`c Read ascii mode, sorted by trajectory (mode=1)`
149	`if (mode.eq.1) then`	157	`if (mode.eq.1) then`
150	`read(fid,*,end=100)`	158	`read(fid,*,end=100)`
151	`do n=1,ntra`	159	`do n=1,ntra`
152	`do i=1,ntim`	160	`do i=1,ntim`
153	`read(fid,*,end=110) (tra(n,i,j),j=1,ncol)`	161	`read(fid,*,end=110) (tra(n,i,j),j=1,ncol)`
154	`enddo`	162	`enddo`
155	`enddo`	163	`enddo`
156		164
157	`c Read ascii mode, sorted by time (mode=2)`	165	`c Read ascii mode, sorted by time (mode=2)`
158	`elseif (mode.eq.2) then`	166	`elseif (mode.eq.2) then`
159	`read(fid,*,end=100)`	167	`read(fid,*,end=100)`
160	`do i=1,ntim`	168	`do i=1,ntim`
161	`do n=1,ntra`	169	`do n=1,ntra`
162	`read(fid,*,end=100) (tra(n,i,j),j=1,ncol)`	170	`read(fid,*,end=100) (tra(n,i,j),j=1,ncol)`
163	`enddo`	171	`enddo`
164	`enddo`	172	`enddo`
165		173
166	`c Read fortran mode (mode=3)`	174	`c Read fortran mode (mode=3)`
167	`elseif (mode.eq.3) then`	175	`elseif (mode.eq.3) then`
168	`read(fid) tra`	176	`read(fid) tra`
169		177
170	`c Read IVE netcdf mode (mode=4)`	178	`c Read IVE netcdf mode (mode=4)`
171	`elseif (mode.eq.4) then`	179	`elseif (mode.eq.4) then`
172	`call gettimes(fid,times,ntimes,ierr)`	180	`call gettimes(fid,times,ntimes,ierr)`
173	`call getvars(fid,nvars,vars,ierr)`	181	`call getvars(fid,nvars,vars,ierr)`
174	`do i=1,ntim`	182	`do i=1,ntim`
175	`do j=1,ncol`	183	`do j=1,ncol`
176	`if (j.eq.1) then`	184	`if (j.eq.1) then`
177	`do n=1,ntra`	185	`do n=1,ntra`
178	`tra(n,i,1)=times(i)`	186	`tra(n,i,1)=times(i)`
179	`enddo`	187	`enddo`
180	`else`	188	`else`
181	`call getdat(fid,vars(j),times(i),0,arr,ierr)`	189	`call getdat(fid,vars(j),times(i),0,arr,ierr)`
182	`do n=1,ntra`	190	`do n=1,ntra`
183	`tra(n,i,j)=arr(n)`	191	`tra(n,i,j)=arr(n)`
184	`enddo`	192	`enddo`
185	`endif`	193	`endif`
186	`enddo`	194	`enddo`
187	`enddo`	195	`enddo`
-		196
-		197	`c Read ascii trajectory with pressure as real`
-		198	`elseif (mode.eq.6) then`
-		199	`read(fid,*,end=100)`
-		200	`do n=1,ntra`
-		201	`do i=1,ntim`
-		202	`read(fid,*,end=110) (tra(n,i,j),j=1,ncol)`
-		203	`enddo`
-		204	`enddo`
188		205
189	`endif`	206	`endif`
190		207
191	`return`	208	`return`
192		209
193	`c End of file has been reached: set negative <fid>`	210	`c End of file has been reached: set negative <fid>`
194	`100 fid=-fid`	211	`100 fid=-fid`
195	`return`	212	`return`
196		213
197	`c Error: incomplete trajectory`	214	`c Error: incomplete trajectory`
198	`110 print*,'<read_tra>: Incomplete trajectory... Stop'`	215	`110 print*,'<read_tra>: Incomplete trajectory... Stop'`
199	`stop`	216	`stop`
200		217
201	`end`	218	`end`
202		219
203		220
204	`c ----------------------------------------------------------------`	221	`c ----------------------------------------------------------------`
205	`c Write a trajectory`	222	`c Write a trajectory`
206	`c ----------------------------------------------------------------`	223	`c ----------------------------------------------------------------`
207		224
208	`subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)`	225	`subroutine write_tra(fid,tra,ntra,ntim,ncol,mode)`
209		226
210	`implicit none`	227	`implicit none`
211		228
212	`c Declaration of subroutine parameters`	229	`c Declaration of subroutine parameters`
213	`integer fid`	230	`integer fid`
214	`integer ntim`	231	`integer ntim`
215	`integer ncol`	232	`integer ncol`
216	`integer ntra`	233	`integer ntra`
217	`real tra(ntra,ntim,ncol)`	234	`real tra(ntra,ntim,ncol)`
218	`integer mode`	235	`integer mode`
219		236
220	`c Auxiliary variables`	237	`c Auxiliary variables`
221	`integer i,j,n`	238	`integer i,j,n`
222	`real arr(ntra)`	239	`real arr(ntra)`
223	`integer ierr`	240	`integer ierr`
224	`real time`	241	`real time`
225	`character*80 vars(ncol+2)`	242	`character*80 vars(ncol+2)`
226	`integer nvars`	243	`integer nvars`
227	`character*20 lonstr,latstr,levstr`	244	`character*20 lonstr,latstr,levstr`
228	`character*80 outstr`	245	`character*80 outstr`
229	`real ref_z(3000),ref_p(3000),ref_t(3000)`	246	`real ref_z(3000),ref_p(3000),ref_t(3000)`
230	`real lev`	247	`real lev`
231	`character*80 path`	248	`character*80 path`
232		249
233	`c Write ascii mode, sorted by trajectory (mode=1)`	250	`c Write ascii mode, sorted by trajectory (mode=1)`
234	`if (mode.eq.1) then`	251	`if (mode.eq.1) then`
235	`do n=1,ntra`	252	`do n=1,ntra`
236	`write(fid,*)`	253	`write(fid,*)`
237	`do i=1,ntim`	254	`do i=1,ntim`
238		255
239	`c Avoid ugly *s or missing space in output`	256	`c Avoid ugly *s or missing space in output`
240	`do j=5,ncol`	257	`do j=5,ncol`
241	`if ( abs(tra(n,i,j)).gt.9999.) then`	258	`if ( abs(tra(n,i,j)).gt.9999.) then`
242	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`	259	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
243	`tra(n,i,j) = -999.99`	260	`tra(n,i,j) = -999.99`
244	`endif`	261	`endif`
245	`enddo`	262	`enddo`
246		263
247	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`	264	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
248	`> (tra(n,i,j),j=1,3), ! time, lon, lat`	265	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
249	`> nint(tra(n,i,4)), ! p`	266	`> nint(tra(n,i,4)), ! p`
250	`> (tra(n,i,j),j=5,ncol) ! fields`	267	`> (tra(n,i,j),j=5,ncol) ! fields`
251	`enddo`	268	`enddo`
252	`enddo`	269	`enddo`
253		270
254	`c Write ascii mode, sorted by time (mode=2)`	271	`c Write ascii mode, sorted by time (mode=2)`
255	`elseif (mode.eq.2) then`	272	`elseif (mode.eq.2) then`
256	`do i=1,ntim`	273	`do i=1,ntim`
257	`write(fid,*)`	274	`write(fid,*)`
258	`do n=1,ntra`	275	`do n=1,ntra`
259		276
260	`c Avoid ugly *s or missing space in output`	277	`c Avoid ugly *s or missing space in output`
261	`do j=5,ncol`	278	`do j=5,ncol`
262	`if ( abs(tra(n,i,j)).gt.9999.) then`	279	`if ( abs(tra(n,i,j)).gt.9999.) then`
263	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`	280	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
264	`tra(n,i,j) = -999.99`	281	`tra(n,i,j) = -999.99`
265	`endif`	282	`endif`
266	`enddo`	283	`enddo`
267		284
268	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`	285	`write(fid,'(1f7.2,f9.2,f8.2,i6,100f10.3)')`
269	`> (tra(n,i,j),j=1,3), ! time, lon, lat`	286	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
270	`> nint(tra(n,i,4)), ! p`	287	`> nint(tra(n,i,4)), ! p`
271	`> (tra(n,i,j),j=5,ncol) ! fields`	288	`> (tra(n,i,j),j=5,ncol) ! fields`
272	`enddo`	289	`enddo`
273	`enddo`	290	`enddo`
274		291
275	`c Write fortran mode (mode=3)`	292	`c Write fortran mode (mode=3)`
276	`elseif (mode.eq.3) then`	293	`elseif (mode.eq.3) then`
277	`write(fid) tra`	294	`write(fid) tra`
278		295
279	`c Write netcdf mode (mode=4)`	296	`c Write netcdf mode (mode=4)`
280	`elseif (mode.eq.4) then`	297	`elseif (mode.eq.4) then`
281	`call getvars(fid,nvars,vars,ierr)`	298	`call getvars(fid,nvars,vars,ierr)`
282	`do i=1,ntim`	299	`do i=1,ntim`
283	`time=tra(1,i,1)`	300	`time=tra(1,i,1)`
284	`do j=2,ncol`	301	`do j=2,ncol`
285	`do n=1,ntra`	302	`do n=1,ntra`
286	`arr(n)=tra(n,i,j)`	303	`arr(n)=tra(n,i,j)`
287	`enddo`	304	`enddo`
288	`call putdat(fid,vars(j),time,0,arr,ierr)`	305	`call putdat(fid,vars(j),time,0,arr,ierr)`
289	`enddo`	306	`enddo`
290	`enddo`	307	`enddo`
291		308
292	`c Write KML mode (mode=5)`	309	`c Write KML mode (mode=5)`
293	`elseif (mode.eq.5) then`	310	`elseif (mode.eq.5) then`
294		311
295	`call getenv('DYN_TOOLS',path)`	312	`call getenv('DYN_TOOLS',path)`
296	`path = trim(path)//'/lagranto.ecmwf/goodies/'`	313	`path = trim(path)//'/lagranto.ecmwf/goodies/'`
297		314
298	`open(fid+1,file=trim(path)//'reformat.refprof')`	315	`open(fid+1,file=trim(path)//'reformat.refprof')`
299		316
300	`do n=1,6`	317	`do n=1,6`
301	`read(fid+1,*)`	318	`read(fid+1,*)`
302	`enddo`	319	`enddo`
303	`do n=1,3000`	320	`do n=1,3000`
304	`read(fid+1,*) ref_z(n),ref_t(n),ref_p(n)`	321	`read(fid+1,*) ref_z(n),ref_t(n),ref_p(n)`
305	`ref_p(n) = 0.01 * ref_p(n)`	322	`ref_p(n) = 0.01 * ref_p(n)`
306	`enddo`	323	`enddo`
307		324
308	`close(fid+1)`	325	`close(fid+1)`
309		326
310	`do n=1,ntra`	327	`do n=1,ntra`
311	`write(fid,"(A)") '<Placemark>'`	328	`write(fid,"(A)") '<Placemark>'`
312	`write(fid,"(A)") '<name>Absolute Extruded</name>'`	329	`write(fid,"(A)") '<name>Absolute Extruded</name>'`
313	`write(fid,"(A)") '<styleUrl>#yellowkLineGreenPoly</styleUrl>'`	330	`write(fid,"(A)") '<styleUrl>#yellowkLineGreenPoly</styleUrl>'`
314	`write(fid,"(A)") '<LineString>'`	331	`write(fid,"(A)") '<LineString>'`
315	`write(fid,"(A)") '<extrude>1</extrude>'`	332	`write(fid,"(A)") '<extrude>1</extrude>'`
316	`write(fid,"(A)") '<tessellate>1</tessellate>'`	333	`write(fid,"(A)") '<tessellate>1</tessellate>'`
317	`write(fid,"(A)") '<altitudeMode>absolute</altitudeMode>'`	334	`write(fid,"(A)") '<altitudeMode>absolute</altitudeMode>'`
318	`write(fid,"(A)") '<coordinates>'`	335	`write(fid,"(A)") '<coordinates>'`
319		336
320	`do i=1,ntim`	337	`do i=1,ntim`
321	`write(lonstr,*) tra(n,i,2)`	338	`write(lonstr,*) tra(n,i,2)`
322	`write(latstr,*) tra(n,i,3)`	339	`write(latstr,*) tra(n,i,3)`
323		340
324	`call binary(lev,tra(n,i,4),ref_z,ref_p)`	341	`call binary(lev,tra(n,i,4),ref_z,ref_p)`
325	`write(levstr,*) lev`	342	`write(levstr,*) lev`
326		343
327	`outstr = trim(adjustl(lonstr))//','//`	344	`outstr = trim(adjustl(lonstr))//','//`
328	`> trim(adjustl(latstr))//','//`	345	`> trim(adjustl(latstr))//','//`
329	`> trim(adjustl(levstr))`	346	`> trim(adjustl(levstr))`
330		347
331	`write(fid,"(A)") outstr`	348	`write(fid,"(A)") outstr`
332		349
333	`enddo`	350	`enddo`
334		351
335	`write(fid,*) '</coordinates>'`	352	`write(fid,*) '</coordinates>'`
336	`write(fid,*) '</LineString>'`	353	`write(fid,*) '</LineString>'`
337	`write(fid,*) '</Placemark>'`	354	`write(fid,*) '</Placemark>'`
338	`enddo`	355	`enddo`
-		356
-		357	`c Write ASCII trajectory with pressure as real`
-		358	`elseif (mode.eq.6) then`
-		359	`do n=1,ntra`
-		360	`write(fid,*)`
-		361	`do i=1,ntim`
339		362
-		363	`c Avoid ugly *s or missing space in output`
-		364	`do j=5,ncol`
-		365	`if ( abs(tra(n,i,j)).gt.9999.) then`
-		366	`print*,'Format problem : ',tra(n,i,j),' -> -999.99'`
-		367	`tra(n,i,j) = -999.99`
-		368	`endif`
-		369	`enddo`
-		370
-		371	`write(fid,'(1f7.2,f9.2,f8.2,f8.2,100f10.3)')`
-		372	`> (tra(n,i,j),j=1,3), ! time, lon, lat`
-		373	`> tra(n,i,4), ! p`
-		374	`> (tra(n,i,j),j=5,ncol) ! fields`
-		375	`enddo`
-		376	`enddo`
340	`endif`	377	`endif`
341		378
342	`end`	379	`end`
343		380
344		381
345	`c ----------------------------------------------------------------`	382	`c ----------------------------------------------------------------`
346	`c Read header from trajectory file`	383	`c Read header from trajectory file`
347	`c ----------------------------------------------------------------`	384	`c ----------------------------------------------------------------`
348		385
349	`subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)`	386	`subroutine read_hea(fid,time,vars,ntra,ntim,ncol,mode)`
350		387
351	`implicit none`	388	`implicit none`
352		389
353	`c Declaration of subroutine parameters`	390	`c Declaration of subroutine parameters`
354	`integer fid`	391	`integer fid`
355	`integer time(6)`	392	`integer time(6)`
356	`integer ntra,ntim,ncol`	393	`integer ntra,ntim,ncol`
357	`character*80 vars(ncol)`	394	`character*80 vars(ncol)`
358	`integer mode`	395	`integer mode`
359		396
360	`c Auxiliary variables`	397	`c Auxiliary variables`
361	`integer i`	398	`integer i`
362	`character ch(500)`	399	`character ch(500)`
363	`character*500 str`	400	`character*500 str`
364	`integer ich(500)`	401	`integer ich(500)`
365	`integer isstr,ileft,iright`	402	`integer isstr,ileft,iright`
366	`character*80 varname`	403	`character*80 varname`
367	`real rtime(6)`	404	`real rtime(6)`
368	`integer ierr`	405	`integer ierr`
369	`integer nvars`	406	`integer nvars`
370	`character*15 str1`	407	`character*15 str1`
371	`character str2`	408	`character str2`
372	`character*13 str3`	409	`character*13 str3`
373	`character*4 str4`	410	`character*4 str4`
374	`character*80 linestr`	411	`character*80 linestr`
375	`integer itmp1,itmp2`	412	`integer itmp1,itmp2`
376		413
377	`c Read ascii format (mode=1,2)`	414	`c Read ascii format (mode=1,2,6)`
378	`if ( (mode.eq.1).or.(mode.eq.2) ) then`	415	`if ( (mode.eq.1).or.(mode.eq.2).or.(mode.eq.6) ) then`
379		416
380	`c Read the time specification (old and new format)`	417	`c Read the time specification (old and new format)`
381	`read(fid,'(a80)') linestr`	418	`read(fid,'(a80)') linestr`
382		419
383	`if ( linestr(1:14).eq.'Reference date' ) then`	420	`if ( linestr(1:14).eq.'Reference date' ) then`
384	`read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)')`	421	`read(linestr,'(a15,i4,i2,i2,a1,i2,i2,a13,i8,a4)')`
385	`> str1,`	422	`> str1,`
386	`> time(1),time(2),time(3),str2,time(4),time(5),`	423	`> time(1),time(2),time(3),str2,time(4),time(5),`
387	`> str3,time(6),str4`	424	`> str3,time(6),str4`
388		425
389	`elseif ( linestr(1:11).eq.'time period' ) then`	426	`elseif ( linestr(1:11).eq.'time period' ) then`
390	`read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)')`	427	`read(linestr,'(a12,i4,i2,i2,a1,i2,a4,i6,a1,i2,a2)')`
391	`> str1,`	428	`> str1,`
392	`> time(1),time(2),time(3),str2,time(4),`	429	`> time(1),time(2),time(3),str2,time(4),`
393	`> str3,itmp1,str3,itmp2,str4`	430	`> str3,itmp1,str3,itmp2,str4`
394	`time(5) = 0`	431	`time(5) = 0`
395	`time(6) = itmp1 * 60 + itmp2`	432	`time(6) = itmp1 * 60 + itmp2`
396		433
397	`endif`	434	`endif`
398		435
399	`c Skip the empty line and read field names`	436	`c Skip the empty line and read field names`
400	`read(fid,*)`	437	`read(fid,*)`
401	`read(fid,'(a500)',end=100) str`	438	`read(fid,'(a500)',end=100) str`
402	`do i=1,500`	439	`do i=1,500`
403	`ch(i)=str(i:i)`	440	`ch(i)=str(i:i)`
404	`enddo`	441	`enddo`
405		442
406	`c Split the input string`	443	`c Split the input string`
407	`isstr=0`	444	`isstr=0`
408	`nvars=0`	445	`nvars=0`
409	`do i=1,500`	446	`do i=1,500`
410	`if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then`	447	`if ( (isstr.eq.0).and.(ch(i).ne.' ') ) then`
411	`isstr=1`	448	`isstr=1`
412	`ileft=i`	449	`ileft=i`
413	`elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then`	450	`elseif ( (isstr.eq.1).and.(ch(i).eq.' ') ) then`
414	`isstr=0`	451	`isstr=0`
415	`iright=i-1`	452	`iright=i-1`
416	`nvars=nvars+1`	453	`nvars=nvars+1`
417	`vars(nvars)=str(ileft:iright)`	454	`vars(nvars)=str(ileft:iright)`
418	`endif`	455	`endif`
419	`enddo`	456	`enddo`
420		457
421	`c Skip the empty line`	458	`c Skip the empty line`
422	`read(fid,*,end=100)`	459	`read(fid,*,end=100)`
423		460
424	`c Read fortran mode (mode=3)`	461	`c Read fortran mode (mode=3)`
425	`elseif (mode.eq.3) then`	462	`elseif (mode.eq.3) then`
426	`read(fid) ntra,ntim,ncol`	463	`read(fid) ntra,ntim,ncol`
427	`read(fid) time`	464	`read(fid) time`
428	`read(fid) vars`	465	`read(fid) vars`
429		466
430	`c Read IVE netcdf mode (mode=4)`	467	`c Read IVE netcdf mode (mode=4)`
431	`elseif (mode.eq.4) then`	468	`elseif (mode.eq.4) then`
432	`call getvars(fid,nvars,vars,ierr)`	469	`call getvars(fid,nvars,vars,ierr)`
433	`varname='BASEDATE'`	470	`varname='BASEDATE'`
434	`call getdat(fid,varname,0.,0,rtime,ierr)`	471	`call getdat(fid,varname,0.,0,rtime,ierr)`
435	`do i=1,6`	472	`do i=1,6`
436	`time(i)=nint(rtime(i))`	473	`time(i)=nint(rtime(i))`
437	`enddo`	474	`enddo`
438	`endif`	475	`endif`
439		476
440	`return`	477	`return`
441		478
442	`c End of file has been reached`	479	`c End of file has been reached`
443	`100 fid=-fid`	480	`100 fid=-fid`
444	`return`	481	`return`
445		482
446	`c Excetion handling`	483	`c Excetion handling`
447	`110 print*,'<read_hea>: Unexspected time format.... Stop'`	484	`110 print*,'<read_hea>: Unexspected time format.... Stop'`
448	`stop`	485	`stop`
449		486
450	`end`	487	`end`
451		488
452		489
453	`c ----------------------------------------------------------------`	490	`c ----------------------------------------------------------------`
454	`c Write header to trajectory file (in ascii mode)`	491	`c Write header to trajectory file (in ascii mode)`
455	`c ----------------------------------------------------------------`	492	`c ----------------------------------------------------------------`
456		493
457	`subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)`	494	`subroutine write_hea(fid,time,vars,ntra,ntim,ncol,mode)`
458		495
459	`implicit none`	496	`implicit none`
460		497
461	`c Declaration of subroutine parameters`	498	`c Declaration of subroutine parameters`
462	`integer fid`	499	`integer fid`
463	`integer time(6)`	500	`integer time(6)`

Subversion Repositories lagranto.ecmwf

(root)/trunk/lib/iotra.f – Rev 15 → 39