c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcmdif/pltloc.f
c------------------------------------------------------------------
c
      subroutine pltloc(utf,hts,ihtsc,l)
c
c Contour ut on x-axis, zp on y-axis if ihtsc <= 0; 
c   if ihtsc > 0, then ht on y-axis
c Location index = l
c
      include 'tigcmloc.h'
      character*56 toplab,fieldlab
      character*80 rec80
      dimension viewport(4),yaxht(kmx),rimx(imx),plt(ntms,kmx),
     +  plt0(ntms,kmx),utf(ntms,kmx,nfplt,nloc),hts(ntms,kmx,nloc),
     +  pltht(ntms,nhtscale)
      data viewport /.12,.86,.26,.91/
c
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',utinc(1))
      call cpsetr('XCM',utinc(ntms))
      call cpsetr('ILX',xmid)
      rily = -.17
      if (gloc(1,l).ne.gmflag.and.gloc(2,l).ne.gmflag) rily = -.34
      call cpsetr('ILY',rily)
      call cpseti('ILP',0)
      call cpsetr('ILS',.016)
      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),
     +           utinc(1),utinc(ntms),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),
     +           utinc(1),utinc(ntms),htscale(1),htscale(nhtscale),1)
      endif
c
c Field loop: 
c
      do 200 ip=1,nftot
        if (ifplt(ip).le.0) goto 200
        if (ip.eq.ixfof2.or.ip.eq.ixhmf2) goto 200
c
c Get field and interpolate to height scale if necessary:
c (Also take log10 if needed)
c
        if (ihtsc.le.0) then 
          do k=izprange(1),izprange(2)
            kk = k-izprange(1)+1
            plt(:,kk) = utf(:,k,ifplt(ip),l)
          enddo
          ny = izprange(2)-izprange(1)+1
          if (logloc.gt.0.and.logterp(ip).gt.0) 
     +      call log10f(plt,ntms*ny,1.e-20,cpspval)
        else
          plt0(:,:) = utf(:,:,ifplt(ip),l)
          call cuthtint(plt0,hts(1,1,l),ntms,kmx,pltht,htscale,
     +      nhtscale,logterp(ip),cpspval,ier,1)
          ny = nhtscale
          if (logloc.gt.0.and.logterp(ip).gt.0) 
     +      call log10f(pltht,ntms*ny,1.e-20,cpspval)
        endif
c
c Contour:
c
        if (icolor.le.0) then
          if (ihtsc.le.0) then
            call contour(plt,ntms,ntms,ny,cint(ip),cmin(ip),cmax(ip))
          else
            call contour(pltht,ntms,ntms,ny,cint(ip),cmin(ip),
     +        cmax(ip))
          endif
        else
          if (ihtsc.le.0) then
            call conclr(plt,ntms,ntms,ny,cint(ip),cmin(ip),cmax(ip))
          else
            call conclr(pltht,ntms,ntms,ny,cint(ip),cmin(ip),cmax(ip))
          endif
        endif
c
c Add axes labels:
c
        isltax = 0
        if (gloc(1,l).ne.gmflag.and.gloc(2,l).ne.gmflag) isltax = 1
        if (ihtsc.le.0) then 
          call labutxy(mtimes,ntms,gcmzp(izprange(1)),ny,'ZP',
     +      0.,isltax,gloc(2,l))
        else
          call labutxy(mtimes,ntms,htscale,nhtscale,'HEIGHT (KM)',
     +      0.,isltax,gloc(2,l))
        endif
c
c Field label at top:
c
        call clearstr(toplab)
        if (logloc.gt.0.and.logterp(ip).gt.0) then
          len = lenstr(flab(ip))
          if (len+6.le.56) then
            write(toplab,"('LOG10 ',a)") flab(ip)(1:len)
          endif
        else
          toplab = flab(ip)
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+.07,0.)
        fieldlab = toplab
c
c Location label (may be at selected lat,lon, or at global means, or
c zonal means at selected latitude):
c
        call clearstr(toplab)
        if (gloc(1,l).eq.gmflag.and.gloc(2,l).eq.gmflag) then
            write(toplab,"('GLOBAL MEANS')")
            xoff = .15
        elseif (gloc(1,l).ne.gmflag.and.gloc(2,l).eq.gmflag) then
          write(toplab,"('LAT= ',f6.2,' LON=ZONAL MEANS')") gloc(1,l)
          xoff = .15
        else
          if (locname(1).eq.'xxxxxxxxxxxxxxxx') then
            write(toplab,"('LAT,LON = ',f6.2,f7.2)") gloc(1,l),
     +        gloc(2,l)
            xoff = .25
          else
            write(toplab,"('LAT,LON = ',f6.2,f7.2,' (',a,')')") 
     +        gloc(1,l),gloc(2,l),locname(l)(1:lenstr(locname(l)))
            xoff = .25
          endif
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(4)+0.03,0.)
        call clearstr(rec80)
        write(rec80,"(a)") toplab(1:lenstr(toplab))
c
c History volume label at bottom:
c
        call clearstr(toplab)
        if (nhvols.eq.1) then
          write(toplab,"('HISTORY=',a)") histvol(ivol)
        else
          write(toplab,"('FIRST HISTORY VOL=',a)") histvol(1)
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-xoff,
     +    .010)
c
c Put height axis on right if pressure on y-axis
c (global or zonal means are already in utf(...,l))
c
        if (ihtsc.le.0) then
          yaxht(:) = 0.
          do k=izprange(1),izprange(2)
            izp = k-izprange(1)+1
            do i=1,ntms
              yaxht(izp) = yaxht(izp) + hts(i,k,l)
            enddo
            yaxht(izp) = yaxht(izp) / float(ntms)
          enddo
          rnd = 10.
          if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5.
          call altyax(ny,yaxht,gcmzp(izprange(1)),rnd,6)
        endif
c
c Wrap it up: 
c (make ascii file if needed)
c
        call frame
        iframe = iframe+1
        if (ihtsc.le.0) then
          if (iwrascii.gt.0) then
            call wrascii(iwrascii,luascii,plt,ntms,ny,utinc,gcmzp,
     +        'UT (HRS)','LN(P0/P)',histvol(ivol),fieldlab,iframe,
     +        rec80,'tigcmloc',dirascii)
          endif
          if (gloc(1,l).eq.gmflag.and.gloc(2,l).eq.gmflag) then
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' (GLOBAL MEANS)  zprange=',2f8.2)") iframe,
     +        labshort(ip),zprange
          elseif (gloc(1,l).ne.gmflag.and.gloc(2,l).eq.gmflag) then
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' lat=',f6.1,' zonal means  zprange=',2f8.2)") 
     +        iframe,labshort(ip),gloc(1,l),zprange
          else
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' lat,lon=',f6.1,',',f7.1,' zprange=',2f8.2)") 
     +        iframe,labshort(ip),gloc(1,l),gloc(2,l),zprange
          endif
        else
          if (iwrascii.gt.0) then
            call wrascii(iwrascii,luascii,pltht,ntms,ny,utinc,htscale,
     +        'UT (HRS)','HEIGHT (KM)',histvol(ivol),fieldlab,iframe,
     +        rec80,'tigcmloc',dirascii)
          endif
          if (gloc(1,l).eq.gmflag.and.gloc(2,l).eq.gmflag) then
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' (GLOBAL MEANS)  htscale=',f6.1,' to ',f6.1)") iframe,
     +        labshort(ip),htscale(1),htscale(nhtscale)
          elseif (gloc(1,l).ne.gmflag.and.gloc(2,l).eq.gmflag) then
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' lat=',f6.1,' zonal means  htscale=',f6.1,' to ',
     +        f6.1)") iframe,labshort(ip),gloc(1,l),htscale(1),
     +        htscale(nhtscale)
          else
            write(6,"('loc frame ',i4,' field ',a8,
     +        ' lat,lon=',f6.1,',',f7.1,' htscale=',f6.1,' to ',
     +        f6.1)") iframe,labshort(ip),gloc(1,l),gloc(2,l),
     +        htscale(1),htscale(nhtscale)
          endif
        endif
c
c End field loop
 200  continue
c
c Plot volume emission rate:
      if (idop.gt.0) call pdop(hts,ihtsc,l)
c
c Plot ve,vu (used in ion vel calculation):
      if (ionvel.eq.2.and.(ifplt(ixui).gt.0.or.ifplt(ixvi).gt.0.or.
     +                     ifplt(ixwi).gt.0))
     +  call pvevu(hts,ihtsc,l)
      
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine pdop(hts,ihtsc,l)
c
c Contour volume emission rate (doppler):
c
      include 'tigcmloc.h'
      include 'doppler.h'
      character*56 toplab,fieldlab
      character*80 rec80
      dimension plt(ntms,kmx),plt0(ntms,kmx),vp(4)
      dimension hts(ntms,kmx,nloc),pltht(ntms,nhtscale),yaxht(kmx)
      data vp /.12,.86,.26,.91/
c
      ip = 5
      xmid = 0.5*(vp(1)+vp(2))
      if (ihtsc.le.0) then 
        do k=izprange(1),izprange(2)
          kk = k-izprange(1)+1
          do it=1,ntms
            plt(it,kk) = ver(it,k,l)
          enddo
        enddo
        ny = izprange(2)-izprange(1)+1
        if (logloc.gt.0.and.log10dop(ip).gt.0) 
     +    call log10f(plt,ntms*ny,1.e-20,cpspval)
      else
        do it=1,ntms
          plt0(it,:) = ver(it,:,l)
        enddo
        call cuthtint(plt0,hts(1,1,l),ntms,kmx,pltht,htscale,
     +    nhtscale,log10dop(ip),cpspval,ier,1)
        ny = nhtscale
        if (logloc.gt.0.and.log10dop(ip).gt.0) 
     +    call log10f(pltht,ntms*ny,1.e-20,cpspval)
      endif
c
c Contour:
c
      if (icolor.le.0) then
        if (ihtsc.le.0) then
          call contour(plt,ntms,ntms,ny,0.,1.,0.)
        else
          call contour(pltht,ntms,ntms,ny,0.,1.,0.)
        endif
      else
        if (ihtsc.le.0) then
          call conclr(plt,ntms,ntms,ny,0.,1.,0.)
        else
          call conclr(pltht,ntms,ntms,ny,0.,1.,0.)
        endif
      endif
c
c Add axes labels:
c
      isltax = 0
      if (gloc(1,l).ne.gmflag.and.gloc(2,l).ne.gmflag) isltax = 1
      if (ihtsc.le.0) then 
        call labutxy(mtimes,ntms,gcmzp(izprange(1)),ny,'ZP',
     +    0.,isltax,gloc(2,l))
      else
        call labutxy(mtimes,ntms,htscale,nhtscale,'HEIGHT (KM)',
     +    0.,isltax,gloc(2,l))
      endif
c
c Field label at top:
c
      call clearstr(toplab)
      if (logloc.gt.0.and.log10dop(ip).gt.0) then
        len = lenstr(dopname(ip))
        if (len+6.le.56) then
          write(toplab,"('LOG10 ',a)") dopname(ip)(1:len)
        endif
      else
        toplab = dopname(ip)
      endif
      call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.07,0.)
      fieldlab = toplab
c
c Location label (may be at selected lat,lon, or at global means, or
c zonal means at selected latitude):
c
      call clearstr(toplab)
      if (locname(1).eq.'xxxxxxxxxxxxxxxx') then
        write(toplab,"('LAT,LON = ',f6.2,f7.2)") gloc(1,l),
     +    gloc(2,l)
        xoff = .25
      else
        write(toplab,"('LAT,LON = ',f6.2,f7.2,' (',a,')')") 
     +    gloc(1,l),gloc(2,l),locname(l)(1:lenstr(locname(l)))
        xoff = .25
      endif
      call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+0.03,0.)
      call clearstr(rec80)
      write(rec80,"(a)") toplab(1:lenstr(toplab))
c
c History volume label at bottom:
c
      call clearstr(toplab)
      if (nhvols.eq.1) then
        write(toplab,"('HISTORY=',a)") histvol(ivol)
      else
        write(toplab,"('FIRST HISTORY VOL=',a)") histvol(1)
      endif
      call wrlab(toplab(1:lenstr(toplab)),xmid,vp(3)-xoff,
     +  .010)
c
c Put height axis on right if pressure on y-axis
c
      if (ihtsc.le.0) then
        yaxht(:) = 0.
        do k=izprange(1),izprange(2)
          izp = k-izprange(1)+1
          do i=1,ntms
            yaxht(izp) = yaxht(izp) + hts(i,k,l)
          enddo
          yaxht(izp) = yaxht(izp) / float(ntms)
        enddo
        rnd = 10.
        if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5.
        call altyax(ny,yaxht,gcmzp(izprange(1)),rnd,6)
      endif
c
c Wrap it up: 
c (make ascii file if needed)
c
      call frame
      iframe = iframe+1
      if (ihtsc.le.0) then
        if (iwrascii.gt.0) then
          call wrascii(iwrascii,luascii,plt,ntms,ny,utinc,gcmzp,
     +      'UT (HRS)','LN(P0/P)',histvol(ivol),fieldlab,iframe,
     +      rec80,'tigcmloc',dirascii)
        endif
        write(6,"('loc frame ',i4,' field ',a,
     +    ' lat,lon=',f6.1,',',f7.1,' zprange=',2f8.2)") 
     +    iframe,dopname(ip)(1:lenstr(dopname(ip))),gloc(1,l),
     +    gloc(2,l),zprange
      else
        if (iwrascii.gt.0) then
          call wrascii(iwrascii,luascii,pltht,ntms,ny,utinc,htscale,
     +      'UT (HRS)','HEIGHT (KM)',histvol(ivol),fieldlab,iframe,
     +      rec80,'tigcmloc',dirascii)
        endif
        write(6,"('loc frame ',i4,' field ',a,
     +    ' lat,lon=',f6.1,',',f7.1,' htscale=',f6.1,' to ',
     +    f6.1)") iframe,dopname(ip)(1:lenstr(dopname(ip))),
     +      gloc(1,l),gloc(2,l),htscale(1),htscale(nhtscale)
      endif
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine pvevu(hts,ihtsc,l)
c
c Contour ve,vu (from mkuivi where ionvel=2):
c   vemag and vumag are saved only if ionvel=2. They are Art's ion drifts
c     perpendicular to mag field where ve is +east and vu is +up
c
      include 'tigcmloc.h'
      character*56 toplab,fieldlab,veuname(2)
      character*80 rec80
      dimension plt(ntms,kmx),plt0(ntms,kmx),vp(4)
      dimension hts(ntms,kmx,nloc),pltht(ntms,nhtscale),yaxht(kmx)
      dimension log10v(2)
      data vp /.12,.86,.26,.91/
      data log10v/0,0/
      data veuname
     + / 'VE (M/S PERPENDICULAR TO MAG FIELD, +EAST)              ',
     +   'VU (M/S PERPENDICULAR TO MAG FIELD, +UP)                '/
c
      xmid = 0.5*(vp(1)+vp(2))
      do ip=1,2	! ve,vu
        if (ihtsc.le.0) then 
          do k=izprange(1),izprange(2)
            kk = k-izprange(1)+1
            if (ip.eq.1) call getvevuk(veutkmx,ntms,kmx,nloc,
     +        plt,k,kk,l)
            if (ip.eq.2) call getvevuk(vuutkmx,ntms,kmx,nloc,
     +        plt,k,kk,l)
          enddo
          ny = izprange(2)-izprange(1)+1
          if (logloc.gt.0.and.log10v(ip).gt.0) 
     +      call log10f(plt,ntms*ny,1.e-20,cpspval)
        else
          do k=1,kmx
            if (ip.eq.1) call getvevuk(veutkmx,ntms,kmx,nloc,
     +        plt0,k,k,l)
            if (ip.eq.2) call getvevuk(vuutkmx,ntms,kmx,nloc,
     +        plt0,k,k,l)
          enddo
          call cuthtint(plt0,hts(1,1,l),ntms,kmx,pltht,htscale,
     +      nhtscale,log10v(ip),cpspval,ier,1)
          ny = nhtscale
          if (logloc.gt.0.and.log10v(ip).gt.0) 
     +      call log10f(pltht,ntms*ny,1.e-20,cpspval)
        endif
c
c Contour:
c
        if (icolor.le.0) then
          if (ihtsc.le.0) then
            call contour(plt,ntms,ntms,ny,0.,1.,0.)
          else
            call contour(pltht,ntms,ntms,ny,0.,1.,0.)
          endif
        else
          if (ihtsc.le.0) then
            call conclr(plt,ntms,ntms,ny,0.,1.,0.)
          else
            call conclr(pltht,ntms,ntms,ny,0.,1.,0.)
          endif
        endif
c
c Add axes labels:
c
        isltax = 0
        if (gloc(1,l).ne.gmflag.and.gloc(2,l).ne.gmflag) isltax = 1
        if (ihtsc.le.0) then 
          call labutxy(mtimes,ntms,gcmzp(izprange(1)),ny,'ZP',
     +      0.,isltax,gloc(2,l))
        else
          call labutxy(mtimes,ntms,htscale,nhtscale,'HEIGHT (KM)',
     +      0.,isltax,gloc(2,l))
        endif
c
c Field label at top:
c
        call clearstr(toplab)
        if (logloc.gt.0.and.log10v(ip).gt.0) then
          len = lenstr(veuname(ip))
          if (len+6.le.56) then
            write(toplab,"('LOG10 ',a)") veuname(ip)(1:len)
          endif
        else
          toplab = veuname(ip)
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+.07,0.)
        fieldlab = toplab
c
c Location label (may be at selected lat,lon, or at global means, or
c zonal means at selected latitude):
c
        call clearstr(toplab)
        if (locname(1).eq.'xxxxxxxxxxxxxxxx') then
          write(toplab,"('LAT,LON = ',f6.2,f7.2)") gloc(1,l),
     +      gloc(2,l)
          xoff = .25
        else
          write(toplab,"('LAT,LON = ',f6.2,f7.2,' (',a,')')") 
     +      gloc(1,l),gloc(2,l),locname(l)(1:lenstr(locname(l)))
          xoff = .25
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,vp(4)+0.03,0.)
        call clearstr(rec80)
        write(rec80,"(a)") toplab(1:lenstr(toplab))
c
c History volume label at bottom:
c
        call clearstr(toplab)
        if (nhvols.eq.1) then
          write(toplab,"('HISTORY=',a)") histvol(ivol)
        else
          write(toplab,"('FIRST HISTORY VOL=',a)") histvol(1)
        endif
        call wrlab(toplab(1:lenstr(toplab)),xmid,vp(3)-xoff,
     +    .010)
c
c Put height axis on right if pressure on y-axis
c
        if (ihtsc.le.0) then
          yaxht(:) = 0.
          do k=izprange(1),izprange(2)
            izp = k-izprange(1)+1
            do i=1,ntms
              yaxht(izp) = yaxht(izp) + hts(i,k,l)
            enddo
            yaxht(izp) = yaxht(izp) / float(ntms)
          enddo
          rnd = 10.
          if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5.
          call altyax(ny,yaxht,gcmzp(izprange(1)),rnd,6)
        endif
c
c Wrap it up: 
c (make ascii file if needed)
c
        call frame
        iframe = iframe+1
        if (ihtsc.le.0) then
          if (iwrascii.gt.0) then
            call wrascii(iwrascii,luascii,plt,ntms,ny,utinc,gcmzp,
     +        'UT (HRS)','LN(P0/P)',histvol(ivol),fieldlab,iframe,
     +        rec80,'tigcmloc',dirascii)
          endif
          write(6,"('loc frame ',i4,' field ',a,
     +      ' lat,lon=',f6.1,',',f7.1,' zprange=',2f8.2)") 
     +      iframe,veuname(ip)(1:lenstr(veuname(ip))),gloc(1,l),
     +      gloc(2,l),zprange
        else
          if (iwrascii.gt.0) then
            call wrascii(iwrascii,luascii,pltht,ntms,ny,utinc,htscale,
     +        'UT (HRS)','HEIGHT (KM)',histvol(ivol),fieldlab,iframe,
     +        rec80,'tigcmloc',dirascii)
          endif
          write(6,"('loc frame ',i4,' field ',a,
     +      ' lat,lon=',f6.1,',',f7.1,' htscale=',f6.1,' to ',
     +      f6.1)") iframe,veuname(ip)(1:lenstr(veuname(ip))),
     +        gloc(1,l),gloc(2,l),htscale(1),htscale(nhtscale)
        endif
      enddo	! ve,vu
      return
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getvevuk(v,ntms,kmx,nloc,plt,kv,kp,l)
      real v(ntms,kmx,nloc),plt(ntms,kmx)
      plt(:,kp) = v(:,kv,l)
      return
      end