NECCTONsimulations/additional_simulation.py at main · Parcels-code/NECCTONsimulations · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
""" Simulation script for NECCTON simulations using Parcels and PlasticParcels.
    This script runs a simulation for plastic particles in the North West European Continental Shelf and Arctic Ocean.
    This script should be run on an HPC cluster, using the simulation.sh submission script.
    The submission script takes 3 arguments, which are passed to this script:
        --startrelease: The first date of particle release (format: YYYY-MM-DD)
        --endrelease: The last date of particle release (format: YYYY-MM-DD)
        --endsimulation: The end date of the simulation (format: YYYY-MM-DD)
    Example usage: sbatch simulation.sh --startrelease 2020-01-01 --endrelease 2020-02-01 --endsimulation 2021-01-01
"""

## Library imports
# Data handling
import numpy as np

from argparse import ArgumentParser
import os

# Lagrangian analysis
import parcels
import plasticparcels as pp

from datetime import date, timedelta

# Import custom functions
from functions import create_full_fieldset, add_additional_fields, create_particleset_from_release_files
from kernels import PolyTEOS10_bsq, Biofouling, StokesDrift, checkThroughSurface, stopExecution, reflectAtSurface, reflectAtBathymetry, checkErrorThroughSurface_2DAdvectionRK2, AdvectionRK2_3D, unbeachingBySamplingAfterwards, belowLatitude

# Parse the arguments
p = ArgumentParser(description="""NECCTON Simulation using Parcels""")
p.add_argument('-startrelease', '--startrelease', default='2020-01-01', help='First release date of the particles')
p.add_argument('-endrelease', '--endrelease', default='2020-02-01', help='Last release date of the particles')
p.add_argument('-endsimulation', '--endsimulation', default='2021-01-01', help='Time the simulation ends')
parsed_args = p.parse_args()

# Simulation date settings
args = {'startrelease': parsed_args.startrelease+'T00:00:00.000000000',
        'endrelease': parsed_args.endrelease+'T00:00:00.000000000',
        'endsimulation': parsed_args.endsimulation+'T00:00:00.000000000',
        'outputdt': 1,              # Daily output [d]
        'dt': 20}                   # Timestep choice in minutes

starttime = np.datetime64(args['startrelease'])
endtime = np.datetime64(args['endsimulation'])
runtime = timedelta(days=int((endtime-starttime).astype('timedelta64[D]').astype('int')))

# Create start date of simulation
start_year, start_month, start_day = [int(dateportion) for dateportion in args['startrelease'][:10].split('-')]
starting_time = date(start_year, start_month, start_day)

# Location to store output data, and the location containing release information
output_dir = '/storage/shared/oceanparcels/output_data/data_Michael/NECCTONsimulations/data/additional_simulations/'
release_folder = '/nethome/denes001/Projects/NECCTONsimulations/release_files_additional/'

# Create fieldset
settings = pp.utils.load_settings('NECCTON_settings.json')

# Set model indices for faster loading
settings['ocean']['indices'] = {'lat': range(1800,3059)}
settings['bgc']['indices'] = {'lat': range(590,1022)}
settings['stokes']['indices'] = {'lat': range(230, 361)}

# Create the simulation settings - used to create the fieldset
settings['simulation'] = {
    'startdate': starting_time,      # Start date of simulation
    'endtime': endtime,              # End time of simulation
    'runtime': runtime,              # Runtime of simulation
    'outputdt': timedelta(days=int(args['outputdt'])),       # Timestep of output
    'dt': timedelta(minutes=int(args['dt'])),                # Timestep of advection
    }

# Create the fieldset
fieldset = create_full_fieldset(settings)

# Add constants to the fieldset
fieldset.add_constant('use_mixing', settings['use_mixing'])
fieldset.add_constant('use_biofouling', settings['use_biofouling'])
fieldset.add_constant('use_stokes', settings['use_stokes'])
fieldset.add_constant('use_wind', settings['use_wind'])
fieldset.add_constant('G', 9.81)  # Gravitational constant [m s-1]
fieldset.add_constant('use_3D', settings['use_3D'])

# Add Stokes and biogeochemistry and unbeaching
fieldset = add_additional_fields(fieldset, settings)

# Create the particle set
pset = create_particleset_from_release_files(args['startrelease'], args['endrelease'], release_folder, settings, fieldset)

# Kernels for simulation
kernels = [PolyTEOS10_bsq,
           AdvectionRK2_3D,
           checkErrorThroughSurface_2DAdvectionRK2,
           StokesDrift,
           Biofouling,
           pp.VerticalMixing,
           reflectAtSurface,
           reflectAtBathymetry,
           checkThroughSurface,
           unbeachingBySamplingAfterwards,
           stopExecution,
           belowLatitude]

# Simulation parameters
dt = settings['simulation']['dt']
outputdt = settings['simulation']['outputdt']
endtime = settings['simulation']['endtime']

# Create the particlefile and run the simulation
startdate_s = settings['simulation']['startdate'].isoformat()

pfilename = output_dir+f'particles_{startdate_s}.zarr'

pfile = pp.ParticleFile(pfilename, pset,
                        settings=settings, outputdt=outputdt, chunks=(len(pset), 7)) #7 obs per chunk = 1 week if daily output
pset.execute(kernels, endtime=endtime, dt=dt, output_file=pfile)

# Write the latest locations at the end of the simulation
pfile.write_latest_locations(pset, time=np.max(pset.time_nextloop))

print("Simulation complete.")