c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timegcm/pltstrm.f
c------------------------------------------------------------------
c
      subroutine pltstrm(fields,utc,ihtsc)
c
c Calculate stream function from VN and zonal means, and contour:
c (adapted from pltlon)
c
      include 'timesproc.h'
      character*56 toplab,fieldlab,saveilt
      character*80 rec80
      dimension plt(jmx,kmx),viewport(4),yaxht(kmx),rimx(imx),
     +  plt0(jmx,kmx),fields(imx,kmx,jmx,nfget),hts(jmx,kmx),
     +  clev(13)
      pointer(ppltht,pltht(1))
      data viewport /.15,.89,.26,.91/
      data clev/-1.E8,-1.E7,-1.E6,-1.E5,-1.E4,-1.E3,0.,1.E3,1.E4,1.E5,
     +         1.E6,1.E7,1.E8/, ncl/13/
      data cint/10000./
      external drawcl,colram
c
c work arrays from ~/libplt/cpwrk.h:
      parameter(lrwrk=2000,liwrk=1000,liama=300000,lcra=10000)          
      common/work/ rwrk(lrwrk),iwrk(liwrk),iama(liama),iiasf(13),
     +  iaia(10),igia(10),xcra(lcra),ycra(lcra)                         
c
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',gcmlat(1))
      call cpsetr('XCM',gcmlat(jmx))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.16)
      call cpseti('ILP',0)
      call cpsetr('ILS',.016)
      if (labels.le.1) call cpsetr('ILS',.02)
      if (ihtsc.le.0) then
        call cpsetr('YC1',gcmzp(izprange(1)))
        call cpsetr('YCN',gcmzp(izprange(2)))
        call set(viewport(1),viewport(2),viewport(3),viewport(4),
     +           -90.,90.,gcmzp(izprange(1)),gcmzp(izprange(2)),1)
      else
        call cpsetr('YC1',htscale(1))
        call cpsetr('YCN',htscale(nhtscale))
        call set(viewport(1),viewport(2),viewport(3),viewport(4),
     +           -90.,90.,htscale(1),htscale(nhtscale),1)
        call hpalloc(ppltht,imx*nhtscale,ier,1)
        if (ier.ne.0) then
          write(6,"('pltstrm: hpalloc error = ',i6,' for pltht')") ier
          stop 'ppltht'
        endif
      endif
      nzp = izprange(2)-izprange(1)+1
      ixlon = ifix(zmflag)
c
c The call getlon stores vn in plt, then calcstrm returns stream func
c in plt0 (vn still in plt). Then if zp on y-axis, transfer desired
c part of stream to plt, otherwise interpolate stream function to
c height scale, storing in pltht. After these calls desired stream
c func is in plt or pltht:
c
      call getlon(fields,plt,ixlon,0,itxv)
      call calcstrm(plt,jmx,kmx,gcmzp,gcmlat,plt0)
      if (ihtsc.le.0) then
        do k=izprange(1),izprange(2)
          kk = k-izprange(1)+1
          plt(:,kk) = plt0(:,k)
        enddo
        ny = izprange(2)-izprange(1)+1
      else
        call getlon(fields,hts,ixlon,0,itxz)
        call cuthtint(plt0,hts,jmx,kmx,pltht,htscale,nhtscale,
     +    0,cpspval,ier,1)
        ny = nhtscale 
      endif
c
c Contour (hardwire contour levels to get log effect):
c (cpcnrc not used for labels=1 because need to hardwire contour levels)
c
      call cpseti('CLS',0)
      call cpseti('NCL',ncl)
      do i=1,ncl
        call cpseti('PAI',i)
        call cpsetr('CLV',clev(i))
        if (clev(i).ge.0.) then
          call cpseti('CLD',65535)
        else
          call cpsetc('CLD -- contour line dash pattern',
     +      '$$''''$$''''')
        endif
c       if (mod(i,3).eq.0) then
c         call cpseti('CLU',3)
c       else
c         call cpseti('CLU',1)
c       endif
        if (icolor.gt.0) call cpseti('CLC',0)
      enddo
      call cpgetc('ILT',saveilt)
      call cpsetc('ILT','MIN $ZMN$, MAX $ZMX$')
      if (ihtsc.le.0) then
        call cprect(plt,jmx,jmx,ny,rwrk,lrwrk,iwrk,liwrk)
        if (icolor.le.0) then
          call cpcldr(plt,rwrk,iwrk)
          call cplbdr(plt,rwrk,iwrk)
        else
          call arinam(iama,liama)
          call cpclam(plt,rwrk,iwrk,iama)
          call arscam(iama,xcra,ycra,lcra,iaia,igia,10,colram)
          call cplbam(plt,rwrk,iwrk,iama)
          call cplbdr(plt,rwrk,iwrk)
          call cpcldm(plt,rwrk,iwrk,iama,drawcl)
        endif
        call labrect(gcmlat,jmx,gcmzp(izprange(1)),ny,'LATITUDE',
     +    'ZP',0.)
c
c Height on y-axis:
      else
        call cprect(pltht,jmx,jmx,ny,rwrk,lrwrk,iwrk,liwrk)
        if (icolor.le.0) then
          call cpcldr(pltht,rwrk,iwrk)
          call cplbdr(pltht,rwrk,iwrk)
        else
          call arinam(iama,liama)
          call cpclam(pltht,rwrk,iwrk,iama)
          call arscam(iama,xcra,ycra,lcra,iaia,igia,10,colram)
          call cplbam(pltht,rwrk,iwrk,iama)
          call cplbdr(pltht,rwrk,iwrk)
          call cpcldm(pltht,rwrk,iwrk,iama,drawcl)
        endif
        call labrect(gcmlat,jmx,htscale,ny,'LATITUDE',
     +    'HEIGHT (KM)',0.)
      endif
      call cpsetc('ILT',saveilt)
c
c Top label with field name, lon/slt, ut, etc:
c
      call clearstr(toplab)
      write(toplab,"('STREAM FUNCTION   UT=',f6.2)") utc
      call clearstr(rec80)
      write(rec80,"(a)") toplab(1:lenstr(toplab))
      call clearstr(fieldlab)
      write(fieldlab,"('STREAM FUNCTION')")
      if (labels.gt.1) then
        call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.05,
     +    0.)
        write(toplab,"('HISTORY=',a)") histvol(ivol)
        call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.15,
     +    .012)
      else
        call wrlablq(toplab(1:lenstr(toplab)),xmid,
     +    viewport(4)+0.05,.02)
      endif
c
c Put height axis on right (use latitudinally averaged heights):
c
      if (ihtsc.gt.0) goto 300
      yaxht(:) = 0.
      do k=izprange(1),izprange(2)
        izp = k-izprange(1)+1
        do j=1,jmx
          rimx(:) = fields(:,k,j,ifget(itxz))
          yaxht(izp) = yaxht(izp)+calcmean(rimx,imx,0,1.e-20,cpspval)
        enddo
        yaxht(izp) = yaxht(izp) / float(jmx)
      enddo 
      rnd = 10.
      if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5.
      call altyax(nzp,yaxht,gcmzp(izprange(1)),rnd,6)
 300  continue
c
c Wrap it up, including ascii data file, if needed:
c
      call frame
      iframe = iframe+1
c
c Pressure on y-axis:
      if (ihtsc.le.0) then
        if (iwrascii.gt.0) then
          call wrascii(iwrascii,luascii,plt,jmx,ny,gcmlat,
     +      gcmzp(izprange(1)),'LATITUDE','LN(P0/P)',
     +      histvol(ivol),fieldlab,iframe,rec80,'timesproc',
     +      dirascii)
        endif
        write(6,"('pltstrm: frame ',i4,' STREAM: zprange=',2f8.2)")
     +    iframe,zprange 
c
c Height on y-axis:
      else
        if (iwrascii.gt.0) then
          call wrascii(iwrascii,luascii,pltht,jmx,ny,gcmlat,
     +      htscale,'LATITUDE','HEIGHT (KM)',histvol(ivol),fieldlab,
     +        iframe,rec80,'timesproc',dirascii)
        endif
        write(6,"('pltstrm: frame ',i4,' STREAM: htscale=',f7.2,
     +    ' to ',f7.2)") iframe,htscale(1),htscale(nhtscale) 
      endif
      if (ihtsc.gt.0) call hpdeallc(ppltht,ier,0)
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine calcstrm(vn,jmx,kmx,gcmzp,gcmlat,strm)
      dimension vn(jmx,kmx),gcmzp(kmx),strm(jmx,kmx),expz(kmx),
     +  gcmlat(jmx),cs(jmx)
c
c Calculate stream function from zonal means of vn
c On input:
c   vn(jmx,kmx) = zonal means of vn
c   jmx,kmx = dimensions (36,45 for times)
c   gcmzp(kmx) = grid pressure levels
c   gcmlat(jmx) = grid geographic latitudes (deg)
c On output:
c   stream function strm(jmx,kmx) is defined
c
      dz = 0.5
      pi = 3.1415926535
      dtr = pi/180.
      rad = 6.37e+8
      p0 = 5.e-4
      grav = 870.
      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)*vn(j,kmx)*cs(j)
        do kk=1,kmx-1
          k = kmx-kk
          strm(j,k) = strm(j,k+1)+0.5*dz*alfa*cs(j)*
     +      (expz(k+1)*vn(j,k+1)+expz(k)*vn(j,k))
        enddo
      enddo
      return
      end