WebSVN – lagranto.um – 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 ) 
           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 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     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         pi180
      parameter    (pi180=3.14159/180.)
      
      real         xp,yp,xq,yq,arg
      character*80 unit
      real         dlon
      
      if ( unit.eq.'km' ) then
 
         arg=sin(pi180*yp)*sin(pi180*yq)+
     >            cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(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.) pause 'bad xa input in splint'
      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)
      IF( ierr /= nf90_NoErr) PRINT *,1,NF90_STRERROR(ierr)
      
c     Define dimensions 
      ierr=nf90_def_dim(ncID,'longitude',nx            , LonDimID )
      IF( ierr /= nf90_NoErr) PRINT *,2,NF90_STRERROR(ierr)
      ierr=nf90_def_dim(ncID,'latitude' ,ny            , LatDimID )
      IF( ierr /= nf90_NoErr) PRINT *,3,NF90_STRERROR(ierr)
      ierr=nf90_def_dim(ncID,'time'     ,nf90_unlimited, TimeDimID)
      IF( ierr /= nf90_NoErr) PRINT *,4,NF90_STRERROR(ierr)
      
c     Define coordinate Variables 
      ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
     >     (/ LonDimID /),varLonID)
      IF( ierr /= nf90_NoErr) PRINT *,5,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
      IF( ierr /= nf90_NoErr) PRINT *,6,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varLonID, "units"      ,"degree_east")
      IF( ierr /= nf90_NoErr) PRINT *,7,NF90_STRERROR(ierr)
      
      ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
     >     (/ LatDimID /),varLatID)
      IF( ierr /= nf90_NoErr) PRINT *,8,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
      IF( ierr /= nf90_NoErr) PRINT *,9,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varLatID, "units"    ,"degree_north")
      IF( ierr /= nf90_NoErr) PRINT *,10,NF90_STRERROR(ierr)
      
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT, 
     >     (/ TimeDimID /), varTimeID)
      IF( ierr /= nf90_NoErr) PRINT *,11,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varTimeID, "axis",            "T")
      IF( ierr /= nf90_NoErr) PRINT *,12,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varTimeID, "calendar", "standard")
      IF( ierr /= nf90_NoErr) PRINT *,13,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
      IF( ierr /= nf90_NoErr) PRINT *,14,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
      IF( ierr /= nf90_NoErr) PRINT *,15,NF90_STRERROR(ierr)
      
c     Write global attributes 
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
      IF( ierr /= nf90_NoErr) PRINT *,16,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',  
     >     'Trajectory Densities')
      IF( ierr /= nf90_NoErr) PRINT *,17,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source', 
     >     'Lagranto Trajectories')
      IF( ierr /= nf90_NoErr) PRINT *,18,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution', 
     >     'ETH Zurich, IACETH')
 
c     Write grid description
      IF( ierr /= nf90_NoErr) PRINT *,19,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'grid',trim(gridtype) )
 
      if ( gridtype.eq.'rotated' ) then
        IF( ierr /= nf90_NoErr) PRINT *,20,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clon',clon )
        IF( ierr /= nf90_NoErr) PRINT *,21,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clat',clat )
        IF( ierr /= nf90_NoErr) PRINT *,22,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )
        IF( ierr /= nf90_NoErr) PRINT *,23,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )
        IF( ierr /= nf90_NoErr) PRINT *,24,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
        IF( ierr /= nf90_NoErr) PRINT *,25,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
        IF( ierr /= nf90_NoErr) PRINT *,26,NF90_STRERROR(ierr)
 
      elseif ( gridtype.eq.'latlon' ) then
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
        IF( ierr /= nf90_NoErr) PRINT *,27,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
        IF( ierr /= nf90_NoErr) PRINT *,28,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dx',dx )
        IF( ierr /= nf90_NoErr) PRINT *,29,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dy',dy )
        IF( ierr /= nf90_NoErr) PRINT *,30,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
        IF( ierr /= nf90_NoErr) PRINT *,31,NF90_STRERROR(ierr)
        ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
        IF( ierr /= nf90_NoErr) PRINT *,32,NF90_STRERROR(ierr)
      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
      
c     Check whether the definition was successful
      ierr = nf90_enddef(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
      
c     Write coordinate data  
      ierr = nf90_put_var(ncID,varLonID ,longitude)
      IF( ierr /= nf90_NoErr) PRINT *,33,NF90_STRERROR(ierr)
      ierr = nf90_put_var(ncID,varLatID ,latitude )
      IF( ierr /= nf90_NoErr) PRINT *,34,NF90_STRERROR(ierr)
      
c     Close netCDF file 
      ierr = nf90_close(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,35,NF90_STRERROR(ierr)
      
 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)
      IF( ierr /= nf90_NoErr) PRINT *,36,NF90_STRERROR(ierr)
 
c     Get the IDs for dimensions
      ierr = nf90_inq_dimid(ncID,'longitude', LonDimID )
      IF( ierr /= nf90_NoErr) PRINT *,37,NF90_STRERROR(ierr)
      ierr = nf90_inq_dimid(ncID,'latitude' , LatDimID )
      IF( ierr /= nf90_NoErr) PRINT *,38,NF90_STRERROR(ierr)
      ierr = nf90_inq_dimid(ncID,'time'     , TimeDimID)
      IF( ierr /= nf90_NoErr) PRINT *,39,NF90_STRERROR(ierr)
 
c     Enter define mode
      ierr = nf90_redef(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,40,NF90_STRERROR(ierr)
 
c     Write definition and add attributes
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
     >                   (/ LonDimID, LatDimID, TimeDimID /),varID)
      IF( ierr /= nf90_NoErr) PRINT *,41,NF90_STRERROR(ierr)
 
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
      IF( ierr /= nf90_NoErr) PRINT *,42,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varID, "units"     , unit     ) 
      IF( ierr /= nf90_NoErr) PRINT *,43,NF90_STRERROR(ierr)
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  ) 
      IF( ierr /= nf90_NoErr) PRINT *,44,NF90_STRERROR(ierr)
 
c     Check whether definition was successful
      ierr = nf90_enddef(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,45,NF90_STRERROR(ierr)
 
c     Close netCDF file 
      ierr = nf90_close(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,46,NF90_STRERROR(ierr)
 
 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)
      IF( ierr /= nf90_NoErr) PRINT *,47,NF90_STRERROR(ierr)
 
c     Get the list of times on the netCDF file
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      IF( ierr /= nf90_NoErr) PRINT *,48,NF90_STRERROR(ierr)
 
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      IF( ierr /= nf90_NoErr) PRINT *,49,NF90_STRERROR(ierr)
 
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      IF( ierr /= nf90_NoErr) PRINT *,50,NF90_STRERROR(ierr)
 
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
      IF( ierr /= nf90_NoErr) PRINT *,51,NF90_STRERROR(ierr)
 
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))
         IF( ierr /= nf90_NoErr) PRINT *,52,NF90_STRERROR(ierr)
      endif
 
c     Close netCDF file 
      ierr = nf90_close(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,53,NF90_STRERROR(ierr)
 
 120  continue
 
c     ---- Write data --------------------------------------------------
 
c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)
      IF( ierr /= nf90_NoErr) PRINT *,54,NF90_STRERROR(ierr)
 
c     Get the varID
      ierr = nf90_inq_varid(ncID,varname, varID )
      IF( ierr /= nf90_NoErr) PRINT *,55,NF90_STRERROR(ierr)
 
c     Get the time index
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      IF( ierr /= nf90_NoErr) PRINT *,56,NF90_STRERROR(ierr)
 
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      IF( ierr /= nf90_NoErr) PRINT *,57,NF90_STRERROR(ierr)
 
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      IF( ierr /= nf90_NoErr) PRINT *,58,NF90_STRERROR(ierr)
 
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
      IF( ierr /= nf90_NoErr) PRINT *,59,NF90_STRERROR(ierr)
      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 /) )
      IF( ierr /= nf90_NoErr) PRINT *,60,NF90_STRERROR(ierr)
 
c     Check whether writing was successful 
      ierr = nf90_close(ncID)
      IF( ierr /= nf90_NoErr) PRINT *,61,NF90_STRERROR(ierr)
 
      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 15	Rev 16
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	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon' ,clon )`	135	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon' ,clon )`
136	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat' ,clat )`	136	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat' ,clat )`
137	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )`	137	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )`
138	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )`	138	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )`
139	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx' ,nx )`	139	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx' ,nx )`
140	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny' ,ny )`	140	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny' ,ny )`
141	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx' ,dx )`	141	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx' ,dx )`
142	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy' ,dy )`	142	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy' ,dy )`
143	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin' ,xmin )`	143	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin' ,xmin )`
144	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin' ,ymin )`	144	`ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin' ,ymin )`
145		145
146	`ierr = nf90_close(ncID)`	146	`ierr = nf90_close(ncID)`
147		147
148	`print,'*** GRID PARAMETERS IMPORTED ',`	148	`print,'*** GRID PARAMETERS IMPORTED ',`
149	`> 'FROM NETCDF FILE!!!! ****'`	149	`> 'FROM NETCDF FILE!!!! ****'`
150	`print*`	150	`print*`
151		151
152	`endif`	152	`endif`
153		153
154	`c Check for consistency`	154	`c Check for consistency`
155	`if ( (step.ne.0).and.(mode.ne.'keep') ) then`	155	`if ( (step.ne.0).and.(mode.ne.'keep') ) then`
156	`print*," ERROR: interpolation is only possible for all",`	156	`print*," ERROR: interpolation is only possible for all",`
157	`> ' time steps... Stop'`	157	`> ' time steps... Stop'`
158	`stop`	158	`stop`
159	`endif`	159	`endif`
160		160
161	`c Set the number of times (just code aesthetics)`	161	`c Set the number of times (just code aesthetics)`
162	`opts=ntime`	162	`opts=ntime`
163		163
164	`c Set grid parameters for rotated grid`	164	`c Set grid parameters for rotated grid`
165	`if ( gridtype.eq.'rotated' ) then`	165	`if ( gridtype.eq.'rotated' ) then`
166	`nx = nlonlat`	166	`nx = nlonlat`
167	`ny = nlonlat`	167	`ny = nlonlat`
168	`dx = dlonlat`	168	`dx = dlonlat`
169	`dy = dlonlat`	169	`dy = dlonlat`
170	`xmin = - real(nlonlat-1)/2. * dx`	170	`xmin = - real(nlonlat-1)/2. * dx`
171	`xmax = + real(nlonlat-1)/2. * dx`	171	`xmax = + real(nlonlat-1)/2. * dx`
172	`ymin = - real(nlonlat-1)/2. * dy`	172	`ymin = - real(nlonlat-1)/2. * dy`
173	`ymax = + real(nlonlat-1)/2. * dy`	173	`ymax = + real(nlonlat-1)/2. * dy`
174	`endif`	174	`endif`
175		175
176	`c Set the flag for dynamic grid adjustment`	176	`c Set the flag for dynamic grid adjustment`
177	`if ( (nx.eq.0).or.(ny.eq.0) ) then`	177	`if ( (nx.eq.0).or.(ny.eq.0) ) then`
178	`dynamic_grid = 1`	178	`dynamic_grid = 1`
179	`else`	179	`else`
180	`dynamic_grid = 0`	180	`dynamic_grid = 0`
181	`endif`	181	`endif`
182		182
183	`c Print status information`	183	`c Print status information`
184	`print*,'---- INPUT PARAMETERS -----------------------------------'`	184	`print*,'---- INPUT PARAMETERS -----------------------------------'`
185	`print*`	185	`print*`
186	`print*,'Input : ',trim(inpfile)`	186	`print*,'Input : ',trim(inpfile)`
187	`print*,'Output : ',trim(outfile)`	187	`print*,'Output : ',trim(outfile)`
188	`print*,'Field : ',trim(field)`	188	`print*,'Field : ',trim(field)`
189	`print*,'Trajectory : ',ntime,nfield,ntra`	189	`print*,'Trajectory : ',ntime,nfield,ntra`
190	`print*,'Grid type : ',trim(gridtype)`	190	`print*,'Grid type : ',trim(gridtype)`
191	`if ( dynamic_grid.eq.1 ) then`	191	`if ( dynamic_grid.eq.1 ) then`
192	`print*,'Grid : dynamic (see below)'`	192	`print*,'Grid : dynamic (see below)'`
193	`elseif ( gridtype.eq.'latlon' ) then`	193	`elseif ( gridtype.eq.'latlon' ) then`
194	`print*,'Grid nlon,nlat : ',nx,ny`	194	`print*,'Grid nlon,nlat : ',nx,ny`
195	`print*,' lonmin,latmin : ',xmin,ymin`	195	`print*,' lonmin,latmin : ',xmin,ymin`
196	`print*,' dlon,dlat : ',dx,dy`	196	`print*,' dlon,dlat : ',dx,dy`
197	`elseif ( gridtype.eq.'rotated' ) then`	197	`elseif ( gridtype.eq.'rotated' ) then`
198	`print*,'Grid clon,clat : ',clon,clat`	198	`print*,'Grid clon,clat : ',clon,clat`
199	`print*,' nlonlat : ',nlonlat`	199	`print*,' nlonlat : ',nlonlat`
200	`print*,' dlonlat : ',dlonlat`	200	`print*,' dlonlat : ',dlonlat`
201	`endif`	201	`endif`
202	`print*,'Filter radius : ',radius,' ',trim(runit)`	202	`print*,'Filter radius : ',radius,' ',trim(runit)`
203	`print*,'Mode : ',trim(mode)`	203	`print*,'Mode : ',trim(mode)`
204	`if ( ( mode.eq.'time' ).or.`	204	`if ( ( mode.eq.'time' ).or.`
205	`> ( mode.eq.'space' ).or.`	205	`> ( mode.eq.'space' ).or.`
206	`> (mode.eq.'grid' ) )`	206	`> (mode.eq.'grid' ) )`
207	`>then`	207	`>then`
208	`print*,'Parameter : ',param`	208	`print*,'Parameter : ',param`
209	`endif`	209	`endif`
210	`if ( step.eq.0 ) then`	210	`if ( step.eq.0 ) then`
211	`print*,'Time step : all'`	211	`print*,'Time step : all'`
212	`elseif (step.gt.0) then`	212	`elseif (step.gt.0) then`
213	`print*,'Time step : ',step`	213	`print*,'Time step : ',step`
214	`endif`	214	`endif`
215	`print*,'Selection file : ',trim(sel_file)`	215	`print*,'Selection file : ',trim(sel_file)`
216	`print*,'Selection format : ',trim(sel_file)`	216	`print*,'Selection format : ',trim(sel_file)`
217	`print*,'Flag <crefile> : ',crefile`	217	`print*,'Flag <crefile> : ',crefile`
218	`print*,'Flag <crevar> : ',crevar`	218	`print*,'Flag <crevar> : ',crevar`
219		219
220	`c Check whether mode is valid`	220	`c Check whether mode is valid`
221	`if ((mode.ne.'keep' ).and.`	221	`if ((mode.ne.'keep' ).and.`
222	`> (mode.ne.'time' ).and.`	222	`> (mode.ne.'time' ).and.`
223	`> (mode.ne.'space' ).and.`	223	`> (mode.ne.'space' ).and.`
224	`> (mode.ne.'grid' ))`	224	`> (mode.ne.'grid' ))`
225	`>then`	225	`>then`
226	`print*,' ERROR: Invalid mode ',trim(mode)`	226	`print*,' ERROR: Invalid mode ',trim(mode)`
227	`stop`	227	`stop`
228	`endif`	228	`endif`
229		229
230	`c Allocate memory for old and new (reparameterised) trajectory`	230	`c Allocate memory for old and new (reparameterised) trajectory`
231	`allocate(olon(ntime),stat=stat)`	231	`allocate(olon(ntime),stat=stat)`
232	`if (stat.ne.0) print,' error allocating array olon *'`	232	`if (stat.ne.0) print,' error allocating array olon *'`
233	`allocate(olat(ntime),stat=stat)`	233	`allocate(olat(ntime),stat=stat)`
234	`if (stat.ne.0) print,' error allocating array olat *'`	234	`if (stat.ne.0) print,' error allocating array olat *'`
235	`allocate(otim(ntime),stat=stat)`	235	`allocate(otim(ntime),stat=stat)`
236	`if (stat.ne.0) print,' error allocating array otim *'`	236	`if (stat.ne.0) print,' error allocating array otim *'`
237	`allocate(nlon(1000*ntime),stat=stat)`	237	`allocate(nlon(1000*ntime),stat=stat)`
238	`if (stat.ne.0) print,' error allocating array nlon *'`	238	`if (stat.ne.0) print,' error allocating array nlon *'`
239	`allocate(nlat(1000*ntime),stat=stat)`	239	`allocate(nlat(1000*ntime),stat=stat)`
240	`if (stat.ne.0) print,' error allocating array nlat *'`	240	`if (stat.ne.0) print,' error allocating array nlat *'`
241	`allocate(ntim(1000*ntime),stat=stat)`	241	`allocate(ntim(1000*ntime),stat=stat)`
242	`if (stat.ne.0) print,' error allocating array ntim *'`	242	`if (stat.ne.0) print,' error allocating array ntim *'`
243	`allocate(odist(ntime),stat=stat)`	243	`allocate(odist(ntime),stat=stat)`
244	`if (stat.ne.0) print,' error allocating array odist *'`	244	`if (stat.ne.0) print,' error allocating array odist *'`
245	`allocate(ndist(1000*ntime),stat=stat)`	245	`allocate(ndist(1000*ntime),stat=stat)`
246	`if (stat.ne.0) print,' error allocating array ndist *'`	246	`if (stat.ne.0) print,' error allocating array ndist *'`
247	`allocate(ofld(ntime),stat=stat)`	247	`allocate(ofld(ntime),stat=stat)`
248	`if (stat.ne.0) print,' error allocating array ofld *'`	248	`if (stat.ne.0) print,' error allocating array ofld *'`
249	`allocate(nfld(1000*ntime),stat=stat)`	249	`allocate(nfld(1000*ntime),stat=stat)`
250	`if (stat.ne.0) print,' error allocating array nfld *'`	250	`if (stat.ne.0) print,' error allocating array nfld *'`
251		251
252	`c Allocate memory for complete trajectory set`	252	`c Allocate memory for complete trajectory set`
253	`allocate(trainp(ntra,ntime,nfield),stat=stat)`	253	`allocate(trainp(ntra,ntime,nfield),stat=stat)`
254	`if (stat.ne.0) print,' error allocating array trainp *'`	254	`if (stat.ne.0) print,' error allocating array trainp *'`
255	`allocate(sel_flag(ntra),stat=stat)`	255	`allocate(sel_flag(ntra),stat=stat)`
256	`if (stat.ne.0) print,' error allocating array sel_flag *'`	256	`if (stat.ne.0) print,' error allocating array sel_flag *'`
257		257
258	`c Allocate memory for auxiliary fields`	258	`c Allocate memory for auxiliary fields`
259	`allocate(traj(nfield),stat=stat)`	259	`allocate(traj(nfield),stat=stat)`
260	`if (stat.ne.0) print,' error allocating array traj *'`	260	`if (stat.ne.0) print,' error allocating array traj *'`
261		261
262	`c Set the format of the input file`	262	`c Set the format of the input file`
263	`call mode_tra(inpmode,inpfile)`	263	`call mode_tra(inpmode,inpfile)`
264	`if (inpmode.eq.-1) inpmode=1`	264	`if (inpmode.eq.-1) inpmode=1`
265		265
266	`c Read the input trajectory file`	266	`c Read the input trajectory file`
267	`call ropen_tra(fid,inpfile,ntra,ntime,nfield,`	267	`call ropen_tra(fid,inpfile,ntra,ntime,nfield,`
268	`> reftime,varsinp,inpmode)`	268	`> reftime,varsinp,inpmode)`
269	`call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)`	269	`call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)`
270	`call close_tra(fid,inpmode)`	270	`call close_tra(fid,inpmode)`
271		271
272	`c Check that first four columns correspond to time,lon,lat,p`	272	`c Check that first four columns correspond to time,lon,lat,p`
273	`if ( (varsinp(1).ne.'time' ).or.`	273	`if ( (varsinp(1).ne.'time' ).or.`
274	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`	274	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`
275	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`	275	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`
276	`> (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p' ) )`	276	`> (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p' ) )`
277	`>then`	277	`>then`
278	`print*,' ERROR: problem with input trajectories ...'`	278	`print*,' ERROR: problem with input trajectories ...'`
279	`stop`	279	`stop`
280	`endif`	280	`endif`
281	`varsinp(1) = 'TIME'`	281	`varsinp(1) = 'TIME'`
282	`varsinp(2) = 'lon'`	282	`varsinp(2) = 'lon'`
283	`varsinp(3) = 'lat'`	283	`varsinp(3) = 'lat'`
284	`varsinp(4) = 'p'`	284	`varsinp(4) = 'p'`
285		285
286	`c Get the index of the field (if needed)`	286	`c Get the index of the field (if needed)`
287	`if ( field.ne.'nil' ) then`	287	`if ( field.ne.'nil' ) then`
288	`ifield = 0`	288	`ifield = 0`
289	`do i=1,nfield`	289	`do i=1,nfield`
290	`if ( varsinp(i).eq.field ) ifield = i`	290	`if ( varsinp(i).eq.field ) ifield = i`
291	`enddo`	291	`enddo`
292	`if ( ifield.eq.0 ) then`	292	`if ( ifield.eq.0 ) then`
293	`print*,' ERROR: field ',trim(field),' not found... Stop'`	293	`print*,' ERROR: field ',trim(field),' not found... Stop'`
294	`stop`	294	`stop`
295	`endif`	295	`endif`
296	`endif`	296	`endif`
297		297
298	`c Write some status information of the input trajectories`	298	`c Write some status information of the input trajectories`
299	`print*`	299	`print*`
300	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`	300	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`
301	`print*`	301	`print*`
302	`print*,' Reference time (year) : ',reftime(1)`	302	`print*,' Reference time (year) : ',reftime(1)`
303	`print*,' (month) : ',reftime(2)`	303	`print*,' (month) : ',reftime(2)`
304	`print*,' (day) : ',reftime(3)`	304	`print*,' (day) : ',reftime(3)`
305	`print*,' (hour) : ',reftime(4)`	305	`print*,' (hour) : ',reftime(4)`
306	`print*,' (min) : ',reftime(5)`	306	`print*,' (min) : ',reftime(5)`
307	`print*,' Time range (min) : ',reftime(6)`	307	`print*,' Time range (min) : ',reftime(6)`
308	`do i=1,nfield`	308	`do i=1,nfield`
309	`if ( i.ne.ifield ) then`	309	`if ( i.ne.ifield ) then`
310	`print*,' Var :',i,trim(varsinp(i))`	310	`print*,' Var :',i,trim(varsinp(i))`
311	`else`	311	`else`
312	`print*,' Var :',i,trim(varsinp(i)),`	312	`print*,' Var :',i,trim(varsinp(i)),`
313	`> ' [ gridding ]'`	313	`> ' [ gridding ]'`
314	`endif`	314	`endif`
315	`enddo`	315	`enddo`
316	`print*,' List of selected times'`	316	`print*,' List of selected times'`
317	`do i=1,ntime`	317	`do i=1,ntime`
318	`if ( (step.eq.0).or.(step.eq.i) ) then`	318	`if ( (step.eq.0).or.(step.eq.i) ) then`
319	`print*,' ',i,' -> ',trainp(1,i,1)`	319	`print*,' ',i,' -> ',trainp(1,i,1)`
320	`endif`	320	`endif`
321	`enddo`	321	`enddo`
322	`print*`	322	`print*`
323		323
324	`c Select flag: all trajectories are selected`	324	`c Select flag: all trajectories are selected`
325	`if ( sel_file.eq.'nil' ) then`	325	`if ( sel_file.eq.'nil' ) then`
326		326
327	`do i=1,ntra`	327	`do i=1,ntra`
328	`sel_flag(i) = 1`	328	`sel_flag(i) = 1`
329	`enddo`	329	`enddo`
330		330
331	`c Select flag: index file`	331	`c Select flag: index file`
332	`elseif ( sel_format.eq.'index' ) then`	332	`elseif ( sel_format.eq.'index' ) then`
333		333
334	`do i=1,ntra`	334	`do i=1,ntra`
335	`sel_flag(i) = 0`	335	`sel_flag(i) = 0`
336	`enddo`	336	`enddo`
337		337
338	`open(10,file=sel_file)`	338	`open(10,file=sel_file)`
339	`142 read(10,*,end=141) ind`	339	`142 read(10,*,end=141) ind`
340	`sel_flag(ind) = 1`	340	`sel_flag(ind) = 1`
341	`goto 142`	341	`goto 142`
342	`141 continue`	342	`141 continue`
343	`close(10)`	343	`close(10)`
344		344
345	`c Select flag: boolean file`	345	`c Select flag: boolean file`
346	`elseif ( sel_format.eq.'boolean' ) then`	346	`elseif ( sel_format.eq.'boolean' ) then`
347		347
348	`open(10,file=sel_file)`	348	`open(10,file=sel_file)`
349	`do i=1,ntra`	349	`do i=1,ntra`
350	`read(10,*) ind`	350	`read(10,*) ind`
351	`if ( ind.eq.1 ) sel_flag(i) = ind`	351	`if ( ind.eq.1 ) sel_flag(i) = ind`
352	`enddo`	352	`enddo`
353	`close(10)`	353	`close(10)`
354		354
355	`endif`	355	`endif`
356		356
357	`c Write status information`	357	`c Write status information`
358	`if ( sel_file.eq.'nil' ) then`	358	`if ( sel_file.eq.'nil' ) then`
359	`print*,' Selected trajectories : all ',ntra`	359	`print*,' Selected trajectories : all ',ntra`
360	`else`	360	`else`
361	`count = 0`	361	`count = 0`
362	`do i=1,ntra`	362	`do i=1,ntra`
363	`if ( sel_flag(i).eq.1 ) count = count + 1`	363	`if ( sel_flag(i).eq.1 ) count = count + 1`
364	`enddo`	364	`enddo`
365	`print*,' #selected trajectories : ',count,`	365	`print*,' #selected trajectories : ',count,`
366	`> ' [ ',real(count)/real(ntra) * 100.,' % ] '`	366	`> ' [ ',real(count)/real(ntra) * 100.,' % ] '`
367	`endif`	367	`endif`
368	`print*`	368	`print*`
369		369
370	`c ---------------------------------------------------------------------`	370	`c ---------------------------------------------------------------------`
371	`c Coordinate transformations and grid adjustment`	371	`c Coordinate transformations and grid adjustment`
372	`c ---------------------------------------------------------------------`	372	`c ---------------------------------------------------------------------`
373		373
374	`c Transform from lat/lon to rotated lat/lon, if requested`	374	`c Transform from lat/lon to rotated lat/lon, if requested`
375	`if ( gridtype.eq.'rotated') then`	375	`if ( gridtype.eq.'rotated') then`
376		376
377	`crot = 0.`	377	`crot = 0.`
378		378
379	`pollon=clon-180.`	379	`pollon=clon-180.`
380	`if (pollon.lt.-180.) pollon=pollon+360.`	380	`if (pollon.lt.-180.) pollon=pollon+360.`
381	`pollat=90.-clat`	381	`pollat=90.-clat`
382	`do i=1,ntra`	382	`do i=1,ntra`
383	`do j=1,ntime`	383	`do j=1,ntime`
384		384
385	`if ( sel_flag(i).eq.1 ) then`	385	`if ( sel_flag(i).eq.1 ) then`
386		386
387	`c Get lat/lon coordinates for trajectory point`	387	`c Get lat/lon coordinates for trajectory point`
388	`lon = trainp(i,j,2)`	388	`lon = trainp(i,j,2)`
389	`lat = trainp(i,j,3)`	389	`lat = trainp(i,j,3)`
390		390
391	`c First Rotation`	391	`c First Rotation`
392	`pollon=clon-180.`	392	`pollon=clon-180.`
393	`if (pollon.lt.-180.) pollon=pollon+360.`	393	`if (pollon.lt.-180.) pollon=pollon+360.`
394	`pollat=90.-clat`	394	`pollat=90.-clat`
395	`rlon0=lmtolms(lat,lon,pollat,pollon)`	395	`rlon0=lmtolms(lat,lon,pollat,pollon)`
396	`rlat0=phtophs(lat,lon,pollat,pollon)`	396	`rlat0=phtophs(lat,lon,pollat,pollon)`
397		397
398	`c Second rotation`	398	`c Second rotation`
399	`pollon=-180.`	399	`pollon=-180.`
400	`pollat=90.+crot`	400	`pollat=90.+crot`
401	`rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)`	401	`rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)`
402	`rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)`	402	`rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)`
403		403
404	`c Get rotated latitude and longitude`	404	`c Get rotated latitude and longitude`
405	`100 if (rlon.lt.xmin) then`	405	`100 if (rlon.lt.xmin) then`
406	`rlon=rlon+360.`	406	`rlon=rlon+360.`
407	`goto 100`	407	`goto 100`
408	`endif`	408	`endif`
409	`102 if (rlon.gt.(xmin+real(nx-1)*dx)) then`	409	`102 if (rlon.gt.(xmin+real(nx-1)*dx)) then`
410	`rlon=rlon-360.`	410	`rlon=rlon-360.`
411	`goto 102`	411	`goto 102`
412	`endif`	412	`endif`
413		413
414	`c Set the new trajectory coordinates`	414	`c Set the new trajectory coordinates`
415	`trainp(i,j,2) = rlon`	415	`trainp(i,j,2) = rlon`
416	`trainp(i,j,3) = rlat`	416	`trainp(i,j,3) = rlat`
417		417
418	`endif`	418	`endif`
419		419
420	`enddo`	420	`enddo`
421	`enddo`	421	`enddo`
422		422
423	`endif`	423	`endif`
424		424
425	`c Dynamic grid adjustment`	425	`c Dynamic grid adjustment`
426	`if ( dynamic_grid.eq.1 ) then`	426	`if ( dynamic_grid.eq.1 ) then`
427		427
428	`c Get the grid parameters`	428	`c Get the grid parameters`
429	`xmin = 180.`	429	`xmin = 180.`
430	`ymin = 90.`	430	`ymin = 90.`
431	`xmax = -180.`	431	`xmax = -180.`
432	`ymax = -90.`	432	`ymax = -90.`
433		433
434	`do i=1,ntra`	434	`do i=1,ntra`
435		435
436	`if ( sel_flag(i).eq.1 ) then`	436	`if ( sel_flag(i).eq.1 ) then`
437		437
438	`if ( step.eq.0 ) then`	438	`if ( step.eq.0 ) then`
439	`do j=1,ntime`	439	`do j=1,ntime`
440	`if ( trainp(i,j,2).lt.xmin) xmin = trainp(i,j,2)`	440	`if ( trainp(i,j,2).lt.xmin) xmin = trainp(i,j,2)`
441	`if ( trainp(i,j,2).gt.xmax) xmax = trainp(i,j,2)`	441	`if ( trainp(i,j,2).gt.xmax) xmax = trainp(i,j,2)`
442	`if ( trainp(i,j,3).lt.ymin) ymin = trainp(i,j,3)`	442	`if ( trainp(i,j,3).lt.ymin) ymin = trainp(i,j,3)`
443	`if ( trainp(i,j,3).gt.ymax) ymax = trainp(i,j,3)`	443	`if ( trainp(i,j,3).gt.ymax) ymax = trainp(i,j,3)`
444	`enddo`	444	`enddo`
445	`else`	445	`else`
446	`if ( trainp(i,step,2).lt.xmin) xmin = trainp(i,step,2)`	446	`if ( trainp(i,step,2).lt.xmin) xmin = trainp(i,step,2)`
447	`if ( trainp(i,step,2).gt.xmax) xmax = trainp(i,step,2)`	447	`if ( trainp(i,step,2).gt.xmax) xmax = trainp(i,step,2)`
448	`if ( trainp(i,step,3).lt.ymin) ymin = trainp(i,step,3)`	448	`if ( trainp(i,step,3).lt.ymin) ymin = trainp(i,step,3)`
449	`if ( trainp(i,step,3).gt.ymax) ymax = trainp(i,step,3)`	449	`if ( trainp(i,step,3).gt.ymax) ymax = trainp(i,step,3)`
450	`endif`	450	`endif`
451		451
452	`endif`	452	`endif`
453		453
454	`enddo`	454	`enddo`
455		455
456	`c Get first guess for "optimal" grid`	456	`c Get first guess for "optimal" grid`
457	`nx = 400`	457	`nx = 400`
458	`ny = 400`	458	`ny = 400`
459	`dx = (xmax - xmin)/real(nx-1)`	459	`dx = (xmax - xmin)/real(nx-1)`
460	`dy = (ymax - ymin)/real(ny-1)`	460	`dy = (ymax - ymin)/real(ny-1)`
461		461
462	`c Make the grid spacing equal in zonal and meridional direction`	462	`c Make the grid spacing equal in zonal and meridional direction`
463	`if ( dx.gt.dy ) then`	463	`if ( dx.gt.dy ) then`
464		464
465	`dy = dx`	465	`dy = dx`
466	`ny = (ymax - ymin)/dy + 1`	466	`ny = (ymax - ymin)/dy + 1`
467	`if (ny.lt.nx/2) ny = nx / 2`	467	`if (ny.lt.nx/2) ny = nx / 2`
468	`if ( real(ny)*dy .ge. 180. ) ny = 180./dy + 1`	468	`if ( real(ny)*dy .ge. 180. ) ny = 180./dy + 1`
469	`ycen = 0.5* (ymin+ymax)`	469	`ycen = 0.5* (ymin+ymax)`
470	`ymin = ycen - 0.5 * real(ny/2) * dy`	470	`ymin = ycen - 0.5 * real(ny/2) * dy`
471	`if (ymin.le.-90.) ymin = -90.`	471	`if (ymin.le.-90.) ymin = -90.`
472		472
473	`else`	473	`else`
474		474
475	`dx = dy`	475	`dx = dy`
476	`nx = (xmax - xmin)/dx + 1`	476	`nx = (xmax - xmin)/dx + 1`
477	`if (nx.lt.ny/2) nx = ny / 2`	477	`if (nx.lt.ny/2) nx = ny / 2`
478	`if ( real(nx)*dx .ge. 360. ) nx = 360./dx + 1`	478	`if ( real(nx)*dx .ge. 360. ) nx = 360./dx + 1`
479	`xcen = 0.5* (xmin+xmax)`	479	`xcen = 0.5* (xmin+xmax)`
480	`xmin = xcen - 0.5 * real(nx/2) * dx`	480	`xmin = xcen - 0.5 * real(nx/2) * dx`
481	`if (xmin.le.-180.) xmin = -180.`	481	`if (xmin.le.-180.) xmin = -180.`
482		482
483	`endif`	483	`endif`
484		484
485	`c Write information`	485	`c Write information`
486	`print*`	486	`print*`
487	`print*,'---- DYNAMIC GRID ADJUSTMENT',`	487	`print*,'---- DYNAMIC GRID ADJUSTMENT',`
488	`> ' ----------------------------'`	488	`> ' ----------------------------'`
489	`print*`	489	`print*`
490	`print*,'Grid nlon,nlat : ',nx,ny`	490	`print*,'Grid nlon,nlat : ',nx,ny`
491	`print*,' lonmin,latmin : ',xmin,ymin`	491	`print*,' lonmin,latmin : ',xmin,ymin`
492	`print*,' dlon,dlat : ',dx,dy`	492	`print*,' dlon,dlat : ',dx,dy`
493	`print*`	493	`print*`
494		494
495	`c Write grid information for rotated grid (if not already done`	495	`c Write grid information for rotated grid (if not already done`
496	`elseif ( gridtype.eq.'rotated') then`	496	`elseif ( gridtype.eq.'rotated') then`
497		497
498	`print*`	498	`print*`
499	`print*,'---- GRID PARAMETERS -------',`	499	`print*,'---- GRID PARAMETERS -------',`
500	`> ' ----------------------------'`	500	`> ' ----------------------------'`

Subversion Repositories lagranto.um

(root)/trunk/density/density.f – Rev 15 → 16