Skip to content

Commit

Permalink
Merge pull request #152 from metno/140-snappy-window-ignores-minutes-…
Browse files Browse the repository at this point in the history 
…in-release-start-time

snappy window ignores minutes in release start time
  • Loading branch information
@heikoklein
heikoklein authored Jan 28, 2025
2 parents d0c0374 + d83ce92 commit 02c0476
Show file tree
Hide file tree
Showing 5 changed files with 71 additions and 68 deletions.
28 changes: 16 additions & 12 deletions src/common/decay.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,15 +48,17 @@ end subroutine decay
!> param: tstep model timestep in seconds
!> NEEDS TO BE RUN BEFORE ::decay
subroutine decayDeps(tstep)
USE iso_fortran_env, only: int16
USE snapfldML, only: depdry, depwet, accdry, accwet, &
total_activity_released, total_activity_lost_domain, total_activity_lost_other
USE snapparML, only: ncomp, run_comp, def_comp
USE snapparML, only: ncomp, nocomp, run_comp, def_comp, output_component
USE iso_fortran_env, only: real64
USE releaseML, only: tpos_bomb

real, intent(in) :: tstep

integer :: m, mm
integer(kind=int16), pointer :: mmo(:)
real :: bomb_decay_rate, current_state, next_state
logical, save :: prepare = .TRUE.
logical, save :: has_bomb_decay = .FALSE.
Expand Down Expand Up @@ -103,17 +105,19 @@ subroutine decayDeps(tstep)
total_steps = total_steps + tstep
end if

do m=1,ncomp
mm = run_comp(m)%to_defined
if(def_comp(mm)%kdecay >= 1) then
depdry(:,:,m) = depdry(:,:,m)*def_comp(mm)%decayrate
depwet(:,:,m) = depwet(:,:,m)*def_comp(mm)%decayrate
accdry(:,:,m) = accdry(:,:,m)*def_comp(mm)%decayrate
accwet(:,:,m) = accwet(:,:,m)*def_comp(mm)%decayrate

total_activity_released(m) = total_activity_released(m)*def_comp(mm)%decayrate
total_activity_lost_domain(m) = total_activity_lost_domain(m)*def_comp(mm)%decayrate
total_activity_lost_other(m) = total_activity_lost_other(m)*def_comp(mm)%decayrate
do m=1,nocomp
mmo => output_component(m)%to_defined
! decay-rates of merged components must be similar, otherwise, decay of output-fields
! does not work well enough
if(def_comp(mmo(1))%kdecay >= 1) then
depdry(:,:,m) = depdry(:,:,m)*def_comp(mmo(1))%decayrate
depwet(:,:,m) = depwet(:,:,m)*def_comp(mmo(1))%decayrate
accdry(:,:,m) = accdry(:,:,m)*def_comp(mmo(1))%decayrate
accwet(:,:,m) = accwet(:,:,m)*def_comp(mmo(1))%decayrate

total_activity_released(m) = total_activity_released(m)*def_comp(mmo(1))%decayrate
total_activity_lost_domain(m) = total_activity_lost_domain(m)*def_comp(mmo(1))%decayrate
total_activity_lost_other(m) = total_activity_lost_other(m)*def_comp(mmo(1))%decayrate
endif
enddo
return
Expand Down
4 changes: 2 additions & 2 deletions src/common/drydep.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ subroutine drydep1(part)
dep = part%scale_rad(1.0 - def_comp(m)%drydeprat)
i = hres_pos(part%x)
j = hres_pos(part%y)
mm = def_comp(m)%to_running
mm = def_comp(m)%to_output
!$OMP atomic
depdry(i,j,mm) = depdry(i,j,mm) + dble(dep)
end if
Expand Down Expand Up @@ -98,7 +98,7 @@ subroutine drydep2(tstep, part)
dep = part%scale_rad(1.0 - deprate)
i = hres_pos(part%x)
j = hres_pos(part%y)
mm = def_comp(m)%to_running
mm = def_comp(m)%to_output
!$OMP atomic
depdry(i,j,mm) = depdry(i,j,mm) + dble(dep)
end if
Expand Down
6 changes: 3 additions & 3 deletions src/common/wetdep.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ subroutine wetdep2(depwet, tstep, part, pextra)
!> uses the precipitation at the particle position
type(extraParticle), intent(in) :: pextra

integer :: m, i, j, mm
integer :: m, i, j, mm, mo
real :: precint, deprate, dep, q

m = part%icomp
Expand All @@ -63,10 +63,10 @@ subroutine wetdep2(depwet, tstep, part, pextra)

i = hres_pos(part%x)
j = hres_pos(part%y)
mm = def_comp(m)%to_running
mo = def_comp(m)%to_output

!$OMP atomic
depwet(i, j, mm) = depwet(i, j, mm) + dble(dep)
depwet(i, j, mo) = depwet(i, j, mo) + dble(dep)
end if
end subroutine wetdep2

Expand Down
70 changes: 44 additions & 26 deletions utils/SnapPy/Snappy/SnapController.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,10 @@
from Snappy.BrowserWidget import BrowserWidget
from Snappy.EcMeteorologyCalculator import EcMeteorologyCalculator
from Snappy.ICONMeteorologyCalculator import ICONMeteorologyCalculator
from Snappy.MeteorologyCalculator import MeteoDataNotAvailableException, MeteorologyCalculator
from Snappy.MeteorologyCalculator import (
MeteoDataNotAvailableException,
MeteorologyCalculator,
)
from Snappy.MailImages import sendPngsFromDir
from Snappy.Resources import Resources, MetModel
from Snappy.SnapInputBomb import SnapInputBomb, ExplosionType
Expand Down Expand Up @@ -107,9 +110,7 @@ def start(self):
self.snap_controller.snapRunning = "running"
debug("started: " + self.snap_controller.snapRunning)
else:
self.snap_controller.write_log(
"starting bsnap_naccident snap.input failed"
)
self.snap_controller.write_log("starting bsnap_naccident snap.input failed")


class SnapController:
Expand All @@ -125,7 +126,9 @@ def __init__(self):

def write_log(self, txt: str):
debug(txt)
self.main.evaluate_javaScript("updateSnapLog({0});".format(json.dumps(html.escape(txt))))
self.main.evaluate_javaScript(
"updateSnapLog({0});".format(json.dumps(html.escape(txt)))
)

def _snap_finished(self):
debug("finished")
Expand Down Expand Up @@ -174,7 +177,7 @@ def _met_finished_run_snap(self):
self.res.getSnapInputMetDefinitions(
self.lastQDict["metmodel"],
self.metcalc.get_meteorology_files(),
**metdefs
**metdefs,
)
)
self._snap_model_run()
Expand Down Expand Up @@ -288,7 +291,7 @@ def _meps25DomainCheck(self, lonf, latf):
return False
return True

def get_bomb_release(self, qDict):
def get_bomb_release(self, qDict, offset_minutes: int):
errors = ""
try:
yld = int(qDict["yield"])
Expand All @@ -300,28 +303,23 @@ def get_bomb_release(self, qDict):
errors += f"unknown explosion_type: {ex}\n"

sib = SnapInputBomb(yld, explosion_type)
sib.minutes = offset_minutes
return (sib.snap_input(), errors)

def get_isotope_release(self, qDict):
def get_isotope_release(self, qDict, offset_minutes: int):
errors = ""
for tag in ("releaseTime", "radius", "lowerHeight", "upperHeight"):
if not re.search(r"\d+", qDict[tag]):
errors += "Cannot interprete {}: {}".format(tag, qDict[tag])

source_tmpl = """
source_term = f"""
MAX.PARTICLES.PER.RELEASE= 2000
TIME.RELEASE.PROFILE.STEPS
RELEASE.HOUR= 0, {releaseTime}
RELEASE.RADIUS.M= {radius}, {radius}
RELEASE.LOWER.M= {lowerHeight}, {lowerHeight}
RELEASE.UPPER.M= {upperHeight}, {upperHeight}
RELEASE.HOUR= {offset_minutes/60:.2f}, {int(qDict["releaseTime"])+offset_minutes/60:.2f}
RELEASE.RADIUS.M= {qDict["radius"]}, {qDict["radius"]}
RELEASE.LOWER.M= {qDict["lowerHeight"]}, {qDict["lowerHeight"]}
RELEASE.UPPER.M= {qDict["upperHeight"]}, {qDict["upperHeight"]}
"""
source_term = source_tmpl.format(
releaseTime=qDict["releaseTime"],
radius=qDict["radius"],
lowerHeight=qDict["lowerHeight"],
upperHeight=qDict["upperHeight"],
)

isotopes = {"relI131": "I131", "relXE133": "Xe133", "relCS137": "Cs137"}
for rel, iso in isotopes.items():
Expand Down Expand Up @@ -353,6 +351,11 @@ def run_snap_query(self, qDict):
if match:
startTime = "{0} {1} {2} {3}".format(*match.group(1, 2, 3, 4))
startDT = datetime.datetime(*tuple(map(int, list(match.group(1, 2, 3, 4)))))
offset_minutes = 0
if match_min := re.search(
r"\d{4}-\d{2}-\d{2}[\+\s]+\d{1,2}:(\d{1,2})", qDict["startTime"]
):
offset_minutes = int(match_min.group(1))
else:
errors += "Cannot interprete startTime: {0}\n".format(qDict["startTime"])

Expand All @@ -378,8 +381,10 @@ def run_snap_query(self, qDict):
except ValueError as ve:
latf = 0.0
lonf = 0.0
errors += "Cannot interprete latitude/longitude: {lat}/{lon}: {ex}\n".format(
lat=lat, lon=lon, ex=ve
errors += (
"Cannot interprete latitude/longitude: {lat}/{lon}: {ex}\n".format(
lat=lat, lon=lon, ex=ve
)
)

if len(errors) > 0:
Expand Down Expand Up @@ -416,9 +421,9 @@ def run_snap_query(self, qDict):
)

if "isBomb" in qDict:
(term, errors) = self.get_bomb_release(qDict)
(term, errors) = self.get_bomb_release(qDict, offset_minutes)
else:
(term, errors) = self.get_isotope_release(qDict)
(term, errors) = self.get_isotope_release(qDict, offset_minutes)
if len(errors) > 0:
debug('updateSnapLog("{0}");'.format(json.dumps("ERRORS:\n\n" + errors)))
self.write_log("ERRORS:\n\n{0}".format(errors))
Expand Down Expand Up @@ -465,12 +470,22 @@ def run_snap_query(self, qDict):
elif qDict["metmodel"] == MetModel.EC0p1Global:
try:
globalRes = EcMeteorologyCalculator.getGlobalMeteoResources()
files = [x[1] for x in sorted(MeteorologyCalculator.findAllGlobalData(globalRes), key=lambda x: x[0])]
files = [
x[1]
for x in sorted(
MeteorologyCalculator.findAllGlobalData(globalRes),
key=lambda x: x[0],
)
]
lat0 = MeteorologyCalculator.getLat0(latf, globalRes.domainHeight)
lon0 = MeteorologyCalculator.getLon0(lonf, globalRes.domainWidth)
with open(os.path.join(self.lastOutputDir, "snap.input"), "a") as fh:
interpol = f"FIMEX.INTERPOLATION=nearest|+proj=latlon +R=6371000 +no_defs|{lon0},{lon0+0.2},...,{lon0+globalRes.domainWidth}|{lat0},{lat0+0.2},...,{lat0+globalRes.domainHeight}|degree\n"
fh.write(self.res.getSnapInputMetDefinitions(qDict["metmodel"], files, interpolation=interpol))
fh.write(
self.res.getSnapInputMetDefinitions(
qDict["metmodel"], files, interpolation=interpol
)
)
self._snap_model_run()
except MeteoDataNotAvailableException as e:
self.write_log("problems finding global EC-met: {}".format(e.args[0]))
Expand Down Expand Up @@ -502,7 +517,10 @@ def run_snap_query(self, qDict):
with open(os.path.join(self.lastOutputDir, "snap.input"), "a") as fh:
fh.write(self.res.getSnapInputMetDefinitions(qDict["metmodel"], files))
self._snap_model_run()
elif qDict["metmodel"] == MetModel.GfsGribFilter or qDict["metmodel"] == MetModel.EC0p1Europe:
elif (
qDict["metmodel"] == MetModel.GfsGribFilter
or qDict["metmodel"] == MetModel.EC0p1Europe
):
files = self.res.getMeteorologyFiles(
qDict["metmodel"], startDT, int(qDict["runTime"]), "best"
)
Expand Down
31 changes: 6 additions & 25 deletions utils/SnapPy/Snappy/SnapInputBomb.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -384,31 +384,13 @@ def snap_input(self) -> str:

lines.append(f"** Explosive yield {self.nuclear_yield}ktonnes")
lines.append("TIME.RELEASE.PROFILE.BOMB")
relradius = []
rellower = []
relupper = []
relstem = []
activity = []
relmins = [0]
if self.minutes > 0:
relmins.append(self.minutes)
relradius.append(0)
rellower.append(0)
relupper.append(0)
relstem.append(0)
activity.append("0.")
relradius.append(self.cloud_radius)
rellower.append(self.cloud_bottom)
relupper.append(self.cloud_top)
relstem.append(self.stem_radius)

lines.append(
f"""
RELEASE.MINUTE= {",".join(map(str, relmins))}
RELEASE.RADIUS.M= {",".join(map(str, relradius))}
RELEASE.LOWER.M= {",".join(map(str, rellower))}
RELEASE.UPPER.M= {",".join(map(str, relupper))}
RELEASE.MUSHROOM.STEM.RADIUS.M= {",".join(map(str, relstem))}
RELEASE.MINUTE= {self.minutes}
RELEASE.RADIUS.M= {self.cloud_radius}
RELEASE.LOWER.M= {self.cloud_bottom}
RELEASE.UPPER.M= {self.cloud_top}
RELEASE.MUSHROOM.STEM.RADIUS.M= {self.stem_radius}
"""
)

Expand Down Expand Up @@ -448,9 +430,8 @@ def snap_input(self) -> str:
)

for i, frac in enumerate(self.size_distribution):
size_activity = activity + [f"{self.activity_after_1hour*frac:.3E}"]
lines.append(
f"RELEASE.BQ/STEP.COMP= {','.join(size_activity)} '{self.component_name(i)}'"
f"RELEASE.BQ/STEP.COMP= {self.activity_after_1hour*frac:.3E} '{self.component_name(i)}'"
)

return "\n".join(lines) + "\n"
Expand Down

0 comments on commit 02c0476

Please sign in to comment.