
      subroutine getglb(fields,plt,ixpl,ht,ip)
      include 'timesproc.h'
      dimension plt0(imx,kmx,jmx),plt(imx,jmx),plt1(imx,jmx),
     +  fields(imx,kmx,jmx,nfget)
      pointer (pfscr,fscr(1)) ! scratch
c
c Return desired global field (at ip) in plt(imx,jmx).
c ixpl = index to desired pressure. If ixpl=0, assume interpolation to
c        constant height surface ht
c
c Field from history:
c
      if (ip.le.nfhist) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ifget(ip))
        else
          call glbhtint(fields(1,1,1,ifget(ip)),
     +      fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt,ht,1,
     +      itimelog(ip),cpspval,ier,0)
        endif
c
c ion drifts:
      elseif (ip.ge.itxui.and.ip.le.itxwi) then
        if (ixpl.gt.0) then
          if (ip.eq.itxui) plt(:,:) = glbui(:,ixpl,:)
          if (ip.eq.itxvi) plt(:,:) = glbvi(:,ixpl,:)
          if (ip.eq.itxwi) plt(:,:) = glbwi(:,ixpl,:)
        else
          if (ip.eq.itxui) 
     +      call glbhtint(glbui,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,
     +        plt,ht,1,itimelog(ip),cpspval,ier,0)
          if (ip.eq.itxvi) 
     +      call glbhtint(glbvi,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,
     +        plt,ht,1,itimelog(ip),cpspval,ier,0)
          if (ip.eq.itxwi) 
     +      call glbhtint(glbwi,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,
     +        plt,ht,1,itimelog(ip),cpspval,ier,0)
        endif
c
c O/O2:
c
      elseif (ip.eq.itxoo2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ifget(itxo1)) / 
     +      fields(:,ixpl,:,ifget(itxo2))
        else
          plt0(:,:,:) = fields(:,:,:,ifget(itxo1)) / 
     +      fields(:,:,:,ifget(itxo2))
          call glbhtint(plt0,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt,
     +      ht,1,itimelog(ip),cpspval,ier,0)
        endif
c
c O/N2:
c
      elseif (ip.eq.itxon2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ifget(itxo1)) / 
     +      fields(:,ixpl,:,ifget(itxn2))
        else
          plt0(:,:,:) = fields(:,:,:,ifget(itxo1)) / 
     +      fields(:,:,:,ifget(itxn2))
          call glbhtint(plt0,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt,
     +      ht,1,itimelog(ip),cpspval,ier,0)
        endif
c
c O2/N2:
c
      elseif (ip.eq.itxo2n2) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ifget(itxo2)) / 
     +      fields(:,ixpl,:,ifget(itxn2))
        else
          plt0(:,:,:) = fields(:,:,:,ifget(itxo2)) / 
     +      fields(:,:,:,ifget(itxn2))
          call glbhtint(plt0,fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt,
     +      ht,1,itimelog(ip),cpspval,ier,0)
        endif
c
c O+O2+N2 (rho):
c
      elseif (ip.eq.itxrho) then
        if (ixpl.gt.0) then
          plt(:,:) = fields(:,ixpl,:,ifget(itxo1)) +
     +               fields(:,ixpl,:,ifget(itxo2)) +
     +               fields(:,ixpl,:,ifget(itxn2))
        else
          call glbhtint(fields(1,1,1,ifget(itxo1)),
     +      fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt,ht,1,itimelog(ip),
     +      cpspval,ier,0)
          call glbhtint(fields(1,1,1,ifget(itxo2)),
     +      fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt1,ht,1,
     +      itimelog(ip),cpspval,ier,0)
          do j=1,jmx
            do i=1,imx
              if (plt1(i,j).ne.cpspval.and.plt(i,j).ne.cpspval) then
                plt(i,j) = plt(i,j) + plt1(i,j)
              else
                plt(i,j) = cpspval
              endif
            enddo
          enddo
          call glbhtint(fields(1,1,1,ifget(itxn2)),
     +      fields(1,1,1,ifget(itxz)),imx,kmx,jmx,plt1,ht,1,
     +      itimelog(ip),cpspval,ier,0)
          do j=1,jmx
            do i=1,imx
              if (plt1(i,j).ne.cpspval.and.plt(i,j).ne.cpspval) then
                plt(i,j) = plt(i,j) + plt1(i,j)
              else
                plt(i,j) = cpspval
              endif
            enddo
          enddo
        endif
      elseif (ip.eq.itxhtot) then
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
        if (ixpl.gt.0) then
          plt(:,:) = 2.*fields(:,ixpl,:,ifget(itxh2o)) +
     +                  fields(:,ixpl,:,ifget(itxh)) +
     +                  fields(:,ixpl,:,ifget(itxoh)) +
     +                  fields(:,ixpl,:,ifget(itxho2)) +
     +               4.*fields(:,ixpl,:,ifget(itxch4)) +
     +               2.*fields(:,ixpl,:,ifget(itxh2))
          plt(:,:) = plt(:,:) / (fields(:,ixpl,:,ifget(itxo1)) +
     +                           fields(:,ixpl,:,ifget(itxo2)) +
     +                           fields(:,ixpl,:,ifget(itxn2)))
        else
          write(6,"('>>> getfield: HTOT not yet available in height')")
        endif
c
c E5577 emissions need t,o2,o,n2:
c
      elseif (ip.eq.ixe5577) then
        if (ixpl.gt.0) then
          do j=1,jmx    
            call mke5577(fields(1,ixpl,j,ifget(itxt)),
     +                   fields(1,ixpl,j,ifget(itxo2)),
     +                   fields(1,ixpl,j,ifget(itxo1)),
     +                   fields(1,ixpl,j,ifget(itxn2)),
     +                   imx,plt(1,j),cpspval)
          enddo
        else
          do j=1,jmx
            call cuthtint(fields(1,1,j,ifget(itxt)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1,ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,2),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo1)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,3),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxn2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,4),ht,1,0,cpspval,ier,0)
            call mke5577(plt1,plt1(1,2),plt1(1,3),plt1(1,4),imx,
     +        plt(1,j),cpspval)
          enddo
        endif
c
c E6300 emissions need t,te,o2,o,n2,o2+,ne:
c
      elseif (ip.eq.ixe6300) then
        if (ixpl.gt.0) then
          do j=1,jmx    
            call mke6300(fields(1,ixpl,j,ifget(itxt)),
     +                   fields(1,ixpl,j,ifget(itxte)),
     +                   fields(1,ixpl,j,ifget(itxo2)),
     +                   fields(1,ixpl,j,ifget(itxo1)),
     +                   fields(1,ixpl,j,ifget(itxn2)),
     +                   fields(1,ixpl,j,ifget(itxo2p)),
     +                   fields(1,ixpl,j,ifget(itxne)),
     +                   imx,plt(1,j),cpspval)
          enddo
        else
          do j=1,jmx
            call cuthtint(fields(1,1,j,ifget(itxt)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1,ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxte)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,2),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,3),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo1)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,4),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxn2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,5),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo2p)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,6),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxne)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,7),ht,1,0,cpspval,ier,0)
            call mke6300(plt1,plt1(1,2),plt1(1,3),plt1(1,4),
     +        plt1(1,5),plt1(1,6),plt1(1,7),imx,plt(1,j),cpspval)
          enddo
        endif
c
c EO200 emissions need t,o2,o,n2:
c
      elseif (ip.eq.ixeo200) then
        if (ixpl.gt.0) then
          do j=1,jmx    
            call mkeo200(fields(1,ixpl,j,ifget(itxt)),
     +                   fields(1,ixpl,j,ifget(itxo2)),
     +                   fields(1,ixpl,j,ifget(itxo1)),
     +                   fields(1,ixpl,j,ifget(itxn2)),
     +                   imx,plt(1,j),cpspval)
          enddo
        else
          do j=1,jmx
            call cuthtint(fields(1,1,j,ifget(itxt)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1,ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,2),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxo1)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,3),ht,1,0,cpspval,ier,0)
            call cuthtint(fields(1,1,j,ifget(itxn2)),
     +                    fields(1,1,j,ifget(itxz)),imx,kmx,
     +                    plt1(1,4),ht,1,0,cpspval,ier,0)
            call mkeo200(plt1,plt1(1,2),plt1(1,3),plt1(1,4),imx,
     +        plt(1,j),cpspval)
          enddo
        endif
      else
        write(6,"('>>> getglb: unknown ip=',i3,' returning zeroes')") 
     +    ip
        plt(:,:) = 0.
      endif
c
c Take log10 if desired:
      if (log10map.gt.0.and.itimelog(ip).gt.0)
     +  call log10f(plt,imx*jmx,1.e-20,cpspval)
c
      return
      end
c
c--------------------------------------------------------------------
c
      subroutine getlon(fields,plt,ixlon,log,ip)
c
c Return lon slice field in plt(jmx,kmx) at longitude index ixlon.
c If ixlon=int(zmflag), return zonal means.
c If log > 0, return log10 of field. Note calcmean finds zonal means
c   and takes log10 of zonal means if necessary. log10f takes log10
c   of selected lon slices.
c
      include 'timesproc.h'
      dimension plt(jmx,kmx),rimx(imx),fields(imx,kmx,jmx,nfget)
      dimension fjmx10(jmx,10)
c
c Field from history:
c
      if (ip.le.nfhist) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ifget(ip))
            enddo
          enddo
          if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(ip))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          enddo
        endif
c
c ion drifts:
      elseif (ip.ge.itxui.and.ip.le.itxwi) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              if (ip.eq.itxui) plt(j,k) = glbui(ixlon,k,j)
              if (ip.eq.itxvi) plt(j,k) = glbvi(ixlon,k,j)
              if (ip.eq.itxwi) plt(j,k) = glbwi(ixlon,k,j)
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              if (ip.eq.itxui) rimx(:) = glbui(:,k,j)
              if (ip.eq.itxvi) rimx(:) = glbvi(:,k,j)
              if (ip.eq.itxwi) rimx(:) = glbwi(:,k,j)
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval,0)
            enddo
          enddo
        endif
c
c O/O2:
c
      elseif (ip.eq.itxoo2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ifget(itxo1)) / 
     +                   fields(ixlon,k,j,ifget(itxo2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(itxo1)) / 
     +          fields(:,k,j,ifget(itxo2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          enddo
        endif
c
c O/N2:
c
      elseif (ip.eq.itxon2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ifget(itxo1)) / 
     +                   fields(ixlon,k,j,ifget(itxn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(itxo1)) / 
     +                  fields(:,k,j,ifget(itxn2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          enddo
        endif
c
c O2/N2:
c
      elseif (ip.eq.itxo2n2) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ifget(itxo2)) / 
     +                   fields(ixlon,k,j,ifget(itxn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(itxo2)) / 
     +                  fields(:,k,j,ifget(itxn2))
              plt(j,k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          enddo
        endif
c
c rho (o+o2+n2):
c
      elseif (ip.eq.itxrho) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = fields(ixlon,k,j,ifget(itxo1)) + 
     +                   fields(ixlon,k,j,ifget(itxo2)) +
     +                   fields(ixlon,k,j,ifget(itxn2))
            enddo
          enddo
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(itxo1))
              plt(j,k) = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxo2))
              plt(j,k) = plt(j,k) + calcmean(rimx,imx-1,0,1.e-20,
     +          cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxn2))
              plt(j,k) = plt(j,k) + calcmean(rimx,imx-1,0,1.e-20,
     +          cpspval,0)
            enddo
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
      elseif (ip.eq.itxhtot) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            do j=1,jmx
              plt(j,k) = 2.*fields(ixlon,k,j,ifget(itxh2o)) +
     +                      fields(ixlon,k,j,ifget(itxh)) +
     +                      fields(ixlon,k,j,ifget(itxoh)) +
     +                      fields(ixlon,k,j,ifget(itxho2)) +
     +                   4.*fields(ixlon,k,j,ifget(itxch4)) +
     +                   2.*fields(ixlon,k,j,ifget(itxh2))
              plt(j,k) = plt(j,k) / (fields(ixlon,k,j,ifget(itxo1)) +
     +                               fields(ixlon,k,j,ifget(itxo2)) +
     +                               fields(ixlon,k,j,ifget(itxn2)))

            enddo
          enddo
c
c Calculate means of each component while building htot:
c
        else
          do k=1,kmx
            do j=1,jmx
              rimx(:) = fields(:,k,j,ifget(itxh2o))
              plt(j,k) = 2.*calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxh))
              plt(j,k) = plt(j,k) + 
     +                   calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxoh))
              plt(j,k) = plt(j,k) + 
     +                   calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxho2))
              plt(j,k) = plt(j,k) + 
     +                   calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxch4))
              plt(j,k) = plt(j,k) + 
     +                   4.*calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxh2))
              plt(j,k) = plt(j,k) + 
     +                   2.*calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxo2))
              tot = calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxo1))
              tot = tot + calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              rimx(:) = fields(:,k,j,ifget(itxn2))
              tot = tot + calcmean(rimx,imx-1,0,1.e-20,cpspval,0)
              plt(j,k) = plt(j,k) / tot 
            enddo
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c E5577 emission:
c
      elseif (ip.eq.ixe5577) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            fjmx10(:,1) = fields(ixlon,k,:,ifget(itxt))
            fjmx10(:,2) = fields(ixlon,k,:,ifget(itxo2))
            fjmx10(:,3) = fields(ixlon,k,:,ifget(itxo1))
            fjmx10(:,4) = fields(ixlon,k,:,ifget(itxn2))
            call mke5577(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        jmx,plt(1,k),cpspval)
          enddo
        else ! zonal means e5577
          do k=1,kmx
            do j=1,jmx
              fjmx10(j,1) = calcmean(fields(1,k,j,ifget(itxt)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,2) = calcmean(fields(1,k,j,ifget(itxo2)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,3) = calcmean(fields(1,k,j,ifget(itxo1)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,4) = calcmean(fields(1,k,j,ifget(itxn2)),
     +          imx-1,0,1.e-20,cpspval,0)
            enddo
            call mke5577(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        jmx,plt(1,k),cpspval)
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c E6300 emission:
c
      elseif (ip.eq.ixe6300) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            fjmx10(:,1) = fields(ixlon,k,:,ifget(itxt))
            fjmx10(:,2) = fields(ixlon,k,:,ifget(itxte))
            fjmx10(:,3) = fields(ixlon,k,:,ifget(itxo2))
            fjmx10(:,4) = fields(ixlon,k,:,ifget(itxo1))
            fjmx10(:,5) = fields(ixlon,k,:,ifget(itxn2))
            fjmx10(:,6) = fields(ixlon,k,:,ifget(itxo2p))
            fjmx10(:,7) = fields(ixlon,k,:,ifget(itxne))
            call mke6300(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        fjmx10(1,5),fjmx10(1,6),fjmx10(1,7),jmx,plt(1,k),cpspval)
          enddo
        else ! zonal means e6300
          do k=1,kmx
            do j=1,jmx
              fjmx10(j,1) = calcmean(fields(1,k,j,ifget(itxt)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,2) = calcmean(fields(1,k,j,ifget(itxte)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,3) = calcmean(fields(1,k,j,ifget(itxo2)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,4) = calcmean(fields(1,k,j,ifget(itxo1)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,5) = calcmean(fields(1,k,j,ifget(itxn2)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,6) = calcmean(fields(1,k,j,ifget(itxo2p)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,7) = calcmean(fields(1,k,j,ifget(itxne)),
     +          imx-1,0,1.e-20,cpspval,0)
            enddo
            call mke6300(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        fjmx10(1,5),fjmx10(1,6),fjmx10(1,7),jmx,plt(1,k),cpspval)
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
c
c EO200 emission:
c
      elseif (ip.eq.ixeo200) then
        if (ixlon.ne.ifix(zmflag)) then
          do k=1,kmx
            fjmx10(:,1) = fields(ixlon,k,:,ifget(itxt))
            fjmx10(:,2) = fields(ixlon,k,:,ifget(itxo2))
            fjmx10(:,3) = fields(ixlon,k,:,ifget(itxo1))
            fjmx10(:,4) = fields(ixlon,k,:,ifget(itxn2))
            call mkeo200(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        jmx,plt(1,k),cpspval)
          enddo
        else ! zonal means eo200
          do k=1,kmx
            do j=1,jmx
              fjmx10(j,1) = calcmean(fields(1,k,j,ifget(itxt)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,2) = calcmean(fields(1,k,j,ifget(itxo2)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,3) = calcmean(fields(1,k,j,ifget(itxo1)),
     +          imx-1,0,1.e-20,cpspval,0)
              fjmx10(j,4) = calcmean(fields(1,k,j,ifget(itxn2)),
     +          imx-1,0,1.e-20,cpspval,0)
            enddo
            call mkeo200(fjmx10,fjmx10(1,2),fjmx10(1,3),fjmx10(1,4),
     +        jmx,plt(1,k),cpspval)
          enddo
        endif
        if (log.gt.0) call log10f(plt,jmx*kmx,1.e-20,cpspval)
      else
        write(6,"('>>> getlon: unknown ip=',i3,' returning zeroes')") 
     +    ip
        plt(:,:) = 0.
      endif
      return
      end
c
c--------------------------------------------------------------------
c
      subroutine getlat(fields,plt,ixlat,log,ip)
      include 'timesproc.h'
      dimension plt(imx,kmx),fields(imx,kmx,jmx,nfget)
c
c Return field ip latitude slice at lat index ixlat. Return log10 of
c field if log > 0.
c
      if (ip.le.nfhist) then
        plt(:,:) = fields(:,:,ixlat,ifget(ip))
      elseif (ip.ge.itxui.and.ip.le.itxwi) then
        if (ip.eq.itxui) plt(:,:) = glbui(:,:,ixlat)
        if (ip.eq.itxvi) plt(:,:) = glbvi(:,:,ixlat)
        if (ip.eq.itxwi) plt(:,:) = glbwi(:,:,ixlat)
      elseif (ip.eq.itxoo2) then
        plt(:,:) = fields(:,:,ixlat,ifget(itxo1)) / 
     +             fields(:,:,ixlat,ifget(itxo2))
      elseif (ip.eq.itxon2) then
        plt(:,:) = fields(:,:,ixlat,ifget(itxo1)) / 
     +             fields(:,:,ixlat,ifget(itxn2))
      elseif (ip.eq.itxo2n2) then
        plt(:,:) = fields(:,:,ixlat,ifget(itxo2)) / 
     +             fields(:,:,ixlat,ifget(itxn2))
      elseif (ip.eq.itxrho) then
        plt(:,:) = fields(:,:,ixlat,ifget(itxo1)) +
     +             fields(:,:,ixlat,ifget(itxo2)) +
     +             fields(:,:,ixlat,ifget(itxn2))
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
      elseif (ip.eq.itxhtot) then
        plt(:,:) = 2.*fields(:,:,ixlat,ifget(itxh2o)) +
     +                fields(:,:,ixlat,ifget(itxh)) +
     +                fields(:,:,ixlat,ifget(itxoh)) +
     +                fields(:,:,ixlat,ifget(itxho2)) +
     +             4.*fields(:,:,ixlat,ifget(itxch4)) +
     +             2.*fields(:,:,ixlat,ifget(itxh2))
        plt(:,:) = plt(:,:) / (fields(:,:,ixlat,ifget(itxo1)) +
     +                         fields(:,:,ixlat,ifget(itxo2)) +
     +                         fields(:,:,ixlat,ifget(itxn2)))
      elseif (ip.eq.ixe5577) then
        do k=1,kmx
          call mke5577(fields(1,k,ixlat,ifget(itxt)),
     +                 fields(1,k,ixlat,ifget(itxo2)),
     +                 fields(1,k,ixlat,ifget(itxo1)),
     +                 fields(1,k,ixlat,ifget(itxn2)),
     +                 imx,plt(1,k),cpspval)
        enddo
      elseif (ip.eq.ixe6300) then
        do k=1,kmx
          call mke6300(fields(1,k,ixlat,ifget(itxt)),
     +                 fields(1,k,ixlat,ifget(itxte)),
     +                 fields(1,k,ixlat,ifget(itxo2)),
     +                 fields(1,k,ixlat,ifget(itxo1)),
     +                 fields(1,k,ixlat,ifget(itxn2)),
     +                 fields(1,k,ixlat,ifget(itxo2p)),
     +                 fields(1,k,ixlat,ifget(itxne)),
     +                 imx,plt(1,k),cpspval)
        enddo
      elseif (ip.eq.ixeo200) then
        do k=1,kmx
          call mkeo200(fields(1,k,ixlat,ifget(itxt)),
     +                 fields(1,k,ixlat,ifget(itxo2)),
     +                 fields(1,k,ixlat,ifget(itxo1)),
     +                 fields(1,k,ixlat,ifget(itxn2)),
     +                 imx,plt(1,k),cpspval)
        enddo
      else
        write(6,"('>>> getlat: unknown ip=',i3,' returning zeroes')") ip
        plt(:,:) = 0.
      endif
      if (log.gt.0) call log10f(plt,imx*kmx,1.e-20,cpspval) 
      return
      end
c
c--------------------------------------------------------------------
c
      subroutine getloc(fields,plt,ilat,ilon,log,ip)
      include 'timesproc.h'
      dimension plt(kmx),glb(imx,jmx),fields(imx,kmx,jmx,nfget),
     +  rimx(imx),plt0(kmx),fkmx10(kmx,10)
c
      ixlat = ilat
      ixlon = ilon
c
c Zonal means at a selected latitude:
c
      if (ixlat.gt.0.and.ixlon.eq.ifix(zmflag)) then
        if (ip.le.nfhist) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(ip))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.ge.itxui.and.ip.le.itxwi) then
          if (ip.eq.itxui) then
            do k=1,kmx
              rimx(:) = glbui(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          elseif (ip.eq.itxvi) then
            do k=1,kmx
              rimx(:) = glbvi(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          elseif (ip.eq.itxwi) then
            do k=1,kmx
              rimx(:) = glbwi(:,k,ixlat)
              plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
            enddo
          endif
        elseif (ip.eq.itxoo2) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxo1)) /
     +                fields(:,k,ixlat,ifget(itxo2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.itxon2) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxo1)) /
     +                fields(:,k,ixlat,ifget(itxn2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.itxo2n2) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxo2)) /
     +                fields(:,k,ixlat,ifget(itxn2))
            plt(k) = calcmean(rimx,imx-1,log,1.e-20,cpspval)
          enddo
        elseif (ip.eq.itxrho) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxo1))
            plt(k) = calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxo2))
            plt0(k) = calcmean(rimx,imx-1,0,1.e-20,cpspval)
          enddo
          plt(:) = plt(:) + plt0(:)
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxn2))
            plt0(k) = calcmean(rimx,imx-1,0,1.e-20,cpspval)
          enddo
          plt(:) = plt(:) + plt0(:)
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
        elseif (ip.eq.itxhtot) then
          do k=1,kmx
            rimx(:) = fields(:,k,ixlat,ifget(itxh2o))
            plt(k) = 2.*calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxh))
            plt(k) = plt(k)+calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxoh))
            plt(k) = plt(k)+calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxho2))
            plt(k) = plt(k)+calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxch4))
            plt(k) = plt(k)+4.*calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxh2))
            plt(k) = plt(k)+2.*calcmean(rimx,imx-1,0,1.e-20,cpspval)
c
            rimx(:) = fields(:,k,ixlat,ifget(itxo2))
            tot = calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxo1))
            tot = tot + calcmean(rimx,imx-1,0,1.e-20,cpspval)
            rimx(:) = fields(:,k,ixlat,ifget(itxn2))
            tot = tot + calcmean(rimx,imx-1,0,1.e-20,cpspval)
            plt(k) = plt(k) / tot
          enddo
c
c zonal means of e5577:
c
        elseif (ip.eq.ixe5577) then
          do k=1,kmx
            ftn = calcmean(fields(1,k,ixlat,ifget(itxt)),imx-1,0,
     +        1.e-20,cpspval)
            fo2 = calcmean(fields(1,k,ixlat,ifget(itxo2)),imx-1,0,
     +        1.e-20,cpspval)
            fo1 = calcmean(fields(1,k,ixlat,ifget(itxo1)),imx-1,0,
     +        1.e-20,cpspval)
            fn2 = calcmean(fields(1,k,ixlat,ifget(itxn2)),imx-1,0,
     +        1.e-20,cpspval)
            call mke5577(ftn,fo2,fo1,fn2,1,plt(k),cpspval)
          enddo
c
c zonal means of e6300:
c
        elseif (ip.eq.ixe6300) then
          do k=1,kmx
            ftn = calcmean(fields(1,k,ixlat,ifget(itxt)),imx-1,0,
     +        1.e-20,cpspval)
            fte = calcmean(fields(1,k,ixlat,ifget(itxte)),imx-1,0,
     +        1.e-20,cpspval)
            fo2 = calcmean(fields(1,k,ixlat,ifget(itxo2)),imx-1,0,
     +        1.e-20,cpspval)
            fo1 = calcmean(fields(1,k,ixlat,ifget(itxo1)),imx-1,0,
     +        1.e-20,cpspval)
            fn2 = calcmean(fields(1,k,ixlat,ifget(itxn2)),imx-1,0,
     +        1.e-20,cpspval)
            fo2p = calcmean(fields(1,k,ixlat,ifget(itxo2p)),imx-1,0,
     +        1.e-20,cpspval)
            fne = calcmean(fields(1,k,ixlat,ifget(itxne)),imx-1,0,
     +        1.e-20,cpspval)
            call mke6300(ftn,fte,fo2,fo1,fn2,fo2p,fne,1,plt(k),cpspval)
          enddo
c
c zonal means of eo200:
c
        elseif (ip.eq.ixeo200) then
          do k=1,kmx
            ftn = calcmean(fields(1,k,ixlat,ifget(itxt)),imx-1,0,
     +        1.e-20,cpspval)
            fo2 = calcmean(fields(1,k,ixlat,ifget(itxo2)),imx-1,0,
     +        1.e-20,cpspval)
            fo1 = calcmean(fields(1,k,ixlat,ifget(itxo1)),imx-1,0,
     +        1.e-20,cpspval)
            fn2 = calcmean(fields(1,k,ixlat,ifget(itxn2)),imx-1,0,
     +        1.e-20,cpspval)
            call mkeo200(ftn,fo2,fo1,fn2,1,plt(k),cpspval)
          enddo
        else
          write(6,"('>>> getloc: unknown ip=',i3,' returning zeroes')") 
     +      ip
          plt(:) = 0.
          return
        endif
c
c Selected lat,lon:
c
      elseif (ixlat.gt.0.and.ixlon.gt.0) then
        if (ip.le.nfhist) then
          plt(:) = fields(ixlon,:,ixlat,ifget(ip))
        elseif (ip.ge.itxui.and.ip.le.itxwi) then
          if (ip.eq.itxui) plt(:) = glbui(ixlon,:,ixlat)
          if (ip.eq.itxvi) plt(:) = glbvi(ixlon,:,ixlat)
          if (ip.eq.itxwi) plt(:) = glbwi(ixlon,:,ixlat)
        elseif (ip.eq.itxoo2) then
          plt(:) = fields(ixlon,:,ixlat,ifget(itxo1)) / 
     +             fields(ixlon,:,ixlat,ifget(itxo2))
        elseif (ip.eq.itxon2) then
          plt(:) = fields(ixlon,:,ixlat,ifget(itxo1)) / 
     +             fields(ixlon,:,ixlat,ifget(itxn2))
        elseif (ip.eq.itxo2n2) then
          plt(:) = fields(ixlon,:,ixlat,ifget(itxo2)) / 
     +             fields(ixlon,:,ixlat,ifget(itxn2))
        elseif (ip.eq.itxrho) then
          plt(:) = fields(ixlon,:,ixlat,ifget(itxo1)) +
     +             fields(ixlon,:,ixlat,ifget(itxo2)) +
     +             fields(ixlon,:,ixlat,ifget(itxn2))
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
        elseif (ip.eq.itxhtot) then
          plt(:) = 2.*fields(ixlon,:,ixlat,ifget(itxh2o)) +
     +                fields(ixlon,:,ixlat,ifget(itxh)) +
     +                fields(ixlon,:,ixlat,ifget(itxoh)) +
     +                fields(ixlon,:,ixlat,ifget(itxho2)) +
     +             4.*fields(ixlon,:,ixlat,ifget(itxch4)) +
     +             2.*fields(ixlon,:,ixlat,ifget(itxh2))
          plt0(:) = fields(ixlon,:,ixlat,ifget(itxo2)) +
     +              fields(ixlon,:,ixlat,ifget(itxo1)) +
     +              fields(ixlon,:,ixlat,ifget(itxn2))
          plt(:) = plt(:) / plt0(:)
c
c E5577 emission:
c
        elseif (ip.eq.ixe5577) then
          fkmx10(:,1) = fields(ixlon,:,ixlat,ifget(itxt)) 
          fkmx10(:,2) = fields(ixlon,:,ixlat,ifget(itxo2)) 
          fkmx10(:,3) = fields(ixlon,:,ixlat,ifget(itxo1)) 
          fkmx10(:,4) = fields(ixlon,:,ixlat,ifget(itxn2)) 
          call mke5577(fkmx10,fkmx10(1,2),fkmx10(1,3),fkmx10(1,4),
     +      kmx,plt,cpspval)
c
c E6300 emission:
c
        elseif (ip.eq.ixe6300) then
          fkmx10(:,1) = fields(ixlon,:,ixlat,ifget(itxt)) 
          fkmx10(:,2) = fields(ixlon,:,ixlat,ifget(itxte)) 
          fkmx10(:,3) = fields(ixlon,:,ixlat,ifget(itxo2)) 
          fkmx10(:,4) = fields(ixlon,:,ixlat,ifget(itxo1)) 
          fkmx10(:,5) = fields(ixlon,:,ixlat,ifget(itxn2)) 
          fkmx10(:,6) = fields(ixlon,:,ixlat,ifget(itxo2p)) 
          fkmx10(:,7) = fields(ixlon,:,ixlat,ifget(itxne)) 
          call mke6300(fkmx10,fkmx10(1,2),fkmx10(1,3),fkmx10(1,4),
     +      fkmx10(1,5),fkmx10(1,6),fkmx10(1,7),kmx,plt,cpspval)
c
c EO200 emission:
c
        elseif (ip.eq.ixeo200) then
          fkmx10(:,1) = fields(ixlon,:,ixlat,ifget(itxt)) 
          fkmx10(:,2) = fields(ixlon,:,ixlat,ifget(itxo2)) 
          fkmx10(:,3) = fields(ixlon,:,ixlat,ifget(itxo1)) 
          fkmx10(:,4) = fields(ixlon,:,ixlat,ifget(itxn2)) 
          call mkeo200(fkmx10,fkmx10(1,2),fkmx10(1,3),fkmx10(1,4),
     +      kmx,plt,cpspval)
        else
          write(6,"('>>> getloc: unknown ip=',i3,' returning zeroes')") 
     +      ip
          plt(:) = 0.
          return
        endif
        if (log.gt.0) call log10f(plt,kmx,1.e-20,cpspval) 
c
c Global means:
c
      elseif (ixlat.eq.0.and.ixlon.eq.0) then
        if (ip.le.nfhist) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(ip))
            plt(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
          enddo
        elseif (ip.ge.itxui.and.ip.le.itxwi) then
          do k=1,kmx
            if (ip.eq.itxui) glb(:,:) = glbui(:,k,:)
            if (ip.eq.itxvi) glb(:,:) = glbvi(:,k,:)
            if (ip.eq.itxwi) glb(:,:) = glbwi(:,k,:)
            plt(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
          enddo
        elseif (ip.eq.itxoo2) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxo1)) / 
     +                 fields(:,k,:,ifget(itxo2))
            plt(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
          enddo
        elseif (ip.eq.itxon2) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxo1)) / 
     +                 fields(:,k,:,ifget(itxn2))
            plt(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
          enddo
        elseif (ip.eq.itxo2n2) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxo2)) / 
     +                 fields(:,k,:,ifget(itxn2))
            plt(k) = fglbm(glb,imx,jmx,gcmlat,dlat,dlon,cpspval)
          enddo
        elseif (ip.eq.itxrho) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxo1))
            plt(k) = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxn2))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
          enddo
c
c Htot = 2*h2o + h + oh + ho2 + 4*ch4 + 2*h2, where all species
c   densities must be in cm-3, then convert to mixing ratios:
c (added 2/94 for Bill Sharp)
c
        elseif (ip.eq.itxhtot) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxh2o))
            plt(k) = 2.*calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxh))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxoh))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxho2))
            plt(k) = plt(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxch4))
            plt(k) = plt(k) + 4.*calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxh2))
            plt(k) = plt(k) + 2.*calcmean(glb,imx*jmx,0,1.e-20,cpspval)
c
            glb(:,:) = fields(:,k,:,ifget(itxo1))
            plt0(k) = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2))
            plt0(k) = plt0(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxn2))
            plt0(k) = plt0(k) + calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            plt(k) = plt(k) / plt0(k)
          enddo
        elseif (ip.eq.ixe5577) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxt))
            gtn = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2))
            go2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo1))
            go1 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxn2))
            gn2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            call mke5577(gtn,go2,go1,gn2,1,plt(k),1,cpspval)
          enddo
        elseif (ip.eq.ixe6300) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxt))
            gtn = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxte))
            gte = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2))
            go2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo1))
            go1 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxn2))
            gn2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2p))
            go2p = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxne))
            gne = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            call mke6300(gtn,gte,go2,go1,gn2,go2p,gne,1,plt(k),1,
     +        cpspval)
          enddo
        elseif (ip.eq.ixeo200) then
          do k=1,kmx
            glb(:,:) = fields(:,k,:,ifget(itxt))
            gtn = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo2))
            go2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxo1))
            go1 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            glb(:,:) = fields(:,k,:,ifget(itxn2))
            gn2 = calcmean(glb,imx*jmx,0,1.e-20,cpspval)
            call mkeo200(gtn,go2,go1,gn2,1,plt(k),1,cpspval)
          enddo
        else
          write(6,"('>>> getloc: unknown ip=',i3,' returning zeroes')") 
     +      ip
          plt(:) = 0.
          return
        endif
      else
        write(6,"('getloc: bad ixlat,ixlon=',2i5)") ixlat,ixlon  
      endif
      return
      end
c
c--------------------------------------------------------------------
c
      function calcmean(rval,nval,log10,rmin,spval)
c
c Return mean of rval(nval), ignoring any values = spval. 
c If log10 flag is set, return log10 of mean, or 
c   spval if rval(i) is < rmin. 
c
      dimension rval(nval)
c
c     write(6,"('in calcmean: nval=',i3,' log10=',i2,' rmin=',
c    +  e12.4,' spval=',e12.4,' rval=',/(6e12.4))") 
c    +  nval,log10,rmin,spval,(rval(i),i=1,nval)
      calcmean = 0.
      ngood = 0
c
c Find the mean, checking for spval:
      do i=1,nval
        if (rval(i).ne.spval) then
          calcmean = calcmean + rval(i)
          ngood = ngood+1
        endif
      enddo
      if (ngood.eq.nval) then				! all good values
        calcmean = calcmean / float(nval)
      elseif (ngood.lt.nval.and.ngood.gt.0) then	! some spval
        calcmean = calcmean / float(ngood)
      else						! all spval
        calcmean = spval
      endif
      if (ngood.ne.nval) then
        write(6,"('Note: calcmean encountered spval: nval=',i3,
     +    ' ngood=',i3,' calcmean = ',e12.4)") nval,ngood,calcmean
      endif
c
c Return mean, with log10 taken if necessary:
      if (log10.gt.0) then
        if (calcmean.eq.spval.or.calcmean.le.rmin) then
          calcmean = spval
        else
          calcmean = alog10(calcmean)
        endif
      endif
      return
      end
