Merge pull request #141 from metno/bomb_dose_fix

fix bomb-dose decay formula to start at H+1

Author: heikoklein

src/common/decay.f90

@@ -45,22 +45,24 @@ end subroutine decay
 !>  Purpose:  Decrease radioactive contents of deposition fields
 !>            due to decay
 !>
+!>     param: tstep model timestep in seconds
 !>     NEEDS TO BE RUN BEFORE ::decay
 subroutine decayDeps(tstep)
   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 iso_fortran_env, only: real64
+  USE releaseML, only: tpos_bomb
 
   real, intent(in) :: tstep
 
   integer :: m, mm
   real :: bomb_decay_rate, current_state, next_state
   logical, save :: prepare = .TRUE.
   logical, save :: has_bomb_decay = .FALSE.
-  !> totalstep in seconds run in decay
+  !> totalstep in seconds run in decay, assumed to be H+1 after release (stable cloud)
   !> start at 1h to satisfy C(t) = C_0 * t^-1.2 (t in [hrs])
-  real(real64), save :: total_steps = 3600.
+  real(real64), save :: total_steps = 0.
 
   if(prepare) then
   !..radioactive decay rate
@@ -78,13 +80,20 @@ subroutine decayDeps(tstep)
     prepare= .FALSE.
   end if
 
-  ! bomb t[h]^-1.2 power-function, 
+  ! bomb t[h]^-1.2 power-function,
   ! see glassstone/dolan: effects of nuclear weapons
   ! converted  to decay-rate factor per tstep
+  ! initial release to be defined at H+1
   if (has_bomb_decay) then
-    current_state = (total_steps/3600.)**(-1.2)
-    next_state = ((total_steps+tstep)/3600.)**(-1.2)
-    bomb_decay_rate = next_state/current_state
+    if (total_steps >= (3600+tpos_bomb)) then
+      ! start decay after 1h after explosion
+      current_state = ((total_steps-tpos_bomb)/3600.)**(-1.2)
+      next_state = ((total_steps+tstep-tpos_bomb)/3600.)**(-1.2)
+      bomb_decay_rate = next_state/current_state
+    else
+      ! no decay before release+1h
+      bomb_decay_rate = 1.
+    end if
     do m=1,ncomp
       mm = run_comp(m)%to_defined
       if (def_comp(mm)%kdecay == 2) then

src/common/release.f90

@@ -99,6 +99,9 @@ module releaseML
 !> max. no. of particles, total in all plumes, preset, can be configured in snap.input
   integer, save, public :: mpart = mpartpre
 
+!> seconds since start of run when bomb-release occurs
+  integer, save, public :: tpos_bomb = 0
+
   contains
 
 subroutine release(istep,nsteph,tf1,tf2,tnow,ierror)

src/common/snap.F90

@@ -169,7 +169,7 @@ PROGRAM bsnap
   USE snapmetML, only: init_meteo_params, met_params
   USE snapparML, only: component, run_comp, &
                        ncomp, def_comp, nparnum, &
-                       time_profile
+                       time_profile, TIME_PROFILE_BOMB
   USE snapposML, only: irelpos, nrelpos, release_positions
   USE snapgrdML, only: modleveldump, ivcoor, &
                        klevel, imslp, itotcomp, gparam, &
@@ -191,7 +191,7 @@ PROGRAM bsnap
   USE decayML, only: decay, decayDeps
   USE posintML, only: posint, posint_init
   USE bldpML, only: bldp
-  USE releaseML, only: release, releases, nrelheight, mprel, &
+  USE releaseML, only: release, releases, tpos_bomb, nrelheight, mprel, &
                        mplume, nplume, iplume, npart, mpart, release_t
   USE init_random_seedML, only: init_random_seed
   USE compheightML, only: compheight
@@ -524,13 +524,24 @@ PROGRAM bsnap
       ncsummary = trim(ncsummary)//". Release "//trim(def_comp(m)%compname) &
                   //" (hour, Bq/s): "
       do i = 1, ntprof
-        write (iulog, *) '  hour,Bq/hour: ', &
-          releases(i)%frelhour, (releases(i)%relbqsec(n, ih)*3600., ih=1, nrelheight)
+        if (time_profile /= TIME_PROFILE_BOMB) then
+          write (iulog, *) '  hour,Bq/hour: ', &
+            releases(i)%frelhour, (releases(i)%relbqsec(n, ih)*3600., ih=1, nrelheight)
+        else
+          write (iulog, *) '  hour,Bq: ', &
+            releases(i)%frelhour, (releases(i)%relbqsec(n, ih), ih=1, nrelheight)
+          if (tpos_bomb == 0) then
+            if (any(releases(i)%relbqsec > 0.)) tpos_bomb = releases(i)%frelhour * 3600
+          end if
+          end if
         write (tempstr, '("(",f5.1,",",ES9.2,")")') &
           releases(i)%frelhour, releases(i)%relbqsec(n, 1)
         ncsummary = trim(ncsummary)//" "//trim(tempstr)
       end do
     end do
+    if (time_profile == TIME_PROFILE_BOMB) then
+      write (iulog, *) 'tpos_bomb:  ', tpos_bomb
+    end if
     write (iulog, *) 'itotcomp:   ', itotcomp
     write (iulog, *) 'blfulmix:   ', blfullmix
     write (error_unit, *) 'Title:      ', trim(nctitle)

src/common/snap.mk

@@ -59,7 +59,7 @@ bldp.o: ../common/bldp.f90 snapdimML.o snaptabML.o snapgrdML.o snapfldML.o ftest
 	${F77} -c ${F77FLAGS} $(INCLUDES) $<
 particleML.o: ../common/particleML.f90
 	${F77} -c ${F77FLAGS} $<
-decay.o: ../common/decay.f90 particleML.o snapfldML.o snapparML.o
+decay.o: ../common/decay.f90 particleML.o snapfldML.o snapparML.o release.o
 	${F77} -c ${F77FLAGS} $(INCLUDES) $<
 compheight.o: ../common/compheight.f90 snapgrdML.o snapfldML.o snaptabML.o snapdimML.o ftest.o
 	${F77} -c ${F77FLAGS} $(INCLUDES) $<