c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcm/dynsatv.f
c------------------------------------------------------------------
c
      subroutine dynsatv(utc,id)
c
c Make sat view projection contours:
c
      include 'tgcmparam.h'
      include 'input.h'
      include 'tigcm.h'
      include 'tgcmhdr.h'
      include 'tiegcmfld.h'
      include 'tigcmfld.h'
      include 'selgrid.h'
      include 'tgcmlab.h'
      include 'cool.h'
      include 'plt.h'
      include 'color.h'
      dimension plotp(imx,jmx)
      character*56 toplab,lab,labgcm1
      character*24 latlab
      character*18 radlab
      save rmintot,rmaxtot
      data rmintot,rmaxtot/1.e36,-1.e36/
c
c Below are essentially dummies for bwcon and clrcon calls:
      parameter (nx=7,ny=7)
      dimension numx(nx),numy(ny)
      data numx/-180,-120,-60,0,60,120,180/
      data numy/-90,-60,-30,0,30,60,90/
      character*8 fmtx,fmty
      character*16 labx,laby
      data fmtx/'(i4)    '/, fmty/'(i3)    '/
      data nfmtx/4/, nfmty/3/, mnrx/3/, mnry/3/
      data laby/' LATITUDE (DEG) '/
      data labx/' LONGITUDE (DEG)'/ 
c
c Default distance from center of earth in earth radii:
      data sa /5.0/, dgrid/15./
c
c Temp:
      data maps/1/, iccont/1/
c
c Set up EZMAP:
c
      proj = 'SV      '  
      iel = 1
      pmapl = .05
      pmapr = .80
      pmapb = .13
      pmapt = .88
      xmid = 0.5*(pmapl+pmapr)
      xc = -180.
      xd = +180.
      yc = -90.
      yd = +90.
      call mappos(pmapl,pmapr,pmapb,pmapt)
      call mapsti('EL',iel)
      call mapstc('OU','CO')
      call mapsti('DO',1)
      call mapstr('SA',sa)
      call mapset('MA',0.,0.,0.,0.)
      call mapsti('G1',1)
      call mapsti('G2',2)
      call mapsti('VS',0)
      call mapstr('GR',dgrid)
      call mapint
c
c Set up conpack:
c (parameters common to all frames by this routine)
c
      rilx = 0.5
      rily = pmapb-0.20
      call cpsetr('ILX - info label x-coord',rilx)
      call cpsetr('ILY - info label y-coord',rily)
      call cpsetr('ORV - out of range value',1.e12)
      call cpseti('MAP - mapping flag',1)
      call cpseti('SET',0)
c
c Loop through desired lat/lon points for center of projection:
c
      do 100 l=1,nloc
c
c input longitude of r12flag means we want longitude at local noon: 
c
        projlat = rloc(1,l)
        projlon = rloc(2,l)
        if (rloc(2,l).eq.r12flag) then
          ixloc(2,l) = ixslt(12.,utc,projlon)
        endif
c
c Tell ezmap lat/lon and projection:
c
        call maproj(proj,projlat,projlon,0.)
        call mapint
        call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg)
c
c Field loop:
        write(6,"(' ')")
        write(6,"('Satellite View projections (dynamo fields):')")
	do 200 ip=1,ntiegcmf
          if (iptiegcm(ip).le.0) goto 200
c
c Selected pressure and heights loop:
c
          do 110 k=1,npls+nhts
            kk = k-npls
c
c At constant pressure surface:
            if (k.le.npls) then
              write(toplab,"('UT = ',f6.2,' ZP = ',f4.1,' DAY = ',
     +          i5,23x)") utc,spls(k),iyd(id)
c
c Define plot array (interpolate to height surface):
              if (ip.le.ntiegcmf-5) then
                plotp(:,:) = dynpnt(:,ispls(k),:,ip) 
              elseif (ip.eq.ixspot) then ! sum potentials
                plotp(:,:) = dynpnt(:,ispls(k),:,ixpot) + hpoten(:,:)
              elseif (ip.gt.ixspot) then
                if (k.gt.1) goto 110
                write(toplab,"('UT = ',f6.2,45x)") utc
              endif
              if (ip.eq.ixhpot) then
                plotp(:,:) = hpoten(:,:)
              elseif (ip.eq.ixhui) then
                plotp(:,:) = hui(:,:)
              elseif (ip.eq.ixhvi) then
                plotp(:,:) = hvi(:,:)
              elseif (ip.eq.ixhwi) then
                plotp(:,:) = hwi(:,:)
              endif
c
c Interpolate to constant height surface:
            else
              write(toplab,"('UT = ',f6.2,' HT = ',f6.1,' DAY = ',
     +          i5,21x)") utc,shts(kk),iyd(id)
              if (ip.le.ntiegcmf-5) then
                call tgcmint(dynpnt(1,1,1,ip),pnt(1,1,1,ixz),gcmlat,jmx,
     +            gcmlon,imx,shts(kk),1,idynlog(ip),2,plotp,imx,jmx,
     +            1,kmx,gcmzp,ier,1,cpspval,0)
              elseif (ip.eq.ixspot) then ! sum potentials
                call tgcmint(dynpnt(1,1,1,ixpot),pnt(1,1,1,ixz),gcmlat,
     +            jmx,gcmlon,imx,shts(kk),1,idynlog(ip),2,plotp,imx,jmx,
     +            1,kmx,gcmzp,ier,1,cpspval,0)
                plotp(:,:) = plotp(:,:) + hpoten(:,:)
              else ! height independent heelis potential
                write(6,"('pltsat skipping heelis because',
     +            ' it is height independent: ip=',i3)") ip
                goto 110
              endif
            endif
c
c Take logs if necessary:
          if (idynlog(ip).gt.0.and.log10pl.gt.0) then
            do 125 j=1,jmx
            do 125 i=1,imx
              if (plotp(i,j).eq.cpspval.or.plotp(i,j).le.1.e-30) 
     +                                                          then
                plotp(i,j) = cpspval
              else
                plotp(i,j) = alog10(plotp(i,j))
              endif
 125        continue
            write(lab,"('LOG10 ',a50)") labtiegcm(ip)(1:50)
          else
            lab = labtiegcm(ip)
          endif
c
c Find min/max at this projection:
c (Used to fix contour levels for movies)
          rmin = 1.e36
          rmax = -1.e36
          do 130 j=1,jmx
          do 130 i=1,imx
            call maptrn(gcmlat(j),gcmlon(i),xpos,ypos)
            if (xpos.eq.1.e12.or.ypos.eq.1.e12) goto 130
            if (plotp(i,j).ne.cpspval) then
              if (plotp(i,j).lt.rmin) rmin = plotp(i,j)
              if (plotp(i,j).gt.rmax) rmax = plotp(i,j)
            endif
 130      continue
          rmintot = amin1(rmintot,rmin)
          rmaxtot = amax1(rmaxtot,rmax)
c
c Contour it, with axes and axis labels:
c
          if (iclrfill.le.0) then
            call bwcon(plotp,imx,imx,jmx,xc,xd,yc,yd,
     +                 cintdyn(ip),cmindyn(ip),cmaxdyn(ip),
     +                 0,nx,numx,mnrx,labx,fmtx,nfmtx,
     +                 ny,numy,mnry,laby,fmty,nfmty)
          else
            lbar = 1
            call clrcon(plotp,imx,imx,jmx,xc,xd,yc,yd,
     +                  cintdyn(ip),cmindyn(ip),cmaxdyn(ip),
     +                  0,nx,numx,mnrx,labx,fmtx,nfmtx,
     +                  ny,numy,mnry,laby,fmty,nfmty,
     +                  lbar,0.83,0.98,pmapb,pmapt,ip)
          endif
c         call box(1)
c
c Change color indices for continental outlines and lat/lon grid:
c
          if (maps.gt.0) then
            call gqplci(ierr,isplci)
            call gqpmci(ierr,ispmci)
            if (iclrfill.le.0) then
              call gsplci(iwhite)
              call gspmci(iwhite)
            else
              call gsplci(iccont)
              call gspmci(iccont)
            endif
c
c Add continental outlines and lat/lon grid and labels
c
            call mapsti('PE',1)
            call mapsti('LA',1)
            call maplot
            call mapgrd
            call maplbl
            call gsplci(isplci)
            call gspmci(ispmci)
          endif
c
c Mark location (center of projection):
c
          szfs = .013
          call maptrn(projlat,projlon,poslat,poslon)
          call plchhq(poslat,poslon,':KRL:G',szfs,0.,0.)
c
c Label location in upper right:
c
          write(latlab,"(' at (',f5.1,',',f6.1,')',6x)")
     +      projlat,projlon
          call getset(vl,vr,vb,vt,wl,wr,wb,wt,lnlg)
          call set(0.,1.,0.,1.,0.,1.,0.,1.,1)
          ypos = vt + 0.135*(vt-vb)
          call plchhq(0.70,ypos,':KRL:G',szfs,0.,0.)
          call plchhq(0.96,ypos,latlab(1:18),szfs,0.,1.)
c
c Label earth radii:
c
          write(radlab,"('(',f4.1,' earth radii)')") sa
          ypos = ypos-0.045*(vt-vb)
          call plchhq(0.96,ypos,radlab,szfs,0.,1.)
          call set(vl,vr,vb,vt,wl,wr,wb,wt,lnlg)
c
c Mark longitude at local noon
c
          rlon12 = 15.*(12.-utc)
          if (rlon12.ge.180.) rlon12 = rlon12-360.
          if (rlon12.lt.-180.) rlon12 = rlon12+360.
c
c For black and white, use stars near each pole, and at equator:
c
          if (iclrfill.le.0) then
            rlat12 = 90.-1.3*dgrid
            call maptrn(rlat12,rlon12,xnoon,ynoon)
            szfs = .025
            if (xnoon.ne.1.e12.and.ynoon.ne.1.e12)
     +        call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.)
            rlat12 = -90.+1.3*dgrid
            call maptrn(rlat12,rlon12,xnoon,ynoon)
            if (xnoon.ne.1.e12.and.ynoon.ne.1.e12)
     +          call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.)
            call maptrn(0.,rlon12,xnoon,ynoon)
            if (xnoon.ne.1.e12.and.ynoon.ne.1.e12)
     +        call plchhq(xnoon,ynoon,':KRL:-',szfs,0.,0.)
          else
c
c For color, draw wide red line along local noon meridian:
c (default line width (LW) is 1000)
c
            call setusv('LW',1900)
            call gsplci(ired)
            call gspmci(ired)
            call setusv('LW',1000)
            call gsplci(iwhite)
            call gspmci(iwhite)
          endif
c
c Add top labels:
c
          if (idynlog(ip).gt.0.and.log10pl.gt.0) then
            write(labgcm1,"('LOG10 ',a50)") labtiegcm(ip)(1:50)
            call top2lab(labgcm1,0.135,toplab,0.09,.017)
          else
            call top2lab(labtiegcm(ip),0.135,toplab,0.09,.017)
          endif
          write(lab,"(12x,a,12x)") runlab
          call botlab(lab,xmid,.05)
c
c Bottom history label:
          if (iclrfill.le.0) then
            write(lab,"(8x,'TIGCM History = ',3a8,8x)")
     +        (histvol(i,ivf),i=1,3)
            call botlab(lab,0.45,.02)
          endif
c
c Add tiegcm ion drift vectors if doing tiegcm potential:
          if (ipltuv.eq.1.and.
     +       (ip.eq.ixpot.or.ip.eq.ixhpot.or.ip.eq.ixspot)) then
            if (ip.eq.ixpot) then
              if (k.le.npls) then
                udat(:,:) = ui(:,ispls(k),:)
                vdat(:,:) = vi(:,ispls(k),:)
              else
c
c Do height interpolation of ui,vi:
c
                call tgcmint(ui,pnt(1,1,1,ixz),gcmlat,jmx,
     +            gcmlon,imx,shts(kk),1,idynlog(ixui),2,udat,imx,jmx,
     +            1,kmx,gcmzp,ier,1,cpspval,0)
                call tgcmint(vi,pnt(1,1,1,ixz),gcmlat,jmx,
     +            gcmlon,imx,shts(kk),1,idynlog(ixvi),2,vdat,imx,jmx,
     +            1,kmx,gcmzp,ier,1,cpspval,0)
              endif
c
c Add heelis ion drift vectors if doing heelis potential
c (height independent):
            elseif (ip.eq.ixhpot) then
              if (k.gt.1) goto 96
              udat(:,:) = hui(:,:)
              vdat(:,:) = hvi(:,:)
c
c Add ion drifts calculated from sum of dynamo and heelis potentials:
c (Cannot yet do this at interpolated heights)
c
            elseif (ip.eq.ixspot) then
              if (k.gt.npls) then
                write(6,"('dynsatv: ion drifts for sum of ',
     +            'potentials at interpolated heights not yet',
     +            ' available')")
                goto 96
              else
                call mkuivis(k,udat,vdat,pnt(1,1,1,ixz),lumag)
              endif
            endif
c
c Check vectors for zeros and plot:
            do 95 j=1,jmx
            do 95 i=1,imx
              if (udat(i,j).eq.0..and.vdat(i,j).eq.0.) then
                udat(i,j) = cpspval
                vdat(i,j) = cpspval
              endif
 95         continue
            uvmax = 0.
            m = 1
            ispvec = 0
            spvec = 1.e36
            iskip = 0
            call pltvec(udat,vdat,uvmax,jmx,iskip,gcmlat(1),
     +        gcmlat(jmx),m,ispvec,spvec)
          endif
 96       continue
c
c Finish up:
c         call box(1)
          call frame
          iframe=iframe+1
          if (k.eq.1.and.(ip.eq.ixfof2.or.ip.eq.ixhmf2)) then
            write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1,
     +        ' lat,lon=',2f8.2)")
     +        iframe,labtiegcm_short(ip),utc,projlat,projlon
          else
            if (k.le.npls) then
              write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1,
     +          ' ZP=',f4.1,' lat,lon=',2f8.2)")
     +          iframe,labtiegcm_short(ip),utc,spls(k),projlat,projlon
            else
              write(6,"('dynsatv frame ',i3,': ',a,' UT=',f4.1,
     +          ' HT=',f6.1,' lat,lon=',2f8.2)")
     +          iframe,labtiegcm_short(ip),utc,shts(kk),projlat,
     +          projlon
            endif
          endif
c
c End pressure/height loop:
 110      continue
c
c End field loop:
 200    continue
c
c End loc loop:
 100  continue
c
      return
      end