

#define skip_interp_details


      subroutine tmp_filename(name)
c--------------------------------------------------------------
c    returns temporary file name
c--------------------------------------------------------------
      implicit none
      character (len=*) :: name
      integer :: i
      real :: x
      logical :: ex=.true.
      do while (ex) 
       do i=1,len(name)
         call random_number(x)
         name(i:i)=char(int(x*25)+65)
       enddo
       inquire(file=name,exist=ex)
      end do
      end subroutine tmp_filename

      subroutine get_commandline_argnr(nr)
c--------------------------------------------------------------
c   returns number of command line arguments
c--------------------------------------------------------------
      integer, intent(out)  :: nr
      integer :: ierr
      character (len=80)   :: arg
      ierr=0; nr=0
 20   call get_commandline_arg(nr,arg,ierr)
      if (ierr==0 .and. arg /= '') then
         nr=nr+1; goto 20
      endif
      nr=nr-1
      end subroutine get_commandline_argnr

      subroutine get_commandline_arg(nr,arg,ierr)
c--------------------------------------------------------------
c     reads command line argument # nr
c     on error ierr is greater than zero
c--------------------------------------------------------------
      integer, intent(in)  :: nr
      integer, intent(out) :: ierr
      character (len=80)   :: arg
      integer :: ilen

      ierr=0
#if defined C90_host || defined T3E_host || defined IRIX_host
      call pxfgetarg(nr,arg,ilen,ierr)
#elif defined SR8000_host
      call getarg(nr+1,arg)
#else
      call getarg(nr,arg)
#endif
      end subroutine get_commandline_arg

      subroutine get_commandline_arg_int(nr,out,ierr)
c--------------------------------------------------------------
c     reads an integer from command line argument # nr
c     return integer out, on error ierr is greater than zero
c--------------------------------------------------------------
      integer, intent(in) :: nr
      integer, intent(out) :: out,ierr
      character (len=80) :: arg
      integer :: ilen
      ierr=0
#if defined C90_host || defined T3E_host || defined IRIX_host
      call pxfgetarg(nr,arg,ilen,ierr)
      if (ierr/=0) return
      read(arg,'(i)',err=9999,end=9999) out
#elif defined ALPHA_host || defined SX5_host
      call getarg(nr,arg)
      read(arg,'(i)',err=9999,end=9999) out
#elif defined SR8000_host
      call getarg(nr+1,arg)
      read(arg,*,err=9999,end=9999) out
#else
      call getarg(nr,arg)
      read(arg,*,err=9999,end=9999) out
#endif
      return
 9999 ierr=1;return
      end subroutine get_commandline_arg_int


      subroutine sub_flush(io)
c-----------------------------------------------------------------------
c     flushes IO unit io to disk
c-----------------------------------------------------------------------
      implicit none
      integer, intent(in) :: io
      integer :: iret
#if defined IRIX_host || defined SX5_host|| defined SR8000_host
      call flush(io,iret)
#elif defined T3E_host
      call flush(io,iret)
#elif defined C90_host
      call flush(io,iret)
#elif defined ALPHA_host
      call flush(io)
#elif defined SUN_host
      call flush(io)
#elif defined AIX_host
      call flush_(io)
#else
      ! do nothing
#endif
      end subroutine sub_flush


      subroutine get_free_iounit (nu)
c-----------------------------------------------------------------------
c     returns the first free IO unit number in nu
c-----------------------------------------------------------------------
      implicit none
      integer nu,n
      logical in_use
      character (len=80) :: name
      do n=7,99
        inquire (n, OPENED=in_use, NAME=name)
#if defined C90_host || defined T3E_host
	if (name .ne. '') then
	  in_use = .true.
	endif
#endif
	if (.not. in_use) then
	  nu = n
	  go to 10
	endif
      enddo
      print *,'Error: exhausted available unit numbers'
      print *,'       Are you forgetting to close units?'
      call halt_stop(' in get_free_iounit')
10    continue
      end subroutine get_free_iounit



      subroutine tolower (string)
      implicit none
      character (len=*) :: string
      character (len=1) :: c
      integer, parameter :: lcshift=32
      integer :: i
      do i=1,len(string)
        c = string(i:i)
        if ('A' .le. c .and. c .le. 'Z') then
          string(i:i) = char(ichar(c) + lcshift)
        end if
      end do
      end subroutine tolower

      subroutine replace_space_zero(name)
      implicit none
      character (len=*) :: name
      integer  :: i
      do i=1,len_trim(name)
          if (name(i:i)==' ')name(i:i)='0'
      enddo
      end subroutine replace_space_zero

      subroutine no_space_chars (string)
      implicit none
      character (len=*) :: string
      integer :: i
      do i=1,len(string)
       if (string(i:i).eq.' ') then
        string(i:len(string)-1)=string(i+1:len(string))
       endif
      enddo 
      end subroutine no_space_chars


      subroutine find_number (number, string)
      implicit none
      character (len=*) :: string
      character (len=1) :: c
      integer :: i,lens,is,ie,number
      do i=len(string),1,-1
        if (string(i:i) .ne. ' ') then
	  lens = i; go to 10
	endif
      enddo
10    is = index (string, '=')
      if (is /= 0) then
        do i=is,lens
          c = string(i:i)
          if (c .ge. '0' .and. c .le. '9') go to 20
        enddo
      else
        print *,' Error: need "=" sign when specifying record length'
	print *,' attributes = ',string
	call halt_stop(' in find_number')
      endif
20    is=i; ie = lens
      do i=is,lens
        c = string(i:i)
        if (c .lt. '0' .or. c .gt. '9') then
	  ie = i-1; go to 30
	endif
      enddo
30    continue
      if (ie .ne. 0) then
        read (string(is:ie), *) number 
      else
        print *,' Error: cannot find record length: string=',string
	call halt_stop(' in find_number')
      endif
      end subroutine find_number


      subroutine getunit (nu, filename, optionlist)
c
c-----------------------------------------------------------------------
c     connect unit number "nu" to "filename" and open the file with
c     the requested options in the string optionlist
c     the options might be:
c     fsr | formatted sequential rewind   : write new formatted file
c     fsa | formatted sequential append   : append to formatted file
c     usr | unformatted sequential rewind : write new unformatted file
c     usa | unformatted sequential append : append to unformatted file
c     ud = recln | unformatted direct words= recln : write to direct access file
# if defined C90_host  || defined T3E_host || defined IRIX_host
c     There are three different assign types which can be specified
c     after the other options:   ieee, t3e, f77  and cray 
#endif
c     Default options are : formatted, sequential, rewind
c-----------------------------------------------------------------------
c
      implicit none
      integer :: nu
      character (len=*) :: filename,optionlist
      character (len=128) :: options,string
      logical :: file_exists = .false., direct = .false.
      character (len=20) :: form, position, access, status
      integer :: len_bytes,ierr,number,len_words
      call get_free_iounit (nu)
c     check for STATUS
      inquire ( FILE=filename, EXIST=file_exists )
      status='old'; if ( .not. file_exists) status='new'
c     prepare options to parse
      options=optionlist
      call tolower (options)
      call no_space_chars (options)
c     Formatted or unformatted
      if     ( index(options,'unformatted').gt.0) then
         form = 'unformatted'
      elseif ( index(options,'formatted').gt.0)   then
         form = 'formatted'
      elseif ( options(1:1).eq.'u') then
         form = 'unformatted'
      elseif ( options(1:1).eq.'f') then
         form = 'formatted'
      else
       print*,' in getunit : '
       print*,' cannot read optionlist (formatted or unformatted)'
       print*,' assuming formatted '
       form = 'formatted'
      endif
c     direct or sequential
      if     ( index(options,'direct').gt.0)      then
        direct = .true.
      elseif ( index(options,'sequential').gt.0)  then
         direct = .false.
      elseif ( options(2:2).eq.'d') then
         direct=.true.
      elseif ( options(2:2).eq.'s') then
         direct=.false.
      else
       print*,' in getunit : '
       print*,' cannot read optionlist (direct or sequential)'
       print*,' assuming sequential '
       direct = .false.
      endif
c     Sequential or direct access
      if (direct) then
        access = 'direct'
	call find_number (len_words, options)
c       convert to bytes assuming 8 bytes per word
        len_bytes = 8*len_words
      else
        access = 'sequential'
        if     ( index(options,'append').gt.0)   then
          position = 'append'
        elseif ( index(options,'rewind').gt.0)   then
          position = 'rewind'
        elseif ( options(3:3).eq.'a') then
          position = 'append'
        elseif ( options(3:3).eq.'r') then
          position = 'rewind'
        else
          print*,' in getunit : '
          print*,' cannot read optionlist (append or rewind)'
          print*,' assuming rewind'
          position = 'rewind'
        endif
      endif
# if defined C90_host  || defined T3E_host || defined IRIX_host
c     parse options for assign statements
       string=' '
       if  ( index(options,'ieee').gt.0) then
	  string = 'assign -F f77 -N ieee f:'//filename
       elseif  ( index(options,'t3e').gt.0) then
	  string = 'assign -N t3e f:'//filename
       elseif  ( index(options,'cray').gt.0) then
	  string = 'assign -N cray f:'//filename
       elseif  ( index(options,'f77').gt.0) then
	  string = 'assign -F f77 f:'//filename
       endif
       if (string /= ' ') then
        call assign (string, ierr)
        if (ierr .ne. 0) then
         print *,' assign failed: ierr=', ierr, ' string=',string
         call assign ('assign -V', ierr)
         call halt_stop(' in getunit')
        endif
       endif
# endif
c      Now open the file 
#  if defined test_iomngr 
       print*,'UNIT    =',nu
       print*,'FILE    = ',filename(1:len_trim(filename))
       print*,'FORM    =',form(1:len_trim(form))
       print*,'ACCESS  =',access(1:len_trim(access))
       print*,'STATUS  =',status(1:len_trim(status))
       print*,'POSITION=',position(1:len_trim(position))
       print*,'RECL    =',len_bytes
#  endif
      if (direct) then
       open (nu, FILE=filename, FORM=form, ACCESS=access
     &,       STATUS=status,RECL=len_bytes)
      else
c    some f77-compilers like GNUs g77 do not like POSITION, 
c    use ACCESS instead
c       open (nu, FILE=filename, FORM=form, ACCESS=position
c     &,       STATUS=status)
       open (nu, FILE=filename, FORM=form, ACCESS=access
     &,       STATUS=status,POSITION=position)
      endif
      end subroutine getunit




      function indp (value, array, ia)
c
c=======================================================================
c
c     indp = index of nearest data point within "array" corresponding to
c            "value".
c
c     inputs:
c
c     value  = arbitrary data...same units as elements in "array"
c     array  = array of data points  (must be monotonically increasing)
c     ia     = dimension of "array"
c
c     output:
c
c     indp =  index of nearest data point to "value"
c             if "value" is outside the domain of "array" then indp = 1
c             or "ia" depending on whether array(1) or array(ia) is
c             closest to "value"
c
c             note: if "array" is dimensioned array(0:ia) in the calling
c                   program, then the returned index should be reduced
c                   by one to account for the zero base.
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c
c     example:
c
c     let model depths be defined by the following:
c     parameter (km=5)
c     dimension z(km)
c     data z /5.0, 10.0, 50.0, 100.0, 250.0/
c
c     k1 = indp (12.5, z, km)
c     k2 = indp (0.0, z, km)
c
c     k1 would be set to 2, & k2 would be set to 1 so that
c     z(k1) would be the nearest data point to 12.5 and z(k2) would
c     be the nearest data point to 0.0
c
c=======================================================================
c
      integer indp
      integer :: stdout=6
      dimension array(ia)
c
      do i=2,ia
        if (array(i) .lt. array(i-1)) then
	  write (stdout,*)
     &   ' => Error: array must be monotonically increasing in "indp"' 
     &,  '           when searching for nearest element to value=',value 
	  write (stdout,*) '           array(i) < array(i-1) for i=',i 
          write (stdout,*) '           array(i) for i=1..ia follows:'
	  do ii=1,ia
	    write (stdout,*) 'i=',ii, ' array(i)=',array(ii)
          enddo
	  stop
	endif
      enddo
      if (value .lt. array(1) .or. value .gt. array(ia)) then
        if (value .lt. array(1))  indp = 1
        if (value .gt. array(ia)) indp = ia
        return
      else
        do i=2,ia
          if (value .le. array(i)) then
            indp = i
            if (array(i)-value .gt. value-array(i-1)) indp = i-1
            go to 101
          endif
        enddo
101     continue
      endif
      end function indp



      subroutine ftc(f, if, jf, xf, yf, c, ic, jc, istart, iend
     &,               jstart, jend, xc, yc, init, work, lenw,spval)
c
c=======================================================================
c
c     "ftc" is a nemonic for "fine to coarse".
c
c     obtain a coarse grid representaion of a fine grid dataset by area
c     averaging grid boxe values on the fine grid which overlay coarse 
c     grid boxes. note: the coarse grid boxes do not have to contain an
c     integral number of fine grid boxes.
c
c
c     inputs:
c
c     f      = data on fine grid
c     if     = inner dimension of "f"
c     jf     = outer dimension of "f"
c     xf     = coordinates for inner dimension of "f" (eg: longitues)
c     yf     = coordinates for outer dimension of "f" (eg: latitues)
c
c     ic     = inner dimension of coarse grid "c"
c     jc     = outer dimension of coarse grid "c"
c     istart = starting index along inner dimension of "c" for which
c              averaged values are desired
c     iend   = ending index along inner dimension of "c" for which
c              averaged values are desired
c     jstart = starting index along outer dimension of "c" for which
c              averaged values are desired
c     jend   = ending index along outer dimension of "c" for which
c              averaged values are desired
c     xc     = coordinates for inner dimension of "c" (eg: longitues)    
c     yc     = coordinates for outer dimension of "c" (eg: latitues)    
c     init   = initialize the averaging factors
c              "init" should be set = 1 on the first call. 
c              "init" <> 1 uses the previously computed factors stored
c              in "work" array.
c     work   = work array of averaging factors when "init" <> 1
c              (previously calculated by "ftc" when "init" = 1)
c     lenw   = size of work array. lenw should be >= 9*max(if,jf)
c     spval  = special value masking out dry points
c       
c     output:
c
c     c  = coarse grid average of "f" defined over
c          ((c(i,j),i=istart,iend),j=jstart,jend)
c     work   = work array of averaging factors when "init" = 1
c
c     restrictions:
c
c     fine and coarse grids are assumed rectangular with "xf" and "xc"
c     having the same units. "yf" and "yc" must also have the same units
c     the coarse domain xc(istart)...xc(iend) must be within 
c     the fine domain xf(1)...xf(if). similarly,
c     yc(jstart)...yc(jend) must be within yf(1)...yf(jf). all 
c     coordinates must be strictly monotonically increasing.
c 
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
!     modified:    J. Dengg, jdengg@awi-bremerhaven.de
!                  Sep-97
c=======================================================================
c
      integer :: stdout=6,indp
      logical error, show_coord
      parameter (len=10000, p5=0.5, c0=0.0)
      dimension iso(0:len), ieo(0:len), jso(0:len), jeo(0:len)
     &,         dx(0:len,2), dy(0:len,2), edgecx(0:len), edgecy(0:len)
     &,         edgefx(0:len), edgefy(0:len)
      dimension f(if,jf), xf(if), yf(jf)
      dimension c(ic,jc), xc(ic), yc(jc)
      dimension work(lenw)
      real,intent(in) :: spval
c
c-----------------------------------------------------------------------
c     initialize weights or use previously calculated weights
c-----------------------------------------------------------------------
c
      if (init .eq. 1) then
        error = .false.
c	write (stdout,*)
c     & '              (initializing the averaging weights)'
c        write (stdout,*) ' '
c
c       test to verify that array sizes do not exceed limits
c
        if (if .gt. len .or. jf .gt. len) then
          i = max(if,jf)
          write (stdout,*) '=>Error: increase "len" in "ftc" to ',i
          stop '=>ftc'        
        endif
	if (lenw .lt. 9*max(if,jf)) then
          write (stdout,*) '=>Error: increase size of "work" array',
     &      ' to at least ',9*max(if,jf),' for calls to "ftc"'
	  error = .true.
        endif
c
c       verify that the "coarse" grid lies within the "fine" grid
c
        if (xf(1) .gt. xc(istart) .or. xf(if) .lt. xc(iend)) then
          write (stdout,*)
     &     '=>Warning: Coarse grid "xc" is outside "fine" grid "xf".'
          if (xf(1) .gt. xc(istart)) then
	    write (stdout,*) '  xc(',istart,')  .lt.  xf(1)'  
	  endif
          if (xf(if) .lt. xc(iend)) then
	    write (stdout,*) '  xc(',iend,')  .gt.  xf(',if,')'  
	  endif
        endif 
        if (yf(1) .gt. yc(jstart) .or. yf(jf) .lt. yc(jend)) then
          write (stdout,*)
     &     '=>Warning: Coarse grid "yc" is outside "fine" grid "yf".'
          if (yf(1) .gt. yc(jstart)) then
	    write (stdout,*) '  yc(',jstart,') .lt. yf(1)' 
	  endif
          if (yf(jf) .lt. yc(jend)) then
	    write (stdout,*) '  yc(',jend,')  .gt.  yf(',jf,')'  
	  endif
        endif
c
c       construct edges of "coarse" grid boxes
c
        do i=1,ic-1
	  edgecx(i) = p5*(xc(i) + xc(i+1))
        enddo
        edgecx(0)  = xc(1) - (edgecx(1) - xc(1))
	edgecx(ic) = xc(ic) + (xc(ic) - edgecx(ic-1))
c
        do j=1,jc-1
	  edgecy(j) = p5*(yc(j) + yc(j+1))
        enddo
        edgecy(0)  = yc(1) - (edgecy(1) - yc(1))
	edgecy(jc) = yc(jc) + (yc(jc) - edgecy(jc-1))
c
c       construct edges of "fine" grid boxes
c
        do i=1,if-1
	  edgefx(i) = p5*(xf(i) + xf(i+1))
        enddo
        edgefx(0)  = xf(1) - (edgefx(1) - xf(1))
	edgefx(if) = xf(if) + (xf(if) - edgefx(if-1))
c
        do j=1,jf-1
	  edgefy(j) = p5*(yf(j) + yf(j+1))
        enddo
        edgefy(0)  = yf(1) - (edgefy(1) - yf(1))
	edgefy(jf) = yf(jf) + (yf(jf) - edgefy(jf-1))
c
c       calculate "dx" and "dy" for the "fine" grid boxes
c
	do i=1,if
	  dx(i,1) = edgefx(i) - edgefx(i-1)       
	  dx(i,2) = dx(i,1)
        enddo
	dx(0,1) = dx(1,1)
	dx(0,2) = dx(1,2)
c
	do j=1,jf
	  dy(j,1) = edgefy(j) - edgefy(j-1)
	  dy(j,2) = dy(j,1)
        enddo
	dy(0,1) = dy(1,1)
	dy(0,2) = dy(1,2)
c
c       modify "dx" and "dy" for possibly partial "fine" grid boxes 
c       near the edges of each coarse grid box.
c       "ii" is the index of the fine grid box which contains the
c       eastern edge of coarse grid box with index "i".
c       dx(ii,1) is the portion of the fine grid box to the west of the
c       edge and dx(ii,2) is the portion to the east. similarly,
c       dy(jj,1) is to the south and dy(jj,2) is to the north of the
c       northern edge of coarse box with index "j".
c        
c       note: edgefx and edgefy are zero based and need the -1 when
c       using "indp" 
c
        do i=0,ic
	  ii   = indp (edgecx(i), edgefx, if+1) - 1
	  frac = abs(edgecx(i) - edgefx(ii))
	  if (edgefx(ii) .lt. edgecx(i)) then
	    ii = ii + 1
	    dx(ii,2) = (edgefx(min(if,ii)) - edgefx(ii-1)) - frac
	    dx(ii,1) = frac
	  else
	    dx(ii,2) = frac
	    dx(ii,1) = (edgefx(ii) - edgefx(max(ii-1,0))) - frac
	  endif
	  ieo(i) = min(if,max(1,ii))
        enddo
        do i=1,ic
	  iso(i) = max(1,ieo(i-1))
        enddo
        iso(0) = ieo(0)
c
        do j=0,jc
	  jj   = indp (edgecy(j), edgefy, jf+1) - 1
	  frac = abs(edgecy(j) - edgefy(jj))
	  if (edgefy(jj) .lt. edgecy(j)) then
	    jj = jj + 1
	    dy(jj,2) = (edgefy(min(jf,jj)) - edgefy(jj-1)) - frac
	    dy(jj,1) = frac
	  else
	    dy(jj,2) = frac
	    dy(jj,1) = (edgefy(jj) - edgefy(max(jj-1,0))) - frac
	  endif
	  jeo(j) = min(jf,max(1,jj))
        enddo
        do j=1,jc
	  jso(j) = max(1,jeo(j-1))
        enddo
        jso(0) = jeo(0)
c
c       store the weights into the "work" array
c
        indx = 1
        do j=0,jc
	  work(indx)   = jso(j)
	  work(indx+1) = jeo(j)
	  indx         = indx + 2
        enddo
c
        do i=0,ic
	  work(indx)   = iso(i)
	  work(indx+1) = ieo(i)
	  indx         = indx + 2
        enddo
c
        do j=0,jf
	  work(indx)   = dy(j,1)
	  work(indx+1) = dy(j,2)
	  indx         = indx + 2
        enddo
c
        do i=0,if
	  work(indx)   = dx(i,1)
	  work(indx+1) = dx(i,2)
	  indx         = indx + 2
        enddo
c
c       verify that coarse grid is coarser than the fine grid
c
        do j=jstart,jend
	  if ((jso(j) .eq. jso(j+1)) .and. (yc(j) .ge. yf(1))
     &      .and. (yc(j) .le. yf(jf))) then
            write (stdout,*)
     &         '=>Warning: "Coarse" grid is finer than "fine" grid'
     &,     ' near yf(',jso(j),') =',yf(jso(j))
     &,     ' (average may not be accurate)'
	  endif
        enddo
c
        do i=istart,iend
	  if ((iso(i) .eq. iso(i+1)) .and. (xc(i) .ge. xf(1))
     &      .and. (xc(i) .le. xf(if))) then
            write (stdout,*)
     &      '=>Warning: "Coarse" grid is finer than "fine" grid'
     &,     ' near xf(',iso(i),') = ',xf(iso(i))
     &,     ' (average may not be accurate)'
	  endif
        enddo
#ifdef skip_interp_details
        show_coord = .false.
c        print '(/a/)'
c     &,   'Remove -Dskip_interp_details to show the indices'
#else
        show_coord = .true.
        print '(/a/)'
     &,   'Use -Dskip_interp_details to not show the following'
#endif
        if (error .or. show_coord) then
          write (stdout,*)
     & ' Indices for averaging fine grid to coarse grid:'
          write (stdout,*)
     & ' (fractional grid boxes are accounted for)'
          write (stdout,8700)
          write (stdout,9000) (m,iso(m),ieo(m),m=istart,iend)
	  write (stdout,*) ' '
          write (stdout,*) ' Coordinates for coarse grid points "xc" ='
          write (stdout,8500) xc
          write (stdout,*) ' Coordinates for fine grid points "xf" ='
          write (stdout,8500) xf
c
          write (stdout,8800)
          write (stdout,9000) (m,jso(m),jeo(m),m=jstart,jend)
	  write (stdout,*) ' '
          write (stdout,*) ' Coordinates for coarse grid points "yc" ='
          write (stdout,8500) yc
          write (stdout,*) ' Coordinates for fine grid points "yf" ='
          write (stdout,8500) yf
        endif
        if (error) stop '=>ftc'
      else
c	write (stdout,*)
c     & '              (using previously initialized averaging weights)'
c        write (stdout,*) ' '
c
c       extract the weights from the "work" array
c
        indx = 1
        do j=0,jc
	  jso(j) = nint(work(indx))   
	  jeo(j) = nint(work(indx+1))
	  indx   = indx + 2
        enddo
c
        do i=0,ic
	  iso(i) = nint(work(indx))
	  ieo(i) = nint(work(indx+1))
	  indx   = indx + 2
        enddo
c
        do j=0,jf
	  dy(j,1) = work(indx)
	  dy(j,2) = work(indx+1)
	  indx    = indx + 2
        enddo
c
        do i=0,if
	  dx(i,1) = work(indx)
	  dx(i,2) = work(indx+1)
	  indx    = indx + 2
        enddo
      endif
c
c-----------------------------------------------------------------------
c     average the "fine" grid to the "coarse" grid
c-----------------------------------------------------------------------
c
      nd = 0              ! initialize counter
      do m=jstart,jend
        do i=istart,iend
          weight = c0
          sum    = c0
          do j=jso(m),jeo(m)
	    indy = 2
	    if (j .eq. jeo(m)) indy = 1
	    wty = dy(j,indy)
            do ii=iso(i),ieo(i)
	      indx = 2
	      if (ii .eq. ieo(i)) indx = 1
              if (f(ii,j).eq.spval) then
                area = c0
              else
                area = dx(ii,indx)*wty
              end if
              weight = weight + area
              sum    = sum + f(ii,j)*area
            enddo
          enddo
          if (weight .ne. c0) then 
            c(i,m) = sum/weight
          else
            c(i,m) = spval
            nd = nd+1
          end if
        enddo
      enddo
#ifndef skip_interp_details
      if (nd.ne.0) write(stdout,'(i8,a,e15.5,a/)')
     &     nd,' points set to special value ', spval, ' during ftc!'
#endif
      return
8500  format (1x,10g11.4)
8700  format (/' Along the 1st dimension, the form is (Coarse grid'
     &,' point "xc": range of fine grid points "xf" to average)'/)
8800  format (/' Along the 2nd dimension, the form is (Coarse grid'
     &,' point "yc": range of fine grid points "yf" to average)'/)
9000  format (5(1x,'(',i4,': ',i4,' to ',i4,')'),/)
      end subroutine ftc




      subroutine ctf(c, ic, jc, xc, yc, f, if, jf, istart, iend
     &,               jstart, jend, xf, yf, init, work, lenw,spval)
c
c=======================================================================
c
c     "ctf" is a nemonic for "coarse to fine".
c
c     obtain a fine grid representaion of a coarse grid dataset by 
c     linear interpolation of grid box values on the coarse grid to grid   
c     boxes on the fine grid.
c
c
c     inputs:
c
c     c      = coarse grid data
c     ic     = inner dimension of coarse grid "c"
c     jc     = outer dimension of coarse grid "c"
c     xc     = coordinates for inner dimension of "c" (eg: longitues)    
c     yc     = coordinates for outer dimension of "c" (eg: latitues)    
c
c     if     = inner dimension of "f"
c     jf     = outer dimension of "f"
c     xf     = coordinates for inner dimension of "f" (eg: longitues)
c     yf     = coordinates for outer dimension of "f" (eg: latitues)
c
c     istart = starting index along inner dimension of "f" for which
c              interpolated values are desired
c     iend   = ending index along inner dimension of "f" for which
c              interpolated values are desired
c     jstart = starting index along outer dimension of "f" for which
c              interpolated values are desired
c     jend   = ending index along outer dimension of "f" for which
c              interpolated values are desired
c     init   = initialize the interpolation factors
c              "init" should be set = 1 on the first call. 
c              "init" <> 1 uses the previously computed factors stored
c              in "work" array.
c     work   = work array of interpolation factors when "init" <> 1
c              (previously calculated by "ctf" when "init" = 1)
c     lenw   = size of work array. lenw should be >= 8*max(if,jf)
c     spval  = special value masking out dry points
c       
c     output:
c
c     f      = interplated data on fine grid defined over
c              ((f(i,j),i=istart,iend),j=jstart,jend)
c     work   = work array of interpolation factors when "init" = 1
c
c     restrictions:
c
c     fine and coarse grids are assumed rectangular with "xf" and "xc"
c     having the same units. "yf" and "yc" must also have the same units
c     the fine domain xf(istart)...xf(iend) must be within 
c     the coarse domain xc(1)...xc(ic). Similarly,   
c     and yc(js)...yc(je) must be within yf(1)...yf(jf). all coordinates
c     must be strictly monotonically increasing.
c 
c
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
!     modified:    J. Dengg, jdengg@awi-bremerhaven.de
!                  Sep-97
c=======================================================================
c
      integer :: stdout=6,indp
      logical error, show_coord
      parameter (len=10000, p5=0.5, c0=0.0)
      dimension indxi(len), indxj(len), dnorth(len), dsouth(len)
     &,         deast(len), dwest(len), width(len), height(len)
      dimension f(if,jf), xf(if), yf(jf)
      dimension c(ic,jc), xc(ic), yc(jc)
      dimension work(lenw)
      real,intent(in) :: spval
c
c-----------------------------------------------------------------------
c     initialize weights or use previously calculated weights
c-----------------------------------------------------------------------
c
      if (init .eq. 1) then
        error = .false.
c
c       test to verify that array sizes do not exceed limits
c
        if (if .gt. len .or. jf .gt. len) then
          i = max(if,jf)
          write (stdout,*) '=>Error: increase "len" in "ctf" to ',i
	  error = .true.
        endif
	if (lenw .lt. 8*max(if,jf)) then
          write (stdout,*) '=>Error: increase size of "work" array',
     &      ' to at least ',8*max(if,jf),' for calls to "ctf"'
	  error = .true.
        endif
c
c       verify that the "fine" grid lies within the "coarse" grid
c
        epsilon = 1.e-5
        xcminus = xc(1) - epsilon*(xc(2)-xc(1))
	xcplus = xc(ic) + epsilon*(xc(ic)-xc(ic-1))
        if (xf(istart) .lt. xcminus .or. xf(iend) .gt. xcplus) then
	  error = .true.
        write (stdout,*)
     &    '=>Warning: "fine" grid outside "coarse" grid in "ctf".'
          if (xf(istart) .lt. xc(1)) then
	    write (stdout,*) '  xf(',istart,')  .lt.  xc(1)'  
	  endif
          if (xc(ic) .lt. xf(iend)) then
	    write (stdout,*) '  xf(',iend,')  .gt.  xc(',ic,')'  
	  endif
        endif 
        ycminus = yc(1) - epsilon*(yc(2)-yc(1))
	ycplus = yc(jc) + epsilon*(yc(jc)-yc(jc-1))
        if (ycminus .gt. yf(jstart) .or. ycplus .lt. yf(jend)) then
	  error = .true.
          write (stdout,*)
     &    '=>Warning: "fine" grid outside "coarse" grid in "ctf".'
          if (yc(1) .gt. yf(jstart)) then
	    write (stdout,*) '  yf(',jstart,') .lt. yc(1)' 
	  endif
          if (yc(jc) .lt. yf(jend)) then
	    write (stdout,*) '  yf(',jend,')  .gt.  yc(',jc,')'  
	  endif
        endif
c
c       find interpolation factors 
c
        indx = 1
        do j=jstart,jend
	  jj = indp (yf(j), yc, jc)
	  if (yc(jj) .gt. yf(j) .or. jj .eq. jc) jj = jj - 1
	  indxj(j) = jj
	  dnorth(j) = yc(jj+1) - yf(j)
	  dsouth(j) = yf(j) - yc(jj)
	  height(j) = yc(jj+1) - yc(jj)
c
c         store into "work" array for future use (when "init" <> 1)
c
          work(indx)   = indxj(j)
	  work(indx+1) = dnorth(j)
	  work(indx+2) = dsouth(j)
	  work(indx+3) = height(j)
	  indx         = indx + 4
        enddo
c
        do i=istart,iend
	  ii = indp (xf(i), xc, ic)
	  if (xc(ii) .gt. xf(i) .or. ii .eq. ic) ii = ii - 1
	  indxi(i) = ii
	  deast(i) = xc(ii+1) - xf(i)
	  dwest(i) = xf(i) - xc(ii)
	  width(i) = xc(ii+1) - xc(ii)
c
c         store into "work" array for future use (when "init" <> 1)
c
          work(indx)   = indxi(i)
	  work(indx+1) = deast(i)
	  work(indx+2) = dwest(i)
	  work(indx+3) = width(i)
	  indx         = indx + 4
        enddo
#ifdef skip_interp_details
        show_coord = .false.
c        print '(/a/)'
c     &,   'Remove -Dskip_interp_details to show the coordinates'
#else
        show_coord = .true.
        print '(/a/)'
     &,   'Use -Dskip_interp_details to not show the following'
#endif
        if (error .or. show_coord) then
	  write (stdout,*) ' '
          write (stdout,*) ' Coordinates for coarse grid points "xc" ='
          write (stdout,8500) xc
	  write (stdout,*) ' Coordinates for fine grid points "xf" ='
          write (stdout,8500) xf

          write (stdout,*) ' Coordinates for coarse grid points "yc" ='
          write (stdout,8500) yc
          write (stdout,*) ' Coordinates for fine grid points "yf" ='
          write (stdout,8500) yf
        endif
	if (error) stop '=>ctf'
      else
c
c       extract previously calculated interpolation weights from "work"
c
        indx = 1
        do j=jstart,jend
          indxj(j)  = nint(work(indx))
	  dnorth(j) = work(indx+1)
	  dsouth(j) = work(indx+2)
	  height(j) = work(indx+3)
	  indx      = indx + 4
        enddo
c
        do i=istart,iend 
          indxi(i) = nint(work(indx))
	  deast(i) = work(indx+1)
	  dwest(i) = work(indx+2)
	  width(i) = work(indx+3)
	  indx     = indx + 4
        enddo
      endif
c
c-----------------------------------------------------------------------
c     interpolate data from "coarse" to "fine" grid
c-----------------------------------------------------------------------
c
      nd = 0
      do jj=jstart,jend
        j = indxj(jj)
        do ii=istart,iend
	  i = indxi(ii)
          if (c(i,j) .eq. spval) then 
             w1=c0; w11=c0
          else
             w1=c(i,j)*deast(ii)*dnorth(jj)
            w11=deast(ii)*dnorth(jj)
          end if 
          if (c(i+1,j) .eq. spval) then
             w2=c0; w22=c0
          else
             w2=c(i+1,j)*dwest(ii)*dnorth(jj)
            w22=dwest(ii)*dnorth(jj)
          end if
          if (c(i,j+1) .eq. spval) then
             w3=c0; w33=c0
          else
             w3=c(i,j+1)*deast(ii)*dsouth(jj)
            w33=deast(ii)*dsouth(jj)
          end if
          if (c(i+1,j+1) .eq. spval) then
             w4=c0; w44=c0
          else
             w4=c(i+1,j+1)*dwest(ii)*dsouth(jj)
            w44=dwest(ii)*dsouth(jj)
          end if
          weight = w11+w22+w33+w44
          if (weight .ne. c0) then
            f(ii,jj)= (w1+w2+w3+w4)/weight   
          else
            f(ii,jj) = spval
            nd = nd+1
          end if
        enddo
      enddo
#ifndef skip_interp_details
      if (nd.ne.0) write(stdout,'(i8,a,e15.5,a/)')
     &     nd,' land points set to special value ', spval,
     &     ' during ctf!'
#endif
c
8500  format (1x,10g11.3)
      end subroutine ctf


      subroutine rotate (glt, gln, phir, thetar, psir, rlt, rln)
c
c=======================================================================
c     subroutine rotate takes a geographic latitude and longitude and 
c     finds the the equivalent latitude and longitude on a rotated grid.
c     when going from a geographic grid to a rotated grid, all of the 
c     defined rotation angles given to rotate by the calling program 
c     are positive, but when going from a rotated grid back to the 
c     geographic grid, the calling program must reverse the angle order 
c     (phir and psir are switched) and all of the angles made negative.
c
c     the first rotation angle phir is defined as a rotation about the
c     original z axis. the second rotation angle thetar is defined as a
c     rotation about the new x axis. the final rotation angle psir is
c     defined as a rotation about the new z axis. these rotation angles
c     are just the Euler angles as defined in "classical mechanics"
c     Goldstein (1951).
c
c     author:   M. Eby            e-mail eby@uvic.ca
c=======================================================================
c
c     g...  = geographic value
c     r...  = rotated value
c     ...lt = latitude (or equivalent spherical coordinate)
c     ...ln = longitude (or equivalent spherical coordinate)
c     ...x  = x coordinate
c     ...y  = y coordinate
c     ...z  = z coordinate
c     psir, thetar, phir = Euler angles defining rotation
c
c     define rad for conversion to radians.
      implicit none
      real (kind=8) glt, gln, phir, thetar, psir, rlt, rln
      real (kind=8) rad,thetas,phis,gx,gy,gz,rx,ry,rz
      real (kind=8) :: c1=1.
c
      rad = acos(-1.)/180.
c
c     convert latitude and longitude to spherical coordinates
      thetas = gln
      if (thetas .gt. 180.) thetas = thetas - 360.
      if (thetas .lt. -180.) thetas = thetas + 360.
      phis = (90. - glt)*rad
      thetas = thetas*rad
c
c     translate point into Cartesian coordinates for rotation.
      gx = sin(phis)*cos(thetas)
      gy = sin(phis)*sin(thetas)
      gz = cos(phis)
c
c     rotate the point (gx, gy, gz) about the z axis by phir then the x
c     axis by thetar and finally about the z axis by psir.
c 
      rx = gx*(cos(psir)*cos(phir) - cos(thetar)*sin(phir)*sin(psir)) +
     &     gy*(cos(psir)*sin(phir) + cos(thetar)*cos(phir)*sin(psir)) +
     &     gz*sin(psir)*sin(thetar)
c
      ry = gx*(-sin(psir)*cos(phir) - cos(thetar)*sin(phir)*cos(psir)) +
     &     gy*(-sin(psir)*sin(phir) + cos(thetar)*cos(phir)*cos(psir)) +
     &     gz*(cos(psir)*sin(thetar))
c
      rz = gx*(sin(thetar)*sin(phir)) + gy*(-sin(thetar)*cos(phir)) +
     &     gz*(cos(thetar))
c
c     convert rotated point back to spherical coordinates
c
c     check for rounding error (arccos(x): abs(x) must be .le. 1)
      rz = min(rz, c1)
      rz = max(rz, -c1)
      rlt = acos(rz)
c     if point is at a pole set rotated longitude equal to initial.
      if (rlt .le. 0. .or. rlt .ge. 180.*rad) then
        rln = thetas
      else
c     if rln lies between -135 and -45 or between 45 and 135 degrees
c     it is more accurate to use an arccos calculation.
        if (abs(rx/sin(rlt)) .lt. cos(45.*rad)) then
          rln = rx/sin(rlt)
c     check for rounding error (arccos(x): abs(x) must be .le. 1)
          rln = min(rln, c1)
          rln = max(rln, -c1)
          rln = acos(rln)
c     arccos will give rln between 0 and 180 degrees.  if the point
c     is negative in y, rln must be equal to negative rln.
          if (ry .lt. 0.) rln = -rln
        else
c     if rln lies between -45 and 45 or between 135 and -135 degrees
c     it is more accurate to use an arcsin calculation.
          rln = ry/sin(rlt)
c     check for rounding error (arcsin(x): abs(x) must be .le. 1)
          rln = min(rln, c1)
          rln = max(rln, -c1)
          rln = asin(rln)
c     arcsin will give rln between -90 and 90 degrees. if the point
c     is negative in x, rln must be equal to 180 degrees minus rln.
          if (rx .lt. 0.) rln = 180.*rad - rln
        endif
      endif
c
c     convert back to degrees of latitude and longitude.
      rlt = 90. - rlt/rad
      rln = rln/rad
      if (rln .gt. 180.) rln = rln - 360.
      if (rln .le. -180.) rln = rln + 360.
      end subroutine rotate






      subroutine rot_intrp_sclr(g_in, xg_in, yg, ig, jg, r, xr, yr
     &,        ir, jr, psir, thetar, phir,spval,verbose)
c
c=======================================================================
c     interpolate scaler data from an geographic data grid to a 
c     rotated model grid
c
c     input
c     psir, thetar, phir = Euler angles defining rotation
c     g  = scaler on geographic data grid
c     xg = longitude of data points on geographic data grid
c     yg = latitude of data points on geographic data grid
c     ig = number of longitudes in on geographic data grid 
c     jg = number of latitudes in on geographic data grid 
c     xr = longitude of points on rotated model grid
c     yr = latitude of points on rotated model grid
c     ir = number of longitudes in rotated model grid 
c     jr = number of latitudes in rotated model grid 
c
c     output
c     r  = scaler on rotated model grid
c
c     internal
c     (rln,rlt) = (longitude,latitude) in rotated coordinates
c     (gln,glt) = (longitude,latitude) in geographic coordinates
c     xg(iw) = point on the geographic grid to the west of (gln,glt)
c     xg(ie) = point on the geographic grid to the east of (gln,glt)
c     yg(js) = point on the geographic grid to the south of (gln,glt)
c     yg(jn) = point on the geographic grid to the north of (gln,glt)
c
c=======================================================================
c
      implicit none
      integer ir,jr,ig,jg
      real g_in(ig,jg), xg_in(ig)
c      real ::  g(ig,jg), xg(ig) 
      real,allocatable ::  g(:,:), xg(:) 
      real yg(jg), r(ir,jr), xr(ir), yr(jr)
      real (kind=8) psir,thetar,phir
      real ::  spval
      logical :: verbose

      real (kind=8) epsln,glt_min,glt_max,gln_min,gln_max,ln_err,lt_err
      integer :: stdout = 6,indp
      integer istrt,i,iend,iw,ie
      integer jstrt,j,jend,js,jn
      real (kind=8) del,wtw,wte,wts,wtn
      real (kind=8) glt,gln,rlt,rln
      real  glt4,gln4
      real (kind=8) wtne,wtse,wtnw,wtsw,wall

      allocate( g(ig,jg), xg(ig) )
c      print*,'ok'

c     make sure that longitudes begin around 0deg

      if (xg_in(ig)>360.) then
       j=1; do i=1,ig
        if (xg_in(i)<=360.)  j=i
       enddo 
       if (verbose) then
         print*,' have to shift geographic data by ',j,' points'
         print*,' assuming cyclic boundary conditions'
       endif
       g(1:ig-(j+1)+1,:)= g_in(j+1:ig,:)
       g(ig-(j+1)+2:ig,:)=g_in(1:j,:) 
       xg(1:ig-(j+1)+1)= xg_in(j+1:ig)-360. 
       xg(ig-(j+1)+2:ig)=xg_in(1:j) 
      else
        xg=xg_in; g=g_in
      endif

c
      epsln = 1.e-10; 
      glt_min = 90.; glt_max = -90.; gln_min = 360.
      gln_max = -360.; ln_err = 0; lt_err = 0
c
c     find longitude points of data within interval [0., 360.]
      istrt = 1; do i=2,ig
        if (xg(i-1) .lt. 0. .and. xg(i) .ge. 0.) istrt = i
      enddo 
      iend = ig; do i=2,ig
        if (xg(i-1) .lt. 360. .and. xg(i) .ge. 360.) iend = i
      enddo
c
c     find latitude points of data within interval [-90., 90.]
      jstrt = 1; do j=2,jg
        if (yg(j-1) .lt. -90. .and. yg(j) .ge. -90.) jstrt = j
      enddo 
      jend = jg; do j=2,jg
        if (yg(j-1) .lt. 90. .and. yg(j) .ge. 90.) jend = j
      enddo
c
c     interpolate data to model grid 
c
      do j=1,jr
        do i=1,ir
          rlt=yr(j);rln=xr(i)
          call rotate (rlt, rln, -psir, -thetar, -phir, glt, gln)
          if (gln .lt. 0.) gln = gln + 360.
          if (gln .ge. 360.) gln = gln - 360.
          glt_min = min(glt,glt_min); glt_max = max(glt,glt_max)
          gln_min = min(gln,gln_min); gln_max = max(gln,gln_max)
          gln4=gln
          iw = indp (gln4, xg(istrt), iend-istrt+1)+istrt-1
          if (xg(iw) .gt. gln) iw = iw - 1
	  ie = iw + 1
	  if (iw .ge. istrt .and. ie .le. iend) then
	    del = xg(ie) - xg(iw); wtw = (xg(ie) - gln)/del
	  else
c     east or west of the last data value. this could be because a
c     cyclic condition is needed or the dataset is too small. in either 
c     case apply a cyclic condition
            ln_err = 1; iw = iend; ie = istrt
	    del = xg(ie) + 360. + epsln - xg(iw) 
            if (xg(ie) .ge. gln) then
  	      wtw = (xg(ie) - gln)/del
            else
	      wtw = (xg(ie) + 360. + epsln - gln)/del
            endif
	  endif
	  wte = 1. - wtw
          glt4=glt
          js = indp (glt4, yg(jstrt), jend-jstrt+1)+jstrt-1
	  if (yg(js) .gt. glt) js = max(js - 1,jstrt)
	  jn = min(js + 1,jend)
	  if (yg(jn) .ne. yg(js) .and. yg(js) .le. glt) then
            wts = (yg(jn) - glt)/(yg(jn) - yg(js))
          else
c     north or south of the last data value. this could be because a
c     pole is not included in the data set or the dataset is too small.
c     in either case extrapolate north or south
            lt_err = 1; wts = 1.
          endif
	  wtn = 1. - wts
          if (g(ie,jn)==spval.or.g(iw,jn)==spval.or.
     &        g(ie,js)==spval.or.g(iw,js)==spval) then
c          care about spvals
c          first approach, leaves many points with sp values
c           r(i,j)=spval
c          second approach
           wtne=wtn*wte;if (g(ie,jn)==spval) wtne=0.
           wtse=wts*wte;if (g(ie,js)==spval) wtse=0.
           wtnw=wtn*wtw;if (g(iw,jn)==spval) wtnw=0.
           wtsw=wts*wtw;if (g(iw,js)==spval) wtsw=0.
           wall = wtne+wtse+wtnw+wtsw
           if (wall /= 0.) then
            r(i,j) = (g(ie,jn)*wtne + g(ie,js)*wtse
     &              + g(iw,jn)*wtnw + g(iw,js)*wtsw)/wall
           else
            r(i,j)=spval
           endif
          else
           r(i,j) = g(ie,jn)*wte*wtn + g(ie,js)*wte*wts
     &            + g(iw,jn)*wtw*wtn + g(iw,js)*wtw*wts
          endif
        enddo
      enddo
c
      if (ln_err .eq. 1 .and. verbose) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: the geographic data set does not extend far   '
     &,   '             enough east or west - a cyclic boundary       '
     &,   '             condition was applied. check if appropriate   '
        write (stdout,'(/,(1x,a,2f8.2))')
     &    '    data required between longitudes:', gln_min, gln_max
     &,   '      data set is between longitudes:', xg(istrt), xg(iend)
c        stop
      endif
c
      if (lt_err .eq. 1 .and. verbose) then
        write (stdout,'(/,(1x,a))')
     &    '==> Warning: the geographic data set does not extend far   '
     &,   '             enough north or south - extrapolation from    '
     &,   '             the nearest data was applied. this may create '
     &,   '             artificial gradients near a geographic pole   ' 
        write (stdout,'(/,(1x,a,2f8.2))')
     &    '    data required between latitudes:', glt_min, glt_max
     &,   '      data set is between latitudes:', yg(jstrt), yg(jend)
c        stop
      endif
      deallocate( g, xg )
      end subroutine rot_intrp_sclr




      subroutine rot_intrp_vctr (g, xg, yg, ig, jg, r, xr, yr
     &,           ir, jr, psir, thetar, phir,spval,verbose)
c
c=======================================================================
c     interpolate vector data from an geographic data grid to a
c     rotated model grid
c
c     input
c     psir, thetar, phir = Euler angles defining rotation
c     g  = vector on geographic data grid
c     xg = longitude of data points on geographic data grid
c     yg = latitude of data points on geographic data grid
c     ig = number of longitudes in geographic data grid 
c     jg = number of latitudes in geographic data grid 
c     xr = longitude of points on rotated model grid
c     yr = latitude of points on rotated model grid
c     ir = number of longitudes in rotated model grid 
c     jr = number of latitudes in rotated model grid 
c
c     output
c     r  = vector on rotated model grid
c
c     internal
c     (rln,rlt) = (longitude,latitude) in rotated coordinates
c     (gln,glt) = (longitude,latitude) in geographic coordinates
c     xg(iw) = point on the geographic grid to the west of (gln,glt)
c     xg(ie) = point on the geographic grid to the east of (gln,glt)
c     yg(js) = point on the geographic grid to the south of (gln,glt)
c     yg(jn) = point on the geographic grid to the north of (gln,glt)
c
c=======================================================================
c
      implicit none
      integer ig,jg,ir,jr
      real g(ig,jg,2), xg(ig), yg(jg), r(ir,jr,2), xr(ir), yr(jr)
      real (kind=8) psir,thetar,phir
      real spval
      logical verbose

      integer :: stdout = 6,i,j
      real (kind=8) rad,vmag,a,glt,gln,angle,rln,rlt
      real (kind=8) :: c1=1.
c
      rad = acos(-1.)/180.

!      print*,' rot_intrp_vec'
c
c     interpolate vector components as scalers on rotated model grid 
c
      call rot_intrp_sclr (g(1,1,1), xg, yg, ig, jg, r(1,1,1), xr, yr
     &,                    ir, jr, psir, thetar, phir,spval,verbose)
      call rot_intrp_sclr (g(1,1,2), xg, yg, ig, jg, r(1,1,2), xr, yr
     &,                    ir, jr, psir, thetar, phir,spval,verbose)
c
c     correct vector direction
c
      do j=1,jr
        do i=1,ir
         if (r(i,j,1)/=spval.and.r(i,j,2)/=spval) then
          vmag = sqrt(r(i,j,1)**2 + r(i,j,2)**2)
          if (vmag >  0.) then 
            a = r(i,j,1)/vmag
            a = min(a, c1)
            a = max(a, -c1)
            a = acos(a)
            if (r(i,j,2) < 0.) a = -a
            rlt=yr(j);rln=xr(i)
            call rotate(rlt, rln, -psir, -thetar, -phir, glt, gln)
            call rotvec(glt, gln, phir, thetar, psir, angle)
            a = a + angle*rad
            r(i,j,1) = vmag*cos(a)
            r(i,j,2) = vmag*sin(a)
          else
	    r(i,j,1) = 0.
	    r(i,j,2) = 0.
          endif
         endif
        enddo
      enddo
!      print*,' rot_intrp_vec done'
      end subroutine rot_intrp_vctr


      subroutine rotvec (glt, gln, phir, thetar, psir, angle)
c
c=======================================================================
c     subroutine rotvec takes a geographic latitude and longitude and 
c     finds the the vector rotation angle angle (in degrees) for a 
c     vector on the rotated grid. when going from the geographic to a 
c     rotated grid, all of the defined rotation angles given to rotvec 
c     by the calling program are positive, but when going from a 
c     rotated grid back to the geographic grid, the calling program 
c     must reverse the angle order (phir and psir are switched) and all 
c     of the angles made negative. if a pole is detected then an angle 
c     of zero is returned.
c
c     rotvec rotates the point defining the head of a very short north
c     or south pointing geographic vector. the angle between this 
c     vector and a similar direction vector defined in the new grid is 
c     calculated using the law of cosines. the accuracy of this 
c     calculation depends on the size of the direction vector (delta) 
c     and the precision of the computation. the smaller the vector the 
c     more accurate the calculation but then the more precision 
c     required. double precision is strongly recommended.
c
c     author:   M. Eby            e-mail eby@uvic.ca
c=======================================================================
c
c     g...   = geographic value
c     r...   = rotated value
c     ...lt  = latitude
c     ...ln  = longitude
c     ...lth = latitude of head of vector
c     ...lnh = longitude of head of vector
c     dst    = distance between heads of vectors
c     delta  = length of vector
c     angle  = angle to rotate vectors back to original orientation
c     psir, thetar, phir = Euler angles defining rotation
c
c     define multiplier rad for conversion to radians.
      implicit none
      real (kind=8) rad,t,rlnh,rlth, delta
      real (kind=8) glt,gln,phir,thetar,psir,rlt
      real (kind=8) rln,glth,dst,angle,glnh
      real (kind=8) :: c1=1.

      rad = acos(-1.)/180.
c
c     define length of direction vector (single precision may require a 
c     longer and thus less accurate vector length).
      delta = 0.001
c
c     if the base is in the north of the geographic grid use a south 
c     pointing vector to avoid any possibility of going over the pole.
      if (glt .ge. 0.) delta = -delta
c
c     find the base of the direction vectors in the rotated grid.
      call rotate (glt, gln, phir, thetar, psir, rlt, rln)
c
c     if base in the rotated grid is near a pole return an angle of zero
      if (abs(rlt) .ge. 90.-abs(delta)) then
        angle = 0.
        return
      endif
c
c     find the head of the geographic grid direction vector in the 
c     rotated grid.
      call rotate (glt+delta, gln, phir, thetar, psir, glth, glnh)
c
c     if the base is in opposite hemispheres switch the vector
c     direction for better accuracy.
      if (glt*rlt .lt. 0) delta = -delta
c
c     find the head of the rotated grid direction vector.
      rlth = rlt + delta
      rlnh = rln
c
c     find the distance between the heads of the direction vectors.
      call dist (glth, glnh, rlth, rlnh, dst)
c
c     find the angle between direction vectors with the law of cosines.
      delta = abs(delta)
      angle = (cos(dst)-cos(delta)**2)/(sin(delta)**2)
      angle = min(angle, c1)
      angle = max(angle, -c1)
      angle = acos(angle)/rad
      t = abs(delta)
c
c     adjust the angle if the direction vectors are opposite.
      if (glt*rlt .lt. 0) angle = 180. - angle
c     determine the sign of the angle by checking the offset longitudes.
      if (glnh+360. .gt. rlnh+360.) angle = -angle
c     change sign if the original direction vector was pointing south.
      if (glt .ge. 0.) angle = -angle
      end subroutine rotvec


      subroutine dist (lat1, lng1, lat2, lng2, dst)
c
c=======================================================================
c     subroutine dist calculates the arc distance between two 
c     points given their latitudes and longitudes
c=======================================================================
c
      implicit none
      real (kind=8) lat1, lng1, lat2, lng2, dst, rad
      real (kind=8) :: c1=1.,c90=90.
c
c     define multiplier rad for conversion to radians.
      rad = acos(-1.)/180.
c
c     check input.
      lat1 = min(lat1, c90)
      lat1 = max(lat1, -c90)
      if (lng1 .lt. 0.) lng1 = lng1 + 360.
      if (lng1 .gt. 360.) lng1 = lng1 - 360.
      lat2 = min(lat2, c90)
      lat2 = max(lat2, -c90)
      if (lng2 .lt. 0.) lng2 = lng2 + 360.
      if (lng2 .gt .360.) lng2 = lng2 - 360.
c
      dst = sin(lat1*rad)*sin(lat2*rad)+cos(lat1*rad)*
     &      cos(lat2*rad)*cos((lng1-lng2)*rad)
      dst = min(dst, c1)
      dst = max(dst, -c1)
      dst = (acos(dst)/rad)
      end subroutine dist



      subroutine tranlon (c, ic, il, jl, cx, fx, ifl, tx)
c
c-----------------------------------------------------------------------
c     translate longitudes  "cx" to "tx" so that tx(i) i=1..ic
c     completely encloses model longitudes  fx(i) i=1..ifl
c     note that "tx" may extend beyond 360 degrees to contain "fx".
c     use same mapping to translate data in "c"
c
c     input:
c     c  = original data array
c     cx = original data longitudes
c     tx = translated data longitudes
c     fx = model longitudes
c
c     output
c     c  = translated data array
c-----------------------------------------------------------------------
c
      implicit none
      integer :: ic,il,jl,ifl
      real ::  c(ic,jl), t(ic), tx(ic), cx(ic), fx(ifl)
      integer :: indp,i,j,iw,im1
c
c-----------------------------------------------------------------------
c     find the index of the 1st model grid point on the data grid
c-----------------------------------------------------------------------
c
      iw = indp (fx(1), cx, ic)
      if (cx(iw) .gt. fx(1)) iw = max(1,iw-1)
c
c-----------------------------------------------------------------------
c     translate data longitudes so that tx(1) = cx(iw), tx(2) = cx(iw+1)
c-----------------------------------------------------------------------
c
      do i=1,ic
        tx(i) = cx(mod(i+iw-2,il) + 1)
        im1   = max(1,i-1)
        if (tx(i) .lt. tx(im1)) tx(i) = tx(i) + 360.0
      enddo
      if (fx(ifl) .gt. tx(ic)) then
        write (6,997) iw, ic, ifl
        write (6,998) 'tx= ',(tx(i),i=1,ic)
        write (6,998) 'fx= ',(fx(i),i=1,ifl)
        stop
      endif
c
c-----------------------------------------------------------------------
c     translate data to match translated longitudes
c-----------------------------------------------------------------------
c
      do j=1,jl
        do i=1,ic
          t(i) = c(mod(i+iw-2,il) + 1,j)
        enddo
        do i=1,ic
          c(i,j) = t(i)
        enddo
      enddo
c
      return
997   format (1x, ' ===>  tx(ic) < fx(ifl) in tranlon. iw=',i6,
     1         ' il=',i6,' ifl=',i6)
998   format (1x,a4,(5x,10e11.4))
      end subroutine tranlon




      subroutine fillgaps(imt, jmt, km, kloc, kmsk, 
     &                    dat, spval, ldoc,stop_on_error)

!------------------------------------------------------------------------
!
!     check a 2-dimensional data field dat(imt,jmt) for data  
!     gaps (marked by special value spval) at wet points (kmsk>kloc, 
!     i.e. where the bottom is deeper than the local k); fill gaps with 
!     arithmetic mean of available points among the 8 nearest horizontal 
!     neighbors, if possible. If no neighbors are available, retain the 
!     gap and repeat the process in the hope that larger gaps will be 
!     filled consecutively from the sides.
!
!     Mask out land/topography with special values after processing 
!     the data field.
!
!     At the moment, this routine is only intended for horizontal 
!     fields! Only maxcount=5000 gaps are allowed, and the routine will 
!     iterate the fill procedure only maxdo=500 times. Any remaining gaps
!     will be reported to the user.
!
!     Variable ldoc determines the level of output documentation:
!     for ldoc.ge.4, the locations of all gaps will be reported. 
!     For ldoc.ge.2, locations of the remaining gaps will be reported.
!
!     Input:  imt, jmt, km, kloc, kmsk, dat, spval, ldoc
!     Output: dat
!
!     J. Dengg, jdengg@awi_bremerhaven.de
!     Nov-97
!
!------------------------------------------------------------------------

      integer :: stdout = 6
      integer,intent(in) :: imt, jmt, km           ! field dimensions
      integer,intent(in) :: kloc                   ! index vertical position

      integer,parameter :: maxcount = 50000
      integer,parameter :: maxdo = 500             ! max. no. iterations

      integer :: igap(maxcount), jgap(maxcount)    ! gap coordinates

      integer,intent(in) :: kmsk(imt,jmt)          ! kmt- or kmu- field
      real,intent(inout) :: dat(imt,jmt)           ! data
      real,intent(in)    :: spval                  ! special value for 
!                                                  !         data gaps
      integer,intent(in) :: ldoc                   ! level of documentation
      logical,intent(in) :: stop_on_error          ! stop if gaps remain

      integer :: ngaps                             ! no. of initial gaps
      integer :: icount                            ! no. of remaining gaps


!     Security Check
!     --------------

      if (kloc > km) then
        write (stdout,*) 'Illegal k-index. ',
     &                   'Error in subroutine fillgaps.' 
        stop '==> fillgaps'
      end if



!     Count gaps
!     ----------

      icount = 0                                 ! initialize counters
        igap = 0
        jgap = 0
       ngaps = 0

      do j=1, jmt
        do i=1, imt
          if (kmsk(i,j)>=kloc .and. dat(i,j)==spval) then
            icount = icount + 1
            igap(icount) = i                     ! register position
            jgap(icount) = j
            if (ldoc.ge.4) write (stdout,*) 'gap #', icount,
     &         'at i=', igap(icount), ', j=', jgap(icount),
     &         ', k=', kloc
          end if
          if (icount>=maxcount) then
            write (stdout,*) 'Sorry, too many data gaps! ',
     &                       'Adjust dimension maxcount !'
            stop '==> fillgaps'
          end if
        end do
      end do

      ngaps = icount
      if (ldoc>=1) write (stdout,*) 'Data gaps: ', ngaps



!     Fill gaps with average of nearest neighbours
!     --------------------------------------------
!           ('maxdo' iterations only)

!        use all 8 adjacent grid points, unless gap is at boundary

      if (ngaps /= 0) then

        do n=1, maxdo                ! repeat process several times
          do ng=1, ngaps

            if (igap(ng)/=0 .and. jgap(ng)/=0) then

              is = max(igap(ng)-1, 1)
              ie = min(igap(ng)+1, imt)
              js = max(jgap(ng)-1, 1)
              je = min(jgap(ng)+1, jmt)
               sum = 0.
              isum = 0

              do j=js, je
                do i=is, ie
                  if (dat(i,j)/=spval) then
                     sum = sum + dat(i,j)
                    isum = isum + 1
                  end if
                end do
              end do

              if (isum/=0) then
                dat(igap(ng),jgap(ng)) = 1.*(sum/isum)
                igap(ng) = 0
                jgap(ng) = 0
                icount = icount-1
              end if

            end if        ! if (igap(ng)/=0 .and. jgap(ng)/=0) 
          end do          ! ng-loop

          if (icount==0) then
            if (ldoc>=1) 
     &      write (stdout,'(a/)') '           all gaps eliminated'
            goto 900
          end if

        end do            ! n-loop

        if (ldoc>=1) 
     &  write (stdout,*) '            ', icount, ' gaps remaining'

        if (ldoc >= 2) then                   ! report remaining gaps 
          do ng=1, ngaps
            if (igap(ng)/=0 .and. jgap(ng)/=0) then
              write (stdout,*) 'i=', igap(ng), ', j=', jgap(ng),
     &                       ', k=', kloc
            end if
          end do
        end if

        if (stop_on_error) stop

      end if              ! if (ngaps /= 0)



!     Assign special value to land/topography points
!     ----------------------------------------------

  900 where (kmsk < kloc) dat = spval

      end subroutine fillgaps


      subroutine fillgaps_poisson(imt, jmt, km, kloc, kmsk, 
     &                    dat, spval, ldoc,stop_on_error)

!------------------------------------------------------------------------
      implicit none

      integer :: stdout = 6
      integer,intent(in) :: imt, jmt, km           ! field dimensions
      integer,intent(in) :: kloc                   ! index vertical position

      integer,intent(in) :: kmsk(imt,jmt)          ! kmt- or kmu- field
      real,intent(inout) :: dat(imt,jmt)           ! data
      real,intent(in)    :: spval                  ! special value for 
!                                                  !         data gaps
      integer,intent(in) :: ldoc                   ! level of documentation
      logical,intent(in) :: stop_on_error          ! stop if gaps remain

      integer :: i,j,itt,max_itts=100
      integer, dimension(imt,jmt) :: ip,im,jp,jm,count
      real :: resid,omega,anorm,anormf,eps_sor=0.1
      integer :: mask(imt,jmt)

c    solve nabla^2 dat = 0 with neumann boundary conditions
c    and dx=dy=1

      where (kmsk < kloc) dat = spval
      mask = 0
      where (dat == spval .and. kmsk >= kloc ) mask=1
      where (mask == 1) dat = 0

      ip=0;im=0;jp=0;jm=0
      do i=2,imt-1
       do j=2,jmt-1
         if (dat(i+1,j) /=spval ) ip(i,j) = 1
         if (dat(i-1,j) /=spval ) im(i,j) = 1
         if (dat(i,j+1) /=spval ) jp(i,j) = 1
         if (dat(i,j-1) /=spval ) jm(i,j) = 1
       enddo
      enddo
      count=ip+jp+im+jm
      where (mask == 0) count = 1
      where (count == 0 .and. mask == 1) count = 1
      where (count == 0) 
        mask = 0
        count = 1
      end where
      
      omega=1.000001
      anormf=0
      anorm=1
      itt = 0
      do while (anorm>anormf*eps_sor .and. itt < max_itts ) 
       anorm=0
       do i=2,imt-1
        do j=2,jmt-1
         resid=dat(i+1,j)*ip(i,j)+dat(i-1,j)*im(i,j)+
     &         dat(i,j+1)*jp(i,j)+dat(i,j-1)*jm(i,j)-count(i,j)*dat(i,j)
         dat(i,j)=dat(i,j)+omega*resid/count(i,j)*mask(i,j)
         anorm=anorm + abs(resid)*mask(i,j)
        enddo
       enddo
       if (anormf==0) anormf=anorm
       print*,' itt=',itt,' anorm = ',anorm
        call sub_flush(6)
       itt=itt+1
      enddo

      if (itt >= max_itts) then
        print*,' solver did not converge'
        call sub_flush(6)
        if (stop_on_error) stop
      endif

      print*,' itt = ',itt
      call sub_flush(6)

      end subroutine fillgaps_poisson



      subroutine fillholes(imt, jmt, km, jloc, kmsk, 
     &                     dat, spval, ldoc, stop_on_error)

!------------------------------------------------------------------------
!
!     Check a 2-dimensional data field dat(imt,km) for data  
!     holes (marked by special value spval) at wet points (kmsk>k, 
!     i.e. where the bottom is deeper than the local k); fill holes  
!     consecutively with the neighbour values from above.
!
!     This leads to labile layered water columns!
!
!     After fillgaps there may be data gaps left, where the topography
!     contain deep holes and the interpolated data set provides no data.
!     A purely horizontal extrapolation does not fill these holes.
!
!     Does not mask out land/topography with special values. This must
!     be done before. It is assumed here that another routine (fillgaps)
!     has seen the data, before you decide to fill data holes in this
!     desperate way.
!
!     At the moment, this routine is only intended for vertical 
!     fields! Before you adapt this routine for horizontal fields
!     be aware that you will get uniform distributed values. 
!
!     Variable ldoc determines the level of output documentation:
!     In any case, if there are remaining holes this number will be reported.
!     for ldoc.ge.1, the locations of remaining holes will be reported.
!     for ldoc.ge.2, the number of initial holes will be reported.
!     for ldoc.ge.3, the locations of initial holes will be reported.
!     for ldoc.ge.4, reports, that all holes has been filled.
!
!     Input:  imt, jmt, km, jloc, kmsk, dat, spval, ldoc
!     Output: dat
!
!     Adapted from fillgaps by
!
!     J. Dengg, jdengg@awi_bremerhaven.de
!     Oct-97
!
!     Christian Dieterich Wed Oct 15 15:42:56 MET DST 1997
!     Modified Thu Oct 16 10:34:47 MET DST 1997
!
!------------------------------------------------------------------------

      integer :: stdout = 6
      integer,intent(in) :: imt, jmt, km           ! field dimensions
      integer,intent(in) :: jloc                   ! index meridional position

      integer,intent(in) :: kmsk(imt,jmt)          ! kmt- or kmu- field
      real,intent(inout) :: dat(imt,km)            ! data
      real,intent(in)    :: spval                  ! special value for 
!                                                  !         data holes
      integer,intent(in) :: ldoc                   ! level of documentation
      logical,intent(in) :: stop_on_error          ! stop if gaps remain

      integer :: nholes                            ! no. of initial holes
      integer :: icount                            ! no. of remaining holes



!     Security Check
!     --------------

      if (jloc > jmt) then
        write (stdout,*) 'Illegal j-index. ',
     &                   'Error in subroutine fillholes.' 
        stop '==> fillholes'
      end if



!     Count and fill holes with nearest neighbour from above
!     ------------------------------------------------------

      icount = 0                                 ! initialize counters
        igap = 0
        kgap = 0
       nholes = 0

      do k=1, km
        do i=1, imt
          if (kmsk(i,jloc)>=k .and. dat(i,k)==spval) then
            icount = icount + 1
            if (ldoc.ge.3) write (stdout,*) 'gap #', icount,
     &      'at i=', i, ', k=', k, ', j=', jloc
            if(k>1)dat(i,k)=dat(i,k-1)
          end if
        end do
      end do

      nholes = icount

      if (ldoc.ge.2) write (stdout,*) 'Data holes: ', nholes

      icount = 0                                 ! initialize counters

      do k=1, km
        do i=1, imt
          if (kmsk(i,jloc)>=k .and. dat(i,k)==spval) then
            icount = icount + 1
            if (ldoc.ge.1) write (stdout,*) 'gap #', icount,
     &      'at i=', i, ', k=', k, ', j=', jloc, ' remaining'
          end if
        end do
      end do
      if (icount==0) then
        if (ldoc.ge.4)
     &  write (stdout,'(a/)') '           all holes eliminated'
        goto 900
      end if

      if (ldoc>=1) 
     &write (stdout,*) '           ', icount, ' holes remaining'
      if (stop_on_error) stop

  900 continue

      end subroutine fillholes




      subroutine read_stamp(stamp, time )
      implicit none
      character*32, intent(in) :: stamp
      real, intent(out) :: time

      character *6  junk6
      character *8  junk8
      character *1  junk1
      integer       year, month, day, hour, min, sec

      integer nmonth
      parameter (nmonth = 12)
      integer daypm(nmonth), msum(nmonth), yrlen
      data daypm /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
      logical leap
      integer i,iyr,nfh,nday

      read  (stamp, '(a6,i2,a1,i2,a1,i4,a8,i2,a1,i2,a1,i2)')
     &       junk6, month, junk1, day, junk1, year, junk8, hour,
     &       junk1, min, junk1, sec
      msum(1) = 0
      do i = 2,nmonth
        msum(i) = msum(i-1) + daypm(i-1)
      enddo
      yrlen = msum(nmonth) + daypm(nmonth)
      if (mod (year, 400) .eq. 0) then
        leap = .true.
      elseif ((mod(year,4) .eq. 0) .and. (mod(year,100) .ne. 0)) then
        leap = .true.
      else
        leap = .false.
      endif
      nday = day + msum(month)
      if (leap .and. (month .gt. 2)) then
        nday = nday + 1
      endif
      iyr = year - 1601
c-----------------------------------------------------------------------
c     make corrections for years before 1601.
c-----------------------------------------------------------------------
      if (iyr .lt. 0) then
        nfh  = 1 + (-iyr)/400
        nday = nday - 146097*nfh
        iyr  = iyr + 400*nfh
      endif
      nday = nday + 365*iyr
c-----------------------------------------------------------------------
c     correct for leap-years between 1601 and iyear-1.
c     A virtual calendar consistent with current leap-year systems
c     is used to extend backwards.  (ie. 1700, 1800, 1900 are not
c     considered leap-year but 2000 is a leap-year)
c-----------------------------------------------------------------------
      iyr = iyr/4
      nday = nday + iyr
      iyr = iyr/25
      nday = nday - iyr
      iyr = iyr/4
      nday = nday + iyr
c-----------------------------------------------------------------------
c     nday is now in days since Dec. 31, 1600.  Convert to days since
c     Dec 31, 1899.
c-----------------------------------------------------------------------
      nday = nday - 109207
      time=float(nday)
     &  + (float(hour) + (float(min)+ float(sec)/60.)/60.  )/24.
      end subroutine read_stamp


