!
! 3/10/94:
! Taken from doptuv.f (this dir), and changed name from redline to doppler.
! This routine takes ht,t,u,v, and ver as input (HM,TM,ZM,VM,VER,
!   dimensioned np (which is usually kmx in calling routines)), and
!   returns scalars doppler u,v,t (VELV,VELZ,ATTPV)
! Note input VER is not necessarilly from 6300A (was assumed to be in
!   doptuv.f)
! Keeping in quadrat at end for use in integrating vol emission rate (ver)
!   to obtain brightness (not used in doppler routine)
! Note: for redline emissions pass in sigma=15867.852
!       for greenline emissions pass  sigma=17924.
!
! 3/22/05:
! Adapted for new glbmean code.
!
!-----------------------------------------------------------------------
      subroutine doppler(np,hm,tm,zm,vm,ver,sigma,
     |  velv,velz,attpv,dattpv,attpz,dattpz)
      implicit none
!
! Args:
      integer,intent(in) :: np
      real,dimension(np),intent(in) :: hm,tm,zm,vm,ver
      real,intent(in) :: sigma
      real,intent(out) :: velv,velz,attpv,dattpv,attpz,dattpz
!
! Local:
      integer,parameter :: jmax=61
      real,dimension(np) :: ht,tt,zt,vt,bt,dl,dop
      real,dimension(jmax) :: as,az,azc,avc,we,aws,asw,wev,wez
      real,dimension(65) :: av
      integer :: i,ii,npn,k,j,jc,ji,iv,iz,ivs,izs
      real,parameter :: 
     |  pi     =3.141592654, 
     |  bb     =3.5811516E-7, 
     |  cspeed =3.0E10,
     |  amass  =16.0
      real :: r,alpha,pi180,b,radalph,rsina,rsinasq,rcosa,twor,sumb,
     |  sumbt,sumbz,sumbv,sumbl,sumbh,sumz,sumzb,doh,dlm,bd,dh,hh,dph,
     |  abt,abz,abv,abl,abh,abzh,sigmav,sigmaz,ddg,dgs,dgss,bsum,
     |  tbt,obt,atop,sigmahz,sigmahv,sigma0z,sigma0v,sigv1,sigv2,
     |  sigz1,sigz2,avz,avv,posz,posv,ddg2
      real :: crov,crb,slope,sb,cr1,crb1,sl1,sb1,r2, ! onepar output
     |        croz,cr

! Diffs w/ old code in last 3-4 decimals points (out of 12):
!     write(6,"('Enter doppler: hm =',/,(4e20.12))") hm
!     write(6,"('Enter doppler: tm =',/,(4e20.12))") tm
!     write(6,"('Enter doppler: zm =',/,(4e20.12))") zm
!     write(6,"('Enter doppler: vm =',/,(4e20.12))") vm
!     write(6,"('Enter doppler: ver=',/,(4e20.12))") ver
!
      r = 6370.0
      alpha = 20.0
      pi180 = pi/180.0
      b = ( bb * sigma)**2/(alog(2.0) * amass)

!     write(6,"('doppler: b=',e20.12)") b

      npn = np-1
      do i = 1,npn ! was do 2
        ii = i + 1
        ht(i) = (hm(i) + hm(ii))/2.0
        tt(i) = (tm(i) + tm(ii))/2.0
        zt(i) = (zm(i) + zm(ii))/2.0
        vt(i) = (vm(i) + vm(ii))/2.0
        bt(i) = (ver(i) + ver(ii))/2.0
      enddo ! i = 1,npn

!     write(6,"('doppler: bt=',/,(4e20.12))") bt
!
      radalph = alpha * pi180
      rsina = r * sin(radalph)
      rsinasq = rsina**2
      rcosa = r * cos(radalph)
      twor = 2.0 * r
      sumb  = 0.0
      sumbt = 0.0
      sumbz = 0.0
      sumbv = 0.0
      sumbl = 0.0
      sumbh = 0.0
      sumz  = 0.0
      sumzb = 0.0
!
      do i = 1,npn ! was do 3
        ii = i + 1
        dh = hm(ii)
        hh = hm(i)
        dph = sqrt( rsinasq + twor*dh + dh**2)
        doh = sqrt( rsinasq + twor*hh + hh**2)
        dlm = (dph + doh)/2.0 - rsina
        dl(i) = dph - doh
        dop(i) = rcosa/(r + ht(i))
        bd = bt(i) * dl(i)
        sumb  = sumb + bd
        sumbt = sumbt + bd*tt(i)
        sumbz = sumbz + bd*dop(i)*zt(i)
        sumbv = sumbv + bd*dop(i)*vt(i)
        sumbl = sumbl + bd*dlm
        sumbh = sumbh + bd*ht(i)
        bd = ( hm(ii) - hm(i)) * bt(i)
        sumz = sumz + bd*ht(i)
        sumzb = sumzb + bd
      enddo ! i = 1,npn

!     write(6,"('doppler: sumb,sumbt,sumbz =',3e20.12)") 
!    |  sumb,sumbt,sumbz
!     write(6,"('doppler: sumbv,sumbl,sumbh=',3e20.12)") 
!    |  sumbv,sumbl,sumbh
!
      abt = sumbt/sumb
      abz = sumbz/sumb
      abv = sumbv/sumb
      abl = sumbl/sumb
!     abh = sumbh/sumb
      abzh = sumz/sumzb
!     ahl = r*(sqrt(1.0 + (abl/r)**2 + 2.0*abl*rsina/r**2) - 1.0)
!     abzc = abz*(r+abzh)/rcosa
!     abvc = abv*(r+abzh)/rcosa
      sumb = sumb/10.
      sumzb = sumzb/10.
      sigmav = sigma * ( abv/cspeed)
      sigmaz = sigma * ( abz/cspeed)
      ddg = sqrt(b*abt*alog(2.0))/10.0

!     write(6,"('doppler: ddg=',e20.12)") ddg

      dgs = ddg*30.0
      jc = jmax/2
      dgss = -dgs
!
      do i=1,jmax ! was do 11
        av(i) = 0.0
        az(i) = 0.0
        as(i) = dgss
        dgss = dgss + ddg
      enddo ! i=1,jmax ! was do 11
!
      bsum = 0.0
      do k = 1,npn ! was do 12
        tbt = b * tt(k)
        obt = -1.0/tbt
        atop = sqrt(-obt/pi)* bt(k) * dl(k)
        bsum = bsum + atop
        sigmahz = sigma * (zt(k) * dop(k)/cspeed)
        sigmahv = sigma * (vt(k) * dop(k)/cspeed)
        sigma0z = sigmahz - sigmaz
        sigma0v = sigmahv - sigmav
        dgss = dgs
        j = jmax + 1
        do i = 1,jc ! was do 13
          ji = j - i
          sigv1 = exp(obt*( -dgss - sigma0v)**2)
          sigv2 = exp(obt*(  dgss - sigma0v)**2)
          sigz1 = exp(obt*( -dgss - sigma0z)**2)
          sigz2 = exp(obt*(  dgss - sigma0z)**2)
          av(i)  = av(i)  + atop * sigv1
          av(ji) = av(ji) + atop * sigv2
          az(i)  = az(i)  + atop * sigz1
          az(ji) = az(ji) + atop * sigz2
          dgss = dgss - ddg
        enddo ! i = 1,jc ! was do 13
        av( jc + 1 ) = av( jc + 1 ) + atop * exp(obt*(sigma0v**2))
        az( jc + 1 ) = az( jc + 1 ) + atop * exp(obt*(sigma0z**2))
      enddo ! k = 1,npn ! was do 12
!
      avz = -1.0e30
      avv = avz
      iv = jc + 1
      iz = iv
      posz = 0.0
!     posv = posz
      do i=1,jmax ! was do 15
        wev(i) = av(i)*0.129
        av(i) = alog(wev(i))
        wez(i) = az(i)*0.129
        az(i) = alog(wez(i))
        we(i) = 1.0
        if (av(i) > avv) then
          avv = av(i)
!         posv = as(i)
          iv = i
        endif
        if (az(i) > avz) then
          avz = az(i)
          posz = as(i)
          iz = i
        endif
      enddo ! i=1,jmax ! was do 15
      ivs = iv - 5
      ddg2 = ddg/2.0
      do i=1,10 ! was do 19
        aws(i) = (av(ivs+1) - av(ivs))/ddg

!       write(6,"('dopper: i=',i3,' ivs=',i3,' av(ivs+1)=',e20.12,
!    |    ' av(ivs)=',e20.12,' ddg=',e20.12,' aws=',e20.12)")
!    |    i,ivs,av(ivs+1),av(ivs),ddg,aws(i)

        asw(i) = as(ivs) + ddg2
        ivs = ivs + 1
      enddo ! i=1,10 ! was do 19

      call onepar(asw,aws,we,10,crov,crb,slope,sb,cr1,crb1,sl1,sb1,r2)

!     write(6,"('doppler after onepar: crov,crb,slope=',3e12.4)")
!    |  crov,crb,slope
!     write(6,"('doppler after onepar: sb,cr1,crb1   =',3e12.4)") 
!    |  sb,cr1,crb1
!     write(6,"('doppler after onepar: sl1,sb1,r2    =',3e12.4)")
!    |  sl1,sb1,r2

      crov = -crov/slope
      velv = (( sigmav+crov)*cspeed/sigma)*(r+abzh)/rcosa
      izs = iz - 5
      do i=1,10 ! was do 20
        aws(i) = ( az(izs + 1 ) - az(izs))/ddg
        asw(i) = as(izs) + ddg2
        izs = izs + 1
      enddo ! i=1,10 ! was do 20
      call onepar(asw,aws,we,10,croz,crb,slope,sb,cr1,crb1,sl1,sb1,r2)

!     write(6,"('doppler after onepar: croz,crb,slope=',3e12.4)")
!    |  croz,crb,slope
!     write(6,"('doppler after onepar: sb,cr1,crb1   =',3e12.4)") 
!    |  sb,cr1,crb1
!     write(6,"('doppler after onepar: sl1,sb1,r2    =',3e12.4)")
!    |  sl1,sb1,r2

      croz = -croz/slope
      velz = (( sigmaz+croz)*cspeed/sigma)*(r+abzh)/rcosa
      do i=1,jmax ! was do 18
        azc(i) = (as(i) + croz)**2
        avc(i) = (as(i) + crov)**2
      enddo ! i=1,jmax ! was do 18
      call onepar(azc,az,wez,jmax,cr,crb,slope,sb,cr1,crb1,sl1,sb1,r2)

!     write(6,"('doppler after onepar: cr,crb,slope=',3e12.4)")
!    |  cr,crb,slope
!     write(6,"('doppler after onepar: sb,cr1,crb1   =',3e12.4)") 
!    |  sb,cr1,crb1
!     write(6,"('doppler after onepar: sl1,sb1,r2    =',3e12.4)")
!    |  sl1,sb1,r2

      attpz = -1.0/(slope * b)
      dattpz = b*sb*attpz**2
      call onepar(avc,av,wev,jmax,cr,crb,slope,sb,cr1,crb1,sl1,sb1,r2)
      attpv = -1.0/(slope * b)
      dattpv = b*sb*attpv**2

      end subroutine doppler
!-----------------------------------------------------------------------
      subroutine onepar(x,y,n,np,cr,crb,sl,sb,cr1,crb1,sl1,sb1,r2)
      implicit none
!
! Args:
      integer,intent(in) :: np
      real,dimension(np),intent(in) :: x,y,n
      real,intent(out) :: cr,crb,sl,sb,cr1,crb1,sl1,sb1,r2
!
! Local:
      real :: sumx,sumxq,sumy,sumyq,sumxy,sumn,sumx2,sumy2,sumxty,
     |  xl,yl
      integer :: i

!
! Diffs w/ old code after around 12th decimal in x and y: 
!     write(6,"('Enter onepar: np=',i3,' x=',/,(4e20.12))") np,x
!     write(6,"('Enter onepar: np=',i3,' y=',/,(4e20.12))") np,y
!     write(6,"('Enter onepar: np=',i3,' n=',/,(4e20.12))") np,n

      sumx   = 0.0
      sumxq  = 0.0
      sumy   = 0.0
      sumyq  = 0.0
      sumxy  = 0.0
      sumn   = 0.0
      sumx2  = 0.0
      sumy2  = 0.0
      sumxty = 0.0
      sl   = 0.0
      cr   = 0.0
!     sx2  = 0.0
!     sy2  = 0.0
!     syx2 = 0.0
      r2   = 0.0
      sb   = 0.0
!     t    = 0.0
      do i = 1,np ! was do 1
        sumx = sumx + x(i) * n(i)

!       write(6,"('onepar: i=',i3,' x=',e20.12,' n=',e12.4,' sumx=',
!    |    e20.12)") i,x(i),n(i),sumx

        sumy = sumy + y(i) * n(i)
        sumxy = sumxy + x(i) * y(i) * n(i)
        sumxq = sumxq + x(i)**2 * n(i)
        sumyq = sumyq + y(i)**2 * n(i)
        sumn = sumn + n(i)
      enddo ! i = 1,np ! was do 1
!
! Diff w/ old code in sumx (because of very small diffs in x on entry)
!     write(6,"('onepar: sumx,y,xy       =',3e12.4)") sumx,sumy,sumxy
!     write(6,"('onepar: sumxq,sumyq,sumn=',3e12.4)") sumxq,sumyq,sumn

      sumxy = sumxy * sumn
      sumxq = sumxq * sumn
      sumxty = sumx * sumy
      sumyq = sumyq * sumn
      sumx2 = sumx**2
      sumy2 = sumy**2
      xl = sumxq - sumx2
      yl = sumyq - sumy2
      if (abs(xl) < 1.0e-25 .or. abs(yl) < 1.0e-25) return
      sl = (sumxy - sumxty) / xl
      sb = sqrt (abs (sumn / xl))
      sl1 = (sumxy - sumxty) / yl
      sb1 = sqrt( abs(sumn/yl))
      xl = xl * sumn
      yl = yl * sumn
!     syx2 = sqrt (abs (sumx2 / xl))
      cr = (sumxq * sumy - sumx * sumxy) / xl
      crb = sumxq/xl
      cr1 = (sumyq*sumx - sumy*sumxy)/ yl
      crb1 = sumyq/yl
      crb1 = sqrt(abs(crb1))
      crb = sqrt( crb)
      r2 = sqrt (abs (sl * sl1))
      end subroutine onepar
C

