Subversion Repositories pvinversion.ecmwf

      PROGRAM def_anomaly

c     ************************************************************************
c     * Define a PV anomaly and set it up for inversion. Additionally, the   *
c     * lateral and upper/lower boundary values for potential temperature    *
c     * and wind are set                                                     *
c     * Michael Sprenger / Spring 2005                                       *
c     ************************************************************************

      implicit none

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

c     Numerical epsilon
      real         eps
      parameter    (eps=0.01)

c     Variables for input Z file  : height level
      character*80 in_zfn,in_varname
      real         in_varmin(4),in_varmax(4),in_stag(4)
      integer      in_vardim(4)
      real         in_mdv
      integer      in_ndim
      integer      in_nx,in_ny,in_nz
      real         in_xmin,in_xmax,in_ymin,in_ymax,in_dx,in_dy
      integer      in_ntimes
      real         in_aklev(1000),in_bklev(1000)
      real         in_aklay(1000),in_bklay(1000)

      real         in_time
      real         in_pollon,in_pollat
      integer      in_idate(5)
      integer      in_nvars
      character*80 in_vnam(100)
      real,allocatable, dimension (:,:,:)    :: in_field
      real,allocatable, dimension (:,:,:)    :: in_filt
      real,allocatable, dimension (:,:,:)    :: in_anom
      real,allocatable, dimension (:,:)      :: in_x,in_y
      real,allocatable, dimension (:,:,:)    :: th_anom
      real,allocatable, dimension (:,:,:)    :: uu_anom,vv_anom

c     Filtering
      real          f_xmin,f_xmax
      real          f_ymin,f_ymax
      real          f_zmin,f_zmax
      integer       nfilt
      real          boundxy,boundz
      real,allocatable, dimension (:,:)      :: tmp

c     Auxiliary variables
      integer      i,j,k,l,n
      integer      cdfid,cstid
      integer      ierr,stat
      character*80 cfn
      real         xpt,ypt,zpt
      integer      isok
      character*80 varname,cdfname
      character*80 parfile
      real         distxy,distz,distz1,distz2
      integer      indx,indy,indz
      real         sum
      integer      count
      character*80 name

c     Externals
      real         kink
      external     kink
      
c     ------------------------------------------------------------------------
c     Preparations
c     ------------------------------------------------------------------------

      print*,'*********************************************************'
      print*,'* def_anomaly                                           *'
      print*,'*********************************************************'

c     Set name of the PV variable
      in_varname='PV'

c     Read parameters
      open(10,file='fort.10')
       read(10,*) in_zfn
       read(10,*) name,f_xmin
       if ( trim(name).ne.'BOX_XMIN') stop
       read(10,*) name,f_xmax
       if ( trim(name).ne.'BOX_XMAX') stop
       read(10,*) name,f_ymin
       if ( trim(name).ne.'BOX_YMIN') stop
       read(10,*) name,f_ymax
       if ( trim(name).ne.'BOX_YMAX') stop
       read(10,*) name,f_zmin
       if ( trim(name).ne.'BOX_ZMIN') stop
       read(10,*) name,f_zmax
       if ( trim(name).ne.'BOX_ZMAX') stop
       read(10,*) name,nfilt
       if ( trim(name).ne.'NFILTER' ) stop
       read(10,*) name,boundxy
       if ( trim(name).ne.'BOUND_XY') stop
       read(10,*) name,boundz
       if ( trim(name).ne.'BOUND_Z')  stop
      close(10)

      print*,'Filter domain, x: ', f_xmin,f_xmax
      print*,'               y: ', f_ymin,f_ymax
      print*,'               z: ', f_zmin,f_zmax
      print*,'Nfilt:            ',nfilt
      print*,'Boundxy,z:        ',boundxy,boundz
      
c     Get grid description for Z file : height levels
      call cdfopn(in_zfn,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 getvars(cdfid,in_nvars,in_vnam,ierr)
      if (ierr.ne.0) goto 998
      call getdef(cdfid,in_varname,in_ndim,in_mdv,in_vardim,
     >            in_varmin,in_varmax,in_stag,ierr)
      if (ierr.ne.0) goto 998
      in_nx  =in_vardim(1)
      in_ny  =in_vardim(2)
      in_nz  =in_vardim(3)
      in_xmin=in_varmin(1)
      in_ymin=in_varmin(2)
      call getlevs(cstid,in_nz,in_aklev,in_bklev,in_aklay,in_bklay,ierr)
      if (ierr.ne.0) goto 998
      call getgrid(cstid,in_dx,in_dy,ierr)
      if (ierr.ne.0) goto 998
      in_xmax=in_xmin+real(in_nx-1)*in_dx
      in_ymax=in_ymin+real(in_ny-1)*in_dy
      call gettimes(cdfid,in_time,in_ntimes,ierr)
      if (ierr.ne.0) goto 998
      call getstart(cstid,in_idate,ierr)
      if (ierr.ne.0) goto 998
      call getpole(cstid,in_pollon,in_pollat,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cstid,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)

c     Allocate memory for grid description
      allocate(in_x(in_nx,in_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_x ***'    
      allocate(in_y(in_nx,in_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_y ***'    

c     Allocate memory for input and output fields
      allocate(in_field(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_field ***'
      allocate(in_anom(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_anom ***'
      allocate(in_filt(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array in_filt ***'
      allocate(th_anom(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array th_anom ***'
      allocate(uu_anom(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array uu_anom ***'
      allocate(vv_anom(in_nx,in_ny,in_nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array vvnom ***'

c     Allocate memory for filtering
      allocate(tmp(in_nx,in_ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp ***'

c     Read variables from file
      call cdfopn(in_zfn,cdfid,ierr)
      if (ierr.ne.0) goto 998

      varname=trim(in_varname)
      print*,'R ',trim(varname),'   ',trim(in_zfn)
      call gettimes(cdfid,in_time,in_ntimes,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,in_time,0,in_field,ierr)
      if (ierr.ne.0) goto 998

      varname='X'
      print*,'R ',trim(varname),'   ',trim(in_zfn)
      call gettimes(cdfid,in_time,in_ntimes,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,in_time,0,in_x,ierr)
      if (ierr.ne.0) goto 998

      varname='Y'
      print*,'R ',trim(varname),'   ',trim(in_zfn)
      call gettimes(cdfid,in_time,in_ntimes,ierr)
      if (ierr.ne.0) goto 998
      call getdat(cdfid,varname,in_time,0,in_y,ierr)
      if (ierr.ne.0) goto 998

      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998

c     ------------------------------------------------------------------------
c     Get the average along the x axis: "zonal" average
c     ------------------------------------------------------------------------

      print*,'A ',trim(in_varname)
      do k=1,in_nz
         do j=1,in_ny

c           Get the mean along x axis
            sum=0.
            count=0
            do i=1,in_nx
               if (abs(in_field(i,j,k)-in_mdv).gt.eps) then
                  sum=sum+in_field(i,j,k)
                  count=count+1
               endif
            enddo
            if (count.ge.1) then
               sum=sum/real(count)
            else
               sum=in_mdv
            endif

c           Get the difference relative to the "zonal" mean
            do i=1,in_nx
               xpt=in_x(i,j)
               ypt=in_y(i,j)
               zpt=in_aklay(k)
               if ( (xpt.ge.f_xmin).and.
     >              (xpt.le.f_xmax).and.
     >              (ypt.ge.f_ymin).and.
     >              (ypt.le.f_ymax).and.
     >              (zpt.ge.f_zmin).and.
     >              (zpt.le.f_zmax)) then
                  in_filt(i,j,k)=sum
               else
                  in_filt(i,j,k)=in_field(i,j,k)
               endif
            enddo
         
         enddo
      enddo

c     ------------------------------------------------------------------------
c     Get the anomaly field
c     ------------------------------------------------------------------------
c 
c     Determine the difference between filtered and unfiltered field
      do i=1,in_nx
         do j=1,in_ny
            do k=1,in_nz
               if ( (abs(in_field(i,j,k)-in_mdv).gt.eps).and.
     >              (abs(in_filt (i,j,k)-in_mdv).gt.eps)) then
                  in_anom(i,j,k)=in_field(i,j,k)-in_filt(i,j,k)
               else
                  in_anom(i,j,k)=in_mdv
               endif
            enddo
         enddo
      enddo

c     Take only the positive part of the difference (mdv is set to 0)
      do i=1,in_nx
         do j=1,in_ny
            do k=1,in_nz
               if ( (abs(in_anom(i,j,k)-in_mdv).lt.eps) ) then
                  in_anom(i,j,k)=0.
               else if ( in_anom(i,j,k).lt.0.) then
                  in_anom(i,j,k)=0.
               endif
            enddo
         enddo
      enddo

c     Smooth transition to zero at boundaries
      do i=1,in_nx
         do j=1,in_ny
            do k=1,in_nz
               if ( (abs(in_anom (i,j,k)-in_mdv).gt.eps) ) then
                  
                  in_anom(i,j,k)=in_anom(i,j,k)*
     >                            kink(f_zmax     -in_aklay(k),boundz )*
     >                            kink(in_aklay(k)-f_zmin     ,boundz )*
     >                            kink(in_x(i,j)  -f_xmin     ,boundxy)*
     >                            kink(f_xmax     -in_x(i,j)  ,boundxy)*
     >                            kink(in_y(i,j)  -f_ymin     ,boundxy)*
     >                            kink(f_ymax     -in_y(i,j)  ,boundxy)

               endif
            enddo
         enddo
      enddo      

c     ------------------------------------------------------------------------
c     Apply a digital filter to the anomaly field
c     ------------------------------------------------------------------------

      print*,'F ',trim(in_varname)
      do k=1,in_nz
         
         do i=1,in_nx
            do j=1,in_ny
               tmp(i,j)=in_anom(i,j,k)
            enddo
         enddo
         
         do n=1,nfilt     
            call filt2d (tmp,tmp,in_nx,in_ny,1.,in_mdv,0,0,0,0)
         enddo
         
         do i=1,in_nx
            do j=1,in_ny
               in_anom(i,j,k)=tmp(i,j)
            enddo
         enddo
         
      enddo

c     ------------------------------------------------------------------------
c     Get the modified PV field (original minus anomaly)
c     ------------------------------------------------------------------------

      do i=1,in_nx
         do j=1,in_ny
            do k=1,in_nz 
               if ( (abs(in_anom (i,j,k)-in_mdv).gt.eps).and.
     >              (abs(in_field(i,j,k)-in_mdv).gt.eps) ) then
                  in_filt(i,j,k)=in_field(i,j,k)-in_anom(i,j,k)
               else
                  in_filt(i,j,k)=in_mdv
               endif
            enddo
         enddo
      enddo

c     -----------------------------------------------------------------
c     Define the boundary values for wind and potential temperature
c     -----------------------------------------------------------------

      do i=1,in_nx
         do j=1,in_ny
            do k=1,in_nz
               th_anom(i,j,k)=0.
               uu_anom(i,j,k)=0.
               vv_anom(i,j,k)=0.
            enddo
         enddo
      enddo
         
c     -----------------------------------------------------------------
c     Write output
c     -----------------------------------------------------------------

c     Open output netcdf file
      call cdfwopn(in_zfn,cdfid,ierr)
      if (ierr.ne.0) goto 997

c     Write the 3d filtered field
      isok=0
      varname=trim(in_varname)//'_FILT'
      print*,'Write ',trim(varname),'  ',trim(in_zfn)
      call check_varok(isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,4,in_mdv,in_vardim,
     >               in_varmin,in_varmax,in_stag,ierr)
         if (ierr.ne.0) goto 997
      endif
      call putdat(cdfid,varname,in_time,0,in_filt,ierr)
      if (ierr.ne.0) goto 997

c     Write the 3d anomaly field
      isok=0
      varname=trim(in_varname)//'_ANOM'
      print*,'Write ',trim(varname),'  ',trim(in_zfn)
      call check_varok(isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,4,in_mdv,in_vardim,
     >               in_varmin,in_varmax,in_stag,ierr)
         if (ierr.ne.0) goto 997
      endif
      call putdat(cdfid,varname,in_time,0,in_anom,ierr)
      if (ierr.ne.0) goto 997

c     Write lower and upper boundary
      isok=0
      varname='TH_ANOM'
      print*,'Write ',trim(varname),'  ',trim(in_zfn)
      call check_varok(isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,4,in_mdv,in_vardim,
     >               in_varmin,in_varmax,in_stag,ierr)
         if (ierr.ne.0) goto 997
      endif
      call putdat(cdfid,varname,in_time,0,th_anom,ierr)
      if (ierr.ne.0) goto 997

c     Write lateral boundary 
      isok=0
      varname='UU_ANOM'
      print*,'Write ',trim(varname),'  ',trim(in_zfn)
      call check_varok(isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,4,in_mdv,in_vardim,
     >               in_varmin,in_varmax,in_stag,ierr)
         if (ierr.ne.0) goto 997
      endif
      call putdat(cdfid,varname,in_time,0,uu_anom,ierr)
      if (ierr.ne.0) goto 997
      isok=0
      varname='VV_ANOM'
      print*,'Write ',trim(varname),'  ',trim(in_zfn)
      call check_varok(isok,varname,in_vnam,in_nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,4,in_mdv,in_vardim,
     >               in_varmin,in_varmax,in_stag,ierr)
         if (ierr.ne.0) goto 997
      endif
      call putdat(cdfid,varname,in_time,0,vv_anom,ierr)
      if (ierr.ne.0) goto 997

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


c     -----------------------------------------------------------------
c     Exception handling and format specs
c     -----------------------------------------------------------------

      stop

 998  print*,'Problems with input netcdf file'
      stop

 997  print*,'Problems with output netcdf file'
      stop

      END


c     ******************************************************************
c     * 2D FILTERING                                                   *
c     ******************************************************************

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     * AUXILIARY SUBROUTINES                                          *
c     ******************************************************************
   
C     -----------------------------------------------------------------
c     Rene's KINK function (for smoothing at bounadries)
C     -----------------------------------------------------------------

      real function kink (x,a)

      implicit none

c     declaration of parameters
      real   x,a

c     parameters
      real   pi
      parameter  (pi=3.1415926535)

      if (x.lt.0.) then
         kink=0.
      elseif (x.gt.a) then
         kink=1.
      else
         kink=(1.+cos(pi*(x-a)/a))/2.
      endif
    
      return
      end    


c     ---------------------------------------------------------
c     Subroutines to write the netcdf output file
c     ---------------------------------------------------------

      subroutine writecdf2D(cdfname,
     >                      varname,arr,time,
     >                      dx,dy,xmin,ymin,nx,ny,
     >                      crefile,crevar)

c     Create and write to the netcdf file <cdfname>. The variable
c     with name <varname> and with time <time> is written. The data
c     are in the two-dimensional array <arr>. The list <dx,dy,xmin,
c     ymin,nx,ny> specifies the output grid. The flags <crefile> and
c     <crevar> determine whether the file and/or the variable should
c     be created. 1: create / 0: not created

      IMPLICIT NONE

c     Declaration of input parameters
      character*80 cdfname,varname
      integer      nx,ny
      real         arr(nx,ny)
      real         dx,dy,xmin,ymin
      real         time
      integer      crefile,crevar

c     Further variables
      real         varmin(4),varmax(4),stag(4)
      integer      ierr,cdfid,ndim,vardim(4)
      character*80 cstname
      real         mdv
      integer      datar(14),stdate(5)
      integer      i

C     Definitions 
      varmin(1)=xmin
      varmin(2)=ymin
      varmin(3)=1050.
      varmax(1)=xmin+real(nx-1)*dx
      varmax(2)=ymin+real(ny-1)*dy
      varmax(3)=1050.
      cstname=trim(cdfname)//'_cst'
      ndim=4
      vardim(1)=nx
      vardim(2)=ny
      vardim(3)=1
      stag(1)=0.
      stag(2)=0.
      stag(3)=0.
      mdv=-999.98999

C     Create the file
      if (crefile.eq.0) then
         call cdfwopn(cdfname,cdfid,ierr)
         if (ierr.ne.0) goto 906
      else if (crefile.eq.1) then
         call crecdf(cdfname,cdfid,
     >        varmin,varmax,ndim,cstname,ierr)
         if (ierr.ne.0) goto 903 

C        Write the constants file
         datar(1)=vardim(1)
         datar(2)=vardim(2)
         datar(3)=int(1000.*varmax(2))
         datar(4)=int(1000.*varmin(1))
         datar(5)=int(1000.*varmin(2))
         datar(6)=int(1000.*varmax(1))
         datar(7)=int(1000.*dx)
         datar(8)=int(1000.*dy)
         datar(9)=1
         datar(10)=1
         datar(11)=0            ! data version
         datar(12)=0            ! cstfile version
         datar(13)=0            ! longitude of pole
         datar(14)=90000        ! latitude of pole     
         do i=1,5
            stdate(i)=0
         enddo
         call wricst(cstname,datar,0.,0.,0.,0.,stdate)
      endif

c     Write the data to the netcdf file, and close the file
      if (crevar.eq.1) then
         call putdef(cdfid,varname,ndim,mdv,
     >        vardim,varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 904
      endif
      call putdat(cdfid,varname,time,0,arr,ierr)
      if (ierr.ne.0) goto 905
      call clscdf(cdfid,ierr)

      return
       
c     Error handling
 903  print*,'*** Problem to create netcdf file ***'
      stop
 904  print*,'*** Problem to write definition ***'
      stop
 905  print*,'*** Problem to write data ***'
      stop
 906  print*,'*** Problem to open netcdf file ***'
      stop

      END


c     ------------------------------------------------------------------
c     Auxiliary routines
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> keeps 
c     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