|
| 1 | +"""CTD_BGC instrument.""" |
| 2 | + |
| 3 | +from dataclasses import dataclass |
| 4 | +from datetime import timedelta |
| 5 | +from pathlib import Path |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +from parcels import FieldSet, JITParticle, ParticleSet, Variable |
| 9 | + |
| 10 | +from ..spacetime import Spacetime |
| 11 | + |
| 12 | + |
| 13 | +@dataclass |
| 14 | +class CTD_BGC: |
| 15 | + """Configuration for a single BGC CTD.""" |
| 16 | + |
| 17 | + spacetime: Spacetime |
| 18 | + min_depth: float |
| 19 | + max_depth: float |
| 20 | + |
| 21 | + |
| 22 | +_CTD_BGCParticle = JITParticle.add_variables( |
| 23 | + [ |
| 24 | + Variable("o2", dtype=np.float32, initial=np.nan), |
| 25 | + Variable("chl", dtype=np.float32, initial=np.nan), |
| 26 | + Variable("raising", dtype=np.int8, initial=0.0), # bool. 0 is False, 1 is True. |
| 27 | + Variable("max_depth", dtype=np.float32), |
| 28 | + Variable("min_depth", dtype=np.float32), |
| 29 | + Variable("winch_speed", dtype=np.float32), |
| 30 | + ] |
| 31 | +) |
| 32 | + |
| 33 | + |
| 34 | +def _sample_o2(particle, fieldset, time): |
| 35 | + particle.o2 = fieldset.o2[time, particle.depth, particle.lat, particle.lon] |
| 36 | + |
| 37 | + |
| 38 | +def _sample_chlorophyll(particle, fieldset, time): |
| 39 | + particle.chl = fieldset.chl[time, particle.depth, particle.lat, particle.lon] |
| 40 | + |
| 41 | + |
| 42 | +def _ctd_bgc_cast(particle, fieldset, time): |
| 43 | + # lowering |
| 44 | + if particle.raising == 0: |
| 45 | + particle_ddepth = -particle.winch_speed * particle.dt |
| 46 | + if particle.depth + particle_ddepth < particle.max_depth: |
| 47 | + particle.raising = 1 |
| 48 | + particle_ddepth = -particle_ddepth |
| 49 | + # raising |
| 50 | + else: |
| 51 | + particle_ddepth = particle.winch_speed * particle.dt |
| 52 | + if particle.depth + particle_ddepth > particle.min_depth: |
| 53 | + particle.delete() |
| 54 | + |
| 55 | + |
| 56 | +def simulate_ctd_bgc( |
| 57 | + fieldset: FieldSet, |
| 58 | + out_path: str | Path, |
| 59 | + ctd_bgcs: list[CTD_BGC], |
| 60 | + outputdt: timedelta, |
| 61 | +) -> None: |
| 62 | + """ |
| 63 | + Use Parcels to simulate a set of BGC CTDs in a fieldset. |
| 64 | +
|
| 65 | + :param fieldset: The fieldset to simulate the BGC CTDs in. |
| 66 | + :param out_path: The path to write the results to. |
| 67 | + :param ctds: A list of BGC CTDs to simulate. |
| 68 | + :param outputdt: Interval which dictates the update frequency of file output during simulation |
| 69 | + :raises ValueError: Whenever provided BGC CTDs, fieldset, are not compatible with this function. |
| 70 | + """ |
| 71 | + WINCH_SPEED = 1.0 # sink and rise speed in m/s |
| 72 | + DT = 10.0 # dt of CTD simulation integrator |
| 73 | + |
| 74 | + if len(ctd_bgcs) == 0: |
| 75 | + print( |
| 76 | + "No BGC CTDs provided. Parcels currently crashes when providing an empty particle set, so no BGC CTD simulation will be done and no files will be created." |
| 77 | + ) |
| 78 | + # TODO when Parcels supports it this check can be removed. |
| 79 | + return |
| 80 | + |
| 81 | + fieldset_starttime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[0]) |
| 82 | + fieldset_endtime = fieldset.time_origin.fulltime(fieldset.U.grid.time_full[-1]) |
| 83 | + |
| 84 | + # deploy time for all ctds should be later than fieldset start time |
| 85 | + if not all( |
| 86 | + [ |
| 87 | + np.datetime64(ctd_bgc.spacetime.time) >= fieldset_starttime |
| 88 | + for ctd_bgc in ctd_bgcs |
| 89 | + ] |
| 90 | + ): |
| 91 | + raise ValueError("BGC CTD deployed before fieldset starts.") |
| 92 | + |
| 93 | + # depth the bgc ctd will go to. shallowest between bgc ctd max depth and bathymetry. |
| 94 | + max_depths = [ |
| 95 | + max( |
| 96 | + ctd_bgc.max_depth, |
| 97 | + fieldset.bathymetry.eval( |
| 98 | + z=0, |
| 99 | + y=ctd_bgc.spacetime.location.lat, |
| 100 | + x=ctd_bgc.spacetime.location.lon, |
| 101 | + time=0, |
| 102 | + ), |
| 103 | + ) |
| 104 | + for ctd_bgc in ctd_bgcs |
| 105 | + ] |
| 106 | + |
| 107 | + # CTD depth can not be too shallow, because kernel would break. |
| 108 | + # This shallow is not useful anyway, no need to support. |
| 109 | + if not all([max_depth <= -DT * WINCH_SPEED for max_depth in max_depths]): |
| 110 | + raise ValueError( |
| 111 | + f"BGC CTD max_depth or bathymetry shallower than maximum {-DT * WINCH_SPEED}" |
| 112 | + ) |
| 113 | + |
| 114 | + # define parcel particles |
| 115 | + ctd_bgc_particleset = ParticleSet( |
| 116 | + fieldset=fieldset, |
| 117 | + pclass=_CTD_BGCParticle, |
| 118 | + lon=[ctd_bgc.spacetime.location.lon for ctd_bgc in ctd_bgcs], |
| 119 | + lat=[ctd_bgc.spacetime.location.lat for ctd_bgc in ctd_bgcs], |
| 120 | + depth=[ctd_bgc.min_depth for ctd_bgc in ctd_bgcs], |
| 121 | + time=[ctd_bgc.spacetime.time for ctd_bgc in ctd_bgcs], |
| 122 | + max_depth=max_depths, |
| 123 | + min_depth=[ctd_bgc.min_depth for ctd_bgc in ctd_bgcs], |
| 124 | + winch_speed=[WINCH_SPEED for _ in ctd_bgcs], |
| 125 | + ) |
| 126 | + |
| 127 | + # define output file for the simulation |
| 128 | + out_file = ctd_bgc_particleset.ParticleFile(name=out_path, outputdt=outputdt) |
| 129 | + |
| 130 | + # execute simulation |
| 131 | + ctd_bgc_particleset.execute( |
| 132 | + [_sample_o2, _sample_chlorophyll, _ctd_bgc_cast], |
| 133 | + endtime=fieldset_endtime, |
| 134 | + dt=DT, |
| 135 | + verbose_progress=False, |
| 136 | + output_file=out_file, |
| 137 | + ) |
| 138 | + |
| 139 | + # there should be no particles left, as they delete themselves when they resurface |
| 140 | + if len(ctd_bgc_particleset.particledata) != 0: |
| 141 | + raise ValueError( |
| 142 | + "Simulation ended before BGC CTD resurfaced. This most likely means the field time dimension did not match the simulation time span." |
| 143 | + ) |
0 commit comments