Skip to content

Rename u arg to step_MOM_dyn_split_RK2 as u_inst #537

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 2 commits into from
Dec 20, 2023
Merged

Rename u arg to step_MOM_dyn_split_RK2 as u_inst #537

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a477554
Rename u arg to step_MOM_dyn_split_RK2 as u_inst
Hallberg-NOAA Dec 15, 2023
bc192e4
Merge branch 'dev/gfdl' into u_inst_argument
marshallward Dec 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 43 additions & 46 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ module MOM_dynamics_split_RK2
use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix
use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units
use MOM_wave_interface, only: wave_parameters_CS, Stokes_PGF
use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF
use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS
use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS
use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member
Expand Down Expand Up @@ -284,16 +284,16 @@ module MOM_dynamics_split_RK2
contains

!> RK2 splitting for time stepping MOM adiabatic dynamics
subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, &
uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, &
MEKE, thickness_diffuse_CSp, pbv, Waves)
subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, forces, &
p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, pbv, Waves)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
target, intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
target, intent(inout) :: u_inst !< Instantaneous zonal velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
target, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
target, intent(inout) :: v_inst !< Instantaneous meridional velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
Expand Down Expand Up @@ -355,9 +355,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1]

real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_OBC ! The starting zonal velocities, which are
! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]
! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! The starting meridional velocities, which are
! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]
! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1]

! GMM, TODO: make these allocatable?
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix
Expand Down Expand Up @@ -423,8 +423,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)

if (CS%debug) then
call MOM_state_chksum("Start predictor ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
call check_redundant("Start predictor u ", u, v, G, unscale=US%L_T_to_m_s)
call MOM_state_chksum("Start predictor ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
call check_redundant("Start predictor u ", u_inst, v_inst, G, unscale=US%L_T_to_m_s)
call check_redundant("Start predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
endif

Expand Down Expand Up @@ -475,7 +475,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil))

call create_group_pass(CS%pass_uv, u, v, G%Domain, halo=max(2,cont_stencil))
call create_group_pass(CS%pass_uv, u_inst, v_inst, G%Domain, halo=max(2,cont_stencil))
call create_group_pass(CS%pass_h, h, G%Domain, halo=max(2,cont_stencil))
call create_group_pass(CS%pass_av_uvh, u_av, v_av, G%Domain, halo=max(2,obc_stencil))
call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(2,obc_stencil))
Expand Down Expand Up @@ -503,7 +503,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
if (Use_Stokes_PGF) then
call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves)

! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv
! will therefore report the sum total PGF and we avoid other
Expand Down Expand Up @@ -576,15 +576,15 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
!$OMP parallel do default(shared)
do k=1,nz
do j=js,je ; do I=Isq,Ieq
up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * u_bc_accel(I,j,k))
up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * u_bc_accel(I,j,k))
enddo ; enddo
do J=Jsq,Jeq ; do i=is,ie
vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * v_bc_accel(i,J,k))
vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * v_bc_accel(i,J,k))
enddo ; enddo
enddo

call enable_averages(dt, Time_local, CS%diag)
call set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp)
call set_viscous_ML(u_inst, v_inst, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp)
call disable_averaging(CS%diag)

if (CS%debug) then
Expand Down Expand Up @@ -628,7 +628,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! u_accel_bt = layer accelerations due to barotropic solver
if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then
call cpu_clock_begin(id_clock_continuity)
call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
call cpu_clock_end(id_clock_continuity)
if (BT_cont_BT_thick) then
Expand All @@ -639,15 +639,15 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
endif

if (CS%BT_use_layer_fluxes) then
uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u ; v_ptr => v
uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u_inst ; v_ptr => v_inst
endif

call cpu_clock_begin(id_clock_btstep)
if (calc_dtbt) call set_dtbt(G, GV, US, CS%barotropic_CSp, eta, CS%pbce)
if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
! This is the predictor step call to btstep.
! The CS%ADp argument here stores the weights for certain integrated diagnostics.
call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
Expand All @@ -661,11 +661,11 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
!$OMP parallel do default(shared)
do k=1,nz
do J=Jsq,Jeq ; do i=is,ie
vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * &
vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt_pred * &
(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
enddo ; enddo
do j=js,je ; do I=Isq,Ieq
up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * &
up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt_pred * &
(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
enddo ; enddo
enddo
Expand All @@ -679,7 +679,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
call MOM_state_chksum("Predictor 1 init", u, v, h, uh, vh, G, GV, US, haloshift=1, &
call MOM_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1, &
symmetric=sym)
call check_redundant("Predictor 1 up", up, vp, G, unscale=US%L_T_to_m_s)
call check_redundant("Predictor 1 uh", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
Expand Down Expand Up @@ -809,7 +809,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
if (Use_Stokes_PGF) then
call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
if (.not.Waves%Passive_Stokes_PGF) then
do k=1,nz
do j=js,je ; do I=Isq,Ieq
Expand Down Expand Up @@ -899,7 +899,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s

if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
! This is the corrector step call to btstep.
call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
Expand All @@ -920,22 +920,22 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
!$OMP parallel do default(shared)
do k=1,nz
do j=js,je ; do I=Isq,Ieq
u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * &
u_inst(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * &
(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
enddo ; enddo
do J=Jsq,Jeq ; do i=is,ie
v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * &
v_inst(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * &
(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
enddo ; enddo
enddo
call cpu_clock_end(id_clock_mom_update)

if (CS%debug) then
call uvchksum("Corrector 1 [uv]", u, v, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, scale=US%L_T_to_m_s)
call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, scale=GV%H_to_MKS)
call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, &
symmetric=sym, scale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, US, haloshift=1)
! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
symmetric=sym)
Expand All @@ -951,26 +951,26 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
do k = 1, nz
do j = js , je
do I = Isq, Ieq
uold(I,j,k) = u(I,j,k)
uold(I,j,k) = u_inst(I,j,k)
enddo
enddo
do J = Jsq, Jeq
do i = is, ie
vold(i,J,k) = v(i,J,k)
vold(i,J,k) = v_inst(i,J,k)
enddo
enddo
enddo
endif

call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
call vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot,waves=waves)

if (CS%fpmix) then
call vertFPmix(u, v, uold, vold, hbl, h, forces, dt, &
call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, &
G, GV, US, CS%vertvisc_CSp, CS%OBC)
call vertvisc(u, v, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
endif

Expand Down Expand Up @@ -1000,7 +1000,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
! h = h + dt * div . uh
! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
call cpu_clock_begin(id_clock_continuity)
call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
call continuity(u_inst, v_inst, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, u_av, v_av)
call cpu_clock_end(id_clock_continuity)
call do_group_pass(CS%pass_h, G%Domain, clock=id_clock_pass)
Expand All @@ -1016,7 +1016,9 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
endif

if (associated(CS%OBC)) then
call radiation_open_bdry_conds(CS%OBC, u, u_old_rad_OBC, v, v_old_rad_OBC, G, GV, US, dt)
!### I suspect that there is a bug here when u_inst is compared with a previous value of u_av
! to estimate the dominant outward group velocity, but a fix is not available yet.
call radiation_open_bdry_conds(CS%OBC, u_inst, u_old_rad_OBC, v_inst, v_old_rad_OBC, G, GV, US, dt)
endif

! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
Expand Down Expand Up @@ -1156,7 +1158,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
if (CS%id_deta_dt > 0) call post_data(CS%id_deta_dt, deta_dt, CS%diag)

if (CS%debug) then
call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
call MOM_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, scale=GV%H_to_MKS)
! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
Expand Down Expand Up @@ -1471,8 +1473,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av

id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
grain=CLOCK_ROUTINE)
id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=CLOCK_ROUTINE)

call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp)
cont_stencil = continuity_stencil(CS%continuity_CSp)
Expand All @@ -1482,17 +1483,14 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, &
CS%SAL_CSp, CS%tides_CSp)
call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
ntrunc, CS%vertvisc_CSp)
call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, ntrunc, CS%vertvisc_CSp)
CS%set_visc_CSp => set_visc
call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, &
activate=is_new_run(restart_CS) )
call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) )

if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp
if (associated(OBC)) then
CS%OBC => OBC
if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, &
activate=is_new_run(restart_CS) )
if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, activate=is_new_run(restart_CS) )
endif
if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp

Expand All @@ -1515,8 +1513,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo

call barotropic_init(u, v, h, CS%eta, Time, G, GV, US, param_file, diag, &
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
CS%SAL_CSp)
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp)

if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
.not. query_initialized(CS%diffv, "diffv", restart_CS)) then
Expand Down