Subversion Repositories pvinversion.ecmwf

      PROGRAM inv_cart

c     ********************************************************************************
c     * CALCULATES FOR A PV AND THETA DISTRIBUTION OTHER PROGNOSTIC FIELDS BY        *
c     * MEANS OF A PV INVERSION                                                      *
c     *                                                                              *
c     * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006            *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Declaration of variables, parameters, externals and common blocks            
c     --------------------------------------------------------------------------------

      implicit none

c     Grid parameters
      integer        nx,ny,nz
      real           xmin,ymin,zmin,xmax,ymax,zmax
      real           dx,dy,dz
      real           deltax,deltay,deltaz
      real           mdv
      real,          allocatable,dimension (:,:) :: coriol

c     Reference state
      real,   allocatable,dimension (:) :: nsqref
      real,   allocatable,dimension (:) :: thetaref
      real,   allocatable,dimension (:) :: rhoref
      real,   allocatable,dimension (:) :: pressref
      real,   allocatable,dimension (:) :: zref


C     Boundary conditions
      real,   allocatable,dimension (:,:) :: thetatop
      real,   allocatable,dimension (:,:) :: thetabot
      real,   allocatable,dimension (:,:) :: ua
      real,   allocatable,dimension (:,:) :: ub
      real,   allocatable,dimension (:,:) :: va
      real,   allocatable,dimension (:,:) :: vb

c     Potentiual vorticity and stream function
      real,   allocatable,dimension (:,:,:) :: psi
      real,   allocatable,dimension (:,:,:) :: pv

c     Auxiliary arrays for numerics
      real,   allocatable,dimension (:,:,:) :: a
      real,   allocatable,dimension (:,:,:) :: b
      real,   allocatable,dimension (:,:,:) :: c

c     Input parameters
      character*80         pvsrcfile
      character*80         referfile

c     Auxiliary variables
      integer             i,j,k
      integer             stat

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

      print*,'********************************************************'
      print*,'* INV_CART                                             *'
      print*,'********************************************************'

c     Read parameter file
      open(10,file='fort.10')
        read(10,*) pvsrcfile
        read(10,*) referfile
      close(10)
      print*
      print*,trim(pvsrcfile)


c     Get lat/lon gid parameters from input file
      call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv,
     >               pvsrcfile)
      print*
      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
      nx=nx-1
      ny=ny-1
      nz=nz-1
      
c     Allocate memory for boundary conditions
      allocate(thetatop(0:nx,0:ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array thetatop ***'
      allocate(thetabot(0:nx,0:ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array thetabot ***'
      allocate(ua(0:nx,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ua ***'
      allocate(ub(0:nx,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ub ***'
      allocate(va(0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array va ***'
      allocate(vb(0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array vb ***'

c     Allocate memory for 3d PV and stream function
      allocate(psi(0:nx,0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array psi ***'
      allocate(pv(0:nx,0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pv ***'

c     Alllocate memory for matrix elements for inversion operator
      allocate(a(0:nx,0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array a ***'
      allocate(b(0:nx,0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array b ***'
      allocate(c(0:nx,0:ny,0:nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array c ***'

c     Allocate memory for reference profile
      allocate(nsqref(0:2*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array nsqref ***'
      allocate(thetaref(0:2*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array thetaref ***'
      allocate(rhoref(0:2*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array rhoref ***'
      allocate(pressref(0:2*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pressref ***'
      allocate(zref(0:2*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zref ***'

c     Allocate memory for Coriolis parameter
      allocate(coriol(0:nx,0:ny),STAT=stat)
      if (stat.ne.0) print*,'error allocating coriol'

c     --------------------------------------------------------------------------------
c     Input
c     --------------------------------------------------------------------------------

c     Read reference profile and grid parameters
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
     >               nx,ny,nz,deltax,deltay,deltaz,coriol,
     >               referfile)
      deltax=1000.*deltax
      deltay=1000.*deltay

c     Read input fields from netcdf 
      call read_inp (pv,thetabot,thetatop,ua,ub,va,vb,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,
     >               pvsrcfile)

c     --------------------------------------------------------------------------------
c     Perform the inversion
c     --------------------------------------------------------------------------------

C     Init matrix elements for inversion
      call matrixel(a,b,c,coriol,nx,ny,nz,nsqref,rhoref,
     >              deltax,deltay,deltaz)

c     Inversion
      call sor(psi,nsqref,thetatop,thetabot,thetaref,rhoref,coriol,
     >         pv,ua,ub,va,vb,a,b,c,nx,ny,nz,deltax,deltay,deltaz)

c     --------------------------------------------------------------------------------
c     Output
c     --------------------------------------------------------------------------------

c     Write output to netcdf
      print*
      call write_out (psi,thetabot,thetatop,ua,ub,va,vb,
     >                nx,ny,nz,deltax,deltay,deltaz,
     >                coriol,thetaref,rhoref,pressref,pvsrcfile)


      END


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

c     --------------------------------------------------------------------------------
c     Output to netcdf file
c     --------------------------------------------------------------------------------

      SUBROUTINE write_out (psi,thetabot,thetatop,ua,ub,va,vb,nx,ny,nz,
     >                      dx,dy,dz,coriol,thetaref,rhoref,pref,
     >                      pvsrcfile)

c     Write the result of the inversion to output netcdf
c     
c        psi               : streamm function as calculated from inversion
c        thetabot,thetatop : potential temperature at lower and upper boundary
c        ua,ub             : Zonal wind at western and eastern boundary
c        va,vb             : Meridional wind at southern and northern boundary
c        nx,ny,nz          : Grid dimension in x, y and z direction
c        dx,dy,dz          : Grid resolution
c        coriol            : Coriolis parameter
c        thetaref          : Reference profile of potential temperature
c        rhoref            : Reference profile of density
c        pref              : Reference profile of pressure
c        pvsrcfile         : Name of the output netcdf file

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 psi(0:nx,0:ny,0:nz)
      real                 thetatop(0:nx,0:ny)
      real                 thetabot(0:nx,0:ny)
      real                 ua(0:nx,0:nz)
      real                 ub(0:nx,0:nz)
      real                 va(0:ny,0:nz)
      real                 vb(0:ny,0:nz)
      character*(30)       pvsrcfile
      real                 dx,dy,dz
      real                 thetaref(0:2*nz)
      real                 rhoref(0:2*nz)
      real                 pref(0:2*nz)
      real                 coriol(0:nx,0:ny)
      
c     Numerical and physical parameters
      real                 eps
      parameter            (eps=0.01)
      real                 g
      parameter            (g=9.81)
      real                 preunit
      parameter            (preunit=0.01)
      real                 kappa
      parameter            (kappa=0.286)

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

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

c     Open output netcdf file
      call cdfwopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 997
      
c     Create the variable if necessary
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 997
      isok=0
      varname='PSI'
      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 997
      endif
      isok=0
      varname='U'
      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 997
      endif
      isok=0
      varname='V'
      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 997
      endif
      isok=0
      varname='TH'
      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 997
      endif
      isok=0
      varname='T'
      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 997
      endif
      isok=0
      varname='P'
      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 997
      endif

c     Write stream function
      varname='PSI'
      call putdat(cdfid,varname,time(1),0,psi,ierr)
      if (stat.ne.0) goto 997
      print*,'W PSI      ',trim(pvsrcfile)

c     Calculate and write velocity U: keep southern and northern boundary
      do k=0,nz
         do i=0,nx
            do j=1,ny-1
               tmp(i,j,k)=(psi(i,j-1,k)-psi(i,j+1,k))/
     >                    (2.*dy*coriol(i,j))
            enddo
            tmp(i,0,k) =ua(i,k)
            tmp(i,ny,k)=ub(i,k)
         enddo
      enddo
      varname='U'
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
      if (stat.ne.0) goto 997
      print*,'W U        ',trim(pvsrcfile)

C     Calculate and write velocity V: keep western and eastern boundary
      do k=0,nz
         do j=0,ny
            do i=1,nx-1
               tmp(i,j,k)=(psi(i+1,j,k)-psi(i-1,j,k))/
     >                    (2.*dx*coriol(i,j))
            enddo
            tmp(0,j,k)=va(j,k)
            tmp(nx,j,k)=vb(j,k)
         enddo
      enddo
      varname='V'
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
      if (stat.ne.0) goto 997
      print*,'W V        ',trim(pvsrcfile)
      
c     Calculate and write potential temperature: keep lower and upper boundary
c     Potential temperature is needed for calculation of temperature: keep it
      do i=0,nx
         do j=0,ny
            th(i,j,0)=thetabot(i,j)
            do k=1,nz-1      
              th(i,j,k)=thetaref(2*k)*
     >                  (psi(i,j,k+1)-psi(i,j,k-1))/(2.*dz*g)

           enddo
           th(i,j,nz)=thetatop(i,j)
         enddo
      enddo    
      varname='TH'
      call putdat(cdfid,varname,time(1),0,th,ierr)
      if (stat.ne.0) goto 997
      print*,'W TH       ',trim(pvsrcfile)

c     Calculate and write pressure: The pressure is directly proportional to the
c     streamfunction. But the streamfunction is determined only up to an additive
c     constant. Shift the streamfunction in such a way that it vanish in the mean
c     on the boundaries. Pressure is needed for calculation of temperature: keep it
      meanpsi=0.
      meancnt=0.
      do i=0,nx
         do j=0,ny
            meanpsi=meanpsi+psi(i,j,0)+psi(i,j,nz)
            meancnt=meancnt+2
         enddo
      enddo
      do i=0,nx
         do k=0,nz
            meanpsi=meanpsi+psi(i,0,k)+psi(i,ny,k)
            meancnt=meancnt+2 
         enddo
      enddo
      do j=0,ny
         do k=0,nz
            meanpsi=meanpsi+psi(0,j,k)+psi(nx,j,k)
            meancnt=meancnt+2 
         enddo
      enddo
      meanpsi=meanpsi/meancnt
      do i=0,nx
         do j=0,ny
            do k=0,nz
               kk=2*k
               pr(i,j,k)=preunit*rhoref(kk)*(psi(i,j,k)-meanpsi)
            enddo
         enddo
      enddo
      varname='P'
      call putdat(cdfid,varname,time(1),0,pr,ierr)
      if (stat.ne.0) goto 997
      print*,'W P        ',trim(pvsrcfile)

c     Calculate and write temperature
      do i=0,nx
         do j=0,ny
            do k=0,nz      
               kk=2*k
               tmp(i,j,k)=((pref(kk)/1000.)**kappa) * 
     >                (th(i,j,k)+kappa*thetaref(kk)*pr(i,j,k)/pref(kk))
           enddo
        enddo
      enddo      
      varname='T'
      call putdat(cdfid,varname,time(1),0,tmp,ierr)
      if (stat.ne.0) goto 997
      print*,'W T        ',trim(pvsrcfile)

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

      return

c     Exception handling
 997  print*,'Problem with output netcdf file... Stop'
      stop

      end


c     --------------------------------------------------------------------------------
c     Read the reference file
c     --------------------------------------------------------------------------------

      SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
     >                     nx,ny,nz,deltax,deltay,deltaz,coriol,
     >                     pvsrcfile) 

c     Read the reference profile from file
c
c         thetaref             : Reference potential temperature (K)
c         pressref             : Reference pressure (Pa)
c         rhoref               : Reference density (kg/m^3)
c         nsqref               : Stratification (s^-1)
c         zref                 : Reference height (m)
c         nx,nny,nz            : Grid dimension in x,y,z direction
c         deltax,deltay,deltaz : Grid spacings used for calculations (m)
c         coriol               : Coriolis parameter (s^-1)
c         pvsrcfile            : Input file

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 nsqref  (0:2*nz)
      real                 thetaref(0:2*nz)
      real                 rhoref  (0:2*nz)
      real                 pressref(0:2*nz)
      real                 zref    (0:2*nz)
      real                 deltax,deltay,deltaz
      real                 coriol  (0:nx,0:ny)
      character*80         pvsrcfile

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

c     Auxiliary variables
      integer              cdfid,stat
      integer              vardim(4)
      real                 misdat
      integer              ndimin
      real                 varmin(4),varmax(4),stag(4)
      integer              i,j,k,nf1
      integer              ntimes
      real                 time(200)  
      character*80         vnam(100),varname
      integer              nvars
      integer              isok,ierr
      real                 x(0:nx,0:ny),y(0:nx,0:ny)
      real                 mean,count

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

c     Open output netcdf file
      call cdfopn(pvsrcfile,cdfid,stat)
      if (stat.ne.0) goto 997
      
c     Create the variable if necessary
      call getvars(cdfid,nvars,vnam,stat)
      if (stat.ne.0) goto 997

c     Read data from netcdf file
      isok=0
      varname='NSQREF'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,nsqref,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='RHOREF'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,rhoref,stat)
      if (stat.ne.0) goto 997
 
      isok=0
      varname='THETAREF'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,thetaref,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='PREREF'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,pressref,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='ZREF'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,zref,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='CORIOL'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,coriol,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='X'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,x,stat)
      if (stat.ne.0) goto 997

      isok=0
      varname='Y'
      print*,'R ',trim(varname),' ',trim(pvsrcfile)
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) goto 997
      call getdat(cdfid,varname,time(1),0,y,stat)
      if (stat.ne.0) goto 997

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

c     Determine the grid spacings <deltax, deltay, deltaz>
      mean=0.
      count=0.
      do i=1,nx
         do j=0,ny
            mean=mean+abs(x(i,j)-x(i-1,j))
            count=count+1.
         enddo
      enddo
      deltax=mean/count

      mean=0.
      count=0.
      do j=1,ny
         do i=0,nx
            mean=mean+abs(y(i,j)-y(i,j-1))
            count=count+1.
         enddo
      enddo
      deltay=mean/count

      mean=0.
      count=0.
      do k=1,nz-1
         mean=mean+abs(zref(k+1)-zref(k-1))
         count=count+1.
      enddo
      deltaz=mean/count

      return

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

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


c     -------------------------------------------------------------------------------
c     Read the input netcdf file
c     --------------------------------------------------------------------------------

      SUBROUTINE read_inp (pv,thetabot,thetatop,ua,ub,va,vb,
     >                     nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,
     >                     pvsrcfile)

c     Read all needed field from netcdf file <pvsrcfile>. The input fields are:
c        pv                : quasigeostrophic potential vorticity
c        thetabot,thetatop : potential temperature at lower and upper boundary
c        ua,ub             : Zonal wind at western and eastern boundary
c        va,vb             : Meridional wind at southern and northern boundary
c     The grid is specified by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed
c     whether the input files are consitent with this grid. The input netcdf file must
c     contain the variables <QGPV,THETA,U,V>. If the netcdf file also contains the fields
c     <DQGPV,DTHETA,DU,DV>, these increments are added to <QGPV,THETA,U,V>.

      implicit   none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 pv(0:nx,0:ny,0:nz)
      real                 thetatop(0:nx,0:ny)
      real                 thetabot(0:nx,0:ny)
      real                 ua(0:nx,0:nz)
      real                 ub(0:nx,0:nz)
      real                 va(0:ny,0:nz)
      real                 vb(0:ny,0:nz)
      character*(30)       pvsrcfile
      real                 dx,dy,dz,xmin,ymin,zmin

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

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                 max_th
      real                 tmp(0:nx,0:ny,0: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='TH'
      call check_varok(isok,varname,vnam,nvars)
      varname='U'
      call check_varok(isok,varname,vnam,nvars)
      varname='V'
      call check_varok(isok,varname,vnam,nvars)
      varname='QGPV'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.ne.4) goto 998

c     Get the grid parameters 
      varname='QGPV'
      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 consitent
      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(nx)-dx).gt.eps).or.
     >     (abs((varmax(2)-varmin(2))/real(ny)-dy).gt.eps).or.
     >     (abs((varmax(3)-varmin(3))/real(nz)-dz).gt.eps) ) then
         print*,'Input grid inconsitency...'
         print*,'xmin : ',xmin,varmin(1)
         print*,'ymin : ',ymin,varmin(2)
         print*,'zmin : ',zmin,varmin(3)
         print*,'dx   : ',dx,(varmax(1)-varmin(1))/real(nx)
         print*,'dy   : ',dy,(varmax(2)-varmin(2))/real(ny)
         print*,'dz   : ',dz,(varmax(3)-varmin(3))/real(nz)
         print*,'nx   : ',nx
         print*,'ny   : ',ny
         print*,'nz   : ',nz
         goto 998
      endif

c     THETA: Load upper and lower boundary values
      varname='TH'
      call getdat(cdfid,varname,time(1),0,tmp,stat)  
      print*,'R TH       ',trim(pvsrcfile)
      if (stat.ne.0) goto 998
      do i=0,nx
         do j=0,ny
            if ( abs(tmp(i,j,0)-misdat).lt.eps ) then
               thetabot(i,j)=0.
            else
               thetabot(i,j)=tmp(i,j,0)     
            endif
            if ( abs(tmp(i,j,nz)-misdat).lt.eps ) then
               thetatop(i,j)=0.
            else
               thetatop(i,j)=tmp(i,j,nz)
            endif
         enddo
      enddo

c     U: Load zonal velocity at southern and northern boundary
      varname='U'
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)
      if (stat.ne.0) goto 998
      call getdat(cdfid,varname,time(1),0,tmp,stat)
      print*,'R U        ',trim(pvsrcfile)
      if (stat.ne.0) goto 998
      do i=0,nx
         do k=0,nz
            if ( abs(tmp(i,0,k)-misdat).lt.eps ) then
               ua(i,k)=0.
            else
               ua(i,k)=tmp(i,0,k)
            endif
            if ( abs(tmp(i,ny,k)-misdat).lt.eps ) then
               ub(i,k)=0.
            else
               ub(i,k)=tmp(i,ny,k)
           endif
         enddo
      enddo

c     Load meridional velocity at western and eastern boundary
      varname='V'
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)
      if (stat.ne.0) goto 998
      call getdat(cdfid,varname,time(1),0,tmp,stat)
      print*,'R V        ',trim(pvsrcfile)
      if (stat.ne.0) goto 998
      do j=0,ny
         do k=0,nz
            if ( abs(tmp(0,j,k)-misdat).lt.eps ) then
               va(j,k)=0.
            else
               va(j,k)=tmp(0,j,k)
            endif
            if ( abs(tmp(nx,j,k)-misdat).lt.eps ) then
               vb(j,k)=0.
            else
               vb(j,k)=tmp(nx,j,k)
            endif
         enddo
      enddo      

c     Load qgPV
      varname='QGPV'
      call getdef(cdfid,varname,ndimin,misdat,vardim,
     >            varmin,varmax,stag,stat)
      if (stat.ne.0) goto 998
      call getdat(cdfid,varname,time(1),0,tmp,stat)
      print*,'R QGPV     ',trim(pvsrcfile)
      if (stat.ne.0) goto 998
      do i=0,nx
         do j=0,ny
            do k=0,nz
               if ( abs(tmp(i,j,k)-misdat).lt.eps ) then
                  pv(i,j,k)=0.
               else
                  pv(i,j,k)=tmp(i,j,k)
               endif
            enddo
         enddo
      enddo

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

      return

c     Exception handling
 998  print*,'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     * INVERSION ROUTINES                                                           *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     SOR algorithm (successive over relaxation)
c     --------------------------------------------------------------------------------

      SUBROUTINE sor(psi,nsq,thetatop,thetabot,thetaref,rhoref,
     >               coriol,pv,ua,ub,va,vb,a,b,c,nx,ny,nz,dx,dy,dz)

c     Solve the qgPV equation by succesive over relaxation (SOR). The subroutine
c     parameters are:
c
c        psi                 : Streamfunction, i.e. result of the PV inversion
c        nsq,rhoref,thetaref : Reference profile
c        thetatop,thetabot   : Upper and lower boundary condition
c        pv                  : quasigeostrophic potential vorticity (qgPV)
c        ua,ub,va,vb         : lateral boundary condition for wind
c        a,b,c               : Matrices for the inversion operator
c        nx,ny,nz,dx,dy,dz   : Grid specification
c        coriol              : Coriolis parameter

      implicit none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 dx,dy,dz
      real                 psi     (0:nx,0:ny,0:nz)
      real                 nsq     (0:2*nz)
      real                 thetatop(0:nx,0:ny)
      real                 thetabot(0:nx,0:ny)
      real                 thetaref(0:2*nz)
      real                 rhoref  (0:2*nz)
      real                 pv      (0:nx,0:ny,0:nz)
      real                 ua      (0:nx,0:nz)
      real                 ub      (0:nx,0:nz)
      real                 va      (0:ny,0:nz)
      real                 vb      (0:ny,0:nz)
      real                 a       (0:nx,0:ny,0:nz)
      real                 b       (0:nx,0:ny,0:nz)
      real                 c       (0:nx,0:ny,0:nz)
      real                 coriol  (0:nx,0:ny)

c     Numerical and physical parameters
      real                 maxspec
      parameter            (maxspec=2.0)
      integer              nofiter
      parameter            (nofiter=500)
      real                 omega
      parameter            (omega=1.81)
     
c     Auxiliary variables
      integer              counter
      integer              i,j,k
      real                 deltasq,psigauge
      real                 specx,specy,specz
      real                 helpx,helpy,helpz

c     Init the output array
      do i=0,nx
         do j=0,ny
            do k=0,nz
               psi(i,j,k)=0.
            enddo
         enddo
      enddo

c     Calculate the spectrum of the matrix 
      i=nx/2
      specx=4.*a(i,0,0)/
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
      specy=2.*(b(i,0,0)+b(i,1,0))/
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
      specz=2.*(c(i,0,0)+c(i,0,1))/
     >      (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1))
      do k=1,nz-2
         do j=1,ny-2

            helpx=4.*a(i,j,k)/
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
            if (helpx.gt.specx) specx=helpx
                       
            helpy=2.*(b(i,j,k)+b(i,j+1,k))/
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
            if (helpy.gt.specy) specy=helpy

            helpz=2.*(c(i,j,k)+c(i,j,k+1))/
     >            (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1))
            if (helpz.gt.specz) specz=helpz

         enddo
      enddo

c     Check whether the dimensions of the grid are sufficient
      print *
      print *, 'Spectrum of the matrix in each direction '
      print *, 'Spectrum = ', specx, specy, specz
      print *
      if ((maxspec*specx.lt.specy).or.(maxspec*specx.lt.specz)) then
         print*,' Nx too small... Stop'
         stop
      endif
      if ((maxspec*specy.lt.specx).or.(maxspec*specy.lt.specz)) then
         print*,'Ny too small... Stop'
         stop
      endif
      if ((maxspec*specz.lt.specx).or.(maxspec*specz.lt.specy)) then
         print*,'Nz too small... Stop'
         stop
      endif

c     Calculate error: control variable for the iteration
      psigauge=0.
      deltasq=0.
      do k=1,nz-1
         do i=1,nx-1
            do j=1,ny-1
               deltasq=deltasq+(-pv(i,j,k)+(
     >              a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)-
     >              2.*psi(i,j,k)) +
     >              b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))-
     >              b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ 
     >              c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))-
     >              c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1))
     >              )/(dx*dy*dz*rhoref(2*k)))**2.
  
            enddo
         enddo
      enddo
      print 102, 'psigauge', psigauge, 'deltasq',
     >           deltasq/(real(nx)*real(ny)*real(nz))

c     Iterations
      do counter=1,nofiter

C        Perform one iteration step
         call psiappsor(omega,pv,psi,nsq,rhoref,thetatop,
     >                  thetabot,thetaref,coriol,ua,ub,va,vb,
     >                  a,b,c,nx,ny,nz,dx,dy,dz)


c        Adjustment
         if (mod(counter,100).eq.0) then
            psigauge=0.
            do i=0,nx
               do j=0,ny
                  if (psi(i,j,0).lt.psigauge) then 
                     psigauge=psi(i,j,0)
                  endif
               enddo
            enddo
            do k=0,nz
               do i=0,nx
                  do j=0,ny
                     psi(i,j,k)=psi(i,j,k)-psigauge
                  enddo
               enddo
            enddo
         endif

c        Calculate error: control variable for the iteration        
         if (mod(counter,nofiter/10).eq.0) then
            deltasq=0.
            do k=1,nz-1
               do i=1,nx-1
                  do j=1,ny-1
                     deltasq=deltasq+(-pv(i,j,k)+(
     >                 a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)-
     >                    2.*psi(i,j,k)) +
     >                 b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))-
     >                 b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ 
     >                 c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))-
     >                 c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1))
     >                    )/(dx*dy*dz*rhoref(2*k)))**2.
                  enddo
               enddo
            enddo
            print 102, 'psigauge', psigauge, 'deltasq',
     >           deltasq/(real(nx)*real(ny)*real(nz))
         endif

      enddo

      return

c     Format specifications
 102  format (a11, ' = ',e10.3,a11, ' = ',e10.3)

      end

c     --------------------------------------------------------------------------------
c     SOR algorithm (successive over relaxation)
c     --------------------------------------------------------------------------------

      subroutine psiappsor(omega,pv,psi,nsq,rhoref,thetatop,
     >                     thetabot,thetaref,coriol,ua,ub,va,vb,
     >                     a,b,c,nx,ny,nz,dx,dy,dz)

c     Perform one relaxation step
c
c        psi                 : Streamfunction, i.e. result of the PV inversion
c        nsq,rhoref,thetaref : Reference profile
c        thetatop,thetabot   : Upper and lower boundary condition
c        pv                  : quasigeostrophic potential vorticity (qgPV)
c        ua,ub,va,vb         : lateral boundary condition for wind
c        a,b,c               : Matrices for the inversion operator
c        nx,ny,nz,dx,dy,dz   : Grid specification
c        nofiter             : Number of iterations
c        omega               : Relaxation parameter
c        coriol              : Coriolis parameter

      implicit none

c     Declaration of subroutine parameters
      integer              nx,ny,nz
      real                 pv(0:nx,0:ny,0:nz)
      real                 psi(0:nx,0:ny,0:nz)
      real                 nsq(0:2*nz)
      real                 rhoref(0:2*nz)
      real                 thetatop(0:nx,0:ny)
      real                 thetabot(0:nx,0:ny)
      real                 thetaref(0:2*nz)
      real                 ua(0:nx,0:nz)
      real                 ub(0:nx,0:nz)
      real                 va(0:ny,0:nz)
      real                 vb(0:ny,0:nz)
      real                 a(0:nx,0:ny,0:nz)
      real                 b(0:nx,0:ny,0:nz)
      real                 c(0:nx,0:ny,0:nz)
      real                 coriol(0:nx,0:ny)
      real                 dx,dy,dz
      real                 omega

c     Numerical and physical parameters
      real                 g
      parameter            (g=9.81)

c     Auxiliary variables
      integer              i,j,k
      real                 dxy,dxz,dyz,dxyz

c     Set the area and volume infinitesimals for integration
      dxy=dx*dy
      dxz=dx*dz
      dyz=dy*dz
      dxyz=dx*dy*dz

c     Inner
      do k=1,nz-1
         do i=1,nx-1
            do j=1,ny-1
               psi(i,j,k)=omega*(-dxyz*
     >                    rhoref(2*k)*pv(i,j,k)+a(i,j,k)*
     >                    (psi(i+1,j,k)+psi(i-1,j,k))+b(i,j,k)*
     >                    psi(i,j+1,k)+b(i,j-1,k)*psi(i,j-1,k)+c(i,j,k)*
     >                    psi(i,j,k+1)+c(i,j,k-1)*psi(i,j,k-1))/
     >                    (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+
     >                    c(i,j,k))+(1.-omega)*psi(i,j,k)
            enddo
         enddo
      enddo

c     ZY plane
      do k=1,nz-1
         do j=1,ny-1
            psi(0,j,k)=omega*(-dyz*
     >           rhoref(2*k)*(dx*pv(0,j,k)+va(j,k))+
     >           a(0,j,k)*psi(1,j,k)+
     >           b(0,j,k)*psi(0,j+1,k)+b(0,j-1,k)*psi(0,j-1,k)+
     >           c(0,j,k)*psi(0,j,k+1)+c(0,j,k-1)*psi(0,j,k-1))/
     >           (a(0,j,k)+b(0,j,k)+b(0,j-1,k)+c(0,j,k-1)+c(0,j,k))
     >           +(1.-omega)*psi(0,j,k)
c     
            psi(nx,j,k)=omega*(-dyz*
     >           rhoref(2*k)*(dx*pv(nx,j,k)-vb(j,k))+
     >           a(nx,j,k)*psi(nx-1,j,k)+
     >           b(nx,j,k)*psi(nx,j+1,k)+b(nx,j-1,k)*psi(nx,j-1,k)+
     >           c(nx,j,k)*psi(nx,j,k+1)+c(nx,j,k-1)*psi(nx,j,k-1))/
     >           (a(nx,j,k)+b(nx,j-1,k)+b(nx,j,k)+c(nx,j,k-1)+c(nx,j,k))
     >           +(1.-omega)*psi(nx,j,k)
         enddo
      enddo

c     ZX plane
      do k=1,nz-1
         do i=1,nx-1
            psi(i,0,k)=omega*(-dxz*
     >           rhoref(2*k)*(dy*pv(i,0,k)-ua(i,k))+
     >           a(i,0,k)*(psi(i+1,0,k)+psi(i-1,0,k))+
     >           b(i,0,k)*psi(i,1,k)+
     >           c(i,0,k)*psi(i,0,k+1)+c(i,0,k-1)*psi(i,0,k-1))/
     >           (2.*a(i,0,k)+b(i,0,k)+c(i,0,k-1)+c(i,0,k))
     >           +(1.-omega)*psi(i,0,k)
c     
            psi(i,ny,k)=omega*(-dxz*
     >           rhoref(2*k)*(dy*pv(i,ny,k)+ub(i,k))+
     >           a(i,ny-1,k)*(psi(i+1,ny,k)+psi(i-1,ny,k))+
     >           b(i,ny-1,k)*psi(i,ny-1,k)+
     >           c(i,ny-1,k)*psi(i,ny,k+1)+c(i,ny-1,k-1)*
     >           psi(i,ny,k-1))/(2.*a(i,ny-1,k)+b(i,ny-1,k)+
     >           c(i,ny-1,k-1)+c(i,ny-1,k))
     >           +(1.-omega)*psi(i,ny,k)
         enddo
      enddo

c     XY plane
      do i=1,nx-1
         do j=1,ny-1
            psi(i,j,0)=omega*(-dxy*rhoref(0)*(
     >           dz*pv(i,j,0)+g*coriol(i,j)*thetabot(i,j)/
     >           (nsq(0)*thetaref(0)))+
     >           a(i,j,0)*(psi(i+1,j,0)+psi(i-1,j,0))+
     >           b(i,j,0)*psi(i,j+1,0)+b(i,j-1,0)*psi(i,j-1,0)+
     >           c(i,j,0)*psi(i,j,1))/
     >           (2.*a(i,j,0)+b(i,j-1,0)+b(i,j,0)+c(i,j,0))
     >           +(1.-omega)*psi(i,j,0)
c     
            psi(i,j,nz)=omega*(-dxy*rhoref(2*nz)*(
     >           dz*pv(i,j,nz)-g*coriol(i,j)*thetatop(i,j)/
     >           (nsq(2*nz)*thetaref(2*nz)))+
     >           a(i,j,nz)*(psi(i+1,j,nz)+psi(i-1,j,nz))+
     >           b(i,j,nz)*psi(i,j+1,nz)+b(i,j-1,nz)*psi(i,j-1,nz)+
     >           c(i,j,nz-1)*psi(i,j,nz-1))/
     >           (2.*a(i,j,nz)+b(i,j-1,nz)+b(i,j,nz)+c(i,j,nz-1))
     >           +(1.-omega)*psi(i,j,nz)
         enddo
      enddo

c     Y edges
      do j=1,ny-1
         psi(0,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(0,j,0)+
     >        dz*va(j,0)+dx*g*coriol(0,j)*thetabot(0,j)/
     >        (nsq(0)*thetaref(0)))+
     >        a(0,j,0)*psi(1,j,0)+
     >        b(0,j,0)*psi(0,j+1,0)+b(0,j-1,0)*psi(0,j-1,0)+
     >        c(0,j,0)*psi(0,j,1))/
     >        (a(0,j,0)+b(0,j-1,0)+b(0,j,0)+c(0,j,0))
     >        +(1.-omega)*psi(0,j,0)
c     
         psi(nx,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(nx,j,0)-
     >        dz*vb(j,0)+dx*g*coriol(nx,j)*thetabot(nx,j)/
     >        (nsq(0)*thetaref(0)))+
     >        a(nx,j,0)*psi(nx-1,j,0)+
     >        b(nx,j,0)*psi(nx,j+1,0)+b(nx,j-1,0)*psi(nx,j-1,0)+
     >        c(nx,j,0)*psi(nx,j,1))/
     >        (a(nx,j,0)+b(nx,j-1,0)+b(nx,j,0)+c(nx,j,0))
     >        +(1.-omega)*psi(nx,j,0)
c
         psi(0,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(0,j,nz)+
     >        dz*va(j,nz)-dx*g*coriol(0,j)*thetatop(0,j)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(0,j,nz)*psi(1,j,nz)+
     >        b(0,j,nz)*psi(0,j+1,nz)+b(0,j-1,nz)*psi(0,j-1,nz)+
     >        c(0,j,nz-1)*psi(0,j,nz-1))/
     >        (a(0,j,nz)+b(0,j-1,nz)+b(0,j,nz)+c(0,j,nz-1))
     >        +(1.-omega)*psi(0,j,nz)
c     
         psi(nx,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(nx,j,nz)-
     >        dz*vb(j,nz)-dx*g*coriol(nx,j)*thetatop(nx,j)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(nx,j,nz)*psi(nx-1,j,nz)+
     >        b(nx,j,nz)*psi(nx,j+1,nz)+b(nx,j-1,nz)*psi(nx,j-1,nz)+
     >        c(nx,j,nz-1)*psi(nx,j,nz-1))/
     >        (a(nx,j,nz)+b(nx,j-1,nz)+b(nx,j,nz)+c(nx,j,nz-1))
     >        +(1.-omega)*psi(nx,j,nz)
      enddo

c     X edges
      do i=1,nx-1
         psi(i,0,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,0,0)-
     >        dz*ua(i,0)+dy*g*coriol(i,0)*thetabot(i,0)/
     >        (nsq(0)*thetaref(0)))+
     >        a(i,0,0)*(psi(i+1,0,0)+psi(i-1,0,0))+
     >        b(i,0,0)*psi(i,1,0)+
     >        c(i,0,0)*psi(i,0,1))/
     >        (2.*a(i,0,0)+b(i,0,0)+c(i,0,0))
     >        +(1.-omega)*psi(i,0,0)
c     
         psi(i,ny,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,ny,0)+
     >        dz*ub(i,0)+dy*g*coriol(i,ny)*thetabot(i,ny)/
     >        (nsq(0)*thetaref(0)))+
     >        a(i,ny,0)*(psi(i+1,ny,0)+psi(i-1,ny,0))+
     >        b(i,ny-1,0)*psi(i,ny-1,0)+
     >        c(i,ny,0)*psi(i,ny,1))/
     >        (2.*a(i,ny,0)+b(i,ny-1,0)+c(i,ny,0))
     >        +(1.-omega)*psi(i,ny,0)
c     
         psi(i,0,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,0,nz)-
     >        dz*ua(i,nz)-dy*g*coriol(i,0)*thetatop(i,0)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(i,0,nz)*(psi(i+1,0,nz)+psi(i-1,0,nz))+
     >        b(i,0,nz)*psi(i,1,nz)+
     >        c(i,0,nz-1)*psi(i,0,nz-1))/
     >        (2.*a(i,0,nz)+b(i,0,nz)+c(i,0,nz-1))
     >        +(1.-omega)*psi(i,0,nz)
c     
         psi(i,ny,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,ny,nz)+
     >        dz*ub(i,nz)-dy*g*coriol(i,ny)*thetatop(i,ny)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(i,ny,nz)*(psi(i+1,ny,nz)+psi(i-1,ny,nz))+
     >        b(i,ny-1,nz)*psi(i,ny-1,nz)+
     >        c(i,ny,nz-1)*psi(i,ny,nz-1))/
     >        (2.*a(i,ny,nz)+b(i,ny-1,nz)+c(i,ny,nz-1))
     >        +(1.-omega)*psi(i,ny,nz)
      enddo

c     Z edges
      do k=1,nz-1
         psi(0,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,0,k)+
     >        dy*va(0,k)-dx*ua(0,k))+
     >        a(0,0,k)*psi(1,0,k)+
     >        b(0,0,k)*psi(0,1,k)+
     >        c(0,0,k)*psi(0,0,k+1)+c(0,0,k-1)*psi(0,0,k-1))/
     >        (a(0,0,k)+b(0,0,k)+c(0,0,k-1)+c(0,0,k))
     >        +(1.-omega)*psi(0,0,k)
c
         psi(nx,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,0,k)-
     >        dy*vb(0,k)-dx*ua(nx,k))+
     >        a(nx,0,k)*psi(nx-1,0,k)+
     >        b(nx,0,k)*psi(nx,1,k)+
     >        c(nx,0,k)*psi(nx,0,k+1)+c(nx,0,k-1)*psi(nx,0,k-1))/
     >        (a(nx,0,k)+b(nx,0,k)+c(nx,0,k-1)+c(nx,0,k))
     >        +(1.-omega)*psi(nx,0,k)
c
         psi(0,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,ny,k)+
     >        dy*va(ny,k)+dx*ub(0,k))+
     >        a(0,ny,k)*psi(1,ny,k)+
     >        b(0,ny-1,k)*psi(0,ny-1,k)+
     >        c(0,ny,k)*psi(0,ny,k+1)+c(0,ny,k-1)*psi(0,ny,k-1))/
     >        (a(0,ny,k)+b(0,ny-1,k)+c(0,ny,k-1)+c(0,ny,k))
     >        +(1.-omega)*psi(0,ny,k)
c
         psi(nx,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,ny,k)-
     >        dy*vb(ny,k)+dx*ub(nx,k))+
     >        a(nx,ny,k)*psi(nx-1,ny,k)+
     >        b(nx,ny-1,k)*psi(nx,ny-1,k)+
     >        c(nx,ny,k)*psi(nx,ny,k+1)+c(nx,ny,k-1)*psi(nx,ny,k-1))/
     >        (a(nx,ny,k)+b(nx,ny-1,k)+c(nx,ny,k-1)+c(nx,ny,k))
     >        +(1.-omega)*psi(nx,ny,k)
      enddo

c     Points
      psi(0,0,0)=omega*(-rhoref(0)*(dxyz*pv(0,0,0)+dyz*va(0,0)-
     >        dxz*ua(0,0)+dxy*g*coriol(0,0)*thetabot(0,0)/
     >        (nsq(0)*thetaref(0)))+
     >        a(0,0,0)*psi(1,0,0)+
     >        b(0,0,0)*psi(0,1,0)+
     >        c(0,0,0)*psi(0,0,1))/
     >        (a(0,0,0)+b(0,0,0)+c(0,0,0))+
     >        (1.-omega)*psi(0,0,0)
c
      psi(nx,0,0)=omega*(-rhoref(0)*(dxyz*pv(nx,0,0)-dyz*vb(0,0)-
     >        dxz*ua(nx,0)+dxy*g*coriol(nx,0)*thetabot(nx,0)/
     >        (nsq(0)*thetaref(0)))+
     >        a(nx,0,0)*psi(nx-1,0,0)+
     >        b(nx,0,0)*psi(nx,1,0)+
     >        c(nx,0,0)*psi(nx,0,1))/
     >        (a(nx,0,0)+b(nx,0,0)+c(nx,0,0))+
     >        (1.-omega)*psi(nx,0,0)
c
      psi(0,ny,0)=omega*(-rhoref(0)*(dxyz*pv(0,ny,0)+dyz*va(ny,0)+
     >        dxz*ub(0,0)+dxy*g*coriol(0,ny)*thetabot(0,ny)/
     >        (nsq(0)*thetaref(0)))+
     >        a(0,ny,0)*psi(1,ny,0)+
     >        b(0,ny-1,0)*psi(0,ny-1,0)+
     >        c(0,ny,0)*psi(0,ny,1))/
     >        (a(0,ny,0)+b(0,ny-1,0)+c(0,ny,0))+
     >        (1.-omega)*psi(0,ny,0)
c
      psi(nx,ny,0)=omega*(-rhoref(0)*(dxyz*pv(nx,ny,0)-dyz*vb(ny,0)+
     >        dxz*ub(nx,0)+dxy*g*coriol(nx,ny)*thetabot(nx,ny)/
     >        (nsq(0)*thetaref(0)))+
     >        a(nx,ny,0)*psi(nx-1,ny,0)+
     >        b(nx,ny-1,0)*psi(nx,ny-1,0)+
     >        c(nx,ny,0)*psi(nx,ny,1))/
     >        (a(nx,ny,0)+b(nx,ny-1,0)+c(nx,ny,0))+
     >        (1.-omega)*psi(nx,ny,0)
c
      psi(0,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,0,nz)+dyz*va(0,nz)-
     >        dxz*ua(0,nz)-dxy*g*coriol(0,0)*thetatop(0,0)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(0,0,nz)*psi(1,0,nz)+
     >        b(0,0,nz)*psi(0,1,nz)+
     >        c(0,0,nz-1)*psi(0,0,nz-1))/
     >       (a(0,0,nz)+b(0,0,nz)+c(0,0,nz-1))+
     >       (1.-omega)*psi(0,0,nz)
c
      psi(nx,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,0,nz)-dyz*vb(0,nz)
     >        -dxz*ua(nx,nz)-dxy*g*coriol(nx,0)*thetatop(nx,0)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(nx,0,nz)*psi(nx-1,0,nz)+
     >        b(nx,0,nz)*psi(nx,1,nz)+
     >        c(nx,0,nz-1)*psi(nx,0,nz-1))/
     >        (a(nx,0,nz)+b(nx,0,nz)+c(nx,0,nz-1))+
     >        (1.-omega)*psi(nx,0,nz)
c
      psi(0,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,ny,nz)+
     >        dyz*va(ny,nz)+dxz*ub(0,nz)-
     >        dxy*g*coriol(0,ny)*thetatop(0,ny)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(0,ny,nz)*psi(1,ny,nz)+
     >        b(0,ny-1,nz)*psi(0,ny-1,nz)+
     >        c(0,ny,nz-1)*psi(0,ny,nz-1))/
     >        (a(0,ny,nz)+b(0,ny-1,nz)+c(0,ny,nz-1))+
     >        (1.-omega)*psi(0,ny,nz)
c
      psi(nx,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,ny,nz)-
     >        dyz*vb(ny,nz)+dxz*ub(nx,nz)-
     >        dxy*g*coriol(nx,ny)*thetatop(nx,ny)/
     >        (nsq(2*nz)*thetaref(2*nz)))+
     >        a(nx,ny,nz)*psi(nx-1,ny,nz)+
     >        b(nx,ny-1,nz)*psi(nx,ny-1,nz)+
     >        c(nx,ny,nz-1)*psi(nx,ny,nz-1))/
     >        (a(nx,ny,nz)+b(nx,ny-1,nz)+c(nx,ny,nz-1))+
     >        (1.-omega)*psi(nx,ny,nz)
c

      end

c     --------------------------------------------------------------------------------
c     Init matrix elements for the inversion
c     --------------------------------------------------------------------------------

      subroutine matrixel(a,b,c,coriol,nx,ny,nz,nsq,rhoref,dx,dy,dz)

c     Define the coefficients for the inversion problem (see page 119ff in Rene's 
c     dissertation).

      implicit none

c     Declaration of subroutine parameters
      integer   nx,nz,ny
      real      a     (0:nx,0:ny,0:nz)
      real      b     (0:nx,0:ny,0:nz)
      real      c     (0:nx,0:ny,0:nz)
      real      coriol(0:nx,0:ny)
      real      nsq   (0:2*nz)
      real      rhoref(0:2*nz)
      real      dx,dy,dz

c     Auxiliary variables
      integer   i,j,k

c     Calculate coefficients
      do i=0,nx
         do j=0,ny
            do k=0,nz
               a(i,j,k)=dy*dz*rhoref(2*k)/(dx*coriol(i,j))
               if (j.lt.ny) then
                  b(i,j,k)=dx*dz*rhoref(2*k)/
     >                     (dy*0.5*(coriol(i,j)+coriol(i,j+1)))
               else
                  b(i,j,k)=0.
               endif
               if (k.lt.nz) then
                  c(i,j,k)=dx*dy*rhoref(2*k+1)*coriol(i,j)/
     >                     (dz*nsq(2*k+1))
               else
                  c(i,j,k)=0.
               endif
            enddo
         enddo
      enddo

      end