From a7f99afa512240daee3751dc2e621e9b5ed09a86 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 30 Oct 2024 12:16:21 -0800 Subject: [PATCH 01/19] Fix for MASK_OUTSIDE_OBCS with MASKING_DEPTH --- src/core/MOM_open_boundary.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 68e1a03669..ff90658460 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -5052,7 +5052,9 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: mask_depth ! The masking depth for ocean points [Z ~> m] + real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside, @@ -5062,6 +5064,12 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & + units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.) + + Dmask = mask_depth + if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth + ! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref allocate(color(G%isd:G%ied, G%jsd:G%jed), source=0.0) @@ -5152,7 +5160,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.) endif - if (color(i,j) == cout) G%bathyT(i,j) = min_depth + if (color(i,j) == cout) G%bathyT(i,j) = Dmask enddo ; enddo if (fatal_error) call MOM_error(FATAL, & "MOM_open_boundary: inconsistent OBC segments.") From ef18f5a214c4da065469c48cda73db67648717d7 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 6 Nov 2024 15:25:49 -0900 Subject: [PATCH 02/19] Fix the tracer advection outside of interior OBCs. - Otherwise, the tracer values just outside the OBC get updated based on fluxes at the OBC and quickly go out of bounds of the equation of state. --- src/tracer/MOM_tracer_advect.F90 | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 0e129b2d03..175006b556 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -648,6 +648,19 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do_i(i+1,j) = .false. + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do_i(i,j) = .false. + endif + endif ; enddo + endif + endif ; endif + ! update tracer concentration from i-flux and save some diagnostics do m=1,ntr @@ -1039,6 +1052,19 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else ; do_i(i,j) = .false. ; endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=is,ie-1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + do_i(i,j+1) = .false. + elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif + endif ; enddo + endif + endif ; endif + ! update tracer and save some diagnostics do m=1,ntr do i=is,ie ; if (do_i(i,j)) then From 01b0dc4edc246c941adc1712530fb39579fd757c Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 7 Nov 2024 16:03:09 -0900 Subject: [PATCH 03/19] A more correct fix to tracer advection - The previous version did the wrong thing at northern boundaries, at a southern corner too. --- src/tracer/MOM_tracer_advect.F90 | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 175006b556..e11ee6e5a5 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -1055,13 +1055,18 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) if (associated(OBC)) then ; if (OBC%OBC_pe) then if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then - do i=is,ie-1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - do_i(i,j+1) = .false. - elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do_i(i,j) = .false. + do i=is,ie + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + do_i(i,j) = .false. + endif endif - endif ; enddo + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif + endif + enddo endif endif ; endif From 4add0e4a9ee8bf8094e75ea628cade72a3d449e8 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 8 Nov 2024 10:09:44 -0900 Subject: [PATCH 04/19] Better fix to eta outside of OBCs issue - It hasn't yet caused a blowup that I know of, but better to prevent any trouble while we're thinking about it. --- src/core/MOM_barotropic.F90 | 40 +++++++++++++++++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index af2beca1fb..0f57335322 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -742,6 +742,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: ioff, joff integer :: l_seg + real :: factor(SZI_(G),SZJ_(G)) ! If non-zero, work on given points. if (.not.CS%module_is_initialized) call MOM_error(FATAL, & "btstep: Module MOM_barotropic must be initialized before it is used.") @@ -2457,17 +2458,52 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m) endif + do j=jsv,jev + do i=isv,iev + factor(i,j) = CS%IareaT(i,j) + enddo + enddo + + ! Update factor so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + do j=jsv,jev + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + factor(i+1,j) = 0.0 + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + factor(i,j) = 0.0 + endif + endif ; enddo + endif + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=isv,iev + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + factor(i,j) = 0.0 + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + factor(i,j) = 0.0 + endif + endif + enddo + endif + enddo + endif ; endif + if (integral_BT_cont) then !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * & + eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + factor(i,j) * & ((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo else !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * factor(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo From 52605fb64378b7aff95e746b715717de70f31db2 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 20 Nov 2024 17:48:32 -0900 Subject: [PATCH 05/19] More tinkering with OBC_in vs CS%OBC --- src/core/MOM.F90 | 20 ++++++++++++ src/core/MOM_open_boundary.F90 | 56 +++++++++++++++++----------------- 2 files changed, 48 insertions(+), 28 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7ee90d746a..059e515051 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2883,6 +2883,26 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! reservoirs are used. call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, & param_file, restart_CSp, use_temperature) + if (turns /= 0) then + if (CS%OBC%radiation_BCs_exist_globally) then + OBC_in%rx_normal => CS%OBC%rx_normal + OBC_in%ry_normal => CS%OBC%ry_normal + endif + if (CS%OBC%oblique_BCs_exist_globally) then + OBC_in%rx_oblique_u => CS%OBC%rx_oblique_u + OBC_in%ry_oblique_u => CS%OBC%ry_oblique_u + OBC_in%rx_oblique_v => CS%OBC%rx_oblique_v + OBC_in%ry_oblique_v => CS%OBC%ry_oblique_v + OBC_in%cff_normal_u => CS%OBC%cff_normal_u + OBC_in%cff_normal_v => CS%OBC%cff_normal_v + endif + if (any(CS%OBC%tracer_x_reservoirs_used)) then + OBC_in%tres_x => CS%OBC%tres_x + endif + if (any(CS%OBC%tracer_y_reservoirs_used)) then + OBC_in%tres_y => CS%OBC%tres_y + endif + endif endif if (present(waves_CSp)) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 69e8cf1796..9543d5b515 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -352,24 +352,24 @@ module MOM_open_boundary type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for !! thickness-based fields on segments type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries - real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of + real, pointer :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of + real, pointer :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation + real, pointer :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation !! rates at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation + real, pointer :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation !! rates at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] - real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] logical :: debug !< If true, write verbose checksums for debugging purposes. real :: silly_h !< A silly value of thickness outside of the domain that can be used to test !! the independence of the OBCs to this external data [Z ~> m]. @@ -1948,15 +1948,15 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) call create_group_pass(OBC%pass_oblique, OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) call do_group_pass(OBC%pass_oblique, G%Domain) endif - if (allocated(OBC%tres_x) .and. allocated(OBC%tres_y)) then + if (associated(OBC%tres_x) .and. associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair) enddo - elseif (allocated(OBC%tres_x)) then + elseif (associated(OBC%tres_x)) then do m=1,OBC%ntr call pass_var(OBC%tres_x(:,:,:,m), G%Domain, position=EAST_FACE) enddo - elseif (allocated(OBC%tres_y)) then + elseif (associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_var(OBC%tres_y(:,:,:,m), G%Domain, position=NORTH_FACE) enddo @@ -2001,16 +2001,16 @@ subroutine open_boundary_dealloc(OBC) if (allocated(OBC%segment)) deallocate(OBC%segment) if (allocated(OBC%segnum_u)) deallocate(OBC%segnum_u) if (allocated(OBC%segnum_v)) deallocate(OBC%segnum_v) - if (allocated(OBC%rx_normal)) deallocate(OBC%rx_normal) - if (allocated(OBC%ry_normal)) deallocate(OBC%ry_normal) - if (allocated(OBC%rx_oblique_u)) deallocate(OBC%rx_oblique_u) - if (allocated(OBC%ry_oblique_u)) deallocate(OBC%ry_oblique_u) - if (allocated(OBC%rx_oblique_v)) deallocate(OBC%rx_oblique_v) - if (allocated(OBC%ry_oblique_v)) deallocate(OBC%ry_oblique_v) - if (allocated(OBC%cff_normal_u)) deallocate(OBC%cff_normal_u) - if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v) - if (allocated(OBC%tres_x)) deallocate(OBC%tres_x) - if (allocated(OBC%tres_y)) deallocate(OBC%tres_y) + if (associated(OBC%rx_normal)) nullify(OBC%rx_normal) + if (associated(OBC%ry_normal)) nullify(OBC%ry_normal) + if (associated(OBC%rx_oblique_u)) nullify(OBC%rx_oblique_u) + if (associated(OBC%ry_oblique_u)) nullify(OBC%ry_oblique_u) + if (associated(OBC%rx_oblique_v)) nullify(OBC%rx_oblique_v) + if (associated(OBC%ry_oblique_v)) nullify(OBC%ry_oblique_v) + if (associated(OBC%cff_normal_u)) nullify(OBC%cff_normal_u) + if (associated(OBC%cff_normal_v)) nullify(OBC%cff_normal_v) + if (associated(OBC%tres_x)) nullify(OBC%tres_x) + if (associated(OBC%tres_y)) nullify(OBC%tres_y) if (associated(OBC%remap_z_CS)) deallocate(OBC%remap_z_CS) if (associated(OBC%remap_h_CS)) deallocate(OBC%remap_h_CS) deallocate(OBC) @@ -3369,7 +3369,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2) endif if (OBC%ntr == 0) return - if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return + if (.not. associated (OBC%tres_x) .or. .not. associated (OBC%tres_y)) return do m=1,OBC%ntr write(var_num,'(I3.3)') m call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, & @@ -5489,7 +5489,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(I,j,k)+ & ((u_L_out+a_out)*Reg%Tr(ntr_id)%t(I+ishift,j,k) - & (u_L_in+a_in)*segment%tr_Reg%Tr(m)%t(I,j,k))) - if (allocated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) + if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif enddo enddo @@ -5529,7 +5529,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(i,J,k) + & ((v_L_out+a_out)*Reg%Tr(ntr_id)%t(i,J+jshift,k) - & (v_L_in+a_in)*segment%tr_Reg%Tr(m)%t(i,J,k))) - if (allocated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + if (associated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif enddo enddo @@ -5605,7 +5605,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:) - if (allocated(OBC%tres_x)) then ; do k=1,nz + if (associated(OBC%tres_x)) then ; do k=1,nz OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif @@ -5672,7 +5672,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:) - if (allocated(OBC%tres_y)) then ; do k=1,nz + if (associated(OBC%tres_y)) then ; do k=1,nz OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif From 677baba59afaed72b8b84cb470f6bd5f7450b826 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Fri, 22 Nov 2024 10:48:08 -0900 Subject: [PATCH 06/19] Another OBC fix, this for the barotropic mode. --- src/core/MOM_open_boundary.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9543d5b515..bf2f7cb2e9 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6117,6 +6117,14 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%buffer_src) endif + if (allocated(segment_in%field(n)%buffer_dst)) then + call allocate_rotated_array(segment_in%field(n)%buffer_dst, & + lbound(segment_in%field(n)%buffer_dst), turns, & + segment%field(n)%buffer_dst) + call rotate_array(segment_in%field(n)%buffer_dst, turns, & + segment%field(n)%buffer_dst) + endif + segment%field(n)%nk_src = segment_in%field(n)%nk_src if (allocated(segment_in%field(n)%dz_src)) then @@ -6130,6 +6138,8 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%value = segment_in%field(n)%value enddo + call rotate_array(segment_in%SSH, turns, segment%SSH) + segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists end subroutine rotate_OBC_segment_data From d60d2eabfe12ef92ea9f51f224f4ede1d598b256 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Sat, 23 Nov 2024 08:06:25 -0900 Subject: [PATCH 07/19] More OBC rotation foxes. Still not out of the woods. --- src/core/MOM_open_boundary.F90 | 43 +++++++++++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index bf2f7cb2e9..8a5b42d16e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6138,7 +6138,48 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%value = segment_in%field(n)%value enddo - call rotate_array(segment_in%SSH, turns, segment%SSH) + if (allocated(segment_in%SSH)) & + call rotate_array(segment_in%SSH, turns, segment%SSH) + if (allocated(segment_in%cg)) & + call rotate_array(segment_in%cg, turns, segment%cg) + if (allocated(segment_in%htot)) & + call rotate_array(segment_in%htot, turns, segment%htot) + if (allocated(segment_in%dztot)) & + call rotate_array(segment_in%dztot, turns, segment%dztot) + if (allocated(segment_in%h)) & + call rotate_array(segment_in%h, turns, segment%h) + if (allocated(segment_in%normal_vel)) & + call rotate_array(segment_in%normal_vel, turns, segment%normal_vel) + if (allocated(segment_in%normal_trans)) & + call rotate_array(segment_in%normal_trans, turns, segment%normal_trans) + if (allocated(segment_in%normal_vel_bt)) & + call rotate_array(segment_in%normal_vel_bt, turns, segment%normal_vel_bt) + if (allocated(segment_in%tangential_vel)) & + call rotate_array(segment_in%tangential_vel, turns, segment%tangential_vel) + if (allocated(segment_in%tangential_grad)) & + call rotate_array(segment_in%tangential_grad, turns, segment%tangential_grad) + if (allocated(segment_in%grad_normal)) & + call rotate_array(segment_in%grad_normal, turns, segment%grad_normal) + if (allocated(segment_in%grad_tan)) & + call rotate_array(segment_in%grad_tan, turns, segment%grad_tan) + if (allocated(segment_in%grad_gradient)) & + call rotate_array(segment_in%grad_gradient, turns, segment%grad_gradient) + if (allocated(segment_in%rx_norm_rad)) & + call rotate_array(segment_in%rx_norm_rad, turns, segment%ry_norm_rad) + if (allocated(segment_in%ry_norm_rad)) & + call rotate_array(segment_in%ry_norm_rad, turns, segment%rx_norm_rad) + if (allocated(segment_in%rx_norm_obl)) & + call rotate_array(segment_in%rx_norm_obl, turns, segment%ry_norm_obl) + if (allocated(segment_in%ry_norm_obl)) & + call rotate_array(segment_in%ry_norm_obl, turns, segment%rx_norm_obl) + if (allocated(segment_in%cff_normal)) & + call rotate_array(segment_in%cff_normal, turns, segment%cff_normal) + if (allocated(segment_in%nudged_normal_vel)) & + call rotate_array(segment_in%nudged_normal_vel, turns, segment%nudged_normal_vel) + if (allocated(segment_in%nudged_tangential_vel)) & + call rotate_array(segment_in%nudged_tangential_vel, turns, segment%nudged_tangential_vel) + if (allocated(segment_in%nudged_tangential_grad)) & + call rotate_array(segment_in%nudged_tangential_grad, turns, segment%nudged_tangential_grad) segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists From 58e7fb6f2f37ad7ecc1fd9a7199fb4a931466300 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 27 Nov 2024 17:26:44 -0900 Subject: [PATCH 08/19] Fixing case on tr_Reg. Still not rototing OBCs correctly. --- src/core/MOM.F90 | 9 +++++++-- src/core/MOM_open_boundary.F90 | 9 +++++---- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 059e515051..0efc18b902 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -113,7 +113,7 @@ module MOM use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields -use MOM_open_boundary, only : open_boundary_setup_vert +use MOM_open_boundary, only : open_boundary_setup_vert, update_OBC_segment_data use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS @@ -3025,8 +3025,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, G, GV, CS%S) endif - if (associated(OBC_in)) & + if (associated(OBC_in)) then call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) + if (CS%OBC%some_need_no_IO_for_data) then + call calc_derived_thermo(CS%tv, CS%h, G, GV, US) + call update_OBC_segment_data(G, GV, US, CS%OBC, CS%tv, CS%h, Time) + endif + endif deallocate(u_in) deallocate(v_in) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 8a5b42d16e..3797f4bec1 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -5463,8 +5463,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(I+ishift,j) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 @@ -5507,8 +5507,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(i,j+jshift) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 @@ -6094,6 +6094,7 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%name = 'Vphase' case ('V') segment%field(n)%name = 'U' + segment%field(n)%value = -segment%field(n)%value case ('Vamp') segment%field(n)%name = 'Uamp' case ('Vphase') From 8ffb6f8084e33f33c9bd00f8686d8bd3897ca503 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Thu, 28 Nov 2024 11:16:36 -0900 Subject: [PATCH 09/19] Looks good for the first few steps, then not so much. --- src/core/MOM_open_boundary.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 3797f4bec1..fb248ad6a7 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6084,6 +6084,14 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%handle = segment_in%field(n)%handle segment%field(n)%dz_handle = segment_in%field(n)%dz_handle + if (allocated(segment_in%field(n)%buffer_dst)) then + call allocate_rotated_array(segment_in%field(n)%buffer_dst, & + lbound(segment_in%field(n)%buffer_dst), turns, & + segment%field(n)%buffer_dst) + call rotate_array(segment_in%field(n)%buffer_dst, turns, & + segment%field(n)%buffer_dst) + endif + if (modulo(turns, 2) /= 0) then select case (segment_in%field(n)%name) case ('U') @@ -6094,7 +6102,7 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%name = 'Vphase' case ('V') segment%field(n)%name = 'U' - segment%field(n)%value = -segment%field(n)%value + segment%field(n)%buffer_dst(:,:,:) = -segment%field(n)%buffer_dst(:,:,:) case ('Vamp') segment%field(n)%name = 'Uamp' case ('Vphase') @@ -6118,14 +6126,6 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%buffer_src) endif - if (allocated(segment_in%field(n)%buffer_dst)) then - call allocate_rotated_array(segment_in%field(n)%buffer_dst, & - lbound(segment_in%field(n)%buffer_dst), turns, & - segment%field(n)%buffer_dst) - call rotate_array(segment_in%field(n)%buffer_dst, turns, & - segment%field(n)%buffer_dst) - endif - segment%field(n)%nk_src = segment_in%field(n)%nk_src if (allocated(segment_in%field(n)%dz_src)) then From 86e026097bf82a6c63dddc679d99a137cfd8e484 Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 4 Dec 2024 12:37:57 -0900 Subject: [PATCH 10/19] Making barotropic fix look symmetric. --- src/core/MOM.F90 | 8 ++------ src/core/MOM_barotropic.F90 | 31 ++++++++++++++++++++++++------- src/core/MOM_open_boundary.F90 | 6 +++--- 3 files changed, 29 insertions(+), 16 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 0efc18b902..db33b523f1 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -114,6 +114,7 @@ module MOM use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields use MOM_open_boundary, only : open_boundary_setup_vert, update_OBC_segment_data +use MOM_open_boundary, only : initialize_segment_data use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS @@ -3025,13 +3026,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, G, GV, CS%S) endif - if (associated(OBC_in)) then + if (associated(OBC_in)) & call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) - if (CS%OBC%some_need_no_IO_for_data) then - call calc_derived_thermo(CS%tv, CS%h, G, GV, US) - call update_OBC_segment_data(G, GV, US, CS%OBC, CS%tv, CS%h, Time) - endif - endif deallocate(u_in) deallocate(v_in) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 0f57335322..bee598b8d6 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2466,16 +2466,33 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! Update factor so that nothing changes outside of the OBC (problem for interior OBCs only) if (associated(OBC)) then ; if (OBC%OBC_pe) then - do j=jsv,jev +! do j=jsv,jev +! if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then +! do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then +! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then +! factor(i+1,j) = 0.0 +! elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then +! factor(i,j) = 0.0 +! endif +! endif ; enddo +! endif + do i=isv,iev if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then - do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - factor(i+1,j) = 0.0 - elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then - factor(i,j) = 0.0 + do j=jsv,jev + if (OBC%segnum_u(I-1,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I-1,j))%direction == OBC_DIRECTION_E) then + factor(i,j) = 0.0 + endif + endif + if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + factor(i,j) = 0.0 + endif endif - endif ; enddo + enddo endif + enddo + do j=jsv,jev if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then do i=isv,iev if (OBC%segnum_v(i,J-1) /= OBC_NONE) then diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index fb248ad6a7..2c553814ae 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6055,13 +6055,13 @@ subroutine rotate_OBC_init(OBC_in, G, GV, US, param_file, tv, restart_CS, OBC) "If true, Temperature and salinity are used as state "//& "variables.", default=.true., do_not_log=.true.) + if (use_temperature) & + call fill_temp_salt_segments(G, GV, US, OBC, tv) + do l = 1, OBC%number_of_segments call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), G%HI%turns) enddo - if (use_temperature) & - call fill_temp_salt_segments(G, GV, US, OBC, tv) - call open_boundary_init(G, GV, US, param_file, OBC, restart_CS) end subroutine rotate_OBC_init From 37fbbbf15bb4590b8132615c8f9538635eeb8b9f Mon Sep 17 00:00:00 2001 From: Kate Hedstrom Date: Wed, 4 Dec 2024 12:37:57 -0900 Subject: [PATCH 11/19] Making barotropic fix look symmetric. --- src/core/MOM.F90 | 3 +-- src/core/MOM_open_boundary.F90 | 44 ++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index db33b523f1..059e515051 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -113,8 +113,7 @@ module MOM use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type use MOM_open_boundary, only : register_temp_salt_segments, update_segment_tracer_reservoirs use MOM_open_boundary, only : open_boundary_register_restarts, remap_OBC_fields -use MOM_open_boundary, only : open_boundary_setup_vert, update_OBC_segment_data -use MOM_open_boundary, only : initialize_segment_data +use MOM_open_boundary, only : open_boundary_setup_vert use MOM_open_boundary, only : rotate_OBC_config, rotate_OBC_init use MOM_porous_barriers, only : porous_widths_layer, porous_widths_interface, porous_barriers_init use MOM_porous_barriers, only : porous_barrier_CS diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 2c553814ae..22322fea82 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6181,6 +6181,50 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) call rotate_array(segment_in%nudged_tangential_vel, turns, segment%nudged_tangential_vel) if (allocated(segment_in%nudged_tangential_grad)) & call rotate_array(segment_in%nudged_tangential_grad, turns, segment%nudged_tangential_grad) + if (associated(segment_in%tr_Reg)) then + do n = 1, segment_in%tr_Reg%ntseg + call rotate_array(segment_in%tr_Reg%tr(n)%tres, turns, segment%tr_Reg%tr(n)%tres) + ! Testing this to see if it works for contant tres values. Probably wrong otherwise. + segment%tr_Reg%Tr(n)%is_initialized=.true. + enddo + endif + + ! Only doing this here for 1 turns, get flow directed in the opposite direction. + ! Currents that are now E_W were N_S and the turns should change their sign. + if (turns == 1 .and. segment%is_E_or_W) then + if (allocated(segment_in%normal_vel)) & + segment%normal_vel(:,:,:) = - segment%normal_vel(:,:,:) + if (allocated(segment_in%normal_trans)) & + segment%normal_trans(:,:,:) = - segment%normal_trans(:,:,:) + if (allocated(segment_in%normal_vel_bt)) & + segment%normal_vel_bt(:,:) = - segment%normal_vel_bt(:,:) + if (allocated(segment_in%tangential_vel)) & + segment%tangential_vel(:,:,:) = - segment%tangential_vel(:,:,:) + if (allocated(segment_in%tangential_grad)) & + segment%tangential_grad(:,:,:) = - segment%tangential_grad(:,:,:) + if (allocated(segment_in%grad_normal)) & + segment%grad_normal(:,:,:) = - segment%grad_normal(:,:,:) + if (allocated(segment_in%grad_tan)) & + segment%grad_tan(:,:,:) = - segment%grad_tan(:,:,:) + if (allocated(segment_in%grad_gradient)) & + segment%grad_gradient(:,:,:) = - segment%grad_gradient(:,:,:) + if (allocated(segment%rx_norm_rad)) & + segment%rx_norm_rad(:,:,:) = - segment%rx_norm_rad(:,:,:) + if (allocated(segment%ry_norm_rad)) & + segment%ry_norm_rad(:,:,:) = - segment%ry_norm_rad(:,:,:) + if (allocated(segment%rx_norm_obl)) & + segment%rx_norm_obl(:,:,:) = - segment%rx_norm_obl(:,:,:) + if (allocated(segment%ry_norm_obl)) & + segment%ry_norm_obl(:,:,:) = - segment%ry_norm_obl(:,:,:) + if (allocated(segment_in%cff_normal)) & + segment%cff_normal(:,:,:) = - segment%cff_normal(:,:,:) + if (allocated(segment_in%nudged_normal_vel)) & + segment%nudged_normal_vel(:,:,:) = - segment%nudged_normal_vel(:,:,:) + if (allocated(segment_in%nudged_tangential_vel)) & + segment%nudged_tangential_vel(:,:,:) = - segment%nudged_tangential_vel(:,:,:) + if (allocated(segment_in%nudged_tangential_grad)) & + segment%nudged_tangential_grad(:,:,:) = - segment%nudged_tangential_grad(:,:,:) + endif segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists From 381a69fc456cc03d1f8251cdb11d18a87810f7a8 Mon Sep 17 00:00:00 2001 From: He Wang <35150900+herrwang0@users.noreply.github.com> Date: Sat, 30 Nov 2024 15:29:57 -0500 Subject: [PATCH 12/19] Bug fix for write_energy with short dt (#749) Fix a bug with subroutine write_energy when using a DT<2. Otherwise, the energy outputs are written at wrong time steps. The reason was that time type divide is essentially a floor. So DT/2 = 0 if DT<2. --- src/diagnostics/MOM_sum_output.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 398241b98c..f5ff19630b 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -21,9 +21,9 @@ module MOM_sum_output use MOM_io, only : attribute_info, set_attribute_info, delete_attribute_info use MOM_io, only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_spatial_means, only : array_global_min_max -use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>) +use MOM_time_manager, only : time_type, get_time, get_date, set_time use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<) +use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<), operator(>) use MOM_time_manager, only : get_calendar_type, time_type_to_real, NO_CALENDAR use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_stocks use MOM_unit_scaling, only : unit_scale_type @@ -489,7 +489,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci CS%write_energy_time = CS%Start_time + CS%energysavedays * & (1 + (day - CS%Start_time) / CS%energysavedays) endif - elseif (day + (dt_force/2) <= CS%write_energy_time) then + elseif (day + (dt_force/2) < CS%write_energy_time) then return ! Do not write this step else ! Determine the next write time before proceeding if (CS%energysave_geometric) then From 4f5044283b078c846ef7aead48ad718e16d69d3e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 30 Nov 2024 05:11:58 -0500 Subject: [PATCH 13/19] Remove extra copy of compute_global_grid_integrals The subroutine compute_global_grid_integrals appeared in both the MOM_state_initialization and MOM_shared_initialization modules, but was only being called from the latter. This commit removes the extra copy in MOM_state_initialization. It also removes some unnecessary parentheses in the copy that is being retained, in part to facilitate the review of this commit. All answers are bitwise identical, and no publicly visible interfaces are altered. --- .../MOM_shared_initialization.F90 | 5 ++--- .../MOM_state_initialization.F90 | 20 ------------------- 2 files changed, 2 insertions(+), 23 deletions(-) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index edf08da3aa..eabd376512 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1327,10 +1327,9 @@ subroutine compute_global_grid_integrals(G, US) G%areaT_global = reproducing_sum(tmpForSumming) if (G%areaT_global == 0.0) & - call MOM_error(FATAL, "compute_global_grid_integrals: "//& - "zero ocean area (check topography?)") + call MOM_error(FATAL, "compute_global_grid_integrals: zero ocean area (check topography?)") - G%IareaT_global = 1.0 / (G%areaT_global) + G%IareaT_global = 1.0 / G%areaT_global end subroutine compute_global_grid_integrals ! ----------------------------------------------------------------------------- diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 769d60d51d..629da43ddc 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2339,26 +2339,6 @@ subroutine set_velocity_depth_max(G) enddo ; enddo end subroutine set_velocity_depth_max -!> Subroutine to pre-compute global integrals of grid quantities for -!! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G, US) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - ! Local variables - real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming ! Masked and unscaled areas for sums [m2] - real :: area_scale ! A conversion factor to prepare for reproducing sums [m2 L-2 ~> 1] - integer :: i,j - - area_scale = US%L_to_m**2 - tmpForSumming(:,:) = 0. - G%areaT_global = 0.0 ; G%IareaT_global = 0.0 - do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) - enddo ; enddo - G%areaT_global = reproducing_sum(tmpForSumming) - G%IareaT_global = 1. / (G%areaT_global) -end subroutine compute_global_grid_integrals - !> This subroutine sets the 4 bottom depths at velocity points to be the !! minimum of the adjacent depths. subroutine set_velocity_depth_min(G) From f32ff090a635f4460045c27fd8ee12259a1ae903 Mon Sep 17 00:00:00 2001 From: "Alan J. Wallcraft" Date: Fri, 22 Nov 2024 13:53:11 +0000 Subject: [PATCH 14/19] Make REMAPPING_USE_OM4_SUBCELLS the default In addition to REMAPPING_USE_OM4_SUBCELLS, for ALE remapping, there are several parameters of the form XXX_REMAPPING_USE_OM4_SUBCELLS, where XXX identifies the target, and they all currently default to True. To simplify setting them all to False, which is recommended, the defaults for the XXX versions is changed to the value of REMAPPING_USE_OM4_SUBCELLS. Answers are only changed if REMAPPING_USE_OM4_SUBCELLS is set to False and the default (now False) is used for one or more of the other parameters. In such cases the original behaviour can be recovered by explicitly setting the other parameters to True. --- src/core/MOM_open_boundary.F90 | 5 ++++- src/diagnostics/MOM_diagnostics.F90 | 5 ++++- src/framework/MOM_diag_mediator.F90 | 4 +++- src/initialization/MOM_state_initialization.F90 | 4 +++- src/initialization/MOM_tracer_initialization_from_Z.F90 | 4 +++- src/parameterizations/lateral/MOM_internal_tides.F90 | 4 +++- .../lateral/MOM_lateral_mixing_coeffs.F90 | 4 +++- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 8 ++++++-- src/tracer/MOM_hor_bnd_diffusion.F90 | 5 ++++- src/tracer/MOM_neutral_diffusion.F90 | 4 +++- 10 files changed, 36 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 22322fea82..a2cd1880fe 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -705,10 +705,13 @@ subroutine open_boundary_config(G, US, param_file, OBC) "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=OBC%om4_remap_via_sub_cells) endif ! OBC%number_of_segments > 0 diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 677fdfe6dc..abce27909b 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1616,10 +1616,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index c28e2e5896..3bb73e4c57 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3280,10 +3280,12 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "DIAG_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for diagnostics. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, & "The vintage of the expressions and order of arithmetic to use for remapping. "//& "Values below 20190101 result in the use of older, less accurate expressions "//& diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 629da43ddc..b566caa531 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2557,10 +2557,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq)) + call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 6e3da385ce..615ce07f0d 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -138,10 +138,12 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index c7101ac6b7..899dcbbbf0 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -3610,10 +3610,12 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "mode speeds are not calculated but are simply reported as 0. This must be "//& "non-negative for the wave_speeds routine to be used.", & units="m s-1", default=0.01, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "UNIFORM_TEST_CG", CS%uniform_test_cg, & "If positive, a uniform group velocity of internal tide for test case", & default=-1., units="m s-1", scale=US%m_s_to_L_T) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 8f388dc263..ac6efe2268 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1834,10 +1834,12 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call wave_speed_init(CS%wave_speed, GV, use_ebt_mode=CS%Resoln_use_ebt, & mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 0dfead633c..c00c72ea3c 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -236,10 +236,12 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& @@ -520,10 +522,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& "Dates before 20190101 give the same answers as the code did in late 2018, "//& diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index 6d8fe881d1..e2718590fb 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -143,10 +143,13 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "HBD_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for horizontal boundary diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping. call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 1aaf7409d2..feb5dde247 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -238,10 +238,12 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "NDIFF_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & From 88a2e4ffe8696166a7197a784e935f6fc76710d4 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 2 Dec 2024 15:51:47 -0500 Subject: [PATCH 15/19] Removed default for mandatory time scale in OBCs Removed two instances of `fail_if_missing=.true., default=0.` which are contradictory: a default value is meaningless if the parameter must be specified. I encountered this when adding the `defaults=` option to `get_param_real_array()`. --- src/core/MOM_open_boundary.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index a2cd1880fe..9f043473cb 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1479,7 +1479,7 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) "Timescales in days for nudging along a segment, "//& "for inflow, then outflow. Setting both to zero should "//& "behave like SIMPLE obcs for the baroclinic velocities.", & - fail_if_missing=.true., default=0., units="days", scale=86400.0*US%s_to_T) + fail_if_missing=.true., units="days", scale=86400.0*US%s_to_T) OBC%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1) OBC%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2) deallocate(tnudge) @@ -1620,7 +1620,7 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) "Timescales in days for nudging along a segment, "//& "for inflow, then outflow. Setting both to zero should "//& "behave like SIMPLE obcs for the baroclinic velocities.", & - fail_if_missing=.true., default=0., units="days", scale=86400.0*US%s_to_T) + fail_if_missing=.true., units="days", scale=86400.0*US%s_to_T) OBC%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1) OBC%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2) deallocate(tnudge) From a89258763f17e9f1d5bc6d954ff280fa3aca47a4 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 25 Nov 2024 20:11:21 -0500 Subject: [PATCH 16/19] Adds a vector of default values to get_param_real_array() The `default=` optional argument to get_param() only provides a uniform value to initialize an array of reals. This commit adds the optional `defaults=` argument that must have the same length as the `values` argument. I've also added a few instances of this optional argument: - by adding the `initialize_thickness_param()` procedure, selected by `THICKNESS_CONFIG = "param"`. The procedure was based on the "uniform" method, and uses the parameter `THICKNESS_INIT_VALUES` which defaults to uniform values derived from `MAXIMUM_DEPTH` - the setting of MLD_EN_VALS in MOM_diabatic_driver.F90 which was previously using a work around to set defaults to 25, 2500, 250000 J/m2. - two vectors of 4 values in user/user_change_diffusivity.F90 There will be some doc file changes, but no answer changes. --- src/framework/MOM_document.F90 | 11 +++- src/framework/MOM_file_parser.F90 | 22 +++++-- .../MOM_state_initialization.F90 | 65 +++++++++++++++++++ .../vertical/MOM_diabatic_driver.F90 | 9 +-- src/user/user_change_diffusivity.F90 | 5 +- 5 files changed, 95 insertions(+), 17 deletions(-) diff --git a/src/framework/MOM_document.F90 b/src/framework/MOM_document.F90 index f32573815f..eceb87d7d4 100644 --- a/src/framework/MOM_document.F90 +++ b/src/framework/MOM_document.F90 @@ -303,14 +303,16 @@ subroutine doc_param_real(doc, varname, desc, units, val, default, debuggingPara end subroutine doc_param_real !> This subroutine handles parameter documentation for arrays of reals. -subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debuggingParam, like_default) +subroutine doc_param_real_array(doc, varname, desc, units, vals, default, defaults, & + debuggingParam, like_default) type(doc_type), pointer :: doc !< A pointer to a structure that controls where the !! documentation occurs and its formatting character(len=*), intent(in) :: varname !< The name of the parameter being documented character(len=*), intent(in) :: desc !< A description of the parameter being documented character(len=*), intent(in) :: units !< The units of the parameter being documented real, intent(in) :: vals(:) !< The array of values to record - real, optional, intent(in) :: default !< The default value of this parameter + real, optional, intent(in) :: default !< A uniform default value of this parameter + real, optional, intent(in) :: defaults(:) !< The element-wise default values of this parameter logical, optional, intent(in) :: debuggingParam !< If present and true, this is a debugging parameter. logical, optional, intent(in) :: like_default !< If present and true, log this parameter as though !! it has the default value, even if there is no default. @@ -334,6 +336,11 @@ subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debugg do i=1,size(vals) ; if (vals(i) /= default) equalsDefault = .false. ; enddo mesg = trim(mesg)//" default = "//trim(real_string(default)) endif + if (present(defaults)) then + equalsDefault = .true. + do i=1,size(vals) ; if (vals(i) /= defaults(i)) equalsDefault = .false. ; enddo + mesg = trim(mesg)//" default = "//trim(real_array_string(defaults)) + endif if (present(like_default)) then ; if (like_default) equalsDefault = .true. ; endif if (mesgHasBeenDocumented(doc, varName, mesg)) return ! Avoid duplicates diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 22d3789ea5..fc496ac1b5 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -1464,7 +1464,7 @@ end subroutine log_param_real !> Log the name and values of an array of real model parameter in documentation files. subroutine log_param_real_array(CS, modulename, varname, value, desc, & - units, default, debuggingParam, like_default, unscale) + units, default, defaults, debuggingParam, like_default, unscale) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1473,7 +1473,8 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, & character(len=*), optional, intent(in) :: desc !< A description of this variable; if not !! present, this parameter is not written to a doc file character(len=*), intent(in) :: units !< The units of this parameter - real, optional, intent(in) :: default !< The default value of the parameter + real, optional, intent(in) :: default !< A uniform default value of the parameter + real, optional, intent(in) :: defaults(:) !< The element-wise defaults of the parameter logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file logical, optional, intent(in) :: like_default !< If present and true, log this parameter as @@ -1498,7 +1499,7 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, & write(myunits(1:240),'(A)') trim(units) if (present(desc)) & - call doc_param(CS%doc, varname, desc, myunits, log_val, default, & + call doc_param(CS%doc, varname, desc, myunits, log_val, default, defaults, & debuggingParam=debuggingParam, like_default=like_default) end subroutine log_param_real_array @@ -1835,7 +1836,7 @@ end subroutine get_param_real !> This subroutine reads the values of an array of real model parameters from a parameter file !! and logs them in documentation files. subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & - default, fail_if_missing, do_not_read, do_not_log, debuggingParam, & + default, defaults, fail_if_missing, do_not_read, do_not_log, debuggingParam, & scale, unscaled) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters @@ -1846,7 +1847,8 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & character(len=*), optional, intent(in) :: desc !< A description of this variable; if not !! present, this parameter is not written to a doc file character(len=*), intent(in) :: units !< The units of this parameter - real, optional, intent(in) :: default !< The default value of the parameter + real, optional, intent(in) :: default !< A uniform default value of the parameter + real, optional, intent(in) :: defaults(:) !< The element-wise defaults of the parameter logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file logical, optional, intent(in) :: do_not_read !< If present and true, do not read a @@ -1865,14 +1867,22 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log + if (present(defaults)) then + if (present(default)) call MOM_error(FATAL, & + "get_param_real_array: Only one of default and defaults can be specified at a time.") + if (size(defaults) /= size(value)) call MOM_error(FATAL, & + "get_param_real_array: The size of defaults nad value are not the same.") + endif + if (do_read) then if (present(default)) value(:) = default + if (present(defaults)) value(:) = defaults(:) call read_param_real_array(CS, varname, value, fail_if_missing) endif if (do_log) then call log_param_real_array(CS, modulename, varname, value, desc, & - units, default, debuggingParam) + units, default, defaults, debuggingParam) endif if (present(unscaled)) unscaled(:) = value(:) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index b566caa531..41b407d6a1 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -278,6 +278,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & " \t uniform - uniform thickness layers evenly distributed \n"//& " \t\t between the surface and MAXIMUM_DEPTH. \n"//& " \t list - read a list of positive interface depths. \n"//& + " \t param - use thicknesses from parameter THICKNESS_INIT_VALUES. \n"//& " \t DOME - use a slope and channel configuration for the \n"//& " \t\t DOME sill-overflow test case. \n"//& " \t ISOMIP - use a configuration for the \n"//& @@ -318,6 +319,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & just_read=just_read) case ("list"); call initialize_thickness_list(dz, depth_tot, G, GV, US, PF, & just_read=just_read) + case ("param"); call initialize_thickness_param(dz, depth_tot, G, GV, US, PF, & + just_read=just_read) case ("DOME"); call DOME_initialize_thickness(dz, depth_tot, G, GV, PF, & just_read=just_read) case ("ISOMIP"); call ISOMIP_initialize_thickness(dz, depth_tot, G, GV, US, PF, tv, & @@ -1011,6 +1014,68 @@ subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_r call callTree_leave(trim(mdl)//'()') end subroutine initialize_thickness_list +!> Initializes thickness based on a run-time parameter with nominal thickness +!! for each layer +subroutine initialize_thickness_param(h, depth_tot, G, GV, US, param_file, just_read) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: just_read !< If true, this call will only read + !! parameters without changing h. + ! Local variables + character(len=40) :: mdl = "initialize_thickness_param" ! This subroutine's name. + real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually + ! negative because it is positive upward. + real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface, + ! positive upward [Z ~> m]. + real :: dz(SZK_(GV)) ! The nominal initial layer thickness [Z ~> m], usually + real :: h0_def(SZK_(GV)) ! Uniform default values for dz [Z ~> m], usually + integer :: i, j, k, is, ie, js, je, nz + + call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90") + if (G%max_depth<=0.) call MOM_error(FATAL, "initialize_thickness_param: "// & + "MAXIMUM_DEPTH has a nonsensical value! Was it set?") + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + h0_def(:) = ( G%max_depth / real(nz) ) * US%Z_to_m + call get_param(param_file, mdl, "THICKNESS_INIT_VALUES", dz, & + "A list of nominal thickness for each layer to initialize with", & + units="m", scale=US%m_to_Z, defaults=h0_def, do_not_log=just_read) + if (just_read) return ! This subroutine has no run-time parameters. + + e0(nz+1) = -G%max_depth + do k=nz, 1, -1 + e0(K) = e0(K+1) + dz(k) + enddo + + do j=js,je ; do i=is,ie + ! This sets the initial thickness (in m) of the layers. The + ! thicknesses are set to insure that: 1. each layer is at least an + ! Angstrom thick, and 2. the interfaces are where they should be + ! based on the resting depths and interface height perturbations, + ! as long at this doesn't interfere with 1. + eta1D(nz+1) = -depth_tot(i,j) + do k=nz,1,-1 + eta1D(K) = e0(K) + if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then + eta1D(K) = eta1D(K+1) + GV%Angstrom_Z + h(i,j,k) = GV%Angstrom_Z + else + h(i,j,k) = eta1D(K) - eta1D(K+1) + endif + enddo + enddo ; enddo + + call callTree_leave(trim(mdl)//'()') +end subroutine initialize_thickness_param + !> Search density space for location of layers (not implemented!) subroutine initialize_thickness_search call MOM_error(FATAL," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED") diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c5297a3cf0..e984d5831c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3256,13 +3256,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m) call get_param(param_file, mdl, "MLD_EN_VALS", CS%MLD_En_vals, & "The energy values used to compute MLDs. If not set (or all set to 0.), the "//& - "default will overwrite to 25., 2500., 250000.", & - units='J/m2', default=0., scale=US%W_m2_to_RZ3_T3*US%s_to_T) - if ((CS%MLD_En_vals(1)==0.).and.(CS%MLD_En_vals(2)==0.).and.(CS%MLD_En_vals(3)==0.)) then - CS%MLD_En_vals = (/ 25.*US%W_m2_to_RZ3_T3*US%s_to_T, & - 2500.*US%W_m2_to_RZ3_T3*US%s_to_T, & - 250000.*US%W_m2_to_RZ3_T3*US%s_to_T /) - endif + "default will overwrite to 25., 2500., 250000.", units='J/m2', & + defaults=(/25., 2500., 250000./), scale=US%W_m2_to_RZ3_T3*US%s_to_T) write(EN1,'(F10.2)') CS%MLD_En_vals(1)*US%RZ3_T3_to_W_m2*US%T_to_s write(EN2,'(F10.2)') CS%MLD_En_vals(2)*US%RZ3_T3_to_W_m2*US%T_to_s write(EN3,'(F10.2)') CS%MLD_En_vals(3)*US%RZ3_T3_to_W_m2*US%T_to_s diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index 9a56c12b9c..1a1881a42b 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -230,14 +230,15 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS) "applied. The four values specify the latitudes at "//& "which the extra diffusivity starts to increase from 0, "//& "hits its full value, starts to decrease again, and is "//& - "back to 0.", units="degrees_N", default=-1.0e9) + "back to 0.", units="degrees_N", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/)) call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", CS%rho_range(:), & "Four successive values that define a range of potential "//& "densities over which the user-given extra diffusivity "//& "is applied. The four values specify the density at "//& "which the extra diffusivity starts to increase from 0, "//& "hits its full value, starts to decrease again, and is "//& - "back to 0.", units="kg m-3", default=-1.0e9, scale=US%kg_m3_to_R) + "back to 0.", units="kg m-3", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/),& + scale=US%kg_m3_to_R) call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", CS%use_abs_lat, & "If true, use the absolute value of latitude when "//& "checking whether a point fits into range of latitudes.", & From 822257e5f6b9f9883b2180c0ca83d9f9db27632a Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sun, 8 Dec 2024 08:06:30 -0500 Subject: [PATCH 17/19] FMS API: Convert real kind of constants Two latent heat constants are imported directly from FMS, which is built independently of MOM6. Previously, it was a safe assumption that both would be built with double precision, but this is no longer the case since FMS now supports both single and double precision. This could cause conflicts with mixed-precision builds. This patch converts the values from FMS-precision to MOM-precision. Single->double should not affect reproducibility since every single-precision number can be exactly represented in double precision. Double->single could affect reproducibility, but this is not an issue since MOM6 does not run in single precision. --- config_src/infra/FMS1/MOM_constants.F90 | 10 +++++++--- config_src/infra/FMS2/MOM_constants.F90 | 10 +++++++--- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/config_src/infra/FMS1/MOM_constants.F90 b/config_src/infra/FMS1/MOM_constants.F90 index 2db177e08c..a632267a7f 100644 --- a/config_src/infra/FMS1/MOM_constants.F90 +++ b/config_src/infra/FMS1/MOM_constants.F90 @@ -3,12 +3,16 @@ module MOM_constants ! This file is part of MOM6. See LICENSE.md for the license. -use constants_mod, only : HLV, HLF +use constants_mod, only : FMS_HLV => HLV +use constants_mod, only : FMS_HLF => HLF implicit none ; private -!> The constant offset for converting temperatures in Kelvin to Celsius real, public, parameter :: CELSIUS_KELVIN_OFFSET = 273.15 -public :: HLV, HLF + !< The constant offset for converting temperatures in Kelvin to Celsius [K] +real, public, parameter :: HLV = real(FMS_HLV, kind=kind(1.0)) + !< Latent heat of vaporization [J kg-1] +real, public, parameter :: HLF = real(FMS_HLF, kind=kind(1.0)) + !< Latent heat of fusion [J kg-1] end module MOM_constants diff --git a/config_src/infra/FMS2/MOM_constants.F90 b/config_src/infra/FMS2/MOM_constants.F90 index 2db177e08c..a632267a7f 100644 --- a/config_src/infra/FMS2/MOM_constants.F90 +++ b/config_src/infra/FMS2/MOM_constants.F90 @@ -3,12 +3,16 @@ module MOM_constants ! This file is part of MOM6. See LICENSE.md for the license. -use constants_mod, only : HLV, HLF +use constants_mod, only : FMS_HLV => HLV +use constants_mod, only : FMS_HLF => HLF implicit none ; private -!> The constant offset for converting temperatures in Kelvin to Celsius real, public, parameter :: CELSIUS_KELVIN_OFFSET = 273.15 -public :: HLV, HLF + !< The constant offset for converting temperatures in Kelvin to Celsius [K] +real, public, parameter :: HLV = real(FMS_HLV, kind=kind(1.0)) + !< Latent heat of vaporization [J kg-1] +real, public, parameter :: HLF = real(FMS_HLF, kind=kind(1.0)) + !< Latent heat of fusion [J kg-1] end module MOM_constants From 8834f68fd0123ae7c8b1e9056b49e516028ed140 Mon Sep 17 00:00:00 2001 From: Chengzhu Xu <135884058+c2xu@users.noreply.github.com> Date: Mon, 9 Dec 2024 17:20:54 -0400 Subject: [PATCH 18/19] Inline harmonic analysis (#744) * Inline harmonic analysis Important bug fix: 1) The Cholesky decomposition was operating on entries below the main diagonal of FtF, whereas in the accumulator of FtF, only entries along and above the main diagonal were calculated. In this revision, I modified HA_accum_FtF so that entries below the main diagonal are accumulated instead. 2) In the accumulator of FtSSH, the first entry for the mean (zero frequency) is moved out of the loop over different tidal constituents, so that it is not accumulated multiple times within a single time step. * Inline harmonic analysis Another bug fix: initial state added back to the mean state. * Inline harmonic analysis Minor update to HA_solver --- src/diagnostics/MOM_harmonic_analysis.F90 | 103 ++++++++++++---------- 1 file changed, 55 insertions(+), 48 deletions(-) diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 index 1e3b9895cb..0013a8ab83 100644 --- a/src/diagnostics/MOM_harmonic_analysis.F90 +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -173,6 +173,7 @@ end subroutine HA_register !> This subroutine accumulates the temporal basis functions in FtF. !! The tidal constituents are those used in MOM_tidal_forcing, plus the mean (of zero frequency). +!! Only the main diagonal and entries below it are calculated, which are needed for Cholesky decomposition. subroutine HA_accum_FtF(Time, CS) type(time_type), intent(in) :: Time !< The current model time type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module @@ -191,27 +192,31 @@ subroutine HA_accum_FtF(Time, CS) nc = CS%nc now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) - ! Accumulate FtF - CS%FtF(1,1) = CS%FtF(1,1) + 1.0 !< For the zero frequency + !< First entry, corresponding to the zero frequency constituent (mean) + CS%FtF(1,1) = CS%FtF(1,1) + 1.0 + do c=1,nc icos = 2*c isin = 2*c+1 cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + + ! First column, corresponding to the zero frequency constituent (mean) CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat - CS%FtF(1,icos) = CS%FtF(icos,1) - CS%FtF(1,isin) = CS%FtF(isin,1) - do cc=c,nc + + do cc=1,c iccos = 2*cc issin = 2*cc+1 ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc)) ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc)) + + ! Interior of the matrix, corresponding to the products of cosine and sine terms CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat CS%FtF(isin,issin) = CS%FtF(isin,issin) + sinomegat * ssinomegat - enddo ! cc=c,nc + enddo ! cc=1,c enddo ! c=1,nc end subroutine HA_accum_FtF @@ -276,14 +281,18 @@ subroutine HA_accum_FtSSH(key, data, Time, G, CS) is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je - ! Accumulate FtF and FtSSH + !< First entry, corresponding to the zero frequency constituent (mean) + do j=js,je ; do i=is,ie + ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) + enddo ; enddo + + !< The remaining entries do c=1,nc icos = 2*c isin = 2*c+1 cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) do j=js,je ; do i=is,ie - ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat enddo ; enddo @@ -315,7 +324,7 @@ subroutine HA_write(ha1, Time, G, CS) ! Local variables real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A] integer :: year, month, day, hour, minute, second - integer :: nc, k, is, ie, js, je + integer :: nc, i, j, k, is, ie, js, je character(len=255) :: filename !< Output file name type(MOM_infra_file) :: cdf !< The file handle for output harmonic constants @@ -348,6 +357,11 @@ subroutine HA_write(ha1, Time, G, CS) call create_MOM_file(cdf, trim(filename), cdf_vars, & 2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G) + ! Add the initial field back to the mean state + do j=js,je ; do i=is,ie + FtSSHw(i,j,1) = FtSSHw(i,j,1) + ha1%ref(i,j) + enddo ; enddo + ! Write data call MOM_write_field(cdf, cdf_fields(1), G%domain, FtSSHw(:,:,1), 0.0) do k=1,nc @@ -362,75 +376,68 @@ subroutine HA_write(ha1, Time, G, CS) end subroutine HA_write -!> This subroutine computes the harmonic constants (stored in FtSSHw) using the dot products of the temporal +!> This subroutine computes the harmonic constants (stored in x) using the dot products of the temporal !! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis !! functions accumulated in FtSSH. The system is solved by Cholesky decomposition, !! -!! FtF * FtSSHw = FtSSH, => FtFw * (FtFw' * FtSSHw) = FtSSH, +!! FtF * x = FtSSH, => L * (L' * x) = FtSSH, => L * y = FtSSH, !! -!! where FtFw is a lower triangular matrix, and the prime denotes matrix transpose. +!! where L is the lower triangular matrix, y = L' * x, and x is the solution vector. !! -subroutine HA_solver(ha1, nc, FtF, FtSSHw) +subroutine HA_solver(ha1, nc, FtF, x) type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field integer, intent(in) :: nc !< Number of harmonic constituents real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim] - real, dimension(:,:,:), allocatable, intent(out) :: FtSSHw !< Work array for Cholesky decomposition [A] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1), & + intent(out) :: x !< Solution vector of harmonic constants [A] ! Local variables - real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim] - real, dimension(:), allocatable :: tmp1 !< Temporary variable for Cholesky decomposition [nondim] - real, dimension(:,:), allocatable :: tmp2 !< Temporary variable for Cholesky decomposition [A] - real, dimension(:,:), allocatable :: FtFw !< Lower triangular matrix for Cholesky decomposition [nondim] - integer :: k, m, n, is, ie, js, je - - is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je - - allocate(tmp1(1:2*nc+1), source=0.0) - allocate(tmp2(is:ie,js:je), source=0.0) - allocate(FtFw(1:2*nc+1,1:2*nc+1), source=0.0) - allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) - - ! Construct FtFw - FtFw(:,:) = 0.0 + real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim] + real, dimension(2*nc+1,2*nc+1) :: L !< Lower triangular matrix of Cholesky decomposition [nondim] + real, dimension(2*nc+1) :: tmp1 !< Inverse of the diagonal entries of L [nondim] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je) :: tmp2 !< 2D temporary array involving FtSSH [A] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1) :: y !< 3D temporary array, i.e., L' * x [A] + integer :: k, m, n + + ! Cholesky decomposition do m=1,2*nc+1 + + ! First, calculate the diagonal entries tmp0 = 0.0 - do k=1,m-1 - tmp0 = tmp0 + FtFw(m,k) * FtFw(m,k) + do k=1,m-1 ! This loop operates along the m-th row + tmp0 = tmp0 + L(m,k) * L(m,k) enddo - FtFw(m,m) = sqrt(FtF(m,m) - tmp0) - tmp1(m) = 1 / FtFw(m,m) - do k=m+1,2*nc+1 + L(m,m) = sqrt(FtF(m,m) - tmp0) ! This is the m-th diagonal entry + + ! Now calculate the off-diagonal entries + tmp1(m) = 1 / L(m,m) + do k=m+1,2*nc+1 ! This loop operates along the column below the m-th diagonal entry tmp0 = 0.0 do n=1,m-1 - tmp0 = tmp0 + FtFw(k,n) * FtFw(m,n) + tmp0 = tmp0 + L(k,n) * L(m,n) enddo - FtFw(k,m) = (FtF(k,m) - tmp0) * tmp1(m) + L(k,m) = (FtF(k,m) - tmp0) * tmp1(m) ! This is the k-th off-diagonal entry below the m-th diagonal entry enddo enddo - ! Solve for (FtFw' * FtSSHw) - FtSSHw(:,:,:) = ha1%FtSSH(:,:,:) + ! Solve for y from L * y = FtSSH do k=1,2*nc+1 tmp2(:,:) = 0.0 do m=1,k-1 - tmp2(:,:) = tmp2(:,:) + FtFw(k,m) * FtSSHw(:,:,m) + tmp2(:,:) = tmp2(:,:) + L(k,m) * y(:,:,m) enddo - FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + y(:,:,k) = (ha1%FtSSH(:,:,k) - tmp2(:,:)) * tmp1(k) enddo - ! Solve for FtSSHw + ! Solve for x from L' * x = y do k=2*nc+1,1,-1 tmp2(:,:) = 0.0 do m=k+1,2*nc+1 - tmp2(:,:) = tmp2(:,:) + FtSSHw(:,:,m) * FtFw(m,k) + tmp2(:,:) = tmp2(:,:) + L(m,k) * x(:,:,m) enddo - FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + x(:,:,k) = (y(:,:,k) - tmp2(:,:)) * tmp1(k) enddo - deallocate(tmp1) - deallocate(tmp2) - deallocate(FtFw) - end subroutine HA_solver !> \namespace harmonic_analysis @@ -441,7 +448,7 @@ end subroutine HA_solver !! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent !! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares, !! -!! (F' * F) * x = F' * SSH_in, +!! (F' * F) * x = F' * SSH_in, => FtF * x = FtSSH, !! !! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by !! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is From ee42fdba03804df9faa12584a8575cbf7f126dbe Mon Sep 17 00:00:00 2001 From: Alan Wallcraft Date: Mon, 9 Dec 2024 17:37:21 -0500 Subject: [PATCH 19/19] Tidal angular frequency has units [rad s-1] (#764) * Tidal angular frequency has units [rad s-1] Tidal frequencies are always angular frequencies to simplify applying sine and cosine. These have MKS units [rad s-1] but they are all currently listed as [s-1]. Updated dOxygen comments for variables, e.g. [T-1 ~> s-1] becomes [rad T-1 ~> rad s-1]. Updated get_param units. e.g. units="s-1" becomes units="rad s-1". No answers are changed, but the logged parameter units are different. There are frequencies in MOM_internal_tides.F90 but these have not been updated because they may be specified incorrectly. They are used as if they are [T-1] but they are calculated as 2PI/period [rad T-1]. real, allocatable, dimension(:) :: frequency !< The frequency of each band [T-1 ~> s-1]. real :: period ! A tidal period read from namelist [T ~> s] ! The periods of the tidal constituents for internal tides raytracing call read_param(param_file, "TIDAL_PERIODS", periods) do fr=1,num_freq period = US%s_to_T*extract_real(periods, " ,", fr, 0.) CS%frequency(fr) = 8.0*atan(1.0)/period enddo All MOM6-examples cases have INTERNAL_TIDES=False and so can't resolve this issue. * fixed too-long line --- src/core/MOM_barotropic.F90 | 6 +++--- src/core/MOM_open_boundary.F90 | 5 +++-- src/parameterizations/lateral/MOM_streaming_filter.F90 | 4 ++-- src/parameterizations/lateral/MOM_tidal_forcing.F90 | 9 +++++---- 4 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index bee598b8d6..aa502bdd4b 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -5323,7 +5323,7 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) character(len=40) :: mdl = "MOM_barotropic" ! This module's name. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB real :: am2, ak1 !< Bandwidth parameters of the M2 and K1 streaming filters [nondim] - real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [T-1 ~> s-1] + real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [rad T-1 ~> rad s-1] isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB @@ -5354,13 +5354,13 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) "Frequency of the M2 tidal constituent. "//& "This is only used if TIDES and TIDE_M2"// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and M2"// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency("M2"), & + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("M2"), & scale=US%T_to_s, do_not_log=.true.) call get_param(param_file, mdl, "TIDE_K1_FREQ", ok1, & "Frequency of the K1 tidal constituent. "//& "This is only used if TIDES and TIDE_K1"// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and K1"// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency("K1"), & + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("K1"), & scale=US%T_to_s, do_not_log=.true.) ALLOC_(CS%ubtav(IsdB:IedB,jsd:jed)) ; CS%ubtav(:,:) = 0.0 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 9f043473cb..b15d28194e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -320,7 +320,8 @@ module MOM_open_boundary logical :: add_tide_constituents = .false. !< If true, add tidal constituents to the boundary elevation !! and velocity. Will be set to true if n_tide_constituents > 0. character(len=2), allocatable, dimension(:) :: tide_names !< Names of tidal constituents to add to the boundary data. - real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents [T-1 ~> s-1]. + real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal + !! constituents [rad T-1 ~> rad s-1]. real, allocatable, dimension(:) :: tide_eq_phases !< Equilibrium phases of chosen tidal constituents [rad]. real, allocatable, dimension(:) :: tide_fn !< Amplitude modulation of boundary tides by nodal cycle [nondim]. real, allocatable, dimension(:) :: tide_un !< Phase modulation of boundary tides by nodal cycle [rad]. @@ -1234,7 +1235,7 @@ subroutine initialize_obc_tides(OBC, US, param_file) "This is only used if TIDES and TIDE_"//trim(OBC%tide_names(c))// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(OBC%tide_names(c))//& " is in OBC_TIDE_CONSTITUENTS.", & - units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s) + units="rad s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s) ! Find equilibrium phase if needed if (OBC%add_eq_phase) then diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index a91f6661f2..f7a4f736c7 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -15,7 +15,7 @@ module MOM_streaming_filter !> The control structure for storing the filter infomation of a particular field type, public :: Filter_CS ; private real :: a, & !< Parameter that determines the bandwidth [nondim] - om, & !< Target frequency of the filter [T-1 ~> s-1] + om, & !< Target frequency of the filter [rad T-1 ~> rad s-1] old_time = -1.0 !< The time of the previous accumulating step [T ~> s] real, allocatable, dimension(:,:) :: s1, & !< Dummy variable [A] u1 !< Filtered data [A] @@ -29,7 +29,7 @@ module MOM_streaming_filter !> This subroutine registers each of the fields to be filtered. subroutine Filt_register(a, om, grid, HI, CS) real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] - real, intent(in) :: om !< Target frequency of the filter [T-1 ~> s-1] + real, intent(in) :: om !< Target frequency of the filter [rad T-1 ~> rad s-1] character(len=*), intent(in) :: grid !< Horizontal grid location: h, u, or v type(hor_index_type), intent(in) :: HI !< Horizontal index type structure type(Filter_CS), intent(out) :: CS !< Control structure for the current field diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 177becf84f..43885cccc3 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -50,7 +50,7 @@ module MOM_tidal_forcing !! and bottom geopotential anomalies [nondim]. integer :: nc !< The number of tidal constituents in use. real, dimension(MAX_CONSTITUENTS) :: & - freq, & !< The frequency of a tidal constituent [T-1 ~> s-1]. + freq, & !< The frequency of a tidal constituent [rad T-1 ~> rad s-1]. phase0, & !< The phase of a tidal constituent at time 0 [rad]. amp, & !< The amplitude of a tidal constituent at time 0 [Z ~> m]. love_no !< The Love number of a tidal constituent at time 0 [nondim]. @@ -151,7 +151,7 @@ end function eq_phase !! Values used here are from previous versions of MOM. function tidal_frequency(constit) character (len=2), intent(in) :: constit !> Constituent to look up - real :: tidal_frequency !> Angular frequency [s-1] + real :: tidal_frequency !> Angular frequency [rad s-1] select case (constit) case ("M2") @@ -246,7 +246,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) phase, & ! The phase of some tidal constituent [radians]. lat_rad, lon_rad ! Latitudes and longitudes of h-points [radians]. real :: deg_to_rad ! A conversion factor from degrees to radians [radian degree-1] - real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [s-1] + real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [rad s-1] real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad] real, dimension(MAX_CONSTITUENTS) :: amp_def ! Default amplitude for each tidal constituent [m] real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim] @@ -480,7 +480,8 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) "Frequency of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(CS%const_name(c))// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c), scale=US%T_to_s) + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=freq_def(c), & + scale=US%T_to_s) call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_AMP", CS%amp(c), & "Amplitude of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &