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 drives the fast microphysics calculations.
6 : !!
7 : !! @author Eric Jensen, Bill McKie
8 : !! @version Sep-1997
9 1418703438 : subroutine microfast(carma, cstate, iz, scale_threshold, rc)
10 :
11 : ! types
12 : use carma_precision_mod
13 : use carma_enums_mod
14 : use carma_constants_mod
15 : use carma_types_mod
16 : use carmastate_mod
17 : use carma_mod
18 :
19 : implicit none
20 :
21 : type(carma_type), intent(in) :: carma !! the carma object
22 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
23 : integer, intent(in) :: iz !! z index
24 : real(kind=f) :: scale_threshold !! Scaling factor for convergence thresholds
25 : integer, intent(inout) :: rc !! return code, negative indicates failure
26 :
27 : ! Local Variables
28 : integer :: ielem ! element index
29 : integer :: ibin ! bin index
30 : integer :: igas ! gas index
31 2837406876 : real(kind=f) :: previous_ice(NGAS) ! total ice at the start of substep
32 2837406876 : real(kind=f) :: previous_liquid(NGAS) ! total liquid at the start of substep
33 2837406876 : real(kind=f) :: previous_supsatl(NGAS) ! supersaturation wrt ice at the start of substep
34 2837406876 : real(kind=f) :: previous_supsati(NGAS) ! supersaturation wrt liquid at the start of substep
35 : real(kind=f) :: supsatold
36 : real(kind=f) :: supsatnew
37 : real(kind=f) :: srat
38 : real(kind=f) :: srat1
39 : real(kind=f) :: srat2
40 : real(kind=f) :: s_threshold
41 :
42 : 1 format(/,'microfast::ERROR - excessive change in supersaturation for ',a,' : iz=',i4,',xc=', &
43 : f7.2,',yc=',f7.2,',srat=',e10.3,',supsatiold=',e10.3,',supsatlold=',e10.3,',supsati=',e10.3, &
44 : ',supsatl=',e10.3,',t=',f6.2)
45 : 2 format('microfast::ERROR - conditions at beginning of the step : gc=',e10.3,',supsati=',e17.10, &
46 : ',supsatl=',e17.10,',t=',f6.2,',d_gc=',e10.3,',d_t=',f6.2)
47 : 3 format(/,'microfast::ERROR - excessive change in supersaturation for ',a,' : iz=',i4,',xc=', &
48 : f7.2,',yc=',f7.2,',supsatiold=',e10.3,',supsatlold=',e10.3,',supsati=',e10.3, &
49 : ',supsatl=',e10.3,',t=',f6.2)
50 :
51 : ! Set production and loss rates to zero.
52 1418703438 : call zeromicro(carma, cstate, iz, rc)
53 1418703438 : if (rc < RC_OK) return
54 :
55 : ! Calculate (implicit) particle loss rates for nucleation, growth,
56 : ! evaporation, melting, etc.
57 1418703438 : if (do_grow) then
58 :
59 : ! Save off the current condensate totals so the gas and latent heating can be
60 : ! figured out in a way that conserves mass and energy.
61 1418703438 : call totalcondensate(carma, cstate, iz, previous_ice, previous_liquid, rc)
62 1418703438 : if (rc < RC_OK) return
63 :
64 4256110314 : do igas = 1, NGAS
65 2837406876 : call supersat(carma, cstate, iz, igas, rc)
66 2837406876 : if (rc < RC_OK) return
67 :
68 2837406876 : previous_supsati(igas) = supsati(iz, igas)
69 4256110314 : previous_supsatl(igas) = supsatl(iz, igas)
70 : end do
71 :
72 : ! Have water vapor and sulfuric acid been defined?
73 1418703438 : if ((igash2o /= 0) .and. (igash2so4 /= 0)) then
74 :
75 : ! Are both gases avaialble?
76 1418703438 : if ((gc(iz, igash2o) > 0._f) .and. (gc(iz,igash2so4) > 0._f)) then
77 :
78 : ! See if any sulfates will form.
79 1418703437 : call sulfnuc(carma, cstate, iz, rc)
80 : endif
81 : end if
82 :
83 1418703438 : call growevapl(carma, cstate, iz, rc)
84 1418703438 : if (rc < RC_OK) return
85 :
86 1418703438 : call actdropl(carma, cstate, iz, rc)
87 1418703438 : if (rc < RC_OK) return
88 :
89 : ! The Koop, Tabazadeh and Mohler routines provide different schemes for aerosol freezing.
90 : ! Only one of these parameterizatons should be active at one time. However, any
91 : ! of these routines can be used in conjunction with heterogenous nucleation of glassy
92 : ! aerosols.
93 1418703438 : call freezaerl_tabazadeh2000(carma, cstate, iz, rc)
94 1418703438 : if (rc < RC_OK) return
95 :
96 1418703438 : call freezaerl_koop2000(carma, cstate, iz, rc)
97 1418703438 : if (rc < RC_OK) return
98 :
99 1418703438 : call freezaerl_mohler2010(carma, cstate, iz, rc)
100 1418703438 : if (rc < RC_OK) return
101 :
102 1418703438 : call freezglaerl_murray2010(carma, cstate, iz, rc)
103 1418703438 : if (rc < RC_OK) return
104 :
105 1418703438 : call hetnucl(carma, cstate, iz, rc)
106 1418703438 : if (rc < RC_OK) return
107 :
108 1418703438 : call freezdropl(carma, cstate, iz, rc)
109 1418703438 : if (rc < RC_OK) return
110 :
111 1418703438 : call melticel(carma, cstate, iz, rc)
112 1418703438 : if (rc < RC_OK) return
113 : endif
114 :
115 1418703438 : if(do_coremasscheck)then
116 0 : call coremasscheck( carma, cstate, iz, .false.,.true.,.true., "Beforegrowth", rc )
117 0 : if (rc < RC_OK) return
118 : end if
119 :
120 : ! Calculate particle production terms and solve for particle
121 : ! concentrations at end of time step.
122 11349627504 : do ielem = 1,NELEM
123 >20996*10^7 : do ibin = 1,NBIN
124 :
125 >19861*10^7 : if( do_grow )then
126 >19861*10^7 : call growp(carma, cstate, iz, ibin, ielem, rc)
127 >19861*10^7 : if (rc < RC_OK) return
128 :
129 >19861*10^7 : call upgxfer(carma, cstate, iz, ibin, ielem, rc)
130 >19861*10^7 : if (rc < RC_OK) return
131 : endif
132 :
133 >19861*10^7 : call psolve(carma, cstate, iz, ibin, ielem, rc)
134 >20854*10^7 : if (rc < RC_OK) return
135 : enddo
136 : enddo
137 :
138 : ! Calculate particle production terms for evaporation;
139 : ! gas loss rates and production terms due to particle nucleation;
140 : ! growth, and evaporation;
141 : ! apply evaporation production terms to particle concentrations;
142 : ! and solve for gas concentrations at end of time step.
143 1418703438 : if (do_grow) then
144 1418703438 : call evapp(carma, cstate, iz, rc)
145 1418703438 : if (rc < RC_OK) return
146 :
147 1418703438 : call downgxfer(carma, cstate, iz, rc)
148 1418703438 : if (rc < RC_OK) return
149 :
150 : ! NOTE: Not needed because changes in gas concentrations and latent
151 : ! heats are now calculated later in gsolve using total condensate.
152 : ! call gasexchange(carma, cstate, iz, rc)
153 : ! if (rc < RC_OK) return
154 :
155 1418703438 : call downgevapply(carma, cstate, iz, rc)
156 1418703438 : if (rc < RC_OK) return
157 :
158 1418703438 : call gsolve(carma, cstate, iz, previous_ice, previous_liquid, scale_threshold, rc)
159 1418703438 : if (rc /=RC_OK) return
160 : endif
161 :
162 1332108521 : call coremasscheck( carma, cstate, iz, .true.,.false.,.false., "AfterGsolve", rc )
163 1332108521 : if (rc < RC_OK) return
164 :
165 : ! Update temperature if thermal processes requested
166 1332108521 : if (do_thermo) then
167 0 : call tsolve(carma, cstate, iz, scale_threshold, rc)
168 0 : if (rc /= RC_OK) return
169 : endif
170 :
171 : ! Update saturation ratios
172 1332108521 : if (do_grow .or. do_thermo) then
173 3996325563 : do igas = 1, NGAS
174 2664217042 : call supersat(carma, cstate, iz, igas, rc)
175 2664217042 : if (rc < RC_OK) return
176 :
177 : ! Check to see how much the supersaturation changed during this step. If it
178 : ! has changed to much, then cause a retry.
179 2664217042 : if (t(iz) >= T0) then
180 1152128432 : supsatold = previous_supsatl(igas)
181 1152128432 : supsatnew = supsatl(iz,igas)
182 : else
183 1512088610 : supsatold = previous_supsati(igas)
184 1512088610 : supsatnew = supsati(iz,igas)
185 : end if
186 :
187 : ! If ds_threshold is positive, then it indicates that the criteria should
188 : ! be based on the percentage change in saturation.
189 3996325563 : if (ds_threshold(igas) > 0._f) then
190 :
191 0 : if (supsatold >= 1.e-4_f) then
192 0 : srat1 = abs(supsatnew / supsatold - 1._f)
193 : else
194 : srat1 = 0._f
195 : end if
196 :
197 0 : if (supsatnew >= 1.e-4_f) then
198 0 : srat2 = abs(supsatold / supsatnew - 1._f)
199 : else
200 : srat2 = 0._f
201 : end if
202 :
203 0 : srat = max(srat1, srat2)
204 :
205 : ! Don't let one substep change the supersaturation by too much.
206 : if (ds_threshold(igas) > 0._f) then
207 : ! if (srat >= ds_threshold(igas)) then
208 0 : if ((srat >= ds_threshold(igas)) .and. (abs(supsatold - supsatnew) > 0.1_f)) then
209 0 : if (do_substep) then
210 0 : if (nretries == maxretries) then
211 0 : if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, &
212 0 : xc, yc, srat, previous_supsati(igas), previous_supsatl(igas), &
213 0 : supsati(iz, igas), supsatl(iz,igas), t(iz)
214 0 : if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz, igas), &
215 0 : supsatlold(iz,igas), told(iz), d_gc(iz, igas), d_t(iz)
216 : end if
217 :
218 0 : rc = RC_WARNING_RETRY
219 : else
220 0 : if (do_print) write(LUNOPRT,1) trim(gasname(igas)), iz, xc, yc, &
221 0 : srat, previous_supsati(igas), previous_supsatl(igas), &
222 0 : supsati(iz, igas), supsatl(iz,igas), t(iz)
223 : end if
224 : end if
225 : end if
226 :
227 :
228 : ! If ds_threshold is negative, then it indicates that the criteria is based
229 : ! upon the supersaturation crossing 0, Indicating a shift from growth to
230 : ! evaporation and a potential overshoot in the result.
231 2664217042 : else if (ds_threshold(igas) < 0._f) then
232 :
233 : ! Adjust the saturation threshold to allow a worse solution if getting a better
234 : ! solution is taking too much time. The particular solution at any individual
235 : ! point is probably not going to affect the overall result by too much.
236 1332108521 : s_threshold = abs(ds_threshold(igas))
237 :
238 1332108521 : if (nretries >= (0.8_f * maxretries)) then
239 5636096 : s_threshold = 4._f * s_threshold
240 1326472425 : else if (nretries >= (0.7_f * maxretries)) then
241 11681906 : s_threshold = 3.5_f * s_threshold
242 1314790519 : else if (nretries >= (0.6_f * maxretries)) then
243 42252495 : s_threshold = 3._f * s_threshold
244 1272538024 : else if (nretries >= (0.5_f * maxretries)) then
245 184273692 : s_threshold = 2.5_f * s_threshold
246 1088264332 : else if (nretries >= (0.4_f * maxretries)) then
247 345062953 : s_threshold = 2._f * s_threshold
248 : end if
249 :
250 : ! If the supersaturation changed signs, then we went from growth to evaporation
251 : ! or vice versa. Don't let the new supersaturation go too far past 0 in one substep.
252 : ! This is to prevent overshooting as growth/evaporation should normally stop when
253 : ! the supersaturation is 0.
254 1332108521 : if (((supsatnew * supsatold) < 0._f) .and. (abs(supsatnew) > s_threshold)) then
255 :
256 0 : if (do_substep) then
257 0 : if (nretries == maxretries) then
258 0 : if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, &
259 0 : xc, yc, previous_supsati(igas), previous_supsatl(igas), &
260 0 : supsati(iz, igas), supsatl(iz,igas), t(iz)
261 0 : if (do_print) write(LUNOPRT,2) gcl(iz,igas), supsatiold(iz, igas), &
262 0 : supsatlold(iz,igas), told(iz), d_gc(iz, igas), d_t(iz)
263 : end if
264 : else
265 0 : if (do_print) write(LUNOPRT,3) trim(gasname(igas)), iz, &
266 0 : xc, yc, previous_supsati(igas), previous_supsatl(igas), &
267 0 : supsati(iz, igas), supsatl(iz,igas), t(iz)
268 : end if
269 :
270 0 : rc = RC_WARNING_RETRY
271 : end if
272 : end if
273 : end do
274 : endif
275 :
276 : ! Return to caller with new particle and gas concentrations.
277 : return
278 1418703438 : end
|
