Skip to content

Commit 9f9843a

Browse files
lizziellizziel
committed
Replace numerical recipe with existing function in shr_spfn_erfc
Signed-off-by: Lizzie Lundgren <[email protected]>
1 parent 1882564 commit 9f9843a

File tree

1 file changed

+76
-99
lines changed

1 file changed

+76
-99
lines changed

‎src/chemistry/geoschem/chemistry.F90‎

Lines changed: 76 additions & 99 deletions
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,6 @@ module chemistry
6565
! Private routines:
6666
!
6767
private :: sect02_mam4
68-
private :: erfc_num_recipes
6968

7069
! Location of valid geoschem_config.yml and species_database.yml
7170
! Use local files in run folder
@@ -4727,33 +4726,6 @@ subroutine chem_emissions( state, cam_in, pbuf )
47274726
ENDDO
47284727

47294728
end subroutine chem_emissions
4730-
!
4731-
! P R E S C R I B E A E R O S O L D I S T R I B U T I O N
4732-
!
4733-
! Based on code from Feng et al., 2021 GMD (WRF-GC v2.0), by Xu Feng et al.
4734-
! in module_diag_aero_size_info.F, originally based from WRF-Chem.
4735-
!
4736-
! Reference:
4737-
! Feng, X., Lin, H., Fu, T.-M., Sulprizio, M. P., Zhuang, J., Jacob, D. J., Tian, H., Ma, Y., Zhang, L., Wang, X., Chen, Q., and Han, Z.: WRF-GC (v2.0): online two-way coupling of WRF (v3.9.1.1) and GEOS-Chem (v12.7.2) for modeling regional atmospheric chemistry–meteorology interactions, Geosci. Model Dev., 14, 3741–3768, https://doi.org/10.5194/gmd-14-3741-2021, 2021.
4738-
!
4739-
4740-
real(8) function erfc_num_recipes( x )
4741-
!
4742-
! from press et al, numerical recipes, 1990, page 164
4743-
!
4744-
implicit none
4745-
real(r8) :: x, erfc_dbl, dum, t, z
4746-
z = abs(x)
4747-
t = 1.0_r8/(1.0_r8 + 0.5_r8*z)
4748-
dum = ( -z*z - 1.26551223_r8 + t*(1.00002368_r8 + t*(0.37409196_r8 + &
4749-
t*(0.09678418_r8 + t*(-0.18628806_r8 + t*(0.27886807_r8 + &
4750-
t*(-1.13520398_r8 + &
4751-
t*(1.48851587_r8 + t*(-0.82215223_r8 + t*0.17087277_r8 )))))))))
4752-
erfc_dbl = t * exp(dum)
4753-
if (x .lt. 0.0_r8) erfc_dbl = 2.0_r8 - erfc_dbl
4754-
erfc_num_recipes = erfc_dbl
4755-
return
4756-
end function erfc_num_recipes
47574729

47584730
! sect02_mam4 is based off sect02_new in WRF-GC, which is based off
47594731
! sect02 in WRF-Chem chem/module_optical_averaging.F.
@@ -4762,76 +4734,81 @@ end function erfc_num_recipes
47624734
! prog calculates mass and number for each section.
47634735
subroutine sect02_mam4(dgnum_um, sigmag, duma, nbin, dlo_sect, dhi_sect, &
47644736
xnum_sect, xmas_sect)
4765-
! INPUT PARAMETERS:
4766-
! dgnum_um *diameter* geometric mean of log-normal distribution [um]
4767-
! sigmag geometric standard deviation of log-normal dist. [unitless]
4768-
! duma 1.0 ?
4769-
! nbin # of target bins (wrf-gc = 4, MAM4 = 3) [count]
4770-
! dlo_sect(nbin) low diameter limit (wrf-gc = 0.0390625) [um]
4771-
! dhi_sect(nbin) high diameter limit (wrf-gc = 10.0) [um]
4772-
4773-
! OUTPUT PARAMETERS:
4774-
! xnum_sect(nbin) aerosol number per bin, ratio of total [unitless]
4775-
! xmas_sect(bin) aerosol mass per bin, ratio of total [unitless]
4776-
4777-
implicit none
4778-
real(8), dimension(nbin), intent(out) :: xnum_sect, xmas_sect
4779-
integer :: n, nbin
4780-
real(8) :: dgnum, dgnum_um, dhi, &
4781-
dlo, duma, dumfrac, &
4782-
dx, sigmag, &
4783-
sx, sxroot2, thi, tlo, x0, x3, &
4784-
xhi, xlo, xmtot, xntot
4785-
real(8), intent(in) :: dlo_sect(nbin), dhi_sect(nbin)
4786-
real(8) :: my_dlo_sect(nbin), my_dhi_sect(nbin)
4787-
real(8) :: pi
4788-
parameter (pi = 3.141592653589_r8)
4789-
4790-
xmtot = duma
4791-
xntot = duma
4792-
4793-
! Compute bins based on number of bins. Originally sect02_new.
4794-
! For MAM4, we prescribe the bin ranges as well.
4795-
! dlo = dlo_um*1.0E-4_r8
4796-
! dhi = dhi_um*1.0E-4_r8
4797-
! xlo = log( dlo )
4798-
! xhi = log( dhi )
4799-
! dx = (xhi - xlo)/nbin
4800-
! do n = 1, nbin
4801-
! dlo_sect(n) = exp( xlo + dx*(n-1) )
4802-
! dhi_sect(n) = exp( xlo + dx*n )
4803-
! end do
4804-
4805-
! dlo_sect and dhi_sect have to be scaled by 1e-4
4806-
! in order to fit parameters in the above calculation, if they are prescribed.
4807-
4808-
my_dlo_sect(:) = dlo_sect(:) * 1.0e-4_r8
4809-
my_dhi_sect(:) = dhi_sect(:) * 1.0e-4_r8
4810-
4811-
dgnum = dgnum_um*1.0E-4_r8
4812-
sx = log( sigmag )
4813-
x0 = log( dgnum )
4814-
x3 = x0 + 3.0_r8*sx*sx
4815-
sxroot2 = sx * sqrt( 2.0_r8 )
4816-
do n = 1, nbin
4817-
xlo = log( my_dlo_sect(n) )
4818-
xhi = log( my_dhi_sect(n) )
4819-
tlo = (xlo - x0)/sxroot2
4820-
thi = (xhi - x0)/sxroot2
4821-
if (tlo .le. 0.0_r8) then
4822-
dumfrac = 0.5_r8*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
4823-
else
4824-
dumfrac = 0.5_r8*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
4825-
end if
4826-
xnum_sect(n) = xntot*dumfrac
4827-
tlo = (xlo - x3)/sxroot2
4828-
thi = (xhi - x3)/sxroot2
4829-
if (tlo .le. 0.0_r8) then
4830-
dumfrac = 0.5_r8*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
4831-
else
4832-
dumfrac = 0.5_r8*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
4833-
endif
4834-
xmas_sect(n) = xmtot*dumfrac
4835-
enddo
4737+
4738+
use shr_spfn_mod, only: erfc => shr_spfn_erfc
4739+
4740+
! INPUT PARAMETERS:
4741+
! dgnum_um *diameter* geometric mean of log-normal distribution [um]
4742+
! sigmag geometric standard deviation of log-normal dist. [unitless]
4743+
! duma 1.0 ?
4744+
! nbin # of target bins (wrf-gc = 4, MAM4 = 3) [count]
4745+
! dlo_sect(nbin) low diameter limit (wrf-gc = 0.0390625) [um]
4746+
! dhi_sect(nbin) high diameter limit (wrf-gc = 10.0) [um]
4747+
4748+
! OUTPUT PARAMETERS:
4749+
! xnum_sect(nbin) aerosol number per bin, ratio of total [unitless]
4750+
! xmas_sect(bin) aerosol mass per bin, ratio of total [unitless]
4751+
4752+
implicit none
4753+
real(8), dimension(nbin), intent(out) :: xnum_sect, xmas_sect
4754+
integer :: n, nbin
4755+
real(8) :: dgnum, dgnum_um, dhi, &
4756+
dlo, duma, dumfrac, &
4757+
dx, sigmag, &
4758+
sx, sxroot2, thi, tlo, x0, x3, &
4759+
xhi, xlo, xmtot, xntot
4760+
real(8), intent(in) :: dlo_sect(nbin), dhi_sect(nbin)
4761+
real(8) :: my_dlo_sect(nbin), my_dhi_sect(nbin)
4762+
real(8) :: pi
4763+
parameter (pi = 3.141592653589_r8)
4764+
4765+
xmtot = duma
4766+
xntot = duma
4767+
4768+
! Compute bins based on number of bins. Originally sect02_new.
4769+
! For MAM4, we prescribe the bin ranges as well.
4770+
! dlo = dlo_um*1.0E-4_r8
4771+
! dhi = dhi_um*1.0E-4_r8
4772+
! xlo = log( dlo )
4773+
! xhi = log( dhi )
4774+
! dx = (xhi - xlo)/nbin
4775+
! do n = 1, nbin
4776+
! dlo_sect(n) = exp( xlo + dx*(n-1) )
4777+
! dhi_sect(n) = exp( xlo + dx*n )
4778+
! end do
4779+
4780+
! dlo_sect and dhi_sect have to be scaled by 1e-4
4781+
! in order to fit parameters in the above calculation, if they are prescribed.
4782+
4783+
my_dlo_sect(:) = dlo_sect(:) * 1.0e-4_r8
4784+
my_dhi_sect(:) = dhi_sect(:) * 1.0e-4_r8
4785+
4786+
dgnum = dgnum_um*1.0E-4_r8
4787+
sx = log( sigmag )
4788+
x0 = log( dgnum )
4789+
x3 = x0 + 3.0_r8*sx*sx
4790+
sxroot2 = sx * sqrt( 2.0_r8 )
4791+
do n = 1, nbin
4792+
xlo = log( my_dlo_sect(n) )
4793+
xhi = log( my_dhi_sect(n) )
4794+
tlo = (xlo - x0)/sxroot2
4795+
thi = (xhi - x0)/sxroot2
4796+
if (tlo .le. 0.0_r8) then
4797+
dumfrac = 0.5_r8*( erfc(-thi) - erfc(-tlo) )
4798+
else
4799+
dumfrac = 0.5_r8*( erfc(tlo) - erfc(thi) )
4800+
end if
4801+
xnum_sect(n) = xntot*dumfrac
4802+
tlo = (xlo - x3)/sxroot2
4803+
thi = (xhi - x3)/sxroot2
4804+
if (tlo .le. 0.0_r8) then
4805+
dumfrac = 0.5_r8*( erfc(-thi) - erfc(-tlo) )
4806+
else
4807+
dumfrac = 0.5_r8*( erfc(tlo) - erfc(thi) )
4808+
endif
4809+
xmas_sect(n) = xmtot*dumfrac
4810+
enddo
4811+
48364812
end subroutine sect02_mam4
4813+
48374814
end module chemistry

0 commit comments

Comments
 (0)