From: Stephen Bougher <[email-masked]>
Date: Fri, 24 Feb 1995 16:49:38 -0700

Dear Ben,

	How is everything with you? Have you had good snow in Jamestown yet?
I have a question about modifying pltlat.f for prints of many fields
at a single latitude slice. This is my alternative to wrascii which is
not detailed enough for me. I can correctly obtain the fields and fill
an array for printing on z-surfaces. However, obtaining fields on
interpolated height surfaces is not straightforward. For instance,
pltht(imx,ny) is what I am assuming should be available from cuthtint.
However, it is already dimensioned in a different way. What is going
on here? Prints on a constant height surface would be very useful to
me. Thanks for your help!

	The pltlat.f file below has been modified for these prints.
	The z-prints work; the height prints do not work.

--Steve Bougher
  (602-621-4900)

c
c------------------------------------------------------------------
c Begin file /home/sting/foster/tigcmproc/pltlat.f
c------------------------------------------------------------------
c
      subroutine pltlat(utc,fields,ihtsc)
c
c Contour latitude slices:
c (lon on x-axis, zp on y-axis if ihtsc <=0, 
c  linear height scale if ihtsc > 0)
c
      include 'vtgcmproc.h'
      character*56 toplab,fieldlab
      character*80 rec80
      dimension plt(imx,kmx),viewport(4),yaxht(kmx),
     +  plt0(imx,kmx),fields(imx,kmx,jmx,nfget),hts(imx,kmx),
     +  RITET(imx,kmx,nftot+2)
c    +  RITET(imx,kmx,nftot+2),pltht(imx,kmx),htscale(kmx)
      pointer(ppltht,pltht(1))
      data viewport /.15,.89,.26,.91/
c
      if (nlat.le.0) then
        write(6,"('pltlat: no selected latitudes -- returning')")
        return
      endif
      xmid = 0.5*(viewport(1)+viewport(2))
      call cpseti('SET',0)
      call cpseti('MAP',0)
      call cpsetr('XC1',gcmlon(1))
      call cpsetr('XCM',gcmlon(imx))
      call cpsetr('ILX',xmid)
      call cpsetr('ILY',-.16)
      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),
     +    gcmlon(1),gcmlon(imx),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),
     +    gcmlon(1),gcmlon(imx),htscale(1),htscale(nhtscale),1)
        call hpalloc(ppltht,imx*nhtscale,ier,1)
        if (ier.ne.0) then
          write(6,"('pltlat: hpalloc error = ',i6,' for pltht')") ier
          stop 'ppltht'
        endif
      endif
c
c Selected latitude loop:
c
      do 100 j=1,nlat
        ixlat = ixfind(gcmlat,jmx,slat(j),dlat)
        if (ixlat.le.0) then
          write(6,"('>>> pltlat: bad latitude=',f10.3,' j=',i3,
     +      ' skipping this lat')") slat(j),j
          goto 100
        endif
        rlat = gcmlat(ixlat)
c
c Field loop: 
c
        write(6,"(' ')")
c
C Index of ritet array
c
        irite=0
        do 200 ip=1,nftot
          if (ifplt(ip).le.0) goto 200
          if (ip.eq.ixunvn) goto 200
	  irite=irite+1
c
c Get field and interpolate to height scale if necessary:
c
          call getlat(fields,plt0,ixlat,ilog(ip),ip)
          if (ihtsc.le.0) then 
            call getlat(fields,plt0,ixlat,ilog(ip),ip)
            call getlat(fields,hts,ixlat,0,ixz)
            do k=izprange(1),izprange(2)
              kk = k-izprange(1)+1
              plt(:,kk) = plt0(:,k)
            enddo
            ny = izprange(2)-izprange(1)+1
          else
            call getlat(fields,plt0,ixlat,0,ip)
            call getlat(fields,hts,ixlat,0,ixz)
            call cuthtint(plt0,hts,imx,kmx,pltht,htscale,nhtscale,
     +        ilog(ip),cpspval,ier,0)
            if (ilog(ip).gt.0) call log10f(pltht,imx*nhtscale,1.e-20,
     +        cpspval)
            ny = nhtscale
          endif
c
c --------------------------------------------------------------------
c FILL PRINTS ARRAY: SINGLE LAT SLICE (I,K) FOR EACH FIELD QUERIED
c I=1,IMX (FULL) and K=1,NY (AS SPECIFIED)
C GENERAL FILL OF RITET ARRAY APPLICABLE TO MANY CASES
c --------ZP ---------------------------------------------------------
c
          if (ihtsc.le.0) then 
	    do i=1,imx
	    do k=1,ny
		  RITET(i,k,nftot)=-10.+(k+izprange(1)-2)*0.5 
		  RITET(i,k,nftot+1)=5.0E-03*EXP(-RITET(i,k,nftot))
		  RITET(i,k,nftot+2)=hts(i,k)
		  RITET(i,k,irite)=plt(i,k)
	    enddo
	    enddo
          else
c ------HTS-----------------------------------------------------------
c           do i=1,imx
c           do k=1,ny
c	          RITET(i,k,nftot+2)=htscale(k)
c	          RITET(i,k,irite)=pltht(i,k)
c           enddo
c           enddo
          endif
c --------------------------------------------------------------------
c
c Contour:
c
          if (icolor.le.0) then
            if (ihtsc.le.0) then
              call contour(plt,imx,imx,ny,cint(ip),cmin(ip),cmax(ip))
            else
              call contour(pltht,imx,imx,ny,cint(ip),cmin(ip),cmax(ip))
            endif
          else
            if (ihtsc.le.0) then
              call conclr(plt,imx,imx,ny,cint(ip),cmin(ip),cmax(ip))
            else
              call conclr(pltht,imx,imx,ny,cint(ip),cmin(ip),cmax(ip))
            endif
          endif
c
c Add axes labels:
c
          if (ihtsc.le.0) then 
            call labrect(gcmlon,imx,gcmzp(izprange(1)),ny,'LONGITUDE',
     +        'ZP',0.)
          else
            call labrect(gcmlon,imx,htscale(1),ny,'LONGITUDE',
     +        'HEIGHT (KM)',0.)
          endif
c
c Field label at top:
c
          call clearstr(toplab)
          if (ilog(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 Lat, ut at top below field label:
c
          call clearstr(toplab)
          write(toplab,"('LAT=',f8.2,' UT=',f6.2)") rlat,utc
          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 label at bottom:
c
          call clearstr(toplab)
          write(toplab,"('HISTORY=',a)") histvol(ivol)
          call wrlab(toplab(1:lenstr(toplab)),xmid,viewport(3)-0.15,
     +      .012)
c
c Put height axis on right if pressure on y-axis
c (use zonally averaged heights):
c
          if (ixfld(ixz).le.0.or.ihtsc.gt.0) goto 300
          yaxht(:) = 0.
          do k=izprange(1),izprange(2)
            izp = k-izprange(1)+1
            do i=1,imx-1
              yaxht(izp) = yaxht(izp) + fields(i,k,ixlat,ixfld(ixz))
            enddo
            yaxht(izp) = yaxht(izp) / float(imx-1)
          enddo 
          rnd = 10.
          if (gcmzp(izprange(2))-gcmzp(izprange(1)).le.5.) rnd = 5.
          call altyax(ny,yaxht,gcmzp(izprange(1)),rnd,6)
c
c Wrap it up (write ascii data file if needed):
c
 300      continue
          call frame
          iframe = iframe+1
          if (ihtsc.le.0) then
            if (iwrascii.gt.0) then
              call wrascii(iwrascii,luascii,plt,imx,ny,gcmlon,
     +          gcmzp(izprange(1)),'LONGITUDE','LN(P0/P)',
     +          histvol(ivol),fieldlab,iframe,rec80,'tigcmproc',
     +          dirascii)
            endif
            write(6,"('pltlat: frame ',i4,' field ',a8,' lat=',f8.2,
     +        ' zprange=',2f8.2)") iframe,labshort(ip),rlat,zprange
          else
            if (iwrascii.gt.0) then
              call wrascii(iwrascii,luascii,pltht,imx,ny,gcmlon,
     +          htscale,'LONGITUDE','HEIGHT (KM)',histvol(ivol),
     +          fieldlab,iframe,rec80,'tigcmproc',dirascii)
            endif
            write(6,"('pltlat: frame ',i4,' field ',a8,' lat=',f8.2,
     +        ' ht scale=',f7.2,' to ',f7.2)") iframe,
     +        labshort(ip),rlat,htscale(1),htscale(nhtscale)
          endif
c
c End field loop: ip=1,nftot  // irite= # of fields to be plotted
 200    continue
c-----------------------------------------------------------------------
c GROUP PRINTS AT EACH QUERIED LAT SLICE : 
c --- FOR EACH LAT SLICE CHOSEN FOR STUDY
c --- EVERY 1-HR INTERVALS
c --- SPANS NY LEVELS (CONSTANT ZP OR HEIGHT)
c --- PRINTS IRITE FIELDS (LESS THAN OR EQUAL TO 14-FIELDS)
c MODIFY 1002 FORMAT FOR EACH REQUIRED PRINT GRP (14-FIELDS MAX)
c-----------------------------------------------------------------------
      if (ihtsc.le.0) then 
C **** PRESSURE LEVEL PRINTS ****
      do 400 i=3,75,3
	 xlt=0.+(i-3.)/3.
         WRITE(20,1000) rlat,xlt
      do 400 k=1,ny
         WRITE(20,1001)RITET(i,k,nftot),RITET(i,k,nftot+1),
     &   RITET(i,k,nftot+2)
         WRITE(20,1002)(RITET(i,k,nprt),nprt=1,7)
         WRITE(20,1002)(RITET(i,k,nprt),nprt=8,14)
 400  continue
      else
c-----------------------------------------------------------------------
C **** ALTITUDE LEVEL PRINTS ****
c     do 500 i=1,imx
c        xlt=0.+(i-3.)/3.
c        WRITE(20,1000) rlat,xlt
c     do 500 k=1,ny
c        WRITE(20,1003)(RITET(i,k,nftot+2)
c        WRITE(20,1002)(RITET(i,k,nprt),nprt=1,7)
c        WRITE(20,1002)(RITET(i,k,nprt),nprt=8,14)
c500  continue
      endif
c-----------------------------------------------------------------------
 1000 FORMAT(//20X,*LAT=*,F10.2,2X,*LT=*,F10.2//)
 1001 FORMAT(1X,F10.2,E12.3,F10.2)
 1002 FORMAT(3X,7E11.3)
 1003 FORMAT(1X,F10.2)
c-----------------------------------------------------------------------
c
c End selected latitude loop
 100  continue
      if (ihtsc.gt.0) call hpdeallc(ppltht,ier,1)
c
      return
      end