-
Notifications
You must be signed in to change notification settings - Fork 7
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #122 from metno/bomb_decay
Bomb decay
- Loading branch information
Showing
22 changed files
with
869 additions
and
341 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,59 @@ | ||
| #! /usr/bin/env python3 | ||
|
|
||
| ''' | ||
| Convert aerosols in a SNAP file to isotopes using the fractional distribution from Tovdal (2002) | ||
| ''' | ||
|
|
||
| import netCDF4 | ||
| from Snappy.BombIsotopeFractions import BombIsotopeFractions | ||
|
|
||
|
|
||
| def snap_add_bomb_isotopes(nc: netCDF4.Dataset): | ||
| ''' | ||
| ncfile: a netcdf-file with Aerosols opened in 'a'-mode | ||
| ''' | ||
| bomb_isotopes = BombIsotopeFractions() | ||
| aerosols = [] | ||
| for var in nc.variables: | ||
| if var.startswith('Aerosol') and var.endswith('acc_concentration'): | ||
| aerosols.append(var[:-18]) | ||
| isos = bomb_isotopes.isotopes() | ||
| hours = nc['time'][:] # snap writes usually hours since start | ||
| for var in ['concentration', 'acc_dry_deposition', 'acc_wet_deposition', 'acc_concentration']: | ||
| basevar = nc[f"{aerosols[0]}_{var}"] | ||
| for iso in isos: | ||
| # declare variables | ||
| name = f"{iso}_{var}" | ||
| if name not in nc.variables: | ||
| nc.createVariable(name, basevar.datatype, basevar.dimensions, zlib=True, | ||
| chunksizes=basevar.chunking()) | ||
| for attr in basevar.ncattrs(): | ||
| nc[name].setncattr(attr, basevar.getncattr(attr)) | ||
| laststepdata = 0 | ||
| for t, hr in enumerate(hours): | ||
| # convert data | ||
| basedata = 0 | ||
| for aero in aerosols: | ||
| basedata += nc[f"{aero}_{var}"][t,:] | ||
| for iso in isos: | ||
| name = f"{iso}_{var}" | ||
| frac = bomb_isotopes.fraction(iso, hr) | ||
| if (var == 'acc_concentration') and t > 1: | ||
| # no decay in dose-equivalent | ||
| nc[name][t,:] = nc[name][t-1,:] + (basedata-laststepdata)*frac | ||
| else: | ||
| nc[name][t,:] = frac*basedata | ||
| laststepdata = basedata | ||
|
|
||
|
|
||
| def main(): | ||
| import argparse | ||
| parser = argparse.ArgumentParser(description="add isotope distribution to a snap.nc with bomb-aerosols") | ||
| parser.add_argument("--nc", help="snap.nc filename", required=True) | ||
|
|
||
| args = parser.parse_args() | ||
| with netCDF4.Dataset(args.nc, 'a') as nc: | ||
| snap_add_bomb_isotopes(args.nc) | ||
|
|
||
| if __name__ == "__main__": | ||
| main() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,76 @@ | ||
| import os | ||
| import csv | ||
| import re | ||
|
|
||
| class BombIsotopeFractions: | ||
| # fractions{isotope} = {0: frac0, 10: frac10, 20: frac20} with 0-40 in hours and frac0-frac40 fractions of total fission-products in Bq | ||
| _fractions = None | ||
| _timesteps = [0,10,20,30,40] | ||
| def __new__(cls): | ||
| '''BombIsotopeFractions singleton data''' | ||
| if cls._fractions is None: | ||
| cls._fractions = dict() | ||
| directory = os.path.join(os.path.dirname(__file__), "resources") | ||
| with open( | ||
| os.path.join(directory, "bomb-isotope-distribution_Tovedal.csv"), mode="r", encoding="UTF-8", newline='' | ||
| ) as fh: | ||
| csvreader = csv.reader(fh, delimiter=',') | ||
| for i in range(2): | ||
| next(csvreader) | ||
| header = next(csvreader) | ||
| offset = 9 | ||
| for i,hrs in enumerate(cls._timesteps): | ||
| if f"t={hrs}" not in header[offset+i]: | ||
| raise Exception(f"error in header for hour {hrs}: {header[offset+i]}") | ||
| for row in csvreader: | ||
| if '-' in row[0]: | ||
| isotope = row[0].replace('-', '') # without - as ususal in snappy | ||
| stepfraction = {} | ||
| for i,hrs in enumerate(cls._timesteps): | ||
| stepfraction[hrs] = float(row[offset+i])/100. | ||
| cls._fractions[isotope] = stepfraction | ||
| obj = object.__new__(cls) | ||
| return obj | ||
|
|
||
| def isotopes(self): | ||
| ''' | ||
| list over isotopes as ['Cs137', 'Cs134', ...] | ||
| ''' | ||
| return self._fractions.keys() | ||
|
|
||
| def fraction(self, isotope: str, hrs: int) -> float: | ||
| ''' | ||
| @param isotope is a isotope name like Cs137 or Cs-137 | ||
| @param hrs since bomb, intra/extrapolated | ||
| return a fraction of the total activity | ||
| ''' | ||
| isotope = isotope.replace("-", "") | ||
| stepfracs = self._fractions[isotope] | ||
| if hrs < 0: | ||
| return 0 | ||
| if hrs > self._timesteps[-1]: | ||
| return stepfracs[self._timesteps[-1]] | ||
| if hrs == self._timesteps[0]: | ||
| return stepfracs[self._timesteps[0]] | ||
|
|
||
| for i, nhr in enumerate(self._timesteps): | ||
| if nhr >= hrs: | ||
| phr = self._timesteps[i-1] | ||
| hfrac = (hrs-phr)/(nhr-phr) | ||
| nfrac = stepfracs[nhr] | ||
| pfrac = stepfracs[phr] | ||
| frac = pfrac + hfrac*(nfrac-pfrac) | ||
| return frac | ||
|
|
||
|
|
||
| if __name__ == "__main__": | ||
| bfracs = BombIsotopeFractions() | ||
| assert('Cs137' in bfracs.isotopes()) | ||
| assert(bfracs.fraction('Cs137',0) == 0.0002/100) | ||
| assert(len(bfracs.isotopes()) > 10) | ||
| for hr in range(0,48): | ||
| tot = 0 | ||
| for iso in bfracs.isotopes(): | ||
| tot += bfracs.fraction(iso, hr) | ||
| assert(tot > .99) | ||
| assert(tot < 1.01) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.