!
      subroutine addiag(tn,o2,ox,he,vn,vc,barm,xnmbar,xnmbari,xnmbarm,z,
     |  zg,n2,lon0,lon1,lev0,lev1,lat0,lat1)
!
! Calculate needed terms vc, barm, xnmbar[i,m], and Z:
!
      use cons_module,only: cs,rmassinv,rmassinv_he,dz,dzgrav,
     |  freq_semidi,dt,p0,boltz,expz,expzmid,expzmid_inv,freq_3m3
      use init_module,only: iter,igetgswm
      use input_module,only: ncep_ncfile,ncep_reanalysis,ecmwf_ncfile
      use ncep_module,only: zncep
      use ncep_rean_module,only: z_ncep
      use ecmwf_module,only: z_ecmwf
      use fields_module,only: tlbc
      use addfld_module,only: addfld
      use lbc,only: z_lbc
      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)
     |  he, ! helium (mmr)
     |  vn  ! meridional wind velocity (cm/s)
!
! Output args:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  vc  ,
     |  barm,
     |  z   ,
     |  zg  ,
     |  n2  ,
     |  xnmbari,
     |  xnmbarm
!
! xnmbar is set by oplus.F, so make it intent(in) here to
! avoid warning from intel compiler:
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(in) :: xnmbar ! see oplus.F
!
! VT vampir tracing:
!
#ifdef VT
#include <VT.inc>
#endif
!
! Local:
      integer :: k,i,j,nlevs,ier,lat
      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)
      complex :: expt,expu
      real :: fmin,fmax
!
#ifdef VT
!     code = 116 ; state = 'addiag' ; activity='ModelCode'
      call vtbegin(116,ier)
#endif
      nlevs = lev1-lev0+1
!
! Latitude scan:
      do j=lat0,lat1
!
! Define 3d N2 (fields.F):
!
        do i=lon0,lon1
          n2(:,i,j) = (1.-o2(:,i,j)-ox(:,i,j)-he(:,i,j))
        enddo
!       write(6,"('addiag: j=',i4,' n2(:,lon0:lon1,j)=',2e12.4)")
!    |    j,minval(n2(:,lon0:lon1,j)),maxval(n2(:,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
!
! 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)+
     |         he(k,i,j)*rmassinv_he+n2(k,i,j)*rmassinv(3))
            xnmbarm(k,i,j)=p0*expz(k)*barm(k,i,j)/
     |        (boltz*tn(k,i,j))
          enddo
        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)
!
! Lower boundary of Z is in z_lbc (lbc.F):
        z(1,lon0:lon1,j) = z_lbc(lon0:lon1,j)
!
! 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
!
! Calculate diagnostic geopotential zg, with varying gravity:
      call calczg(tn,o2,ox,n2,he,z,zg,lon0,lon1,lev0,lev1,lat0,lat1)
!
#ifdef VT
!     code = 116 ; state = 'addiag' ; activity='ModelCode'
      call vtend(116,ier)
#endif
      end subroutine addiag
!-----------------------------------------------------------------------
      subroutine calczg(tn,o2,ox,n2,he,z,zg,lon0,lon1,lev0,lev1,
     |  lat0,lat1)
!
! Given geopotential z (calculated with the model constant gravity),
!   calculate geopotential zg with varying gravity. This is taken from
!   tgcmproc_f90, routines calchts and glatf in proclat.F.
! ZG will be put on secondary histories, along with the regular Z.
!
      use params_module,only: dz,glat
      use cons_module,only: boltz,avo
      use addfld_module,only: addfld
!
! 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)
     |  he, ! helium (mmr)
     |  n2, ! molecular nitrogen (mmr)
     |  z   ! geopotential calculated with constant gravity (from addiag)
      real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),
     |  intent(out) ::
     |  zg  ! output geopotential calculated with varying gravity
!
! Local:
      integer :: i,j,k
      real :: g0,r0,c2
      real,dimension(lev0:lev1) :: xmas,g
      real,parameter :: dgtr=1.74533E-2
!
! Latitude scan:
      do j=lat0,lat1
        c2 = cos(2.*dgtr*glat(j))
        g0 = 980.616*(1.-.0026373*c2)
        r0 = 2.*g0/(3.085462e-6 + 2.27e-9*c2) ! effective earth radius
!
! Longitude scan:
        do i=lon0,lon1
          g(1)=g0*(r0/(r0+0.5*(z(1,i,j)+z(2,i,j))))**2
          xmas(:) = 1./(ox(:,i,j)/16.+o2(:,i,j)/32.+he(:,i,j)/4.+
     |                  n2(:,i,j)/28.)/avo
!
! Levels:
          zg(lev0,i,j) = z(lev0,i,j)
          do k=lev0+1,lev1-1
            zg(k,i,j) = zg(k-1,i,j) + boltz*dz*tn(k-1,i,j) /
     |        (xmas(k-1)*g(k-1))
            g(k)=g0*(r0/(r0+0.5*(zg(k,i,j)+z(k+1,i,j))))**2
          enddo ! k=lev0+1,lev1-1
          zg(lev1,i,j) = 1.5*zg(lev1-1,i,j)-0.5*zg(lev1-2,i,j)
        enddo ! i=lon0,lon1
!
! Save ZG to secondary histories:
! PLEASE DO NOT COMMENT THIS OUT -- ZG IS A MANDATORY FIELD ON SECH HISTORIES
!
        call addfld('ZG','Geometric Height ZG',
     |    'cm',zg(:,lon0:lon1,j),'ilev',lev0,lev1,'lon',lon0,lon1,j)
      enddo ! j=lat0,lat1
      end subroutine calczg
