Subversion Repositories pvinversion.ecmwf

      PROGRAM z2s

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

      implicit none

c     ---------------------------------------------------------------
c     Declaration of parameters
c     ---------------------------------------------------------------

      real      tzero
      parameter (tzero=273.15)
      real      kappa
      parameter (kappa=0.6078)
      real      rd
      parameter (rd=287.)

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,oro
      real,allocatable, dimension (:,:,:)   :: out1,out2
      real,allocatable, dimension (:,:,:)   :: inp
      real,allocatable,dimension (:)        :: nsqref
      real,allocatable,dimension (:)        :: thetaref
      real,allocatable,dimension (:)        :: tref
      real,allocatable,dimension (:)        :: rhoref
      real,allocatable,dimension (:)        :: pressref
      real,allocatable,dimension (:)        :: zref
      real         deltax,deltay,deltaz
      integer      nvars
      character*80 vnam(100),varname
      integer      isok
      integer      cdfid,cstid
      character*3  mode

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 (:,:)   :: tmp
      character*80 vnam2(100)
      integer      nvars2

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

      print*,'*********************************************************'
      print*,'* qvec_analysis                                         *'
      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     Decide whether to enter ANO or MOD/ORG mode
      mode = ofn(1:3)
      if ( (mode.ne.'ANO').and.
     >     (mode.ne.'MOD').and.
     >     (mode.ne.'ORG') )
     >then
         print*,'Unsupported mode ',trim(mode)
         stop
      endif

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.'GEO' ).and.
     >     (fieldname.ne.'AGEO' ).and.
     >     (fieldname.ne.'DIV_UV' ).and.
     >     (fieldname.ne.'QVEC' ).and.
     >     (fieldname.ne.'DIV_Q' ) ) then
         print*,'Variable not supported ',trim(fieldname)
         stop
      endif
      
c     Set dependencies
      if (fieldname.eq.'GEO') then
         nfields=2
         field_nam(1)='T'
         field_nam(2)='P'
      elseif (fieldname.eq.'AGEO') then
         nfields=4
         field_nam(1)='U'
         field_nam(2)='UG'
         field_nam(3)='V'
         field_nam(4)='VG'
      else if (fieldname.eq.'DIV_UV') then
         nfields=2
         field_nam(1)='U'
         field_nam(2)='V'     
      else if (fieldname.eq.'QVEC') then
         nfields=3
         field_nam(1)='UG'
         field_nam(2)='VG'
         field_nam(3)='TH'
      else if (fieldname.eq.'DIV_Q') then
         nfields=2
         field_nam(1)='QX'
         field_nam(2)='QY'     
      endif

c     Allocate memory
      allocate(field_dat(nfields,nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating field_dat'
      allocate(out1(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating out1'
      allocate(out2(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'error allocating out2'
      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'
      allocate(oro(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating oro'

C     Allocate memory for reference profile
      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(tref    (0:2*nz),STAT=stat)
      if (stat.ne.0) print*,'error allocating tref'

c     Allocate auxiliary fields
      allocate(tmp(nx,ny),stat=stat)
      if (stat.ne.0) print*,'error allocating tmp'
      
c     Read reference profile and grid parameters
      call read_ref (nsqref,rhoref,thetaref,pressref,zref,
     >               nx,ny,nz,deltax,deltay,deltaz,f2,oro,gri)

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     Calculate reference profile of temperature
      print*,'C T_ref = TH_ref * ( P_ref/1000)^kappa'
      do i=0,2*nz
         tref(i) = thetaref(i) * ( pressref(i)/1000. ) ** kappa
      enddo

c     Init height levels
      do i=1,nx
         do j=1,ny
            do k=1,nz
               z3(i,j,k)=zref(2*k-1)
            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     Add reference profile if ANO file is provided
      if ( mode.eq.'ANO' ) then

          print*,'C T -> T_ano + T_ref'
          n=0
          do i=1,nfields
              if (field_nam(i).eq.'T') 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) = field_dat(n,i,j,k) + tref(2*k-1)
              enddo
            enddo
           enddo
          endif

          print*,'C TH -> TH_ano + TH_ref'
          n=0
          do i=1,nfields
              if (field_nam(i).eq.'TH') 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) = field_dat(n,i,j,k)+thetaref(2*k-1)
              enddo
            enddo
           enddo
          endif

          print*,'C P -> P_ano + P_ref'
          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) = field_dat(n,i,j,k)+pressref(2*k-1)
              enddo
            enddo
           enddo
          endif

      endif

c     Change units: P (hPa -> Pa), T(K -> degC)
      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

      n=0
      do i=1,nfields
         if (field_nam(i).eq.'T') n=i
      enddo
      if (n.ne.0) then
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  if ( field_dat(n,i,j,1).gt.100. ) then
                    field_dat(n,i,j,k)=field_dat(n,i,j,k) - tzero
                  endif
               enddo
            enddo
         enddo
      endif

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

c     Geostrophic wind (GEO)
      if (fieldname.eq.'GEO') then

         print*,'C ',trim(fieldname)
         field_nam(1)='RHO'
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  if ( mode.eq.'ANO' ) then
                    field_dat(1,i,j,k)=rhoref(2*k-1)
                  else
                    field_dat(1,i,j,k)=
     >               1./rd*field_dat(2,i,j,k)/(field_dat(1,i,j,k)+tzero)
                  endif
                enddo
            enddo
         enddo
         call calc_geo (out1,out2,                           ! (Ug,Vg)
     >                  field_dat(1,:,:,:),                  ! RHO
     >                  field_dat(2,:,:,:),                  ! P
     >                  z3,x2,y2,f2,                         ! Z,X,Y,CORIOL
     >                  nx,ny,nz,mdv)

c     Ageostrophic wind (AGEO)
      elseif (fieldname.eq.'AGEO') then

         print*,'C ',trim(fieldname)
         call calc_ageo (out1,out2,                           ! (Ua,Va)
     >                   field_dat(1,:,:,:),                  ! U
     >                   field_dat(2,:,:,:),                  ! UG
     >                   field_dat(3,:,:,:),                  ! V
     >                   field_dat(4,:,:,:),                  ! VG
     >                   nx,ny,nz,mdv)

c     Divergence of wind (DIV_UV)
      else if (fieldname.eq.'DIV_UV') then

         print*,'C ',trim(fieldname)
         call calc_div_uv (out1,                             ! Div(U,V))
     >                     field_dat(1,:,:,:),               ! U
     >                     field_dat(2,:,:,:),               ! V
     >                     z3,x2,y2,f2,                      ! Z,X,Y,CORIOL
     >                     nx,ny,nz,mdv)

c     Q vector (QVEC)
      else if (fieldname.eq.'QVEC') then

         print*,'C ',trim(fieldname)
         call calc_qvec (out1,out2,                           ! (Qx,Qy)
     >                   field_dat(1,:,:,:),                  ! UG
     >                   field_dat(2,:,:,:),                  ! VG
     >                   field_dat(3,:,:,:),                  ! TH
     >                   z3,x2,y2,f2,                         ! Z,X,Y,CORIOL
     >                   nx,ny,nz,mdv)

c     Divergence of Q vector (DIV_Q)
      else if (fieldname.eq.'DIV_Q') then
         
         print*,'C ',trim(fieldname)
         call calc_div_q  (out1,                             ! Div(U,V))
     >                     field_dat(1,:,:,:),               ! QX
     >                     field_dat(2,:,:,:),               ! QY
     >                     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
                tmp(i,j)=out1(i,j,k)
             enddo
          enddo
          do n=1,nfilt     
             call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
          enddo
          do i=1,nx
             do j=1,ny
                out1(i,j,k)=tmp(i,j)
             enddo
          enddo

          do i=1,nx
             do j=1,ny
                tmp(i,j)=out2(i,j,k)
             enddo
          enddo
          do n=1,nfilt     
             call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0)
          enddo
          do i=1,nx
             do j=1,ny
                out2(i,j,k)=tmp(i,j)
             enddo
          enddo
          
       enddo

c     ----------------------------------------------------------------
c     Save result onto netcdf file
c     ----------------------------------------------------------------

c     Open output file
      call cdfwopn(ofn,cdfid,ierr)
      if (ierr.ne.0) goto 998

c     Save geostrophic wind
      if (fieldname.eq.'GEO') then
         isok=0
         varname='UG'
         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,out1,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)
                  isok=0
         varname='VG'
         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,out2,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)

c     Save ageostrophic wind
      elseif (fieldname.eq.'AGEO') then
         isok=0
         varname='UA'
         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,out1,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)
                  isok=0
         varname='VA'
         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,out2,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)

c     Save divergence of wind field
      else if (fieldname.eq.'DIV_UV') then
         isok=0
         varname='DIV_UV'
         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,out1,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)

c     Save components of Q vector
      else if (fieldname.eq.'QVEC') then
         isok=0
         varname='QX'
         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,out1,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)
                  isok=0
         varname='QY'
         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,out2,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)

c     Save divergence of wind field
      else if (fieldname.eq.'DIV_Q') then
         isok=0
         varname='DIV_Q'
         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,out1,ierr)
         if (ierr.ne.0) goto 998
         print*,'W ',trim(varname),' ',trim(ofn)

      endif

c     Close output file
      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 geostrophic wind components
c     ----------------------------------------------------------------

      subroutine calc_geo (ug,vg,rho,p,
     >                     z3,x2,y2,f2,nx,ny,nz,mdv) 

c     Calculate the geostrophic wind components (ug,vg) if the temperature 
c     (t) and the pressure (p) are given. The grid and the missing data
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 

      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      ug (nx,ny,nz)
      real      vg (nx,ny,nz)
      real      rho(nx,ny,nz)
      real      p  (nx,ny,nz)
      real      z3 (nx,ny,nz)
      real      x2 (nx,ny)
      real      y2 (nx,ny)
      real      f2 (nx,ny)
      real      mdv

c     Physical parameters and numerical constants
      real      eps
      parameter (eps=0.01)
      real       g
      parameter  (g=9.80616)


c     Auxiliray variables
      integer   i,j,k
      real      dpdx(nx,ny,nz)
      real      dpdy(nx,ny,nz)

c     Calculate horizontal derivatives of pressure
      call deriv(dpdx,p,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv(dpdy,p,'y',z3,x2,y2,nx,ny,nz,mdv)

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

               if ((abs(rho (i,j,k)-mdv).gt.eps).and.
     >             (abs(dpdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dpdy(i,j,k)-mdv).gt.eps)) then

                  ug(i,j,k)=-1./(rho(i,j,k)*f2(i,j))*dpdy(i,j,k)
                  vg(i,j,k)= 1./(rho(i,j,k)*f2(i,j))*dpdx(i,j,k)
                  
               endif

            enddo
         enddo
      enddo

      end

c     ----------------------------------------------------------------
c     Calculate ageostrophic wind components
c     ----------------------------------------------------------------

      subroutine calc_ageo (ua,va,u,ug,v,vg,nx,ny,nz,mdv)

c     Calculate the geostrophic wind components (ug,vg) if the temperature
c     (t) and the pressure (p) are given. The grid and the missing data
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv)

      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      ua (nx,ny,nz)
      real      va (nx,ny,nz)
      real      u  (nx,ny,nz)
      real      ug (nx,ny,nz)
      real      v  (nx,ny,nz)
      real      vg (nx,ny,nz)
      real      mdv

c     Parameters
      real      eps
      parameter (eps=0.01)

c     Auxiliray variables
      integer   i,j,k

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

               if ((abs(u (i,j,k)-mdv).gt.eps).and.
     >             (abs(ug(i,j,k)-mdv).gt.eps).and.
     >             (abs(v (i,j,k)-mdv).gt.eps).and.
     >             (abs(vg(i,j,k)-mdv).gt.eps)) then

                  ua(i,j,k)=u(i,j,k) - ug(i,j,k)
                  va(i,j,k)=v(i,j,k) - vg(i,j,k)

               endif

            enddo
         enddo
      enddo

      end

c     ----------------------------------------------------------------
c     Calculate divergence of wind field
c     ----------------------------------------------------------------

      subroutine calc_div_uv (div_uv,u,v,
     >                        z3,x2,y2,f2,nx,ny,nz,mdv) 

c     Calculate the divergence (div_uv) of the horizontal wind field if+
c     the wind components (u,v) are given. The grid and the missing data
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 

      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      div_uv(nx,ny,nz)
      real      u     (nx,ny,nz)
      real      v     (nx,ny,nz)
      real      z3    (nx,ny,nz)
      real      x2    (nx,ny)
      real      y2    (nx,ny)
      real      f2    (nx,ny)
      real      mdv

c     Physical parameters and numerical constants
      real      eps
      parameter (eps=0.01)

c     Auxiliray variables
      integer   i,j,k
      real      rho
      real      dudx(nx,ny,nz)
      real      dvdy(nx,ny,nz)

c     Calculate horizontal derivatives of pressure
      call deriv(dudx,u,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv(dvdy,v,'y',z3,x2,y2,nx,ny,nz,mdv)

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

               if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dvdy(i,j,k)-mdv).gt.eps)) then

                  div_uv(i,j,k)=dudx(i,j,k)+dvdy(i,j,k)
                  
               endif

            enddo
         enddo
      enddo

      end

c     ----------------------------------------------------------------
c     Calculate Q vector
c     ----------------------------------------------------------------

      subroutine calc_qvec (qx3,qy3,
     >                     th3,u3,v3,
     >                     z3,x2,y2,f2,nx,ny,nz,mdv) 

c     Calculate teh Q vector components <qx3> and <qy3>, as well as the divergence
c     of the Q vector <divq3> on the model grid. The grid is specified in the horizontal
c     by <xmin,ymin,dx,dy,nx,ny>. The number of vertical levels is <nz>. The input field
c     are: potential temperature <th3>, horizontal wind <u3> and <v3>, pressure <p3>.
c     The calculation follows the one described in "Weather Analysis, Dusan Djuric"

      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      qx3(nx,ny,nz)
      real      qy3(nx,ny,nz)
      real      th3(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      scale1,scale2
      parameter (scale1=1.E10,scale2=1.E14)
      real      eps
      parameter (eps=0.01)
      real      g
      parameter (g=9.80616)
      real      tzero
      parameter (tzero=273.16)

c     Auxiliary variables
      real      dudx(nx,ny,nz)
      real      dudy(nx,ny,nz)
      real      dvdx(nx,ny,nz)
      real      dvdy(nx,ny,nz)
      real      dtdx(nx,ny,nz)
      real      dtdy(nx,ny,nz)
      integer   i,j,k

c     Needed derivatives
      call deriv (dudx, u3,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv (dvdy, v3,'y',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)

c     Calculate vector components of the Q vector
      do i=1,nx
         do j=1,ny
            do k=1,nz

c              Evaluate Q vector formula with missing data check
               if ((abs(dudx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dudy(i,j,k)-mdv).gt.eps).and.
     >             (abs(dvdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dvdy(i,j,k)-mdv).gt.eps).and.
     >             (abs(dtdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dtdy(i,j,k)-mdv).gt.eps)) then

                  qx3(i,j,k) = -g/tzero * (dudx(i,j,k)*dtdx(i,j,k)
     >                                     +dvdx(i,j,k)*dtdy(i,j,k))
                  qy3(i,j,k) = -g/tzero * (dudy(i,j,k)*dtdx(i,j,k)
     >                                     +dvdy(i,j,k)*dtdy(i,j,k))

               else
                  qx3(i,j,k)=mdv
                  qy3(i,j,k)=mdv
               endif

            enddo
         enddo
      enddo

c     Scale the output
      do i=1,nx
         do j=1,ny
            do k=1,nz
               if (abs(qx3(i,j,k)-mdv).gt.eps) then
                  qx3(i,j,k)=scale1*qx3(i,j,k)
               endif
               if (abs(qy3(i,j,k)-mdv).gt.eps) then
                  qy3(i,j,k)=scale1*qy3(i,j,k)
               endif
            enddo
         enddo
      enddo

      end

c     ----------------------------------------------------------------
c     Calculate divergence of wind field
c     ----------------------------------------------------------------

      subroutine calc_div_q (div_q,qx,qy,
     >                       z3,x2,y2,f2,nx,ny,nz,mdv) 

c     Calculate the divergence (div_q) of the Q vector field if
c     the components (qx,qy) are given. The grid and the missing data
c     value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) 

      implicit none

c     Declaration of subroutine parameters
      integer   nx,ny,nz
      real      div_q(nx,ny,nz)
      real      qx   (nx,ny,nz)
      real      qy   (nx,ny,nz)
      real      z3   (nx,ny,nz)
      real      x2   (nx,ny)
      real      y2   (nx,ny)
      real      f2   (nx,ny)
      real      mdv

c     Physical parameters and numerical constants
      real      eps
      parameter (eps=0.01)

c     Auxiliray variables
      integer   i,j,k
      real      rho
      real      dqxdx(nx,ny,nz)
      real      dqydy(nx,ny,nz)

c     Calculate horizontal derivatives of pressure
      call deriv(dqxdx,qx,'x',z3,x2,y2,nx,ny,nz,mdv)
      call deriv(dqydy,qy,'y',z3,x2,y2,nx,ny,nz,mdv)

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

               if ((abs(dqxdx(i,j,k)-mdv).gt.eps).and.
     >             (abs(dqydy(i,j,k)-mdv).gt.eps)) then

                  div_q(i,j,k)=dqxdx(i,j,k)+dqydy(i,j,k)
                  
               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       deltay
      parameter  (deltay=111.1775E3)
      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

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