!
      subroutine minor(tn,o2,ox,w,difkk,fin,fnm,fhd,flbc,fubc,
     |  sloss,sprod,small,rmass,phix,alfax,xy,ilbc,iubc,fout,fout_nm,
     |  lev0,lev1,lon0,lon1,lat0,lat1,idebug)
!
! Advances minor species by one time step.
!
! Lower boundary:
!   If ilbc == 0, then
!     flbc(1) = A, flbc(2) = B, flbc(3) = C, where A*DPSX/DZ + B*PSX + C = 0.
!   If ilbc == 1,      (not the normal case)
!     flbc(3) = total upward number flux of F at lower boundary.
!
! Upper boundary:
!   If iubc == 0, then
!     fubc = diffusive upward number flux of F at upper boundary.
!   If iubc == 1, then
!     fubc = total upward number flux of F at upper boundary.
!   If iubc == 2, then
!     fubc = value of PSX at upper boundary.
!
! Sources:
!   sloss = Sx/n(x), where sx is portion of number density
!           source proportional to n(x), the minor species number
!           density     (level k+1/2)
!   sprod = Portion of number density source independent
!           of n(x).    (level k+1/2)
!
      use params_module,only: nlonp4,dz
      use lbc,only: b,fb
      use cons_module,only: rmassinv_o2,rmassinv_ox,rmassinv_n2,p0,
     |  expz,expzmid,expzmid_inv,boltz,rmass_o2,rmass_ox,rmass_n2,
     |  dtx2inv,grav,avo,dtsmooth,dtsmooth_div2,shapiro
      use input_module,only: difhor
      use fields_module,only: tlbc
      use addfld_module,only: addfld
#ifdef MPI
      use mpi_module,only: mp_bndlons_f3d, mp_periodic_f3d
#endif
      implicit none
!
! Input Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1,iubc,ilbc,
     |  idebug
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,o2,ox,w,difkk,fin,fnm,fhd
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(in)::
     |  sloss,  ! portion of number density source proportional to minor sp s1,s14
     |  sprod   ! portion of number density source independent of minor sp  s2,s15
      real,dimension(lon0:lon1),intent(in) ::
     |  flbc(lon0:lon1,3,lat0:lat1), ! 3 components of lower boundary
     |  fubc(lon0:lon1,lat0:lat1)    ! upper boundary
      real,intent(in) :: rmass,phix(3),alfax,small(lev1-lev0+1),xy
!
! Output Args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) :: fout,fout_nm
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1,lat,nlevs,lonbeg,lonend
      real :: phi(2,3),tau,t00,alfa12,alfa21,alfax1,alfax2
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  hadvec,    ! horizontal diffusion (output of routine advec)
     |  barm,      ! mbar (probably same as 4d barm, which could be passed in) 
     |  barmi,     ! mbar at interfaces ! s11
     |  s0mbar,    ! src0*mass/mbar     ! s15
     |  dmdz,      ! dm/dz              ! s10
     |  oxi,       ! o2 at interfaces   ! s9
     |  o2i,       ! ox at interfaces   ! s8
     |  oxdz,      ! dox/dz             ! s7
     |  o2dz,      ! do2/dz             ! s6
     |  tni,       ! tn at interfaces   ! s5
     |  phialfa_ox,                     ! s1
     |  alfa_o2,                        ! s2
     |  alfa_ox,                        ! s3
     |  phialfa_o2,                     ! s4
     |  thdiff,    ! thermal diffusion  ! s12
     |  ax,                             ! s11
     |  pcoef, qcoef, rcoef, fcoef      ! s1,2,3,4
      real,dimension(lon0:lon1) :: 
     |  lbc_o2, lbc_ox, lbcmbar, dfactor,
     |  flbc1,  flbc2,  flbc3,   fubc1
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) ::
     |  ftm1_smooth  ! time-smoothed field at time n-1 (s10)
!     real,dimension(lev0:lev1,lon0-2:lon1+2,lat0:lat1) ::
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2) ::
     |  ftm1_jdif ! lat diffs (includes lon boundaries) (s8)
      integer :: jm2,jm1,j0,jp1,jp2 ! lat-2,lat-1,lat,lat+1,lat+2
! 
      phi(:,1)=(/0.  ,0.673/)
      phi(:,2)=(/1.35,0.   /)
      phi(:,3)=(/1.11,0.769/)
      tau=1.86E+3
      t00=273.
!
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
      nlevs = nk
!
      lonbeg = lon0
      if (lon0==1) lonbeg = 3
      lonend = lon1
      if (lon1==nlonp4) lonend = lon1-2 ! nlonp4==nlonp4
!
! Check inputs:
!
      if (idebug > 0) then
!       write(6,"(/,'minor: iubc,ilbc=',2i3)") iubc,ilbc
!       write(6,"('minor: rmass=',e12.4)") rmass
!       write(6,"('minor: phix =',3e12.4)") phix
!       write(6,"('minor: alfax=',3e12.4)") alfax
!       write(6,"('minor: xy   =',3e12.4)") xy

        do lat=lat0,lat1
          call addfld('TN_MNR' ,' ',' ',tn (:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('O2_MNR' ,' ',' ',o2 (:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('OX_MNR' ,' ',' ',ox (:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('SPROD'   ,' ',' ',sprod(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('SLOSS'   ,' ',' ',sloss(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FIN_MNR',' ',' ',fin(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FNM_MNR',' ',' ',fnm(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('PDH_MNR',' ',' ',fhd(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
 
!         write(6,"('minor: lat=',i3,' flbc(1)=',/,(6e12.4))") 
!    |      lat,flbc(:,1,lat)
!         write(6,"('minor: lat=',i3,' flbc(2)' ,/,(6e12.4))") 
!    |      lat,flbc(:,2,lat)
!         write(6,"('minor: lat=',i3,' flbc(3)=',/,(6e12.4))") 
!    |      lat,flbc(:,3,lat)
!         write(6,"('minor: lat=',i3,' fubc ='  ,/,(6e12.4))") 
!    |      lat,fubc(:,lat)
!         write(6,"('minor: lat=',i3,' small='  ,/,(6e12.4))") 
!    |      lat,small

        enddo ! lat=lat0,lat1
      endif ! idebug > 0
!
! Latitude scan over task subdomain:
! real,dimension(lev0:lev1,lon0-2:lon1+2,lat0:lat1) ::
!   ftm1_jdif ! lat diffs (includes lon boundaries) (s8)

      ftm1_jdif = 0. ! init
      do lat=lat0,lat1
        jm2 = lat-2
        jm1 = lat-1
        j0  = lat
        jp1 = lat+1
        jp2 = lat+2
!
! Shapiro time smoother:
!
! real,dimension(lev0:lev1,lon0-2:lon1+2,lat0:lat1) ::
! ftm1_jdif = 4th order diffs in latitude at time n-1:
! f = f(j)-shapiro * ( f(j+2)+f(j-2) - 4.*(f(j+1)+f(j-1)) + 6.*f(j) )
!
        do i=lon0,lon1
          do k=lev0,lev1
            ftm1_jdif(k,i,lat) = fnm(k,i,j0) - shapiro *
     |           (fnm(k,i,jp2) + fnm(k,i,jm2)  -
     |        4.*(fnm(k,i,jp1) + fnm(k,i,jm1)) +
     |        6.*fnm(k,i,j0))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1
      enddo ! lat=lat0,lat1
!------------------------- End first latitude scan ---------------------
#ifdef MPI
!
! Get subdomain boundary longitudes for ftm1_jdif, to complete shapiro
! time smoother.
!
      call mp_periodic_f3d(ftm1_jdif(:,lon0:lon1,lat0:lat1),
     |  lev0,lev1,lon0,lon1,lat0,lat1,1)
      call mp_bndlons_f3d(ftm1_jdif,nlevs,lon0,lon1,lat0,lat1,1,0)
#endif
!----------------------- Begin second latitude scan --------------------
      do lat=lat0,lat1
!
! ftm1_smooth = zonally-smoothed field at time n-1:
! f = f(i)-shapiro * ( f(i+2)+f(i-2) - 4.*(f(i+1)+f(i-1)) + 6.*f(i) )
        do i=lonbeg,lonend
          do k=lev0,lev1-1
            ftm1_smooth(k,i,lat) = ftm1_jdif(k,i,lat) - shapiro *
     |           (ftm1_jdif(k,i+2,lat) + ftm1_jdif(k,i-2,lat)  -
     |        4.*(ftm1_jdif(k,i+1,lat) + ftm1_jdif(k,i-1,lat)) +
     |        6.* ftm1_jdif(k,i,lat))
          enddo ! k=lev0,lev1-1
        enddo ! i=lonbeg,lonend
!
! Set periodic points to zero:
        if (lon0==1) ftm1_smooth(:,lon0:lon0+1,lat) = 0.
        if (lon1==nlonp4) ftm1_smooth(:,lon1-1:lon1,lat) = 0.
!
! Advec input:
! real,dimension(lev0:lev1,lon0-2:lon1+2,lat-2:lat+2),intent(in) :: f
!
        call advec(fin(:,:,lat-2:lat+2),hadvec,lev0,lev1,lon0,lon1,lat)
        if (idebug > 0) 
     |    call addfld('HADVEC',' ',' ',hadvec(lev0:lev1-1,:),
     |      'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! Lower boundary:
        do i=lon0,lon1
          lbc_o2(i) = b(i,1,1)*o2(lev0,i,lat)+                    ! t4
     |                b(i,1,2)*ox(lev0,i,lat)+fb(i,1,lat)
          lbc_ox(i) = b(i,2,1)*o2(lev0,i,lat)+                    ! t5
     |                b(i,2,2)*ox(lev0,i,lat)+fb(i,2,lat)
          lbcmbar(i) = 1./(lbc_o2(i)*rmassinv_o2+                 ! t6
     |                     lbc_ox(i)*rmassinv_ox+
     |                     (1.-lbc_o2(i)-lbc_ox(i))*rmassinv_n2)
        enddo ! i=lon0,lon1
!
! Set mean molecular weight (could probably be passed in from 4d barm):
        do i=lon0,lon1
          do k=lev0,lev1
            barm(k,i) = 1./(o2(k,i,lat)*rmassinv_o2+         ! s12
     |                      ox(k,i,lat)*rmassinv_ox+
     |        (1.-o2(k,i,lat)-ox(k,i,lat))*rmassinv_n2)
            s0mbar(k,i) = sprod(k,i,lat)*rmass*boltz*tn(k,i,lat)/ ! s15
     |        (p0*expz(k)*barm(k,i))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1
!
! Lower boundary:
        do i=lon0,lon1
          barmi(lev0,i) = 0.5*(lbcmbar(i)+barm(lev0,i))     ! s11
          dmdz(lev0,i)  = (barm(lev0,i)-lbcmbar(i))/dz      ! s10
          oxi (lev0,i)  = 0.5*(lbc_ox(i)+ox(lev0,i,lat))    ! s9
          o2i (lev0,i)  = 0.5*(lbc_o2(i)+o2(lev0,i,lat))    ! s8
          oxdz(lev0,i)  = (ox(lev0,i,lat)-lbc_ox(i))/dz     ! s7
          o2dz(lev0,i)  = (o2(lev0,i,lat)-lbc_o2(i))/dz     ! s6
!
! Levels 2->top
          do k=lev0+1,lev1
            barmi(k,i) = 0.5*(barm(k,i)+barm(k-1,i))        ! s11
            dmdz(k,i) = (barm(k,i)-barm(k-1,i))/dz          ! s10
            oxi(k,i) = 0.5*(ox(k,i,lat)+ox(k-1,i,lat))      ! s9
            o2i(k,i) = 0.5*(o2(k,i,lat)+o2(k-1,i,lat))      ! s8
            oxdz(k,i) = (ox(k,i,lat)-ox(k-1,i,lat))/dz      ! s7
            o2dz(k,i) = (o2(k,i,lat)-o2(k-1,i,lat))/dz      ! s6
          enddo ! k=lev0+1,lev1
        enddo ! i=lon0,lon1

        if (idebug > 0) then
          call addfld('BARMX',' ',' ',barm ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('BARMI',' ',' ',barmi,
     |      'ilev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('DMDZ' ,' ',' ',dmdz ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('OXI'  ,' ',' ',oxi  ,
     |      'ilev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('O2I'  ,' ',' ',o2i  ,
     |      'ilev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('OXDZ' ,' ',' ',oxdz ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('O2DZ' ,' ',' ',o2dz ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('S0MBAR',' ',' ',s0mbar,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
        endif
!
! Tn at interfaces (s5):
! (Note bottom boundary of tn is stored in top slot)
! 5/4/06 btf: bottom boundary of TN is in tlbc:
        do i=lon0,lon1
!         tni(lev0,i) = tn(lev1  ,i,lat) ! bottom boundary in top slot (s5)
          tni(lev0,i) = tlbc(i,lat)
          tni(lev1,i) = tn(lev1-1,i,lat) ! top boundary same as lev1-1
          do k=lev0+1,lev1-1
            tni(k,i) = 0.5*(tn(k,i,lat)+tn(k-1,i,lat))
          enddo
        enddo

! Evaluate alfa12, alfa21, alfam1, alfam2
        alfa12=phi(1,2)-phi(1,3)
        alfa21=phi(2,1)-phi(2,3)
        alfax1=phix(1)-phix(3)
        alfax2=phix(2)-phix(3)
!
        do i=lon0,lon1
          do k=lev0,lev1
            phialfa_ox(k,i) = -(phi(1,3)+alfa12*oxi(k,i))                ! s1
            alfa_o2(k,i) = alfa12*o2i(k,i)                               ! s2
            alfa_ox(k,i) = alfa21*oxi(k,i)                               ! s3
            phialfa_o2(k,i) = -(phi(2,3)+alfa21*o2i(k,i))                ! s4
            thdiff(k,i) = ((alfax1*phialfa_o2(k,i)-alfax2*alfa_ox(k,i))* ! s12
     |        (o2dz(k,i)-(1.-(rmass_o2+dmdz(k,i))/barmi(k,i))*o2i(k,i))+ 
     |        (alfax2*phialfa_ox(k,i)-alfax1*alfa_o2(k,i))*
     |        (oxdz(k,i)-(1.-(rmass_ox+dmdz(k,i))/barmi(k,i))*oxi(k,i)))
     |        /(phialfa_ox(k,i)*phialfa_o2(k,i)-alfa_o2(k,i)*
     |        alfa_ox(k,i))+1.-(rmass+dmdz(k,i))/barmi(k,i)
            dmdz(k,i) = dmdz(k,i)/barmi(k,i)                             ! s10
            ax(k,i) = -barmi(k,i)/(tau*rmass_n2)*(t00/tni(k,i))**0.25/   ! s11
     |        (phix(3)+alfax1*o2i(k,i)+alfax2*oxi(k,i))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

        if (idebug > 0) then
          call addfld('PHIALFOX',' ',' ',phialfa_ox,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('ALFAO2'  ,' ',' ',alfa_o2   ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('ALFAOX'  ,' ',' ',alfa_ox   ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('PHIALFO2',' ',' ',phialfa_o2,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('THDIFF'  ,' ',' ',thdiff    ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('DMDZ'    ,' ',' ',dmdz      ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('AX'      ,' ',' ',ax        ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
        endif

        do i=lon0,lon1
          do k=lev0+1,lev1-1
            thdiff(k,i) = thdiff(k,i)-alfax*(tni(k+1,i)-tni(k-1,i))/     ! s12
     |        (2.*dz*tni(k,i))
          enddo ! k=lev0+1,lev1-1
!
! Lower boundary:
          thdiff(lev0,i) = thdiff(lev0,i)-alfax*
     |      (tni(lev0+1,i)-tni(lev0,i))/(dz*tni(lev0,i))
!
! Upper boundary:
          thdiff(lev1,i) = thdiff(lev1,i)-alfax*
     |      (tni(lev1,i)-tni(lev1-1,i))/(dz*tni(lev1,i))
        enddo ! i=lon0,lon1

        if (idebug > 0)
     |    call addfld('THDIFF1'  ,' ',' ',thdiff    ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)

        if (difhor == 1) then
          call dfact(dfactor,lon0,lon1,lat)
        else
          dfactor(:) = 1.
        endif
!
! Shapiro smoother: 

!
! Set coefficients for tridiagonal solver:
        do i=lon0,lon1
          do k=lev0,lev1-1
!
! Note difkk (mgw.F) will match tgcm24 only if both are set zero. 
! This is because mgw gravity wave parameterization is new in timegcm1.
!
            pcoef(k,i) = ax(k,i)/dz*(1./dz+0.5*thdiff(k,i))-expz(k)*     ! s1
     |        (expzmid_inv*difkk(k,i,lat)*dfactor(i)*(1./dz-0.5*
     |        dmdz(k,i))+0.25*(w(k,i,lat)+w(k+1,i,lat)))/dz
! 
            rcoef(k,i) = ax(k+1,i)/dz*(1./dz-0.5*thdiff(k+1,i))-expz(k)* ! s3
     |        (expzmid*difkk(k+1,i,lat)*dfactor(i)*(1./dz+0.5*
     |        dmdz(k+1,i))-0.25*(w(k,i,lat)+w(k+1,i,lat)))/dz
!
            qcoef(k,i) = -(ax(k,i)/dz*(1./dz-0.5*thdiff(k,i))+ax(k+1,i)/ ! s2
     |        dz*(1./dz+0.5*thdiff(k+1,i)))+expz(k)*((expzmid_inv*
     |        difkk(k,i,lat)*(1./dz+0.5*dmdz(k,i))+expzmid*
     |        difkk(k+1,i,lat)*(1./dz-0.5*dmdz(k+1,i)))*dfactor(i)/
     |        dz-sloss(k,i,lat)+dtx2inv)
!
            fcoef(k,i) = expz(k)*(ftm1_smooth(k,i,lat)*dtx2inv-          ! s4
     |        hadvec(k,i)+fhd(k,i,lat)+s0mbar(k,i))
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
        
        if (idebug > 0) then
!
! Small diamond diffs with tgcm24 at lbc of S0MBAR and FCOEF.
          call addfld('AX'     ,' ',' ',ax     ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('THDIFF' ,' ',' ',thdiff ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('DIFKK'  ,' ',' ',difkk(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('OMEGA'  ,' ',' ',w(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('SLOSS'  ,' ',' ',sloss(:,i0:i1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('DMDZ'   ,' ',' ',dmdz   ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('HADVEC' ,' ',' ',hadvec(lev0:lev1-1,:),
     |      'lev',lev0,lev1-1,'lon',i0,i1,lat)
          call addfld('S0MBAR' ,' ',' ',s0mbar ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FSMOOTH',' ',' ',
     |      ftm1_smooth(lev0:lev1-1,i0:i1,lat),
     |      'lev',lev0,lev1-1,'lon',i0,i1,lat)
          call addfld('PCOEF' ,' ',' ',pcoef ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('RCOEF' ,' ',' ',rcoef ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('QCOEF' ,' ',' ',qcoef ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FCOEF' ,' ',' ',fcoef ,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
        endif
!
! If ilbc==1, set flux boundary condition at bottom (normally, this 
!   is called with ilbc==0).
        if (ilbc==1) then ! (normally not true)
          do i=lon0,lon1
            flbc1(i) = 1.                                           ! t1
            flbc2(i) = dmdz(lev0,i)-w(lev0,i,lat)/difkk(lev0,i,lat) ! t2
            flbc3(i) = flbc(i,3,lat)*grav*rmass/(difkk(lev0,i,lat)* ! t3
     |        p0*expz(1)*expzmid_inv*avo) 
          enddo ! i=lon0,lon1
        else
!
! Use input flbc(lon0:lon1,3,lat0:lat1) ! 3 components of lower boundary
          do i=lon0,lon1
            flbc1(i) = flbc(i,1,lat)
            flbc2(i) = flbc(i,2,lat)
            flbc3(i) = flbc(i,3,lat)
          enddo ! i=lon0,lon1
        endif

!       if (idebug > 0) then
!         write(6,"('minor: ilbc=',i3,' lat=',i3,' flbc1=',/,(6e12.4))") 
!    |      ilbc,lat,flbc1
!         write(6,"('minor: ilbc=',i3,' lat=',i3,' flbc2=',/,(6e12.4))") 
!    |      ilbc,lat,flbc2
!         write(6,"('minor: ilbc=',i3,' lat=',i3,' flbc3=',/,(6e12.4))") 
!    |      ilbc,lat,flbc3
!       endif

!
! Modify ex if iubc==1:
        if (iubc==1) then
          do i=lon0,lon1
            fubc1(i) = expz(lev1-1)*expzmid*w(lev1,i,lat)/ax(lev1,i) ! t7
            thdiff(lev1,i) = thdiff(lev1,i)-fubc1(i)                 ! s12
          enddo !i=lon0,lon1
        endif
!
! Lower boundary:
        do i=lon0,lon1
          qcoef(lev0,i) = qcoef(lev0,i)+pcoef(lev0,i)*(flbc1(i)+     ! s2
     |      0.5*flbc2(i)*dz)/(flbc1(i)-0.5*flbc2(i)*dz)
          fcoef(lev0,i) = fcoef(lev0,i)-pcoef(lev0,i)*flbc3(i)*dz/   ! s4
     |      (flbc1(i)-0.5*flbc2(i)*dz)
          pcoef(lev0,i) = 0.                                         ! s1
        enddo ! i=lon0,lon1
!
! Upper boundary (this routine is normally called with iubc==1):
        if (iubc==0.or.iubc==1) then ! upward number flux is in fubc(i,lat)
          do i=lon0,lon1
            pcoef(lev1,i) = 1.+0.5*dz*thdiff(lev1,i)
            qcoef(lev1,i) = pcoef(lev1,i)-2.
            rcoef(lev1,i) = 0.
            fcoef(lev1,i) = -grav*rmass*fubc(i,lat)*dz/(p0*ax(lev1,i)*
     |                      avo)
          enddo
        else                         ! value at upper boundary is in fubc
          do i=lon0,lon1
            pcoef(lev1,i) = 0.5
            qcoef(lev1,i) = 0.5
            rcoef(lev1,i) = 0.
            fcoef(lev1,i) = fubc(i,lat)
          enddo
        endif

        if (idebug > 0) then
          call addfld('PCOEF' ,' ',' ',pcoef,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('QCOEF' ,' ',' ',qcoef,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('RCOEF' ,' ',' ',rcoef,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FCOEF' ,' ',' ',fcoef,
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
        endif 
!
! Solve tridiagonal system:
!
        call trsolv(pcoef,qcoef,rcoef,fcoef,fout(:,lon0:lon1,lat),
     |    lev0,lev1,lev0,lev1,lon0,lon1,nlonp4,lat,0)
!
        if (idebug > 0)
     |    call addfld('TRSOLV',' ',' ',fout(:,lon0:lon1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
!
! End latitude scan:
      enddo ! lat=lat0,lat1
!
! Filter fout (all latitudes):
      call filter_minor(fout(:,lon0:lon1,lat0:lat1),
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Resume latitude scan:
      do lat=lat0,lat1

        if (idebug > 0)
     |    call addfld('FILTER',' ',' ',fout(:,lon0:lon1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Note: Periodic points for fout and fout_nm are set here in tgcm24.
!       Am skipping this for now, as it may not be necessary.

!
! Time smoothing:
        do i=lon0,lon1
          do k=lev0,lev1
            fout_nm(k,i,lat) = dtsmooth*fin(k,i,lat)+dtsmooth_div2*
     |        (fnm(k,i,lat)+fout(k,i,lat))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1
! 
! Insure positive fout:
        do i=lon0,lon1
          do k=lev0+1,lev1-1
            if (fout(k,i,lat) < small(k)*xy) 
     |        fout(k,i,lat) = small(k)*xy
            if (fout_nm(k,i,lat) < small(k)*xy) 
     |        fout_nm(k,i,lat) = small(k)*xy
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

        if (idebug > 0) then
          call addfld('FOUT'   ,' ',' ',fout(:,lon0:lon1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FOUT_NM',' ',' ',fout_nm(:,lon0:lon1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
          call addfld('FNM'    ,' ',' ',fnm(:,lon0:lon1,lat),
     |      'lev',lev0,lev1,'lon',i0,i1,lat)
        endif

      enddo ! lat=lat0,lat1

      end subroutine minor
!-----------------------------------------------------------------------
      subroutine minor2(tn,o2,ox,w,difkk,fin,fnm,fhd,flbc,fubc,small,
     |  rmass,phix,alfax,xy,ilbc,iubc,fout,fout_nm,
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Advances minor species by one time step (Currently called for hox only).
!
! Lower boundary (for hox, ilbc==0):
!   If ilbc == 0, then
!     flbc(1) = A, flbc(2) = B, flbc(3) = C, where A*DPSX/DZ + B*PSX + C = 0.
!   If ilbc == 1,      (not the normal case)
!     flbc(3) = total upward number flux of F at lower boundary.
!
! Upper boundary (for hox, iubc==1):
!   If iubc == 0, then
!     fubc = diffusive upward number flux of F at upper boundary.
!   If iubc == 1, then
!     fubc = total upward number flux of F at upper boundary.
!   If iubc == 2, then
!     fubc = value of PSX at upper boundary.
!
! Normally (when called from comp_hox, normally ilbc==0 and iubc==1) 
!
      use params_module,only: nlonp4,dz
      use lbc,only: b,fb
      use cons_module,only: rmassinv_o2,rmassinv_ox,rmassinv_n2,p0,
     |  expz,expzmid,expzmid_inv,boltz,rmass_o2,rmass_ox,rmass_n2,
     |  dtx2inv,grav,avo,dtsmooth,dtsmooth_div2
      use input_module,only: difhor
      use fields_module,only: tlbc
      use addfld_module,only: addfld
      implicit none
!
! Input Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1,iubc,ilbc
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in)::
     |  tn,o2,ox,w,difkk,fin,fhd
      real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(inout) ::
     |  fnm
      real,dimension(lon0:lon1),intent(in) ::
     |  flbc(lon0:lon1,3,lat0:lat1), ! 3 components of lower boundary
     |  fubc(lon0:lon1,lat0:lat1)    ! upper boundary
      real,intent(in) :: rmass,phix(3),alfax,small(lev1-lev0+1),xy
!
! Output Args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) :: fout
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(inout) :: fout_nm
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1,lat
      real :: phi(2,3),tau,t00,alfa12,alfa21,alfax1,alfax2
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  hadvec,    ! horizontal diffusion (output of routine advec)
     |  source0,   ! portion of number density source proportional to minor sp
     |  sourcex,   ! portion of number density source independent of minor sp
     |  barm,      ! mbar (probably same as 4d barm, which could be passed in) 
     |  barmi,     ! mbar at interfaces ! s11
     |  s0mbar,    ! source0*mass/mbar  ! s15
     |  dmdz,      ! dm/dz              ! s10
     |  oxi,       ! o2 at interfaces   ! s9
     |  o2i,       ! ox at interfaces   ! s8
     |  oxdz,      ! dox/dz             ! s7
     |  o2dz,      ! do2/dz             ! s6
     |  tni,       ! tn at interfaces   ! s5
     |  phialfa_ox,                     ! s1
     |  alfa_o2,                        ! s2
     |  alfa_ox,                        ! s3
     |  phialfa_o2,                     ! s4
     |  thdiff,    ! thermal diffusion  ! s12
     |  ax,                             ! s11
     |  pcoef, qcoef, rcoef, fcoef      ! s1,2,3,4
      real,dimension(lon0:lon1) :: 
     |  lbc_o2, lbc_ox, lbcmbar, dfactor,
     |  flbc1,  flbc2,  flbc3,   fubc1
! 
      i0=lon0 ; i1=lon1 ; nk=lev1-lev0+1 ; nkm1=nk-1
!
! Check inputs:
!
!     write(6,"(/,'minor2: iubc,ilbc=',2i3)") iubc,ilbc
!     write(6,"('minor2: rmass=',e12.4)") rmass
!     write(6,"('minor2: phix =',3e12.4)") phix
!     write(6,"('minor2: alfax=',3e12.4)") alfax
!     write(6,"('minor2: xy   =',3e12.4)") xy

!     do lat=lat0,lat1
!       call addfld('TN_MNR' ,' ',' ',tn (:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('O2_MNR' ,' ',' ',o2 (:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('OX_MNR' ,' ',' ',ox (:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FIN_MNR',' ',' ',fin(:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FNM_MNR',' ',' ',fnm(:,:,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('PDH_MNR',' ',' ',fhd(:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
!       write(6,"('minor2: lat=',i3,' flbc(1)=',/,(6e12.4))") 
!    |    lat,flbc(:,1,lat)
!       write(6,"('minor2: lat=',i3,' flbc(2)' ,/,(6e12.4))") 
!    |    lat,flbc(:,2,lat)
!       write(6,"('minor2: lat=',i3,' flbc(3)=',/,(6e12.4))") 
!    |    lat,flbc(:,3,lat)
!       write(6,"('minor2: lat=',i3,' fubc ='  ,/,(6e12.4))") 
!    |    lat,fubc(:,lat)
!       write(6,"('minor2: lat=',i3,' small='  ,/,(6e12.4))") 
!    |    lat,small
!     enddo ! lat=lat0,lat1

      phi(:,1)=(/0.  ,0.673/)
      phi(:,2)=(/1.35,0.   /)
      phi(:,3)=(/1.11,0.769/)
      tau=1.86E+3
      t00=273.
!
! Advec input:
! real,dimension(lev0:lev1,lon0-2:lon1+2,lat-2:lat+2),intent(in) :: f
!
      do lat=lat0,lat1
        call advec(fin(:,:,lat-2:lat+2),hadvec,lev0,lev1,lon0,lon1,lat)
!       call addfld('HADVEC',' ',' ',hadvec,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! In tgcm24/cmphox, s1 and s2 are set to zero before calling minor2.
! In minor2, they are transferred to s14 and s15.
!
        source0(:,:) = 0.   ! s1,s14
        sourcex(:,:) = 0.   ! s2,s15
!
! Lower boundary:
        do i=lon0,lon1
          lbc_o2(i) = b(i,1,1)*o2(lev0,i,lat)+                    ! t4
     |                b(i,1,2)*ox(lev0,i,lat)+fb(i,1,lat)
          lbc_ox(i) = b(i,2,1)*o2(lev0,i,lat)+                    ! t5
     |                b(i,2,2)*ox(lev0,i,lat)+fb(i,2,lat)
          lbcmbar(i) = 1./(lbc_o2(i)*rmassinv_o2+                 ! t6
     |                     lbc_ox(i)*rmassinv_ox+
     |                     (1.-lbc_o2(i)-lbc_ox(i))*rmassinv_n2)
        enddo ! i=lon0,lon1
!
! Set mean molecular weight (could probably be passed in from 4d barm):
        do i=lon0,lon1
          do k=lev0,lev1
            barm(k,i) = 1./(o2(k,i,lat)*rmassinv_o2+         ! s12
     |                      ox(k,i,lat)*rmassinv_ox+
     |        (1.-o2(k,i,lat)-ox(k,i,lat))*rmassinv_n2)
            s0mbar(k,i) = source0(k,i)*rmass*boltz*tn(k,i,lat)/ ! s15
     |        (p0*expz(k)*barm(k,i))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1
!
! Lower boundary:
        do i=lon0,lon1
          barmi(lev0,i) = 0.5*(lbcmbar(i)+barm(lev0,i))     ! s11
          dmdz(lev0,i)  = (barm(lev0,i)-lbcmbar(i))/dz      ! s10
          oxi (lev0,i)  = 0.5*(lbc_ox(i)+ox(lev0,i,lat))    ! s9
          o2i (lev0,i)  = 0.5*(lbc_o2(i)+o2(lev0,i,lat))    ! s8
          oxdz(lev0,i)  = (ox(lev0,i,lat)-lbc_ox(i))/dz     ! s7
          o2dz(lev0,i)  = (o2(lev0,i,lat)-lbc_o2(i))/dz     ! s6
!
! Levels 2->top
          do k=lev0+1,lev1
            barmi(k,i) = 0.5*(barm(k,i)+barm(k-1,i))        ! s11
            dmdz(k,i) = (barm(k,i)-barm(k-1,i))/dz          ! s10
            oxi(k,i) = 0.5*(ox(k,i,lat)+ox(k-1,i,lat))      ! s9
            o2i(k,i) = 0.5*(o2(k,i,lat)+o2(k-1,i,lat))      ! s8
            oxdz(k,i) = (ox(k,i,lat)-ox(k-1,i,lat))/dz      ! s7
            o2dz(k,i) = (o2(k,i,lat)-o2(k-1,i,lat))/dz      ! s6
          enddo ! k=lev0+1,lev1
        enddo ! i=lon0,lon1

!       call addfld('BARMX',' ',' ',barm ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('BARMI',' ',' ',barmi,
!    |    'ilev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DMDZ' ,' ',' ',dmdz ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('OXI'  ,' ',' ',oxi  ,
!    |    'ilev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('O2I'  ,' ',' ',o2i  ,
!    |    'ilev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('OXDZ' ,' ',' ',oxdz ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('O2DZ' ,' ',' ',o2dz ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('S0MBAR',' ',' ',s0mbar ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

!
! Tn at interfaces (s5):
! (Note bottom boundary of tn is stored in top slot)
! 5/15/06 btf: bottom boundary of TN is in tlbc:
        do i=lon0,lon1
!         tni(lev0,i) = tn(lev1  ,i,lat) ! bottom boundary in top slot (s5)
          tni(lev0,i) = tlbc(i,lat)
          tni(lev1,i) = tn(lev1-1,i,lat) ! top boundary same as lev1-1
          do k=lev0+1,lev1-1
            tni(k,i) = 0.5*(tn(k,i,lat)+tn(k-1,i,lat))
          enddo
        enddo

! Evaluate alfa12, alfa21, alfam1, alfam2
        alfa12=phi(1,2)-phi(1,3)
        alfa21=phi(2,1)-phi(2,3)
        alfax1=phix(1)-phix(3)
        alfax2=phix(2)-phix(3)
!
        do i=lon0,lon1
          do k=lev0,lev1
            phialfa_ox(k,i) = -(phi(1,3)+alfa12*oxi(k,i))                ! s1
            alfa_o2(k,i) = alfa12*o2i(k,i)                               ! s2
            alfa_ox(k,i) = alfa21*oxi(k,i)                               ! s3
            phialfa_o2(k,i) = -(phi(2,3)+alfa21*o2i(k,i))                ! s4
            thdiff(k,i) = ((alfax1*phialfa_o2(k,i)-alfax2*alfa_ox(k,i))* ! s12
     |        (o2dz(k,i)-(1.-(rmass_o2+dmdz(k,i))/barmi(k,i))*o2i(k,i))+ 
     |        (alfax2*phialfa_ox(k,i)-alfax1*alfa_o2(k,i))*
     |        (oxdz(k,i)-(1.-(rmass_ox+dmdz(k,i))/barmi(k,i))*oxi(k,i)))
     |        /(phialfa_ox(k,i)*phialfa_o2(k,i)-alfa_o2(k,i)*
     |        alfa_ox(k,i))+1.-(rmass+dmdz(k,i))/barmi(k,i)
            dmdz(k,i) = dmdz(k,i)/barmi(k,i)                             ! s10
            ax(k,i) = -barmi(k,i)/(tau*rmass_n2)*(t00/tni(k,i))**0.25/   ! s11
     |        (phix(3)+alfax1*o2i(k,i)+alfax2*oxi(k,i))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

!       call addfld('PHIALFOX',' ',' ',phialfa_ox,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('ALFAO2'  ,' ',' ',alfa_o2   ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('ALFAOX'  ,' ',' ',alfa_ox   ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('PHIALFO2',' ',' ',phialfa_o2,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('THDIFF'  ,' ',' ',thdiff    ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DMDZ'    ,' ',' ',dmdz      ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('AX'      ,' ',' ',ax        ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

        do i=lon0,lon1
          do k=lev0+1,lev1-1
            thdiff(k,i) = thdiff(k,i)-alfax*(tni(k+1,i)-tni(k-1,i))/     ! s12
     |        (2.*dz*tni(k,i))
          enddo ! k=lev0+1,lev1-1
!
! Lower boundary:
          thdiff(lev0,i) = thdiff(lev0,i)-alfax*
     |      (tni(lev0+1,i)-tni(lev0,i))/(dz*tni(lev0,i))
!
! Upper boundary:
          thdiff(lev1,i) = thdiff(lev1,i)-alfax*
     |      (tni(lev1,i)-tni(lev1-1,i))/(dz*tni(lev1,i))
        enddo ! i=lon0,lon1

!       call addfld('THDIFF1'  ,' ',' ',thdiff    ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

        if (difhor == 1) then
          call dfact(dfactor,lon0,lon1,lat)
        else
          dfactor(:) = 1.
        endif
!
! Set coefficients for tridiagonal solver:
        do i=lon0,lon1
          do k=lev0,lev1-1
!
! Note difkk (mgw.F) will match tgcm24 only if both are set zero. 
! This is because mgw gravity wave parameterization is new in timegcm1.
!
            pcoef(k,i) = ax(k,i)/dz*(1./dz+0.5*thdiff(k,i))-expz(k)*     ! s1
     |        (expzmid_inv*difkk(k,i,lat)*dfactor(i)*(1./dz-0.5*
     |        dmdz(k,i))+0.25*(w(k,i,lat)+w(k+1,i,lat)))/dz
! 
            rcoef(k,i) = ax(k+1,i)/dz*(1./dz-0.5*thdiff(k+1,i))-expz(k)* ! s3
     |        (expzmid*difkk(k+1,i,lat)*dfactor(i)*(1./dz+0.5*
     |        dmdz(k+1,i))-0.25*(w(k,i,lat)+w(k+1,i,lat)))/dz
!
            qcoef(k,i) = -(ax(k,i)/dz*(1./dz-0.5*thdiff(k,i))+ax(k+1,i)/ ! s2
     |        dz*(1./dz+0.5*thdiff(k+1,i)))+expz(k)*((expzmid_inv*
     |        difkk(k,i,lat)*(1./dz+0.5*dmdz(k,i))+expzmid*
     |        difkk(k+1,i,lat)*(1./dz-0.5*dmdz(k+1,i)))*dfactor(i)/
     |        dz-source0(k,i)+dtx2inv)
!
            fcoef(k,i) = expz(k)*(fnm(k,i,lat)*dtx2inv-hadvec(k,i)+      ! s4
     |        fhd(k,i,lat)+s0mbar(k,i))
          enddo ! k=lev0,lev1-1
        enddo ! i=lon0,lon1
        
!       call addfld('AX'     ,' ',' ',ax     ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('THDIFF' ,' ',' ',thdiff ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DIFKK'  ,' ',' ',difkk(:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('OMEGA'  ,' ',' ',w(:,i0:i1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('SOURCE0',' ',' ',source0,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('DMDZ'   ,' ',' ',dmdz   ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('HADVEC' ,' ',' ',hadvec(lev0:lev1-1,:),
!    |    'lev',lev0,lev1-1,'lon',i0,i1,lat)
!       call addfld('S0MBAR' ,' ',' ',s0mbar ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

!       call addfld('PCOEF' ,' ',' ',pcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('RCOEF' ,' ',' ',rcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('QCOEF' ,' ',' ',qcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FCOEF' ,' ',' ',fcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

!
! If ilbc==1, set flux boundary condition at bottom (normally, this 
!   is called with ilbc==0).
        if (ilbc==1) then ! (normally not true)
          do i=lon0,lon1
            flbc1(i) = 1.                                           ! t1
            flbc2(i) = dmdz(lev0,i)-w(lev0,i,lat)/difkk(lev0,i,lat) ! t2
            flbc3(i) = flbc(i,3,lat)*grav*rmass/(difkk(lev0,i,lat)* ! t3
     |        p0*expz(1)*expzmid_inv*avo) 
          enddo ! i=lon0,lon1
        else
!
! Use input flbc(lon0:lon1,3,lat0:lat1) ! 3 components of lower boundary
          do i=lon0,lon1
            flbc1(i) = flbc(i,1,lat)
            flbc2(i) = flbc(i,2,lat)
            flbc3(i) = flbc(i,3,lat)
          enddo ! i=lon0,lon1
        endif
!
! Modify ex if iubc==1:
        if (iubc==1) then
          do i=lon0,lon1
            fubc1(i) = expz(lev1-1)*expzmid*w(lev1,i,lat)/ax(lev1,i) ! t7
            thdiff(lev1,i) = thdiff(lev1,i)-fubc1(i)                 ! s12
          enddo !i=lon0,lon1
        endif
!
! Lower boundary:
        do i=lon0,lon1
          qcoef(lev0,i) = qcoef(lev0,i)+pcoef(lev0,i)*(flbc1(i)+
     |      0.5*flbc2(i)*dz)/(flbc1(i)-0.5*flbc2(i)*dz)
          fcoef(lev0,i) = fcoef(lev0,i)-pcoef(lev0,i)*flbc3(i)*dz/
     |      (flbc1(i)-0.5*flbc2(i)*dz)
          pcoef(lev0,i) = 0.
        enddo ! i=lon0,lon1
!
! Upper boundary (this routine is normally called with iubc==1):
        if (iubc==0.or.iubc==1) then ! upward number flux is in fubc(i,lat)
          do i=lon0,lon1
            pcoef(lev1,i) = 1.+0.5*dz*thdiff(lev1,i)
            qcoef(lev1,i) = pcoef(lev1,i)-2.
            rcoef(lev1,i) = 0.
            fcoef(lev1,i) = -grav*rmass*fubc(i,lat)*dz/(p0*ax(lev1,i)*
     |                      avo)
          enddo
        else                         ! value at upper boundary is in fubc
          do i=lon0,lon1
            pcoef(lev1,i) = 0.5
            qcoef(lev1,i) = 0.5
            rcoef(lev1,i) = 0.
            fcoef(lev1,i) = fubc(i,lat)
          enddo
        endif

!       call addfld('PCOEF' ,' ',' ',pcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('QCOEF' ,' ',' ',qcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('RCOEF' ,' ',' ',rcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FCOEF' ,' ',' ',fcoef ,
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Solve tridiagonal system:
!
        call trsolv(pcoef,qcoef,rcoef,fcoef,fout(:,lon0:lon1,lat),
     |    lev0,lev1,lev0,lev1,lon0,lon1,nlonp4,lat,0)

!       call addfld('TRSOLV',' ',' ',fout(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Note: Periodic points for fout and fout_nm are set here in tgcm24.
!       Am skipping this for now, as it may not be necessary.

! 
! Insure positive fout:
        do i=lon0,lon1
          do k=lev0,lev1
            if (fout(k,i,lat) < small(k)*xy) 
     |        fout(k,i,lat) = small(k)*xy
            if (fout_nm(k,i,lat) < small(k)*xy) 
     |        fout_nm(k,i,lat) = small(k)*xy
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1
!
! End latitude scan:
      enddo ! lat=lat0,lat1
!
! Filter fout:
      call filter_minor(fout(:,lon0:lon1,lat0:lat1),
     |  lev0,lev1,lon0,lon1,lat0,lat1)
!
! Resume latitude scan:
      do lat=lat0,lat1

!       call addfld('FILTER2',' ',' ',fout(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Time smoothing:
        do i=lon0,lon1
          do k=lev0,lev1
            fnm(k,i,lat) = fout_nm(k,i,lat)    
            fout_nm(k,i,lat) = dtsmooth*fin(k,i,lat)+dtsmooth_div2*
     |        (fnm(k,i,lat)+fout(k,i,lat))
          enddo ! k=lev0,lev1
        enddo ! i=lon0,lon1

!       call addfld('FOUT'   ,' ',' ',fout(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FOUT_NM',' ',' ',fout_nm(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)
!       call addfld('FNM'    ,' ',' ',fnm(:,lon0:lon1,lat),
!    |    'lev',lev0,lev1,'lon',i0,i1,lat)

      enddo ! lat=lat0,lat1

      end subroutine minor2
!-----------------------------------------------------------------------
      subroutine filter_minor(fout,lev0,lev1,lon0,lon1,lat0,lat1)
!
! Filter minor species (called from minor after trsolv)
!
      use params_module,only: nlonp4,nlon,nlat
      use filter_module,only: filter2
      use addfld_module,only: addfld
!
#ifdef MPI
      use mpi_module,only: mp_gatherlons_f3d,mp_scatterlons_f3d,mytidi
      implicit none
#else
      implicit none
      integer :: mytidi=0
#endif
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1
      real,intent(inout) :: fout(lev0:lev1,lon0:lon1,lat0:lat1)
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: i,j,ier
      real :: f_ik(nlonp4,lev0:lev1),f_kij(lev0:lev1,nlonp4,lat0:lat1)
      real :: fmin,fmax
!
#ifdef VT
!     code = 123 ; state = 'filter_minor' ; activity='Filtering'
      call vtbegin(123,ier)
#endif
!
! Define lons in f_ki from current task:
      f_kij = 0.
      do j=lat0,lat1
        do i=lon0,lon1
          f_kij(:,i,j) = fout(:,i,j)
        enddo
      enddo
!
#ifdef MPI
!
! Gather longitudes into tasks in first longitude column of task table
!   (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0
!   gather lons from other tasks in that row). This includes all latitudes.
!
      call mp_gatherlons_f3d(f_kij,lev0,lev1,lon0,lon1,lat0,lat1,1)
#endif
!
! Only leftmost tasks at each j-row of tasks does the global filtering:
      if (mytidi==0) then
!
! Loop through subdomain latitudes:
        latscan: do j=lat0,lat1
!
! Define 2d array with all longitudes for filter2:
          do i=1,nlonp4
            f_ik(i,:) = f_kij(:,i,j)
          enddo ! i=1,nlonp4
!
! Do the filtering (requires fft):
          call filter2(f_ik,lev0,lev1,j)
!
! Return filtered array to f_kij:
          do i=1,nlonp4
            f_kij(:,i,j) = f_ik(i,:)
          enddo ! i=1,nlonp4
        enddo latscan ! j=lat0,lat1
      endif ! mytidi==0
#ifdef MPI
!
! Now leftmost task at each j-row must redistribute filtered data
! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude):
!
      call mp_scatterlons_f3d(f_kij,lev0,lev1,lon0,lon1,lat0,lat1,1)
#endif
!
! Return filtered array to fout at current task longitudes:
      do j=lat0,lat1
        do i=lon0,lon1
          fout(:,i,j) = f_kij(:,i,j)
        enddo
      enddo ! j=lat0,lat1
!
#ifdef VT
!     code = 123 ; state = 'filter_minor' ; activity='Filtering'
      call vtend(123,ier)
#endif
      end subroutine filter_minor