WebSVN – lagranto.ecmwf – 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,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                                    pollon,pollat
      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
 
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.'rotated') 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
      close(10)
 
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 )
           if ( gridtype.eq.'rotated' ) then
             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  )
           endif
           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 rotated grid
      if ( gridtype.eq.'rotated' ) 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.'rotated' ) then
         print*,'Grid   clon,clat     : ',clon,clat
         print*,'       nlonlat       : ',nlonlat
         print*,'       dlonlat       : ',dlonlat
      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     Bring trajectories into grid interval
      print*,xmin,ymin,dx,dy,nx,ny
      if ( ( abs(xmin+180.).lt.eps).and.
     >     ( abs(ymin+90. ).lt.eps).and.
     >     ( abs(dx-1.    ).lt.eps).and.
     >     ( abs(dy-1.    ).lt.eps).and.
     >     ( nx.eq.360    ).and.
     >     ( ny.eq.180    ) )
     >then
        print*,'HALLO'
        do i=1,ntra
          do j=1,ntime
            if ( trainp(i,j,2).lt.-180. ) then
                trainp(i,j,2) = trainp(i,j,2) + 360.
            endif
            if ( trainp(i,j,2).gt.180. ) then
                trainp(i,j,2) = trainp(i,j,2) - 360.
            endif
         enddo
        enddo
      endif
 
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.'ppos' ).and.(varsinp(4).ne.'p'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      varsinp(1) = 'TIME'
      varsinp(2) = 'lon'
      varsinp(3) = 'lat'
      varsinp(4) = 'p'
 
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     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.'rotated') 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.'rotated') 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)
      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)
      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)
      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)
 
         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=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)
      if (h.eq.0.) then
         print*,'bad xa input in splint'
         stop
      endif
      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)
 
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
 
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) ) 
      if ( gridtype.eq.'rotated' ) then
         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 )
      endif
      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 )
 
 
 
 
 
 
 
 
      
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  
      do i = 1,nx
         longitude(i) = xmin + real(i-1) * dx
      enddo
      do i = 1,ny
         latitude(i)  = ymin + real(i-1) * dy
      enddo
      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
      print*,trim(varname),' defined on ',trim(cdfname)
 
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
 
c     ------------------------------------------------------------------
c     Compute Cos/Sin of an argument in Degree instead of Radian
c     ------------------------------------------------------------------
 
      real function cosd(arg)
 
      real,intent(IN) :: arg
      real,parameter :: grad2rad=3.1415926/180.
      cosd=cos(arg*grad2rad)
      return
      end
 
      real function sind(arg)
 
      real,intent(IN) :: arg
      real,parameter :: grad2rad=3.1415926/180.
      sind=sin(arg*grad2rad)
      return
      end
 
 

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

Subversion Repositories lagranto.ecmwf

(root)/trunk/density/density.f – Rev 29 → 39