FORALL (j=1:nlat,i=1:nlon+1) t31(i,j) = 0
DO k = 1,nlev1
FORALL (j=1:nlat,i=1:nlon+1) t31(i,j) = t31(i,j)+(p(i,j,k)
. +p(i1,j,k)p(i,j,k+1)p(i1,j,k+1))*(u(i1,j,k)+u(i1,j,k+1))
ENDDO
FORALL (j=1:nlat,i=1:nlon+1) t31(i,j) = 0.25*t31(i,j)
FORALL (j=1:nlat+1,i=1:nlon) t35(i,j) = 0
DO k = 1,nlev1
FORALL (j=1:nlat+1,i=1:nlon) t35(i,j) = t35(i,j)+(p(i,j,k)
. +p(i,j1,k)p(i,j,k+1)p(i,j1,k+1))*(v(i,j1,k)+v(i,j1,k+1))
ENDDO
FORALL (j=1:nlat+1,i=1:nlon) t35(i,j) = 0.25*t35(i,j)
FORALL (j=1:nlat,i=1:nlon) pst(i,j) = 1.56961230576048E7*(
. (1/dlam*(hyu(i,j)*t31(i+1,j)hyu(i1,j)*t31(i,j))+1/dtheta
. *(hxv(i,j)*t35(i,j+1)hxv(i,j1)*t35(i,j)))/(hxt(i,j)*hyt(i,j))
. )
