The code can be found in the PML GitLab. Adjust the makefile to point to your netCDF library paths and compile the code with make. The output file is called ptraj.

Note: the code is not parallelised, so runs slowly and on a single core only.

The model configuration is controlled by the file casename_run.dat. An example is given below.

!=========================================================================================!
! INPUT FILE FOR PARAMETERS CONTROLLING EXECUTION OF OFF-LINE PARTICLE TRACKING           !
! DESCRIPTION OF VARIABLES AND SUGGESTED PARAMETERS CAN BE FOUND AT BOTTOM                !
!                                                                                         !
!        FORMAT:                                                                          !
!       1.) VARIABLE  = VALUE  (EQUAL SIGN MUST BE USED)                                  !
!       2.) FLOATING POINT VARIABLES MUST CONTAIN A PERIOD "." EX: 1.3, 2.,etc            !
!       3.) BLANK LINES ARE IGNORED AS ARE LINES BEGINNING WITH ! (F90 COMMENT)           !
!       4.) COMMENTS CAN FOLLOW VALUES IF MARKED BY !                                     !
!       5.) ORDER OF VARIABLES IS NOT IMPORTANT                                           !
!       6.) FOR MULTIPLE VALUE VARIABLES FIRST ENTRY IS NUMBER OF VARIABLES               !
!           TO FOLLOW (OR 0 IF NONE)                                                      !
!       7.) DO NOT USE COMMAS TO SEPARATE VARIABLES                                       !
!       8.) DO NOT EXCEED EIGHTY CHARACTERS PER LINE                                      !
!       9.) FOR LINE CONTINUATION ADD \\ TO END OF LINE TO FORCE CONTINUE                 !
!           TO NEXT LINE.  MAXIMUM 4 CONTINUATIONS                                        !
!       10.) TRUE = T, FALSE = F                                                          !
!                                                                                         !    
!=========================================================================================!

!=========Parameters Controlling Time Integration==========================================

DTI      = 120.0 ! Time step for scheme resolution. mod(INSTP/DTI) should be 0
INSTP    = 3600  ! INPUT time step (s).
DTOUT    = 1     ! output time step in hour
TDRIFT   = 720 !480 !720 ! Drift duration (hour)

!=========Parameters Controlling Starting Tracking Date====================================

YEARLAG   = 1995 ! Year of tracking start
MONTHLAG  = 1    ! Month of tracking start
DAYLAG    = 1    ! Day of tracking start
HOURLAG   = 0    ! Hour of tracking start

!=========Parameters Controlling Input/Output LOCATION=====================================

INPDIR   = INPDIR
GEOAREA  = tst
OUTDIR   = OUTDIR_tst
INFOFILE = screen
LAGINI   = particle_tst   ! Input file with particle location
METRICS  = tst            ! Prefix for the grid metrics netCDF file (saved in OUTDIR)

!=========Parameters Controlling SIGMA/CARTESIAN===========================================

F_DEPTH  = F              ! Cartesian depth is kept constant if T.
P_SIGMA  = F              ! Input particle depth in sigma if T.
OUT_SIGMA  = F            ! Output particle depth in sigma if T.

!=========Parameters Controlling Random Walk===============================================

IRW   = 0   ! IRW (0-w/o rw; 1-hor rw; 2-vert rw; 3:hor+vert rw)
DHOR  = 20. ! DHOR (horizontal diff coeff m^2/s, miller 10^6m^2/day~11.57m^2/s)
DTRW  = 6.  ! time step (s) for RW (from the visser's criterion). mod(DTI/DTRW) should be 0

Parameters of note are YEARLAG, MONTHLAG, DAYLAG and HOURLAG which determine when the particle tracking should begin. In spite of YEARLAG being defined in this control file, the source code currently has hardcoded values for the year, so you will need to edit offlag.f90 around lines 45 and 46 to use the correct year for your model run.

LAGINI is the file containing the initial particle positions.

GEOAREA is the equivalent of casename in the FVCOM manual. Your .dat file must be called GEOAREA_run.dat (i.e. casename_run.dat).

Set the INPDIR and OUTDIR values to reflect the input and output directories you have from your FVCOM model run. Output files must be in a directory for each year. The naming of the model output files must be casename_####.nc where #### is a number starting from 0001.

    INPDIR/YEARLAG/casename_0001.nc
    INPDIR/YEARLAG/casename_0002.nc
    ...
    INPDIR/YEARLAG/casename_000n.nc

A script to chop up monthly model output files into day long individual files can be found here.

INFOFILE is either screen to output run time information to standard output (i.e. to the terminal) or a file name.

If you are using the medusa-v3.1.0-LAG branch of the offline Lagrangian code, you must add an extra parameter to the Parameters Controlling Input/Output LOCATION section:

    !=========Parameters Controlling Input/Output LOCATION=====================================
    
    INPDIR   = INPDIR
    GEOAREA  = tst
    OUTDIR   = OUTDIR_tst
    INFOFILE = screen
    LAGINI   = particle_tst   ! Input file with particle location
    METRICS  = tst            ! Prefix for the grid metrics netCDF file (saved in OUTDIR)

The METRICS parameter controls the file name for the grid metrics. This netCDF file contains the grid metrics which are necesssary for advecting the particles. This can be a time consuming operation if the grid is complicated, so saving those outputs for successive runs can be a significant time saving.

The way it works is if the file does not exist in the input directory, the grid metrics are calculated and saved to a netCDF file named tst_metrics.nc in the example above (change METRICS to change the prefix). Once the file has been generated on first run, successive runs can use it if it is moved into the input directory (INPDIR/YEARLAG).

The offline code requires model output which contains the following variables:

|Variable name|Long name|Dimensions|Units| |x |nodal x-coordinate |node |metres| |y |nodal y-coordinate |node |metres| |nv |nodes surrounding element |three, nele |--| |h |Bathymetry |node |metres| |siglev |Sigma Levels |siglev, node |--| |a1u |a1u |three, nele |--| |a2u |a2u |three, nele |--| |aw0 |aw0 |three, nele |--| |awx |awx |three, nele |--| |awy |awy |three, nele |--| |time |time |time |days since 1858-11-17 00:00:00| |zeta |Water Surface Elevation |time, node |metres| |u |Eastward Water Velocity |time, siglay, nele |ms^-1| |v |Northward Water Velocity |time, siglay, nele |ms^-1| |omega |Vertical Sigma Coordinate Velocity |time, siglev, node |s^-1| |km |Turbulent Eddy Viscosity For Momentum |time, siglev, node |m² s^-1|

Note: there is no support for spherical coordinates, so you will have to either run your original model in cartesian, or enable FLAG_6 = -DPROJ when compiling FVCOM (and subsequently set the PROJECTION_REFERENCE in your casename_run.nml file).

Note: your FVCOM model output must be in day-long chunks. The particle tracking will only read the first 24 hours of any results file it's given. A script to chop up monthly model output files into day long individual files can be found here.

The offline code requires only a single file to define the release sites for the particles. An example is given below:

3
1 510876 6555485 5.0
2 500064 6563972 0.0
3 530350 6561007 10.0

The first line defines the number of particles to track (this is read in with a DO loop in offlag.f90, so if you have a long list of positions but only wish to use the first few, set the header number to the number you wish to read in, leaving the remaining points defined in the list; the code will ignore them).

The definition of each point is ID XPOS YPOS DEPTH with depth being positive down. Set the name of this file to match the LAGINI value in the casename_run.dat file.