diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 96fa98cbf3..b361cd7a82 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -1331,7 +1331,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & local_open_BC = .false. if (present(OBC)) then ; if (associated(OBC)) then - local_open_BC = OBC%open_u_BCs_exist_globally + local_open_BC = OBC%open_v_BCs_exist_globally endif ; endif do i=ish,ieh ; if (do_I(i)) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f35748dd4a..0573886c3e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -118,6 +118,7 @@ module MOM_open_boundary logical :: nudged_grad !< Optional supplement to nudge normal gradient of tangential velocity. logical :: specified !< Boundary normal velocity fixed to external value. logical :: specified_tan !< Boundary tangential velocity fixed to external value. + logical :: specified_grad !< Boundary gradient of 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 any external OBC fields are needed. @@ -436,6 +437,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) OBC%segment(l)%nudged_grad = .false. OBC%segment(l)%specified = .false. OBC%segment(l)%specified_tan = .false. + OBC%segment(l)%specified_grad = .false. OBC%segment(l)%open = .false. OBC%segment(l)%gradient = .false. OBC%segment(l)%values_needed = .false. @@ -941,6 +943,10 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) OBC%segment%u_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then OBC%segment(l_seg)%specified_tan = .true. + OBC%segment%v_values_needed = .true. + elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then + OBC%segment(l_seg)%specified_grad = .true. + OBC%segment%g_values_needed = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") @@ -1078,6 +1084,10 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) OBC%segment%v_values_needed = .true. elseif (trim(action_str(a_loop)) == 'SIMPLE_TAN') then OBC%segment(l_seg)%specified_tan = .true. + OBC%segment%u_values_needed = .true. + elseif (trim(action_str(a_loop)) == 'SIMPLE_GRAD') then + OBC%segment(l_seg)%specified_grad = .true. + OBC%segment%g_values_needed = .true. else call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc: "//& "String '"//trim(action_str(a_loop))//"' not understood.") @@ -2774,7 +2784,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt) enddo ; enddo endif if (segment%nudged_grad) then - do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB ! dhdt gets set to 0 on inflow in oblique case if (ry_tang_rad(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in @@ -3211,7 +3221,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0 endif if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. & - segment%oblique_grad) then + segment%oblique_grad .or. segment%specified_grad) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -3254,7 +3264,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0 endif if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. & - segment%oblique_grad) then + segment%oblique_grad .or. segment%specified_grad) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -3761,8 +3771,9 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) endif endif - if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then - if (segment%field(m)%fid>0) then ! calculate external BT velocity and transport if needed + if (segment%field(m)%fid>0) then + ! calculate external BT velocity and transport if needed + if (trim(segment%field(m)%name) == 'U' .or. trim(segment%field(m)%name) == 'V') then if (trim(segment%field(m)%name) == 'U' .and. segment%is_E_or_W) then I=is_obc do j=js_obc+1,je_obc @@ -3809,23 +3820,27 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) 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) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k) - enddo + endif + 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) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k) + if (associated(segment%nudged_tangential_grad)) & + segment%nudged_tangential_grad(I,J,:) = segment%tangential_grad(I,J,:) 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) = US%T_to_s*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) = US%T_to_s*segment%field(m)%buffer_dst(I,J,k) + if (associated(segment%nudged_tangential_grad)) & + segment%nudged_tangential_grad(I,J,:) = segment%tangential_grad(I,J,:) enddo - endif + enddo endif endif diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index b16b3a341c..3e07ac839e 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -49,6 +49,7 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth ) ! Local variables integer :: i, j real :: x, y, delta, dblen, dbfrac + logical :: dbrotate call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, & 'Lateral Length scale for dumbbell.',& @@ -56,20 +57,35 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth ) call get_param(param_file, mdl,"DUMBBELL_FRACTION",dbfrac, & 'Meridional fraction for narrow part of dumbbell.',& units='nondim', default=0.5, do_not_log=.false.) + call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & + 'Logical for rotation of dumbbell domain.',& + units='nondim', default=.false., do_not_log=.false.) if (G%x_axis_units == 'm') then dblen=dblen*1.e3 endif - do j=G%jsc,G%jec ; do i=G%isc,G%iec - ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) ) / dblen - y = ( G%geoLatT(i,j) ) / G%len_lat - D(i,j) = G%max_depth - if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then - D(i,j) = 0.0 - endif - enddo ; enddo + if (dbrotate) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! Compute normalized zonal coordinates (x,y=0 at center of domain) + x = ( G%geoLonT(i,j) ) / G%len_lon + y = ( G%geoLatT(i,j) ) / dblen + D(i,j) = G%max_depth + if ((y>=-0.25 .and. y<=0.25) .and. (x <= -0.5*dbfrac .or. x >= 0.5*dbfrac)) then + D(i,j) = 0.0 + endif + enddo ; enddo + else + do j=G%jsc,G%jec ; do i=G%isc,G%iec + ! Compute normalized zonal coordinates (x,y=0 at center of domain) + x = ( G%geoLonT(i,j) ) / dblen + y = ( G%geoLatT(i,j) ) / G%len_lat + D(i,j) = G%max_depth + if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then + D(i,j) = 0.0 + endif + enddo ; enddo + endif end subroutine dumbbell_initialize_topography @@ -209,6 +225,7 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file real :: x, y, dblen real :: T_ref, T_Light, T_Dense, S_ref, S_Light, S_Dense, a1, frac_dense, k_frac, res_rat logical :: just_read ! If true, just read parameters but set nothing. + logical :: dbrotate ! If true, rotate the domain. character(len=20) :: verticalCoordinate, density_profile is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -230,6 +247,9 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, & 'Lateral Length scale for dumbbell ',& units='k', default=600., do_not_log=just_read) + call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & + 'Logical for rotation of dumbbell domain.',& + units='nondim', default=.false., do_not_log=just_read) if (G%x_axis_units == 'm') then dblen=dblen*1.e3 @@ -238,7 +258,12 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file do j=G%jsc,G%jec do i=G%isc,G%iec ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) ) / dblen + if (dbrotate) then + ! This is really y in the rotated case + x = ( G%geoLatT(i,j) ) / dblen + else + x = ( G%geoLonT(i,j) ) / dblen + endif do k=1,nz T(i,j,k)=T_surf enddo @@ -278,9 +303,13 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, integer :: i, j, k, nz real :: x, zi, zmid, dist, min_thickness, dblen real :: mld, S_ref, S_range, S_dense, T_ref, sill_height + logical :: dbrotate ! If true, rotate the domain. call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, & 'Lateral Length scale for dumbbell ',& units='k', default=600., do_not_log=.true.) + call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & + 'Logical for rotation of dumbbell domain.',& + units='nondim', default=.false., do_not_log=.true.) if (G%x_axis_units == 'm') then dblen=dblen*1.e3 @@ -307,7 +336,12 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, do i = G%isc,G%iec if (G%mask2dT(i,j) > 0.) then ! nondimensional x position - x = (G%geoLonT(i,j) ) / dblen + if (dbrotate) then + ! This is really y in the rotated case + x = ( G%geoLatT(i,j) ) / dblen + else + x = ( G%geoLonT(i,j) ) / dblen + endif if (x > 0.25 .or. x < -0.25) then ! scale restoring by depth into sponge Idamp(i,j) = 1. / sponge_time_scale @@ -339,18 +373,23 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, do j=G%jsc,G%jec ; do i=G%isc,G%iec ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) ) / dblen - if (x>=0.25 ) then - do k=1,nz - S(i,j,k)=S_ref + 0.5*S_range - enddo - endif - if (x<=-0.25 ) then - do k=1,nz - S(i,j,k)=S_ref - 0.5*S_range - enddo - endif - enddo ; enddo + if (dbrotate) then + ! This is really y in the rotated case + x = ( G%geoLatT(i,j) ) / dblen + else + x = ( G%geoLonT(i,j) ) / dblen + endif + if (x>=0.25 ) then + do k=1,nz + S(i,j,k)=S_ref + 0.5*S_range + enddo + endif + if (x<=-0.25 ) then + do k=1,nz + S(i,j,k)=S_ref - 0.5*S_range + enddo + endif + enddo ; enddo endif if (associated(tv%S)) call set_up_ALE_sponge_field(S, G, tv%S, ACSp) diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index d6d6dea11a..2d19cce6dd 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -189,6 +189,7 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) real :: S_surf, S_range real :: x, y integer :: i, j + logical :: dbrotate ! If true, rotate the domain. #include "version_variable.h" character(len=40) :: mdl = "dumbbell_surface_forcing" ! This module's name. @@ -224,6 +225,9 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) call get_param(param_file, mdl, "DUMBBELL_SLP_PERIOD", CS%slp_period, & "Periodicity of SLP forcing in reservoirs.", & units="days", default = 1.0) + call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, & + 'Logical for rotation of dumbbell domain.',& + units='nondim', default=.false., do_not_log=.true.) call get_param(param_file, mdl,"INITIAL_SSS", S_surf, & "Initial surface salinity", units="1e-3", default=34.0, do_not_log=.true.) call get_param(param_file, mdl,"INITIAL_S_RANGE", S_range, & @@ -250,8 +254,14 @@ subroutine dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS) do j=G%jsc,G%jec do i=G%isc,G%iec ! Compute normalized zonal coordinates (x,y=0 at center of domain) - x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 - y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 +! x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 +! y = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 + if (dbrotate) then + ! This is really y in the rotated case + x = ( G%geoLatT(i,j) - G%south_lat ) / G%len_lat - 0.5 + else + x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon - 0.5 + endif CS%forcing_mask(i,j)=0 CS%S_restore(i,j) = S_surf if ((x>0.25)) then