!
      subroutine addfsech(name,long_name,units,f2d,londim,levdim,
     |  nlev,lat)
!
! Add a field to secondary histories.
! This routine is user-callable from inside the main latitude loop.
! It adds a 2d latitude slice of a diagnostic field to fsech(i)%data,
! which will be written to the secondary histories.
!
      use fields_module,only: fsech,fprog,fsechmag
      use hist_module,only: isechist
      implicit none
      include "params.h"
!
! Args:
      character(len=*),intent(in) :: name,long_name,units
      integer,intent(in) :: londim,levdim,nlev,lat
      real,intent(in) :: f2d(londim,levdim)
!
! Local:
      integer :: mxf,ix,ixmag,k,nlev_used
      real :: fmin,fmax
!
! External:
      integer,external :: strloc
!
! Return silently if not writing secondary histories:
      if (isechist==0) return
!
! Check that secondary histories are being written:
      if (all(len_trim(fsech%name)==0).and.
     |    all(len_trim(fsechmag%name)==0)) then
        write(6,"('>>> addfsech: not writing secondary histories --',
     |    ' returning.')")
        return
      endif
!
! Check that field is not a prognostic:
      mxf = size(fprog)
      ix = strloc(fprog%name,mxf,name) 
      if (ix > 0) then
        write(6,"('>>> addfsech: field ',a,' is a prognostic',
     |    ' -- returning.')") trim(name)
        return
      endif
!
! Check that field was requested on secondary histories:
! (i.e., was on either SECFLDS or SECFMAG inputs)
! Check geographic fields:
      mxf = size(fsech)
      ixmag = 0
      ix = strloc(fsech%name,mxf,name) 
      if (ix == 0) then ! was not in SECFLDS (not geographic)
!
! Check magnetic fields:
        mxf = size(fsechmag)
        ixmag = strloc(fsechmag%name,mxf,name) 
        if (ixmag==0) then
          if (lat==1) write(6,"('>>> addfsech: field ',a,' not found ',
     |      'in requested secondary history field names.')") trim(name)
          return
        endif
      endif
!
! Call for geographic or magnetic:
      if (ix > 0) then
        call addfsech_geo(name,long_name,units,f2d,ix,londim,levdim,
     |    nlev,lat)
      else
        call addfsech_mag(name,long_name,units,f2d,ixmag,londim,levdim,
     |    nlev,lat)
      endif
      return
      end subroutine addfsech
!-----------------------------------------------------------------------
      subroutine addfsech_geo(name,long_name,units,f2d,ix,londim,levdim,
     |    nlev,lat)
!
! Add a field on the geographic grid to secondary histories.
! This adds a 2d latitude slice of a diagnostic field to fsech(i)%data,
! which will be written to the secondary histories.
!
      use fields_module,only: fsech,fprog
      implicit none
      include "params.h"
!
! Args:
      character(len=*),intent(in) :: name,long_name,units
      integer,intent(in) :: londim,levdim,nlev,lat,ix
      real,intent(in) :: f2d(londim,levdim)
!
! Local:
      integer :: k,nlev_used
      real :: fmin,fmax
!
! Check dimensions:
      if (londim /= zimxp) then
        write(6,"('>>> addfsech_geo: londim must equal zimxp: londim=',i4,
     |    ' zimxp=',i4)") londim,zimxp
        stop 'addfsech_geo'
      endif
      if (levdim /= zkmxp) then
        write(6,"('>>> addfsech_geo: levdim must equal zkmxp: levdim=',i4,
     |    ' zkmxp=',i4)") levdim,zkmxp
        stop 'addfsech_geo'
      endif
      if (lat < 0 .or. lat > zjmx) then
        write(6,"('>>> addfsech_geo: bad lat=',i4,' (zjmx=',i3,')')")
     |    lat,zjmx
        stop 'addfsech_geo'
      endif
!
! nlev must be <= levdim:
      nlev_used = nlev
      if (nlev > levdim) then 
        write(6,"('>>> addfsech_geo: nlev must be <= levdim: nlev=',i4,
     |    ' levdim=',i4,' Will use nlev==levdim')") nlev,levdim
        nlev_used = levdim
      endif
!
! Define fsech(ix)%data for current latitude:
!           allocate(fsech(nfsech)%data(zimxp,zjmx,zkmxp),stat=ier)
! (if nlev < levdim, pad top levels nlev+1,levdim with spval.
!
      if (nlev_used == levdim) then
        fsech(ix)%data(:,lat,:) = f2d(:,:) 
      else ! nlev < levdim
        do k=1,nlev_used
          fsech(ix)%data(:,lat,k) = f2d(:,k)
        enddo
        do k=nlev_used+1,levdim
          fsech(ix)%data(:,lat,k) = spval
        enddo
      endif
!
      fsech(ix)%long_name = long_name
      fsech(ix)%units = units
!
!     call fminmax(f2d(1:zimxp,1:nlev),zimxp*nlev,fmin,fmax)
!     write(6,"('addfsech_geo: defined fsech diag field ',a,' at latitude ',
!    |  i3,' fmin,max=',2e12.4)") trim(fsech(ix)%name),lat,fmin,fmax
!
      end subroutine addfsech_geo
!-----------------------------------------------------------------------
      subroutine addfsech_mag(name,long_name,units,f2d,ix,londim,levdim,
     |    nlev,lat)
!
! Add a field on the magnetic grid field to secondary histories.
! This adds a 2d latitude slice of a diagnostic field to fsechmag(i)%data,
! which will be written to the secondary histories.
!
      use fields_module,only: fsechmag
      implicit none
      include "params.h"
!
! Args:
      character(len=*),intent(in) :: name,long_name,units
      integer,intent(in) :: londim,levdim,nlev,lat,ix
      real,intent(in) :: f2d(londim,levdim)
!
! Local:
      integer :: k,nlev_used
      real :: fmin,fmax
!
! Check dimensions:
      if (londim /= imaxmp) then
        write(6,"('>>> addfsech_mag: londim must equal imaxmp: ',
     |    'londim=',i4,' imaxmp=',i4)") londim,imaxmp
        stop 'addfsech_geo'
      endif
      if (levdim /= zkmxp+3) then
        write(6,"('>>> addfsech_mag: levdim must equal zkmxp+3: ',
     |    'levdim=',i4,' zkmxp+3=',i4)") levdim,zkmxp+3
        stop 'addfsech_mag'
      endif
      if (lat < 0 .or. lat > jmaxm) then
        write(6,"('>>> addfsech_mag: bad lat=',i4,' (jmaxm=',i3,')')")
     |    lat,jmaxm
        stop 'addfsech_mag'
      endif
!
! nlev must be <= levdim:
      nlev_used = nlev
      if (nlev > levdim) then 
        write(6,"('>>> addfsech_mag: nlev must be <= levdim: nlev=',i4,
     |    ' levdim=',i4,' Will use nlev==levdim')") nlev,levdim
        nlev_used = levdim
      endif
!
! Define fsechmag(ix)%data for current latitude:
!         allocate(fsechmag(nfsech)%data(imaxmp,jmaxm,zkmxp+3),
! (if nlev < levdim, pad top levels nlev+1,levdim with spval.
!
      if (nlev_used == levdim) then
        fsechmag(ix)%data(:,lat,:) = f2d(:,:) 
      else ! nlev < levdim
        do k=1,nlev_used
          fsechmag(ix)%data(:,lat,k) = f2d(:,k)
        enddo
        do k=nlev_used+1,levdim
          fsechmag(ix)%data(:,lat,k) = spval
        enddo
      endif
!
      fsechmag(ix)%long_name = long_name
      fsechmag(ix)%units = units
!
!     call fminmax(f2d(1:imaxmp,1:nlev),imaxmp*nlev,fmin,fmax)
!     write(6,"('addfsech_mag: defined fsechmag diag field ',a,
!    |  ' at latitude ',i3,' fmin,max=',2e12.4)") 
!    |  trim(fsechmag(ix)%name),lat,fmin,fmax
!
      end subroutine addfsech_mag