import numpy as np

def calculate3dpressbnd(p0,ps,hyai,hybi):
#----------------------
#   
    nlev = hyai.shape[0] # the value is actually nlev+1)
    nlat, nlon = ps.shape
#
    pressbnd = np.zeros( ( nlev  , nlat, nlon ) )

    for ilev in range(nlev) :

        pressbnd[ilev,:,:] = hyai[ilev] * p0 + hybi[ilev] * ps[:,:] # boundary pressure
#
    return pressbnd