WebSVN – lagranto.wrf – Diff – /trunk/goodies/wrfmap.f


      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.)
 
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
      real,allocatable,dimension (:,:)  :: connect1 
      integer      connectval1 
      real,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 
 
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
 
          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)
 
      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
 
          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)
 
      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
 
          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)
 
      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
 
          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)
 
      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=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
         if (arg.lt.-1.) arg=-1.
         if (arg.gt.1.) arg=1.
         sdis=re*acos(arg)
         
      elseif ( unit.eq.'deg' ) then
 
         dlon = xp-xq
         if ( dlon.gt. 180. ) dlon = dlon - 360.
         if ( dlon.lt.-180. ) dlon = dlon + 360.
         sdis = sqrt( dlon**2 + (yp-yq)**2 )
 
      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
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Rev 8	Rev 11
1	`PROGRAM wrfmap`	1	`PROGRAM wrfmap`
2		2
3	`c **********************************************************************`	3	`c **********************************************************************`
4	`c * Transform between lon/lat and x/y coordinates *`	4	`c * Transform between lon/lat and x/y coordinates *`
5	`c * Michael Sprenger / Spring 2013 *`	5	`c * Michael Sprenger / Spring 2013 *`
6	`c **********************************************************************`	6	`c **********************************************************************`
7		7
8	`implicit none`	8	`implicit none`
9		9
10	`c ------------------------------------------------------------`	10	`c ------------------------------------------------------------`
11	`c Declaration of parameters and variables`	11	`c Declaration of parameters and variables`
12	`c ------------------------------------------------------------`	12	`c ------------------------------------------------------------`
13		13
14	`c Input parameters`	14	`c Input parameters`
15	`character*80 mode`	15	`character*80 mode`
16	`character*80 inpfile`	16	`character*80 inpfile`
17	`character*80 outfile`	17	`character*80 outfile`
18	`integer ntra,ntim,ncol`	18	`integer ntra,ntim,ncol`
19	`character*80 anglemode`	19	`character*80 anglemode`
20		20
21	`c Fixed parameters`	21	`c Fixed parameters`
22	`character*80 mapfile ! netCDF file with map transformation`	22	`character*80 mapfile ! netCDF file with map transformation`
23	`parameter (mapfile ='wrfmap.nc')`	23	`parameter (mapfile ='wrfmap.nc')`
24	`real mdv ! missing data value`	24	`real mdv ! missing data value`
25	`parameter (mdv=-999.)`	25	`parameter (mdv=-999.)`
26		26
27	`c Numerical grid`	27	`c Numerical grid`
28	`real xmin,xmax,ymin,ymax ! Domain size`	28	`real xmin,xmax,ymin,ymax ! Domain size`
29	`real dx,dy ! Horizontal resolution`	29	`real dx,dy ! Horizontal resolution`
30	`integer nx,ny,nz ! Grid dimensions`	30	`integer nx,ny,nz ! Grid dimensions`
31	`real,allocatable,dimension (:,:) :: lon,lat ! lon/lat on numerical grid`	31	`real,allocatable,dimension (:,:) :: lon,lat ! lon/lat on numerical grid`
32	`real,allocatable,dimension (:,:) :: mpx,mpy ! map scale factors in x/y direction`	32	`real,allocatable,dimension (:,:) :: mpx,mpy ! map scale factors in x/y direction`
33		33
34		34
35	`c Geographical grid`	35	`c Geographical grid`
36	`real latmin,latmax,lonmin,lonmax ! Geographical domain`	36	`real latmin,latmax,lonmin,lonmax ! Geographical domain`
37	`real dlon,dlat ! Minimum spacing in geographical space`	37	`real dlon,dlat ! Minimum spacing in geographical space`
38	`integer nlon,nlat,nlev ! Geographical grid dimension`	38	`integer nlon,nlat,nlev ! Geographical grid dimension`
39	`real,allocatable,dimension (:,:) :: x,y ! x/y on geographical grid`	39	`real,allocatable,dimension (:,:) :: x,y ! x/y on geographical grid`
40		40
41	`c Vertical grid`	41	`c Vertical grid`
42	`real,allocatable,dimension (:,:,:) :: p3 ! 3D pressure`	42	`real,allocatable,dimension (:,:,:) :: p3 ! 3D pressure`
43	`real,allocatable,dimension (:,:) :: ps ! surface pressure`	43	`real,allocatable,dimension (:,:) :: ps ! surface pressure`
44	`real,allocatable,dimension (:,:,:) :: z3 ! 3D geopotential height`	44	`real,allocatable,dimension (:,:,:) :: z3 ! 3D geopotential height`
45	`real,allocatable,dimension (:,:) :: zb ! surface geopotential height`	45	`real,allocatable,dimension (:,:) :: zb ! surface geopotential height`
46		46
47	`c Transformation between numerical and geographical grid`	47	`c Transformation between numerical and geographical grid`
48	`real,allocatable,dimension (:,:) :: connect1`	48	`real,allocatable,dimension (:,:) :: connect1`
49	`integer connectval1`	49	`integer connectval1`
50	`real,allocatable,dimension (:,:) :: connect2`	50	`real,allocatable,dimension (:,:) :: connect2`
51	`integer connectval2`	51	`integer connectval2`
52	`real,allocatable,dimension (:,:) :: xc,yc`	52	`real,allocatable,dimension (:,:) :: xc,yc`
53	`real radius`	53	`real radius`
54	`real xval,yval`	54	`real xval,yval`
55		55
56	`c Trajectories`	56	`c Trajectories`
57	`real,allocatable, dimension (:,:,:) :: tra ! Input start coordinates (ntra,ntim,ncol)`	57	`real,allocatable, dimension (:,:,:) :: tra ! Input start coordinates (ntra,ntim,ncol)`
58	`integer refdate(6) ! Reference date`	58	`integer refdate(6) ! Reference date`
59	`character*80 vars(200) ! Field names`	59	`character*80 vars(200) ! Field names`
60	`real,allocatable, dimension (:) :: xx0,yy0,zz0 ! Position of air parcels`	60	`real,allocatable, dimension (:) :: xx0,yy0,zz0 ! Position of air parcels`
61	`integer inpmode`	61	`integer inpmode`
62	`integer outmode`	62	`integer outmode`
63	`integer npoints`	63	`integer npoints`
64		64
65	`c Auxiliary variables`	65	`c Auxiliary variables`
66	`integer fid`	66	`integer fid`
67	`integer stat`	67	`integer stat`
68	`character*80 varname,cdfname`	68	`character*80 varname,cdfname`
69	`integer i,j`	69	`integer i,j`
70	`character*80 radunit`	70	`character*80 radunit`
71	`real rdummy`	71	`real rdummy`
72	`real rid,rjd,rkd`	72	`real rid,rjd,rkd`
73	`real xpos,ypos,zpos,ppos,lonpos,latpos`	73	`real xpos,ypos,zpos,ppos,lonpos,latpos`
74	`real lon1,lat1,lon2,lat2`	74	`real lon1,lat1,lon2,lat2`
-		75	`integer nfilter`
75		76
76	`c Externals`	77	`c Externals`
77	`real int_index3`	78	`real int_index3`
78	`external int_index3`	79	`external int_index3`
79	`real sdis`	80	`real sdis`
80	`external sdis`	81	`external sdis`
81		82
82	`c ------------------------------------------------------------`	83	`c ------------------------------------------------------------`
83	`c Preparations`	84	`c Preparations`
84	`c ------------------------------------------------------------`	85	`c ------------------------------------------------------------`
85		86
86	`c Read parameters`	87	`c Read parameters`
87	`open(10,file='wrfmap.param')`	88	`open(10,file='wrfmap.param')`
88		89
89	`read(10,*) mode`	90	`read(10,*) mode`
90		91
91	`if ( mode.eq.'-create' ) then ! create mapping file`	92	`if ( mode.eq.'-create' ) then ! create mapping file`
92	`read(10,*) inpfile`	93	`read(10,*) inpfile`
93	`read(10,*) anglemode`	94	`read(10,*) anglemode`
94	`endif`	95	`endif`
95		96
96	`if ( mode.eq.'-ll2xy' ) then ! convert lon/lat -> x/y`	97	`if ( mode.eq.'-ll2xy' ) then ! convert lon/lat -> x/y`
97	`read(10,*) inpfile`	98	`read(10,*) inpfile`
98	`read(10,*) outfile`	99	`read(10,*) outfile`
99	`read(10,*) ntra,ntim,ncol`	100	`read(10,*) ntra,ntim,ncol`
100	`endif`	101	`endif`
101		102
102	`if ( mode.eq.'-ll2xy.single' ) then ! convert single lon/lat -> x/y`	103	`if ( mode.eq.'-ll2xy.single' ) then ! convert single lon/lat -> x/y`
103	`read(10,*) lonpos,latpos`	104	`read(10,*) lonpos,latpos`
104	`endif`	105	`endif`
105		106
106	`if ( mode.eq.'-xy2ll' ) then ! convert x/y -> lon/lat`	107	`if ( mode.eq.'-xy2ll' ) then ! convert x/y -> lon/lat`
107	`read(10,*) inpfile`	108	`read(10,*) inpfile`
108	`read(10,*) outfile`	109	`read(10,*) outfile`
109	`read(10,*) ntra,ntim,ncol`	110	`read(10,*) ntra,ntim,ncol`
110	`endif`	111	`endif`
111		112
112	`if ( mode.eq.'-xy2ll.single' ) then ! convert single x/y -> lon/lat`	113	`if ( mode.eq.'-xy2ll.single' ) then ! convert single x/y -> lon/lat`
113	`read(10,*) xpos,ypos`	114	`read(10,*) xpos,ypos`
114	`endif`	115	`endif`
115		116
116	`if ( mode.eq.'-p2z' ) then ! convert x/y/p -> x/y/z`	117	`if ( mode.eq.'-p2z' ) then ! convert x/y/p -> x/y/z`
117	`read(10,*) inpfile`	118	`read(10,*) inpfile`
118	`read(10,*) outfile`	119	`read(10,*) outfile`
119	`read(10,*) ntra,ntim,ncol`	120	`read(10,*) ntra,ntim,ncol`
120	`endif`	121	`endif`
121		122
122	`if ( mode.eq.'-p2z.single' ) then ! convert single x/y/p -> x/y/z`	123	`if ( mode.eq.'-p2z.single' ) then ! convert single x/y/p -> x/y/z`
123	`read(10,*) xpos,ypos,ppos`	124	`read(10,*) xpos,ypos,ppos`
124	`endif`	125	`endif`
125		126
126	`if ( mode.eq.'-z2p' ) then ! convert x/y/z -> x/y/p`	127	`if ( mode.eq.'-z2p' ) then ! convert x/y/z -> x/y/p`
127	`read(10,*) inpfile`	128	`read(10,*) inpfile`
128	`read(10,*) outfile`	129	`read(10,*) outfile`
129	`read(10,*) ntra,ntim,ncol`	130	`read(10,*) ntra,ntim,ncol`
130	`endif`	131	`endif`
131		132
132	`if ( mode.eq.'-z2p.single' ) then ! convert single x/y/z -> x/y/p`	133	`if ( mode.eq.'-z2p.single' ) then ! convert single x/y/z -> x/y/p`
133	`read(10,*) xpos,ypos,zpos`	134	`read(10,*) xpos,ypos,zpos`
134	`endif`	135	`endif`
135		136
-		137	`if ( mode.eq.'-mapscale' ) then ! calculate mapping scale factors`
-		138	`read(10,*) nfilter`
-		139	`endif`
-		140
136	`close(10)`	141	`close(10)`
137		142
138	`c Read the mapping file - except for mode '-create'`	143	`c Read the mapping file - except for mode '-create'`
139	`if ( mode.ne.'-create' ) then`	144	`if ( mode.ne.'-create' ) then`
140		145
141	`c Get dimensions`	146	`c Get dimensions`
142	`cdfname = mapfile`	147	`cdfname = mapfile`
143	`varname = 'Z'`	148	`varname = 'Z'`
144	`nx = 1`	149	`nx = 1`
145	`ny = 1`	150	`ny = 1`
146	`nz = 1`	151	`nz = 1`
147	`call readcdf2D(cdfname,varname,rdummy,0.,`	152	`call readcdf2D(cdfname,varname,rdummy,0.,`
148	`> dx,dy,xmin,ymin,nx,ny,nz)`	153	`> dx,dy,xmin,ymin,nx,ny,nz)`
149		154
150	`cdfname = mapfile`	155	`cdfname = mapfile`
151	`varname = 'X'`	156	`varname = 'X'`
152	`nlon = 1`	157	`nlon = 1`
153	`nlat = 1`	158	`nlat = 1`
154	`nlev = 1`	159	`nlev = 1`
155	`call readcdf2D(cdfname,varname,rdummy,0.,`	160	`call readcdf2D(cdfname,varname,rdummy,0.,`
156	`> dlon,dlat,lonmin,latmin,nlon,nlat,nlev)`	161	`> dlon,dlat,lonmin,latmin,nlon,nlat,nlev)`
157		162
158	`c Allocate memory`	163	`c Allocate memory`
159	`allocate(lon(nx,ny),stat=stat)`	164	`allocate(lon(nx,ny),stat=stat)`
160	`if (stat.ne.0) print,' error allocating array lon *'`	165	`if (stat.ne.0) print,' error allocating array lon *'`
161	`allocate(lat(nx,ny),stat=stat)`	166	`allocate(lat(nx,ny),stat=stat)`
162	`if (stat.ne.0) print,' error allocating array lat *'`	167	`if (stat.ne.0) print,' error allocating array lat *'`
163	`allocate(x(nlon,nlat),stat=stat)`	168	`allocate(x(nlon,nlat),stat=stat)`
164	`if (stat.ne.0) print,' error allocating array x *'`	169	`if (stat.ne.0) print,' error allocating array x *'`
165	`allocate(y(nlon,nlat),stat=stat)`	170	`allocate(y(nlon,nlat),stat=stat)`
166	`if (stat.ne.0) print,' error allocating array y *'`	171	`if (stat.ne.0) print,' error allocating array y *'`
167	`allocate(p3(nx,ny,nz),stat=stat)`	172	`allocate(p3(nx,ny,nz),stat=stat)`
168	`if (stat.ne.0) print,' error allocating array p3 *'`	173	`if (stat.ne.0) print,' error allocating array p3 *'`
169	`allocate(z3(nx,ny,nz),stat=stat)`	174	`allocate(z3(nx,ny,nz),stat=stat)`
170	`if (stat.ne.0) print,' error allocating array z3 *'`	175	`if (stat.ne.0) print,' error allocating array z3 *'`
171	`allocate(zb(nx,ny),stat=stat)`	176	`allocate(zb(nx,ny),stat=stat)`
172	`if (stat.ne.0) print,' error allocating array zb *'`	177	`if (stat.ne.0) print,' error allocating array zb *'`
173	`allocate(ps(nx,ny),stat=stat)`	178	`allocate(ps(nx,ny),stat=stat)`
174	`if (stat.ne.0) print,' error allocating array ps *'`	179	`if (stat.ne.0) print,' error allocating array ps *'`
175		180
176	`c Read data`	181	`c Read data`
177	`cdfname = mapfile`	182	`cdfname = mapfile`
178	`varname = 'LON'`	183	`varname = 'LON'`
179	`call readcdf2D(cdfname,varname,lon,0.,`	184	`call readcdf2D(cdfname,varname,lon,0.,`
180	`> dx,dy,xmin,ymin,nx,ny,1)`	185	`> dx,dy,xmin,ymin,nx,ny,1)`
181	`cdfname = mapfile`	186	`cdfname = mapfile`
182	`varname = 'LAT'`	187	`varname = 'LAT'`
183	`call readcdf2D(cdfname,varname,lat,0.,`	188	`call readcdf2D(cdfname,varname,lat,0.,`
184	`> dx,dy,xmin,ymin,nx,ny,1)`	189	`> dx,dy,xmin,ymin,nx,ny,1)`
185	`cdfname = mapfile`	190	`cdfname = mapfile`
186	`varname = 'X'`	191	`varname = 'X'`
187	`call readcdf2D(cdfname,varname,x,0.,`	192	`call readcdf2D(cdfname,varname,x,0.,`
188	`> dlon,dlat,lonmin,latmin,nlon,nlat,1)`	193	`> dlon,dlat,lonmin,latmin,nlon,nlat,1)`
189	`cdfname = mapfile`	194	`cdfname = mapfile`
190	`varname = 'Y'`	195	`varname = 'Y'`
191	`call readcdf2D(cdfname,varname,y,0.,`	196	`call readcdf2D(cdfname,varname,y,0.,`
192	`> dlon,dlat,lonmin,latmin,nlon,nlat,1)`	197	`> dlon,dlat,lonmin,latmin,nlon,nlat,1)`
193	`cdfname = mapfile`	198	`cdfname = mapfile`
194	`varname = 'Z'`	199	`varname = 'Z'`
195	`call readcdf2D(cdfname,varname,z3,0.,`	200	`call readcdf2D(cdfname,varname,z3,0.,`
196	`> dx,dy,xmin,ymin,nx,ny,nz)`	201	`> dx,dy,xmin,ymin,nx,ny,nz)`
197	`cdfname = mapfile`	202	`cdfname = mapfile`
198	`varname = 'P'`	203	`varname = 'P'`
199	`call readcdf2D(cdfname,varname,p3,0.,`	204	`call readcdf2D(cdfname,varname,p3,0.,`
200	`> dx,dy,xmin,ymin,nx,ny,nz)`	205	`> dx,dy,xmin,ymin,nx,ny,nz)`
201	`cdfname = mapfile`	206	`cdfname = mapfile`
202	`varname = 'TOPO'`	207	`varname = 'TOPO'`
203	`call readcdf2D(cdfname,varname,zb,0.,`	208	`call readcdf2D(cdfname,varname,zb,0.,`
204	`> dx,dy,xmin,ymin,nx,ny,1)`	209	`> dx,dy,xmin,ymin,nx,ny,1)`
205		210
206	`c Set surface prtessure to lowest 3d pressure level`	211	`c Set surface prtessure to lowest 3d pressure level`
207	`do i=1,nx`	212	`do i=1,nx`
208	`do j=1,ny`	213	`do j=1,ny`
209	`ps(i,j) = p3(i,j,1)`	214	`ps(i,j) = p3(i,j,1)`
210	`enddo`	215	`enddo`
211	`enddo`	216	`enddo`
212		217
213	`endif`	218	`endif`
214		219
215	`c ------------------------------------------------------------`	220	`c ------------------------------------------------------------`
216	`c Create mapping file`	221	`c Create mapping file`
217	`c ------------------------------------------------------------`	222	`c ------------------------------------------------------------`
218		223
219	`if ( mode.eq.'-create' ) then`	224	`if ( mode.eq.'-create' ) then`
220		225
221	`c Open mapping file`	226	`c Open mapping file`
222	`call input_open(fid,inpfile)`	227	`call input_open(fid,inpfile)`
223		228
224	`c Get grid description`	229	`c Get grid description`
225	`call input_grid_wrf(fid,`	230	`call input_grid_wrf(fid,`
226	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`	231	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`
227		232
228	`c Write grid information to screen`	233	`c Write grid information to screen`
229	`print*,' xmin,xmax = ',xmin,xmax`	234	`print*,' xmin,xmax = ',xmin,xmax`
230	`print*,' ymin,ymax = ',ymin,ymax`	235	`print*,' ymin,ymax = ',ymin,ymax`
231	`print*,' dx,dy = ',dx,dy`	236	`print*,' dx,dy = ',dx,dy`
232	`print*,' nx,ny,nz = ',nx,ny,nz`	237	`print*,' nx,ny,nz = ',nx,ny,nz`
233		238
234	`c Transform grid specifier into grid coordinates`	239	`c Transform grid specifier into grid coordinates`
235	`xmin = 1.`	240	`xmin = 1.`
236	`ymin = 1.`	241	`ymin = 1.`
237	`dx = 1.`	242	`dx = 1.`
238	`dy = 1.`	243	`dy = 1.`
239		244
240	`c Allocate memory for lon/lat on numerical grid`	245	`c Allocate memory for lon/lat on numerical grid`
241	`allocate(lon(nx,ny),stat=stat)`	246	`allocate(lon(nx,ny),stat=stat)`
242	`if (stat.ne.0) print,' error allocating array lon *'`	247	`if (stat.ne.0) print,' error allocating array lon *'`
243	`allocate(lat(nx,ny),stat=stat)`	248	`allocate(lat(nx,ny),stat=stat)`
244	`if (stat.ne.0) print,' error allocating array lat *'`	249	`if (stat.ne.0) print,' error allocating array lat *'`
245	`allocate(p3(nx,ny,nz),stat=stat)`	250	`allocate(p3(nx,ny,nz),stat=stat)`
246	`if (stat.ne.0) print,' error allocating array p3 *'`	251	`if (stat.ne.0) print,' error allocating array p3 *'`
247	`allocate(z3(nx,ny,nz),stat=stat)`	252	`allocate(z3(nx,ny,nz),stat=stat)`
248	`if (stat.ne.0) print,' error allocating array z3 *'`	253	`if (stat.ne.0) print,' error allocating array z3 *'`
249	`allocate(zb(nx,ny),stat=stat)`	254	`allocate(zb(nx,ny),stat=stat)`
250	`if (stat.ne.0) print,' error allocating array zb *'`	255	`if (stat.ne.0) print,' error allocating array zb *'`
251		256
252	`c Read lon/lat on the model grid`	257	`c Read lon/lat on the model grid`
253	`varname='XLONG'`	258	`varname='XLONG'`
254	`call input_var_wrf(fid,varname,lon,nx,ny,1)`	259	`call input_var_wrf(fid,varname,lon,nx,ny,1)`
255	`varname='XLAT'`	260	`varname='XLAT'`
256	`call input_var_wrf(fid,varname,lat,nx,ny,1)`	261	`call input_var_wrf(fid,varname,lat,nx,ny,1)`
257	`varname='pressure'`	262	`varname='pressure'`
258	`call input_var_wrf(fid,varname,p3 ,nx,ny,nz)`	263	`call input_var_wrf(fid,varname,p3 ,nx,ny,nz)`
259	`varname='geopot'`	264	`varname='geopot'`
260	`call input_var_wrf(fid,varname,z3 ,nx,ny,nz)`	265	`call input_var_wrf(fid,varname,z3 ,nx,ny,nz)`
261	`varname='geopots'`	266	`varname='geopots'`
262	`call input_var_wrf(fid,varname,zb ,nx,ny,1)`	267	`call input_var_wrf(fid,varname,zb ,nx,ny,1)`
263		268
264	`c Transform longitude depending on <anglemode>`	269	`c Transform longitude depending on <anglemode>`
265	`if ( anglemode.eq.'greenwich' ) then`	270	`if ( anglemode.eq.'greenwich' ) then`
266	`do i=1,nx`	271	`do i=1,nx`
267	`do j=1,ny`	272	`do j=1,ny`
268	`if ( lon(i,j).lt.0. ) lon(i,j) = lon(i,j) + 360.`	273	`if ( lon(i,j).lt.0. ) lon(i,j) = lon(i,j) + 360.`
269	`enddo`	274	`enddo`
270	`enddo`	275	`enddo`
271	`elseif ( anglemode.eq.'dateline' ) then`	276	`elseif ( anglemode.eq.'dateline' ) then`
272	`do i=1,nx`	277	`do i=1,nx`
273	`do j=1,ny`	278	`do j=1,ny`
274	`if ( lon(i,j).gt.180. ) lon(i,j) = lon(i,j) - 360.`	279	`if ( lon(i,j).gt.180. ) lon(i,j) = lon(i,j) - 360.`
275	`enddo`	280	`enddo`
276	`enddo`	281	`enddo`
277	`else`	282	`else`
278	`print*,' ERROR: unsupported angle mode... Stop'`	283	`print*,' ERROR: unsupported angle mode... Stop'`
279	`stop`	284	`stop`
280	`endif`	285	`endif`
281		286
282	`c Close input file file`	287	`c Close input file file`
283	`call input_close(fid)`	288	`call input_close(fid)`
284		289
285	`c Check for 'date line' and for pole`	290	`c Check for 'date line' and for pole`
286	`do i=2,nx`	291	`do i=2,nx`
287	`do j=1,ny`	292	`do j=1,ny`
288	`if ( abs( lat(i,j) ).gt.90. ) then`	293	`if ( abs( lat(i,j) ).gt.90. ) then`
289	`print*,' ERROR: mapping over pole not supported... Stop'`	294	`print*,' ERROR: mapping over pole not supported... Stop'`
290	`stop`	295	`stop`
291	`endif`	296	`endif`
292	`if ( abs( lon(i,j)-lon(i-1,j) ).gt.180. ) then`	297	`if ( abs( lon(i,j)-lon(i-1,j) ).gt.180. ) then`
293	`print*,' ERROR: mapping over date line not supported... Stop'`	298	`print*,' ERROR: mapping over date line not supported... Stop'`
294	`stop`	299	`stop`
295	`endif`	300	`endif`
296	`enddo`	301	`enddo`
297	`enddo`	302	`enddo`
298		303
299	`c Determine the extreme coordinates`	304	`c Determine the extreme coordinates`
300	`latmin = lat(1,1)`	305	`latmin = lat(1,1)`
301	`latmax = lat(1,1)`	306	`latmax = lat(1,1)`
302	`lonmin = lon(1,1)`	307	`lonmin = lon(1,1)`
303	`lonmax = lon(1,1)`	308	`lonmax = lon(1,1)`
304	`do i=1,nx`	309	`do i=1,nx`
305	`do j=1,ny`	310	`do j=1,ny`
306	`if ( lat(i,j).lt.latmin ) latmin = lat(i,j)`	311	`if ( lat(i,j).lt.latmin ) latmin = lat(i,j)`
307	`if ( lat(i,j).gt.latmax ) latmax = lat(i,j)`	312	`if ( lat(i,j).gt.latmax ) latmax = lat(i,j)`
308	`if ( lon(i,j).lt.lonmin ) lonmin = lon(i,j)`	313	`if ( lon(i,j).lt.lonmin ) lonmin = lon(i,j)`
309	`if ( lon(i,j).gt.lonmax ) lonmax = lon(i,j)`	314	`if ( lon(i,j).gt.lonmax ) lonmax = lon(i,j)`
310	`enddo`	315	`enddo`
311	`enddo`	316	`enddo`
312		317
313	`print*,' lonmin,max = ',lonmin,lonmax`	318	`print*,' lonmin,max = ',lonmin,lonmax`
314	`print*,' latmin,max = ',latmin,latmax`	319	`print*,' latmin,max = ',latmin,latmax`
315		320
316	`c Determine the extreme dlon/dlat spacing`	321	`c Determine the extreme dlon/dlat spacing`
317	`dlon = abs( lon(2,1)-lon(1,1) )`	322	`dlon = abs( lon(2,1)-lon(1,1) )`
318	`do i=2,nx`	323	`do i=2,nx`
319	`do j=1,ny`	324	`do j=1,ny`
320	`if ( abs( lon(i,j)-lon(i-1,j) ).lt.dlon ) then`	325	`if ( abs( lon(i,j)-lon(i-1,j) ).lt.dlon ) then`
321	`dlon = abs( lon(i,j)-lon(i-1,j) )`	326	`dlon = abs( lon(i,j)-lon(i-1,j) )`
322	`endif`	327	`endif`
323	`enddo`	328	`enddo`
324	`enddo`	329	`enddo`
325	`dlat = abs( lat(1,2)-lat(1,1) )`	330	`dlat = abs( lat(1,2)-lat(1,1) )`
326	`do i=1,nx`	331	`do i=1,nx`
327	`do j=2,ny`	332	`do j=2,ny`
328	`if ( abs( lat(i,j)-lat(i,j-1) ).lt.dlat ) then`	333	`if ( abs( lat(i,j)-lat(i,j-1) ).lt.dlat ) then`
329	`dlat = abs( lat(i,j)-lat(i,j-1) )`	334	`dlat = abs( lat(i,j)-lat(i,j-1) )`
330	`endif`	335	`endif`
331	`enddo`	336	`enddo`
332	`enddo`	337	`enddo`
333		338
334	`c Set dimension of geographical grid`	339	`c Set dimension of geographical grid`
335	`nlon = ceiling( (lonmax-lonmin) / dlon + 1. )`	340	`nlon = ceiling( (lonmax-lonmin) / dlon + 1. )`
336	`nlat = ceiling( (latmax-latmin) / dlat + 1. )`	341	`nlat = ceiling( (latmax-latmin) / dlat + 1. )`
337	`lonmax = lonmin + real(nlon-1) * dlon`	342	`lonmax = lonmin + real(nlon-1) * dlon`
338	`latmax = latmin + real(nlat-1) * dlat`	343	`latmax = latmin + real(nlat-1) * dlat`
339		344
340	`print*,' dlon,dlat = ',dlon,dlat`	345	`print*,' dlon,dlat = ',dlon,dlat`
341	`print*,' nlon,nlat = ',nlon,nlat`	346	`print*,' nlon,nlat = ',nlon,nlat`
342		347
343	`c Allocate memory for x/y on geographical grid`	348	`c Allocate memory for x/y on geographical grid`
344	`allocate(x(nlon,nlat),stat=stat)`	349	`allocate(x(nlon,nlat),stat=stat)`
345	`if (stat.ne.0) print,' error allocating array x *'`	350	`if (stat.ne.0) print,' error allocating array x *'`
346	`allocate(y(nlon,nlat),stat=stat)`	351	`allocate(y(nlon,nlat),stat=stat)`
347	`if (stat.ne.0) print,' error allocating array y *'`	352	`if (stat.ne.0) print,' error allocating array y *'`
348	`allocate(connect1(nlon,nlat),stat=stat)`	353	`allocate(connect1(nlon,nlat),stat=stat)`
349	`if (stat.ne.0) print,' error allocating array connect1*'`	354	`if (stat.ne.0) print,' error allocating array connect1*'`
350	`allocate(connect2(nlon,nlat),stat=stat)`	355	`allocate(connect2(nlon,nlat),stat=stat)`
351	`if (stat.ne.0) print,' error allocating array connect2*'`	356	`if (stat.ne.0) print,' error allocating array connect2*'`
352	`allocate(xc(nlon,nlat),stat=stat)`	357	`allocate(xc(nlon,nlat),stat=stat)`
353	`if (stat.ne.0) print,' error allocating array xc *'`	358	`if (stat.ne.0) print,' error allocating array xc *'`
354	`allocate(yc(nlon,nlat),stat=stat)`	359	`allocate(yc(nlon,nlat),stat=stat)`
355	`if (stat.ne.0) print,' error allocating array yc *'`	360	`if (stat.ne.0) print,' error allocating array yc *'`
356		361
357	`c Init the mapping arrays`	362	`c Init the mapping arrays`
358	`connectval1 = 0`	363	`connectval1 = 0`
359	`connectval2 = 0`	364	`connectval2 = 0`
360	`do i=1,nlon`	365	`do i=1,nlon`
361	`do j=1,nlat`	366	`do j=1,nlat`
362	`x(i,j) = 0.`	367	`x(i,j) = 0.`
363	`y(i,j) = 0.`	368	`y(i,j) = 0.`
364	`xc(i,j) = 0.`	369	`xc(i,j) = 0.`
365	`yc(i,j) = 0.`	370	`yc(i,j) = 0.`
366	`connect1(i,j) = 1`	371	`connect1(i,j) = 1`
367	`connect2(i,j) = 1`	372	`connect2(i,j) = 1`
368	`enddo`	373	`enddo`
369	`enddo`	374	`enddo`
370		375
371	`c Get the good radius for gridding - avoiding gaps`	376	`c Get the good radius for gridding - avoiding gaps`
372	`radunit = 'deg'`	377	`radunit = 'deg'`
373	`radius = 0.`	378	`radius = 0.`
374	`do i=2,nx`	379	`do i=2,nx`
375	`do j=2,ny`	380	`do j=2,ny`
376	`if ( abs( lat(i,j)-lat(i,j-1) ).gt.radius ) then`	381	`if ( abs( lat(i,j)-lat(i,j-1) ).gt.radius ) then`
377	`radius = abs( lat(i,j)-lat(i,j-1) )`	382	`radius = abs( lat(i,j)-lat(i,j-1) )`
378	`endif`	383	`endif`
379	`if ( abs( lon(i,j)-lon(i-1,j) ).gt.radius ) then`	384	`if ( abs( lon(i,j)-lon(i-1,j) ).gt.radius ) then`
380	`radius = abs( lon(i,j)-lon(i-1,j) )`	385	`radius = abs( lon(i,j)-lon(i-1,j) )`
381	`endif`	386	`endif`
382	`enddo`	387	`enddo`
383	`enddo`	388	`enddo`
384	`radius = 3. * radius`	389	`radius = 3. * radius`
385		390
386	`c Do the mapping`	391	`c Do the mapping`
387	`do i=1,nx`	392	`do i=1,nx`
388	`do j=1,ny`	393	`do j=1,ny`
389		394
390	`connectval1 = connectval1 + 1`	395	`connectval1 = connectval1 + 1`
391	`call gridding1(lat(i,j),lon(i,j),real(i),radius,radunit,`	396	`call gridding1(lat(i,j),lon(i,j),real(i),radius,radunit,`
392	`> connect1,connectval1,`	397	`> connect1,connectval1,`
393	`> x,xc,nlon,nlat,lonmin,latmin,dlon,dlat)`	398	`> x,xc,nlon,nlat,lonmin,latmin,dlon,dlat)`
394		399
395	`connectval2 = connectval2 + 1`	400	`connectval2 = connectval2 + 1`
396	`call gridding1(lat(i,j),lon(i,j),real(j),radius,radunit,`	401	`call gridding1(lat(i,j),lon(i,j),real(j),radius,radunit,`
397	`> connect2,connectval2,`	402	`> connect2,connectval2,`
398	`> y,yc,nlon,nlat,lonmin,latmin,dlon,dlat)`	403	`> y,yc,nlon,nlat,lonmin,latmin,dlon,dlat)`
399		404
400		405
401	`enddo`	406	`enddo`
402	`enddo`	407	`enddo`
403		408
404	`c Normalize output and set miising data flag`	409	`c Normalize output and set miising data flag`
405	`do i=1,nlon`	410	`do i=1,nlon`
406	`do j=1,nlat`	411	`do j=1,nlat`
407	`if ( xc(i,j).gt.0. ) then`	412	`if ( xc(i,j).gt.0. ) then`
408	`x(i,j) = x(i,j)/xc(i,j)`	413	`x(i,j) = x(i,j)/xc(i,j)`
409	`else`	414	`else`
410	`x(i,j) = mdv`	415	`x(i,j) = mdv`
411	`endif`	416	`endif`
412	`if ( yc(i,j).gt.0. ) then`	417	`if ( yc(i,j).gt.0. ) then`
413	`y(i,j) = y(i,j)/yc(i,j)`	418	`y(i,j) = y(i,j)/yc(i,j)`
414	`else`	419	`else`
415	`y(i,j) = mdv`	420	`y(i,j) = mdv`
416	`endif`	421	`endif`
417	`enddo`	422	`enddo`
418	`enddo`	423	`enddo`
419		424
420	`c Write data to netCDF`	425	`c Write data to netCDF`
421	`cdfname = mapfile`	426	`cdfname = mapfile`
422	`varname = 'X'`	427	`varname = 'X'`
423	`call writecdf2D(cdfname,varname,x,0.,`	428	`call writecdf2D(cdfname,varname,x,0.,`
424	`> dlon,dlat,lonmin,latmin,nlon,nlat,1,1,1)`	429	`> dlon,dlat,lonmin,latmin,nlon,nlat,1,1,1)`
425		430
426	`cdfname = mapfile`	431	`cdfname = mapfile`
427	`varname = 'Y'`	432	`varname = 'Y'`
428	`call writecdf2D(cdfname,varname,y,0.,`	433	`call writecdf2D(cdfname,varname,y,0.,`
429	`> dlon,dlat,lonmin,latmin,nlon,nlat,1,0,1)`	434	`> dlon,dlat,lonmin,latmin,nlon,nlat,1,0,1)`
430		435
431	`cdfname = mapfile`	436	`cdfname = mapfile`
432	`varname = 'LON'`	437	`varname = 'LON'`
433	`call writecdf2D(cdfname,varname,lon,0.,`	438	`call writecdf2D(cdfname,varname,lon,0.,`
434	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`	439	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`
435		440
436	`cdfname = mapfile`	441	`cdfname = mapfile`
437	`varname = 'LAT'`	442	`varname = 'LAT'`
438	`call writecdf2D(cdfname,varname,lat,0.,`	443	`call writecdf2D(cdfname,varname,lat,0.,`
439	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`	444	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`
440		445
441	`cdfname = mapfile`	446	`cdfname = mapfile`
442	`varname = 'Z'`	447	`varname = 'Z'`
443	`call writecdf2D(cdfname,varname,z3,0.,`	448	`call writecdf2D(cdfname,varname,z3,0.,`
444	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`	449	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`
445		450
446	`cdfname = mapfile`	451	`cdfname = mapfile`
447	`varname = 'P'`	452	`varname = 'P'`
448	`call writecdf2D(cdfname,varname,p3,0.,`	453	`call writecdf2D(cdfname,varname,p3,0.,`
449	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`	454	`> dx,dy,xmin,ymin,nx,ny,nz,0,1)`
450		455
451	`cdfname = mapfile`	456	`cdfname = mapfile`
452	`varname = 'TOPO'`	457	`varname = 'TOPO'`
453	`call writecdf2D(cdfname,varname,zb,0.,`	458	`call writecdf2D(cdfname,varname,zb,0.,`
454	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`	459	`> dx,dy,xmin,ymin,nx,ny,1,0,1)`
455		460
456		461
457	`endif`	462	`endif`
458		463
459		464
460	`c ------------------------------------------------------------`	465	`c ------------------------------------------------------------`
461	`c Convert a lat/lon/z list to a x/y/z list`	466	`c Convert a lat/lon/z list to a x/y/z list`
462	`c ------------------------------------------------------------`	467	`c ------------------------------------------------------------`
463		468
464	`if ( mode.eq.'-ll2xy' ) then`	469	`if ( mode.eq.'-ll2xy' ) then`
465		470
466	`c Allocate memory for input and output lists`	471	`c Allocate memory for input and output lists`
467	`allocate(tra(ntra,ntim,ncol),stat=stat)`	472	`allocate(tra(ntra,ntim,ncol),stat=stat)`
468	`if (stat.ne.0) print,' error allocating array tra *'`	473	`if (stat.ne.0) print,' error allocating array tra *'`
469	`allocate(xx0(ntra*ntim),stat=stat)`	474	`allocate(xx0(ntra*ntim),stat=stat)`
470	`if (stat.ne.0) print,' error allocating array xx0 *'`	475	`if (stat.ne.0) print,' error allocating array xx0 *'`
471	`allocate(yy0(ntra*ntim),stat=stat)`	476	`allocate(yy0(ntra*ntim),stat=stat)`
472	`if (stat.ne.0) print,' error allocating array yy0 *'`	477	`if (stat.ne.0) print,' error allocating array yy0 *'`
473	`allocate(zz0(ntra*ntim),stat=stat)`	478	`allocate(zz0(ntra*ntim),stat=stat)`
474	`if (stat.ne.0) print,' error allocating array zz0 *'`	479	`if (stat.ne.0) print,' error allocating array zz0 *'`
475		480
476	`c Get format of input and output file`	481	`c Get format of input and output file`
477	`call mode_tra(inpmode,inpfile)`	482	`call mode_tra(inpmode,inpfile)`
478	`call mode_tra(outmode,outfile)`	483	`call mode_tra(outmode,outfile)`
479	`if ( outmode.eq.-1 ) outmode = 1`	484	`if ( outmode.eq.-1 ) outmode = inpmode`
480		485
481	`c Read coordinates from file - Format (lon,lat,lev)`	486	`c Read coordinates from file - Format (lon,lat,lev)`
482	`if (inpmode.eq.-1) then`	487	`if (inpmode.eq.-1) then`
483	`fid = 10`	488	`fid = 10`
484	`npoints = 0`	489	`npoints = 0`
485	`open(fid,file=inpfile)`	490	`open(fid,file=inpfile)`
486	`do i=1,ntra`	491	`do i=1,ntra`
487	`npoints = npoints + 1`	492	`npoints = npoints + 1`
488	`read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)`	493	`read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)`
489	`enddo`	494	`enddo`
490	`close(fid)`	495	`close(fid)`
491		496
492	`c Read coordinates from trajectory file`	497	`c Read coordinates from trajectory file`
493	`else`	498	`else`
494	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)`	499	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)`
495	`call read_tra (fid,tra,ntra,ntim,ncol,inpmode)`	500	`call read_tra (fid,tra,ntra,ntim,ncol,inpmode)`

Subversion Repositories lagranto.wrf

(root)/trunk/goodies/wrfmap.f – Rev 8 → 11