diff --git a/.gitignore b/.gitignore
index a9be715e..f5caf08b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -82,36 +82,18 @@ Snappy.egg-info
/src/naccident/snap_batch_copy.f90
/src/naccident/snap_batch_copy.F90
-/src/nbomb/*.inc
-/src/nbomb/*.mod
-/src/nbomb/bsnap_nbomb
-/src/nbomb/create_bomb_input
-/src/nbomb/snap_batch_copy.f90
-/src/nbomb/snap_batch_copy.F90
-
/src/test/testDateCalc
/src/test/test_array_utils
/src/test/find_parameters_fi_test
/src/test/snap_batch_copy.F90
-/src/traj/*.inc
-/src/traj/*.mod
-/src/traj/bsnap_traj
-/src/traj/create_traj_input
-/src/traj/snap_batch_copy.f90
-/src/traj/snap_batch_copy.F90
-
-/src/volcano/*.mod
-/src/volcano/*.inc
-/src/volcano/snapargos.inc
-/src/volcano/bsnap_volcano
-/src/volcano/snap_batch_copy.f90
-/src/volcano/snap_batch_copy.F90
src/test/snap_testdata
-
docs/html
docs/latex
bsnap_*
+
+utils/SnapPy/Snappy.egg-info
+src/test/snap_testdata
diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 00000000..60cd4325
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,37 @@
+FROM ubuntu:22.04 AS builder
+
+ENV SNAP_FIMEX_VERSION=1.9
+
+RUN apt-get update -y && \
+ apt-get install -y software-properties-common && \
+ add-apt-repository -y ppa:met-norway/fimex && \
+ apt-get update -y && \
+ apt-get install -y libnetcdff-dev gfortran libfimex-$SNAP_FIMEX_VERSION-dev make && \
+ apt-get remove -y software-properties-common && \
+ apt-get -y autoremove && \
+ rm -rf /var/lib/apt/lists/*
+
+WORKDIR /snap
+ARG VERSION="latest"
+ADD src .
+RUN ln --symbolic --force gcc_pkgconfig.mk current.mk
+ENV SNAP_USE_OMP=1
+RUN make clean && make
+
+FROM ubuntu:22.04
+
+ENV SNAP_FIMEX_VERSION=1.9
+
+RUN apt-get update -y && \
+ apt-get install -y software-properties-common && \
+ add-apt-repository -y ppa:met-norway/fimex && \
+ apt-get update -y && \
+ apt-get install -y libnetcdff7 libfimex-$SNAP_FIMEX_VERSION-0 tini && \
+ apt-get remove -y software-properties-common && \
+ apt-get -y autoremove && \
+ rm -rf /var/lib/apt/lists/*
+
+WORKDIR /snap
+COPY --from=builder /snap/naccident/bsnap_naccident /snap/bsnap
+
+ENTRYPOINT ["/usr/bin/env", "TINI_VERBOSITY=0", "/usr/bin/tini", "-g", "--", "/snap/bsnap"]
diff --git a/src/common/allocateFields.f90 b/src/common/allocateFields.f90
index 906916e9..11d09ad7 100644
--- a/src/common/allocateFields.f90
+++ b/src/common/allocateFields.f90
@@ -24,10 +24,14 @@ module allocateFieldsML
depdry, depwet, accprec, avgprec, avghbl, precip, &
pmsl1, pmsl2, field1, field2, field3, field4, field3d1, xm, ym, &
garea, field_hr1, field_hr2, field_hr3, hbl_hr, &
+ precip3d, cw3d, &
max_column_scratch, max_column_concentration, &
aircraft_doserate, aircraft_doserate_scratch, t1_abs, t2_abs, &
- aircraft_doserate_threshold_height, &
- total_activity_released, total_activity_lost_domain, total_activity_lost_other
+ aircraft_doserate_threshold_height, vd_dep, &
+ xflux, yflux, hflux, t2m, z0, leaf_area_index, &
+ roa, ustar, monin_l, raero, vs, rs, &
+ total_activity_released, total_activity_lost_domain, total_activity_lost_other, &
+ wscav, cloud_cover
USE snapfilML, only: idata, fdata
USE snapgrdML, only: ahalf, bhalf, vhalf, alevel, blevel, vlevel, imodlevel, &
compute_column_max_conc, compute_aircraft_doserate, aircraft_doserate_threshold
@@ -44,6 +48,7 @@ subroutine allocateFields
USE snapdimML, only: nx, ny, nk, output_resolution_factor, ldata, maxsiz
USE snapparML, only: ncomp, nocomp, iparnum
USE releaseML, only: mplume, iplume, plume_release, mpart
+ USE snapmetML, only: met_params
logical, save :: FirstCall = .TRUE.
integer :: AllocateStatus
@@ -147,6 +152,14 @@ subroutine allocateFields
ALLOCATE ( precip(nx,ny), STAT = AllocateStatus)
IF (AllocateStatus /= 0) ERROR STOP errmsg
+ if (met_params%use_3d_precip) then
+ ALLOCATE(precip3d(nx,ny,nk), cw3d(nx,ny,nk), STAT=AllocateStatus)
+ if (AllocateStatus /= 0) ERROR STOP errmsg
+ ALLOCATE(wscav(nx,ny,nk,ncomp),STAT=AllocateStatus)
+ if (AllocateStatus /= 0) ERROR STOP errmsg
+ wscav(:,:,:,:) = 0.0
+ allocate(cloud_cover(nx,ny,nk))
+ endif
! the calculation-fields
ALLOCATE ( avghbl(nx,ny), STAT = AllocateStatus)
@@ -223,6 +236,18 @@ subroutine allocateFields
total_activity_lost_domain(:) = 0.0
total_activity_lost_other(:) = 0.0
+ ! Dry deposition fields
+ block
+ use drydepml, only: requires_extra_fields_to_be_read
+ 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, leaf_area_index, mold=ps2)
+ allocate(roa(nx, ny))
+ allocate(ustar, monin_l, raero, vs, rs, mold=roa)
+ endif
+ end block
+
end subroutine allocateFields
@@ -276,6 +301,10 @@ subroutine deAllocateFields
DEALLOCATE ( pmsl2)
DEALLOCATE ( precip)
+ if (allocated(precip3d)) deallocate(precip3d)
+ if (allocated(cw3d)) deallocate(cw3d)
+ if (allocated(wscav)) deallocate(wscav)
+ if (allocated(cloud_cover)) deallocate(cloud_cover)
DEALLOCATE ( avghbl )
DEALLOCATE ( avgprec )
@@ -316,6 +345,11 @@ subroutine deAllocateFields
DEALLOCATE ( plume_release )
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, leaf_area_index)
+ deallocate(roa, ustar, monin_l, raero, vs, rs)
+ endif
end subroutine deAllocateFields
end module allocateFieldsML
diff --git a/src/common/drydep.f90 b/src/common/drydep.f90
index a7182ea4..7be401b9 100644
--- a/src/common/drydep.f90
+++ b/src/common/drydep.f90
@@ -15,14 +15,132 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
-module drydep
+module drydepml
+ use ISO_FORTRAN_ENV, only: real64, int8
implicit none
private
- public drydep1, drydep2
+ public :: drydep, gravitational_settling, preprocess_landfraction, unload, &
+ requires_extra_fields_to_be_read, drydep_precompute
+
+ integer, parameter, public :: DRYDEP_SCHEME_UNDEFINED = 0
+ integer, parameter, public :: DRYDEP_SCHEME_OLD = 1
+ integer, parameter, public :: DRYDEP_SCHEME_NEW = 2
+ integer, parameter, public :: DRYDEP_SCHEME_EMEP = 3
+ integer, parameter, public :: DRYDEP_SCHEME_ZHANG = 4
+ integer, parameter, public :: DRYDEP_SCHEME_EMERSON = 5
+
+ !> The active dry deposition scheme
+ integer, save, public :: drydep_scheme = DRYDEP_SCHEME_UNDEFINED
+
+
+ real(real64), parameter :: R = 287.05
+ real(real64), parameter :: grav = 9.8
+ real(real64), parameter :: CP = 1005.0
+ real(real64), parameter :: pi = 2.0*asin(1.0)
+
+ !> Kinematic viscosity of air, m2 s-1 at +15 C
+ real(real64), parameter :: ny = 1.5e-5
+ !> Mean free path of air molecules [m]
+ real(real64), parameter :: lambda = 0.065e-6
+
+ !> Boltzmanns constant, J K-1
+ real(real64), parameter :: bolzc = 1.380649e-23
+
+ character(len=256), save, public :: largest_landfraction_file = "not set"
+
+ integer(int8), save, public, allocatable :: classnr(:, :)
contains
+subroutine drydep(tstep, part)
+ use snapfldML, only: vd_dep
+ use particleML, only: Particle
+ real, intent(in) :: tstep
+ type(particle), intent(inout) :: part
+
+
+ if (drydep_scheme == DRYDEP_SCHEME_OLD) call drydep1(part)
+ if (drydep_scheme == DRYDEP_SCHEME_NEW) call drydep2(tstep, part)
+ if (drydep_scheme == DRYDEP_SCHEME_EMEP .or. &
+ drydep_scheme == DRYDEP_SCHEME_ZHANG .or. &
+ drydep_scheme == DRYDEP_SCHEME_EMERSON) call drydep_nonconstant_vd(tstep, vd_dep, part)
+end subroutine
+
+pure logical function requires_extra_fields_to_be_read()
+ requires_extra_fields_to_be_read = ( &
+ (drydep_scheme == DRYDEP_SCHEME_ZHANG).or. &
+ (drydep_scheme == DRYDEP_SCHEME_EMERSON).or. &
+ (drydep_scheme == DRYDEP_SCHEME_EMEP))
+end function
+
+!> Precompute dry deposition constants for (x, y)
+!> At every time step dry deposition is computed
+!> by lookup into the vd_dep matrix
+pure subroutine drydep_precompute(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, leaf_area_index, diam, density, classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ use iso_fortran_env, only: real64, int8
+ real, intent(in) :: surface_pressure(:,:) !> [hPa]
+ real, intent(in) :: t2m(:,:)
+ real, intent(in) :: yflux(:,:), xflux(:,:)
+ real, intent(in) :: z0(:,:), hflux(:,:)
+ real, intent(in) :: leaf_area_index(:,:)
+ real, intent(in) :: diam
+ real, intent(in) :: density
+ integer(int8), intent(in) :: classnr(:,:)
+ real, intent(out) :: vd_dep(:,:)
+ real(real64), intent(out) :: roa(:,:), monin_obukhov_length(:,:), raero(:,:), vs(:,:), ustar(:,:), rs(:,:)
+
+ select case(drydep_scheme)
+ case (DRYDEP_SCHEME_EMEP)
+ call drydep_emep_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, leaf_area_index, diam, density, classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ case (DRYDEP_SCHEME_ZHANG)
+ call drydep_zhang_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, real(diam, kind=real64), real(density, kind=real64), classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ case (DRYDEP_SCHEME_EMERSON)
+ call drydep_emerson_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, real(diam, kind=real64), real(density, kind=real64), classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ case default
+ error stop "Precomputation should not be called for this dry deposition scheme"
+ end select
+end subroutine
+
+!> Stoer landfraction values
+subroutine preprocess_landfraction(values)
+ use iso_fortran_env, only: real32, error_unit
+ real(real32), intent(in) :: values(:,:)
+
+ write(error_unit,*) "We do not currently check the landclasses programatically"
+ write(error_unit,*) "The classes must be:"
+ write(error_unit,*) " 11: Sea"
+ write(error_unit,*) " 12: Inland water"
+ write(error_unit,*) " 13: Tundra/desert"
+ write(error_unit,*) " 14: Ice and ice sheets"
+ write(error_unit,*) " 15: Urban"
+ write(error_unit,*) " 16: Crops"
+ write(error_unit,*) " 17: Grass"
+ write(error_unit,*) " 18: Wetlands"
+ write(error_unit,*) " 19: Evergreen needleleaf"
+ write(error_unit,*) " 20: Deciduous broadleaf"
+ write(error_unit,*) " 21: Mixed forest"
+ write(error_unit,*) " 22: Shrubs and interrupted woodlands"
+
+ if (allocated(classnr)) then
+ error stop "preprocess_landfraction is to be called once only"
+ endif
+ allocate(classnr(size(values,dim=1),size(values,dim=2)))
+ classnr(:,:) = nint(values)
+end subroutine
+
+subroutine unload()
+ if (allocated(classnr)) deallocate(classnr)
+end subroutine
+
!> Purpose: Compute dry deposition for each particle and each component
!> and store depositions in nearest gridpoint in a field
!>
@@ -76,7 +194,9 @@ subroutine drydep2(tstep, part)
integer :: m,i,j,mm
real :: deprate, dep
- real, parameter :: h = 30.0
+ real, parameter :: h = 30.0 ! [m]
+ real, parameter :: vd_gas = 0.008 ! [m/s]
+ real, parameter :: vd_particles = 0.002 ! [m/s]
m = part%icomp
!#### 30m = surface-layer (deposition-layer); sigma(hybrid)=0.996 ~ 30m
@@ -85,13 +205,13 @@ subroutine drydep2(tstep, part)
if (def_comp(m)%radiusmym <= 0.05) then
! gas
- deprate = 1.0 - exp(-tstep*(0.008)/h)
+ deprate = 1.0 - exp(-tstep*(vd_gas)/h)
else if (def_comp(m)%radiusmym <= 10.0) then
! particle 0.05=10
- deprate = 1.0 - exp(-tstep*(0.002+part%grv)/h)
+ ! particle r>=10.0
+ deprate = 1.0 - exp(-tstep*(vd_particles + part%grv)/h)
! complete deposition when particle hits ground
if (part%z == vlevel(1)) deprate = 1.
endif
@@ -103,4 +223,312 @@ subroutine drydep2(tstep, part)
depdry(i,j,mm) = depdry(i,j,mm) + dble(dep)
end if
end subroutine drydep2
-end module drydep
+
+pure function aerodynres(L, ustar, z0) result(raero)
+ real(real64), intent(in) :: L, ustar, z0
+ real(real64) :: raero
+
+ real(real64), parameter :: ka = 0.4
+ real(real64), parameter :: z = 30 ! Assumed height of surface/constant flux layer
+
+ real(real64) :: fac, fi
+
+ if (L < 0) then
+ fac = 0.598 + 0.390 * log(-z / L) - 0.09 * (log(-z / L)) ** 2
+ fi = exp(fac)
+ else if (L > 0) then
+ fi = -5 * (z/L)
+ else
+ fi = 0
+ endif
+
+ raero = (1 / (ka * ustar)) * (log(z/z0) - fi)
+end function
+
+pure elemental function gravitational_settling(roa, diam, ro_part) result(vs)
+ real(real64), intent(in) :: roa
+ !> Diameter in m
+ real(real64), intent(in) :: diam
+ !> Density in kg m-3
+ real(real64), intent(in) :: ro_part
+
+ real(real64) :: vs
+
+ real(real64) :: my ! Dynamic visocity of air [kg/(m s)]
+
+ real(real64) :: fac1, cslip
+
+ my = ny * roa
+ fac1 = -0.55 * diam / lambda
+ cslip = 1 + 2 * lambda / diam * ( 1.257 + 0.4 * exp(fac1) )
+
+ vs = ro_part * diam * diam * grav * cslip / (18*my)
+end function
+
+!> Dry deposition velocity given by
+!> Simpson et al. 2012, The EMEP MSC-W chemical transport model - technical description
+!> https://doi.org/10.5194/acp-12-7825-2012
+pure elemental subroutine drydep_emep_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, leaf_area_index, diam, density, classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ !> In hPa
+ real, intent(in) :: surface_pressure
+ real, intent(in) :: t2m
+ real, intent(in) :: yflux, xflux
+ real, intent(in) :: z0, hflux
+ real, intent(in) :: leaf_area_index
+ real, intent(in) :: diam
+ real, intent(in) :: density
+ integer(int8), intent(in) :: classnr
+ real, intent(out) :: vd_dep
+
+ real(real64), intent(out) :: roa
+ real(real64), intent(out) :: ustar
+ real(real64), intent(out) :: monin_obukhov_length
+ real(real64) :: SAI
+
+ real(real64), parameter :: k = 0.4
+ ! real, parameter :: a1 = 0.002
+ real(real64) :: a1
+ real(real64), parameter :: a2 = 300
+ real(real64) :: a1sai
+ real(real64), intent(out) :: vs
+ real(real64), intent(out) :: raero
+ real(real64), intent(out) :: rs
+ real(real64) :: fac
+
+
+ roa = surface_pressure / (t2m * R)
+ vs = real(gravitational_settling(real(roa, real64), real(diam, real64), real(density, real64)))
+
+ ustar = hypot(yflux, xflux) / sqrt(roa)
+ monin_obukhov_length = - roa * CP * t2m * (ustar**3) / (k * grav * hflux)
+
+ SAI = leaf_area_index + 1
+ if (classnr >= 19 .and. classnr <= 21) then
+ a1sai = 0.008 * SAI / 10
+ else
+ a1sai = 0.0
+ endif
+
+ a1 = max(a1sai, 0.002)
+
+ raero = aerodynres(real(monin_obukhov_length, real64), real(ustar, real64), real(z0, real64))
+
+ if (monin_obukhov_length > -25.0 .and. monin_obukhov_length < 0.0) then
+ monin_obukhov_length = -25.0
+ endif
+ if (monin_obukhov_length > 0) then
+ rs = 1.0 / (ustar * a1)
+ else
+ fac = (-a2 / monin_obukhov_length) ** ( 2.0 / 3.0 )
+ rs = 1.0 / (ustar * a1 * (1 + fac))
+ endif
+
+ vd_dep = 1.0 / (rs + raero) + vs
+end subroutine
+
+pure real(kind=real64) function lookup_A(classnr)
+ use ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
+ integer(int8), intent(in) :: classnr
+
+ select case(classnr)
+ ! Sea, Inland_water, Tundra_and_desert, Ice_and_ice_sheets
+ case(11,12,13,14)
+ lookup_A = IEEE_VALUE(lookup_A, IEEE_QUIET_NAN)
+ ! Urban, Wetlands, Shurbs_and_interrupted_woodlands
+ case(15,18,22)
+ lookup_A = 10.0
+ ! Crops, Grass, Evergreen_needleleaf
+ case(16,17,19)
+ lookup_A = 2.0
+ ! Decidious_broadleaf, Mixed_forest
+ case(20,21)
+ lookup_A = 5.0
+ case default
+ lookup_A = IEEE_VALUE(lookup_A, IEEE_QUIET_NAN)
+ end select
+
+end function
+
+!> Dry deposition velocites based on
+!> Zhang et al. 2001, A size-segregated particle dry deposition scheme for an atmospheric aerosol module
+!> https://doi.org/10.1016/S1352-2310(00)00326-5
+pure elemental subroutine drydep_zhang_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, diam, density, classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ use ieee_arithmetic, only: ieee_is_nan
+ !> In hPa
+ real, intent(in) :: surface_pressure
+ real, intent(in) :: t2m
+ real, intent(in) :: yflux, xflux
+ real, intent(in) :: z0, hflux
+ real(real64), intent(in) :: diam
+ real(real64), intent(in) :: density
+ integer(int8), intent(in) :: classnr
+ real, intent(out) :: vd_dep
+ real(real64), intent(out) :: roa, monin_obukhov_length, raero, vs, ustar, rs
+
+ !> Aerial factor for interception (table 3 Zhang et al. (2001)), corresponding to evergreen
+ !> needleleaf trees (i.e. close to maximum deposition)
+ real(real64) :: A
+ real(real64), parameter :: k = 0.4
+
+ real(real64) :: fac1, cslip, bdiff, my, sc, EB, EIM, EIN, stokes
+ real(real64) :: Apar
+
+ roa = surface_pressure / (t2m * R)
+ vs = gravitational_settling(roa, diam, density)
+
+ ustar = hypot(yflux, xflux) / 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
+
+ fac1 = -0.55 * diam / lambda
+ ! Cunningham slip factor
+ cslip = 1 + 2 * lambda / diam * (1.257 + 0.4*exp(fac1))
+
+ ! Brownian diffusion
+ ! Brownian diffusivity of air (see equation 19 of Giardiana and Buffa, 2018)
+ ! bdiff=2.83e-11 # Browian diffusion coefficient for 1 um particle (see Brereton, 2014)
+ bdiff = bolzc * t2m * cslip / (3 * pi * my * diam)
+
+ sc = ny / bdiff
+ ! A range og 0.5-0.58 dependening on the surface is given, 0.54=grass
+ EB = sc ** (-0.54)
+
+ Apar = lookup_A(classnr)
+ A = Apar * 1e-3
+ ! Impaction
+ if (.not. ieee_is_nan(A)) then
+ ! Stokes number for vegetated surfaces (Zhang (2001)
+ stokes = vs * ustar / (grav * A)
+ else
+ stokes = vs * ustar * ustar / (grav * ny)
+ endif
+
+ EIM = (stokes / (0.8 + stokes)) ** 2
+
+ ! Interception
+ if (ieee_is_nan(A)) then
+ ! No interception over water surfaces
+ EIN = 0.0
+ else
+ EIN = 0.5 * (diam / A) ** 2
+ endif
+
+ rs = 1.0 / (3.0 * ustar * (EB + EIM + EIN))
+
+ vd_dep = 1.0 / (raero + rs) + vs
+end subroutine
+
+!> Dry deposition velocites based on
+!> Emerson et al. 2020, Revisiting particle dry deposition and its role in radiative effect estimates
+!> https://doi.org/10.1073/pnas.2014761117
+pure elemental subroutine drydep_emerson_vd(surface_pressure, t2m, yflux, xflux, z0, &
+ hflux, diam, density, classnr, vd_dep, &
+ roa, ustar, monin_obukhov_length, raero, vs, rs)
+ use ieee_arithmetic, only: ieee_is_nan
+ !> In hPa
+ real, intent(in) :: surface_pressure
+ real, intent(in) :: t2m
+ real, intent(in) :: yflux, xflux
+ real, intent(in) :: z0, hflux
+ real(real64), intent(in) :: diam
+ real(real64), intent(in) :: density
+ integer(int8), intent(in) :: classnr
+ real, intent(out) :: vd_dep
+ real(real64), intent(out) :: roa, monin_obukhov_length, raero, vs, ustar, rs
+
+ !> Aerial factor for interception (table 3 Zhang et al. (2001)), corresponding to evergreen
+ !> needleleaf trees (i.e. close to maximum deposition)
+ real(real64) :: A
+ real(real64), parameter :: k = 0.4
+
+ real(real64) :: fac1, cslip, bdiff, my, sc, EB, EIM, EIN, stokes
+ real(real64) :: Apar
+
+ roa = surface_pressure / (t2m * R)
+ vs = gravitational_settling(roa, diam, density)
+
+ ustar = hypot(yflux, xflux) / 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
+
+ fac1 = -0.55 * diam / lambda
+ ! Cunningham slip factor
+ cslip = 1 + 2 * lambda / diam * (1.257 + 0.4*exp(fac1))
+
+ ! Brownian diffusion
+ ! Brownian diffusivity of air (see equation 19 of Giardiana and Buffa, 2018)
+ ! bdiff=2.83e-11 # Browian diffusion coefficient for 1 um particle (see Brereton, 2014)
+ bdiff = bolzc * t2m * cslip / (3 * pi * my * diam)
+
+ sc = ny / bdiff
+ EB = 0.2 * sc ** (-2.0 / 3.0)
+
+ Apar = lookup_A(classnr)
+ A = Apar * 1e-3
+ ! Impaction
+ if (.not. ieee_is_nan(A)) then
+ ! Stokes number for vegetated surfaces (Zhang (2001)
+ stokes = vs * ustar / (grav * A)
+ else
+ stokes = vs * ustar * ustar / (grav * ny)
+ endif
+ EIM = 0.4 * (stokes / (0.8 + stokes)) ** 1.7
+
+ ! Interception
+ if (ieee_is_nan(A)) then
+ ! No interception over water surfaces
+ EIN = 0.0
+ else
+ EIN = 2.5 * (diam / A) ** 0.8
+ endif
+
+ rs = 1.0 / (3.0 * ustar * (EB + EIM + EIN))
+ vd_dep = 1.0 / (raero + rs) + vs
+end subroutine
+
+!> Apply precomputed dry deposition velocities
+subroutine drydep_nonconstant_vd(tstep, vd, part)
+ use particleML, only: Particle
+ use snapfldML, only: depdry
+ use snapparML, only: def_comp
+ use snapdimML, only: hres_pos
+
+ !> Timestep of the simulation in seconds
+ real, intent(in) :: tstep
+ !> Particle which observes deposition
+ type(Particle), intent(inout) :: part
+ !> Deposition velocity
+ real, intent(in) :: vd(:,:, :)
+
+ real, parameter :: h = 30.0
+
+ integer :: m, mm, i, j
+ real :: dep, deprate_m1
+
+ m = part%icomp
+ if (def_comp(m)%kdrydep == 1 .and. part%z > 0.996) then
+
+ mm = def_comp(m)%to_running
+
+ i = nint(part%x)
+ j = nint(part%y)
+
+ deprate_m1 = exp(-tstep*vd(i,j,mm)/h)
+
+ i = hres_pos(part%x)
+ j = hres_pos(part%y)
+ dep = part%scale_rad(deprate_m1)
+ !$OMP atomic
+ depdry(i,j,mm) = depdry(i,j,mm) + dble(dep)
+ end if
+end subroutine
+
+end module drydepml
diff --git a/src/common/fldout_nc.f90 b/src/common/fldout_nc.f90
index d35e16e3..6809650b 100644
--- a/src/common/fldout_nc.f90
+++ b/src/common/fldout_nc.f90
@@ -51,11 +51,16 @@ module fldout_ncML
integer :: ic
integer :: icml
integer :: conc_column = -1
+ !> Dry deposition velocity
+ integer :: vd = -1
+ !> Wet scavenging rate
+ integer :: wetscavrate = -1
end type
!> Variables in a file
type :: common_var
integer :: accum_prc
+ integer :: instant_prc = -1
integer :: prc
integer :: mslp
integer :: icblt
@@ -77,6 +82,21 @@ module fldout_ncML
integer :: aircraft_doserate_threshold_height
integer :: components
type(component_var) :: comp(mcomp)
+ integer :: xflux = -1
+ integer :: yflux = -1
+ integer :: hflux = -1
+ integer :: z0 = -1
+ integer :: t2m = -1
+ integer :: lai = -1
+ integer :: roa = -1
+ integer :: ustar = -1
+ integer :: monin_l = -1
+ integer :: raero = -1
+ integer :: vs = -1
+ integer :: rs = -1
+ integer :: ps_vd = -1
+ integer :: landfraction = -1
+ integer :: wetdeprate = -1
end type
!> dimensions used in a file
@@ -100,25 +120,28 @@ module fldout_ncML
character(len=256), save, public :: massbalance_filename = ""
integer, save, allocatable :: massbalance_file
+ logical, save, public :: output_wetdeprate = .false.
+
contains
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, modlevel_is_average, imslp, precipitation_in_output, &
+ USE snapgrdML, only: imodlevel, 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, &
- field1, &
+ aircraft_doserate_threshold, output_column, &
+ output_column, output_vd, output_vd_debug
+ USE snapfldML, only: field_hr1, field_hr2, field_hr3, hbl_hr
+ USE snapfldML, only: field1, &
depdry, depwet, &
avgbq1, avgbq2, garea, pmsl1, pmsl2, hbl1, hbl2, &
accdry, accwet, avgprec, concen, ps1, ps2, avghbl, &
concacc, accprec, max_column_concentration, aircraft_doserate, &
aircraft_doserate_threshold_height, &
- total_activity_released, total_activity_lost_domain, total_activity_lost_other
- USE snapparML, only: time_profile, nocomp, output_component, def_comp, &
- TIME_PROFILE_BOMB
+ total_activity_released, total_activity_lost_domain, total_activity_lost_other, &
+ vd_dep, precip3d
+ USE snapparML, only: output_component, def_comp, nocomp
USE snapdebug, only: iulog, idebug
USE ftestML, only: ftest
USE snapdimML, only: nx, ny, output_resolution_factor, hres_field, hres_pos
@@ -141,7 +164,6 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
integer :: nptot1,nptot2
real(real64) :: bqtot1,bqtot2
- real(real64) :: dblscale
integer :: i,j,m,n
integer(int16), pointer, dimension(:) :: mm
@@ -200,6 +222,12 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
call hres_field(field1, field_hr1)
call check(nf90_put_var(iunit, varid%accum_prc, start=ipos, count=isize, &
values=field_hr1), "set_accum_prc")
+
+ if (allocated(precip3d)) then
+ field1(:,:) = precip3d(:,:,2)
+ call check(nf90_put_var(iunit, varid%instant_prc, start=ipos, count=isize, &
+ values=field1), "set_accum_prc")
+ endif
end if
!..mslp (if switched on)
@@ -321,6 +349,12 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
call check(nf90_put_var(iunit, varid%comp(m)%idd, start=ipos, count=isize, &
values=field_hr1), "set_idd(m)")
+
+ if (output_vd) then
+ call hres_field(vd_dep(:,:,m), field_hr2)
+ call check(nf90_put_var(iunit, varid%comp(m)%vd, start=ipos, count=shape(field_hr2), &
+ values=field_hr2), "dry_deposition_velocity(m)")
+ endif
end if
!..accumulated dry deposition
@@ -352,6 +386,24 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
values=field_hr1), "set_accwd(m)")
end if
+ if (output_component(m)%has_wetdep) then
+ block
+ use snapfldML, only: wscav
+ use snapdimML, only: nk
+ integer :: k
+ if (output_wetdeprate) then
+ do k=2,nk
+ call hres_field(wscav(:,:,k,m), field_hr1)
+ call check(nf90_put_var(iunit, varid%comp(m)%wetscavrate, start=[1,1,k-1,ihrs_pos], &
+ count=[nx*output_resolution_factor,ny*output_resolution_factor,1,1], values=field_hr1), "wscavrate")
+ end do
+ endif
+ end block
+ endif
+
+ !..instant part of Bq in boundary layer
+ if(idebug == 1) call ftest('pbq', field_hr3, contains_undef=.true.)
+
!..average part of Bq in boundary layer
scale=100.
@@ -540,6 +592,39 @@ subroutine fldout_nc(filename, itime,tf1,tf2,tnow, &
[nx*output_resolution_factor, ny*output_resolution_factor, 1, 1], rt1, rt2)
endif
+ if (output_vd_debug) then
+ block
+ use snapfldml, only: t2m, xflux, yflux, z0, hflux, leaf_area_index, &
+ roa, ustar, monin_l, raero, vs, rs, 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(z0, field_hr1)
+ call check(nf90_put_var(iunit, varid%z0, start=ipos, count=isize, values=field_hr1))
+ call hres_field(hflux, field_hr1)
+ call check(nf90_put_var(iunit, varid%hflux, start=ipos, count=isize, values=field_hr1))
+ call hres_field(leaf_area_index, field_hr1)
+ call check(nf90_put_var(iunit, varid%lai, start=ipos, count=isize, values=field_hr1))
+ call hres_field(roa, field_hr1)
+ call check(nf90_put_var(iunit, varid%roa, start=ipos, count=isize, values=field_hr1))
+ call hres_field(ustar, field_hr1)
+ call check(nf90_put_var(iunit, varid%ustar, start=ipos, count=isize, values=field_hr1))
+ call hres_field(monin_l, field_hr1)
+ call check(nf90_put_var(iunit, varid%monin_l, start=ipos, count=isize, values=field_hr1))
+ call hres_field(raero, field_hr1)
+ call check(nf90_put_var(iunit, varid%raero, start=ipos, count=isize, values=field_hr1))
+ call hres_field(vs, field_hr1)
+ call check(nf90_put_var(iunit, varid%vs, start=ipos, count=isize, values=field_hr1))
+ call hres_field(rs, field_hr1)
+ call check(nf90_put_var(iunit, varid%rs, start=ipos, count=isize, values=field_hr1))
+ end block
+ endif
+
! reset fields
do m=1,nocomp
if (output_component(m)%has_drydep) then
@@ -576,7 +661,7 @@ subroutine write_ml_fields(iunit, varid, average, ipos_in, isize, rt1, rt2)
type(Particle) :: part
real :: avg, total, dh
integer :: ivlvl
- integer :: i, j, k, loop, m, maxage, n, npl
+ integer :: i, j, k, m, maxage, n, npl
integer :: ipos(4)
logical :: inactivated_ ! dummy param
@@ -671,7 +756,7 @@ subroutine write_ml_fields(iunit, varid, average, ipos_in, isize, rt1, rt2)
end subroutine
-subroutine nc_declare(iunit, dimids, varid, varname, units, stdname, chunksize)
+subroutine nc_declare(iunit, dimids, varid, varname, units, stdname, chunksize, datatype)
USE snapdebug, only: iulog
integer, intent(in) :: iunit
integer, intent(out) :: varid
@@ -681,13 +766,21 @@ subroutine nc_declare(iunit, dimids, varid, varname, units, stdname, chunksize)
character(len=*), intent(in) :: varname
character(len=*), intent(in), optional :: stdname
integer, intent(in), optional :: chunksize(:)
+ integer, intent(in), optional :: datatype
+
+ integer :: datatype_internal
write(iulog,"('declaring ' (a) ' ' (a) )",advance="NO") varname, units
if (present(stdname)) write(iulog, "(' ' (a))",advance="NO") trim(stdname)
write(iulog,'()',advance="YES")
+ if (present(datatype)) then
+ datatype_internal = datatype
+ else
+ datatype_internal = NF90_FLOAT
+ endif
call check(nf90_def_var(iunit, TRIM(varname), &
- NF90_FLOAT, dimids, varid), "def_"//varname)
+ datatype_internal, dimids, varid), "def_"//varname)
if (present(chunksize)) then
call check(nf90_def_var_chunking(iunit, varid, NF90_CHUNKED, chunksize))
call check(nf90_def_var_deflate(iunit, varid, &
@@ -987,22 +1080,24 @@ subroutine initialize_output(filename, itime, ierror)
USE snapfilML, only: ncsummary, nctitle, simulation_start
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
+ aircraft_doserate_threshold, &
+ output_vd, output_column, output_vd_debug
USE snapfldML, only: &
garea, &
xm, ym, &
nhfout
USE snapparML, only: nocomp, output_component, def_comp
USE ftestML, only: ftest
- USE snapdimML, only: nx, ny, nk, output_resolution_factor
+ USE snapdimML, only: nx, ny, nk, output_resolution_factor, hres_field
USE particleML, only: Particle
+ use snapmetML, only: met_params
character(len=*), intent(in) :: filename
type(datetime_t), intent(in) :: itime
integer, intent(out) :: ierror
integer :: iunit
- integer :: m, mm
+ integer :: m
character(len=256) :: string
integer :: dimids2d(2), dimids3d(3), dimids4d(4)
integer :: chksz3d(3), chksz4d(4)
@@ -1035,7 +1130,8 @@ subroutine initialize_output(filename, itime, ierror)
call nc_set_projection(iunit, dimid%x, dimid%y, &
igtype, gparam, garea, xm, ym, simulation_start)
- if (imodlevel) then
+ if (imodlevel .or. output_wetdeprate .or. &
+ (precipitation_in_output .and. met_params%use_3d_precip)) then
call nc_set_vtrans(iunit, dimid%k, varid%k, varid%ap, varid%b)
endif
@@ -1074,6 +1170,12 @@ subroutine initialize_output(filename, itime, ierror)
call nc_declare(iunit, dimids3d, varid%prc, &
"lwe_precipitation_rate", units="mm/("//itoa(nhfout)//"hr)", &
stdname="lwe_precipitation_rate", chunksize=chksz3d)
+
+ if (met_params%use_3d_precip) then
+ call nc_declare(iunit, dimids3d, varid%instant_prc, &
+ "lwe_instantaneous_precipitation_rate_ml", units="mm/hr", &
+ stdname="lwe_instantaneous_precipitation_rate_ml", chunksize=chksz3d)
+ endif
endif
call nc_declare(iunit, dimids3d, varid%ihbl, &
@@ -1100,6 +1202,11 @@ subroutine initialize_output(filename, itime, ierror)
endif
endif
+ if (output_wetdeprate) then
+ call nc_declare(iunit, dimids3d, varid%wetdeprate, &
+ "wetdeprate", units="m/s", chunksize=chksz3d)
+ endif
+
call check(nf90_def_var(iunit, "components", NF90_CHAR, [dimid%maxcompname, dimid%nocomp], varid%components))
do m=1,nocomp
@@ -1126,6 +1233,7 @@ subroutine initialize_output(filename, itime, ierror)
trim(output_component(m)%name)//"_acc_dry_deposition", &
units="Bq/m2", chunksize=chksz3d)
end if
+
if (output_component(m)%has_wetdep) then
call nc_declare(iunit, dimids3d, varid%comp(m)%iwd, &
trim(output_component(m)%name)//"_wet_deposition", &
@@ -1133,7 +1241,12 @@ subroutine initialize_output(filename, itime, ierror)
call nc_declare(iunit, dimids3d, varid%comp(m)%accwd, &
trim(output_component(m)%name)//"_acc_wet_deposition", &
units="Bq/m2", chunksize=chksz3d)
+ if (output_wetdeprate) then
+ call nc_declare(iunit, dimids4d, varid%comp(m)%wetscavrate, &
+ trim(output_component(m)%name)//"_wetdeprate", units="1/s", chunksize=chksz4d)
+ endif
end if
+
if (imodlevel) then
if (modleveldump > 0.) then
string = trim(output_component(m)%name)//"_concentration_dump_ml"
@@ -1153,6 +1266,71 @@ subroutine initialize_output(filename, itime, ierror)
trim(output_component(m)%name)//"_column_concentration", &
units="Bq/m2", chunksize=chksz3d)
endif
+
+ if (output_component(m)%has_drydep .and. output_vd) then
+ call nc_declare(iunit, dimids3d, varid%comp(m)%vd, &
+ trim(output_component(m)%name)//"_dry_deposition_velocity", &
+ units="m/s", chunksize=chksz3d)
+ endif
+
+ if (output_vd_debug) then
+ block
+ use snapmetml, only: downward_momentum_flux_units, surface_heat_flux_units, &
+ leaf_area_index_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, &
+ chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%hflux, &
+ "hflux", units=surface_heat_flux_units, &
+ chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%z0, &
+ "z0", units=surface_roughness_length_units, &
+ chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%lai, &
+ "lai", units=leaf_area_index_units, &
+ chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%t2m, &
+ "t2m", units=temp_units, &
+ chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%roa, &
+ "roa", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%ustar, &
+ "ustar", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%monin_l, &
+ "monin_l", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%raero, &
+ "raero", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%vs, &
+ "vs", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%rs, &
+ "rs", units="??", chunksize=chksz3d)
+ call nc_declare(iunit, dimids3d, varid%ps_vd, &
+ "ps_vd", units="hPa", chunksize=chksz3d)
+ end block
+ endif
+
+ block
+ use iso_fortran_env, only: int8
+ use drydepml, only: largest_landfraction_file, classnr
+ integer(kind=int8), allocatable :: classnr_hr(:,:)
+ if (largest_landfraction_file /= "not set") then
+ call nc_declare(iunit, dimids2d, varid%landfraction, &
+ "largest_land_fraction", units="1", datatype=NF90_BYTE)
+ call hres_field(classnr, classnr_hr)
+ call check(nf90_put_var(iunit, varid%landfraction, start=[1, 1], count=shape(classnr_hr), &
+ values=classnr_hr), "Put landfraction")
+ call check(nf90_put_att(iunit, varid%landfraction, "flag_values", &
+ [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]))
+ call check(nf90_put_att(iunit, varid%landfraction, "flag_meanings", &
+ "sea inland_water tundra_and_desert ice_and_ice_sheets urban crops" // &
+ " grass wetlands evergreen_needleleaf deciduous_broadleaf" // &
+ " mixed_forest shrubs_and_interrupted_woodlands"))
+ endif
+ end block
+
end do
if (itotcomp == 1) then
call nc_declare(iunit, dimids3d, varid%icblt, &
@@ -1205,7 +1383,7 @@ subroutine unload()
end subroutine
subroutine get_varids(iunit, varid, ierror)
- USE snapparML, only: nocomp, output_component, def_comp
+ USE snapparML, only: nocomp, output_component, nocomp
USE snapgrdML, only: imodlevel, modleveldump
integer, intent(in) :: iunit
type(common_var), intent(out) :: varid
@@ -1227,6 +1405,8 @@ subroutine get_varids(iunit, varid, ierror)
! Optional
ierror = nf90_inq_varid(iunit, "precipitation_amount_acc", varid%accum_prc)
if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "lwe_instantaneous_precipitation_rate_ml", varid%instant_prc)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
ierror = nf90_inq_varid(iunit, "lwe_precipitation_rate", varid%prc)
if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
ierror = nf90_inq_varid(iunit, "air_pressure_at_sea_level", varid%mslp)
@@ -1260,6 +1440,33 @@ subroutine get_varids(iunit, varid, ierror)
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)
+ 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
+ ierror = nf90_inq_varid(iunit, "lai", varid%lai)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "z0", varid%z0)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "t2m", varid%t2m)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "roa", varid%roa)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "ustar", varid%ustar)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "monin_l", varid%monin_l)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "raero", varid%raero)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "vs", varid%vs)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "rs", varid%rs)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ ierror = nf90_inq_varid(iunit, "ps_vd", varid%ps_vd)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+
do m=1,nocomp
varname = trim(output_component(m)%name) // "_concentration"
ierror = nf90_inq_varid(iunit, varname, varid%comp(m)%ic)
@@ -1285,6 +1492,10 @@ subroutine get_varids(iunit, varid, ierror)
ierror = nf90_inq_varid(iunit, varname, varid%comp(m)%accdd)
if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+ varname = trim(output_component(m)%name) // "_dry_deposition_velocity"
+ ierror = nf90_inq_varid(iunit, varname, varid%comp(m)%vd)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+
varname = trim(output_component(m)%name) // "_wet_deposition"
ierror = nf90_inq_varid(iunit, varname, varid%comp(m)%iwd)
if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
@@ -1293,6 +1504,11 @@ subroutine get_varids(iunit, varid, ierror)
ierror = nf90_inq_varid(iunit, varname, varid%comp(m)%accwd)
if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+
+ ierror = nf90_inq_varid(iunit, trim(output_component(m)%name)//"_wetdeprate", varid%comp(m)%wetscavrate)
+ if (ierror /= NF90_NOERR .and. .not. ierror == NF90_ENOTVAR) return
+
+
if (imodlevel) then
if (modleveldump > 0.) then
varname = trim(output_component(m)%name) // "_concentration_dump_ml"
@@ -1315,8 +1531,7 @@ subroutine get_varids(iunit, varid, ierror)
!> accumulate model level fields only
subroutine accumulate_ml_field()
- USE snapgrdML, only: imodlevel, modlevel_is_average, &
- ivlayer
+ USE snapgrdML, only: imodlevel, ivlayer
USE snapfldML, only: ml_bq
USE snapparML, only: def_comp
USE snapdimml, only: hres_pos
@@ -1351,7 +1566,7 @@ subroutine accumulate_fields(tf1, tf2, tnow, tstep, nsteph)
ps1, ps2, t1_abs, t2_abs, aircraft_doserate_scratch, aircraft_doserate, &
aircraft_doserate_threshold_height
USE snapdimml, only: nx, ny, nk, output_resolution_factor, hres_pos, lres_pos
- USE snapparML, only: nocomp, def_comp, output_component
+ USE snapparML, only: nocomp, def_comp
USE ftestML, only: ftest
USE releaseML, only: npart
USE particleML, only: pdata, Particle
diff --git a/src/common/posint.f90 b/src/common/posint.f90
index 9ff5f839..9117d448 100644
--- a/src/common/posint.f90
+++ b/src/common/posint.f90
@@ -19,45 +19,12 @@ module posintML
implicit none
private
- public posint, posint_init
-
- real, save :: vminprec = -1.
- real, parameter :: precmin = 0.01
+ public :: posint
contains
-!> initialisation function for ::posint
-subroutine posint_init()
- USE snapgrdML, only: alevel, blevel, vlevel
- USE snapdimML, only: nk
- USE snapdebug, only: iulog
-
- integer :: k
- real :: p1,p2
- real, parameter :: plim = 550.0
-
- if(vminprec < 0.) then
- p2 = 1000.
- k = 1
- do while (p2 > plim .AND. k < nk)
- k = k+1
- p1 = p2
- p2 = alevel(k) + blevel(k)*1000.
- end do
- if(k > 1) then
- vminprec = vlevel(k-1) &
- + (vlevel(k)-vlevel(k-1))*(p1-plim)/(p1-p2)
- else
- vminprec = vlevel(nk)
- end if
- write(iulog,*) 'POSINT. precmin,vminprec: ',precmin,vminprec
- end if
-end subroutine
-
!> Purpose: Interpolation of boundary layer top and height
!> and precipitation intensity to particle positions
-!>
-!> must call ::posint_init first
subroutine posint(part,tf1,tf2,tnow,pextra)
USE particleML, only: Particle, extraParticle
USE snapgrdML, only: gparam
@@ -133,13 +100,5 @@ subroutine posint(part,tf1,tf2,tnow,pextra)
pextra%rmy=rmy/dygrid
pextra%prc=pr
- !..reset precipitation to zero if pressure less than approx. 550 hPa.
- !..and if less than a minimum precipitation intensity (mm/hour)
- if(part%z < vminprec .OR. &
- pextra%prc < precmin) then
- pextra%prc=0.
- endif
-
- return
end subroutine posint
end module posintML
diff --git a/src/common/readfield_fi.f90 b/src/common/readfield_fi.f90
index ba735c0e..a26d4183 100644
--- a/src/common/readfield_fi.f90
+++ b/src/common/readfield_fi.f90
@@ -32,11 +32,12 @@ module readfield_fiML
implicit none
private
- public readfield_fi, fi_checkload, check, fimex_open
+ public :: readfield_fi, fi_checkload, check, fimex_open, read_largest_landfraction
!> @brief load and check an array from a source
interface fi_checkload
- module procedure fi_checkload1d, fi_checkload2d, fi_checkload3d
+ module procedure :: fi_checkload1d, fi_checkload2d, fi_checkload3d, &
+ fi_checkload2d_64, fi_checkload3d_64
end interface
contains
@@ -293,28 +294,34 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
!..mean sea level pressure, not used in computations,
!..(only for output to results file)
if (imslp /= 0) then
- if (.NOT. met_params%mslpv == '') then
+ if (met_params%mslpv /= '') then
+ call fi_checkload(fio, met_params%mslpv, pressure_units, pmsl2(:, :), nt=timepos, nr=nr, nz=1)
+ else if (met_params%psv /= '') then
+ call fi_checkload(fio, met_params%psv, pressure_units, pmsl2(:, :), nt=timepos, nr=nr, nz=1)
+ else
write (iulog, *) 'Mslp not found. Not important.'
imslp = 0
- else
- call fi_checkload(fio, met_params%mslpv, pressure_units, pmsl2(:, :), nt=timepos, nr=nr, nz=1)
- end if
+ endif
end if
+ if (first_time_read) then
+ call compute_vertical_coords(alev, blev, ptop)
+ endif
+
if (met_params%need_precipitation) then
call read_precipitation(fio, nhdiff_precip, timepos, timeposm1)
else
precip = 0.0
endif
+ call read_drydep_required_fields(fio, timepos, timeposm1, nr)
+
call check(fio%close(), "close fio")
! first time initialized data
if (first_time_read) then
first_time_read = .false.
- call compute_vertical_coords(alev, blev, ptop)
-
!..compute map ratio
call mapfield(1, 0, igtype, gparam, nx, ny, xm, ym, &
xm, & ! Ignored when icori = 0
@@ -448,10 +455,10 @@ subroutine readfield_fi(istep, backward, itimei, ihr1, ihr2, itimefi, ierror)
end do
end if
- return
end subroutine readfield_fi
-!> read precipitation
+!> read precipitation fields and precompute
+!> for wet deposition
subroutine read_precipitation(fio, nhdiff, timepos, timeposm1)
use iso_fortran_env, only: error_unit
use snapdebug, only: iulog
@@ -459,6 +466,7 @@ subroutine read_precipitation(fio, nhdiff, timepos, timeposm1)
precip_rate_units, precip_units_ => precip_units, precip_units_fallback
use snapfldML, only: field1, field2, field3, field4, precip, &
enspos, precip
+ use wetdepML, only: requires_extra_precip_fields, wetdep_precompute
!> open netcdf file
TYPE(FimexIO), intent(inout) :: fio
@@ -478,6 +486,74 @@ subroutine read_precipitation(fio, nhdiff, timepos, timeposm1)
nr = enspos + 1
if (enspos <= 0) nr = 1
+ if (met_params%use_3d_precip) then
+ block ! read_precip
+ use snaptabML, only: g
+ use snapfldML, only: ps2, precip3d, cw3d, cloud_cover
+ use snapgrdML, only: ahalf, bhalf, klevel
+ use snapdimML, only: nx, ny, nk
+ real(real64), allocatable :: rain_in_air(:,:), graupel_in_air(:,:), snow_in_air(:,:)
+ real(real64), allocatable :: cloud_water(:,:), cloud_ice(:,:)
+ real(real64), allocatable :: pdiff(:,:)
+
+ integer :: ilevel, k
+ character(len=*), parameter :: mass_fraction_units = "kg/kg"
+
+ allocate(rain_in_air(nx,ny),graupel_in_air(nx,ny),snow_in_air(nx,ny),pdiff(nx,ny))
+ allocate(cloud_water(nx,ny),cloud_ice(nx,ny))
+
+ precip3d(:,:,:) = 0.0
+ cw3d(:,:,:) = 0.0
+
+ do k=nk,2
+ ilevel = klevel(k)
+ call fi_checkload(fio, "mass_fraction_of_rain_in_air_ml", mass_fraction_units, &
+ rain_in_air, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, "mass_fraction_of_graupel_in_air_ml", mass_fraction_units, &
+ graupel_in_air, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, "mass_fraction_of_snow_in_air_ml", mass_fraction_units, &
+ snow_in_air, nt=timepos, nz=ilevel, nr=nr)
+
+ where (rain_in_air < 0.0)
+ rain_in_air = 0.0
+ end where
+ where (graupel_in_air < 0.0)
+ graupel_in_air = 0.0
+ end where
+ where (snow_in_air < 0.0)
+ snow_in_air = 0.0
+ end where
+
+ pdiff(:,:) = 100*( (ahalf(k-1) - ahalf(k)) + (bhalf(k-1) - bhalf(k))*ps2 )
+
+ precip3d(:,:,k) = rain_in_air + graupel_in_air + snow_in_air
+ precip3d(:,:,k) = precip3d(:,:,k) * pdiff / g
+
+ call fi_checkload(fio, "mass_fraction_of_cloud_condensed_water_in_air_ml", mass_fraction_units, &
+ cloud_water, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, "mass_fraction_of_cloud_ice_in_air_ml", mass_fraction_units, &
+ cloud_ice, nt=timepos, nz=ilevel, nr=nr)
+
+ where (cloud_water < 0.0)
+ cloud_water = 0.0
+ end where
+ where (cloud_ice < 0.0)
+ cloud_ice = 0.0
+ end where
+ cw3d(:,:,k) = cloud_water + cloud_ice
+ cw3d(:,:,k) = cw3d(:,:,k) * pdiff / g
+
+ call fi_checkload(fio, "cloud_area_fraction_ml", "%", &
+ cloud_cover(:,:,k), nt=timepos, nz=ilevel, nr=nr)
+ enddo
+
+ ! Accumulate precipitation from top down
+ do k=nk-1,2,-1
+ precip3d(:,:,k) = precip3d(:,:,k) + precip3d(:,:,k+1)
+ enddo
+ end block
+ end if
+
if (met_params%precaccumv /= '') then
!..precipitation between input time 't1' and 't2'
if (timepos == 1) then
@@ -534,8 +610,182 @@ subroutine read_precipitation(fio, nhdiff, timepos, timeposm1)
where (precip < 0.0)
precip = 0.0
end where
+
+ if (requires_extra_precip_fields()) then
+ if ((met_params%mass_fraction_rain_in_air /= "") .and. &
+ (met_params%mass_fraction_cloud_condensed_water_in_air /= "")) then
+ call read_extra_precipitation_fields(fio, timepos)
+ else if (met_params%mass_fraction_cloud_condensed_water_in_air /= "") then
+ call read_extra_precipitation_fields_infer_3d_precip(fio, timepos)
+ else
+ error stop "Can not read extra 3D precipition for this meteorology"
+ endif
+ endif
+
+ call wetdep_precompute()
end subroutine read_precipitation
+ subroutine read_extra_precipitation_fields(fio, timepos)
+ use iso_fortran_env, only: error_unit
+ use snaptabML, only: g
+ use snapfldML, only: ps2, precip3d, cw3d, cloud_cover, enspos
+ use snapgrdML, only: ahalf, bhalf, klevel
+ use snapdimML, only: nx, ny, nk
+ use snapmetML, only: mass_fraction_units, cloud_fraction_units, met_params
+!> open netcdf file
+ TYPE(FimexIO), intent(inout) :: fio
+!> timestep in file
+ integer, intent(in) :: timepos
+
+ real(real64), allocatable :: rain_in_air(:,:), graupel_in_air(:,:), snow_in_air(:,:)
+ real(real64), allocatable :: cloud_water(:,:), cloud_ice(:,:)
+ real(real64), allocatable :: pdiff(:,:)
+
+ integer :: ilevel, k, nr
+
+!.. get the correct ensemble/realization position, nr starting with 1, enspos starting with 0
+ nr = enspos + 1
+ if (enspos <= 0) nr = 1
+
+ allocate(rain_in_air(nx,ny),graupel_in_air(nx,ny),snow_in_air(nx,ny),pdiff(nx,ny))
+ allocate(cloud_water(nx,ny),cloud_ice(nx,ny))
+
+ precip3d(:,:,:) = 0.0
+ cw3d(:,:,:) = 0.0
+
+ do k=nk,2,-1
+ ilevel = klevel(k)
+ call fi_checkload(fio, met_params%mass_fraction_rain_in_air, mass_fraction_units, &
+ rain_in_air, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, met_params%mass_fraction_graupel_in_air, mass_fraction_units, &
+ graupel_in_air, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, met_params%mass_fraction_snow_in_air, mass_fraction_units, &
+ snow_in_air, nt=timepos, nz=ilevel, nr=nr)
+
+ where (rain_in_air < 0.0)
+ rain_in_air = 0.0
+ end where
+ where (graupel_in_air < 0.0)
+ graupel_in_air = 0.0
+ end where
+ where (snow_in_air < 0.0)
+ snow_in_air = 0.0
+ end where
+
+ pdiff(:,:) = 100*( (ahalf(k-1) - ahalf(k)) + (bhalf(k-1) - bhalf(k))*ps2 )
+
+ precip3d(:,:,k) = rain_in_air + graupel_in_air + snow_in_air
+ precip3d(:,:,k) = precip3d(:,:,k) * pdiff / g
+
+ call fi_checkload(fio, met_params%mass_fraction_cloud_condensed_water_in_air, mass_fraction_units, &
+ cloud_water, nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, met_params%mass_fraction_cloud_ice_in_air, mass_fraction_units, &
+ cloud_ice, nt=timepos, nz=ilevel, nr=nr)
+
+ where (cloud_water < 0.0)
+ cloud_water = 0.0
+ end where
+ where (cloud_ice < 0.0)
+ cloud_ice = 0.0
+ end where
+ cw3d(:,:,k) = cloud_water + cloud_ice
+ cw3d(:,:,k) = cw3d(:,:,k) * pdiff / g
+
+ call fi_checkload(fio, met_params%cloud_fraction, cloud_fraction_units, &
+ cloud_cover(:,:,k), nt=timepos, nz=ilevel, nr=nr)
+ enddo
+ end subroutine
+
+ !> Read and convert fields from ecemep input to what we need for 3D precip
+ subroutine read_extra_precipitation_fields_infer_3d_precip(fio, timepos)
+ use iso_fortran_env, only: error_unit
+ use snaptabML, only: g
+ use snapfldML, only: ps2, precip3d, cw3d, cloud_cover, enspos, precip
+ use snapgrdML, only: ahalf, bhalf, klevel
+ use snapdimML, only: nx, ny, nk
+ use snapmetML, only: mass_fraction_units, cloud_fraction_units, met_params
+!> open netcdf file
+ TYPE(FimexIO), intent(inout) :: fio
+!> timestep in file
+ integer, intent(in) :: timepos
+
+ real(real64), allocatable :: normaliser(:,:)
+ real(real64), allocatable :: pdiff(:,:)
+ real(real64), allocatable :: cloud_water(:,:), cloud_ice(:,:)
+
+ integer :: ilevel, k, nr, i, j
+
+!.. get the correct ensemble/realization position, nr starting with 1, enspos starting with 0
+ nr = enspos + 1
+ if (enspos <= 0) nr = 1
+
+ allocate(pdiff(nx,ny))
+ allocate(normaliser(nx,ny))
+ allocate(cloud_water(nx,ny), cloud_ice(nx,ny))
+
+ precip3d(:,:,:) = 0.0
+ cw3d(:,:,:) = 0.0
+
+ normaliser(:,:) = 0.0
+
+ do k=nk,2,-1
+ ilevel = klevel(k)
+
+ pdiff(:,:) = 100*( (ahalf(k-1) - ahalf(k)) + (bhalf(k-1) - bhalf(k))*ps2 )
+ call fi_checkload(fio, met_params%mass_fraction_cloud_condensed_water_in_air, mass_fraction_units, &
+ cloud_water(:,:), nt=timepos, nz=ilevel, nr=nr)
+ call fi_checkload(fio, met_params%mass_fraction_cloud_ice_in_air, mass_fraction_units, &
+ cloud_ice(:,:), nt=timepos, nz=ilevel, nr=nr)
+
+ cw3d(:,:,k) = (abs(cloud_water(:,:)) + abs(cloud_ice(:,:))) * pdiff / g
+
+ ! Use cloud water to assign precipitation at model levels
+ normaliser(:,:) = normaliser + cw3d(:,:,k)
+ precip3d(:,:,k) = precip * cw3d(:,:,k)
+
+ call fi_checkload(fio, met_params%cloud_fraction, cloud_fraction_units, &
+ cloud_cover(:,:,k), nt=timepos, nz=ilevel, nr=nr)
+ enddo
+
+ block
+ use snapgrdML, only: ivlevel
+ integer :: klimit
+ klimit = ivlevel(nint(0.67*10000.0))
+ do k=nk,2,-1
+ do j=1,ny
+ do i=1,nx
+ if (normaliser(i,j) > 0.0) then
+ precip3d(i,j,k) = precip3d(i,j,k) / normaliser(i,j)
+ elseif (k == klimit) then
+ ! Put all precip at level closest to 0.67, analoguous
+ ! with the old formulation of the precipitation
+ precip3d(i,j,k) = precip(i,j)
+ endif
+ enddo
+ enddo
+ enddo
+ end block
+
+ ! block ! Debug
+ ! use snapfldML, only: garea
+ ! real :: total_precip
+ ! real :: total_precip2
+ ! real, allocatable :: tmp(:,:)
+
+ ! allocate(tmp(nx,ny))
+
+ ! tmp(:,:) = precip * garea
+ ! total_precip = sum(tmp)
+ ! tmp(:,:) = sum(precip3d, dim=3)
+ ! tmp(:,:) = tmp * garea
+ ! total_precip2 = sum(tmp)
+
+ ! write(*,*) "PRECIP LOSS: ", total_precip2 - total_precip
+ ! write(*,*) "Precip from 2D fields: ", total_precip
+ ! write(*,*) "Precip from 3d fields: ", total_precip2
+ ! end block
+ end subroutine
+
subroutine check(status, errmsg)
integer, intent(in) :: status
character(len=*), intent(in), optional :: errmsg
@@ -591,6 +841,21 @@ subroutine fi_checkload2d(fio, varname, units, field, nt, nz, nr, ierror)
deallocate(zfield)
end subroutine fi_checkload2d
+ subroutine fi_checkload2d_64(fio, varname, units, field, nt, nz, nr, ierror)
+ TYPE(FimexIO), intent(inout) :: fio
+ character(len=*), intent(in) :: varname, units
+ integer, intent(in), optional :: nt, nz, nr
+ real(real64), intent(out) :: field(:, :)
+ integer, intent(out), optional :: ierror
+
+ real(kind=real64), dimension(:), allocatable, target :: zfield
+
+ call fi_checkload_intern(fio, varname, units, zfield, nt, nz, nr, ierror)
+
+ field(:,:) = reshape(zfield, shape(field))
+ deallocate(zfield)
+ end subroutine fi_checkload2d_64
+
subroutine fi_checkload3d(fio, varname, units, field, nt, nz, nr, ierror)
TYPE(FimexIO), intent(inout) :: fio
character(len=*), intent(in) :: varname, units
@@ -606,6 +871,21 @@ subroutine fi_checkload3d(fio, varname, units, field, nt, nz, nr, ierror)
deallocate(zfield)
end subroutine fi_checkload3d
+ subroutine fi_checkload3d_64(fio, varname, units, field, nt, nz, nr, ierror)
+ TYPE(FimexIO), intent(inout) :: fio
+ character(len=*), intent(in) :: varname, units
+ integer, intent(in), optional :: nt, nz, nr
+ real(real64), intent(out) :: field(:, :, :)
+ integer, intent(out), optional :: ierror
+
+ real(kind=real64), dimension(:), allocatable, target :: zfield
+
+ call fi_checkload_intern(fio, varname, units, zfield, nt, nz, nr, ierror)
+
+ field(:,:,:) = reshape(zfield, shape(field))
+ deallocate(zfield)
+ end subroutine fi_checkload3d_64
+
!> internal implementation, allocating the zfield
subroutine fi_checkload_intern(fio, varname, units, zfield, nt, nz, nr, ierror)
USE Fimex, ONLY: FimexIO, AXIS_GeoX, AXIS_GeoY, AXIS_Lon, AXIS_Lat, AXIS_GeoZ, &
@@ -689,4 +969,114 @@ 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_drydep_required_fields(fio, timepos, timeposm1, nr)
+ 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, leaf_area_index_units
+ use drydepml, only: drydep_precompute, 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, leaf_area_index, t2m, &
+ roa, ustar, monin_l, raero, vs, rs
+ type(FimexIO), intent(inout) :: fio
+ integer, intent(in) :: timepos, timeposm1
+ integer, intent(in) :: nr
+
+ real, allocatable :: tmp1(:, :), tmp2(:, :)
+ integer :: i, mm
+ real(real64) :: diam, dens
+
+ 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)
+ 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
+
+ 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
+ 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)
+ else if (timepos == 13) then
+ ! Weird AROME data is invalid at t=12
+ 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+1, nr=nr)
+ hflux(:,:) = (tmp2 - tmp1)/2
+ else if (timepos == 14) then
+ ! Weird AROME data is invalid at t=13
+ call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, tmp1(:, :), nt=timeposm1-1, nr=nr)
+ call fi_checkload(fio, met_params%hflux, surface_heat_flux_units, tmp2(:, :), nt=timepos, nr=nr)
+ hflux(:,:) = (tmp2 - tmp1)/2
+ 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
+ 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)
+ call fi_checkload(fio, met_params%leaf_area_index, leaf_area_index_units, leaf_area_index(:, :), nt=timepos, nr=nr)
+ where (ieee_is_nan(leaf_area_index))
+ leaf_area_index = 0.0
+ endwhere
+ call fi_checkload(fio, met_params%t2m, temp_units, t2m(:, :), nt=timepos, nr=nr)
+
+ do i=1,ncomp
+ mm = run_comp(i)%to_defined
+
+ if (def_comp(mm)%kdrydep == 1) then
+ diam = 2*def_comp(mm)%radiusmym*1e-6
+ dens = def_comp(mm)%densitygcm3*1e3
+ call drydep_precompute(ps2*100, t2m, yflux, xflux, z0, &
+ hflux, leaf_area_index, real(diam), real(dens), classnr, vd_dep(:, :, i), &
+ roa, ustar, monin_l, raero, vs, rs)
+ endif
+ end do
+ end subroutine
+
+ subroutine read_largest_landfraction(inputfile)
+ use snapdimML, only: nx, ny
+ use drydepml, only: preprocess_landfraction
+ use fimex, only: INTERPOL_NEAREST_NEIGHBOR
+ use ISO_C_BINDING, only: C_INT
+ character(len=*), intent(in) :: inputfile
+
+ type(fimexIO) :: fio, fio_intern
+
+ real(kind=real32), allocatable :: arr(:,:)
+
+ if (fint%method >= 0) then
+ call check(fio_intern%open(inputfile, "", "nc4"), &
+ "can't open largest landfraction file")
+ call check(fio%interpolate(fio_intern, INTERPOL_NEAREST_NEIGHBOR, fint%proj, fint%x_axis, &
+ fint%y_axis, fint%unit_is_degree), &
+ "Can't interpolate largest landfraction file")
+ else
+ call check(fio%open(inputfile, "", "nc4"), "Can't open largest landraction file")
+ endif
+
+ allocate(arr(nx, ny))
+ call fi_checkload(fio, "Main_Nature_Cover", "1", arr)
+
+ call preprocess_landfraction(arr)
+ end subroutine
+
end module readfield_fiML
diff --git a/src/common/readfield_nc.f90 b/src/common/readfield_nc.f90
index 94a4fb65..4e0705e5 100644
--- a/src/common/readfield_nc.f90
+++ b/src/common/readfield_nc.f90
@@ -359,11 +359,13 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
!..mean sea level pressure, not used in computations,
!..(only for output to results file)
if(imslp /= 0) then
- if ( .NOT. met_params%mslpv == '') then
+ if (met_params%mslpv /= '') then
+ call nfcheckload(ncid, met_params%mslpv, start3d, count3d, pmsl2(:,:))
+ else if (met_params%psv /= '') then
+ call nfcheckload(ncid, met_params%psv, start3d, count3d, pmsl2(:,:))
+ else
write(iulog,*) 'Mslp not found. Not important.'
imslp=0
- else
- call nfcheckload(ncid, met_params%mslpv, start3d, count3d, pmsl2(:,:))
end if
end if
@@ -514,7 +516,7 @@ subroutine readfield_nc(istep, backward, itimei, ihr1, ihr2, &
end do
end if
- return
+ call read_additional_drydep_fields()
end subroutine readfield_nc
!> Reads `units` attribute of the precipitation variable
@@ -899,7 +901,6 @@ subroutine nfcheckload3d(ncid, varname, start, length, field, return_status)
end if
end subroutine nfcheckload3d
-
subroutine compute_vertical_coords(alev, blev, ptop)
use iso_fortran_env, only: error_unit
use snapgrdML, only: alevel, blevel, vlevel, ivcoor, klevel, &
@@ -968,4 +969,80 @@ subroutine compute_vertical_coords(alev, blev, ptop)
vhalf(nk) = vlevel(nk)
end subroutine
+subroutine read_additional_drydep_fields()
+ use drydepml, only: requires_extra_fields_to_be_read
+ if (requires_extra_fields_to_be_read()) then
+ error stop "Reading of extra dry deposition fields is not implemented for netCDF"
+ endif
+end subroutine
+
+ subroutine compute_vertical_levels(alev, blev, ptop)
+ use iso_fortran_env, only: error_unit
+ use snapgrdML, only: alevel, blevel, ahalf, bhalf, vlevel, vhalf, klevel, &
+ ivcoor
+ use snapdimML, only: nk
+ use snapmetML, only: met_params
+ use snaptabML, only: standard_atmosphere
+
+ real, intent(in) :: alev(:), blev(:)
+ real, intent(in) :: ptop
+
+ integer :: k
+
+ do k = 2, nk
+ alevel(k) = alev(k)
+ blevel(k) = blev(k)
+ end do
+
+
+ if (ivcoor == 2) then
+ !..sigma levels (norlam)
+ do k = 2, nk
+ alevel(k) = ptop*(1.-blevel(k))
+ end do
+ end if
+
+ !..surface
+ alevel(1) = 0.
+ blevel(1) = 1.
+
+ if (ivcoor == 2) then
+ !..sigma levels ... vlevel=sigma
+ vlevel(:) = blevel
+ elseif (ivcoor == 10) then
+ !..eta (hybrid) levels ... vlevel=eta (eta as defined in Hirlam)
+ vlevel(:) = alevel / standard_atmosphere + blevel
+ else
+ write (error_unit, *) 'PROGRAM ERROR. ivcoor= ', ivcoor
+ error stop 255
+ end if
+
+ !..half levels where height is found,
+ !..alevel and blevel are in the middle of each layer
+ ahalf(1) = alevel(1)
+ bhalf(1) = blevel(1)
+ vhalf(1) = vlevel(1)
+ !..check if subselection of levels
+ do k = 2, nk - 1
+ if (klevel(k + 1) /= klevel(k) - 1) then
+ met_params%manual_level_selection = .TRUE.
+ endif
+ end do
+ do k = 2, nk - 1
+ if (.NOT. met_params%manual_level_selection) then
+ ahalf(k) = alevel(k) + (alevel(k) - ahalf(k - 1))
+ bhalf(k) = blevel(k) + (blevel(k) - bhalf(k - 1))
+ vhalf(k) = ahalf(k)/standard_atmosphere + bhalf(k)
+ else
+ ahalf(k) = (alevel(k) + alevel(k + 1))*0.5
+ bhalf(k) = (blevel(k) + blevel(k + 1))*0.5
+ vhalf(k) = ahalf(k)/standard_atmosphere + bhalf(k)
+ end if
+ end do
+
+ ! Top half level is set to zero pressure
+ ahalf(nk) = 0.0
+ bhalf(nk) = 0.0
+ vhalf(nk) = 0.0
+ end subroutine
end module readfield_ncML
diff --git a/src/common/releasefile.f90 b/src/common/releasefile.f90
index b8dab7e0..dee1698d 100644
--- a/src/common/releasefile.f90
+++ b/src/common/releasefile.f90
@@ -23,6 +23,19 @@ module releasefileML
contains
+function should_skip_line(line)
+ character(len=*), intent(in) :: line
+ logical :: should_skip_line
+
+ if (line == "") then
+ should_skip_line = .true.
+ else if (line(1:1) == "*") then
+ should_skip_line = .true.
+ else
+ should_skip_line = .false.
+ endif
+end function
+
!> reading of input-files with hourly data (hours since run-start)
!> comment-rows start with #
!> hour height[upper_in_m] component release[kg/s]
@@ -42,7 +55,7 @@ subroutine releasefile(filename, release1)
integer :: ifd, ios, iexit, nlines
integer :: i, AllocateStatus
real :: hour, lasthour
- integer :: height
+ real :: height
integer :: ihour, iheight, icmp
real :: rel_s
character(32) :: comp
@@ -65,8 +78,7 @@ subroutine releasefile(filename, release1)
error stop 1
endif
-! header row
- nlines=0
+ nlines = 0
lasthour = -1
ihour = 0
inputlines: do
@@ -78,20 +90,14 @@ subroutine releasefile(filename, release1)
goto 11
endif
if (debugrelfile) write (error_unit,*) 'cinput (',nlines,'):',cinput
+ if (should_skip_line(cinput)) cycle inputlines
if (cinput == "end") exit inputlines
- if (cinput(1:1) == '*') cycle inputlines
read(cinput, *, err=12) hour, height, comp, rel_s
- if (hour < lasthour) then
- write (error_unit,*) 'hour must increase monotonic: ', &
- hour, ' < ', lasthour
- goto 12
- end if
if (hour > lasthour) then
! add new release timestep
lasthour = hour
- ihour = ihour + 1
if (.not.allocated(releases)) then
allocate(releases(1))
else
@@ -101,10 +107,24 @@ subroutine releasefile(filename, release1)
releases(1:size(tmp_release)) = tmp_release
deallocate(tmp_release)
endif
+ ihour = size(releases)
releases(ihour)%frelhour = hour
! make sure all initial release are undefined
releases(ihour)%relbqsec(:,:) = -1
- end if
+ else
+ ihour = 0
+ do i=1,size(releases)
+ if (releases(i)%frelhour == hour) then
+ ihour = i
+ exit
+ endif
+ enddo
+ if (ihour == 0) then
+ write(error_unit,*) "Unknown release hour: ", hour
+ goto 12
+ endif
+ endif
+
! find the component
icmp = 0
diff --git a/src/common/snap.F90 b/src/common/snap.F90
index 8ab7e541..66611d7e 100644
--- a/src/common/snap.F90
+++ b/src/common/snap.F90
@@ -32,10 +32,14 @@
! RANDOM.WALK.OFF
! BOUNDARY.LAYER.FULL.MIX.OFF
! BOUNDARY.LAYER.FULL.MIX.ON
-! DRY.DEPOSITION.OLD .......................... (default)
-! DRY.DEPOSITION.NEW
-! WET.DEPOSITION.OLD .......................... (default)
-! WET.DEPOSITION.NEW
+! DRY.DEPOSITION.OLD * deprecated
+! DRY.DEPOSITION.NEW * deprecated
+! DRY.DEPOSITION.SCHEME = emerson / new
+! DRY.DEPOSITION.SAVE
+! DRY.DEPOSITION.LARGEST_LANDFRACTION_FILE = "landclasses.nc"
+! WET.DEPOSITION.NEW ! deprecated
+! WET.DEPOSITION.SCHEME = Bartnicki ! (default) bartnicki-takemura
+! WET.DEPOSITION.SAVE ! Outputs wet scavenging rate (for 3D output only)
! TIME.STEP= 900.
! TIME.RELEASE.PROFILE.CONSTANT
! TIME.RELEASE.PROFILE.BOMB
@@ -74,7 +78,6 @@
! WET.DEP.OFF
! DRY.DEP.HEIGHT= 44.
! DRY.DEP.RATIO= 0.1
-! WET.DEP.RATIO= 0.2
! RADIOACTIVE.DECAY.ON
! RADIOACTIVE.DECAY.BOMB
! RADIOACTIVE.DECAY.OFF
@@ -119,8 +122,7 @@
! GRID.GPARAM = 3,-46.400002,-36.400002,0.10800000,0.10800000, 0.0000000, 65.000000
! * emep 1x1 deg lat lon
! * GRID.GPARAM = 2, -179.,-89.5,1.,1., 0., 0.
-! GRID.RUN= 88,1814, 1,1,1
-! * Norlam (sigma levels)
+! GRID.RUN ! deprecated
! DATA.SIGMA.LEVELS
! * Hirlam (eta levels)
! DATA.ETA.LEVELS
@@ -129,7 +131,7 @@
! * wind.surface = wind.10m (one 0 level first in list)
! LEVELS.INPUT= 14, 0,31,30,29,28,27,26,25,23,21,17,13,7,1
! FIELD.TYPE=felt|netcdf
-! * number of steps to skip in the beginning of a file
+! * number of steps to skip in the beginning of a file (obs: step 0 is not special in case of prognosis)
! FIELD.SPINUPSTEPS= 1
! FIELD.INPUT= arklam.dat
! FIELD.INPUT= feltlam.dat
@@ -139,7 +141,7 @@
! * increase output-resolution to factor*input_resolution
! FIELD.OUTPUT_RESOLUTION_FACTOR= 1
! FIELD.OUTTYPE=netcdf
-! FIELD.DAILY.OUTPUT.ON
+! FIELD.DAILY.OUTPUT.ON ! Output a file per day
! FIELD.USE_MODEL_WIND_INSTEAD_OF_10M= [.false.]/.true
! OUTPUT.COLUMN_MAX_CONC.ENABLE
! OUTPUT.COLUMN_MAX_CONC.DISABLE
@@ -186,12 +188,26 @@ PROGRAM bsnap
USE checkdomainML, only: check_in_domain
USE rwalkML, only: rwalk, rwalk_init
USE milibML, only: xyconvert
- use snapfldML, only: depwet, total_activity_lost_domain
+ use snapfldML, only: total_activity_lost_domain
USE forwrdML, only: forwrd, forwrd_init
- USE wetdep, only: wetdep2, wetdep2_init
- USE drydep, only: drydep1, drydep2
+ USE wetdepML, only: wetdep, wetdep_scheme, wetdep_scheme_t, &
+ WETDEP_SUBCLOUD_SCHEME_UNDEFINED, WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, WETDEP_INCLOUD_SCHEME_TAKEMURA, &
+ WETDEP_INCLOUD_SCHEME_UNDEFINED, &
+ operator(==), operator(/=), &
+ wetdep_init => init, wetdep_deinit => deinit
+#if defined(SNAP_EXPERIMENTAL)
+ USE wetdepML, only: WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL, &
+ wet_deposition_conventional_params => conventional_params, &
+ wet_deposition_RATM => RATM_params, &
+ WETDEP_INCLOUD_SCHEME_ROSELLE
+#endif
+ USE drydepml, only: drydep, drydep_scheme, &
+ DRYDEP_SCHEME_OLD, DRYDEP_SCHEME_NEW, DRYDEP_SCHEME_EMEP, &
+ DRYDEP_SCHEME_ZHANG, DRYDEP_SCHEME_EMERSON, DRYDEP_SCHEME_UNDEFINED, &
+ largest_landfraction_file, drydep_unload => unload
USE decayML, only: decay, decayDeps
- USE posintML, only: posint, posint_init
+ USE posintML, only: posint
USE bldpML, only: bldp
USE releaseML, only: release, releases, tpos_bomb, nrelheight, mprel, &
mplume, nplume, iplume, npart, mpart, release_t
@@ -239,7 +255,7 @@ PROGRAM bsnap
integer :: k, ierror, i, n
integer, allocatable :: klevel_manual(:)
integer :: ih
- integer :: idrydep = 0, wetdep_version = 0, idecay
+ integer :: idecay
integer :: ntimefo
integer :: nsteph, nstep, nstepr
integer :: ihread, isteph, lstepr, iendrel, istep, nhleft
@@ -247,6 +263,7 @@ PROGRAM bsnap
integer :: ihdiff, ifldout, idailyout = 0, ihour, split_particle_after_step, split_particle_hours
integer :: date_time(8)
logical :: warning = .false.
+ integer :: npartmax
!> tstep: timestep in seconds
real :: tstep = 900
real :: rmlimit = -1.0, rnhrel, tf1, tf2, tnow, tnext
@@ -256,13 +273,8 @@ PROGRAM bsnap
integer :: ntprof
type(duration_t) :: dur
logical :: out_of_domain
-! ipcount(mdefcomp, nk)
-! integer, dimension(:,:), allocatable:: ipcount
-! npcount(nk)
-! integer, dimension(:), allocatable:: npcount
-! b_start
+
real :: mhmin, mhmax ! minimum and maximum of mixing height
-! b_end
!> Information for reading from a releasefile
type(release_t) :: release1
@@ -486,8 +498,11 @@ PROGRAM bsnap
write (iulog, *) 'mprel: ', mprel
write (iulog, *) 'ifltim: ', ifltim
write (iulog, *) 'irwalk: ', use_random_walk
- write (iulog, *) 'idrydep: ', idrydep
- write (iulog, *) 'wetdep_version: ', wetdep_version
+ write (iulog, *) 'drydep_scheme: ', drydep_scheme
+ write (iulog, *) 'wetdep_scheme: subcloud scheme: ', wetdep_scheme%subcloud%description
+ write (iulog, *) 'wetdep_scheme: incloud scheme: ', wetdep_scheme%incloud%description
+ write (iulog, *) 'wetdep_scheme: use vertical: ', wetdep_scheme%use_vertical
+ write (iulog, *) 'wetdep_scheme: use cloudfraction: ', wetdep_scheme%use_cloudfraction
write (iulog, *) 'idecay: ', idecay
write (iulog, *) 'rmlimit: ', rmlimit
write (iulog, *) 'ndefcomp:', size(def_comp)
@@ -508,7 +523,6 @@ PROGRAM bsnap
write (iulog, *) ' drydephgt: ', def_comp(m)%drydephgt
write (iulog, *) ' drydeprat: ', def_comp(m)%drydeprat
write (iulog, *) ' kwetdep: ', def_comp(m)%kwetdep
- write (iulog, *) ' wetdeprat: ', def_comp(m)%wetdeprat
write (iulog, *) ' kdecay: ', def_comp(m)%kdecay
write (iulog, *) ' halftime: ', def_comp(m)%halftime
write (iulog, *) ' decayrate: ', def_comp(m)%decayrate
@@ -560,6 +574,7 @@ PROGRAM bsnap
mhmax = -10.0
! b_end
+
! reset readfield_nc (eventually, traj will rerun this loop)
if (ftype == "netcdf") then
call readfield_nc(-1, nhrun < 0, time_start, nhfmin, nhfmax, &
@@ -627,6 +642,7 @@ PROGRAM bsnap
! start time loop
itimei = time_start
+ npartmax = 0
time_loop: do istep = 0, nstep
call timeloop_timer%start()
write (iulog, *) 'istep,nplume,npart: ', istep, nplume, npart
@@ -683,6 +699,7 @@ PROGRAM bsnap
!..release one plume of particles
call release(istep, nsteph, tf1, tf2, tnow, ierror)
+ npartmax = max(npartmax, npart)
if (ierror == 0) then
lstepr = istep
@@ -701,9 +718,7 @@ PROGRAM bsnap
if (idecay == 1) call decayDeps(tstep)
! prepare particle functions once before loop
if (init) then
- ! setting particle-number to 0 means init
- call posint_init()
- if (wetdep_version == 2) call wetdep2_init(tstep)
+ call wetdep_init(tstep)
call forwrd_init()
if (use_random_walk) call rwalk_init(tstep)
init = .FALSE.
@@ -729,12 +744,11 @@ PROGRAM bsnap
!..radioactive decay
if (idecay == 1) call decay(pdata(np))
- !..dry deposition (1=old, 2=new version)
- if (idrydep == 1) call drydep1(pdata(np))
- if (idrydep == 2) call drydep2(tstep, pdata(np))
+ !..dry deposition
+ call drydep(tstep, pdata(np))
- !..wet deposition (1=old, 2=new version)
- if (wetdep_version == 2) call wetdep2(depwet, tstep, pdata(np), pextra)
+ !..wet deposition
+ call wetdep(tstep, pdata(np), pextra)
!..move all particles forward, save u and v to pextra
call forwrd(tf1, tf2, tnow, tstep, pdata(np), pextra)
@@ -763,6 +777,7 @@ PROGRAM bsnap
if (split_particle_after_step > 0) then
call split_particles(split_particle_after_step)
end if
+ npartmax = max(npartmax, npart)
!$OMP PARALLEL DO REDUCTION(max : mhmax) REDUCTION(min : mhmin)
do n = 1, npart
@@ -847,6 +862,8 @@ PROGRAM bsnap
call snap_error_exit(iulog)
end if
+ write(error_unit, *) 'npart: Used a maximum of ', npartmax, ' out of available ', mpart
+ write(output_unit, *) 'npart: Used a maximum of ', npartmax, ' out of available ', mpart
! b_240311
write (error_unit, *)
@@ -866,6 +883,8 @@ PROGRAM bsnap
call fldout_unload()
! deallocate all fields
+ call wetdep_deinit()
+ call drydep_unload()
CALL deAllocateFields()
close (iulog)
@@ -909,7 +928,7 @@ subroutine read_inputfile(snapinput_unit)
TIME_PROFILE_UNDEFINED
use snapfimexML, only: parse_interpolator
use snapgrdml, only: compute_column_max_conc, compute_aircraft_doserate, aircraft_doserate_threshold, &
- output_column
+ output_column, output_vd, output_vd_debug
use init_random_seedML, only: extra_seed
use fldout_ncML, only: surface_layer_is_lowest_level, surface_height_m
@@ -1109,34 +1128,200 @@ subroutine read_inputfile(snapinput_unit)
endif
case ('dry.deposition.old')
!..dry.deposition.old
- if (idrydep /= 0 .AND. idrydep /= 1) goto 12
- idrydep = 1
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_OLD) goto 12
+ write(error_unit,*) "dry.deposition.old is deprecated, use dry.deposition.scheme=old"
+ warning = .true.
+ drydep_scheme = DRYDEP_SCHEME_OLD
case ('dry.deposition.new')
!..dry.deposition.new
- if (idrydep /= 0 .AND. idrydep /= 2) goto 12
- idrydep = 2
- case ('wet.deposition.old')
- write (error_unit, *) "This option is deprecated and removed"
- goto 12
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_NEW) goto 12
+ write(error_unit,*) "dry.deposition.new is deprecated, use dry.deposition.scheme=new"
+ warning = .true.
+ drydep_scheme = DRYDEP_SCHEME_NEW
+ case ('dry.deposition.scheme')
+ if (.not. has_value) goto 12
+ call to_lowercase(cinput(pname_start:pname_end))
+ select case (cinput(pname_start:pname_end))
+ case ('old')
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_OLD) goto 12
+ drydep_scheme = DRYDEP_SCHEME_OLD
+ case ('new')
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_NEW) goto 12
+ drydep_scheme = DRYDEP_SCHEME_NEW
+#if defined(SNAP_EXPERIMENTAL)
+ case ('emep')
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_EMEP) goto 12
+ drydep_scheme = DRYDEP_SCHEME_EMEP
+ case ('zhang')
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_ZHANG) goto 12
+ drydep_scheme = DRYDEP_SCHEME_ZHANG
+#endif
+ case ('emerson')
+ if (drydep_scheme /= 0 .AND. drydep_scheme /= DRYDEP_SCHEME_EMERSON) goto 12
+ drydep_scheme = DRYDEP_SCHEME_EMERSON
+ case default
+ write(error_unit, *) "Scheme ", cinput(pname_start:pname_end), " is unknown"
+ goto 12
+ end select
+ case ('dry.deposition.save')
+ if (drydep_scheme /= DRYDEP_SCHEME_EMEP .and. &
+ drydep_scheme /= DRYDEP_SCHEME_ZHANG .and. &
+ drydep_scheme /= DRYDEP_SCHEME_EMERSON) then
+ write(error_unit, *) "The drydep scheme is not set to a compatible value, ignoring"
+ else
+ output_vd = .true.
+ endif
+ if (has_value) goto 12
+ case ('dry.deposition.save.debug')
+ if (.not.output_vd) then
+ write(error_unit, *) "dry.deposition.save must be set"
+ goto 12
+ endif
+ warning = .true.
+ write(error_unit, *) "This option is a hack, only use for a single component"
+ output_vd_debug = .true.
+ case ('dry.deposition.largest_landfraction_file')
+ if (.not.has_value) goto 12
+ largest_landfraction_file = cinput(pname_start:pname_end)
case ('wet.deposition.new')
!..wet.deposition.new
- write (error_unit, *) "Deprecated, please use wet.deposition.version = 2"
+ write (error_unit, *) "Deprecated, please use wet.deposition.scheme = Bartnicki"
warning = .true.
- if (wetdep_version /= 0) then
- write (error_unit, *) "already set"
+ if (wetdep_scheme%subcloud /= WETDEP_SUBCLOUD_SCHEME_UNDEFINED .or. &
+ wetdep_scheme%incloud /= WETDEP_INCLOUD_SCHEME_UNDEFINED) then
+ write (error_unit, *) "wet deposition already set"
goto 12
endif
- wetdep_version = 2
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .false., .false.)
case ('wet.deposition.version')
- if (wetdep_version /= 0) then
- write (error_unit, *) "already set"
+ write (error_unit, *) "Deprecated, please use wet.deposition.scheme = Bartnicki"
+ warning = .true.
+ if (wetdep_scheme%subcloud /= WETDEP_SUBCLOUD_SCHEME_UNDEFINED .or. &
+ wetdep_scheme%incloud /= WETDEP_INCLOUD_SCHEME_UNDEFINED) then
+ write (error_unit, *) "wet deposition already set"
goto 12
endif
if (.not. has_value) then
write (error_unit, *) "expected a keyword"
goto 12
endif
- read (cinput(pname_start:pname_end), *, err=12) wetdep_version
+ block
+ integer :: scheme_number
+ read (cinput(pname_start:pname_end), *, err=12) scheme_number
+ if (scheme_number /= 2) goto 12
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .false., &
+ .false.)
+ end block
+#if defined(SNAP_EXPERIMENTAL)
+ case ('wet.deposition.conventional.a')
+ if (wet_deposition_conventional_params%A /= 0.0) then
+ write(error_unit,*) "wet deposition parameter already set"
+ goto 12
+ endif
+ if (.not.has_value) goto 12
+ read(cinput(pname_start:pname_end), *) wet_deposition_conventional_params%A
+ case ('wet.deposition.conventional.b')
+ if (wet_deposition_conventional_params%B /= 0.0) then
+ write(error_unit,*) "wet deposition parameter already set"
+ goto 12
+ endif
+ if (.not.has_value) goto 12
+ read(cinput(pname_start:pname_end), *) wet_deposition_conventional_params%B
+#endif
+ case ('wet.deposition.scheme')
+ if (.not.has_value) goto 12
+ if (wetdep_scheme%subcloud /= WETDEP_SUBCLOUD_SCHEME_UNDEFINED .or. &
+ wetdep_scheme%incloud /= WETDEP_INCLOUD_SCHEME_UNDEFINED) then
+ write (error_unit, *) "wet deposition already set"
+ goto 12
+ endif
+ call to_lowercase(cinput(pname_start:pname_end))
+ select case(cinput(pname_start:pname_end))
+ case("bartnicki")
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .false., .false. &
+ )
+ case("bartnicki-takemura")
+ met_params%use_3d_precip = .true.
+ met_params%use_ccf = .true.
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_TAKEMURA, &
+ .true., .true. &
+ )
+ case("bartnicki-vertical")
+ met_params%use_3d_precip = .true.
+ met_params%use_ccf = .true.
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .true., .true. &
+ )
+#if defined(SNAP_EXPERIMENTAL)
+ case("conventional")
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .false., .false. &
+ )
+ case("bartnicki-roselle")
+ met_params%use_3d_precip = .true.
+ met_params%use_ccf = .true.
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_ROSELLE, &
+ .true., .true. &
+ )
+ case("ratm-roselle")
+ if (wet_deposition_conventional_params%A /= 0.0 .or. &
+ wet_deposition_conventional_params%B /= 0.0) then
+ write(error_unit,*) "wet deposition parameter already set"
+ goto 12
+ endif
+ wet_deposition_conventional_params = wet_deposition_RATM
+ met_params%use_3d_precip = .true.
+ met_params%use_ccf = .true.
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL, &
+ WETDEP_INCLOUD_SCHEME_ROSELLE, &
+ .true., .true. &
+ )
+ case("ratm-takemura")
+ if (wet_deposition_conventional_params%A /= 0.0 .or. &
+ wet_deposition_conventional_params%B /= 0.0) then
+ write(error_unit,*) "wet deposition parameter already set"
+ goto 12
+ endif
+ wet_deposition_conventional_params = wet_deposition_RATM
+ met_params%use_3d_precip = .true.
+ met_params%use_ccf = .true.
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL, &
+ WETDEP_INCLOUD_SCHEME_TAKEMURA, &
+ .true., .true. &
+ )
+#endif
+ case default
+ write(error_unit,*) "Unknown scheme ", cinput(pname_start:pname_end)
+ goto 12
+ end select
+ case ('wet.deposition.save')
+ if (has_value) goto 12
+ if (.not.wetdep_scheme%use_vertical) then
+ write(error_unit,*) "wet.depositon.save is only allowed when 3D wetdep scheme is used"
+ endif
+ block
+ use fldout_ncml, only: output_wetdeprate
+ output_wetdeprate = .true.
+ end block
case ('time.step')
!..time.step=
if (.not. has_value) goto 12
@@ -1357,11 +1542,8 @@ subroutine read_inputfile(snapinput_unit)
if (d_comp%drydeprat >= 0.) goto 12
read (cinput(pname_start:pname_end), *, err=12) d_comp%drydeprat
case ('wet.dep.ratio')
- !..wet.dep.ratio=
- if (.not. has_value) goto 12
- if (.not. associated(d_comp)) goto 12
- if (d_comp%wetdeprat >= 0.) goto 12
- read (cinput(pname_start:pname_end), *, err=12) d_comp%wetdeprat
+ write (error_unit, *) "wet.dep.ratio is no longer used"
+ warning = .true.
case ('radioactive.decay.on')
!..radioactive.decay.on
if (.not. associated(d_comp)) goto 12
@@ -1806,10 +1988,11 @@ subroutine conform_input(ierror)
use find_parameter, only: detect_gridparams, get_klevel
#if defined(FIMEX)
use find_parameters_fi, only: detect_gridparams_fi
+ use readfield_fiML, only: read_largest_landfraction
#endif
integer, intent(out) :: ierror
- integer :: i1, i2
+ integer :: i1
logical :: error_release_profile
@@ -1944,10 +2127,6 @@ subroutine conform_input(ierror)
end if
do m = 1, size(def_comp) - 1
- if (def_comp(m)%idcomp < 1) then
- write (error_unit, *) 'Component has no field identification: ', &
- trim(def_comp(m)%compname)
- end if
do i = m + 1, size(def_comp)
if (def_comp(m)%compname == def_comp(i)%compname) then
write (error_unit, *) 'Component defined more than once: ', &
@@ -2035,22 +2214,25 @@ subroutine conform_input(ierror)
end if
end do
- if (idrydep == 0) idrydep = 1
- if (wetdep_version == 0) then ! Set default wetdep version
- wetdep_version = 2
- endif
- if (wetdep_version /= 2) then
- write (error_unit, *) "Unknown wet deposition version"
- ierror = 1
+ if (drydep_scheme == DRYDEP_SCHEME_UNDEFINED) drydep_scheme = DRYDEP_SCHEME_OLD
+
+ ! Set default wetdep schemes
+ if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_UNDEFINED .and. &
+ wetdep_scheme%incloud == WETDEP_INCLOUD_SCHEME_UNDEFINED) then
+ wetdep_scheme = wetdep_scheme_t( &
+ WETDEP_SUBCLOUD_SCHEME_BARTNICKI, &
+ WETDEP_INCLOUD_SCHEME_NONE, &
+ .false., .false. &
+ )
endif
+
i1 = 0
- i2 = 0
idecay = 0
do n = 1, ncomp
m = run_comp(n)%to_defined
if (m == 0) cycle
- if (idrydep == 1 .AND. def_comp(m)%kdrydep == 1) then
+ if (drydep_scheme == DRYDEP_SCHEME_OLD .AND. def_comp(m)%kdrydep == 1) then
if (def_comp(m)%drydeprat > 0. .AND. def_comp(m)%drydephgt > 0.) then
i1 = i1 + 1
else
@@ -2058,7 +2240,7 @@ subroutine conform_input(ierror)
def_comp(m)%drydeprat, def_comp(m)%drydephgt
ierror = 1
end if
- elseif (idrydep == 2 .AND. def_comp(m)%kdrydep == 1) then
+ elseif (drydep_scheme == DRYDEP_SCHEME_NEW .AND. def_comp(m)%kdrydep == 1) then
if (def_comp(m)%grav_type == 1 .AND. def_comp(m)%gravityms > 0.) then
i1 = i1 + 1
elseif (def_comp(m)%grav_type == 2) then
@@ -2068,12 +2250,21 @@ subroutine conform_input(ierror)
def_comp(m)%gravityms
ierror = 1
end if
+ elseif (((drydep_scheme == DRYDEP_SCHEME_EMEP) .or. &
+ (drydep_scheme == DRYDEP_SCHEME_EMERSON) .or. &
+ (drydep_scheme == DRYDEP_SCHEME_ZHANG)) .and. &
+ def_comp(m)%kdrydep == 1) then
+ ! Check if component has the necessary definitions to compute
+ ! the dry deposition
+ if (.true.) then
+ i1 = i1 + 1
+ else
+ write (error_unit, *) 'Dry deposition error'
+ end if
end if
- if (wetdep_version == 2 .AND. def_comp(m)%kwetdep == 1) then
- if (def_comp(m)%radiusmym > 0.) then
- i2 = i2 + 1
- else
+ if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_BARTNICKI .AND. def_comp(m)%kwetdep == 1) then
+ if (def_comp(m)%radiusmym <= 0.) then
write (error_unit, *) 'Wet deposition error. radius: ', &
def_comp(m)%radiusmym
ierror = 1
@@ -2090,7 +2281,14 @@ subroutine conform_input(ierror)
end if
end do
- if (i1 == 0) idrydep = 0
+ if (i1 == 0) drydep_scheme = DRYDEP_SCHEME_UNDEFINED
+ if (drydep_scheme /= DRYDEP_SCHEME_UNDEFINED .and. largest_landfraction_file /= "not set") then
+#if defined(FIMEX)
+ call read_largest_landfraction(largest_landfraction_file)
+#else
+ error stop "Reading of largest landfraction requires fimex support"
+#endif
+ endif
if (itotcomp == 1 .AND. ncomp == 1) itotcomp = 0
diff --git a/src/common/snap.mk b/src/common/snap.mk
index 5284db66..f82cf0b2 100644
--- a/src/common/snap.mk
+++ b/src/common/snap.mk
@@ -36,8 +36,8 @@ snap.o: ../common/snap.F90 $(MODELOBJ)
#--------------------------------
fimex.o: ../common/fimex.f90
- ${F77} -c $(F77FLAGS) $(INCLUDES) $<
-allocateFields.o: ../common/allocateFields.f90 particleML.o snapparML.o snapfldML.o snapfilML.o snapgrdML.o release.o snapdimML.o
+ ${F77} -c $(F77FLAGS) $(INCLUDES) -fno-module-private $<
+allocateFields.o: ../common/allocateFields.f90 particleML.o snapparML.o snapfldML.o snapfilML.o snapgrdML.o release.o snapdimML.o snapmetML.o
${F77} -c $(F77FLAGS) $(INCLUDES) $<
array_utils.o: ../common/array_utils.f90
${F77} -c $(F77FLAGS) $(INCLUDES) $<
@@ -77,9 +77,9 @@ epinterp.o: ../common/epinterp.f90
${F77} -c ${F77FLAGS} $(INCLUDES) $<
om2edot.o: ../common/om2edot.f90 snapgrdML.o snapfldML.o snapdimML.o snapdebugML.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
-readfield_nc.o: ../common/readfield_nc.f90 particleML.o snapfilML.o snapgrdML.o snapmetML.o snaptabML.o snapdebugML.o snapdimML.o om2edot.o ftest.o milibML.o snapfldML.o datetime.o
+readfield_nc.o: ../common/readfield_nc.f90 particleML.o snapfilML.o snapgrdML.o snapmetML.o snaptabML.o snapdebugML.o snapdimML.o om2edot.o ftest.o milibML.o snapfldML.o datetime.o drydep.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
-readfield_fi.o: ../common/readfield_fi.f90 snapfimexML.o particleML.o snapfilML.o snapgrdML.o snapmetML.o snaptabML.o snapdebugML.o snapdimML.o om2edot.o ftest.o milibML.o fimex.o datetime.o readfield_nc.o utils.o
+readfield_fi.o: ../common/readfield_fi.f90 snapfimexML.o particleML.o snapfilML.o snapgrdML.o snapmetML.o snaptabML.o snapdebugML.o snapdimML.o om2edot.o ftest.o milibML.o fimex.o datetime.o readfield_nc.o utils.o drydep.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
filesort_nc.o: ../common/filesort_nc.f90 dateCalc.o snapfilML.o snapdimML.o snapgrdML.o snapfldML.o snapmetML.o snapdebugML.o readfield_nc.o datetime.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
@@ -113,7 +113,7 @@ split_particles.o: ../common/split_particles.f90 snapparML.o release.o snapdebug
${F77} -c ${F77FLAGS} $<
utils.o: ../common/utils.f90
${F77} -c ${F77FLAGS} $(INCLUDES) $<
-vgravtables.o: ../common/vgravtables.f90 snapparML.o snapdimML.o
+vgravtables.o: ../common/vgravtables.f90 snapparML.o snapdimML.o drydep.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
wetdep.o: ../common/wetdep.f90 particleML.o snapgrdML.o snapfldML.o snapparML.o snaptabML.o snapdimML.o snapdebugML.o
${F77} -c ${F77FLAGS} $(INCLUDES) $<
diff --git a/src/common/snapdimML.f90 b/src/common/snapdimML.f90
index 3a96dea0..00115407 100644
--- a/src/common/snapdimML.f90
+++ b/src/common/snapdimML.f90
@@ -60,10 +60,15 @@ module snapdimML
!> low_resolution input grid position
public :: lres_pos
+!> translate a field in normal resolution to high output resolution
+ interface hres_field
+ module procedure :: hres_field_real32, hres_field_int8, hres_field_real64_real32
+ end interface
+
contains
!> translate a field in normal resolution to high output resolution
- subroutine hres_field(field, field_hres, bilinear)
+ subroutine hres_field_real32(field, field_hres, bilinear)
USE iso_fortran_env, only: real32
real(kind=real32), intent(in) :: field(:,:)
real(kind=real32), intent(inout) :: field_hres(:,:)
@@ -111,11 +116,85 @@ subroutine hres_field(field, field_hres, bilinear)
end do
end do
end if
- end subroutine hres_field
+ end subroutine
+
+!> translate a field in normal resolution to high output resolution
+ subroutine hres_field_real64_real32(field, field_hres, bilinear)
+ USE iso_fortran_env, only: real32, real64
+ real(kind=real64), intent(in) :: field(:,:)
+ real(kind=real32), intent(inout) :: field_hres(:,:)
+ logical, optional, intent(in) :: bilinear
+ logical :: is_bilinear
+ integer :: i, j, k, l, or_2
+ real dd, dx, dy, c1, c2, c3, c4
+
+ is_bilinear = .false.
+ if (present(bilinear)) is_bilinear = bilinear
+ if (output_resolution_factor == 1) is_bilinear = .false.
+
+ ! nearest neighbor, even if bilinear, to fix the borders
+ do j = 1, ny
+ do l = 1, output_resolution_factor
+ do i = 1, nx
+ do k = 1, output_resolution_factor
+ field_hres(output_resolution_factor*(i-1)+k, output_resolution_factor*(j-1)+l) = field(i,j)
+ end do
+ end do
+ end do
+ end do
+
+ if (is_bilinear) then
+ or_2 = int(output_resolution_factor / 2.)
+ dd = 1./output_resolution_factor
+ do j = 1, ny-1
+ do l = 1, output_resolution_factor
+ do i = 1, nx-1
+ do k = 1, output_resolution_factor
+ ! this will partly extrapolate to left/bottom
+ ! but otherwise too many corner-cases needed
+ dx=(k-1-or_2)*dd
+ dy=(l-1-or_2)*dd
+ c1=(1.-dy)*(1.-dx)
+ c2=(1.-dy)*dx
+ c3=dy*(1.-dx)
+ c4=dy*dx
+ ! go here from -5 to 4 (for case resolution-factor 10)
+ field_hres(output_resolution_factor*i+k-or_2, output_resolution_factor*j+l-or_2) = &
+ c1 * field(i, j) + c2 * field(i+1, j) + &
+ c3 * field(i, j+1) + c4 * field(i+1, j+1)
+ end do
+ end do
+ end do
+ end do
+ end if
+ end subroutine
+
+ subroutine hres_field_int8(field, field_hres)
+ USE iso_fortran_env, only: int8
+ integer(kind=int8), intent(in) :: field(:,:)
+ integer(kind=int8), allocatable, intent(out) :: field_hres(:,:)
+
+ integer :: newshape(2)
+
+ integer :: i, j, k, l
+
+ newshape(:) = shape(field)*output_resolution_factor
+ allocate(field_hres(newshape(1), newshape(2)))
+
+ do j = 1, ny
+ do l = 1, output_resolution_factor
+ do i = 1, nx
+ do k = 1, output_resolution_factor
+ field_hres(output_resolution_factor*(i-1)+k, output_resolution_factor*(j-1)+l) = field(i,j)
+ end do
+ end do
+ end do
+ end do
+ end subroutine
!> translate a x or y position in the input-grid to the
!> high_resolution output grid position
- integer function hres_pos(lres_pos)
+ pure integer function hres_pos(lres_pos)
USE iso_fortran_env, only: real64
real(kind=real64), intent(in) :: lres_pos
! convert to 0.5-starting position (cell 1 from [0.5,1.5[, extend to new range, convert to 1-start
@@ -123,8 +202,8 @@ integer function hres_pos(lres_pos)
end function hres_pos
!> translate a x or y position in the output-grid to the
-!> lwo_resolution input grid position
- integer function lres_pos(hres_pos)
+!> low_resolution input grid position
+ pure integer function lres_pos(hres_pos)
USE iso_fortran_env, only: real64
integer, intent(in) :: hres_pos
! convert to 0-starting positions, extend to new range, convert to 1-start
diff --git a/src/common/snapfldML.f90 b/src/common/snapfldML.f90
index 46c39c3c..970d05b2 100644
--- a/src/common/snapfldML.f90
+++ b/src/common/snapfldML.f90
@@ -53,6 +53,16 @@ module snapfldML
!> hourly precipitation intensity (mm/hour)
real(kind=real32), allocatable, save, public :: precip(:,:)
+!> instant precipitation intensity in three dimensions (mm/hour)
+ real(kind=real32), allocatable, save, public :: precip3d(:,:,:)
+!> Cloud water content (mm)
+ real(kind=real32), allocatable, save, public :: cw3d(:,:,:)
+
+!> Cloud cover fraction
+ real(kind=real32), allocatable, save, public :: cloud_cover(:,:,:)
+
+!> wet scavenging rate
+ real(kind=real32), allocatable, save, public :: wscav(:,:,:,:)
!> surface pressure (time step 1)
real(kind=real32), allocatable, save, public :: ps1(:,:)
@@ -162,4 +172,9 @@ module snapfldML
! > Activity lost through exiting rmlimit, maxage
real(kind=real64), allocatable, save, public :: total_activity_lost_other(:)
+ !> Deposition velocity on the grid per species
+ real, allocatable, save, public :: vd_dep(:, :, :)
+
+ real, allocatable, save, public :: xflux(:, :), yflux(:, :), hflux(:, :), z0(:, :), leaf_area_index(:, :), t2m(:, :)
+ real(real64), allocatable, save, public :: roa(:,:), ustar(:,:), monin_l(:,:), raero(:,:), vs(:,:), rs(:,:)
end module snapfldML
diff --git a/src/common/snapgrdML.f90 b/src/common/snapgrdML.f90
index f21948a7..5bd1a204 100644
--- a/src/common/snapgrdML.f90
+++ b/src/common/snapgrdML.f90
@@ -126,4 +126,10 @@ module snapgrdML
!> Output concentration in each column
logical, save, public :: output_column = .false.
+
+!> Save dry deposition velocity in output file
+ logical, save, public :: output_vd = .false.
+!> Extra debug flags for dry deposition
+ logical, save, public :: output_vd_debug = .false.
+
end module snapgrdML
diff --git a/src/common/snapmetML.f90 b/src/common/snapmetML.f90
index b4f4026e..388a80c1 100644
--- a/src/common/snapmetML.f90
+++ b/src/common/snapmetML.f90
@@ -45,6 +45,13 @@ module snapmetML
character(len=80) :: precconvrt = ''
character(len=80) :: total_column_rain = ''
+ character(len=80) :: t2m = ''
+ character(len=80) :: yflux = ''
+ character(len=80) :: xflux = ''
+ character(len=80) :: hflux = ''
+ character(len=80) :: z0 = ''
+ character(len=80) :: leaf_area_index = ''
+
! flags when reading the data
logical :: temp_is_abs = .false.
logical :: has_dummy_dim = .false.
@@ -54,6 +61,20 @@ module snapmetML
!> Use lowest level in #xwindv/#ywindv in place of
!> #xwind10mv/#ywind10mv
logical :: use_model_wind_for_10m = .false.
+ !> Wet deposition: Use 3D precip
+ logical :: use_3d_precip = .false.
+ !> Wet deposition: Use cloud cover fraction
+ logical :: use_ccf = .false.
+
+ !> Cloud water (3D)
+ character(len=80) :: mass_fraction_rain_in_air = ''
+ character(len=80) :: mass_fraction_graupel_in_air = ''
+ character(len=80) :: mass_fraction_snow_in_air = ''
+ !> Precip (3D)
+ character(len=80) :: mass_fraction_cloud_condensed_water_in_air = ''
+ character(len=80) :: mass_fraction_cloud_ice_in_air = ''
+ !> Cloud fraction (3D)
+ character(len=80) :: cloud_fraction = ''
end type
type(met_params_t), save, public :: met_params = met_params_t()
@@ -70,8 +91,15 @@ module snapmetML
character(len=*), parameter, public :: precip_units_fallback = 'mm'
character(len=*), parameter, public :: temp_units = 'K'
+ 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 :: leaf_area_index_units = '1'
+
+ character(len=*), parameter, public :: cloud_fraction_units = '%'
+ character(len=*), parameter, public :: mass_fraction_units = 'kg/kg'
- public init_meteo_params, requires_precip_deaccumulation
+ public :: init_meteo_params, requires_precip_deaccumulation
contains
@@ -200,6 +228,18 @@ subroutine init_meteo_params(nctype, ierr)
met_params%precaccumv = 'precipitation_amount_acc'
met_params%precstrativrt = ''
met_params%precconvrt = ''
+
+ met_params%t2m = 'air_temperature_2m'
+ met_params%xflux = 'downward_northward_momentum_flux_in_air'
+ met_params%yflux = 'downward_eastward_momentum_flux_in_air'
+ met_params%z0 = 'surface_roughness_length'
+ met_params%hflux = 'integral_of_surface_downward_sensible_heat_flux_wrt_time'
+ met_params%leaf_area_index = 'leaf_area_index'
+
+ met_params%mass_fraction_cloud_condensed_water_in_air = "cloudwater"
+ met_params%mass_fraction_cloud_ice_in_air = "cloudice"
+
+ met_params%cloud_fraction = "3D_cloudcover"
!..get grid parameters from field identification
case('dmi_eps')
met_params%has_dummy_dim = .true.
@@ -262,6 +302,15 @@ subroutine init_meteo_params(nctype, ierr)
!.. non accumulated precipitation rates in m/s
met_params%precstrativrt = 'large_scale_precipitations'
met_params%precconvrt = 'convective_precipitations'
+
+ 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"
+
+ met_params%mass_fraction_cloud_condensed_water_in_air = "mass_fraction_of_cloud_condensed_water_in_air_ml"
+ met_params%mass_fraction_cloud_ice_in_air = "mass_fraction_of_cloud_ice_in_air_ml"
+
+ met_params%cloud_fraction = "cloud_area_fraction_ml"
!..get grid parameters from field identification
! set as long as sortfield still is called
case('gfs_grib_filter_fimex')
diff --git a/src/common/snapparML.f90 b/src/common/snapparML.f90
index 5b66d935..9f700dc9 100644
--- a/src/common/snapparML.f90
+++ b/src/common/snapparML.f90
@@ -67,8 +67,6 @@ module snapparML
!> radioactive decay (rate)
real :: decayrate = -1.0
-!> wet deposition rate
- real :: wetdeprat = -1.0
!> for each component: 0=wet deposition off 1=wet dep. on
integer :: kwetdep = -1
diff --git a/src/common/snaptabML.f90 b/src/common/snaptabML.f90
index 900528ba..47683a90 100644
--- a/src/common/snaptabML.f90
+++ b/src/common/snaptabML.f90
@@ -31,7 +31,7 @@ module snaptabML
real, parameter, public :: g=9.81, r=287.0, cp=1004.0
real, parameter, public :: rcp = r/cp
-
+!> Units of hPa
real, parameter, public :: standard_atmosphere = 1013.25
!> Height of surface layer
diff --git a/src/common/vgravtables.f90 b/src/common/vgravtables.f90
index c0bc4acf..36f686f4 100644
--- a/src/common/vgravtables.f90
+++ b/src/common/vgravtables.f90
@@ -16,6 +16,7 @@
! along with this program. If not, see .
module vgravtablesML
+ use drydepml, only: drydep_scheme, DRYDEP_SCHEME_EMEP, DRYDEP_SCHEME_ZHANG, DRYDEP_SCHEME_EMERSON
implicit none
private
@@ -26,7 +27,7 @@ module vgravtablesML
!> table of gravity in m/s
!>
-!> (temperature as first index, pressure second)
+!> vgtable(temperature,pressure,component)
real, save, allocatable, public :: vgtable(:,:,:)
real, parameter, public :: tincrvg = 200.0/float(numtempvg - 1)
@@ -34,21 +35,31 @@ module vgravtablesML
real, parameter, public :: pincrvg = 1200./float(numpresvg-1)
real, parameter, public :: pbasevg = 0. - pincrvg
- public vgravtables
+ public :: vgravtables
contains
!> program for calculating the gravitational settling velocities
-!> for small and large particles (outside the Stokes low)
+!> for small and large particles (outside the Stokes law)
subroutine vgravtables
+ USE ISO_FORTRAN_ENV, only: real64
USE snapparML, only: ncomp, run_comp, def_comp
-
- real :: t ! absolute temperature (K)
- real :: dp ! particle size um
- real :: rp ! particle density in g/m3
- real :: vg ! gravitational setling in cm/s (Stokes low)
- real :: vgmod ! gravitational setling in cm/s after iteration
- real :: p ! atmospheric pressure
+ USE drydepml, only: gravitational_settling
+
+ real, parameter :: R = 287.05
+ !> absolute temperature (K)
+ real :: t
+ !> particle size
+ real :: diam_part
+ !> particle density
+ real :: rho_part
+ !> gravitational setling
+ real :: vg
+ !> gravitational setling after iteration
+ real :: vgmod
+ !> atmospheric pressure
+ real :: p
+ real(real64) :: roa
integer :: n,m,ip,it
!---------------------------------------
@@ -62,27 +73,46 @@ subroutine vgravtables
if (def_comp(m)%grav_type /= 2) then ! Not using gravity table
cycle do_comp
endif
- ! radius to diameter
- dp= 2. * def_comp(m)%radiusmym
- rp= def_comp(m)%densitygcm3
- do ip=1,numpresvg
+ select case(drydep_scheme)
+ case (DRYDEP_SCHEME_EMEP,DRYDEP_SCHEME_ZHANG,DRYDEP_SCHEME_EMERSON)
+ ! expected kg/m3
+ rho_part = def_comp(m)%densitygcm3 / 1000.0
+ ! expected m
+ diam_part = 2 * def_comp(m)%radiusmym / 1e6
+
+ do ip=1,numpresvg
+ ! Expecting pascal
+ p = max(1.0, pbasevg + ip*pincrvg) / 100.0
+ do it=1,numtempvg
+ t = tbasevg + it*tincrvg
+ roa = p / (real(t, kind=real64) * R)
+ vgtable(it, ip, n) = gravitational_settling(roa, real(diam_part, kind=real64), real(rho_part, kind=real64))
+ end do
+ end do
+ case default
+ ! radius to diameter
+ diam_part = 2 * def_comp(m)%radiusmym
+ rho_part = def_comp(m)%densitygcm3
+
+ do ip=1,numpresvg
- p= pbasevg + ip*pincrvg
- if(p < 1.) p=1.
+ p= pbasevg + ip*pincrvg
+ if(p < 1.) p=1.
- do it=1,numtempvg
+ do it=1,numtempvg
- t= tbasevg + it*tincrvg
+ t= tbasevg + it*tincrvg
- vg=vgrav(dp,rp,p,t)
- call iter(vgmod,vg,dp,rp,p,t)
+ vg=vgrav(diam_part,rho_part,p,t)
+ call iter(vgmod,vg,diam_part,rho_part,p,t)
- !..table in unit m/s (cm/s computed)
- vgtable(it,ip,n)= vgmod*0.01
+ !..table in unit m/s (cm/s computed)
+ vgtable(it,ip,n)= vgmod*0.01
+ end do
end do
- end do
+ end select
end do do_comp
end subroutine vgravtables
diff --git a/src/common/wetdep.f90 b/src/common/wetdep.f90
index 62f5423e..a425a348 100644
--- a/src/common/wetdep.f90
+++ b/src/common/wetdep.f90
@@ -15,14 +15,151 @@
! You should have received a copy of the GNU General Public License
! along with this program. If not, see .
-module wetdep
+module wetdepml
+ use iso_fortran_env, only: real64
implicit none
private
- public wetdep2, wetdep2_init
+ real, save :: vminprec = -1.
+ real, parameter :: precmin = 0.01
+
+ public :: wetdep, init, deinit, requires_extra_precip_fields, &
+ wetdep_precompute
+ public :: operator(==), operator(/=)
+
+ type, public :: wetdep_subcloud_scheme_t
+ integer, private :: scheme
+ character(len=32), public :: description
+ end type
+
+ type(wetdep_subcloud_scheme_t), parameter, public :: WETDEP_SUBCLOUD_SCHEME_UNDEFINED = &
+ wetdep_subcloud_scheme_t(0, "Not defined")
+ type(wetdep_subcloud_scheme_t), parameter, public :: WETDEP_SUBCLOUD_SCHEME_NONE = &
+ wetdep_subcloud_scheme_t(1, "No scheme (skip)")
+ type(wetdep_subcloud_scheme_t), parameter, public :: WETDEP_SUBCLOUD_SCHEME_BARTNICKI = &
+ wetdep_subcloud_scheme_t(2, "Bartnicki")
+ type(wetdep_subcloud_scheme_t), parameter, public :: WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL = &
+ wetdep_subcloud_scheme_t(3, "Conventional")
+
+ type, public :: wetdep_incloud_scheme_t
+ integer, private :: scheme
+ character(len=32), public :: description
+ end type
+
+ type(wetdep_incloud_scheme_t), parameter, public :: WETDEP_INCLOUD_SCHEME_UNDEFINED = &
+ wetdep_incloud_scheme_t(0, "Not defined")
+ type(wetdep_incloud_scheme_t), parameter, public :: WETDEP_INCLOUD_SCHEME_NONE = &
+ wetdep_incloud_scheme_t(1, "No scheme (skip)")
+ type(wetdep_incloud_scheme_t), parameter, public :: WETDEP_INCLOUD_SCHEME_TAKEMURA = &
+ wetdep_incloud_scheme_t(2, "Takemura")
+ type(wetdep_incloud_scheme_t), parameter, public :: WETDEP_INCLOUD_SCHEME_ROSELLE = &
+ wetdep_incloud_scheme_t(3, "Roselle")
+
+ type, public :: wetdep_scheme_t
+ type(wetdep_subcloud_scheme_t) :: subcloud
+ type(wetdep_incloud_scheme_t) :: incloud
+ !> Use 3D precip and precomputed parameters
+ logical :: use_vertical
+ !> Whether to use cloud fraction correction for belowcloud schemes
+ logical :: use_cloudfraction
+ end type
+
+ interface operator (==)
+ module procedure :: equal_subcloud_scheme, equal_incloud_scheme
+ end interface
+
+ interface operator (/=)
+ module procedure :: not_equal_subcloud_scheme, not_equal_incloud_scheme
+ end interface
+
+
+ type, public :: conventional_params_t
+ real(real64) :: A
+ real(real64) :: B
+ end type
+
+ type(conventional_params_t), save, public :: conventional_params = conventional_params_t(0.0, 0.0)
+ type(conventional_params_t), parameter, public :: RATM_params = conventional_params_t(2.98e-5, 0.75)
+ real(real64), allocatable, save, public :: conventional_deprate_m1(:,:)
+
+ type(wetdep_scheme_t), save, public :: wetdep_scheme = &
+ wetdep_scheme_t(WETDEP_SUBCLOUD_SCHEME_UNDEFINED, WETDEP_INCLOUD_SCHEME_UNDEFINED, .false., .false.)
contains
+ logical pure function equal_subcloud_scheme(this, other) result(eq)
+ type(wetdep_subcloud_scheme_t), intent(in) :: this, other
+ eq = this%scheme == other%scheme
+ end function
+
+ logical pure function not_equal_subcloud_scheme(this, other) result(eq)
+ type(wetdep_subcloud_scheme_t), intent(in) :: this, other
+ eq = .not. (this == other)
+ end function
+
+ logical pure function equal_incloud_scheme(this, other) result(eq)
+ type(wetdep_incloud_scheme_t), intent(in) :: this, other
+ eq = this%scheme == other%scheme
+ end function
+
+ logical pure function not_equal_incloud_scheme(this, other) result(eq)
+ type(wetdep_incloud_scheme_t), intent(in) :: this, other
+ eq = .not. (this == other)
+ end function
+
+ subroutine init(tstep)
+ use snapdimML, only: nx, ny
+ real, intent(in) :: tstep
+ if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL .and. &
+ wetdep_scheme%use_vertical) then
+ allocate(conventional_deprate_m1(nx,ny))
+ endif
+
+ if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_BARTNICKI.and. &
+ .not.wetdep_scheme%use_vertical) then
+ call wetdep2_init(tstep)
+ endif
+
+ end subroutine
+
+ subroutine deinit()
+ if (allocated(conventional_deprate_m1)) deallocate(conventional_deprate_m1)
+ end subroutine
+
+ pure logical function requires_extra_precip_fields()
+ requires_extra_precip_fields = wetdep_scheme%use_vertical
+ end function
+
+ !> Move activity from particle to wet deposition field
+ !> depending on the precipitation at the place of the
+ !> particle.
+ subroutine wetdep(tstep, part, pextra)
+ use iso_fortran_env, only: real64
+ USE particleML, only: particle, extraParticle
+ use snapparML, only: def_comp
+ use snapfldml, only: depwet
+
+ real, intent(in) :: tstep
+ type(particle), intent(inout) :: part
+ type(extraParticle), intent(in) :: pextra
+
+ if (def_comp(part%icomp)%kwetdep == 1) then
+ if (wetdep_scheme%use_vertical) then
+ block
+ use snapfldML, only: wscav, depwet
+ call wetdep_using_precomputed_wscav(part, wscav, depwet, tstep)
+ end block
+ else
+ if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_BARTNICKI) then
+ call wetdep2(depwet, tstep, part, pextra)
+ elseif (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL) then
+ call wetdep_conventional(depwet, part, tstep)
+ endif
+ endif
+ endif
+ end subroutine
+
+
!> Purpose: Compute wet deposition for each particle and each component
!> and store depositions in nearest gridpoint in a field
!>
@@ -48,8 +185,11 @@ subroutine wetdep2(depwet, tstep, part, pextra)
real :: precint, deprate, dep, q
m = part%icomp
- if (def_comp(m)%kwetdep == 1 .AND. pextra%prc > 0.0 &
- .AND. part%z > 0.67) then
+
+ !..reset precipitation to zero if pressure less than approx. 550 hPa.
+ !..and if less than a minimum precipitation intensity (mm/hour)
+ if (def_comp(m)%kwetdep == 1 .AND. pextra%prc > precmin &
+ .AND. part%z > 0.67 .AND. part%z > vminprec) then
!..find particles with wet deposition and
!..reset precipitation to zero if not wet deposition
precint = pextra%prc
@@ -71,10 +211,11 @@ subroutine wetdep2(depwet, tstep, part, pextra)
end subroutine wetdep2
!> Initialisation routine for ::wetdep2
+!> Setting vminprec
+!> Output diagnostics
subroutine wetdep2_init(tstep)
! initalization
USE snapdebug, only: iulog
- USE particleML, only: Particle, extraParticle
USE snapparML, only: ncomp
USE snapparML, only: run_comp, def_comp
@@ -109,6 +250,33 @@ subroutine wetdep2_init(tstep)
1010 format(1x, f5.1, ':', 12f7.4)
end do
write (iulog, *) '-------------------------------------------------'
+
+ block
+ USE snapgrdML, only: alevel, blevel, vlevel
+ USE snapdimML, only: nk
+ USE snapdebug, only: iulog
+
+ integer :: k
+ real :: p1,p2
+ real, parameter :: plim = 550.0
+
+ if(vminprec < 0.) then
+ p2 = 1000.
+ k = 1
+ do while (p2 > plim .AND. k < nk)
+ k = k+1
+ p1 = p2
+ p2 = alevel(k) + blevel(k)*1000.
+ end do
+ if(k > 1) then
+ vminprec = vlevel(k-1) &
+ + (vlevel(k)-vlevel(k-1))*(p1-plim)/(p1-p2)
+ else
+ vminprec = vlevel(nk)
+ end if
+ write(iulog,*) 'POSINT. precmin,vminprec: ',precmin,vminprec
+ end if
+ end block
end subroutine
pure real function wet_deposition_constant(rm) result(depconst)
@@ -123,21 +291,25 @@ pure real function wet_deposition_constant(rm) result(depconst)
depconst = b0 + b1*rm + b2*rm*rm + b3*rm*rm*rm
end function
- pure real function wet_deposition_rate(radius, q, depconst, tstep) result(deprate)
+ pure elemental real function wet_deposition_rate_bartnicki_imm(radius, q, depconst, no_use_convective_rain) result(rkw)
!> radius in micrometer
real, intent(in) :: radius
!> precipitation intensity in mm/h
real, intent(in) :: q
!> deposition constant
real, intent(in) :: depconst
- !> length of a timestep in hours
- real, intent(in) :: tstep
+ logical, optional, intent(in) :: no_use_convective_rain
real, parameter :: a0 = 8.4e-5
real, parameter :: a1 = 2.7e-4
real, parameter :: a2 = -3.618e-6
+ logical :: use_convective
- real :: rkw
+ if (present(no_use_convective_rain)) then
+ use_convective = .not.no_use_convective_rain
+ else
+ use_convective = .true.
+ endif
rkw = 0
if (radius > 0.05 .AND. radius <= 1.4) then
@@ -149,12 +321,344 @@ pure real function wet_deposition_rate(radius, q, depconst, tstep) result(deprat
if (radius > 10.0) then
rkw = a1*q + a2*q*q
endif
- if (q > 7.0) then ! convective
+ if (use_convective .and. q > 7.0) then ! convective
rkw = 3.36e-4*q**0.79
endif
if (radius <= 0.1) then ! gas
rkw = 1.12e-4*q**0.79
endif
+ end function
+
+ pure elemental real function wet_deposition_rate(radius, q, depconst, tstep) result(deprate)
+ !> radius in micrometer
+ real, intent(in) :: radius
+ !> precipitation intensity in mm/h
+ real, intent(in) :: q
+ !> deposition constant
+ real, intent(in) :: depconst
+ !> length of a timestep in seconds
+ real, intent(in) :: tstep
+
+ real :: rkw
+
+ rkw = wet_deposition_rate_bartnicki_imm(radius, q, depconst)
+
deprate = 1.0 - exp(-tstep*rkw)
end function
-end module wetdep
+
+ subroutine wetdep_conventional_compute(precip)
+ real, intent(in) :: precip(:,:)
+
+ conventional_deprate_m1(:,:) = conventional_params%A * (precip ** conventional_params%B)
+ end subroutine
+
+ subroutine wetdep_conventional(depwet, part, tstep)
+ USE iso_fortran_env, only: real64
+ USE particleML, only: Particle, extraParticle
+ USE snapparML, only: def_comp
+ USE snapdimML, only: hres_pos
+
+!> Field which ret deposition gets added to
+ real(real64), intent(inout) :: depwet(:, :, :)
+!> particle
+ type(Particle), intent(inout) :: part
+ real, intent(in) :: tstep
+
+ integer :: m, i, j, mm
+ real :: dep, deprate
+
+ m = part%icomp
+
+ i = nint(part%x)
+ j = nint(part%y)
+
+ deprate = exp(-tstep * conventional_deprate_m1(i, j))
+
+ i = hres_pos(part%x)
+ j = hres_pos(part%y)
+ dep = part%scale_rad(deprate)
+
+ mm = def_comp(m)%to_running
+
+ !$OMP atomic
+ depwet(i, j, mm) = depwet(i, j, mm) + dble(dep)
+ end subroutine
+
+ !> Aerosol rainout process also known as GCM-type wet deposition process
+ subroutine wetdep_incloud_takemura(lambda, q, cloud_water, cloud_fraction)
+ real, intent(out) :: lambda(:,:)
+ !> vertical flux of hydrometeors
+ real, intent(in) :: q(:,:)
+ !> Fraction of aerosol mass in cloud water to total aerosol mass in the grid
+ !> or the absorbtion coefficient
+ !> Usually very high
+ real, parameter :: f_inc = 1.0
+ !> cloud water
+ real, intent(in) :: cloud_water(:,:)
+ !> Cloud fraction
+ real, intent(in) :: cloud_fraction(:,:)
+
+
+ where (q > 0.0)
+ lambda = q/(q + cloud_water) * f_inc * cloud_fraction / 3600.0
+ elsewhere
+ lambda = 0.0
+ end where
+ end subroutine
+
+ !> Aerosol rainout process from Roselle and Binkowski 1999
+ subroutine wetdep_incloud_roselle(lambda, q, cloud_water, cloud_fraction)
+ real, intent(out) :: lambda(:,:)
+ !> vertical flux of hydrometeors
+ real, intent(in) :: q(:,:)
+ !> Fraction of aerosol mass in cloud water to total aerosol mass in the grid
+ !> or the absorbtion coefficient
+ !> Usually very high
+ real, parameter :: f_inc = 1.0
+ !> cloud water
+ real, intent(in) :: cloud_water(:,:)
+ !> Cloud fraction
+ real, intent(in) :: cloud_fraction(:,:)
+
+ !> Characteristic time of cloud (equal to timestep?)
+ real, parameter :: tau_cloud = 1.0
+
+ real, allocatable :: tau_washout(:,:)
+
+ allocate(tau_washout, mold=lambda)
+
+ where (q > 0.0)
+ tau_washout = cloud_water / q
+ lambda = (1.0 - exp(-tau_cloud / tau_washout)) / tau_cloud * f_inc * cloud_fraction / 3600.0
+ elsewhere
+ lambda = 0.0
+ endwhere
+ end subroutine
+
+ subroutine wet_deposition_rate_bartnicki(wscav, radius, precip, ccf, use_ccf)
+ real, intent(out) :: wscav(:,:)
+ real, intent(in) :: precip(:,:)
+ real, intent(in) :: radius
+ real, intent(in) :: ccf(:,:)
+ logical, intent(in) :: use_ccf
+
+ real :: depconst
+ integer :: nx, ny
+
+ nx = size(precip,1)
+ ny = size(precip,2)
+
+ depconst = wet_deposition_constant(radius)
+ if (.not.use_ccf) then
+ wscav(:,:) = wet_deposition_rate_bartnicki_imm(radius, precip, depconst, no_use_convective_rain=.true.)
+ else
+ block
+ integer :: i, j
+ real :: precip_scaled
+
+ do j=1,ny
+ do i=1,nx
+ if (precip(i,j) <= 0.0) then
+ wscav(i,j) = 0.0
+ cycle
+ endif
+
+ if (ccf(i,j) > 0.0) then
+ ! Scale up precip intensity
+ precip_scaled = precip(i,j) / ccf(i,j)
+ else
+ precip_scaled = precip(i,j)
+ endif
+
+ wscav(i,j) = wet_deposition_rate_bartnicki_imm(radius, precip_scaled, depconst, no_use_convective_rain=.true.)
+
+ if (ccf(i,j) > 0.0) then
+ ! Scale down efficiency
+ wscav(i,j) = wscav(i,j) * ccf(i,j)
+ endif
+ enddo
+ enddo
+ end block
+ endif
+ end subroutine
+
+
+ subroutine wet_deposition_rate_ratm(wscav, precip, ccf, use_ccf)
+ real, intent(out) :: wscav(:,:)
+ real, intent(in) :: precip(:,:)
+ real, intent(in) :: ccf(:,:)
+ logical, intent(in) :: use_ccf
+
+ integer :: i, j
+ real adj_precip
+ integer :: nx, ny
+
+ nx = size(precip,1)
+ ny = size(precip,2)
+
+ do j=1,ny
+ do i=1,nx
+ if (use_ccf .and. (ccf(i,j) > 0.0)) then
+ adj_precip = precip(i,j) / ccf(i,j)
+ wscav(i,j) = ccf(i,j) * conventional_params%A * adj_precip ** conventional_params%B
+ else
+ wscav(i,j) = conventional_params%A * (precip(i,j) ** conventional_params%B)
+ endif
+ enddo
+ enddo
+ end subroutine
+
+ !> Should be called every input timestep to
+ !> prepare the scavenging rates
+ subroutine wetdep_precompute()
+ use snapparML, only: ncomp, run_comp
+ use snapfldML, only: wscav, cw3d, precip3d, cloud_cover, precip
+
+ integer :: i
+
+ do i=1,ncomp
+ if (.not.run_comp(i)%defined%kdrydep == 1) cycle
+ if (wetdep_scheme%use_vertical) then
+ if (.not.(allocated(precip3d).and.allocated(cw3d).and.allocated(wscav))) then
+ error stop "Some wetdep/precip fields not allocated"
+ endif
+ call prepare_wetdep_3d(wscav(:,:,:,i), run_comp(i)%defined%radiusmym, precip3d, cw3d, cloud_cover)
+ else if (wetdep_scheme%subcloud == WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL) then
+ call wetdep_conventional_compute(precip)
+ endif
+ enddo
+ end subroutine
+
+
+ !> Precompute wet scavenging coefficients
+ subroutine prepare_wetdep_3d(wscav, radius, precip, cw, ccf)
+ !> Wet scavenging coefficient [1/s]
+ real, intent(out) :: wscav(:,:,:)
+ !> Precipitation intensity
+ real, intent(in) :: precip(:,:,:)
+ !> Radius of particle
+ real, intent(in) :: radius
+ !> Cloud water
+ real, intent(in) :: cw(:,:,:)
+ !> Cloud cover fraction
+ real, intent(in) :: ccf(:,:,:)
+
+ real, allocatable :: wscav_tmp(:,:,:)
+ real, allocatable :: accum_precip(:,:), accum_ccf(:,:)
+
+ integer :: nk, k, nx, ny
+
+ nx = size(wscav,1)
+ ny = size(wscav,2)
+ nk = size(wscav,3)
+
+ allocate(wscav_tmp, mold=wscav)
+ allocate(accum_precip(nx,ny), accum_ccf(nx,ny))
+
+ accum_precip(:,:) = 0.0
+ accum_ccf(:,:) = 0.0
+ do k=nk,1,-1
+ ! Accumulated precipitation in the column
+ accum_precip(:,:) = accum_precip(:,:) + precip(:,:,k)
+ accum_ccf(:,:) = accum_ccf(:,:) + ccf(:,:,k)
+ where (accum_ccf >= 1.0)
+ accum_ccf = 1.0
+ endwhere
+
+ ! Subcloud
+ select case (wetdep_scheme%subcloud%scheme)
+ case (WETDEP_SUBCLOUD_SCHEME_BARTNICKI%scheme)
+ call wet_deposition_rate_bartnicki(wscav(:,:,k), radius, accum_precip(:,:), &
+ accum_ccf(:,:), use_ccf=wetdep_scheme%use_cloudfraction)
+ case (WETDEP_SUBCLOUD_SCHEME_CONVENTIONAL%scheme)
+ call wet_deposition_rate_ratm(wscav(:,:,k), accum_precip(:,:), &
+ accum_ccf(:,:), use_ccf=wetdep_scheme%use_cloudfraction)
+ case (WETDEP_SUBCLOUD_SCHEME_NONE%scheme)
+ wscav(:,:,k) = 0.0
+ case default
+ error stop wetdep_scheme%subcloud%description
+ end select
+
+ ! Incloud
+ select case (wetdep_scheme%incloud%scheme)
+ case (WETDEP_INCLOUD_SCHEME_NONE%scheme)
+ wscav_tmp(:,:,k) = 0.0
+ case (WETDEP_INCLOUD_SCHEME_ROSELLE%scheme)
+ call wetdep_incloud_roselle(wscav_tmp(:,:,k), precip(:,:,k), cw(:,:,k), ccf(:,:,k))
+ case (WETDEP_INCLOUD_SCHEME_TAKEMURA%scheme)
+ call wetdep_incloud_takemura(wscav_tmp(:,:,k), precip(:,:,k), cw(:,:,k), ccf(:,:,k))
+ case default
+ error stop wetdep_scheme%incloud%description
+ end select
+ wscav(:,:,k) = max(wscav(:,:,k), wscav_tmp(:,:,k))
+ end do
+ end subroutine
+
+ subroutine wetdep_using_precomputed_wscav(part, wscav, dep, tstep)
+ use iso_fortran_env, only: real64
+ use particleml, only: particle
+ use snapparML, only: def_comp
+ use snapgrdML, only: ivlevel
+ use snapdimML, only: hres_pos
+ type(particle), intent(inout) :: part
+ real, intent(in) :: wscav(:,:,:,:)
+ real(real64), intent(inout) :: dep(:,:,:)
+ real, intent(in) :: tstep
+
+ real :: radlost, rkw
+ integer :: ivlvl, i, j, k, mm
+
+ ivlvl = part%z * 10000.0
+ k = ivlevel(ivlvl)
+ i = nint(part%x)
+ j = nint(part%y)
+
+ mm = def_comp(part%icomp)%to_running
+
+ rkw = wscav(i,j,k,mm)
+
+ radlost = part%scale_rad(exp(-rkw*tstep))
+
+ i = hres_pos(part%x)
+ j = hres_pos(part%y)
+ !$OMP atomic
+ dep(i, j, mm) = dep(i, j, mm) + real(radlost, kind=real64)
+ end subroutine
+
+
+ !> Laakso et al. 2003, Ultrafine particle scavenging coefficients
+ !> Valid for particle diameter in between 10 nm < diameter < 510 nm
+ !> Valid for precipitation intensity between 0-20 mm/h
+ elemental subroutine laakso_preciprate(dia_p, p, lambda)
+ !> Diameter of particle [m]
+ real, intent(in) :: dia_p
+ !> Precipitation intensity [mm/hr]
+ real, intent(in) :: p
+ !> Scavenging rate [1/s]
+ real, intent(out) :: lambda
+
+ ! [1/s]
+ real, parameter :: lambda0 = 1.0
+ ! [m]
+ real, parameter :: dia_p0 = 1.0
+ ! [mm/hr]
+ real, parameter :: p0 = 1.0
+
+ real, parameter :: a = 274.35758
+ real, parameter :: b = 332839.59273
+ real, parameter :: c = 226656.57259
+ real, parameter :: d = 58005.91340
+ real, parameter :: e = 6588.38582
+ real, parameter :: f = 0.244984
+
+ real(real64) :: dp, log10_l_l0
+
+ dp = log10(dia_p/ dia_p0)
+
+ log10_l_l0 = a + b*dp**-4 + c*dp**-3 + d * dp**-2 + e * dp**-1 + f*(p/p0)**0.5
+
+ lambda = lambda0 * 10.0**log10_l_l0
+
+ end subroutine
+
+end module wetdepml
diff --git a/src/gcc_pkgconfig.mk b/src/gcc_pkgconfig.mk
index a04e7e9c..1a61d790 100644
--- a/src/gcc_pkgconfig.mk
+++ b/src/gcc_pkgconfig.mk
@@ -5,13 +5,16 @@ F77 = gfortran
F77FPEFLAGS= -ffpe-trap=invalid,zero,overflow
F77BOUNDFLAGS= -fbounds-check
-F77FLAGS=-DVERSION=\"$(VERSION)\" -O2 -ftree-vectorize -fno-math-errno -g -mavx2 -mfma -Wall -Wextra -fimplicit-none -fmodule-private -Wno-conversion
+F77FLAGS=-DVERSION=\"$(VERSION)\" -O2 -ftree-vectorize -fno-math-errno -g -mavx2 -mfma -Wall -Wextra -fimplicit-none -fmodule-private -Wno-conversion -Wno-compare-reals
ifdef SNAP_DEBUG_CHECKS
F77FLAGS+=$(F77BOUNDFLAGS) $(F77FPEFLAGS)
endif
ifdef SNAP_BOUND_CHECKS
F77FLAGS+=$(F77BOUNDFLAGS)
endif
+ifdef SNAP_USE_OMP
+ F77FLAGS+=-fopenmp
+endif
# optional versioned fimex
FIMEX = fimex
diff --git a/src/test/.gitignore b/src/test/.gitignore
index 53e3a496..3198dd14 100644
--- a/src/test/.gitignore
+++ b/src/test/.gitignore
@@ -2,3 +2,5 @@ bsnap_naccident
data/snap.log
data/snap.nc
data/snap_ecemep.nc
+data/snap_meps_fimex.nc
+data/snap_meps_fimex.log
diff --git a/src/test/Makefile b/src/test/Makefile
index dab383ef..0198979d 100644
--- a/src/test/Makefile
+++ b/src/test/Makefile
@@ -13,8 +13,8 @@ include ../common/snap.mk
TESTFILES = snap_testdata/meteo20200415_00_ringhals.nc \
snap_testdata/meteo20200415_01_ringhals.nc \
snap_testdata/meps_det_2_5km_20200316T00Z_ringhals.nc \
- snap_testdata/snap_ecemep_expected4.nc \
- snap_testdata/snap_meps_interpolated_expected3.nc
+ snap_testdata/snap_ecemep_expected5.nc \
+ snap_testdata/snap_meps_interpolated_expected4.nc
install:
diff --git a/src/test/data/snap.input_meps_fimex b/src/test/data/snap.input_meps_fimex
index 34459184..7ad37f31 100644
--- a/src/test/data/snap.input_meps_fimex
+++ b/src/test/data/snap.input_meps_fimex
@@ -110,8 +110,8 @@ FIELD_TIME.FORECAST
** OUTPUT FILES
*
FIELD.OUTTYPE=netcdf
-FIELD.OUTPUT= snap.nc
-LOG.FILE= snap.log
+FIELD.OUTPUT= snap_meps_fimex.nc
+LOG.FILE= snap_meps_fimex.log
*
* DEBUG.OFF
DEBUG.ON
diff --git a/src/test/snap_forward_test.py b/src/test/snap_forward_test.py
index cfea8ed0..0220de4a 100644
--- a/src/test/snap_forward_test.py
+++ b/src/test/snap_forward_test.py
@@ -16,7 +16,7 @@ class SnapEcEMEPForwardTestCase(SnapTestCase):
snap = "../bsnap_naccident"
datadir = os.path.join(os.path.dirname(os.path.realpath(__file__)), "data")
- snapExpected = "snap_testdata/snap_ecemep_expected4.nc"
+ snapExpected = "snap_testdata/snap_ecemep_expected5.nc"
def setUp(self):
pass
@@ -87,7 +87,7 @@ class SnapMEPSForwardTestCase(SnapTestCase):
snap = "../bsnap_naccident"
datadir = os.path.join(os.path.dirname(os.path.realpath(__file__)), "data")
- snapExpected = "snap_testdata/snap_meps_interpolated_expected3.nc"
+ snapExpected = "snap_testdata/snap_meps_interpolated_expected4.nc"
def setUp(self):
pass