Skip to content

Commit

Permalink
Pure radiation2D is unstable, grr.
Browse files Browse the repository at this point in the history
  • Loading branch information
@kshedstrom
kshedstrom committed Sep 9, 2016
1 parent 6c2849c commit 1c00163
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 20 deletions.
47 changes: 29 additions & 18 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,9 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg)
if (Je_obc>Js_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E
if (Je_obc<Js_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_W
OBC%OBC_segment_list(l_seg)%radiation = .true.
! This is a total hack for the tangential flow! - KSH
! But it wasn't on and I still get the right answer??
OBC%OBC_segment_list(l_seg)%gradient = .true.
if (Je_obc>Js_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA
if (Je_obc<Js_obc) OBC%apply_OBC_u_flather_west = .true. ! This line will not be needed soon - AJA
elseif (trim(action_str) == 'RADIATION2D') then
Expand Down Expand Up @@ -369,8 +372,8 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg)
if (Ie_obc<Is_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_N
OBC%OBC_segment_list(l_seg)%radiation = .true.
OBC%OBC_segment_list(l_seg)%radiation2D = .true.
if (Ie_obc>Is_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not bee needed soon - AJA
if (Ie_obc<Is_obc) OBC%apply_OBC_v_flather_north = .true. ! This line will not bee needed soon - AJA
if (Ie_obc>Is_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA
if (Ie_obc<Is_obc) OBC%apply_OBC_v_flather_north = .true. ! This line will not be needed soon - AJA
elseif (trim(action_str) == 'SIMPLE') then
this_kind = OBC_SIMPLE
OBC%OBC_segment_list(l_seg)%specified = .true.
Expand Down Expand Up @@ -748,8 +751,10 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
cff = max(dhdx*dhdx + dhdy*dhdy, eps)
Cx = dhdt*dhdx
Cy = min(cff,max(dhdt*dhdy,-cff))
u_new(I,j,k) = (cff*u_old(I,j,k) + Cx*u_new(I-1,j,k) - &
max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J)) / (cff + Cx)
! Turning off radiation2D part
Cy = 0
u_new(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I-1,j,k)) - &
(max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J))) / (cff + Cx)
elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then
dhdt = u_old(I-1,j,k)-u_new(I-1,j,k) !old-new
dhdx = u_new(I-1,j,k)-u_new(I-2,j,k) !in new time backward sasha for I-1
Expand Down Expand Up @@ -786,8 +791,10 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
cff = max(dhdx*dhdx + dhdy*dhdy, eps)
Cx = dhdt*dhdx
Cy = min(cff,max(dhdt*dhdy,-cff))
u_new(I,j,k) = (cff*u_old(I,j,k) + Cx*u_new(I+1,j,k) - &
max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J)) / (cff + Cx)
! Turning off radiation2D part
Cy = 0
u_new(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I+1,j,k)) - &
(max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J))) / (cff + Cx)
elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then
dhdt = u_old(I+1,j,k)-u_new(I+1,j,k) !old-new
dhdx = u_new(I+1,j,k)-u_new(I+2,j,k) !in new time backward sasha for I+1
Expand Down Expand Up @@ -828,8 +835,8 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) &
Cx = min(cff, max(dhdt*dhdx, -cff))
Cy = dhdt*dhdy
u_new(I,j,k) = (cff*u_old(I,j,k) + Cy*u_new(I,j-1,k) - &
max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j))/(cff + Cy)
u_new(I,j,k) = ((cff*u_old(I,j,k) + Cy*u_new(I,j-1,k)) - &
(max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j)))/(cff + Cy)
endif
endif

Expand All @@ -854,8 +861,8 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) &
Cx = min(cff, max(dhdt*dhdx, -cff))
Cy = dhdt*dhdy
u_new(I,j,k) = (cff*u_old(I,j,k) + Cy*u_new(I,j+1,k) - &
max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j))/(cff + Cy)
u_new(I,j,k) = ((cff*u_old(I,j,k) + Cy*u_new(I,j+1,k)) - &
(max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j)))/(cff + Cy)
endif
endif
endif ; enddo ; enddo ; enddo
Expand All @@ -878,8 +885,10 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
cff = max(dhdx*dhdx + dhdy*dhdy, eps)
Cy = dhdt*dhdy
Cx = min(cff,max(dhdt*dhdx,-cff))
v_new(i,J,k) = (cff*v_old(i,J,k) + Cy*v_new(i,J-1,k) - &
max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J)) / (cff + Cy)
! Turning off radiation2D part
Cx = 0
v_new(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J-1,k)) - &
(max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J))) / (cff + Cy)
elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then
dhdt = v_old(i,J-1,k)-v_new(i,J-1,k) !old-new
dhdy = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1
Expand Down Expand Up @@ -916,8 +925,10 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
cff = max(dhdx*dhdx + dhdy*dhdy, eps)
Cy = dhdt*dhdy
Cx = min(cff,max(dhdt*dhdx,-cff))
v_new(i,J,k) = (cff*v_old(i,J,k) + Cy*v_new(i,J+1,k) - &
max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J)) / (cff + Cy)
! Turning off radiation2D part
Cx = 0
v_new(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J+1,k)) - &
(max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J))) / (cff + Cy)
elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then
dhdt = v_old(i,J+1,k)-v_new(i,J+1,k) !old-new
dhdy = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J+1
Expand Down Expand Up @@ -958,8 +969,8 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
Cy = 0.0
if (OBC%OBC_segment_list(OBC%OBC_segment_v(I,j))%radiation2D) &
Cy = min(cff, max(dhdt*dhdy, -cff))
v_new(i,J,k) = (cff*v_old(i,J,k) + Cx*v_new(i-1,J,k) - &
max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i,j+1))/(cff + Cx)
v_new(i,J,k) = ((cff*v_old(i,J,k) + Cx*v_new(i-1,J,k)) - &
(max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i,j+1)))/(cff + Cx)
endif
endif
if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_W) then
Expand All @@ -983,8 +994,8 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, &
Cy = 0.0
if (OBC%OBC_segment_list(OBC%OBC_segment_v(I,j))%radiation2D) &
Cy = min(cff, max(dhdt*dhdy, -cff))
v_new(i,J,k) = (cff*v_old(i,J,k) + Cx*v_new(i+1,J,k) - &
max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i+1,j))/(cff + Cx)
v_new(i,J,k) = ((cff*v_old(i,J,k) + Cx*v_new(i+1,J,k)) - &
(max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i+1,j)))/(cff + Cx)
endif
endif
endif ; enddo ; enddo ; enddo
Expand Down
4 changes: 2 additions & 2 deletions src/user/soliton_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ subroutine soliton_initialize_thickness(h, G)
call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness")

x0 = 2.0*G%len_lon/3.0
y0 = 0.5*G%len_lat
y0 = 0.0
val1 = 0.395
val2 = 0.771*(val1*val1)

Expand Down Expand Up @@ -76,7 +76,7 @@ subroutine soliton_initialize_velocity(u, v, h, G)
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

x0 = 2.0*G%len_lon/3.0
y0 = 0.5*G%len_lat
y0 = 0.0
val1 = 0.395
val2 = 0.771*(val1*val1)

Expand Down

0 comments on commit 1c00163

Please sign in to comment.