From 0ba6e6958a0479ca50188d6f97cee48d4269b5e3 Mon Sep 17 00:00:00 2001 From: Mike Manyin Date: Fri, 9 Jan 2026 17:25:50 -0500 Subject: [PATCH 1/2] Added a version of chem_settling that does not involve w_c --- Shared/Chem_Shared/Chem_SettlingMod.F90 | 303 ++++++++++++++++++++++++ 1 file changed, 303 insertions(+) diff --git a/Shared/Chem_Shared/Chem_SettlingMod.F90 b/Shared/Chem_Shared/Chem_SettlingMod.F90 index c0595921..37073c9f 100644 --- a/Shared/Chem_Shared/Chem_SettlingMod.F90 +++ b/Shared/Chem_Shared/Chem_SettlingMod.F90 @@ -30,6 +30,7 @@ module Chem_SettlingMod ! PUBLIC Chem_Settling + PUBLIC Chem_Settling2 PUBLIC Chem_SettlingSimple PUBLIC Chem_CalcVsettle @@ -350,6 +351,308 @@ subroutine Chem_Settling ( i1, i2, j1, j2, km, nbeg, nend, nbins, flag, & end subroutine Chem_Settling +! This is Chem_Settling with just 1 3D dataset, and no w_c +! nbeg, nend, nbins - removed +! delp, rh, data3d - added + subroutine Chem_Settling2 ( i1, i2, j1, j2, km, flag, & + radiusInp, rhopInp, cdt, delp, rh, tmpu, rhoa, & + hsurf, hghte, data3d, fluxout, rc, & + vsettleOut, correctionMaring ) + +! !USES: + + implicit NONE + +! !INPUT PARAMETERS: + + integer, intent(in) :: i1, i2, j1, j2, km + integer, intent(in) :: flag ! flag to control particle swelling (see note) + real, intent(in) :: cdt + real :: radiusInp, rhopInp + real, pointer, dimension(:,:) :: hsurf + real, pointer, dimension(:,:,:) :: delp, rh + real, pointer, dimension(:,:,:) :: tmpu, rhoa, hghte + +! !OUTPUT PARAMETERS: + + real, pointer, dimension(:,:,:) :: data3d ! 3D Chemical tracer fields, + type(Chem_Array), pointer :: fluxout ! Mass lost by settling + ! to surface, kg/m2/s + integer, intent(out) :: rc ! Error return code: + ! 0 - all is well + ! 1 - +! Optionally output the settling velocity calculated + type(Chem_Array), pointer, optional :: vsettleOut + +! Optionally correct the settling velocity following Maring et al, 2003 + logical, optional, intent(in) :: correctionMaring + + character(len=*), parameter :: myname = 'Settling' + +! !DESCRIPTION: Gravitational settling of aerosol between vertical +! layers. Assumes input radius in [m] and density (rhop) +! in [kg m-3]. If flag is set, use the Fitzgerald 1975 (flag = 1) +! or Gerber 1985 (flag = 2) parameterization to update the +! particle radius for the calculation (local variables radius +! and rhop). +! +! !REVISION HISTORY: +! +! 17Sep2004 Colarco Strengthen sedimentation flux out at surface +! by setting removal to be valid from middle of +! surface layer +! 06Nov2003 Colarco Based on Ginoux +! 23Jan2003 da Silva Standardization +! +!EOP +!------------------------------------------------------------------------- + +! !Local Variables + integer :: i, j, k, iit + real, parameter :: rhow = 1000. ! Density of water [kg m-3] + real :: vsettle(i1:i2,j1:j2,km) ! fall speed [m s-1] + real*8 :: cmass_before(i1:i2,j1:j2) + real*8 :: cmass_after(i1:i2,j1:j2) + real :: diff_coef ! Brownian diffusion coefficient [m2 s-1] +! real*8 :: qdel(i1:i2,j1:j2), qsrc(i1:i2,j1:j2), dp(i1:i2,j1:j2), dpm1(i1:i2,j1:j2) + +real*8 :: qdel, qsrc, dp, dpm1 + + real :: dz(i1:i2,j1:j2,km) ! layer thickness [m] + real*8 :: dzd(i1:i2,j1:j2,km), vsd(i1:i2,j1:j2,km), qa(i1:i2,j1:j2,km) + +! The following parameters relate to the swelling of seasalt like particles +! following Fitzgerald, Journal of Applied Meteorology, 1975. + real, parameter :: epsilon = 1. ! soluble fraction of deliqeuscing particle + real, parameter :: alphaNaCl = 1.35 + real :: alpha, alpha1, alpharat, beta, theta, f1, f2 + +! parameter from Gerber 1985 (units require radius in cm, see rcm) + real :: rcm + real, parameter :: c1=0.7674, c2=3.079, c3=2.573e-11, c4=-1.424 +! parameters for ammonium sulfate + real, parameter :: SU_c1=0.4809, SU_c2=3.082, SU_c3=3.110e-11, SU_c4=-1.428 + + +! parameters from Maring et al, 2003 + real, parameter :: v_upwardMaring = 0.33e-2 ! upward velocity, [m s-1] + real, parameter :: diameterMaring = 7.30e-6 ! particle diameter, [m] + +! + real :: sat, rrat + real :: radius, rhop ! particle radius and density passed to + ! fall velocity calculation + real :: minTime, qmin, qmax + integer :: nSubSteps, dk, ijl + real*8 :: dt_settle, g + + rc = 0 + + ijl = ( i2 - i1 + 1 ) * ( j2 - j1 + 1 ) + g = grav + +! Handle the fact that hghte may be in the range [1,km+1] or [0,km] +! ----------------------------------------------------------------- + dk = lbound(hghte,3) - 1 ! This is either 0 or 1 + +! Layer thickness from hydrostatic equation + k = km + dz(:,:,k) = hghte(:,:,k+dk)-hsurf(:,:) + do k = km-1, 1, -1 + dz(:,:,k) = hghte(:,:,k+dk) - hghte(:,:,k+dk+1) + enddo + dzd = dz + + qa = data3d + + radius = radiusInp + rhop = rhopInp + +! Reset a (large) minimum time to cross a grid cell in settling + minTime = cdt + + if( associated(fluxout%data2d) ) fluxout%data2d(i1:i2,j1:j2) = 0.0 + cmass_before(:,:) = 0.d0 + cmass_after(:,:) = 0.d0 + +! If radius le 0 then get out of loop + if(radius .gt. 0.) then + + do k = 1, km + do j = j1, j2 + do i = i1, i2 + +! Find the column dry mass before sedimentation + cmass_before(i,j) = cmass_before(i,j) + qa(i,j,k)/g * delp(i,j,k) + +! Adjust the particle size for relative humidity effects + sat = max(rh(i,j,k),tiny(1.0)) ! to avoid zero FPE + +! Fitzgerald + if(flag .eq. 1 .and. sat .ge. 0.80) then +! parameterization blows up for RH > 0.995, so set that as max +! rh needs to be scaled 0 - 1 + sat = min(0.995,sat) +! Calculate the alpha and beta parameters for the wet particle +! relative to amonium sulfate + beta = exp( (0.00077*sat) / (1.009-sat) ) + if(sat .le. 0.97) then + theta = 1.058 + else + theta = 1.058 - (0.0155*(sat-0.97)) /(1.02-sat**1.4) + endif + alpha1 = 1.2*exp( (0.066*sat) / (theta-sat) ) + f1 = 10.2 - 23.7*sat + 14.5*sat**2. + f2 = -6.7 + 15.5*sat - 9.2*sat**2. + alpharat = 1. - f1*(1.-epsilon) - f2*(1.-epsilon**2.) + alpha = alphaNaCl * (alpha1*alpharat) +! radius is the radius of the wet particle + radius = alpha * radiusInp**beta + rrat = (radiusInp/radius)**3. + rhop = rrat*rhopInp + (1.-rrat)*rhow + elseif(flag .eq. 2) then ! Gerber + sat = min(0.995,sat) + rcm = radiusInp*100. + radius = 0.01 * ( c1*rcm**c2 / (c3*rcm**c4-alog10(sat)) & + + rcm**3.)**(1./3.) + rrat = (radiusInp/radius)**3. + rhop = rrat*rhopInp + (1.-rrat)*rhow + elseif(flag .eq. 3) then +! Gerber parameterization for Ammonium Sulfate + sat = min(0.995,sat) + rcm = radiusInp*100. + radius = 0.01 * ( SU_c1*rcm**SU_c2 / (SU_c3*rcm**SU_c4-alog10(sat)) & + + rcm**3.)**(1./3.) + rrat = (radiusInp/radius)**3. + rhop = rrat*rhopInp + (1.-rrat)*rhow + elseif(flag .eq. 4) then +! Petters and Kreidenweis (ACP2007) parameterization + sat = min(0.99,sat) + radius = (radiusInp**3 * (1+1.19*sat/(1-sat)))**(1./3.) + rrat = (radiusInp/radius)**3 + rhop = rrat*rhopInp + (1.-rrat)*rhow + endif + +! Calculate the settling velocity + call Chem_CalcVsettle(radius, rhop, rhoa(i,j,k), & + tmpu(i,j,k), diff_coef, vsettle(i,j,k)) + end do + end do + end do + + if(present(correctionMaring)) then + if ((correctionMaring) .and. (radiusInp .le. (0.5*diameterMaring))) then + vsettle = max(1.0e-9, vsettle - v_upwardMaring) + endif + endif + + vsd = vsettle + + if(present(vsettleOut)) then + vsettleOut%data3d = vsettle + endif + +! Determine global max/min time to cross grid cell +! call pmaxmin ( 'Chem_Settling: dt', dz(i1:i2,j1:j2,1:km)/vsettle(i1:i2,j1:j2,1:km), & +! qmin, qmax, ijl, km, 0. ) +! minTime = min(minTime,qmin) + + +! Loop over sub-timestep + do j = j1, j2 + do i = i1, i2 + ! Determine global min time to cross grid cell + qmin = minval(dz(i,j,:)/vsettle(i,j,:)) + minTime = min(cdt,qmin) + ! Now, how many iterations do we need to do? + if ( minTime < 0 ) then + nSubSteps = 0 + else if(minTime .ge. cdt) then + nSubSteps = 1 + dt_settle = cdt + else + nSubSteps = cdt/minTime+1 + dt_settle = cdt/nSubSteps + endif + + do iit = 1, nSubSteps +! Try a simple forward Euler scheme + qdel = qa(i,j,1)*dt_settle*vsd(i,j,1)/dzd(i,j,1) + qa(i,j,1) = qa(i,j,1) - qdel + + do k = 2, km + dp = delp(i,j,k) + dpm1 = delp(i,j,k-1) + qsrc = qdel * dpm1 / dp + qdel = qa(i,j,k)*dt_settle*vsd(i,j,k)/dzd(i,j,k) + qa(i,j,k) = qa(i,j,k) - qdel + qsrc + end do + end do !iit + end do + end do + + + + +#if 0 +! Now, how many iterations do we need to do? + if ( minTime < 0 ) then + nSubSteps = 0 + call mpout_log(myname,'no Settling because minTime = ', minTime ) + else if(minTime .ge. cdt) then + nSubSteps = 1 + dt_settle = cdt + else + nSubSteps = cdt/minTime+1 + dt_settle = cdt/nSubSteps + endif + +! Loop over sub-timestep + do iit = 1, nSubSteps + +! Try a simple forward Euler scheme + + qdel = qa(i1:i2,j1:j2,1)*dt_settle*vsd(i1:i2,j1:j2,1)/dzd(i1:i2,j1:j2,1) + qa(i1:i2,j1:j2,1) = qa(i1:i2,j1:j2,1) - qdel + + do k = 2, km + dp = delp(i1:i2,j1:j2,k) + dpm1 = delp(i1:i2,j1:j2,k-1) + qsrc = qdel * dpm1 / dp + qdel = qa(i1:i2,j1:j2,k)*dt_settle*vsd(i1:i2,j1:j2,k)/dzd(i1:i2,j1:j2,k) + qa(i1:i2,j1:j2,k) = qa(i1:i2,j1:j2,k) - qdel + qsrc + enddo + +! An alternative accumulator approach to computing the outgoing flux +! if( associated(fluxout%data2d) ) then +! fluxout%data2d = fluxout%data2d + qdel * pdog/grav / dt_settle +! endif + + end do ! iit +#endif + + +! Find the column dry mass after sedimentation and thus the loss flux + do k = 1, km + do j = j1, j2 + do i = i1, i2 + cmass_after(i,j) = cmass_after(i,j) + qa(i,j,k)/ g * delp(i,j,k) + enddo + enddo + enddo + + if( associated(fluxout%data2d) ) then + fluxout%data2d(i1:i2,j1:j2) & + = (cmass_before(i1:i2,j1:j2) - cmass_after(i1:i2,j1:j2))/cdt + endif + + data3d = qa + + endif ! radius .gt. 0 + + end subroutine Chem_Settling2 + + From fbaf98eca66f11193d96f6222ed222050ce73704 Mon Sep 17 00:00:00 2001 From: Mike Manyin Date: Fri, 9 Jan 2026 17:55:17 -0500 Subject: [PATCH 2/2] edit CHANGELOG --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index bee63cb2..57326183 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added + +- Added Chem_Settling2 routine, a version of Chem_Settling that does not require w_c (Chem_Bundle) + ### Removed ### Changed ### Fixed