Skip to content

Commit 16e92dd

Browse files
mathomp4mathomp4
authored
Merge pull request #365 from zhanglikate/bugfix-ufs-wetremoval
2 parents fcedfca + e723d19 commit 16e92dd

File tree

6 files changed

+144
-26
lines changed

6 files changed

+144
-26
lines changed

‎CHANGELOG.md‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
3434
- fwet removed from children GridCompMod and placed in respective instance RC files
3535
- Updated settling routine and calls to allow settling velocity diagnostics in output field
3636
- Updated pressure lids in instance.rc files for use in GCMv12
37+
3738
### Fixed
3839

3940
- In SU an incorrect (older) optics table was specified in AMIP/AMIP.20C sub-directories.
@@ -44,6 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
4445
- Removed erroneous/extraneous friendly attributes to internal state for DU and NI
4546
when in data_driven mode
4647
- typo in filepath for BR optics file
48+
- Bug fix for "WetRemovalUFS"
4749

4850
### Added
4951

@@ -52,6 +54,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
5254
coupling of GMI OH, NO3, H2O2 to sulfur chemistry mechanism
5355
- Added a callback to allow a chemistry module to call for optical properties needed
5456
for photolysis calculation; currently used for GMI with CloudJ
57+
- Implementation of the "WetRemovalUFS" option for sulfate
5558

5659
## [v2.4.3] - 2025-07-21
5760

@@ -66,7 +69,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
6669
### Added
6770

6871
- Added option for choosing SettlingSolver
69-
- `gocart' - Default GOCART Settling scheme | 'ufs' - New Settling scheme
72+
- 'gocart' - Default GOCART Settling scheme | 'ufs' - New Settling scheme
7073

7174
## [v2.4.2] - 2025-06-12
7275

‎ESMF/GOCART2G_GridComp/SU2G_GridComp/AMIP.20C/SU2G_instance_SU.rc‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#
2-
# Resource file for Sulfer parameters.
2+
# Resource file for Sulfur parameters.
33
#
44

55
aerosol_radBands_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opticsBands_SU.v1_6.RRTMG.nc
@@ -23,7 +23,15 @@ aircraft_fuel_emission_factor: 1.0000
2323
# Scavenging efficiency per bin [km-1]
2424
fscav: 0.0 0.0 0.4 0.4
2525

26-
#Dummy wet removal efficiency for the GOCART scheme. This is not used for sulfate as of GOCART v2.5 and should be revisited.
26+
# Wet Removal Scheme Option | gocart, ufs
27+
wet_removal_scheme: gocart
28+
29+
# Rainout efficiency for sulfate wet_removal_scheme == ufs
30+
fwet_ice: 1.0
31+
fwet_snow: 1.0
32+
fwet_rain: 1.0
33+
34+
#Dummy wet removal efficiency for wet_removal_scheme == gocart. This is not used for sulfate as of GOCART v2.5 and should be revisited.
2735
fwet: 1.0 1.0 1.0 1.0
2836

2937
# Dry particle radius [um], used for settling
@@ -53,3 +61,7 @@ sigma: -1 -1 2.03 -1
5361
# -------------------------------------------------------------------------------------
5462
using_GMI: .false.
5563
disable_emissions: .false.
64+
65+
# SettlingSolver options
66+
# Options: 'gocart' or 'ufs'
67+
settling_scheme: gocart

‎ESMF/GOCART2G_GridComp/SU2G_GridComp/AMIP/SU2G_instance_SU.rc‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#
2-
# Resource file for Sulfer parameters.
2+
# Resource file for Sulfur parameters.
33
#
44

55
aerosol_radBands_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opticsBands_SU.v1_6.RRTMG.nc
@@ -23,7 +23,15 @@ aircraft_fuel_emission_factor: 1.0000
2323
# Scavenging efficiency per bin [km-1]
2424
fscav: 0.0 0.0 0.4 0.4
2525

26-
#Dummy wet removal efficiency for the GOCART scheme. This is not used for sulfate as of GOCART v2.5 and should be revisited.
26+
# Wet Removal Scheme Option | gocart, ufs
27+
wet_removal_scheme: gocart
28+
29+
# Rainout efficiency for sulfate wet_removal_scheme == ufs
30+
fwet_ice: 1.0
31+
fwet_snow: 1.0
32+
fwet_rain: 1.0
33+
34+
#Dummy wet removal efficiency for wet_removal_scheme == gocart. This is not used for sulfate as of GOCART v2.5 and should be revisited.
2735
fwet: 1.0 1.0 1.0 1.0
2836

2937
# Dry particle radius [um], used for settling
@@ -56,4 +64,4 @@ disable_emissions: .false.
5664

5765
# SettlingSolver options
5866
# Options: 'gocart' or 'ufs'
59-
settling_scheme: 'gocart'
67+
settling_scheme: gocart

‎ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_GridCompMod.F90‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -370,6 +370,7 @@ subroutine SetServices ( GC, RC )
370370
#include "SU2G_Import___.h"
371371
#include "SU2G_Internal___.h"
372372
if (MAPL_AM_I_ROOT()) then
373+
write (*,*) trim(Iam)//": Wet removal scheme is "//trim(self%wet_removal_scheme)
373374
write (*,*) trim(Iam)//": Settling scheme is "//trim(self%settling_scheme)
374375
end if
375376
end if
@@ -1172,6 +1173,7 @@ subroutine Run2 (GC, import, export, clock, RC)
11721173
type(ThreadWorkspace), pointer :: workspace
11731174
integer :: thread
11741175
integer :: i1, j1, i2, j2, km
1176+
real, dimension(3) :: rainout_eff
11751177
real, target, allocatable, dimension(:,:,:) :: RH20,RH80
11761178
real, pointer, dimension(:,:) :: flux_ptr
11771179
integer :: settling_opt
@@ -1311,12 +1313,33 @@ subroutine Run2 (GC, import, export, clock, RC)
13111313
__RC__)
13121314

13131315
KIN = .true.
1316+
select case (self%wet_removal_scheme)
1317+
case ('gocart')
13141318
call SU_Wet_Removal ( self%km, self%nbins, self%klid, self%cdt, kin, MAPL_GRAV, MAPL_AIRMW, &
13151319
delp, fMassSO4, fMassSO2, &
13161320
h2o2_init, ple, airdens, cn_prcp, ncn_prcp, pfl_lsan, pfi_lsan, t, &
13171321
nDMS, nSO2, nSO4, nMSA, DMS, SO2, SO4, dummyMSA, &
13181322
SUWT, SUPSO4, SUPSO4WT, PSO4, PSO4WET, __RC__ )
13191323

1324+
case ('ufs')
1325+
1326+
do n = 1, self%nbins
1327+
rainout_eff = 0.0
1328+
rainout_eff(1) = self%fwet_ice(n)
1329+
rainout_eff(2) = self%fwet_snow(n)
1330+
rainout_eff(3) = self%fwet_rain(n)
1331+
1332+
call MAPL_VarSpecGet(InternalSpec(n), SHORT_NAME=short_name, __RC__)
1333+
call MAPL_GetPointer(internal, NAME=short_name, ptr=int_ptr, __RC__)
1334+
call WetRemovalUFS (self%km, self%klid, n, self%cdt, 'sulfate', &
1335+
KIN, MAPL_GRAV, self%radius(n), rainout_eff, self%washout_tuning, &
1336+
self%wet_radius_thr, int_ptr, ple, t, airdens, pfl_lsan, pfi_lsan, SUWT, __RC__)
1337+
enddo
1338+
1339+
case default
1340+
_ASSERT_RC(.false.,'Unsupported wet removal scheme: '//trim(self%wet_removal_scheme),ESMF_RC_NOT_IMPL)
1341+
end select
1342+
13201343
! Certain variables are multiplied by 1.0e-9 to convert from nanometers to meters
13211344
call SU_Compute_Diags ( self%km, self%klid, self%rhop(nSO4), &
13221345
MAPL_GRAV, MAPL_PI, nSO4, self%diag_Mie, &

‎ESMF/GOCART2G_GridComp/SU2G_GridComp/SU2G_instance_SU.rc‎

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
#
2-
# Resource file for Sulfer parameters.
2+
# Resource file for Sulfur parameters.
33
#
44

55
aerosol_radBands_optics_file: ExtData/chemistry/AerosolOptics/v0.0.0/x/opticsBands_SU.v1_6.RRTMG.nc
@@ -20,10 +20,18 @@ so4_anthropogenic_fraction: 0.03
2020
# Aircraft emission factor: convert input unit to kg SO2
2121
aircraft_fuel_emission_factor: 1.0000
2222

23-
# Scavenging efficiency per bin [km-1]
23+
# Scavenging efficiency per bin [km-1]
2424
fscav: 0.0 0.0 0.4 0.4
2525

26-
#Dummy wet removal efficiency for the GOCART scheme. This is not used for sulfate as of GOCART v2.5 and should be revisited.
26+
# Wet Removal Scheme Option | gocart, ufs
27+
wet_removal_scheme: gocart
28+
29+
# Rainout efficiency for sulfate wet_removal_scheme == ufs
30+
fwet_ice: 1.0
31+
fwet_snow: 1.0
32+
fwet_rain: 1.0
33+
34+
#Dummy wet removal efficiency for wet_removal_scheme == gocart. This is not used for sulfate as of GOCART v2.5 and should be revisited.
2735
fwet: 1.0 1.0 1.0 1.0
2836

2937
# Dry particle radius [um], used for settling
@@ -53,4 +61,4 @@ disable_emissions: .false.
5361

5462
# SettlingSolver options
5563
# Options: 'gocart' or 'ufs'
56-
settling_scheme: 'gocart'
64+
settling_scheme: gocart

‎Process_Library/GOCART2G_Process.F90‎

Lines changed: 80 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -3783,11 +3783,34 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
37833783
do j = jl, ju
37843784
do i = il, iu
37853785
! -- compute column quantities
3786-
do k = ktop, kbot
3786+
! ---- TOP LAYER (k = 1) ----
3787+
k = ktop
3788+
km1 = k ! prevents k=0 indexing
3789+
3790+
! layer pressure thickness (use k and k+1)
3791+
delp = ple(i,j,k+1) - ple(i,j,k)
3792+
dpog(k) = delp / grav
3793+
delz = dpog(k) / rhoa(i,j,k)
3794+
delz_cm(k) = delz * m_to_cm
3795+
3796+
! no precipitation entering from above
3797+
qq(k) = 0.0
3798+
pdwn(k) = 0.0
3799+
3800+
! initialize aerosol column mass and wet-deposition accumulator
3801+
conc(k) = aerosol(i,j,k) * dpog(k)
3802+
dconc(k) = 0.0
3803+
! ---- REMAINING LAYERS (k = 2 .. km) ----
3804+
do k = ktop+1, kbot
37873805
km1 = k - 1
37883806

37893807
! -- initialize auxiliary arrays
3790-
delp = ple(i,j,k) - ple(i,j,km1)
3808+
if (k < kbot) then
3809+
delp = ple(i,j,k+1) - ple(i,j,k)
3810+
else
3811+
delp = ple(i,j,k) - ple(i,j,k-1) ! bottom-most layer
3812+
end if
3813+
37913814
dpog(k) = delp / grav
37923815
delz = dpog(k) / rhoa(i,j,k) ! thickness of layer [m]
37933816
delz_cm(k) = delz * m_to_cm ! thickness of layer [cm]
@@ -3816,7 +3839,11 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
38163839
dconc(k) = zero
38173840

38183841
! -- compute mixing ratio of saturated water vapour over ice
3819-
pres = 0.5 * ( ple(i,j,km1) + ple(i,j,k) )
3842+
if (k < kbot) then
3843+
pres = 0.5 * ( ple(i,j,k) + ple(i,j,k+1) )
3844+
else
3845+
pres = 0.5 * ( ple(i,j,k-1) + ple(i,j,k) )
3846+
end if
38203847
c_h2o(k) = 10._dp ** (-2663.5_dp / tmpu(i,j,k) + 12.537_dp ) / pres
38213848

38223849
! -- estimate cloud ice and liquid water content
@@ -3836,18 +3863,25 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
38363863
if (qq(k) > qq_thr) then
38373864
! -- compute rainout rate
38383865
k_rain = k_min + qq(k) / cwc
3839-
f = qq(k) / ( k_rain * cwc )
3866+
! f = qq(k) / ( k_rain * cwc )
3867+
3868+
! -- Safety check: ensure denominator is positive
3869+
if (k_rain * cwc > 1.0e-12) then
3870+
f = qq(k) / ( k_rain * cwc )
3871+
f = min(one, max(zero, f)) ! Bound f between 0 and 1
3872+
else
3873+
f = zero
3874+
end if
38403875

38413876
call rainout( kin, rainout_eff, f, k_rain, dt, tmpu(i,j,k), delz_cm(k), &
38423877
pdwn(k), c_h2o(k), cldice(k), cldliq(k), spc, lossfrac )
38433878

38443879
! -- compute and apply effective loss fraction
38453880
wetloss = lossfrac * conc(k)
3881+
wetloss = max(zero, min(wetloss, conc(k)))
38463882
conc(k) = conc(k) - wetloss
38473883
dconc(k) = wetloss
38483884

3849-
! -- add to total column deposition flux
3850-
dconc(k) = wetloss
38513885
end if
38523886

38533887
! -- middle layers
@@ -3860,6 +3894,7 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
38603894
if (qq(k) > qq_thr) then
38613895
k_rain = k_min + qq(k) / cwc
38623896
f_prime = qq(k) / ( k_rain * cwc )
3897+
f_prime = min(one, max(zero, f_prime))
38633898
end if
38643899

38653900
! -- account for precipitation flux
@@ -3882,6 +3917,13 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
38823917
wetloss = lossfrac * conc(k)
38833918
conc(k) = conc(k) - wetloss
38843919
dconc(k) = dconc(km1) + wetloss
3920+
! If precipitation evaporates significantly before reaching next level
3921+
if (k < kbot .and. pdwn(k+1) < 0.5 * pdwn(k) .and. pdwn(k) > pdwn_thr) then
3922+
! Partial re-evaporation anticipated
3923+
alpha = (pdwn(k) - pdwn(k+1)) / pdwn(k)
3924+
alpha = min(one, max(zero, alpha))
3925+
! This will be handled in next level, just note it here if needed
3926+
end if
38853927
end if
38863928
if ( f_washout > zero ) then
38873929
if ( f_rainout > zero ) then
@@ -3897,11 +3939,21 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
38973939
if ( kin ) then
38983940
! -- adjust loss fraction for aerosols
38993941
lossfrac = lossfrac * f_washout / f
3900-
3901-
alpha = abs( qq(k) ) * delz_cm(k) / pdwn(km1)
3902-
alpha = min( one, alpha )
3903-
gain = 0.5 * alpha * dconc(km1)
3904-
wetloss = conc(k) * lossfrac - gain
3942+
gain = zero
3943+
if ( f_rainout > zero ) then
3944+
! Washout from precipitation entering from top
3945+
if (pdwn(km1) > pdwn_thr) then
3946+
alpha = abs( qq(k) ) * delz_cm(k) / pdwn(km1)
3947+
alpha = min( one, alpha )
3948+
gain = 0.5 * alpha * dconc(km1)
3949+
end if
3950+
else
3951+
! Washout from precipitation leaving through bottom
3952+
! No re-evaporation gain in this case (precipitation already present)
3953+
gain = zero
3954+
end if
3955+
wetloss = conc(k) * lossfrac - gain
3956+
wetloss = max(zero, wetloss) ! Prevent negative loss
39053957
! -- skip sulfate
39063958
else
39073959
washed = f_washout * conc(k) + dconc(km1)
@@ -3915,16 +3967,28 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
39153967
end if
39163968
end if
39173969
else
3918-
! -- complete resuspension of rainout + washout from level above
3919-
conc(k) = conc(k) + dconc(km1)
3920-
end if
3921-
3970+
!-- complete resuspension of rainout + washout from level above
3971+
! conc(k) = conc(k) + dconc(km1)
3972+
! Check for actual evaporation: negative formation rate OR decreasing flux
3973+
! Re-evaporation only if precip enters the layer and flux decreases
3974+
if ( (pdwn(km1) > 1.0e-20) .and. (pdwn(k) < pdwn(km1)) ) then
3975+
alpha = max(zero, min(one, (pdwn(km1) - pdwn(k)) / pdwn(km1)))
3976+
conc(k) = conc(k) + alpha * dconc(km1)
3977+
dconc(k) = (one - alpha) * dconc(km1)
3978+
else
3979+
! No evaporation here: just carry scavenged mass downward
3980+
dconc(k) = dconc(km1)
3981+
! conc(k) unchanged
3982+
end if
3983+
end if
3984+
39223985
ftop = f
39233986

39243987
end do
39253988

39263989
! -- surface level
39273990
k = kbot
3991+
km1=k-1
39283992
if (pdwn(km1) > pdwn_thr) then
39293993
f = ftop
39303994
if ( f > zero ) then
@@ -3948,7 +4012,7 @@ subroutine WetRemovalUFS( km, klid, bin_ind, cdt, aero_type, kin, grav, radius,
39484012
aerosol(i,j,k) = conc(k) / dpog(k)
39494013
end do
39504014

3951-
if (associated(fluxout)) fluxout(i,j,bin_ind) = sum(dconc) / dt
4015+
if (associated(fluxout)) fluxout(i,j,bin_ind) = dconc(kbot) / dt
39524016

39534017
end do
39544018
end do

0 commit comments

Comments
 (0)