!      
      module hme_module
!
! This software is part of the NCAR TIE-GCM.  Use is governed by the 
! Open Source Academic Research License Agreement contained in the file 
! tiegcmlicense.txt.
!      
! Read HME boundary perturbations T, U, V, W, RHO and interpolate
!   to model grid and time, and calculates Z.
!
! HME tidal components: There is only one file with all the components
! included
!
      use params_module,only: nlon,nlat,nlonp4,nlonp2
      use mpi_module,only: lon0,lon1,lat0,lat1
      use nchist_module,only:nc_open,nc_close,handle_ncerr
      use addfld_module,only: addfld
      implicit none
!
#ifdef SUN
#include <netcdf3.inc>
#else
#include <netcdf.inc>
#endif 
! 
! HME boundary data at model grid, output by this module,
!   to be used by dt, duv, addiag, etc:
!
! 11/1/05: later change these from global to allocated subdomains:
!   real,dimension(:,:) :: ! will be (lon0:lon1,lat0,lat1)
!
      real,dimension(nlonp4,nlat),save ::
     |  hme_z, 
     |  hme_t, 
     |  hme_u, 
     |  hme_v, 
     |  hme_w, 
     |  hme_rho 
! Private module data, read by sub rdhme:
      integer,parameter,private :: nmonth= 1, nhour = 24
      real,dimension(nlonp4,nlat,nmonth,nhour),private,save ::
     |  w_all,
     |  rho_all,
     |  t_all,
     |  u_all,
     |  v_all 
!
      contains 
!-----------------------------------------------------------------------
      subroutine gethme
!
! Module driver to read nc files, and do time interpolations.
!
! Files provided by user via namelist read:
      use input_module,only: 
     |  hme_ncfile
      use init_module,only: istep,secs,iday,iyear
!
! Driver for obtaining HME data, called once per timestep from advance.
!
! Args:
!
! Local:
      integer :: iprint
!
      iprint = 0
      if (istep==1) iprint = 1
!
! Get hme data:
      if (istep==1) call rdhme(hme_ncfile,' ')
      call mkhme(iyear,iday,int(secs),iprint,' ')
!
      end subroutine gethme
!-----------------------------------------------------------------------
      subroutine rdhme(ncfile,type)
      use input_module,only: mxlen_filename
      implicit none
!
! Args:
      character(len=*),intent(in) :: ncfile,type
!
! Local:
      integer :: ncid,istat
      character(len=mxlen_filename) :: dskfile
      integer :: nlon_rd, nlat_rd, nmonth_rd, nhour_rd
      integer :: id_nmonth, id_nhour, id_nlon, id_nlat
      integer :: idv_w, idv_t, idv_u, idv_v, idv_rho
      character(len=240) :: char240
      real,dimension(nlonp4,nlat,nmonth,nhour) :: t,u,v,w,rho
      real,dimension(nlon,nlat,nmonth,nhour) :: tmp
!
      dskfile = ' '
      call getfile(ncfile,dskfile)
      write(6,"(/,72('-'))")
      write(6,"('Reading HME data file ',a)") trim(ncfile)
!
! Open netcdf file:
      call nc_open(ncid,dskfile,'OLD','READ')
      if (ncid==0) then
        write(6,"(/,'>>> rdhme: error opening netcdf HME ',
     |    'file ',a,' dskfile ',a)") trim(ncfile),trim(dskfile)
        call shutdown('rdhme')
      endif
!
! Check nmonth dimension:
      istat = nf_inq_dimid(ncid,'nmonth',id_nmonth)
      istat = nf_inq_dimlen(ncid,id_nmonth,nmonth_rd)
      if (istat /= NF_NOERR) then
        write(char240,"('rdhme: Error getting nmonth dimension from ',
     |    'file ',a)") trim(ncfile)
        call handle_ncerr(istat,char240)
      endif
      if (nmonth_rd /= nmonth) then
        write(6,"(/,'>>> rdhme: nmonth_rd=',i4,' not equal to nmonth=',
     |    i4)") nmonth_rd,nmonth
        write(6,"('hme data file: ',a)") trim(ncfile)
        call shutdown('rdhme')
      endif
!
! Get nhour (time) dimension:
      istat = nf_inq_dimid(ncid,'time',id_nhour)
      istat = nf_inq_dimlen(ncid,id_nhour,nhour_rd)
      if (istat /= NF_NOERR) then
        write(char240,"('rdhme: Error getting time dimension from ',
     |    'file ',a)") trim(ncfile)
        call handle_ncerr(istat,char240)
      endif
      if (nhour_rd /= nhour) then
        write(6,"(/,'>>> rdhme: nhour_rd=',i4,' not equal to nhour=',
     |    i4)") nhour_rd,nhour
        write(6,"('hme data file: ',a)") trim(ncfile)
        call shutdown('rdhme')
      endif
!
! Get nlon dimension:
      istat = nf_inq_dimid(ncid,'lon',id_nlon)
      istat = nf_inq_dimlen(ncid,id_nlon,nlon_rd)
      if (istat /= NF_NOERR) then
        write(char240,"('rdhme: Error getting nlon dimension from ',
     |    'file ',a)") trim(ncfile)
        call handle_ncerr(istat,char240)
      endif
      if (nlon_rd /= nlon) then
        write(6,"(/,'>>> rdhme: nlon_rd=',i4,' not equal to nlon=',
     |    i4)") nlon_rd,nlon
        write(6,"('hme data file: ',a)") trim(ncfile)
        call shutdown('rdhme')
      endif
!
! Get nlat dimension:
      istat = nf_inq_dimid(ncid,'lat',id_nlat)
      istat = nf_inq_dimlen(ncid,id_nlat,nlat_rd)
      if (istat /= NF_NOERR) then
        write(char240,"('rdhme: Error getting nlat dimension from ',
     |    'file ',a)") trim(ncfile)
        call handle_ncerr(istat,char240)
      endif
      if (nlat_rd /= nlat) then
        write(6,"(/,'>>> rdhme: nlat_rd=',i4,' not equal to nlat=',
     |    i4)") nlat_rd,nlat
        write(6,"('hme data file: ',a)") trim(ncfile)
        call shutdown('rdhme')
      endif
!
! Get W vertical velocity [m/s]:
!
      istat = nf_inq_varid(ncid,'WN',idv_w)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting WN var id')
      w = 0. ! init
      istat = nf_get_var_double(ncid,idv_w,tmp)
      w(3:nlon+2,:,:,:) = tmp(:,:,:,:)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting variable W')
!     write(6,"('rdhme ',a,': w min,max=',2e12.4)") 
!    |  type,minval(w),maxval(w)
!
! Get TN perturbation [deg K]:
!
      istat = nf_inq_varid(ncid,'TN',idv_t)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting TN var id')
      t = 0.
      istat = nf_get_var_double(ncid,idv_t,tmp)
      t(3:nlon+2,:,:,:) = tmp(:,:,:,:)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting variable TN')
!     write(6,"('rdhme ',a,': tn min,max=',2e12.4)") 
!    |  type,minval(t),maxval(t)
!
! Get UN perturbation [m/s]:
!
      istat = nf_inq_varid(ncid,'UN',idv_u)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting UN var id')
      u = 0.
      istat = nf_get_var_double(ncid,idv_u,tmp)
      u(3:nlon+2,:,:,:) = tmp(:,:,:,:)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting variable UN')
      u = u*100. 	! convert to cm/s
!     write(6,"('rdhme ',a,': un min,max=',2e12.4)") 
!    |  type,minval(u),maxval(u)
!
! Get VN perturbation [m/s]:
!
      istat = nf_inq_varid(ncid,'VN',idv_v)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting VN var id')
      v = 0.
      istat = nf_get_var_double(ncid,idv_v,tmp)
      v(3:nlon+2,:,:,:) = tmp(:,:,:,:)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting variable VN')
      v = v*100. 	! convert to cm/s
!     write(6,"('rdhme ',a,': vn min,max=',2e12.4)") 
!    |  type,minval(v),maxval(v)
!
! Get DN perturbation [%]:
!
      istat = nf_inq_varid(ncid,'DN',idv_rho)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting DN var id')
      rho = 0.
      istat = nf_get_var_double(ncid,idv_rho,tmp)
      rho(3:nlon+2,:,:,:) = tmp(:,:,:,:)
      if (istat /= NF_NOERR) call handle_ncerr(istat,
     |  'rdhme: Error getting variable VN')
      rho = rho*0.01 	! convert from [%] to [-]
!     write(6,"('rdhme ',a,': rho min,max=',2e12.4)") 
!    |  type,minval(rho),maxval(rho)
!
! Transfer to private module data (whole-array assignments):
      w_all = w
      rho_all = rho
      t_all = t
      u_all = u
      v_all = v
!
! Close the file:
      call nc_close(ncid)
      write(6,"('Completed read from HME data file ',a)") trim(ncfile)
      write(6,"(/,72('-'))")
      end subroutine rdhme
!-----------------------------------------------------------------------
      subroutine mkhme(iyear,iday,isecs,iprint,type)
!
! Use data read from hme file to return perturbation in Z
! and T at current model date and time. It is assumed that the hme file
! provides data for a whole year with one day of data for each month
! with hourly values (0 UT to 23 UT).
! The data is linearly interpolated first to the model UT for the
! corresponding months (current and next month) and then linearly
! interpolated to the modelday
!
! am 4/2012: I have only one month therefore assume day is the one for
!   that month w/o interpolating to the date, interpolation to hour will
!   still be done
!
      use hist_module,only: modeltime
#ifdef MPI
      use mpi_module,only: mp_periodic_f2d
#endif
      implicit none
!
! Args:
      integer,intent(in) :: iyear,iday,isecs,iprint
      character(len=*),intent(in) :: type
!
! Local:
      real :: difsec,difday,sec,secint,dayint
      integer :: i,j,ndmon_nl(13),ndmon_lp(13),nonleap,
     |   mon_cur,mon_nxt,ihr_cur,ihr_nxt,nointp=0,
     |   lonbeg,lonend,i0,i1,j0,j1,iday_tmp
      real,dimension(lon0:lon1,lat0:lat1) :: t,u,v,rho,
     |  rho_curmo, t_curmo, u_curmo, v_curmo,
     |  rho_nxtmo, t_nxtmo, u_nxtmo, v_nxtmo,
     |  rho_curmo_nxthr, t_curmo_nxthr, u_curmo_nxthr, v_curmo_nxthr,
     |  rho_nxtmo_nxthr, t_nxtmo_nxthr, u_nxtmo_nxthr, v_nxtmo_nxthr
      real :: fhme(lon0:lon1,lat0:lat1,4) ! for mp calls
!   
! Data: GSWM data given at the middle of each month 
!       -> day of year for this data assuming non-leap year
!         J  F  M  A   M   J   J   A   S   O   N   D   J
      data ndmon_nl                                ! non_leap
     |  /15,46,74,105,135,166,196,227,258,288,319,349,380/ 
! 
! External:
      real,external :: finterp
!
! For addfld calls:
      i0 = lon0 ; i1 = lon1
      j0 = lat0 ; j1 = lat1
! 
! Get month of model run
      do i = 1,13
        if(iday.le.ndmon_nl(i)) goto 10
      enddo
 10   mon_nxt = i		! next month
      if(mon_nxt == 13 ) mon_nxt = 1
      mon_cur = i-1		! current month
      if(mon_cur == 0 ) mon_cur = 12
!
! am 4/2012 for the initial testing we have only one month
! therefore mon_cur =1
      mon_cur = 1    
      mon_nxt = 1 
      iday_tmp = 15  ! this is since the subroutine thinks we are in January, and we do not want
                     ! to do the interpolation in doy   
! 
! Get hour of model run (model hours from 0 to 23 )
      ihr_cur = modeltime(2)	! current hour
      ihr_nxt = modeltime(2)+1  ! next hour
      if(ihr_nxt == 24 ) ihr_nxt = 1
!
! Subdomains at current month, current hour:
      rho_curmo = hmedat(type,mon_cur,ihr_cur+1,'rho')
      t_curmo = hmedat(type,mon_cur,ihr_cur+1,'t')
      u_curmo = hmedat(type,mon_cur,ihr_cur+1,'u')
      v_curmo = hmedat(type,mon_cur,ihr_cur+1,'v')
!     call addfld('z_curmo','z_curmo',' ',z_curmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('t_curmo','t_curmo',' ',t_curmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('u_curmo','u_curmo',' ',u_curmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('v_curmo','v_curmo',' ',v_curmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!
! Subdomains at next month, current hour:
      rho_nxtmo = hmedat(type,mon_nxt,ihr_cur+1,'rho')
      t_nxtmo = hmedat(type,mon_nxt,ihr_cur+1,'t')
      u_nxtmo = hmedat(type,mon_nxt,ihr_cur+1,'u')
      v_nxtmo = hmedat(type,mon_nxt,ihr_cur+1,'v')
!     call addfld('z_nxtmo','z_nxtmo',' ',z_nxtmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('t_nxtmo','z_nxtmo',' ',t_nxtmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('u_nxtmo','z_nxtmo',' ',u_nxtmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!     call addfld('v_nxtmo','z_nxtmo',' ',v_nxtmo,
!    |  'lon',i0,i1,'lat',j0,j1,0)
!
! Interpolate to month:
      if (isecs /= 0) then
        difsec = 60.*60. ! difference in sec between ihr_cur and ihr_nxt
        secint = float(modeltime(3)*60 + modeltime(4))  ! interpolation time [sec]
!
! Current month, next hour:
        rho = rho_curmo ; t = t_curmo ; u = u_curmo ; v = v_curmo
        rho_curmo_nxthr = hmedat(type,mon_cur,ihr_nxt+1,'rho')
        t_curmo_nxthr = hmedat(type,mon_cur,ihr_nxt+1,'t')
        u_curmo_nxthr = hmedat(type,mon_cur,ihr_nxt+1,'u')
        v_curmo_nxthr = hmedat(type,mon_cur,ihr_nxt+1,'v')
        call timeinterp(rho,rho_curmo_nxthr,difsec,secint,rho_curmo,1)
        call timeinterp(t, t_curmo_nxthr, difsec, secint, t_curmo, 1)
        call timeinterp(u, u_curmo_nxthr, difsec, secint, u_curmo, 1)
        call timeinterp(v, v_curmo_nxthr, difsec, secint, v_curmo, 1)
!
! Interpolate to next month:
        rho = rho_nxtmo ; t = t_nxtmo ; u = u_nxtmo ; v = v_nxtmo
!
! Next month, next hour:
        rho_nxtmo_nxthr = hmedat(type,mon_nxt,ihr_nxt+1,'rho')
        t_nxtmo_nxthr = hmedat(type,mon_nxt,ihr_nxt+1,'t')
        u_nxtmo_nxthr = hmedat(type,mon_nxt,ihr_nxt+1,'u')
        v_nxtmo_nxthr = hmedat(type,mon_nxt,ihr_nxt+1,'v')

        call timeinterp(rho,rho_nxtmo_nxthr,difsec,secint,rho_nxtmo,1)
        call timeinterp(t, t_nxtmo_nxthr, difsec, secint, t_nxtmo,1)
        call timeinterp(u, u_nxtmo_nxthr, difsec, secint, u_nxtmo,1)
        call timeinterp(v, v_nxtmo_nxthr, difsec, secint, v_nxtmo,1)

      endif ! interpolate to month or not
!
! Check if interpolation to ut is necessary
! am 4/2012 change iday to iday_tmp since we have only one month of data
      nointp = 0
      if(iday_tmp.eq.ndmon_nl(mon_cur)) then ! same day as cur. month
         nointp = 1
         goto 20
      endif
      if(iday_tmp.eq.ndmon_nl(mon_nxt)) then ! same day as next month
         nointp = 2
         goto 20
      endif
!
! If ut interpolation is necessary, calculate time differences
      if(mon_cur /= 12) then                           ! not December
        difday = ndmon_nl(mon_nxt)-ndmon_nl(mon_cur)   ! difference in days
        dayint = iday - ndmon_nl(mon_cur) ! difference to interpolation day
      else                                             ! December wrap around
        difday = ndmon_nl(mon_cur+1)-ndmon_nl(mon_cur) ! difference in days
        if(iday.lt.ndmon_nl(mon_nxt)) then ! difference to interpolation day
           dayint = 365. - ndmon_nl(mon_cur)+  iday
        else
          dayint = iday - ndmon_nl(mon_cur)
        endif
      endif
 20   continue          ! if no interpolation is necessary (nointp /= 0)
!
! Interpolate to ut if necessary to local rho,t,u,v:
      select case (nointp)
        case (0)                  ! interpolate
          call timeinterp(rho_curmo,rho_nxtmo,difday,dayint,rho,1)
          call timeinterp(t_curmo,t_nxtmo,difday,dayint,t,1)
          call timeinterp(u_curmo,u_nxtmo,difday,dayint,u,1)
          call timeinterp(v_curmo,v_nxtmo,difday,dayint,v,1)
        case (1)                  ! no interp (same day as current month)
          rho = rho_curmo
          t = t_curmo
          u = u_curmo
          v = v_curmo
        case (2)                  ! no interp (same day as next month)
          rho = rho_nxtmo
          t = t_nxtmo
          u = u_nxtmo
          v = v_nxtmo
        case default
          write(6,"(/,'>>> mkhme: unknown nointp=',i4)") nointp
          call shutdown('mkhme')
      end select
!
! Transfer to module data according to type:
      hme_rho(lon0:lon1,lat0:lat1) = rho
      hme_t(lon0:lon1,lat0:lat1) = t
      hme_u(lon0:lon1,lat0:lat1) = u
      hme_v(lon0:lon1,lat0:lat1) = v
!
! Do mpi periodic points exchange for gswm with f2d(:)
! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
!
#ifdef MPI
           fhme(:,:,1) = hme_rho(lon0:lon1,lat0:lat1)
           fhme(:,:,2) = hme_t(lon0:lon1,lat0:lat1)
           fhme(:,:,3) = hme_u(lon0:lon1,lat0:lat1)
           fhme(:,:,4) = hme_v(lon0:lon1,lat0:lat1)
           call mp_periodic_f2d(fhme,lon0,lon1,lat0,lat1,4)
           hme_rho(lon0:lon1,lat0:lat1) = fhme(:,:,1)
           hme_t(lon0:lon1,lat0:lat1) = fhme(:,:,2)
           hme_u(lon0:lon1,lat0:lat1) = fhme(:,:,3)
           hme_v(lon0:lon1,lat0:lat1) = fhme(:,:,4)
#else
          call set_periodic_f2d(hme_rho)
          call set_periodic_f2d(hme_t)
          call set_periodic_f2d(hme_u)
          call set_periodic_f2d(hme_v)
#endif

!       call addfld('mi_di_z','GSWM migrating diurnal Z','cm',
!    |    gswm_mi_di_z(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0)
!       call addfld('mi_di_t','GSWM migrating diurnal TN','K',
!    |    gswm_mi_di_t(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0)
!       call addfld('mi_di_u','GSWM migrating diurnal UN','cm/s',
!    |    gswm_mi_di_u(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0)
!       call addfld('mi_di_v','GSWM migrating diurnal VN','cm/s',
!    |    gswm_mi_di_v(i0:i1,j0:j1),'lon',i0,i1,'lat',j0,j1,0)
      end subroutine mkhme
!-----------------------------------------------------------------------
      function hmedat(type,mon,ihr,fname)
!
! Return subdomain array from read arrays according to requested field, 
!   month, and hour.
! The 4 global fields t, u,v,rho are private module data, read by 
!   sub rdhme at beginning of model run.
!
! Args:
      character(len=*),intent(in) :: type,fname
      integer,intent(in) :: mon,ihr
!
! Function output dimension:
      real :: hmedat(lon0:lon1,lat0:lat1)
!
! Field type must be t,u,v, or z:
      if (fname/='t'.and.fname/='u'.and.fname/='v'.and.fname/='rho')then
        write(6,"(/,'>>> hmedat: unknown fname=',a)") fname
        write(6,"('Must be rho, u, v, or t')")
        call shutdown('hmedat')
      endif
      
      if (fname=='t') hmedat(:,:)= t_all(lon0:lon1,lat0:lat1,mon,ihr)
      if (fname=='u') hmedat(:,:)= u_all(lon0:lon1,lat0:lat1,mon,ihr)
      if (fname=='v') hmedat(:,:)= v_all(lon0:lon1,lat0:lat1,mon,ihr)
      if (fname=='rho')hmedat(:,:)= rho_all(lon0:lon1,lat0:lat1,mon,ihr)
!
      end function hmedat
!-----------------------------------------------------------------------
      subroutine timeinterp(d1,d2,difd1d2,difint,fout,
     |  iprnt)
!
! Interpolate from fields d1,d2  linearly to time difint
!   (d1 must be at 0 unit time), returning fout
! On input:
!   d1,d2   = the field (d1 at 0 unit, d2 at difd1d2 units)
!   difd1d2 = time difference [unit] between d1 & d2
!   difint  = time of interpolation [unit] (counted from d1=0 time)
! In output:
!   fout is defined at difint
!
! Args:
      integer,intent(in) :: iprnt
      real,intent(in)  :: difd1d2,difint,d1(lon0:lon1,lat0:lat1),
     |   d2(lon0:lon1,lat0:lat1)
      real,intent(out) :: fout(lon0:lon1,lat0:lat1)
!
! Local:
      integer :: i,j
      real :: frac0,difd1d2_inv
!
      fout = 0. ! initialize output      
!
! Interpolation:
! from data_cur/d1 to data_nxt/d2 the difference is in time difd1d2
!      
      difd1d2_inv = 1./difd1d2
      frac0 = difint*difd1d2_inv	! x(int)/[x(d2)-x(d1)]
 
! linear interpolation: special case x(d1) = 0    
!  fout = (d2-d1)*frac0 + d1 - (d2-d1)*difd1d2_inv*x(d1)
! (d2-d1)*difd1d2_inv*x(d1) = 0 since x(d1) = 0
      
      do i = lon0,lon1
        do j = lat0,lat1
          fout(i,j) = (d2(i,j)-d1(i,j))*frac0 + d1(i,j)
	enddo
      enddo
      
      end subroutine timeinterp
!-----------------------------------------------------------------------
      subroutine set_periodic_f2d(f)
!
! Set periodic points for all f2d fields (serial or non-mpi only):
!
      real,intent(inout) :: f(nlonp4,nlat)

! lons 1,2 <- nlonp4-3,nlonp4-2 and nlonp4-1,nlonp4 <- 3,4
      f(1:2,:)              = f(nlonp4-3:nlonp4-2,:)
      f(nlonp4-1:nlonp4,:)  = f(3:4,:)
!      
      end subroutine set_periodic_f2d
!-----------------------------------------------------------------------
      end module hme_module
!-----------------------------------------------------------------------
      subroutine cal_hme_z(istep,itp,tbgrd,zprt,tprt,rhoprt)
      
      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(nlonp4,nlat),  ! tn pert HME [k]
     |    rhoprt(nlonp4,nlat) ! rho_pertr/rho_bgrd  HME [-]
      real,intent(out):: zprt(nlonp4,nlat)
      
      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_hme_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,*) 'HME1',i,j,boltz/grav/barm_lb(i,j)*tprt(i,j)
!          write(6,*) 'HME2',i,j,
!     |  	boltz/grav/barm_lb(i,j)*tbgrd(i,j)*
!     |         (rhoprt(i,j))
!          write(6,*) 'HME3',i,j,
!     |  	boltz/grav/barm_lb(i,j)*tbgrd(i,j)*
!     |         (tprt(i,j)/tbgrd(i,j))
        enddo
      enddo
!     write(6,"('cal_hme_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_hme_z