!
      module get_flds
      use get_hist,only: gethist
      use fields,only: field,flds_known,mxfknown
      use mk_derived
      use den_convert
      use epflux
!
! qbary module not available on IBM (because fitpack is not available):
#ifndef AIX
      use qbary
#endif
      implicit none
      contains
!-------------------------------------------------------------------
      subroutine getflds(vols,nvols,mxdiskvols,ivol,luin,mtime,tmpdir,
     +  f,nf,h,modelhts,iden,ionvel,iscntr,iprhist)
      use proc_lat
      use proc_lat_old
      use proc,only: gcmlat
      use hist,only: history
      use netcdf_module,only: f4d
!
! Get history at model time mtime from vols(nvols) (see module get_hist)
! If history is found, define requested fields f(nf).
!
! Args:
      character(len=*), intent(in) :: vols(nvols),tmpdir
      integer,intent(in) :: 
     +  nvols,		! number of history volumes
     +  mxdiskvols,	! max number of hist vols allowed on the disk
     +  nf,		! number of fields requested in f
     |  iscntr          ! iscntr == 1 if this is a diffs control case
      integer,intent(inout) :: 
     +  mtime(3),	! model time (day,hr,min)
     +  luin		! logical unit number to use
!
! On input, start search at file vols(ivol). If ivol<0, start at
!   vols(-ivol), but force file open (see gethist).
! On output, return ivol as index to vols(ivol) where history was
!   found. If i/o error or not found, return ivol=0 (see gethist).
!
      integer,intent(inout) :: ivol,ionvel
      type(field),intent(inout) :: f(nf)
      type(history),intent(out) :: h
      integer, intent(in) :: modelhts,iden,iprhist
!
! Locals:
      integer :: k ! temporary for wrphi
      integer :: i,ix,j,n,ier,ixz,ixui,ixvi,ixwi,ixphi,ixuivi,nlu
      integer :: ixoh,ixh,ixho2,ixt,ixu,ixv,ixw
      real :: rj,fmin,fmax
      real,allocatable :: flatslice(:,:,:)	! fields in lat slice
      character(len=8),allocatable :: fproclat(:) ! list of fields to be
						  ! processed by proclat
!
      logical :: isvtgcm
      interface
        function vecsum(u,v,id1,id2,spv)
          implicit none
          integer,intent(in) :: id1,id2
          real,intent(in) :: u(id1,id2),v(id1,id2) 
          real:: vecsum(id1,id2)			! result variable
          real,intent(in),optional :: spv
        end function vecsum
      end interface
!
! Get history (gethist will position file before lat slices):
!
      call gethist(h,vols,nvols,mxdiskvols,ivol,luin,mtime,tmpdir,
     |  iprhist)
      if (ivol <= 0) then
        write(6,"('getflds: history ',i3,':',i2,':',i2,
     +    ' not found')") mtime
        return
      endif
!
! Set up for a mars tgcm history:
!
      if (h%ismtgcm) call setmtgcm(f,nf,h,iprhist)
!
! Set up for a venus tgcm history:
! (make vtgcm CO2 a derived field and rotate ut by 12 hrs)
! (vtgcm N2 remains a derived field)
!
      isvtgcm = .false.
      if (trim(h%version)=='VTGCM'.or.h%version=='vtgcm') then
        isvtgcm = .true.
        call setvtgcm(f,nf,h,iprhist)
      endif
!
! Set up for jupiter tgcm history:
!
      if (h%isjtgcm) call setjtgcm(f,nf,h,iprhist)
!
! Set vertical pressure levels of history fields:
! (vertical coords of height-independent and height-only fields were
!  set by fset_known)
!
      call setflev(f,nf,h)
!
! Allocate space for list of history field names to be processed
! by proclat:
      allocate(fproclat(h%nflds),stat=ier)
      if (ier /= 0) call allocerr(ier,"getflds allocate fproclat")
!
! Allocate f(ix)%data, where fields f(ix) are available.
! (if f(ix) is derived, allocate it if all f(ix)%fneed are available,
!  including nested derived fields to arbitrary depth)
! Names of fields that need to be processed by proclat are returned
!   by allocfdat in fproclat(h%nflds).
! Allocfdat returns number of fields for which data was allocated.
!   (these fields may or may not have been requested)
!
      n = allocfdat(f,nf,h,fproclat,iprhist)
!
! Report to stdout fields that will be processed from the history:
      if (iprhist > 0) then
        n = count(fproclat/=' ')
        if (n > 0) then
          write(6,"(/'The following ',i3,' fields will be processed',
     +      ' from the history:')") n
          write(6,"(8a8)") fproclat(1:n)
        endif
      endif
!
! Return if no requested fields were allocated:
      n = 0
      do j=1,nf
        if (f(j)%requested.and.associated(f(j)%data)) n= n+1
      enddo
      if (n <= 0) then
        write(6,"(/'>>> tgcmproc: no fields to process:')")
        write(6,"( '    All requested fields are either not ',
     +    'available on the history,',/,'    or fields required to',
     +    ' calculated them are not available.'/)")
        ivol = -1 
        return
      endif
!
! Report to stdout the fields that were allocated:
      n = 0
      do j=1,nf
        if (associated(f(j)%data)) n = n+1
      enddo
      if (iprhist > 0) then
        write(6,"(/'The following ',i3,' field structures have been ',
     +    'allocated:')") n
        write(6,"('  (only requested fields will be plotted/output)')")
        write(6,"('Field     Requested   Known   Derived    Vtype',
     |    '    Default zptype')")
        write(6,"(72('-'))")
        do j=1,nf
          if (associated(f(j)%data))
     |      write(6,"(a8,2x,l1,11x,l1,7x,l1,10x,a,7x,a)") f(j)%fname8,
     |        f(j)%requested,any(flds_known%fname8==f(j)%fname8),
     |        f(j)%derived,trim(f(j)%vtype),trim(f(j)%zptype)
        enddo
      endif
!
! Turn off ionvel if necessary.
! 3/29/03 bf: May need to set ionvel = -1 here if ion vels have been 
!   requested but are not allocated or available.
! (recall that if ionvel==0, ion vels are read from the history,
!  otherwise they are derived, see fset_known.f)
!

!
! Allocate fields for latitude slices: 
!
      allocate(flatslice(h%nlon,h%nzp,h%nflds),stat=ier)
      if (ier /= 0) call allocerr(ier,"getflds allocate lat slice")
      if (iprhist > 0)
     +  write(6,"(/'Allocated history lat slice: nlon,nzp,nflds=',
     +    i3,',',i3,',',i3)") h%nlon,h%nzp,h%nflds
!
! Read and process latitude slices:
! (transfer of h%lu to n for read statement is temporary for sgi 
!  compiler. The patch for this should be available soon (12/97))
! (note h%nlat==18 for vtgcm)
!
      if (h%lu > 0 .and. h%ncid > 0) then
        write(6,"(/,'>>> WARNING getflds: both h%lu and h%ncid',
     |    ' > 0: h%lu=',i3,' h%ncid=',i3)") h%lu,h%ncid
      endif
      if (h%lu <=0 .and. h%ncid <= 0) then
        write(6,"(/,'>>> WARNING getflds: both h%lu and h%ncid',
     |    ' <= 0: h%lu=',i3,' h%ncid=',i3)") h%lu,h%ncid
      endif
      nlu = h%lu
      jloop: do j=1,h%nlat
!
! Transfer data read from netcdf file by nc_rdhist (nchist_mod.F): 
! flatslice(h%nlon,h%nzp,h%nflds) <= f4d(h%nlon-2,h%nlat,h%nzp,h%nflds)
!
        do n=1,h%nflds
          do i=1,h%nlon-2
            flatslice(i,:,n) = f4d(i,j,:,n)
          enddo
          flatslice(h%nlon-1,:,n) = flatslice(1,:,n)
          flatslice(h%nlon,  :,n) = flatslice(2,:,n)
        enddo

!       write(6,"(' ')")
!       do n=1,h%nflds
!         if (n <= 16) then
!           call fminmax(flatslice(1,1,n),h%nlon*h%nzp,fmin,fmax)
!           write(6,"('getflds: j=',i3,' n=',i3,' h%nflds=',i3,
!    |        ' h%fnames(n)=',a,' fmin,max=',2e15.6)") j,n,h%nflds,
!    |        h%fnames(n),fmin,fmax
!         endif
!       enddo
!
! Convert flatslice from midpoint to interface levels, change some units,
! and calculate heights if required. Proclat changes flatslice:
!
        if (h%isnew) then
          call proclat(flatslice(1:h%nlon-1,:,:),fproclat,f,nf,h,
     +      h%nlon-1,h%nzp,j,modelhts,0)
        else
          call proclat_old(flatslice(1:h%nlon-1,:,:),fproclat,h,
     +      h%nlon-1,h%nzp,j,modelhts,0)
        endif
!
! Allocate and calculate barm and pkt for density conversions 
! (see denconv.f) (note h%nzp==35 for vtgcm)
!
        call mkdenparms(flatslice(1:h%nlon-1,:,:),h,h%nlon-1,h%nzp,j)
!
! Define underived (history) fields needed in flds from flatslice:
! (just assign f(ix)%data(:,lat,:) <- flatslice(:,:,ixh))
! Then convert underived density species according to iden.
!
        call mkunderived(f,nf,flatslice(1:h%nlon-1,:,:),h,j)
        if (iden > 0) then
          call denconvert(f,nf,flatslice(1:h%nlon-1,:,:),h,h%nlon-1,
     +      h%nzp,j,iden)
        endif
!
! Make derived fields from flatslice:
! subroutine mkderived(f,nf,flat,h,imx,kmx,lat,iden,ionvel)
!   real,intent(inout) :: flat(imx,kmx,h%nflds)
!
        call mkderived(f,nf,flatslice(1:h%nlon-1,:,:),h,
     +    h%nlon-1,h%nzp,j,iden,ionvel)
!
! Change zptypes if requested (f%zptype_req was set from namelist grid_levels
! by sub setflev, called above):
!
        call set_zptype(f,nf,h,j)
      enddo jloop ! j=1,h%nlat
!
! 1/21/08 btf: Save global geopotential at midpoints and interfaces:
      call mkgeopot(f,nf,iscntr)
!
! Deallocate f4d:
      if (allocated(f4d)) deallocate(f4d)
!
! Set fields for vtgcm (n2 from external file and mirrored hemispheres):
! 9/06 btf: vtgcm histories now have global latitudes, so no need to
!           mirror hemispheres.
!     if (h%isvtgcm) call setvtgcmf(f,nf,h)
!
! Finish ion drifts (necessary 3d components were saved by savephi,
! called from mkderived in above lat loop)
!
      if (ionvel > 0) then
        ixui   = ixfindc(f%fname8,nf,'UI      ')
        ixvi   = ixfindc(f%fname8,nf,'VI      ')
        ixwi   = ixfindc(f%fname8,nf,'WI      ')
        ixuivi = ixfindc(f%fname8,nf,'UIVI    ')
        ixphi  = ixfindc(h%fnames,h%nflds,'POTEN   ') ! not used
        ixz    = ixfindc(h%fnames,h%nflds,'Z       ')
        do j=1,h%nlat
          if (f(ixui)%derived.or.f(ixvi)%derived.or.
     |        f(ixwi)%derived) then
          call mkdrifts(
     |      f(ixui)%data(:,j,:),f(ixui)%derived,
     |      f(ixvi)%data(:,j,:),f(ixvi)%derived,
     |      f(ixwi)%data(:,j,:),f(ixwi)%derived,
     |      h,h%nlon-1,h%nzp,h%nlat,j,ionvel)
!
          endif
! Calculate vector sum for UIVI:
          if (ixuivi>0) f(ixuivi)%data(:,j,:) = 
     +      vecsum(f(ixui)%data(:,j,:),f(ixvi)%data(:,j,:),
     +        size(f(ixui)%data,1),size(f(ixui)%data,3))
        enddo
      endif
!
! Calculate Eliasson-Palm fluxes:
! (using 3d quantities saved by save_epv, called from mkderived) 
!
      if (mkepv) call getepv(f,nf,h)
!
! Calculate qbary:
#if defined(AIX) || defined(SUN) || defined(LINUX)
#else
      if (mkqbary) call calc_qbary(f,nf,h,h%nlon-1,h%nzp,h%nlat)
#endif
!
! Release local space:
      deallocate(flatslice)
      deallocate(fproclat)
      deallocate(barm)
      deallocate(pkt)
      return
      end subroutine getflds
!-------------------------------------------------------------------
      integer function allocfdat(f,nf,h,fproc,iprhist)
!
! Allocate data components of requested field structures f(nf).
! Allocate only if either the field is available on the history,
!   or it is a derived field and all of its dependencies are
!   available on the history.
! Provide list of field names that will be required from the history
!   in fproc(h%nflds) (i.e, fproc/=' ').
! On input:
!   f(nf) = requested field structures
!   h%fnames(h%nflds) = fields available on the current history
! On output:
!   fproc(h%nflds) is list of field names that need to be processed
!     from the history (by proclat).
!   Return value allocfdat is number of fields in f for which the
!     data component was allocated. 
!
! Following this routine field loops can check that both f(i)%requested
!   and associated(f(i)%data) are true before processing/plotting, etc.
!
      use proc
      use flist
      implicit none 
!
! Args:
      integer,intent(in) :: nf,iprhist
      type(field),intent(inout) :: f(nf) ! inout so can change O+ to OP
      type(history),intent(in) :: h
      character(len=*),intent(out) :: fproc(h%nflds)
!
! Locals:
      integer :: i,ix,nfdep,nfproc,ier
      logical avail
      character(len=len(fproc)) :: fdep(mxfproc)
! 
      fproc = '        '	! init
      allocfdat = 0
      nfproc = 0
!
! Loop over requested fields:
      floop: do ix=1,nf
!
! 8/3/07 btf: Use nullify func to define the pointer (as empty). 
! If this is not done, linux pgf90 built code will fail in a 
! deallocate after testing w/ the associated function because
! its results are undefined if given an undefined pointer. 
!
        nullify(f(ix)%data)
        avail = any(h%fnames==f(ix)%fname8)
        if (.not.avail.and.(f(ix)%fname8)=='O+      ') then
          avail = any(h%fnames=='OP      ')
          if (avail.and.iprhist > 0) then
            i = ixfindc(h%fnames,h%nflds,'OP      ')
            f(ix)%fname8 = 'OP      '
            write(6,"('Note: will use OP field ',i2,' instead of O+')")i
          endif
        endif
        if (.not.avail.and.(f(ix)%fname8)=='O2+     ') then
          avail = any(h%fnames=='O2P     ')
          if (avail.and.iprhist > 0) then
            i = ixfindc(h%fnames,h%nflds,'O2P     ')
            f(ix)%fname8 = 'O2P     '
            write(6,"('Note: will use O2P field ',i2,
     |        ' instead of O2+')") i
          endif
        endif
!
! Field is on the history:
        if (avail) then	
!         write(6,"('allocfdat: Field ',a,' is available on the ',
!    |      'history: ix=',i3,' derived=',l1)") f(ix)%fname8,ix,
!    |      f(ix)%derived
!
! If data already associated, deallocate before allocating.
! IBM crashes after several histories if this check is not made.
! (was ok on SGI).
!
! 12/5/05 btf: hao linux (crumb) stops here on first unknown field:
!              "0: DEALLOCATE: memory at 31563248 not allocated"
!              Even tho associated(f(ix)%data) returns true.
!              
          if (associated(f(ix)%data)) then

!           write(6,"('allocfdat deallocating history field ',a,' ix=',
!    |        i3,' associated(f(ix)%data=',l1)") f(ix)%fname8,ix,
!    |        associated(f(ix)%data)

            deallocate(f(ix)%data)
          endif
          allocate(f(ix)%data(h%nlon-1,h%nlat,f(ix)%nlev),stat=ier)
          if (ier /= 0) then
            write(6,"(/,'>>> allocfdat: Error allocating history field '
     |        ,a,' h%nlon-1=',i3,' h%nlat=',i3,' f(ix)%nlev=',i3)") 
     |        f(ix)%fname8,h%nlon-1,h%nlat,f(ix)%nlev
            call allocerr(ier,"allocfdat allocating fields%data") 
          else
!           write(6,"('Allocated history field ',a,' h%nlon-1=',i3,
!    |        ' h%nlat=',i3,' ix=',i3,' f(ix)%nlev=',i3)") 
!    |        f(ix)%fname8,h%nlon-1,h%nlat,ix,f(ix)%nlev
          endif

          if (f(ix)%derived) then
            write(6,"('allocfdat: derived field ',a,' was found on',
     |        ' the history (ix=',i4,'), so setting it NOT derived')")
     |        f(ix)%fname8,ix
            f(ix)%derived = .false.
!
! If WN is on the history, it will be cm/s, not m/s (the units should
! be read from the history, but its inconvenient to do that here):
            if (f(ix)%fname8=='WN      ') f(ix)%units = 'cm/s'
          endif

          allocfdat = allocfdat+1

          if (.not.any(fproc==f(ix)%fname8)) then
            nfproc = nfproc+1
            fproc(nfproc) = f(ix)%fname8
          endif
!
! If history field is family species, process components as well
! (components will be needed for density conversion if iden>0)
          if (associated(f(ix)%fcomponents)) then
            do i=1,size(f(ix)%fcomponents)
              if (.not.any(fproc==f(ix)%fcomponents(i))) then
                nfproc = nfproc+1
                fproc(nfproc) = f(ix)%fcomponents(i)
              endif
            enddo
          endif
          cycle floop
!
! Field not on history and not derived:
        elseif (.not.f(ix)%derived) then
          if (iprhist > 0)
     +      write(6,"('WARNING: field ',a,' not available ',
     +        'on this history.')") f(ix)%fname8
!
! 8/3/07 btf: Use nullify func to define the pointer (as empty). 
! If this is not done, linux pgf90 built code will fail in a 
! deallocate after testing w/ the associated function because
! its results are undefined if given an undefined pointer. 
!
!         nullify(f(ix)%data)
!
! Derived field (not on history) -- check for dependencies:
! (getdep returns list of dependency names in fdep)
        else
          call getdep(f(ix),fdep,1,0)
          nfdep = shiftblanks(fdep,len(fdep),size(fdep)) ! number of dependencies

!         write(6,"('allocfdep after getdep: ix=',i3,' f(ix)%fname8=',a,
!    |      ' nfdep=',i3,' fdep(1:nfdep)=',/,(8a8))") ix,f(ix)%fname8,
!    |      nfdep,fdep(1:nfdep)
!
! Do not allocate data if any dependencies are not available:
!
          do i=1,nfdep
            if (.not.any(h%fnames==fdep(i))) then
!
! 11/97 Add: If checking dependencies for OH-V/B: if HO2 is the only
!   unavailable dependency for OH-V/B, allow HO2 to be calculated
!   w/ NO (if NO is available on the history and the other dependencies
!   (see mkho2 in util.f))
!
              if ((trim(f(ix)%type)=='OH-VIB'.or.
     +            trim(f(ix)%type)=='OH-BAND').and.trim(fdep(i))=='HO2') 
     +          then
                if (.not.any(h%fnames=='NO      ')) then
                  write(6,"('WARNING: ',a,' not allocated ',
     +              ' (HO2 and NO not available for OH-V/B fields)')")
     +              f(ix)%fname8
                  cycle floop  
                else
!                 write(6,"('Note: will use mkho2 with NO to calculate',
!    +              ' ',a,' because HO2 is not on the history.')")
!    +              f(ix)%fname8
                  fdep(i) = 'NO      '
                endif
!
! O2P is an alias for O2+ 
              elseif (trim(fdep(i))=='O2+') then 
                if (any(h%fnames=='O2P     ')) then
                  if (iprhist > 0)
     |              write(6,"('Note: using O2P alias for ',a)") fdep(i)
                  fdep(i) = 'O2P     '
                  cycle
                endif              
!
! OP is an alias for O+ 
              elseif (trim(fdep(i))=='O+') then 
                if (any(h%fnames=='OP      ')) then
                  if (iprhist > 0)
     |              write(6,"('Note: using OP alias for ',a)") fdep(i)
                  fdep(i) = 'OP      '
                  cycle
                endif              
              elseif (trim(fdep(i))=='W') then 
                if (any(h%fnames=='OMEGA   ')) then
                  write(6,"('Note: using OMEGA instead of W as ',
     |              'dependency for ',a)") f(ix)%fname8
                  fdep(i) = 'OMEGA   '
                  cycle
                else
                  write(6,"('>>> WARNING: cannot find W or OMEGA as ',
     |              'dependency for ',a)") f(ix)%fname8
                endif
              else
                write(6,"('WARNING: ',a,' not allocated ',
     |            '(dependency ',a,' not available)')") 
     |            f(ix)%fname8,fdep(i)
                cycle floop  
              endif
            endif
          enddo
!
! All dependencies available, so allocate the derived field:
          if (associated(f(ix)%data)) then
!           write(6,"('allocfdat deallocating derived field ',a)") 
!    |        f(ix)%fname8
            deallocate(f(ix)%data)
          endif
          allocate(f(ix)%data(h%nlon-1,h%nlat,f(ix)%nlev),stat=ier)
          if (ier /= 0) then
            write(6,"('>>> allocfdat: Error allocating derived ',
     |        'field ',a)") f(ix)%fname8
            call allocerr(ier,"allocfdat allocating fields%data") 
          endif
          allocfdat = allocfdat+1
!
! Add dependencies to fproc list:
          do i=1,nfdep
            if (.not.any(fproc==fdep(i))) then
              nfproc = nfproc+1
              fproc(nfproc) = fdep(i)
            endif
          enddo     
          nfproc = shiftblanks(fproc,len(fproc),size(fproc))
        endif	! available and/or derived
      enddo floop
!
! Report to stdout:
      if (allocfdat>0.and.iprhist>0) then
        write(6,"('Allocated f(ix)%data for ',i3,' fields:')") 
     +    allocfdat
        do ix=1,nf
          if (associated(f(ix)%data)) write(6,"(a)",advance='NO') 
     +      f(ix)%fname8
        enddo 
        write(6,"(' ')")
      endif
      end function allocfdat
!-------------------------------------------------------------------
      subroutine mkunderived(f,nf,flat,h,lat)
!
! Define underived (history) fields in f(nf) from flat
! (flat has been converted to intervace levels, with non-density unit
!  changes).
!
! Args:
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
      integer,intent(in) :: nf,lat
      real,intent(in) :: flat(h%nlon-1,h%nzp,h%nflds) 
!
! Locals:
      integer ix,ixh,k
      real :: fmin,fmax
!
      fields_loop: do ix=1,nf
!       ixh = ixfindc(h%fnames,h%nflds,f(ix)%fname8)
!       if (f(ix)%derived.or..not.associated(f(ix)%data).or.
!    |    ixh==0) cycle fields_loop
        if (f(ix)%derived.or..not.associated(f(ix)%data)) 
     +    cycle fields_loop
        ixh = ixfindc(h%fnames,h%nflds,f(ix)%fname8)
        if (ixh==0.and.f(ix)%fname8=='O+      ')
     |    ixh = ixfindc(h%fnames,h%nflds,'OP      ')
        if (ixh==0.and.f(ix)%fname8=='O2+     ')
     |    ixh = ixfindc(h%fnames,h%nflds,'O2P     ')
        if (ixh == 0) then
          write(6,"('WARNING mkunderived: could not find ',
     +      'underived field ',a,' in h%fnames.')")
     +      f(ix)%fname8         
          cycle fields_loop
        endif
        f(ix)%data(:,lat,:) = flat(:,:,ixh)

!       write(6,"('mkunderived: field ',a)") h%fnames(ixh)

      enddo fields_loop
      end subroutine mkunderived
!-------------------------------------------------------------------
      subroutine setflev(f,nf,h)
      use input,only: v5d_zprange,cdf_zprange,grid_levels
      use mk_v5d,only: v5d
      use proc,only: spval,mxfproc,gcmlev,gcmilev
!
! Set vertical pressure coords of history fields (vert coords of 
! known non-pressure fields (ht-indep and height only fields) were 
! set by fset_known).
!
! 3/06 btf: Vertical components of fields read from the history
!   were set by sub nc_rdfld (nchist_mod.F).
!
! Args:
      integer,intent(in) :: nf
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
!
! Locals:
      integer :: i,ix,ixf,ier,k,istat
      real :: dlev
      real,allocatable :: zp(:)
      character(len=80) :: char80
      integer,external :: ixfind
!
! Vertical coordinate components of fields on the history (known and unknown)
! were set by sub nc_rdfld (nchist_mod.F) while reading the history. Here we 
! set vertical coord for fields requested but not on the history, i.e., derived 
! fields), according to the f%vtype set by fset_known.
!
      do ix=1,nf
        if (f(ix)%vtype=='ZP      '.and..not.associated(f(ix)%lev)) 
     |    then
!
! For old histories, all fields are on interfaces:
          if (.not.h%isnew) f(ix)%zptype = 'INTERFACES'
!
          f(ix)%nlev = h%nzp
          allocate(f(ix)%lev(h%nzp),stat=istat)
          if (istat /= 0) then
            write(char80,"('setflev: error allocating lev for field ',
     |        a,': h%nzp=',i3)") f(ix)%fname8,h%nzp
            call allocerr(ier,char80)
          endif
          if (trim(f(ix)%zptype)=='MIDPOINTS') then
            f(ix)%lev = h%lev
          else
            f(ix)%lev = h%ilev
          endif
          f(ix)%dlev = 
     |      (f(ix)%lev(h%nzp)-f(ix)%lev(1))/float(h%nzp-1)

!         write(6,"('setflev: ix=',i3,' derived field ',a,' zptype=',a,
!    |      ' nlev=',i3,' lev=',f7.2,' to ',f7.2,' by ',f7.2)")
!    |      ix,f(ix)%fname8,f(ix)%zptype,f(ix)%nlev,f(ix)%lev(1),
!    |      f(ix)%lev(f(ix)%nlev),f(ix)%dlev

        endif
      enddo
!
! Transfer user requested zptype from grid_levels(2,mxflds) to zptype_req:
! character(len=80)      :: grid_levels(2,mxfproc)=" "  ! midpoints or interfaces

      f%zptype_req = f%zptype ! default
!     write(6,"('setflev: grid_levels(1,1)=',a,' h%isnew=',l1)") 
!    |  grid_levels(1,1),h%isnew
!
! grid_levels(1,1) == INTERFACES -> all fields on interfaces:
      if (trim(grid_levels(1,1))=='INTERFACES') then    ! user requests all i-faces
        f%zptype_req = 'INTERFACES'
        write(6,"('setflev: set all zptype_req to INTERFACES.')")
!
! grid_levels(1,1) == MIDPOINTS -> all fields on midpoints:
      elseif (trim(grid_levels(1,1))=='MIDPOINTS') then ! user requests all m-points
        f%zptype_req = 'MIDPOINTS'
        write(6,"('setflev: set all zptype_req to MIDPOINTS.')")
!
! grid_levels(1,1) == HISTORY -> all fields as on history (T,U,V on midpoints)
      elseif (trim(grid_levels(1,1))=='HISTORY') then
        if (h%isnew) then
          do ix=1,nf
            ixf = ixfindc(flds_known%fname8,mxfknown,f(ix)%fname8)
            if (ixf > 0) then
              f(ix)%zptype_req = flds_known(ixf)%zptype
            endif
          enddo
        else ! fields on "old" histories are on interfaces
          f(:)%zptype_req = 'INTERFACES'
        endif
!       write(6,"('setflev: set zptype_req to flds_known%zptype')")
!
! grid_levels(1,1) == HIST+TUV -> as on the history, except T,U,V on interfaces
      elseif (trim(grid_levels(1,1))=='HIST+TUV') then
        if (h%isnew) then
          do ix=1,nf
            ixf = ixfindc(flds_known%fname8,mxfknown,f(ix)%fname8)
            if (ixf > 0) then
              f(ix)%zptype_req = flds_known(ixf)%zptype
              if (trim(f(ix)%fname8)=='TN'.or.
     |            trim(f(ix)%fname8)=='UN'.or.
     |            trim(f(ix)%fname8)=='VN') then
                f(ix)%zptype_req = 'INTERFACES'
              endif
            endif
          enddo
        else ! fields on "old" histories are on interfaces
          f(:)%zptype_req = 'INTERFACES'
        endif
!
! User specified individual fields and zptypes:
      else 
        do i=1,mxfproc
          if (len_trim(grid_levels(1,i))>0.and.
     |        len_trim(grid_levels(2,i))>0) then
            ixf = ixfindc(f%fname8,nf,grid_levels(1,i))
            if (ixf > 0) then
              f(ixf)%zptype_req = grid_levels(2,i)
!             write(6,"('setflev: ixf=',i3,' field ',a,' zptype=',a,
!    |          ' zptype_req=',a)") ixf,f(ixf)%fname8,
!    |          trim(f(ixf)%zptype),trim(f(ixf)%zptype_req)
            endif ! found grid_levels(1,i) name in f
          endif
        enddo ! i=1,mxfproc
      endif
!
! Validate cdf_zprange if given: 
      if (any(cdf_zprange /= spval)) then
        if (cdf_zprange(1) < h%zpib) then
          write(6,"(/,'>>> Note setflev: cdf_zprange(1)=',e12.4,
     |      ' is below bottom interface boundary of model (',f8.2,')')")
     |      cdf_zprange(1),h%zpib
          write(6,"('Setting cdf_zprange(1) = ',e12.4,
     |      ' (bottom boundary)')") h%zpib
          cdf_zprange(1) = h%zpib
        endif
        if (cdf_zprange(2) > h%zpit) then
          write(6,"(/,'>>> Note setflev: cdf_zprange(2)=',e12.4,
     |      ' is above top interface boundary of model (',f8.2,')')")
     |      cdf_zprange(2),h%zpit
          write(6,"('Setting cdf_zprange(2) = ',e12.4,
     |      ' (top boundary)')") h%zpit
          cdf_zprange(2) = h%zpit
        endif
!
! Set cdf_zprange bot,top to nearest vertical coord:
        dlev = (h%zpit-h%zpib)/float(h%nzp-1)
        allocate(zp(h%nzp),stat=ier)
        if (ier /= 0) call allocerr(ier,"setflev allocate zp")
        do ix=1,h%nzp
          zp(ix) = h%zpib+(ix-1)*dlev
        enddo
        ix = ixfind(zp,h%nzp,cdf_zprange(1),dlev) 
        if (ix<=0) then 
          write(6,"('>>> setflev: error finding cdf_zprange(1)=',
     |      f8.2,' in zp=',/(9f8.2))") cdf_zprange(1),zp
        else
          cdf_zprange(1) = zp(ix)
        endif
        ix = ixfind(zp,h%nzp,cdf_zprange(2),dlev) 
        if (ix<=0) then 
          write(6,"('>>> setflev: error finding cdf_zprange(2)=',
     |      f8.2,' in zp=',/(9f8.2))") cdf_zprange(2),zp
        else
          cdf_zprange(2) = zp(ix)
        endif
        deallocate(zp)
!       write(6,"('setflev: cdf_zprange=',2f8.2)") cdf_zprange
      endif
!
! Zp range for v5d files is optionally user input:
! (if not provided by user, set to full vertical range of history)
! (v5d is vis5d type -- see module mk_v5d and input)
      if (len_trim(v5d%sendv5d)>0.or.len_trim(v5d%sendmsv5d)>0) then
        if (v5d_zprange(1)==spval.or.v5d_zprange(2)==spval) then
          v5d%zprange(1) = h%zpb
          v5d%zprange(2) = h%zpt
        else	! v5d_zprange was set -- verify its within range:
#ifdef SUN
          write(6,"('>>> setflev: v5d not available on the Sun.')")
          stop 'v5d on Sun'
#else
          if (v5d_zprange(1) < h%zpb) then
            write(6,"('>>> setflev: bad v5d_zprange bottom=',f8.2,
     +        ' will raise to h%zpb=',f8.2)") v5d_zprange(1),h%zpb
            v5d%zprange(1) = h%zpb
          endif
          if (v5d_zprange(2) > h%zpt) then
            write(6,"('>>> setflev: bad v5d_zprange top=',f8.2,
     +        ' will lower to h%zpt=',f8.2)") v5d_zprange(2),h%zpt
            v5d%zprange(2) = h%zpt
          endif
!
! Set v5d_zprange bot,top to nearest vertical coord:
          dlev = (h%zpt-h%zpb)/float(h%nzp-1)
          allocate(zp(h%nzp),stat=ier)
          if (ier /= 0) call allocerr(ier,"setflev allocate zp")
          do ix=1,h%nzp
            zp(ix) = h%zpb+(ix-1)*dlev
          enddo
          ix = ixfind(zp,h%nzp,v5d_zprange(1),dlev) 
          if (ix<=0) then 
            write(6,"('>>> setflev: error finding v5d_zprange(1)=',
     +        f8.2,' in zp=',/(9f8.2))") v5d_zprange(1),zp
          else
            v5d%zprange(1) = zp(ix)
          endif
          ix = ixfind(zp,h%nzp,v5d_zprange(2),dlev) 
          if (ix<=0) then 
            write(6,"('>>> setflev: error finding v5d_zprange(2)=',
     +        f8.2,' in zp=',/(9f8.2))") v5d_zprange(2),zp
          else
            v5d%zprange(2) = zp(ix)
          endif
          deallocate(zp)
#endif
        endif
!       write(6,"('setflev: v5d%zprange=',2f8.2)") v5d%zprange
      endif
      end subroutine setflev
!-------------------------------------------------------------------
      subroutine setjtgcm(f,nf,h,iprhist)
      use proc,only: p0
!
! Current history is from jtgcm (jupiter tgcm) -- make some necessary
!   changes:
! If JTGCM, then H2 is a derived field: h2 = 1.-h-he, so is 
!   dependent on h and he.
! For jtgcm, RHO=H+HE+H2
! Reset p0 for jupiter, and rotate ut by 12 hrs for jtgcm histories.
!
! Args:
      integer,intent(in) :: nf,iprhist
      type(field),intent(inout) :: f(nf)
      type(history),intent(inout) :: h
!
! Locals:
      integer :: n,j,ier
!
      if (iprhist > 0) then
        write(6,"(/72('-'))")
        write(6,"('This is a jtgcm (jupiter) history:')")
      endif
!
! Make jtgcm H2 a derived field:
!
      n = ixfindc(flds_known%fname8,mxfknown,'H2      ')
      allocate(flds_known(n)%fneed(2),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for jtgcm H2')
      write(flds_known(n)%fneed(1),"('H       ')") 
      write(flds_known(n)%fneed(2),"('HE      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'H2      ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     |      write(6,"('  Made H2 a derived field ',
     |        'for jtgcm, dependent on h and h2.')")
        endif
      endif
!
! For jtgcm, rho = h+he+h2 (where h2=1-h-he):
      n = ixfindc(flds_known%fname8,mxfknown,'RHO     ')
      write(flds_known(n)%fname56,"('TOTAL DENSITY (H+HE+H2)')")
      if (associated(flds_known(n)%fneed)) 
     +  deallocate(flds_known(n)%fneed)
      allocate(flds_known(n)%fneed(2),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for jtgcm rho')
      write(flds_known(n)%fneed(1),"('H       ')") 
      write(flds_known(n)%fneed(2),"('HE      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'RHO     ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     |      write(6,"('  Changed RHO dependencies to H, HE')")
        endif
      endif
!
! Set p0 for jtgcm histories:
      p0 = 4.5
      if (iprhist > 0) write(6,"('  Set p0=',e10.2,
     +  ' for jtgcm history.')") p0
!
! Do 12-hr rotation only if this is an old cray-blocked jtgcm file:
!   (jtgcm15 and later files are netcdf)
!
! 8/20/02: jtgcm15 histories do NOT have 12-hr offset.
!          See ~foster/tgcm/jtgcm15/cray2nc.
! 9/25/02: Do 12-hr rotation on old cray-blocked jtgcm histories:
!
! 11/28/05 btf: cannot read cray-blocked histories anymore.
!
      if (trim(h%version)=='JTGCM_CR') then
        write(6,"('>>> setjtgcm: cannot read JTGCM_CR cray-blocked',
     |    'histories.')")
        stop 'jtgcm'

!       if (iprhist > 0)
!    |    write(6,"('  Rotating old cray-blocked jtgcm history ',
!    |      'ut by 12 hours: mtime=',3i4)") h%mtime
!       h%ut = float(mod(h%mtime(2)+12,24))
!       if (iprhist > 0) write(6,"('  Rotated ut = ',f9.2)") h%ut
!       if (iprhist > 0) write(6,"(72('-')/)")
      endif
!
      end subroutine setjtgcm
!-------------------------------------------------------------------
      subroutine setmtgcm(f,nf,h,iprhist)
      use proc,only: p0
!
! Current history is from mtgcm (mars tgcm) -- make some necessary
!   changes:
! If MTGCM, then CO2 is a derived field: co2 = 1.-o-co-n2, so is
!   dependent on o, co, and n2.
!   (for tgcm, CO2 is a non-derived history field, see fset_known)
! N2 is written on mtgcm histories (is not derived, as in tgcm)
! For mtgcm, RHO=CO2+CO+N2+O
! Reset p0 for mars, and rotate ut by 12 hrs for mtgcm histories.
!
! Args:
      integer,intent(in) :: nf,iprhist
      type(field),intent(inout) :: f(nf)
      type(history),intent(inout) :: h
!
! Locals:
      integer :: n,j,ier
!
      if (iprhist > 0) then
        write(6,"(/72('-'))")
        write(6,"('This is an mtgcm (mars) history:')")
      endif
!
! Make mars CO2 a derived field:
!
      n = ixfindc(flds_known%fname8,mxfknown,'CO2     ')
      allocate(flds_known(n)%fneed(3),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for mtgcm co2')
      write(flds_known(n)%fneed(1),"('O1      ')") 
      write(flds_known(n)%fneed(2),"('N2      ')") 
      write(flds_known(n)%fneed(3),"('CO      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'CO2     ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     +      write(6,"('  Made CO2 a derived field ',
     +        'for mtgcm, dependent on o, co, and n2.')")
        endif
      endif
!
! Make mars N2 a non-derived field (is on mtgcm histories):
      n = ixfindc(flds_known%fname8,mxfknown,'N2      ')
      flds_known(n)%derived = .false. 
      if (associated(flds_known(n)%fneed)) 
     +  deallocate(flds_known(n)%fneed)
      j = n
      n = ixfindc(f%fname8,nf,'N2      ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     +      write(6,"('  Made N2 an underived field for mtgcm ',
     +        'history.')")
        endif
      endif
!
! For [mv]tgcm, rho = co2+n2+co+o:
      n = ixfindc(flds_known%fname8,mxfknown,'RHO     ')
      write(flds_known(n)%fname56,"('TOTAL DENSITY (CO2+CO+N2+O)')")
      if (associated(flds_known(n)%fneed)) 
     +  deallocate(flds_known(n)%fneed)
      allocate(flds_known(n)%fneed(4),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for mtgcm rho')
      write(flds_known(n)%fneed(1),"('CO2     ')") 
      write(flds_known(n)%fneed(2),"('CO      ')") 
      write(flds_known(n)%fneed(3),"('N2      ')") 
      write(flds_known(n)%fneed(4),"('O1      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'RHO     ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     +      write(6,"('  Changed RHO dependencies to CO2,CO,N2,O')")
        endif
      endif
!
! Set p0 for mtgcm histories:
      p0 = 1.2e-3
      if (iprhist > 0) write(6,"('  Set p0=',e10.2,
     +  ' for mtgcm history.')") p0
!
! mtgcm histories have 12-hr offset between mtime and ut:
!
! 2/15/07 btf: New CVS mtgcm no longer requires ut rotation.
!   (however, this will incorrectly not rotate old mtgcm16 histories)
! 3/6/07 btf: Rotate 12 hours only if mtgcm15 or mtgcm16 history
!   (mtgcm1.0 and later histories will not be rotated)
!
      write(6,"('  Mars model version = ',a)") trim(h%version)
      if (trim(h%version)=='mtgcm15'.or.trim(h%version)=='mtgcm16')
     |  then
        if (iprhist > 0) then
          write(6,"('  Rotating mtgcm ut by 12 hours: mtime=',
     |      3i4,', rotated ut=')",advance='no') h%mtime
        endif
        h%ut = float(mod(h%mtime(2)+12,24))
        if (iprhist > 0) write(6,"(f9.2)") h%ut
      else
        if (iprhist > 0) then
          write(6,"('  NOT rotating ut by 12 hours.')")
          write(6,"('  mtime=',3i4,' ut=',f9.2)") h%mtime,h%ut
        endif
      endif
      if (iprhist > 0) write(6,"(72('-')/)")
      end subroutine setmtgcm
!-------------------------------------------------------------------
      subroutine setvtgcm(f,nf,h,iprhist)
      use proc,only: p0
!
! Current history is from vtgcm (venus tgcm) -- make some necessary
!   changes:
! If VTGCM, then CO2 is a derived field: co2 = 1.-o-co-n2, so is
!   dependent on o, co, and n2.
!   (for tgcm, CO2 is a non-derived history field, see fset_known)
! Rotate ut by 12 hrs for mtgcm histories.
!
! Args:
      integer,intent(in) :: nf,iprhist
      type(field),intent(inout) :: f(nf)
      type(history),intent(inout) :: h
!
! Locals:
      integer :: n,j,ier
!
      if (iprhist > 0) then
        write(6,"(/72('-'))")
        write(6,"('This is an vtgcm (venus) history:')")
      endif
!
! Make venus CO2 a derived field:
!
      n = ixfindc(flds_known%fname8,mxfknown,'CO2     ')
      allocate(flds_known(n)%fneed(3),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for vtgcm co2')
      write(flds_known(n)%fneed(1),"('O1      ')") 
      write(flds_known(n)%fneed(2),"('N2      ')") 
      write(flds_known(n)%fneed(3),"('CO      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'CO2     ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     +      write(6,"('  Made CO2 a derived field ',
     +        'for vtgcm, dependent on o, co, and n2.')")
        endif
      endif
!
! 10/25/06 btf: If this is a new history, vtgcm N2 is not derived.
! (it should be on both secondary and primary new histories) 
! If this is an old history, will need to use vtgcm_n2, which was 
! read from external file, see sub getvtgcm in gethist.F.
! (see also mkderived.F)
!
      if (h%isnew) then
        n = ixfindc(flds_known%fname8,mxfknown,'N2      ')
        flds_known(n)%derived = .FALSE.
        n = ixfindc(f%fname8,nf,'N2      ')
        if (n > 0) then
          f(n)%derived = .FALSE.
          write(6,"('setvtgcm: h%isnew=',l1,' n=',i3,' field ',a,
     |      ': made underived..')") h%isnew,n,trim(f(n)%fname8)
        endif
      endif
!
! For [mv]tgcm, rho = co2+n2+co+o, so fix label:
      n = ixfindc(flds_known%fname8,mxfknown,'RHO     ')
      write(flds_known(n)%fname56,"('TOTAL DENSITY (CO2+CO+N2+O)')")
      if (associated(flds_known(n)%fneed)) 
     +  deallocate(flds_known(n)%fneed)
      allocate(flds_known(n)%fneed(4),stat=ier)
      if (ier /= 0) 
     +  call allocerr(ier,'allocating fneed for vtgcm rho')
      write(flds_known(n)%fneed(1),"('CO2     ')") 
      write(flds_known(n)%fneed(2),"('CO      ')") 
      write(flds_known(n)%fneed(3),"('N2      ')") 
      write(flds_known(n)%fneed(4),"('O1      ')") 
      flds_known(n)%derived = .TRUE.
      j = n
      n = ixfindc(f%fname8,nf,'RHO     ')
      if (n > 0) then
        ier = 0
        if (f(n)%requested) ier = 1
        f(n) = flds_known(j)	! struct assign
        if (ier==1) then
          f(n)%requested = .true.
          if (iprhist > 0)
     +      write(6,"('  Changed RHO dependencies to CO2,CO,N2,O')")
        endif
      endif
!
! Set p0 for vtgcm histories:
      p0 = 5.e-3
      if (iprhist > 0) write(6,"('  Set p0=',e10.2,
     +  ' for vtgcm history.')") p0
!
! vtgcm histories have 12-hr offset between mtime and ut:
!     if (iprhist > 0)
!    +  write(6,"('  Rotating vtgcm ut by 12 hours: mtime=',
!    +    3i4,', rotated ut=')",advance='no') h%mtime
!     h%ut = float(mod(h%mtime(2)+12,24))
!     if (iprhist > 0) write(6,"(f9.2)") h%ut
      if (iprhist > 0) write(6,"('NOT rotating vtgcm history: mtime=',
     |  3i4)") h%mtime
      if (iprhist > 0) write(6,"(72('-')/)")
      end subroutine setvtgcm
!-------------------------------------------------------------------
      subroutine setvtgcmf(f,nf,h)
!
! For vtgcm: h%nlat==18
! f%data has been defined by proclat for northern hem only at j=1->18,
!   so these need to be transferred to f%data(j=19->36) (north indices).
! These are mirrored across the equator for southern hem, with sign
!   change for vn, so:
!
! f north (j=19->36) gets f(j=1->18) from proclat, then
! f south (j=18->1 ) gets f(j=1->18), with sign change if vn
!
! Args:
      integer,intent(in) :: nf
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
!
! Locals:
      integer :: j,ixf
      real :: xfct
!
      do ixf = 1,nf
        xfct = 1.
        if (f(ixf)%fname8=='VN      ') xfct = -1.
!
! Transfer processed j slices to northern indices (19-36 <- 1-18)
        do j=1,h%nlat
          f(ixf)%data(:,j+h%nlat,:) = f(ixf)%data(:,j,:)
        enddo
!
! Mirror across equator, defining southern indices (18-1 <- 19-36)
        do j=1,h%nlat
          f(ixf)%data(:,h%nlat-j+1,:) = xfct*f(ixf)%data(:,j+h%nlat,:)
        enddo
      enddo
      end subroutine setvtgcmf
!-------------------------------------------------------------------
      subroutine getepv(f,nf,h)
!
! Calculate Eliasson-Palm fluxes:
! (necessary 3d and zm fields were saved by save_epv, called by 
!  mkderived in getflds lat loop)
!
! Args:
      integer,intent(in) :: nf
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
!
! Locals:
      integer :: j,ier,ixz,ixepvy,ixepvz,ixepvdiv,ixepvyz,ixepvyzmag
      real,allocatable :: fepvz(:,:)
      interface
        function vecsum(u,v,id1,id2,spv)
          implicit none
          integer,intent(in) :: id1,id2
          real,intent(in) :: u(id1,id2),v(id1,id2) 
          real:: vecsum(id1,id2)			! result variable
          real,intent(in),optional :: spv
        end function vecsum
      end interface
!
        ixepvy   = ixfindc(f%fname8,nf,'EPVY    ')
        ixepvz   = ixfindc(f%fname8,nf,'EPVZ    ')
        ixepvdiv = ixfindc(f%fname8,nf,'EPVDIV  ')
        ixepvyz  = ixfindc(f%fname8,nf,'EPVYZ   ')
        ixepvyzmag  = ixfindc(f%fname8,nf,'EPVYZMAG')
        if (ixepvy > 0) then
          if (.not.associated(f(ixepvy)%data)) ixepvy = 0
        endif
        if (ixepvz > 0) then
          if (.not.associated(f(ixepvz)%data)) ixepvz = 0
        endif
        if (ixepvdiv > 0) then
          if (.not.associated(f(ixepvdiv)%data)) ixepvdiv = 0
        endif
        if (ixepvyz > 0) then
          if (.not.associated(f(ixepvyz)%data)) ixepvyz = 0
          if (ixepvyz > 0.and.(ixepvy==0.or.ixepvz==0)) ixepvyz = 0
        endif
        if (ixepvyzmag > 0) then
          if (.not.associated(f(ixepvyzmag)%data)) ixepvyzmag = 0
          if (ixepvyzmag > 0.and.(ixepvy==0.or.ixepvz==0)) ixepvyzmag=0
        endif
        ixz = ixfindc(f%fname8,nf,'Z       ')
!
! epvy3d was saved by save_epv (is not dependent on adjacent lats):
!
        if (ixepvy > 0) then
          do j=1,h%nlat
            f(ixepvy)%data(:,j,:) = epvy3d(:,:,j)
          enddo
        endif
!
! epvdiv is dependent on epvz. fepvz stores epvz for epvdiv if epvz 
! was not requested: 
!
        if (ixepvdiv > 0 .and. ixepvz == 0) then
          allocate(fepvz(h%nlon-1,h%nzp),stat=ier)
          if (ier /= 0) 
     +      call allocerr(ier,"getflds allocating fepvz") 
        endif
!
! epvz and epvdiv are dependent on adjacent latitudes, so they are
! 1st determined from 2 to nlat-1, then boundaries are extrapolated 
! below.
!
        do j=2,h%nlat-1
!
! Calculate epvz:
!     subroutine epfluxz(t,u,v,w,fout,nlon,nlev,glat,lat)
!
          if (ixepvz > 0) then
            call epfluxz(t3d(:,:,j),u3d(:,:,j),v3d(:,:,j),
     +        w3d(:,:,j),f(ixepvz)%data(:,j,:),h%nlon-1,h%nzp,
     +        gcmlat(j),j)
          endif
!
! Calculate epvdiv (need epvz):
!     subroutine epfluxdiv(t,u,v,w,epvz,fout,nlon,nlev,glat,lat)
! (If epvz was requested, it was obtained above and is now used
!  in epfluxdiv call. If epvz was not requested, it is obtained
!  with call to epfluxz and stored in fepvz. Fepvz is then used
!  in the call to epfluxdiv)
!
          if (ixepvdiv > 0) then
            if (ixepvz > 0) then
              call epfluxdiv(t3d(:,:,j),u3d(:,:,j),v3d(:,:,j),
     +          w3d(:,:,j),f(ixepvz)%data(:,j,:),
     +          f(ixepvdiv)%data(:,j,:),h%nlon-1,h%nzp,gcmlat(j),j)
            else ! get epvz into fepvz, and pass it to epfluxdiv
              call epfluxz(t3d(:,:,j),u3d(:,:,j),v3d(:,:,j),
     +          w3d(:,:,j),fepvz,h%nlon-1,h%nzp,gcmlat(j),j)
              call epfluxdiv(t3d(:,:,j),u3d(:,:,j),v3d(:,:,j),
     +          w3d(:,:,j),fepvz,f(ixepvdiv)%data(:,j,:),h%nlon-1,
     +          h%nzp,gcmlat(j),j)
            endif
          endif
!
! Get vector magnitudes epvy+epvz:
! (EPVYZ    -> vector arrow plots of EPVY+EPVZ)
! (EPVYZMAG -> contour vector magnitudes EPVY+EPVZ)
!
          if (ixepvyz > 0) then
            f(ixepvyz)%data(:,j,:) = vecsum(f(ixepvy)%data(:,j,:),
     +        f(ixepvz)%data(:,j,:),size(f(ixepvy)%data,1),
     +        size(f(ixepvy)%data,3))
          endif
          if (ixepvyzmag > 0) then
            f(ixepvyzmag)%data(:,j,:) = vecsum(f(ixepvy)%data(:,j,:),
     +        f(ixepvz)%data(:,j,:),size(f(ixepvy)%data,1),
     +        size(f(ixepvy)%data,3))
          endif
        enddo		! j=2,h%nlat-1
!
! Extrapolate latitude boundaries of epvz and epvdiv (they are 
! dependent on adjacent lats, see subs epfluxz and epfluxdiv):
        if (ixepvz > 0) then
          f(ixepvz)%data(:,1,:) = 2.*f(ixepvz)%data(:,2,:)-
     +                               f(ixepvz)%data(:,3,:)
          f(ixepvz)%data(:,nlat,:) = 2.*f(ixepvz)%data(:,nlat-1,:)-
     +                               f(ixepvz)%data(:,nlat-2,:)
        endif 
        if (ixepvdiv > 0) then
          f(ixepvdiv)%data(:,1,:) = 2.*f(ixepvdiv)%data(:,2,:)-
     +                                 f(ixepvdiv)%data(:,3,:)
          f(ixepvdiv)%data(:,nlat,:) = 2.*f(ixepvdiv)%data(:,nlat-1,:)-
     +                                    f(ixepvdiv)%data(:,nlat-2,:)
        endif
        if (ixepvyz > 0) then ! vectors epvy+epvz (vectors plot)
          f(ixepvyz)%data(:,1,:) = vecsum(f(ixepvy)%data(:,1,:),
     +      f(ixepvz)%data(:,1,:),size(f(ixepvy)%data,1),
     +      size(f(ixepvy)%data,3))
          f(ixepvyz)%data(:,nlat,:) = vecsum(f(ixepvy)%data(:,nlat,:),
     +      f(ixepvz)%data(:,nlat,:),size(f(ixepvy)%data,1),
     +      size(f(ixepvy)%data,3))
        endif
        if (ixepvyzmag > 0) then ! vectors epvy+epvz (contour magnitudes)
          f(ixepvyzmag)%data(:,1,:) = vecsum(f(ixepvy)%data(:,1,:),
     +      f(ixepvz)%data(:,1,:),size(f(ixepvy)%data,1),
     +      size(f(ixepvy)%data,3))
          f(ixepvyzmag)%data(:,nlat,:)=vecsum(f(ixepvy)%data(:,nlat,:),
     +      f(ixepvz)%data(:,nlat,:),size(f(ixepvy)%data,1),
     +      size(f(ixepvy)%data,3))
        endif
        if (allocated(fepvz)) deallocate(fepvz) ! local
        call cleanepv	! release space in module epflux
      end subroutine getepv
!-------------------------------------------------------------------
      subroutine print_jslice(f,name,id1,id2,lat)
      implicit none
!
! Print values of 2d array to stdout (typically a j-slice at latitude 
! index lat):
!
! Args:
      integer,intent(in) :: id1,id2,lat 
      real,intent(in) :: f(id1,id2)
      character(len=*),intent(in) :: name 
!
! Locals:
      integer :: k
!
      write(6,"(/'Print 2d field ',a,' (id1=',i4,' id2=',i4,') lat=',
     +  i3)") name,id1,id2,lat
      do k=1,id2
        write(6,"('Field ',a,'(',i4,' ,1-',i4,')=')") name,k,id2
!       write(6,"(5e15.7)") f(k,:)
        write(6,"(4e23.15)") f(k,:)
      enddo
      end subroutine print_jslice
!-----------------------------------------------------------------------
      subroutine set_zptype(f,nf,h,lat)
      use hist,only: history	! use only history type-def
      use fields,only: field	! use only field type-def
      implicit none
!
! Args:
      integer,intent(in) :: nf,lat
      type(field),intent(inout) :: f(nf)
      type(history),intent(in) :: h
!
! Local:
      integer :: ix,k
!
! For old histories, all fields are plotted on interfaces:
      if (.not.h%isnew) then
        do ix=1,nf
          if (trim(f(ix)%vtype)/='HEIGHT') then
            f(ix)%zptype = 'INTERFACES'
            f(ix)%lev = h%ilev 
          endif
        enddo
      endif
!
! Switch zptype if requested: 
      floop: do ix=1,nf

!       if (lat==1)
!    |  write(6,"('set_zptype: lat=',i3,' ix=',i3,' field ',a,
!    |    ' associated(data)=',l1,' vtype=',a,' zptype=',a,
!    |    ' zptype_req=',a)") lat,ix,f(ix)%fname8,
!    |    associated(f(ix)%data),trim(f(ix)%vtype),trim(f(ix)%zptype),
!    |    trim(f(ix)%zptype_req)

        if (.not.associated(f(ix)%data)) cycle floop
        if (trim(f(ix)%vtype)=='HEIGHT') cycle floop
        if (len_trim(f(ix)%zptype_req)>0) then
          if (f(ix)%zptype(1:3) /= f(ix)%zptype_req(1:3)) then

            if (lat==1) 
     |      write(6,"('set_zptype: field ',a,': switching zptype ',
     |        'from ',a,' to ',a)") f(ix)%fname8,f(ix)%zptype,
     |        f(ix)%zptype_req
!
! Interpolate data:
            do k=f(ix)%nlev-1,2,-1
              f(ix)%data(:,lat,k) = 0.5*(f(ix)%data(:,lat,k)+
     |          f(ix)%data(:,lat,k-1))
            enddo
!
! Set structure components (last lat only):
            if (lat==h%nlat) then
              if (f(ix)%zptype_req(1:3)=='INT') then
                f(ix)%zptype = 'INTERFACES'
                f(ix)%lev = h%ilev 
!               write(6,"('set_zptype converted ',a,
!    |            ' from MIDPOINTS to INTERFACES')") f(ix)%fname8
              elseif (f(ix)%zptype_req(1:3)=='MID') then
                f(ix)%zptype = 'MIDPOINTS'
                f(ix)%lev = h%lev 
!               write(6,"('set_zptype converted ',a,
!    |            ' from INTERFACES to MIDPOINTS')") f(ix)%fname8
              else
                write(6,"('>>> set_zptype field ',a,': unknown ',
     |            'zptype_req=',a)") f(ix)%fname8,f(ix)%zptype_req
              endif
            endif ! lat==1
          endif ! zptype /= zptype_req
        endif ! zptype_req is non-blank

!       write(6,"('set_zptype: field ',a,' zptype=',a,' lev=',/,
!    |    (9f8.3))") f(ix)%fname8,f(ix)%zptype,f(ix)%lev

      enddo floop ! ix=1,nf
      end subroutine set_zptype
!-----------------------------------------------------------------------
      subroutine mkgeopot(f,nf,iscntr)
      use hist,only: history
      use fields,only: field,z_ifaces,z_midpts,zcntr_ifaces,zcntr_midpts
      use proc,only: spval
!
! Args:
      integer,intent(in) :: nf         ! number of field structs
      integer,intent(in) :: iscntr     ! iscntr==1 if control case
      type(field),intent(in) :: f(nf)  ! field structs
!
! Local:
      integer :: ixz,idim1,idim2,idim3,ier,k,nlev
      real,allocatable :: zifaces(:,:,:),zmidpts(:,:,:)
!
! External:
      integer,external :: ixfindc
!
      ixz = ixfindc(f%fname8,nf,'Z       ')
      if (ixz <= 0) then
        write(6,"('>>> mkgeopot: need geopotential in fields: ',
     |    ' field names = ',(8a8))") f%fname8
        call shutdown('geopotential')
      endif
!
! Get field dimensions (will be lon,lat,lev):
      idim1 = size(f(ixz)%data(:,1,1))
      idim2 = size(f(ixz)%data(1,:,1))
      idim3 = size(f(ixz)%data(1,1,:))
!
! Allocate local space:
      allocate(zifaces(idim1,idim2,idim3),stat=ier) 
      if (ier /= 0) 
     |  call allocerr(ier,"mkgeopot: error allocating zifaces")
      allocate(zmidpts(idim1,idim2,idim3),stat=ier) 
      if (ier /= 0) 
     |  call allocerr(ier,"mkgeopot: error allocating zmidpts")
      nlev = f(ixz)%nlev ! should be same as idim3
!
! Note that f(ixz) has already been processed by proclat, including
!   converted from cm to km, and possibly recalculated by sub calchts.
! Z is always calculated on interfaces by the model, but f(ixz) might
!   have been interpolated to midpoints (if user requested midpoints via
!   namelist read parameter grid_levels).
! (note midpoints are 1/2 dzp above interfaces)
!
! f(ixz) is at interfaces:
      if (f(ixz)%zptype == 'INTERFACES') then
        zifaces(:,:,:) = f(ixz)%data(:,:,:)
!
! Interpolate to midpoints (extrapolate to top level):
        do k=1,nlev-1
          zmidpts(:,:,k) = 0.5*(f(ixz)%data(:,:,k)+
     |      f(ixz)%data(:,:,k+1))
        enddo
        zmidpts(:,:,nlev) = 
     |    1.5*zmidpts(:,:,nlev-1)-0.5*zmidpts(:,:,nlev-2)
!
! f(ixz) is at midpoints:
      elseif (f(ixz)%zptype == 'MIDPOINTS') then
        zmidpts(:,:,:) = f(ixz)%data(:,:,:)
!
! Interpolate to interfaces (extrapolate to bottom level):
        do k=2,nlev
          zifaces(:,:,k) = 0.5*(f(ixz)%data(:,:,k)+
     |      f(ixz)%data(:,:,k-1))
        enddo
        zifaces(:,:,1) = 1.5*zifaces(:,:,2)-0.5*zifaces(:,:,3)
      else
        write(6,"('>>> mkgeopot: unknown zptype for Z: ',a)")
     |    f(ixz)%zptype
        call shutdown('mkgeopot')
      endif
!
! Allocate and assign to module data (fields.F) according to
! whether processing control or perturbed fields (if not a
! diffs run, then iscntr==0)
!
! iscntr == 0 -> is either not a diffs run, or perturbed fields of 
!   a diffs run are being processed:
!
      if (iscntr == 0) then 
        if (allocated(z_ifaces)) deallocate(z_ifaces)
        if (allocated(z_midpts)) deallocate(z_midpts)
        allocate(z_ifaces(idim1,idim2,idim3),stat=ier) 
        if (ier /= 0) 
     |    call allocerr(ier,"mkgeopot: error allocating z_ifaces")
        allocate(z_midpts(idim1,idim2,idim3),stat=ier) 
        if (ier /= 0) 
     |    call allocerr(ier,"mkgeopot: error allocating z_midpts")
        z_ifaces = zifaces
        z_midpts = zmidpts
!
! iscntr == 1 -> processing control fields of a diffs run:
!
      else          
        if (allocated(zcntr_ifaces)) deallocate(zcntr_ifaces)
        if (allocated(zcntr_midpts)) deallocate(zcntr_midpts)
        allocate(zcntr_ifaces(idim1,idim2,idim3),stat=ier) 
        if (ier /= 0) 
     |    call allocerr(ier,"mkgeopot: error allocating zcntr_ifaces")
        allocate(zcntr_midpts(idim1,idim2,idim3),stat=ier) 
        if (ier /= 0) 
     |    call allocerr(ier,"mkgeopot: error allocating zcntr_midpts")
        zcntr_ifaces = zifaces
        zcntr_midpts = zmidpts
      endif
!
! Release local memory:
      deallocate(zifaces)
      deallocate(zmidpts)
      end subroutine mkgeopot
!-----------------------------------------------------------------------
      end module get_flds