Line data Source code
1 : module mo_imp_sol
2 : use shr_kind_mod, only : r8 => shr_kind_r8
3 : use chem_mods, only : clscnt4, gas_pcnst, clsmap
4 : use cam_logfile, only : iulog
5 : implicit none
6 : private
7 : public :: imp_slv_inti, imp_sol
8 : save
9 : real(r8), parameter :: rel_err = 1.e-3_r8
10 : real(r8), parameter :: high_rel_err = 1.e-4_r8
11 : !-----------------------------------------------------------------------
12 : ! Newton-Raphson iteration limits
13 : !-----------------------------------------------------------------------
14 : integer, parameter :: itermax = 11
15 : integer, parameter :: cut_limit = 5
16 : real(r8), parameter :: small = 1.e-40_r8
17 : real(r8) :: epsilon(clscnt4)
18 : logical :: factor(itermax)
19 : contains
20 0 : subroutine imp_slv_inti
21 : !-----------------------------------------------------------------------
22 : ! ... Initialize the implict solver
23 : !-----------------------------------------------------------------------
24 : use mo_chem_utls, only : get_spc_ndx
25 : implicit none
26 : !-----------------------------------------------------------------------
27 : ! ... Local variables
28 : !-----------------------------------------------------------------------
29 : integer :: m, ox_ndx, o3a_ndx
30 : real(r8) :: eps(gas_pcnst)
31 0 : factor(:) = .true.
32 : eps(:) = rel_err
33 0 : ox_ndx = get_spc_ndx( 'OX' )
34 0 : if( ox_ndx < 1 ) then
35 0 : ox_ndx = get_spc_ndx( 'O3' )
36 : end if
37 0 : if( ox_ndx > 0 ) then
38 0 : eps(ox_ndx) = high_rel_err
39 : end if
40 0 : m = get_spc_ndx( 'NO' )
41 0 : if( m > 0 ) then
42 0 : eps(m) = high_rel_err
43 : end if
44 0 : m = get_spc_ndx( 'NO2' )
45 0 : if( m > 0 ) then
46 0 : eps(m) = high_rel_err
47 : end if
48 0 : m = get_spc_ndx( 'NO3' )
49 0 : if( m > 0 ) then
50 0 : eps(m) = high_rel_err
51 : end if
52 0 : m = get_spc_ndx( 'HNO3' )
53 0 : if( m > 0 ) then
54 0 : eps(m) = high_rel_err
55 : end if
56 0 : m = get_spc_ndx( 'HO2NO2' )
57 0 : if( m > 0 ) then
58 0 : eps(m) = high_rel_err
59 : end if
60 0 : m = get_spc_ndx( 'N2O5' )
61 0 : if( m > 0 ) then
62 0 : eps(m) = high_rel_err
63 : end if
64 0 : m = get_spc_ndx( 'OH' )
65 0 : if( m > 0 ) then
66 0 : eps(m) = high_rel_err
67 : end if
68 0 : m = get_spc_ndx( 'HO2' )
69 0 : if( m > 0 ) then
70 0 : eps(m) = high_rel_err
71 : end if
72 0 : o3a_ndx = get_spc_ndx( 'O3A' )
73 0 : if( o3a_ndx > 0 ) then
74 0 : eps(m) = high_rel_err
75 : end if
76 0 : m = get_spc_ndx( 'XNO' )
77 0 : if( m > 0 ) then
78 0 : eps(m) = high_rel_err
79 : end if
80 0 : m = get_spc_ndx( 'XNO2' )
81 0 : if( m > 0 ) then
82 0 : eps(m) = high_rel_err
83 : end if
84 0 : m = get_spc_ndx( 'XNO3' )
85 0 : if( m > 0 ) then
86 0 : eps(m) = high_rel_err
87 : end if
88 0 : m = get_spc_ndx( 'XHNO3' )
89 0 : if( m > 0 ) then
90 0 : eps(m) = high_rel_err
91 : end if
92 0 : m = get_spc_ndx( 'XHO2NO2' )
93 0 : if( m > 0 ) then
94 0 : eps(m) = high_rel_err
95 : end if
96 0 : m = get_spc_ndx( 'XNO2NO3' )
97 0 : if( m > 0 ) then
98 0 : eps(m) = high_rel_err
99 : end if
100 0 : m = get_spc_ndx( 'NO2XNO3' )
101 0 : if( m > 0 ) then
102 0 : eps(m) = high_rel_err
103 : end if
104 : do m = 1,clscnt4
105 : epsilon(m) = eps(clsmap(m,4))
106 : end do
107 0 : end subroutine imp_slv_inti
108 0 : subroutine imp_sol( base_sol, reaction_rates, het_rates, extfrc, delt, &
109 0 : ncol,nlev, lchnk, prod_out, loss_out )
110 : !-----------------------------------------------------------------------
111 : ! ... imp_sol advances the volumetric mixing ratio
112 : ! forward one time step via the fully implicit euler scheme.
113 : ! this source is meant for small l1 cache machines such as
114 : ! the intel pentium and itanium cpus
115 : !-----------------------------------------------------------------------
116 : use chem_mods, only : rxntot, extcnt, nzcnt, permute, cls_rxt_cnt
117 : use mo_tracname, only : solsym
118 : use mo_lin_matrix, only : linmat
119 : use mo_nln_matrix, only : nlnmat
120 : use mo_lu_factor, only : lu_fac
121 : use mo_lu_solve, only : lu_slv
122 : use mo_prod_loss, only : imp_prod_loss
123 : use mo_indprd, only : indprd
124 : use time_manager, only : get_nstep
125 : use perf_mod, only : t_startf, t_stopf
126 : implicit none
127 : !-----------------------------------------------------------------------
128 : ! ... dummy args
129 : !-----------------------------------------------------------------------
130 : integer, intent(in) :: ncol ! columns in chunck
131 : integer, intent(in) :: nlev
132 : integer, intent(in) :: lchnk ! chunk id
133 : real(r8), intent(in) :: delt ! time step (s)
134 : real(r8), intent(in) :: reaction_rates(ncol,nlev,max(1,rxntot)) ! rxt rates (1/cm^3/s)
135 : real(r8), intent(in) :: extfrc(ncol,nlev,max(1,extcnt)) ! external in-situ forcing (1/cm^3/s)
136 : real(r8), intent(in) :: het_rates(ncol,nlev,max(1,gas_pcnst)) ! washout rates (1/s)
137 : real(r8), intent(inout) :: base_sol(ncol,nlev,gas_pcnst) ! species mixing ratios (vmr)
138 : real(r8), intent(out) :: prod_out(ncol,nlev,max(1,clscnt4))
139 : real(r8), intent(out) :: loss_out(ncol,nlev,max(1,clscnt4))
140 : !-----------------------------------------------------------------------
141 : ! ... local variables
142 : !-----------------------------------------------------------------------
143 : integer :: nr_iter, &
144 : lev, &
145 : i, &
146 : j, &
147 : k, l, &
148 : m
149 : integer :: fail_cnt, cut_cnt, stp_con_cnt
150 : integer :: nstep
151 : real(r8) :: interval_done, dt, dti
152 : real(r8) :: max_delta(max(1,clscnt4))
153 : real(r8) :: sys_jac(max(1,nzcnt))
154 : real(r8) :: lin_jac(max(1,nzcnt))
155 : real(r8), dimension(max(1,clscnt4)) :: &
156 : solution, &
157 : forcing, &
158 : iter_invariant, &
159 : prod, &
160 : loss
161 : real(r8) :: lrxt(max(1,rxntot))
162 : real(r8) :: lsol(max(1,gas_pcnst))
163 : real(r8) :: lhet(max(1,gas_pcnst))
164 : real(r8), dimension(ncol,nlev,max(1,clscnt4)) :: &
165 0 : ind_prd
166 : logical :: convergence
167 : logical :: frc_mask, iter_conv
168 : logical :: converged(max(1,clscnt4))
169 0 : solution(:) = 0._r8
170 : !-----------------------------------------------------------------------
171 : ! ... class independent forcing
172 : !-----------------------------------------------------------------------
173 0 : if( cls_rxt_cnt(1,4) > 0 .or. extcnt > 0 ) then
174 : call indprd( 4, ind_prd, clscnt4, base_sol, extfrc, &
175 0 : reaction_rates, ncol )
176 : else
177 0 : do m = 1,max(1,clscnt4)
178 0 : ind_prd(:,:,m) = 0._r8
179 : end do
180 : end if
181 0 : level_loop : do lev = 1,nlev
182 0 : column_loop : do i = 1,ncol
183 : !-----------------------------------------------------------------------
184 : ! ... transfer from base to local work arrays
185 : !-----------------------------------------------------------------------
186 : do m = 1,rxntot
187 : lrxt(m) = reaction_rates(i,lev,m)
188 : end do
189 : if( gas_pcnst > 0 ) then
190 : do m = 1,gas_pcnst
191 : lhet(m) = het_rates(i,lev,m)
192 : end do
193 : end if
194 : !-----------------------------------------------------------------------
195 : ! ... time step loop
196 : !-----------------------------------------------------------------------
197 0 : dt = delt
198 0 : cut_cnt = 0
199 0 : fail_cnt = 0
200 0 : stp_con_cnt = 0
201 0 : interval_done = 0._r8
202 : time_step_loop : do
203 0 : dti = 1._r8 / dt
204 : !-----------------------------------------------------------------------
205 : ! ... transfer from base to local work arrays
206 : !-----------------------------------------------------------------------
207 : do m = 1,gas_pcnst
208 : lsol(m) = base_sol(i,lev,m)
209 : end do
210 : !-----------------------------------------------------------------------
211 : ! ... transfer from base to class array
212 : !-----------------------------------------------------------------------
213 : do k = 1,clscnt4
214 : j = clsmap(k,4)
215 : m = permute(k,4)
216 : solution(m) = lsol(j)
217 : end do
218 : !-----------------------------------------------------------------------
219 : ! ... set the iteration invariant part of the function f(y)
220 : !-----------------------------------------------------------------------
221 0 : if( cls_rxt_cnt(1,4) > 0 .or. extcnt > 0 ) then
222 : do m = 1,clscnt4
223 : iter_invariant(m) = dti * solution(m) + ind_prd(i,lev,m)
224 : end do
225 : else
226 : do m = 1,clscnt4
227 : iter_invariant(m) = dti * solution(m)
228 : end do
229 : end if
230 : !-----------------------------------------------------------------------
231 : ! ... the linear component
232 : !-----------------------------------------------------------------------
233 0 : if( cls_rxt_cnt(2,4) > 0 ) then
234 0 : call t_startf( 'lin_mat' )
235 0 : call linmat( lin_jac, lsol, lrxt, lhet )
236 0 : call t_stopf( 'lin_mat' )
237 : end if
238 : !=======================================================================
239 : ! the newton-raphson iteration for f(y) = 0
240 : !=======================================================================
241 0 : iter_loop : do nr_iter = 1,itermax
242 : !-----------------------------------------------------------------------
243 : ! ... the non-linear component
244 : !-----------------------------------------------------------------------
245 0 : if( factor(nr_iter) ) then
246 0 : call t_startf( 'nln_mat' )
247 0 : call nlnmat( sys_jac, lsol, lrxt, lin_jac, dti )
248 0 : call t_stopf( 'nln_mat' )
249 : !-----------------------------------------------------------------------
250 : ! ... factor the "system" matrix
251 : !-----------------------------------------------------------------------
252 0 : call t_startf( 'lu_fac' )
253 0 : call lu_fac( sys_jac )
254 0 : call t_stopf( 'lu_fac' )
255 : end if
256 : !-----------------------------------------------------------------------
257 : ! ... form f(y)
258 : !-----------------------------------------------------------------------
259 0 : call t_startf( 'prod_loss' )
260 0 : call imp_prod_loss( prod, loss, lsol, lrxt, lhet )
261 0 : call t_stopf( 'prod_loss' )
262 : do m = 1,clscnt4
263 : forcing(m) = solution(m)*dti - (iter_invariant(m) + prod(m) - loss(m))
264 : end do
265 : !-----------------------------------------------------------------------
266 : ! ... solve for the mixing ratio at t(n+1)
267 : !-----------------------------------------------------------------------
268 0 : call t_startf( 'lu_slv' )
269 0 : call lu_slv( sys_jac, forcing )
270 0 : call t_stopf( 'lu_slv' )
271 : do m = 1,clscnt4
272 : solution(m) = solution(m) + forcing(m)
273 : end do
274 : !-----------------------------------------------------------------------
275 : ! ... convergence measures
276 : !-----------------------------------------------------------------------
277 : if( nr_iter > 1 ) then
278 : do k = 1,clscnt4
279 : m = permute(k,4)
280 : if( abs(solution(m)) > 1.e-20_r8 ) then
281 : max_delta(k) = abs( forcing(m)/solution(m) )
282 : else
283 : max_delta(k) = 0._r8
284 : end if
285 : end do
286 : end if
287 : !-----------------------------------------------------------------------
288 : ! ... limit iterate
289 : !-----------------------------------------------------------------------
290 0 : where( solution(:) < 0._r8 )
291 : solution(:) = 0._r8
292 : endwhere
293 : !-----------------------------------------------------------------------
294 : ! ... transfer latest solution back to work array
295 : !-----------------------------------------------------------------------
296 : do k = 1,clscnt4
297 : j = clsmap(k,4)
298 : m = permute(k,4)
299 : lsol(j) = solution(m)
300 : end do
301 : !-----------------------------------------------------------------------
302 : ! ... check for convergence
303 : !-----------------------------------------------------------------------
304 0 : converged(:) = .true.
305 0 : if( nr_iter > 1 ) then
306 : do k = 1,clscnt4
307 : m = permute(k,4)
308 : frc_mask = abs( forcing(m) ) > small
309 : if( frc_mask ) then
310 : converged(k) = abs(forcing(m)) <= epsilon(k)*abs(solution(m))
311 : else
312 : converged(k) = .true.
313 : end if
314 : end do
315 0 : convergence = all( converged(:) )
316 0 : if( convergence ) then
317 : exit
318 : end if
319 : end if
320 : end do iter_loop
321 : !-----------------------------------------------------------------------
322 : ! ... check for newton-raphson convergence
323 : !-----------------------------------------------------------------------
324 0 : if( .not. convergence ) then
325 : !-----------------------------------------------------------------------
326 : ! ... non-convergence
327 : !-----------------------------------------------------------------------
328 0 : fail_cnt = fail_cnt + 1
329 0 : nstep = get_nstep()
330 : write(iulog,'('' imp_sol: Time step '',1p,e21.13,'' failed to converge @ (lchnk,lev,col,nstep) = '',4i6)') &
331 0 : dt,lchnk,lev,i,nstep
332 0 : stp_con_cnt = 0
333 0 : if( cut_cnt < cut_limit ) then
334 0 : cut_cnt = cut_cnt + 1
335 0 : if( cut_cnt < cut_limit ) then
336 0 : dt = .5_r8 * dt
337 : else
338 0 : dt = .1_r8 * dt
339 : end if
340 : cycle time_step_loop
341 : else
342 : write(iulog,'('' imp_sol: Failed to converge @ (lchnk,lev,col,nstep,dt,time) = '',4i6,1p,2e21.13)') &
343 0 : lchnk,lev,i,nstep,dt,interval_done+dt
344 : do m = 1,clscnt4
345 0 : if( .not. converged(m) ) then
346 : write(iulog,'(1x,a8,1x,1pe10.3)') solsym(clsmap(m,4)), max_delta(m)
347 : end if
348 : end do
349 : end if
350 : end if
351 : !-----------------------------------------------------------------------
352 : ! ... check for interval done
353 : !-----------------------------------------------------------------------
354 0 : interval_done = interval_done + dt
355 0 : if( abs( delt - interval_done ) <= .0001_r8 ) then
356 0 : if( fail_cnt > 0 ) then
357 0 : write(iulog,*) 'imp_sol : @ (lchnk,lev,col) = ',lchnk,lev,i,' failed ',fail_cnt,' times'
358 : end if
359 : exit time_step_loop
360 : else
361 : !-----------------------------------------------------------------------
362 : ! ... transfer latest solution back to base array
363 : !-----------------------------------------------------------------------
364 0 : if( convergence ) then
365 0 : stp_con_cnt = stp_con_cnt + 1
366 : end if
367 : do m = 1,gas_pcnst
368 : base_sol(i,lev,m) = lsol(m)
369 : end do
370 0 : if( stp_con_cnt >= 2 ) then
371 0 : dt = 2._r8*dt
372 0 : stp_con_cnt = 0
373 : end if
374 0 : dt = min( dt,delt-interval_done )
375 : ! write(iulog,'('' imp_sol: New time step '',1p,e21.13)') dt
376 : end if
377 : end do time_step_loop
378 : !-----------------------------------------------------------------------
379 : ! ... Transfer latest solution back to base array
380 : !-----------------------------------------------------------------------
381 0 : cls_loop: do k = 1,clscnt4
382 : j = clsmap(k,4)
383 : m = permute(k,4)
384 : base_sol(i,lev,j) = solution(m)
385 : ! output diagnostics
386 : prod_out(i,lev,k) = prod(k) + ind_prd(i,lev,k)
387 : loss_out(i,lev,k) = loss(k)
388 : end do cls_loop
389 : end do column_loop
390 : end do level_loop
391 0 : end subroutine imp_sol
392 : end module mo_imp_sol
|
