!
! Utility procedures for tgcmproc:
! This is NOT a module, just a collection of externals.
!
!-------------------------------------------------------------------
      subroutine fminmax(f,n,rmin,rmax)
!
! Return min and max of array f(n), excluding any values == spval:
!
      use proc, only: spval
      implicit none
      integer, intent(in)  :: n
      real,    intent(in)  :: f(n)
      real,    intent(out) :: rmin,rmax
      integer i
!
!     rmin = spval
!     rmax = -spval
!
! On dataproc with 4-byte reals, huge(x) = 0.34028235E+39
      rmin = huge(rmin)
      rmax = -rmin
      do i=1,n
        if (f(i).ne.spval.and.f(i).gt.rmax) rmax = f(i)         
        if (f(i).ne.spval.and.f(i).lt.rmin) rmin = f(i)         
      enddo
      return
      end subroutine fminmax
!-------------------------------------------------------------------
      subroutine min_to_mtime(mins,iday,ihr,min)
!
! Given total time in minutes (min), return iday,ihr,min:
!
      implicit none
      integer,intent(in)  :: mins
      integer,intent(out) :: iday,ihr,min
      integer :: minutes
!
      minutes = mins
      iday = minutes/1440
      minutes = minutes - iday*1440
      ihr = minutes/60
      min = minutes - ihr*60
      return 
      end subroutine min_to_mtime
!-------------------------------------------------------------------
      integer function mtime_to_mins(mtime)
      implicit none
!
! Given model time mtime (day,hr,min), return equivalent time
! in minutes (includes day):
!
! Arg:
      integer,intent(in)  :: mtime(3)
!
      mtime_to_mins = mtime(1)*24*60+mtime(2)*60+mtime(3)
      end function mtime_to_mins
!-------------------------------------------------------------------
      logical function int_is_str(int,str)
      implicit none
!
! Return true if integer i is string str, false otherwise:
!
! Args:
      integer,intent(in) :: int
      character(len=*),intent(in) :: str
!
! Locals:
      character(len=80) str80
      integer :: i,lstr,lstr80
!
      int_is_str = .false.
      str80 = ' '
      write(str80,"(a)") int
      lstr = len_trim(str)
      lstr80 = len_trim(str80)
      if (lstr.ne.lstr80) return
      if (str(1:lstr).eq.str80(1:lstr80)) int_is_str = .true.
      return
      end function int_is_str
!-------------------------------------------------------------------
      logical function float_is_str8(f,str)
      implicit none
!
! Return true if float f is string str, false otherwise:
!
! Args:
      real(kind=8),intent(in) :: f
      character(len=*),intent(in) :: str
!
! Locals:
      character(len=80) str80
      integer :: lstr,lstr80
!
      float_is_str8 = .false.
      str80 = ' '
      write(str80,"(a)") f
      lstr = len_trim(str)
      lstr80 = len_trim(str80)
      if (lstr /= lstr80) return
      if (str(1:lstr)==str80(1:lstr80)) float_is_str8 = .true.
      return
      end function float_is_str8
!-------------------------------------------------------------------
      logical function float_is_str(f,str)
      implicit none
!
! Return true if float f is string str, false otherwise:
!
! Args:
      real,intent(in) :: f
      character(len=*),intent(in) :: str
!
! Locals:
      character(len=80) str80
      integer :: lstr,lstr80
!
      float_is_str = .false.
      str80 = ' '
      write(str80,"(a)") f
      lstr = len_trim(str)
      lstr80 = len_trim(str80)
      if (lstr /= lstr80) return
      if (str(1:lstr)==str80(1:lstr80)) float_is_str = .true.
      return
      end function float_is_str
!-------------------------------------------------------------------
      integer function nextlu()
      implicit none
      logical isopen
      integer lu
      integer,save :: lureq(99-7+1)=0	! lu's given out so far
!
! Return an unopened fortan logical unit number (not 5 or 6):
! 9/97: Do not return a previously returned lu.
!
      do lu=7,99
        inquire(lu,opened=isopen)
        if (.not.isopen) then
          if (lureq(lu)<=0) then	! do not use previously returned lu
            nextlu = lu
            lureq(lu) = lu 
            return
          endif
        endif
      enddo
      write(6,"(/'>>> nextlu: all logical units apparently in use')")
      nextlu = 0
      stop 'nextlu'
      end function nextlu
!-------------------------------------------------------------------
      integer function ixfindc(strarr,nstr,searchstr)
      implicit none
!
! Given string array strarr, return index to strarr which contains 
!   string searchstr, or return 0 if searchstr not found,
!   searchstr is 0 length, or strarr is 0 length.
!
! Args:
      character(len=*),intent(in) :: strarr(nstr),searchstr
      integer,intent(in) :: nstr ! extent of strarr to search
!
! Locals:
      integer :: lenele,	! length of a strarr element
     |           lsearch,	! length of the search string
     |           i		! loop index
!
      ixfindc = 0
      if (nstr == 0) then
        write(6,"('WARNING ixfindc: nstr=0 (length string array)')")
        return
      endif
      lsearch = len_trim(searchstr)
      do i=1,nstr
        lenele = len_trim(strarr(i))
        if (lenele > 0) then
          if (strarr(i)(1:lenele) == searchstr(1:lsearch)) then
            ixfindc = i
            return 
          endif
        endif
      enddo
      return
      end function ixfindc
!-------------------------------------------------------------------
      subroutine allocerr(ier,msg)
      implicit none
!
! Args:
      character(len=*) :: msg
      integer,intent(in) :: ier
!
      write(6,"(/'>>> An error has occurred from an allocate ',
     +  'statement')")
      write(6,"('ier=',i8)") ier
      if (len_trim(msg) > 0) write(6,"(a)") msg
      stop 'allocerr'
      end subroutine allocerr
!-------------------------------------------------------------------
      function vecsum(u,v,id1,id2,spv)
!
! Return vector sum sqrt(u**2+v**2), where u and v are 2d arrays
!   of the same shape.
! This is an array-valued function which returns the vector sum
!   in an array which is the same 2d shape as u and v.
!
      implicit none
!
      integer,intent(in) :: id1,id2
      real,intent(in) :: u(id1,id2),v(id1,id2) 
      real:: vecsum(id1,id2)			! result variable
      real,intent(in),optional :: spv
!
      integer :: i,ii
!
      if (.not.present(spv)) then
        vecsum = sqrt(u**2 + v**2)
      else
        do i=1,id2
          do ii=1,id1
            if (u(ii,i)/=spv.and.v(ii,i)/=spv) then
              vecsum(ii,i) = sqrt(u(ii,i)**2+v(ii,i)**2)
            else
              vecsum(ii,i) = spv
            endif
          enddo
        enddo
      endif
      return
      end function vecsum
!-------------------------------------------------------------------
      subroutine hnmf2(ht,ne,hmf2,nmf2,fof2,name,imx,kmx)
!
! Find the NmF2, foF2, and hmF2 based on the electron profile.
! btf: This is taken from an idl procedure provided by Liying.
! fof2 added April, 2012.
!
      implicit none
!
! Args:
      integer,intent(in) :: imx,kmx
      real,intent(in),dimension(imx,kmx) :: ht,ne
      real,intent(out),dimension(imx) :: hmf2,nmf2,fof2
      character(len=*),intent(in) :: name
!
! Local:
      integer :: k,kx,i
      real :: h(3),p(3),h12,h22,h32,deltx,atx,ax,btx,bx,ctx,cx
!
      do i=1,imx
        kx = 0
        kloop: do k = kmx-1,2,-1
          if (ne(i,k) >= ne(i,k-1).and.ne(i,k) >= ne(i,k+1)) then
            kx = k
            exit kloop
          endif
        enddo kloop
        if (kx==0) then
          write(6,"('>>> hnmf2: could not find kx -- ',
     |      'setting hmf2=nmf2=fof2=0')") 
          hmf2(:) = 0. ; nmf2(:) = 0. ; fof2(:) = 0.
          return
        endif
        h = (/ht(i,kx-1),ht(i,kx),ht(i,kx+1)/)
        p = (/ne(i,kx-1),ne(i,kx),ne(i,kx+1)/)
        h12 = h(1)*h(1)
        h22 = h(2)*h(2)
        h32 = h(3)*h(3)
        deltx=h12*h(2)+h22*h(3)+h32*h(1)-h32*h(2)-h12*h(3)-h22*h(1)
        atx=p(1)*h(2)+p(2)*h(3)+p(3)*h(1)-h(2)*p(3)-h(3)*p(1)-h(1)*p(2)
!
! Periodic points are not set in ht and ne, so they are zero
! at i=1,2. Setting hmf2and nmf2 zero here is ok, since they
! are saved on histories at i=3,nlon. Without this check, code 
! will stop here w/ divide by zero deltx when debug flags are set.
!
!       if (deltx == 0.) then
!         hmf2(:) = 0. ; nmf2(:) = 0. ; fof2(:) = 0.
!         cycle
!       endif
        ax=atx/deltx
        btx=h12*p(2)+h22*p(3)+h32*p(1)-h32*p(2)-h12*p(3)-h22*p(1)
        bx=btx/deltx
        ctx=h12*h(2)*p(3)+h22*h(3)*p(1)+h32*h(1)*p(2)-h32*h(2)*p(1)-
     |    h12*h(3)*p(2)-h22*h(1)*p(3)
        cx=ctx/deltx
        hmf2(i)=-(bx/(2.*ax))
        nmf2(i)=-((bx*bx-4.*ax*cx)/(4.*ax))
        fof2(i)=sqrt(nmf2(i)*1.e6/1.24e10)
      enddo ! i=lon0,lon1
!     if (name=='FOF2') then
!       write(6,"('hnmf2: fof2 min,max=',2e12.4)") 
!    |    minval(fof2),maxval(fof2)
!     elseif (name=='HMF2') then
!       write(6,"('hnmf2: hmf2 min,max=',2e12.4)") 
!    |    minval(hmf2),maxval(hmf2)
!     elseif (name=='NMF2') then
!       write(6,"('hnmf2: nmf2 min,max=',2e12.4)") 
!    |    minval(nmf2),maxval(nmf2)
!     endif
      end subroutine hnmf2
!-------------------------------------------------------------------
      real function fslt(slt,ut,glon,iflag)
      implicit none
!
! If iflag = 1 -> return slt, given ut and glon (slt=ut+glon/15.)
! If iflag = 2 -> return ut, given slt and glon (ut=slt-glon/15.)
! If iflag = 3 -> return glon, given slt and ut (glon=(slt-ut)*15.)
! (slt and ut always in decimal hours, glon in -180 -> 180)
!
! Args:
      real :: slt,ut,glon
      integer,intent(in) :: iflag
!
      if (iflag==1) then
        fslt = ut + glon/15.
        if (fslt < 0.) fslt = fslt+24.
        if (fslt >= 24.) fslt = fslt-24.
      elseif (iflag==2) then
        fslt = slt - glon/15.  
        if (fslt < 0.) fslt = fslt+24.
        if (fslt > 24.) fslt = fslt-24.
      elseif (iflag==3) then
        fslt = (slt-ut) * 15.
        if (fslt >= 180.) fslt = fslt-360.
        if (fslt < -180.) fslt = fslt+360.
      else
        write(6,"('>>> fslt: bad iflag=',i3)") iflag
        stop 'fslt'
      endif
      return
      end
!-------------------------------------------------------------------
      real function quadrat(f,hts,npts)
!
! Integrate f over hts (hts are input in km):
!
      use proc, only: spval
      implicit none
!
! Args:
      integer,intent(in) :: npts
      real,intent(in) :: f(npts),hts(npts)
!
! Locals:
      real :: htcm(npts)
      integer :: k
!
      htcm = hts*1.e5	! whole array op
      quadrat = 0.
      do k=1,npts-1
        if (f(k) /= spval.and.f(k+1) /= spval) then

!         write(6,"('quadrat: k=',i3,' f(k)=',e12.4,' f(k+1)=',e12.4,
!    +      ' htcm(k)=',f10.1,' htcm(k+1)=',f10.1)")
!    +      k,f(k),f(k+1),htcm(k),htcm(k+1)

          quadrat = quadrat+0.5*(f(k)+f(k+1))*(htcm(k+1)-htcm(k)) 
        endif
      enddo
      return                                                            
      end                                                               
!-------------------------------------------------------------------
      subroutine glb_integ(f3d,fhts,fglb,imx,kmx,jmx)
      implicit none
!
! Do global height integration of f3d, return result in fglb:
! (input fhts are in km)
!
! Args:
      integer,intent(in) :: imx,kmx,jmx
      real,intent(in) :: f3d(imx,jmx,kmx),fhts(imx,jmx,kmx)
      real,intent(out) :: fglb(imx,jmx)
!
! Locals:
      integer :: i,j
      real,external :: quadrat
!
      do j=1,jmx
        do i=1,imx
          fglb(i,j) = quadrat(f3d(i,j,:),fhts(i,j,:),kmx)
        enddo
      enddo
      return
      end
!-------------------------------------------------------------------
      integer function ixfind8(z,nz,val,del)
      implicit none
!
! Find nearest index to val in z(nz), with increment del
! If not found, return -1
!
! Args:
      integer,intent(in) :: nz
      real,intent(in) :: z(nz),del
      real(kind=8),intent(in) :: val
!
! Locals:
      integer :: n
!
      do n=1,nz
        if (val.ge.z(n)-0.5*del.and.val.le.z(n)+0.5*del) then
          ixfind8 = n
          return
        endif
      enddo
      ixfind8 = -1
      return
      end
!-------------------------------------------------------------------
      integer function ixfind(z,nz,val,del)
      implicit none
!
! Find nearest index to val in z(nz), with increment del
! If not found, return -1
!
! Args:
      integer,intent(in) :: nz
      real,intent(in) :: z(nz),val,del
!
! Locals:
      integer :: n
!
      do n=1,nz
        if (val.ge.z(n)-0.5*del.and.val.le.z(n)+0.5*del) then
          ixfind = n
          return
        endif
      enddo
      ixfind = -1
      return
      end
!-------------------------------------------------------------------
      real function rnd(val,rndval)
      real,intent(in) :: val,rndval
      integer :: ival
      ival = (val + (rndval-.0001)) / maxval((/rndval,1./))
      rnd = ival * maxval((/rndval,1./))
      return
      end 
!-------------------------------------------------------------------
      subroutine vecterp(fy,f,ny,fylin,flin,nlin,spval,iprnt)
      implicit none
!
! Args:
      integer,intent(in) :: ny,nlin,iprnt
      real,intent(in) :: fy(ny),f(ny),fylin(nlin)
      real,intent(out) :: flin(nlin)
      real,intent(in) :: spval
!
! Locals
      integer :: inc,k,kbot,nk,k1,k0,ier,i
!      
      if (ny.lt.2) then
        write(6,"('>>> vecterp: ny must be >= 2: ny=',i3)") ny
        return
      endif
!
! Determine if fy array increases or decreases:
!
      inc = 1
      do k=2,ny
        if (fy(k).ne.spval.and.fy(k-1).ne.spval) then
          if (fy(k)-fy(k-1).lt.0.) inc = -1
          goto 100
        endif
      enddo
      if (iprnt.gt.0)
     +  write(6,"('>>> warning vecterp: could not determine ',
     +    'if fy is increasing or decreasing...')")
 100  continue
!
! Use only that part of fy .ne. spval.
! kbot = first non-spval in fy:
!
      kbot = 0
      do k=1,ny
        if (fy(k).ne.spval) then
          kbot = k
          goto 101
        endif 
      enddo
      if (iprnt.gt.0) then
        write(6,"('>>> warning vecterp: fy all spval?')")
        do k=1,nlin
          flin(k) = spval
        enddo
        return
      endif
 101  continue
!
! Assume that all spval's are either at beginning and/or at end
! (no isolated spval's in middle of array)
! nk = number of fy values to use:
!
      nk = 0
      do k=kbot,ny
        if (fy(k).ne.spval) nk = nk+1
      enddo
!
!     subroutine bracket(x,xx,nx,inc,n1,n2,ier)
! Bracket x in xx(nx), returning lower index in n1 and upper index in n2
!
      do k=1,nlin
        call bracket(fylin(k),fy(kbot),nk,inc,k0,k1,ier)
        if (ier.gt.0) then
          if (iprnt.gt.0) 
     +    write(6,"('>>> vecterp: err from bracket: k=',i2,' nlin=',i3,
     +      ' fylin(k)=',
     +      e12.4,' ny=',i2,' inc=',i2,' fy=',/(6e12.4))") 
     +      k,nlin,fylin(k),ny,inc,(fy(i),i=1,ny)
          flin(k) = spval
        else
          if (fy(k1)-fy(k0).gt.1.e-20.or.fy(k1)-fy(k0).lt.-1.e-20) then
            flin(k) = f(k0) + (f(k1)-f(k0))*(fylin(k)-fy(k0))/
     +        (fy(k1)-fy(k0))
          else
            flin(k) = spval
          endif
        endif
      enddo
      return
      end
!-------------------------------------------------------------------
      subroutine tail(path,file)
      implicit none
!
! Given path, return tail part (i.e., file name)
!
! Args:
      character(len=*),intent(in) :: path
      character(len=*),intent(out) :: file
!
! Locals:
      integer :: lenpath,islash,ix,i
!
      lenpath = len_trim(path)
      if (lenpath <= 0) then
        write(6,"('tail: lenpath=',i3,' (path is all blanks)',
     +    ' -- returning a single blank')") lenpath
        file = ' '
        return
      endif
      islash = 0
      ix = 1
      do i=1,lenpath
        if (path(i:i).eq.'/') then
          islash = 1
        else
          if (islash.gt.0) then
            ix = i
            islash = 0
          endif
        endif
      enddo
      file = path(ix:lenpath)
      return
      end
!-------------------------------------------------------------------
      real function convlon(glon,max)
      implicit none
!
! Given a longitude glon:
! If max = 180, return longitude equal to glon but in -180 to +180 system
! If max = 360, return longitude equal to glon but in 0 to 360 system
!
      real,intent(in) :: glon
      integer,intent(in) :: max
      convlon = glon
      if (max==180) then
        if (convlon >  180.) convlon = convlon-360.
        if (convlon < -180.) convlon = convlon+360.
      elseif (max==360) then
        if (convlon.lt.0.) convlon = convlon + 360.
      endif
      return
      end
!-------------------------------------------------------------------
      integer function ipat(ibts)
      implicit none
!
! Return decimal integer converted from 16-bit pattern
! provided as int ibts(16) (for use in dash patterns, etc):
!
      integer,intent(in) :: ibts(16)
      integer :: i
      ipat = 0
      do i=1,16
        ipat = ior(ishft(ipat,1),ibts(i))
      enddo
      return
      end
!-------------------------------------------------------------------
      subroutine log10f(x,n,rmin,spv)
      implicit none
!
! Given a field x, n vals long, change field to log10 of itself
! unless a value is <= rmin. If a value is <= rmin, set that
! value to spval. (if a value is already spval, leave it)
!
      integer,intent(in) :: n
      real,intent(inout) :: x(n)
      real,intent(in) :: spv,rmin
      integer :: i
      do i=1,n
        if (x(i) /= spv) then
          if (x(i) <= rmin) then
            x(i) = spv
          else
            x(i) = log10(x(i))
          endif
        endif
      enddo
      return
      end
!-------------------------------------------------------------------
      character(len=*) function mkflnm(path,checkopen)
      implicit none
!
! Given full path, return string same as path except with '/' replaced
! by '.', except when '/' occurs in 1st position, in which case it is
! removed from output string. 
! (this is used to make a disk file name from an absolute path)
!
! Args:
      character(len=*),intent(in) :: path
      logical,intent(in) :: checkopen
!
! Locals:
      integer :: lenpath,indx=1
      logical :: exists,isopen
      character(len=2) :: strnum
!
      mkflnm = path
      lenpath = len_trim(path)
      if (lenpath.eq.0) then
        write(6,"('>>> WARNING mkflnm: zero length path')") 
        return
      endif
      do
        indx = index(mkflnm,'/')
        if (indx.eq.0) exit
        if (indx.eq.1) then
          mkflnm(indx:indx) = ' '
        else
          mkflnm(indx:indx) = '.'
        endif
      enddo
      mkflnm = trim(adjustl(mkflnm))
!
! If checkopen is set and file is already open, make new names 
!   until file either does not exist or exists but is not open: 
!
      if (checkopen) then
        lenpath = len_trim(mkflnm)
        inquire(file=mkflnm(1:lenpath),exist=exists,opened=isopen)
        if (.not.exists.or..not.isopen) return
        if (lenpath+2 > len(mkflnm)) then
          write(6,"('>>> WARNING mkflnm: name too long -- ',
     +      'need 2 more chars. len(mkflnm)=',i3)") len(mkflnm)
          return
        endif
        do indx=0,99
          if (.not.exists.or..not.isopen) then
            write(6,"(8x,'will use file name ',a)") trim(mkflnm)
            return
          endif
          write(6,"(/'mkflnm: file ',a,' already open. Making new',
     +      ' file name...')") trim(mkflnm)
          write(strnum,"(i2.2)") indx 
          mkflnm = mkflnm(1:lenpath)//trim(adjustl(strnum))
          inquire(file=trim(mkflnm),exist=exists,opened=isopen)
        enddo
        write(6,"('>>> WARNING mkflnm: could not find an unopened',
     +    ' file name: mkflnm=',a)") trim(mkflnm)
      endif	! check open
      end function mkflnm
!-------------------------------------------------------------------
      character(len=16) function float_to_str(f)
      implicit none
      real,intent(in) :: f
      write(float_to_str,"(a)") f
      return
      end
!-------------------------------------------------------------------
      character(len=16) function float8_to_str(f)
      implicit none
      real(kind=8),intent(in) :: f
      write(float8_to_str,"(a)") f
      return
      end
!-------------------------------------------------------------------
      real function fmean(f,n,spval,iprnt)
      implicit none
!
! Return mean of f(n), where f(i) /= spval:
!
! Args:
      integer,intent(in) :: n,iprnt
      real,intent(in) :: f(n),spval
!
! Locals:
      integer :: i,nok
!
      fmean = 0.
      nok = 0	! number of values /= spval
      do i=1,n
        if (f(i) /= spval) then
          fmean = fmean+f(i)
          nok = nok+1
        endif
      enddo
      if (nok==n) then			! all good values (no spval)
        fmean = fmean / float(n)
      elseif (nok<n.and.nok>0) then	! some spval
        fmean = fmean / float(nok)
      else				! all spval
        fmean = spval
      endif
      if (nok /= n.and.iprnt > 0)
     +  write(6,"('Note: fmean encountered spval: n=',i3,
     +    ' nok=',i3,' fmean=',e12.4)") n,nok,fmean
      return
      end function fmean
!-------------------------------------------------------------------
      logical function isslt(f,slt)
      use proc,only: sltflag
      implicit none
!
! Return true if first 2 or 3 chars of float f are "SLT" or "LT",
!   (or "slt" or "lt"), otherwise return false.
! If true, then return slt, which is in the rest of f, i.e., the
!   requested local time. If false, set slt = -1.
!
! Args:
      real(kind=8),intent(in) :: f
      real,intent(out) :: slt
!
! Locals:
      character(len=8) sltstr
!
      isslt = .false.
      slt = -1.
!
! Check for sltflag (AIX) (see sub checkslt in fixnamelist.f):
      if (f >= sltflag .and. f <= sltflag+24.) then
        isslt = .true.
        slt = f-sltflag
      endif
!
! Check for slt encoded in real (SGI):
      sltstr=' '
      write(sltstr,"(a)") f
      if (sltstr(1:2).eq.'LT'.or.sltstr(1:3).eq.'SLT'.or.
     +    sltstr(1:2).eq.'lt'.or.sltstr(1:3).eq.'slt') then
        isslt = .true.
        if (sltstr(1:2).eq.'LT'.or.sltstr(1:2).eq.'lt')
     +    read(sltstr(3:8),"(f6.2)") slt
        if (sltstr(1:3).eq.'SLT'.or.sltstr(1:3).eq.'slt')
     +    read(sltstr(4:8),"(f5.2)") slt
      endif
      return
      end
!-------------------------------------------------------------------
      real function fglbm(f,imx,jmx,gcmlat,dlat,dlon,spv)
      implicit none
!
! Return global mean of f(imx,jmx), weighted by cos(lat):
!
! Args:
      integer,intent(in) :: imx,jmx
      real,intent(in) :: f(imx,jmx),gcmlat(jmx),dlat,dlon,spv
!
! Locals:
! 9/11/03 bf: Changed re=63711.e+5 to 6371.e+5
!
      real,parameter :: pi=3.14156, dtr=pi/180., re=6371.e+5
      integer,save :: ncalls=0
      real,save :: sud,sum
      real :: abb
      integer :: i,j
!
      ncalls = ncalls+1
      if (ncalls==1) then
        sud = re*re*dlon*dlat*dtr*dtr
        sum = 0.
        do i=1,imx
          do j=1,jmx
            sum = sum + sud*cos(gcmlat(j)*dtr)
          enddo
        enddo
      endif
!
! Should find more efficient way to check for spval:
      fglbm = 0.
      jloop: do j=1,jmx 
        abb = cos(gcmlat(j)*dtr) * sud / sum
        do i=1,imx
          if (f(i,j)==spv) then
            fglbm = spv
            exit jloop
          else
            fglbm = fglbm + f(i,j) * abb
          endif
        enddo
      enddo jloop
      return
      end
!-------------------------------------------------------------------
      subroutine wrdat(iwr,lu,f,nx,ny,xx,yy,xlab,ylab,fieldlab,
     +  infolab,histlab,iframe,proclab,senddat)
      use input,only: flnm_dat
      implicit none
!
! Write ascii data file containing field f(nx,ny)
!
! On input:
!   lu = logical unit connected to file (opened by input.F)
!   f(nx,ny) = data field to write (nx or ny can be 1 for 1d fields)
!   xx(nx) = x-axis values
!   yy(ny) = y-axis values
!   xlab,ylab = x,y axis labels
!   fieldlab = field name
!   infolab = information re grid slice, ut, etc.
!   histlab = history label
!   iframe = frame number
!   proclab = used to create file name (will be "proclab".dat)
!   senddat = machine:path to which individual frames are transferred (iwr=1 only)
! On output:
!   4-line header, xx, yy, and f, are written to lu (6e13.5 format for data)
!   Single-frame file(s) are sent to senddat if iwr=1
!   Output file(s) may be read and plotted from IDL code in 
!     ~foster/tgcmvis/pltdat
!
! 11/29/07 btf:
! Fixed this so flnm_dat is always overwritten at beginning of program
!   if it pre-exists.
! The ascii data file flnm_dat is now opened by sub input (if the user
!   set flnm_dat in namelist read), and closed at end of main routine 
!   tgcmproc.F. 
! Formerly, the iwr flag was 0, 1, or 2 for no write, write, or append,
!   but now is either 0 or 1, never 2 (see input.F).
!
! Args:
      integer,intent(in) :: iwr,lu,nx,ny,iframe
      real,intent(in) :: f(nx,ny),xx(nx),yy(ny)
      character(len=*),intent(in) :: infolab,proclab,histlab,fieldlab,
     +  senddat,xlab,ylab
!
! Locals:
      character(len=80) :: rec80,fname
      logical :: isopen,exists
      integer :: i
!
      if (iwr.le.0) return
c
c Field name (1st of 4-line header):
c
      rec80 = ' '
      if (len_trim(fieldlab).le.80) then
        write(rec80,"(a)") trim(fieldlab)
      else
        write(rec80,"(a)") fieldlab(1:80)
      endif
      write(lu,"(a)") rec80
c
c Info lab (2nd of 4-line header):
c
      rec80 = ' '
      if (len_trim(infolab).le.80) then
        write(rec80,"(a)") trim(infolab)
      else
        write(rec80,"(a)") infolab(1:80)
      endif
      write(lu,"(a)") rec80
c
c History volume (3rd of 4-line header):
c
      rec80 = ' '
      if (len_trim(histlab).le.80) then
        write(rec80,"(a)") trim(histlab)
      else
        write(rec80,"(a)") histlab(1:80)
      endif
      write(lu,"(a)") rec80
c
c Frame number (4th of 4-line header):
c
      rec80 = ' '
      write(rec80,"('Frame ',i4,' ')") iframe
      write(lu,"(a)") rec80
c
c X-axis, y-axis, and the field values:
c
      write(lu,"('nx=',i4,' xx  ',a)") nx,xlab
      if (nx.gt.1) write(lu,"(6e13.5)") xx
      write(lu,"('ny=',i4,' yy  ',a)") ny,ylab
      if (ny.gt.1) write(lu,"(6e13.5)") yy
      if (nx.gt.1.and.ny.gt.1) then
        write(lu,"(6e13.5)") f
      elseif (nx.eq.1.and.ny.gt.1) then
        write(lu,"(6e13.5)") (f(1,i),i=1,ny)
      elseif (ny.eq.1.and.nx.gt.1) then
        write(lu,"(6e13.5)") (f(i,1),i=1,nx)
      endif
      return
      end
!-------------------------------------------------------------------
      subroutine wrxdr(flnm,f,nx,ny,xx,yy,xlab,ylab,lab1,
     +  lab2,lab3,lab4,mtime,iclose)
      implicit none
!
! Write xdr data file containing field f(nx,ny)
!
! On input:
!   flnm = xdr file name
!   f(nx,ny) = data field to write (nx or ny can be 1 for 1d fields)
!   xx(nx) = x-axis values
!   yy(ny) = y-axis values
!   xlab,ylab = x,y axis labels
!   lab1,lab2,lab3,lab4 = information labels
!   mtime(3) = model time (day,hr,min) (set 0,0,0 for multi-ut plots)
! On output:
!   4-line header, xx, yy, and f, are written to lu (6e13.5 format for data)
! btf 4/8/97: added mtime (see also wrxdrc.c)
!
! Args:
      integer,intent(in) :: nx,ny,iclose,mtime(3)
      character(len=*),intent(in) :: flnm,lab1,lab2,lab3,lab4,xlab,ylab
      real,intent(in) :: f(nx,ny),xx(nx),yy(ny)
!
! Locals:
      character(len=80) :: hdr1,hdr2,hdr3,hdr4,xlabel,ylabel
      integer :: lenlab,istat
!
! Externals:
#if defined(__sgi)
      integer,external :: wrxdrc
#elif defined(AIX)
      integer,external :: wrxdrc_
#endif
!
      if (iclose > 0) goto 100
!
! lab1:
!
      hdr1 = ' '
      lenlab = len_trim(lab1)
      if (lenlab.le.80) then
        write(hdr1,"(a)") lab1(1:lenlab)
      else
        write(hdr1,"(a)") lab1(1:80)
      endif
!
! lab2:
!
      hdr2 = ' '
      lenlab = len_trim(lab2)
      if (lenlab.le.80) then
        write(hdr2,"(a)") lab2(1:lenlab)
      else
        write(hdr2,"(a)") lab2(1:80)
      endif
!
! lab3:
!
      hdr3 = ' '
      lenlab = len_trim(lab3)
      if (lenlab.le.80) then
        write(hdr3,"(a)") lab3(1:lenlab)
      else
        write(hdr3,"(a)") lab3(1:80)
      endif
!
! lab4:
!
      hdr4 = ' '
      lenlab = len_trim(lab4)
      if (lenlab.le.80) then
        write(hdr4,"(a)") lab4(1:lenlab)
      else
        write(hdr4,"(a)") lab4(1:80)
      endif
!
! x and y labels:
!
      xlabel = ' '
      lenlab = len_trim(xlab)
      if (lenlab.le.80) then
        write(xlabel,"(a)") xlab(1:lenlab)
      else
        write(xlabel,"(a)") xlab(1:80)
      endif
!
      ylabel = ' '
      lenlab = len_trim(ylab)
      if (lenlab.le.80) then
        write(ylabel,"(a)") ylab(1:lenlab)
      else
        write(ylabel,"(a)") ylab(1:80)
      endif
!
! int WRXDR(_fcd flnm, _fcd hdr1, _fcd hdr2, _fcd hdr3, _fcd hdr4,
!	  float *f2d, short *nx, short *ny, float *xx, float *yy,
!	  _fcd xlab, _fcd ylab, short *iclose)
!
 100  continue
#if defined(__sgi)
      istat = wrxdrc(flnm(1:len_trim(flnm)),hdr1,hdr2,hdr3,hdr4,
     +  f,nx,ny,xx,yy,xlabel,ylabel,mtime,iclose)
#elif defined(AIX)
      istat = wrxdrc_(flnm(1:len_trim(flnm)),hdr1,hdr2,hdr3,hdr4,
     +  f,nx,ny,xx,yy,xlabel,ylabel,mtime,iclose)
#endif
      return
      end
!-------------------------------------------------------------------
      function rotslt(forig,slt,ut,gcmlon,nlon,dlon,spval)
      implicit none
!
      real :: rotslt(nlon)	! function result
!
! Return rotslt(nlon) from forig(nlon), rotated such that local time
!   slt is in the center. 
! On input, forig(nlon) is ordered by gcmlon(nlon), and is unchanged.
!
! Args:
      integer,intent(in) :: nlon
      real,intent(in) :: forig(nlon),slt,ut,gcmlon(nlon),dlon,spval
!
! Locals
      integer :: i,ii,isltlon
      real :: sltlon,dum
      real,external :: fslt
      integer,external :: ixfind
!
      rotslt(:) = spval
      if (slt>=24..or.slt<0.) then
        write(6,"('>>> rotslt: bad slt=',f8.3)") slt
        return
      endif
      sltlon = fslt(slt,ut,dum,3)		! center longitude
      isltlon = ixfind(gcmlon,nlon,sltlon,dlon)	! index to center lon
      if (isltlon<=0) then
        write(6,"('>>> rotslt: could not find sltlon=',f8.2,
     +    ' in gcmlon.')")
        return
      endif
      ii = isltlon-nlon/2+1
      if (ii > nlon) ii=ii-nlon
      if (ii < 1) ii=ii+nlon
      do i=1,nlon
        if (ii > nlon) ii=ii-nlon+1
        if (ii < 1) ii=ii+nlon
        rotslt(i) = forig(ii)
        ii = ii+1
      enddo
      return
      end function rotslt
!-------------------------------------------------------------------
      function mkpmb(zp,nzp,p0,zprange,idimmb,spval,v5dhts)
!
! Given array of zp pressures zp(nzp), standard pressure p0,
!   and desired range zprange(2) (bottom,top), return array
!   of equivalent millibar values (mkpmb(k)==spval where zp(k)
!   is outside zprange).
! To convert from zp to mb use: mb = p0*exp(-zp)*1.e-3
! Optionally return v5dhts(idimmb) as vis5d-style heights
!   corresponding to the mb pressures (h=-7.2*ln(mb/1012.5))
!
      implicit none
!
      real :: mkpmb(idimmb)	! function result
!
! Args:
      integer,intent(in) :: nzp,idimmb
      real,intent(in) :: zp(nzp),p0,zprange(2),spval
      real,optional,intent(out) :: v5dhts(idimmb)
!
! Locals:
      integer :: k,kk
!
      mkpmb = 0.	! init
      kk = 0
      do k=1,nzp
        if (zp(k) >= zprange(1) .and. zp(k) <= zprange(2)) then
          kk = kk+1
          if (kk>idimmb) write(6,"('>>> mkpmb WARNING: kk=',i3,
     +      ' idimmb=',i3)") kk,idimmb
          mkpmb(kk) = p0*exp(-zp(k))*1.e-3 
!
! optionally return v5d type heights corresponding to each mb pressure:
          if (present(v5dhts)) v5dhts(kk) = -7.2*log(mkpmb(kk)/1012.5)
        endif
      enddo
      return
      end
!-------------------------------------------------------------------
      subroutine ddaytomt(decday,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal days (decday):
c
      md = int(decday)
      mh = (decday-float(md))*24.
      mm = (decday*(24.*60.)-float(md)*(24.*60.)-float(mh)*60.)+.5
      if (mm.eq.60) then
        mm = 0
        mh = mh+1
      endif
      if (md > 366) md = md-366
      return
      end
!-------------------------------------------------------------------
      subroutine dhrstomt(dechrs,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal hours (dechrs):
c
      rminday = 24.*60.
      md = dechrs/24.
      mh = dechrs-float(md)*24.
      if (mh.eq.24) then
        mh = 0
        md = md+1
      endif
      mm = dechrs*60.-float(md)*rminday-float(mh)*60.+.5
      return
      end
!-------------------------------------------------------------------
      subroutine dmintomt(decmin,md,mh,mm)
c
c Return model time (md,mh,mm), given time in decimal minutes (decmin):
c
      rminday = 24.*60.
      md = decmin/rminday+1.e-10
      mh = (decmin-float(md)*rminday)/60.
      mm = decmin-float(md)*rminday-float(mh)*60.
      if (mm.eq.60) mm = 0
      if (mh.eq.24) mh = 0
      return
      end
!-------------------------------------------------------------------
      function dhrs2ut(dechrs)
c
c Return ut, given time in dechrs (dechrs includes model day)
c
      call dhrstomt(dechrs,md,mh,mm)
      dhrs2ut = float(mh)+float(mm)/60.
      return
      end
!-------------------------------------------------------------------
      real function dut2dslt(dechrs,rlon)
c
c Convert decimal ut (includes model day) to decimal slt (including 
c model day)
c
      call dhrstomt(dechrs,md,mh,mm)
      ut = float(mh)+float(mm)/60.
      slt = fslt(0.,ut,rlon,1)
      dut2dslt = float(md)*24.+slt 
      return
      end
!-------------------------------------------------------------------
      real function cmt2dhrs(md,mh,mm)
        integer,intent(in) :: md,mh,mm
        cmt2dhrs = float(md)*24.+float(mh)+float(mm)/60.
      end function cmt2dhrs
!-------------------------------------------------------------------
      real function cmt2dday(md,mh,mm)
        integer,intent(in) :: md,mh,mm
        cmt2dday = float(md)+float(mh)/24.+float(mm)/(24.*60.)
      end function cmt2dday
!-------------------------------------------------------------------
      real function cmt2dmin(md,mh,mm) 
        integer,intent(in) :: md,mh,mm
        cmt2dmin = float(md)*24.*60.+float(mh)*60.+float(mm)
      end function cmt2dmin
!-------------------------------------------------------------------
      integer function addstrele(strarr,mx,newstr)
!
! Add newstr in 1st blank element of strarr(mx)
! (do not add if newstr alread in strarr)
! (return 1 if successful, otherwise 0)
!
! Args:
      integer,intent(in) :: mx
      character(len=*),intent(inout) :: strarr(mx)
      character(len=*),intent(in) :: newstr 
!
! Locals:
      integer :: i
!
      addstrele = 0
!
! Check that newstr will fit into strarr size:
      if (len(newstr)>len(strarr(1))) then
        write(6,"('>>> WARNING addstrele: len(newstr)=',i3,' > ',
     +    ' len(strarr(1)=',i3)") len(newstr),len(strarr(1))
        return
      endif
!
! Do not add redundant newstr:
      do i=1,mx
        if (trim(strarr(i))==trim(newstr)) return
      enddo
!
! Search and add if blank element of strarr is available:
      do i=1,mx
        if (len_trim(strarr(i))==0) then
          write(strarr(i),"(a)") trim(newstr)
          addstrele = 1
          return
        endif
      enddo
      return
      end function addstrele
!-------------------------------------------------------------------
      integer function ixohband(fname,hi,lo)
!
! Given a field name that ends in a 2-digit oh emission band
! (e.g., 62, 98, etc), return index to oh-band field array, as
! returned by getoh. Also return hi and lo as 1st and 2nd digits.
! E.g., if fname=='OH-98', then ixohband will be 10, hi=9, lo=8
!
! Array returned by getoh is dimensioned (nohalt,60) for 60
!   emission bands, as follows (only 39 are valid): 
!
! 0  1-0 2-1 3-2 4-3 5-4 6-5 7-6 8-7 9-8
! 0   0  2-0 3-1 4-2 5-3 6-4 7-5 8-6 9-7
! 0   0   0  3-0 4-1 5-2 6-3 7-4 8-5 9-6
! 0   0   0   0  4-0 5-1 6-2 7-3 8-4 9-5
! 0   0   0   0   0  5-0 6-1 7-2 8-3 9-4
! 0   0   0   0   0   0  6-0 7-1 8-2 9-3
!
      character(len=*),intent(in) :: fname
      integer,intent(out) :: hi,lo
      integer :: iband,len
!
      len = len_trim(fname)
      read(fname(len-1:len),"(i2)",err=100) iband
      hi = iband/10
      lo = iband-hi*10
      ixohband = (hi-lo-1)*10+hi+1
      return
 100  continue
      write(6,"('>>> WARNING ixohband: error reading emission ',
     +  'band from fname=',a)") fname
      ixohband = 0
      end function ixohband
!-------------------------------------------------------------------
      real function getsza(iday,slt,glat,glon)
      implicit none
      integer,intent(in) :: iday
      real,intent(in) :: slt,glat,glon
      real :: dtor,rtod,pi,c1,dle,dles,dlec,glatr,glonr,seczi
c
c Return solar zenith angle using input day of year iday (1-366),
c solar local time slt, and geographic glat,glon:
c
      dtor=atan(1.)/45.
      rtod=1./dtor
      pi=180./rtod
      c1=23.5*dtor
      dle=atan(tan(c1)*sin(2.*pi*(iday-80.)/365.))
      dles=sin(dle)
      dlec=cos(dle)
      glatr=glat*dtor
      glonr=glon*dtor
      seczi=1./(dles*sin(glatr)+cos(glatr)*cos(pi*(slt-12.)/12.)*dlec)  
      getsza=acos(1./seczi)*rtod
      return
      end
!-------------------------------------------------------------------
      subroutine mkho2(tn,fo2,fo1,fn2,fh,fo3,foh,fno,fho2_out,kmx)
      implicit none
!
! Return fho2_out(kmx), given tn and number densities of o2,o1,h,o3,oh,no
! (used by getoh when ho2 is not available on histories)
! (this is the same way the model calculates ho2)
!
      integer,intent(in) :: kmx
      real,dimension(kmx),intent(in) :: tn,fo2,fo1,fh,fo3,foh,fno,fn2
      real,intent(out) :: fho2_out(kmx)
!
      real :: r36,r29,r26,r34,b10	! reaction rates
      real :: rho
      integer :: k
!
      do k=1,kmx
        r36 = 5.7e-32 * (300./tn(k))**1.6
        r29 = 1.6e-12 * exp(-940./tn(k))
        r26 = 3.0e-11 * exp(200./tn(k))
        r34 = 1.1e-14 * exp(-500./tn(k))
        b10 = 3.5e-12 * exp(250./tn(k))
        fho2_out(k) = (r36*fo2(k)*fh(k)*(fo2(k)+fo1(k)+fn2(k))+
     +    r29*fo3(k)*foh(k)) / (r26*fo1(k)+r34*fo3(k)+b10*fno(k))
      enddo
      end subroutine mkho2
!-------------------------------------------------------------------
      real function bilin(f,xa,ya,nx,ny,x,y,spv,iprnt)
      implicit none
!
! Bilinear interpolation of f(nx,ny)
!
! Args:
      integer,intent(in) :: nx,ny,iprnt
      real,intent(in) :: f(nx,ny),xa(nx),ya(ny),x,y,spv 
!
! Locals:
      integer :: i,ix0,ix1,iy0,iy1
      real :: fracx,fracy
!
! Bracket requested x value:
      ix0 = 0 ; ix1 = 1
      xloop1: do i=1,nx
        if (xa(i)==x) then 
          ix0 = i ; ix1 = i ; cycle xloop1
        endif
      enddo xloop1
      if (ix0==0) then
        xloop2: do i=1,nx-1
          if ((x>xa(i).and.x<xa(i+1)).or.(x<xa(i).and.x>xa(i+1))) then
            ix0 = i ; ix1 = i+1 ; cycle xloop2
          endif
        enddo xloop2
      endif
      if (ix0==0) then
        if (iprnt > 0) 
     +    write(6,"('>>> WARNING bilin: could not find x=',e12.4)") x
        bilin = spv
        return 
      endif
      if (abs(xa(ix1)-xa(ix0)) <= 1.e-20) then
        fracx = 1.
      else
        fracx = (x-xa(ix0)) / abs(xa(ix1)-xa(ix0))
      endif
!
! Bracket requested y value:
      iy0 = 0 ; iy1 = 1
      yloop1: do i=1,ny
        if (ya(i)==y) then 
          iy0 = i ; iy1 = i ; cycle yloop1
        endif
      enddo yloop1
      if (iy0==0) then
        yloop2: do i=1,ny-1
          if ((y>ya(i).and.y<ya(i+1)).or.(y<ya(i).and.y>ya(i+1))) then
            iy0 = i ; iy1 = i+1 ; cycle yloop2
          endif
        enddo yloop2
      endif
      if (iy0==0) then
        if (iprnt > 0) 
     +    write(6,"('>>> WARNING bilin: could not find y=',e12.4)") y
        bilin = spv
        return 
      endif
      if (abs(ya(iy1)-ya(iy0)) <= 1.e-20) then
        fracy = 1.
      else
        fracy = (y-ya(iy0)) / abs(ya(iy1)-ya(iy0))
      endif
!
! Do bilinear interpolation:
      bilin = fracx*(fracy*f(ix1,iy1)+(1.-fracy)*fracx*f(ix1,iy0))+
     +  (1.-fracx)*(fracy*f(ix0,iy1)+(1.-fracy)*(1.-fracx)*f(ix0,iy0))
      end function bilin
!-------------------------------------------------------------------
      subroutine mkemislab(fname,ie5577,ie6300,labout)
      use fields,only: ne5577,e5577lab
      implicit none
!
! Args:
      character(len=8),intent(in) :: fname
      integer,intent(in) :: ie5577(ne5577),ie6300 
      character(len=*),intent(out) :: labout
!
! Locals:
      integer :: i,lenout
!
      lenout = len(labout)
      labout = ' '
      if (trim(fname)=='E5577') then
        e5577loop: do i=1,ne5577
          if (ie5577(i) > 0) then
            if (len_trim(labout)+len_trim(e5577lab(i))+1 > lenout) then
              write(6,"('>>> mkemislab: labout too short: lenout=',i3)")
     +          lenout
              exit e5577loop
            endif
            labout = trim(labout)//trim(e5577lab(i))//'+'
          endif
        enddo e5577loop
        i = len_trim(labout) 
        if (i==0) then
          write(6,"('>>> mkemislab: no e5577 components? ie5577=',
     +      5i3)") ie5577
          return
        endif
        labout(i:i) = ' '
      elseif (trim(fname)=='E6300') then
        if (ie6300 > 0) then
          write(labout,"('[SR63 ADDED]')")
        else
          write(labout,"('[SR63 NOT ADDED]')")
        endif
      endif
      end subroutine mkemislab
!-------------------------------------------------------------------
      subroutine bracket(x,xx,nx,inc,n1,n2,ier)
      implicit none
!
! Bracket x in xx(nx), returning lower index in n1 and upper index in n2
! If inc > 0 -> array increases from bottom to top,
! If inc <= 0 -> array increases from top to bottom
!
      integer,intent(in) :: nx,inc
      integer,intent(out) :: n1,n2,ier
      real,intent(in) :: x,xx(nx)
      integer :: i
!
! Array increases from bottom to top:
      ier = 0
      if (inc.gt.0) then
        if (x.lt.xx(1)) then
          n1 = 1
          n2 = 2
          ier = 1
          return
        endif
        if (x.gt.xx(nx)) then
          n1 = nx-1
          n2 = nx
          ier = 2
          return
        endif
        do i=1,nx-1
          if (x.ge.xx(i).and.x.le.xx(i+1)) then
            n1 = i
            n2 = i+1
            return
          endif
        enddo
!
! Array increases from top to bottom:
      else
        if (x.gt.xx(1)) then
          n1 = 1
          n2 = 2
          ier = 3
          return
        endif
        if (x.lt.xx(nx)) then
          n1 = nx-1
          n2 = nx
          ier = 4
          return
        endif
        do i=1,nx-1
          if (x.le.xx(i).and.x.ge.xx(i+1)) then
            n1 = i
            n2 = i+1
            return
          endif
        enddo
      endif
      return
      end subroutine bracket
!-------------------------------------------------------------------
      subroutine mkdiffs(pert,cntr,diffs,n,difftype)
      use proc, only: spval
!
! Take diffs of pert(n)-cntr(n) according to difftype (raw or percent),
!   return result in diffs(n). If (pert(i)==spval.or.cntr(i)==spval,
!   return diffs(i)=spval.
! A local array (diffs_out) is used to store the diffs before returning 
!   them in the output array provided (diffs) so that this routine can
!   be called when pert or cntr is the same array as diffs.
!   (e.g. call mkdiffs(a,b,a,...) (i.e., a = a-b), or 
!         call mkdiffs(b,a,a,...) (i.e., a = b-a))
!
! Args:
      integer,intent(in) :: n
      real,intent(in) :: pert(n),cntr(n)
      character(len=*),intent(in) :: difftype
      real,intent(out) :: diffs(n)
!
! Locals:
      real :: diffs_out(n)
!
! Init:
      diffs_out = spval
!
! Take raw diffs:
      if (trim(difftype)=='RAW') then
        where(pert/=spval.and.cntr/=spval) diffs_out = pert - cntr
!
! Take percent diffs:
      elseif (trim(difftype)=='PERCENT') then
        where(pert/=spval.and.cntr/=spval.and.cntr>1.e-20)
     +    diffs_out = (pert - cntr) / cntr * 100.
!
! Unknown diff type:
      else
        write(6,"('>>> WARNING mkdiffs: unknown difftype=',a)") 
     +    difftype
      endif
      diffs = diffs_out
      end subroutine mkdiffs
!-------------------------------------------------------------------
      subroutine calcstrm(vw,jmx,kmx,gcmzp,gcmlat,strm)
      use proc,only: dlev,pi,p0
      implicit none
!
! Calculate stream function from zonal means of vn or w
!
! On input:
!   vw(jmx,kmx) = zonal means of vn or w
!   jmx,kmx = latitude and pressure dimensions
!   gcmzp(kmx) = grid pressure levels
!   gcmlat(jmx) = grid geographic latitudes (deg)
! On output:
!   stream function strm(jmx,kmx) is defined
!
! Args:
      integer,intent(in) :: jmx,kmx
      real,intent(in) :: vw(jmx,kmx),gcmzp(kmx),gcmlat(jmx)
      real,intent(out) :: strm(jmx,kmx)
!
! Locals:
      real :: expz(kmx),cs(jmx)
      real :: dtr,alfa
      integer :: k,j,kk
      real,parameter :: grav=870., rad=6.37e+8
!
! Begin exec:
      dtr = pi/180.
      alfa = 2.*pi*rad*p0/grav
      do k=1,kmx
        expz(k) = exp(-gcmzp(k))
      enddo
      do j=1,jmx
        cs(j) = cos(dtr*gcmlat(j))
      enddo
      do j=1,jmx
        strm(j,kmx) = alfa*expz(kmx)*vw(j,kmx)*cs(j)
        do kk=1,kmx-1
          k = kmx-kk
          strm(j,k) = strm(j,k+1)+0.5*dlev*alfa*cs(j)*
     +      (expz(k+1)*vw(j,k+1)+expz(k)*vw(j,k))
        enddo
      enddo
      end subroutine calcstrm
!-------------------------------------------------------------------
      integer function fseries(f,mxf,fmin,fmax)
      use proc, only: spval
      implicit none
!
! Validate each element of f(mxf) according to fmin,fmax.
! If "short form" has been used (e.g., x,'to',y,'by',z), then
! expand to the full series and validate. Return number of
! valid elements of final expanded series.
!
! Args:
      real,intent(in) :: fmin,fmax
      integer,intent(in) :: mxf
      real(kind=8),intent(inout) :: f(mxf)
!
! Locals:
      logical isshort
      real :: del,f1,f2,ftmp
      integer :: i
!
! Externals:
      logical,external :: float_is_str8
!
! Determine if "short" form is used:
! short form: f(1),'to',f(3),'by',f(5)
!
      isshort = .false.
      if ((float_is_str8(f(2),'to').or.
     |     float_is_str8(f(2),'TO')).and.
     |    (float_is_str8(f(4),'by').or.
     |     float_is_str8(f(4),'BY'))) isshort = .true.
!
      fseries = 0
!
! Short form is in use -- validate and expand series:
!
      if (isshort) then
        if (f(1) < fmin .or. f(1) > fmax) then
          write(6,"('>>> fseries: bad 1st value: f(1)=',e12.4,
     |      ' fmin,fmax=',2e12.4)") f(1),fmin,fmax
          return
        endif
        if (f(3) < fmin .or. f(3) > fmax) then
          write(6,"('>>> fseries: bad last value: f(3)=',e12.4,
     |      ' fmin,fmax=',2e12.4)") f(3),fmin,fmax
          return
        endif
        f1 = f(1)   ! 1st value
        f2 = f(3)   ! last value
        del = f(5)  ! delta
        fseries = 0
        f = spval   ! reset f
        loop: do i=1,mxf
          ftmp = f1+(i-1)*del
          if (ftmp < fmin .or. ftmp > fmax .or. ftmp > f2) exit loop
          f(i) = ftmp
          fseries = fseries+1 
        enddo loop
!
! Short form not in use -- validate given values:
!
      else
        do i=1,mxf
          if (f(i) /= spval) then
            if (f(i) < fmin .or. f(i) > fmax) then
              write(6,"('>>> WARNING fseries: value ',e12.4,
     |          ' out of range: fmin,max=',2e12.4)") f(i),fmin,fmax
              f(i) = spval
            else
              fseries = fseries+1
            endif            
          endif
        enddo
      endif ! short form
!     write(6,"('fseries returning: fseries=',i3,' f=',
!    |  /(6e12.4))") fseries,f
      end function fseries
!-------------------------------------------------------------------
      integer function nunique_r(r,n,spval)
!
! Given real r(n), return number of unique non-spval elements in r:
!
! Args:
      integer,intent(in) :: n
      real(kind=8),intent(in) :: r(n)
      real,intent(in) :: spval
!
! Local:
      integer :: i
      real(kind=8) :: rtmp(n)
!
      nunique_r = 0
      rtmp = spval
      do i=1,n
        if (r(i) /= spval) then
          if (.not.any(rtmp==r(i))) then
            nunique_r = nunique_r+1
            rtmp(nunique_r) = r(i)
          endif
        endif
      enddo
      end function nunique_r
!-------------------------------------------------------------------
      integer function rmcomments(lu,comcharin,echo)
      implicit none
!
! Read input lines from unit lu. If current line contains the comment
!   character comcharin, strip the line from position of comchar to end,
!   and write any remaining line to a new unit. If no comment in current 
!   line, write entire line to new unit. 
! Return new unit, rewound (e.g., ready to be read by namelist).
! If echo > 0, echo output lines to stdout.
! If comcharin is ' ', then default comment char is ';'
!
! Args:
      integer,intent(in) :: lu,echo
      character(len=1),intent(in) :: comcharin
! Local:
      character(len=1) :: comchar
      logical isopen
      integer :: i,lens,ios,compos,nline,nlcfields
      character(len=10000) :: line
      character(len=64) :: newcfields(30)
! Externals:
      integer,external :: nextlu
!
      if (lu <= 0) then
        write(6,"('>>> rmcomments: bad input lu=',i5)") lu
        rmcomments = -1
        return
      endif
      if (len_trim(comcharin) > 0) then
        comchar = comcharin
      else
        comchar = ';'
        write(6,"('rmcomments: using default semicolon as ',
     +    'comment character.')")
      endif
      inquire(unit=lu,opened=isopen)
      if (.not.isopen) then
        open(unit=lu,iostat=ios)
        if (ios /= 0) then
          write(6,"('>>> WARNING rmcomments: error opening input',
     +      ' file with unit lu=',i2,' ios=',i5)") lu,ios
          rmcomments = -1
          return
        endif
      endif
      rmcomments = nextlu()
!     rewind lu ! not allowed to rewind unit 5 on sgi
      nline = 0
      read_loop: do
        line = ' '
        read(lu,"(a)",iostat=ios) line
        if (ios > 0) then
          write(6,"('>>> rmcomments: error reading from input',
     +      ' unit lu=',i3,' at line ',i5)") lu,nline
          return
        endif
        if (ios < 0) exit read_loop	! eof
        nline = nline+1
!
! Remove line if it has only "E" in column 1 (this was an
! old "Echo" directive from f77/cray namelist):
!
        if (line(1:1)=='E'.and.trim(line)=='E') cycle read_loop
!
! Use only non-commented part of line:
!
        compos = index(line,comchar)
        if (compos == 1) cycle read_loop
        if (compos > 0) line = line(1:compos-1)
        if (len_trim(adjustl(line))==0) cycle read_loop
!
! Write to new unit:
        write(rmcomments,"(a)") trim(line)
        if (echo > 0) write(6,"(a)") line(1:len_trim(line))
      enddo read_loop  
      rewind rmcomments
      return      
      end function rmcomments
!-------------------------------------------------------------------
      subroutine mkdiskflnm(msspath,diskname)
      implicit none
!
! Given a mss path (msspath), construct and return a local
! file name (diskname) which is the same as msspath except
! that slashes ('/') in the msspath are replaced by dots ('.')
! in the diskname. (The initial slash in msspath is ignored)
!
! For example, if on input msspath = '/FOSTER/dir1/file', then
! diskname would be returned as 'FOSTER.dir1.file'
!
! Args:
      character(len=*),intent(in) :: msspath
      character(len=*),intent(out) :: diskname
!
! Local:
      integer :: lmsspath,i,ii
!
      lmsspath = len_trim(msspath)
      if (lmsspath==0) then
        write(6,"('WARNING mkdiskflnm: zero length msspath.',
     +    ' Returning diskname: ''disk.file''')")
        diskname = 'disk.file'
        return
      endif
      diskname = ' '
      do i=1,len(diskname)
!
! Do not replace leading slash of msspath.
! If there are no occurrences of '/' in the msspath (i.e., its
!   not really an mss path, but just a file name), return the
!   msspath unchanged.
!                 
        ii = i+1
        if (index(msspath,'/')==0) ii = i
        if (ii > lmsspath) exit
        if (msspath(ii:ii)=='/') then
          diskname(i:i)='.'
        else
          diskname(i:i) = msspath(ii:ii)
        endif
      enddo
      if (len_trim(diskname)==0) then
        write(6,"('WARNING mkdiskflnm: zero length diskname output.')")
        write(6,"('  msspath=',a)") trim(msspath)
        write(6,"('Returning diskname: ''disk.file''')")
        diskname = 'disk.file'
      endif
!     write(6,"('mkdiskflnm returning:',/,'  msspath=',a,
!    +  /,'  diskname=',a)") trim(msspath),trim(diskname)
      end subroutine mkdiskflnm
!-------------------------------------------------------------------
      integer function isystem(command)
      implicit none
!
! Execute command to the shell. UNICOS and SGI use ishell(),
! AIX uses system().
!
      character(len=*),intent(in) :: command
#if defined(UNICOS) || defined(SGI)
      integer,external :: ishell
      isystem = ishell(trim(command))
#elif defined (AIX)
      integer,external :: system
      isystem = system(trim(command)//"\0")
#elif defined(OSF1)
      integer,external :: system
      isystem = system(trim(command))
#elif defined(SUN) || defined(LINUX)
      integer,external :: system
      isystem = system(trim(command))
#endif
      end function isystem
!-------------------------------------------------------------------
      integer function iunlink(file,iprint)
      implicit none
      character(len=*),intent(in) :: file
      integer,intent(in) :: iprint
      integer,external :: unlink
!
      iunlink = 1 ! to avoid ifort warning
#if defined(UNICOS) || defined(SGI) || defined(SUN) || defined(OSF1)
      iunlink = unlink(trim(file))
      if (iunlink.eq.0) then
        if (iprint > 0)
     |    write(6,"('Unlinked file ',a)") trim(file)
      else
        if (iprint > 0)
     |    write(6,"('Note: unlink of ',a,' failed',
     |      ' (possibly non-existant file).')") trim(file)
      endif
#elif AIX
      iunlink = unlink(trim(file)//"\0")
      if (iunlink.eq.0) then
        if (iprint > 0)
     |    write(6,"('Unlinked file ',a)") trim(file)
      else
        if (iprint > 0)
     |    write(6,"('Note: unlink of ',a,' failed',
     |      ' (possibly non-existant file).')") trim(file)
      endif
#endif
      end function iunlink
!-------------------------------------------------------------------
      integer function ilink(from,to,iprint)
      implicit none
!
! Args:
      character(len=*),intent(in) :: from,to
      integer,intent(in) :: iprint
!
! External:
      integer,external :: link
!
      ilink = 1 ! to avoid ifort warning
#if defined(UNICOS) || defined(SGI) || defined(SUN) || defined(OSF1)
      ilink = link(trim(from),trim(to))
      if (ilink.eq.0) then
        if (iprint > 0)
     |  write(6,"('Linked file ',a,' to ',a)") trim(from),trim(to)
      else
        write(6,"(/,'>>> WARNING: link command of ',a,' to ',a,
     |    ' failed: ierror=',i3,/)") trim(from),trim(to),ilink
!       stop 'ilink'
      endif
#elif defined(AIX)
!
! Note: link will fail if executing from different file system
!       than tmpdir (e.g., executing from /home but tmpdir is
!       in /ptmp)
!
      ilink = link(trim(from)//"\0",trim(to)//"\0")
      if (ilink.eq.0) then
        if (iprint > 0)
     |    write(6,"('Linked file ',a,' to ',a)") trim(from),trim(to)
      else
        write(6,"('>>> WARNING: link command of ',a,' to ',a,
     |    ' failed: ierror=',i3)") trim(from),trim(to),ilink
      endif
#endif
      end function ilink
!-------------------------------------------------------------------
      subroutine datetime(curdate,curtime)
!
! Return character*8 values for current date and time.
! (sub date_and_time is an f90 intrinsic)
!
      implicit none
!
! Args:
      character(len=*),intent(out) :: curdate,curtime
!
! Local:
      character(len=8) :: date
      character(len=10) :: time
      character(len=5) :: zone
      integer :: values(8)
!
      curdate = ' '
      curtime = ' '
      call date_and_time(date,time,zone,values)
!
!     write(6,"('datetime: date=',a,' time=',a,' zone=',a)")
!    |  date,time,zone
!     write(6,"('datetime: values=',8i8)") values
!
      curdate(1:2) = date(5:6)
      curdate(3:3) = '/'
      curdate(4:5) = date(7:8)
      curdate(6:6) = '/'
      curdate(7:8) = date(3:4)
!
      curtime(1:2) = time(1:2)
      curtime(3:3) = ':'
      curtime(4:5) = time(3:4)
      curtime(6:6) = ':'
      curtime(7:8) = time(5:6)
!
      end subroutine datetime
!-------------------------------------------------------------------
      subroutine handle_ncerr(istat,msg)
      implicit none
      include 'netcdf.inc'
!
! Handle a netcdf lib error:
!
      integer,intent(in) :: istat
      character(len=*),intent(in) :: msg
!
      write(6,"(/72('-'))")
      write(6,"('>>> Error from netcdf library:')")
      write(6,"(a)") trim(msg)
      write(6,"('istat=',i5)") istat
      write(6,"(a)") nf_strerror(istat)
      write(6,"(72('-')/)")
      return
      end subroutine handle_ncerr
!-------------------------------------------------------------------
      subroutine scpfile(lu,local,remot)
!
! Send file local to remot via scp.
! remot is machine:path specification.
! If lu is > 0 and <= 99 it is closed before doing the scp.
!
      implicit none
!
! Args:
      integer,intent(in) :: lu
      character(len=*),intent(in) :: local,remot
!
! Locals:
      character(len=124) :: command = ' '
      character(len=80) :: msserr = ' '
      logical :: isopen,exists
      integer :: istat
!
! Externals:
      integer,external :: isystem
!
      if (lu.gt.0.and.lu.le.99) then
        inquire(lu,opened=isopen)
        if (isopen) then
          write(6,"('scpfile: closing unit ',i2)") lu
          close(lu)
        endif
      endif
      inquire(file=local,exist=exists)
      if (.not.exists) then
        write(6,"('scpfile: cannot find file ',a,' -- no scp done')")
     +    local
        return
      endif
      command = ' '
!
! 5/6/04 btf:
! The -i option to scp specifies a private key file for 
! unattended scp. Using scp without -i does not work, even 
! if ~/.ssh/id_dsa exists and works from the command line. 
! It does work if the -i option is given and the arg to -i exists.
!
      write(command,"('scp -i ~/id_dsa ',a,' ',a)") 
     |  trim(local),trim(remot)

      write(6,"(/,'Tgcmproc using scp to send ',a,' to ',a)")
     |  trim(local),trim(remot)
      write(6,"('Command = ',a)") trim(command)
      command = trim(command)
      istat = isystem(command)

      if (istat.eq.0) then
        write(6,"('Copied file ',a,' to remote ',a)") 
     +    trim(local),trim(remot)
      else
        write(6,"(/72('>'))")
        write(6,"('AN ERROR HAS OCCURRED USING SCP TO TRANSFER ')")
        write(6,"('LOCAL FILE ',a,' TO REMOTE FILE (OR DIR) ',a)")
     +    trim(local),trim(remot)
        write(6,"(' ')")
        write(6,"('The command that failed was as follows:',/a)")
     +    trim(command)
        write(6,"('Note that for this unattended scp to work, you ',
     |    'must have a private key file in ~/id_dsa')")
        write(6,"(' ')")
      endif
      return
      end
!-----------------------------------------------------------------------
      subroutine getcwd(cwd)
!
! Return current working directory in cwd:
!
      implicit none
      character(len=*),intent(out) :: cwd
      integer,external :: isystem,nextlu
      integer :: istat,lu,ier
!
      istat = isystem('pwd > cwd')
      if (istat /= 0) write(6,"('>>> WARNING getcwd: error return ',
     |  'isystem(''pwd > cwd''): istat=',i8)") istat
      lu = nextlu()
      open(lu,file='cwd',status='old')
      read(lu,"(a)") cwd
      close(lu) 
      end subroutine getcwd
!-----------------------------------------------------------------------
      subroutine check_nans(f,id1,id2,id3,name,n_total,ispval,spval,
     |  iprint,ifatal)
!
! Check for existence of +/-INF and NaN's in field f(id1,id2,id3).
!   If ispval > 0 -> replace any INF or NaNs with spval
!   If iprint==1  -> print warnings only if INF or NaNs are found
!   If iprint==2  -> always print number of INF and NaNs found
!   If ifatal > 0 -> stop program when first INF or NaNs are found
! Note: Can distinguish between +/-INF (not really NaNs), but cannot 
!       distinguish between types of actual NaNs (+/-NaNQ and NaNS).
! IBM only. See pp318-319 User's Guide Version 8 XL Fortran for AIX
!
      implicit none
!
! Args:
      integer,intent(in) :: id1,id2,id3,iprint,ifatal,ispval
      integer,intent(out) :: n_total ! total number of +/-INF+NaNs
      real,intent(inout) :: f(id1,id2,id3)
      real,intent(in) :: spval
      character(len=*),intent(in) :: name 
!
! Local:
      real :: plus_inf,minus_inf,plus_nanq,minus_nanq,sig_nan
!
! For double precision 8-byte reals (-qrealsize=8):
!     data plus_inf   /z'7ff0000000000000'/  ! INF   (overflow)
!     data minus_inf  /z'fff0000000000000'/  ! -INF  (underflow)
!     data plus_nanq  /z'7ff8000000000000'/  ! NaNQ  (plus quiet NaN)
!     data minus_nanq /z'fff8000000000000'/  ! -NaNQ (minus quiet NaN)
!     data sig_nan    /z'7ff0000000000001'/  ! NaNS  (signalling NaN)
!
! For single precision (4-byte) reals:
      data plus_inf   /z'7f800000'/  ! INF   (overflow)
      data minus_inf  /z'ff800000'/  ! -INF  (underflow)
      data plus_nanq  /z'7fc00000'/  ! NaNQ  (plus quiet NaN)
      data minus_nanq /z'ffc00000'/  ! -NaNQ (minus quiet NaN)
      data sig_nan    /z'7f800001'/  ! NaNS  (signalling NaN)
!
      integer :: i1,i2,i3
      integer :: 
     |  n_plus_inf,   ! number of INF
     |  n_minus_inf,  ! number of -INF
     |  n_nan         ! total number of NaNs (+/-NaNQ and NaNS)
!
! Init:
      n_plus_inf = 0
      n_minus_inf = 0
      n_nan = 0
      n_total = 0
!
! Scan array:
      do i3=1,id3
        do i2=1,id2
!
! +/-INF are detected by simple comparison:
          n_plus_inf   = n_plus_inf   + count(f(:,i2,i3)==plus_inf) 
          n_minus_inf  = n_minus_inf  + count(f(:,i2,i3)==minus_inf) 
!
! NaNs (NaNQ or NaNS) are detected by (a /= a):
          n_nan        = n_nan        + count(f(:,i2,i3)/=f(:,i2,i3))
          n_total = n_plus_inf+n_minus_inf+n_nan
!
!         write(6,"('i3=',i3,' i2=',i3,' n_plus_inf=',i8,' n_minus_inf='
!    |      ,i8,' n_nan=',i8,' n_total=',i8)") i3,i2,n_plus_inf,
!    |      n_minus_inf,n_nan,n_total
!
! Fatal when first INF or NaN is found:
          if (ifatal > 0 .and. n_total > 0) then
            write(6,"(/,'>>> FATAL: Found INF and/or NaNs in field ',
     |        a)") name
            write(6,"('  Dimensions id1,id2,id3=',3i4)") id1,id2,id3
            write(6,"('  First INF or NaN found at i2=',i4,', i3=',i4)") 
     |        i2,i3
            write(6,"('  n_plus_inf   = ',i6)") n_plus_inf
            write(6,"('  n_minus_inf  = ',i6)") n_minus_inf
            write(6,"('  n_nan (NaNS or NaNQ) = ',i6)") n_nan
            write(6,"('  data(:,',i3,',',i3,') = ',/,(6e12.4))") 
     |        i2,i3,f(:,i2,i3)
            stop 'check_nans'
          endif ! ifatal > 0
!
! Replace any INF or NaNs with spval:
          if (ispval > 0 .and. n_total > 0) then
            do i1=1,id1
              if (f(i1,i2,i3)==plus_inf.or.f(i1,i2,i3)==minus_inf.or.
     |            f(i1,i2,i3)/=f(i1,i2,i3)) f(i1,i2,i3) = spval
            enddo
          endif 
        enddo ! i2=1,id2
      enddo ! i3=1,id3
!
! Print level 1 (print warnings only if INF or NaNs are found):
      if (iprint==1) then
        if (n_plus_inf > 0) write(6,"('>>> WARNING: found ',
     |    i6,' INF values in field ',a,' (id1,2,3=',3i4,')')") 
     |    n_plus_inf,name,id1,id2,id3
        if (n_minus_inf > 0) write(6,"('>>> WARNING: found ',
     |    i6,' -INF values in field ',a,' (id1,2,3=',3i4,')')") 
     |    n_minus_inf,name,id1,id2,id3
        if (n_nan > 0) write(6,"('>>> WARNING: found ',i6,
     |    ' NaNS or NaNQ values in field ',a,' (id1,2,3=',3i4,')')") 
     |    n_nan,name,id1,id2,id3
!       if (ispval > 0 .and. n_total > 0) 
!    |    write(6,"('>>> Replaced ',i8,' values with spval ',e12.4)")
!    |      n_total,spval
!
! Print level 2 (always print number of nans found):
      elseif (iprint==2) then 
        write(6,"('Checking for INF and NaNs in field ',a,' id1,2,3=',
     |    3i4)") name,id1,id2,id3
        print *,'  n_plus_inf   (',plus_inf,  ') = ',n_plus_inf
        print *,'  n_minus_inf  (',minus_inf, ') = ',n_minus_inf
        print *,'  n_nan        (',plus_nanq,'+',sig_nan,') = ',n_nan
        print *,'  n_total      (total INF+NaNs) = ',n_total
!       if (ispval > 0)
!    |  print *,'  Replaced ',n_total,' values with spval ',spval
      endif
      end subroutine check_nans
!
!-----------------------------------------------------------------------
      subroutine smoother(f,nf,npts,spv,iprint)
      implicit none
!
! Args:
      integer,intent(in) :: nf,npts,iprint
      real,intent(inout) :: f(nf)
      real,intent(in) :: spv
!
! Local:
      integer :: i,nspv
!
      nspv = 0
!
! 3-point smoother:
      if (npts==3) then
        do i=2,nf-1
          if (f(i-1)/=spv.and.f(i)/=spv.and.f(i+1)/=spv) then
            f(i) = (f(i-1) + 2.*f(i) + f(i+1)) / 4.
          else
            nspv = nspv+1
          endif
        enddo ! i=2,nf-1
      elseif (npts==5) then
!
! 5-point smoother:
        do i=3,nf-2
          if (f(i-2)/=spv.and.f(i-1)/=spv.and.f(i)/=spv.and.
     |        f(i+1)/=spv.and.f(i+2)/=spv) then
            f(i) = (f(i-2) + f(i-1) + f(i) + f(i+1) + f(i+2)) / 5.
          else
            nspv = nspv+1
          endif
        enddo ! i=2,nf-1
!
! Only 3 and 5-point smoothers currently available:
      else
        write(6,"('>>> WARNING smoother: bad npts=',i4)") npts
      endif 
!
      if (nspv > 0 .and. iprint > 0) then
        write(6,"('Smoother: nspv=',i5)") nspv
      endif
      end subroutine smoother
!-----------------------------------------------------------------------
      recursive subroutine expand_path(path)
!
! Expand any environment variables in path, return expanded path.
!   If '$' is found, then an env var is defined as that part of path
!   following the '$' until the next slash. If no '$' is found,
!   return without changing path.
! Example:
!   Here, "TGCMDATA" is considered the env var:
!   path = '$TGCMDATA/dir1/file'
! This routine is recursive, so multiple env vars can be used, e.g.:
!   path = '$MYDIR/$MYSUBDIR/file.nc'
!
      implicit none
!
! Args:
      character(len=*),intent(inout) :: path
!
! Local:
      character(len=224) :: path_out,envvar_value
      character(len=80) :: envvar_name
      integer :: i,idollar,islash
!
      if (len_trim(path)==0) then
        write(6,"('>>> WARNING expand_path: path is empty.')")
        return
      endif
!
!     write(6,"('Enter expand_path: path=',a)") trim(path)
!
      idollar = index(path,'$')
      if (idollar <= 0) return ! no env var in path
!
! Env var is between idollar and next slash 
! (or end of path if there is no slash after idollar):
!
      islash = index(path(idollar+1:len_trim(path)),'/')
      if (islash <= 0) islash = len_trim(path)
      islash = islash+idollar
      envvar_name = path(idollar+1:islash-1)
!
! Get value of env var (getenv is f90 intrinsic):
      call getenv(trim(envvar_name),envvar_value)
      if (len_trim(envvar_value) <= 0) then
        write(6,"('>>> WARNING expand_path: error retrieving ',
     |    'value for env var ',a)") trim(envvar_name)
        return
      endif
!
! Put together expanded output path:
      if (idollar > 1) then
        if (islash < len_trim(path)) then
          path_out = path(1:idollar-1)//trim(envvar_value)//
     |      path(islash:len_trim(path))
        else
          path_out = path(1:idollar-1)//trim(envvar_value)
        endif
      else     ! idollar == 1
        if (islash < len_trim(path)) then
          path_out = trim(envvar_value)//path(islash:len_trim(path))
        else
          path_out = trim(envvar_value)
        endif
      endif
!
! Return new path, and make recursive call for more env vars:
      path = trim(path_out)
      call expand_path(path) ! expand next env var
!
      end subroutine expand_path
!-------------------------------------------------------------------
      integer function isymlink(from,to,iprint)
      implicit none
!
! Args:
      character(len=*),intent(in) :: from,to
      integer,intent(in) :: iprint
!
! Local:
      logical :: exists
      integer :: istat
      integer,parameter :: maxlen=1024
      character(len=maxlen) :: command
!
! External:
      integer,external :: isystem
!
! Check for existence of "from" file:
      inquire(file=trim(from),exist=exists)
!     write(6,"('isymlink: Linking from=',a,' to=',a,' exists=',l1)")
!    |  trim(from),trim(to),exists
!
! Make symbolic link:
      write(command,"('ln -s -f ',a,' ',a)") trim(from),trim(to)
      isymlink = isystem(command)
      if (isymlink == 0) then
        if (iprint > 0)
     |    write(6,"('Made symbolic link of ',a,' to ',a)")
     |      trim(from),trim(to)
      else
        write(6,"('>>> WARNING: symbolic link of ',a,' to ',a,
     |    ' failed: ierror=',i3,/)") trim(from),trim(to),isymlink
!       call shutdown('isymlink')
      endif
      end function isymlink
!-----------------------------------------------------------------------
      integer function mssrcp(opts,src,dest)
      implicit none
!
! Build and execute msrcp command to the shell. Assume "mss:" prefixes
! either src (for a get) or dest (for a put).
!
! Args:
      character(len=*),intent(in) :: opts,src,dest
!     
! Local:
      character(len=1024) :: cmd
      integer :: ier
!
! External:
      integer,external :: isystem
!
      if (src(1:4) /= "mss:" .and. dest(1:4) /= "mss:") then
        write(6,"(/,'>>> mssrcp: need mss: prefix on either src or ',
     |    'dest')")
        write(6,"('src=',a)") trim(src)
        write(6,"('dest=',a)") trim(dest)
        call shutdown('mssrcp')
      endif
!
! Only master task does the msrcp:
      if (src(1:4) == "mss:") then
        write(6,"('mssrcp: Obtaining file ',a,' from mss.')")
     |    src(5:len_trim(src))
      else
        write(6,"('mssrcp: Disposing file ',a,' to mss.')")
     |    dest(5:len_trim(dest))
      endif
      cmd = ' '         
      cmd = "msrcp "//trim(opts)//' '//trim(src)//' '//trim(dest)
      write(6,"('mssrcp: Executing ',a)") trim(cmd)
      mssrcp = isystem(cmd)
      write(6,"(a)") trim(cmd)
      end function mssrcp
!-----------------------------------------------------------------------
      logical function is_mspath(path)
!
! Return true if input path begins with "/LOGNAME/" or "/TGCM/",
! otherwise return false.
!
      implicit none
!
! Args:
      character(len=*),intent(in) :: path
!      
! Local:
      character(len=64) :: logname,lognamedir
!       
! External:
      character(len=64),external :: trcase ! (util.F)
!
      is_mspath = .false.
      if (path(1:1) /= '/') return
!
! Get logname in upper case:
      logname = ' '
      call getenv('LOGNAME',logname) ! get logname from env var
      if (len_trim(logname)==0) then
        write(6,"('>>> WARNING is_mspath: Cannot get LOGNAME ',
     |    'environment variable.')")
        return
      elseif (len_trim(logname) > 62) then
        write(6,"('>>> WARNING is_mspath: please lengthen logname',
     |    ' (must be at least 62: len_trim(logname)=',i3)")
     |    len_trim(logname)
        return
      endif
      logname = trcase(logname) ! translate from lower to upper case
!
      lognamedir = '/'//trim(logname)//'/'
      if (path(1:len_trim(lognamedir))==lognamedir) is_mspath = .true.
      if (is_mspath) return
!
! Try /TGCM/:
      lognamedir = '/'//'TGCM'//'/'
      if (path(1:len_trim(lognamedir))==lognamedir) is_mspath = .true.
!
      end function is_mspath
!-----------------------------------------------------------------------
      subroutine check_mspath(mspathin,mspathout)
!
! Check validity of mspathin as an mss path 
! (must begin with /LOGNAME or /TGCM)
!
      implicit none
!
! Args:
      character(len=*),intent(in) :: mspathin
      character(len=*),intent(out) :: mspathout
!
! Local:
      integer :: i,ii,istat,ipos,oldlen,newlen
      character(len=32) :: logname
      character(len=len(logname)),external :: trcase
!
! Init:
!     write(6,"('Enter check_mspath: mspathin=',a)") trim(mspathin)
      mspathout = mspathin
!
! Mss path must begin with a slash:
      if (mspathin(1:1) /= '/') then
        if (len(mspathout) < len_trim(mspathin)+1) then
          write(6,"('>>> WARNING check_mspath: length of mspathout',
     |      ' too small: ',i4, ' (must be at least ',i4,')')") 
     |      len(mspathout),len_trim(mspathin)+1
        endif
        mspathout = ' '       ! reinit
        mspathout(1:1) = '/'
        do i=2,len_trim(mspathin)+1
          mspathout(i:i) = mspathin(i-1:i-1)    
        enddo
!       write(6,"('check_mspath: added slash to front of mspath: ',a)")
!    |    trim(mspathout)
      endif
!
! Get logname in upper case:
      logname = ' '
      call getenv('LOGNAME',logname) ! get logname from env var
      if (len_trim(logname)==0) then
        write(6,"('>>> WARNING check_mspath: Cannot get LOGNAME ',
     |    'environment variable.')")
        return
      endif
      logname = trcase(logname) ! translate from lower to upper case
!
! Mss path must begin with /LOGNAME:
      if (mspathout(2:len_trim(logname)+1) /= trim(logname).and.
     |    mspathout(2:5) /= 'TGCM') then
        if (len(mspathout) < len_trim(mspathout)+len_trim(logname)+1) 
     |    then
          write(6,"('>>> WARNING check_mspath: length of mspathout',
     |      ' not long enough to add /',a)") trim(logname)
          return
        endif 
        oldlen = len_trim(mspathout)
        newlen = oldlen+len_trim(logname)+1 ! +1 for leading slash
        do i=newlen,newlen-oldlen+1,-1
          ii = oldlen-newlen+i
          mspathout(i:i) = mspathout(ii:ii) ! copy existing part
        enddo
        mspathout(1:len_trim(logname)+1) = '/'//trim(logname) ! add /LOGNAME
!       write(6,"('check_mspath: added /',a,' to front of mspath: ',a)")
!    |    trim(logname),trim(mspathout)
      endif
!     write(6,"('check_mspath returning: mspathout=',a)") 
!    |  trim(mspathout)
      end subroutine check_mspath
!-----------------------------------------------------------------------
      character(len=*) function trcase(string)
!
! Translate case of input string string, i.e., return string like string,
! but with lower case character converted to upper case, and vice-versa.
!
      implicit none
!
! Args:
      character(len=*),intent(in) :: string
!
! Local:
      character(len=*),parameter :: ucase = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
      character(len=*),parameter :: lcase = "abcdefghijklmnopqrstuvwxyz"
      integer :: i,lpos,upos
!
      trcase = ' '
      do i=1,len_trim(string)
        trcase(i:i) = string(i:i)
        lpos = index(lcase,trcase(i:i))
        upos = index(ucase,trcase(i:i))
        if (lpos > 0) then
          trcase(i:i) = ucase(lpos:lpos) ! lower to upper
        elseif (upos > 0) then
          trcase(i:i) = lcase(upos:upos) ! upper to lower
        endif
      enddo
      end function trcase
!-----------------------------------------------------------------------
      subroutine shutdown(msg)
!
! An fatal error has occurred -- shut down the processor.
!
      implicit none
!
! Args:
      character(len=*) :: msg
!
! Local:
      integer :: ier
      character(len=80) :: errorcode
!
      write(6,"(/,28('>'),' PROCESSOR SHUTDOWN ',28('<'))")
      write(6,"('Shutdown: stop message: ',a)") trim(msg)
      stop 'shutdown'
      end subroutine shutdown
!-----------------------------------------------------------------------
      subroutine rotate_periodic(fin,fout,nf,npts)
      implicit none
!
! Rotate periodic array fin(nf) npts indices.  Return rotated array
! in fout(nf)
!
! Args:
      integer,intent(in) :: nf,npts
      real,intent(in) :: fin(nf) 
      real,intent(out) :: fout(nf) 
!
! Local:
      integer :: i,ii
      real :: ftmp(nf) 
!
      if (nf <= 0) then
        write(6,"('>>> rotate_periodic: nf must be >= 0: nf=',i5)") nf
        stop 'rotate_periodic'
      endif
      if (npts < 0 .or. npts > nf) then
        write(6,"('>>> rotate_periodic: npts must be > 0 and <= nf: ',
     |    ' nf=',i5,' npts=',i5)") nf,npts
        stop 'rotate_periodic'
      endif
      fout = fin ! init to input field
      if (npts==0) return
      do i=1,nf
        ii = mod(i+npts,nf)
        if (ii==0) ii=nf
!       write(6,"('i=',i3,' ii=',i3,' nf=',i3)") i,ii,nf
        fout(i) = fin(ii)
      enddo ! i=1,nf
      end subroutine rotate_periodic
!-----------------------------------------------------------------------
      subroutine plch(xpos,ypos,chrs,size,angd,cntr)
      use input,only: label_quality
      implicit none
!
! Args:
      real,intent(in) :: xpos,ypos,size,angd,cntr
      character(len=*),intent(in) :: chrs

      if (label_quality(1:3)=='low') then
        call plchlq(xpos,ypos,chrs,size,angd,cntr)
      elseif (label_quality(1:3)=='med') then
        call plchmq(xpos,ypos,chrs,size,angd,cntr)
      elseif (label_quality(1:3)=='hig') then
        call plchhq(xpos,ypos,chrs,size,angd,cntr)
      else
        write(6,"('>>> Warning plch: unknown label_quality=',a)")
     |    label_quality
      endif
      end subroutine plch
