diff --git a/src/common/decay.f90 b/src/common/decay.f90 index 81ef176f..c19a2192 100644 --- a/src/common/decay.f90 +++ b/src/common/decay.f90 @@ -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. @@ -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 diff --git a/src/common/drydep.f90 b/src/common/drydep.f90 index 037540a5..a7182ea4 100644 --- a/src/common/drydep.f90 +++ b/src/common/drydep.f90 @@ -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 @@ -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 diff --git a/src/common/wetdep.f90 b/src/common/wetdep.f90 index ea70bb7c..62f5423e 100644 --- a/src/common/wetdep.f90 +++ b/src/common/wetdep.f90 @@ -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 @@ -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 diff --git a/utils/SnapPy/Snappy/SnapController.py b/utils/SnapPy/Snappy/SnapController.py index e3e10db8..efe72eaf 100644 --- a/utils/SnapPy/Snappy/SnapController.py +++ b/utils/SnapPy/Snappy/SnapController.py @@ -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 @@ -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: @@ -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") @@ -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() @@ -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"]) @@ -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(): @@ -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"]) @@ -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: @@ -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)) @@ -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])) @@ -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" ) diff --git a/utils/SnapPy/Snappy/SnapInputBomb.py b/utils/SnapPy/Snappy/SnapInputBomb.py index 893ee627..b5028502 100644 --- a/utils/SnapPy/Snappy/SnapInputBomb.py +++ b/utils/SnapPy/Snappy/SnapInputBomb.py @@ -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} """ ) @@ -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"