allow model-level output

src/common/allocateFields.f90

@@ -20,7 +20,7 @@ module allocateFieldsML
   USE snapparML, only: ncomp, iparnum
   USE snapfldML, only: u1, u2, v1, v2, w1, w2, bl1, bl2, t1, t2, &
       ps1, ps2, hbl1, hbl2, hlevel1, hlevel2, hlayer1, hlayer2, &
-      concacc, avgbq1, avgbq2, instmlbq, accwet, accdry, concen, &
+      concacc, avgbq1, avgbq2, ml_bq, accwet, accdry, concen, &
       depdry, depwet, accprec, avgprec, avghbl, precip, &
       pmsl1, pmsl2, field1, field2, field3, field4, field3d1, xm, ym, &
       garea, field_hr1, field_hr2, field_hr3, hbl_hr, &
@@ -181,7 +181,7 @@ subroutine allocateFields
   IF (AllocateStatus /= 0) ERROR STOP errmsg
 
   if (imodlevel) then
-    ALLOCATE ( instmlbq(nxhr,nyhr,nk-1,nocomp), STAT = AllocateStatus)
+    ALLOCATE ( ml_bq(nxhr,nyhr,nk-1,nocomp), STAT = AllocateStatus)
     IF (AllocateStatus /= 0) ERROR STOP errmsg
   endif
 
@@ -290,8 +290,8 @@ subroutine deAllocateFields
   DEALLOCATE ( avgbq1 )
   DEALLOCATE ( avgbq2 )
 
-  if (allocated(instmlbq)) then
-    deallocate(instmlbq)
+  if (allocated(ml_bq)) then
+    deallocate(ml_bq)
   endif
 
   if (allocated(max_column_concentration)) then

src/common/fldout_nc.f90

@@ -106,7 +106,7 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
     ierror)
   USE iso_fortran_env, only: int16
   USE array_utils, only: is_in_array
-  USE snapgrdML, only: imodlevel, imslp, precipitation_in_output, &
+  USE snapgrdML, only: imodlevel, modlevel_is_average, imslp, precipitation_in_output, &
       itotcomp, compute_column_max_conc, compute_aircraft_doserate, &
       aircraft_doserate_threshold, output_column
   USE snapfldML, only: field_hr1, field_hr2, field_hr3, hbl_hr, &
@@ -536,7 +536,7 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
 
 !..model level fields...................................................
   if (imodlevel) then
-    call write_ml_fields(iunit, varid, [1, 1, -1, ihrs_pos], &
+    call write_ml_fields(iunit, varid, average, [1, 1, -1, ihrs_pos], &
         [nx*output_resolution_factor, ny*output_resolution_factor, 1, 1], rt1, rt2)
   endif
 
@@ -555,19 +555,20 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
 end subroutine fldout_nc
 
 
-subroutine write_ml_fields(iunit, varid, ipos_in, isize, rt1, rt2)
+subroutine write_ml_fields(iunit, varid, average, ipos_in, isize, rt1, rt2)
   USE releaseML, only: nplume, iplume
   USE particleML, only: pdata, Particle
   USE snapparML, only: def_comp, nocomp
   USE snapfldML, only: field_hr1, field_hr2, &
-      hlayer1, hlayer2, garea, instmlbq
+      hlayer1, hlayer2, garea, ml_bq
   USE ftestML, only: ftest
   USE snapdebug, only: idebug
-  USE snapgrdML, only: itotcomp, modleveldump, ivlayer
+  USE snapgrdML, only: itotcomp, modlevel_is_average, modleveldump, ivlayer
   USE snapdimML, only: nx,ny,nk,output_resolution_factor, hres_pos
 
   integer, intent(in) :: iunit
   type(common_var), intent(in) :: varid
+  real, intent(in) :: average
   integer, intent(in) :: ipos_in(4)
   integer, intent(in) :: isize(4)
   real, intent(in) :: rt1, rt2
@@ -589,34 +590,38 @@ subroutine write_ml_fields(iunit, varid, ipos_in, isize, rt1, rt2)
   total=0.
   maxage=0
 
-  instmlbq = 0.0
-
-  do npl = 1, nplume
-    do n = iplume(npl)%start, iplume(npl)%end
-      part = pdata(n)
-      i = hres_pos(part%x)
-      j = hres_pos(part%y)
-      ivlvl = part%z*10000.
-      k = ivlayer(ivlvl)
-      m = def_comp(part%icomp)%to_output
-    !..in each sigma/eta (input model) layer
-      if (modleveldump > 0) then
-      !.. dump and remove old particles, don't touch  new ones
-        if (iplume(npl)%ageInSteps >= nint(modleveldump)) then
-          maxage = max(maxage, int(iplume(npl)%ageInSteps, kind(maxage)))
-          instmlbq(i,j,k,m) = instmlbq(i,j,k,m) + part%rad()
-          total = total + part%rad()
-
-          inactivated_ =  part%inactivate()
-          pdata(n) = part
-        end if
-      else
-        instmlbq(i,j,k,m)=instmlbq(i,j,k,m)+pdata(n)%rad()
-      endif
+  ml_bq = 0.0
+  if (modlevel_is_average) then
+    avg = average
+  else
+    avg = 1.0
+    do npl = 1, nplume
+      do n = iplume(npl)%start, iplume(npl)%end
+        part = pdata(n)
+        i = hres_pos(part%x)
+        j = hres_pos(part%y)
+        ivlvl = part%z*10000.
+        k = ivlayer(ivlvl)
+        m = def_comp(part%icomp)%to_output
+      !..in each sigma/eta (input model) layer
+        if (modleveldump > 0) then
+        !.. dump and remove old particles, don't touch  new ones
+          if (iplume(npl)%ageInSteps >= nint(modleveldump)) then
+            maxage = max(maxage, int(iplume(npl)%ageInSteps, kind(maxage)))
+            ml_bq(i,j,k,m) = ml_bq(i,j,k,m) + part%rad()
+            total = total + part%rad()
+
+            inactivated_ =  part%inactivate()
+            pdata(n) = part
+          end if
+        else
+          ml_bq(i,j,k,m)=ml_bq(i,j,k,m)+pdata(n)%rad()
+        endif
+      end do
     end do
-  end do
-  if (modleveldump > 0) then
-    write (error_unit,*) "dumped; maxage, total", maxage, total
+    if (modleveldump > 0) then
+      write (error_unit,*) "dumped; maxage, total", maxage, total
+    endif
   endif
 
   do k=1,nk-1
@@ -627,32 +632,40 @@ subroutine write_ml_fields(iunit, varid, ipos_in, isize, rt1, rt2)
       end do
     end do
     do m=1,nocomp
-      instmlbq(:,:,k,m) = instmlbq(:,:,k,m)/field_hr2
+      ml_bq(:,:,k,m) = ml_bq(:,:,k,m)/(field_hr2*avg)
     end do
   end do
 
-!..average concentration in each layer for each type
+!..concentration in each layer for each type
   do m=1,nocomp
     do k=1,nk-1
-      field_hr1(:,:) = cscale*instmlbq(:,:,k,m)
-      if(idebug == 1) call ftest('instconcml', field_hr1)
+      field_hr1(:,:) = cscale*ml_bq(:,:,k,m)
+      if (modlevel_is_average) then
+        if(idebug == 1) call ftest('tinstconcml', field_hr1)
+      else
+        if(idebug == 1) call ftest('tavgconcml', field_hr1)
+      end if
 
       ipos(3) = k
       call check(nf90_put_var(iunit, varid%comp(m)%icml, start=ipos, &
           count=isize, values=field_hr1), "icml(m)")
     end do
   end do
 
-!..total average concentration in each layer
+!..total concentration in each layer
   if(nocomp > 1 .AND. itotcomp == 1) then
     do m=2,nocomp
       do k=1,nk-1
-        instmlbq(:,:,k,1) = instmlbq(:,:,k,1) + instmlbq(:,:,k,m)
+        ml_bq(:,:,k,1) = ml_bq(:,:,k,1) + ml_bq(:,:,k,m)
       end do
     end do
     do k=1,nk-1
-      field_hr1(:,:) = cscale*instmlbq(:,:,k,1)
-      if(idebug == 1) call ftest('tinstconcml', field_hr1)
+      field_hr1(:,:) = cscale*ml_bq(:,:,k,1)
+      if (modlevel_is_average) then
+        if(idebug == 1) call ftest('tinstconcml', field_hr1)
+      else
+        if(idebug == 1) call ftest('tavgconcml', field_hr1)
+      end if
     end do
   end if
 end subroutine
@@ -972,7 +985,7 @@ end subroutine nc_set_projection
 
 subroutine initialize_output(filename, itime, ierror)
   USE snapfilML, only: ncsummary, nctitle, simulation_start
-  USE snapgrdML, only: gparam, igtype, imodlevel, imslp, precipitation_in_output, &
+  USE snapgrdML, only: gparam, igtype, imodlevel, modlevel_is_average, imslp, precipitation_in_output, &
       itotcomp, modleveldump, compute_column_max_conc, compute_aircraft_doserate, &
       aircraft_doserate_threshold, output_column
   USE snapfldML, only:  &
@@ -1129,6 +1142,11 @@ subroutine initialize_output(filename, itime, ierror)
         endif
         call nc_declare(iunit, dimids4d, varid%comp(m)%icml, &
           string, units="Bq/m3", chunksize=chksz4d)
+        if (modlevel_is_average) then
+          call check(nf90_put_att(iunit, varid%comp(m)%icml, "cell_method", "time: mean"))
+        else
+          call check(nf90_put_att(iunit, varid%comp(m)%icml, "cell_method", "time: point"))
+        end if
       end if
       if (output_column) then
         call nc_declare(iunit, dimids3d, varid%comp(m)%conc_column, &
@@ -1297,13 +1315,13 @@ subroutine get_varids(iunit, varid, ierror)
 
 !> accumulation for average fields
 subroutine accumulate_fields(tf1, tf2, tnow, tstep, nsteph)
-  USE snapgrdML, only: imodlevel, &
+  USE snapgrdML, only: imodlevel, modlevel_is_average, &
       ivlayer, compute_column_max_conc, compute_aircraft_doserate, &
       alevel, blevel, aircraft_doserate_threshold
   USE snapfldML, only:  &
       avgbq1, avgbq2, hlayer1, hlayer2, hbl1, hbl2, &
       avgprec, concen, avghbl, &
-      concacc, precip, &
+      concacc, precip, ml_bq, &
       max_column_scratch, max_column_concentration, garea, &
       ps1, ps2, t1_abs, t2_abs, aircraft_doserate_scratch, aircraft_doserate, &
       aircraft_doserate_threshold_height
@@ -1332,6 +1350,10 @@ subroutine accumulate_fields(tf1, tf2, tnow, tstep, nsteph)
     avgbq1 = 0.0
     avgbq2 = 0.0
 
+    if (imodlevel .and. modlevel_is_average) then
+      ml_bq = 0.0
+    end if
+
     if (compute_column_max_conc) then
       max_column_concentration = 0.0
     endif
@@ -1423,6 +1445,19 @@ subroutine accumulate_fields(tf1, tf2, tnow, tstep, nsteph)
 
   end block
 
+  if(imodlevel .and. modlevel_is_average) then
+    do n=1,npart
+      part = pdata(n)
+      i = hres_pos(part%x)
+      j = hres_pos(part%y)
+      ivlvl = part%z*10000.
+      k = ivlayer(ivlvl)
+      m = def_comp(part%icomp)%to_running
+    !..in each sigma/eta (input model) layer
+      ml_bq(i,j,k,m) = ml_bq(i,j,k,m) + part%rad()
+    end do
+  end if
+
   if (compute_column_max_conc) then
     max_column_scratch = 0.0
     do n=1,npart

src/common/snap.F90

@@ -106,8 +106,8 @@
 ! *
 ! * default is 30 meters, can set specific value by setting the parameter
 ! SURFACE.LAYER.CONCENTRATION.AT.HEIGHT (= 30)
+! * MODEL.LEVEL.FIELDS can be ON/INST/AVG/OFF (where ON == INST)
 ! MODEL.LEVEL.FIELDS.ON
-! MODEL.LEVEL.FIELDS.OFF
 ! * only write particles to model-level, which are at least DUMPTIME in h
 ! * old. If DUMPTIME > 0, the dumped particles will be removed from model
 ! MODEL.LEVEL.FIELDS.DUMPTIME= 4.0
@@ -175,7 +175,7 @@ PROGRAM bsnap
   USE snapposML, only: irelpos, nrelpos, release_positions
   USE snapgrdML, only: modleveldump, ivcoor, &
                        klevel, imslp, itotcomp, gparam, &
-                       igtype, imodlevel, precipitation_in_output
+                       igtype, imodlevel, modlevel_is_average, precipitation_in_output
   USE snaptabML, only: tabcon
   USE particleML, only: pdata, extraParticle
   USE allocateFieldsML, only: allocateFields, deallocateFields
@@ -1515,6 +1515,14 @@ subroutine read_inputfile(snapinput_unit)
       case ('model.level.fields.on')
         !..model.level.fields.on
         imodlevel = .true.
+      case ('model.level.fields.inst')
+        !..model.level.fields.off
+        imodlevel = .true.
+        modlevel_is_average = .false.
+      case ('model.level.fields.avg')
+        !..model.level.fields.off
+        imodlevel = .true.
+        modlevel_is_average = .true.
       case ('model.level.fields.off')
         !..model.level.fields.off
         imodlevel = .false.

src/common/snapfldML.f90

@@ -137,7 +137,7 @@ module snapfldML
 !>
 !> only used if imodlevel
 !> This is not weighted by area
-  real(kind=real64), allocatable, save, public :: instmlbq(:,:,:,:)
+  real(kind=real64), allocatable, save, public :: ml_bq(:,:,:,:)
 
 !> Scratch space for finding max column concentration
 !>

src/common/snapgrdML.f90

@@ -40,6 +40,8 @@ module snapgrdML
     integer, save, public :: imslp = 0
 !> * output of concentrations in model levels
     logical, save, public :: imodlevel = .false.
+!> * output of concentrations in model levels average/instantaneous
+    logical, save, public :: modlevel_is_average = .false.
 !> remove particles from model after at least that many steps
     real, save, public :: modleveldump = 0.0
 !> * 0=not output of total of all components (e.g. when each component is released with a mass unit)