Feat: emerson drydep for ecemep

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, 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, 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, hflux, t2m, z0)
     deallocate(ustar, raero)
   endif
 

src/common/drydep.f90

@@ -77,24 +77,23 @@ 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
 
   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

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,16 +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, 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(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(surface_stress, field_hr1)
+      call check(nf90_put_var(iunit, varid%surface_stress, 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)
@@ -1205,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, &
@@ -1420,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

src/common/readfield_fi.f90

@@ -893,19 +893,40 @@ 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(real32), intent(out) :: field(:, :)
+
+    real, allocatable :: tmp1(:, :), tmp2(:, :)
+    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
     USE snapmetML, only: met_params, &
       temp_units, downward_momentum_flux_units, surface_roughness_length_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 ftestML, only: ftest
     use snapdebug, only: idebug
 
     use snapparML, only: ncomp, run_comp, def_comp
-    use snapfldML, only: ps2, vd_dep, xflux, yflux, hflux, z0, t2m, &
+    use snapfldML, only: ps2, vd_dep, surface_stress, hflux, z0, t2m, &
       ustar, raero, my
     use snaptimers, only: metcalc_timer
 
@@ -915,50 +936,58 @@ subroutine read_drydep_required_fields(fio, timepos, timeposm1, nr, itimefi)
     integer, intent(in) :: nr
     type(datetime_t), intent(in) :: itimefi
 
-    real, allocatable :: tmp1(:, :), tmp2(:, :)
+    real, allocatable :: xflux(:, :), yflux(:, :)
     integer :: i, mm
 
     if (.not.requires_extra_fields_to_be_read()) then
       return
     endif
 
-    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)
+    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
+      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
+        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
 
-      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
+        ! 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
-    ! 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
+      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, 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)
     endif
-    ! TODO: Normalise by difference between intervals
-    hflux(:,:) = -hflux / 3600 ! Follow conventions for up/down
+    hflux(:,:) = -hflux ! Follow conventions for up/down
 
-
-    call fi_checkload(fio, met_params%z0, surface_roughness_length_units, z0(:, :), nt=timepos, nr=nr)
+    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 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

src/common/readfield_nc.f90

@@ -1174,12 +1174,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
-    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, hflux, z0, t2m, vd_dep, ustar, &
       ps2, raero, my, enspos
     use drydepml, only: drydep_precompute_meteo, drydep_precompute_particle, &
       requires_extra_fields_to_be_read, classnr
@@ -1191,16 +1214,14 @@ subroutine read_drydep_required_fields(ncid, timepos, timeposm1, itimefi)
 
 
     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
@@ -1210,39 +1231,49 @@ 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
 
-    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(:, :))
 
-    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

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(:, :), hflux(:, :), z0(:, :), t2m(:, :)
   real(real64), allocatable, save, public :: ustar(:,:), raero(:,:), my(:,:)
 end module snapfldML