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 <[email protected]>

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

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

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() !<Lagrangian particles
 end type MOM_control_struct
 
 public initialize_MOM, finish_MOM_initialization, MOM_end
@@ -1097,6 +1101,13 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
 
   endif ! -------------------------------------------------- end SPLIT
 
+   if (CS%do_dynamics) then!run particles whether or not stepping is split
+     if (CS%use_particles) then
+       call particles_run(CS%particles, Time_local, CS%u, CS%v, CS%h, CS%tv) ! Run the particles model
+     endif
+   endif
+
+
   if (CS%thickness_diffuse .and. .not.CS%thickness_diffuse_first) then
     call cpu_clock_begin(id_clock_thick_diff)
 
@@ -2092,6 +2103,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
        "If true, enables the ice shelf model.", default=.false.)
   endif
 
+  call get_param(param_file, "MOM", "USE_PARTICLES", CS%use_particles, &
+                 "If true, use the particles package.", default=.false.)
+
   CS%ensemble_ocean=.false.
   call get_param(param_file, "MOM", "ENSEMBLE_OCEAN", CS%ensemble_ocean, &
                  "If False, The model is being run in serial mode as a single realization. "//&
@@ -2870,6 +2884,11 @@ subroutine finish_MOM_initialization(Time, dirs, CS, restart_CSp)
   call fix_restart_scaling(GV)
   call fix_restart_unit_scaling(US)
 
+
+  if (CS%use_particles) then
+    call particles_init(CS%particles, G, CS%Time, CS%dt_therm, CS%u, CS%v)
+  endif
+
   ! Write initial conditions
   if (CS%write_IC) then
     allocate(restart_CSp_tmp)
@@ -3562,6 +3581,10 @@ end subroutine get_ocean_stocks
 subroutine MOM_end(CS)
   type(MOM_control_struct), intent(inout) :: CS   !< MOM control structure
 
+  if (CS%use_particles) then
+    call particles_save_restart(CS%particles)
+  endif
+
   call MOM_sum_output_end(CS%sum_output_CSp)
 
   if (CS%use_ALE_algorithm) call ALE_end(CS%ALE_CSp)
@@ -3592,6 +3615,11 @@ subroutine MOM_end(CS)
     call end_dyn_unsplit(CS%dyn_unsplit_CSp)
   endif
 
+  if (CS%use_particles) then
+     call particles_end(CS%particles)
+     deallocate(CS%particles)
+  endif
+
   call thickness_diffuse_end(CS%thickness_diffuse_CSp, CS%CDp)
   deallocate(CS%thickness_diffuse_CSp)