Subversion Repositories lagranto.wrf

      PROGRAM wrfmap

c     **********************************************************************
c     * Transform between lon/lat and x/y coordinates                      *
c     * Michael Sprenger / Spring 2013                                     *
c     **********************************************************************

      implicit none
     
c     ------------------------------------------------------------
c     Declaration of parameters and variables
c     ------------------------------------------------------------

c     Input parameters
      character*80 mode
      character*80 inpfile
      character*80 outfile
      integer      ntra,ntim,ncol
      character*80 anglemode

c     Fixed parameters
      character*80 mapfile                          ! netCDF file with map transformation
      parameter    (mapfile ='wrfmap.nc')
      real         mdv                              ! missing data value
      parameter    (mdv=-999.)
      real         eps                              ! Numerical epsilon
      parameter    (eps= 0.001)

c     Numerical grid
      real         xmin,xmax,ymin,ymax              ! Domain size
      real         dx,dy                            ! Horizontal resolution
      integer      nx,ny,nz                         ! Grid dimensions
      real,allocatable,dimension (:,:)  :: lon,lat  ! lon/lat on numerical grid
      real,allocatable,dimension (:,:)  :: mpx,mpy  ! map scale factors in x/y direction


c     Geographical grid
      real         latmin,latmax,lonmin,lonmax      ! Geographical domain
      real         dlon,dlat                        ! Minimum spacing in geographical space
      integer      nlon,nlat,nlev                        ! Geographical grid dimension
      real,allocatable,dimension (:,:)   :: x,y      ! x/y on geographical grid

c     Vertical grid
      real,allocatable,dimension (:,:,:)   :: p3    ! 3D pressure
      real,allocatable,dimension (:,:)     :: ps    ! surface pressure
      real,allocatable,dimension (:,:,:)   :: z3    ! 3D geopotential height
      real,allocatable,dimension (:,:)     :: zb    ! surface geopotential height

c     Transformation between numerical and geographical grid
      integer,allocatable,dimension (:,:)  :: connect1 
      integer      connectval1 
      integer,allocatable,dimension (:,:)  :: connect2 
      integer      connectval2
      real,allocatable,dimension (:,:)  :: xc,yc
      real         radius
      real         xval,yval

c     Trajectories
      real,allocatable, dimension (:,:,:) :: tra             ! Input start coordinates (ntra,ntim,ncol)
      integer                                refdate(6)      ! Reference date
      character*80                           vars(200)       ! Field names
      real,allocatable, dimension (:)     :: xx0,yy0,zz0     ! Position of air parcels
      integer                                inpmode
      integer                                outmode
      integer                                npoints

c     Auxiliary variables
      integer      fid
      integer      stat
      character*80 varname,cdfname
      integer      i,j
      character*80 radunit
      real         rdummy
      real         rid,rjd,rkd
      real         xpos,ypos,zpos,ppos,lonpos,latpos
      real         lon1,lat1,lon2,lat2
      integer      nfilter

c     Externals
      real         int_index3
      external     int_index3
      real         sdis
      external     sdis

c     ------------------------------------------------------------
c     Preparations
c     ------------------------------------------------------------

c     Read parameters
      open(10,file='wrfmap.param')

       read(10,*) mode

       if ( mode.eq.'-create' ) then         ! create mapping file
          read(10,*) inpfile
          read(10,*) anglemode
       endif

       if ( mode.eq.'-ll2xy'  ) then         ! convert lon/lat -> x/y
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif

       if ( mode.eq.'-ll2xy.single'  ) then   ! convert single lon/lat -> x/y
          read(10,*) lonpos,latpos
       endif

       if ( mode.eq.'-xy2ll'  ) then         ! convert x/y -> lon/lat
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif

       if ( mode.eq.'-xy2ll.single'  ) then   ! convert single x/y -> lon/lat
          read(10,*) xpos,ypos
       endif

       if ( mode.eq.'-p2z'  ) then            ! convert x/y/p -> x/y/z
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif

       if ( mode.eq.'-p2z.single'  ) then     ! convert single x/y/p -> x/y/z
          read(10,*) xpos,ypos,ppos
       endif

       if ( mode.eq.'-z2p'  ) then            ! convert x/y/z -> x/y/p
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif

       if ( mode.eq.'-z2p.single'  ) then     ! convert single x/y/z -> x/y/p
          read(10,*) xpos,ypos,zpos
       endif

       if ( mode.eq.'-mapscale'  ) then       ! calculate mapping scale factors
          read(10,*) nfilter
       endif

      close(10)

c     Read the mapping file - except for mode '-create'
      if ( mode.ne.'-create' ) then

c        Get dimensions 
         cdfname = mapfile
         varname = 'Z'
         nx      = 1
         ny      = 1
         nz      = 1
         call readcdf2D(cdfname,varname,rdummy,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)

         cdfname = mapfile
         varname = 'X'
         nlon    = 1
         nlat    = 1
         nlev    = 1
         call readcdf2D(cdfname,varname,rdummy,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,nlev)

c        Allocate memory
         allocate(lon(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array lon ***'
         allocate(lat(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array lat ***'
         allocate(x(nlon,nlat),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array x   ***'
         allocate(y(nlon,nlat),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array y   ***'
         allocate(p3(nx,ny,nz),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array p3  ***'
         allocate(z3(nx,ny,nz),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array z3  ***'
         allocate(zb(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array zb  ***'
         allocate(ps(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array ps  ***'

c        Read data
         cdfname = mapfile
         varname = 'LON'
         call readcdf2D(cdfname,varname,lon,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)
         cdfname = mapfile
         varname = 'LAT'
         call readcdf2D(cdfname,varname,lat,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)
         cdfname = mapfile
         varname = 'X'
         call readcdf2D(cdfname,varname,x,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,1)
         cdfname = mapfile
         varname = 'Y'
         call readcdf2D(cdfname,varname,y,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,1)
         cdfname = mapfile
         varname = 'Z'
         call readcdf2D(cdfname,varname,z3,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)
         cdfname = mapfile
         varname = 'P'
         call readcdf2D(cdfname,varname,p3,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)
         cdfname = mapfile
         varname = 'TOPO'
         call readcdf2D(cdfname,varname,zb,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)

c            Set surface prtessure to lowest 3d pressure level
         do i=1,nx
            do j=1,ny
               ps(i,j) = p3(i,j,1)
            enddo
         enddo

      endif

c     ------------------------------------------------------------
c     Create mapping file
c     ------------------------------------------------------------
 
      if ( mode.eq.'-create' ) then

c     Open mapping file
      call input_open(fid,inpfile)

c     Get grid description
      call input_grid_wrf(fid,
     >                    xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)

c     Write grid information to screen
      print*,' xmin,xmax  = ',xmin,xmax
      print*,' ymin,ymax  = ',ymin,ymax
      print*,' dx,dy      = ',dx,dy
      print*,' nx,ny,nz   = ',nx,ny,nz 
      print*,' anglemode  = ',trim(anglemode)

c     Transform grid specifier into grid coordinates
      xmin = 1.
      ymin = 1.
      dx   = 1.
      dy   = 1.

c     Allocate memory for lon/lat on numerical grid
      allocate(lon(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon ***'
      allocate(lat(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat ***'
      allocate(p3(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3  ***'
      allocate(z3(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array z3  ***'
      allocate(zb(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zb  ***'

c     Read lon/lat on the model grid
      varname='XLONG'
      call input_var_wrf(fid,varname,lon,nx,ny,1)
      varname='XLAT'
      call input_var_wrf(fid,varname,lat,nx,ny,1)
      varname='pressure'
      call input_var_wrf(fid,varname,p3 ,nx,ny,nz)
      varname='geopot'
      call input_var_wrf(fid,varname,z3 ,nx,ny,nz)
      varname='geopots'
      call input_var_wrf(fid,varname,zb ,nx,ny,1)

c     Transform longitude depending on <anglemode>
      if ( anglemode.eq.'greenwich' ) then
         do i=1,nx
           do j=1,ny
              if ( lon(i,j).lt.0. ) lon(i,j) = lon(i,j) + 360.
            enddo
         enddo
      elseif ( anglemode.eq.'dateline' ) then
         do i=1,nx
            do j=1,ny
               if ( lon(i,j).gt.180. ) lon(i,j) = lon(i,j) - 360.
            enddo
         enddo
      else
         print*,' ERROR: unsupported angle mode... Stop'
         stop
      endif

c     Close input file file
      call input_close(fid)

c     Check for 'date line' and for pole
      do i=2,nx
      do j=1,ny
        if ( abs( lat(i,j) ).gt.90. ) then
           print*,' ERROR: mapping over pole not supported... Stop'
           stop
        endif
        if ( abs( lon(i,j)-lon(i-1,j) ).gt.180. ) then
           print*,' ERROR: mapping over date line not supported... Stop'
           stop
        endif
      enddo
      enddo

c     Determine the extreme coordinates
      latmin  = lat(1,1)
      latmax  = lat(1,1)
      lonmin  = lon(1,1)
      lonmax  = lon(1,1)
      do i=1,nx
            do j=1,ny
          if ( lat(i,j).lt.latmin ) latmin = lat(i,j)
          if ( lat(i,j).gt.latmax ) latmax = lat(i,j)
          if ( lon(i,j).lt.lonmin ) lonmin = lon(i,j)
          if ( lon(i,j).gt.lonmax ) lonmax = lon(i,j)
        enddo
      enddo

      print*,' lonmin,max = ',lonmin,lonmax
      print*,' latmin,max = ',latmin,latmax

c     Determine the extreme dlon/dlat spacing
      dlon = abs( lon(2,1)-lon(1,1) )
      do i=2,nx
            do j=1,ny
           if ( abs( lon(i,j)-lon(i-1,j) ).lt.dlon ) then
              dlon = abs( lon(i,j)-lon(i-1,j) )
           endif
         enddo
      enddo
      dlat = abs( lat(1,2)-lat(1,1) )
      do i=1,nx
            do j=2,ny
           if ( abs( lat(i,j)-lat(i,j-1) ).lt.dlat ) then
              dlat = abs( lat(i,j)-lat(i,j-1) )
           endif
         enddo
      enddo

c     Set dimension of geographical grid 
      nlon   = ceiling( (lonmax-lonmin) / dlon + 1. )
      nlat   = ceiling( (latmax-latmin) / dlat + 1. )
      lonmax = lonmin + real(nlon-1) * dlon
      latmax = latmin + real(nlat-1) * dlat

      print*,' dlon,dlat  = ',dlon,dlat
      print*,' nlon,nlat  = ',nlon,nlat

c     Allocate memory for x/y on geographical grid
      allocate(x(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array x       ***'
      allocate(y(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array y       ***'
      allocate(connect1(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect1***'
      allocate(connect2(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect2***'
      allocate(xc(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xc      ***'
      allocate(yc(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yc      ***'

c     Init the mapping arrays
      connectval1 = 0
      connectval2 = 0
      do i=1,nlon
        do j=1,nlat
           x(i,j)        = 0.
           y(i,j)        = 0.
           xc(i,j)       = 0.
           yc(i,j)       = 0.
           connect1(i,j) = 1
           connect2(i,j) = 1
        enddo
      enddo

c     Get the good radius for gridding - avoiding gaps
      radunit = 'deg'
      radius  = 0.
      do i=2,nx
        do j=2,ny
           if ( abs( lat(i,j)-lat(i,j-1) ).gt.radius ) then
              radius = abs( lat(i,j)-lat(i,j-1) )
           endif
           if ( abs( lon(i,j)-lon(i-1,j) ).gt.radius ) then
              radius = abs( lon(i,j)-lon(i-1,j) )
           endif
        enddo
      enddo
      radius = 3. * radius 

c     Do the mapping
      do i=1,nx
        do j=1,ny

          connectval1 = connectval1 + 1
          call gridding1(lat(i,j),lon(i,j),real(i),radius,radunit,
     >                   connect1,connectval1,
     >                   x,xc,nlon,nlat,lonmin,latmin,dlon,dlat)

          connectval2 = connectval2 + 1
          call gridding1(lat(i,j),lon(i,j),real(j),radius,radunit,
     >                   connect2,connectval2,
     >                   y,yc,nlon,nlat,lonmin,latmin,dlon,dlat)


        enddo
      enddo

c     Normalize output and set miising data flag
      do i=1,nlon
         do j=1,nlat
            if ( xc(i,j).gt.0. ) then
               x(i,j) = x(i,j)/xc(i,j)
            else
               x(i,j) = mdv
            endif
            if ( yc(i,j).gt.0. ) then
               y(i,j) = y(i,j)/yc(i,j)
            else
               y(i,j) = mdv
            endif
         enddo
      enddo

c     Write data to netCDF
      cdfname = mapfile
      varname = 'X'
      call writecdf2D(cdfname,varname,x,0.,
     >                dlon,dlat,lonmin,latmin,nlon,nlat,1,1,1)

      cdfname = mapfile
      varname = 'Y'
      call writecdf2D(cdfname,varname,y,0.,
     >                dlon,dlat,lonmin,latmin,nlon,nlat,1,0,1)

      cdfname = mapfile
      varname = 'LON'
      call writecdf2D(cdfname,varname,lon,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)

      cdfname = mapfile
      varname = 'LAT'
      call writecdf2D(cdfname,varname,lat,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)

      cdfname = mapfile
      varname = 'Z'
      call writecdf2D(cdfname,varname,z3,0.,
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)

      cdfname = mapfile
      varname = 'P'
      call writecdf2D(cdfname,varname,p3,0.,
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
      cdfname = mapfile
      varname = 'TOPO'
      call writecdf2D(cdfname,varname,zb,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)


      endif


c     ------------------------------------------------------------
c     Convert a lat/lon/z list to a x/y/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-ll2xy' ) then

c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'

c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1 ) outmode = inpmode

c     Read coordinates from file - Format (lon,lat,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)

c     Read coordinates from trajectory file 
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)

         if ( (vars(2).ne.'lon').and.(vars(3).ne.'lat') ) then
            print*,' ERROR: input must be in lon/lat format... Stop'       
            stop
         endif

         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2) 
               yy0(npoints) = tra(i,j,3) 
               zz0(npoints) = tra(i,j,4) 
            enddo
         enddo
      endif 

c     Transform coordinates
      do i=1,npoints

          if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >         ( abs(yy0(i)-mdv).gt.eps )  )
     >    then
             rid=(xx0(i)-lonmin)/dlon+1.
             rjd=(yy0(i)-latmin)/dlat+1.
             xx0(i) = int_index3(x,nlon,nlat,1,rid,rjd,1.,mdv)
             yy0(i) = int_index3(y,nlon,nlat,1,rid,rjd,1.,mdv)
          else
             xx0(i) = mdv
             yy0(i) = mdv
          endif

      enddo

c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)

c     Write output to trajectory
      else         
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif

      endif

c     ------------------------------------------------------------
c     Convert a x/y/z list to a lon/lat/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-xy2ll' ) then

c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'

c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode

c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)

c     Read coordinates from trajectory file 
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)

         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'       
            stop
         endif

         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2) 
               yy0(npoints) = tra(i,j,3) 
               zz0(npoints) = tra(i,j,4) 
            enddo
         enddo
      endif 

c     Transform coordinates
      do i=1,npoints

          if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >         ( abs(yy0(i)-mdv).gt.eps )  )
     >    then
              rid=(xx0(i)-xmin)/dx+1.
              rjd=(yy0(i)-ymin)/dy+1.
              xx0(i) = int_index3(lon,nx,ny,1,rid,rjd,1.,mdv)
              yy0(i) = int_index3(lat,nx,ny,1,rid,rjd,1.,mdv)
          else
              xx0(i) = mdv
              yy0(i) = mdv
          endif

      enddo

c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)

c     Write output to trajectory
      else         
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'lon'
         vars(3) = 'lat'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif

      endif

c     ------------------------------------------------------------
c     Convert a single lat/lon/z list to a x/y/z list
c     ------------------------------------------------------------

      if ( mode.eq.'-ll2xy.single' ) then

c     Transform coordinates
      rid  = (lonpos-lonmin)/dlon+1.
      rjd  = (latpos-latmin)/dlat+1.
      xpos = int_index3(x,nlon,nlat,1,rid,rjd,1.,mdv)
      ypos = int_index3(y,nlon,nlat,1,rid,rjd,1.,mdv)

c         Write result to screen
      print*,xpos,ypos

      endif

c     ------------------------------------------------------------
c     Convert a single x/y/z list to a lon/lat/z list
c     ------------------------------------------------------------

      if ( mode.eq.'-xy2ll.single' ) then

c     Transform coordinates
      rid    = (xpos-xmin)/dx+1.
      rjd    = (ypos-ymin)/dy+1.
      lonpos = int_index3(lon,nx,ny,1,rid,rjd,1.,mdv)
      latpos = int_index3(lat,nx,ny,1,rid,rjd,1.,mdv)

c         Write result to screen
      print*,lonpos,latpos

      endif

c     ------------------------------------------------------------
c     Calculate numerically the maps scale factor
c     ------------------------------------------------------------

      if ( mode.eq.'-mapscale' ) then

c         Allocate memory for map scale factors
      allocate(mpx(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array mpx       ***'
      allocate(mpy(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array mpy       ***'

c     Get map scale for interior points
      do i=2,nx-1
         do j=2,ny-1

c           Map scale in x direction
            lon1     = lon(i-1,j)
            lat1     = lat(i-1,j)
            lon2     = lon(i+1,j)
            lat2     = lat(i+1,j)
            radunit  = 'km.haversine'
            mpx(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)

c           Map scale in y direction
            lon1     = lon(i,j-1)
            lat1     = lat(i,j-1)
            lon2     = lon(i,j+1)
            lat2     = lat(i,j+1)
            radunit  = 'km.haversine'
            mpy(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)

         enddo
      enddo

c     Copy map scale for boundary line points
      do i=2,nx-1
        mpx(i, 1) = mpx(i,   2)
        mpx(i,ny) = mpx(i,ny-1)
        mpy(i, 1) = mpy(i,   2)
        mpy(i,ny) = mpy(i,ny-1)
      enddo
      do j=2,ny-1
        mpx(1, j) = mpx(2,   j)
        mpx(nx,j) = mpx(nx-1,j)
        mpy(1, j) = mpy(2,   j)
        mpy(nx,j) = mpy(nx-1,j)
      enddo

c     Copy map scale factor for boundary corner points
      mpx( 1, 1) = mpx(   2,   2)
      mpy( 1, 1) = mpy(   2,   2)
      mpx(nx, 1) = mpx(nx-1,   2)
      mpy(nx, 1) = mpy(nx-1,   2)
      mpx(nx,ny) = mpx(nx-1,ny-1)
      mpy(nx,ny) = mpy(nx-1,ny-1)
      mpx( 1,ny) = mpx(   2,ny-1)
      mpy( 1,ny) = mpy(   2,ny-1)

c     Apply a diffuive filter
      do i=1,nfilter
         call filt2d (mpx,mpx,nx,ny,1.0,-999.,0,0,0,0)
         call filt2d (mpy,mpy,nx,ny,1.0,-999.,0,0,0,0)
      enddo

c     Write map factors to mapping file
      cdfname = mapfile
      varname = 'MAPSCALE_X'
      call writecdf2D(cdfname,varname,mpx,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
      cdfname = mapfile
      varname = 'MAPSCALE_Y'
      call writecdf2D(cdfname,varname,mpy,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)


      endif

c     ------------------------------------------------------------
c     Convert from pressure to height
c     ------------------------------------------------------------

      if ( mode.eq.'-p2z' ) then

c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'

c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode

c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)

c     Read coordinates from trajectory file
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)

         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'
            stop
         endif

         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2)
               yy0(npoints) = tra(i,j,3)
               zz0(npoints) = tra(i,j,4)
            enddo
         enddo
      endif

c     Transform coordinates
      do i=1,npoints

         if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >        ( abs(yy0(i)-mdv).gt.eps ).and.
     >        ( abs(zz0(i)-mdv).gt.eps ) )
     >    then
             call get_index3 (rid,rjd,rkd,xx0(i),yy0(i),-zz0(i),1,
     >                        -p3,-ps,nx,ny,nz,xmin,ymin,dx,dy)

             zz0(i) = int_index3(z3,nx,ny,nz,rid,rjd,rkd,mdv)
          else
             xx0(i) = mdv
             yy0(i) = mdv
             zz0(i) = mdv
          endif

      enddo

c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)

c     Write output to trajectory
      else
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif

      endif

c     ------------------------------------------------------------
c     Convert single x/y/p to x/y/z
c     ------------------------------------------------------------

      if ( mode.eq.'-p2z.single' ) then

c     Transform coordinates - change from descending to ascending (minus sign)
      call get_index3 (rid,rjd,rkd,xpos,ypos,-ppos,1,
     >                 -p3,-ps,nx,ny,nz,xmin,ymin,dx,dy)

      zpos = int_index3(z3,nx,ny,nz,rid,rjd,rkd,mdv)

c         Write result to screen
      print*,xpos,ypos,zpos

      endif

c     ------------------------------------------------------------
c     Convert from height to pressure
c     ------------------------------------------------------------

      if ( mode.eq.'-z2p' ) then

c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'

c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode

c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)

c     Read coordinates from trajectory file
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)

         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'
            stop
         endif

         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2)
               yy0(npoints) = tra(i,j,3)
               zz0(npoints) = tra(i,j,4)
            enddo
         enddo
      endif

c     Transform coordinates
      do i=1,npoints

         if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >        ( abs(yy0(i)-mdv).gt.eps ).and.
     >        ( abs(zz0(i)-mdv).gt.eps ) )
     >    then
               call get_index3 (rid,rjd,rkd,xx0(i),yy0(i),zz0(i),1,
     >                          z3,zb,nx,ny,nz,xmin,ymin,dx,dy)

               zz0(i) = int_index3(p3,nx,ny,nz,rid,rjd,rkd,mdv)
          else
               xx0(i) = mdv
               yy0(i) = mdv
               zz0(i) = mdv
          endif

      enddo

c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)

c     Write output to trajectory
      else
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'p'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif

      endif

c     ------------------------------------------------------------
c     Convert single x/y/z to x/y/p
c     ------------------------------------------------------------

      if ( mode.eq.'-z2p.single' ) then

c     Transform coordinates - change from descending to ascending (minus sign)
      call get_index3 (rid,rjd,rkd,xpos,ypos,zpos,1,
     >                 z3,zb,nx,ny,nz,xmin,ymin,dx,dy)

      ppos = int_index3(p3,nx,ny,nz,rid,rjd,rkd,mdv)

c         Write result to screen
      print*,xpos,ypos,ppos

      endif

      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,outc,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),outc(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
      character ch

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
            outc(tab_x(i),tab_y(i))=outc(tab_x(i),tab_y(i))+
     >                             1.*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
               outc(nx,tab_y(i))=outc(nx,tab_y(i))+
     >                             1.*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,parameter ::       re=6371.229
      real,parameter ::       rad2deg=180./3.14159265
      real,parameter ::       deg2rad=3.141592654/180.

      real         xp,yp,xq,yq,arg
      character*80 unit
      real         dlon,dlat,alpha,rlat1,ralt2

      if ( unit.eq.'km' ) then

         arg=sin(yp*deg2rad)*sin(yq*deg2rad)+
     >           cos(yp*deg2rad)*cos(yq*deg2rad)*cos((xp-xq)*deg2rad)
         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 )

      elseif ( unit.eq.'km.haversine' ) then

        dlat   =  (yp - yq)*deg2rad
        dlon   =  (xp - xq)*deg2rad
        rlat1   =  yp*deg2rad
        rlat2   =  yq*deg2rad

        alpha  = ( sin(0.5*dlat)**2 ) +
     >           ( sin(0.5*dlon)**2 )*cos(rlat1)*cos(rlat2)

        sdis = 4. * re * atan2( sqrt(alpha),1. + sqrt( 1. - alpha ) )

      endif  

      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=1.-0.65*r/radius 
      else
          weight=0.
      endif

      end


c     ********************************************************************
c     * INPUT / OUTPUT SUBROUTINES                                       *
c     ********************************************************************

c     --------------------------------------------------------------------
c     Subroutines to write the netcdf output file
c     --------------------------------------------------------------------

      subroutine writecdf2D(cdfname,
     >                      varname,arr,time,
     >                      dx,dy,xmin,ymin,nx,ny,nz,
     >                      crefile,crevar)

      IMPLICIT NONE

c     Declaration of input parameters
      character*80 cdfname,varname
      integer nx,ny,nz
      real arr(nx,ny,nz)
      real dx,dy,xmin,ymin
      real time
      integer crefile,crevar

c     Further variables
      real varmin(4),varmax(4),stag(4)
      integer ierr,cdfid,ndim,vardim(4)
      character*80 cstname
      real mdv
      integer i
      integer nvars
      character*80 vars(300)

C     Definitions
      varmin(1)=xmin
      varmin(2)=ymin
      varmin(3)=1050.
      varmax(1)=xmin+real(nx-1)*dx
      varmax(2)=ymin+real(ny-1)*dy
      varmax(3)=1050.
      cstname='no_constants_file'
      ndim=4
      vardim(1)=nx
      vardim(2)=ny
      vardim(3)=nz
      stag(1)=0.
      stag(2)=0.
      stag(3)=0.
      mdv=-999.98999

C     Create the file
      if (crefile.eq.0) then
         call cdfwopn(cdfname,cdfid,ierr)
         if (ierr.ne.0) goto 906
      else if (crefile.eq.1) then
         call crecdf(cdfname,cdfid,
     >        varmin,varmax,ndim,cstname,ierr)
         if (ierr.ne.0) goto 903
      endif

c     Check whether the variable is already on the file
      call getvars(cdfid,nvars,vars,ierr)
      if (ierr.ne.0) goto 906
      do i=1,nvars
         if ( vars(i).eq.varname ) crevar = 0
      enddo

c     Write the data to the netcdf file, and close the file
      if (crevar.eq.1) then
         call putdef(cdfid,varname,ndim,mdv,
     >        vardim,varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 904
      endif
      call putdat(cdfid,varname,time,0,arr,ierr)
      if (ierr.ne.0) goto 905
      call clscdf(cdfid,ierr)

      return

c     Error handling
 903  print*,'*** Problem to create netcdf file ***'
      stop
 904  print*,'*** Problem to write definition ***'
      stop
 905  print*,'*** Problem to write data ***'
      stop
 906  print*,'*** Problem to open netcdf file ***'
      stop

      end


c     -----------------------------------------------
c     Subroutine to read a netcdf output file
c     -----------------------------------------------

      subroutine readcdf2D(cdfname,varname,arr,time,
     >                     dx,dy,xmin,ymin,nx,ny,nz)

      implicit none

c     Declaration of input parameters
      character*80 cdfname,varname
      integer      nx,ny,nz
      real         arr(nx,ny,nz)
      real         dx,dy,xmin,ymin
      real         time

c     Internal variables for the netcdf file
      real         varmin(4),varmax(4),stag(4)
      integer      ierr,cdfid,ndim,vardim(4)
      real         mdv
      real         delx,dely
      real         times(500)
      integer      ntimes,flag

c     Numerical epsilon
      real         eps
      parameter    (eps=0.01)

c     Auxiliary variables
      integer      i,j,k

c     Initialize the ouput array
      do i=1,nx
         do j=1,ny
            do k=1,nz
               arr(i,j,k)=0.
            enddo
         enddo
      enddo

c     Try to open the netcdf file and to read the variable
      call cdfopn(cdfname,cdfid,ierr)
      if (ierr.ne.0) goto 801

c     Check whether the specifieed time is available
      call gettimes(cdfid,times,ntimes,ierr)
      if (ierr.ne.0) goto 801
      flag=0
      do i=1,ntimes
          if (abs(times(i)-time).lt.eps) flag=1
      enddo
      if (flag.eq.0) then
         print*,'No such time on cdf',time
         stop
      endif

c     Get the variable definition, and the grid
      call getdef(cdfid,varname,ndim,mdv,vardim,varmin,varmax,
     >     stag,ierr)
      if (ierr.ne.0) goto 801

c     Set parameters or read parameters - depending on nx,ny
      if ( (nx.eq.1).and.(ny.eq.1) ) then
         dx   = (varmax(1)-varmin(1))/real(vardim(1)-1)
         dy   = (varmax(2)-varmin(2))/real(vardim(2)-1)
         nx   = vardim(1)
         ny   = vardim(2)
         nz   = vardim(3)
         xmin = varmin(1)
         ymin = varmin(2)
      else
         call getdat(cdfid,varname,time,0,arr,ierr)
         if (ierr.ne.0) goto 801
      endif

c     Close data file
      call clscdf(cdfid,ierr)

      return
      
c     Exception handling
 801  print*,'Problems in reading netcdf file'
      stop
      
      end


c     -----------------------------------------------
c     Diffusive filter
c     -----------------------------------------------

      subroutine filt2d (a,af,nx,ny,fil,misdat,
     &                   iperx,ipery,ispol,inpol)

c     Apply a conservative diffusion operator onto the 2d field a,
c     with full missing data checking.
c
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
c     af    real   out  filtered array, dimensioned (nx,ny), can be
c                       equivalenced with array a in the calling routine
c     f1    real        workarray, dimensioned (nx+1,ny)
c     f2    real        workarray, dimensioned (nx,ny+1)
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
c                       corresponds to one application of the linear filter.
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
c                       the corresponding value is not available. The
c                       misdat-checking can be switched off with with misdat=0.
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
c
c     Christoph Schaer, 1993

c     argument declaration
      integer     nx,ny
      real        a(nx,ny),af(nx,ny),fil,misdat
      integer     iperx,ipery,inpol,ispol

c     local variable declaration
      integer     i,j,is
      real        fh
      real        f1(nx+1,ny),f2(nx,ny+1)

c     compute constant fh
      fh=0.125*fil

c     compute fluxes in x-direction
      if (misdat.eq.0.) then
        do j=1,ny
        do i=2,nx
          f1(i,j)=a(i-1,j)-a(i,j)
        enddo
        enddo
      else
        do j=1,ny
        do i=2,nx
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
            f1(i,j)=0.
          else
            f1(i,j)=a(i-1,j)-a(i,j)
          endif
        enddo
        enddo
      endif
      if (iperx.eq.1) then
c       do periodic boundaries in the x-direction
        do j=1,ny
          f1(1,j)=f1(nx,j)
          f1(nx+1,j)=f1(2,j)
        enddo
      else
c       set boundary-fluxes to zero
        do j=1,ny
          f1(1,j)=0.
          f1(nx+1,j)=0.
        enddo
      endif

c     compute fluxes in y-direction
      if (misdat.eq.0.) then
        do j=2,ny
        do i=1,nx
          f2(i,j)=a(i,j-1)-a(i,j)
        enddo
        enddo
      else
        do j=2,ny
        do i=1,nx
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
            f2(i,j)=0.
          else
            f2(i,j)=a(i,j-1)-a(i,j)
          endif
        enddo
        enddo
      endif
c     set boundary-fluxes to zero
      do i=1,nx
        f2(i,1)=0.
        f2(i,ny+1)=0.
      enddo
      if (ipery.eq.1) then
c       do periodic boundaries in the x-direction
        do i=1,nx
          f2(i,1)=f2(i,ny)
          f2(i,ny+1)=f2(i,2)
        enddo
      endif
      if (iperx.eq.1) then
        if (ispol.eq.1) then
c         do south-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
          enddo
        endif
        if (inpol.eq.1) then
c         do north-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
          enddo
        endif
      endif

c     compute flux-convergence -> filter
      if (misdat.eq.0.) then
        do j=1,ny
        do i=1,nx
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
        enddo
        enddo
      else
        do j=1,ny
        do i=1,nx
          if (a(i,j).eq.misdat) then
            af(i,j)=misdat
          else
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
          endif
        enddo
        enddo
      endif
      end