c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timegcm/mksrc/getforbes.f
c------------------------------------------------------------------
c
      subroutine getforbe(ieqnx,zpbot,gcmlat,jmx,tng,tnlb,htlb,unlb,
     +  iplt,iframe)
c
      parameter(nflat=19)
      dimension tng(jmx),tnlb(jmx),htlb(jmx),unlb(jmx),gcmlat(jmx)
      dimension rit(nflat),riu(nflat),riz(nflat),ril(nflat)
      common/atm/zt(nflat*2-1),zu(nflat*2-1),y(nflat*2-1)
      character*40 toplab
      dimension viewport(4)
      data viewport /.11,.89,.21,.86/
c
c Either equinox or solstice:
c
      month = 12
      if (ieqnx.gt.0) month = 3
c
c Convert zpbot (desired lower boundary) from -ln(p/p0) to mb:
c
      pbot = exp(-zpbot) * 5.e-7   ! mb press at zp -13 (.221)
c
c Call Forbes (it uses msis90):
c
      write(6,"(' ')")
      write(6,"('getforbe: month=',i3,' zpbot=',f5.1,
     +  ' pbot(mb)=',f10.4)") month,zpbot,pbot
      call lbtgcm(month,pbot)
      kt=0
      do k=1,nflat*2-1,2
	kt=kt+1
	ril(kt)=90.-float(k-1)*5.  ! latitude
	riz(kt)=y(k)/1000.         ! height
	rit(kt)=zt(k)              ! zonal mean tn
	riu(kt)=zu(k)              ! zonal mean winds
      enddo
c
c Interpolate to gcm latitude grid:
c
      do j=1,jmx
        call bracket(gcmlat(j),ril,nflat,0,j1,j2)
        f1 = (gcmlat(j)-ril(j1)) / (ril(j2)-ril(j1))
        tnlb(j) = f1*rit(j2) + (1.-f1)*rit(j1)
        htlb(j) = f1*riz(j2) + (1.-f1)*riz(j1)
        unlb(j) = f1*riu(j2) + (1.-f1)*riu(j1)
c       write(6,"(' j=',i3,' gcmlat=',f6.2,' j1 j2=',2i3,/
c    +    ' rit(j1,2)=',2f10.3,' tnlb(j)=',f10.3,/
c    +    ' riz(j1,2)=',2f10.3,' htlb(j)=',f10.3,/
c    +    ' riu(j1,2)=',2f10.3,' unlb(j)=',f10.3)")
c    +    j,gcmlat(j),j1,j2,rit(j1),rit(j2),tnlb(j),
c    +    riz(j1),riz(j2),htlb(j),riu(j1),riu(j2),unlb(j)
      enddo
c
c Un must be zero at poles:
c
      unlb(1) = 0.
      unlb(jmx) = 0.
c
c Make line plots:
c height:
c
      if (iplt.gt.0) then
      call aggeti('SET.',iset)
      call agseti('SET.',1)
c
      if (ieqnx.gt.0) then
        write(toplab,"('ZATMOS HEIGHT (EQUINOX)',17x)")
      else
        write(toplab,"('ZATMOS HEIGHT (SOLSTICE)',16x)")
      endif
      call displa(1,1,1)
      call anotat('LATITUDE (GCM GRID)','HEIGHT (KM)',0,1,-1,L)
      call ezxy(gcmlat,htlb,jmx,toplab)
      iframe = iframe+1
      write(6,"('getforbes frame ',i3,': made line plot of htlb')")
     +  iframe
c
c tn:
c
      if (ieqnx.gt.0) then
        write(toplab,"('ZATMOS TEMPERATURE (EQUINOX)',12x)")
      else
        write(toplab,"('ZATMOS TEMPERATURE (SOLSTICE)',11x)")
      endif
      call displa(1,1,1)
      call anotat('LATITUDE (GCM GRID)','TN (K)',0,1,-1,L)
      call ezxy(gcmlat,tnlb,jmx,toplab)
      iframe = iframe+1
      write(6,"('getforbes frame ',i3,': made line plot of tnlb')")
     +  iframe
c
c un:
c
      if (ieqnx.gt.0) then
        write(toplab,"('ZATMOS ZONAL WIND (EQUINOX)',13x)")
      else
        write(toplab,"('ZATMOS ZONAL WIND (SOLSTICE)',12x)")
      endif
      call displa(1,1,1)
      call anotat('LATITUDE (GCM GRID)','UN (M/S)',0,1,-1,L)
      call ezxy(gcmlat,unlb,jmx,toplab)
      iframe = iframe+1
      write(6,"('getforbes frame ',i3,': made line plot of unlb')")
     +  iframe
c
      call agseti('SET.',iset)
      endif ! iplt
c
      return
      end