subroutine calczg(tn,o2,o1,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
      use init_module,only: glat,istep
      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)
     |  o1, ! atomic oxygen (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,n2
      real,parameter :: dgtr=1.74533E-2
!
! Latitude scan:
! 1/20/10 btf: changed float(j) to glat(j) in cos of c2 calculation
      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
          n2(:) = (1.-o2(:,i,j)-o1(:,i,j))
          xmas(:) = 1./(o1(:,i,j)/16.+o2(:,i,j)/32.+n2(:)/28.)/avo
!         write(6,"('calczg: j=',i3,' i=',i3,' xmas=',/,(6e12.4))")
!    |      j,i,xmas
!
! 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