From 49a8a482d28571a653c1adb83b1dab81a41522d0 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 26 Apr 2018 21:50:45 -0800 Subject: [PATCH 01/10] Start of oblique tweaks. --- src/core/MOM_open_boundary.F90 | 139 +++++++++++++++++---------------- 1 file changed, 72 insertions(+), 67 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 35b442c816..b08be83b15 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -143,6 +143,8 @@ module MOM_open_boundary !! segment (m s-1) real, pointer, dimension(:,:,:) :: rx_normal=>NULL() !< The rx_old_u value for radiation coeff !! for normal velocity + real, pointer, dimension(:,:,:) :: ry_normal=>NULL() !< The tangential value for radiation coeff + !! for normal velocity real, pointer, dimension(:,:,:) :: nudged_normal_vel=>NULL() !< The layer velocity normal to the OB segment !! that values should be nudged towards (m s-1). real, pointer, dimension(:,:,:) :: nudged_tangential_vel=>NULL() !< The layer velocity tangential to the OB segment @@ -1411,7 +1413,24 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) do k=1,G%ke J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied - segment%rx_normal(i,J,k) = OBC%ry_normal(i,J,k) + segment%ry_normal(i,J,k) = OBC%ry_normal(i,J,k) + enddo + enddo + endif + if (segment%is_E_or_W .and. segment%oblique) then + do k=1,G%ke + I=segment%HI%IsdB + do j=segment%HI%jsd,segment%HI%jed + segment%rx_normal(I,j,k) = OBC%rx_normal(I,j,k) + segment%ry_normal(I,j,k) = OBC%ry_normal(I,j,k) + enddo + enddo + elseif (segment%is_N_or_S .and. segment%oblique) then + do k=1,G%ke + J=segment%HI%JsdB + do i=segment%HI%isd,segment%HI%ied + segment%rx_normal(i,J,k) = OBC%rx_normal(i,J,k) + segment%ry_normal(i,J,k) = OBC%ry_normal(i,J,k) enddo enddo endif @@ -1444,15 +1463,15 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) elseif (segment%oblique) 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 - if (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) > 0.0) then - dhdy = segment%grad_normal(J-1,1,k) - elseif (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) == 0.0) then - dhdy = 0.0 - else - dhdy = segment%grad_normal(J,1,k) - endif + if (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) > 0.0) then + ry_new = min(segment%grad_normal(J-1,1,k), rx_max) + elseif (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) == 0.0) then + ry_new = 0.0 + else + ry_new = min(segment%grad_normal(J,1,k), rx_max) + endif if (dhdt*dhdx < 0.0) dhdt = 0.0 - Cx = min(dhdt*dhdx,rx_max) ! default to normal radiation + Cx = dhdt*dhdx cff = max(dhdx*dhdx + dhdy*dhdy, eps) Cy = min(cff,max(dhdt*dhdy,-cff)) segment%normal_vel(I,j,k) = ((cff*u_new(I,j,k) + Cx*u_new(I-1,j,k)) - & @@ -1494,23 +1513,17 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) elseif (segment%oblique) 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 forward sasha for I+1 -! if (segment%oblique) then - if (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) > 0.0) then - dhdy = segment%grad_normal(J-1,1,k) - elseif (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) == 0.0) then - dhdy = 0.0 - else - dhdy = segment%grad_normal(J,1,k) - endif -! endif + if (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) > 0.0) then + dhdy = segment%grad_normal(J-1,1,k) + elseif (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) == 0.0) then + dhdy = 0.0 + else + dhdy = segment%grad_normal(J,1,k) + endif if (dhdt*dhdx < 0.0) dhdt = 0.0 - Cx = min(dhdt*dhdx,rx_max) ! default to normal flow only -! Cy = 0. - cff = max(dhdx*dhdx, eps) -! if (segment%oblique) then - cff = max(dhdx*dhdx + dhdy*dhdy, eps) - Cy = min(cff,max(dhdt*dhdy,-cff)) -! endif + Cx = dhdt*dhdx + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cy = min(cff,max(dhdt*dhdy,-cff)) segment%normal_vel(I,j,k) = ((cff*u_new(I,j,k) + Cx*u_new(I+1,j,k)) - & (max(Cy,0.0)*segment%grad_normal(J-1,2,k) + min(Cy,0.0)*segment%grad_normal(J,2,k))) / (cff + Cx) elseif (segment%gradient) then @@ -1538,35 +1551,30 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) dhdy = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1 ry_new = 0.0 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max) - ry_avg = (1.0-gamma_v)*segment%rx_normal(I,j,k) + gamma_v*ry_new - segment%rx_normal(i,J,k) = ry_avg + ry_avg = (1.0-gamma_v)*segment%ry_normal(I,j,k) + gamma_v*ry_new + segment%ry_normal(i,J,k) = ry_avg ! The new boundary value is interpolated between future interior ! value, v_new(J-1) and past boundary value but with barotropic ! accelerations, v_new(J). segment%normal_vel(i,J,k) = (v_new(i,J,k) + ry_avg*v_new(i,J-1,k)) / (1.0+ry_avg) ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues ! implemented as a work-around to limitations in restart capability - OBC%ry_normal(i,J,k) = segment%rx_normal(i,J,k) + OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) elseif (segment%oblique) 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 -! if (segment%oblique) then - if (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) > 0.0) then - dhdx = segment%grad_normal(I-1,1,k) - elseif (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) == 0.0) then - dhdx = 0.0 - else - dhdx = segment%grad_normal(I,1,k) - endif -! endif + segment%ry_normal(i,J,k) = ry_avg + if (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) > 0.0) then + dhdx = segment%grad_normal(I-1,1,k) + elseif (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) == 0.0) then + dhdx = 0.0 + else + dhdx = segment%grad_normal(I,1,k) + endif if (dhdt*dhdy < 0.0) dhdt = 0.0 - Cy = min(dhdt*dhdy,rx_max) ! default to normal flow only -! Cx = 0 - cff = max(dhdy*dhdy, eps) -! if (segment%oblique) then - cff = max(dhdx*dhdx + dhdy*dhdy, eps) - Cx = min(cff,max(dhdt*dhdx,-cff)) -! endif + Cy = dhdt*dhdy + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = min(cff,max(dhdt*dhdx,-cff)) segment%normal_vel(i,J,k) = ((cff*v_new(i,J,k) + Cy*v_new(i,J-1,k)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) elseif (segment%gradient) then @@ -1595,35 +1603,29 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) dhdy = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J-1 ry_new = 0.0 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max) - ry_avg = (1.0-gamma_v)*segment%rx_normal(I,j,k) + gamma_v*ry_new - segment%rx_normal(i,J,k) = ry_avg + ry_avg = (1.0-gamma_v)*segment%ry_normal(I,j,k) + gamma_v*ry_new + segment%ry_normal(i,J,k) = ry_avg ! The new boundary value is interpolated between future interior ! value, v_new(J+1) and past boundary value but with barotropic ! accelerations, v_new(J). segment%normal_vel(i,J,k) = (v_new(i,J,k) + ry_avg*v_new(i,J+1,k)) / (1.0+ry_avg) ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues ! implemented as a work-around to limitations in restart capability - OBC%ry_normal(i,J,k) = segment%rx_normal(i,J,k) + OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) elseif (segment%oblique) 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 -! if (segment%oblique) then - if (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) > 0.0) then - dhdx = segment%grad_normal(I-1,1,k) - elseif (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) == 0.0) then - dhdx = 0.0 - else - dhdx = segment%grad_normal(I,1,k) - endif -! endif + if (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) > 0.0) then + dhdx = segment%grad_normal(I-1,1,k) + elseif (dhdt*(segment%grad_normal(I,1,k) + segment%grad_normal(I-1,1,k)) == 0.0) then + dhdx = 0.0 + else + dhdx = segment%grad_normal(I,1,k) + endif if (dhdt*dhdy < 0.0) dhdt = 0.0 - Cy = min(dhdt*dhdy,rx_max) ! default to normal flow only -! Cx = 0 - cff = max(dhdy*dhdy, eps) -! if (segment%oblique) then - cff = max(dhdx*dhdx + dhdy*dhdy, eps) - Cx = min(cff,max(dhdt*dhdx,-cff)) -! endif + Cy = dhdt*dhdy + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = min(cff,max(dhdt*dhdx,-cff)) segment%normal_vel(i,J,k) = ((cff*v_new(i,J,k) + Cy*v_new(i,J+1,k)) - & (max(Cx,0.0)*segment%grad_normal(I-1,2,k) + min(Cx,0.0)*segment%grad_normal(I,2,k))) / (cff + Cy) elseif (segment%gradient) then @@ -1903,7 +1905,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%eta(IsdB:IedB,jsd:jed)); segment%eta(:,:)=0.0 allocate(segment%normal_trans_bt(IsdB:IedB,jsd:jed)); segment%normal_trans_bt(:,:)=0.0 if (segment%radiation) then - allocate(segment%rx_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%rx_normal(:,:,:)=0.0 + allocate(segment%rx_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%rx_normal(:,:,:)=0.0 endif allocate(segment%normal_vel(IsdB:IedB,jsd:jed,OBC%ke)); segment%normal_vel(:,:,:)=0.0 allocate(segment%normal_vel_bt(IsdB:IedB,jsd:jed)); segment%normal_vel_bt(:,:)=0.0 @@ -1913,7 +1915,8 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif if (segment%oblique) then - allocate(segment%grad_normal(JsdB:JedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 + allocate(segment%grad_normal(JsdB:JedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 + allocate(segment%ry_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%ry_normal(:,:,:)=0.0 endif endif @@ -1925,7 +1928,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%eta(isd:ied,JsdB:JedB)); segment%eta(:,:)=0.0 allocate(segment%normal_trans_bt(isd:ied,JsdB:JedB)); segment%normal_trans_bt(:,:)=0.0 if (segment%radiation) then - allocate(segment%rx_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%rx_normal(:,:,:)=0.0 + allocate(segment%ry_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%ry_normal(:,:,:)=0.0 endif allocate(segment%normal_vel(isd:ied,JsdB:JedB,OBC%ke)); segment%normal_vel(:,:,:)=0.0 allocate(segment%normal_vel_bt(isd:ied,JsdB:JedB)); segment%normal_vel_bt(:,:)=0.0 @@ -1935,7 +1938,8 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif if (segment%oblique) then - allocate(segment%grad_normal(IsdB:IedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 + allocate(segment%grad_normal(IsdB:IedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 + allocate(segment%rx_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%rx_normal(:,:,:)=0.0 endif endif @@ -1956,6 +1960,7 @@ subroutine deallocate_OBC_segment_data(OBC, segment) if (associated (segment%eta)) deallocate(segment%eta) if (associated (segment%normal_trans_bt)) deallocate(segment%normal_trans_bt) if (associated (segment%rx_normal)) deallocate(segment%rx_normal) + if (associated (segment%ry_normal)) deallocate(segment%ry_normal) if (associated (segment%normal_vel)) deallocate(segment%normal_vel) if (associated (segment%normal_vel_bt)) deallocate(segment%normal_vel_bt) if (associated (segment%normal_trans)) deallocate(segment%normal_trans) From 8ca33b0eef418d705a5657bffd8dc4ebca83009f Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 27 Apr 2018 11:20:49 -0800 Subject: [PATCH 02/10] Fix allocate bug on road to time-filtering OBLIQUE. --- src/core/MOM_open_boundary.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index b08be83b15..2317fbc146 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1916,6 +1916,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) endif if (segment%oblique) then allocate(segment%grad_normal(JsdB:JedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 + allocate(segment%rx_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%rx_normal(:,:,:)=0.0 allocate(segment%ry_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%ry_normal(:,:,:)=0.0 endif endif @@ -1940,6 +1941,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%oblique) then allocate(segment%grad_normal(IsdB:IedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 allocate(segment%rx_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%rx_normal(:,:,:)=0.0 + allocate(segment%ry_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%ry_normal(:,:,:)=0.0 endif endif @@ -3024,7 +3026,7 @@ subroutine open_boundary_register_restarts(HI, GV, OBC_CS,restart_CSp) ! This implementation uses 3D arrays solely for restarts. We need ! to be able to add 2D ( x,z or y,z ) data to restarts to avoid using ! so much memory and disk space. *** - if (OBC_CS%radiation_BCs_exist_globally) then + if (OBC_CS%radiation_BCs_exist_globally .or. OBC_CS%oblique_BCs_exist_globally) then allocate(OBC_CS%rx_normal(HI%isdB:HI%iedB,HI%jsd:HI%jed,GV%ke)) OBC_CS%rx_normal(:,:,:) = 0.0 vd = var_desc("rx_normal","m s-1", "Normal Phase Speed for EW OBCs",'u','L') From d6bf3e13d676ba39cd26b74fdcd3e16de9957379 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 1 May 2018 06:54:23 -0800 Subject: [PATCH 03/10] Added OBC_COMPUTED_VORTICITY option. - Helps with dumbbell cases, uses external tangential flow. --- src/core/MOM_CoriolisAdv.F90 | 14 ++++++ src/core/MOM_open_boundary.F90 | 68 ++++++++++++++++++++++++------ src/user/Kelvin_initialization.F90 | 5 +++ 3 files changed, 75 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index c32352acb1..44449bbd2d 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -290,6 +290,13 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) if (OBC%freeslip_vorticity) then ; do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB dudy(I,J) = 0. enddo ; endif + if (OBC%computed_vorticity) then ; do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB + if (OBC%segment(n)%direction == OBC_DIRECTION_N) then + dudy(I,J) = 2.0*(OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))*G%dxCu(I,j) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_S) + dudy(I,J) = 2.0*(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%dxCu(I,j+1) + endif + enddo ; endif ! Project thicknesses across OBC points with a no-gradient condition. do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied) @@ -316,6 +323,13 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) if (OBC%freeslip_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB dvdx(I,J) = 0. enddo ; endif + if (OBC%computed_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB + if (OBC%segment(n)%direction == OBC_DIRECTION_E) then + dvdx(I,J) = 2.0*(OBC%segment(n)%tangential_vel(I,J,k) - v(i,J,k))*G%dyCv(i,J) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_W) + dvdx(I,J) = 2.0*(v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%dyCv(i+1,J) + endif + enddo ; endif ! Project thicknesses across OBC points with a no-gradient condition. do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 2317fbc146..b017aadd55 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -107,7 +107,9 @@ module MOM_open_boundary !! If False, a gradient condition is applied. logical :: oblique !< Oblique waves supported at radiation boundary. logical :: nudged !< Optional supplement to radiation boundary. - logical :: specified !< Boundary fixed to external value. + logical :: nudged_tan !< Optional supplement to nudge tangential velocity. + logical :: specified !< Boundary normal velocity fixed to external value. + logical :: specified_tan !< Boundary tangential velocity fixed to external value. logical :: open !< Boundary is open for continuity solver. logical :: gradient !< Zero gradient at boundary. logical :: values_needed !< Whether or not external OBC fields are needed. @@ -132,6 +134,8 @@ module MOM_open_boundary real, pointer, dimension(:,:,:) :: h=>NULL() !< The cell thickness (m) at OBC-points. real, pointer, dimension(:,:,:) :: normal_vel=>NULL() !< The layer velocity normal to the OB !! segment (m s-1). + real, pointer, dimension(:,:,:) :: tangential_vel=>NULL() !< The layer velocity tangential to the + !! OB segment (m s-1). real, pointer, dimension(:,:,:) :: normal_trans=>NULL() !< The layer transport normal to the OB !! segment (m3 s-1). real, pointer, dimension(:,:) :: normal_vel_bt=>NULL() !< The barotropic velocity normal to @@ -188,6 +192,8 @@ module MOM_open_boundary logical :: zero_vorticity = .false. !< If True, sets relative vorticity to zero on open boundaries. logical :: freeslip_vorticity = .false. !< If True, sets normal gradient of tangential velocity to zero !! in the relative vorticity on open boundaries. + logical :: computed_vorticity = .false. !< If True, uses external data for tangential velocity + !! in the relative vorticity on open boundaries. logical :: zero_strain = .false. !< If True, sets strain to zero on open boundaries. logical :: freeslip_strain = .false. !< If True, sets normal gradient of tangential velocity to zero !! in the strain on open boundaries. @@ -299,9 +305,14 @@ subroutine open_boundary_config(G, param_file, OBC) "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the relative vorticity on open boundaries. This cannot\n"// & "be true if OBC_ZERO_VORTICITY is True.", default=.false.) - if (OBC%zero_vorticity .and. OBC%freeslip_vorticity) call MOM_error(FATAL, & - "MOM_open_boundary.F90, open_boundary_config: "//& - "Only one of OBC_ZERO_VORTICITY and OBC_FREESLIP_VORTICITY can be True at once.") + call get_param(param_file, mdl, "OBC_COMPUTED_VORTICITY", OBC%computed_vorticity, & + "If true, uses the external values of tangential velocity\n"// & + "in the relative vorticity on open boundaries. This cannot\n"// & + "be true if OBC_ZERO_VORTICITY or OBC_FREESLIP_VORTICITY is True.", default=.false.) + if (OBC%zero_vorticity .and. OBC%freeslip_vorticity .and. OBC%computed_vorticity) & + call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& + "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY and OBC_COMPUTED_VORTICITY\n"//& + "can be True at once.") call get_param(param_file, mdl, "OBC_ZERO_STRAIN", OBC%zero_strain, & "If true, sets the strain used in the stress tensor to zero on open boundaries.", & default=.false.) @@ -344,7 +355,9 @@ subroutine open_boundary_config(G, param_file, OBC) OBC%segment(l)%radiation = .false. OBC%segment(l)%oblique = .false. OBC%segment(l)%nudged = .false. + OBC%segment(l)%nudged_tan = .false. OBC%segment(l)%specified = .false. + OBC%segment(l)%specified_tan = .false. OBC%segment(l)%open = .false. OBC%segment(l)%gradient = .false. OBC%segment(l)%values_needed = .false. @@ -770,6 +783,8 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg, PF) elseif (trim(action_str(a_loop)) == 'SIMPLE') then OBC%segment(l_seg)%specified = .true. OBC%specified_u_BCs_exist_globally = .true. ! This avoids deallocation + elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then + OBC%segment(l_seg)%specified_tan = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") @@ -871,6 +886,8 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg, PF) elseif (trim(action_str(a_loop)) == 'SIMPLE') then OBC%segment(l_seg)%specified = .true. OBC%specified_v_BCs_exist_globally = .true. ! This avoids deallocation + elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then + OBC%segment(l_seg)%specified_tan = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") @@ -1464,11 +1481,11 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) 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 if (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) > 0.0) then - ry_new = min(segment%grad_normal(J-1,1,k), rx_max) + ry_new = segment%grad_normal(J-1,1,k) elseif (dhdt*(segment%grad_normal(J,1,k) + segment%grad_normal(J-1,1,k)) == 0.0) then ry_new = 0.0 else - ry_new = min(segment%grad_normal(J,1,k), rx_max) + ry_new = segment%grad_normal(J,1,k) endif if (dhdt*dhdx < 0.0) dhdt = 0.0 Cx = dhdt*dhdx @@ -1912,6 +1929,11 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%normal_trans(IsdB:IedB,jsd:jed,OBC%ke)); segment%normal_trans(:,:,:)=0.0 if (segment%nudged) then allocate(segment%nudged_normal_vel(IsdB:IedB,jsd:jed,OBC%ke)); segment%nudged_normal_vel(:,:,:)=0.0 + endif + if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan) then + allocate(segment%tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_vel(:,:,:)=0.0 + endif + if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif if (segment%oblique) then @@ -1936,6 +1958,11 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%normal_trans(isd:ied,JsdB:JedB,OBC%ke)); segment%normal_trans(:,:,:)=0.0 if (segment%nudged) then allocate(segment%nudged_normal_vel(isd:ied,JsdB:JedB,OBC%ke)); segment%nudged_normal_vel(:,:,:)=0.0 + endif + if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan) then + allocate(segment%tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_vel(:,:,:)=0.0 + endif + if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif if (segment%oblique) then @@ -1967,6 +1994,7 @@ subroutine deallocate_OBC_segment_data(OBC, segment) if (associated (segment%normal_vel_bt)) deallocate(segment%normal_vel_bt) if (associated (segment%normal_trans)) deallocate(segment%normal_trans) if (associated (segment%nudged_normal_vel)) deallocate(segment%nudged_normal_vel) + if (associated (segment%tangential_vel)) deallocate(segment%tangential_vel) if (associated (segment%nudged_tangential_vel)) deallocate(segment%nudged_tangential_vel) if (associated (segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg) @@ -2288,8 +2316,8 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (segment%direction == OBC_DIRECTION_E) ishift=0 I=is_obc if (segment%field(m)%name == 'V') then - ! Only do q points within the segment - do J=js_obc+1,je_obc-1 + ! Do q points for the whole segment + do J=js_obc,je_obc ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer @@ -2331,8 +2359,8 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (segment%direction == OBC_DIRECTION_N) jshift=0 J=js_obc if (segment%field(m)%name == 'U') then - ! Only do q points within the segment - do I=is_obc+1,ie_obc-1 + ! Do q points for the whole segment + do I=is_obc,ie_obc segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer if (G%mask2dCv(i,J)>0. .and. G%mask2dCv(i+1,J)>0.) then ! Using the h remapping approach @@ -2342,13 +2370,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & G%ke, h_stack, segment%field(m)%buffer_dst(I,J,:)) - else if (G%mask2dCv(i,J)>0.) then + elseif (G%mask2dCv(i,J)>0.) then h_stack(:) = h(i,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & G%ke, h_stack, segment%field(m)%buffer_dst(I,J,:)) - else if (G%mask2dCv(i+1,J)>0.) then + elseif (G%mask2dCv(i+1,J)>0.) then h_stack(:) = h(i+1,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & @@ -2434,6 +2462,22 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) G%dxCv(i,J)) if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,J,:) = segment%normal_vel(i,J,:) enddo + elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_E_or_W .and. associated(segment%tangential_vel)) then + I=is_obc + do J=js_obc+1,je_obc-1 + do k=1,G%ke + segment%tangential_vel(I,J,k) = segment%field(m)%buffer_dst(I,J,k) + enddo + if (associated(segment%nudged_tangential_vel)) segment%nudged_tangential_vel(I,J,:) = segment%tangential_vel(I,J,:) + enddo + elseif (trim(segment%field(m)%name) == 'U' .and. segment%is_N_or_S .and. associated(segment%tangential_vel)) then + J=js_obc + do I=is_obc+1,ie_obc-1 + do k=1,G%ke + segment%tangential_vel(I,J,k) = segment%field(m)%buffer_dst(I,J,k) + enddo + if (associated(segment%nudged_tangential_vel)) segment%nudged_tangential_vel(I,J,:) = segment%tangential_vel(I,J,:) + enddo endif endif endif diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 3b249864e4..63d61bea35 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -211,6 +211,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time) if (segment%direction == OBC_DIRECTION_W) then IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB jsd = segment%HI%jsd ; jed = segment%HI%jed + JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB do j=jsd,jed ; do I=IsdB,IedB x1 = 1000. * G%geoLonCu(I,j) y1 = 1000. * G%geoLatCu(I,j) @@ -242,6 +243,10 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time) endif endif enddo ; enddo +! if (allocated(segment%tangential_vel)) then +! do J=JsdB,JedB ; do I=IsdB,IedB +! enddo ; enddo +! endif else isd = segment%HI%isd ; ied = segment%HI%ied JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB From 17c4ffa8f967c210d992e3991e001f0d11f90ac1 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Tue, 1 May 2018 11:30:59 -0800 Subject: [PATCH 04/10] Added OBC_COMPUTED_STRAIN option. - slight improvement to dumbbell. --- src/core/MOM_open_boundary.F90 | 21 ++++++++++++++----- .../lateral/MOM_hor_visc.F90 | 14 ++++++++++++- 2 files changed, 29 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index b017aadd55..1a4a60024c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -197,6 +197,8 @@ module MOM_open_boundary logical :: zero_strain = .false. !< If True, sets strain to zero on open boundaries. logical :: freeslip_strain = .false. !< If True, sets normal gradient of tangential velocity to zero !! in the strain on open boundaries. + logical :: computed_strain = .false. !< If True, uses external data for tangential velocity to compute + !! normal gradient in the strain on open boundaries. logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for !! use in the biharmonic viscosity term. logical :: brushcutter_mode = .false. !< If True, read data on supergrid. @@ -309,7 +311,9 @@ subroutine open_boundary_config(G, param_file, OBC) "If true, uses the external values of tangential velocity\n"// & "in the relative vorticity on open boundaries. This cannot\n"// & "be true if OBC_ZERO_VORTICITY or OBC_FREESLIP_VORTICITY is True.", default=.false.) - if (OBC%zero_vorticity .and. OBC%freeslip_vorticity .and. OBC%computed_vorticity) & + if ((OBC%zero_vorticity .and. OBC%freeslip_vorticity) .or. & + (OBC%zero_vorticity .and. OBC%computed_vorticity) .or. & + (OBC%freeslip_vorticity .and. OBC%computed_vorticity)) & call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY and OBC_COMPUTED_VORTICITY\n"//& "can be True at once.") @@ -319,10 +323,17 @@ subroutine open_boundary_config(G, param_file, OBC) call get_param(param_file, mdl, "OBC_FREESLIP_STRAIN", OBC%freeslip_strain, & "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & - "be true if OBC_ZERO_STRAIN is True.", default=.false.) - if (OBC%zero_strain .and. OBC%freeslip_strain) call MOM_error(FATAL, & - "MOM_open_boundary.F90, open_boundary_config: "//& - "Only one of OBC_ZERO_STRAIN and OBC_FREESLIP_STRAIN can be True at once.") + "be true if OBC_ZERO_STRAIN or OBC_COMPUTED_STRAIN is True.", default=.false.) + call get_param(param_file, mdl, "OBC_COMPUTED_STRAIN", OBC%computed_strain, & + "If true, sets the normal gradient of tangential velocity to\n"// & + "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & + "be true if OBC_ZERO_STRAIN or OBC_FREESLIP_STRAIN is True.", default=.false.) + if ((OBC%zero_strain .and. OBC%freeslip_strain) .or. & + (OBC%zero_strain .and. OBC%computed_strain) .or. & + (OBC%freeslip_strain .and. OBC%computed_strain)) & + call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& + "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN and OBC_COMPUTED_STRAIN\n"//& + "can be True at once.") call get_param(param_file, mdl, "OBC_ZERO_BIHARMONIC", OBC%zero_biharmonic, & "If true, zeros the Laplacian of flow on open boundaries in the biharmonic\n"//& "viscosity term.", default=.false.) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 0e02cefba2..85bb438f8f 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -443,13 +443,19 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, ! thicknesses on open boundaries. if (apply_OBC) then ; do n=1,OBC%number_of_segments J = OBC%segment(n)%HI%JsdB ; I = OBC%segment(n)%HI%IsdB - if (OBC%zero_strain .or. OBC%freeslip_strain) then + if (OBC%zero_strain .or. OBC%freeslip_strain .or. OBC%computed_strain) then if (OBC%segment(n)%is_N_or_S .and. (J >= js-2) .and. (J <= Jeq+1)) then do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB if (OBC%zero_strain) then dvdx(I,J) = 0. ; dudy(I,J) = 0. elseif (OBC%freeslip_strain) then dudy(I,J) = 0. + elseif (OBC%computed_strain) then + if (OBC%segment(n)%direction == OBC_DIRECTION_N) then + dudy(I,J) = CS%DX_dyBu(I,J)*(OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))*G%IdxCu(I,j) + else + dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdxCu(I,j+1) + endif endif enddo elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-2) .and. (I <= Ieq+1)) then @@ -458,6 +464,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, dvdx(I,J) = 0. ; dudy(I,J) = 0. elseif (OBC%freeslip_strain) then dvdx(I,J) = 0. + elseif (OBC%computed_strain) then + if (OBC%segment(n)%direction == OBC_DIRECTION_E) then + dvdx(I,J) = 2.0*CS%DY_dxBu(I,J)*(OBC%segment(n)%tangential_vel(I,J,k) - v(i,J,k))*G%IdyCv(i,J) + else + dvdx(I,J) = 2.0*CS%DY_dxBu(I,J)*(v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdyCv(i+1,J) + endif endif enddo endif From 9ca0b159b9c90a3dc29237517b1cc55b03923098 Mon Sep 17 00:00:00 2001 From: matthew harrison Date: Wed, 2 May 2018 13:27:45 -0400 Subject: [PATCH 05/10] ensure that new_sponges logical is set properly --- .../vertical/MOM_ALE_sponge.F90 | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index e42601d51b..9eac0285e8 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -21,7 +21,7 @@ module MOM_ALE_sponge use MOM_time_manager, only : time_type, init_external_field, get_external_field_size, time_interp_external_init use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer - +use mpp_mod, only : mpp_pe ! GMM - Planned extension: Support for time varying sponge targets. implicit none ; private @@ -166,6 +166,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ "PPM reconstruction will also be used within boundary cells.", & default=.false., do_not_log=.true.) + CS%new_sponges = .false. CS%nz = G%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed @@ -178,13 +179,10 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ CS%num_col = CS%num_col + 1 enddo ; enddo - if (CS%num_col > 0) then - allocate(CS%Iresttime_col(CS%num_col)) ; CS%Iresttime_col = 0.0 allocate(CS%col_i(CS%num_col)) ; CS%col_i = 0 allocate(CS%col_j(CS%num_col)) ; CS%col_j = 0 - ! pass indices, restoring time to the CS structure col = 1 do j=G%jsc,G%jec ; do i=G%isc,G%iec @@ -194,16 +192,12 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ col = col +1 endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data CS%nz_data = nz_data allocate(CS%Ref_h%p(CS%nz_data,CS%num_col)) do col=1,CS%num_col ; do K=1,CS%nz_data CS%Ref_h%p(K,col) = data_h(CS%col_i(col),CS%col_j(col),K) enddo; enddo - CS%new_sponges = .false. - - endif total_sponge_cols = CS%num_col @@ -354,6 +348,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) "PPM reconstruction will also be used within boundary cells.", & default=.false., do_not_log=.true.) + CS%new_sponges = .true. CS%nz = G%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed @@ -368,11 +363,9 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) if (CS%num_col > 0) then - allocate(CS%Iresttime_col(CS%num_col)) ; CS%Iresttime_col = 0.0 allocate(CS%col_i(CS%num_col)) ; CS%col_i = 0 allocate(CS%col_j(CS%num_col)) ; CS%col_j = 0 - ! pass indices, restoring time to the CS structure col = 1 do j=G%jsc,G%jec ; do i=G%isc,G%iec @@ -382,9 +375,6 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) col = col +1 endif enddo ; enddo - - CS%new_sponges = .true. - endif total_sponge_cols = CS%num_col @@ -821,7 +811,8 @@ subroutine apply_ALE_sponge(h, dt, G, CS, Time) deallocate(sp_val, mask_z) enddo - else + else + print *,'CS%nz_data= ',mpp_pe(),CS%nz_data nz_data = CS%nz_data endif From 968597521a676ad8cf8f83109e74fe6a6b8c089a Mon Sep 17 00:00:00 2001 From: matthew harrison Date: Wed, 2 May 2018 13:34:58 -0400 Subject: [PATCH 06/10] cleanup --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 9eac0285e8..99ca06bb66 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -21,7 +21,6 @@ module MOM_ALE_sponge use MOM_time_manager, only : time_type, init_external_field, get_external_field_size, time_interp_external_init use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer -use mpp_mod, only : mpp_pe ! GMM - Planned extension: Support for time varying sponge targets. implicit none ; private @@ -811,8 +810,7 @@ subroutine apply_ALE_sponge(h, dt, G, CS, Time) deallocate(sp_val, mask_z) enddo - else - print *,'CS%nz_data= ',mpp_pe(),CS%nz_data + else nz_data = CS%nz_data endif From 75b6ee0369f521a625836938f5a1886b79279a5a Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 2 May 2018 09:38:39 -0800 Subject: [PATCH 07/10] Fix indexing? --- src/core/MOM_open_boundary.F90 | 73 ++++++++++++++++++++++++---------- 1 file changed, 51 insertions(+), 22 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 1a4a60024c..7f9e1ecaab 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -136,6 +136,8 @@ module MOM_open_boundary !! segment (m s-1). real, pointer, dimension(:,:,:) :: tangential_vel=>NULL() !< The layer velocity tangential to the !! OB segment (m s-1). + real, pointer, dimension(:,:,:) :: tangential_grad=>NULL() !< The gradient of the velocity tangential + !! to the OB segment (m s-1). real, pointer, dimension(:,:,:) :: normal_trans=>NULL() !< The layer transport normal to the OB !! segment (m3 s-1). real, pointer, dimension(:,:) :: normal_vel_bt=>NULL() !< The barotropic velocity normal to @@ -194,11 +196,15 @@ module MOM_open_boundary !! in the relative vorticity on open boundaries. logical :: computed_vorticity = .false. !< If True, uses external data for tangential velocity !! in the relative vorticity on open boundaries. + logical :: imported_vorticity = .false. !< If True, uses external data for tangential velocity + !! gradients in the relative vorticity on open boundaries. logical :: zero_strain = .false. !< If True, sets strain to zero on open boundaries. logical :: freeslip_strain = .false. !< If True, sets normal gradient of tangential velocity to zero !! in the strain on open boundaries. logical :: computed_strain = .false. !< If True, uses external data for tangential velocity to compute !! normal gradient in the strain on open boundaries. + logical :: imported_strain = .false. !< If True, uses external data for tangential velocity gradients + !! to compute strain on open boundaries. logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for !! use in the biharmonic viscosity term. logical :: brushcutter_mode = .false. !< If True, read data on supergrid. @@ -306,34 +312,48 @@ subroutine open_boundary_config(G, param_file, OBC) call get_param(param_file, mdl, "OBC_FREESLIP_VORTICITY", OBC%freeslip_vorticity, & "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the relative vorticity on open boundaries. This cannot\n"// & - "be true if OBC_ZERO_VORTICITY is True.", default=.false.) + "be true if another OBC_XXX_VORTICITY option is True.", default=.false.) call get_param(param_file, mdl, "OBC_COMPUTED_VORTICITY", OBC%computed_vorticity, & "If true, uses the external values of tangential velocity\n"// & "in the relative vorticity on open boundaries. This cannot\n"// & - "be true if OBC_ZERO_VORTICITY or OBC_FREESLIP_VORTICITY is True.", default=.false.) + "be true if another OBC_XXX_VORTICITY option is True.", default=.false.) + call get_param(param_file, mdl, "OBC_IMPORTED_VORTICITY", OBC%imported_vorticity, & + "If true, uses the external values of tangential velocity\n"// & + "in the relative vorticity on open boundaries. This cannot\n"// & + "be true if another OBC_XXX_VORTICITY option is True.", default=.false.) if ((OBC%zero_vorticity .and. OBC%freeslip_vorticity) .or. & (OBC%zero_vorticity .and. OBC%computed_vorticity) .or. & - (OBC%freeslip_vorticity .and. OBC%computed_vorticity)) & + (OBC%zero_vorticity .and. OBC%imported_vorticity) .or. & + (OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or. & + (OBC%freeslip_vorticity .and. OBC%imported_vorticity) .or. & + (OBC%computed_vorticity .and. OBC%imported_vorticity)) & call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& - "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY and OBC_COMPUTED_VORTICITY\n"//& - "can be True at once.") + "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//& + "and OBC_IMPORTED_VORTICITY can be True at once.") call get_param(param_file, mdl, "OBC_ZERO_STRAIN", OBC%zero_strain, & "If true, sets the strain used in the stress tensor to zero on open boundaries.", & default=.false.) call get_param(param_file, mdl, "OBC_FREESLIP_STRAIN", OBC%freeslip_strain, & "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & - "be true if OBC_ZERO_STRAIN or OBC_COMPUTED_STRAIN is True.", default=.false.) + "be true if another OBC_XXX_STRAIN option is True.", default=.false.) call get_param(param_file, mdl, "OBC_COMPUTED_STRAIN", OBC%computed_strain, & "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & - "be true if OBC_ZERO_STRAIN or OBC_FREESLIP_STRAIN is True.", default=.false.) + "be true if another OBC_XXX_STRAIN option is True.", default=.false.) + call get_param(param_file, mdl, "OBC_IMPORTED_STRAIN", OBC%imported_strain, & + "If true, sets the normal gradient of tangential velocity to\n"// & + "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & + "be true if another OBC_XXX_STRAIN option is True.", default=.false.) if ((OBC%zero_strain .and. OBC%freeslip_strain) .or. & (OBC%zero_strain .and. OBC%computed_strain) .or. & - (OBC%freeslip_strain .and. OBC%computed_strain)) & + (OBC%zero_strain .and. OBC%imported_strain) .or. & + (OBC%freeslip_strain .and. OBC%computed_strain) .or. & + (OBC%freeslip_strain .and. OBC%imported_strain) .or. & + (OBC%computed_strain .and. OBC%imported_strain)) & call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& - "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN and OBC_COMPUTED_STRAIN\n"//& - "can be True at once.") + "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN\n"//& + "and OBC_IMPORTED_STRAIN can be True at once.") call get_param(param_file, mdl, "OBC_ZERO_BIHARMONIC", OBC%zero_biharmonic, & "If true, zeros the Laplacian of flow on open boundaries in the biharmonic\n"//& "viscosity term.", default=.false.) @@ -1941,12 +1961,16 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged) then allocate(segment%nudged_normal_vel(IsdB:IedB,jsd:jed,OBC%ke)); segment%nudged_normal_vel(:,:,:)=0.0 endif - if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan) then + if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan .or. & + OBC%computed_strain) then allocate(segment%tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_vel(:,:,:)=0.0 endif if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif + if (OBC%imported_vorticity .or. OBC%imported_strain) then + allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 + endif if (segment%oblique) then allocate(segment%grad_normal(JsdB:JedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 allocate(segment%rx_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%rx_normal(:,:,:)=0.0 @@ -1970,12 +1994,16 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged) then allocate(segment%nudged_normal_vel(isd:ied,JsdB:JedB,OBC%ke)); segment%nudged_normal_vel(:,:,:)=0.0 endif - if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan) then + if (OBC%computed_vorticity .or. segment%nudged_tan .or. segment%specified_tan .or. & + OBC%computed_strain) then allocate(segment%tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_vel(:,:,:)=0.0 endif if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif + if (OBC%imported_vorticity .or. OBC%imported_strain) then + allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 + endif if (segment%oblique) then allocate(segment%grad_normal(IsdB:IedB,2,OBC%ke)); segment%grad_normal(:,:,:) = 0.0 allocate(segment%rx_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%rx_normal(:,:,:)=0.0 @@ -2007,6 +2035,7 @@ subroutine deallocate_OBC_segment_data(OBC, segment) if (associated (segment%nudged_normal_vel)) deallocate(segment%nudged_normal_vel) if (associated (segment%tangential_vel)) deallocate(segment%tangential_vel) if (associated (segment%nudged_tangential_vel)) deallocate(segment%nudged_tangential_vel) + if (associated (segment%tangential_grad)) deallocate(segment%tangential_grad) if (associated (segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg) @@ -2140,10 +2169,10 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) ni_seg = segment%ie_obc-segment%is_obc+1 nj_seg = segment%je_obc-segment%js_obc+1 - is_obc = max(segment%is_obc,isd-1) - ie_obc = min(segment%ie_obc,ied) - js_obc = max(segment%js_obc,jsd-1) - je_obc = min(segment%je_obc,jed) + is_obc = max(segment%is_obc,isd+1) + ie_obc = min(segment%ie_obc,ied-1) + js_obc = max(segment%js_obc,jsd+1) + je_obc = min(segment%je_obc,jed-1) ! Calculate auxiliary fields at staggered locations. ! Segment indices are on q points: @@ -2200,7 +2229,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) if (siz(3) /= segment%field(m)%nk_src) call MOM_error(FATAL,'nk_src inconsistency') if (segment%field(m)%nk_src > 1) then if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,G%ke)) @@ -2210,7 +2239,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) segment%field(m)%bt_vel(:,:)=0.0 endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,G%ke)) @@ -2222,7 +2251,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) endif else if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1)) else allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,1)) @@ -2232,7 +2261,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) segment%field(m)%bt_vel(:,:)=0.0 endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,1)) else allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,1)) @@ -2263,13 +2292,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) call time_interp_external(segment%field(m)%fid,Time, tmp_buffer) if (OBC%brushcutter_mode) then if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)+1:2*(je_obc+G%jdg_offset)+1:2,:) else segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)+1:2*(je_obc+G%jdg_offset):2,:) endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)+1:2*(ie_obc+G%idg_offset)+1:2,1,:) else segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)+1:2*(ie_obc+G%idg_offset):2,1,:) From 85ea88e95112fac91e1752d27e35b536b0e9add0 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 2 May 2018 10:40:39 -0800 Subject: [PATCH 08/10] Better fix to array-bounds trouble --- src/core/MOM_open_boundary.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 7f9e1ecaab..32b9f6d197 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -2169,10 +2169,10 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) ni_seg = segment%ie_obc-segment%is_obc+1 nj_seg = segment%je_obc-segment%js_obc+1 - is_obc = max(segment%is_obc,isd+1) - ie_obc = min(segment%ie_obc,ied-1) - js_obc = max(segment%js_obc,jsd+1) - je_obc = min(segment%je_obc,jed-1) + is_obc = max(segment%is_obc,isd-1) + ie_obc = min(segment%ie_obc,ied) + js_obc = max(segment%js_obc,jsd-1) + je_obc = min(segment%je_obc,jed) ! Calculate auxiliary fields at staggered locations. ! Segment indices are on q points: @@ -2357,7 +2357,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) I=is_obc if (segment%field(m)%name == 'V') then ! Do q points for the whole segment - do J=js_obc,je_obc + do J=max(js_obc,jsd),min(je_obc,jed-1) ! Using the h remapping approach ! Pretty sure we need to check for source/target grid consistency here segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer @@ -2400,7 +2400,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) J=js_obc if (segment%field(m)%name == 'U') then ! Do q points for the whole segment - do I=is_obc,ie_obc + do I=max(is_obc,isd),min(ie_obc,ied-1) segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer if (G%mask2dCv(i,J)>0. .and. G%mask2dCv(i+1,J)>0.) then ! Using the h remapping approach From 8eba1c94c6753d4ec52278a5b824692fb1ebc88b Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 2 May 2018 11:22:57 -0800 Subject: [PATCH 09/10] Fix MATT's trailing space. --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 99ca06bb66..cff8b6e372 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -347,7 +347,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) "PPM reconstruction will also be used within boundary cells.", & default=.false., do_not_log=.true.) - CS%new_sponges = .true. + CS%new_sponges = .true. CS%nz = G%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed From 7d1a50cc1eb6fa2d009203458b4bbefe432e286e Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 2 May 2018 19:13:29 -0800 Subject: [PATCH 10/10] Adding "specified" option to vorticity and strain. --- src/core/MOM_CoriolisAdv.F90 | 16 ++++- src/core/MOM_open_boundary.F90 | 62 ++++++++++++------- .../lateral/MOM_hor_visc.F90 | 16 ++++- 3 files changed, 69 insertions(+), 25 deletions(-) diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index 44449bbd2d..99dfdaa568 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -297,6 +297,13 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) dudy(I,J) = 2.0*(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%dxCu(I,j+1) endif enddo ; endif + if (OBC%specified_vorticity) then ; do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB + if (OBC%segment(n)%direction == OBC_DIRECTION_N) then + dudy(I,J) = OBC%segment(n)%tangential_grad(I,J,k)*G%dxCu(I,j)*G%dyBu(I,J) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_S) + dudy(I,J) = OBC%segment(n)%tangential_grad(I,J,k)*G%dxCu(I,j+1)*G%dyBu(I,J) + endif + enddo ; endif ! Project thicknesses across OBC points with a no-gradient condition. do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied) @@ -323,13 +330,20 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, CS) if (OBC%freeslip_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB dvdx(I,J) = 0. enddo ; endif - if (OBC%computed_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB + if (OBC%specified_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB if (OBC%segment(n)%direction == OBC_DIRECTION_E) then dvdx(I,J) = 2.0*(OBC%segment(n)%tangential_vel(I,J,k) - v(i,J,k))*G%dyCv(i,J) else ! (OBC%segment(n)%direction == OBC_DIRECTION_W) dvdx(I,J) = 2.0*(v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%dyCv(i+1,J) endif enddo ; endif + if (OBC%specified_vorticity) then ; do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB + if (OBC%segment(n)%direction == OBC_DIRECTION_E) then + dvdx(I,J) = OBC%segment(n)%tangential_grad(I,J,k)*G%dyCv(i,J)*G%dxBu(I,J) + else ! (OBC%segment(n)%direction == OBC_DIRECTION_W) + dvdx(I,J) = OBC%segment(n)%tangential_grad(I,J,k)*G%dyCv(i+1,J)*G%dxBu(I,J) + endif + enddo ; endif ! Project thicknesses across OBC points with a no-gradient condition. do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 32b9f6d197..fc78abf62f 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -196,14 +196,14 @@ module MOM_open_boundary !! in the relative vorticity on open boundaries. logical :: computed_vorticity = .false. !< If True, uses external data for tangential velocity !! in the relative vorticity on open boundaries. - logical :: imported_vorticity = .false. !< If True, uses external data for tangential velocity + logical :: specified_vorticity = .false. !< If True, uses external data for tangential velocity !! gradients in the relative vorticity on open boundaries. logical :: zero_strain = .false. !< If True, sets strain to zero on open boundaries. logical :: freeslip_strain = .false. !< If True, sets normal gradient of tangential velocity to zero !! in the strain on open boundaries. logical :: computed_strain = .false. !< If True, uses external data for tangential velocity to compute !! normal gradient in the strain on open boundaries. - logical :: imported_strain = .false. !< If True, uses external data for tangential velocity gradients + logical :: specified_strain = .false. !< If True, uses external data for tangential velocity gradients !! to compute strain on open boundaries. logical :: zero_biharmonic = .false. !< If True, zeros the Laplacian of flow on open boundaries for !! use in the biharmonic viscosity term. @@ -317,16 +317,16 @@ subroutine open_boundary_config(G, param_file, OBC) "If true, uses the external values of tangential velocity\n"// & "in the relative vorticity on open boundaries. This cannot\n"// & "be true if another OBC_XXX_VORTICITY option is True.", default=.false.) - call get_param(param_file, mdl, "OBC_IMPORTED_VORTICITY", OBC%imported_vorticity, & + call get_param(param_file, mdl, "OBC_SPECIFIED_VORTICITY", OBC%specified_vorticity, & "If true, uses the external values of tangential velocity\n"// & "in the relative vorticity on open boundaries. This cannot\n"// & "be true if another OBC_XXX_VORTICITY option is True.", default=.false.) if ((OBC%zero_vorticity .and. OBC%freeslip_vorticity) .or. & (OBC%zero_vorticity .and. OBC%computed_vorticity) .or. & - (OBC%zero_vorticity .and. OBC%imported_vorticity) .or. & + (OBC%zero_vorticity .and. OBC%specified_vorticity) .or. & (OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or. & - (OBC%freeslip_vorticity .and. OBC%imported_vorticity) .or. & - (OBC%computed_vorticity .and. OBC%imported_vorticity)) & + (OBC%freeslip_vorticity .and. OBC%specified_vorticity) .or. & + (OBC%computed_vorticity .and. OBC%specified_vorticity)) & call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& "Only one of OBC_ZERO_VORTICITY, OBC_FREESLIP_VORTICITY, OBC_COMPUTED_VORTICITY\n"//& "and OBC_IMPORTED_VORTICITY can be True at once.") @@ -341,16 +341,16 @@ subroutine open_boundary_config(G, param_file, OBC) "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & "be true if another OBC_XXX_STRAIN option is True.", default=.false.) - call get_param(param_file, mdl, "OBC_IMPORTED_STRAIN", OBC%imported_strain, & + call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN", OBC%specified_strain, & "If true, sets the normal gradient of tangential velocity to\n"// & "zero in the strain use in the stress tensor on open boundaries. This cannot\n"// & "be true if another OBC_XXX_STRAIN option is True.", default=.false.) if ((OBC%zero_strain .and. OBC%freeslip_strain) .or. & (OBC%zero_strain .and. OBC%computed_strain) .or. & - (OBC%zero_strain .and. OBC%imported_strain) .or. & + (OBC%zero_strain .and. OBC%specified_strain) .or. & (OBC%freeslip_strain .and. OBC%computed_strain) .or. & - (OBC%freeslip_strain .and. OBC%imported_strain) .or. & - (OBC%computed_strain .and. OBC%imported_strain)) & + (OBC%freeslip_strain .and. OBC%specified_strain) .or. & + (OBC%computed_strain .and. OBC%specified_strain)) & call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config:\n"//& "Only one of OBC_ZERO_STRAIN, OBC_FREESLIP_STRAIN, OBC_COMPUTED_STRAIN\n"//& "and OBC_IMPORTED_STRAIN can be True at once.") @@ -1968,7 +1968,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif - if (OBC%imported_vorticity .or. OBC%imported_strain) then + if (OBC%specified_vorticity .or. OBC%specified_strain) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -2001,7 +2001,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged_tan) then allocate(segment%nudged_tangential_vel(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_vel(:,:,:)=0.0 endif - if (OBC%imported_vorticity .or. OBC%imported_strain) then + if (OBC%specified_vorticity .or. OBC%specified_strain) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -2306,13 +2306,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) endif else if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+G%jdg_offset+1:je_obc+G%jdg_offset+1,:) else segment%field(m)%buffer_src(is_obc,:,:)=tmp_buffer(1,js_obc+G%jdg_offset+1:je_obc+G%jdg_offset,:) endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+G%idg_offset+1:ie_obc+G%idg_offset+1,1,:) else segment%field(m)%buffer_src(:,js_obc,:)=tmp_buffer(is_obc+G%idg_offset+1:ie_obc+G%idg_offset,1,:) @@ -2323,13 +2323,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer) if (OBC%brushcutter_mode) then if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)+1:2*(je_obc+G%jdg_offset)+1:2,:) else segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,2*(js_obc+G%jdg_offset)+1:2*(je_obc+G%jdg_offset):2,:) endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)+1:2*(ie_obc+G%idg_offset)+1:2,1,:) else segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(2*(is_obc+G%idg_offset)+1:2*(ie_obc+G%idg_offset):2,1,:) @@ -2337,13 +2337,13 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) endif else if (segment%is_E_or_W) then - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+G%jdg_offset+1:je_obc+G%jdg_offset+1,:) else segment%field(m)%dz_src(is_obc,:,:)=tmp_buffer(1,js_obc+G%jdg_offset+1:je_obc+G%jdg_offset,:) endif else - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+G%idg_offset+1:ie_obc+G%idg_offset+1,1,:) else segment%field(m)%dz_src(:,js_obc,:)=tmp_buffer(is_obc+G%idg_offset+1:ie_obc+G%idg_offset,1,:) @@ -2355,7 +2355,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) ishift=1 if (segment%direction == OBC_DIRECTION_E) ishift=0 I=is_obc - if (segment%field(m)%name == 'V') then + if (segment%field(m)%name == 'V' .or. segment%field(m)%name == 'DVDX') then ! Do q points for the whole segment do J=max(js_obc,jsd),min(je_obc,jed-1) ! Using the h remapping approach @@ -2398,7 +2398,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) jshift=1 if (segment%direction == OBC_DIRECTION_N) jshift=0 J=js_obc - if (segment%field(m)%name == 'U') then + if (segment%field(m)%name == 'U' .or. segment%field(m)%name == 'DUDY') then ! Do q points for the whole segment do I=max(is_obc,isd),min(ie_obc,ied-1) segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer @@ -2451,6 +2451,8 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) else if (segment%field(m)%name == 'U') then allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc:ie_obc,js_obc+1:je_obc)) + else if (segment%field(m)%name == 'DVDX') then + allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc+1:je_obc,G%ke)) endif @@ -2461,6 +2463,8 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) else if (segment%field(m)%name == 'V') then allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,G%ke)) allocate(segment%field(m)%bt_vel(is_obc+1:ie_obc,js_obc:je_obc)) + else if (segment%field(m)%name == 'DUDY') then + allocate(segment%field(m)%buffer_dst(is_obc:ie_obc,js_obc:je_obc,G%ke)) else allocate(segment%field(m)%buffer_dst(is_obc+1:ie_obc,js_obc:je_obc,G%ke)) endif @@ -2504,7 +2508,7 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) enddo elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_E_or_W .and. associated(segment%tangential_vel)) then I=is_obc - do J=js_obc+1,je_obc-1 + do J=js_obc,je_obc do k=1,G%ke segment%tangential_vel(I,J,k) = segment%field(m)%buffer_dst(I,J,k) enddo @@ -2512,12 +2516,26 @@ subroutine update_OBC_segment_data(G, GV, OBC, tv, h, Time) enddo elseif (trim(segment%field(m)%name) == 'U' .and. segment%is_N_or_S .and. associated(segment%tangential_vel)) then J=js_obc - do I=is_obc+1,ie_obc-1 + do I=is_obc,ie_obc do k=1,G%ke segment%tangential_vel(I,J,k) = segment%field(m)%buffer_dst(I,J,k) enddo if (associated(segment%nudged_tangential_vel)) segment%nudged_tangential_vel(I,J,:) = segment%tangential_vel(I,J,:) enddo + elseif (trim(segment%field(m)%name) == 'DVDX' .and. segment%is_E_or_W .and. associated(segment%tangential_grad)) then + I=is_obc + do J=js_obc,je_obc + do k=1,G%ke + segment%tangential_grad(I,J,k) = segment%field(m)%buffer_dst(I,J,k) + enddo + enddo + elseif (trim(segment%field(m)%name) == 'DUDY' .and. segment%is_N_or_S .and. associated(segment%tangential_grad)) then + J=js_obc + do I=is_obc,ie_obc + do k=1,G%ke + segment%tangential_grad(I,J,k) = segment%field(m)%buffer_dst(I,J,k) + enddo + enddo endif endif endif diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 85bb438f8f..5a7dbf7208 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -452,9 +452,15 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, dudy(I,J) = 0. elseif (OBC%computed_strain) then if (OBC%segment(n)%direction == OBC_DIRECTION_N) then - dudy(I,J) = CS%DX_dyBu(I,J)*(OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))*G%IdxCu(I,j) + dudy(I,J) = 2.0*CS%DX_dyBu(I,J)*(OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))*G%IdxCu(I,j) else - dudy(I,J) = CS%DX_dyBu(I,J)*(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdxCu(I,j+1) + dudy(I,J) = 2.0*CS%DX_dyBu(I,J)*(u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdxCu(I,j+1) + endif + elseif (OBC%specified_strain) then + if (OBC%segment(n)%direction == OBC_DIRECTION_N) then + dudy(I,J) = CS%DX_dyBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdxCu(I,j)*G%dxBu(I,J) + else + dudy(I,J) = CS%DX_dyBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdxCu(I,j+1)*G%dxBu(I,J) endif endif enddo @@ -470,6 +476,12 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, else dvdx(I,J) = 2.0*CS%DY_dxBu(I,J)*(v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdyCv(i+1,J) endif + elseif (OBC%specified_strain) then + if (OBC%segment(n)%direction == OBC_DIRECTION_E) then + dvdx(I,J) = CS%DY_dxBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdyCv(i,J)*G%dxBu(I,J) + else + dvdx(I,J) = CS%DY_dxBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdyCv(i+1,J)*G%dxBu(I,J) + endif endif enddo endif