I have derived the implicit upstream radiation condition from Marchesiello et al.(2001). I found it should have an implicit assumption dx=dy (or om=on in ROMS). If not, the variables r_x (Cx in ROMS) and r_y (Cy in ROMS) should be
Code: Select all
dy*dy*dudt*dudx*alpha dx*dx*dudt*dudy*beta
r_x= - --------------------------------- r_y= - ---------------------------------
dy*dy*dudx*dudx+dx*dx*dudy*dudy dy*dy*dudx*dudx+dx*dx*dudy*dudy
x(b-1,j)-x(b-2,j) y(b-1,j)-y(b-1,j-1)
alpha= ------------------- beta=---------------------
x(b,j)-x(b-1,j) y(b,j)-y(b,j-1)
Code: Select all
dxzxm=GRID(ng)%om_r(Istr+1,j))
aphazxm=dxzxm/GRID(ng)%om_r(Istr,j)
IF ((dUdt*dUdx).lt.0.0_r8) dUdt=0.0_r8
IF ((dUdt*(grad(Istr+1,j )+ &
& grad(Istr+1,j+1))).gt.0.0_r8) THEN
dUde=grad(Istr+1,j )
dyzxm=GRID(ng)%on_p(Istr+1,j))
betazxm=dyzxm/GRID(ng)%on_p(Istr,j)
ELSE
dUde=grad(Istr+1,j+1)
dyzxm=GRID(ng)%on_p(Istr+1,j+1))
betazxm=dyzxm/GRID(ng)%on_p(Istr,j+1)
END IF
!! cff=MAX(dUdx*dUdx+dUde*dUde,eps)
!! Cx=dUdt*dUdx
cff=MAX(dyzxm*dyzxm*dUdx*dUdx+dxzxm*dxzxm*dUde*dUde,eps)
Cx=dUdt*dUdx*dyzxm*dyzxm*aphazxm
# ifdef RADIATION_2D
Ce=MIN(cff,MAX(dUdt*dUde,-cff))*dxzxm*dxzxm*betazxm
# else
Ce=0.0_r8
# endif
By the way, this is also the reason I have posted on the title "question about the caculating of om_p, om_u, on_p, on_u" in viewtopic.php?t=4805