-
Notifications
You must be signed in to change notification settings - Fork 9
/
freezglaerl_murray2010.F90
136 lines (114 loc) · 6.07 KB
/
freezglaerl_murray2010.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
! Include shortname defintions, so that the F77 code does not have to be modified to
! reference the CARMA structure.
#include "carma_globaer.h"
!! This routine evaluates particle loss rates due to nucleation <rnuclg>:
!! heterogeneous nucleation of glassy aerosols only,.
!!
!! The parameterization of glass aerosols is described in Murray et al.
!! [Nature Geosciences, 2010], and is based upon measurements of the nucleation of
!! citric acid aerosols at cold temperatures.
!!
!! NOTE: This implementation assumes that the aerosol being nucleated is the total
!! aerosol population and not just the fraction of aerosols that are glassy. To
!! account for homogenous freezing of the aerosol population, the routine freezaerl
!! also needs to be called and the overall nucleation rate is the sum of
!! the rates for homogeneous freezing and for heterogenous nucleation.
!!
!! The parameter fglass is the fraction of the total aerosol population that will be
!! in a glassy state for T <= 212K.
!!
!! The loss rates for all particle elements in a particle group are equal.
!!
!! @author Chuck Bardeen, Eric Jensen
!! @version Apr-2010
subroutine freezglaerl_murray2010(carma, cstate, iz, rc)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
type(carma_type), intent(in) :: carma !! the carma object
type(carmastate_type), intent(inout) :: cstate !! the carma state object
integer, intent(in) :: iz !! z index
integer, intent(inout) :: rc !! return code, negative indicates failure
! Local declarations
! Define parameters needed for freezing nucleation calculations.
real(kind=f), parameter :: kice1 = 7.7211e-5_f ! Fit constant from Murray et al.
real(kind=f), parameter :: kice2 = 9.2688e-3_f ! Fit constant from Murray et al.
real(kind=f), parameter :: ssmin = 0.21_f ! Minimum supersaturation for nucleation
real(kind=f), parameter :: ssmax = 0.7_f ! Maximum supersaturation for nucleation
real(kind=f), parameter :: tglass = 212._f ! Maximum temperature for glassy state
real(kind=f), parameter :: fglass = 0.5_f ! Fraction of aerosols that can become glassy
integer :: igas ! gas index
integer :: igroup ! group index
integer :: ibin ! bin index
integer :: iepart ! element for condensing group index
integer :: inuc ! nucleating element index
integer :: ienucto ! index of target nucleation element
integer :: ignucto ! index of target nucleation group
integer :: inucto ! index of target nucleation bin
real(kind=f) :: dfice ! difference in fraction of aerosol nucleated
real(kind=f) :: ssi, ssiold
! Assume success.
rc = RC_OK
! Loop over particle groups.
do igroup = 1,NGROUP
igas = inucgas(igroup) ! condensing gas
iepart = ienconc( igroup ) ! particle number density element
if (igas /= 0) then
! Calculate nucleation loss rates. Do not allow nucleation into
! an evaporating bin.
do inuc = 1, nnuc2elem(iepart)
ienucto = inuc2elem(inuc,iepart)
if (ienucto /= 0) then
ignucto = igelem(ienucto)
! Only compute nucleation rate for glassy aerosol freezing.
if ((iand(inucproc(iepart,ienucto), I_AF_MURRAY_2010) /= 0)) then
! Is it cold enough for aerosols to be in a glassy state.
if (t(iz) <= tglass) then
! Loop over particle bins. Loop from largest to smallest for
! evaluation of index of smallest bin nucleated during time step <inucstep>.
do ibin = NBIN, 1, -1
! Bypass calculation if few particles are present or if it isn't cold enough
! for the aerosols to be present in a glassy state.
if (pconmax(iz,igroup) > FEW_PC) then
! Murray et al. [2010] doesn't really give a nucleation rate. Instead it gives
! a fraction of glassy aerosol particles that have been nucleated as a function
! of ice supersaturation.
!
! Since CARMA really wants to work with rates, use the difference in relative
! humidity and the length of the timestep to come up with an approximation to
! a nucleation rate.
! The supersaturation must be greater than .21 for heterogeneous nucleation to
! commence. The fraction of glassy aerosol nucleated is:
!
! fice = 7.7211e-5 * RHi(%) - 9.2688e-3 for 121 % < RHi < 170 %
!
! To get a pseudo production rate, use
!
! rnuclg = (fice(RHi) - fice(RHi_old)) / dtime
!
ssi = supsati(iz,igas)
ssiold = supsatiold(iz,igas)
if ((ssi >= ssmin) .and. (ssi > ssiold)) then
dfice = kice1 * (1._f + min(ssmax, ssi)) * 100._f - kice2
if (ssiold >= ssmin) then
dfice = dfice - (kice1 * (1._f + min(ssmax, ssiold)) * 100._f - kice2)
endif
! Add the rate of heterogenous freezing to the rate of homogeneous
rnuclg(ibin,igroup,ignucto) = rnuclg(ibin,igroup,ignucto) + fglass * dfice / dtime
endif
endif
enddo
endif
endif
endif
enddo
endif
enddo
! Return to caller with particle loss rates due to nucleation evaluated.
return
end