Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added a new inorganic nucleation subroutine to tomas_mod.F90. #2528

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Added a new inorganic nucleation subroutine to tomas_mod.F90. #2528

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
184 changes: 180 additions & 4 deletions GeosCore/tomas_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,8 @@ MODULE TOMAS_MOD
REAL(fp), SAVE, ALLOCATABLE :: ECSCALE30(:)
REAL(fp), SAVE, ALLOCATABLE :: ECSCALE100(:)

INTEGER :: bin_nuc = 1, tern_nuc = 1 ! Switches for nucleation type.
INTEGER :: bin_nuc = 0, tern_nuc = 0 ! Switches for nucleation type.
INTEGER :: dunn_nuc = 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If these are intended to be constant values, consider declaring these with the PARAMETER attribute (i.e. INTEGER, PARAMETER :: bin_nuc =0, etc. That will provide some extra protection in case someone else tries to modify the value elsewhere.

INTEGER :: act_nuc = 0 ! in BL
INTEGER :: ion_nuc = 0 ! 1 for modgil, 2 for Yu
! (Yu currently broken, JRP 202101)
Expand Down Expand Up @@ -462,9 +463,15 @@ SUBROUTINE AEROPHYS( Input_Opt, State_Chm, State_Grid, State_Met, &
if(tern_nuc == 1) then
write(*,*)' Nucleation: Ternary (Napari ', &
'et al. 2002) and Binary (Vehkamaki et al. 2002)'
else
endif
if (bin_nuc == 1) then
write(*,*)' Nucleation: Binary (Vehkamaki et al. 2002)'
endif

if(dunn_nuc == 1) then
write(*,*)' Nucleation: Inorganic nucleation (Dunne ', &
'et al. 2016)'
endif

firsttime = .false.
endif
Expand Down Expand Up @@ -1768,10 +1775,12 @@ SUBROUTINE getH2SO4conc(Nk, Mk, H2SO4rate, CS, NH3conc, gasConc, &
elseif (bin_nuc.eq.1)then
max_H2SO4conc=1.0e+11_fp*boxvol/1000.e+0_fp*98.e+0_fp/6.022e+23_fp
else
max_H2SO4conc = 1.0e+100_fp
!max_H2SO4conc = 1.0e+100_fp
max_H2SO4conc = 1.0e+11*boxvol/1000.e+0_fp*98.e+0_fp/6.022e+23_fp ! SamO changed this (was 1.0e+100_fp)
endif
else
max_H2SO4conc = 1.0e+100_fp
!max_H2SO4conc = 1.0e+100_fp
max_H2SO4conc = 1.0e+11*boxvol/1000.e+0_fp*98.e+0_fp/6.022e+23_fp ! SamO Changed this
endif

! Checks for when condensation sink is very small
Expand Down Expand Up @@ -1938,6 +1947,7 @@ SUBROUTINE getNucRate(Nk, Mk, Gci,fn,mnuc,nflg, ionrate,surf_area, &
! !LOCAL VARIABLES:
!
double precision nh3ppt ! gas phase ammonia in pptv
double precision nh3moleccm3 ! gas phase ammonia in molec/cm3
double precision h2so4 ! gas phase h2so4 in molec cc-1
double precision gtime ! time to grow to first size bin [s]
double precision ltc, ltc1, ltc2 ! coagulation loss rates [s-1]
Expand All @@ -1964,6 +1974,11 @@ SUBROUTINE getNucRate(Nk, Mk, Gci,fn,mnuc,nflg, ionrate,surf_area, &
double precision h1,h2,h3,h4,h5,h6
double precision dum1,dum2,dum3,dum4 ! dummy variables
double precision rhin,tempin ! rel hum in

DOUBLE PRECISION Jbn2, Jtn2, Jbi2, Jti2 ! SamO
DOUBLE PRECISION Mair ! SamO
DOUBLE PRECISION TOG_LVOC ! SamO
DOUBLE PRECISION fntemp1,fntemp2

real(fp) mydummy
!
Expand All @@ -1979,8 +1994,13 @@ SUBROUTINE getNucRate(Nk, Mk, Gci,fn,mnuc,nflg, ionrate,surf_area, &
h2so4 = Gci(srtso4)/boxvol*1000.e+0_fp/98.e+0_fp*6.022e+23_fp
nh3ppt = Gci(srtnh4)/17.e+0_fp/(boxmass/29.e+0_fp)*1e+12_fp* &
PRES/101325.*273./TEMPTMS ! corrected for pressure (because this should be concentration)

nh3moleccm3 = Gci(srtnh4)/boxvol*1000.e+0_fp/17e+0_fp*6.022e+23_fp ! Changed by SamO
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For consistency, please use the MW of 17.04 for NH3 as defined in the species_database.yml file here:

https://github.com/samuelod-atmos/geos-chem/blob/aedbeb6b4fb2384828131425d831bcd95364d6d0/run/shared/species_database.yml#L3399-L3418

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also please use the AIRMW parameter from Headers/physconstants.F90 for the MW of air. This is defined as 28.9644 g/mol.

https://github.com/samuelod-atmos/geos-chem/blob/main/Headers/physconstants.F90


Mair = 2.69E19*273.15/TEMPTMS*PRES/101325.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use _fp to prevent loss of precision:

Mair = 2.69E19_fp * 273.15_fp / TEMPTMS * PRES / 101325.0_fp

fn = 0.e+0_fp
fntemp1 = 0.e+0_fp !SamO
fntemp2 = 0.e+0_fp !SamO
rnuc = 0.e+0_fp

!print*,'h2so4',h2so4,'nh3ppt',nh3ppt
Expand Down Expand Up @@ -2054,6 +2074,16 @@ SUBROUTINE getNucRate(Nk, Mk, Gci,fn,mnuc,nflg, ionrate,surf_area, &
fn=0.
rnuc=1E-9
nflg=.true.

elseif (nh3moleccm3.gt.0.1.and.dunn_nuc.eq.1) then
tempin=dble(TEMPTMS)
call dunne_inorg_nucl(tempin,ionrate,h2so4,nh3moleccm3, &
Mair, fntemp2,rnuc,Jbn2,Jtn2,Jbi2,Jti2)
! print*,'Jbn2=', Jbn2,'Jtn2=', Jtn2,'Jbi2=', Jbi2,'Jti2=', Jti2
fn = fn + fntemp2
rnuc=0.85d0 ! [nm]
nflg=.true.

else
nflg=.false.
endif
Expand Down Expand Up @@ -2730,6 +2760,7 @@ SUBROUTINE NUCLEATION(Nki,Mki,Gci,Nkf,Mkf,Gcf,fn,fn1,totsulf, &
double precision dt
double precision fn ! nucleation rate of clusters cm-3 s-1
double precision fn1 ! formation rate of particles to first size bin cm-3 s-1
DOUBLE PRECISION fntemp1,fntemp2 !SamO

LOGICAL PDBG ! Signal print for debug
integer lev ! layer of model
Expand All @@ -2746,6 +2777,7 @@ SUBROUTINE NUCLEATION(Nki,Mki,Gci,Nkf,Mkf,Gcf,fn,fn1,totsulf, &
! !LOCAL VARIABLES:
!
double precision nh3ppt ! gas phase ammonia in pptv
double precision nh3moleccm3 ! gas phase ammonia in molec/cm3
double precision h2so4 ! gas phase h2so4 in molec cc-1
double precision rnuc ! critical nucleation radius [nm]
double precision gtime ! time to grow to first size bin [s]
Expand Down Expand Up @@ -2773,6 +2805,9 @@ SUBROUTINE NUCLEATION(Nki,Mki,Gci,Nkf,Mkf,Gcf,fn,fn1,totsulf, &
double precision h1,h2,h3,h4,h5,h6
double precision dum1,dum2,dum3,dum4 ! dummy variables
double precision rhin,tempin ! rel hum in

DOUBLE PRECISION Jbn2, Jtn2, Jbi2, Jti2 ! SamO
DOUBLE PRECISION Mair ! SamO

LOGICAL ERRORSWITCH
!
Expand All @@ -2787,12 +2822,16 @@ SUBROUTINE NUCLEATION(Nki,Mki,Gci,Nkf,Mkf,Gcf,fn,fn1,totsulf, &

errorswitch = .false.

Mair = 2.69E19*273.15/TEMPTMS*PRES/101325.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use _fp to prevent loss of precision:

Mair = 2.69E19_fp * 273.15_fp / TEMPTMS * PRES / 101325.0_fp

h2so4 = Gci(srtso4)/boxvol*1000.e+0_fp/98.e+0_fp*6.022e+23_fp
nh3ppt = Gci(srtnh4)/17.e+0_fp/(boxmass/29.e+0_fp)*1e+12_fp* &
PRES/101325.*273./TEMPTMS ! corrected for pressure (because this should be concentration)
nh3moleccm3 = Gci(srtnh4)/boxvol*1000.e+0_fp/17e+0_fp*6.022e+23_fp ! Changed by SamO
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use 17.04 for the MW of NH3 (see above comment)


fn = 0.e+0_fp
fn1 = 0.e+0_fp
fntemp1 = 0.e+0_fp !SamO
fntemp2 = 0.e+0_fp !SamO
rnuc = 0.e+0_fp
gtime = 0.e+0_fp
nuc_bin = 1 ! added by Pengfei Liu,initialize nuc_bin value
Expand Down Expand Up @@ -2859,6 +2898,18 @@ SUBROUTINE NUCLEATION(Nki,Mki,Gci,Nkf,Mkf,Gcf,fn,fn1,totsulf, &
fn=0.
rnuc=1E-9
endif

if (nh3moleccm3.gt.0.1.and.dunn_nuc.eq.1) then
tempin=dble(TEMPTMS)
call dunne_inorg_nucl(tempin,ionrate,h2so4,nh3moleccm3, &
Mair, fntemp2,rnuc,Jbn2,Jtn2,Jbi2,Jti2)
! print*,'Jbn2=', Jbn2,'Jtn2=', Jtn2,'Jbi2=', Jbi2,'Jti2=', Jti2
fn = fn + fntemp2
rnuc=0.85d0 ! [nm]
endif



if((act_nuc.eq.1).and.(lev.le.7))then
call bl_nucl(h2so4,fn,rnuc)
endif
Expand Down Expand Up @@ -3014,6 +3065,131 @@ END SUBROUTINE NUCLEATION
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
! SUBROUTINE dunne_inorg_nucl

! WRITTEN BY Jeff Pierce, July 2019

! This subroutine calculates the binary and ternary (neutral and ion-mediated)
! nucleation rats from Dunne et al. 2016

! Dunne et al., Science, DOI: 10.1126/science.aaf2649, 2016
!

SUBROUTINE dunne_inorg_nucl(tempi,fioni,cnai,nh3i,Mairi,fn,rnuc, &
Jbn,Jtn,Jbi,Jti)

IMPLICIT NONE

!-----INPUTS------------------------------------------------------------

double precision tempi ! temperature of air [Ka]
double precision fioni ! formation rate of ions [pairs cm-3 s-1[
double precision cnai ! concentration of gas phase sulfuric acid [molec cm-3]
double precision nh3i ! concentration of ammonia [molec cm-3]
double precision Mairi ! concentration of air [molec cm-3]

!-----OUTPUTS-----------------------------------------------------------

double precision fn ! nucleation rate [cm-3 s-1]
double precision rnuc ! critical cluster radius [nm]

!-----INCLUDE FILES-----------------------------------------------------

!-----ARGUMENT DECLARATIONS---------------------------------------------

!-----VARIABLE DECLARATIONS---------------------------------------------

double precision temp ! temperature of air [K]
double precision nh3 ! concentration of gas phase ammonia [molec cm-3]
double precision cna ! concentration of gas phase sulfuric acid [molec cm-3]
double precision fion, Mair
! Simple temperature dependence (comment out this or complex)
!double precision pbn, ubn, vbn, ptn, utn, vtn, pAn, an ! parameters for neutral nuc
!double precision pbi, ubi, vbi, pti, uti, vti, pAi, ai ! parameters for ion nuc
! Complex temperature dependence (comment out this or simple)
double precision pbn, ubn, vbn, wbn, ptn, utn, vtn, wtn, pAn, an ! parameters for neutral nuc
double precision pbi, ubi, vbi, wbi, pti, uti, vti, wti, pAi, ai ! parameters for ion nuc
double precision kbn, ktn, kbi, kti ! rate constants for 4 nucleation types
double precision ionc ! steady-state ion conc. [cm-3]
double precision alpha_ion ! ion recombination coef [cm3/ion/s]
double precision ffn, ffi ! intermediat calculation for ternary
double precision Jbn, Jtn, Jbi, Jti ! nucleation rates from 4 mechanisms
integer i ! counter

!-----ADJUSTABLE PARAMETERS---------------------------------------------

! Simple temperature dependence (comment out this or complex)
!parameter(pbn=3.62, ubn=46.3, vbn=245., ptn=2.82, utn=41.2)
!parameter(vtn=252., pAn=6.76, an=1.3d-4)
!parameter(pbi=2.73, ubi=24.1, vbi=166., pti=2.86, uti=18.3)
!parameter(vti=208., pAi=5.00, ai=5.0d-7)
! Complex temperature dependence (comment out this or simple)
parameter(pbn=3.95, ubn=9.70, vbn=12.6, wbn=-0.00707, ptn=2.89)
parameter(utn=182., vtn=1.20, wtn=-4.19, pAn=8.00, an=1.6d-6)
parameter(pbi=3.37, ubi=-11.5, vbi=25.5, wbi=0.181, pti=3.14)
parameter(uti=-23.8, vti=37.0, wti=0.227, pAi=3.07, ai=0.00485)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Each of the parameters should be followed by d0 to prevent a loss of precision



!-----CODE--------------------------------------------------------------
!nh3i = 1e10 !SamO
temp=tempi !SamO
!temp=278.0
nh3=nh3i*1E-6 ! I think they want it in units of 1E6 molec cm-3
cna=cnai*1E-6
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use 1d-6 to prevent loss of precision

!cna=1E7*1E-6 !SamO
fion=fioni
Mair=Mairi
!Mair=2.4116E+019 !SamO
!fion = 75.0 !SamO

! CALCULATE ION CONCENTRATION
alpha_ion = 6d-8*sqrt(300./temp) + 6d-26*Mair*(300./temp)**4 ! Need Mair in air molec per cm3
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use 300.0d0 to prevent a loss of precision.

ionc = sqrt(fion/alpha_ion) ! assume that ion-ion recombination dominates from conversion with svensmakr, need to verify

! CALCULATE k VALUES
! Simple temperature dependence (comment out this or complex)
!kbn = exp(ubn - vbn*(temp/1000.))
!ktn = exp(utn - vtn*(temp/1000.))
!kbi = exp(ubi - vbi*(temp/1000.))
!kti = exp(uti - vti*(temp/1000.))
! Complex temperature dependence (comment out this or simple)
kbn = exp(ubn - exp(vbn*(temp/1000. - wbn)))
ktn = exp(utn - exp(vtn*(temp/1000. - wtn)))
kbi = exp(ubi - exp(vbi*(temp/1000. - wbi)))
kti = exp(uti - exp(vti*(temp/1000. - wti)))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use 1000.0d0 to prevent a loss of precision.



! CALCULATE f VALUES
if (nh3.gt.1d-10)then
ffn = nh3*cna**ptn/(an + (cna**ptn)/(nh3**pAn))
ffi = nh3*cna**pti/(ai + (cna**pti)/(nh3**pAi))
else
ffn = 0.
ffi = 0.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use 0.0d0 to prevent a loss of precision

endif

! CALCULATE NUCLEATION RATES
Jbn = kbn*cna**pbn
Jtn = ktn*ffn
Jbi = kbi*ionc*cna**pbi
Jti = kti*ionc*ffi
! print*,'Jbn=', Jbn,'Jtn=', Jtn,'Jbi=', Jbi,'Jti=', Jti

fn = Jbn + Jtn + Jbi + Jti ! sum them
!print*,'fn=',fn,cna,fion,nh3

fn = fn*1000.d0 !SamO - added this scalar

! cluster radius
rnuc=0.85d0 ! [nm]

RETURN
END SUBROUTINE dunne_inorg_nucl
!EOC
!------------------------------------------------------------------------------
! GEOS-Chem Global Chemical Transport Model !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: ezcond
!
Expand Down