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
      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
 
          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=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
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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

Subversion Repositories lagranto.wrf

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