Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : !! This routine calculates new average particle densities.
6 : !!
7 : !! The particle mass density can change at each time-step due to
8 : !! changes in the core mass fraction.
9 : !!
10 : !! For particles that are hydrophilic and whose particle size changes based
11 : !! upon the relative humidity, and wet radius and density are also calculated.
12 : !! For particles that do not swell, the wet and dry radius and densities are
13 : !! the same.
14 : !!
15 : !! @author Chuck Bardeen Eric Jensen
16 : !! @ version May-2009; Oct-1995
17 : !!
18 : !! @see wetr
19 1050624 : subroutine rhopart(carma, cstate, rc)
20 :
21 : ! types
22 : use carma_precision_mod
23 : use carma_enums_mod
24 : use carma_constants_mod
25 : use carma_types_mod
26 : use carmastate_mod
27 : use carma_mod
28 : use sulfate_utils
29 : use wetr
30 :
31 : implicit none
32 :
33 : type(carma_type), intent(in) :: carma !! the carma object
34 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
35 : integer, intent(inout) :: rc !! return code, negative indicates failure
36 :
37 : ! Local declarations
38 : integer :: iz !! z index
39 : integer :: igroup !! group index
40 : integer :: ibin !! bin index
41 : integer :: iepart !! element in group containing the particle concentration
42 : integer :: jcore
43 : integer :: iecore
44 2101248 : real(kind=f) :: vcore(NBIN)
45 2101248 : real(kind=f) :: mcore(NBIN)
46 : real(kind=f) :: r_ratio
47 : real(kind=f) :: h2o_mass
48 :
49 : 1 format(/,'rhopart::WARNING - core mass > total mass, truncating : iz=',i4,',igroup=',&
50 : i4,',ibin=',i4,',total mass=',e10.3,',core mass=',e10.3,',using rhop=',f9.4)
51 :
52 : ! Calculate hygroscopicity parameter, kappa, for mixed aerosols (Yu et al., JAMES, 2015)
53 1050624 : call hygroscopicity(carma, cstate, rc)
54 :
55 : ! Calculate average particle mass density for each group
56 3151872 : do igroup = 1,NGROUP
57 :
58 : ! Define particle # concentration element index for current group
59 2101248 : iepart = ienconc(igroup) ! element of particle number concentration
60 :
61 150239232 : do iz = 1, NZ
62 :
63 : ! If there are no cores, than the density of the particle is just the density
64 : ! of the element.
65 147087360 : if (ncore(igroup) < 1) then
66 1544417280 : rhop(iz,:,igroup) = rhoelem(:,iepart)
67 :
68 : ! Otherwise, the density changes depending on the amount of core and volatile
69 : ! components.
70 : else ! ncore(igroup) >= 1
71 :
72 : ! Calculate volume of cores and the mass of shell material
73 : ! <vcore> is the volume of core material and <rmshell> is the
74 : ! mass of shell material.
75 1544417280 : vcore(:) = 0._f
76 1544417280 : mcore(:) = 0._f
77 :
78 441262080 : do jcore = 1,ncore(igroup)
79 367718400 : iecore = icorelem(jcore,igroup) ! core element
80 :
81 7722086400 : mcore(:) = mcore(:) + pc(iz,:,iecore)
82 7795630080 : vcore(:) = vcore(:) + pc(iz,:,iecore) / rhoelem(:,iecore)
83 : end do ! jcore = 1,ncore(igroup)
84 :
85 : ! Calculate average density
86 1544417280 : do ibin = 1,NBIN
87 :
88 : ! If there is no core, the the density is that of the volatile element.
89 1544417280 : if (mcore(ibin) == 0._f) then
90 13363 : rhop(iz,ibin,igroup) = rhoelem(ibin,iepart)
91 : else! mcore(ibin) /= 0._f
92 :
93 : ! Since core mass and particle number (i.e. total mass) are advected separately,
94 : ! numerical diffusion during advection can cause problems where the core mass
95 : ! becomes greater than the total mass. To prevent adevction errors from making the
96 : ! group inconsistent, we will truncate core mass if it is larger than the total
97 : ! mass.
98 1470860237 : if (mcore(ibin) > (rmass(ibin,igroup) * pc(iz,ibin,iepart))) then
99 :
100 : ! Calculate the density.
101 5513213 : rhop(iz,ibin,igroup) = mcore(ibin) / vcore(ibin)
102 :
103 : ! NOTE: This error happens a lot, so this error message is commented out
104 : ! by default.
105 : ! if (do_print) write(LUNOPRT,1) iz, igroup, ibin, pc(iz,ibin,iepart)*rmass(ibin,igroup), &
106 : ! mcore(ibin), rhop(iz,ibin,igroup)
107 : ! rc = RC_WARNING
108 :
109 : ! Repair total mass.
110 5513213 : pc(iz,ibin,iepart) = mcore(ibin) / rmass(ibin,igroup)
111 : else
112 0 : rhop(iz,ibin,igroup) = (rmass(ibin,igroup) * pc(iz,ibin,iepart)) / &
113 1465347024 : ((pc(iz,ibin,iepart)*rmass(ibin,igroup) - mcore(ibin))/rhoelem(ibin,iepart) + vcore(ibin))
114 : end if ! mcore(ibin) > (rmass(ibin,igroup) * pc(iz,ibin,iepart))
115 : end if ! mcore(ibin) == 0._f
116 : end do ! ibin = 1,NBIN
117 : end if ! ncore(igroup) < 1
118 :
119 : ! If these particles are hygroscopic and grow in response to the relative
120 : ! humidity, then caclulate a wet radius and wet density. Otherwise the wet
121 : ! and dry radius are the same.
122 :
123 : ! Determine the weight percent of sulfate, and store it for later use.
124 147087360 : if (irhswell(igroup) == I_WTPCT_H2SO4 .or. irhswell(igroup) == I_PETTERS) then
125 147087360 : h2o_mass = gc(iz, igash2o) / zmet(iz)
126 : end if
127 :
128 : ! Loop over particle size bins.
129 3090935808 : do ibin = 1,NBIN
130 :
131 : ! If humidity affects the particle, then determine the equilbirium
132 : ! radius and density based upon the relative humidity.
133 3088834560 : if (irhswell(igroup) == I_WTPCT_H2SO4) then
134 :
135 : ! rlow
136 0 : call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &
137 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
138 0 : h2o_vp=pvapl(iz, igash2o), temp=t(iz))
139 1470873600 : if (rc < 0) return
140 1470873600 :
141 : ! rup
142 : call getwetr(carma, igroup, relhum(iz), rup(ibin,igroup), rup_wet(iz,ibin,igroup), &
143 1470873600 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
144 : h2o_vp=pvapl(iz, igash2o), temp=t(iz))
145 0 : if (rc < 0) return
146 2941747200 :
147 1470873600 : ! r
148 : call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
149 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
150 1470873600 : h2o_vp=pvapl(iz, igash2o), temp=t(iz))
151 : if (rc < 0) return
152 0 :
153 2941747200 : else if (irhswell(igroup) == I_PETTERS) then
154 1470873600 :
155 : call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &
156 1470873600 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc,h2o_mass=h2o_mass, &
157 : h2o_vp=pvapl(iz, igash2o), temp=t(iz), kappa=kappahygro(iz,ibin,igroup))
158 0 : if (rc < 0) return
159 0 :
160 1470873600 : ! rup
161 1470873600 : call getwetr(carma, igroup, relhum(iz), rup(ibin,igroup), rup_wet(iz,ibin,igroup), &
162 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
163 : h2o_vp=pvapl(iz, igash2o),temp=t(iz), kappa=kappahygro(iz,ibin,igroup))
164 1470873600 : if (rc < 0) return
165 :
166 0 : ! r
167 2941747200 : call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
168 1470873600 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc, h2o_mass=h2o_mass, &
169 : h2o_vp=pvapl(iz, igash2o),temp=t(iz), kappa=kappahygro(iz,ibin,igroup))
170 : if (rc < 0) return
171 1470873600 :
172 : else ! I_GERBER and I_FITZGERALD
173 0 :
174 2941747200 : call getwetr(carma, igroup, relhum(iz), rlow(ibin,igroup), rlow_wet(iz,ibin,igroup), &
175 1470873600 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc)
176 : if (rc < 0) return
177 :
178 : ! rup
179 0 : call getwetr(carma, igroup, relhum(iz), rup(ibin,igroup), rup_wet(iz,ibin,igroup), &
180 0 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc)
181 0 : if (rc < 0) return
182 :
183 : ! r
184 0 : call getwetr(carma, igroup, relhum(iz), r(ibin,igroup), r_wet(iz,ibin,igroup), &
185 : rhop(iz,ibin,igroup), rhop_wet(iz,ibin,igroup), rc)
186 0 : if (rc < 0) return
187 0 :
188 : end if
189 : end do ! ibin = 1,NBIN
190 0 : end do ! iz = 1, NZ
191 : end do ! igroup = 1,NGROUP
192 0 :
193 0 : ! Return to caller with new particle number densities.
194 : return
195 : end subroutine rhopart
|
