#include "options.inc"


c-----------------------------------------------------------------------
! Time filtering a la Killworth
! to data where it is apropriate
! which means taux,tauy and ustar and maybe obc psi
! Notes:
! program uses LAPACK routines to invert tridiagonal matrix
! a calendar with non equal months but equal years is assumed
c-----------------------------------------------------------------------

      program killfilt
      implicit none

      call killfilt_it('forcing_folded_wind.cdf','taux')
      call killfilt_it('forcing_folded_wind.cdf','tauy')
      call killfilt_it('forcing_folded_wind.cdf','ustar')
      call cut_negative_values('forcing_folded_wind.cdf','ustar')
      call killfilt_it('forcing_folded_obc.cdf','psi_obc_south')


      call killfilt_it('forcing_folded.cdf','tauxBB')
      call killfilt_it('forcing_folded.cdf','tauyBB')
      call killfilt_it('forcing_folded.cdf','hflxBB')
      call killfilt_it('forcing_folded.cdf','sstm_BB')

      call killfilt_it('forcing_folded_sflx.cdf','sst_flux')
      call killfilt_it('forcing_folded_sflx.cdf','sss_flux')

      end program killfilt
 

      subroutine killfilt_it(fname,vname)
      implicit none
      character (len=*) :: fname, vname
#include "netcdf.inc"
      integer :: iret,ncid,vid,dimid(5),dim(5),ndim,n,xtype
      integer :: start(5),count(5)
      real (kind=4), allocatable :: buf(:,:,:)
      real (kind=4) :: spval
      character (len=80) :: s

      print*,' applying time filter a la Killworth et al to ',
     &  'variable ',vname(1:len_trim(vname)),' in file ',
     &   fname(1:len_trim(fname))

      iret=nf_open(fname,NF_WRITE,ncid)
      if (iret/=0) then
        print*,' ERROR: cannot open file'
        stop
      endif
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      iret=nf_inq_varid(ncid,vname,vid)
      if (iret/=0) then
        print*,' ERROR: cannot find variable'
        call ncclos (ncid, iret)
        return
      endif
      iret=NF_INQ_VARNDIMS(ncid,vid, ndim)
      if (ndim>3) then
        print*,' more than 3 dimensions are not allowed so far'
        print*,' but found ',ndim
        call ncclos (ncid, iret)
        stop
      endif
      if (ndim==1) then
        print*,' less than 2 dimensions are not allowed so far'
        print*,' but found ',ndim
        call ncclos (ncid, iret)
        stop
      endif
      iret=NF_INQ_VARTYPE(ncid,vid, xtype)
      if (xtype /= ncfloat) then
        print*,' cannot handle other types than ncfloat=',ncfloat
        print*,' but found ',xtype
        call ncclos (ncid, iret)
        stop
      endif

      iret = nf_get_att_real(ncid, vid, 'missing_value', spval)
      if (iret/=0) then
        print*,' ERROR: cannot find attribute missing value'
        print*,' for variable'
        call ncclos (ncid, iret)
        stop
      endif

      iret = nf_get_att_text(ncid, vid, 'killworth_filtered', s)
      if (iret==0) then
        print*,' ERROR: variable was already filtered'
        call ncclos (ncid, iret)
        stop
      endif

      iret=NF_INQ_VARDIMID (ncid,vid,dimid)
      dim=1
      do n=1,ndim
       iret=NF_INQ_DIMLEN(ncid, dimid(n), dim(n))
      enddo
      print*,' dimensions : ',dim(1:ndim)

      allocate( buf(dim(1),dim(2),dim(3)) )
      start=1; count=dim;
      iret= nf_get_vara_real(ncid,vid,start,count,buf)
      if (ndim==3) then
         call prep_killfilt(dim(ndim),dim(1),dim(2),buf,buf,spval)
      elseif (ndim==2) then
       call prep_killfilt(dim(ndim),dim(1),1,buf,buf,spval)
      else
       stop
      endif
      iret= nf_put_vara_real(ncid,vid,start,count,buf)

      deallocate(buf)

      iret = nf_redef(ncid)
      iret= nf_put_att_text(ncid,vid ,'killworth_filtered',4,'yeah')
      iret = nf_enddef(ncid)

      call ncclos (ncid, iret)

      end subroutine


      subroutine cut_negative_values(fname,vname)
      implicit none
      character (len=*) :: fname, vname
#include "netcdf.inc"
      integer :: iret,ncid,vid,dimid(5),dim(5),ndim,n,xtype
      integer :: start(5),count(5)
      real (kind=4), allocatable :: buf(:,:,:)
      real (kind=4):: spval
      character (len=80) :: s

      print*,' setting negative values to zero for ',
     &  'variable ',vname(1:len_trim(vname)),' in file ',
     &   fname(1:len_trim(fname))

      iret=nf_open(fname,NF_WRITE,ncid)
      if (iret/=0) then
        print*,' ERROR: cannot open file'
        stop
      endif
      iret=nf_set_fill(ncid, NF_NOFILL, iret)

      iret=nf_inq_varid(ncid,vname,vid)
      if (iret/=0) then
        print*,' ERROR: cannot find variable'
        call ncclos (ncid, iret)
        return
      endif
      iret=NF_INQ_VARNDIMS(ncid,vid, ndim)
      if (ndim>3) then
        print*,' more than 3 dimensions are not allowed so far'
        print*,' but found ',ndim
        call ncclos (ncid, iret)
        stop
      endif
      if (ndim==1) then
        print*,' less than 2 dimensions are not allowed so far'
        print*,' but found ',ndim
        call ncclos (ncid, iret)
        stop
      endif
      iret=NF_INQ_VARTYPE(ncid,vid, xtype)
      if (xtype /= ncfloat) then
        print*,' cannot handle other types than ncfloat=',ncfloat
        print*,' but found ',xtype
        call ncclos (ncid, iret)
        stop
      endif

      iret = nf_get_att_real(ncid, vid, 'missing_value', spval)
      if (iret/=0) then
        print*,' ERROR: cannot find attribute missing value'
        print*,' for variable'
        call ncclos (ncid, iret)
        stop
      endif

      iret = nf_get_att_text(ncid, vid, 'negative_values_cutted', s)
      if (iret==0) then
        print*,' ERROR: neg. values of variable already cutted'
        call ncclos (ncid, iret)
        stop
      endif

      iret=NF_INQ_VARDIMID (ncid,vid,dimid)
      dim=1
      do n=1,ndim
       iret=NF_INQ_DIMLEN(ncid, dimid(n), dim(n))
      enddo
      print*,' dimensions : ',dim(1:ndim)

      allocate( buf(dim(1),dim(2),dim(3)) )
      start=1; count=dim;
      iret= nf_get_vara_real(ncid,vid,start,count,buf)
      where( buf /=spval .and. buf <0.) buf=0.
      iret= nf_put_vara_real(ncid,vid,start,count,buf)

      deallocate(buf)

      iret = nf_redef(ncid)
      iret= nf_put_att_text(ncid,vid ,'negative_values_cutted',4,'yeah')
      iret = nf_enddef(ncid)

      call ncclos (ncid, iret)

      end subroutine


      subroutine prep_killfilt(itm,idim, jdim, rawdata, outdata, 
     &                         spval)
      implicit none
      integer :: itm,idim,jdim
      real (kind=4) :: rawdata(idim,jdim,itm), outdata(idim,jdim,itm)
      real (kind=4) :: spval
      real (kind=4) :: datold(itm),datnew
      real (kind=4), save, allocatable, dimension(:,:) :: pkf       
      logical,save :: first=.true.
      integer, save :: itm_save 
      integer i,j,m,l

      if (first .or. itm/=itm_save) then
        if (.not. first ) deallocate(pkf)
        allocate(pkf(itm,itm))
        call prep_killfilt_init (itm,pkf)
        first = .false.
        itm_save = itm
      end if

      do j=1, jdim
        do i=1, idim
          do m=1,itm
           datold(m) = rawdata(i,j,m)
          enddo
          do m=1, itm
            datnew = 0.
            do l=1, itm
              if (datold(l) .ne. spval) then
                datnew = datnew + pkf(m,l)*datold(l)
              else
                datnew = spval
              end if
            end do
            outdata(i,j,m) = datnew
          end do
        end do
      end do
      end subroutine prep_killfilt


      subroutine prep_killfilt_init (itm,pkf)
      implicit none
      real (kind=4), parameter :: r8 = 0.125,  p75 = 0.75
      real (kind=4), parameter :: c0 = 0., c1 = 1.
      integer :: itm
      integer :: indx (itm)
      real (kind=4) :: fkp (itm,itm)  ! interpolation matrix, automatic array
      real (kind=4) :: pkf (itm,itm)   !       inverse matrix
      real (kind=4) :: month(itm)
      real (kind=4) :: mmonth(12)
      data  mmonth/31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./
      integer :: lwork ,m,n,i,info

      print*,' initializing time interpolation matrix'
      lwork=itm*itm
      do m=1,12
       do n=m,itm,12
        month(n)=mmonth(m)
       enddo
      enddo
      fkp = c0
      do i=2,itm
       fkp(i,i-1)=month(i)/( 4.0*(month(i-1)+month(i)))
      enddo
      fkp(1,itm)=month(1)/( 4.0*(month(itm)+month(1)))
      do i=1,itm-1
       fkp(i,i+1)=month(i)/( 4.0*(month(i)+month(i+1) ))
      enddo
      fkp(itm,1)=month(itm)/( 4.0*(month(itm)+month(1)))
      do  i=2,itm-1
       fkp(i,i)=1.0-fkp(i,i-1)-fkp(i,i+1)
      enddo
      fkp(1,1)=1.0-fkp(1,itm)-fkp(1,2)
      fkp(itm,itm)=1.0-fkp(itm,itm-1)-fkp(itm,1)
      print*,' inverting matrix using lapack '
      call sgetrf(itm,itm,fkp,itm,indx,info)
      call sgetri(itm,fkp,itm,indx,pkf,lwork,info)
      print*,' done'
      pkf = fkp
      end  subroutine prep_killfilt_init


