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 manages the calculations that update state variables
6 : !! of the model with new values at the current simulation time. It supports
7 : !! a retry mechanism, so the the number of steps can be increased dynamically
8 : !! if the fast microphysics was not able to generate a valid solution. The
9 : !! validity of the solution is control by the convergence thresholds
10 : !! (dgc_threshold, dt_threshold and ds_threshold)
11 : !!
12 : !! NOTE: For cloud models, this routine may get called multiple times, once for
13 : !! in-cloud calculations and again for clear sky.
14 : !!
15 : !! @author Bardeen
16 : !! @version Jan 2012
17 1050624 : subroutine newstate_calc(carma, cstate, scale_threshold, rc)
18 :
19 : ! types
20 : use carma_precision_mod
21 : use carma_enums_mod
22 : use carma_constants_mod
23 : use carma_types_mod
24 : use carmastate_mod
25 : use carma_mod
26 :
27 : implicit none
28 :
29 : type(carma_type), intent(in) :: carma !! the carma object
30 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
31 : real(kind=f), intent(in) :: scale_threshold(NZ) !! Scaling factor for convergence thresholds
32 : integer, intent(inout) :: rc !! return code, negative indicates failure
33 :
34 2101248 : real(kind=f) :: sedlayer(NBIN,NELEM)
35 2101248 : real(kind=f) :: pcd_last(NBIN,NELEM)
36 : integer :: kb
37 : integer :: ke
38 : integer :: idk
39 : integer :: iz
40 : integer :: isubstep
41 : integer :: igroup
42 : integer :: igas
43 : integer :: ielem
44 : integer :: ibin
45 : integer :: ntsubsteps
46 : logical :: takeSteps
47 : real(kind=f) :: fraction ! Fraction of dT, dgc and pdc to be added in a substep.
48 :
49 : 1 format(/,'newstate::ERROR - Substep failed, maximum retries execeed. : iz=',i4,',isubstep=',i12, &
50 : ',ntsubsteps=',i12,',nretries=',F9.0)
51 :
52 :
53 : ! Redetermine the maximum particle values.
54 1050624 : if ((do_vtran) .or. do_incloud) then
55 74594304 : do iz = 1, NZ
56 73543680 : call maxconc(carma, cstate, iz, rc)
57 74594304 : if (rc < RC_OK) return
58 : end do
59 : end if
60 :
61 : ! Calculate changes in particle concentrations due to microphysical
62 : ! processes, part 1. (potentially slower microphysical calcs)
63 : ! All spatial points are handled by one call to this routine.
64 1050624 : if (do_coag) then
65 1050624 : call microslow(carma, cstate, rc)
66 1050624 : if (rc < RC_OK) return
67 : endif
68 :
69 1050624 : call fixcorecol(carma, cstate, rc)
70 :
71 : ! If there is any microsphysics that happens on a faster time scale,
72 : ! then check to see if the time step needs to be subdivided and then
73 : ! perform the fast microphysical calculations.
74 1050624 : if (do_grow) then
75 :
76 : ! Set vertical loop index to increment downwards
77 : ! (for substepping of sedimentation)
78 1050624 : if (igridv .eq. I_CART) then
79 0 : kb = NZ
80 0 : ke = 1
81 0 : idk = -1
82 : else
83 1050624 : kb = 1
84 1050624 : ke = NZ
85 1050624 : idk = 1
86 : endif
87 :
88 : ! Initialize sedimentation source to zero at top of model
89 155492352 : dpc_sed(:,:) = 0._f
90 :
91 : ! Save the results from the slow operations, since we might need to retry the
92 : ! fast operations
93 10451607552 : pcl(:,:,:) = pc(:,:,:)
94 :
95 1050624 : if (do_substep) then
96 3151872 : do igas = 1,NGAS
97 150239232 : gcl(:,igas) = gc(:,igas)
98 : end do
99 74594304 : told(:) = t(:)
100 : endif
101 :
102 :
103 74594304 : do iz = kb,ke,idk
104 :
105 : ! Compute or specify number of sub-timestep intervals for current spatial point
106 : ! (Could be same for all spatial pts, or could vary as a function of location)
107 73543680 : ntsubsteps = minsubsteps
108 :
109 73543680 : call nsubsteps(carma, cstate, iz, dtime_orig, ntsubsteps, rc)
110 73543680 : if (rc < RC_OK) return
111 :
112 : ! Grab sedimentation source for entire step for this layer
113 : ! and set accumlated source for underlying layer to zero
114 10884464640 : sedlayer(:,:) = dpc_sed(:,:)
115 :
116 : ! Do sub-timestepping for current spatial grid point, and allow for
117 : ! retrying should this level of substepping not be enough to keep the
118 : ! gas concentration from going negative.
119 73543680 : nretries = 0._f
120 73543680 : takeSteps = .true.
121 :
122 233682277 : do while (takeSteps)
123 :
124 : ! Compute sub-timestep time interval for current spatial grid point
125 160138597 : dtime = dtime_orig / ntsubsteps
126 :
127 : ! Don't retry unless requested.
128 160138597 : takeSteps = .false.
129 :
130 : ! Reset the amount that has been collected to sedimented down to the
131 : ! layer below.
132 23700512356 : dpc_sed(:,:) = 0._f
133 :
134 : ! Reset the total nucleation for the step.
135 23700512356 : pc_nucl(iz,:,:) = 0._f
136 :
137 : ! Remember the amount of detrained particles.
138 160138597 : if (do_detrain) then
139 0 : pcd_last(:,:) = pcd(iz,:,:)
140 : end if
141 :
142 : ! Reset average heating rates.
143 160138597 : rlheat(iz) = 0._f
144 160138597 : partheat(iz) = 0._f
145 :
146 1565790798 : do isubstep = 1,ntsubsteps
147 :
148 : ! If substepping, then increment the gas concentration and the temperature by
149 : ! an amount for one substep.
150 1418703438 : if (do_substep) then
151 :
152 : ! Since we don't really know how the gas and temperature changes arrived during the
153 : ! step, we can try different assumptions for how the gas and temperature are add to
154 : ! the values from the previous substep.
155 :
156 : ! Linear increment for substepping.
157 1418703438 : fraction = 1._f / ntsubsteps
158 :
159 4256110314 : do igas = 1,NGAS
160 4256110314 : gc(iz,igas) = gc(iz,igas) + d_gc(iz,igas) * fraction
161 : enddo
162 :
163 1418703438 : t(iz) = t(iz) + d_t(iz) * fraction
164 :
165 :
166 : ! Detrainment puts the full gridbox amount into the incloud portion.
167 1418703438 : if (do_detrain) then
168 0 : pc(iz,:,:) = pc(iz,:,:) + pcd_last(:,:) * fraction
169 0 : pcd(iz,:,:) = pcd(iz,:,:) - pcd_last(:,:) * fraction
170 : end if
171 : endif
172 :
173 :
174 : ! Redetermine maximum particle concentrations.
175 1418703438 : call maxconc(carma, cstate, iz, rc)
176 1418703438 : if (rc < RC_OK) return
177 :
178 1418703438 : call coremasscheck( carma, cstate, iz, .true.,.false.,.false., "BeforeMicrofast", rc )
179 1418703438 : if (rc < RC_OK) return
180 :
181 : ! Calculate changes in particle concentrations for current spatial point
182 : ! due to microphysical processes, part 2. (faster microphysical calcs)
183 1418703438 : call microfast(carma, cstate, iz, scale_threshold(iz), rc)
184 1418703438 : if (rc < RC_OK) return
185 :
186 1418703438 : call coremasscheck( carma, cstate, iz, .true.,.false.,.false., "AfterMicrofast", rc )
187 1418703438 : if (rc < RC_OK) return
188 :
189 : ! If there was a retry warning message and substepping is enabled, then retry
190 : ! the operation with more substepping.
191 2910950556 : if (rc == RC_WARNING_RETRY) then
192 86594917 : if (do_substep) then
193 :
194 : ! Only retry for so long ...
195 86594917 : nretries = nretries + 1
196 :
197 86594917 : if (nretries > maxretries) then
198 0 : if (do_print) write(LUNOPRT,1) iz, isubstep, ntsubsteps, nretries - 1._f
199 0 : rc = RC_ERROR
200 : exit
201 : end if
202 :
203 : ! Try twice the substeps
204 : !
205 : ! NOTE: We are going to rely upon retries, so don't clutter the log
206 : ! with retry print statements. They slow down the run.
207 86594917 : ntsubsteps = ntsubsteps * 2
208 :
209 : ! if (do_print) write(LUNOPRT,*) "newstate::WARNING - Substep failed, retrying with ", ntsubsteps, " substeps."
210 :
211 : ! Reset the state to the beginning of the step
212 12816047716 : pc(iz,:,:) = pcl(iz,:,:)
213 12816047716 : pcd(iz,:,:) = pcd_last(:,:)
214 86594917 : t(iz) = told(iz)
215 259784751 : do igas = 1,NGAS
216 173189834 : gc(iz,igas) = gcl(iz,igas)
217 :
218 : ! Now that we have reset the gas concentration, we need to recalculate the supersaturation.
219 173189834 : call supersat(carma, cstate, iz, igas, rc)
220 259784751 : if (rc < RC_OK) return
221 : end do
222 :
223 86594917 : rc = RC_OK
224 86594917 : takeSteps = .true.
225 86594917 : exit
226 :
227 :
228 : ! If substepping is not enabled, than the retry warning should be treated as an error.
229 : else
230 :
231 0 : if (do_print) write(LUNOPRT,*) "newstate::ERROR - Step failed, suggest enabling substepping."
232 0 : rc = RC_ERROR
233 : exit
234 : end if
235 : end if
236 : end do
237 : end do
238 :
239 :
240 : ! Keep track of substepping and retry statistics for performance tuning.
241 73543680 : max_nsubstep = max(max_nsubstep, ntsubsteps)
242 73543680 : max_nretry = max(max_nretry, nretries)
243 :
244 73543680 : nstep = nstep + 1._f
245 73543680 : nsubstep = nsubstep + ntsubsteps
246 73543680 : nretry = nretry + nretries
247 :
248 74594304 : if (do_substep) zsubsteps(iz) = ntsubsteps
249 : end do
250 :
251 : ! Restore normal timestep
252 1050624 : dtime = dtime_orig
253 :
254 : else
255 :
256 : ! If there is no reason to substep, but substepping was enabled, get the gas and
257 : ! temperature back to their final states.
258 0 : if (do_substep) then
259 :
260 0 : do igas = 1,NGAS
261 0 : gc(:,igas) = gc(:,igas) + d_gc(:,igas)
262 : enddo
263 :
264 0 : t(:) = t(:) + d_t(:)
265 : end if
266 :
267 : ! Do the detrainment, if it was being done in the growth loop.
268 0 : if (do_detrain) then
269 0 : pc(:,:,:) = pc(:,:,:) + pcd(:,:,:)
270 :
271 : ! Remove the ice from the detrained ice, so that total ice will be conserved.
272 0 : pcd(:,:,:) = 0._f
273 : end if
274 : end if
275 :
276 : ! Calculate average heating rates.
277 1050624 : if (do_grow) then
278 74594304 : rlheat(:) = rlheat(:) / dtime
279 74594304 : partheat(:) = partheat(:) / dtime
280 : end if
281 :
282 : ! Return to caller with new state computed
283 : return
284 1050624 : end
|
