!
      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,nfprog,len_name
      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 :: i,ix,ixmag,k,nlev_used,istat
      real :: fmin,fmax
      character(len=len_name) :: fprog_names(nfprog),
     |                           fsech_names(mxfsech)
!
! External:
      integer,external :: strloc
!
! Return silently if not writing secondary histories:
      if (isechist==0) return
!
! Check that secondary histories are being written:
!
      i = 0
      do ix=1,mxfsech
        i = i+len_trim(fsech(ix)%name)
        i = i+len_trim(fsechmag(ix)%name)
      enddo
      if (i==0) then
        write(6,"('>>> addfsech: not writing secondary histories --',
     |    ' returning.')")
        return
      endif
!
! Check that field is not a prognostic:
      fprog_names(:) = fprog(:)%name

      write(6,"('addfsech lat=',i2,' nfprog =',i3,' loc(fprog_names)=',i8,
     |  ' len(fprog_names)=',i3,' name=',a,' loc(name)=',i8,
     |  ' len(name)=',i3)") lat,nfprog,loc(fprog_names),
     |  len(fprog_names),trim(name),loc(name),len(name)

      ix = strloc(fprog_names,nfprog,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:
      ixmag = 0
      fsech_names(:) = fsech(:)%name

      write(6,"('addfsech lat=',i2,' mxfsech=',i3,' loc(geo names)=',
     |  i8,' len(fsech_names)=',i3,' name=',a,' loc(name)=',i8,
     |  ' len(name)=',i3)") lat,mxfsech,loc(fsech_names),
     |  len(fsech_names),trim(name),loc(name),len(name)

      ix = strloc(fsech_names,mxfsech,name) 
      if (ix == 0) then ! was not in SECFLDS (not geographic)
!
! Check magnetic fields:
        fsech_names(:) = fsechmag(:)%name

      write(6,"('addfsech lat=',i2,' mxfsech=',i3,' loc(mag names)=',
     |  i8,' len(fsech_names)=',i3,' name=',a,' loc(name)=',i8,
     |  ' len(name)=',i3)") lat,mxfsech,loc(fsech_names),
     |  len(fsech_names),trim(name),loc(name),len(name)

        ixmag = strloc(fsech_names,mxfsech,name) 
        if (ixmag==0) then
          if (lat==2) 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
!
! fsech was allocated in set_fsech (flds_mod.f):
!   allocate(fsech(nfsech)%data(zimxp,zjmx,zkmxp),stat=ier)
! Double check that this pointer is ok:
!
      if (.not.associated(fsech(ix)%data)) then
        write(6,"('>>> addfsech_geo: fsech(ix)%data data not ',
     |    'associated. ix=',i3)") ix 
        stop 'addfsech_geo'
      endif
!
! Define fsech(ix)%data for current latitude:
! (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
!
!     if (lat==2) then
!       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
!     endif
!
      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
!
! fsechmag was allocated in set_fsechmag (flds_mod.f):
!   allocate(fsechmag(nfsech)%data(imaxmp,jmaxm,zkmxp+3),
! Double check that this pointer is ok:
!
      if (.not.associated(fsechmag(ix)%data)) then
        write(6,"('>>> addfsech_mag: fsechmag(ix)%data data not ',
     |    'associated. ix=',i3)") ix 
        stop 'addfsech_mag'
      endif
!
! Define fsechmag(ix)%data for current latitude:
! (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
!
!     if (lat==2) then
!       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
!     endif
!
      end subroutine addfsech_mag