
      program test
c
c     parameter (cpspval=1.e36,nlat=36,nzp=37)
      parameter (cpspval=1.e36,nlat=36,nzp=19)
      dimension vb2(nlat,nzp),gcmlat(nlat),fub2(nlat,nzp),
     +  gcmzp(nzp),d(8),t(2)
      dimension iasf(13),viewport(4)
      character*56 toplab
      data toplab
     + /'                                                        '/
      data iasf/13*1/
      data viewport /.11,.89,.21,.86/

      call opngks
      call gsclip(0)
      call gsasf(iasf)
      call gsfais(1)
      call cpsetup(cpspval)
c
c Timegcm zp and lat grid:
c
      do k=1,nzp
        gcmzp(k) = -13. + (k-1)*0.5
      enddo
      do j=1,nlat
        gcmlat(j) = -87.5+(j-1)*5.
      enddo
      iyd = 76355
      sec = 0.
      glong = 0.
      slt = 0.
      f107a = 72.
      f107 = 67.
      ap = 4.
c
c  Input:	z   =  altitude in km
c	     xlat   =  latitude in degrees
c  Output:
c              vb   =  prevailing northward wind (m/s)
c             fub   =  eastward acceleration (m/s2)
c
      do k=1,nzp
c
c Find altitudes corresponding to zp (using msis90):
c (December solstice case only) and call wvparam at those alts:
c
        zpmb = exp(-gcmzp(k)) * 5.e-7
        write(6,"('k=',i2,' zp=',f5.1,' zpmb=',e12.4)") k,gcmzp(k),zpmb
        do j=1,nlat
          call ghp6(iyd,sec,alt,gcmlat(j),glong,slt,f107a,f107,ap,d,t,
     +      zpmb)
c         write(6,"('gcmlat=',f8.2,' gcmzp=',f5.1,' zpmb=',e12.4,
c    +      ' alt=',f10.2,' rho=',e12.4)") 
c    +      gcmlat(j),gcmzp(k),zpmb,alt,d(6)
          call wvparam(alt,gcmlat(j),vb,fub)
          vb2(j,k) = vb
          fub2(j,k) = fub
        enddo
      enddo
c
c     subroutine contour(zz,idimx,nx,ny)
c
      call set(viewport(1),viewport(2),viewport(3),viewport(4),
     +         -90.,90.,gcmzp(1),gcmzp(nzp),1)
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',gcmlat(1))
      call cpsetr('XCM',gcmlat(nlat))
      call cpsetr('YC1',gcmzp(1))
      call cpsetr('YCN',gcmzp(nzp))
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.20)
c
c Contour northward wind:
c
      call contour(vb2,nlat,nlat,nzp)
      write(toplab,"('WAVEDRAG NORTHWARD WIND (M/S)  ',
     +  'DECEMBER SOLSTICE',8x)") 
      call labrect(gcmlat,nlat,gcmzp,nzp,'LATITUDE (DEG)',
     +  'ZP',0.)
      call wrlab(toplab(1:lenstr(toplab)),xmid,0.9,0.)
      call frame
      write(6,"('wavedrag: completed frame: ',a)") 
     +  toplab(1:lenstr(toplab))
c
c Contour eastward acceleration:
c
      call contour(fub2,nlat,nlat,nzp)
      write(toplab,"('EASTWARD ACCELERATION (M/S2)  ',
     +  'DECEMBER SOLSTICE',9x)") 
      call labrect(gcmlat,nlat,gcmzp,nzp,'LATITUDE (DEG)',
     +  'ZP',0.)
      call wrlab(toplab(1:lenstr(toplab)),xmid,0.9,0.)
      call frame
      write(6,"('wavedrag: completed frame: ',a)") 
     +  toplab(1:lenstr(toplab))

      call clsgks
c
      stop
      end
