WebSVN – lagranto.wrf – Diff – /branches/met.nps.monterey/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
 
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
 
      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 = 1
 
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=1
 
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     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=1
 
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=1
 
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
 
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     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
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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

Subversion Repositories lagranto.wrf

(root)/branches/met.nps.monterey/goodies/wrfmap.f – Rev 8 → 10