WebSVN – lagranto.icon – Diff – /trunk/density/density.f


      PROGRAM density
 
      use netcdf
 
      implicit none
 
c     ---------------------------------------------------------------------
c     Declaration of variables
c     ---------------------------------------------------------------------
      
c     Parameter and working arrays
      real                                    radius
      character*80                            runit
      integer                                 nx,ny
      integer                                 nlonlat
      real                                    dlonlat
      real                                    xmin,ymin,dx,dy
      real                                    clon,clat
      integer                                 ntime,nfield,ntra
      character*80                            inpfile
      character*80                            outfile
      character*80                            mode
      real                                    param
      integer                                 opts,npts
      integer                                 step
      character*80                            gridtype
      character*80                            field
      integer                                 crefile,crevar
      real                                    pollon,pollat,polgam
      real,allocatable,    dimension (:,:) :: cnt,res,fld,area
      real,allocatable,    dimension (:)   :: traj
      real,allocatable,    dimension (:)   :: olon,olat,otim,ofld
      real,allocatable,    dimension (:)   :: nlon,nlat,ntim,nfld
 
c     Output format
      character*80                            outformat
 
c     Physical and mathematical constants
      real                                    pi180
      parameter                               (pi180=3.14159/180.)
      real                                    deltay
      parameter                               (deltay=111.)
      real                                    eps
      parameter                               (eps=0.001)
 
c     Input trajectories (see iotra.f)
      integer                                 inpmode
      real,allocatable, dimension (:,:,:) ::  trainp     
      integer                                 reftime(6)      
      character*80                            varsinp(100)   
      integer,allocatable, dimension (:) ::   sel_flag
      character*80                            sel_file
      character*80                            sel_format
 
c     Auxiliary variables
      character*80                            cdfname,varname
      integer                                 i,j,k
      integer                                 stat
      integer,allocatable, dimension (:,:) :: connect0
      integer                                 connectval0
      integer,allocatable, dimension (:,:) :: connect1
      integer                                 connectval1
      integer,allocatable, dimension (:,:) :: connect2
      integer                                 connectval2
      real                                    slat
      integer                                 ipre
      real                                    addvalue
      real                                    xmax,ymax
      real ,allocatable, dimension (:)  ::    odist,ndist
      real                                    dt
      integer                                 fid
      integer                                 dynamic_grid
      real                                    ycen,xcen
      integer                                 indx,indy
      character*80                            unit
      real                                    rlon0,rlat0,rlon,rlat
      real                                    lon,lat
      real                                    crot
      integer                                 count
      character*80                            longname, varunit
      real                                    time
      integer                                 ind
      integer                                 ifield
      real                                    hhmm,frac
      integer                                 ierr,ncID
      integer                                 rotflag
 
c     External functions
      real         lmstolm,lmtolms
      real         phstoph,phtophs
      external     lmstolm,lmtolms,phstoph,phtophs
      
      real         sdis
      external     sdis
 
c     ---------------------------------------------------------------------
c     Preparations
c     ---------------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM DENSITY ***'
      print*
 
c     Read input parameters
      open(10,file='density.param')
       read(10,*) inpfile
       read(10,*) outfile
       read(10,*) field
       read(10,*) ntime,nfield,ntra
       read(10,*) gridtype
       if ( gridtype.eq.'latlon' ) then 
          read(10,*) nx,ny,xmin,ymin,dx,dy
       elseif ( gridtype.eq.'centered') then
          read(10,*) clon,clat,nlonlat,dlonlat
       else
          print*,' ERROR: unsupported grid type ',trim(gridtype)
          stop
       endif
       read(10,*) radius,runit
       read(10,*) mode
       read(10,*) param
       read(10,*) step
       read(10,*) sel_file
       read(10,*) sel_format
       read(10,*) crefile
       read(10,*) crevar
       read(10,*) pollon,pollat
      close(10)
 
c	  Set the flag for trajectory rotation
      if ( abs(pollat-90.).gt.eps ) then
          rotflag = 1
      else
          rotflag = 0
      endif
 
c     Get the grid parameters if <crefile=0>
      if ( crefile.eq.0 ) then
 
           ierr = nf90_open  (trim(outfile), NF90_NOWRITE  , ncID)
 
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'grid'   ,gridtype ) 
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon'   ,clon     )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat'   ,clat     )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat  )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat  )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx'     ,nx       )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny'     ,ny       )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx'     ,dx       )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy'     ,dy       )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin'   ,xmin     )
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin'   ,ymin     )
 
           ierr = nf90_close(ncID)
 
           print*,'**** GRID PARAMETERS IMPORTED ',
     >            'FROM NETCDF FILE!!!! ****'
           print*
           
      endif
 
c     Check for consistency
      if ( (step.ne.0).and.(mode.ne.'keep') ) then
         print*," ERROR: interpolation is only possible for all",
     >                   ' time steps... Stop'
         stop
      endif
 
c     Set the number of times (just code aesthetics)
      opts=ntime
 
c     Set grid parameters for centered grid
      if ( gridtype.eq.'centered' ) then
         nx   = nlonlat
         ny   = nlonlat
         dx   = dlonlat
         dy   = dlonlat
         xmin = - real(nlonlat-1)/2. * dx
         xmax = + real(nlonlat-1)/2. * dx
         ymin = - real(nlonlat-1)/2. * dy
         ymax = + real(nlonlat-1)/2. * dy
      endif
      
c     Set the flag for dynamic grid adjustment
      if ( (nx.eq.0).or.(ny.eq.0) ) then
         dynamic_grid = 1
      else
         dynamic_grid = 0
      endif
 
c     Print status information
      print*,'---- INPUT PARAMETERS -----------------------------------'
      print* 
      print*,'Input                : ',trim(inpfile)
      print*,'Output               : ',trim(outfile)
      print*,'Field                : ',trim(field)
      print*,'Trajectory           : ',ntime,nfield,ntra
      print*,'Grid type            : ',trim(gridtype)
      if ( dynamic_grid.eq.1 ) then
         print*,'Grid                 : dynamic (see below)'
      elseif ( gridtype.eq.'latlon' ) then
         print*,'Grid   nlon,nlat     : ',nx,ny
         print*,'       lonmin,latmin : ',xmin,ymin
         print*,'       dlon,dlat     : ',dx,dy
      elseif ( gridtype.eq.'centered' ) then
         print*,'Grid   clon,clat     : ',clon,clat
         print*,'       nlonlat       : ',nlonlat
         print*,'       dlonlat       : ',dlonlat
      endif
      if ( rotflag.eq.0 ) then
          print*,'Pole lon/lat         : ',pollon,pollat
      else
          print*,'Pole lon/lat         : ',pollon,pollat,' (ROTATION)'
      endif
      print*,'Filter radius        : ',radius,' ',trim(runit)
      print*,'Mode                 : ',trim(mode)
      if ( ( mode.eq.'time'  ).or.
     >     ( mode.eq.'space' ).or.
     >     ( mode.eq.'grid' ) )
     >then
         print*,'Parameter            : ',param
      endif
      if ( step.eq.0 ) then
         print*,'Time step            : all'
      elseif (step.gt.0) then
         print*,'Time step            : ',step
      endif
      print*,'Selection file       : ',trim(sel_file)
      print*,'Selection format     : ',trim(sel_file)
      print*,'Flag <crefile>       : ',crefile
      print*,'Flag <crevar>        : ',crevar
 
c     Check whether mode is valid
      if ((mode.ne.'keep'  ).and.
     >    (mode.ne.'time'  ).and.
     >    (mode.ne.'space' ).and.
     >    (mode.ne.'grid'  ))  
     >then
         print*,' ERROR: Invalid mode ',trim(mode)
         stop
      endif
 
c     Allocate memory for old and new (reparameterised) trajectory
      allocate(olon(ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array olon ***'
      allocate(olat(ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array olat ***'
      allocate(otim(ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array otim ***'
      allocate(nlon(1000*ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array nlon ***'
      allocate(nlat(1000*ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array nlat ***'
      allocate(ntim(1000*ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ntim ***'
      allocate(odist(ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array odist ***'
      allocate(ndist(1000*ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ndist ***'
      allocate(ofld(ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ofld ***'
      allocate(nfld(1000*ntime),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array nfld ***'
 
c     Allocate memory for complete trajectory set
      allocate(trainp(ntra,ntime,nfield),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp ***'
      allocate(sel_flag(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array sel_flag ***'
 
c     Allocate memory for auxiliary fields
      allocate(traj(nfield),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array traj ***'
 
c     Set the format of the input file
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1
 
c     Read the input trajectory file
      call ropen_tra(fid,inpfile,ntra,ntime,nfield,
     >                   reftime,varsinp,inpmode)
      call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)
      call close_tra(fid,inpmode)
 
c     Check that first four columns correspond to time,lon,lat,p
      if ( (varsinp(1).ne.'time' ).or.
     >     (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
     >     (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
     >     (varsinp(4).ne.'zpos' ).and.(varsinp(4).ne.'z'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      varsinp(1) = 'TIME'
      varsinp(2) = 'lon'
      varsinp(3) = 'lat'
      varsinp(4) = 'z'
 
c     Get the index of the field (if needed)
      if ( field.ne.'nil' ) then
         ifield = 0
         do i=1,nfield
            if ( varsinp(i).eq.field ) ifield = i
         enddo
         if ( ifield.eq.0 ) then
            print*,' ERROR: field ',trim(field),' not found... Stop'
            stop
         endif
      endif
 
c     Rotate coordinates - if pollat<>90
      do i=1,ntra
        do j=1,ntime
 
            lon = trainp(i,j,2)
            lat = trainp(i,j,3)
 
            if ( rotflag.eq.1 ) then
               rlon = lmtolms(lat,lon,pollat,pollon)
               rlat = phtophs(lat,lon,pollat,pollon)
            else
               rlon = lon
               rlat = lat
            endif
 
            trainp(i,j,2) = rlon
            trainp(i,j,3) = rlat
 
        enddo
      enddo
 
c     Write some status information of the input trajectories
      print*
      print*,'---- INPUT TRAJECTORIES ---------------------------------'
      print*
      print*,' Reference time (year)  : ',reftime(1)
      print*,'                (month) : ',reftime(2)
      print*,'                (day)   : ',reftime(3)
      print*,'                (hour)  : ',reftime(4)
      print*,'                (min)   : ',reftime(5)
      print*,' Time range (min)       : ',reftime(6)
      do i=1,nfield
         if ( i.ne.ifield ) then
            print*,' Var                    :',i,trim(varsinp(i))
         else
            print*,' Var                    :',i,trim(varsinp(i)),
     >                                        '       [ gridding ]'
         endif
      enddo
      print*,' List of selected times'
      do i=1,ntime
         if ( (step.eq.0).or.(step.eq.i) ) then
            print*,'     ',i,'  -> ',trainp(1,i,1)
         endif
      enddo
      print*
 
c     Select flag: all trajectories are selected
      if ( sel_file.eq.'nil' ) then
 
         do i=1,ntra
            sel_flag(i) = 1
         enddo
 
c     Select flag: index file
      elseif ( sel_format.eq.'index' ) then
 
         do i=1,ntra
            sel_flag(i) = 0
         enddo
         
         open(10,file=sel_file)
 142      read(10,*,end=141) ind
          sel_flag(ind) = 1
          goto 142 
 141     continue
         close(10)
 
c     Select flag: boolean file
      elseif ( sel_format.eq.'boolean' ) then
       
         open(10,file=sel_file)
          do i=1,ntra
            read(10,*) ind
            if ( ind.eq.1 ) sel_flag(i) = ind
          enddo
         close(10)
         
      endif
 
c     Write status information
      if ( sel_file.eq.'nil' ) then
          print*,' Selected trajectories  : all ',ntra          
       else
          count = 0
          do i=1,ntra
             if ( sel_flag(i).eq.1 ) count = count + 1
          enddo
          print*,' #selected trajectories : ',count,
     >            ' [ ',real(count)/real(ntra) * 100.,' % ] '
       endif
       print*
 
c     ---------------------------------------------------------------------
c     Coordinate transformations and grid adjustment
c     ---------------------------------------------------------------------
 
c     Transform from lat/lon to rotated lat/lon, if requested
      if ( gridtype.eq.'centered') then
 
         crot = 0.
 
         pollon=clon-180.
         if (pollon.lt.-180.) pollon=pollon+360.
         pollat=90.-clat
         do i=1,ntra
            do j=1,ntime
               
               if ( sel_flag(i).eq.1 ) then
 
c                Get lat/lon coordinates for trajectory point
                 lon = trainp(i,j,2)
                 lat = trainp(i,j,3)
 
c                First Rotation
                 pollon=clon-180.
                 if (pollon.lt.-180.) pollon=pollon+360.
                 pollat=90.-clat
                 rlon0=lmtolms(lat,lon,pollat,pollon)
                 rlat0=phtophs(lat,lon,pollat,pollon)            
 
c                Second rotation
                 pollon=-180.
                 pollat=90.+crot
                 rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)
                 rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)   
 
c                Get rotated latitude and longitude
 100             if (rlon.lt.xmin) then
                  rlon=rlon+360.
                  goto 100
                 endif
 102             if (rlon.gt.(xmin+real(nx-1)*dx)) then
                  rlon=rlon-360.
                  goto 102
                 endif
 
c                Set the new trajectory coordinates
                 trainp(i,j,2) = rlon
                 trainp(i,j,3) = rlat
 
              endif
 
            enddo
         enddo
            
      endif
 
c     Dynamic grid adjustment
      if ( dynamic_grid.eq.1 ) then
 
c        Get the grid parameters
         xmin =  180.
         ymin =   90.
         xmax = -180.
         ymax =  -90.
 
         do i=1,ntra
 
            if ( sel_flag(i).eq.1 ) then
 
              if ( step.eq.0 ) then
               do j=1,ntime
                  if ( trainp(i,j,2).lt.xmin) xmin =  trainp(i,j,2)
                  if ( trainp(i,j,2).gt.xmax) xmax =  trainp(i,j,2)
                  if ( trainp(i,j,3).lt.ymin) ymin =  trainp(i,j,3)
                  if ( trainp(i,j,3).gt.ymax) ymax =  trainp(i,j,3)
               enddo
              else
                if ( trainp(i,step,2).lt.xmin) xmin =  trainp(i,step,2)
                if ( trainp(i,step,2).gt.xmax) xmax =  trainp(i,step,2)
                if ( trainp(i,step,3).lt.ymin) ymin =  trainp(i,step,3)
                if ( trainp(i,step,3).gt.ymax) ymax =  trainp(i,step,3)
              endif
            
            endif
 
         enddo
 
c        Get first guess for "optimal" grid
         nx = 400
         ny = 400
         dx = (xmax - xmin)/real(nx-1)
         dy = (ymax - ymin)/real(ny-1)
 
c        Make the grid spacing equal in zonal and meridional direction
         if ( dx.gt.dy ) then
            
            dy = dx
            ny = (ymax - ymin)/dy + 1
            if (ny.lt.nx/2)              ny = nx / 2
            if ( real(ny)*dy .ge. 180. ) ny = 180./dy + 1
            ycen = 0.5* (ymin+ymax)
            ymin = ycen - 0.5 * real(ny/2) * dy
            if (ymin.le.-90.) ymin = -90.
 
         else
               
            dx = dy
            nx = (xmax - xmin)/dx + 1
            if (nx.lt.ny/2)              nx = ny / 2
            if ( real(nx)*dx .ge. 360. ) nx = 360./dx + 1
            xcen = 0.5* (xmin+xmax)
            xmin = xcen - 0.5 * real(nx/2) * dx
            if (xmin.le.-180.) xmin = -180.
 
         endif
            
c        Write information
         print*
         print*,'---- DYNAMIC GRID ADJUSTMENT',
     >          ' ----------------------------'  
         print*
         print*,'Grid   nlon,nlat     : ',nx,ny
         print*,'       lonmin,latmin : ',xmin,ymin
         print*,'       dlon,dlat     : ',dx,dy
         print*
 
c     Write grid information for rotated grid (if not already done
      elseif ( gridtype.eq.'centered') then
         
         print*
         print*,'---- GRID PARAMETERS -------',
     >          ' ----------------------------'  
         print*
         print*,'Grid   nlon,nlat     : ',nx,ny
         print*,'       lonmin,latmin : ',xmin,ymin
         print*,'       dlon,dlat     : ',dx,dy
         print*
       
 
      endif
 
c     Set the grid boundaries
      xmax=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy
 
c     Allocate memory for output array and auxiliary gridding array 
      allocate(cnt(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array cnt  ***'
      allocate(res(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array res  ***'
      allocate(fld(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array fld  ***'
      allocate(area(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array area ***'
 
      allocate(connect0(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect0 ***'
      allocate(connect1(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect1 ***'
      allocate(connect2(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect2 ***'
 
 
c     Init the output array
      do i=1,nx
         do j=1,ny
            connect0(i,j) = 0
            connect1(i,j) = 0
            connect2(i,j) = 0
            cnt(i,j)      = 0.
            res(i,j)      = 0.
            fld(i,j)      = 0.
         enddo
      enddo  
 
c     ---------------------------------------------------------------------
c     Gridding
c     ---------------------------------------------------------------------
 
c     Write some status information 
      print*,'---- GRIDDING -------------------------------------------'
      print*
 
c     Loop over all entries of sampling table
      connectval0 = 0
      connectval1 = 0
      connectval2 = 0
      count       = 0
 
      do i=1,ntra
 
         if (mod(i,100).eq.0) print*,i,' of ',ntra
 
c        Skip all trajectories which are not selected
         if ( sel_flag(i).eq.0 ) goto 300
 
c        ------- Read a complete trajectory ---------------------------
         do j=1,ntime
            otim(j) = trainp(i,j,1)
            olon(j) = trainp(i,j,2)
            olat(j) = trainp(i,j,3)
            if ( field.ne.'nil' ) then
               ofld(j) =trainp(i,j,ifield)
            endif
         enddo
 
c        -------- Convert hh.m time into fractional time --------------
         do j=1,ntime
            hhmm    = otim(j)
            call hhmm2frac (hhmm,frac)
            otim(j) = frac
         enddo
 
c        -------- Interpolation ---------------------------------------
 
c        Keep the trajectory points as they are
         if ( ( mode.eq.'keep').and.(step.eq.0) ) then
            npts=opts
            do j=1,opts
               ntim(j)=otim(j)
               nlon(j)=olon(j)
               nlat(j)=olat(j)
               if ( field.ne.'nil' ) then
                  nfld(j)=ofld(j)
               endif
            enddo
 
c        Select a single time step
         elseif ( ( mode.eq.'keep').and.(step.gt.0) ) then
            npts    = 1
            ntim(1) = otim(step)
            nlon(1) = olon(step)
            nlat(1) = olat(step)
            if ( field.ne.'nil' ) then
               nfld(1) = ofld(step)
            endif
 
c        Perform a reparameterisation in time
         else if ( (mode.eq.'time').and.(step.eq.0) ) then
 
c           Get the new number of trajectory points
            npts=nint(abs(otim(opts)-otim(1))/param)+1
 
c           Handle date line problem
            do j=2,opts
               if ( (olon(j-1)-olon(j)).gt.180. ) then
                  olon(j) = olon(j) + 360.
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
                  olon(j) = olon(j) - 360.
               endif
            enddo
 
c           Cubic spline fitting
            call curvefit(otim,olon,opts,ntim,nlon,npts)
            call curvefit(otim,olat,opts,ntim,nlat,npts)
            if ( field.ne.'nil' ) then
               call curvefit(otim,ofld,opts,ntim,nfld,npts)
            endif
 
c           Reverse date line handling
            do j=1,npts
               if ( nlon(j).gt.xmax ) then
                  nlon(j) = nlon(j) -360.
               else if ( nlon(j).lt.xmin ) then
                  nlon(j) = nlon(j) +360.
               endif
            enddo
 
c        Perform a reparameterisation with equally spaced gridpoint
         elseif ( (mode.eq.'space').and.(step.eq.0) ) then
            
c           Calculate the distance and spacing
            odist(1) = 0.
            unit     = 'km'
            do j=2,ntime
               odist(j)=odist(j-1) + 
     >                  sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
            enddo
            
c           Determine the new number of trajectory points
            npts=nint(odist(ntime)/param)+1
            if (npts.eq.0) then
               npts=1.
            endif
        
c           Handle date line problem
            do j=2,opts
               if ( (olon(j-1)-olon(j)).gt.180. ) then
                  olon(j) = olon(j) + 360.
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
                  olon(j) = olon(j) - 360.
               endif
            enddo
                  
c           Cubic spline fitting
            call curvefit(odist,olon,opts,ndist,nlon,npts)
            call curvefit(odist,olat,opts,ndist,nlat,npts)
            call curvefit(odist,otim,opts,ndist,ntim,npts)
            if ( field.ne.'nil' ) then
               call curvefit(odist,ofld,opts,ndist,nfld,npts)
            endif
 
c           Reverse date line handling
            do j=1,npts
               if ( nlon(j).gt.xmax ) then
                  nlon(j) = nlon(j) -360.
               else if ( nlon(j).lt.xmin ) then
                  nlon(j) = nlon(j) +360.
               endif
            enddo
 
c        Perform a reparameterisation with equally spaced gridpoint
         elseif ( (mode.eq.'grid').and.(step.eq.0) ) then
            
c           Calculate the distance and spacing
            odist(1) = 0.
            unit     = 'deg'
            do j=2,ntime
               odist(j)=odist(j-1) + 
     >                  sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
            enddo
            
c           Determine the new number of trajectory points
            npts=nint(odist(ntime)/param)+1
            if (npts.eq.0) then
               npts=1.
            endif
        
c           Handle date line problem
            do j=2,opts
               if ( (olon(j-1)-olon(j)).gt.180. ) then
                  olon(j) = olon(j) + 360.
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
                  olon(j) = olon(j) - 360.
               endif
            enddo
                  
c           Cubic spline fitting
            call curvefit(odist,olon,opts,ndist,nlon,npts)
            call curvefit(odist,olat,opts,ndist,nlat,npts)
            call curvefit(odist,otim,opts,ndist,ntim,npts)
            if ( field.ne.'nil' ) then
               call curvefit(odist,ofld,opts,ndist,nfld,npts)
            endif
 
c           Reverse date line handling
            do j=1,npts
               if ( nlon(j).gt.xmax ) then
                  nlon(j) = nlon(j) -360.
               else if ( nlon(j).lt.xmin ) then
                  nlon(j) = nlon(j) +360.
               endif
            enddo
 
         endif
 
c        -------- Do the gridding -------------------------------------
 
c        Gridding of trajectory
         do j=1,npts
 
c           Check whether point is in data domain
	    if ( (nlon(j).gt.xmin).and.(nlon(j).lt.xmax).and.
     >           (nlat(j).gt.ymin).and.(nlat(j).lt.ymax))
     >      then
 
c              Increase counter for gridded points
               count = count + 1
 
c              ----------------- Gridding: simple count -----------------
               connectval0 = connectval0+1
               
               addvalue    = 1.
               
               call  gridding1
     >              (nlat(j),nlon(j),addvalue,
     >               radius,runit,connect0,connectval0,
     >               cnt,nx,ny,xmin,ymin,dx,dy)
 
c              ----------------- Gridding: residence time ---------------
               connectval1 = connectval1+1
               
               if ( ntime.eq.1 ) then
                  addvalue = 0.
               elseif ( j.eq.1 )  then
                  addvalue=abs(ntim(2)-ntim(1))
               else
                  addvalue=abs(ntim(j)-ntim(j-1))
               endif
               
               call  gridding1
     >              (nlat(j),nlon(j),addvalue,
     >               radius,runit,connect1,connectval1,
     >               res,nx,ny,xmin,ymin,dx,dy)
 
 
c              --------------- Gridding: field -------------------------
               if ( field.ne.'nil' ) then
 
                   connectval2 = connectval2+1
               
                   addvalue    = nfld(j)
               
                   call  gridding1
     >                  (nlat(j),nlon(j),addvalue,
     >                  radius,runit,connect2,connectval2,
     >                  fld,nx,ny,xmin,ymin,dx,dy)
 
               endif
 
	    endif
 
         enddo
 
c        Exit point for loop over all trajectories
 300     continue
 
      enddo
 
c     Write status information
      print*
      print*,' # gridded points       : ',count
 
c     ---------------------------------------------------------------------
c     Unit conversions and output to netCDF file
c     ---------------------------------------------------------------------
 
c     Write some status information 
      print*
      print*,'---- WRITE OUTPUT ---------------------------------------'
      print*
 
c     Area (in km^2)
      do i=1,nx	         
         do j=1,ny	
            slat=ymin+real(j-1)*dy
            if (abs(abs(slat)-90.).gt.eps) then
               area(i,j) = dy*dx*cos(pi180*slat)*deltay**2
            else
               area(i,j) = 0.
            endif
         enddo
      enddo
 
c     Normalise gridded field
      if ( field.ne.'nil' ) then
         do i=1,nx
            do j=1,ny
               if ( cnt(i,j).gt.0. ) then
                  fld(i,j) = fld(i,j) / cnt(i,j)
               endif
            enddo
         enddo
      endif
 
c     Set the time for the output netCDF files - if a composite is
c     calculatd, then the time is set to 
      if ( step.eq.0 ) then
         time = -999.
         print*,'   ... COMPOSITE OVER ALL TRAJECTORY TIMES (-999)'
         print*
      else
         time = trainp(1,step,1)
      endif
 
c     Write output to CF netCDF
      cdfname  = outfile
      
      varname  = 'COUNT'
      longname = 'trajectory counts'
      varunit  = 'counts per grid point'
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
     >       clon,clat,nlonlat,dlonlat,cnt,time,dx,dy,xmin,ymin,nx,
     >       ny,crefile,crefile,1,pollon,pollat)
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
     >     ' [ time = ',time,' ]'    
 
      varname  = 'RESIDENCE'
      longname = 'residence time'
      varunit  = 'hours per grid point'
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
     >       clon,clat,nlonlat,dlonlat,res,time,dx,dy,xmin,ymin,nx,
     >       ny,0,crefile,1,pollon,pollat)
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
     >     ' [ time = ',time,' ]'    
 
      varname  = 'AREA'
      longname = 'area corresponding to grid points'
      varunit  = 'square kilometers'
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
     >       clon,clat,nlonlat,dlonlat,area,time,dx,dy,xmin,ymin,nx,
     >       ny,0,crefile,1,pollon,pollat)
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
     >     ' [ time = ',time,' ]'    
      
      if ( field.ne.'nil' ) then
         varname  = field
         longname = field
         varunit  = 'as on trajectory file'
         call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
     >       clon,clat,nlonlat,dlonlat,fld,time,dx,dy,xmin,ymin,nx,
     >       ny,0,crevar,1,pollon,pollat)
 
         write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
     >        '    ... ',trim(varname),' -> ',trim(cdfname),
     >        ' [ time = ',time,' ]'    
      endif
 
c     Write status information
      print*
      print*,'              *** END OF PROGRAM DENSITY **'
      print*,'========================================================='
 
      end
 
c     ********************************************************************
c     * GRIDDING SUBROUTINES                                             *
c     ********************************************************************
 
c     ---------------------------------------------------------------------
c     Gridding of one single data point (smoothing in km, deg, gridp)
c     ---------------------------------------------------------------------
 
      subroutine gridding1 (lat,lon,addval,radius,unit,
     >                      connect,connectval,
     >                      out,nx,ny,xmin,ymin,dx,dy)
 
      implicit none
 
c     Declaration of subroutine parameters
      real         lat,lon
      integer      nx,ny
      real         xmin,ymin,dx,dy
      real         out(nx,ny)
      real         radius
      character*80 unit
      integer      connectval
      integer      connect(nx,ny)
      real         addval
 
c     Auxiliary variables
      integer   i,j,k
      integer   mu,md,nr,nl,n,m
      integer   stackx(nx*ny),stacky(nx*ny)
      integer   tab_x(nx*ny),tab_y(nx*ny)
      real      tab_r(nx*ny)
      integer   sp
      real      lat2,lon2
      real      dist,sum
      real      xmax
      integer   periodic
      integer   test
 
c     Numerical epsilon
      real      eps
      parameter (eps=0.01)
 
c     Externals
      real      sdis,weight
      external  sdis,weight
 
c     Check whether lat/lon point is valid
      xmax=xmin+real(nx-1)*dx
      if (lon.lt.xmin-eps) lon=lon+360.
      if (lon.gt.xmax+eps) lon=lon-360.
      if (abs(lat-90).lt.eps) lat=90.
      if (abs(lat+90).lt.eps) lat=-90.
      if ((abs(lat).gt.(90.+eps)).or.
     >    (lon.lt.xmin-eps).or.(lon.gt.xmax+eps)) then
         print*,'Invalid lat/lon point ',lat,lon
         return
      endif
 
c     Set flag for periodic domain
      if (abs(xmax-xmin-360.).lt.eps) then
         periodic=1
      else if (abs(xmax-xmin-360+dx).lt.eps) then
         periodic=2
      else
         periodic=0
      endif
 
c     Get indices of one coarse grid point within search radius
      i=nint((lon-xmin)/dx)+1
      if ((i.eq.nx).and.(periodic.eq.1)) i=1
      j=nint((lat-ymin)/dy)+1
      lat2=ymin+real(j-1)*dy
      lon2=xmin+real(i-1)*dx
      dist=sdis(lon,lat,lon2,lat2,unit)
      if (dist.gt.radius) then
         print*,'1: Search radius is too small...'
         stop
      endif
 
c     Get connected points
      k=0
      stackx(1)=i
      stacky(1)=j
      sp=1
      do while (sp.ne.0) 
         
c        Get an element from stack
         n=stackx(sp)
         m=stacky(sp)
         sp=sp-1
                  
c        Get distance from reference point
         lat2=ymin+real(m-1)*dy
         lon2=xmin+real(n-1)*dx
         dist=sdis(lon,lat,lon2,lat2,unit)
 
c        Check whether distance is smaller than search radius: connected
         if (dist.lt.radius) then
 
c           Make entry in filter mask
            k=k+1
            tab_x(k)=n
            tab_y(k)=m
            tab_r(k)=weight(dist,radius)
 
c           Mark this point as visited
            connect(n,m)=connectval
                     
c           Get coordinates of neighbouring points
            nr=n+1
            if ((nr.gt.nx)  .and.(periodic.eq.0)) nr=nx
            if ((nr.gt.nx-1).and.(periodic.eq.1)) nr=1
            if ((nr.gt.nx)  .and.(periodic.eq.2)) nr=1
            nl=n-1
            if ((nl.lt.1).and.(periodic.eq.0)) nl=1
            if ((nl.lt.1).and.(periodic.eq.1)) nl=nx-1
            if ((nl.lt.1).and.(periodic.eq.2)) nl=nx
            mu=m+1
            if (mu.gt.ny) mu=ny
            md=m-1
            if (md.lt.1) md=1
 
c           Update stack
            if (connect(nr,m).ne.connectval) then
               connect(nr,m)=connectval
               sp=sp+1
               stackx(sp)=nr
               stacky(sp)=m
            endif
            if (connect(nl,m).ne.connectval) then
               connect(nl,m)=connectval
               sp=sp+1
               stackx(sp)=nl
               stacky(sp)=m
            endif
            if (connect(n,mu).ne.connectval) then
               connect(n,mu)=connectval
               sp=sp+1
               stackx(sp)=n
               stacky(sp)=mu
            endif
            if (connect(n,md).ne.connectval) then
               connect(n,md)=connectval
               sp=sp+1
               stackx(sp)=n
               stacky(sp)=md
            endif
         endif
         
      end do
 
      if (k.ge.1) then
         sum=0.
         do i=1,k
            sum=sum+tab_r(i)
         enddo
         do i=1,k
            out(tab_x(i),tab_y(i))=out(tab_x(i),tab_y(i))+
     >                             addval*tab_r(i)/sum
 
            if ((tab_x(i).eq.1).and.(periodic.eq.1)) then
               out(nx,tab_y(i))=out(nx,tab_y(i))+
     >                             addval*tab_r(i)/sum
            endif
         enddo
      else
         print*,'2: Search radius is too small...'
         stop
      endif
 
      end
 
 
c     ----------------------------------------------------------------------
c     Get spherical distance between lat/lon points
c     ----------------------------------------------------------------------
            
      real function sdis(xp,yp,xq,yq,unit)
 
c     Calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
 
      real         re
      parameter    (re=6370.)
      real         xp,yp,xq,yq,arg
      character*80 unit
      real         dlon
 
      if ( unit.eq.'km' ) then
 
         arg=sin(yp*pi180)*sin(yq*pi180)+cos(yp*pi180)*cos(yq*pi180)*
     >       cos((xp-xq)*pi180)
c         arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
         if (arg.lt.-1.) arg=-1.
         if (arg.gt.1.) arg=1.
         sdis=re*acos(arg)
         
      elseif ( unit.eq.'deg' ) then
 
         dlon = xp-xq
         if ( dlon.gt. 180. ) dlon = dlon - 360.
         if ( dlon.lt.-180. ) dlon = dlon + 360.
         sdis = sqrt( dlon**2 + (yp-yq)**2 )
 
      endif
      
 
c     Quick and dirty trick to avoid zero distances
      if (sdis.eq.0.) sdis=0.1
 
      end
 
c     ----------------------------------------------------------------------
c     Weight function for the filter mask
c     ----------------------------------------------------------------------
 
      real function weight (r,radius)
 
c     Attribute to each distanc r its corresponding weight in the filter mask
 
      implicit none
 
c     Declaration of subroutine parameters
      real r
      real radius
 
c     Simple 0/1 mask
      if (r.lt.radius) then
         weight=exp(-r/radius)
      else
         weight=0.
      endif
 
      end
 
 
c     ********************************************************************
c     * REPARAMETERIZATION SUBROUTINES                                   *
c     ********************************************************************
 
c     -------------------------------------------------------------
c     Interpolation of the trajectory with a natural cubic spline
c     -------------------------------------------------------------
 
      SUBROUTINE curvefit (time,lon,n,
     >                     sptime,splon,spn)
 
c     Given the curve <time,lon> with <n> data points, fit a
c     cubic spline to this curve. The new curve is returned in 
c     <sptime,splon,spn> with <spn> data points. The parameter
c     <spn> specifies on entry the number of spline interpolated points
c     along the curve.
      
      implicit none
 
c     Declaration of subroutine parameters
      integer n
      real time(n),lon(n)
      integer spn
      real sptime(spn),splon(spn)
 
c     Auxiliary variables
      real y2ax(n)
      real dt
      real s
      integer i
      real order
 
c     Determine whether the input array is ascending or descending
      if (time(1).gt.time(n)) then
         order=-1.
      else
         order= 1.
      endif
 
c     Bring the time array into ascending order
      do i=1,n
         time(i)=order*time(i)
      enddo
 
c     Prepare the (natural) cubic spline interpolation
      call spline (time,lon,n,1.e30,1.e30,y2ax)
      dt=(time(n)-time(1))/real(spn-1)
      do i=1,spn
         sptime(i)=time(1)+real(i-1)*dt
      enddo
      
c     Do the spline interpolation
      do i=1,spn
         call splint(time,lon,y2ax,n,sptime(i),s)
         splon(i)=s
      enddo
 
c     Change the time arrays back
      do i=1,spn
         sptime(i)=order*sptime(i)
      enddo
      do i=1,n
         time(i)=order*time(i)
      enddo
 
      return
      end
 
c     -------------------------------------------------------------
c     Basic routines for spline interpolation (Numerical Recipes)
c     -------------------------------------------------------------
 
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
 
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
c      if (h.eq.0.) pause 'bad xa input in splint'
      if (h.eq.0.) then
        write(*,*) 'bad xa input in splint'
        read(*,*)
      end if
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END
 
c     ********************************************************************
c     * INPUT / OUTPUT SUBROUTINES                                       *
c     ********************************************************************
 
 
c     --------------------------------------------------------------------
c     Subroutines to write the CF netcdf output file
c     --------------------------------------------------------------------
 
      subroutine writecdf2D_cf 
     >          (cdfname,varname,longname,unit,gridtype,clon,clat,
     >           nlonlat,dlonlat,arr,time,dx,dy,xmin,ymin,nx,ny,
     >           crefile,crevar,cretime,pollon,pollat)
 
c     Create and write to the CF netcdf file <cdfname>. The variable
c     with name <varname> and with time <time> is written. The data
c     are in the two-dimensional array <arr>. The list <dx,dy,xmin,
c     ymin,nx,ny> specifies the output grid. The flags <crefile> and
c     <crevar> determine whether the file and/or the variable should
c     be created; correspondingly for the unlimited dimension <time>
c     with the flag <cretime>.
 
      USE netcdf
 
      IMPLICIT NONE
 
c     Declaration of input parameters
      character*80 cdfname
      character*80 varname,longname,unit
      integer      nx,ny
      real         arr(nx,ny)
      real         dx,dy,xmin,ymin
      real         time
      integer      crefile,crevar,cretime
      character*80 gridtype
      real         clon,clat
      integer      nlonlat
      real         dlonlat
      real         pollon,pollat
 
c     Local variables
      integer      ierr
      integer      ncID
      integer      LonDimId,    varLonID
      integer      LatDimID,    varLatID
      integer      TimeDimID,   varTimeID
      real         longitude(nx)
      real         latitude (ny)
      real         timeindex
      integer      i
      integer      nvars,varids(100)
      integer      ndims,dimids(100)
      real         timelist(1000)
      integer      ntimes
      integer      ind
      integer      varID
 
c     Quick an dirty solution for fieldname conflict
      if ( varname.eq.'time' ) varname = 'TIME'
 
c     Initially set error to indicate no errors.
      ierr = 0
 
c     ---- Create the netCDF - skip if <crefile=0> ----------------------
      if ( crefile.ne.1 ) goto 100
 
c     Create the file 
      ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
      
c     Define dimensions 
      ierr=nf90_def_dim(ncID,'longitude',nx            , LonDimID )
      ierr=nf90_def_dim(ncID,'latitude' ,ny            , LatDimID )
      ierr=nf90_def_dim(ncID,'time'     ,nf90_unlimited, TimeDimID)
      
c     Define coordinate Variables 
      ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
     >     (/ LonDimID /),varLonID)
      ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
      ierr = nf90_put_att(ncID, varLonID, "units"      ,"degree_east")
      
      ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
     >     (/ LatDimID /),varLatID)
      ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
      ierr = nf90_put_att(ncID, varLatID, "units"    ,"degree_north")
      
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT, 
     >     (/ TimeDimID /), varTimeID)
      ierr = nf90_put_att(ncID, varTimeID, "axis",            "T")
      ierr = nf90_put_att(ncID, varTimeID, "calendar", "standard")
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
      
c     Write global attributes 
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',  
     >     'Trajectory Densities')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source', 
     >     'Lagranto Trajectories')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution', 
     >     'ETH Zurich, IACETH')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'grid',trim(gridtype) ) 
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clon',clon )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clat',clat )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dx',dx )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dy',dy )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'pollon',pollon )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'pollat',pollat )
 
c     Write coordinate data 
      do i = 1,nx+1
         longitude(i) = xmin + real(i-1) * dx
      enddo
      do i = 1,ny+1
         latitude(i)  = ymin + real(i-1) * dy
      enddo
      
c     Check whether the definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ', 
     >        'finish definition mode.'
         stop
      endif
      
c     Write coordinate data  
      ierr = nf90_put_var(ncID,varLonID ,longitude)
      ierr = nf90_put_var(ncID,varLatID ,latitude )
      
c     Close netCDF file 
      ierr = nf90_close(ncID)
      
 100  continue
 
c     ---- Define a new variable - skip if <crevar=0> -----------------------
 
      if ( crevar.ne.1 ) goto 110
      
c     Open the file for read(write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
 
c     Get the IDs for dimensions
      ierr = nf90_inq_dimid(ncID,'longitude', LonDimID )
      ierr = nf90_inq_dimid(ncID,'latitude' , LatDimID )
      ierr = nf90_inq_dimid(ncID,'time'     , TimeDimID)
 
c     Enter define mode
      ierr = nf90_redef(ncID)
 
c     Write definition and add attributes
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
     >                   (/ LonDimID, LatDimID, TimeDimID /),varID)
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
      ierr = nf90_put_att(ncID, varID, "units"     , unit     ) 
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  ) 
 
c     Check whether definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ', 
     >           'finish definition mode.'
         stop
      endif
 
c     Close netCDF file 
      ierr = nf90_close(ncID)
 
 110  continue
 
c     ---- Create a new time (unlimited dimension) - skip if <cretime=0> ------
 
      if ( cretime.ne.1 ) goto 120
 
c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE, ncID)
      
c     Get the list of times on the netCDF file
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      if ( ierr.ne.0 ) then
         print*,'Time dimension is not defined on ',trim(cdfname),
     >          ' .... Stop'
         stop
      endif
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      if ( ierr.ne.0 ) then
         print*,'Variable time is not defined on ',trim(cdfname),
     >          ' ... Stop'
         stop
      endif
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
 
c     Decide whether a new time must be written
      ind = 0
      do i=1,ntimes
         if ( time.eq.timelist(i) ) ind = i
      enddo
 
c     Extend the time list if required 
      if ( ind.eq.0 ) then
         ntimes           = ntimes + 1
         timelist(ntimes) = time
         ierr = nf90_put_var(ncID,varTimeID,timelist(1:ntimes))
      endif
 
c     Close netCDF file 
      ierr = nf90_close(ncID)
 
 120  continue
 
c     ---- Write data --------------------------------------------------
 
c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)
 
c     Get the varID
      ierr = nf90_inq_varid(ncID,varname, varID )
      if (ierr.ne.0) then
         print*,'Variable ',trim(varname),' is not defined on ',
     >          trim(cdfname)
         stop
      endif
 
c     Get the time index
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      if ( ierr.ne.0 ) then
         print*,'Time dimension is not defined on ',trim(cdfname),
     >          ' .... Stop'
         stop
      endif
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      if ( ierr.ne.0 ) then
         print*,'Variable time is not defined on ',trim(cdfname),
     >          ' ... Stop'
         stop
      endif
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
      ind = 0
      do i=1,ntimes
         if ( time.eq.timelist(i) ) ind = i
      enddo
      if (ind.eq.0) then
         print*,'Time',time,' is not defined on the netCDF file',
     >          trim(cdfname),' ... Stop'
         stop
      endif
 
c     Write data block      
      ierr = nf90_put_var(ncID,varID,arr,
     >                    start = (/ 1, 1, ind /), 
     >                    count = (/ nx, ny, 1 /) )
 
c     Check whether writing was successful 
      ierr = nf90_close(ncID)
      if (ierr.ne.0) then
         write(*,*) trim(nf90_strerror(ierr))
         write(*,*) 'An error occurred while attempting to ', 
     >              'close the netcdf file.'
         write(*,*) 'in clscdf_CF'
      endif
 
      end
 
 
c     ********************************************************************************
c     * Transformation routine: LMSTOLM and PHSTOPH from library gm2em               *
c     ********************************************************************************
 
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHIS   = ZPIR18*PHIS
      ZLAMS   = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS   = ZPIR18*ZLAMS
 
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) -
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) +
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMSTOLM =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMSTOLAM =  90.0
            ELSE
              LMSTOLAM = -90.0
            ENDIF
      ELSE
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END
 
 
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      SINPOL = SIN(ZPIR18*POLPHI)
      COSPOL = COS(ZPIR18*POLPHI)
      ZPHIS  = ZPIR18*PHIS
      ZLAMS  = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS  = ZPIR18*ZLAMS
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
 
      PHSTOPH = ZRPI18*ASIN(ARG)
 
      RETURN
      END
 
 
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL =     ZPIR18*POLLAM
      ZPHI    =     ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
 
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMTOLMS =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMTOLMS =  90.0
            ELSE
              LMTOLMS = -90.0
            ENDIF
      ELSE
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END
 
 
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHI    = ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 
      RETURN
      END

Rev 3	Rev 5
1	`PROGRAM density`	1	`PROGRAM density`
2		2
3	`use netcdf`	3	`use netcdf`
4		4
5	`implicit none`	5	`implicit none`
6		6
7	`c ---------------------------------------------------------------------`	7	`c ---------------------------------------------------------------------`
8	`c Declaration of variables`	8	`c Declaration of variables`
9	`c ---------------------------------------------------------------------`	9	`c ---------------------------------------------------------------------`
10		10
11	`c Parameter and working arrays`	11	`c Parameter and working arrays`
12	`real radius`	12	`real radius`
13	`character*80 runit`	13	`character*80 runit`
14	`integer nx,ny`	14	`integer nx,ny`
15	`integer nlonlat`	15	`integer nlonlat`
16	`real dlonlat`	16	`real dlonlat`
17	`real xmin,ymin,dx,dy`	17	`real xmin,ymin,dx,dy`
18	`real clon,clat`	18	`real clon,clat`
19	`integer ntime,nfield,ntra`	19	`integer ntime,nfield,ntra`
20	`character*80 inpfile`	20	`character*80 inpfile`
21	`character*80 outfile`	21	`character*80 outfile`
22	`character*80 mode`	22	`character*80 mode`
23	`real param`	23	`real param`
24	`integer opts,npts`	24	`integer opts,npts`
25	`integer step`	25	`integer step`
26	`character*80 gridtype`	26	`character*80 gridtype`
27	`character*80 field`	27	`character*80 field`
28	`integer crefile,crevar`	28	`integer crefile,crevar`
29	`real pollon,pollat,polgam`	29	`real pollon,pollat,polgam`
30	`real,allocatable, dimension (:,:) :: cnt,res,fld,area`	30	`real,allocatable, dimension (:,:) :: cnt,res,fld,area`
31	`real,allocatable, dimension (:) :: traj`	31	`real,allocatable, dimension (:) :: traj`
32	`real,allocatable, dimension (:) :: olon,olat,otim,ofld`	32	`real,allocatable, dimension (:) :: olon,olat,otim,ofld`
33	`real,allocatable, dimension (:) :: nlon,nlat,ntim,nfld`	33	`real,allocatable, dimension (:) :: nlon,nlat,ntim,nfld`
34		34
35	`c Output format`	35	`c Output format`
36	`character*80 outformat`	36	`character*80 outformat`
37		37
38	`c Physical and mathematical constants`	38	`c Physical and mathematical constants`
39	`real pi180`	39	`real pi180`
40	`parameter (pi180=3.14159/180.)`	40	`parameter (pi180=3.14159/180.)`
41	`real deltay`	41	`real deltay`
42	`parameter (deltay=111.)`	42	`parameter (deltay=111.)`
43	`real eps`	43	`real eps`
44	`parameter (eps=0.001)`	44	`parameter (eps=0.001)`
45		45
46	`c Input trajectories (see iotra.f)`	46	`c Input trajectories (see iotra.f)`
47	`integer inpmode`	47	`integer inpmode`
48	`real,allocatable, dimension (:,:,:) :: trainp`	48	`real,allocatable, dimension (:,:,:) :: trainp`
49	`integer reftime(6)`	49	`integer reftime(6)`
50	`character*80 varsinp(100)`	50	`character*80 varsinp(100)`
51	`integer,allocatable, dimension (:) :: sel_flag`	51	`integer,allocatable, dimension (:) :: sel_flag`
52	`character*80 sel_file`	52	`character*80 sel_file`
53	`character*80 sel_format`	53	`character*80 sel_format`
54		54
55	`c Auxiliary variables`	55	`c Auxiliary variables`
56	`character*80 cdfname,varname`	56	`character*80 cdfname,varname`
57	`integer i,j,k`	57	`integer i,j,k`
58	`integer stat`	58	`integer stat`
59	`integer,allocatable, dimension (:,:) :: connect0`	59	`integer,allocatable, dimension (:,:) :: connect0`
60	`integer connectval0`	60	`integer connectval0`
61	`integer,allocatable, dimension (:,:) :: connect1`	61	`integer,allocatable, dimension (:,:) :: connect1`
62	`integer connectval1`	62	`integer connectval1`
63	`integer,allocatable, dimension (:,:) :: connect2`	63	`integer,allocatable, dimension (:,:) :: connect2`
64	`integer connectval2`	64	`integer connectval2`
65	`real slat`	65	`real slat`
66	`integer ipre`	66	`integer ipre`
67	`real addvalue`	67	`real addvalue`
68	`real xmax,ymax`	68	`real xmax,ymax`
69	`real ,allocatable, dimension (:) :: odist,ndist`	69	`real ,allocatable, dimension (:) :: odist,ndist`
70	`real dt`	70	`real dt`
71	`integer fid`	71	`integer fid`
72	`integer dynamic_grid`	72	`integer dynamic_grid`
73	`real ycen,xcen`	73	`real ycen,xcen`
74	`integer indx,indy`	74	`integer indx,indy`
75	`character*80 unit`	75	`character*80 unit`
76	`real rlon0,rlat0,rlon,rlat`	76	`real rlon0,rlat0,rlon,rlat`
77	`real lon,lat`	77	`real lon,lat`
78	`real crot`	78	`real crot`
79	`integer count`	79	`integer count`
80	`character*80 longname, varunit`	80	`character*80 longname, varunit`
81	`real time`	81	`real time`
82	`integer ind`	82	`integer ind`
83	`integer ifield`	83	`integer ifield`
84	`real hhmm,frac`	84	`real hhmm,frac`
85	`integer ierr,ncID`	85	`integer ierr,ncID`
86	`integer rotflag`	86	`integer rotflag`
87		87
88	`c External functions`	88	`c External functions`
89	`real lmstolm,lmtolms`	89	`real lmstolm,lmtolms`
90	`real phstoph,phtophs`	90	`real phstoph,phtophs`
91	`external lmstolm,lmtolms,phstoph,phtophs`	91	`external lmstolm,lmtolms,phstoph,phtophs`
92		92
93	`real sdis`	93	`real sdis`
94	`external sdis`	94	`external sdis`
95		95
96	`c ---------------------------------------------------------------------`	96	`c ---------------------------------------------------------------------`
97	`c Preparations`	97	`c Preparations`
98	`c ---------------------------------------------------------------------`	98	`c ---------------------------------------------------------------------`
99		99
100	`c Write start message`	100	`c Write start message`
101	`print*,'========================================================='`	101	`print*,'========================================================='`
102	`print,' START OF PROGRAM DENSITY *'`	102	`print,' START OF PROGRAM DENSITY *'`
103	`print*`	103	`print*`
104		104
105	`c Read input parameters`	105	`c Read input parameters`
106	`open(10,file='density.param')`	106	`open(10,file='density.param')`
107	`read(10,*) inpfile`	107	`read(10,*) inpfile`
108	`read(10,*) outfile`	108	`read(10,*) outfile`
109	`read(10,*) field`	109	`read(10,*) field`
110	`read(10,*) ntime,nfield,ntra`	110	`read(10,*) ntime,nfield,ntra`
111	`read(10,*) gridtype`	111	`read(10,*) gridtype`
112	`if ( gridtype.eq.'latlon' ) then`	112	`if ( gridtype.eq.'latlon' ) then`
113	`read(10,*) nx,ny,xmin,ymin,dx,dy`	113	`read(10,*) nx,ny,xmin,ymin,dx,dy`
114	`elseif ( gridtype.eq.'centered') then`	114	`elseif ( gridtype.eq.'centered') then`
115	`read(10,*) clon,clat,nlonlat,dlonlat`	115	`read(10,*) clon,clat,nlonlat,dlonlat`
116	`else`	116	`else`
117	`print*,' ERROR: unsupported grid type ',trim(gridtype)`	117	`print*,' ERROR: unsupported grid type ',trim(gridtype)`
118	`stop`	118	`stop`
119	`endif`	119	`endif`
120	`read(10,*) radius,runit`	120	`read(10,*) radius,runit`
121	`read(10,*) mode`	121	`read(10,*) mode`
122	`read(10,*) param`	122	`read(10,*) param`
123	`read(10,*) step`	123	`read(10,*) step`
124	`read(10,*) sel_file`	124	`read(10,*) sel_file`
125	`read(10,*) sel_format`	125	`read(10,*) sel_format`
126	`read(10,*) crefile`	126	`read(10,*) crefile`
127	`read(10,*) crevar`	127	`read(10,*) crevar`
128	`read(10,*) pollon,pollat`	128	`read(10,*) pollon,pollat`
129	`close(10)`	129	`close(10)`
130		130
131	`c Set the flag for trajectory rotation`	131	`c Set the flag for trajectory rotation`
132	`if ( abs(pollat-90.).gt.eps ) then`	132	`if ( abs(pollat-90.).gt.eps ) then`
133	`rotflag = 1`	133	`rotflag = 1`
134	`else`	134	`else`
135	`rotflag = 0`	135	`rotflag = 0`
136	`endif`	136	`endif`
137		137
138	`c Get the grid parameters if <crefile=0>`	138	`c Get the grid parameters if <crefile=0>`
139	`if ( crefile.eq.0 ) then`	139	`if ( crefile.eq.0 ) then`
140		140
141	`ierr = nf90_open (trim(outfile), NF90_NOWRITE , ncID)`	141	`ierr = nf90_open (trim(outfile), NF90_NOWRITE , ncID)`
142		142
143	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'grid' ,gridtype )`	143	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'grid' ,gridtype )`
144	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon' ,clon )`	144	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon' ,clon )`
145	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat' ,clat )`	145	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat' ,clat )`
146	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )`	146	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )`
147	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )`	147	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )`
148	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx' ,nx )`	148	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx' ,nx )`
149	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny' ,ny )`	149	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny' ,ny )`
150	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx' ,dx )`	150	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx' ,dx )`
151	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy' ,dy )`	151	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy' ,dy )`
152	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin' ,xmin )`	152	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin' ,xmin )`
153	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin' ,ymin )`	153	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin' ,ymin )`
154		154
155	`ierr = nf90_close(ncID)`	155	`ierr = nf90_close(ncID)`
156		156
157	`print,'*** GRID PARAMETERS IMPORTED ',`	157	`print,'*** GRID PARAMETERS IMPORTED ',`
158	`> 'FROM NETCDF FILE!!!! ****'`	158	`> 'FROM NETCDF FILE!!!! ****'`
159	`print*`	159	`print*`
160		160
161	`endif`	161	`endif`
162		162
163	`c Check for consistency`	163	`c Check for consistency`
164	`if ( (step.ne.0).and.(mode.ne.'keep') ) then`	164	`if ( (step.ne.0).and.(mode.ne.'keep') ) then`
165	`print*," ERROR: interpolation is only possible for all",`	165	`print*," ERROR: interpolation is only possible for all",`
166	`> ' time steps... Stop'`	166	`> ' time steps... Stop'`
167	`stop`	167	`stop`
168	`endif`	168	`endif`
169		169
170	`c Set the number of times (just code aesthetics)`	170	`c Set the number of times (just code aesthetics)`
171	`opts=ntime`	171	`opts=ntime`
172		172
173	`c Set grid parameters for centered grid`	173	`c Set grid parameters for centered grid`
174	`if ( gridtype.eq.'centered' ) then`	174	`if ( gridtype.eq.'centered' ) then`
175	`nx = nlonlat`	175	`nx = nlonlat`
176	`ny = nlonlat`	176	`ny = nlonlat`
177	`dx = dlonlat`	177	`dx = dlonlat`
178	`dy = dlonlat`	178	`dy = dlonlat`
179	`xmin = - real(nlonlat-1)/2. * dx`	179	`xmin = - real(nlonlat-1)/2. * dx`
180	`xmax = + real(nlonlat-1)/2. * dx`	180	`xmax = + real(nlonlat-1)/2. * dx`
181	`ymin = - real(nlonlat-1)/2. * dy`	181	`ymin = - real(nlonlat-1)/2. * dy`
182	`ymax = + real(nlonlat-1)/2. * dy`	182	`ymax = + real(nlonlat-1)/2. * dy`
183	`endif`	183	`endif`
184		184
185	`c Set the flag for dynamic grid adjustment`	185	`c Set the flag for dynamic grid adjustment`
186	`if ( (nx.eq.0).or.(ny.eq.0) ) then`	186	`if ( (nx.eq.0).or.(ny.eq.0) ) then`
187	`dynamic_grid = 1`	187	`dynamic_grid = 1`
188	`else`	188	`else`
189	`dynamic_grid = 0`	189	`dynamic_grid = 0`
190	`endif`	190	`endif`
191		191
192	`c Print status information`	192	`c Print status information`
193	`print*,'---- INPUT PARAMETERS -----------------------------------'`	193	`print*,'---- INPUT PARAMETERS -----------------------------------'`
194	`print*`	194	`print*`
195	`print*,'Input : ',trim(inpfile)`	195	`print*,'Input : ',trim(inpfile)`
196	`print*,'Output : ',trim(outfile)`	196	`print*,'Output : ',trim(outfile)`
197	`print*,'Field : ',trim(field)`	197	`print*,'Field : ',trim(field)`
198	`print*,'Trajectory : ',ntime,nfield,ntra`	198	`print*,'Trajectory : ',ntime,nfield,ntra`
199	`print*,'Grid type : ',trim(gridtype)`	199	`print*,'Grid type : ',trim(gridtype)`
200	`if ( dynamic_grid.eq.1 ) then`	200	`if ( dynamic_grid.eq.1 ) then`
201	`print*,'Grid : dynamic (see below)'`	201	`print*,'Grid : dynamic (see below)'`
202	`elseif ( gridtype.eq.'latlon' ) then`	202	`elseif ( gridtype.eq.'latlon' ) then`
203	`print*,'Grid nlon,nlat : ',nx,ny`	203	`print*,'Grid nlon,nlat : ',nx,ny`
204	`print*,' lonmin,latmin : ',xmin,ymin`	204	`print*,' lonmin,latmin : ',xmin,ymin`
205	`print*,' dlon,dlat : ',dx,dy`	205	`print*,' dlon,dlat : ',dx,dy`
206	`elseif ( gridtype.eq.'centered' ) then`	206	`elseif ( gridtype.eq.'centered' ) then`
207	`print*,'Grid clon,clat : ',clon,clat`	207	`print*,'Grid clon,clat : ',clon,clat`
208	`print*,' nlonlat : ',nlonlat`	208	`print*,' nlonlat : ',nlonlat`
209	`print*,' dlonlat : ',dlonlat`	209	`print*,' dlonlat : ',dlonlat`
210	`endif`	210	`endif`
211	`if ( rotflag.eq.0 ) then`	211	`if ( rotflag.eq.0 ) then`
212	`print*,'Pole lon/lat : ',pollon,pollat`	212	`print*,'Pole lon/lat : ',pollon,pollat`
213	`else`	213	`else`
214	`print*,'Pole lon/lat : ',pollon,pollat,' (ROTATION)'`	214	`print*,'Pole lon/lat : ',pollon,pollat,' (ROTATION)'`
215	`endif`	215	`endif`
216	`print*,'Filter radius : ',radius,' ',trim(runit)`	216	`print*,'Filter radius : ',radius,' ',trim(runit)`
217	`print*,'Mode : ',trim(mode)`	217	`print*,'Mode : ',trim(mode)`
218	`if ( ( mode.eq.'time' ).or.`	218	`if ( ( mode.eq.'time' ).or.`
219	`> ( mode.eq.'space' ).or.`	219	`> ( mode.eq.'space' ).or.`
220	`> ( mode.eq.'grid' ) )`	220	`> ( mode.eq.'grid' ) )`
221	`>then`	221	`>then`
222	`print*,'Parameter : ',param`	222	`print*,'Parameter : ',param`
223	`endif`	223	`endif`
224	`if ( step.eq.0 ) then`	224	`if ( step.eq.0 ) then`
225	`print*,'Time step : all'`	225	`print*,'Time step : all'`
226	`elseif (step.gt.0) then`	226	`elseif (step.gt.0) then`
227	`print*,'Time step : ',step`	227	`print*,'Time step : ',step`
228	`endif`	228	`endif`
229	`print*,'Selection file : ',trim(sel_file)`	229	`print*,'Selection file : ',trim(sel_file)`
230	`print*,'Selection format : ',trim(sel_file)`	230	`print*,'Selection format : ',trim(sel_file)`
231	`print*,'Flag <crefile> : ',crefile`	231	`print*,'Flag <crefile> : ',crefile`
232	`print*,'Flag <crevar> : ',crevar`	232	`print*,'Flag <crevar> : ',crevar`
233		233
234	`c Check whether mode is valid`	234	`c Check whether mode is valid`
235	`if ((mode.ne.'keep' ).and.`	235	`if ((mode.ne.'keep' ).and.`
236	`> (mode.ne.'time' ).and.`	236	`> (mode.ne.'time' ).and.`
237	`> (mode.ne.'space' ).and.`	237	`> (mode.ne.'space' ).and.`
238	`> (mode.ne.'grid' ))`	238	`> (mode.ne.'grid' ))`
239	`>then`	239	`>then`
240	`print*,' ERROR: Invalid mode ',trim(mode)`	240	`print*,' ERROR: Invalid mode ',trim(mode)`
241	`stop`	241	`stop`
242	`endif`	242	`endif`
243		243
244	`c Allocate memory for old and new (reparameterised) trajectory`	244	`c Allocate memory for old and new (reparameterised) trajectory`
245	`allocate(olon(ntime),stat=stat)`	245	`allocate(olon(ntime),stat=stat)`
246	`if (stat.ne.0) print,' error allocating array olon *'`	246	`if (stat.ne.0) print,' error allocating array olon *'`
247	`allocate(olat(ntime),stat=stat)`	247	`allocate(olat(ntime),stat=stat)`
248	`if (stat.ne.0) print,' error allocating array olat *'`	248	`if (stat.ne.0) print,' error allocating array olat *'`
249	`allocate(otim(ntime),stat=stat)`	249	`allocate(otim(ntime),stat=stat)`
250	`if (stat.ne.0) print,' error allocating array otim *'`	250	`if (stat.ne.0) print,' error allocating array otim *'`
251	`allocate(nlon(1000*ntime),stat=stat)`	251	`allocate(nlon(1000*ntime),stat=stat)`
252	`if (stat.ne.0) print,' error allocating array nlon *'`	252	`if (stat.ne.0) print,' error allocating array nlon *'`
253	`allocate(nlat(1000*ntime),stat=stat)`	253	`allocate(nlat(1000*ntime),stat=stat)`
254	`if (stat.ne.0) print,' error allocating array nlat *'`	254	`if (stat.ne.0) print,' error allocating array nlat *'`
255	`allocate(ntim(1000*ntime),stat=stat)`	255	`allocate(ntim(1000*ntime),stat=stat)`
256	`if (stat.ne.0) print,' error allocating array ntim *'`	256	`if (stat.ne.0) print,' error allocating array ntim *'`
257	`allocate(odist(ntime),stat=stat)`	257	`allocate(odist(ntime),stat=stat)`
258	`if (stat.ne.0) print,' error allocating array odist *'`	258	`if (stat.ne.0) print,' error allocating array odist *'`
259	`allocate(ndist(1000*ntime),stat=stat)`	259	`allocate(ndist(1000*ntime),stat=stat)`
260	`if (stat.ne.0) print,' error allocating array ndist *'`	260	`if (stat.ne.0) print,' error allocating array ndist *'`
261	`allocate(ofld(ntime),stat=stat)`	261	`allocate(ofld(ntime),stat=stat)`
262	`if (stat.ne.0) print,' error allocating array ofld *'`	262	`if (stat.ne.0) print,' error allocating array ofld *'`
263	`allocate(nfld(1000*ntime),stat=stat)`	263	`allocate(nfld(1000*ntime),stat=stat)`
264	`if (stat.ne.0) print,' error allocating array nfld *'`	264	`if (stat.ne.0) print,' error allocating array nfld *'`
265		265
266	`c Allocate memory for complete trajectory set`	266	`c Allocate memory for complete trajectory set`
267	`allocate(trainp(ntra,ntime,nfield),stat=stat)`	267	`allocate(trainp(ntra,ntime,nfield),stat=stat)`
268	`if (stat.ne.0) print,' error allocating array trainp *'`	268	`if (stat.ne.0) print,' error allocating array trainp *'`
269	`allocate(sel_flag(ntra),stat=stat)`	269	`allocate(sel_flag(ntra),stat=stat)`
270	`if (stat.ne.0) print,' error allocating array sel_flag *'`	270	`if (stat.ne.0) print,' error allocating array sel_flag *'`
271		271
272	`c Allocate memory for auxiliary fields`	272	`c Allocate memory for auxiliary fields`
273	`allocate(traj(nfield),stat=stat)`	273	`allocate(traj(nfield),stat=stat)`
274	`if (stat.ne.0) print,' error allocating array traj *'`	274	`if (stat.ne.0) print,' error allocating array traj *'`
275		275
276	`c Set the format of the input file`	276	`c Set the format of the input file`
277	`call mode_tra(inpmode,inpfile)`	277	`call mode_tra(inpmode,inpfile)`
278	`if (inpmode.eq.-1) inpmode=1`	278	`if (inpmode.eq.-1) inpmode=1`
279		279
280	`c Read the input trajectory file`	280	`c Read the input trajectory file`
281	`call ropen_tra(fid,inpfile,ntra,ntime,nfield,`	281	`call ropen_tra(fid,inpfile,ntra,ntime,nfield,`
282	`> reftime,varsinp,inpmode)`	282	`> reftime,varsinp,inpmode)`
283	`call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)`	283	`call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)`
284	`call close_tra(fid,inpmode)`	284	`call close_tra(fid,inpmode)`
285		285
286	`c Check that first four columns correspond to time,lon,lat,p`	286	`c Check that first four columns correspond to time,lon,lat,p`
287	`if ( (varsinp(1).ne.'time' ).or.`	287	`if ( (varsinp(1).ne.'time' ).or.`
288	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`	288	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`
289	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`	289	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`
290	`> (varsinp(4).ne.'zpos' ).and.(varsinp(4).ne.'z' ) )`	290	`> (varsinp(4).ne.'zpos' ).and.(varsinp(4).ne.'z' ) )`
291	`>then`	291	`>then`
292	`print*,' ERROR: problem with input trajectories ...'`	292	`print*,' ERROR: problem with input trajectories ...'`
293	`stop`	293	`stop`
294	`endif`	294	`endif`
295	`varsinp(1) = 'TIME'`	295	`varsinp(1) = 'TIME'`
296	`varsinp(2) = 'lon'`	296	`varsinp(2) = 'lon'`
297	`varsinp(3) = 'lat'`	297	`varsinp(3) = 'lat'`
298	`varsinp(4) = 'z'`	298	`varsinp(4) = 'z'`
299		299
300	`c Get the index of the field (if needed)`	300	`c Get the index of the field (if needed)`
301	`if ( field.ne.'nil' ) then`	301	`if ( field.ne.'nil' ) then`
302	`ifield = 0`	302	`ifield = 0`
303	`do i=1,nfield`	303	`do i=1,nfield`
304	`if ( varsinp(i).eq.field ) ifield = i`	304	`if ( varsinp(i).eq.field ) ifield = i`
305	`enddo`	305	`enddo`
306	`if ( ifield.eq.0 ) then`	306	`if ( ifield.eq.0 ) then`
307	`print*,' ERROR: field ',trim(field),' not found... Stop'`	307	`print*,' ERROR: field ',trim(field),' not found... Stop'`
308	`stop`	308	`stop`
309	`endif`	309	`endif`
310	`endif`	310	`endif`
311		311
312	`c Rotate coordinates - if pollat<>90`	312	`c Rotate coordinates - if pollat<>90`
313	`do i=1,ntra`	313	`do i=1,ntra`
314	`do j=1,ntime`	314	`do j=1,ntime`
315		315
316	`lon = trainp(i,j,2)`	316	`lon = trainp(i,j,2)`
317	`lat = trainp(i,j,3)`	317	`lat = trainp(i,j,3)`
318		318
319	`if ( rotflag.eq.1 ) then`	319	`if ( rotflag.eq.1 ) then`
320	`rlon = lmtolms(lat,lon,pollat,pollon)`	320	`rlon = lmtolms(lat,lon,pollat,pollon)`
321	`rlat = phtophs(lat,lon,pollat,pollon)`	321	`rlat = phtophs(lat,lon,pollat,pollon)`
322	`else`	322	`else`
323	`rlon = lon`	323	`rlon = lon`
324	`rlat = lat`	324	`rlat = lat`
325	`endif`	325	`endif`
326		326
327	`trainp(i,j,2) = rlon`	327	`trainp(i,j,2) = rlon`
328	`trainp(i,j,3) = rlat`	328	`trainp(i,j,3) = rlat`
329		329
330	`enddo`	330	`enddo`
331	`enddo`	331	`enddo`
332		332
333	`c Write some status information of the input trajectories`	333	`c Write some status information of the input trajectories`
334	`print*`	334	`print*`
335	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`	335	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`
336	`print*`	336	`print*`
337	`print*,' Reference time (year) : ',reftime(1)`	337	`print*,' Reference time (year) : ',reftime(1)`
338	`print*,' (month) : ',reftime(2)`	338	`print*,' (month) : ',reftime(2)`
339	`print*,' (day) : ',reftime(3)`	339	`print*,' (day) : ',reftime(3)`
340	`print*,' (hour) : ',reftime(4)`	340	`print*,' (hour) : ',reftime(4)`
341	`print*,' (min) : ',reftime(5)`	341	`print*,' (min) : ',reftime(5)`
342	`print*,' Time range (min) : ',reftime(6)`	342	`print*,' Time range (min) : ',reftime(6)`
343	`do i=1,nfield`	343	`do i=1,nfield`
344	`if ( i.ne.ifield ) then`	344	`if ( i.ne.ifield ) then`
345	`print*,' Var :',i,trim(varsinp(i))`	345	`print*,' Var :',i,trim(varsinp(i))`
346	`else`	346	`else`
347	`print*,' Var :',i,trim(varsinp(i)),`	347	`print*,' Var :',i,trim(varsinp(i)),`
348	`> ' [ gridding ]'`	348	`> ' [ gridding ]'`
349	`endif`	349	`endif`
350	`enddo`	350	`enddo`
351	`print*,' List of selected times'`	351	`print*,' List of selected times'`
352	`do i=1,ntime`	352	`do i=1,ntime`
353	`if ( (step.eq.0).or.(step.eq.i) ) then`	353	`if ( (step.eq.0).or.(step.eq.i) ) then`
354	`print*,' ',i,' -> ',trainp(1,i,1)`	354	`print*,' ',i,' -> ',trainp(1,i,1)`
355	`endif`	355	`endif`
356	`enddo`	356	`enddo`
357	`print*`	357	`print*`
358		358
359	`c Select flag: all trajectories are selected`	359	`c Select flag: all trajectories are selected`
360	`if ( sel_file.eq.'nil' ) then`	360	`if ( sel_file.eq.'nil' ) then`
361		361
362	`do i=1,ntra`	362	`do i=1,ntra`
363	`sel_flag(i) = 1`	363	`sel_flag(i) = 1`
364	`enddo`	364	`enddo`
365		365
366	`c Select flag: index file`	366	`c Select flag: index file`
367	`elseif ( sel_format.eq.'index' ) then`	367	`elseif ( sel_format.eq.'index' ) then`
368		368
369	`do i=1,ntra`	369	`do i=1,ntra`
370	`sel_flag(i) = 0`	370	`sel_flag(i) = 0`
371	`enddo`	371	`enddo`
372		372
373	`open(10,file=sel_file)`	373	`open(10,file=sel_file)`
374	`142 read(10,*,end=141) ind`	374	`142 read(10,*,end=141) ind`
375	`sel_flag(ind) = 1`	375	`sel_flag(ind) = 1`
376	`goto 142`	376	`goto 142`
377	`141 continue`	377	`141 continue`
378	`close(10)`	378	`close(10)`
379		379
380	`c Select flag: boolean file`	380	`c Select flag: boolean file`
381	`elseif ( sel_format.eq.'boolean' ) then`	381	`elseif ( sel_format.eq.'boolean' ) then`
382		382
383	`open(10,file=sel_file)`	383	`open(10,file=sel_file)`
384	`do i=1,ntra`	384	`do i=1,ntra`
385	`read(10,*) ind`	385	`read(10,*) ind`
386	`if ( ind.eq.1 ) sel_flag(i) = ind`	386	`if ( ind.eq.1 ) sel_flag(i) = ind`
387	`enddo`	387	`enddo`
388	`close(10)`	388	`close(10)`
389		389
390	`endif`	390	`endif`
391		391
392	`c Write status information`	392	`c Write status information`
393	`if ( sel_file.eq.'nil' ) then`	393	`if ( sel_file.eq.'nil' ) then`
394	`print*,' Selected trajectories : all ',ntra`	394	`print*,' Selected trajectories : all ',ntra`
395	`else`	395	`else`
396	`count = 0`	396	`count = 0`
397	`do i=1,ntra`	397	`do i=1,ntra`
398	`if ( sel_flag(i).eq.1 ) count = count + 1`	398	`if ( sel_flag(i).eq.1 ) count = count + 1`
399	`enddo`	399	`enddo`
400	`print*,' #selected trajectories : ',count,`	400	`print*,' #selected trajectories : ',count,`
401	`> ' [ ',real(count)/real(ntra) * 100.,' % ] '`	401	`> ' [ ',real(count)/real(ntra) * 100.,' % ] '`
402	`endif`	402	`endif`
403	`print*`	403	`print*`
404		404
405	`c ---------------------------------------------------------------------`	405	`c ---------------------------------------------------------------------`
406	`c Coordinate transformations and grid adjustment`	406	`c Coordinate transformations and grid adjustment`
407	`c ---------------------------------------------------------------------`	407	`c ---------------------------------------------------------------------`
408		408
409	`c Transform from lat/lon to rotated lat/lon, if requested`	409	`c Transform from lat/lon to rotated lat/lon, if requested`
410	`if ( gridtype.eq.'centered') then`	410	`if ( gridtype.eq.'centered') then`
411		411
412	`crot = 0.`	412	`crot = 0.`
413		413
414	`pollon=clon-180.`	414	`pollon=clon-180.`
415	`if (pollon.lt.-180.) pollon=pollon+360.`	415	`if (pollon.lt.-180.) pollon=pollon+360.`
416	`pollat=90.-clat`	416	`pollat=90.-clat`
417	`do i=1,ntra`	417	`do i=1,ntra`
418	`do j=1,ntime`	418	`do j=1,ntime`
419		419
420	`if ( sel_flag(i).eq.1 ) then`	420	`if ( sel_flag(i).eq.1 ) then`
421		421
422	`c Get lat/lon coordinates for trajectory point`	422	`c Get lat/lon coordinates for trajectory point`
423	`lon = trainp(i,j,2)`	423	`lon = trainp(i,j,2)`
424	`lat = trainp(i,j,3)`	424	`lat = trainp(i,j,3)`
425		425
426	`c First Rotation`	426	`c First Rotation`
427	`pollon=clon-180.`	427	`pollon=clon-180.`
428	`if (pollon.lt.-180.) pollon=pollon+360.`	428	`if (pollon.lt.-180.) pollon=pollon+360.`
429	`pollat=90.-clat`	429	`pollat=90.-clat`
430	`rlon0=lmtolms(lat,lon,pollat,pollon)`	430	`rlon0=lmtolms(lat,lon,pollat,pollon)`
431	`rlat0=phtophs(lat,lon,pollat,pollon)`	431	`rlat0=phtophs(lat,lon,pollat,pollon)`
432		432
433	`c Second rotation`	433	`c Second rotation`
434	`pollon=-180.`	434	`pollon=-180.`
435	`pollat=90.+crot`	435	`pollat=90.+crot`
436	`rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)`	436	`rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)`
437	`rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)`	437	`rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)`
438		438
439	`c Get rotated latitude and longitude`	439	`c Get rotated latitude and longitude`
440	`100 if (rlon.lt.xmin) then`	440	`100 if (rlon.lt.xmin) then`
441	`rlon=rlon+360.`	441	`rlon=rlon+360.`
442	`goto 100`	442	`goto 100`
443	`endif`	443	`endif`
444	`102 if (rlon.gt.(xmin+real(nx-1)*dx)) then`	444	`102 if (rlon.gt.(xmin+real(nx-1)*dx)) then`
445	`rlon=rlon-360.`	445	`rlon=rlon-360.`
446	`goto 102`	446	`goto 102`
447	`endif`	447	`endif`
448		448
449	`c Set the new trajectory coordinates`	449	`c Set the new trajectory coordinates`
450	`trainp(i,j,2) = rlon`	450	`trainp(i,j,2) = rlon`
451	`trainp(i,j,3) = rlat`	451	`trainp(i,j,3) = rlat`
452		452
453	`endif`	453	`endif`
454		454
455	`enddo`	455	`enddo`
456	`enddo`	456	`enddo`
457		457
458	`endif`	458	`endif`
459		459
460	`c Dynamic grid adjustment`	460	`c Dynamic grid adjustment`
461	`if ( dynamic_grid.eq.1 ) then`	461	`if ( dynamic_grid.eq.1 ) then`
462		462
463	`c Get the grid parameters`	463	`c Get the grid parameters`
464	`xmin = 180.`	464	`xmin = 180.`
465	`ymin = 90.`	465	`ymin = 90.`
466	`xmax = -180.`	466	`xmax = -180.`
467	`ymax = -90.`	467	`ymax = -90.`
468		468
469	`do i=1,ntra`	469	`do i=1,ntra`
470		470
471	`if ( sel_flag(i).eq.1 ) then`	471	`if ( sel_flag(i).eq.1 ) then`
472		472
473	`if ( step.eq.0 ) then`	473	`if ( step.eq.0 ) then`
474	`do j=1,ntime`	474	`do j=1,ntime`
475	`if ( trainp(i,j,2).lt.xmin) xmin = trainp(i,j,2)`	475	`if ( trainp(i,j,2).lt.xmin) xmin = trainp(i,j,2)`
476	`if ( trainp(i,j,2).gt.xmax) xmax = trainp(i,j,2)`	476	`if ( trainp(i,j,2).gt.xmax) xmax = trainp(i,j,2)`
477	`if ( trainp(i,j,3).lt.ymin) ymin = trainp(i,j,3)`	477	`if ( trainp(i,j,3).lt.ymin) ymin = trainp(i,j,3)`
478	`if ( trainp(i,j,3).gt.ymax) ymax = trainp(i,j,3)`	478	`if ( trainp(i,j,3).gt.ymax) ymax = trainp(i,j,3)`
479	`enddo`	479	`enddo`
480	`else`	480	`else`
481	`if ( trainp(i,step,2).lt.xmin) xmin = trainp(i,step,2)`	481	`if ( trainp(i,step,2).lt.xmin) xmin = trainp(i,step,2)`
482	`if ( trainp(i,step,2).gt.xmax) xmax = trainp(i,step,2)`	482	`if ( trainp(i,step,2).gt.xmax) xmax = trainp(i,step,2)`
483	`if ( trainp(i,step,3).lt.ymin) ymin = trainp(i,step,3)`	483	`if ( trainp(i,step,3).lt.ymin) ymin = trainp(i,step,3)`
484	`if ( trainp(i,step,3).gt.ymax) ymax = trainp(i,step,3)`	484	`if ( trainp(i,step,3).gt.ymax) ymax = trainp(i,step,3)`
485	`endif`	485	`endif`
486		486
487	`endif`	487	`endif`
488		488
489	`enddo`	489	`enddo`
490		490
491	`c Get first guess for "optimal" grid`	491	`c Get first guess for "optimal" grid`
492	`nx = 400`	492	`nx = 400`
493	`ny = 400`	493	`ny = 400`
494	`dx = (xmax - xmin)/real(nx-1)`	494	`dx = (xmax - xmin)/real(nx-1)`
495	`dy = (ymax - ymin)/real(ny-1)`	495	`dy = (ymax - ymin)/real(ny-1)`
496		496
497	`c Make the grid spacing equal in zonal and meridional direction`	497	`c Make the grid spacing equal in zonal and meridional direction`
498	`if ( dx.gt.dy ) then`	498	`if ( dx.gt.dy ) then`
499		499
500	`dy = dx`	500	`dy = dx`

Subversion Repositories lagranto.icon

(root)/trunk/density/density.f – Rev 3 → 5