Skip to content

Commit

Permalink
Merge branch 'NOAA-GFDL:dev/gfdl' into c2xu/nodal_modulation
Browse files Browse the repository at this point in the history
  • Loading branch information
@c2xu
c2xu committed Sep 14, 2024
2 parents 211cab5 + 5fc90eb commit 241f7b7
Show file tree
Hide file tree
Showing 15 changed files with 2,002 additions and 579 deletions.
13 changes: 9 additions & 4 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2890,7 +2890,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
endif

if (.not. CS%adiabatic) then
call register_diabatic_restarts(G, US, param_file, CS%int_tide_CSp, restart_CSp)
call register_diabatic_restarts(G, GV, US, param_file, CS%int_tide_CSp, restart_CSp)
endif

call callTree_waypoint("restart registration complete (initialize_MOM)")
Expand Down Expand Up @@ -2955,20 +2955,20 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
! These arrays are not initialized in most solo cases, but are needed
! when using an ice shelf. Passing the ice shelf diagnostics CS from MOM
! for legacy reasons. The actual ice shelf diag CS is internal to the ice shelf
call initialize_ice_shelf(param_file, G_in, Time, ice_shelf_CSp, diag_ptr, &
call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr, &
Time_init, dirs%output_directory, calve_ice_shelf_bergs=point_calving)
allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0)
allocate(mass_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0)
allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0)
allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0)
call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf)
call ice_shelf_query(ice_shelf_CSp, G, CS%frac_shelf_h, CS%mass_shelf)
! MOM_initialize_state is using the unrotated metric
call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in)
call rotate_array(CS%mass_shelf, -turns, mass_shelf_in)
call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, &
frac_shelf_h=frac_shelf_in, mass_shelf = mass_shelf_in)
frac_shelf_h=frac_shelf_in, mass_shelf=mass_shelf_in)
else
call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
Expand Down Expand Up @@ -4173,9 +4173,14 @@ subroutine save_MOM_restart(CS, directory, time, G, time_stamped, filename, &
logical, optional, intent(in) :: write_IC
!< If present and true, initial conditions are being written

logical :: showCallTree
showCallTree = callTree_showQuery()

if (showCallTree) call callTree_waypoint("About to call save_restart (step_MOM)")
call save_restart(directory, time, G, CS%restart_CS, &
time_stamped=time_stamped, filename=filename, GV=GV, &
num_rest_files=num_rest_files, write_IC=write_IC)
if (showCallTree) call callTree_waypoint("Done with call to save_restart (step_MOM)")

if (CS%use_particles) call particles_save_restart(CS%particles, CS%h, directory, time, time_stamped)
end subroutine save_MOM_restart
Expand Down
32 changes: 29 additions & 3 deletions src/core/MOM_PressureForce_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,9 @@ module MOM_PressureForce_FV
!! for the finite volume pressure gradient calculation.
!! By the default (1) is for a piecewise linear method

logical :: use_SSH_in_Z0p !< If true, adjust the height at which the pressure used in the
!! equation of state is 0 to account for the displacement of the sea
!! surface including adjustments for atmospheric or sea-ice pressure.
logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
integer :: tides_answer_date !< Recover old answers with tides in Boussinesq mode
integer :: id_e_tide = -1 !< Diagnostic identifier
Expand Down Expand Up @@ -522,6 +525,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! [Z ~> m].
e_tide_sal, & ! The bottom geopotential anomaly due to harmonic self-attraction and loading
! specific to tides [Z ~> m].
Z_0p, & ! The height at which the pressure used in the equation of state is 0 [Z ~> m]
SSH, & ! The sea surface height anomaly, in depth units [Z ~> m].
dM ! The barotropic adjustment to the Montgomery potential to
! account for a reduced gravity model [L2 T-2 ~> m2 s-2].
Expand Down Expand Up @@ -577,6 +581,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
! in roundoff and can be neglected [H ~> m].
real :: I_Rho0 ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1].
real :: G_Rho0 ! G_Earth / Rho0 in [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1].
real :: I_g_rho ! The inverse of the density times the gravitational acceleration [Z T2 L-2 R-1 ~> m Pa-1]
real :: rho_ref ! The reference density [R ~> kg m-3].
real :: dz_neglect ! A minimal thickness [Z ~> m], like e.
real :: H_to_RL2_T2 ! A factor to convert from thickness units (H) to pressure
Expand Down Expand Up @@ -753,6 +758,21 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
enddo ; enddo
endif

if (CS%use_SSH_in_Z0p .and. use_p_atm) then
I_g_rho = 1.0 / (CS%rho0*GV%g_Earth)
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = e(i,j,1) + p_atm(i,j) * I_g_rho
enddo ; enddo
elseif (CS%use_SSH_in_Z0p) then
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = e(i,j,1)
enddo ; enddo
else
do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
Z_0p(i,j) = G%Z_ref
enddo ; enddo
endif

do k=1,nz
! Calculate 4 integrals through the layer that are required in the
! subsequent calculation.
Expand All @@ -769,19 +789,19 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
intx_dpa(:,:,k), inty_dpa(:,:,k), &
MassWghtInterp=CS%MassWghtInterp, &
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=Z_0p)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa(:,:,k), intz_dpa(:,:,k), &
intx_dpa(:,:,k), inty_dpa(:,:,k), &
MassWghtInterp=CS%MassWghtInterp, Z_0p=G%Z_ref)
MassWghtInterp=CS%MassWghtInterp, Z_0p=Z_0p)
endif
else
call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), &
rho_ref, CS%Rho0, GV%g_Earth, G%HI, tv%eqn_of_state, US, dpa(:,:,k), &
intz_dpa(:,:,k), intx_dpa(:,:,k), inty_dpa(:,:,k), G%bathyT, dz_neglect, &
CS%MassWghtInterp, Z_0p=G%Z_ref)
CS%MassWghtInterp, Z_0p=Z_0p)
endif
if (GV%Z_to_H /= 1.0) then
!$OMP parallel do default(shared)
Expand Down Expand Up @@ -1042,6 +1062,12 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, SAL_CSp,
endif
call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
"If true, calculate self-attraction and loading.", default=CS%tides)
call get_param(param_file, mdl, "SSH_IN_EOS_PRESSURE_FOR_PGF", CS%use_SSH_in_Z0p, &
"If true, include contributions from the sea surface height in the height-based "//&
"pressure used in the equation of state calculations for the Boussinesq pressure "//&
"gradient forces, including adjustments for atmospheric or sea-ice pressure.", &
default=.false., do_not_log=.not.GV%Boussinesq)

call get_param(param_file, "MOM", "USE_REGRIDDING", use_ALE, &
"If True, use the ALE algorithm (regridding/remapping). "//&
"If False, use the layered isopycnal algorithm.", default=.false. )
Expand Down
70 changes: 70 additions & 0 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ module MOM_barotropic
use MOM_restart, only : query_initialized, MOM_restart_CS
use MOM_self_attr_load, only : scalar_SAL_sensitivity
use MOM_self_attr_load, only : SAL_CS
use MOM_streaming_filter, only : Filt_register, Filt_accum, Filter_CS
use MOM_tidal_forcing, only : tidal_frequency
use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-)
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : BT_cont_type, alloc_bt_cont_type
Expand Down Expand Up @@ -248,6 +250,10 @@ module MOM_barotropic
logical :: linearized_BT_PV !< If true, the PV and interface thicknesses used
!! in the barotropic Coriolis calculation is time
!! invariant and linearized.
logical :: use_filter_m2 !< If true, apply streaming band-pass filter for detecting
!! instantaneous tidal signals.
logical :: use_filter_k1 !< If true, apply streaming band-pass filter for detecting
!! instantaneous tidal signals.
logical :: use_wide_halos !< If true, use wide halos and march in during the
!! barotropic time stepping for efficiency.
logical :: clip_velocity !< If true, limit any velocity components that are
Expand Down Expand Up @@ -291,6 +297,10 @@ module MOM_barotropic
type(hor_index_type), pointer :: debug_BT_HI => NULL() !< debugging copy of horizontal index_type
type(SAL_CS), pointer :: SAL_CSp => NULL() !< Control structure for SAL
type(harmonic_analysis_CS), pointer :: HA_CSp => NULL() !< Control structure for harmonic analysis
type(Filter_CS) :: Filt_CS_um2, & !< Control structures for the M2 streaming filter
Filt_CS_vm2, & !< Control structures for the M2 streaming filter
Filt_CS_uk1, & !< Control structures for the K1 streaming filter
Filt_CS_vk1 !< Control structures for the K1 streaming filter
logical :: module_is_initialized = .false. !< If true, module has been initialized

integer :: isdw !< The lower i-memory limit for the wide halo arrays.
Expand Down Expand Up @@ -598,6 +608,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
DCor_v, & ! An averaged total thickness at v points [H ~> m or kg m-2].
Datv ! Basin depth at v-velocity grid points times the x-grid
! spacing [H L ~> m2 or kg m-1].
real, dimension(:,:), pointer :: um2, uk1, vm2, vk1
! M2 and K1 velocities from the output of streaming filters [m s-1]
real, target, dimension(SZIW_(CS),SZJW_(CS)) :: &
eta, & ! The barotropic free surface height anomaly or column mass
! anomaly [H ~> m or kg m-2]
Expand Down Expand Up @@ -1586,6 +1598,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
endif ; enddo ; enddo
endif

! Here is an example of how the filter equations are time stepped to determine the M2 and K1 velocities.
! The filters are initialized and registered in subroutine barotropic_init.
if (CS%use_filter_m2) then
call Filt_accum(ubt, um2, CS%Time, US, CS%Filt_CS_um2)
call Filt_accum(vbt, vm2, CS%Time, US, CS%Filt_CS_vm2)
endif
if (CS%use_filter_k1) then
call Filt_accum(ubt, uk1, CS%Time, US, CS%Filt_CS_uk1)
call Filt_accum(vbt, vk1, CS%Time, US, CS%Filt_CS_vk1)
endif

! Zero out the arrays for various time-averaged quantities.
if (find_etaav) then
!$OMP do
Expand Down Expand Up @@ -5247,6 +5270,8 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
type(vardesc) :: vd(3)
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]

isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed
IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB
Expand All @@ -5259,6 +5284,33 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
"sum(u dh_dt) while also correcting for truncation errors.", &
default=.false., do_not_log=.true.)

call get_param(param_file, mdl, "STREAMING_FILTER_M2", CS%use_filter_m2, &
"If true, turn on streaming band-pass filter for detecting "//&
"instantaneous tidal signals.", default=.false.)
call get_param(param_file, mdl, "STREAMING_FILTER_K1", CS%use_filter_k1, &
"If true, turn on streaming band-pass filter for detecting "//&
"instantaneous tidal signals.", default=.false.)
call get_param(param_file, mdl, "FILTER_ALPHA_M2", am2, &
"Bandwidth parameter of the streaming filter targeting the M2 frequency. "//&
"Must be positive. To turn off filtering, set FILTER_ALPHA_M2 <= 0.0.", &
default=0.0, units="nondim", do_not_log=.not.CS%use_filter_m2)
call get_param(param_file, mdl, "FILTER_ALPHA_K1", ak1, &
"Bandwidth parameter of the streaming filter targeting the K1 frequency. "//&
"Must be positive. To turn off filtering, set FILTER_ALPHA_K1 <= 0.0.", &
default=0.0, units="nondim", do_not_log=.not.CS%use_filter_k1)
call get_param(param_file, mdl, "TIDE_M2_FREQ", om2, &
"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"), &
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"), &
scale=US%T_to_s, do_not_log=.true.)

ALLOC_(CS%ubtav(IsdB:IedB,jsd:jed)) ; CS%ubtav(:,:) = 0.0
ALLOC_(CS%vbtav(isd:ied,JsdB:JedB)) ; CS%vbtav(:,:) = 0.0
if (CS%gradual_BT_ICs) then
Expand Down Expand Up @@ -5287,6 +5339,24 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS)
call register_restart_field(CS%dtbt, "DTBT", .false., restart_CS, &
longname="Barotropic timestep", units="seconds", conversion=US%T_to_s)

! Initialize and register streaming filters
if (CS%use_filter_m2) then
if (am2 > 0.0 .and. om2 > 0.0) then
call Filt_register(am2, om2, 'u', HI, CS%Filt_CS_um2)
call Filt_register(am2, om2, 'v', HI, CS%Filt_CS_vm2)
else
CS%use_filter_m2 = .false.
endif
endif
if (CS%use_filter_k1) then
if (ak1 > 0.0 .and. ok1 > 0.0) then
call Filt_register(ak1, ok1, 'u', HI, CS%Filt_CS_uk1)
call Filt_register(ak1, ok1, 'v', HI, CS%Filt_CS_vk1)
else
CS%use_filter_k1 = .false.
endif
endif

end subroutine register_barotropic_restarts

!> \namespace mom_barotropic
Expand Down
Loading

0 comments on commit 241f7b7

Please sign in to comment.