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.)
      real         eps                              ! Numerical epsilon
      parameter    (eps= 0.001)
 
c     Numerical grid
      real         xmin,xmax,ymin,ymax              ! Domain size
      real         dx,dy                            ! Horizontal resolution
      integer      nx,ny,nz                         ! Grid dimensions
      real,allocatable,dimension (:,:)  :: lon,lat  ! lon/lat on numerical grid
      real,allocatable,dimension (:,:)  :: mpx,mpy  ! map scale factors in x/y direction
 
 
c     Geographical grid
      real         latmin,latmax,lonmin,lonmax      ! Geographical domain
      real         dlon,dlat                        ! Minimum spacing in geographical space
      integer      nlon,nlat,nlev                        ! Geographical grid dimension
      real,allocatable,dimension (:,:)   :: x,y      ! x/y on geographical grid
 
c     Vertical grid
      real,allocatable,dimension (:,:,:)   :: p3    ! 3D pressure
      real,allocatable,dimension (:,:)     :: ps    ! surface pressure
      real,allocatable,dimension (:,:,:)   :: z3    ! 3D geopotential height
      real,allocatable,dimension (:,:)     :: zb    ! surface geopotential height
 
c     Transformation between numerical and geographical grid
      integer,allocatable,dimension (:,:)  :: connect1 
      integer      connectval1 
      integer,allocatable,dimension (:,:)  :: connect2 
      integer      connectval2
      real,allocatable,dimension (:,:)  :: xc,yc
      real         radius
      real         xval,yval
 
c     Trajectories
      real,allocatable, dimension (:,:,:) :: tra             ! Input start coordinates (ntra,ntim,ncol)
      integer                                refdate(6)      ! Reference date
      character*80                           vars(200)       ! Field names
      real,allocatable, dimension (:)     :: xx0,yy0,zz0     ! Position of air parcels
      integer                                inpmode
      integer                                outmode
      integer                                npoints
 
c     Auxiliary variables
      integer      fid
      integer      stat
      character*80 varname,cdfname
      integer      i,j
      character*80 radunit
      real         rdummy
      real         rid,rjd,rkd
      real         xpos,ypos,zpos,ppos,lonpos,latpos
      real         lon1,lat1,lon2,lat2
      integer      nfilter
 
c     Externals
      real         int_index3
      external     int_index3
      real         sdis
      external     sdis
 
c     ------------------------------------------------------------
c     Preparations
c     ------------------------------------------------------------
 
c     Read parameters
      open(10,file='wrfmap.param')
 
       read(10,*) mode
 
       if ( mode.eq.'-create' ) then         ! create mapping file
          read(10,*) inpfile
          read(10,*) anglemode
       endif
 
       if ( mode.eq.'-ll2xy'  ) then         ! convert lon/lat -> x/y
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif
 
       if ( mode.eq.'-ll2xy.single'  ) then   ! convert single lon/lat -> x/y
          read(10,*) lonpos,latpos
       endif
 
       if ( mode.eq.'-xy2ll'  ) then         ! convert x/y -> lon/lat
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif
 
       if ( mode.eq.'-xy2ll.single'  ) then   ! convert single x/y -> lon/lat
          read(10,*) xpos,ypos
       endif
 
       if ( mode.eq.'-p2z'  ) then            ! convert x/y/p -> x/y/z
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif
 
       if ( mode.eq.'-p2z.single'  ) then     ! convert single x/y/p -> x/y/z
          read(10,*) xpos,ypos,ppos
       endif
 
       if ( mode.eq.'-z2p'  ) then            ! convert x/y/z -> x/y/p
          read(10,*) inpfile
          read(10,*) outfile
          read(10,*) ntra,ntim,ncol
       endif
 
       if ( mode.eq.'-z2p.single'  ) then     ! convert single x/y/z -> x/y/p
          read(10,*) xpos,ypos,zpos
       endif
 
       if ( mode.eq.'-mapscale'  ) then       ! calculate mapping scale factors
          read(10,*) nfilter
       endif
 
      close(10)
 
c     Read the mapping file - except for mode '-create'
      if ( mode.ne.'-create' ) then
 
c        Get dimensions 
         cdfname = mapfile
         varname = 'Z'
         nx      = 1
         ny      = 1
         nz      = 1
         call readcdf2D(cdfname,varname,rdummy,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)
 
         cdfname = mapfile
         varname = 'X'
         nlon    = 1
         nlat    = 1
         nlev    = 1
         call readcdf2D(cdfname,varname,rdummy,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,nlev)
 
c        Allocate memory
         allocate(lon(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array lon ***'
         allocate(lat(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array lat ***'
         allocate(x(nlon,nlat),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array x   ***'
         allocate(y(nlon,nlat),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array y   ***'
         allocate(p3(nx,ny,nz),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array p3  ***'
         allocate(z3(nx,ny,nz),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array z3  ***'
         allocate(zb(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array zb  ***'
         allocate(ps(nx,ny),stat=stat)
         if (stat.ne.0) print*,'*** error allocating array ps  ***'
 
c        Read data
         cdfname = mapfile
         varname = 'LON'
         call readcdf2D(cdfname,varname,lon,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)
         cdfname = mapfile
         varname = 'LAT'
         call readcdf2D(cdfname,varname,lat,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)
         cdfname = mapfile
         varname = 'X'
         call readcdf2D(cdfname,varname,x,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,1)
         cdfname = mapfile
         varname = 'Y'
         call readcdf2D(cdfname,varname,y,0.,
     >                  dlon,dlat,lonmin,latmin,nlon,nlat,1)
         cdfname = mapfile
         varname = 'Z'
         call readcdf2D(cdfname,varname,z3,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)
         cdfname = mapfile
         varname = 'P'
         call readcdf2D(cdfname,varname,p3,0.,
     >                  dx,dy,xmin,ymin,nx,ny,nz)
         cdfname = mapfile
         varname = 'TOPO'
         call readcdf2D(cdfname,varname,zb,0.,
     >                  dx,dy,xmin,ymin,nx,ny,1)
 
c	     Set surface prtessure to lowest 3d pressure level
         do i=1,nx
            do j=1,ny
               ps(i,j) = p3(i,j,1)
            enddo
         enddo
 
      endif
 
c     ------------------------------------------------------------
c     Create mapping file
c     ------------------------------------------------------------
 
      if ( mode.eq.'-create' ) then
 
c     Open mapping file
      call input_open(fid,inpfile)
 
c     Get grid description
      call input_grid_wrf(fid,
     >		          xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)
 
c     Write grid information to screen
      print*,' xmin,xmax  = ',xmin,xmax
      print*,' ymin,ymax  = ',ymin,ymax
      print*,' dx,dy      = ',dx,dy
      print*,' nx,ny,nz   = ',nx,ny,nz 
      print*,' anglemode  = ',trim(anglemode)
 
c     Transform grid specifier into grid coordinates
      xmin = 1.
      ymin = 1.
      dx   = 1.
      dy   = 1.
 
c     Allocate memory for lon/lat on numerical grid
      allocate(lon(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon ***'
      allocate(lat(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat ***'
      allocate(p3(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3  ***'
      allocate(z3(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array z3  ***'
      allocate(zb(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zb  ***'
 
c     Read lon/lat on the model grid
      varname='XLONG'
      call input_var_wrf(fid,varname,lon,nx,ny,1)
      varname='XLAT'
      call input_var_wrf(fid,varname,lat,nx,ny,1)
      varname='pressure'
      call input_var_wrf(fid,varname,p3 ,nx,ny,nz)
      varname='geopot'
      call input_var_wrf(fid,varname,z3 ,nx,ny,nz)
      varname='geopots'
      call input_var_wrf(fid,varname,zb ,nx,ny,1)
 
c     Transform longitude depending on <anglemode>
      if ( anglemode.eq.'greenwich' ) then
         do i=1,nx
           do j=1,ny
              if ( lon(i,j).lt.0. ) lon(i,j) = lon(i,j) + 360.
            enddo
         enddo
      elseif ( anglemode.eq.'dateline' ) then
         do i=1,nx
            do j=1,ny
               if ( lon(i,j).gt.180. ) lon(i,j) = lon(i,j) - 360.
            enddo
         enddo
      else
         print*,' ERROR: unsupported angle mode... Stop'
         stop
      endif
 
c     Close input file file
      call input_close(fid)
 
c     Check for 'date line' and for pole
      do i=2,nx
      do j=1,ny
        if ( abs( lat(i,j) ).gt.90. ) then
           print*,' ERROR: mapping over pole not supported... Stop'
           stop
        endif
        if ( abs( lon(i,j)-lon(i-1,j) ).gt.180. ) then
           print*,' ERROR: mapping over date line not supported... Stop'
           stop
        endif
      enddo
      enddo
 
c     Determine the extreme coordinates
      latmin  = lat(1,1)
      latmax  = lat(1,1)
      lonmin  = lon(1,1)
      lonmax  = lon(1,1)
      do i=1,nx
	    do j=1,ny
          if ( lat(i,j).lt.latmin ) latmin = lat(i,j)
          if ( lat(i,j).gt.latmax ) latmax = lat(i,j)
          if ( lon(i,j).lt.lonmin ) lonmin = lon(i,j)
          if ( lon(i,j).gt.lonmax ) lonmax = lon(i,j)
        enddo
      enddo
 
      print*,' lonmin,max = ',lonmin,lonmax
      print*,' latmin,max = ',latmin,latmax
 
c     Determine the extreme dlon/dlat spacing
      dlon = abs( lon(2,1)-lon(1,1) )
      do i=2,nx
	    do j=1,ny
           if ( abs( lon(i,j)-lon(i-1,j) ).lt.dlon ) then
              dlon = abs( lon(i,j)-lon(i-1,j) )
           endif
         enddo
      enddo
      dlat = abs( lat(1,2)-lat(1,1) )
      do i=1,nx
	    do j=2,ny
           if ( abs( lat(i,j)-lat(i,j-1) ).lt.dlat ) then
              dlat = abs( lat(i,j)-lat(i,j-1) )
           endif
         enddo
      enddo
 
c     Set dimension of geographical grid 
      nlon   = ceiling( (lonmax-lonmin) / dlon + 1. )
      nlat   = ceiling( (latmax-latmin) / dlat + 1. )
      lonmax = lonmin + real(nlon-1) * dlon
      latmax = latmin + real(nlat-1) * dlat
 
      print*,' dlon,dlat  = ',dlon,dlat
      print*,' nlon,nlat  = ',nlon,nlat
 
c     Allocate memory for x/y on geographical grid
      allocate(x(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array x       ***'
      allocate(y(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array y       ***'
      allocate(connect1(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect1***'
      allocate(connect2(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array connect2***'
      allocate(xc(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xc      ***'
      allocate(yc(nlon,nlat),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yc      ***'
 
c     Init the mapping arrays
      connectval1 = 0
      connectval2 = 0
      do i=1,nlon
	do j=1,nlat
           x(i,j)        = 0.
           y(i,j)        = 0.
           xc(i,j)       = 0.
           yc(i,j)       = 0.
           connect1(i,j) = 1
           connect2(i,j) = 1
        enddo
      enddo
 
c     Get the good radius for gridding - avoiding gaps
      radunit = 'deg'
      radius  = 0.
      do i=2,nx
        do j=2,ny
           if ( abs( lat(i,j)-lat(i,j-1) ).gt.radius ) then
              radius = abs( lat(i,j)-lat(i,j-1) )
           endif
           if ( abs( lon(i,j)-lon(i-1,j) ).gt.radius ) then
              radius = abs( lon(i,j)-lon(i-1,j) )
           endif
        enddo
      enddo
      radius = 3. * radius 
 
c     Do the mapping
      do i=1,nx
	do j=1,ny
 
          connectval1 = connectval1 + 1
          call gridding1(lat(i,j),lon(i,j),real(i),radius,radunit,
     >                   connect1,connectval1,
     >                   x,xc,nlon,nlat,lonmin,latmin,dlon,dlat)
 
          connectval2 = connectval2 + 1
          call gridding1(lat(i,j),lon(i,j),real(j),radius,radunit,
     >                   connect2,connectval2,
     >                   y,yc,nlon,nlat,lonmin,latmin,dlon,dlat)
 
 
	enddo
      enddo
 
c     Normalize output and set miising data flag
      do i=1,nlon
         do j=1,nlat
            if ( xc(i,j).gt.0. ) then
               x(i,j) = x(i,j)/xc(i,j)
            else
               x(i,j) = mdv
            endif
            if ( yc(i,j).gt.0. ) then
               y(i,j) = y(i,j)/yc(i,j)
            else
               y(i,j) = mdv
            endif
         enddo
      enddo
 
c     Write data to netCDF
      cdfname = mapfile
      varname = 'X'
      call writecdf2D(cdfname,varname,x,0.,
     >                dlon,dlat,lonmin,latmin,nlon,nlat,1,1,1)
 
      cdfname = mapfile
      varname = 'Y'
      call writecdf2D(cdfname,varname,y,0.,
     >                dlon,dlat,lonmin,latmin,nlon,nlat,1,0,1)
 
      cdfname = mapfile
      varname = 'LON'
      call writecdf2D(cdfname,varname,lon,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 
      cdfname = mapfile
      varname = 'LAT'
      call writecdf2D(cdfname,varname,lat,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 
      cdfname = mapfile
      varname = 'Z'
      call writecdf2D(cdfname,varname,z3,0.,
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
 
      cdfname = mapfile
      varname = 'P'
      call writecdf2D(cdfname,varname,p3,0.,
     >                dx,dy,xmin,ymin,nx,ny,nz,0,1)
      cdfname = mapfile
      varname = 'TOPO'
      call writecdf2D(cdfname,varname,zb,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 
 
      endif
 
 
c     ------------------------------------------------------------
c     Convert a lat/lon/z list to a x/y/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-ll2xy' ) then
 
c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 
c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1 ) outmode = inpmode
 
c     Read coordinates from file - Format (lon,lat,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)
 
c     Read coordinates from trajectory file 
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)
 
         if ( (vars(2).ne.'lon').and.(vars(3).ne.'lat') ) then
            print*,' ERROR: input must be in lon/lat format... Stop'       
            stop
         endif
 
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2) 
               yy0(npoints) = tra(i,j,3) 
               zz0(npoints) = tra(i,j,4) 
            enddo
         enddo
      endif 
 
c     Transform coordinates
      do i=1,npoints
 
          if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >         ( abs(yy0(i)-mdv).gt.eps )  )
     >    then
             rid=(xx0(i)-lonmin)/dlon+1.
             rjd=(yy0(i)-latmin)/dlat+1.
 
             xx0(i) = int_index3(x,nlon,nlat,1,rid,rjd,1.,mdv)
             yy0(i) = int_index3(y,nlon,nlat,1,rid,rjd,1.,mdv)
          else
             xx0(i) = mdv
             yy0(i) = mdv
          endif
 
      enddo
 
c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)
 
c     Write output to trajectory
      else         
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif
 
      endif
 
c     ------------------------------------------------------------
c     Convert a x/y/z list to a lon/lat/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-xy2ll' ) then
 
c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 
c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode
 
c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)
 
c     Read coordinates from trajectory file 
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)
 
         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'       
            stop
         endif
 
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2) 
               yy0(npoints) = tra(i,j,3) 
               zz0(npoints) = tra(i,j,4) 
            enddo
         enddo
      endif 
 
c     Transform coordinates
      do i=1,npoints
 
          if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >         ( abs(yy0(i)-mdv).gt.eps )  )
     >    then
              rid=(xx0(i)-xmin)/dx+1.
              rjd=(yy0(i)-ymin)/dy+1.
 
              xx0(i) = int_index3(lon,nx,ny,1,rid,rjd,1.,mdv)
              yy0(i) = int_index3(lat,nx,ny,1,rid,rjd,1.,mdv)
          else
              xx0(i) = mdv
              yy0(i) = mdv
          endif
 
      enddo
 
c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)
 
c     Write output to trajectory
      else         
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'lon'
         vars(3) = 'lat'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif
 
      endif
 
c     ------------------------------------------------------------
c     Convert a single lat/lon/z list to a x/y/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-ll2xy.single' ) then
 
c     Transform coordinates
      rid  = (lonpos-lonmin)/dlon+1.
      rjd  = (latpos-latmin)/dlat+1.
      xpos = int_index3(x,nlon,nlat,1,rid,rjd,1.,mdv)
      ypos = int_index3(y,nlon,nlat,1,rid,rjd,1.,mdv)
 
c	  Write result to screen
      print*,xpos,ypos
 
      endif
 
c     ------------------------------------------------------------
c     Convert a single x/y/z list to a lon/lat/z list
c     ------------------------------------------------------------
 
      if ( mode.eq.'-xy2ll.single' ) then
 
c     Transform coordinates
      rid    = (xpos-xmin)/dx+1.
      rjd    = (ypos-ymin)/dy+1.
      lonpos = int_index3(lon,nx,ny,1,rid,rjd,1.,mdv)
      latpos = int_index3(lat,nx,ny,1,rid,rjd,1.,mdv)
 
c	  Write result to screen
      print*,lonpos,latpos
 
      endif
 
c     ------------------------------------------------------------
c     Calculate numerically the maps scale factor
c     ------------------------------------------------------------
 
      if ( mode.eq.'-mapscale' ) then
 
c	  Allocate memory for map scale factors
      allocate(mpx(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array mpx       ***'
      allocate(mpy(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array mpy       ***'
 
c     Get map scale for interior points
      do i=2,nx-1
         do j=2,ny-1
 
c           Map scale in x direction
            lon1     = lon(i-1,j)
            lat1     = lat(i-1,j)
            lon2     = lon(i+1,j)
            lat2     = lat(i+1,j)
            radunit  = 'km.haversine'
            mpx(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)
 
c           Map scale in y direction
            lon1     = lon(i,j-1)
            lat1     = lat(i,j-1)
            lon2     = lon(i,j+1)
            lat2     = lat(i,j+1)
            radunit  = 'km.haversine'
            mpy(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)
 
         enddo
      enddo
 
c     Copy map scale for boundary line points
      do i=2,nx-1
        mpx(i, 1) = mpx(i,   2)
        mpx(i,ny) = mpx(i,ny-1)
        mpy(i, 1) = mpy(i,   2)
        mpy(i,ny) = mpy(i,ny-1)
      enddo
      do j=2,ny-1
        mpx(1, j) = mpx(2,   j)
        mpx(nx,j) = mpx(nx-1,j)
        mpy(1, j) = mpy(2,   j)
        mpy(nx,j) = mpy(nx-1,j)
      enddo
 
c     Copy map scale factor for boundary corner points
      mpx( 1, 1) = mpx(   2,   2)
      mpy( 1, 1) = mpy(   2,   2)
      mpx(nx, 1) = mpx(nx-1,   2)
      mpy(nx, 1) = mpy(nx-1,   2)
      mpx(nx,ny) = mpx(nx-1,ny-1)
      mpy(nx,ny) = mpy(nx-1,ny-1)
      mpx( 1,ny) = mpx(   2,ny-1)
      mpy( 1,ny) = mpy(   2,ny-1)
 
c     Apply a diffuive filter
      do i=1,nfilter
         call filt2d (mpx,mpx,nx,ny,1.0,-999.,0,0,0,0)
         call filt2d (mpy,mpy,nx,ny,1.0,-999.,0,0,0,0)
      enddo
 
c     Write map factors to mapping file
      cdfname = mapfile
      varname = 'MAPSCALE_X'
      call writecdf2D(cdfname,varname,mpx,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
      cdfname = mapfile
      varname = 'MAPSCALE_Y'
      call writecdf2D(cdfname,varname,mpy,0.,
     >                dx,dy,xmin,ymin,nx,ny,1,0,1)
 
 
      endif
 
c     ------------------------------------------------------------
c     Convert from pressure to height
c     ------------------------------------------------------------
 
      if ( mode.eq.'-p2z' ) then
 
c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 
c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode
 
c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)
 
c     Read coordinates from trajectory file
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)
 
         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'
            stop
         endif
 
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2)
               yy0(npoints) = tra(i,j,3)
               zz0(npoints) = tra(i,j,4)
            enddo
         enddo
      endif
 
c     Transform coordinates
      do i=1,npoints
 
         if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >        ( abs(yy0(i)-mdv).gt.eps ).and.
     >        ( abs(zz0(i)-mdv).gt.eps ) )
     >    then
             call get_index3 (rid,rjd,rkd,xx0(i),yy0(i),-zz0(i),1,
     >                        -p3,-ps,nx,ny,nz,xmin,ymin,dx,dy)
 
             zz0(i) = int_index3(z3,nx,ny,nz,rid,rjd,rkd,mdv)
          else
             xx0(i) = mdv
             yy0(i) = mdv
             zz0(i) = mdv
          endif
 
      enddo
 
c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)
 
c     Write output to trajectory
      else
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'z'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif
 
      endif
 
c     ------------------------------------------------------------
c     Convert single x/y/p to x/y/z
c     ------------------------------------------------------------
 
      if ( mode.eq.'-p2z.single' ) then
 
c     Transform coordinates - change from descending to ascending (minus sign)
      call get_index3 (rid,rjd,rkd,xpos,ypos,-ppos,1,
     >                 -p3,-ps,nx,ny,nz,xmin,ymin,dx,dy)
 
      zpos = int_index3(z3,nx,ny,nz,rid,rjd,rkd,mdv)
 
c	  Write result to screen
      print*,xpos,ypos,zpos
 
      endif
 
c     ------------------------------------------------------------
c     Convert from height to pressure
c     ------------------------------------------------------------
 
      if ( mode.eq.'-z2p' ) then
 
c     Allocate memory for input and output lists
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra         ***'
      allocate(xx0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0         ***'
      allocate(yy0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0         ***'
      allocate(zz0(ntra*ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zz0         ***'
 
c     Get format of input and output file
      call mode_tra(inpmode,inpfile)
      call mode_tra(outmode,outfile)
      if ( outmode.eq.-1) outmode=inpmode
 
c     Read coordinates from file - Format (x,y,lev)
      if (inpmode.eq.-1) then
         fid = 10
         npoints = 0
         open(fid,file=inpfile)
          do i=1,ntra
             npoints = npoints + 1
             read(fid,*) xx0(npoints),yy0(npoints),zz0(npoints)
          enddo
         close(fid)
 
c     Read coordinates from trajectory file
      else
         call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
         call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
         call close_tra(fid,inpmode)
 
         if ( (vars(2).ne.'x').and.(vars(3).ne.'y') ) then
            print*,' ERROR: input must be in x/y format... Stop'
            stop
         endif
 
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               xx0(npoints) = tra(i,j,2)
               yy0(npoints) = tra(i,j,3)
               zz0(npoints) = tra(i,j,4)
            enddo
         enddo
      endif
 
c     Transform coordinates
      do i=1,npoints
 
         if ( ( abs(xx0(i)-mdv).gt.eps ).and.
     >        ( abs(yy0(i)-mdv).gt.eps ).and.
     >        ( abs(zz0(i)-mdv).gt.eps ) )
     >    then
               call get_index3 (rid,rjd,rkd,xx0(i),yy0(i),zz0(i),1,
     >                          z3,zb,nx,ny,nz,xmin,ymin,dx,dy)
 
               zz0(i) = int_index3(p3,nx,ny,nz,rid,rjd,rkd,mdv)
          else
               xx0(i) = mdv
               yy0(i) = mdv
               zz0(i) = mdv
          endif
 
      enddo
 
c     Write output to list
      if ( outmode.eq.-1 ) then
         fid = 10
         open(fid,file=outfile)
         do i=1,ntra
            write(fid,'(f9.2,f8.2,i6)') xx0(i),yy0(i),nint(zz0(i))
         enddo
         close(fid)
 
c     Write output to trajectory
      else
         npoints = 0
         do i=1,ntra
            do j=1,ntim
               npoints = npoints + 1
               tra(i,j,2) = xx0(npoints)
               tra(i,j,3) = yy0(npoints)
               tra(i,j,4) = zz0(npoints)
            enddo
         enddo
         if ( ncol.eq.3 ) ncol = 4
         vars(1) = 'time'
         vars(2) = 'x'
         vars(3) = 'y'
         vars(4) = 'p'
         call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
         call write_tra(fid,tra,ntra,ntim,ncol,outmode)
         call close_tra(fid,outmode)
      endif
 
      endif
 
c     ------------------------------------------------------------
c     Convert single x/y/z to x/y/p
c     ------------------------------------------------------------
 
      if ( mode.eq.'-z2p.single' ) then
 
c     Transform coordinates - change from descending to ascending (minus sign)
      call get_index3 (rid,rjd,rkd,xpos,ypos,zpos,1,
     >                 z3,zb,nx,ny,nz,xmin,ymin,dx,dy)
 
      ppos = int_index3(p3,nx,ny,nz,rid,rjd,rkd,mdv)
 
c	  Write result to screen
      print*,xpos,ypos,ppos
 
      endif
 
      end
 
 
c     ********************************************************************
c     * GRIDDING SUBROUTINES                                             *
c     ********************************************************************
 
c     ---------------------------------------------------------------------
c     Gridding of one single data point (smoothing in km, deg, gridp)
c     ---------------------------------------------------------------------
 
      subroutine gridding1 (lat,lon,addval,radius,unit,
     >                      connect,connectval,
     >                      out,outc,nx,ny,xmin,ymin,dx,dy)
 
      implicit none
 
c     Declaration of subroutine parameters
      real         lat,lon
      integer      nx,ny
      real         xmin,ymin,dx,dy
      real         out(nx,ny),outc(nx,ny)
      real         radius
      character*80 unit
      integer      connectval
      integer      connect(nx,ny)
      real         addval
 
c     Auxiliary variables
      integer   i,j,k
      integer   mu,md,nr,nl,n,m
      integer   stackx(nx*ny),stacky(nx*ny)
      integer   tab_x(nx*ny),tab_y(nx*ny)
      real      tab_r(nx*ny)
      integer   sp
      real      lat2,lon2
      real      dist,sum
      real      xmax
      integer   periodic
      integer   test
      character ch
 
c     Numerical epsilon
      real      eps
      parameter (eps=0.01)
 
c     Externals
      real      sdis,weight
      external  sdis,weight
 
c     Check whether lat/lon point is valid
      xmax=xmin+real(nx-1)*dx
      if (lon.lt.xmin-eps) lon=lon+360.
      if (lon.gt.xmax+eps) lon=lon-360.
      if (abs(lat-90).lt.eps) lat=90.
      if (abs(lat+90).lt.eps) lat=-90.
      if ((abs(lat).gt.(90.+eps)).or.
     >    (lon.lt.xmin-eps).or.(lon.gt.xmax+eps)) then
         print*,'Invalid lat/lon point ',lat,lon
         return
      endif
 
c     Set flag for periodic domain
      if (abs(xmax-xmin-360.).lt.eps) then
         periodic=1
      else if (abs(xmax-xmin-360+dx).lt.eps) then
         periodic=2
      else
         periodic=0
      endif
 
c     Get indices of one coarse grid point within search radius
      i=nint((lon-xmin)/dx)+1
      if ((i.eq.nx).and.(periodic.eq.1)) i=1
      j=nint((lat-ymin)/dy)+1
      lat2=ymin+real(j-1)*dy
      lon2=xmin+real(i-1)*dx
      dist=sdis(lon,lat,lon2,lat2,unit)
      if (dist.gt.radius) then
         print*,'1: Search radius is too small...'
         stop
      endif
 
c     Get connected points
      k=0
      stackx(1)=i
      stacky(1)=j
      sp    = 1
      do while (sp.ne.0) 
         
c        Get an element from stack
         n=stackx(sp)
         m=stacky(sp)
         sp=sp-1
                  
c        Get distance from reference point
         lat2=ymin+real(m-1)*dy
         lon2=xmin+real(n-1)*dx
         dist=sdis(lon,lat,lon2,lat2,unit)
 
c        Check whether distance is smaller than search radius: connected
         if (dist.lt.radius) then
 
c           Make entry in filter mask
            k=k+1
            tab_x(k)=n
            tab_y(k)=m
            tab_r(k)=weight(dist,radius)
 
c           Mark this point as visited
            connect(n,m)=connectval
                     
c           Get coordinates of neighbouring points
            nr=n+1
            if ((nr.gt.nx)  .and.(periodic.eq.0)) nr=nx
            if ((nr.gt.nx-1).and.(periodic.eq.1)) nr=1
            if ((nr.gt.nx)  .and.(periodic.eq.2)) nr=1
            nl=n-1
            if ((nl.lt.1).and.(periodic.eq.0)) nl=1
            if ((nl.lt.1).and.(periodic.eq.1)) nl=nx-1
            if ((nl.lt.1).and.(periodic.eq.2)) nl=nx
            mu=m+1
            if (mu.gt.ny) mu=ny
            md=m-1
            if (md.lt.1) md=1
 
c           Update stack
            if (connect(nr,m).ne.connectval) then
               connect(nr,m)=connectval
               sp=sp+1
               stackx(sp)=nr
               stacky(sp)=m
            endif
            if (connect(nl,m).ne.connectval) then
               connect(nl,m)=connectval
               sp=sp+1
               stackx(sp)=nl
               stacky(sp)=m
            endif
            if (connect(n,mu).ne.connectval) then
               connect(n,mu)=connectval
               sp=sp+1
               stackx(sp)=n
               stacky(sp)=mu
            endif
            if (connect(n,md).ne.connectval) then
               connect(n,md)=connectval
               sp=sp+1
               stackx(sp)=n
               stacky(sp)=md
            endif
         endif     
 
      end do
 
      if (k.ge.1) then
         sum=0.
         do i=1,k
            sum=sum+tab_r(i)
         enddo
         do i=1,k
            out (tab_x(i),tab_y(i))=out(tab_x(i),tab_y(i))+
     >                             addval*tab_r(i)/sum
            outc(tab_x(i),tab_y(i))=outc(tab_x(i),tab_y(i))+
     >                             1.*tab_r(i)/sum
 
            if ((tab_x(i).eq.1).and.(periodic.eq.1)) then
               out(nx,tab_y(i))=out(nx,tab_y(i))+
     >                             addval*tab_r(i)/sum
               outc(nx,tab_y(i))=outc(nx,tab_y(i))+
     >                             1.*tab_r(i)/sum
 
            endif
         enddo
      else
         print*,'2: Search radius is too small...'
         stop
      endif
 
      end
 
 
c     ----------------------------------------------------------------------
c     Get spherical distance between lat/lon points
c     ----------------------------------------------------------------------
        
      real function sdis(xp,yp,xq,yq,unit)
 
c     Calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
 
      real,parameter ::       re=6371.229
      real,parameter ::       rad2deg=180./3.14159265
      real,parameter ::       deg2rad=3.141592654/180.
 
      real         xp,yp,xq,yq,arg
      character*80 unit
      real         dlon,dlat,alpha,rlat1,ralt2
 
      if ( unit.eq.'km' ) then
 
         arg=sin(yp*deg2rad)*sin(yq*deg2rad)+
     >           cos(yp*deg2rad)*cos(yq*deg2rad)*cos((xp-xq)*deg2rad)
         if (arg.lt.-1.) arg=-1.
         if (arg.gt.1.) arg=1.
         sdis=re*acos(arg)
         
      elseif ( unit.eq.'deg' ) then
 
         dlon = xp-xq
         if ( dlon.gt. 180. ) dlon = dlon - 360.
         if ( dlon.lt.-180. ) dlon = dlon + 360.
         sdis = sqrt( dlon**2 + (yp-yq)**2 )
 
      elseif ( unit.eq.'km.haversine' ) then
 
        dlat   =  (yp - yq)*deg2rad
        dlon   =  (xp - xq)*deg2rad
        rlat1   =  yp*deg2rad
        rlat2   =  yq*deg2rad
 
        alpha  = ( sin(0.5*dlat)**2 ) +
     >           ( sin(0.5*dlon)**2 )*cos(rlat1)*cos(rlat2)
 
        sdis = 4. * re * atan2( sqrt(alpha),1. + sqrt( 1. - alpha ) )
 
      endif  
 
      end
 
c     ----------------------------------------------------------------------
c     Weight function for the filter mask
c     ----------------------------------------------------------------------
 
      real function weight (r,radius)
 
c     Attribute to each distanc r its corresponding weight in the filter mask
 
      implicit none
 
c     Declaration of subroutine parameters
      real r
      real radius
 
c     Simple 0/1 mask
      if (r.lt.radius) then
          weight=1.-0.65*r/radius 
      else
          weight=0.
      endif
 
      end
 
 
c     ********************************************************************
c     * INPUT / OUTPUT SUBROUTINES                                       *
c     ********************************************************************
 
c     --------------------------------------------------------------------
c     Subroutines to write the netcdf output file
c     --------------------------------------------------------------------
 
      subroutine writecdf2D(cdfname,
     >                      varname,arr,time,
     >                      dx,dy,xmin,ymin,nx,ny,nz,
     >                      crefile,crevar)
 
      IMPLICIT NONE
 
c     Declaration of input parameters
      character*80 cdfname,varname
      integer nx,ny,nz
      real arr(nx,ny,nz)
      real dx,dy,xmin,ymin
      real time
      integer crefile,crevar
 
c     Further variables
      real varmin(4),varmax(4),stag(4)
      integer ierr,cdfid,ndim,vardim(4)
      character*80 cstname
      real mdv
      integer i
      integer nvars
      character*80 vars(300)
 
C     Definitions
      varmin(1)=xmin
      varmin(2)=ymin
      varmin(3)=1050.
      varmax(1)=xmin+real(nx-1)*dx
      varmax(2)=ymin+real(ny-1)*dy
      varmax(3)=1050.
      cstname='no_constants_file'
      ndim=4
      vardim(1)=nx
      vardim(2)=ny
      vardim(3)=nz
      stag(1)=0.
      stag(2)=0.
      stag(3)=0.
      mdv=-999.98999
 
C     Create the file
      if (crefile.eq.0) then
         call cdfwopn(cdfname,cdfid,ierr)
         if (ierr.ne.0) goto 906
      else if (crefile.eq.1) then
         call crecdf(cdfname,cdfid,
     >        varmin,varmax,ndim,cstname,ierr)
         if (ierr.ne.0) goto 903
      endif
 
c     Check whether the variable is already on the file
      call getvars(cdfid,nvars,vars,ierr)
      if (ierr.ne.0) goto 906
      do i=1,nvars
         if ( vars(i).eq.varname ) crevar = 0
      enddo
 
c     Write the data to the netcdf file, and close the file
      if (crevar.eq.1) then
         call putdef(cdfid,varname,ndim,mdv,
     >        vardim,varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 904
      endif
      call putdat(cdfid,varname,time,0,arr,ierr)
      if (ierr.ne.0) goto 905
      call clscdf(cdfid,ierr)
 
      return
 
c     Error handling
 903  print*,'*** Problem to create netcdf file ***'
      stop
 904  print*,'*** Problem to write definition ***'
      stop
 905  print*,'*** Problem to write data ***'
      stop
 906  print*,'*** Problem to open netcdf file ***'
      stop
 
      end
 
 
c     -----------------------------------------------
c     Subroutine to read a netcdf output file
c     -----------------------------------------------
 
      subroutine readcdf2D(cdfname,varname,arr,time,
     >                     dx,dy,xmin,ymin,nx,ny,nz)
 
      implicit none
 
c     Declaration of input parameters
      character*80 cdfname,varname
      integer      nx,ny,nz
      real         arr(nx,ny,nz)
      real         dx,dy,xmin,ymin
      real         time
 
c     Internal variables for the netcdf file
      real         varmin(4),varmax(4),stag(4)
      integer      ierr,cdfid,ndim,vardim(4)
      real         mdv
      real         delx,dely
      real         times(500)
      integer      ntimes,flag
 
c     Numerical epsilon
      real         eps
      parameter    (eps=0.01)
 
c     Auxiliary variables
      integer      i,j,k
 
c     Initialize the ouput array
      do i=1,nx
         do j=1,ny
            do k=1,nz
               arr(i,j,k)=0.
            enddo
         enddo
      enddo
 
c     Try to open the netcdf file and to read the variable
      call cdfopn(cdfname,cdfid,ierr)
      if (ierr.ne.0) goto 801
 
c     Check whether the specifieed time is available
      call gettimes(cdfid,times,ntimes,ierr)
      if (ierr.ne.0) goto 801
      flag=0
      do i=1,ntimes
          if (abs(times(i)-time).lt.eps) flag=1
      enddo
      if (flag.eq.0) then
         print*,'No such time on cdf',time
         stop
      endif
 
c     Get the variable definition, and the grid
      call getdef(cdfid,varname,ndim,mdv,vardim,varmin,varmax,
     >     stag,ierr)
      if (ierr.ne.0) goto 801
 
c     Set parameters or read parameters - depending on nx,ny
      if ( (nx.eq.1).and.(ny.eq.1) ) then
         dx   = (varmax(1)-varmin(1))/real(vardim(1)-1)
         dy   = (varmax(2)-varmin(2))/real(vardim(2)-1)
         nx   = vardim(1)
         ny   = vardim(2)
         nz   = vardim(3)
         xmin = varmin(1)
         ymin = varmin(2)
      else
         call getdat(cdfid,varname,time,0,arr,ierr)
         if (ierr.ne.0) goto 801
      endif
 
c     Close data file
      call clscdf(cdfid,ierr)
 
      return
      
c     Exception handling
 801  print*,'Problems in reading netcdf file'
      stop
      
      end
 
 
c     -----------------------------------------------
c     Diffusive filter
c     -----------------------------------------------
 
      subroutine filt2d (a,af,nx,ny,fil,misdat,
     &                   iperx,ipery,ispol,inpol)
 
c     Apply a conservative diffusion operator onto the 2d field a,
c     with full missing data checking.
c
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
c     af    real   out  filtered array, dimensioned (nx,ny), can be
c                       equivalenced with array a in the calling routine
c     f1    real        workarray, dimensioned (nx+1,ny)
c     f2    real        workarray, dimensioned (nx,ny+1)
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
c                       corresponds to one application of the linear filter.
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
c                       the corresponding value is not available. The
c                       misdat-checking can be switched off with with misdat=0.
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
c
c     Christoph Schaer, 1993
 
c     argument declaration
      integer     nx,ny
      real        a(nx,ny),af(nx,ny),fil,misdat
      integer     iperx,ipery,inpol,ispol
 
c     local variable declaration
      integer     i,j,is
      real        fh
      real        f1(nx+1,ny),f2(nx,ny+1)
 
c     compute constant fh
      fh=0.125*fil
 
c     compute fluxes in x-direction
      if (misdat.eq.0.) then
        do j=1,ny
        do i=2,nx
          f1(i,j)=a(i-1,j)-a(i,j)
        enddo
        enddo
      else
        do j=1,ny
        do i=2,nx
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
            f1(i,j)=0.
          else
            f1(i,j)=a(i-1,j)-a(i,j)
          endif
        enddo
        enddo
      endif
      if (iperx.eq.1) then
c       do periodic boundaries in the x-direction
        do j=1,ny
          f1(1,j)=f1(nx,j)
          f1(nx+1,j)=f1(2,j)
        enddo
      else
c       set boundary-fluxes to zero
        do j=1,ny
          f1(1,j)=0.
          f1(nx+1,j)=0.
        enddo
      endif
 
c     compute fluxes in y-direction
      if (misdat.eq.0.) then
        do j=2,ny
        do i=1,nx
          f2(i,j)=a(i,j-1)-a(i,j)
        enddo
        enddo
      else
        do j=2,ny
        do i=1,nx
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
            f2(i,j)=0.
          else
            f2(i,j)=a(i,j-1)-a(i,j)
          endif
        enddo
        enddo
      endif
c     set boundary-fluxes to zero
      do i=1,nx
        f2(i,1)=0.
        f2(i,ny+1)=0.
      enddo
      if (ipery.eq.1) then
c       do periodic boundaries in the x-direction
        do i=1,nx
          f2(i,1)=f2(i,ny)
          f2(i,ny+1)=f2(i,2)
        enddo
      endif
      if (iperx.eq.1) then
        if (ispol.eq.1) then
c         do south-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
          enddo
        endif
        if (inpol.eq.1) then
c         do north-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
          enddo
        endif
      endif
 
c     compute flux-convergence -> filter
      if (misdat.eq.0.) then
        do j=1,ny
        do i=1,nx
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
        enddo
        enddo
      else
        do j=1,ny
        do i=1,nx
          if (a(i,j).eq.misdat) then
            af(i,j)=misdat
          else
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
          endif
        enddo
        enddo
      endif
      end
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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

Subversion Repositories lagranto.wrf

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