!
! 9/98: For SGI: surfn1 and surfd are not in SGI fitpack lib
! (it does have surf1 and surf2)
!
      module qbary
      use proc,only: spval,dlon,p0,gcmlat,nlat
      use hist,only: history
      use fields,only: field
      implicit none
      real,parameter,private ::
     +  g=9.8,			! gravity
     +  cp=1004.,		! specific heat
     +  ts=300.,		! surface temperature
     +  pi=3.1415,		! pi
     +  r=287.,			! scale height reference constant
     +  ps=1000.,               ! 
     +  dtr=pi/180.,            ! degree to radian
     +  omega=7.292e-5,
     +  re=6370e3
!
! Zonal means are saved by save_qbary:
      real,allocatable ::
     +  tzm(:,:),uzm(:,:),      ! zonal mean t,u (jmx,kmx)
     +  rhozm(:,:),             ! zonal mean rho (jmx,kmx) (kg/m3)
     +  pottn(:,:)              ! potential temperature (jmx,kmx)
      logical mkqbary
#if defined(SUN) || defined(LINUX)
      end module qbary
#else
      contains
!-------------------------------------------------------------------
      subroutine save_qbary(flat,rho,h,imx,kmx,jmx,glat,lat)
!
! Save 2d zonal means (jmx,kmx) for later use by calc_qbary.
! (This routine called for each latitude at index lat)
!
! Args:
      integer,intent(in) :: imx,kmx,jmx,lat
      type(history),intent(in) :: h
      real,intent(in) :: flat(imx,kmx,h%nflds),rho(imx,kmx),glat
!
! Locals:
      integer,save :: ixt,ixu
      integer :: ier,k
      real :: pmb,zp,dlev
!
! Externals:
      real,external :: fmean
      integer,external :: ixfindc
!
! Assume 1st entry if tzm not allocated:
      if (.not.allocated(tzm)) then
!
! Save indices to needed fields in flat:
        ixt  = ixfindc(h%fnames,h%nflds,'TN      ')
        ixu  = ixfindc(h%fnames,h%nflds,'UN      ')
        if (ixt==0.or.ixu==0) then
          write(6,"('>>> WARNING save_qbary: need t and u for ',
     +      'qbary: ixt,u=',2i4)") ixt,ixu
          stop 'save_qbary'
        endif
!
! Allocate zonal mean fields:
        allocate(tzm(jmx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"save_qbary allocating tzm")
        allocate(uzm(jmx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"save_qbary allocating uzm")
        allocate(rhozm(jmx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"save_qbary allocating rhozm")
        allocate(pottn(jmx,kmx),stat=ier)
        if (ier /= 0) call allocerr(ier,"save_qbary allocating pottn")
      endif ! allocated(tzm)
!
! Save 2d zonal means t and u (jmx,kmx), for current lat:
! (Also save log potential temperature pottn)
!
      dlev = (h%zpt-h%zpb)/float(h%nzp-1)
      do k=1,kmx
        tzm(lat,k) = fmean(flat(:,k,ixt),imx,spval,0)
        uzm(lat,k) = fmean(flat(:,k,ixu),imx,spval,0)
        rhozm(lat,k) = fmean(rho(:,k),imx,spval,0)
        zp = h%zpb+(k-1)*dlev
        pmb = p0*exp(-zp)*1.e-3 ! pressure in mbars
        pottn(lat,k) = log(tzm(lat,k)*(ps/pmb)**(r/cp))
      enddo
      end subroutine save_qbary
!-------------------------------------------------------------------
      subroutine calc_qbary(f,nf,h,imx,kmx,jmx)
!
! Calculate qbary, return in f(ixqbary)%data:
! This routine modified from Scott Palo's "PROGRAM COMP_QBARY_TIMEGCM" 
!   and related subroutines COMPBVFREQ, COMPPOTTEMP, COMPPVGRAD.
! This routine called only once, after save_qbary has been called
!   at all latitudes.
!
! Args:
      integer,intent(in) :: 
     +  nf,    ! number of fields
     +  imx,   ! number of longitudes
     +  kmx,   ! number of pressures
     +  jmx    ! number of latitudes
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
!
! Locals:
      integer :: ixqbary,ixqh,ixqv,nwork,ier,i,k,j
      real,parameter :: sigma=1.
      real :: rlat(jmx),zp,dlev
      real :: zpart(jmx,kmx,3)    ! partial derivatives output by surfn1
      real :: zz(jmx,kmx),zx(jmx,kmx),zy(jmx,kmx),zxx(jmx,kmx),
     +  zxy(jmx,kmx),zyy(jmx,kmx) ! outputs from surfd
      real :: bvf(jmx,kmx)        ! brunt-viasala frequency
      real :: qh(jmx,kmx),qv(jmx,kmx),beta(jmx,kmx)
      real,allocatable :: 
     +  work(:),                  ! work array for surfn1
     +  zpht(:)                   ! height from zp(kmx)
!
! Externals:
      integer,external :: ixfindc
!
! Index to qbary in f:
      ixqbary   = ixfindc(f%fname8,nf,'QBARY   ')
      ixqh      = ixfindc(f%fname8,nf,'QH      ')
      ixqv      = ixfindc(f%fname8,nf,'QV      ')
      if (ixqbary+ixqh+ixqv == 0) return
!
! Latitudes in radians (array op):
      rlat = gcmlat*dtr
!
! Calculate height in meters from zp and reference scale height:
      allocate(zpht(kmx),stat=ier)
      if (ier /= 0) call allocerr(ier,"calc_qbary allocating zpht")
      dlev = (h%zpt-h%zpb)/float(h%nzp-1)
      do k=1,kmx
        zp = h%zpb+(k-1)*dlev ! log pressure, as in tgcm vertical grid
        zpht(k) = zp*r*ts/g
      enddo
!
! Work array for surfn1:
      nwork = max(jmx,kmx)
      allocate(work(nwork),stat=ier) 
      if (ier /= 0) call allocerr(ier,"calc_qbary allocating work") 
!
! Surfn1 and surfd are in fitpack (/usr/local/lib on the crays)
! Surfn1 here returns partial deriviatives of pottn in zpart:
! 
      call surfn1(jmx,kmx,rlat,zpht,pottn,jmx,zpart,work,sigma,ier)
      if (ier /= 0) then
        write(6,"('>>> calc_qbary: error return from surfn1 for ',
     +    'pottn: ier=',i2)") ier
      endif
!
! Compute Brunt-Viasala Frequency (bvf) from vertical gradient of
! log potential temperature:
!
      do k=1,kmx
        do j=1,jmx
          call surfd(rlat(j),zpht(k),zz(j,k),zx(j,k),zy(j,k),zxx(j,k),
     +      zxy(j,k),zyy(j,k),jmx,kmx,rlat,zpht,pottn,jmx,zpart,sigma)
          bvf(j,k) = tzm(j,k)/ts*g*zy(j,k)
        enddo
      enddo
!
! Derivatives of zonal mean zonal wind:
      call surfn1(jmx,kmx,rlat,zpht,uzm,jmx,zpart,work,sigma,ier)
      if (ier /= 0) then
        write(6,"('>>> calc_qbary: error return from surfn1 for ',
     +    'uzm: ier=',i2)") ier
      endif
      do k=1,kmx
        do j=1,jmx
      write(6,"(/,'>>> qbary not available on the Sun.')")
      stop 'qbary on Sun'
          call surfd(rlat(j),zpht(k),zz(j,k),zx(j,k),zy(j,k),
     +      zxx(j,k),zxy(j,k),zyy(j,k),jmx,kmx,rlat,zpht,uzm,
     +      jmx,zpart,sigma) 
!
! qh = horizontal contribution to pv.
! qv = vertical contribution to pv.
!
          qh(j,k) = (-1./re**2)*((-1./cos(rlat(j))**2)*uzm(j,k)-
     +      tan(rlat(j))*zx(j,k)+zxx(j,k))
          qv(j,k) = (rhozm(j,k)/bvf(j,k))*zy(j,k)
        enddo
      enddo
!
! Derivatives of intermediate qv:
      write(6,"(/,'>>> qbary not available on the Sun.')")
      stop 'qbary on Sun'
      call surfn1(jmx,kmx,rlat,zpht,qv,jmx,zpart,work,sigma,ier)
      if (ier /= 0) then
        write(6,"('>>> calc_qbary: error return from surfn1 for ',
     +    'qv: ier=',i2)") ier
      endif
      do k=1,kmx
        do j=1,jmx
          call surfd(rlat(j),zpht(k),zz(j,k),zx(j,k),zy(j,k),
     +      zxx(j,k),zxy(j,k),zyy(j,k),jmx,kmx,rlat,zpht,qv,
     +      jmx,zpart,sigma) 
          qv(j,k)=((-1.*(2.*omega*sin(rlat(j)))**2)/rhozm(j,k))*zy(j,k)
          beta(j,k) = (2.*omega/re)*cos(rlat(j)) ! redundant in vertical
        enddo
      enddo
!
! qbary (redundant in longitude):
! (note f%data(i,:,:), beta, qh, and qv are all (jmx,kmx))
      if (ixqbary > 0) then
        do i=1,imx
          f(ixqbary)%data(i,:,:) = beta+qh+qv
        enddo
      endif
!
! Save qh and qv if requested (also redundant in longitude):
      if (ixqh > 0) then
        do i=1,imx
          f(ixqh)%data(i,:,:) = qh
        enddo
      endif
      if (ixqv > 0) then
        do i=1,imx
          f(ixqv)%data(i,:,:) = qv
        enddo
      endif
!
! Release space (work and zpht are local, rest are module):
      deallocate(work)
      deallocate(zpht)
      deallocate(tzm)
      deallocate(uzm)
      deallocate(rhozm)
      deallocate(pottn)
      end subroutine calc_qbary
      end module qbary
#endif
