#include "options.inc"


      module nc_util_module
      implicit none
      integer ::  memory = 1 !  0 high memory, 1 low memory
      logical :: verbose = .false.
!      logical :: verbose = .true.
      end module nc_util_module


      subroutine mult_variable(ncin,name,fx)
c=======================================================================
c   multiply variable name in open netcdf file ncin with 
c   constant factor fx
c=======================================================================
      use nc_util_module
      implicit none
      integer :: ncin
      character (len=*) :: name
      integer, parameter :: ndmax=100
      integer, dimension(ndmax)  :: dims,start,count,start2,count2
      integer :: iret,ndims,xtype,varin,n
      real , allocatable :: buf(:)
      integer :: buflen,i1,i2
      real :: spval,fx
#include "netcdf.inc"

      if (verbose) print*,' multiply variable ',name(1:len_trim(name)),
     &   ' with ',fx
      iret = NF_INQ_VARID(NCin,name,varin) 
      if (iret/=0) goto 20
      iret = NF_INQ_VARNDIMS(NCin, VARin, ndims)
      if (iret/=0) goto 20
      if (ndims>ndmax) then
        print*,' ERROR: number of dims exceed ',ndmax
        stop
      endif
      iret = NF_INQ_VARTYPE (NCin, VARin, xtype)
      iret = NF_INQ_VARdimid (NCin, VARin,dims)
      buflen=1
      if (memory==0 .or. ndims < 3) then
       do n=1,ndims
        start(n)=1
        iret=NF_INQ_DIMlen (NCin, DIMs(n), count(n))
        buflen=buflen*count(n)
       enddo
       allocate(buf(buflen))
       iret= nf_get_vara_real(ncin,varin,start,count,buf)
       iret = nf_get_att_real(ncin,varin, 'missing_value', spval) 
       if (iret==0) then
        if (verbose) print*,' found special value of ',spval
        where(buf/=spval) buf =buf*fx
       else
         buf =buf*fx
       endif
       iret= nf_put_vara_real(ncin,varin,start,count,buf)
       deallocate(buf)
      elseif (memory == 1) then

       do n=1,ndims
        start2(n)=1
        iret=NF_INQ_DIMlen (NCin, DIMs(n), count2(n))
       enddo
       do n=1,ndims-1
         buflen=buflen*count2(n)
       enddo
       allocate(buf(buflen))
       do i1=1,count2(ndims)
         if (verbose) print'(a,i5,a,i5)','  i=',i1,'/',count2(ndims)
         start=start2; count=count2
         start(ndims)=i1; count(ndims)=1
         iret= nf_get_vara_real(ncin,varin,start,count,buf)
         iret = nf_get_att_real(ncin,varin, 'missing_value', spval) 
         if (iret==0) then
           where(buf/=spval) buf =buf*fx
         else
           buf =buf*fx
         endif
         iret= nf_put_vara_real(ncin,varin,start,count,buf)
       enddo
       deallocate(buf)

      else
        print*,' ERROR: memory flag >1 '
        stop
      endif
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine mult_variable


      subroutine add_to_variable(ncin,ncout,name)
c=======================================================================
c   add variable name in open netCDF file ncin to same
c   variable in open file ncout and store results in ncout
c   both file must be in data mode
c=======================================================================
      use nc_util_module
      implicit none
      integer :: ncin,ncout
      character (len=*) :: name
      integer, parameter :: ndmax=100
      integer, dimension(ndmax)  :: dims,start,count,start2,count2
      integer :: iret,ndims,xtype,varin,varout,n
      real , allocatable :: buf(:),buf2(:)
      integer :: buflen,i1,i2
      real :: spval
#include "netcdf.inc"

      if (verbose) print*,' adding to variable ',name(1:len_trim(name))
      iret = NF_INQ_VARID(NCin,name,varin) 
      if (iret/=0) goto 20
      iret = NF_INQ_VARNDIMS(NCin, VARin, ndims)
      if (ndims>ndmax) then
        print*,' ERROR: number of dims exceed ',ndmax
        stop
      endif
      iret = NF_INQ_VARTYPE (NCin, VARin, xtype)
      iret = NF_INQ_VARdimid (NCin, VARin,dims)
      iret = NF_INQ_VARID(NCout,name,varout) 
      if (iret/=0) goto 20
      buflen=1
      if (memory==0 .or. ndims < 3) then
       do n=1,ndims
        start(n)=1
        iret=NF_INQ_DIMlen (NCin, DIMs(n), count(n))
        buflen=buflen*count(n)
       enddo
       allocate(buf(buflen),buf2(buflen))
       iret= nf_get_vara_real(ncin,varin,start,count,buf)
       iret= nf_get_vara_real(ncout,varout,start,count,buf2)
       iret = nf_get_att_real(ncin,varin, 'missing_value', spval) 
       if (iret==0) then
        if (verbose) print*,' found special value of ',spval
        where(buf/=spval) buf =buf+buf2
       else
        buf =buf+buf2
       endif
       iret= nf_put_vara_real(ncout,varout,start,count,buf)
       deallocate(buf,buf2)
      elseif (memory == 1) then

       do n=1,ndims
        start2(n)=1
        iret=NF_INQ_DIMlen (NCin, DIMs(n), count2(n))
       enddo
       do n=1,ndims-1
         buflen=buflen*count2(n)
       enddo
       allocate(buf(buflen),buf2(buflen))
       do i1=1,count2(ndims)
         if (verbose) print'(a,i5,a,i5)','  i=',i1,'/',count2(ndims)
         start=start2; count=count2
         start(ndims)=i1; count(ndims)=1
         iret= nf_get_vara_real(ncin,varin,start,count,buf)
         iret= nf_get_vara_real(ncout,varout,start,count,buf2)
         iret = nf_get_att_real(ncin,varin, 'missing_value', spval) 
         if (iret==0) then
          where(buf/=spval) buf =buf+buf2
         else
           buf =buf+buf2
         endif
         iret= nf_put_vara_real(ncout,varout,start,count,buf)
       enddo
       deallocate(buf,buf2)
      else
        print*,' ERROR: memory flag >1 '
        stop
      endif
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine add_to_variable



      recursive subroutine copy_variable(ncin,ncout,name,expand,nlen)
c=======================================================================
c   define variable name in ncin in ncout.
c   copy also values. both files must be in definition mode.
c   if expand = true then expand last dimension to nlen
c=======================================================================
      use nc_util_module
      implicit none
      integer  :: ncin,ncout
      character (len=*) :: name
      logical :: expand
      integer :: nlen
      integer :: iret,varin,xtype,ndims,natts,varout,n
      integer, parameter :: ndmax=100
      integer, dimension(ndmax)  :: dimin,dimout,start,count
      real , allocatable :: buf(:)
      integer :: buflen
      character (len=80) :: dimname
#include "netcdf.inc"
c    inquire
      iret=nf_sync(ncout)
      iret = NF_INQ_VARID(NCin,name,varin) 
      if (iret/=0) goto 20
      iret = NF_INQ_VARNDIMS(NCin, VARin, ndims)
      if (ndims>ndmax) then
        print*,' ERROR: number of dims exceed ',ndmax
        stop
      endif
      iret = NF_INQ_VAR (NCin, VARin,name,xtype, ndims,dimin, natts)
      if (verbose) then 
         print'(a,a,a,i3,a,i3,a,i3,a)',
     &  ' copying variable ',name(1:len_trim(name)),
     & ' (# of dims :',ndims,' # of attr :',natts,' type :',xtype,')'
         if (expand) print*,' expanding last dimension to ',nlen
      endif
c   define dimensions if needed, expand last dimension if requested
      do n=1,ndims
        dimname=''
        iret=NF_INQ_DIMNAME (NCin, DIMIn(n), dimname)
        iret=NF_INQ_DIMID(NCout,dimname,dimout(n)) 
        if (iret/=0) then 
          if (n==ndims) then
           call copy_dimension(ncin,ncout,dimname,expand,nlen)
          else
           call copy_dimension(ncin,ncout,dimname,.false.,0)
          endif
        endif
      enddo
c   define variable
      if (verbose) print*,' variable ',name(1:len_trim(name)),
     &            ' resides on the following dimensions'
      do n=1,ndims
        dimname=''
        iret=NF_INQ_DIMNAME (NCin, DIMIn(n), dimname)
        if (verbose) print*,dimname(1:len_trim(dimname))
        iret=NF_INQ_DIMID(NCout,dimname,dimout(n)) 
      enddo
      iret=NF_INQ_VARID(NCout,name,varout) 
      if (iret==0) return
      if (verbose) print*,' defining variable ',name(1:len_trim(name)),
     &   ' in output file'
      varout = ncvdef (ncout,name, xtype,ndims,dimout,iret)
c   define attributes
      do n=1,natts
       dimname=''
       iret= NF_INQ_ATTNAME (NCin, varin, n, dimname)
       iret= NF_COPY_ATT(NCIN, varin, dimname, NCOUT, varout)
      enddo
c   copy values of variable 
      iret = nf_enddef(ncout)
      buflen=1; start=1
      do n=1,ndims
        iret=NF_INQ_DIMlen (NCin, DIMIn(n), count(n))
        buflen=buflen*count(n)
      enddo
      if (xtype == ncdouble) buflen=buflen*2
      allocate(buf(buflen))
      if (xtype == nf_int .or. xtype == nf_float  .or.
     &    xtype == nf_byte .or. xtype == nf_short  .or.
     &    xtype == nf_char  ) then
       iret= nf_get_vara_real(ncin,varin,start,count,buf)
       iret= nf_put_vara_real(ncout,varout,start,count,buf)
      elseif (xtype == nf_double) then
       iret= nf_get_vara_double(ncin,varin,start,count,buf)
       iret= nf_put_vara_double(ncout,varout,start,count,buf)
      else
       print*,' do not know how to handle NCtype ',xtype
       stop
      endif
      deallocate(buf)
      iret = nf_redef(ncout)
      iret=nf_sync(ncout)
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine copy_variable

      recursive subroutine copy_variable_novalues
&                         (ncin,ncout,name,expand,nlen)
c=======================================================================
c   define variable name in ncin in ncout.
c   but do not copy values of variable. both files must be in definition mode.
c   if expand = true then expand last dimension to nlen
c=======================================================================
      use nc_util_module
      implicit none
      integer  :: ncin,ncout
      character (len=*) :: name
      logical :: expand
      integer :: nlen
      integer :: iret,varin,xtype,ndims,natts,varout,n
      integer, parameter :: ndmax=100
      integer, dimension(ndmax)  :: dimin,dimout,start,count
      real , allocatable :: buf(:)
      integer :: buflen
      character (len=80) :: dimname
#include "netcdf.inc"
c    inquire
      iret=nf_sync(ncout)
      iret = NF_INQ_VARID(NCin,name,varin) 
      if (iret/=0) goto 20
      iret = NF_INQ_VARNDIMS(NCin, VARin, ndims)
      if (ndims>ndmax) then
        print*,' ERROR: number of dims exceed ',ndmax
        stop
      endif
      iret = NF_INQ_VAR (NCin, VARin,name,xtype, ndims,dimin, natts)
      if (verbose) then 
         print'(a,a,a,i3,a,i3,a,i3,a)',
     &  ' copying variable ',name(1:len_trim(name)),
     & ' (# of dims :',ndims,' # of attr :',natts,' type :',xtype,')'
         if (expand) print*,' expanding last dimension to ',nlen
      endif
c   define dimensions if needed, expand last dimension if requested
      do n=1,ndims
        dimname=''
        iret=NF_INQ_DIMNAME (NCin, DIMIn(n), dimname)
        iret=NF_INQ_DIMID(NCout,dimname,dimout(n)) 
        if (iret/=0) then 
          if (n==ndims) then
           call copy_dimension(ncin,ncout,dimname,expand,nlen)
          else
           call copy_dimension(ncin,ncout,dimname,.false.,0)
          endif
        endif
      enddo
c   define variable
      if (verbose) print*,' variable ',name(1:len_trim(name)),
     &            ' resides on the following dimensions'
      do n=1,ndims
        dimname=''
        iret=NF_INQ_DIMNAME (NCin, DIMIn(n), dimname)
        if (verbose) print*,dimname(1:len_trim(dimname))
        iret=NF_INQ_DIMID(NCout,dimname,dimout(n)) 
      enddo
      iret=NF_INQ_VARID(NCout,name,varout) 
      if (iret==0) return
      if (verbose) print*,' defining variable ',name(1:len_trim(name)),
     &   ' in output file'
      varout = ncvdef (ncout,name, xtype,ndims,dimout,iret)
c   define attributes
      do n=1,natts
       dimname=''
       iret= NF_INQ_ATTNAME (NCin, varin, n, dimname)
       iret= NF_COPY_ATT(NCIN, varin, dimname, NCOUT, varout)
      enddo
      iret=nf_sync(ncout)
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine copy_variable_novalues


      
      recursive subroutine copy_dimension(ncin,ncout,name,expand,nlen)
c=======================================================================
c   copy dimension name from ncin to ncout
c   both files must be in definition mode
c   if expand = true then dimension to nlen
c=======================================================================
      use nc_util_module
      implicit none
      integer :: ncin,ncout
      character (len=*) :: name
      logical :: expand
      integer :: nlen
      integer :: iret,dimin,dimout,len,varin,varout,n,natts
      character (len=80) :: attname,s
#include "netcdf.inc"
      iret=nf_sync(ncout)
      iret=NF_INQ_DIMID(NCin,name,dimin) 
      if (iret/=0) goto 20
      iret= NF_INQ_DIMLEN (NCin, DIMin, len)
      if (verbose) then
       print*,' copying dimension ',name(1:len_trim(name)),
     &      ' (length :',len,')'
       if (expand) print*,' expanding dimension to ',nlen
      endif
      if (expand) len=nlen
      dimout = ncddef(ncout,name,len,iret)
      if (iret/=0) goto 20
      iret=NF_INQ_VARID(NCin,name,varin) 
      if (iret==0) then 
        iret=NF_INQ_VARID(NCout,name,varout) 
        if (iret/=0) call copy_variable(ncin,ncout,name,expand,nlen)       
c     care about pointers to other axis
        iret=NF_INQ_VARnatts(ncin,varin,natts) 
        do n=1,natts
         attname=''; iret= NF_INQ_ATTNAME (NCin, varin, n, attname)
         if (attname == 'edges' .or. attname == 'EDGES' ) then
          s='';iret= nf_get_att_text(ncin,varin,attname,s)
          if (verbose) print*,' found edges definition :', 
     &           s(1:len_trim(s)),', defining as well'
          iret=NF_INQ_varid(NCin,s,varin) 
          if (iret==0) then 
            call copy_variable(ncin,ncout,s,.false.,0)
          else
            if (verbose) print*,' cannot find variable ',
     &                   s(1:len_trim(s)),' in file'
          endif
         endif
        enddo

      endif
      iret=nf_sync(ncout)
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine copy_dimension




      subroutine copy_file_attributes(ncin,ncout)
c=======================================================================
c   copy global attributes from ncin to ncout
c   both files must be in definition mode.
c=======================================================================
      use nc_util_module
      implicit none
      integer  :: ncin,ncout
      integer :: iret,natts,n
      character (len=80) :: name
#include "netcdf.inc"
c    inquire
      iret = NF_INQ_natts(NCin,natts) 
      if (iret/=0) goto 20
c   define attributes
      do n=1,natts
       name=''
       iret= NF_INQ_ATTNAME (NCin, nf_global, n, name)
       if (verbose) print*,'copying attribute : ',name
       iret= NF_COPY_ATT(NCIN, nf_global, name, NCOUT, nf_global)
      enddo
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine copy_file_attributes





      subroutine search_nr_variables(ncin,nr)
c=======================================================================
c=======================================================================
      use nc_util_module
      implicit none
      integer :: ncin,nr,iret
      integer :: nd,n,m,nr2,i
      character (len=80) :: name,dname
#include "netcdf.inc"

      iret = NF_INQ_nvars(NCin,nr) 
      iret = NF_INQ_ndims(NCin,nd) 
      if (iret/=0) goto 20
c   search for axis variables
      nr2=nr
      do m=1,nr
        iret = NF_INQ_varname(NCin,m,name) 
        i=len_trim(name)
        do n=1,nd
         iret = NF_INQ_dimname(NCin,n,dname) 
         if (name(1:i) == dname(1:i) ) nr2=nr2-1
        enddo
      enddo
      nr=max(0,nr2)
      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
      end subroutine search_nr_variables




      subroutine search_variables(ncin,nri,name)
c=======================================================================
c=======================================================================
      use nc_util_module
      implicit none
      integer :: ncin,nr,iret,m,nd,nr2,nri,nn,n
      character (len=*) :: name(nri)
      character (len=80) :: vname,dname
      integer :: i
#include "netcdf.inc"
      iret = NF_INQ_nvars(NCin,nr) 
      iret = NF_INQ_ndims(NCin,nd) 
      nr2=nr
      if (iret/=0) goto 20
c   search for axis variables
      nn=1
      do m=1,nr
        iret = NF_INQ_varname(NCin,m,vname) 
        name(nn)=vname
        i=len_trim(vname)
        do n=1,nd
         iret = NF_INQ_dimname(NCin,n,dname) 
         if (vname(1:i) == dname(1:i) ) then
          nr2=nr2-1
          nn=nn-1
         endif
        enddo
        nn=nn+1
        if (nn>nri) goto 30
      enddo
      nr=max(0,nr2)


      return
 20   print*,'ERROR: netcdf lib says:',nf_strerror(iret)
      stop
 30   print*,'ERROR: variables count wrong'
      stop
      end subroutine search_variables
