-
Notifications
You must be signed in to change notification settings - Fork 9
/
microfast.F90
278 lines (227 loc) · 10.4 KB
/
microfast.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
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
! 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 drives the fast microphysics calculations.
!!
!! @author Eric Jensen, Bill McKie
!! @version Sep-1997
subroutine microfast(carma, cstate, iz, scale_threshold, 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
real(kind=f) :: scale_threshold !! Scaling factor for convergence thresholds
integer, intent(inout) :: rc !! return code, negative indicates failure
! Local Variables
integer :: ielem ! element index
integer :: ibin ! bin index
integer :: igas ! gas index
real(kind=f) :: previous_ice(NGAS) ! total ice at the start of substep
real(kind=f) :: previous_liquid(NGAS) ! total liquid at the start of substep
real(kind=f) :: previous_supsatl(NGAS) ! supersaturation wrt ice at the start of substep
real(kind=f) :: previous_supsati(NGAS) ! supersaturation wrt liquid at the start of substep
real(kind=f) :: supsatold
real(kind=f) :: supsatnew
real(kind=f) :: srat
real(kind=f) :: srat1
real(kind=f) :: srat2
real(kind=f) :: s_threshold
1 format(/,'microfast::ERROR - excessive change in supersaturation for ',a,' : iz=',i4,',xc=', &
f7.2,',yc=',f7.2,',srat=',e10.3,',supsatiold=',e10.3,',supsatlold=',e10.3,',supsati=',e10.3, &
',supsatl=',e10.3,',t=',f6.2)
2 format('microfast::ERROR - conditions at beginning of the step : gc=',e10.3,',supsati=',e17.10, &
',supsatl=',e17.10,',t=',f6.2,',d_gc=',e10.3,',d_t=',f6.2)
3 format(/,'microfast::ERROR - excessive change in supersaturation for ',a,' : iz=',i4,',xc=', &
f7.2,',yc=',f7.2,',supsatiold=',e10.3,',supsatlold=',e10.3,',supsati=',e10.3, &
',supsatl=',e10.3,',t=',f6.2)
! Set production and loss rates to zero.
call zeromicro(carma, cstate, iz, rc)
if (rc < RC_OK) return
! Calculate (implicit) particle loss rates for nucleation, growth,
! evaporation, melting, etc.
if (do_grow) then
! Save off the current condensate totals so the gas and latent heating can be
! figured out in a way that conserves mass and energy.
call totalcondensate(carma, cstate, iz, previous_ice, previous_liquid, rc)
if (rc < RC_OK) return
do igas = 1, NGAS
call supersat(carma, cstate, iz, igas, rc)
if (rc < RC_OK) return
previous_supsati(igas) = supsati(iz, igas)
previous_supsatl(igas) = supsatl(iz, igas)
end do
! Have water vapor and sulfuric acid been defined?
if ((igash2o /= 0) .and. (igash2so4 /= 0)) then
! Are both gases avaialble?
if ((gc(iz, igash2o) > 0._f) .and. (gc(iz,igash2so4) > 0._f)) then
! See if any sulfates will form.
call sulfnuc(carma, cstate, iz, rc)
endif
end if
call growevapl(carma, cstate, iz, rc)
if (rc < RC_OK) return
call actdropl(carma, cstate, iz, rc)
if (rc < RC_OK) return
! The Koop, Tabazadeh and Mohler routines provide different schemes for aerosol freezing.
! Only one of these parameterizatons should be active at one time. However, any
! of these routines can be used in conjunction with heterogenous nucleation of glassy
! aerosols.
call freezaerl_tabazadeh2000(carma, cstate, iz, rc)
if (rc < RC_OK) return
call freezaerl_koop2000(carma, cstate, iz, rc)
if (rc < RC_OK) return
call freezaerl_mohler2010(carma, cstate, iz, rc)
if (rc < RC_OK) return
call freezglaerl_murray2010(carma, cstate, iz, rc)
if (rc < RC_OK) return
call hetnucl(carma, cstate, iz, rc)
if (rc < RC_OK) return
call freezdropl(carma, cstate, iz, rc)
if (rc < RC_OK) return
call melticel(carma, cstate, iz, rc)
if (rc < RC_OK) return
endif
if(do_coremasscheck)then
call coremasscheck( carma, cstate, iz, .false.,.true.,.true., "Beforegrowth", rc )
if (rc < RC_OK) return
end if
! Calculate particle production terms and solve for particle
! concentrations at end of time step.
do ielem = 1,NELEM
do ibin = 1,NBIN
if( do_grow )then
call growp(carma, cstate, iz, ibin, ielem, rc)
if (rc < RC_OK) return
call upgxfer(carma, cstate, iz, ibin, ielem, rc)
if (rc < RC_OK) return
endif
call psolve(carma, cstate, iz, ibin, ielem, rc)
if (rc < RC_OK) return
enddo
enddo
! Calculate particle production terms for evaporation;
! gas loss rates and production terms due to particle nucleation;
! growth, and evaporation;
! apply evaporation production terms to particle concentrations;
! and solve for gas concentrations at end of time step.
if (do_grow) then
call evapp(carma, cstate, iz, rc)
if (rc < RC_OK) return
call downgxfer(carma, cstate, iz, rc)
if (rc < RC_OK) return
! NOTE: Not needed because changes in gas concentrations and latent
! heats are now calculated later in gsolve using total condensate.
! call gasexchange(carma, cstate, iz, rc)
! if (rc < RC_OK) return
call downgevapply(carma, cstate, iz, rc)
if (rc < RC_OK) return
call gsolve(carma, cstate, iz, previous_ice, previous_liquid, scale_threshold, rc)
if (rc /=RC_OK) return
endif
call coremasscheck( carma, cstate, iz, .true.,.false.,.false., "AfterGsolve", rc )
if (rc < RC_OK) return
! Update temperature if thermal processes requested
if (do_thermo) then
call tsolve(carma, cstate, iz, scale_threshold, rc)
if (rc /= RC_OK) return
endif
! Update saturation ratios
if (do_grow .or. do_thermo) then
do igas = 1, NGAS
call supersat(carma, cstate, iz, igas, rc)
if (rc < RC_OK) return
! Check to see how much the supersaturation changed during this step. If it
! has changed to much, then cause a retry.
if (t(iz) >= T0) then
supsatold = previous_supsatl(igas)
supsatnew = supsatl(iz,igas)
else
supsatold = previous_supsati(igas)
supsatnew = supsati(iz,igas)
end if
! If ds_threshold is positive, then it indicates that the criteria should
! be based on the percentage change in saturation.
if (ds_threshold(igas) > 0._f) then
if (supsatold >= 1.e-4_f) then
srat1 = abs(supsatnew / supsatold - 1._f)
else
srat1 = 0._f
end if
if (supsatnew >= 1.e-4_f) then
srat2 = abs(supsatold / supsatnew - 1._f)
else
srat2 = 0._f
end if
srat = max(srat1, srat2)
! Don't let one substep change the supersaturation by too much.
if (ds_threshold(igas) > 0._f) then
! if (srat >= ds_threshold(igas)) then
if ((srat >= ds_threshold(igas)) .and. (abs(supsatold - supsatnew) > 0.1_f)) then
if (do_substep) then
if (nretries == maxretries) then
if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, &
xc, yc, srat, previous_supsati(igas), previous_supsatl(igas), &
supsati(iz, igas), supsatl(iz,igas), t(iz)
if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz, igas), &
supsatlold(iz,igas), told(iz), d_gc(iz, igas), d_t(iz)
end if
rc = RC_WARNING_RETRY
else
if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, xc, yc, &
srat, previous_supsati(igas), previous_supsatl(igas), &
supsati(iz, igas), supsatl(iz,igas), t(iz)
end if
end if
end if
! If ds_threshold is negative, then it indicates that the criteria is based
! upon the supersaturation crossing 0, Indicating a shift from growth to
! evaporation and a potential overshoot in the result.
else if (ds_threshold(igas) < 0._f) then
! Adjust the saturation threshold to allow a worse solution if getting a better
! solution is taking too much time. The particular solution at any individual
! point is probably not going to affect the overall result by too much.
s_threshold = abs(ds_threshold(igas))
if (nretries >= (0.8_f * maxretries)) then
s_threshold = 4._f * s_threshold
else if (nretries >= (0.7_f * maxretries)) then
s_threshold = 3.5_f * s_threshold
else if (nretries >= (0.6_f * maxretries)) then
s_threshold = 3._f * s_threshold
else if (nretries >= (0.5_f * maxretries)) then
s_threshold = 2.5_f * s_threshold
else if (nretries >= (0.4_f * maxretries)) then
s_threshold = 2._f * s_threshold
end if
! If the supersaturation changed signs, then we went from growth to evaporation
! or vice versa. Don't let the new supersaturation go too far past 0 in one substep.
! This is to prevent overshooting as growth/evaporation should normally stop when
! the supersaturation is 0.
if (((supsatnew * supsatold) < 0._f) .and. (abs(supsatnew) > s_threshold)) then
if (do_substep) then
if (nretries == maxretries) then
if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, &
xc, yc, previous_supsati(igas), previous_supsatl(igas), &
supsati(iz, igas), supsatl(iz,igas), t(iz)
if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz, igas), &
supsatlold(iz,igas), told(iz), d_gc(iz, igas), d_t(iz)
end if
else
if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, &
xc, yc, previous_supsati(igas), previous_supsatl(igas), &
supsati(iz, igas), supsatl(iz,igas), t(iz)
end if
rc = RC_WARNING_RETRY
end if
end if
end do
endif
! Return to caller with new particle and gas concentrations.
return
end