Subversion Repositories pvinversion.ecmwf

Rev

Blame | Last modification | View Log | Download | RSS feed 

      PROGRAM pv_to_qgpv

c     ********************************************************************************
c     * TRANSFORM ERTEL'S PV TO QUASI-GEOSTROPHIC PV                                 *
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     Input and output file
      character*80   pvsrcfile
      character*80   referfile
      character*80   anomafile

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

c     Numerical and physical parameters
      real           pi180                           ! Pi/180
      parameter      (pi180=3.141592654/180.)
      real           rerd                            ! Earth's radius
      parameter      (rerd=6.371229e6)
      real           eps                             ! Numerical epsilon
      parameter      (eps=0.01)
      real           scale                           ! Scale for PV unit
      parameter      (scale=1e6)
      real           minagl                          ! No PV and qgPV below this height AGL
      parameter      (minagl=1000.)

c     Reference state and grid parameters
      real,   allocatable,dimension (:)   :: nsqref
      real,   allocatable,dimension (:)   :: thetaref
      real,   allocatable,dimension (:)   :: rhoref
      real,   allocatable,dimension (:)   :: pressref
      real,   allocatable,dimension (:)   :: zref
      real,   allocatable,dimension (:,:) :: coriol
      real,   allocatable,dimension (:,:) :: oro
      real    deltax,deltay,deltaz

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

c     Auxiliary variables
      real         zpos
      integer      i,j,k
      integer      stat
      character*80 varname
      integer      istep
      real         mean,rmsq,min,max
      integer      step
      real,allocatable,dimension (:,:) :: tmp
      
c     --------------------------------------------------------------------------------
c     Input
c     --------------------------------------------------------------------------------

      print*,'********************************************************'
      print*,'* PV_TO_QGPV                                           *'
      print*,'********************************************************'


c     Read parameter file
      open(10,file='fort.10')
       read(10,*) pvsrcfile
       read(10,*) referfile
       read(10,*) anomafile
      close(10) 
      print*
      print*,trim(pvsrcfile)
      print*,trim(referfile)
      print*,trim(anomafile)
      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 reference profile and grid parameters
      allocate(rhoref  (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating rhoref'
      allocate(pressref(0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating pressref'
      allocate(thetaref(0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating thetaref'
      allocate(nsqref  (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating nsqref'
      allocate(zref    (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating zref'
      allocate(coriol  (0:nx,0:ny),STAT=stat)
      if (stat.ne.0) print*,'error allocating coriol'
      allocate(oro     (0:nx,0:ny),STAT=stat)
      if (stat.ne.0) print*,'error allocating oro'
      
         
c     Allocate memory for calculation of qgPV and Ertel's 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(pv  (0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating pv'      
      allocate(qgpv(0:nx,0:ny,0:nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating qgpv'

c     Allocate memory for temporary array
      allocate(tmp(0:nx,0:ny),STAT=stat)
      if (stat.ne.0) print*,'error allocating tmp'

c     --------------------------------------------------------------------------------
c     Calculate the qgPV from Ertel's PV and put it onto file
c     --------------------------------------------------------------------------------

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_AIM'
      call read_inp (pv2,varname,pvsrcfile,
     >               nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv)
      
c     Read reference profile and grid parameters
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
     >               nx,ny,nz,deltax,deltay,deltaz,coriol,oro,
     >               referfile)

c     If the PV is negative, set it to zero
      do i=0,nx
         do j=0,ny
            do k=0,nz
               if (pv1(i,j,k).lt.0.) pv1(i,j,k)=0.
               if (pv2(i,j,k).lt.0.) pv2(i,j,k)=0.
            enddo
         enddo
      enddo

c     Get the difference of Ertel's PV and set all missing values to 0
      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
                  pv(i,j,k)=pv1(i,j,k)-pv2(i,j,k)
               else
                  pv(i,j,k)=0.
               endif
            enddo
         enddo
      enddo

c     Calculate qgPV
      call epv_to_qgpv (qgpv,pv,
     >                  rhoref,pressref,nsqref,thetaref,
     >                  nx,ny,nz,mdv)

c     Set values on the boundaries to zero
      do i=0,nx
         do j=0,ny
            qgpv(i,j, 0)=0.
            qgpv(i,j,nz)=0.
         enddo
      enddo
      do i=0,nx
         do k=0,nz
            qgpv(i, 0,k)=0.
            qgpv(i,ny,k)=0.
         enddo
      enddo
      do j=0,ny
         do k=0,nz
            qgpv( 0,j,k)=0.
            qgpv(nx,j,k)=0.
         enddo
      enddo

c     Set all values to zero which are too near to the surface
      do i=0,nx
         do j=0,ny
            do k=0,nz
               zpos=zmin+real(k)*dz
               if (zpos.lt.(oro(i,j)+minagl)) then
                  pv(i,j,k)=0.
                  qgpv(i,j,k)=0.
               endif
            enddo
         enddo
      enddo

c     Write result to netcdf file
      varname='QGPV'
      call write_inp (qgpv,varname,anomafile,nx,ny,nz)
      varname='PV'
      call write_inp (pv,varname,anomafile,nx,ny,nz)

c     Write some info
      print*
      print*,'PV -> qgPV:     k    z     min     max    mean      rmsq'
      step=nz/10
      if (step.lt.1) step=1
      do k=0,nz,step
         do i=0,nx
            do j=0,ny
               tmp(i,j)=qgpv(i,j,k)
            enddo
         enddo
         call calc_error(min,max,mean,rmsq,tmp,nx,ny)
         write(*,'(8x,i3,f10.1,4f10.2)')
     >             k,zmin+real(k)*dz,scale*min,scale*max,
     >                               scale*mean,scale*rmsq
      enddo
      print*         

c     --------------------------------------------------------------------------------
c     Format specifications
c     --------------------------------------------------------------------------------

 111  format (5f20.9)
 106  format (2f20.9)
          
      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     Read refernece profile from netcdf
c     --------------------------------------------------------------------------------

      SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref,
     >                     nx,ny,nz,deltax,deltay,deltaz,coriol,oro,
     >                     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         oro                  : Height of orography (m)
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)
      real                 oro     (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='ORO'
      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,oro,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)-x(i-1))
            count=count+1.
         enddo
      enddo
      deltax=mean/count

      mean=0.
      count=0.
      do j=1,ny
         do i=0,nx
            mean=mean+abs(y(j)-y(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     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


c     ********************************************************************************
c     * CALCULATION                                                                  *
c     ********************************************************************************

c     --------------------------------------------------------------------------------
c     Calculate qgPV from Ertels's PV 
c     --------------------------------------------------------------------------------               

      subroutine epv_to_qgpv (qgpv,pv,
     >                        rhoref,pressref,nsqref,thetaref,
     >                        nx,ny,nz,mdv)

c     Calculate the qgPV from Ertel's PV according to equation 2.11 p16, Thesis 
c     from Rene Fehlmann.

      implicit none
      
c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      qgpv(0:nx,0:ny,0:nz),pv(0:nx,0:ny,0:nz)
      real      rhoref  (0:2*nz)
      real      nsqref  (0:2*nz)
      real      thetaref(0:2*nz)
      real      pressref(0:2*nz)
      real      mdv

c     Numerical epsilon
      real       g
      parameter  (g=9.81)
      real       eps
      parameter  (eps=0.01)
      real       scale
      parameter  (scale=1e6)

c     Auxiliary variables
      integer i,j,k
      integer kk

c     Calculation
      do i=0,nx
         do j=0,ny
            do k=0,nz

               kk=2*k
               
               if (( abs(rhoref(kk)  -mdv).gt.eps).and. 
     >             ( abs(thetaref(kk)-mdv).gt.eps).and. 
     >             ( abs(nsqref(kk)  -mdv).gt.eps).and. 
     >             ( abs(pv(i,j,k)   -mdv).gt.eps)) then

                  qgpv(i,j,k)=rhoref(kk)*g*pv(i,j,k)/thetaref(kk)/ 
     >                        nsqref(kk)/scale
                     
               else
                  qgpv(i,j,k)=0.
               endif

            enddo
         enddo
      enddo

      end

c     --------------------------------------------------------------------------------
c     Calculate error statistics
c     --------------------------------------------------------------------------------               

      subroutine calc_error (min,max,mean,rmsq,tmp,nx,ny)

c     Calculate the error statistics for the two-dimensional error field <tmp>. The
c     following error measures are calculated: the minimum <min>, the maximum <max>,
c     the mean <mean>, the root-mean square <rmsq>

      implicit none

c     Declaration of subroutine parameters
      integer nx,ny
      real    tmp(0:nx,0:ny)
      real    mean,rmsq
      real    min,max

c     Auxiliary variables
      integer i,j
      real    sum
      integer cnt

c     Calculate the minimum and maximum
      min=tmp(0,0)
      max=tmp(0,0)
      do i=0,nx
         do j=0,ny
            if (tmp(i,j).lt.min) min=tmp(i,j)
            if (tmp(i,j).gt.max) max=tmp(i,j)
         enddo
      enddo

c     Calculate the mean
      sum=0.
      cnt=0
      do i=0,nx
         do j=0,ny
            cnt=cnt+1
            sum=sum+tmp(i,j)
         enddo
      enddo
      if (cnt.ge.1) then
         mean=sum/real(cnt)
      else
         mean=0.
      endif

c     Calculate rmsq
      rmsq=0.
      cnt=0
      do i=0,nx
         do j=0,ny
            cnt=cnt+1
            rmsq=rmsq+(tmp(i,j)-mean)**2
         enddo
      enddo
      if (cnt.ge.1) then
         rmsq=1./real(cnt)*sqrt(rmsq)
      else
         rmsq=0.
      endif

      end