From 3906a24b4ee8a222eac18edc1c204d2f0e604f3c Mon Sep 17 00:00:00 2001
From: lucieb92 <lucie.bakels@univie.ac.at>
Date: Tue, 27 Feb 2024 16:38:28 +0100
Subject: [PATCH] Bugfix in wetdeposition by Massimo Cassiani
Former-commit-id: 24c1c5773d63d9641e395a8699765073ffb74646
---
src/par_mod.f90 | 7 ++++--
src/wetdepo_mod.f90 | 37 ++++++++++++++++++++++++++-----
tests/pathnames_nests | 3 ++-
tests/run_default_options_test.sh | 3 ++-
tests/run_nests_test.sh | 1 +
5 files changed, 42 insertions(+), 9 deletions(-)
diff --git a/src/par_mod.f90 b/src/par_mod.f90
index 09cedf7b..17d2e75e 100644
--- a/src/par_mod.f90
+++ b/src/par_mod.f90
@@ -80,7 +80,7 @@ module par_mod
!real,parameter :: d_trop=50., d_strat=0.1
real :: d_trop=50., d_strat=0.1, fturbmeso=0.16 ! turbulence factors can change for different runs
real,parameter :: rho_water=1000. !ZHG 2015 [kg/m3]
- real,parameter :: ratio_incloud=6.2 !ZHG MAR2016
+ real,parameter :: ratio_incloud=0.1 !MC 2024
real,parameter :: wet_a=1.e-5, wet_b=0.8 !AT
! karman Karman's constant
@@ -92,7 +92,7 @@ module par_mod
! yield the scales for the mesoscale wind velocity fluctuations
! d_trop [m2/s] Turbulent diffusivity for horiz components in the troposphere
! d_strat [m2/s] Turbulent diffusivity for vertical component in the stratosphere
- ! ratio_incloud ZHG MAR2016
+ ! ratio_incloud MC 2024, dimensionless ratio that should be <= 1
! wet_a, wet_b for wetscav=wet_a*prec**wet_b if no cloud found, but precipitation occurs
@@ -306,4 +306,7 @@ module par_mod
logical,parameter :: lpartoutputperfield=.false.
+ ! Temporary parameter to switch off the gridfaction calculation in the wetdeposition
+ logical,parameter :: lgridfraction=.false.
+
end module par_mod
diff --git a/src/wetdepo_mod.f90 b/src/wetdepo_mod.f90
index 7f696a92..5765d74d 100644
--- a/src/wetdepo_mod.f90
+++ b/src/wetdepo_mod.f90
@@ -325,6 +325,8 @@ subroutine get_wetscav(itime,jpart,ks,gridfract,wetscav)
! Interpolate cloud information
call interpol_rain(itime,kz,lsp,convp,cc,t_particle,cl,icbot,ictop,icmv)
+ ! cc = total cloud cover
+ ! cl = ctwc
! If total precipitation is less than precsub=0.01 mm/h - no scavenging
! Note: original PS version (in order avoid step at 0.01)
@@ -678,7 +680,8 @@ subroutine get_wetscav_incld_aerosol(ks,gridfract,prec,cl,cc,t_particle,wetscav)
! scavenging coefficient based on Hertel et al 1995 -
! using si (S_i in paper) for both gas and aerosol
- !SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
+ ! wetscav = ratio_incloud*si*prec/3.6E6/cloud_thickness
+ ! cloud_thickness cancels out since cl is computed without
wetscav=ratio_incloud*si* prec/3.6E6 ! mm/h -> m/s
end subroutine get_wetscav_incld_aerosol
@@ -710,7 +713,8 @@ subroutine get_wetscav_incld_gas(ks,gridfract,prec,cl,cc,t_particle,wetscav)
! scavenging coefficient based on Hertel et al 1995 -
! using si (S_i in paper) for both gas and aerosol
- !SEC wetscav fix, the cloud height is no longer needed, it gives wrong results
+ ! wetscav = ratio_incloud*si*prec/3.6E6/cloud_thickness
+ ! cloud_thickness cancels out since cl is computed without
wetscav=ratio_incloud*si* prec/3.6E6 ! mm/h -> m/s
end subroutine get_wetscav_incld_gas
@@ -751,7 +755,12 @@ subroutine get_gridfract(lsp,convp,cc,gridfract)
real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
integer :: i, j
-
+
+ if (.not. lgridfraction) then
+ gridfract=1.0
+ return
+ endif
+
if (lsp.gt.20.) then
i=5
else if (lsp.gt.8.) then
@@ -789,10 +798,23 @@ subroutine get_cloud_liquid(gridfract,prec,cc,cl)
real, intent(in) :: prec,cc,gridfract ! precipitation in sub-grid cell
real, intent(inout) :: cl ! scavenging coefficient
!ZHG 2015 use cloud liquid & ice water (CLWC+CIWC) from ECMWF
+
+ ! MC -- Integrated water content:
+ ! CTWC = SUM(CLWC * rho_water * cloud_thickness) [kg/kg * kg/m3 * m]
+ ! -> Average water content: cl = CTWC/rho_water/cloud_thickness [m3(water)/m3(cloud)]
+
+ ! Note that cloud_thickness is not included, since this will cancel out when
+ ! computing the wetscavenging: Wetscav=ratio_incloud*S_i*(prec/3.6E6)/cloud_thickness
+ ! S_i=1/cl
! Mother grid
if (ngrid.le.0) then
if (lcw) then
- cl = cl*(gridfract/cc)
+ if (lgridfraction) then
+ cl=cl/rho_water*(gridfract/cc)
+ else
+ cl=cl/rho_water ! Grythe et al. eq (1), no cloud_thickness since this will cancel out later
+ endif
+ ! cl = cl*(gridfract/cc)
! A.Plach 2021 cl should not become too small
! cl=max(0.2*prec**0.36, cl*(gridfract/cc))
else ! no cloud water available, use parameterisation for cloud water [m2/m3]
@@ -806,7 +828,12 @@ subroutine get_cloud_liquid(gridfract,prec,cc,cl)
endif
else
if (lcw_nest(ngrid)) then
- cl = cl*(gridfract/cc)
+ if (lgridfraction) then
+ cl=cl/rho_water*(gridfract/cc)
+ else
+ cl=cl/rho_water ! Grythe et al. eq (1), no cloud_thickness since this will cancel out later
+ endif
+ !cl = cl*(gridfract/cc)
else
cl=0.2*prec**0.36
endif
diff --git a/tests/pathnames_nests b/tests/pathnames_nests
index 864327fa..7ee7831a 100644
--- a/tests/pathnames_nests
+++ b/tests/pathnames_nests
@@ -2,4 +2,5 @@
./output/
./default_winds/
./default_winds/AVAILABLE_glob
-./default_winds/
\ No newline at end of file
+./default_winds/
+============================================
\ No newline at end of file
diff --git a/tests/run_default_options_test.sh b/tests/run_default_options_test.sh
index 6abd64b0..c2919614 100644
--- a/tests/run_default_options_test.sh
+++ b/tests/run_default_options_test.sh
@@ -283,7 +283,7 @@ rm -rf ./current ./output/*
# sed -i "/LDIRECT=/c\ LDIRECT= -1," ./current/COMMAND
# sed -i "/IOUTPUTFOREACHRELEASE=/c\ IOUTPUTFOREACHRELEASE= 1," ./current/COMMAND
# sed -i "/IOUT=/c\ IOUT= 1," ./current/COMMAND
-# sed -i "/IETIME=/c\ LOUTSTEP= 010000," ./current/COMMAND
+# sed -i "/IBTIME=/c\ IBTIME= 020000," ./current/COMMAND
# sed -i "/LOUTSTEP=/c\ LOUTSTEP= 3600," ./current/COMMAND
# sed -i "/LOUTAVER=/c\ LOUTAVER= 3600," ./current/COMMAND
# cp -rf part_ic.nc output/
@@ -292,6 +292,7 @@ rm -rf ./current ./output/*
# STATUS=$((STATUS + $?))
# TESTSRUN=$((TESTSRUN + 1))
# # clean up
+# rm -rf ./current ./output/*
#
#IFLUX
cp -rf ./default_options ./current
diff --git a/tests/run_nests_test.sh b/tests/run_nests_test.sh
index d3512dbc..11be1af3 100644
--- a/tests/run_nests_test.sh
+++ b/tests/run_nests_test.sh
@@ -67,3 +67,4 @@ echo "[$MM] Tests failed: $STATUS / $TESTSRUN"
#
# Return collective error status
#
+exit $STATUS
\ No newline at end of file
--
GitLab