Subversion Repositories pvinversion.ecmwf

      PROGRAM modify_anomaly

c     ******************************************************************************
c     * Read the modified and unmodified PV from the R file and                  *
c     * perform some non-standard modifications: (a) Change of amplitude,          *
c     * (b) Stretching and shrinking along an axis, (c) rotation, or               *
c     * (d) change in position.     
c     *                                                                            
c     # modify_anomaly.sh R_20060116_18 shifty=-10,rot=10,stretchy=1.5,rot=-10,shifty=10
c     # modify_anomaly.sh R_20060116_18 fac=2
c     # modify_anomaly.sh R_20060116_18 cex=-30,cey=-30,rot=20
c     *                                                                            *
c     * Michael Sprenger / Spring, Summer 2007                                     *
c     ******************************************************************************

      implicit none

c     -----------------------------------------------------------------------------
c     Declaration of parameters and variables
c     -----------------------------------------------------------------------------

c     Input/output file and command string
      character*80   pvsrcfile
      character*80   commandstr

c     Grid parameters
      integer        nx,ny,nz
      real           xmin,ymin,zmin,xmax,ymax,zmax
      real           dx,dy,dz
      real           mdv

c     3d fields for calculation of qgPV and Ertel's PV
      real,allocatable,dimension (:,:,:) :: pv1,pv2,ano

c     Numerical epsilon
      real           eps
      parameter      (eps=0.01)

c     Parameters for the transformations
      real           cex,cey              ! Centre point for rotation
      real           angle                ! Rotation angle

c     Auxiliary variables
      real           zpos
      integer        i,j,k
      integer        stat
      character*80   varname
      integer        n
      real           par
      character*80   com

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

      print*,'********************************************************'
      print*,'* MODIFY_ANOMALY                                       *'
      print*,'********************************************************'

c     Read parameter file
      open(10,file='fort.10')
       read(10,*) pvsrcfile
       read(10,*) commandstr
      close(10)
      print*
      print*,'Input file : ',trim(pvsrcfile)
      print*,'Command    : ',trim(commandstr)
      print*

c     Get lat/lon gid parameters from input file
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >               pvsrcfile)
      print*,'Read_Dim: nx,ny,nz         = ',nx,ny,nz
      print*,'          dx,dy,dz         = ',dx,dy,dz
      print*,'          xmin,ymin,zmin   = ',xmin,ymin,zmin
      print*,'          mdv              = ',mdv
      print*

c     Count from 0, not from 1: consistent with <inv_cart.f>.
      nx=nx-1
      ny=ny-1
      nz=nz-1

c     Allocate memory for modified and non-modified PV
      allocate(pv1 (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating pv1'
      allocate(pv2 (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating pv2'
      allocate(ano  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating ano'   
   
c     Read data from file
      varname='PV'
      call read_inp (pv1,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      varname='PV_FILT'
      call read_inp (pv2,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)

c     Define the anomaly
      do i=0,nx
         do j=0,ny
            do k=0,nz
               if ( (abs(pv1(i,j,k)-mdv).gt.eps).and.
     >              (abs(pv2(i,j,k)-mdv).gt.eps) ) then
                  ano(i,j,k)=pv1(i,j,k)-pv2(i,j,k)
               else
                  ano(i,j,k)=0.
               endif
            enddo
         enddo
      enddo


c     -------------------------------------------------------------------------------
c     Modifications
c     -------------------------------------------------------------------------------

c     Set the default values for parameters
      cex   = 0.
      cey   = 0.
      angle = 0.

c     Set the counter for the command string
      n=1

c     Loop over all commands
 100  continue

c     Extract new command/parameter pair; exit if no new command
      call next_command(commandstr,n,com,par)
      if (com.eq.'nil') goto 200
      print*,trim(com),par

c     Multiply the anomaly by a constant factor
      if (com.eq.'fac') then
         call mod_factor (ano,par,nx,ny,nz,mdv)
      endif

c     Set the centre point (needed for rotations)
      if (com.eq.'cex') then
         cex=par
      endif
      if (com.eq.'cey') then
         cey=par
      endif

c     Rotation around the centrre point
      if (com.eq.'rot') then
         call mod_rotation (ano,par,cex,cey,
     >                       nx,ny,nz,xmin,ymin,dx,dy,mdv)
      endif

c     Shift in x or in y direction
      if (com.eq.'shiftx') then
         call  mod_shift (ano,par,0.,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
      endif
      if (com.eq.'shifty') then
         call  mod_shift (ano,0.,par,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
      endif      

c     Stretch/shrink in x or in y direction
      if (com.eq.'stretchx') then
         call  mod_stretch (ano,par,1.,
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
      endif
      if (com.eq.'stretchy') then
         call  mod_stretch (ano,1.,par,
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
      endif      


c     Goto next command/parameter pair
      goto 100

c     All commands handled
 200  continue

c     --------------------------------------------------------------------------------
c     Write modified fields
c     --------------------------------------------------------------------------------

c     Build the modified PV distribution from the anomaly
      do i=0,nx
         do j=0,ny
            do k=0,nz
               if ( (abs(pv1(i,j,k)-mdv).gt.eps).and.
     >              (abs(ano(i,j,k)-mdv).gt.eps) ) then
                  pv2(i,j,k)=pv1(i,j,k)-ano(i,j,k)
               else
                  pv2(i,j,k)=mdv
               endif
            enddo
         enddo
      enddo

c     Write result to netcdf file
      varname='PV_FILT'
      call write_inp (pv2,varname,pvsrcfile,nx,ny,nz)

c     Write the modified anomaly also to the file
      varname='PV_ANOM'
      call write_inp (ano,varname,pvsrcfile,nx,ny,nz)

      end

c     ********************************************************************************
c     * Command handling and transformations                                         *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Extract next command from command string
c     --------------------------------------------------------------------------------
      
      subroutine next_command (commandstr,n,com,par)
      
c     Given the command string <commandstr>, extract the next command/parameter pair
c     <com,par> starting at position <n> of the command string.

      implicit none
      
c     Declaration of subroutine parameters
      character*80 commandstr
      character*80 com
      real         par
      integer      n

c     Auxiliary variables
      integer      i,j,k

c     Check whether end of command line reached
      if (n.ge.80) then
         com='nil'
         par=0.
         goto 120
      endif

c     Set indices to next next command/parameter pair
      i=n
      j=n
 100  if ( (commandstr(j:j).ne.'=').and.(j.lt.80) ) then
        j=j+1
        goto 100
      endif
      k=j+1
 110  if ( (commandstr(k:k).ne.',' ).and.
     >     (commandstr(k:k).ne.';' ).and.
     >     (commandstr(k:k).ne.' ' ).and.
     >     (k              .lt.80  ) ) then
        k=k+1
        goto 110
      endif

c     Check whether this is a valid command/parameter pair
      if ( ((j-1).lt.i).or.((k-1).lt.(j+1)) ) then
         com='nil'
         par=0.
         goto 120
      endif

c     Extract tzhe command and the parameter
      com=commandstr(i:j-1)
      read(commandstr(j+1:k-1),*) par

c     Set the counter to the next command/parameter pair
      n=k+1

c     Exit point
 120  continue

      end

c     --------------------------------------------------------------------------------
c     Multiply the anomaly by a factor
c     --------------------------------------------------------------------------------
 
      subroutine mod_factor (field,factor,nx,ny,nz,mdv)
      
c     Multiply the anomaly <field>  by a constant factor <factor>. The grid and the
c     missing data value are given by <nx,ny,nz,mdv>.
      
      implicit none
 
c     Declaration of subroutine parameters
      integer    nx,ny,nz
      real       field(0:nx,0:ny,0:nz)
      real       mdv
      real       factor
      
c     Parameters
      real       eps
      parameter  (eps=0.01)

c     Auxiliary variables
      integer    i,j,k

c     Do the transformation
      do i=0,nx
         do j=0,ny
            do k=0,nz
               if ( (abs(field(i,j,k)-mdv).gt.eps) ) then 
                  field(i,j,k)=factor*field(i,j,k)
               endif
            enddo
         enddo
      enddo
      
      end

c     --------------------------------------------------------------------------------
c     Rotate the anomaly
c     --------------------------------------------------------------------------------

      subroutine mod_rotation (field,angle,cex,cey,
     >                         nx,ny,nz,xmin,ymin,dx,dy,mdv)

c     Rotate the anomaly <field> in the horizontal by the angle <angle>. The centre
c     for the rotation is given by <cex,cey>, expressed in rotated longitude, latitude.
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
      
      
      implicit none
 
c     Declaration of subroutine parameters
      integer    nx,ny,nz
      real       field(0:nx,0:ny,0:nz)
      real       xmin,ymin,dx,dy
      real       mdv
      real       angle
      real       cex,cey
      
c     Parameters
      real       eps
      parameter  (eps=0.01)
      real       pi180
      parameter  (pi180=3.14159/180.)

c     Auxiliary variables
      integer    i,j,k
      real       lon,lat,rlon,rlat
      real       xmax,ymax
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)

c     Externals 
      real       int2d
      external   int2d

c     Maximal grid extension
      xmax=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy

c     Do the Transformation
      do k=0,nz

c        Copy level to 2d array
         do i=0,nx
            do j=0,ny
               tmp1(i,j)=field(i,j,k)
            enddo
         enddo

c        Rotate each grid point
         do i=0,nx
            do j=0,ny
               
c              Lon/lat coordinates of grid point
               lon=xmin+real(i)*dx-cex
               lat=ymin+real(j)*dy-cey

c              Rotation
               rlon = cex + lon*cos(angle*pi180) + lat*sin(angle*pi180)
               rlat = cey - lon*sin(angle*pi180) + lat*cos(angle*pi180)

c              Do the interpolation
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
     >         then
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
               else
                  tmp2(i,j)=0.
               endif
               
            enddo
         enddo
         
c        Copy 2d array to level
         do i=0,nx
            do j=0,ny
               field(i,j,k)=tmp2(i,j)
            enddo
         enddo
         
      enddo

      end

c     --------------------------------------------------------------------------------
c     Shift the anomaly in x and y direction
c     --------------------------------------------------------------------------------

      subroutine mod_shift (field,shiftx,shifty,
     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)

c     Shift the anomaly <field> in the horizontal by the vector <shiftx,shifty>. 
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
      
      
      implicit none
 
c     Declaration of subroutine parameters
      integer    nx,ny,nz
      real       field(0:nx,0:ny,0:nz)
      real       xmin,ymin,dx,dy
      real       mdv
      real       shiftx,shifty
      
c     Parameters
      real       eps
      parameter  (eps=0.01)

c     Auxiliary variables
      integer    i,j,k
      real       lon,lat,rlon,rlat
      real       xmax,ymax
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)

c     Externals 
      real       int2d
      external   int2d

c     Maximal grid extension
      xmax=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy

c     Do the Transformation
      do k=0,nz

c        Copy level to 2d array
         do i=0,nx
            do j=0,ny
               tmp1(i,j)=field(i,j,k)
            enddo
         enddo

c        Rotate each grid point
         do i=0,nx
            do j=0,ny
               
c              Lon/lat coordinates of grid point
               lon=xmin+real(i)*dx
               lat=ymin+real(j)*dy

c              shifted coordinates
               rlon = lon - shiftx
               rlat = lat - shifty

c              Do the interpolation
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
     >         then
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
               else
                  tmp2(i,j)=0.
               endif
               
            enddo
         enddo
         
c        Copy 2d array to level
         do i=0,nx
            do j=0,ny
               field(i,j,k)=tmp2(i,j)
            enddo
         enddo
         
      enddo

      end

c     --------------------------------------------------------------------------------
c     Stretch/shrink the anomaly in x and y direction
c     --------------------------------------------------------------------------------

      subroutine mod_stretch (field,stretchx,stretchy,
     >                        nx,ny,nz,xmin,ymin,dx,dy,mdv)

c     Stretch the anomaly <field> in the horizontal by the fatcors <stretchx,stretchy>. 
c     The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>.
      
      
      implicit none
 
c     Declaration of subroutine parameters
      integer    nx,ny,nz
      real       field(0:nx,0:ny,0:nz)
      real       xmin,ymin,dx,dy
      real       mdv
      real       stretchx,stretchy
      
c     Parameters
      real       eps
      parameter  (eps=0.01)

c     Auxiliary variables
      integer    i,j,k
      real       lon,lat,rlon,rlat
      real       xmax,ymax
      real       tmp1(0:nx,0:ny),tmp2(0:nx,0:ny)

c     Externals 
      real       int2d
      external   int2d

c     Maximal grid extension
      xmax=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy

c     Do the Transformation
      do k=0,nz

c        Copy level to 2d array
         do i=0,nx
            do j=0,ny
               tmp1(i,j)=field(i,j,k)
            enddo
         enddo

c        Rotate each grid point
         do i=0,nx
            do j=0,ny
               
c              Lon/lat coordinates of grid point
               lon=xmin+real(i)*dx
               lat=ymin+real(j)*dy

c              shifted coordinates
               rlon = 1./stretchx * lon
               rlat = 1./stretchy * lat

c              Do the interpolation
               if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and.
     >              (rlat.gt.ymin).and.(rlat.lt.ymax) ) 
     >         then
                  tmp2(i,j)=int2d(tmp1,rlon,rlat,
     >                            dx,dy,xmin,ymin,nx,ny,mdv)
               else
                  tmp2(i,j)=0.
               endif
               
            enddo
         enddo
         
c        Copy 2d array to level
         do i=0,nx
            do j=0,ny
               field(i,j,k)=tmp2(i,j)
            enddo
         enddo
         
      enddo

      end

c     ------------------------------------------------------------------------
c     Two-dimensional interpolation         
c     ------------------------------------------------------------------------

      real function int2d(ar,lon,lat,
     >                    dx,dy,xmin,ymin,nx,ny,mdv)

c     Interpolate the field <ar(nx,ny)> to the position <lon,lat>. The
c     grid is specified by <dx,dy,xmin,ymin,nx,ny>.

      implicit none

c     Declaration of subroutine paramters
      integer   nx,ny
      real      ar(0:nx,0:ny)
      real      lon,lat
      real      dx,dy,xmin,ymin
      real      mdv

c     Parameters
      real      eps
      parameter (eps=0.01)

c     Auxiliary variables
      real      ri,rj
      integer   i,j,ir,ju     
      real      frac0i,frac0j,frac1i,frac1j
      real      tmp,nor

c     Get index
      ri=(lon-xmin)/dx
      rj=(lat-ymin)/dy
      i=int(ri)
      j=int(rj)
      if ((i.lt.0).or.(i.gt.nx).or.
     >    (j.lt.0).or.(j.gt.ny)) then
         print*,'lat/lon interpolation not possible....'
         stop
      endif
      
c     Get inidices of left and upper neighbours
      ir=i+1
      if (ir.gt.nx) ir=nx
      ju=j+1
      if (ju.gt.ny) ju=ny

c     Get the weights for the bilinear interpolation
      frac0i=ri-float(i)
      frac0j=rj-float(j)
      frac1i=1.-frac0i
      frac1j=1.-frac0j

c     Bilinear interpolation with missing data check
      tmp=0.
      nor=0.
      if ( (abs(ar(i ,j )-mdv).gt.eps) ) then
         tmp = tmp + ar(i ,j ) * frac1i * frac1j
         nor = nor + frac1i * frac1j
      endif
      if ( (abs(ar(i ,ju)-mdv).gt.eps) ) then
         tmp = tmp + ar(i ,ju) * frac1i * frac0j
         nor = nor + frac1i * frac0j
      endif
      if ( (abs(ar(ir,j )-mdv).gt.eps) ) then
         tmp = tmp + ar(ir,j ) * frac0i * frac1j
         nor = nor + frac0i * frac1j
      endif
      if ( (abs(ar(ir,ju)-mdv).gt.eps) ) then
         tmp = tmp + ar(ir,ju) * frac0i * frac0j
         nor = nor + frac0i * frac0j
      endif

c     Return result
      int2d=tmp/nor

      end

c     ********************************************************************************
c     * NETCDF INPUT AND OUTPUT                                                      *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Write input field to netcdf
c     --------------------------------------------------------------------------------

      SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz)

c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
c     files are consitent with this grid. 

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 field (0:nx,0:ny,0:nz)
      character*80         fieldname
      character*80         pvsrcfile

c     Auxiliary variables
      integer              cdfid,stat
      integer              vardim(4)
      real                 misdat
      real                 varmin(4),varmax(4),stag(4)
      integer              ndimin,outid,i,j,k
      real                 tmp(0:nx,0:ny,0:nz)
      integer              ntimes
      real                 time(200)       
      integer              nvars
      character*80         vnam(100),varname
      integer              isok

c     Get grid parameters from PV
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname='PV'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)    
      if (stat.ne.0) goto 998
      time(1)=0.
      call gettimes(cdfid,time,ntimes,stat)  
      if (stat.ne.0) goto 998
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 998

c     Save variables (write definition, if necessary)
      call cdfwopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname=fieldname
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndimin,misdat,vardim,
     >               varmin,varmax,stag,stat)
         if (stat.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time(1),0,field,stat)
      print*,'W ',trim(varname),' ',trim(pvsrcfile)
      if (stat.ne.0) goto 998

c     Close input netcdf file
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 998
      
      return

c     Exception handling
 998  print*,'Write_Inp: Problem with input netcdf file... Stop'
      stop

      end


c     --------------------------------------------------------------------------------
c     Read input fields for reference profile
c     --------------------------------------------------------------------------------

      SUBROUTINE read_inp (field,fieldname,pvsrcfile,
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)

c     Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified 
c     by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input 
c     files are consitent with this grid. The missing data value is set to <mdv>.

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 field(0:nx,0:ny,0:nz)
      character*80         fieldname
      character*80         pvsrcfile
      real                 dx,dy,dz,xmin,ymin,zmin
      real                 mdv

c     Numerical and physical parameters
      real                 eps
      parameter            (eps=0.01)

c     Auxiliary variables
      integer              cdfid,stat,cdfid99
      integer              vardim(4)
      real                 misdat
      real                 varmin(4),varmax(4),stag(4)
      integer              ndimin,outid,i,j,k
      real                 tmp(nx,ny,nz)
      integer              ntimes
      real                 time(200)       
      integer              nvars
      character*80         vnam(100),varname
      integer              isok

c     Open the input netcdf file
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 998

c     Check whether needed variables are on file
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 998
      isok=0
      varname=trim(fieldname)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998

c     Get the grid parameters from theta     
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)    
      if (stat.ne.0) goto 998
      time(1)=0.
      call gettimes(cdfid,time,ntimes,stat)  
      if (stat.ne.0) goto 998

c     Check whether grid parameters are consistent
      if ( (vardim(1).ne.(nx+1)).or.
     >     (vardim(2).ne.(ny+1)).or.
     >     (vardim(3).ne.(nz+1)).or.
     >     (abs(varmin(1)-xmin).gt.eps).or.
     >     (abs(varmin(2)-ymin).gt.eps).or.
     >     (abs(varmin(3)-zmin).gt.eps).or.
     >     (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or.
     >     (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or.
     >     (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) 
     >then
         print*,'Input grid inconsitency...'
         print*,'  Nx      = ',vardim(1),nx+1
         print*,'  Ny      = ',vardim(2),ny+1
         print*,'  Nz      = ',vardim(3),nz+1
         print*,'  Varminx = ',varmin(1),xmin
         print*,'  Varminy = ',varmin(2),ymin
         print*,'  Varminz = ',varmin(3),zmin
         print*,'  Dx      = ',(varmax(1)-varmin(1))/real(nx-1),dx
         print*,'  Dy      = ',(varmax(2)-varmin(2))/real(ny-1),dy
         print*,'  Dz      = ',(varmax(3)-varmin(3))/real(nz-1),dz
         goto 998
      endif

c     Load variables
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)
      if (stat.ne.0) goto 998
      call getdat(cdfid,varname,time(1),0,field,stat)
      print*, 'R ',trim(varname),' ',trim(pvsrcfile)
      if (stat.ne.0) goto 998

c     Close input netcdf file
      call clscdf(cdfid,stat)
      if (stat.ne.0) goto 998

c     Set missing data value to <mdv>
      do i=1,nx
         do j=1,ny
            do k=1,nz
               if (abs(field(i,j,k)-misdat).lt.eps) then
                  field(i,j,k)=mdv
               endif
            enddo
         enddo
      enddo

      return

c     Exception handling
 998  print*,'Read_Inp: Problem with input netcdf file... Stop'
      stop

      end

c     --------------------------------------------------------------------------------
c     Check whether variable is found on netcdf file
c     --------------------------------------------------------------------------------

      subroutine check_varok (isok,varname,varlist,nvars)

c     Check whether the variable <varname> is in the list <varlist(nvars)>. If this is 
C     the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.

      implicit none

c     Declaraion of subroutine parameters
      integer      isok
      integer      nvars
      character*80 varname
      character*80 varlist(nvars)

c     Auxiliary variables
      integer      i

c     Main
      do i=1,nvars
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
      enddo

      end

c     --------------------------------------------------------------------------------
c     Get grid parameters
c     --------------------------------------------------------------------------------

      subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >                     pvsrcfile)

c     Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>.
c     The grid parameters are
c        nx,ny,nz                : Number of grid points in x, y and z direction
c        xmin,ymin,zmin          : Minimum domain coordinates in x, y and z direction
c        xmax,ymax,zmax          : Maximal domain coordinates in x, y and z direction
c        dx,dy,dz                : Horizontal and vertical resolution
c     Additionally, it is checked whether the vertical grid is equally spaced. If ok,
c     the grid paramters are transformed from lon/lat to distance (in meters)

      implicit none

c     Declaration of subroutine parameters
      character*80   pvsrcfile
      integer        nx,ny,nz
      real           dx,dy,dz
      real           xmin,ymin,zmin,xmax,ymax,zmax
      real           mdv

c     Numerical epsilon and other physical/geoemtrical parameters
      real           eps
      parameter      (eps=0.01)

c     Auxiliary variables
      integer        cdfid,cstid
      integer        ierr
      character*80   vnam(100),varname
      integer        nvars
      integer        isok
      integer        vardim(4)
      real           misdat
      real           varmin(4),varmax(4),stag(4)
      real           aklev(1000),bklev(1000),aklay(1000),bklay(1000)
      real           dh
      character*80   csn
      integer        ndim
      integer        i

c     Get all grid parameters
      call cdfopn(pvsrcfile,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,ierr)
      if (ierr.ne.0) goto 998
      isok=0
      varname='PV'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 998
      call getcfn(cdfid,csn,ierr)
      if (ierr.ne.0) goto 998
      call cdfopn(csn,cstid,ierr)
      if (ierr.ne.0) goto 998
      call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax,
     >            stag,ierr)
      if (ierr.ne.0) goto 998
      nx=vardim(1)
      ny=vardim(2)
      nz=vardim(3)
      xmin=varmin(1)
      ymin=varmin(2)
      zmin=varmin(3)
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
      if (ierr.ne.0) goto 998
      call getgrid(cstid,dx,dy,ierr)
      if (ierr.ne.0) goto 998
      xmax=varmax(1)
      ymax=varmax(2)
      zmax=varmax(3)
      dz=(zmax-zmin)/real(nz-1)
      call clscdf(cstid,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998

c     Check whether the grid is equallay spaced in the vertical
      do i=1,nz-1
         dh=aklev(i+1)-aklev(i)
         if (abs(dh-dz).gt.eps) then
            print*,'Aklev: Vertical grid must be equally spaced... Stop'
            print*,(aklev(i),i=1,nz)
            stop
         endif
         dh=aklay(i+1)-aklay(i)
         if (abs(dh-dz).gt.eps) then
            print*,'Aklay: Vertical grid must be equally spaced... Stop'
            print*,(aklay(i),i=1,nz)
            stop
         endif
      enddo

c     Set missing data value
      mdv=misdat

      return

c     Exception handling
 998  print*,'Read_Dim: Problem with input netcdf file... Stop'
      stop

      end