!
      subroutine addiag(tn,o2,ox,vn,vc,barm,xnmbar,xnmbari,xnmbarm,z,
     |  lon0,lon1,lev0,lev1,lat0,lat1)
!
! Calculate needed terms vc, barm, xnmbar[i,m], and Z:
!
      use cons_module,only: cs,rmassinv,dz,dzgrav,freq_semidi,dt,p0,
     |  boltz,expz,expzmid,expzmid_inv,freq_3m3
      use bndry_module,only: zb,zb2,zba,bnd,bnd2,bnda,ci,
     |  zb3m3,bnd3m3,lbc_gswm_addiag,zbplanet,zbkelvin
      use init_module,only: iter,istep,igetgswm
      use input_module,only: ncep_ncfile,ncep_reanalysis,ecmwf_ncfile
      use ncep_module,only: zncep
      use ncep_reanalysis_module,only: z_ncep
      use ecmwf_module,only: z_ecmwf
      use fields_module,only: tlbc
      use addfld_module,only: addfld
      implicit none
!
! Input args:
      integer,intent(in) :: lon0,lon1,lev0,lev1,lat0,lat1
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: 
     |  tn, ! neutral temperature (deg K)
     |  o2, ! molecular oxygen (mmr)
     |  ox, ! atomic oxygen (mmr)
     |  vn  ! meridional wind velocity (cm/s)
!
! Output args:
!     real,dimension(lev0:lev1,lon0  :lon1  ,lat0  :lat1  ),
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  vc  ,
     |  barm,
     |  z   ,
     |  xnmbar , 
     |  xnmbari,
     |  xnmbarm
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: k,i,j,nlevs,ier
      real :: 
     |  barm1(lon0:lon1),
     |  tni  (lev0:lev1,lon0:lon1), ! tn at interfaces
     |  expzi(lev0:lev1,lon0:lon1), ! e(-z) at interfaces
     |  w1   (lev0:lev1,lon0:lon1),
     |  zlbc_diag(lev0:lev1,lon0:lon1)
      complex :: expt,expu
      real :: fmin,fmax
!
#ifdef VT
!     code = 116 ; state = 'addiag' ; activity='ModelCode'
      call vtbegin(116,ier)
#endif
      nlevs = lev1-lev0+1
!
! Check for gswm (see bndry.F):
      if (igetgswm > 0) then
        call lbc_gswm_addiag(z(1,lon0:lon1,lat0:lat1),
     |    lon0,lon1,lat0,lat1)

!       do j=lat0,lat1
!         write(6,"('addiag after lbc_gswm_addiag: j=',i3,
!    |      ' z(1,:,j)=',/,(6e12.4))") j,z(1,:,j)
!       enddo

      endif
!
! Latitude scan:
      do j=lat0,lat1

        do k=lev0,lev1
          zlbc_diag(k,lon0:lon1) = z(1,lon0:lon1,j) ! redundant in vertical
        enddo
!       call addfld('ZGSWM',' ',' ',zlbc_diag,
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! vc = cos(phi)*v
!
        do i=lon0,lon1
          do k=lev0,lev1
            vc(k,i,j) = cs(j)*vn(k,i,j)
          enddo
        enddo

!       call addfld('diag_vn',' ',' ',vn(:,lon0:lon1,j),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!       call addfld('diag_vc',' ',' ',vc(:,lon0:lon1,j),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! barm = mean molecular weight (k+1/2):
!
        do i=lon0,lon1
          do k=lev0,lev1
            barm(k,i,j) = 1./
     |        (o2(k,i,j)*rmassinv(1)+ox(k,i,j)*rmassinv(2)+
     |        (1.-o2(k,i,j)-ox(k,i,j))*rmassinv(3))
            xnmbarm(k,i,j)=p0*expz(k)*barm(k,i,j)/
     |        (boltz*tn(k,i,j))
          enddo
!
! 5/4/06 btf: set boundaries of xnmbarm:
!         xnmbarm(lev0,i,j)=p0*expz(lev0)*barm(lev0,i,j)/
!    |      (boltz*tlbc(i,j))
!         xnmbarm(lev1,i,j)=p0*expz(lev1)*barm(lev1,i,j)/
!    |      (boltz*tn(lev1-1,i,j))
        enddo
!
! barm1 = barm(k=0) (linear extrapolation)
!
        do i=lon0,lon1
          barm1(i) = 1.5*barm(1,i,j)-0.5*barm(2,i,j)
        enddo
!
! barm(k) = 0.5*(barm(k+1/2)+barm(k-1/2)), k = kmaxp1,2,1
!
        do i=lon0,lon1
          do k=lev1,lev0+1,-1
            barm(k,i,j) = 0.5*(barm(k,i,j)+barm(k-1,i,j))
          enddo
        enddo
!
! barm(1) = barm1
!
        do i=lon0,lon1
          barm(lev0,i,j) = barm1(i)
        enddo

!       call addfld('barm',' ',' ',barm(:,lon0:lon1,j),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! xnmbar is output by sub oplus (oplus.F)
!
! xnmbar = p0*e(-z)*barm/kT at midpoints (used in conversion from mmr to cm3).
! (set by oplus)
!       do i=lon0,lon1
!         do k=lev0,lev1-1
!           xnmbar(k,i,j)=p0*expz(k)*.5*(barm(k,i,j)+barm(k+1,i,j))
!    |        /(boltz*tn(k,i,j))
!         enddo
!       enddo
!
! xnmbari = p0*e(-z)*barm/kT at interfaces (used by qrj and qinite):

        do i=lon0,lon1
!         tni(1,i) = tn(lev1,i,j) ! tn bottom boundary is stored in top slot
          tni(lev0,i) = tlbc(i,j)    ! Lower boundary is in tlbc
          expzi(1,i) = expzmid_inv*expz(1)
          do k=lev0+1,lev1-1
            tni(k,i) = .5*(tn(k-1,i,j)+tn(k,i,j))
            expzi(k,i) = expzmid_inv*expz(k)
          enddo
          tni(lev1,i) = tn(lev1-1,i,j) ! nlevp1 <- nlev
          expzi(lev1,i) = expzmid*expz(lev1-1)
          do k=lev0,lev1
            xnmbari(k,i,j) = p0*expzi(k,i)*barm(k,i,j)/
     |        (boltz*tni(k,i))
          enddo
        enddo
!       call addfld('xnmbari',' ',' ',xnmbari(:,lon0:lon1,j),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! w1 = barm
        do i=lon0,lon1
          do k=lev0,lev1-1
            w1(k,i) = (barm(k,i,j)+barm(k+1,i,j))*0.5
          enddo
        enddo
!       call addfld('w1a',' ',' ',w1,'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! w1 = tn/w1 (old model comment: s1=s2/s1=(t+t0)/m)
        do i=lon0,lon1
          do k=lev0,lev1-1
            w1(k,i) = tn(k,i,j)/w1(k,i)
          enddo
        enddo
!       call addfld('diag_tn',' ',' ',tn(:,lon0:lon1,j),
!    |    'lev',lev0,lev1,'lon',lon0,lon1,j)
!       call addfld('w1b',' ',' ',w1,'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! w1=(ds*r/g)*w1
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            w1(k,i) = (dz/dzgrav) * w1(k,i)
          enddo
        enddo
!       call addfld('w1c',' ',' ',w1,'lev',lev0,lev1,'lon',lon0,lon1,j)
!
! Cannot run with both GSWM and NCEP (this should have been checked
!   in input, but just in case):
! 3/25/04 btf: gswm+ncep is allowed.
!       if (igetswm > 0 .and. len_trim(ncep_ncfile) > 0) then
!         write(6,"(/,'>>> ADDIAG: Cannot run with both GSWM and NCEP:',
!    |      ' gswm=',l1,' ncep_ncfile=',a)") gswm,ncep_ncfile
!         call shutdown('addiag')
!       endif
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
! Neither GSWM nor NCEP nor ECMWF:
        if (igetgswm <= 0 .and. 
     |    len_trim(ncep_ncfile)     <= 0.and.
     |    len_trim(ncep_reanalysis) <= 0.and.
     |    len_trim(ecmwf_ncfile)    <= 0) then
!
! z(1) = zb
          expt=cexp(ci*freq_semidi*dt*iter)
          do i=lon0,lon1
            z(1,i,j) = real(zb(j)*bnd(i)*expt) 
          enddo
!
! Add in effect of annual tidal component to lbc
!
          expt = 1.
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+real(zba(j)*bnda(i)*expt)
          enddo
!
! Add in effect of (1,1) tidal component to lbc
!
          expt = cexp(ci*.5*freq_semidi*dt*iter)
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+real(zb2(j)*bnd2(i)*expt)
          enddo
!
! Add planetary and kelvin waves:
!
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+zbplanet(i,j)
            z(1,i,j) = z(1,i,j)+zbkelvin(i,j)
          enddo
!
! Add (3,-3) 2-day wave (unconditional on ncep and gswm)
!
          expt = cexp(ci*freq_3m3*dt*iter)
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+real(zb3m3(j)*bnd3m3(i)*expt)
          enddo
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
! NCEP or ECMWF run (no GSWM):
        elseif (igetgswm <= 0 .and. (len_trim(ncep_ncfile) > 0.or.
     |    len_trim(ncep_reanalysis) > 0.or.len_trim(ecmwf_ncfile) > 0)) 
     |    then
!
! Add ncep or ecmwf boundaries:
!
! Add in effect of (1,1) and (2,2) tidal component to lbc
!
          expu = cexp(ci*freq_semidi*dt*iter)
          expt = cexp(ci*.5*freq_semidi*dt*iter)
          do i = lon0,lon1
            if (len_trim(ncep_ncfile) > 0) then
              z(1,i,j) = zncep(i,j,2)+real(zb2(j)*bnd2(i)*expt)+
     |          real(zb(j)*bnd(i)*expu)
            elseif (len_trim(ncep_reanalysis) > 0) then
              z(1,i,j) = z_ncep(i,j)+real(zb2(j)*bnd2(i)*expt)+
     |          real(zb(j)*bnd(i)*expu)
            elseif (len_trim(ecmwf_ncfile) > 0) then
              z(1,i,j) = z_ecmwf(i,j)+real(zb2(j)*bnd2(i)*expt)+
     |          real(zb(j)*bnd(i)*expu)
            endif
          enddo
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
! GSWM run (no NCEP or ECMWF) (note lbc_gswm_addiag was called above):
        elseif (igetgswm > 0 .and. 
     |          len_trim(ncep_ncfile)     <= 0.and.
     |          len_trim(ncep_reanalysis) <= 0.and.
     |          len_trim(ecmwf_ncfile)    <= 0) then
!
! Add in effect of annual tidal component to lbc
!
          expt = 1.
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+real(zba(j)*bnda(i)*expt)
          enddo
!
! Add planetary and kelvin waves:
!
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+zbplanet(i,j)
            z(1,i,j) = z(1,i,j)+zbkelvin(i,j)
          enddo
!
! Add effect of (3,-3) 2-day wave (unconditional on ncep and gswm)
!
          expt = cexp(ci*freq_3m3*dt*iter)
          do i=lon0,lon1
            z(1,i,j) = z(1,i,j)+real(zb3m3(j)*bnd3m3(i)*expt)
          enddo

!         write(6,"('addiag: j=',i3,' z(1,:,j)=',/,(6e12.4))")
!    |      j,z(1,:,j)

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
! GSWM and NCEP run:
        elseif (igetgswm > 0 .and. (len_trim(ncep_ncfile) > 0.or.
     |    len_trim(ncep_reanalysis) > 0.or.len_trim(ecmwf_ncfile) > 0)) 
     |    then
          do i = lon0,lon1
            if (len_trim(ncep_ncfile) > 0) then
              z(1,i,j) = z(1,i,j)+zncep(i,j,2)
            elseif (len_trim(ncep_reanalysis) > 0) then
              z(1,i,j) = z(1,i,j)+z_ncep(i,j)
            elseif (len_trim(ecmwf_ncfile) > 0) then
              z(1,i,j) = z(1,i,j)+z_ecmwf(i,j)
            endif
          enddo
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
! End GSWM or NCEP:
        endif
!
! z(k+1)=s1(k)+z(k)
!
        do i=lon0,lon1
          do k=lev0,lev1-1
            z(k+1,i,j) = w1(k,i)+z(k,i,j)
          enddo
        enddo

!       do k=lev0,lev1
!         write(6,"('addiag: j=',i3,' k=',i3,' z(k,:,j)=',/,(6e12.4))")
!    |      j,k,z(k,:,j)
!       enddo

!       call addfld('ADIAG_Z','Geopotential from addiag',
!    |    'cm',z(:,lon0:lon1,j),'ilev',lev0,lev1,'lon',lon0,lon1,j)

      enddo ! j=lat0,lat1
!
#ifdef VT
!     code = 116 ; state = 'addiag' ; activity='ModelCode'
      call vtend(116,ier)
#endif
      end subroutine addiag
