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