From 80a1597173b1ff3dc66cc9d500c14712bda35e83 Mon Sep 17 00:00:00 2001 From: Spencer Jones <41342785+cspencerjones@users.noreply.github.com> Date: Mon, 11 Oct 2021 18:12:45 -0500 Subject: [PATCH] Particle API (#1504) * first draft API based on Luyu's code * fixed various errors * Code for particles in MOM.F90 * moved particles_run into dynamics step * added particles_end * Fixed particle time * Fixed some documentation * Further documentation edits * converted pointers to allocatables in particles_gridded * Remove trailing space * Further doxygen tweaks * another trailing space * removed set_time Co-authored-by: Cory Spencer Jones --- .../external/drifters/MOM_particles.F90 | 59 +++++++ .../external/drifters/MOM_particles_types.F90 | 161 ++++++++++++++++++ src/core/MOM.F90 | 28 +++ 3 files changed, 248 insertions(+) create mode 100644 config_src/external/drifters/MOM_particles.F90 create mode 100644 config_src/external/drifters/MOM_particles_types.F90 diff --git a/config_src/external/drifters/MOM_particles.F90 b/config_src/external/drifters/MOM_particles.F90 new file mode 100644 index 0000000000..135f5d284c --- /dev/null +++ b/config_src/external/drifters/MOM_particles.F90 @@ -0,0 +1,59 @@ +!> A set of dummy interfaces for compiling the MOM6 drifters code +module MOM_particles_mod + +use MOM_grid, only : ocean_grid_type +use MOM_time_manager, only : time_type, get_date, operator(-) +use MOM_variables, only : thermo_var_ptrs + + +use particles_types_mod, only: particles, particles_gridded + +public particles_run, particles_init, particles_save_restart, particles_end + +contains + +!> Initializes particles container "parts" +subroutine particles_init(parts, Grid, Time, dt, u, v) + ! Arguments + type(particles), pointer, intent(out) :: parts !< Container for all types and memory + type(ocean_grid_type), target, intent(in) :: Grid !< Grid type from parent model + type(time_type), intent(in) :: Time !< Time type from parent model + real, intent(in) :: dt !< particle timestep in seconds + real, dimension(:,:,:),intent(in) :: u !< Zonal velocity field + real, dimension(:,:,:),intent(in) :: v !< Meridional velocity field + +end subroutine particles_init + +!> The main driver the steps updates particles +subroutine particles_run(parts, time, uo, vo, ho, tv, stagger) + ! Arguments + type(particles), pointer :: parts !< Container for all types and memory + type(time_type), intent(in) :: time !< Model time + real, dimension(:,:,:),intent(in) :: uo !< Ocean zonal velocity (m/s) + real, dimension(:,:,:),intent(in) :: vo !< Ocean meridional velocity (m/s) + 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 + integer, optional, intent(in) :: stagger !< Flag for whether velocities are staggered + +end subroutine particles_run + + +!>Save particle locations (and sometimes other vars) to restart file +subroutine particles_save_restart(parts,temp,salt) +! Arguments +type(particles), pointer :: parts !< Container for all types and memory +real,dimension(:,:,:),optional,intent(in) :: temp !< Optional container for temperature +real,dimension(:,:,:),optional,intent(in) :: salt !< Optional container for salinity + +end subroutine particles_save_restart + +!> Deallocate all memory and disassociated pointer +subroutine particles_end(parts,temp,salt) +! Arguments +type(particles), pointer :: parts !< Container for all types and memory +real,dimension(:,:,:),optional,intent(in) :: temp !< Optional container for temperature +real,dimension(:,:,:),optional,intent(in) :: salt !< Optional container for salinity + +end subroutine particles_end + +end module MOM_particles_mod diff --git a/config_src/external/drifters/MOM_particles_types.F90 b/config_src/external/drifters/MOM_particles_types.F90 new file mode 100644 index 0000000000..b7bc01acb9 --- /dev/null +++ b/config_src/external/drifters/MOM_particles_types.F90 @@ -0,0 +1,161 @@ +!> Dummy data structures and methods for drifters package +module particles_types_mod + +use MOM_grid, only : ocean_grid_type +use mpp_domains_mod, only: domain2D + + +!> Container for gridded fields +type :: particles_gridded + type(domain2D), pointer :: domain !< MPP parallel domain + integer :: halo !< Nominal halo width + integer :: isc !< Start i-index of computational domain + integer :: iec !< End i-index of computational domain + integer :: jsc !< Start j-index of computational domain + integer :: jec !< End j-index of computational domain + integer :: isd !< Start i-index of data domain + integer :: ied !< End i-index of data domain + integer :: jsd !< Start j-index of data domain + integer :: jed !< End j-index of data domain + integer :: isg !< Start i-index of global domain + integer :: ieg !< End i-index of global domain + integer :: jsg !< Start j-index of global domain + integer :: jeg !< End j-index of global domain + integer :: is_offset=0 !< add to i to recover global i-index + integer :: js_offset=0 !< add to j to recover global j-index + integer :: my_pe !< MPI PE index + integer :: pe_N !< MPI PE index of PE to the north + integer :: pe_S !< MPI PE index of PE to the south + integer :: pe_E !< MPI PE index of PE to the east + integer :: pe_W !< MPI PE index of PE to the west + logical :: grid_is_latlon !< Flag to say whether the coordinate is in lat-lon degrees, or meters + logical :: grid_is_regular !< Flag to say whether point in cell can be found assuming regular Cartesian grid + real :: Lx !< Length of the domain in x direction + real, dimension(:,:), allocatable :: lon !< Longitude of cell corners (degree E) + real, dimension(:,:), allocatable :: lat !< Latitude of cell corners (degree N) + real, dimension(:,:), allocatable :: lonc !< Longitude of cell centers (degree E) + real, dimension(:,:), allocatable :: latc !< Latitude of cell centers (degree N) + real, dimension(:,:), allocatable :: dx !< Length of cell edge (m) + real, dimension(:,:), allocatable :: dy !< Length of cell edge (m) + real, dimension(:,:), allocatable :: area !< Area of cell (m^2) + real, dimension(:,:), allocatable :: msk !< Ocean-land mask (1=ocean) + real, dimension(:,:), allocatable :: cos !< Cosine from rotation matrix to lat-lon coords + real, dimension(:,:), allocatable :: sin !< Sine from rotation matrix to lat-lon coords + real, dimension(:,:), allocatable :: ocean_depth !< Depth of ocean (m) + real, dimension(:,:), allocatable :: uo !< Ocean zonal flow (m/s) + real, dimension(:,:), allocatable :: vo !< Ocean meridional flow (m/s) + real, dimension(:,:), allocatable :: tmp !< Temporary work space + real, dimension(:,:), allocatable :: tmpc !< Temporary work space + real, dimension(:,:), allocatable :: parity_x !< X component of vector point from i,j to i+1,j+1 + real, dimension(:,:), allocatable :: parity_y !< Y component of vector point from i,j to i+1,j+1 + ! (For detecting tri-polar fold) + integer, dimension(:,:), allocatable :: particle_counter_grd !< Counts particles created for naming purposes + !>@{ + !! Diagnostic handle + integer :: id_uo=-1, id_vo=-1, id_unused=-1 + integer :: id_count=-1, id_chksum=-1 + !>@} + +end type particles_gridded + + +!>xyt is a data structure containing particle position and velocity fields. +type :: xyt + real :: lon !< Longitude of particle (degree N or unit of grid coordinate) + real :: lat !< Latitude of particle (degree N or unit of grid coordinate) + real :: day !< Day of this record (days) + real :: lat_old !< Previous latitude + real :: lon_old !< Previous longitude + real :: uvel !< Zonal velocity of particle (m/s) + real :: vvel !< Meridional velocity of particle (m/s) + real :: uvel_old !< Previous zonal velocity component (m/s) + real :: vvel_old !< Previous meridional velocity component (m/s) + integer :: year !< Year of this record + integer :: particle_num !< Current particle number + integer(kind=8) :: id = -1 !< Particle Identifier + type(xyt), pointer :: next=>null() !< Pointer to the next position in the list +end type xyt + +!>particle types are data structures describing a tracked particle +type :: particle + type(particle), pointer :: prev=>null() !< Previous link in list + type(particle), pointer :: next=>null() !< Next link in list +! State variables (specific to the particles, needed for restarts) + real :: lon !< Longitude of particle (degree N or unit of grid coordinate) + real :: lat !< Latitude of particle (degree E or unit of grid coordinate) + real :: depth !< Depth of particle + real :: uvel !< Zonal velocity of particle (m/s) + real :: vvel !< Meridional velocity of particle (m/s) + real :: lon_old !< previous lon (degrees) + real :: lat_old !< previous lat (degrees) + real :: uvel_old !< previous uvel + real :: vvel_old !< previous vvel + real :: start_lon !< starting longitude where particle was created + real :: start_lat !< starting latitude where particle was created + real :: start_day !< origination position (degrees) and day + integer :: start_year !< origination year + real :: halo_part !< equal to zero for particles on the computational domain, and 1 for particles on the halo + integer(kind=8) :: id !< particle identifier + integer(kind=8) :: drifter_num !< particle identifier + integer :: ine !< nearest i-index in NE direction (for convenience) + integer :: jne !< nearest j-index in NE direction (for convenience) + real :: xi !< non-dimensional x-coordinate within current cell (0..1) + real :: yj !< non-dimensional y-coordinate within current cell (0..1) + real :: uo !< zonal ocean velocity + real :: vo !< meridional ocean velocity + !< by the particle (m/s) + type(xyt), pointer :: trajectory=>null() !< Trajectory for this particle +end type particle + + +!>A buffer structure for message passing +type :: buffer + integer :: size=0 !< Size of buffer + real, dimension(:,:), pointer :: data !< Buffer memory +end type buffer + +!> A wrapper for the particle linked list (since an array of pointers is not allowed) +type :: linked_list + type(particle), pointer :: first=>null() !< Pointer to the beginning of a linked list of parts +end type linked_list + + +!> A grand data structure for the particles in the local MOM domain +type :: particles !; private + type(particles_gridded) :: grd !< Container with all gridded data + type(linked_list), dimension(:,:), allocatable :: list !< Linked list of particles + type(xyt), pointer :: trajectories=>null() !< A linked list for detached segments of trajectories + real :: dt !< Time-step between particle calls + integer :: current_year !< Current year (years) + real :: current_yearday !< Current year-day, 1.00-365.99, (days) + integer :: traj_sample_hrs !< Period between sampling for trajectories (hours) + integer :: traj_write_hrs !< Period between writing of trajectories (hours) + integer :: verbose_hrs !< Period between terminal status reports (hours) + !>@{ + !! Handles for clocks + integer :: clock, clock_mom, clock_the, clock_int, clock_cal, clock_com, clock_ini, clock_ior, clock_iow, clock_dia + integer :: clock_trw, clock_trp + !>@} + logical :: restarted=.false. !< Indicate whether we read state from a restart or not + logical :: Runge_not_Verlet=.True. !< True=Runge-Kutta, False=Verlet. + logical :: ignore_missing_restart_parts=.False. !< True allows the model to ignore particles missing in the restart. + logical :: halo_debugging=.False. !< Use for debugging halos (remove when its working) + logical :: save_short_traj=.false. !< True saves only lon,lat,time,id in particle_trajectory.nc + logical :: ignore_traj=.False. !< If true, then model does not write trajectory data at all + logical :: use_new_predictive_corrective =.False. !< Flag to use Bob's predictive corrective particle scheme + !Added by Alon + integer(kind=8) :: debug_particle_with_id = -1 !< If positive, monitors a part with this id + type(buffer), pointer :: obuffer_n=>null() !< Buffer for outgoing parts to the north + type(buffer), pointer :: ibuffer_n=>null() !< Buffer for incoming parts from the north + type(buffer), pointer :: obuffer_s=>null() !< Buffer for outgoing parts to the south + type(buffer), pointer :: ibuffer_s=>null() !< Buffer for incoming parts from the south + type(buffer), pointer :: obuffer_e=>null() !< Buffer for outgoing parts to the east + type(buffer), pointer :: ibuffer_e=>null() !< Buffer for incoming parts from the east + type(buffer), pointer :: obuffer_w=>null() !< Buffer for outgoing parts to the west + type(buffer), pointer :: ibuffer_w=>null() !< Buffer for incoming parts from the west + type(buffer), pointer :: obuffer_io=>null() !< Buffer for outgoing parts during i/o + type(buffer), pointer :: ibuffer_io=>null() !< Buffer for incoming parts during i/o +end type particles + + +end module particles_types_mod diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 1f1492f2e5..73a7cd58a7 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -150,6 +150,7 @@ module MOM use MOM_offline_main, only : offline_advection_layer, offline_transport_end use MOM_ALE, only : ale_offline_tracer_final, ALE_main_offline use MOM_ice_shelf, only : ice_shelf_CS, ice_shelf_query, initialize_ice_shelf +use MOM_particles_mod, only : particles, particles_init, particles_run, particles_save_restart, particles_end implicit none ; private @@ -329,6 +330,8 @@ module MOM logical :: answers_2018 !< If true, use expressions for the surface properties that recover !! the answers from the end of 2018. Otherwise, use more appropriate !! expressions that differ at roundoff for non-Boussinsq cases. + logical :: use_particles !< Turns on the particles package + character(len=10) :: particle_type !< Particle types include: surface(default), profiling and sail drone. type(MOM_diag_IDs) :: IDs !< Handles used for diagnostics. type(transport_diag_IDs) :: transport_IDs !< Handles used for transport diagnostics. @@ -396,6 +399,7 @@ module MOM type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling !! ensemble model state vectors and data assimilation !! increments and priors + type(particles), pointer :: particles => NULL() !