From 14f2c97aa9ee2f5d508aa281c8b2fc4a0747bcbd Mon Sep 17 00:00:00 2001 From: Cory Spencer Jones Date: Mon, 23 Dec 2024 15:59:40 -0600 Subject: [PATCH] Update particles_run call to allow accurate particle advection using u, v or uh, vh The drifters package is able to run in 2 different modes: one where the particles are advected using u and v, and one where they are advected using uh and vh. When u and v are used (i.e. the drifters are advected using the resolved velocities only, with no sub-grd component), the particles package needs the layer thickness that is consistent with these velocities, so particles_run needs to be called before any subgrid scale parameterizations are applied. This was not implemented correctly in my last PR, and is corrected here. The particles package also needs the correct timestep for particle advection. This is added here. --- .../external/drifters/MOM_particles.F90 | 4 +-- src/core/MOM.F90 | 26 +++++++++++++------ 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/config_src/external/drifters/MOM_particles.F90 b/config_src/external/drifters/MOM_particles.F90 index 95470e6510..1c41170582 100644 --- a/config_src/external/drifters/MOM_particles.F90 +++ b/config_src/external/drifters/MOM_particles.F90 @@ -28,7 +28,7 @@ subroutine particles_init(parts, Grid, Time, dt, u, v, h) end subroutine particles_init !> The main driver the steps updates particles -subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) +subroutine particles_run(parts, time, uo, vo, ho, tv, dt_adv, use_uh) ! Arguments type(particles), pointer :: parts !< Container for all types and memory type(time_type), intent(in) :: time !< Model time @@ -40,8 +40,8 @@ subroutine particles_run(parts, time, uo, vo, ho, tv, use_uh, stagger) !! that are used to advect tracers [H L2 ~> m3 or kg] real, dimension(:,:,:), intent(in) :: ho !< Ocean layer thickness [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< structure containing pointers to available thermodynamic fields + real, intent(in) :: dt_adv !< timestep for advecting particles [s] logical :: use_uh !< Flag for whether u and v are weighted by thickness - integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered end subroutine particles_run diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index e5d1f46086..8cfab91cfc 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1283,6 +1283,18 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif ! -------------------------------------------------- end SPLIT + if (CS%use_particles .and. CS%do_dynamics .and. (.not. CS%use_uh_particles)) then + if (CS%thickness_diffuse_first) call MOM_error(WARNING,"particles_run: "//& + "Thickness_diffuse_first is true and use_uh_particles is false. "//& + "This is usually a bad combination.") + !Run particles using unweighted velocity + call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & + CS%tv, dt, CS%use_uh_particles) + call particles_to_z_space(CS%particles, h) + endif + + + ! Update the model's current to reflect wind-wave growth if (Waves%Stokes_DDT .and. (.not.Waves%Passive_Stokes_DDT)) then do J=jsq,jeq ; do i=is,ie @@ -1368,23 +1380,21 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & endif call disable_averaging(CS%diag) + ! Advance the dynamics time by dt. + CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt + if (CS%use_particles .and. CS%do_dynamics .and. CS%use_uh_particles) then !Run particles using thickness-weighted velocity call particles_run(CS%particles, Time_local, CS%uhtr, CS%vhtr, CS%h, & - CS%tv, CS%use_uh_particles) - elseif (CS%use_particles .and. CS%do_dynamics) then - !Run particles using unweighted velocity - call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, & - CS%tv, CS%use_uh_particles) + CS%tv, CS%t_dyn_rel_adv, CS%use_uh_particles) endif - - ! Advance the dynamics time by dt. - CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt CS%n_dyn_steps_in_adv = CS%n_dyn_steps_in_adv + 1 if (CS%alternate_first_direction) then call set_first_direction(G, MODULO(G%first_direction+1,2)) CS%first_dir_restart = real(G%first_direction) + elseif (CS%use_particles .and. CS%do_dynamics .and. (.not.CS%use_uh_particles)) then + call particles_to_k_space(CS%particles, h) endif CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt if (abs(CS%t_dyn_rel_thermo) < 1e-6*dt) CS%t_dyn_rel_thermo = 0.0