


      subroutine density(idm,jdm,tempp,salnp,den,kmt,klev,refp)
      integer idm,jdm,klev
      real tempp(idm,jdm),salnp(idm,jdm),ref_p
      integer kmt(idm,jdm)
      real den(idm,jdm)
!===================================================================
! This program is used to calculate the density of each ocean grid
! point. Now the default is to calculate the potential density of
! the sea water -- sigma-0.
! The inputs are the dimensions of i, j, k direction, temperature
! and salinity, ocean grid mask (kmt). The outputs are density
! (sigma), temperature and salionity.
!
! Coded by Aixue Hu, May 9, 2002.
!
! input:
!	idm -- i dimension of the input array
!	jdm -- j dimension of the input array
!	kmax -- k dimension of the input array
!	tempp -- input potential temperature
!	salnp -- input salinity (g/kg, or psu)
!	kmt -- ocean depth file
!	sigma -- output potential density
!	ref_p -- reference pressure (decibars)
!===================================================================

      integer i,j,k,lev
      real sigma(idm,jdm)


! --- calculate the potential density

      do j=1,jdm
      do i=1,idm
       if(kmt(i,j).ge.1.and.klev.le.kmt(i,j)) then
             call svannew(salnp(i,j),tempp(i,j),ref_p,sig0)
             den(i,j) = sig0
       else
             den(i,j) = 1.e30
       end if
      enddo
      enddo

      return
      end
