From 48e13974300e7e2e05ab9130013da119a463b843 Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Thu, 13 Nov 2025 14:11:22 +0000
Subject: [PATCH 1/6] Enable emerson drydep scheme for ecemep meteorology

- Allow scalar reynold stress
- Allow non-accumulated vector fluxes
---
 src/common/allocateFields.f90 |  6 +--
 src/common/drydep.f90         |  9 ++--
 src/common/fldout_nc.f90      |  5 +-
 src/common/readfield_fi.f90   | 86 +++++++++++++++++++++++++----------
 src/common/readfield_nc.f90   |  7 ++-
 src/common/snapfldML.f90      |  2 +-
 src/common/snapmetML.f90      | 16 ++++++-
 7 files changed, 93 insertions(+), 38 deletions(-)

diff --git a/src/common/allocateFields.f90 b/src/common/allocateFields.f90
index 7af7e589..599a988f 100644
--- a/src/common/allocateFields.f90
+++ b/src/common/allocateFields.f90
@@ -28,7 +28,7 @@ module allocateFieldsML
       max_column_scratch, max_column_concentration, &
       aircraft_doserate, aircraft_doserate_scratch, t1_abs, t2_abs, &
       aircraft_doserate_threshold_height, vd_dep, &
-      xflux, yflux, hflux, t2m, z0, &
+      surface_stress, xflux, yflux, hflux, t2m, z0, &
       ustar, raero, my, &
       total_activity_released, total_activity_lost_domain, total_activity_lost_other, &
       wscav, cloud_cover
@@ -242,7 +242,7 @@ subroutine allocateFields
     if (requires_extra_fields_to_be_read()) then
       allocate(vd_dep(nx,ny,ncomp), STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
-      allocate(xflux, yflux, hflux, t2m, z0, mold=ps2, STAT=AllocateStatus)
+      allocate(surface_stress, xflux, yflux, hflux, t2m, z0, mold=ps2, STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
       allocate(raero(nx, ny), STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
@@ -350,7 +350,7 @@ subroutine deAllocateFields
   DEALLOCATE( total_activity_released, total_activity_lost_domain, total_activity_lost_other )
   if (allocated(vd_dep)) then
     deallocate(vd_dep)
-    deallocate(xflux, yflux, hflux, t2m, z0)
+    deallocate(surface_stress, xflux, yflux, hflux, t2m, z0)
     deallocate(ustar, raero)
   endif
 
diff --git a/src/common/drydep.f90 b/src/common/drydep.f90
index 31755389..040db462 100644
--- a/src/common/drydep.f90
+++ b/src/common/drydep.f90
@@ -77,16 +77,15 @@ pure logical function requires_extra_fields_to_be_read()
 
 !> Precompute dry deposition meteorology
 !> returns roa, ustar, monin_obukhov_length, raero, my
-elemental subroutine drydep_precompute_meteo(surface_pressure, t2m, yflux, xflux, z0, hflux, &
+elemental subroutine drydep_precompute_meteo(surface_pressure, t2m, surface_stress, z0, hflux, &
                                              ustar, raero, my)
   use iso_fortran_env, only: real64
   use datetime, only: datetime_t
   real, intent(in) :: surface_pressure !> [Pa]
   real, intent(in) :: t2m !> [K]
-  real, intent(in) :: yflux !> [N/m^2]
-  real, intent(in) :: xflux !> [N/m^2]
+  real, intent(in) :: surface_stress !> [N/m^2]
   real, intent(in) :: z0 !> [m]
-  real, intent(in) :: hflux !> [W s/m^2]
+  real, intent(in) :: hflux !> [W/m^2]
   real(real64), intent(out) :: raero, ustar, my
 
   real(real64), parameter :: k = 0.4
@@ -94,7 +93,7 @@ elemental subroutine drydep_precompute_meteo(surface_pressure, t2m, yflux, xflux
   real(real64) :: roa, monin_obukhov_length
 
   roa = surface_pressure / (t2m * R)
-  ustar = hypot(yflux, xflux) / sqrt(roa)
+  ustar = surface_stress / sqrt(roa)
   monin_obukhov_length = - roa * CP * t2m * (ustar**3) / (k * grav * hflux)
   raero = aerodynres(monin_obukhov_length, ustar, real(z0, real64))
   my = ny * roa
diff --git a/src/common/fldout_nc.f90 b/src/common/fldout_nc.f90
index cbd4b87c..747f56ed 100644
--- a/src/common/fldout_nc.f90
+++ b/src/common/fldout_nc.f90
@@ -82,6 +82,7 @@ module fldout_ncML
     integer :: aircraft_doserate_threshold_height
     integer :: components
     type(component_var) :: comp(mcomp)
+    integer :: surface_stress = -1
     integer :: xflux = -1
     integer :: yflux = -1
     integer :: hflux = -1
@@ -596,12 +597,14 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
 
   if (output_vd_debug) then
     block
-      use snapfldml, only: t2m, xflux, yflux, z0, hflux, &
+      use snapfldml, only: t2m, surface_stress, xflux, yflux, z0, hflux, &
         ustar, raero, ps2
       call hres_field(ps2, field_hr1)
       call check(nf90_put_var(iunit, varid%ps_vd, start=ipos, count=isize, values=field_hr1))
       call hres_field(t2m, field_hr1)
       call check(nf90_put_var(iunit, varid%t2m, start=ipos, count=isize, values=field_hr1))
+      call hres_field(surface_stress, field_hr1)
+      call check(nf90_put_var(iunit, varid%surface_stress, start=ipos, count=isize, values=field_hr1))
       call hres_field(xflux, field_hr1)
       call check(nf90_put_var(iunit, varid%xflux, start=ipos, count=isize, values=field_hr1))
       call hres_field(yflux, field_hr1)
diff --git a/src/common/readfield_fi.f90 b/src/common/readfield_fi.f90
index 0178c0ab..35816af8 100644
--- a/src/common/readfield_fi.f90
+++ b/src/common/readfield_fi.f90
@@ -896,11 +896,12 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     USE iso_fortran_env, only: real64
     USE snapmetML, only: met_params, &
       temp_units, downward_momentum_flux_units, surface_roughness_length_units, &
-      surface_heat_flux_units
+      accumulated_surface_heat_flux_units, surface_heat_flux_units
     use drydepml, only: drydep_precompute_meteo, drydep_precompute_particle, &
       requires_extra_fields_to_be_read, classnr
     use snapparML, only: ncomp, run_comp, def_comp
-    use snapfldML, only: ps2, vd_dep, xflux, yflux, hflux, z0, t2m, &
+    ! TODO1: make xflux and yflux temporary variables in readfield
+    use snapfldML, only: ps2, vd_dep, surface_stress, xflux, yflux, hflux, z0, t2m, &
       ustar, raero, my
     use snaptimers, only: metcalc_timer
 
@@ -919,41 +920,76 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
 
     allocate(tmp1, tmp2, MOLD=ps2)
 
-
-    ! Fluxes are integrated: Deaccumulate
-    if (timepos == 1) then
-      call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, xflux(:, :), nt=timepos, nr=nr)
-      call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, yflux(:, :), nt=timepos, nr=nr)
+    ! TODO: refactor
+    if (met_params%surface_stress /= "") then
+      ! Load surface_stress
+      call fi_checkload(fio, met_params%surface_stress, downward_momentum_flux_units, surface_stress(:, :), nt=timepos, nr=nr)
+    else if (met_params%xflux == "" .OR. met_params%yflux == "") then
+      ! Either surface_stress or xflux and yflux needs to be defined
+      error stop "Assertion error: Either surface_stress or xflux and yflux needs to be defined in the meteorological parameters"
     else
-      call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-      call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-      xflux(:,:) = tmp2 - tmp1
+      ! Load and combine surface stress/momentum flux components
+      if (met_params%xflux_is_accumulated) then
+        if (timepos == 1) then
+          ! Note: Arome files have the wrong units for downward_momentum_flux_units, missing a unit of time
+          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, xflux(:, :), nt=timepos, nr=nr)
+        else
+          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
+          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
+          xflux(:,:) = tmp2 - tmp1
+        endif
+        ! TODO: Normalise by difference between intervals
+        xflux(:,:) =  xflux / 3600
+      else
+        call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, xflux(:, :), nt=timepos, nr=nr)
+      endif
 
-      call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-      call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-      yflux(:,:) = tmp2 - tmp1
+      if (met_params%yflux_is_accumulated) then
+        if (timepos == 1) then
+          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, yflux(:, :), nt=timepos, nr=nr)
+        else
+          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
+          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
+          yflux(:,:) = tmp2 - tmp1
+        endif
+        ! TODO: Normalise by difference between intervals
+        yflux(:,:) =  yflux / 3600
+      else
+        call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, yflux(:, :), nt=timepos, nr=nr)
+      endif
+      surface_stress = hypot(yflux, xflux)
     endif
-    ! TODO: Normalise by difference between intervals
-    xflux(:,:) =  xflux / 3600
-    yflux(:,:) =  yflux / 3600
 
-    if (timepos == 1) then
-      call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, hflux(:, :), nt=timepos, nr=nr)
+    if (met_params%hflux_is_accumulated) then
+      if (timepos == 1) then
+        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, hflux(:, :), nt=timepos, nr=nr)
+      else
+        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
+        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, tmp2(:, :), nt=timepos, nr=nr)
+        hflux(:,:) = tmp2 - tmp1
+      endif
+      ! TODO: Normalise by difference between intervals
+      hflux(:,:) = -hflux / 3600 ! Follow conventions for up/down
     else
-      call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-      call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-      hflux(:,:) = tmp2 - tmp1
+      call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, hflux(:, :), nt=timepos, nr=nr)
+      hflux(:,:) = -hflux ! Follow conventions for up/down
     endif
-    ! TODO: Normalise by difference between intervals
-    hflux(:,:) = -hflux / 3600 ! Follow conventions for up/down
 
 
-    call fi_checkload(fio, met_params%z0, surface_roughness_length_units, z0(:, :), nt=timepos, nr=nr)
+    deallocate(tmp1)
+    deallocate(tmp2)
+
+    if (met_params%z0 == "") then
+      write(error_unit,*) "WARNING!!!! Surface roughness set to 1 in liu of loaded value."
+      z0(:,:) = 1
+    else
+      call fi_checkload(fio, met_params%z0, surface_roughness_length_units, z0(:, :), nt=timepos, nr=nr)
+    endif
 
     call fi_checkload(fio, met_params%t2m, temp_units, t2m(:, :), nt=timepos, nr=nr)
 
     call metcalc_timer%start()
-    call drydep_precompute_meteo(ps2*100., t2m, yflux, xflux, z0, hflux, &
+    call drydep_precompute_meteo(ps2*100., t2m, surface_stress, z0, hflux, &
       ustar, raero, my)
     !$OMP PARALLEL DO PRIVATE(i,mm)
     do i=1,ncomp
diff --git a/src/common/readfield_nc.f90 b/src/common/readfield_nc.f90
index b2fa77c1..81066811 100644
--- a/src/common/readfield_nc.f90
+++ b/src/common/readfield_nc.f90
@@ -1177,7 +1177,9 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
     USE iso_fortran_env, only: real64
     use datetime, only: datetime_t
     use snapmetML, only: met_params
-    use snapfldML, only: xflux, yflux, hflux, z0, t2m, vd_dep, ustar, &
+    ! TODO1: make xflux and yflux temporary variables in readfield
+    ! TODO2: apply changes in readfield_fi to readfield_nc
+    use snapfldML, only: surface_stress, xflux, yflux, hflux, z0, t2m, vd_dep, ustar, &
       ps2, raero, my, enspos
     use drydepml, only: classnr, requires_extra_fields_to_be_read, drydep_precompute_meteo, drydep_precompute_particle
     use snapdimML, only: nx, ny
@@ -1221,6 +1223,7 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
     ! TODO: Normalise by difference between intervals
     xflux(:,:) =  xflux / 3600
     yflux(:,:) =  yflux / 3600
+    surface_stress = hypot(yflux, xflux)
 
     if (timepos == 1) then
       call nfcheckload(ncid, met_params%hflux, start, count, hflux(:,:))
@@ -1236,7 +1239,7 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
 
     call nfcheckload(ncid, met_params%t2m, start, count, t2m(:, :))
 
-    call drydep_precompute_meteo(ps2*100., t2m, yflux, xflux, z0, hflux, &
+    call drydep_precompute_meteo(ps2*100., t2m, surface_stress, z0, hflux, &
       ustar, raero, my)
     do i=1,ncomp
       mm = run_comp(i)%to_defined
diff --git a/src/common/snapfldML.f90 b/src/common/snapfldML.f90
index e9e96248..9b24bd29 100644
--- a/src/common/snapfldML.f90
+++ b/src/common/snapfldML.f90
@@ -175,6 +175,6 @@ module snapfldML
   !> Deposition velocity on the grid per species
   real, allocatable, save, public :: vd_dep(:, :, :)
 
-  real, allocatable, save, public :: xflux(:, :), yflux(:, :), hflux(:, :), z0(:, :), t2m(:, :)
+  real, allocatable, save, public :: surface_stress(:, :), xflux(:, :), yflux(:, :), hflux(:, :), z0(:, :), t2m(:, :)
   real(real64), allocatable, save, public :: ustar(:,:), raero(:,:), my(:,:)
 end module snapfldML
diff --git a/src/common/snapmetML.f90 b/src/common/snapmetML.f90
index 2ec9d3a8..c339f3ff 100644
--- a/src/common/snapmetML.f90
+++ b/src/common/snapmetML.f90
@@ -48,10 +48,14 @@ module snapmetML
     character(len=80) :: t2m = ''
     character(len=80) :: yflux = ''
     character(len=80) :: xflux = ''
+    character(len=80) :: surface_stress = ''
     character(len=80) :: hflux = ''
     character(len=80) :: z0 = ''
 
     ! flags when reading the data
+    logical :: xflux_is_accumulated = .false.
+    logical :: yflux_is_accumulated = .false.
+    logical :: hflux_is_accumulated = .false.
     logical :: temp_is_abs = .false.
     logical :: has_dummy_dim = .false.
     logical :: manual_level_selection = .false.
@@ -92,7 +96,8 @@ module snapmetML
 
   character(len=*), parameter, public :: downward_momentum_flux_units = 'N/m^2'
   character(len=*), parameter, public :: surface_roughness_length_units = 'm'
-  character(len=*), parameter, public :: surface_heat_flux_units = 'W s/m^2'
+  character(len=*), parameter, public :: surface_heat_flux_units = 'W/m^2'
+  character(len=*), parameter, public :: accumulated_surface_heat_flux_units = 'W s/m^2'
 
   character(len=*), parameter, public :: cloud_fraction_units = '%'
   character(len=*), parameter, public :: mass_fraction_units = 'kg/kg'
@@ -238,10 +243,13 @@ subroutine init_meteo_params(nctype, ierr)
 
       met_params%t2m = 'air_temperature_2m'
       met_params%xflux = 'downward_northward_momentum_flux_in_air'
+      met_params%xflux_is_accumulated = .true.
       met_params%yflux = 'downward_eastward_momentum_flux_in_air'
+      met_params%yflux_is_accumulated = .true.
       ! met_params%z0 = 'surface_roughness_length'
       met_params%z0 = "SFX_Z0"
       met_params%hflux = 'integral_of_surface_downward_sensible_heat_flux_wrt_time'
+      met_params%hflux_is_accumulated = .true.
 
       met_params%mass_fraction_rain_in_air = "mass_fraction_of_rain_in_air_ml"
       met_params%mass_fraction_graupel_in_air = "mass_fraction_of_graupel_in_air_ml"
@@ -314,6 +322,12 @@ subroutine init_meteo_params(nctype, ierr)
       met_params%precstrativrt = 'large_scale_precipitations'
       met_params%precconvrt = 'convective_precipitations'
 
+      met_params%t2m = 'temperature_2m'
+      met_params%surface_stress = 'surface_stress'
+      ! met_params%z0 = 'surface_roughness_length'
+      met_params%z0 = ''
+      met_params%hflux = 'surface_flux_sensible_heat'
+
       met_params%mass_fraction_rain_in_air = "mass_fraction_of_rain_in_air_ml"
       met_params%mass_fraction_graupel_in_air = "mass_fraction_of_graupel_in_air_ml"
       met_params%mass_fraction_snow_in_air = "mass_fraction_of_snow_in_air_ml"

From 704cf258ebea35b9ca5e8342213a4506dc57e69c Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Thu, 13 Nov 2025 14:22:45 +0000
Subject: [PATCH 2/6] Fix misnamed units

---
 src/common/readfield_fi.f90 | 2 +-
 src/common/snapmetML.f90    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/common/readfield_fi.f90 b/src/common/readfield_fi.f90
index 35816af8..51bc21c0 100644
--- a/src/common/readfield_fi.f90
+++ b/src/common/readfield_fi.f90
@@ -896,7 +896,7 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     USE iso_fortran_env, only: real64
     USE snapmetML, only: met_params, &
       temp_units, downward_momentum_flux_units, surface_roughness_length_units, &
-      accumulated_surface_heat_flux_units, surface_heat_flux_units
+      accum_surface_heat_flux_units, surface_heat_flux_units
     use drydepml, only: drydep_precompute_meteo, drydep_precompute_particle, &
       requires_extra_fields_to_be_read, classnr
     use snapparML, only: ncomp, run_comp, def_comp
diff --git a/src/common/snapmetML.f90 b/src/common/snapmetML.f90
index c339f3ff..b3dfb49a 100644
--- a/src/common/snapmetML.f90
+++ b/src/common/snapmetML.f90
@@ -97,7 +97,7 @@ module snapmetML
   character(len=*), parameter, public :: downward_momentum_flux_units = 'N/m^2'
   character(len=*), parameter, public :: surface_roughness_length_units = 'm'
   character(len=*), parameter, public :: surface_heat_flux_units = 'W/m^2'
-  character(len=*), parameter, public :: accumulated_surface_heat_flux_units = 'W s/m^2'
+  character(len=*), parameter, public :: accum_surface_heat_flux_units = 'W s/m^2'
 
   character(len=*), parameter, public :: cloud_fraction_units = '%'
   character(len=*), parameter, public :: mass_fraction_units = 'kg/kg'

From b9f4d7f954af7f7a5a94a30a3b4190826616ce50 Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Thu, 13 Nov 2025 15:35:36 +0000
Subject: [PATCH 3/6] Move accumulation code into subroutine

---
 src/common/readfield_fi.f90 | 67 +++++++++++++++++--------------------
 1 file changed, 30 insertions(+), 37 deletions(-)

diff --git a/src/common/readfield_fi.f90 b/src/common/readfield_fi.f90
index 51bc21c0..90fb59ed 100644
--- a/src/common/readfield_fi.f90
+++ b/src/common/readfield_fi.f90
@@ -891,6 +891,28 @@ subroutine fi_checkload_intern(fio, varname, units, zfield, nt, nz, nr, ierror)
     write (iulog, *) "reading "//trim(varname)//", min, max: ", minval(zfield), maxval(zfield)
   end subroutine fi_checkload_intern
 
+
+  subroutine read_accumulated_field(fio, timepos, timeposm1, varname, units, field, nr)
+    TYPE(FimexIO), intent(inout) :: fio
+    character(len=*), intent(in) :: varname, units
+    integer, intent(in) :: timepos, timeposm1, nr
+    real, allocatable :: tmp1(:, :), tmp2(:, :)
+    real(real32), intent(out) :: field(:, :)
+
+    allocate(tmp1, tmp2, MOLD=field)
+
+    if (timepos == 1) then
+      call fi_checkload(fio, varname, units, field(:, :), nt=timepos, nr=nr)
+    else
+      call fi_checkload(fio, varname, units, tmp1(:, :), nt=timeposm1, nr=nr)
+      call fi_checkload(fio, varname, units, tmp2(:, :), nt=timepos, nr=nr)
+      field(:,:) = tmp2 - tmp1
+    endif
+    ! TODO: Normalise by difference between intervals
+    field(:,:) =  field / 3600
+  end subroutine read_accumulated_field
+
+
   subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     USE ieee_arithmetic, only: ieee_is_nan
     USE iso_fortran_env, only: real64
@@ -911,15 +933,12 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     integer, intent(in) :: nr
     type(datetime_t), intent(in) :: itimefi
 
-    real, allocatable :: tmp1(:, :), tmp2(:, :)
     integer :: i, mm
 
     if (.not.requires_extra_fields_to_be_read()) then
       return
     endif
 
-    allocate(tmp1, tmp2, MOLD=ps2)
-
     ! TODO: refactor
     if (met_params%surface_stress /= "") then
       ! Load surface_stress
@@ -930,57 +949,31 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     else
       ! Load and combine surface stress/momentum flux components
       if (met_params%xflux_is_accumulated) then
-        if (timepos == 1) then
-          ! Note: Arome files have the wrong units for downward_momentum_flux_units, missing a unit of time
-          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, xflux(:, :), nt=timepos, nr=nr)
-        else
-          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-          call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-          xflux(:,:) = tmp2 - tmp1
-        endif
-        ! TODO: Normalise by difference between intervals
-        xflux(:,:) =  xflux / 3600
+        ! Note: Arome files have the wrong units for downward_momentum_flux_units, missing a unit of time
+        call read_accumulated_field(fio, timepos, timeposm1, met_params%xflux, downward_momentum_flux_units, xflux(:, :),  nr=nr)
       else
         call fi_checkload(fio, met_params%xflux, downward_momentum_flux_units, xflux(:, :), nt=timepos, nr=nr)
       endif
 
       if (met_params%yflux_is_accumulated) then
-        if (timepos == 1) then
-          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, yflux(:, :), nt=timepos, nr=nr)
-        else
-          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-          call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-          yflux(:,:) = tmp2 - tmp1
-        endif
-        ! TODO: Normalise by difference between intervals
-        yflux(:,:) =  yflux / 3600
+        ! Note: Arome files have the wrong units for downward_momentum_flux_units, missing a unit of time
+        call read_accumulated_field(fio, timepos, timeposm1, met_params%yflux, downward_momentum_flux_units, yflux(:, :),  nr=nr)
       else
         call fi_checkload(fio, met_params%yflux, downward_momentum_flux_units, yflux(:, :), nt=timepos, nr=nr)
       endif
+
       surface_stress = hypot(yflux, xflux)
     endif
 
     if (met_params%hflux_is_accumulated) then
-      if (timepos == 1) then
-        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, hflux(:, :), nt=timepos, nr=nr)
-      else
-        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, tmp1(:, :), nt=timeposm1, nr=nr)
-        call fi_checkload(fio, met_params%hflux, accum_surface_heat_flux_units, tmp2(:, :), nt=timepos, nr=nr)
-        hflux(:,:) = tmp2 - tmp1
-      endif
-      ! TODO: Normalise by difference between intervals
-      hflux(:,:) = -hflux / 3600 ! Follow conventions for up/down
+      call read_accumulated_field(fio, timepos, timeposm1, met_params%hflux, accum_surface_heat_flux_units, hflux(:, :), nr=nr)
     else
       call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, hflux(:, :), nt=timepos, nr=nr)
-      hflux(:,:) = -hflux ! Follow conventions for up/down
     endif
-
-
-    deallocate(tmp1)
-    deallocate(tmp2)
+    hflux(:,:) = -hflux ! Follow conventions for up/down
 
     if (met_params%z0 == "") then
-      write(error_unit,*) "WARNING!!!! Surface roughness set to 1 in liu of loaded value."
+      write(error_unit,*) "WARNING!!!! Surface roughness set to 1 in lieu of loaded value."
       z0(:,:) = 1
     else
       call fi_checkload(fio, met_params%z0, surface_roughness_length_units, z0(:, :), nt=timepos, nr=nr)

From c5299a648a94226b3600e26fc9cfe3251b40b0e5 Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Thu, 13 Nov 2025 16:31:10 +0000
Subject: [PATCH 4/6] Move xflux and yflux to local scope

---
 src/common/allocateFields.f90 |  6 +++---
 src/common/fldout_nc.f90      | 18 ++++--------------
 src/common/readfield_fi.f90   | 10 +++++-----
 src/common/snapfldML.f90      |  2 +-
 4 files changed, 13 insertions(+), 23 deletions(-)

diff --git a/src/common/allocateFields.f90 b/src/common/allocateFields.f90
index 599a988f..93b26b9d 100644
--- a/src/common/allocateFields.f90
+++ b/src/common/allocateFields.f90
@@ -28,7 +28,7 @@ module allocateFieldsML
       max_column_scratch, max_column_concentration, &
       aircraft_doserate, aircraft_doserate_scratch, t1_abs, t2_abs, &
       aircraft_doserate_threshold_height, vd_dep, &
-      surface_stress, xflux, yflux, hflux, t2m, z0, &
+      surface_stress, hflux, t2m, z0, &
       ustar, raero, my, &
       total_activity_released, total_activity_lost_domain, total_activity_lost_other, &
       wscav, cloud_cover
@@ -242,7 +242,7 @@ subroutine allocateFields
     if (requires_extra_fields_to_be_read()) then
       allocate(vd_dep(nx,ny,ncomp), STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
-      allocate(surface_stress, xflux, yflux, hflux, t2m, z0, mold=ps2, STAT=AllocateStatus)
+      allocate(surface_stress, hflux, t2m, z0, mold=ps2, STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
       allocate(raero(nx, ny), STAT=AllocateStatus)
       if (AllocateStatus /= 0) ERROR STOP errmsg
@@ -350,7 +350,7 @@ subroutine deAllocateFields
   DEALLOCATE( total_activity_released, total_activity_lost_domain, total_activity_lost_other )
   if (allocated(vd_dep)) then
     deallocate(vd_dep)
-    deallocate(surface_stress, xflux, yflux, hflux, t2m, z0)
+    deallocate(surface_stress, hflux, t2m, z0)
     deallocate(ustar, raero)
   endif
 
diff --git a/src/common/fldout_nc.f90 b/src/common/fldout_nc.f90
index 747f56ed..a6109d6f 100644
--- a/src/common/fldout_nc.f90
+++ b/src/common/fldout_nc.f90
@@ -597,7 +597,7 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
 
   if (output_vd_debug) then
     block
-      use snapfldml, only: t2m, surface_stress, xflux, yflux, z0, hflux, &
+      use snapfldml, only: t2m, surface_stress, z0, hflux, &
         ustar, raero, ps2
       call hres_field(ps2, field_hr1)
       call check(nf90_put_var(iunit, varid%ps_vd, start=ipos, count=isize, values=field_hr1))
@@ -605,10 +605,6 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
       call check(nf90_put_var(iunit, varid%t2m, start=ipos, count=isize, values=field_hr1))
       call hres_field(surface_stress, field_hr1)
       call check(nf90_put_var(iunit, varid%surface_stress, start=ipos, count=isize, values=field_hr1))
-      call hres_field(xflux, field_hr1)
-      call check(nf90_put_var(iunit, varid%xflux, start=ipos, count=isize, values=field_hr1))
-      call hres_field(yflux, field_hr1)
-      call check(nf90_put_var(iunit, varid%yflux, start=ipos, count=isize, values=field_hr1))
       call hres_field(z0, field_hr1)
       call check(nf90_put_var(iunit, varid%z0, start=ipos, count=isize, values=field_hr1))
       call hres_field(hflux, field_hr1)
@@ -1208,11 +1204,8 @@ subroutine initialize_output(filename, itime, ierror)
       block
         use snapmetml, only: downward_momentum_flux_units, surface_heat_flux_units, &
           surface_roughness_length_units, temp_units
-      call nc_declare(iunit, dimids3d, varid%xflux, &
-        "xflux", units=downward_momentum_flux_units, &
-        chunksize=chksz3d)
-      call nc_declare(iunit, dimids3d, varid%yflux, &
-        "yflux", units=downward_momentum_flux_units, &
+      call nc_declare(iunit, dimids3d, varid%surface_stress, &
+        "surface_stress", units=downward_momentum_flux_units, &
         chunksize=chksz3d)
       call nc_declare(iunit, dimids3d, varid%hflux, &
         "hflux", units=surface_heat_flux_units, &
@@ -1423,10 +1416,7 @@ subroutine get_varids(iunit, varid, ierror)
   if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
   ierror = nf90_inq_varid(iunit, "aircraft_doserate_threshold_height", varid%aircraft_doserate_threshold_height)
   if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
-
-  ierror = nf90_inq_varid(iunit, "xflux", varid%xflux)
-  if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
-  ierror = nf90_inq_varid(iunit, "yflux", varid%yflux)
+  ierror = nf90_inq_varid(iunit, "surface_stress", varid%surface_stress)
   if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
   ierror = nf90_inq_varid(iunit, "hflux", varid%hflux)
   if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
diff --git a/src/common/readfield_fi.f90 b/src/common/readfield_fi.f90
index 90fb59ed..2a6dca05 100644
--- a/src/common/readfield_fi.f90
+++ b/src/common/readfield_fi.f90
@@ -896,9 +896,9 @@ subroutine read_accumulated_field(fio, timepos, timeposm1, varname, units, field
     TYPE(FimexIO), intent(inout) :: fio
     character(len=*), intent(in) :: varname, units
     integer, intent(in) :: timepos, timeposm1, nr
-    real, allocatable :: tmp1(:, :), tmp2(:, :)
     real(real32), intent(out) :: field(:, :)
 
+    real, allocatable :: tmp1(:, :), tmp2(:, :)
     allocate(tmp1, tmp2, MOLD=field)
 
     if (timepos == 1) then
@@ -915,15 +915,13 @@ end subroutine read_accumulated_field
 
   subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     USE ieee_arithmetic, only: ieee_is_nan
-    USE iso_fortran_env, only: real64
     USE snapmetML, only: met_params, &
       temp_units, downward_momentum_flux_units, surface_roughness_length_units, &
       accum_surface_heat_flux_units, surface_heat_flux_units
     use drydepml, only: drydep_precompute_meteo, drydep_precompute_particle, &
       requires_extra_fields_to_be_read, classnr
     use snapparML, only: ncomp, run_comp, def_comp
-    ! TODO1: make xflux and yflux temporary variables in readfield
-    use snapfldML, only: ps2, vd_dep, surface_stress, xflux, yflux, hflux, z0, t2m, &
+    use snapfldML, only: ps2, vd_dep, surface_stress, hflux, z0, t2m, &
       ustar, raero, my
     use snaptimers, only: metcalc_timer
 
@@ -933,13 +931,13 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     integer, intent(in) :: nr
     type(datetime_t), intent(in) :: itimefi
 
+    real, allocatable :: xflux(:, :), yflux(:, :)
     integer :: i, mm
 
     if (.not.requires_extra_fields_to_be_read()) then
       return
     endif
 
-    ! TODO: refactor
     if (met_params%surface_stress /= "") then
       ! Load surface_stress
       call fi_checkload(fio, met_params%surface_stress, downward_momentum_flux_units, surface_stress(:, :), nt=timepos, nr=nr)
@@ -947,6 +945,7 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
       ! Either surface_stress or xflux and yflux needs to be defined
       error stop "Assertion error: Either surface_stress or xflux and yflux needs to be defined in the meteorological parameters"
     else
+      allocate(xflux, yflux, MOLD=ps2)
       ! Load and combine surface stress/momentum flux components
       if (met_params%xflux_is_accumulated) then
         ! Note: Arome files have the wrong units for downward_momentum_flux_units, missing a unit of time
@@ -973,6 +972,7 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     hflux(:,:) = -hflux ! Follow conventions for up/down
 
     if (met_params%z0 == "") then
+      ! TODO3: Load z0 from land use data if not defined here
       write(error_unit,*) "WARNING!!!! Surface roughness set to 1 in lieu of loaded value."
       z0(:,:) = 1
     else
diff --git a/src/common/snapfldML.f90 b/src/common/snapfldML.f90
index 9b24bd29..4c5a7fb1 100644
--- a/src/common/snapfldML.f90
+++ b/src/common/snapfldML.f90
@@ -175,6 +175,6 @@ module snapfldML
   !> Deposition velocity on the grid per species
   real, allocatable, save, public :: vd_dep(:, :, :)
 
-  real, allocatable, save, public :: surface_stress(:, :), xflux(:, :), yflux(:, :), hflux(:, :), z0(:, :), t2m(:, :)
+  real, allocatable, save, public :: surface_stress(:, :), hflux(:, :), z0(:, :), t2m(:, :)
   real(real64), allocatable, save, public :: ustar(:,:), raero(:,:), my(:,:)
 end module snapfldML

From 877bdc5f1e187b2bdaf6d55049fb8c4bc07b6759 Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Thu, 13 Nov 2025 16:31:56 +0000
Subject: [PATCH 5/6] Update readfield_nc to match readfield_fi

---
 src/common/readfield_nc.f90 | 86 ++++++++++++++++++++++++-------------
 1 file changed, 57 insertions(+), 29 deletions(-)

diff --git a/src/common/readfield_nc.f90 b/src/common/readfield_nc.f90
index 81066811..c923c12e 100644
--- a/src/common/readfield_nc.f90
+++ b/src/common/readfield_nc.f90
@@ -1172,14 +1172,35 @@ subroutine compute_vertical_coords(alev, blev, ptop)
   vhalf(nk) = vlevel(nk)
 end subroutine
 
+
+  subroutine read_accumulated_field(ncid, varname, timepos, start, startm1, count, field)
+    use iso_fortran_env, only: real32
+    integer, intent(in) :: ncid, timepos, start(7), startm1(7), count(7)
+    character(len=*), intent(in) :: varname
+    real(real32), intent(out) :: field(:, :)
+
+    real, allocatable :: tmp1(:, :), tmp2(:, :)
+    allocate(tmp1, tmp2, MOLD=field)
+
+    if (timepos == 1) then
+      call nfcheckload(ncid, varname, start, count, field(:, :))
+    else
+      call nfcheckload(ncid, varname, start, count, tmp1(:, :))
+      call nfcheckload(ncid, varname, startm1, count, tmp2(:, :))
+      field(:,:) = tmp2 - tmp1
+    endif
+    ! TODO: Normalise by difference between intervals
+    field(:,:) =  field / 3600
+  end subroutine read_accumulated_field
+
   subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
     USE ieee_arithmetic, only: ieee_is_nan
-    USE iso_fortran_env, only: real64
+    USE iso_fortran_env, only: real64, error_unit
     use datetime, only: datetime_t
     use snapmetML, only: met_params
     ! TODO1: make xflux and yflux temporary variables in readfield
     ! TODO2: apply changes in readfield_fi to readfield_nc
-    use snapfldML, only: surface_stress, xflux, yflux, hflux, z0, t2m, vd_dep, ustar, &
+    use snapfldML, only: surface_stress, hflux, z0, t2m, vd_dep, ustar, &
       ps2, raero, my, enspos
     use drydepml, only: classnr, requires_extra_fields_to_be_read, drydep_precompute_meteo, drydep_precompute_particle
     use snapdimML, only: nx, ny
@@ -1187,16 +1208,14 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
     use vgravtablesML, only: vgrav
 
     integer, intent(in) :: ncid
-    integer, intent(in) :: timepos
-    integer, intent(in) :: timeposm1
+    integer, intent(in) :: timepos, timeposm1
     type(datetime_t), intent(in) :: itimefi
 
     integer :: start(7), startm1(7)
     integer :: count(7)
+    real, allocatable :: xflux(:, :), yflux(:, :)
     integer :: i, mm
 
-    real, allocatable :: tmp1(:, :), tmp2(:, :)
-
     if (.not.requires_extra_fields_to_be_read()) then
       return
     endif
@@ -1206,36 +1225,45 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
     call calc_2d_start_length(startm1, count, nx, ny, -1, &
           enspos, timeposm1, met_params%has_dummy_dim)
 
-    allocate(tmp1(nx,ny), tmp2(nx,ny))
 
-    ! Fluxes are integrated: Deaccumulate
-    if (timepos == 1) then
-      call nfcheckload(ncid, met_params%xflux, start, count, xflux(:,:))
-      call nfcheckload(ncid, met_params%yflux, start, count, yflux(:,:))
+    if (met_params%surface_stress /= "") then
+      ! Load surface_stress
+      call nfcheckload(ncid, met_params%surface_stress, start, count, surface_stress(:, :))
+    else if (met_params%xflux == "" .OR. met_params%yflux == "") then
+      ! Either surface_stress or xflux and yflux needs to be defined
+      error stop "Assertion error: Either surface_stress or xflux and yflux needs to be defined in the meteorological parameters"
     else
-      call nfcheckload(ncid, met_params%xflux, start, count, tmp1(:,:))
-      call nfcheckload(ncid, met_params%xflux, startm1, count, tmp2(:,:))
-      xflux(:,:) = tmp2 - tmp1
-      call nfcheckload(ncid, met_params%yflux, start, count, tmp1(:,:))
-      call nfcheckload(ncid, met_params%yflux, startm1, count, tmp2(:,:))
-      yflux(:,:) = tmp2 - tmp1
+      allocate(xflux, yflux, MOLD=ps2)
+      ! Load and combine surface stress/momentum flux components
+      if (met_params%xflux_is_accumulated) then
+        call read_accumulated_field(ncid, met_params%xflux, timepos, start, startm1, count, xflux(:,:))
+      else
+        call nfcheckload(ncid, met_params%xflux, start, count, xflux(:,:))
+      endif
+
+      if (met_params%yflux_is_accumulated) then
+        call read_accumulated_field(ncid, met_params%yflux, timepos, start, startm1, count, yflux(:,:))
+      else
+        call nfcheckload(ncid, met_params%yflux, start, count, yflux(:,:))
+      endif
+
+      surface_stress = hypot(yflux, xflux)
     endif
-    ! TODO: Normalise by difference between intervals
-    xflux(:,:) =  xflux / 3600
-    yflux(:,:) =  yflux / 3600
-    surface_stress = hypot(yflux, xflux)
 
-    if (timepos == 1) then
-      call nfcheckload(ncid, met_params%hflux, start, count, hflux(:,:))
+    if (met_params%hflux_is_accumulated) then
+      call read_accumulated_field(ncid, met_params%hflux, timepos, start, startm1, count, hflux(:,:))
     else
-      call nfcheckload(ncid, met_params%hflux, start, count, tmp1(:,:))
-      call nfcheckload(ncid, met_params%hflux, start, count, tmp2(:,:))
-      hflux(:,:) = tmp2 - tmp1
+      call nfcheckload(ncid, met_params%hflux, start, count, hflux(:,:))
     endif
-    ! TODO: Normalise by difference between intervals
-    hflux(:,:) = -hflux / 3600 ! Follow conventions for up/down
+    hflux(:,:) = -hflux ! Follow conventions for up/down
 
-    call nfcheckload(ncid, met_params%z0, start, count, z0(:, :))
+    if (met_params%z0 == "") then
+      ! TODO3: Load z0 from land use data if not defined here
+      write(error_unit,*) "WARNING!!!! Surface roughness set to 1 in lieu of loaded value."
+      z0(:,:) = 1
+    else
+      call nfcheckload(ncid, met_params%z0, start, count, z0(:, :))
+    endif
 
     call nfcheckload(ncid, met_params%t2m, start, count, t2m(:, :))
 

From c092105a5cd769cf6b3a28b517ab60640fa4ee01 Mon Sep 17 00:00:00 2001
From: Halvard Sutterud <halvard.sutterud@gmail.com>
Date: Fri, 21 Nov 2025 10:31:22 +0000
Subject: [PATCH 6/6] Combine xflux and yflux in drydep test

---
 src/test/testDryDep.f90 | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/test/testDryDep.f90 b/src/test/testDryDep.f90
index c76731d0..c48baf91 100644
--- a/src/test/testDryDep.f90
+++ b/src/test/testDryDep.f90
@@ -13,7 +13,7 @@ program testDryDep
   type(datetime_t) :: itimefi
   integer :: i
   integer(kind=1) :: classnr
-  real(4) :: ps2, t2m, yflux, xflux, z0, hflux, &
+  real(4) :: ps2, t2m, yflux, xflux, surface_stress, z0, hflux, &
         vd_dep
   real(kind=real64) :: ustar, raero, my
 
@@ -55,10 +55,11 @@ program testDryDep
   t2m = 280.0 ! K
   yflux = 1.0 ! N /m2/hr
   xflux = 1.0 ! N /m2/hr
+  surface_stress = HYPOT(yflux, xflux)
   z0 = 0.1 ! m
   hflux = 100.0 ! W hr/m2
   classnr = 21 ! class number 11=sea, 21=mixed forest
-  call drydep_precompute_meteo(ps2*100., t2m, yflux, xflux, z0, hflux, &
+  call drydep_precompute_meteo(ps2*100., t2m, surface_stress, z0, hflux, &
     ustar, raero, my)
   ! Test dry-dep velocity for Cs137
   call drydep_precompute_particle(ps2*100., t2m, &