!-----------------------------------------------------------------------
      subroutine cal_ctmt_z(istep,itp,tbgrd,zprt,tprt,rhoprt)

! cal_ctmt_z based on Maute's cal_hme_z and called from lbc.F
!  with 2 extra arguements istep,itp in tuzvz_lbc called from advance.F
! ilev = interface points for Z,ZG,DEN,NE,UI_ExB,VI_ExB,WI_ExB
! lev(=ilev) midpoints for UN,VN,WN,TN,TI,TE,O1,O2,NO
      
      use mpi_module,only: lon0,lon1,lat0,lat1
      use params_module,only: nlonp4,nlat
      use cons_module,only: 
     |   grav,           ! 870. [cm/s2] accel due to gravity
     |   boltz,          ! 1.38E-16 [ergs/K] boltzman's constant 
     |   rmassinv,
     !   avo             ! 6.023e23 [#/mol] Avogadro constant 
      use fields_module,only: o2,o1,barm
      implicit none
!
! calculates the zpert 
!   zprt = -R/g0*tprt - R/g0*tbgrd*(pprt/p0 + tprt/tbgrd)
! with R = kb/m [cm2/s2]  kb Boltzman constant; m mass
!
      integer,intent(in):: istep,  ! step
     |          itp               ! index of previous step
      real,intent(in) :: 
     |    tbgrd(nlonp4,nlat), ! t_lbc
     |    tprt(lon0:lon1,lat0:lat1),  ! tn pert CTMT [K]
     |    rhoprt(lon0:lon1,lat0:lat1) ! rho_pertr/rho_bgrd  CTMT [-]
      real,intent(out):: zprt(lon0:lon1,lat0:lat1)
      
      integer :: i,j,k
      real :: barm_lb(lon0:lon1,lat0:lat1),
     |     barm_tmp(2)
      
!     if(istep == 1) then
! 
! barm = mean molecular weight (k+1/2): [g/mol]
! barm_lb [g]
!
        do j=lat0,lat1
          do i=lon0,lon1
            do k=1,2
              barm_tmp(k) = 1./
     |  	(o2(k,i,j,itp)*rmassinv(1)+o1(k,i,j,itp)*rmassinv(2)+
     |  	(1.-o2(k,i,j,itp)-o1(k,i,j,itp))*rmassinv(3))
            enddo
            barm_lb(i,j) = 1.5*barm_tmp(1)-0.5*barm_tmp(2)
            barm_lb(i,j) = barm_lb(i,j)/avo
          enddo
        enddo
!     else  ! not first timestep
!
! barm1 = barm(k=0) (linear extrapolation)
!
!       do j=lat0,lat1
!         do i=lon0,lon1
!           barm_lb(i,j) = 1.5*barm(1,i,j,itp)-0.5*barm(2,i,j,itp)
!           barm_lb(i,j) = barm_lb(i,j)/avo
!         enddo
!       enddo
!     endif ! end first timestep
!    
! kb/m/g*Tn  [ergs/K /g / (cm/s2) *K] = [cm]
      do j=lat0,lat1
        do i=lon0,lon1
          zprt(i,j) = - boltz/grav/barm_lb(i,j)*tprt(i,j)-
     |  	boltz/grav/barm_lb(i,j)*tbgrd(i,j)*
     |         (rhoprt(i,j)+tprt(i,j)/tbgrd(i,j))

!         write(6,"('cal_ctmt_z: istep=',i3,' j=',i3,' i=',i3,
!    |      ' barm_lb=',e12.4,' tprt=',e12.4,' tbgrd=',e12.4,
!    |      ' rhoprt=',e12.4,' zprt=',e12.4)") istep,j,i,
!    |      barm_lb(i,j),tprt(i,j),tbgrd(i,j),rhoprt(i,j),zprt(i,j)

!          write(6,*) 'CTMT1',i,j,boltz/grav/barm_lb(i,j)*tprt(i,j)
!          write(6,*) 'CTMT2',i,j,
!     |  	boltz/grav/barm_lb(i,j)*tbgrd(i,j)*
!     |         (rhoprt(i,j))
!          write(6,*) 'CTMT3',i,j,
!     |  	boltz/grav/barm_lb(i,j)*tbgrd(i,j)*
!     |         (tprt(i,j)/tbgrd(i,j))
        enddo
      enddo
!      write(6,"('cal_cmtm_z: istep=',i3,' barm_lb min,max=',2e12.4,
!     |  ' zprt min,max=',2e12.4)") istep,minval(barm_lb),
!     |  maxval(barm_lb),minval(zprt),maxval(zprt)
        
      end subroutine cal_ctmt_z