diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index c32352acb1..99dfdaa568 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -290,6 +290,20 @@ 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 + 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) @@ -316,6 +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%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 35b442c816..fc78abf62f 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,10 @@ 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(:,:,:) :: 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 @@ -143,6 +149,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 @@ -186,9 +194,17 @@ 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 :: 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 :: 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. logical :: brushcutter_mode = .false. !< If True, read data on supergrid. @@ -296,20 +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.) - 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.") + "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 another OBC_XXX_VORTICITY option is True.", default=.false.) + 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%specified_vorticity) .or. & + (OBC%freeslip_vorticity .and. OBC%computed_vorticity) .or. & + (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.") 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 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 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 another OBC_XXX_STRAIN option is True.", default=.false.) + 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%specified_strain) .or. & + (OBC%freeslip_strain .and. OBC%computed_strain) .or. & + (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.") 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.) @@ -342,7 +386,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. @@ -768,6 +814,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.") @@ -869,6 +917,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.") @@ -1411,7 +1461,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 +1511,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 = 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 = segment%grad_normal(J,1,k) + 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 +1561,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 +1599,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 +1651,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,17 +1953,28 @@ 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 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 .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%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 - 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%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 @@ -1925,17 +1986,28 @@ 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 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 .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%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 - 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 + allocate(segment%ry_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%ry_normal(:,:,:)=0.0 endif endif @@ -1956,11 +2028,14 @@ 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) 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) @@ -2154,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)) @@ -2164,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)) @@ -2176,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)) @@ -2186,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)) @@ -2217,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,:) @@ -2231,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,:) @@ -2248,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,:) @@ -2262,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,:) @@ -2280,9 +2355,9 @@ 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 - ! Only do q points within the segment - do J=js_obc+1,je_obc-1 + 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 ! 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 @@ -2323,9 +2398,9 @@ 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 - ! Only do q points within the segment - do I=is_obc+1,ie_obc-1 + 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 if (G%mask2dCv(i,J)>0. .and. G%mask2dCv(i+1,J)>0.) then ! Using the h remapping approach @@ -2335,13 +2410,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,:), & @@ -2376,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 @@ -2386,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 @@ -2427,6 +2506,36 @@ 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,je_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) == 'U' .and. segment%is_N_or_S .and. associated(segment%tangential_vel)) then + J=js_obc + 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 @@ -3019,7 +3128,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') diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 0e02cefba2..5a7dbf7208 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -443,13 +443,25 @@ 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) = 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) = 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 elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-2) .and. (I <= Ieq+1)) then @@ -458,6 +470,18 @@ 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 + 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 diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index e42601d51b..cff8b6e372 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 - ! GMM - Planned extension: Support for time varying sponge targets. implicit none ; private @@ -166,6 +165,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 +178,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 +191,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 +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%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 +362,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 +374,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 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