Subversion Repositories pvinversion.ecmwf

      PROGRAM z2s

c     Calculate secondary fields on z levels
c     Michael Sprenger / Summer 2006

      implicit none

c     ---------------------------------------------------------------
c     Declaration of variables
c     ---------------------------------------------------------------

c     Variables for output Z file  : height level
      character*80 cfn
      real         varmin(4),varmax(4),stag(4)
      integer      vardim(4)
      real         mdv
      integer      ndim
      integer      nx,ny,nz
      real         xmin,xmax,ymin,ymax,dx,dy
      integer      ntimes
      real         aklev(1000),bklev(1000)
      real         aklay(1000),bklay(1000)
      real         time
      real         pollon,pollat
      integer      idate(5)
      integer      nfields
      character*80 field_nam(100)
      real,allocatable, dimension (:,:,:,:) :: field_dat
      real,allocatable, dimension (:,:,:)   :: z3
      real,allocatable, dimension (:,:)     :: x2,y2,f2
      real,allocatable, dimension (:,:,:)   :: out
      real,allocatable, dimension (:,:,:)   :: inp
      integer      nvars
      character*80 vnam(100),varname
      integer      isok
      integer      cdfid,cstid

c     Parameter file
      character*80 fieldname
      integer      nfilt
      character*80 ofn,gri

c     Auxiliary variables
      integer      ierr,stat
      integer      i,j,k,n
      real,allocatable, dimension (:,:)   :: out2
      character*80 vnam2(100)
      integer      nvars2

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

      print*,'*********************************************************'
      print*,'* z2s                                                   *'
      print*,'*********************************************************'

c     Read parameter file
      open(10,file='fort.10')
       read(10,*) fieldname
       read(10,*) ofn
       read(10,*) gri
       read(10,*) nfilt
      close(10)

c     Get grid description for Z file : height level
      call cdfopn(ofn,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getcfn(cdfid,cfn,ierr)
      if (ierr.ne.0) goto 998
      call cdfopn(cfn,cstid,ierr)
      if (ierr.ne.0) goto 998
      call getdef(cdfid,'T',ndim,mdv,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)
      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=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy
      call gettimes(cdfid,time,ntimes,ierr)
      if (ierr.ne.0) goto 998
      call getstart(cstid,idate,ierr)
      if (ierr.ne.0) goto 998
      call getpole(cstid,pollon,pollat,ierr)
      if (ierr.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cstid,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998

c     Get a list of all variables on the GRID file
      if (trim(gri).ne.trim(ofn)) then
         call cdfopn(gri,cdfid,ierr)
         if (ierr.ne.0) goto 998
         call getvars(cdfid,nvars2,vnam2,ierr)
         if (ierr.ne.0) goto 998
         do i=1,nvars2
            vnam(nvars+i)=vnam2(i)
         enddo
         nvars=nvars+nvars2
      endif

c     Check whether calculation of <fieldname> is supported
      if ( (fieldname.ne.'TH' ).and.
     >     (fieldname.ne.'NSQ')  .and.
     >     (fieldname.ne.'PV' )   .and.
     >     (fieldname.ne.'RHO')) then
         print*,'Variable not supported ',trim(fieldname)
         stop
      endif
      
c     Set dependencies
      if (fieldname.eq.'TH') then
         nfields=2
         field_nam(1)='T'
         field_nam(2)='P'
      else if (fieldname.eq.'RHO') then
         nfields=2
         field_nam(1)='T'
         field_nam(2)='P'
      else if (fieldname.eq.'NSQ') then
         nfields=2
         field_nam(1)='T'
         field_nam(2)='Q'
      else if (fieldname.eq.'PV') then
         nfields=4
         field_nam(1)='T'
         field_nam(2)='P'
         field_nam(3)='U'
         field_nam(4)='V'
      endif

c     Allocate memory
      allocate(field_dat(nfields,nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating field_dat'
      allocate(out(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating out'
      allocate(inp(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating inp'
      allocate(z3(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating z3'
      allocate(x2(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating x2'
      allocate(y2(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating y2'
      allocate(f2(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating f2'

c     Allocate auxiliary fields
      allocate(out2(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating out2'
      
c     Read X grid
      varname='X'
      isok=0
      call check_varok (isok,varname,vnam,nvars)
      if (isok.eq.0) then
         print*,'Variable ',trim(varname),' missing... Stop'
         stop
      endif
      call cdfopn(gri,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,time,0,x2,ierr)
      if (ierr.ne.0) goto 998
      print*,'R ',trim(varname),' ',trim(gri)
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998

c     Read Y grid
      varname='Y'
      isok=0
      call check_varok (isok,varname,vnam,nvars)
      if (isok.eq.0) then
         print*,'Variable ',trim(varname),' missing... Stop'
         stop
      endif
      call cdfopn(gri,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,time,0,y2,ierr)
      if (ierr.ne.0) goto 998
      print*,'R ',trim(varname),' ',trim(gri)
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998
      
c     Read Coriolis parametzer
      varname='CORIOL'
      isok=0
      call check_varok (isok,varname,vnam,nvars)
      if (isok.eq.0) then
         print*,'Variable ',trim(varname),' missing... Stop'
         stop
      endif
      call cdfopn(gri,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,time,0,f2,ierr)
      if (ierr.ne.0) goto 998
      print*,'R ',trim(varname),' ',trim(gri)
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998      

c     Init height levels
      do i=1,nx
         do j=1,ny
            do k=1,nz
               z3(i,j,k)=aklay(k)
            enddo
         enddo
      enddo

c     Load needed variables
      do n=1,nfields

c        Check whether variable is available on file
         varname=field_nam(n)
         isok=0
         call check_varok (isok,varname,vnam,nvars)

c        Variable is available on file
         if (isok.eq.1) then

            call cdfopn(ofn,cdfid,ierr)
            if (ierr.ne.0) goto 998
            call getdat(cdfid,varname,time,0,inp,ierr)
            if (ierr.ne.0) goto 998
            print*,'R ',trim(varname),' ',trim(ofn)
            call clscdf(cdfid,ierr)
            if (ierr.ne.0) goto 998
            
            do i=1,nx
               do j=1,ny
                  do k=1,nz
                     field_dat(n,i,j,k)=inp(i,j,k)
                  enddo
               enddo
            enddo
            
         else            
            print*,'Variable missing : ',trim(varname)
            stop
         endif

      enddo

c     Change unit of pressure from hPa to Pa
      n=0
      do i=1,nfields
         if (field_nam(i).eq.'P') n=i
      enddo
      if (n.ne.0) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  field_dat(n,i,j,k)=100.*field_dat(n,i,j,k)
               enddo
            enddo
         enddo
      endif

c     ----------------------------------------------------------------
c     Calculations
c     ----------------------------------------------------------------

c     Call to the defining routines
      if (fieldname.eq.'RHO') then

         print*,'C ',trim(fieldname)
         call calc_rho (out,                                 ! RHO
     >                  field_dat(1,:,:,:),                  ! T
     >                  field_dat(2,:,:,:),                  ! P
     >                  nx,ny,nz,mdv)

      else if (fieldname.eq.'TH') then

         print*,'C ',trim(fieldname)
         call calc_theta (out,                               ! TH
     >                    field_dat(1,:,:,:),                ! T
     >                    field_dat(2,:,:,:),                ! P
     >                    nx,ny,nz,mdv)

       else if (fieldname.eq.'NSQ') then

          print*,'C ',trim(fieldname)
          call calc_nsq (out,                                 ! NSQ
     >                   field_dat(1,:,:,:),                  ! T
     >                   field_dat(2,:,:,:),                  ! Q
     >                   z3,                                  ! Z      
     >                   nx,ny,nz,mdv)
      
       else if (fieldname.eq.'PV') then

          print*,'C ',trim(fieldname)
          call calc_pv (out,                                  ! PV
     >                  field_dat(1,:,:,:),                   ! T
     >                  field_dat(2,:,:,:),                   ! P
     >                  field_dat(3,:,:,:),                   ! U
     >                  field_dat(4,:,:,:),                   ! V
     >                  z3,x2,y2,f2,                          ! Z,X,Y,CORIOL
     >                  nx,ny,nz,mdv)

       endif

c      Horizontal filtering of the resulting fields
       print*,'F ',trim(fieldname)
       do k=1,nz
          
          do i=1,nx
             do j=1,ny
                out2(i,j)=out(i,j,k)
             enddo
          enddo
          do n=1,nfilt     
             call filt2d (out2,out2,nx,ny,1.,mdv,0,0,0,0)
          enddo

          do i=1,nx
             do j=1,ny
                out(i,j,k)=out2(i,j)
             enddo
          enddo
          
       enddo

c     ----------------------------------------------------------------
c     Save result onto netcdf file
c     ----------------------------------------------------------------
       
      call cdfwopn(ofn,cdfid,ierr)
      if (ierr.ne.0) goto 998
      isok=0
      varname=fieldname
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,out,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(ofn)
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998


c     ----------------------------------------------------------------
c     Exception handling
c     ----------------------------------------------------------------
   
      stop

 998  print*,'Problem with netcdf file'
      stop

      end

       

c     ****************************************************************
c     * SUBROUTINE SECTION: AUXILIARY ROUTINES                       *
c     ****************************************************************     

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)>. 
c     If this is the case, <isok> is incremented by 1. Otherwise <isok> 
c     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     * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS               *
c     ****************************************************************

c     -----------------------------------------------------------------
c     Calculate density
c     -----------------------------------------------------------------

      subroutine calc_rho (rho3,t3,p3,
     >                     nx,ny,nz,mdv)

c     Calculate the density <rho3> (in kg/m^3) from temperature <t3> 
c     (in deg C) and pressure <p3> (in hPa).

      implicit none
      
c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      mdv
      real      rho3(nx,ny,nz)
      real      t3  (nx,ny,nz)
      real      p3  (nx,ny,nz)

c     Declaration of physical constants
      real      eps 
      parameter (eps=0.01)
      real      rd
      parameter (rd=287.)
      real      tzero
      parameter (tzero=273.15)

c     Auxiliary variables
      integer   i,j,k
      
c     Calculation
      do i=1,nx
         do j=1,ny
            do k=1,nz
            
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
     >             (abs(p3(i,j,k)-mdv).gt.eps)) then
               
                  rho3(i,j,k)=1./rd*p3(i,j,k)/(t3(i,j,k)+tzero)

               else
                  rho3(i,j,k)=mdv
               endif
               
            enddo
         enddo
      enddo

      end


c     -----------------------------------------------------------------
c     Calculate potential temperature
c     -----------------------------------------------------------------

      subroutine calc_theta (th3,t3,p3,
     >                       nx,ny,nz,mdv)

c     Calculate potential temperature <th3> (in K) from temperature <t3> 
c     (in deg C) and pressure <p3> (in hPa). 
    
      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      mdv
      real      th3(nx,ny,nz)
      real      t3 (nx,ny,nz)
      real      p3 (nx,ny,nz)
      
c     Declaration of physical constants
      real      rdcp,tzero,p0
      parameter (rdcp=0.286)
      parameter (tzero=273.15)
      parameter (p0=100000.)
      real      eps
      parameter (eps=0.01)
      
c     Auxiliary variables
      integer  i,j,k

c     Calculation
      do i=1,nx
         do j=1,ny
            do k=1,nz
            
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
     >             (abs(p3(i,j,k)-mdv).gt.eps)) then
               
                  th3(i,j,k)=(t3(i,j,k)+tzero)*( (p0/p3(i,j,k))**rdcp )

               else
                  th3(i,j,k)=mdv
               endif
               
            enddo
         enddo
      enddo

      end

c     -----------------------------------------------------------------
c     Calculate stratification
c     -----------------------------------------------------------------

      subroutine calc_nsq (nsq3,t3,q3,
     >                     z3,nx,ny,nz,mdv)

c     Calculate stratification <nsq3> on the model grid. The grid is
c     specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number
c     of vertical levels is <nz>. The input field are: temperature <t3>, 
c     specific humidity <q3>, height <z3> and horizontal grid <x2>, <y2>. 

      implicit none
      
c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      nsq3(nx,ny,nz)
      real      t3  (nx,ny,nz)
      real      q3  (nx,ny,nz)
      real      z3  (nx,ny,nz)
      real      x2  (nx,ny)
      real      y2  (nx,ny)
      real      mdv

c     Physical and numerical parameters
      real      scale
      parameter (scale=1.)
      real      g
      parameter (g=9.80665)
      real      eps
      parameter (eps=0.01)
      real      tzero
      parameter (tzero=273.15)
      real      kappa
      parameter (kappa=0.6078)
      real      cp
      parameter (cp=1005.)
      real      zerodiv
      parameter (zerodiv=0.0000000001)
      real      R
      parameter (R=287.)

c     Auxiliary variables
      real      tv   (nx,ny,nz)
      real      dtvdz(nx,ny,nz)
      integer   i,j,k
      real      scaledz

c     Calculate 3d virtual temperature
      do k=1,nz
         do i=1,nx
            do j=1,ny
               if ((abs(t3(i,j,k)-mdv).gt.eps).and.
     >             (abs(q3(i,j,k)-mdv).gt.eps)) 
     >         then
                  tv(i,j,k) = (t3(i,j,k)+tzero)*(1.+kappa*q3(i,j,k))
               else
                  tv(i,j,k) = mdv
               endif
            enddo
         enddo
      enddo

c     Vertical derivative of virtual temperature
      call  deriv (dtvdz,tv,'z',z3,x2,y2,nx,ny,nz,mdv)

c     Stratification
      do i=1,nx
         do j=1,ny
            do k=1,nz
               if ((abs(dtvdz(i,j,k)-mdv).gt.eps).and.
     >             (abs(tv   (i,j,k)-mdv).gt.eps)) 
     >         then
                  nsq3(i,j,k) = g/tv(i,j,k) * (dtvdz(i,j,k) + g/cp)  
               else
                  nsq3(i,j,k) = mdv
               endif
           enddo
        enddo
      enddo

      end
          
c     -----------------------------------------------------------------
c     Calculate potential vorticity
c     -----------------------------------------------------------------

      subroutine calc_pv (pv3,t3,p3,u3,v3,
     >                    z3,x2,y2,f2,nx,ny,nz,mdv)

c     Calculate potential vorticity <pv3> on the model grid. The grid is
c     specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number
c     of vertical levels is <nz>. The input field are: potential temperature
c     <th3>, horizontal wind <u3> and <v3>, density <rho3>, height <z3>. 
c     The terms including the vertical velocity are neglected in the calculation 
c     of the PV. 

      implicit none
      
c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      pv3 (nx,ny,nz)
      real      t3  (nx,ny,nz)
      real      p3  (nx,ny,nz)
      real      u3  (nx,ny,nz)
      real      v3  (nx,ny,nz)
      real      z3  (nx,ny,nz)
      real      x2  (nx,ny)
      real      y2  (nx,ny)
      real      f2  (nx,ny)

      real      mdv

c     Physical and numerical parameters
      real      scale
      parameter (scale=1.E6)
      real      omega
      parameter (omega=7.292E-5)
      real      pi180
      parameter (pi180=3.141592654/180.)
      real      eps
      parameter (eps=0.01)
      
c     Auxiliary variables
      real      dtdz(nx,ny,nz)
      real      dtdx(nx,ny,nz)
      real      dtdy(nx,ny,nz)
      real      dudz(nx,ny,nz)
      real      dvdz(nx,ny,nz)
      real      dvdx(nx,ny,nz)
      real      dudy(nx,ny,nz)
      real      rho3(nx,ny,nz)
      real      th3 (nx,ny,nz)

      integer   i,j,k

c     Calculate density and potential temperature
      call calc_rho   (rho3,t3,p3,nx,ny,nz,mdv)
      call calc_theta (th3 ,t3,p3,nx,ny,nz,mdv)

c     Calculate needed derivatives
      call deriv (dudz, u3,'z',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dvdz, v3,'z',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dtdz,th3,'z',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv)     
      call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)     

c     Calculate potential vorticity
      do j=1,ny
         do i=1,nx
            do k=1,nz

c              Evaluate PV formula with missing data check
               if ((abs(dtdz(i,j,k)-mdv).gt.eps).and.
     >             (abs(dudz(i,j,k)-mdv).gt.eps).and.
     >             (abs(dtdy(i,j,k)-mdv).gt.eps).and.
     >             (abs(dvdz(i,j,k)-mdv).gt.eps).and.
     >             (abs(rho3(i,j,k)-mdv).gt.eps).and.              
     >             (abs(dtdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dvdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dudy(i,j,k)-mdv).gt.eps)) then

                  pv3(i,j,k)=scale*1./rho3(i,j,k)*
     >                   ((dvdx(i,j,k)-dudy(i,j,k)+f2(i,j))*dtdz(i,j,k)
     >                                         +dudz(i,j,k)*dtdy(i,j,k)
     >                                         -dvdz(i,j,k)*dtdx(i,j,k))
               else
                  pv3(i,j,k)=mdv
               endif

            enddo
         enddo
      enddo

      end


c     ****************************************************************
c     * SUBROUTINE SECTION: GRID HANDLING                            *
c     ****************************************************************

c     -----------------------------------------------------------------
c     Horizontal and vertical derivatives for 3d fields
c     -----------------------------------------------------------------

      subroutine deriv (df,f,direction,
     >                  z3,x2,y2,nx,ny,nz,mdv)

c     Calculate horizontal and vertical derivatives of the 3d field <f>.
c     The direction of the derivative is specified in <direction> 
c         'x','y'          : Horizontal derivative in x and y direction
c         'p','z','t','m'  : Vertical derivative (pressure, height, theta, model)
c     The 3d field <z3> specifies the isosurfaces along which the horizontal 
c     derivatives are calculated or the levels for the vertical derivatives.

      implicit none

c     Input and output parameters
      integer    nx,ny,nz
      real       df (nx,ny,nz)
      real       f  (nx,ny,nz)
      real       z3 (nx,ny,nz)
      real       x2 (nx,ny)
      real       y2 (nx,ny)
      character  direction
      real       mdv
      
c     Numerical and physical parameters
      real       pi180
      parameter  (pi180=3.141592654/180.)
      real       zerodiv
      parameter  (zerodiv=0.00000001)
      real       eps
      parameter  (eps=0.01) 

c     Auxiliary variables
      integer    i,j,k
      real       vmin,vmax
      real       scale,lat
      real       vu,vl,vuvl,vlvu
      integer    o,u,w,e,n,s

c     Vertical derivative
      if ((direction.eq.'z').or.
     >    (direction.eq.'th').or.
     >    (direction.eq.'p').or.
     >    (direction.eq.'m').and.
     >    (nz.gt.1)) then

         do i=1,nx
            do j=1,ny
               do k=1,nz
                  
                  o=k+1
                  if (o.gt.nz) o=nz
                  u=k-1
                  if (u.lt.1) u=1                  

                  if ((abs(f(i,j,o)-mdv).gt.eps).and.
     >                (abs(f(i,j,u)-mdv).gt.eps).and.
     >                (abs(f(i,j,k)-mdv).gt.eps)) then

                     vu = z3(i,j,k)-z3(i,j,o)
                     vl = z3(i,j,u)-z3(i,j,k)
                     vuvl = vu/(vl+zerodiv)
                     vlvu = 1./(vuvl+zerodiv)

                     df(i,j,k) = 1./(vu+vl)
     >                           * (vuvl*(f(i,j,u)-f(i,j,k)) 
     >                           +  vlvu*(f(i,j,k)-f(i,j,o))) 

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo
              
c     Horizontal derivative in the y direction: 3d
      elseif (direction.eq.'y') then
         
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  
                  s=j-1
                  if (s.lt.1) s=1
                  n=j+1
                  if (n.gt.ny) n=ny

                  if ((abs(f(i,n,k)-mdv).gt.eps).and.
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
     >                (abs(f(i,s,k)-mdv).gt.eps)) then  
 
                     vu = 1000.*(y2(i,j)-y2(i,n))
                     vl = 1000.*(y2(i,s)-y2(i,j))
                     vuvl = vu/(vl+zerodiv)
                     vlvu = 1./(vuvl+zerodiv)
 
                     df(i,j,k) =  1./(vu+vl)
     >                            * (vuvl*(f(i,s,k)-f(i,j,k))
     >                            +  vlvu*(f(i,j,k)-f(i,n,k)))

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo

c     Horizontal derivative in the x direction: 3d
      elseif (direction.eq.'x') then

         do i=1,nx
            do j=1,ny
               do k=1,nz
                  
                  w=i-1
                  if (w.lt.1) w=1
                  e=i+1
                  if (e.gt.nx) e=nx

                  if ((abs(f(w,j,k)-mdv).gt.eps).and.
     >                (abs(f(i,j,k)-mdv).gt.eps).and.
     >                (abs(f(e,j,k)-mdv).gt.eps)) then  
 
                     vu = 1000.*(x2(i,j)-x2(e,j))
                     vl = 1000.*(x2(w,j)-x2(i,j))
                     vuvl = vu/(vl+zerodiv)
                     vlvu = 1./(vuvl+zerodiv)
 
                     df(i,j,k) =  1./(vu+vl)
     >                            * (vuvl*(f(w,j,k)-f(i,j,k))
     >                            +  vlvu*(f(i,j,k)-f(e,j,k)))

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo
            
c     Undefined direction for derivative
      else
         
         print*,'Invalid direction of derivative... Stop'
         stop
      
      endif

      end

c     -----------------------------------------------------------------
c     Horizontal filter
c     -----------------------------------------------------------------

      subroutine filt2d (a,af,nx,ny,fil,misdat,
     &                   iperx,ipery,ispol,inpol)

c     Apply a conservative diffusion operator onto the 2d field a,
c     with full missing data checking.
c
c     a     real   inp  array to be filtered, dimensioned (nx,ny)
c     af    real   out  filtered array, dimensioned (nx,ny), can be
c                       equivalenced with array a in the calling routine
c     f1    real        workarray, dimensioned (nx+1,ny)
c     f2    real        workarray, dimensioned (nx,ny+1)
c     fil   real   inp  filter-coeff., 0<afil<=1. Maximum filtering with afil=1
c                       corresponds to one application of the linear filter.
c     misdat real  inp  missing-data value, a(i,j)=misdat indicates that
c                       the corresponding value is not available. The
c                       misdat-checking can be switched off with with misdat=0.
c     iperx int    inp  periodic boundaries in the x-direction (1=yes,0=no)
c     ipery int    inp  periodic boundaries in the y-direction (1=yes,0=no)
c     inpol int    inp  northpole at j=ny  (1=yes,0=no)
c     ispol int    inp  southpole at j=1   (1=yes,0=no)
c
c     Christoph Schaer, 1993

c     argument declaration
      integer     nx,ny
      real        a(nx,ny),af(nx,ny),fil,misdat
      integer     iperx,ipery,inpol,ispol

c     local variable declaration
      integer     i,j,is
      real        fh
      real        f1(nx+1,ny),f2(nx,ny+1)

c     compute constant fh
      fh=0.125*fil

c     compute fluxes in x-direction
      if (misdat.eq.0.) then
        do j=1,ny
        do i=2,nx
          f1(i,j)=a(i-1,j)-a(i,j)
        enddo
        enddo
      else
        do j=1,ny
        do i=2,nx
          if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then
            f1(i,j)=0.
          else
            f1(i,j)=a(i-1,j)-a(i,j)
          endif
        enddo
        enddo
      endif
      if (iperx.eq.1) then
c       do periodic boundaries in the x-direction
        do j=1,ny
          f1(1,j)=f1(nx,j)
          f1(nx+1,j)=f1(2,j)
        enddo
      else
c       set boundary-fluxes to zero
        do j=1,ny
          f1(1,j)=0.
          f1(nx+1,j)=0.
        enddo
      endif

c     compute fluxes in y-direction
      if (misdat.eq.0.) then
        do j=2,ny
        do i=1,nx
          f2(i,j)=a(i,j-1)-a(i,j)
        enddo
        enddo
      else
        do j=2,ny
        do i=1,nx
          if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then
            f2(i,j)=0.
          else
            f2(i,j)=a(i,j-1)-a(i,j)
          endif
        enddo
        enddo
      endif
c     set boundary-fluxes to zero
      do i=1,nx
        f2(i,1)=0.
        f2(i,ny+1)=0.
      enddo
      if (ipery.eq.1) then
c       do periodic boundaries in the x-direction
        do i=1,nx
          f2(i,1)=f2(i,ny)
          f2(i,ny+1)=f2(i,2)
        enddo
      endif
      if (iperx.eq.1) then
        if (ispol.eq.1) then
c         do south-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,1)=-f2(mod(i-1+is,nx)+1,2)
          enddo
        endif
        if (inpol.eq.1) then
c         do north-pole
          is=(nx-1)/2
          do i=1,nx
            f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny)
          enddo
        endif
      endif

c     compute flux-convergence -> filter
      if (misdat.eq.0.) then
        do j=1,ny
        do i=1,nx
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
        enddo
        enddo
      else
        do j=1,ny
        do i=1,nx
          if (a(i,j).eq.misdat) then
            af(i,j)=misdat
          else
            af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1))
          endif
        enddo
        enddo
      endif
      end