forked from ESCOMP/CARMA_base
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgrowevapl.F90
256 lines (215 loc) · 9.06 KB
/
growevapl.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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
! 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 evaluate particle loss rates due to condensational
!! growth and evaporation for all condensing gases.
!!
!! The loss rates for each group are <growlg> and <evaplg>.
!!
!! Units are [s^-1].
!!
!! @author Andy Ackerman
!! @version Dec-1995
subroutine growevapl(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
integer :: igroup
integer :: iepart
integer :: igas
integer :: ibin
integer :: isol
integer :: nother
integer :: ieoth_rel
integer :: ieoth_abs
integer :: jother
real(kind=f) :: argsol
real(kind=f) :: othermtot
real(kind=f) :: condm
real(kind=f) :: akas
real(kind=f) :: expon
real(kind=f) :: g0
real(kind=f) :: g1
real(kind=f) :: g2
real(kind=f) :: ss
real(kind=f) :: pvap
real(kind=f) :: dpc
real(kind=f) :: dpc1
real(kind=f) :: dpcm1
real(kind=f) :: rat1
real(kind=f) :: rat2
real(kind=f) :: rat3
real(kind=f) :: rat4
real(kind=f) :: ratt1
real(kind=f) :: ratt2
real(kind=f) :: ratt3
real(kind=f) :: den1
real(kind=f) :: test1
real(kind=f) :: test2
real(kind=f) :: x
integer :: ieother(NELEM)
real(kind=f) :: otherm(NELEM)
real(kind=f) :: dela(NBIN)
real(kind=f) :: delma(NBIN)
real(kind=f) :: aju(NBIN)
real(kind=f) :: ar(NBIN)
real(kind=f) :: al(NBIN)
real(kind=f) :: a6(NBIN)
real(kind=f) :: dmdt(NBIN)
real(kind=f) :: growlg_max
do igroup = 1,NGROUP
! element of particle number concentration
iepart = ienconc(igroup)
! condensing gas
igas = igrowgas(iepart)
if (igas .ne. 0) then
! Only valid for condensing liquid water and sulfric acid currently.
if ((igas /= igash2o) .and. (igas .ne. igash2so4)) then
if (do_print) write(LUNOPRT,*) 'growevapl::ERROR - Invalid gas (', igas, ').'
rc = -1
return
endif
! Treat condensation of gas <igas> to/from particle group <igroup>.
!
! Bypass calculation if few particles are present
if( pconmax(iz,igroup) .gt. FEW_PC )then
do ibin = 1,NBIN-1
! Determine the growth rate (dmdt). This calculation may take into account
! radiative effects on the particle which can affect the growth rates.
call pheat(carma, cstate, iz, igroup, iepart, ibin, igas, dmdt(ibin), rc)
enddo ! ibin = 1,NBIN-1
! Now calculate condensation/evaporation production and loss rates.
! Use Piecewise Polynomial Method [Colela and Woodard, J. Comp. Phys.,
! 54, 174-201, 1984]
!
! First, use cubic fits to estimate concentration values at bin
! boundaries
do ibin = 2,NBIN-1
dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
dpc1 = pc(iz,ibin+1,iepart) / dm(ibin+1,igroup)
dpcm1 = pc(iz,ibin-1,iepart) / dm(ibin-1,igroup)
ratt1 = pratt(1,ibin,igroup)
ratt2 = pratt(2,ibin,igroup)
ratt3 = pratt(3,ibin,igroup)
dela(ibin) = ratt1 * ( ratt2*(dpc1-dpc) + ratt3*(dpc-dpcm1) )
delma(ibin) = 0._f
if( (dpc1-dpc)*(dpc-dpcm1) .gt. 0._f ) &
delma(ibin) = min( abs(dela(ibin)), 2._f*abs(dpc-dpc1), &
2._f*abs(dpc-dpcm1) ) * sign(1._f, dela(ibin))
enddo ! ibin = 2,NBIN-2
do ibin = 2,NBIN-2
dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
dpc1 = pc(iz,ibin+1,iepart) / dm(ibin+1,igroup)
dpcm1 = pc(iz,ibin-1,iepart) / dm(ibin-1,igroup)
rat1 = prat(1,ibin,igroup)
rat2 = prat(2,ibin,igroup)
rat3 = prat(3,ibin,igroup)
rat4 = prat(4,ibin,igroup)
den1 = pden1(ibin,igroup)
! <aju(ibin)> is the estimate for concentration (dn/dm) at bin
! boundary <ibin>+1/2.
aju(ibin) = dpc + rat1*(dpc1-dpc) + 1._f/den1 * &
( rat2*(rat3-rat4)*(dpc1-dpc) - &
dm(ibin,igroup)*rat3*delma(ibin+1) + &
dm(ibin+1,igroup)*rat4*delma(ibin) )
enddo ! ibin = 2,NBIN-2
! Now construct polynomial functions in each bin
do ibin = 3,NBIN-2
al(ibin) = aju(ibin-1)
ar(ibin) = aju(ibin)
enddo
! Use linear functions in first two and last two bins
if( NBIN .gt. 1 )then
ibin = NBIN
ar(2) = aju(2)
al(2) = pc(iz,1,iepart)/dm(1,igroup) + &
palr(1,igroup) * &
(pc(iz,2,iepart)/dm(2,igroup)- &
pc(iz,1,iepart)/dm(1,igroup))
ar(1) = al(2)
al(1) = pc(iz,1,iepart)/dm(1,igroup) + &
palr(2,igroup) * &
(pc(iz,2,iepart)/dm(2,igroup)- &
pc(iz,1,iepart)/dm(1,igroup))
al(ibin-1) = aju(ibin-2)
ar(ibin-1) = pc(iz,ibin-1,iepart)/dm(ibin-1,igroup) + &
palr(3,igroup) * &
(pc(iz,ibin,iepart)/dm(ibin,igroup)- &
pc(iz,ibin-1,iepart)/dm(ibin-1,igroup))
al(ibin) = ar(ibin-1)
ar(ibin) = pc(iz,ibin-1,iepart)/dm(ibin-1,igroup) + &
palr(4,igroup) * &
(pc(iz,ibin,iepart)/dm(ibin,igroup)- &
pc(iz,ibin-1,iepart)/dm(ibin-1,igroup))
endif
! Next, ensure that polynomial functions do not deviate beyond the
! range [<al(ibin)>,<ar(ibin)>]
do ibin = 1,NBIN
dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
if( (ar(ibin)-dpc)*(dpc-al(ibin)) .le. 0._f )then
al(ibin) = dpc
ar(ibin) = dpc
endif
test1 = (ar(ibin)-al(ibin))*(dpc - 0.5_f*(al(ibin)+ar(ibin)))
test2 = 1._f/6._f*(ar(ibin)-al(ibin))**2
if( test1 .gt. test2 )then
al(ibin) = 3._f*dpc - 2._f*ar(ibin)
elseif( test1 .lt. -test2 )then
ar(ibin) = 3._f*dpc - 2._f*al(ibin)
endif
enddo
! Lastly, calculate fluxes across each bin boundary.
!
! Use upwind advection when courant number > 1.
do ibin = 1,NBIN
dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
dela(ibin) = ar(ibin) - al(ibin)
a6(ibin) = 6._f * ( dpc - 0.5_f*(ar(ibin)+al(ibin)) )
enddo
do ibin = 1,NBIN-1
if( dmdt(ibin) .gt. 0._f .and. &
pc(iz,ibin,iepart) .gt. SMALL_PC )then
x = dmdt(ibin)*dtime/dm(ibin,igroup)
if( x .lt. 1._f )then
growlg(ibin,igroup) = dmdt(ibin)/pc(iz,ibin,iepart) &
* ( ar(ibin) - 0.5*dela(ibin)*x + &
(x/2._f - x**2/3._f)*a6(ibin) )
else
growlg(ibin,igroup) = dmdt(ibin) / dm(ibin,igroup)
endif
elseif( dmdt(ibin) .lt. 0._f .and. &
pc(iz,ibin+1,iepart) .gt. SMALL_PC )then
x = -dmdt(ibin)*dtime/dm(ibin+1,igroup)
if( x .lt. 1._f )then
evaplg(ibin+1,igroup) = -dmdt(ibin)/ &
pc(iz,ibin+1,iepart) &
* ( al(ibin+1) + 0.5_f*dela(ibin+1)*x + &
(x/2._f - (x**2)/3._f)*a6(ibin+1) )
else
evaplg(ibin+1,igroup) = -dmdt(ibin) / dm(ibin+1,igroup)
endif
! Boundary conditions: for evaporation out of first bin (with cores),
! use evaporation rate from second bin.
! if( ibin .eq. 1 .and. ncore(igroup) .gt. 0 )then
if( ibin .eq. 1)then
evaplg(1,igroup) = -dmdt(1) / dm(1,igroup)
endif
endif
enddo ! ibin = 1,NBIN-1
endif ! (pconmax .gt. FEW_PC)
endif ! (igas = igrowgas(ielem)) .ne. 0
enddo ! igroup = 1,NGROUP
! Return to caller with particle loss rates for growth and evaporation
! evaluated.
return
end