#include "dims.h"
!
      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,modeltime
      use init_module,only: istep
      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,iseries_p,iseries_s
      real :: fmin,fmax
      character(len=len_name) :: fprog_names(nfprog),
     |                           fsech_names(mxfsech)
      logical :: time2write,wrprim,wrsech,save_prim,save_sech,
     |  newseries_prim,newseries_sech
!
! External:
      integer,external :: strloc
      logical,external :: wrhist
!
! Return silently if not writing secondary histories:
      if (isechist==0) return
!
! Save only if we are writing a secondary history this step:
      time2write = wrhist(istep,modeltime,
     |  wrprim, save_prim, newseries_prim, iseries_p,
     |  wrsech, save_sech, newseries_sech, iseries_s) 
      if (.not.wrsech) return
!
!     write(6,"('addfsech: lat=',i2,' istep=',i6,' modeltime=',4i4,
!    |  ' name=',a)") lat,istep,modeltime,name
!
! Check that field is not a prognostic:
      fprog_names(:) = fprog(:)%name

!     write(6,"('addfsech lat=',i2,' nfprog=',i2,' name=',a,
!    |  ' loc(fprog_names)=',i8,' len(fprog_names)=',i2,
!    |  ' loc(name)=',i8)") lat,nfprog,name,loc(fprog_names),
!    |  len(fprog_names),loc(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 
      ix = strloc(fsech_names,mxfsech,name) 
      if (ix == 0) then ! was not in SECFLDS (not geographic)
!
! Check magnetic fields:
        fsech_names(:) = fsechmag(:)%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
!       write(6,"('Calling addfsech_geo for field ',a,' at lat=',i3)")
!    |    name,lat
        call addfsech_geo(name,long_name,units,f2d,ix,londim,levdim,
     |    nlev,lat)
      else
!       write(6,"('Calling addfsech_mag for field ',a,' at lat=',i3)")
!    |    name,lat
        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 field ',a,
!    |    ' at lat ',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 field ',a,
!    |    ' at lat ',i3,' fmin,max=',2e12.4)") 
!    |    trim(fsechmag(ix)%name),lat,fmin,fmax
!     endif
!
      end subroutine addfsech_mag