import numpy as np

import calculate3dpressure

def calculate3ddensity(t,ps,p0,hyam,hybm) :
#---------------------
#
    Rair=287.058
#
    nlev, nlat, nlon = t.shape
#
    density  = np.zeros( ( nlev , nlat, nlon ) )
#
    pressure = calculate3dpressure.calculate3dpressure(ps=ps,p0=p0,hyam=hyam,hybm=hybm)
# 
    for ilev in range(nlev) :
#
        density[ilev,:,:] = pressure[ilev,:,:] / Rair / t[ilev,:,:]
#
    return density