c    @(#)bqgofk.f	1.3
      real function bqgofk(ip,k)
c
c
c
c leonard walstad
c
c
c
c-----------------------------------------------------------------------
c
c
      include'mladv.inc'
c
c
      real mltlst
c
c-----------------------------------------------------------------------
c
      common /threed/ hz(k0), scsig, sigz(k01), eigval(k0), ainv(k0,k0),
     .                amat(k0,k0)

      common/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +     topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +     prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)

      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva
c
      common /extpar/ rlat0,rlng0,t0,v0,dhor,ht,time0,r0,f0star
c
      dimension dbdz(k01)
c
      data icall/0/
c
c_______________________________________________________________________
c
c
      if(icall.ne.0)go to 99
        icall=1
        iplast=-99
        mltlst=-9999.9
        scale=f0star*v0*dhor/(ht*scsig*sigz(1))
        do 1 kqg=1,kz+1
           dbdz(kqg)=(f0star*dhor/ht)**2./(scsig*sigz(kqg))
  1      continue
c
c
c_______________________________________________________________________
c
c
c
   99 continue
      if(ip.ne.iplast)go to 101
      if(mltime.eq.mltlst)go to 110
  101   continue
        delta=topden(ip) + artop
c
c
  110 continue

      bqgofk = buoysur + scale * delta

      z=zml(k)

      if (ifmqg .eq. 0) then
         bqgofk = bqgofk + z * dbdz(1)
         goto 3
        endif

      klev = 1
      zlast = 0.0
  2   continue
      if (z .lt. -zqg(klev)) then
         bqgofk = bqgofk + (-zqg(klev) - zlast) * dbdz(klev)
         zlast = -zqg(klev)
         klev = klev + 1
         goto 2
       endif

      bqgofk = bqgofk + (z - zlast) * dbdz(klev) 
c
c
  3   continue

      mltlst=mltime
      iplast=ip
      end
