Line data Source code
1 : !=======================================================================
2 : !
3 : ! Viscous-plastic sea ice dynamics model
4 : ! Computes ice velocity and deformation
5 : !
6 : ! See:
7 : !
8 : ! Lemieux, J.‐F., Tremblay, B., Thomas, S., Sedláček, J., and Mysak, L. A. (2008),
9 : ! Using the preconditioned Generalized Minimum RESidual (GMRES) method to solve
10 : ! the sea‐ice momentum equation, J. Geophys. Res., 113, C10004, doi:10.1029/2007JC004680.
11 : !
12 : ! Hibler, W. D., and Ackley, S. F. (1983), Numerical simulation of the Weddell Sea pack ice,
13 : ! J. Geophys. Res., 88( C5), 2873– 2887, doi:10.1029/JC088iC05p02873.
14 : !
15 : ! Y. Saad. A Flexible Inner-Outer Preconditioned GMRES Algorithm. SIAM J. Sci. Comput.,
16 : ! 14(2):461–469, 1993. URL: https://doi.org/10.1137/0914028, doi:10.1137/0914028.
17 : !
18 : ! C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995.
19 : ! (https://www.siam.org/books/textbooks/fr16_book.pdf)
20 : !
21 : ! Y. Saad, Iterative Methods for Sparse Linear Systems. SIAM, 2003.
22 : ! (http://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf)
23 : !
24 : ! Walker, H. F., & Ni, P. (2011). Anderson Acceleration for Fixed-Point Iterations.
25 : ! SIAM Journal on Numerical Analysis, 49(4), 1715–1735. https://doi.org/10.1137/10078356X
26 : !
27 : ! Fang, H., & Saad, Y. (2009). Two classes of multisecant methods for nonlinear acceleration.
28 : ! Numerical Linear Algebra with Applications, 16(3), 197–221. https://doi.org/10.1002/nla.617
29 : !
30 : ! Birken, P. (2015) Termination criteria for inexact fixed‐point schemes.
31 : ! Numer. Linear Algebra Appl., 22: 702– 716. doi: 10.1002/nla.1982.
32 : !
33 : ! authors: JF Lemieux, ECCC, Philppe Blain, ECCC
34 : !
35 :
36 : module ice_dyn_vp
37 :
38 : use ice_kinds_mod
39 : use ice_blocks, only: nx_block, ny_block
40 : use ice_boundary, only: ice_halo
41 : use ice_communicate, only: my_task, master_task, get_num_procs
42 : use ice_constants, only: field_loc_center, field_loc_NEcorner, &
43 : field_type_scalar, field_type_vector
44 : use ice_constants, only: c0, p027, p055, p111, p166, &
45 : p222, p25, p333, p5, c1
46 : use ice_domain, only: nblocks, distrb_info
47 : use ice_domain_size, only: max_blocks
48 : use ice_dyn_shared, only: dyn_prep1, dyn_prep2, dyn_finish, &
49 : cosw, sinw, fcor_blk, uvel_init, vvel_init, & ! LCOV_EXCL_LINE
50 : seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & ! LCOV_EXCL_LINE
51 : seabed_stress, Ktens, stack_fields, unstack_fields, fld2, fld3, fld4
52 : use ice_fileunits, only: nu_diag
53 : use ice_flux, only: fmU
54 : use ice_global_reductions, only: global_sum
55 : use ice_grid, only: dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, uarear
56 : use ice_exit, only: abort_ice
57 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
58 : use icepack_intfc, only: icepack_ice_strength, icepack_query_parameters
59 :
60 : implicit none
61 : private
62 : public :: implicit_solver, init_vp
63 :
64 : ! namelist parameters
65 :
66 : integer (kind=int_kind), public :: &
67 : maxits_nonlin , & ! max nb of iteration for nonlinear solver ! LCOV_EXCL_LINE
68 : dim_fgmres , & ! size of fgmres Krylov subspace ! LCOV_EXCL_LINE
69 : dim_pgmres , & ! size of pgmres Krylov subspace ! LCOV_EXCL_LINE
70 : maxits_fgmres , & ! max nb of iteration for fgmres ! LCOV_EXCL_LINE
71 : maxits_pgmres , & ! max nb of iteration for pgmres ! LCOV_EXCL_LINE
72 : fpfunc_andacc , & ! fixed point function for Anderson acceleration: ! LCOV_EXCL_LINE
73 : ! 1: g(x) = FMGRES(A(x),b(x)), 2: g(x) = x - A(x)x + b(x)
74 : dim_andacc , & ! size of Anderson minimization matrix (number of saved previous residuals)
75 : start_andacc ! acceleration delay factor (acceleration starts at this iteration)
76 :
77 : logical (kind=log_kind), public :: &
78 : monitor_nonlin , & ! print nonlinear residual norm ! LCOV_EXCL_LINE
79 : monitor_fgmres , & ! print fgmres residual norm ! LCOV_EXCL_LINE
80 : monitor_pgmres , & ! print pgmres residual norm ! LCOV_EXCL_LINE
81 : use_mean_vrel ! use mean of previous 2 iterates to compute vrel (see Hibler and Ackley 1983)
82 :
83 : real (kind=dbl_kind), public :: &
84 : reltol_nonlin , & ! nonlinear stopping criterion: reltol_nonlin*res(k=0) ! LCOV_EXCL_LINE
85 : reltol_fgmres , & ! fgmres stopping criterion: reltol_fgmres*res(k) ! LCOV_EXCL_LINE
86 : reltol_pgmres , & ! pgmres stopping criterion: reltol_pgmres*res(k) ! LCOV_EXCL_LINE
87 : damping_andacc , & ! damping factor for Anderson acceleration ! LCOV_EXCL_LINE
88 : reltol_andacc ! relative tolerance for Anderson acceleration
89 :
90 : character (len=char_len), public :: &
91 : precond , & ! preconditioner for fgmres: 'ident' (identity), 'diag' (diagonal), ! LCOV_EXCL_LINE
92 : ! 'pgmres' (Jacobi-preconditioned GMRES)
93 : algo_nonlin , & ! nonlinear algorithm: 'picard' (Picard iteration), 'anderson' (Anderson acceleration)
94 : ortho_type ! type of orthogonalization for FGMRES ('cgs' or 'mgs')
95 :
96 : ! module variables
97 :
98 : integer (kind=int_kind), allocatable :: &
99 : icellT(:) , & ! no. of cells where iceTmask = .true. ! LCOV_EXCL_LINE
100 : icellU(:) ! no. of cells where iceUmask = .true.
101 :
102 : integer (kind=int_kind), allocatable :: &
103 : indxTi(:,:) , & ! compressed index in i-direction ! LCOV_EXCL_LINE
104 : indxTj(:,:) , & ! compressed index in j-direction ! LCOV_EXCL_LINE
105 : indxUi(:,:) , & ! compressed index in i-direction ! LCOV_EXCL_LINE
106 : indxUj(:,:) ! compressed index in j-direction
107 :
108 : !=======================================================================
109 :
110 : contains
111 :
112 : !=======================================================================
113 :
114 : ! Initialize parameters and variables needed for the vp dynamics
115 : ! author: Philippe Blain, ECCC
116 :
117 0 : subroutine init_vp
118 :
119 : use ice_blocks, only: get_block, block
120 : use ice_boundary, only: ice_HaloUpdate
121 : use ice_constants, only: c1, &
122 : field_loc_center, field_type_scalar
123 : use ice_domain, only: blocks_ice, halo_info
124 : use ice_calendar, only: dt_dyn
125 : use ice_dyn_shared, only: init_dyn_shared
126 : ! use ice_grid, only: tarea
127 :
128 : ! local variables
129 :
130 : integer (kind=int_kind) :: &
131 : i, j, iblk, & ! LCOV_EXCL_LINE
132 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
133 :
134 : type (block) :: &
135 : this_block ! block information for current block
136 :
137 0 : call init_dyn_shared(dt_dyn)
138 :
139 : ! Initialize module variables
140 0 : allocate(icellT(max_blocks), icellU(max_blocks))
141 0 : allocate(indxTi(nx_block*ny_block, max_blocks), &
142 : indxTj(nx_block*ny_block, max_blocks), & ! LCOV_EXCL_LINE
143 : indxUi(nx_block*ny_block, max_blocks), & ! LCOV_EXCL_LINE
144 0 : indxUj(nx_block*ny_block, max_blocks))
145 :
146 0 : end subroutine init_vp
147 :
148 : !=======================================================================
149 :
150 : ! Viscous-plastic dynamics driver
151 : !
152 : #ifdef CICE_IN_NEMO
153 : ! Wind stress is set during this routine from the values supplied
154 : ! via NEMO (unless calc_strair is true). These values are supplied
155 : ! rotated on u grid and multiplied by aice. strairxT = 0 in this
156 : ! case so operations in dyn_prep1 are pointless but carried out to
157 : ! minimise code changes.
158 : #endif
159 : !
160 : ! author: JF Lemieux, A. Qaddouri and F. Dupont ECCC
161 :
162 0 : subroutine implicit_solver (dt)
163 :
164 : use ice_arrays_column, only: Cdn_ocn
165 : use ice_boundary, only: ice_HaloMask, ice_HaloUpdate, &
166 : ice_HaloDestroy, ice_HaloUpdate_stress
167 : use ice_blocks, only: block, get_block, nx_block, ny_block
168 : use ice_domain, only: blocks_ice, halo_info, maskhalo_dyn
169 : use ice_domain_size, only: max_blocks, ncat
170 : use ice_dyn_shared, only: deformations, iceTmask, iceUmask
171 : use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
172 : strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, & ! LCOV_EXCL_LINE
173 : strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, & ! LCOV_EXCL_LINE
174 : strax, stray, & ! LCOV_EXCL_LINE
175 : TbU, hwater, & ! LCOV_EXCL_LINE
176 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
177 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
178 : stress12_1, stress12_2, stress12_3, stress12_4
179 : use ice_grid, only: tmask, umask, dxT, dyT, cxp, cyp, cxm, cym, &
180 : tarear, grid_type, grid_average_X2Y, & ! LCOV_EXCL_LINE
181 : grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
182 : use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
183 : aice_init, aice0, aicen, vicen, strength
184 : use ice_timers, only: timer_dynamics, timer_bound, &
185 : ice_timer_start, ice_timer_stop
186 :
187 : real (kind=dbl_kind), intent(in) :: &
188 : dt ! time step
189 :
190 : ! local variables
191 :
192 : integer (kind=int_kind) :: &
193 : ntot , & ! size of problem for Anderson ! LCOV_EXCL_LINE
194 : iblk , & ! block index ! LCOV_EXCL_LINE
195 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
196 : i, j, ij
197 :
198 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
199 : uocnU , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
200 : vocnU , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
201 : ss_tltxU , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
202 : ss_tltyU , & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
203 : cdn_ocnU , & ! ocn drag coefficient ! LCOV_EXCL_LINE
204 : tmass , & ! total mass of ice and snow (kg/m^2) ! LCOV_EXCL_LINE
205 : waterxU , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
206 : wateryU , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
207 : forcexU , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
208 : forceyU , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
209 : bxfix , & ! part of bx that is constant during Picard ! LCOV_EXCL_LINE
210 : byfix , & ! part of by that is constant during Picard ! LCOV_EXCL_LINE
211 : Cb , & ! seabed stress coefficient ! LCOV_EXCL_LINE
212 : fpresx , & ! fixed point residual vector, x components: fx = uvel - uprev_k ! LCOV_EXCL_LINE
213 : fpresy , & ! fixed point residual vector, y components: fy = vvel - vprev_k ! LCOV_EXCL_LINE
214 : umass , & ! total mass of ice and snow (u grid) ! LCOV_EXCL_LINE
215 0 : umassdti ! mass of U-cell/dte (kg/m^2 s)
216 :
217 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4):: &
218 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
219 : etax2 , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
220 0 : rep_prs ! replacement pressure
221 :
222 : logical (kind=log_kind) :: calc_strair
223 :
224 : integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: &
225 0 : halomask ! generic halo mask
226 :
227 : type (ice_halo) :: &
228 : halo_info_mask ! ghost cell update info for masked halo
229 :
230 : type (block) :: &
231 : this_block ! block information for current block
232 :
233 : real (kind=dbl_kind), allocatable :: &
234 0 : sol(:) ! solution vector
235 :
236 : character(len=*), parameter :: subname = '(implicit_solver)'
237 :
238 0 : call ice_timer_start(timer_dynamics) ! dynamics
239 :
240 : !-----------------------------------------------------------------
241 : ! Initialize
242 : !-----------------------------------------------------------------
243 :
244 : ! This call is needed only if dt changes during runtime.
245 : ! call set_evp_parameters (dt)
246 :
247 : !-----------------------------------------------------------------
248 : ! boundary updates
249 : ! commented out because the ghost cells are freshly
250 : ! updated after cleanup_itd
251 : !-----------------------------------------------------------------
252 :
253 : ! call ice_timer_start(timer_bound)
254 : ! call ice_HaloUpdate (aice, halo_info, &
255 : ! field_loc_center, field_type_scalar)
256 : ! call ice_HaloUpdate (vice, halo_info, &
257 : ! field_loc_center, field_type_scalar)
258 : ! call ice_HaloUpdate (vsno, halo_info, &
259 : ! field_loc_center, field_type_scalar)
260 : ! call ice_timer_stop(timer_bound)
261 :
262 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
263 0 : do iblk = 1, nblocks
264 :
265 0 : do j = 1, ny_block
266 0 : do i = 1, nx_block
267 0 : rdg_conv (i,j,iblk) = c0
268 0 : rdg_shear(i,j,iblk) = c0
269 0 : divu (i,j,iblk) = c0
270 0 : shear(i,j,iblk) = c0
271 : enddo
272 : enddo
273 :
274 : !-----------------------------------------------------------------
275 : ! preparation for dynamics
276 : !-----------------------------------------------------------------
277 :
278 0 : this_block = get_block(blocks_ice(iblk),iblk)
279 0 : ilo = this_block%ilo
280 0 : ihi = this_block%ihi
281 0 : jlo = this_block%jlo
282 0 : jhi = this_block%jhi
283 :
284 : call dyn_prep1 (nx_block, ny_block, &
285 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
286 : aice (:,:,iblk), vice (:,:,iblk), & ! LCOV_EXCL_LINE
287 : vsno (:,:,iblk), tmask (:,:,iblk), & ! LCOV_EXCL_LINE
288 0 : tmass (:,:,iblk), iceTmask(:,:,iblk))
289 :
290 : enddo ! iblk
291 : !$OMP END PARALLEL DO
292 :
293 0 : call ice_timer_start(timer_bound)
294 : call ice_HaloUpdate (iceTmask, halo_info, &
295 0 : field_loc_center, field_type_scalar)
296 0 : call ice_timer_stop(timer_bound)
297 :
298 : !-----------------------------------------------------------------
299 : ! convert fields from T to U grid
300 : !-----------------------------------------------------------------
301 :
302 0 : call stack_fields(tmass, aice_init, cdn_ocn, fld3)
303 : call ice_HaloUpdate (fld3, halo_info, &
304 0 : field_loc_center, field_type_scalar)
305 0 : call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
306 : call ice_HaloUpdate (fld4, halo_info, &
307 0 : field_loc_center, field_type_vector)
308 0 : call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
309 0 : call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)
310 :
311 0 : call grid_average_X2Y('S',tmass , 'T' , umass , 'U')
312 0 : call grid_average_X2Y('S',aice_init, 'T' , aiU , 'U')
313 0 : call grid_average_X2Y('S',cdn_ocn , 'T' , cdn_ocnU, 'U')
314 0 : call grid_average_X2Y('S',uocn , grid_ocn_dynu, uocnU , 'U')
315 0 : call grid_average_X2Y('S',vocn , grid_ocn_dynv, vocnU , 'U')
316 0 : call grid_average_X2Y('S',ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
317 0 : call grid_average_X2Y('S',ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')
318 :
319 : !----------------------------------------------------------------
320 : ! Set wind stress to values supplied via NEMO or other forcing
321 : ! This wind stress is rotated on u grid and multiplied by aice
322 : !----------------------------------------------------------------
323 0 : call icepack_query_parameters(calc_strair_out=calc_strair)
324 0 : call icepack_warnings_flush(nu_diag)
325 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
326 0 : file=__FILE__, line=__LINE__)
327 :
328 0 : if (.not. calc_strair) then
329 0 : call grid_average_X2Y('F', strax, grid_atm_dynu, strairxU, 'U')
330 0 : call grid_average_X2Y('F', stray, grid_atm_dynv, strairyU, 'U')
331 : else
332 : call ice_HaloUpdate (strairxT, halo_info, &
333 0 : field_loc_center, field_type_vector)
334 : call ice_HaloUpdate (strairyT, halo_info, &
335 0 : field_loc_center, field_type_vector)
336 0 : call grid_average_X2Y('F',strairxT,'T',strairxU,'U')
337 0 : call grid_average_X2Y('F',strairyT,'T',strairyU,'U')
338 : endif
339 :
340 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,ij,i,j)
341 0 : do iblk = 1, nblocks
342 :
343 : !-----------------------------------------------------------------
344 : ! more preparation for dynamics
345 : !-----------------------------------------------------------------
346 :
347 0 : this_block = get_block(blocks_ice(iblk),iblk)
348 0 : ilo = this_block%ilo
349 0 : ihi = this_block%ihi
350 0 : jlo = this_block%jlo
351 0 : jhi = this_block%jhi
352 :
353 : call dyn_prep2 (nx_block, ny_block, &
354 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
355 : icellT(iblk), icellU(iblk), & ! LCOV_EXCL_LINE
356 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
357 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
358 : aiU (:,:,iblk), umass (:,:,iblk), & ! LCOV_EXCL_LINE
359 : umassdti (:,:,iblk), fcor_blk (:,:,iblk), & ! LCOV_EXCL_LINE
360 : umask (:,:,iblk), & ! LCOV_EXCL_LINE
361 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
362 : strairxU (:,:,iblk), strairyU (:,:,iblk), & ! LCOV_EXCL_LINE
363 : ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & ! LCOV_EXCL_LINE
364 : iceTmask (:,:,iblk), iceUmask (:,:,iblk), & ! LCOV_EXCL_LINE
365 : fmU (:,:,iblk), dt, & ! LCOV_EXCL_LINE
366 : strtltxU (:,:,iblk), strtltyU (:,:,iblk), & ! LCOV_EXCL_LINE
367 : strocnxU (:,:,iblk), strocnyU (:,:,iblk), & ! LCOV_EXCL_LINE
368 : strintxU (:,:,iblk), strintyU (:,:,iblk), & ! LCOV_EXCL_LINE
369 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
370 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
371 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
372 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
373 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
374 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
375 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
376 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
377 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
378 : uvel_init (:,:,iblk), vvel_init (:,:,iblk), & ! LCOV_EXCL_LINE
379 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
380 0 : TbU (:,:,iblk))
381 :
382 : call calc_bfix (nx_block , ny_block , &
383 : icellU(iblk) , & ! LCOV_EXCL_LINE
384 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
385 : umassdti (:,:,iblk), & ! LCOV_EXCL_LINE
386 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
387 : uvel_init (:,:,iblk), vvel_init (:,:,iblk), & ! LCOV_EXCL_LINE
388 0 : bxfix (:,:,iblk), byfix (:,:,iblk))
389 :
390 : !-----------------------------------------------------------------
391 : ! ice strength
392 : !-----------------------------------------------------------------
393 :
394 0 : strength(:,:,iblk) = c0 ! initialize
395 0 : do ij = 1, icellT(iblk)
396 0 : i = indxTi(ij, iblk)
397 0 : j = indxTj(ij, iblk)
398 : call icepack_ice_strength (ncat, &
399 : aice (i,j, iblk), & ! LCOV_EXCL_LINE
400 : vice (i,j, iblk), & ! LCOV_EXCL_LINE
401 : aice0 (i,j, iblk), & ! LCOV_EXCL_LINE
402 : aicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
403 : vicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
404 0 : strength(i,j, iblk))
405 : enddo ! ij
406 :
407 : enddo ! iblk
408 : !$OMP END PARALLEL DO
409 :
410 0 : call icepack_warnings_flush(nu_diag)
411 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
412 0 : file=__FILE__, line=__LINE__)
413 :
414 0 : call ice_timer_start(timer_bound)
415 : call ice_HaloUpdate (strength, halo_info, &
416 0 : field_loc_center, field_type_scalar)
417 : ! velocities may have changed in dyn_prep2
418 0 : call stack_fields(uvel, vvel, fld2)
419 : call ice_HaloUpdate (fld2, halo_info, &
420 0 : field_loc_NEcorner, field_type_vector)
421 0 : call unstack_fields(fld2, uvel, vvel)
422 0 : call ice_timer_stop(timer_bound)
423 :
424 0 : if (maskhalo_dyn) then
425 0 : call ice_timer_start(timer_bound)
426 0 : halomask = 0
427 0 : where (iceUmask) halomask = 1
428 0 : call ice_HaloUpdate (halomask, halo_info, &
429 0 : field_loc_center, field_type_scalar)
430 0 : call ice_timer_stop(timer_bound)
431 0 : call ice_HaloMask(halo_info_mask, halo_info, halomask)
432 : endif
433 :
434 : !-----------------------------------------------------------------
435 : ! seabed stress factor TbU (TbU is part of Cb coefficient)
436 : !-----------------------------------------------------------------
437 0 : if (seabed_stress) then
438 0 : if ( seabed_stress_method == 'LKD' ) then
439 0 : !$OMP PARALLEL DO PRIVATE(iblk)
440 0 : do iblk = 1, nblocks
441 : call seabed_stress_factor_LKD (nx_block, ny_block, &
442 : icellU (iblk), & ! LCOV_EXCL_LINE
443 : indxUi(:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
444 : vice(:,:,iblk), aice(:,:,iblk), & ! LCOV_EXCL_LINE
445 0 : hwater(:,:,iblk), TbU(:,:,iblk))
446 : enddo
447 : !$OMP END PARALLEL DO
448 :
449 0 : elseif ( seabed_stress_method == 'probabilistic' ) then
450 0 : !$OMP PARALLEL DO PRIVATE(iblk)
451 0 : do iblk = 1, nblocks
452 :
453 : call seabed_stress_factor_prob (nx_block, ny_block, &
454 : icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & ! LCOV_EXCL_LINE
455 : icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
456 : aicen(:,:,:,iblk), vicen(:,:,:,iblk), & ! LCOV_EXCL_LINE
457 0 : hwater(:,:,iblk), TbU(:,:,iblk))
458 : enddo
459 : !$OMP END PARALLEL DO
460 :
461 : endif
462 : endif
463 :
464 :
465 : !-----------------------------------------------------------------
466 : ! calc size of problem (ntot) and allocate solution vector
467 : !-----------------------------------------------------------------
468 :
469 0 : ntot = 0
470 0 : do iblk = 1, nblocks
471 0 : ntot = ntot + icellU(iblk)
472 : enddo
473 0 : ntot = 2 * ntot ! times 2 because of u and v
474 :
475 0 : allocate(sol(ntot))
476 :
477 : !-----------------------------------------------------------------
478 : ! Start of nonlinear iteration
479 : !-----------------------------------------------------------------
480 : call anderson_solver (icellT , icellU , &
481 : indxTi , indxTj , & ! LCOV_EXCL_LINE
482 : indxUi , indxUj , & ! LCOV_EXCL_LINE
483 : aiU , ntot , & ! LCOV_EXCL_LINE
484 : uocnU , vocnU , & ! LCOV_EXCL_LINE
485 : waterxU , wateryU, & ! LCOV_EXCL_LINE
486 : bxfix , byfix , & ! LCOV_EXCL_LINE
487 : umassdti, sol , & ! LCOV_EXCL_LINE
488 : fpresx , fpresy , & ! LCOV_EXCL_LINE
489 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
490 : rep_prs , cdn_ocnU,& ! LCOV_EXCL_LINE
491 0 : Cb, halo_info_mask)
492 : !-----------------------------------------------------------------
493 : ! End of nonlinear iteration
494 : !-----------------------------------------------------------------
495 :
496 0 : deallocate(sol)
497 :
498 0 : if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)
499 :
500 : !-----------------------------------------------------------------
501 : ! Compute stresses
502 : !-----------------------------------------------------------------
503 0 : !$OMP PARALLEL DO PRIVATE(iblk)
504 0 : do iblk = 1, nblocks
505 : call stress_vp (nx_block , ny_block , &
506 : icellT(iblk) , & ! LCOV_EXCL_LINE
507 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
508 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
509 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
510 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
511 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
512 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
513 : rep_prs (:,:,iblk,:), & ! LCOV_EXCL_LINE
514 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
515 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
516 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
517 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
518 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
519 0 : stress12_3(:,:,iblk), stress12_4(:,:,iblk))
520 : enddo ! iblk
521 : !$OMP END PARALLEL DO
522 :
523 : !-----------------------------------------------------------------
524 : ! Compute deformations
525 : !-----------------------------------------------------------------
526 0 : !$OMP PARALLEL DO PRIVATE(iblk)
527 0 : do iblk = 1, nblocks
528 : call deformations (nx_block , ny_block , &
529 : icellT(iblk) , & ! LCOV_EXCL_LINE
530 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
531 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
532 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
533 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
534 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
535 : tarear (:,:,iblk), & ! LCOV_EXCL_LINE
536 : shear (:,:,iblk), divu (:,:,iblk), & ! LCOV_EXCL_LINE
537 0 : rdg_conv (:,:,iblk), rdg_shear (:,:,iblk))
538 : enddo
539 : !$OMP END PARALLEL DO
540 :
541 : !-----------------------------------------------------------------
542 : ! Compute seabed stress (diagnostic)
543 : !-----------------------------------------------------------------
544 0 : if (seabed_stress) then
545 0 : !$OMP PARALLEL DO PRIVATE(iblk)
546 0 : do iblk = 1, nblocks
547 : call calc_seabed_stress (nx_block , ny_block , &
548 : icellU(iblk) , & ! LCOV_EXCL_LINE
549 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
550 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
551 : Cb (:,:,iblk), & ! LCOV_EXCL_LINE
552 0 : taubxU (:,:,iblk), taubyU (:,:,iblk))
553 : enddo
554 : !$OMP END PARALLEL DO
555 : endif
556 :
557 : ! Force symmetry across the tripole seam
558 0 : if (trim(grid_type) == 'tripole') then
559 0 : if (maskhalo_dyn) then
560 : !-------------------------------------------------------
561 : ! set halomask to zero because ice_HaloMask always keeps
562 : ! local copies AND tripole zipper communication
563 : !-------------------------------------------------------
564 0 : halomask = 0
565 0 : call ice_HaloMask(halo_info_mask, halo_info, halomask)
566 :
567 : call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info_mask, &
568 0 : field_loc_center, field_type_scalar)
569 : call ice_HaloUpdate_stress(stressp_3, stressp_1, halo_info_mask, &
570 0 : field_loc_center, field_type_scalar)
571 : call ice_HaloUpdate_stress(stressp_2, stressp_4, halo_info_mask, &
572 0 : field_loc_center, field_type_scalar)
573 : call ice_HaloUpdate_stress(stressp_4, stressp_2, halo_info_mask, &
574 0 : field_loc_center, field_type_scalar)
575 :
576 : call ice_HaloUpdate_stress(stressm_1, stressm_3, halo_info_mask, &
577 0 : field_loc_center, field_type_scalar)
578 : call ice_HaloUpdate_stress(stressm_3, stressm_1, halo_info_mask, &
579 0 : field_loc_center, field_type_scalar)
580 : call ice_HaloUpdate_stress(stressm_2, stressm_4, halo_info_mask, &
581 0 : field_loc_center, field_type_scalar)
582 : call ice_HaloUpdate_stress(stressm_4, stressm_2, halo_info_mask, &
583 0 : field_loc_center, field_type_scalar)
584 :
585 : call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info_mask, &
586 0 : field_loc_center, field_type_scalar)
587 : call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info_mask, &
588 0 : field_loc_center, field_type_scalar)
589 : call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info_mask, &
590 0 : field_loc_center, field_type_scalar)
591 : call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info_mask, &
592 0 : field_loc_center, field_type_scalar)
593 :
594 0 : call ice_HaloDestroy(halo_info_mask)
595 : else
596 : call ice_HaloUpdate_stress(stressp_1, stressp_3, halo_info, &
597 0 : field_loc_center, field_type_scalar)
598 : call ice_HaloUpdate_stress(stressp_3, stressp_1, halo_info, &
599 0 : field_loc_center, field_type_scalar)
600 : call ice_HaloUpdate_stress(stressp_2, stressp_4, halo_info, &
601 0 : field_loc_center, field_type_scalar)
602 : call ice_HaloUpdate_stress(stressp_4, stressp_2, halo_info, &
603 0 : field_loc_center, field_type_scalar)
604 :
605 : call ice_HaloUpdate_stress(stressm_1, stressm_3, halo_info, &
606 0 : field_loc_center, field_type_scalar)
607 : call ice_HaloUpdate_stress(stressm_3, stressm_1, halo_info, &
608 0 : field_loc_center, field_type_scalar)
609 : call ice_HaloUpdate_stress(stressm_2, stressm_4, halo_info, &
610 0 : field_loc_center, field_type_scalar)
611 : call ice_HaloUpdate_stress(stressm_4, stressm_2, halo_info, &
612 0 : field_loc_center, field_type_scalar)
613 :
614 : call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info, &
615 0 : field_loc_center, field_type_scalar)
616 : call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info, &
617 0 : field_loc_center, field_type_scalar)
618 : call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info, &
619 0 : field_loc_center, field_type_scalar)
620 : call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info, &
621 0 : field_loc_center, field_type_scalar)
622 : endif ! maskhalo
623 : endif ! tripole
624 :
625 : !-----------------------------------------------------------------
626 : ! ice-ocean stress
627 : !-----------------------------------------------------------------
628 :
629 0 : !$OMP PARALLEL DO PRIVATE(iblk)
630 0 : do iblk = 1, nblocks
631 :
632 : call dyn_finish &
633 : (nx_block, ny_block, & ! LCOV_EXCL_LINE
634 : icellU (iblk), Cdn_ocnU(:,:,iblk), & ! LCOV_EXCL_LINE
635 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
636 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
637 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
638 : aiU (:,:,iblk), fmU (:,:,iblk), & ! LCOV_EXCL_LINE
639 : ! strintxU(:,:,iblk), strintyU(:,:,iblk), & ! LCOV_EXCL_LINE
640 : ! strairxU(:,:,iblk), strairyU(:,:,iblk), & ! LCOV_EXCL_LINE
641 0 : strocnxU(:,:,iblk), strocnyU(:,:,iblk))
642 :
643 : enddo
644 : !$OMP END PARALLEL DO
645 :
646 : ! shift velocity components from CD grid locations (N, E) to B grid location (U) for transport
647 : ! commented out in order to focus on EVP for now within the cdgrid
648 : ! should be used when routine is ready
649 : ! if (grid_ice == 'CD' .or. grid_ice == 'C') then
650 : ! call grid_average_X2Y('E2US',uvelE,uvel)
651 : ! call grid_average_X2Y('N2US',vvelN,vvel)
652 : ! endif
653 : !end comment out
654 0 : call ice_timer_stop(timer_dynamics) ! dynamics
655 :
656 0 : end subroutine implicit_solver
657 :
658 : !=======================================================================
659 :
660 : ! Solve the nonlinear equation F(u,v) = 0, where
661 : ! F(u,v) := A(u,v) * (u,v) - b(u,v)
662 : ! using Anderson acceleration (accelerated fixed point (Picard) iteration)
663 : !
664 : ! author: JF Lemieux, A. Qaddouri, F. Dupont and P. Blain ECCC
665 : !
666 : ! Anderson algorithm adadpted from:
667 : ! H. F. Walker, “Anderson Acceleration: Algorithms and Implementations”
668 : ! [Online]. Available: https://users.wpi.edu/~walker/Papers/anderson_accn_algs_imps.pdf
669 :
670 0 : subroutine anderson_solver (icellT , icellU , &
671 : indxTi , indxTj , & ! LCOV_EXCL_LINE
672 : indxUi , indxUj , & ! LCOV_EXCL_LINE
673 : aiU , ntot , & ! LCOV_EXCL_LINE
674 : uocn , vocn , & ! LCOV_EXCL_LINE
675 : waterxU , wateryU, & ! LCOV_EXCL_LINE
676 : bxfix , byfix , & ! LCOV_EXCL_LINE
677 : umassdti, sol , & ! LCOV_EXCL_LINE
678 : fpresx , fpresy , & ! LCOV_EXCL_LINE
679 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
680 : rep_prs , cdn_ocn, & ! LCOV_EXCL_LINE
681 0 : Cb, halo_info_mask)
682 :
683 : use ice_blocks, only: nx_block, ny_block
684 : use ice_boundary, only: ice_HaloUpdate
685 : use ice_constants, only: c1
686 : use ice_domain, only: maskhalo_dyn, halo_info
687 : use ice_domain_size, only: max_blocks
688 : use ice_flux, only: fmU, TbU
689 : use ice_grid, only: dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
690 : uarear
691 : use ice_dyn_shared, only: DminTarea
692 : use ice_state, only: uvel, vvel, strength
693 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
694 :
695 : integer (kind=int_kind), intent(in) :: &
696 : ntot ! size of problem for Anderson
697 :
698 : integer (kind=int_kind), dimension(max_blocks), intent(in) :: &
699 : icellT , & ! no. of cells where iceTmask = .true. ! LCOV_EXCL_LINE
700 : icellU ! no. of cells where iceUmask = .true.
701 :
702 : integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: &
703 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
704 : indxTj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
705 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
706 : indxUj ! compressed index in j-direction
707 :
708 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
709 : aiU , & ! ice fraction on u-grid ! LCOV_EXCL_LINE
710 : uocn , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
711 : vocn , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
712 : cdn_ocn , & ! ocn drag coefficient ! LCOV_EXCL_LINE
713 : waterxU , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
714 : wateryU , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
715 : bxfix , & ! part of bx that is constant during Picard ! LCOV_EXCL_LINE
716 : byfix , & ! part of by that is constant during Picard ! LCOV_EXCL_LINE
717 : umassdti ! mass of U-cell/dte (kg/m^2 s)
718 :
719 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(out) :: &
720 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
721 : etax2 , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
722 : rep_prs ! replacement pressure
723 :
724 : type (ice_halo), intent(in) :: &
725 : halo_info_mask ! ghost cell update info for masked halo
726 :
727 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
728 : fpresx , & ! fixed point residual vector, x components: fx = uvel - uprev_k ! LCOV_EXCL_LINE
729 : fpresy , & ! fixed point residual vector, y components: fy = vvel - vprev_k ! LCOV_EXCL_LINE
730 : Cb ! seabed stress coefficient
731 :
732 : real (kind=dbl_kind), dimension (ntot), intent(inout) :: &
733 : sol ! current approximate solution
734 :
735 : ! local variables
736 :
737 : integer (kind=int_kind) :: &
738 : it_nl , & ! nonlinear loop iteration index ! LCOV_EXCL_LINE
739 : res_num , & ! current number of stored residuals ! LCOV_EXCL_LINE
740 : j , & ! iteration index for QR update ! LCOV_EXCL_LINE
741 : iblk , & ! block index ! LCOV_EXCL_LINE
742 : nbiter ! number of FGMRES iterations performed
743 :
744 : integer (kind=int_kind), parameter :: &
745 : inc = 1 ! increment value for BLAS calls
746 :
747 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
748 : uprev_k , & ! uvel at previous Picard iteration ! LCOV_EXCL_LINE
749 : vprev_k , & ! vvel at previous Picard iteration ! LCOV_EXCL_LINE
750 : ulin , & ! uvel to linearize vrel ! LCOV_EXCL_LINE
751 : vlin , & ! vvel to linearize vrel ! LCOV_EXCL_LINE
752 : vrel , & ! coeff for tauw ! LCOV_EXCL_LINE
753 : bx , & ! b vector ! LCOV_EXCL_LINE
754 : by , & ! b vector ! LCOV_EXCL_LINE
755 : diagx , & ! Diagonal (x component) of the matrix A ! LCOV_EXCL_LINE
756 : diagy , & ! Diagonal (y component) of the matrix A ! LCOV_EXCL_LINE
757 : Au , & ! matvec, Fx = bx - Au ! LCOV_EXCL_LINE
758 : Av , & ! matvec, Fy = by - Av ! LCOV_EXCL_LINE
759 : Fx , & ! x residual vector, Fx = bx - Au ! LCOV_EXCL_LINE
760 : Fy , & ! y residual vector, Fy = by - Av ! LCOV_EXCL_LINE
761 : solx , & ! solution of FGMRES (x components) ! LCOV_EXCL_LINE
762 0 : soly ! solution of FGMRES (y components)
763 :
764 : real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
765 : stress_Pr, & ! x,y-derivatives of the replacement pressure ! LCOV_EXCL_LINE
766 0 : diag_rheo ! contributions of the rhelogy term to the diagonal
767 :
768 : real (kind=dbl_kind), dimension (ntot) :: &
769 : res , & ! current residual ! LCOV_EXCL_LINE
770 : res_old , & ! previous residual ! LCOV_EXCL_LINE
771 : res_diff , & ! difference between current and previous residuals ! LCOV_EXCL_LINE
772 : fpfunc , & ! current value of fixed point function ! LCOV_EXCL_LINE
773 : fpfunc_old , & ! previous value of fixed point function ! LCOV_EXCL_LINE
774 : tmp ! temporary vector for BLAS calls
775 :
776 : real (kind=dbl_kind), dimension(ntot,dim_andacc) :: &
777 : Q , & ! Q factor for QR factorization of F (residuals) matrix ! LCOV_EXCL_LINE
778 : G_diff ! Matrix containing the differences of g(x) (fixed point function) evaluations
779 :
780 : real (kind=dbl_kind), dimension(dim_andacc,dim_andacc) :: &
781 : R ! R factor for QR factorization of F (residuals) matrix
782 :
783 : real (kind=dbl_kind), dimension(dim_andacc) :: &
784 : rhs_tri , & ! right hand side vector for matrix-vector product ! LCOV_EXCL_LINE
785 : coeffs ! coeffs used to combine previous solutions
786 :
787 : real (kind=dbl_kind) :: &
788 : ! tol , & ! tolerance for fixed point convergence: reltol_andacc * (initial fixed point residual norm) [unused for now] ! LCOV_EXCL_LINE
789 : tol_nl , & ! tolerance for nonlinear convergence: reltol_nonlin * (initial nonlinear residual norm) ! LCOV_EXCL_LINE
790 : fpres_norm , & ! norm of current fixed point residual : f(x) = g(x) - x ! LCOV_EXCL_LINE
791 : prog_norm , & ! norm of difference between current and previous solution ! LCOV_EXCL_LINE
792 0 : nlres_norm ! norm of current nonlinear residual : F(x) = A(x)x -b(x)
793 :
794 : #ifdef USE_LAPACK
795 : real (kind=dbl_kind) :: &
796 : ddot, dnrm2 ! external BLAS functions
797 : #endif
798 :
799 : character(len=*), parameter :: subname = '(anderson_solver)'
800 :
801 : ! Initialization
802 0 : res_num = 0
803 :
804 0 : !$OMP PARALLEL DO PRIVATE(iblk)
805 0 : do iblk = 1, nblocks
806 0 : uprev_k(:,:,iblk) = uvel(:,:,iblk)
807 0 : vprev_k(:,:,iblk) = vvel(:,:,iblk)
808 : enddo
809 : !$OMP END PARALLEL DO
810 :
811 : ! Start iterations
812 0 : do it_nl = 0, maxits_nonlin ! nonlinear iteration loop
813 : ! Compute quantities needed for computing PDE residual and g(x) (fixed point map)
814 : !-----------------------------------------------------------------
815 : ! Calc zetax2, etax2, dPr/dx, dPr/dy, Cb and vrel = f(uprev_k, vprev_k)
816 : !-----------------------------------------------------------------
817 0 : !$OMP PARALLEL DO PRIVATE(iblk,stress_Pr)
818 0 : do iblk = 1, nblocks
819 :
820 0 : if (use_mean_vrel) then
821 0 : ulin(:,:,iblk) = p5*uprev_k(:,:,iblk) + p5*uvel(:,:,iblk)
822 0 : vlin(:,:,iblk) = p5*vprev_k(:,:,iblk) + p5*vvel(:,:,iblk)
823 : else
824 0 : ulin(:,:,iblk) = uvel(:,:,iblk)
825 0 : vlin(:,:,iblk) = vvel(:,:,iblk)
826 : endif
827 0 : uprev_k(:,:,iblk) = uvel(:,:,iblk)
828 0 : vprev_k(:,:,iblk) = vvel(:,:,iblk)
829 :
830 : call calc_zeta_dPr (nx_block , ny_block , &
831 : icellT (iblk), & ! LCOV_EXCL_LINE
832 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
833 : uprev_k (:,:,iblk), vprev_k (:,:,iblk), & ! LCOV_EXCL_LINE
834 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
835 : dxhy (:,:,iblk), dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
836 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
837 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
838 : DminTarea (:,:,iblk),strength (:,:,iblk),& ! LCOV_EXCL_LINE
839 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:),& ! LCOV_EXCL_LINE
840 0 : rep_prs(:,:,iblk,:), stress_Pr (:,:,:))
841 :
842 : call calc_vrel_Cb (nx_block , ny_block , &
843 : icellU (iblk), Cdn_ocn (:,:,iblk), & ! LCOV_EXCL_LINE
844 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
845 : aiU (:,:,iblk), TbU (:,:,iblk), & ! LCOV_EXCL_LINE
846 : uocn (:,:,iblk), vocn (:,:,iblk), & ! LCOV_EXCL_LINE
847 : ulin (:,:,iblk), vlin (:,:,iblk), & ! LCOV_EXCL_LINE
848 0 : vrel (:,:,iblk), Cb (:,:,iblk))
849 :
850 : ! prepare b vector (RHS)
851 : call calc_bvec (nx_block , ny_block , &
852 : icellU (iblk), & ! LCOV_EXCL_LINE
853 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
854 : stress_Pr (:,:,:), uarear (:,:,iblk), & ! LCOV_EXCL_LINE
855 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
856 : bxfix (:,:,iblk), byfix (:,:,iblk), & ! LCOV_EXCL_LINE
857 : bx (:,:,iblk), by (:,:,iblk), & ! LCOV_EXCL_LINE
858 0 : vrel (:,:,iblk))
859 :
860 : ! Compute nonlinear residual norm (PDE residual)
861 : call matvec (nx_block , ny_block , &
862 : icellU (iblk) , icellT (iblk), & ! LCOV_EXCL_LINE
863 : indxUi (:,iblk) , indxUj (:,iblk), & ! LCOV_EXCL_LINE
864 : indxTi (:,iblk) , indxTj (:,iblk), & ! LCOV_EXCL_LINE
865 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
866 : dxhy (:,:,iblk) , dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
867 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
868 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
869 : uprev_k (:,:,iblk) , vprev_k (:,:,iblk), & ! LCOV_EXCL_LINE
870 : vrel (:,:,iblk) , Cb (:,:,iblk), & ! LCOV_EXCL_LINE
871 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
872 : umassdti (:,:,iblk) , fmU (:,:,iblk), & ! LCOV_EXCL_LINE
873 : uarear (:,:,iblk) , & ! LCOV_EXCL_LINE
874 0 : Au (:,:,iblk) , Av (:,:,iblk))
875 : call residual_vec (nx_block , ny_block , &
876 : icellU (iblk), & ! LCOV_EXCL_LINE
877 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
878 : bx (:,:,iblk), by (:,:,iblk), & ! LCOV_EXCL_LINE
879 : Au (:,:,iblk), Av (:,:,iblk), & ! LCOV_EXCL_LINE
880 0 : Fx (:,:,iblk), Fy (:,:,iblk))
881 : enddo
882 : !$OMP END PARALLEL DO
883 : nlres_norm = global_norm(nx_block, ny_block, &
884 : icellU , & ! LCOV_EXCL_LINE
885 : indxUi , indxUj , & ! LCOV_EXCL_LINE
886 0 : Fx , Fy )
887 0 : if (my_task == master_task .and. monitor_nonlin) then
888 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_nonlin: iter_nonlin= ", it_nl, &
889 0 : " nonlin_res_L2norm= ", nlres_norm
890 : endif
891 : ! Compute relative tolerance at first iteration
892 0 : if (it_nl == 0) then
893 0 : tol_nl = reltol_nonlin*nlres_norm
894 : endif
895 :
896 : ! Check for nonlinear convergence
897 0 : if (nlres_norm < tol_nl) then
898 0 : exit
899 : endif
900 :
901 : ! Put initial guess for FGMRES in solx,soly and sol (needed for anderson)
902 0 : solx = uprev_k
903 0 : soly = vprev_k
904 : call arrays_to_vec (nx_block , ny_block , &
905 : nblocks , max_blocks , & ! LCOV_EXCL_LINE
906 : icellU (:), ntot , & ! LCOV_EXCL_LINE
907 : indxUi (:,:), indxUj (:,:), & ! LCOV_EXCL_LINE
908 : uprev_k (:,:,:), vprev_k (:,:,:), & ! LCOV_EXCL_LINE
909 0 : sol (:))
910 :
911 : ! Compute fixed point map g(x)
912 0 : if (fpfunc_andacc == 1) then
913 : ! g_1(x) = FGMRES(A(x), b(x))
914 :
915 : ! Prepare diagonal for preconditioner
916 0 : if (precond == 'diag' .or. precond == 'pgmres') then
917 0 : !$OMP PARALLEL DO PRIVATE(iblk,diag_rheo)
918 0 : do iblk = 1, nblocks
919 : ! first compute diagonal contributions due to rheology term
920 : call formDiag_step1 (nx_block , ny_block , &
921 : icellU (iblk) , & ! LCOV_EXCL_LINE
922 : indxUi (:,iblk) , indxUj(:,iblk), & ! LCOV_EXCL_LINE
923 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
924 : dxhy (:,:,iblk) , dyhx(:,:,iblk), & ! LCOV_EXCL_LINE
925 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
926 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
927 : zetax2 (:,:,iblk,:), etax2(:,:,iblk,:), & ! LCOV_EXCL_LINE
928 0 : diag_rheo(:,:,:))
929 : ! second compute the full diagonal
930 : call formDiag_step2 (nx_block , ny_block , &
931 : icellU (iblk), & ! LCOV_EXCL_LINE
932 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
933 : diag_rheo (:,:,:), vrel (:,:,iblk), & ! LCOV_EXCL_LINE
934 : umassdti (:,:,iblk), & ! LCOV_EXCL_LINE
935 : uarear (:,:,iblk), Cb (:,:,iblk), & ! LCOV_EXCL_LINE
936 0 : diagx (:,:,iblk), diagy (:,:,iblk))
937 : enddo
938 : !$OMP END PARALLEL DO
939 : endif
940 :
941 : ! FGMRES linear solver
942 : call fgmres (zetax2 , etax2 , &
943 : Cb , vrel , & ! LCOV_EXCL_LINE
944 : umassdti , & ! LCOV_EXCL_LINE
945 : halo_info_mask, & ! LCOV_EXCL_LINE
946 : bx , by , & ! LCOV_EXCL_LINE
947 : diagx , diagy , & ! LCOV_EXCL_LINE
948 : reltol_fgmres , dim_fgmres, & ! LCOV_EXCL_LINE
949 : maxits_fgmres , & ! LCOV_EXCL_LINE
950 : solx , soly , & ! LCOV_EXCL_LINE
951 0 : nbiter)
952 : ! Put FGMRES solution solx,soly in fpfunc vector (needed for Anderson)
953 : call arrays_to_vec (nx_block , ny_block , &
954 : nblocks , max_blocks , & ! LCOV_EXCL_LINE
955 : icellU (:), ntot , & ! LCOV_EXCL_LINE
956 : indxUi (:,:), indxUj (:,:), & ! LCOV_EXCL_LINE
957 : solx (:,:,:), soly (:,:,:), & ! LCOV_EXCL_LINE
958 0 : fpfunc (:))
959 0 : elseif (fpfunc_andacc == 2) then
960 : ! g_2(x) = x - A(x)x + b(x) = x - F(x)
961 : call abort_ice(error_message=subname // " Fixed point function g_2(x) not yet implemented (fpfunc_andacc = 2)" , &
962 0 : file=__FILE__, line=__LINE__)
963 : endif
964 :
965 : ! Compute fixed point residual f(x) = g(x) - x
966 0 : res = fpfunc - sol
967 : #ifdef USE_LAPACK
968 : fpres_norm = global_sum(dnrm2(size(res), res, inc)**2, distrb_info)
969 : #else
970 : call vec_to_arrays (nx_block , ny_block , &
971 : nblocks , max_blocks , & ! LCOV_EXCL_LINE
972 : icellU (:), ntot , & ! LCOV_EXCL_LINE
973 : indxUi (:,:), indxUj(:,:) , & ! LCOV_EXCL_LINE
974 : res (:), & ! LCOV_EXCL_LINE
975 0 : fpresx (:,:,:), fpresy (:,:,:))
976 : fpres_norm = global_norm(nx_block, ny_block, &
977 : icellU , & ! LCOV_EXCL_LINE
978 : indxUi , indxUj , & ! LCOV_EXCL_LINE
979 0 : fpresx , fpresy )
980 : #endif
981 0 : if (my_task == master_task .and. monitor_nonlin) then
982 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_nonlin: iter_nonlin= ", it_nl, &
983 0 : " fixed_point_res_L2norm= ", fpres_norm
984 : endif
985 :
986 : ! Not used for now (only nonlinear residual is checked)
987 : ! ! Store initial residual norm
988 : ! if (it_nl == 0) then
989 : ! tol = reltol_andacc*fpres_norm
990 : ! endif
991 : !
992 : ! ! Check residual
993 : ! if (fpres_norm < tol) then
994 : ! exit
995 : ! endif
996 :
997 0 : if (dim_andacc == 0 .or. it_nl < start_andacc) then
998 : ! Simple fixed point (Picard) iteration in this case
999 0 : sol = fpfunc
1000 : else
1001 : #ifdef USE_LAPACK
1002 : ! Begin Anderson acceleration
1003 : if (get_num_procs() > 1) then
1004 : ! Anderson solver is not yet parallelized; abort
1005 : if (my_task == master_task) then
1006 : call abort_ice(error_message=subname // " Anderson solver (algo_nonlin = 'anderson') is not yet parallelized, and nprocs > 1 " , &
1007 : file=__FILE__, line=__LINE__)
1008 : endif
1009 : endif
1010 : if (it_nl > start_andacc) then
1011 : ! Update residual difference vector
1012 : res_diff = res - res_old
1013 : ! Update fixed point function difference matrix
1014 : if (res_num < dim_andacc) then
1015 : ! Add column
1016 : G_diff(:,res_num+1) = fpfunc - fpfunc_old
1017 : else
1018 : ! Delete first column and add column
1019 : G_diff(:,1:res_num-1) = G_diff(:,2:res_num)
1020 : G_diff(:,res_num) = fpfunc - fpfunc_old
1021 : endif
1022 : res_num = res_num + 1
1023 : endif
1024 : res_old = res
1025 : fpfunc_old = fpfunc
1026 : if (res_num == 0) then
1027 : sol = fpfunc
1028 : else
1029 : if (res_num == 1) then
1030 : ! Initialize QR factorization
1031 : R(1,1) = dnrm2(size(res_diff), res_diff, inc)
1032 : Q(:,1) = res_diff/R(1,1)
1033 : else
1034 : if (res_num > dim_andacc) then
1035 : ! Update factorization since 1st column was deleted
1036 : call qr_delete(Q,R)
1037 : res_num = res_num - 1
1038 : endif
1039 : ! Update QR factorization for new column
1040 : do j = 1, res_num - 1
1041 : R(j,res_num) = ddot(ntot, Q(:,j), inc, res_diff, inc)
1042 : res_diff = res_diff - R(j,res_num) * Q(:,j)
1043 : enddo
1044 : R(res_num, res_num) = dnrm2(size(res_diff) ,res_diff, inc)
1045 : Q(:,res_num) = res_diff / R(res_num, res_num)
1046 : endif
1047 : ! TODO: here, drop more columns to improve conditioning
1048 : ! if (droptol) then
1049 :
1050 : ! endif
1051 : ! Solve least square problem for coefficients
1052 : ! 1. Compute rhs_tri = Q^T * res
1053 : call dgemv ('t', size(Q,1), res_num, c1, Q(:,1:res_num), size(Q,1), res, inc, c0, rhs_tri, inc)
1054 : ! 2. Solve R*coeffs = rhs_tri, put result in rhs_tri
1055 : call dtrsv ('u', 'n', 'n', res_num, R(1:res_num,1:res_num), res_num, rhs_tri, inc)
1056 : coeffs = rhs_tri
1057 : ! Update approximate solution: x = fpfunc - G_diff*coeffs, put result in fpfunc
1058 : call dgemv ('n', size(G_diff,1), res_num, -c1, G_diff(:,1:res_num), size(G_diff,1), coeffs, inc, c1, fpfunc, inc)
1059 : sol = fpfunc
1060 : ! Apply damping
1061 : if (damping_andacc > 0 .and. damping_andacc /= 1) then
1062 : ! x = x - (1-beta) (res - Q*R*coeffs)
1063 :
1064 : ! tmp = R*coeffs
1065 : call dgemv ('n', res_num, res_num, c1, R(1:res_num,1:res_num), res_num, coeffs, inc, c0, tmp, inc)
1066 : ! res = res - Q*tmp
1067 : call dgemv ('n', size(Q,1), res_num, -c1, Q(:,1:res_num), size(Q,1), tmp, inc, c1, res, inc)
1068 : ! x = x - (1-beta)*res
1069 : sol = sol - (1-damping_andacc)*res
1070 : endif
1071 : endif
1072 : #else
1073 : ! Anderson solver is not usable without LAPACK; abort
1074 : call abort_ice(error_message=subname // " CICE was not compiled with LAPACK, "// &
1075 : "and Anderson solver was chosen (algo_nonlin = 'anderson')" , & ! LCOV_EXCL_LINE
1076 0 : file=__FILE__, line=__LINE__)
1077 : #endif
1078 : endif
1079 :
1080 : !-----------------------------------------------------------------------
1081 : ! Put vector sol in uvel and vvel arrays
1082 : !-----------------------------------------------------------------------
1083 : call vec_to_arrays (nx_block , ny_block , &
1084 : nblocks , max_blocks , & ! LCOV_EXCL_LINE
1085 : icellU (:), ntot , & ! LCOV_EXCL_LINE
1086 : indxUi (:,:), indxUj (:,:), & ! LCOV_EXCL_LINE
1087 : sol (:), & ! LCOV_EXCL_LINE
1088 0 : uvel (:,:,:), vvel (:,:,:))
1089 :
1090 : ! Do halo update so that halo cells contain up to date info for advection
1091 0 : call stack_fields(uvel, vvel, fld2)
1092 0 : call ice_timer_start(timer_bound)
1093 0 : if (maskhalo_dyn) then
1094 : call ice_HaloUpdate (fld2, halo_info_mask, &
1095 0 : field_loc_NEcorner, field_type_vector)
1096 : else
1097 : call ice_HaloUpdate (fld2, halo_info, &
1098 0 : field_loc_NEcorner, field_type_vector)
1099 : endif
1100 0 : call ice_timer_stop(timer_bound)
1101 0 : call unstack_fields(fld2, uvel, vvel)
1102 :
1103 : ! Compute "progress" residual norm
1104 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1105 0 : do iblk = 1, nblocks
1106 0 : fpresx(:,:,iblk) = uvel(:,:,iblk) - uprev_k(:,:,iblk)
1107 0 : fpresy(:,:,iblk) = vvel(:,:,iblk) - vprev_k(:,:,iblk)
1108 : enddo
1109 : !$OMP END PARALLEL DO
1110 : prog_norm = global_norm(nx_block, ny_block, &
1111 : icellU , & ! LCOV_EXCL_LINE
1112 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1113 0 : fpresx , fpresy )
1114 0 : if (my_task == master_task .and. monitor_nonlin) then
1115 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_nonlin: iter_nonlin= ", it_nl, &
1116 0 : " progress_res_L2norm= ", prog_norm
1117 : endif
1118 :
1119 : enddo ! nonlinear iteration loop
1120 :
1121 0 : end subroutine anderson_solver
1122 :
1123 : !=======================================================================
1124 :
1125 : ! Computes the viscosities and dPr/dx, dPr/dy
1126 :
1127 0 : subroutine calc_zeta_dPr (nx_block, ny_block, &
1128 : icellT , & ! LCOV_EXCL_LINE
1129 : indxTi , indxTj , & ! LCOV_EXCL_LINE
1130 : uvel , vvel , & ! LCOV_EXCL_LINE
1131 : dxT , dyT , & ! LCOV_EXCL_LINE
1132 : dxhy , dyhx , & ! LCOV_EXCL_LINE
1133 : cxp , cyp , & ! LCOV_EXCL_LINE
1134 : cxm , cym , & ! LCOV_EXCL_LINE
1135 : DminTarea,strength, & ! LCOV_EXCL_LINE
1136 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
1137 0 : rep_prs , stPr)
1138 :
1139 : use ice_dyn_shared, only: strain_rates, visc_replpress, &
1140 : capping
1141 :
1142 : integer (kind=int_kind), intent(in) :: &
1143 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1144 : icellT ! no. of cells where iceTmask = .true.
1145 :
1146 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1147 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1148 : indxTj ! compressed index in j-direction
1149 :
1150 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1151 : strength , & ! ice strength (N/m) ! LCOV_EXCL_LINE
1152 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1153 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1154 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1155 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1156 : dxhy , & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
1157 : dyhx , & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
1158 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1159 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1160 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1161 : cxm , & ! 0.5*HTN - 1.5*HTS ! LCOV_EXCL_LINE
1162 : DminTarea ! deltaminVP*tarea
1163 :
1164 : real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(out) :: &
1165 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
1166 : etax2 , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
1167 : rep_prs ! replacement pressure
1168 :
1169 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: &
1170 : stPr ! stress combinations from replacement pressure
1171 :
1172 : ! local variables
1173 :
1174 : integer (kind=int_kind) :: &
1175 : i, j, ij
1176 :
1177 : real (kind=dbl_kind) :: &
1178 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1179 : tensionne, tensionnw, tensionse, tensionsw , & ! tension ! LCOV_EXCL_LINE
1180 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1181 : Deltane, Deltanw, Deltase, Deltasw , & ! Delt ! LCOV_EXCL_LINE
1182 : ssigpn, ssigps, ssigpe, ssigpw, ssigp1, ssigp2, & ! LCOV_EXCL_LINE
1183 : csigpne, csigpnw, csigpsw, csigpse , & ! LCOV_EXCL_LINE
1184 : stressp_1, stressp_2, stressp_3, stressp_4 , & ! LCOV_EXCL_LINE
1185 0 : strp_tmp
1186 :
1187 : character(len=*), parameter :: subname = '(calc_zeta_dPr)'
1188 :
1189 : ! Initialize stPr, zetax2 and etax2 to zero
1190 : ! (for cells where iceTmask is false)
1191 0 : stPr = c0
1192 0 : zetax2 = c0
1193 0 : etax2 = c0
1194 :
1195 0 : do ij = 1, icellT
1196 0 : i = indxTi(ij)
1197 0 : j = indxTj(ij)
1198 :
1199 : !-----------------------------------------------------------------
1200 : ! strain rates
1201 : ! NOTE these are actually strain rates * area (m^2/s)
1202 : !-----------------------------------------------------------------
1203 : call strain_rates (nx_block , ny_block , &
1204 : i , j , & ! LCOV_EXCL_LINE
1205 : uvel , vvel , & ! LCOV_EXCL_LINE
1206 : dxT , dyT , & ! LCOV_EXCL_LINE
1207 : cxp , cyp , & ! LCOV_EXCL_LINE
1208 : cxm , cym , & ! LCOV_EXCL_LINE
1209 : divune , divunw , & ! LCOV_EXCL_LINE
1210 : divuse , divusw , & ! LCOV_EXCL_LINE
1211 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1212 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1213 : shearne , shearnw , & ! LCOV_EXCL_LINE
1214 : shearse , shearsw , & ! LCOV_EXCL_LINE
1215 : Deltane , Deltanw , & ! LCOV_EXCL_LINE
1216 0 : Deltase , Deltasw)
1217 :
1218 : !-----------------------------------------------------------------
1219 : ! viscosities and replacement pressure
1220 : !-----------------------------------------------------------------
1221 :
1222 0 : call visc_replpress (strength(i,j) , DminTarea(i,j) , &
1223 : Deltane , zetax2 (i,j,1), & ! LCOV_EXCL_LINE
1224 : etax2 (i,j,1), rep_prs (i,j,1), & ! LCOV_EXCL_LINE
1225 0 : capping)
1226 0 : call visc_replpress (strength(i,j) , DminTarea(i,j) , &
1227 : Deltanw , zetax2 (i,j,2), & ! LCOV_EXCL_LINE
1228 : etax2 (i,j,2), rep_prs (i,j,2), & ! LCOV_EXCL_LINE
1229 0 : capping)
1230 0 : call visc_replpress (strength(i,j) , DminTarea(i,j) , &
1231 : Deltasw , zetax2 (i,j,3), & ! LCOV_EXCL_LINE
1232 : etax2 (i,j,3), rep_prs (i,j,3), & ! LCOV_EXCL_LINE
1233 0 : capping)
1234 0 : call visc_replpress (strength(i,j) , DminTarea(i,j) , &
1235 : Deltase , zetax2 (i,j,4), & ! LCOV_EXCL_LINE
1236 : etax2 (i,j,4), rep_prs (i,j,4), & ! LCOV_EXCL_LINE
1237 0 : capping)
1238 :
1239 : !-----------------------------------------------------------------
1240 : ! the stresses ! kg/s^2
1241 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
1242 : !-----------------------------------------------------------------
1243 :
1244 0 : stressp_1 = -rep_prs(i,j,1)
1245 0 : stressp_2 = -rep_prs(i,j,2)
1246 0 : stressp_3 = -rep_prs(i,j,3)
1247 0 : stressp_4 = -rep_prs(i,j,4)
1248 :
1249 : !-----------------------------------------------------------------
1250 : ! combinations of the Pr related stresses for the momentum equation ! kg/s^2
1251 : !-----------------------------------------------------------------
1252 :
1253 0 : ssigpn = stressp_1 + stressp_2
1254 0 : ssigps = stressp_3 + stressp_4
1255 0 : ssigpe = stressp_1 + stressp_4
1256 0 : ssigpw = stressp_2 + stressp_3
1257 0 : ssigp1 =(stressp_1 + stressp_3)*p055
1258 0 : ssigp2 =(stressp_2 + stressp_4)*p055
1259 :
1260 0 : csigpne = p111*stressp_1 + ssigp2 + p027*stressp_3
1261 0 : csigpnw = p111*stressp_2 + ssigp1 + p027*stressp_4
1262 0 : csigpsw = p111*stressp_3 + ssigp2 + p027*stressp_1
1263 0 : csigpse = p111*stressp_4 + ssigp1 + p027*stressp_2
1264 :
1265 : !-----------------------------------------------------------------
1266 : ! for dF/dx (u momentum)
1267 : !-----------------------------------------------------------------
1268 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
1269 :
1270 : ! northeast (i,j)
1271 0 : stPr(i,j,1) = -strp_tmp &
1272 0 : + dxhy(i,j)*(-csigpne)
1273 :
1274 : ! northwest (i+1,j)
1275 0 : stPr(i,j,2) = strp_tmp &
1276 0 : + dxhy(i,j)*(-csigpnw)
1277 :
1278 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
1279 :
1280 : ! southeast (i,j+1)
1281 0 : stPr(i,j,3) = -strp_tmp &
1282 0 : + dxhy(i,j)*(-csigpse)
1283 :
1284 : ! southwest (i+1,j+1)
1285 0 : stPr(i,j,4) = strp_tmp &
1286 0 : + dxhy(i,j)*(-csigpsw)
1287 :
1288 : !-----------------------------------------------------------------
1289 : ! for dF/dy (v momentum)
1290 : !-----------------------------------------------------------------
1291 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
1292 :
1293 : ! northeast (i,j)
1294 0 : stPr(i,j,5) = -strp_tmp &
1295 0 : - dyhx(i,j)*(csigpne)
1296 :
1297 : ! southeast (i,j+1)
1298 0 : stPr(i,j,6) = strp_tmp &
1299 0 : - dyhx(i,j)*(csigpse)
1300 :
1301 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
1302 :
1303 : ! northwest (i+1,j)
1304 0 : stPr(i,j,7) = -strp_tmp &
1305 0 : - dyhx(i,j)*(csigpnw)
1306 :
1307 : ! southwest (i+1,j+1)
1308 0 : stPr(i,j,8) = strp_tmp &
1309 0 : - dyhx(i,j)*(csigpsw)
1310 :
1311 : enddo ! ij
1312 :
1313 0 : end subroutine calc_zeta_dPr
1314 :
1315 : !=======================================================================
1316 :
1317 : ! Computes the VP stresses (as diagnostic)
1318 :
1319 : ! Lemieux, J.-F., and Dupont, F. (2020), On the calculation of normalized
1320 : ! viscous-plastic sea ice stresses, Geosci. Model Dev., 13, 1763–1769,
1321 :
1322 0 : subroutine stress_vp (nx_block , ny_block , &
1323 : icellT , & ! LCOV_EXCL_LINE
1324 : indxTi , indxTj , & ! LCOV_EXCL_LINE
1325 : uvel , vvel , & ! LCOV_EXCL_LINE
1326 : dxT , dyT , & ! LCOV_EXCL_LINE
1327 : cxp , cyp , & ! LCOV_EXCL_LINE
1328 : cxm , cym , & ! LCOV_EXCL_LINE
1329 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
1330 : rep_prs , & ! LCOV_EXCL_LINE
1331 : stressp_1 , stressp_2 , & ! LCOV_EXCL_LINE
1332 : stressp_3 , stressp_4 , & ! LCOV_EXCL_LINE
1333 : stressm_1 , stressm_2 , & ! LCOV_EXCL_LINE
1334 : stressm_3 , stressm_4 , & ! LCOV_EXCL_LINE
1335 : stress12_1, stress12_2, & ! LCOV_EXCL_LINE
1336 0 : stress12_3, stress12_4)
1337 :
1338 : use ice_dyn_shared, only: strain_rates
1339 :
1340 : integer (kind=int_kind), intent(in) :: &
1341 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1342 : icellT ! no. of cells where iceTmask = .true.
1343 :
1344 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1345 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1346 : indxTj ! compressed index in j-direction
1347 :
1348 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1349 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1350 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1351 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1352 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1353 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1354 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1355 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1356 : cxm ! 0.5*HTN - 1.5*HTS
1357 :
1358 : real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: &
1359 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
1360 : etax2 , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
1361 : rep_prs
1362 :
1363 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1364 : stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1365 : stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1366 : stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12
1367 :
1368 : ! local variables
1369 :
1370 : integer (kind=int_kind) :: &
1371 : i, j, ij
1372 :
1373 : real (kind=dbl_kind) :: &
1374 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1375 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1376 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1377 0 : Deltane, Deltanw, Deltase, Deltasw ! Delt
1378 :
1379 : character(len=*), parameter :: subname = '(stress_vp)'
1380 :
1381 0 : do ij = 1, icellT
1382 0 : i = indxTi(ij)
1383 0 : j = indxTj(ij)
1384 :
1385 : !-----------------------------------------------------------------
1386 : ! strain rates
1387 : ! NOTE these are actually strain rates * area (m^2/s)
1388 : !-----------------------------------------------------------------
1389 : call strain_rates (nx_block , ny_block , &
1390 : i , j , & ! LCOV_EXCL_LINE
1391 : uvel , vvel , & ! LCOV_EXCL_LINE
1392 : dxT , dyT , & ! LCOV_EXCL_LINE
1393 : cxp , cyp , & ! LCOV_EXCL_LINE
1394 : cxm , cym , & ! LCOV_EXCL_LINE
1395 : divune , divunw , & ! LCOV_EXCL_LINE
1396 : divuse , divusw , & ! LCOV_EXCL_LINE
1397 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1398 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1399 : shearne , shearnw , & ! LCOV_EXCL_LINE
1400 : shearse , shearsw , & ! LCOV_EXCL_LINE
1401 : Deltane , Deltanw , & ! LCOV_EXCL_LINE
1402 0 : Deltase , Deltasw)
1403 :
1404 : !-----------------------------------------------------------------
1405 : ! the stresses ! kg/s^2
1406 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
1407 : !-----------------------------------------------------------------
1408 :
1409 0 : stressp_1(i,j) = zetax2(i,j,1)*divune - rep_prs(i,j,1)
1410 0 : stressp_2(i,j) = zetax2(i,j,2)*divunw - rep_prs(i,j,2)
1411 0 : stressp_3(i,j) = zetax2(i,j,3)*divusw - rep_prs(i,j,3)
1412 0 : stressp_4(i,j) = zetax2(i,j,4)*divuse - rep_prs(i,j,4)
1413 :
1414 0 : stressm_1(i,j) = etax2(i,j,1)*tensionne
1415 0 : stressm_2(i,j) = etax2(i,j,2)*tensionnw
1416 0 : stressm_3(i,j) = etax2(i,j,3)*tensionsw
1417 0 : stressm_4(i,j) = etax2(i,j,4)*tensionse
1418 :
1419 0 : stress12_1(i,j) = etax2(i,j,1)*shearne*p5
1420 0 : stress12_2(i,j) = etax2(i,j,2)*shearnw*p5
1421 0 : stress12_3(i,j) = etax2(i,j,3)*shearsw*p5
1422 0 : stress12_4(i,j) = etax2(i,j,4)*shearse*p5
1423 :
1424 : enddo ! ij
1425 :
1426 0 : end subroutine stress_vp
1427 :
1428 : !=======================================================================
1429 :
1430 : ! Compute vrel and seabed stress coefficients
1431 :
1432 0 : subroutine calc_vrel_Cb (nx_block, ny_block, &
1433 : icellU , Cw , & ! LCOV_EXCL_LINE
1434 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1435 : aiU , TbU , & ! LCOV_EXCL_LINE
1436 : uocn , vocn , & ! LCOV_EXCL_LINE
1437 : uvel , vvel , & ! LCOV_EXCL_LINE
1438 0 : vrel , Cb)
1439 :
1440 : use ice_dyn_shared, only: u0 ! residual velocity for seabed stress (m/s)
1441 :
1442 : integer (kind=int_kind), intent(in) :: &
1443 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1444 : icellU ! total count when iceUmask = .true.
1445 :
1446 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1447 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1448 : indxUj ! compressed index in j-direction
1449 :
1450 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1451 : TbU, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
1452 : aiU , & ! ice fraction on u-grid ! LCOV_EXCL_LINE
1453 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
1454 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
1455 : Cw ! ocean-ice neutral drag coefficient
1456 :
1457 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1458 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1459 : vvel ! y-component of velocity (m/s)
1460 :
1461 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1462 : vrel , & ! coeff for tauw ! LCOV_EXCL_LINE
1463 : Cb ! seabed stress coeff
1464 :
1465 : ! local variables
1466 :
1467 : integer (kind=int_kind) :: &
1468 : i, j, ij
1469 :
1470 : real (kind=dbl_kind) :: &
1471 0 : rhow !
1472 :
1473 : character(len=*), parameter :: subname = '(calc_vrel_Cb)'
1474 :
1475 0 : call icepack_query_parameters(rhow_out=rhow)
1476 0 : call icepack_warnings_flush(nu_diag)
1477 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1478 0 : file=__FILE__, line=__LINE__)
1479 :
1480 0 : do ij = 1, icellU
1481 0 : i = indxUi(ij)
1482 0 : j = indxUj(ij)
1483 :
1484 : ! (magnitude of relative ocean current)*rhow*drag*aice
1485 0 : vrel(i,j) = aiU(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uvel(i,j))**2 + &
1486 0 : (vocn(i,j) - vvel(i,j))**2) ! m/s
1487 :
1488 0 : Cb(i,j) = TbU(i,j) / (sqrt(uvel(i,j)**2 + vvel(i,j)**2) + u0) ! for seabed stress
1489 : enddo ! ij
1490 :
1491 0 : end subroutine calc_vrel_Cb
1492 :
1493 : !=======================================================================
1494 :
1495 : ! Compute seabed stress (diagnostic)
1496 :
1497 0 : subroutine calc_seabed_stress (nx_block, ny_block, &
1498 : icellU , & ! LCOV_EXCL_LINE
1499 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1500 : uvel , vvel , & ! LCOV_EXCL_LINE
1501 : Cb , & ! LCOV_EXCL_LINE
1502 0 : taubxU , taubyU)
1503 :
1504 : integer (kind=int_kind), intent(in) :: &
1505 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1506 : icellU ! total count when iceUmask = .true.
1507 :
1508 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1509 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1510 : indxUj ! compressed index in j-direction
1511 :
1512 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1513 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1514 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1515 : Cb ! seabed stress coefficient
1516 :
1517 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1518 : taubxU , & ! seabed stress, x-direction (N/m^2) ! LCOV_EXCL_LINE
1519 : taubyU ! seabed stress, y-direction (N/m^2)
1520 :
1521 : ! local variables
1522 :
1523 : integer (kind=int_kind) :: &
1524 : i, j, ij
1525 :
1526 : character(len=*), parameter :: subname = '(calc_seabed_stress)'
1527 :
1528 0 : do ij = 1, icellU
1529 0 : i = indxUi(ij)
1530 0 : j = indxUj(ij)
1531 :
1532 0 : taubxU(i,j) = -uvel(i,j)*Cb(i,j)
1533 0 : taubyU(i,j) = -vvel(i,j)*Cb(i,j)
1534 : enddo ! ij
1535 :
1536 0 : end subroutine calc_seabed_stress
1537 :
1538 : !=======================================================================
1539 :
1540 : ! Computes the matrix vector product A(u,v) * (u,v)
1541 : ! Au = A(u,v)_[x] * uvel (x components of A(u,v) * (u,v))
1542 : ! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v))
1543 :
1544 0 : subroutine matvec (nx_block, ny_block, &
1545 : icellU , icellT , & ! LCOV_EXCL_LINE
1546 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1547 : indxTi , indxTj , & ! LCOV_EXCL_LINE
1548 : dxT , dyT , & ! LCOV_EXCL_LINE
1549 : dxhy , dyhx , & ! LCOV_EXCL_LINE
1550 : cxp , cyp , & ! LCOV_EXCL_LINE
1551 : cxm , cym , & ! LCOV_EXCL_LINE
1552 : uvel , vvel , & ! LCOV_EXCL_LINE
1553 : vrel , Cb , & ! LCOV_EXCL_LINE
1554 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
1555 : umassdti, fmU , & ! LCOV_EXCL_LINE
1556 : uarear , & ! LCOV_EXCL_LINE
1557 0 : Au , Av)
1558 :
1559 : use ice_dyn_shared, only: strain_rates
1560 :
1561 : integer (kind=int_kind), intent(in) :: &
1562 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1563 : icellU, & ! total count when iceUmask = .true. ! LCOV_EXCL_LINE
1564 : icellT ! total count when iceTmask = .true.
1565 :
1566 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1567 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1568 : indxUj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
1569 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1570 : indxTj ! compressed index in j-direction
1571 :
1572 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1573 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1574 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1575 : dxhy , & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
1576 : dyhx , & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
1577 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1578 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1579 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1580 : cxm ! 0.5*HTN - 1.5*HTS
1581 :
1582 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1583 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1584 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1585 : vrel , & ! coefficient for tauw ! LCOV_EXCL_LINE
1586 : Cb , & ! coefficient for seabed stress ! LCOV_EXCL_LINE
1587 : umassdti, & ! mass of U-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
1588 : fmU , & ! Coriolis param. * mass in U-cell (kg/s) ! LCOV_EXCL_LINE
1589 : uarear ! 1/uarea
1590 :
1591 : real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: &
1592 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
1593 : etax2 ! etax2 = 2*eta (shear viscosity)
1594 :
1595 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1596 : Au , & ! matvec, Fx = bx - Au (N/m^2) ! LCOV_EXCL_LINE
1597 : Av ! matvec, Fy = by - Av (N/m^2)
1598 :
1599 : ! local variables
1600 :
1601 : integer (kind=int_kind) :: &
1602 : i, j, ij
1603 :
1604 : real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
1605 0 : str
1606 :
1607 : real (kind=dbl_kind) :: &
1608 : ccaimp,ccb , & ! intermediate variables ! LCOV_EXCL_LINE
1609 0 : strintx, strinty ! divergence of the internal stress tensor
1610 :
1611 : real (kind=dbl_kind) :: &
1612 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1613 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1614 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1615 : Deltane, Deltanw, Deltase, Deltasw , & ! Delt ! LCOV_EXCL_LINE
1616 : ssigpn, ssigps, ssigpe, ssigpw , & ! LCOV_EXCL_LINE
1617 : ssigmn, ssigms, ssigme, ssigmw , & ! LCOV_EXCL_LINE
1618 : ssig12n, ssig12s, ssig12e, ssig12w , & ! LCOV_EXCL_LINE
1619 : ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122, & ! LCOV_EXCL_LINE
1620 : csigpne, csigpnw, csigpse, csigpsw , & ! LCOV_EXCL_LINE
1621 : csigmne, csigmnw, csigmse, csigmsw , & ! LCOV_EXCL_LINE
1622 : csig12ne, csig12nw, csig12se, csig12sw , & ! LCOV_EXCL_LINE
1623 : str12ew, str12we, str12ns, str12sn , & ! LCOV_EXCL_LINE
1624 0 : strp_tmp, strm_tmp
1625 :
1626 : real (kind=dbl_kind) :: &
1627 : stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 (without Pr) ! LCOV_EXCL_LINE
1628 : stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1629 0 : stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12
1630 :
1631 : character(len=*), parameter :: subname = '(matvec)'
1632 :
1633 : !-----------------------------------------------------------------
1634 : ! Initialize
1635 : !-----------------------------------------------------------------
1636 :
1637 0 : str(:,:,:) = c0
1638 :
1639 0 : do ij = 1, icellT
1640 0 : i = indxTi(ij)
1641 0 : j = indxTj(ij)
1642 :
1643 : !-----------------------------------------------------------------
1644 : ! strain rates
1645 : ! NOTE these are actually strain rates * area (m^2/s)
1646 : !-----------------------------------------------------------------
1647 : call strain_rates (nx_block , ny_block , &
1648 : i , j , & ! LCOV_EXCL_LINE
1649 : uvel , vvel , & ! LCOV_EXCL_LINE
1650 : dxT , dyT , & ! LCOV_EXCL_LINE
1651 : cxp , cyp , & ! LCOV_EXCL_LINE
1652 : cxm , cym , & ! LCOV_EXCL_LINE
1653 : divune , divunw , & ! LCOV_EXCL_LINE
1654 : divuse , divusw , & ! LCOV_EXCL_LINE
1655 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1656 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1657 : shearne , shearnw , & ! LCOV_EXCL_LINE
1658 : shearse , shearsw , & ! LCOV_EXCL_LINE
1659 : Deltane , Deltanw , & ! LCOV_EXCL_LINE
1660 0 : Deltase , Deltasw)
1661 :
1662 : !-----------------------------------------------------------------
1663 : ! the stresses ! kg/s^2
1664 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
1665 : ! NOTE: commented part of stressp is from the replacement pressure Pr
1666 : !-----------------------------------------------------------------
1667 :
1668 0 : stressp_1 = zetax2(i,j,1)*divune! - Deltane*(c1-Ktens))
1669 0 : stressp_2 = zetax2(i,j,2)*divunw! - Deltanw*(c1-Ktens))
1670 0 : stressp_3 = zetax2(i,j,3)*divusw! - Deltasw*(c1-Ktens))
1671 0 : stressp_4 = zetax2(i,j,4)*divuse! - Deltase*(c1-Ktens))
1672 :
1673 0 : stressm_1 = etax2(i,j,1)*tensionne
1674 0 : stressm_2 = etax2(i,j,2)*tensionnw
1675 0 : stressm_3 = etax2(i,j,3)*tensionsw
1676 0 : stressm_4 = etax2(i,j,4)*tensionse
1677 :
1678 0 : stress12_1 = etax2(i,j,1)*shearne*p5
1679 0 : stress12_2 = etax2(i,j,2)*shearnw*p5
1680 0 : stress12_3 = etax2(i,j,3)*shearsw*p5
1681 0 : stress12_4 = etax2(i,j,4)*shearse*p5
1682 :
1683 : !-----------------------------------------------------------------
1684 : ! combinations of the stresses for the momentum equation ! kg/s^2
1685 : !-----------------------------------------------------------------
1686 :
1687 0 : ssigpn = stressp_1 + stressp_2
1688 0 : ssigps = stressp_3 + stressp_4
1689 0 : ssigpe = stressp_1 + stressp_4
1690 0 : ssigpw = stressp_2 + stressp_3
1691 0 : ssigp1 =(stressp_1 + stressp_3)*p055
1692 0 : ssigp2 =(stressp_2 + stressp_4)*p055
1693 :
1694 0 : ssigmn = stressm_1 + stressm_2
1695 0 : ssigms = stressm_3 + stressm_4
1696 0 : ssigme = stressm_1 + stressm_4
1697 0 : ssigmw = stressm_2 + stressm_3
1698 0 : ssigm1 =(stressm_1 + stressm_3)*p055
1699 0 : ssigm2 =(stressm_2 + stressm_4)*p055
1700 :
1701 0 : ssig12n = stress12_1 + stress12_2
1702 0 : ssig12s = stress12_3 + stress12_4
1703 0 : ssig12e = stress12_1 + stress12_4
1704 0 : ssig12w = stress12_2 + stress12_3
1705 0 : ssig121 =(stress12_1 + stress12_3)*p111
1706 0 : ssig122 =(stress12_2 + stress12_4)*p111
1707 :
1708 0 : csigpne = p111*stressp_1 + ssigp2 + p027*stressp_3
1709 0 : csigpnw = p111*stressp_2 + ssigp1 + p027*stressp_4
1710 0 : csigpsw = p111*stressp_3 + ssigp2 + p027*stressp_1
1711 0 : csigpse = p111*stressp_4 + ssigp1 + p027*stressp_2
1712 :
1713 0 : csigmne = p111*stressm_1 + ssigm2 + p027*stressm_3
1714 0 : csigmnw = p111*stressm_2 + ssigm1 + p027*stressm_4
1715 0 : csigmsw = p111*stressm_3 + ssigm2 + p027*stressm_1
1716 0 : csigmse = p111*stressm_4 + ssigm1 + p027*stressm_2
1717 :
1718 : csig12ne = p222*stress12_1 + ssig122 &
1719 0 : + p055*stress12_3
1720 : csig12nw = p222*stress12_2 + ssig121 &
1721 0 : + p055*stress12_4
1722 : csig12sw = p222*stress12_3 + ssig122 &
1723 0 : + p055*stress12_1
1724 : csig12se = p222*stress12_4 + ssig121 &
1725 0 : + p055*stress12_2
1726 :
1727 0 : str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
1728 0 : str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
1729 0 : str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
1730 0 : str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)
1731 :
1732 : !-----------------------------------------------------------------
1733 : ! for dF/dx (u momentum)
1734 : !-----------------------------------------------------------------
1735 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
1736 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)
1737 :
1738 : ! northeast (i,j)
1739 0 : str(i,j,1) = -strp_tmp - strm_tmp - str12ew &
1740 0 : + dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne
1741 :
1742 : ! northwest (i+1,j)
1743 0 : str(i,j,2) = strp_tmp + strm_tmp - str12we &
1744 0 : + dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw
1745 :
1746 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
1747 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)
1748 :
1749 : ! southeast (i,j+1)
1750 0 : str(i,j,3) = -strp_tmp - strm_tmp + str12ew &
1751 0 : + dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se
1752 :
1753 : ! southwest (i+1,j+1)
1754 0 : str(i,j,4) = strp_tmp + strm_tmp + str12we &
1755 0 : + dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw
1756 :
1757 : !-----------------------------------------------------------------
1758 : ! for dF/dy (v momentum)
1759 : !-----------------------------------------------------------------
1760 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
1761 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)
1762 :
1763 : ! northeast (i,j)
1764 0 : str(i,j,5) = -strp_tmp + strm_tmp - str12ns &
1765 0 : - dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne
1766 :
1767 : ! southeast (i,j+1)
1768 0 : str(i,j,6) = strp_tmp - strm_tmp - str12sn &
1769 0 : - dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se
1770 :
1771 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
1772 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)
1773 :
1774 : ! northwest (i+1,j)
1775 0 : str(i,j,7) = -strp_tmp + strm_tmp + str12ns &
1776 0 : - dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw
1777 :
1778 : ! southwest (i+1,j+1)
1779 0 : str(i,j,8) = strp_tmp - strm_tmp + str12sn &
1780 0 : - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw
1781 :
1782 : enddo ! ij - icellT
1783 :
1784 : !-----------------------------------------------------------------
1785 : ! Form Au and Av
1786 : !-----------------------------------------------------------------
1787 :
1788 0 : do ij = 1, icellU
1789 0 : i = indxUi(ij)
1790 0 : j = indxUj(ij)
1791 :
1792 0 : ccaimp = umassdti(i,j) + vrel(i,j) * cosw + Cb(i,j) ! kg/m^2 s
1793 :
1794 0 : ccb = fmU(i,j) + sign(c1,fmU(i,j)) * vrel(i,j) * sinw ! kg/m^2 s
1795 :
1796 : ! divergence of the internal stress tensor
1797 0 : strintx = uarear(i,j)* &
1798 0 : (str(i,j,1) + str(i+1,j,2) + str(i,j+1,3) + str(i+1,j+1,4))
1799 0 : strinty = uarear(i,j)* &
1800 0 : (str(i,j,5) + str(i,j+1,6) + str(i+1,j,7) + str(i+1,j+1,8))
1801 :
1802 0 : Au(i,j) = ccaimp*uvel(i,j) - ccb*vvel(i,j) - strintx
1803 0 : Av(i,j) = ccaimp*vvel(i,j) + ccb*uvel(i,j) - strinty
1804 : enddo ! ij - icellU
1805 :
1806 0 : end subroutine matvec
1807 :
1808 : !=======================================================================
1809 :
1810 : ! Compute the constant component of b(u,v) i.e. the part of b(u,v) that
1811 : ! does not depend on (u,v) and thus do not change during the nonlinear iteration
1812 :
1813 0 : subroutine calc_bfix (nx_block , ny_block , &
1814 : icellU , & ! LCOV_EXCL_LINE
1815 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1816 : umassdti , & ! LCOV_EXCL_LINE
1817 : forcexU , forceyU , & ! LCOV_EXCL_LINE
1818 : uvel_init, vvel_init, & ! LCOV_EXCL_LINE
1819 0 : bxfix , byfix)
1820 :
1821 : integer (kind=int_kind), intent(in) :: &
1822 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1823 : icellU ! no. of cells where iceUmask = .true.
1824 :
1825 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1826 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1827 : indxUj ! compressed index in j-direction
1828 :
1829 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1830 : uvel_init,& ! x-component of velocity (m/s), beginning of time step ! LCOV_EXCL_LINE
1831 : vvel_init,& ! y-component of velocity (m/s), beginning of time step ! LCOV_EXCL_LINE
1832 : umassdti, & ! mass of U-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
1833 : forcexU , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
1834 : forceyU ! work array: combined atm stress and ocn tilt, y
1835 :
1836 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1837 : bxfix , & ! bx = taux + bxfix ! LCOV_EXCL_LINE
1838 : byfix ! by = tauy + byfix
1839 :
1840 : ! local variables
1841 :
1842 : integer (kind=int_kind) :: &
1843 : i, j, ij
1844 :
1845 : character(len=*), parameter :: subname = '(calc_bfix)'
1846 :
1847 0 : do ij = 1, icellU
1848 0 : i = indxUi(ij)
1849 0 : j = indxUj(ij)
1850 :
1851 0 : bxfix(i,j) = umassdti(i,j)*uvel_init(i,j) + forcexU(i,j)
1852 0 : byfix(i,j) = umassdti(i,j)*vvel_init(i,j) + forceyU(i,j)
1853 : enddo
1854 :
1855 0 : end subroutine calc_bfix
1856 :
1857 : !=======================================================================
1858 :
1859 : ! Compute the vector b(u,v), i.e. the part of the nonlinear function F(u,v)
1860 : ! that cannot be written as A(u,v)*(u,v), where A(u,v) is a matrix with entries
1861 : ! depending on (u,v)
1862 :
1863 0 : subroutine calc_bvec (nx_block, ny_block, &
1864 : icellU , & ! LCOV_EXCL_LINE
1865 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1866 : stPr , uarear , & ! LCOV_EXCL_LINE
1867 : waterxU , wateryU , & ! LCOV_EXCL_LINE
1868 : bxfix , byfix , & ! LCOV_EXCL_LINE
1869 : bx , by , & ! LCOV_EXCL_LINE
1870 0 : vrel)
1871 :
1872 : integer (kind=int_kind), intent(in) :: &
1873 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1874 : icellU ! total count when iceUmask = .true.
1875 :
1876 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1877 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1878 : indxUj ! compressed index in j-direction
1879 :
1880 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1881 : uarear , & ! 1/uarea ! LCOV_EXCL_LINE
1882 : waterxU , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
1883 : wateryU , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
1884 : bxfix , & ! bx = taux + bxfix ! LCOV_EXCL_LINE
1885 : byfix , & ! by = tauy + byfix ! LCOV_EXCL_LINE
1886 : vrel ! relative ice-ocean velocity
1887 :
1888 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(in) :: &
1889 : stPr
1890 :
1891 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
1892 : bx , & ! b vector, bx = taux + bxfix (N/m^2) ! LCOV_EXCL_LINE
1893 : by ! b vector, by = tauy + byfix (N/m^2)
1894 :
1895 : ! local variables
1896 :
1897 : integer (kind=int_kind) :: &
1898 : i, j, ij
1899 :
1900 : real (kind=dbl_kind) :: &
1901 : taux, tauy , & ! part of ocean stress term ! LCOV_EXCL_LINE
1902 : strintx, strinty , & ! divergence of the internal stress tensor (only Pr contributions) ! LCOV_EXCL_LINE
1903 0 : rhow !
1904 :
1905 : character(len=*), parameter :: subname = '(calc_bvec)'
1906 :
1907 : !-----------------------------------------------------------------
1908 : ! calc b vector
1909 : !-----------------------------------------------------------------
1910 :
1911 0 : call icepack_query_parameters(rhow_out=rhow)
1912 0 : call icepack_warnings_flush(nu_diag)
1913 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1914 0 : file=__FILE__, line=__LINE__)
1915 :
1916 0 : do ij = 1, icellU
1917 0 : i = indxUi(ij)
1918 0 : j = indxUj(ij)
1919 :
1920 : ! ice/ocean stress
1921 0 : taux = vrel(i,j)*waterxU(i,j) ! NOTE this is not the entire
1922 0 : tauy = vrel(i,j)*wateryU(i,j) ! ocn stress term
1923 :
1924 : ! divergence of the internal stress tensor (only Pr part, i.e. dPr/dx, dPr/dy)
1925 0 : strintx = uarear(i,j)* &
1926 0 : (stPr(i,j,1) + stPr(i+1,j,2) + stPr(i,j+1,3) + stPr(i+1,j+1,4))
1927 0 : strinty = uarear(i,j)* &
1928 0 : (stPr(i,j,5) + stPr(i,j+1,6) + stPr(i+1,j,7) + stPr(i+1,j+1,8))
1929 :
1930 0 : bx(i,j) = bxfix(i,j) + taux + strintx
1931 0 : by(i,j) = byfix(i,j) + tauy + strinty
1932 : enddo ! ij
1933 :
1934 0 : end subroutine calc_bvec
1935 :
1936 : !=======================================================================
1937 :
1938 : ! Compute the non linear residual F(u,v) = b(u,v) - A(u,v) * (u,v),
1939 : ! with Au, Av precomputed as
1940 : ! Au = A(u,v)_[x] * uvel (x components of A(u,v) * (u,v))
1941 : ! Av = A(u,v)_[y] * vvel (y components of A(u,v) * (u,v))
1942 :
1943 0 : subroutine residual_vec (nx_block , ny_block, &
1944 : icellU , & ! LCOV_EXCL_LINE
1945 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1946 : bx , by , & ! LCOV_EXCL_LINE
1947 : Au , Av , & ! LCOV_EXCL_LINE
1948 0 : Fx , Fy )
1949 :
1950 : integer (kind=int_kind), intent(in) :: &
1951 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1952 : icellU ! total count when iceUmask = .true.
1953 :
1954 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1955 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1956 : indxUj ! compressed index in j-direction
1957 :
1958 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1959 : bx , & ! b vector, bx = taux + bxfix (N/m^2) ! LCOV_EXCL_LINE
1960 : by , & ! b vector, by = tauy + byfix (N/m^2) ! LCOV_EXCL_LINE
1961 : Au , & ! matvec, Fx = bx - Au (N/m^2) ! LCOV_EXCL_LINE
1962 : Av ! matvec, Fy = by - Av (N/m^2)
1963 :
1964 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1965 : Fx , & ! x residual vector, Fx = bx - Au (N/m^2) ! LCOV_EXCL_LINE
1966 : Fy ! y residual vector, Fy = by - Av (N/m^2)
1967 :
1968 : ! local variables
1969 :
1970 : integer (kind=int_kind) :: &
1971 : i, j, ij
1972 :
1973 : character(len=*), parameter :: subname = '(residual_vec)'
1974 :
1975 : !-----------------------------------------------------------------
1976 : ! compute residual
1977 : !-----------------------------------------------------------------
1978 :
1979 0 : do ij = 1, icellU
1980 0 : i = indxUi(ij)
1981 0 : j = indxUj(ij)
1982 :
1983 0 : Fx(i,j) = bx(i,j) - Au(i,j)
1984 0 : Fy(i,j) = by(i,j) - Av(i,j)
1985 : enddo ! ij
1986 :
1987 0 : end subroutine residual_vec
1988 :
1989 : !=======================================================================
1990 :
1991 : ! Form the diagonal of the matrix A(u,v) (first part of the computation)
1992 : ! Part 1: compute the contributions to the diagonal from the rheology term
1993 :
1994 0 : subroutine formDiag_step1 (nx_block, ny_block, &
1995 : icellU , & ! LCOV_EXCL_LINE
1996 : indxUi , indxUj , & ! LCOV_EXCL_LINE
1997 : dxT , dyT , & ! LCOV_EXCL_LINE
1998 : dxhy , dyhx , & ! LCOV_EXCL_LINE
1999 : cxp , cyp , & ! LCOV_EXCL_LINE
2000 : cxm , cym , & ! LCOV_EXCL_LINE
2001 : zetax2 , etax2 , & ! LCOV_EXCL_LINE
2002 0 : Drheo)
2003 :
2004 : integer (kind=int_kind), intent(in) :: &
2005 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2006 : icellU ! no. of cells where iceUmask = .true.
2007 :
2008 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2009 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2010 : indxUj ! compressed index in j-direction
2011 :
2012 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2013 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
2014 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
2015 : dxhy , & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
2016 : dyhx , & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
2017 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
2018 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
2019 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
2020 : cxm ! 0.5*HTN - 1.5*HTS
2021 :
2022 : real (kind=dbl_kind), dimension(nx_block,ny_block,4), intent(in) :: &
2023 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
2024 : etax2 ! etax2 = 2*eta (shear viscosity)
2025 :
2026 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: &
2027 : Drheo ! intermediate value for diagonal components of matrix A associated
2028 : ! with rheology term
2029 :
2030 : ! local variables
2031 :
2032 : integer (kind=int_kind) :: &
2033 : i, j, ij, iu, ju, di, dj, cc
2034 :
2035 : real (kind=dbl_kind) :: &
2036 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
2037 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
2038 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
2039 : uij,ui1j,uij1,ui1j1,vij,vi1j,vij1,vi1j1 , & ! == c0 or c1 ! LCOV_EXCL_LINE
2040 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
2041 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
2042 : stress12_1,stress12_2,stress12_3,stress12_4,& ! LCOV_EXCL_LINE
2043 : ssigpn, ssigps, ssigpe, ssigpw , & ! LCOV_EXCL_LINE
2044 : ssigmn, ssigms, ssigme, ssigmw , & ! LCOV_EXCL_LINE
2045 : ssig12n, ssig12s, ssig12e, ssig12w , & ! LCOV_EXCL_LINE
2046 : ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122, & ! LCOV_EXCL_LINE
2047 : csigpne, csigpnw, csigpse, csigpsw , & ! LCOV_EXCL_LINE
2048 : csigmne, csigmnw, csigmse, csigmsw , & ! LCOV_EXCL_LINE
2049 : csig12ne, csig12nw, csig12se, csig12sw , & ! LCOV_EXCL_LINE
2050 : str12ew, str12we, str12ns, str12sn , & ! LCOV_EXCL_LINE
2051 0 : strp_tmp, strm_tmp
2052 :
2053 : character(len=*), parameter :: subname = '(formDiag_step1)'
2054 :
2055 : !-----------------------------------------------------------------
2056 : ! Initialize
2057 : !-----------------------------------------------------------------
2058 :
2059 0 : Drheo(:,:,:) = c0
2060 :
2061 : ! Be careful: Drheo contains 4 terms for u and 4 terms for v.
2062 : ! These 8 terms come from the surrounding T cells but are all
2063 : ! refrerenced to the i,j (u point) :
2064 :
2065 : ! Drheo(i,j,1) corresponds to str(i,j,1)
2066 : ! Drheo(i,j,2) corresponds to str(i+1,j,2)
2067 : ! Drheo(i,j,3) corresponds to str(i,j+1,3)
2068 : ! Drheo(i,j,4) corresponds to str(i+1,j+1,4))
2069 : ! Drheo(i,j,5) corresponds to str(i,j,5)
2070 : ! Drheo(i,j,6) corresponds to str(i,j+1,6)
2071 : ! Drheo(i,j,7) corresponds to str(i+1,j,7)
2072 : ! Drheo(i,j,8) corresponds to str(i+1,j+1,8))
2073 :
2074 0 : do cc = 1, 8 ! 4 for u and 4 for v
2075 :
2076 0 : if (cc == 1) then ! u comp, T cell i,j
2077 0 : uij = c1
2078 0 : ui1j = c0
2079 0 : uij1 = c0
2080 0 : ui1j1 = c0
2081 0 : vij = c0
2082 0 : vi1j = c0
2083 0 : vij1 = c0
2084 0 : vi1j1 = c0
2085 0 : di = 0
2086 0 : dj = 0
2087 0 : elseif (cc == 2) then ! u comp, T cell i+1,j
2088 0 : uij = c0
2089 0 : ui1j = c1
2090 0 : uij1 = c0
2091 0 : ui1j1 = c0
2092 0 : vij = c0
2093 0 : vi1j = c0
2094 0 : vij1 = c0
2095 0 : vi1j1 = c0
2096 0 : di = 1
2097 0 : dj = 0
2098 0 : elseif (cc == 3) then ! u comp, T cell i,j+1
2099 0 : uij = c0
2100 0 : ui1j = c0
2101 0 : uij1 = c1
2102 0 : ui1j1 = c0
2103 0 : vij = c0
2104 0 : vi1j = c0
2105 0 : vij1 = c0
2106 0 : vi1j1 = c0
2107 0 : di = 0
2108 0 : dj = 1
2109 0 : elseif (cc == 4) then ! u comp, T cell i+1,j+1
2110 0 : uij = c0
2111 0 : ui1j = c0
2112 0 : uij1 = c0
2113 0 : ui1j1 = c1
2114 0 : vij = c0
2115 0 : vi1j = c0
2116 0 : vij1 = c0
2117 0 : vi1j1 = c0
2118 0 : di = 1
2119 0 : dj = 1
2120 0 : elseif (cc == 5) then ! v comp, T cell i,j
2121 0 : uij = c0
2122 0 : ui1j = c0
2123 0 : uij1 = c0
2124 0 : ui1j1 = c0
2125 0 : vij = c1
2126 0 : vi1j = c0
2127 0 : vij1 = c0
2128 0 : vi1j1 = c0
2129 0 : di = 0
2130 0 : dj = 0
2131 0 : elseif (cc == 6) then ! v comp, T cell i,j+1
2132 0 : uij = c0
2133 0 : ui1j = c0
2134 0 : uij1 = c0
2135 0 : ui1j1 = c0
2136 0 : vij = c0
2137 0 : vi1j = c0
2138 0 : vij1 = c1
2139 0 : vi1j1 = c0
2140 0 : di = 0
2141 0 : dj = 1
2142 0 : elseif (cc == 7) then ! v comp, T cell i+1,j
2143 0 : uij = c0
2144 0 : ui1j = c0
2145 0 : uij1 = c0
2146 0 : ui1j1 = c0
2147 0 : vij = c0
2148 0 : vi1j = c1
2149 0 : vij1 = c0
2150 0 : vi1j1 = c0
2151 0 : di = 1
2152 0 : dj = 0
2153 0 : elseif (cc == 8) then ! v comp, T cell i+1,j+1
2154 0 : uij = c0
2155 0 : ui1j = c0
2156 0 : uij1 = c0
2157 0 : ui1j1 = c0
2158 0 : vij = c0
2159 0 : vi1j = c0
2160 0 : vij1 = c0
2161 0 : vi1j1 = c1
2162 0 : di = 1
2163 0 : dj = 1
2164 : endif
2165 :
2166 0 : do ij = 1, icellU
2167 :
2168 0 : iu = indxUi(ij)
2169 0 : ju = indxUj(ij)
2170 0 : i = iu + di
2171 0 : j = ju + dj
2172 :
2173 : !-----------------------------------------------------------------
2174 : ! strain rates
2175 : ! NOTE these are actually strain rates * area (m^2/s)
2176 : !-----------------------------------------------------------------
2177 : ! divergence = e_11 + e_22
2178 0 : divune = cyp(i,j)*uij - dyT(i,j)*ui1j &
2179 0 : + cxp(i,j)*vij - dxT(i,j)*vij1
2180 0 : divunw = cym(i,j)*ui1j + dyT(i,j)*uij &
2181 0 : + cxp(i,j)*vi1j - dxT(i,j)*vi1j1
2182 0 : divusw = cym(i,j)*ui1j1 + dyT(i,j)*uij1 &
2183 0 : + cxm(i,j)*vi1j1 + dxT(i,j)*vi1j
2184 0 : divuse = cyp(i,j)*uij1 - dyT(i,j)*ui1j1 &
2185 0 : + cxm(i,j)*vij1 + dxT(i,j)*vij
2186 :
2187 : ! tension strain rate = e_11 - e_22
2188 0 : tensionne = -cym(i,j)*uij - dyT(i,j)*ui1j &
2189 0 : + cxm(i,j)*vij + dxT(i,j)*vij1
2190 0 : tensionnw = -cyp(i,j)*ui1j + dyT(i,j)*uij &
2191 0 : + cxm(i,j)*vi1j + dxT(i,j)*vi1j1
2192 0 : tensionsw = -cyp(i,j)*ui1j1 + dyT(i,j)*uij1 &
2193 0 : + cxp(i,j)*vi1j1 - dxT(i,j)*vi1j
2194 0 : tensionse = -cym(i,j)*uij1 - dyT(i,j)*ui1j1 &
2195 0 : + cxp(i,j)*vij1 - dxT(i,j)*vij
2196 :
2197 : ! shearing strain rate = 2*e_12
2198 0 : shearne = -cym(i,j)*vij - dyT(i,j)*vi1j &
2199 0 : - cxm(i,j)*uij - dxT(i,j)*uij1
2200 0 : shearnw = -cyp(i,j)*vi1j + dyT(i,j)*vij &
2201 0 : - cxm(i,j)*ui1j - dxT(i,j)*ui1j1
2202 0 : shearsw = -cyp(i,j)*vi1j1 + dyT(i,j)*vij1 &
2203 0 : - cxp(i,j)*ui1j1 + dxT(i,j)*ui1j
2204 0 : shearse = -cym(i,j)*vij1 - dyT(i,j)*vi1j1 &
2205 0 : - cxp(i,j)*uij1 + dxT(i,j)*uij
2206 :
2207 : !-----------------------------------------------------------------
2208 : ! the stresses ! kg/s^2
2209 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
2210 : !-----------------------------------------------------------------
2211 :
2212 0 : stressp_1 = zetax2(i,j,1)*divune
2213 0 : stressp_2 = zetax2(i,j,2)*divunw
2214 0 : stressp_3 = zetax2(i,j,3)*divusw
2215 0 : stressp_4 = zetax2(i,j,4)*divuse
2216 :
2217 0 : stressm_1 = etax2(i,j,1)*tensionne
2218 0 : stressm_2 = etax2(i,j,2)*tensionnw
2219 0 : stressm_3 = etax2(i,j,3)*tensionsw
2220 0 : stressm_4 = etax2(i,j,4)*tensionse
2221 :
2222 0 : stress12_1 = etax2(i,j,1)*shearne*p5
2223 0 : stress12_2 = etax2(i,j,2)*shearnw*p5
2224 0 : stress12_3 = etax2(i,j,3)*shearsw*p5
2225 0 : stress12_4 = etax2(i,j,4)*shearse*p5
2226 :
2227 : !-----------------------------------------------------------------
2228 : ! combinations of the stresses for the momentum equation ! kg/s^2
2229 : !-----------------------------------------------------------------
2230 :
2231 0 : ssigpn = stressp_1 + stressp_2
2232 0 : ssigps = stressp_3 + stressp_4
2233 0 : ssigpe = stressp_1 + stressp_4
2234 0 : ssigpw = stressp_2 + stressp_3
2235 0 : ssigp1 =(stressp_1 + stressp_3)*p055
2236 0 : ssigp2 =(stressp_2 + stressp_4)*p055
2237 :
2238 0 : ssigmn = stressm_1 + stressm_2
2239 0 : ssigms = stressm_3 + stressm_4
2240 0 : ssigme = stressm_1 + stressm_4
2241 0 : ssigmw = stressm_2 + stressm_3
2242 0 : ssigm1 =(stressm_1 + stressm_3)*p055
2243 0 : ssigm2 =(stressm_2 + stressm_4)*p055
2244 :
2245 0 : ssig12n = stress12_1 + stress12_2
2246 0 : ssig12s = stress12_3 + stress12_4
2247 0 : ssig12e = stress12_1 + stress12_4
2248 0 : ssig12w = stress12_2 + stress12_3
2249 0 : ssig121 =(stress12_1 + stress12_3)*p111
2250 0 : ssig122 =(stress12_2 + stress12_4)*p111
2251 :
2252 0 : csigpne = p111*stressp_1 + ssigp2 + p027*stressp_3
2253 0 : csigpnw = p111*stressp_2 + ssigp1 + p027*stressp_4
2254 0 : csigpsw = p111*stressp_3 + ssigp2 + p027*stressp_1
2255 0 : csigpse = p111*stressp_4 + ssigp1 + p027*stressp_2
2256 :
2257 0 : csigmne = p111*stressm_1 + ssigm2 + p027*stressm_3
2258 0 : csigmnw = p111*stressm_2 + ssigm1 + p027*stressm_4
2259 0 : csigmsw = p111*stressm_3 + ssigm2 + p027*stressm_1
2260 0 : csigmse = p111*stressm_4 + ssigm1 + p027*stressm_2
2261 :
2262 : csig12ne = p222*stress12_1 + ssig122 &
2263 0 : + p055*stress12_3
2264 : csig12nw = p222*stress12_2 + ssig121 &
2265 0 : + p055*stress12_4
2266 : csig12sw = p222*stress12_3 + ssig122 &
2267 0 : + p055*stress12_1
2268 : csig12se = p222*stress12_4 + ssig121 &
2269 0 : + p055*stress12_2
2270 :
2271 0 : str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
2272 0 : str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
2273 0 : str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
2274 0 : str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)
2275 :
2276 : !-----------------------------------------------------------------
2277 : ! for dF/dx (u momentum)
2278 : !-----------------------------------------------------------------
2279 :
2280 0 : if (cc == 1) then ! T cell i,j
2281 :
2282 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
2283 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)
2284 :
2285 : ! northeast (i,j)
2286 0 : Drheo(iu,ju,1) = -strp_tmp - strm_tmp - str12ew &
2287 0 : + dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne
2288 :
2289 0 : elseif (cc == 2) then ! T cell i+1,j
2290 :
2291 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
2292 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)
2293 :
2294 : ! northwest (i+1,j)
2295 0 : Drheo(iu,ju,2) = strp_tmp + strm_tmp - str12we &
2296 0 : + dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw
2297 :
2298 0 : elseif (cc == 3) then ! T cell i,j+1
2299 :
2300 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
2301 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)
2302 :
2303 : ! southeast (i,j+1)
2304 0 : Drheo(iu,ju,3) = -strp_tmp - strm_tmp + str12ew &
2305 0 : + dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se
2306 :
2307 0 : elseif (cc == 4) then ! T cell i+1,j+1
2308 :
2309 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
2310 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)
2311 :
2312 : ! southwest (i+1,j+1)
2313 0 : Drheo(iu,ju,4) = strp_tmp + strm_tmp + str12we &
2314 0 : + dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw
2315 :
2316 : !-----------------------------------------------------------------
2317 : ! for dF/dy (v momentum)
2318 : !-----------------------------------------------------------------
2319 :
2320 0 : elseif (cc == 5) then ! T cell i,j
2321 :
2322 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
2323 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)
2324 :
2325 : ! northeast (i,j)
2326 0 : Drheo(iu,ju,5) = -strp_tmp + strm_tmp - str12ns &
2327 0 : - dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne
2328 :
2329 0 : elseif (cc == 6) then ! T cell i,j+1
2330 :
2331 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
2332 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)
2333 :
2334 : ! southeast (i,j+1)
2335 0 : Drheo(iu,ju,6) = strp_tmp - strm_tmp - str12sn &
2336 0 : - dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se
2337 :
2338 0 : elseif (cc == 7) then ! T cell i,j+1
2339 :
2340 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
2341 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)
2342 :
2343 : ! northwest (i+1,j)
2344 0 : Drheo(iu,ju,7) = -strp_tmp + strm_tmp + str12ns &
2345 0 : - dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw
2346 :
2347 0 : elseif (cc == 8) then ! T cell i+1,j+1
2348 :
2349 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
2350 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)
2351 :
2352 : ! southwest (i+1,j+1)
2353 0 : Drheo(iu,ju,8) = strp_tmp - strm_tmp + str12sn &
2354 0 : - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw
2355 :
2356 : endif
2357 :
2358 : enddo ! ij
2359 :
2360 : enddo ! cc
2361 :
2362 0 : end subroutine formDiag_step1
2363 :
2364 : !=======================================================================
2365 :
2366 : ! Form the diagonal of the matrix A(u,v) (second part of the computation)
2367 : ! Part 2: compute diagonal
2368 :
2369 0 : subroutine formDiag_step2 (nx_block, ny_block, &
2370 : icellU , & ! LCOV_EXCL_LINE
2371 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2372 : Drheo , vrel , & ! LCOV_EXCL_LINE
2373 : umassdti, & ! LCOV_EXCL_LINE
2374 : uarear , Cb , & ! LCOV_EXCL_LINE
2375 0 : diagx , diagy)
2376 :
2377 : integer (kind=int_kind), intent(in) :: &
2378 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2379 : icellU ! total count when iceUmask = .true.
2380 :
2381 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2382 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2383 : indxUj ! compressed index in j-direction
2384 :
2385 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2386 : vrel, & ! coefficient for tauw ! LCOV_EXCL_LINE
2387 : Cb, & ! coefficient for seabed stress ! LCOV_EXCL_LINE
2388 : umassdti, & ! mass of U-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
2389 : uarear ! 1/uarea
2390 :
2391 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(in) :: &
2392 : Drheo
2393 :
2394 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
2395 : diagx , & ! Diagonal (x component) of the matrix A ! LCOV_EXCL_LINE
2396 : diagy ! Diagonal (y component) of the matrix A
2397 :
2398 : ! local variables
2399 :
2400 : integer (kind=int_kind) :: &
2401 : i, j, ij
2402 :
2403 : real (kind=dbl_kind) :: &
2404 : ccaimp , & ! intermediate variables ! LCOV_EXCL_LINE
2405 0 : strintx, strinty ! diagonal contributions to the divergence
2406 :
2407 : character(len=*), parameter :: subname = '(formDiag_step2)'
2408 :
2409 : !-----------------------------------------------------------------
2410 : ! integrate the momentum equation
2411 : !-----------------------------------------------------------------
2412 :
2413 0 : strintx = c0
2414 0 : strinty = c0
2415 :
2416 : ! Be careful: Drheo contains 4 terms for u and 4 terms for v.
2417 : ! These 8 terms come from the surrounding T cells but are all
2418 : ! refrerenced to the i,j (u point) :
2419 :
2420 : ! Drheo(i,j,1) corresponds to str(i,j,1)
2421 : ! Drheo(i,j,2) corresponds to str(i+1,j,2)
2422 : ! Drheo(i,j,3) corresponds to str(i,j+1,3)
2423 : ! Drheo(i,j,4) corresponds to str(i+1,j+1,4))
2424 : ! Drheo(i,j,5) corresponds to str(i,j,5)
2425 : ! Drheo(i,j,6) corresponds to str(i,j+1,6)
2426 : ! Drheo(i,j,7) corresponds to str(i+1,j,7)
2427 : ! Drheo(i,j,8) corresponds to str(i+1,j+1,8))
2428 :
2429 0 : do ij = 1, icellU
2430 0 : i = indxUi(ij)
2431 0 : j = indxUj(ij)
2432 :
2433 0 : ccaimp = umassdti(i,j) + vrel(i,j) * cosw + Cb(i,j) ! kg/m^2 s
2434 :
2435 0 : strintx = uarear(i,j)* &
2436 0 : (Drheo(i,j,1) + Drheo(i,j,2) + Drheo(i,j,3) + Drheo(i,j,4))
2437 0 : strinty = uarear(i,j)* &
2438 0 : (Drheo(i,j,5) + Drheo(i,j,6) + Drheo(i,j,7) + Drheo(i,j,8))
2439 :
2440 0 : diagx(i,j) = ccaimp - strintx
2441 0 : diagy(i,j) = ccaimp - strinty
2442 : enddo ! ij
2443 :
2444 0 : end subroutine formDiag_step2
2445 :
2446 : !=======================================================================
2447 :
2448 : ! Compute global l^2 norm of a vector field (field_x, field_y)
2449 :
2450 0 : function global_norm (nx_block, ny_block, &
2451 : icellU , & ! LCOV_EXCL_LINE
2452 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2453 : field_x , field_y ) & ! LCOV_EXCL_LINE
2454 0 : result(norm)
2455 :
2456 : use ice_domain, only: distrb_info
2457 : use ice_domain_size, only: max_blocks
2458 :
2459 : integer (kind=int_kind), intent(in) :: &
2460 : nx_block, ny_block ! block dimensions
2461 :
2462 : integer (kind=int_kind), dimension (max_blocks), intent(in) :: &
2463 : icellU ! total count when iceumask is true
2464 :
2465 : integer (kind=int_kind), dimension (nx_block*ny_block,max_blocks), intent(in) :: &
2466 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2467 : indxUj ! compressed index in j-direction
2468 :
2469 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
2470 : field_x , & ! x-component of vector field ! LCOV_EXCL_LINE
2471 : field_y ! y-component of vector field
2472 :
2473 : real (kind=dbl_kind) :: &
2474 : norm ! l^2 norm of vector field
2475 :
2476 : ! local variables
2477 :
2478 : integer (kind=int_kind) :: &
2479 : i, j, ij, iblk
2480 :
2481 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
2482 : squared ! temporary array for summed squared components
2483 :
2484 : character(len=*), parameter :: subname = '(global_norm)'
2485 :
2486 : norm = sqrt(global_dot_product (nx_block , ny_block , &
2487 : icellU , & ! LCOV_EXCL_LINE
2488 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2489 : field_x , field_y , & ! LCOV_EXCL_LINE
2490 0 : field_x , field_y ))
2491 :
2492 0 : end function global_norm
2493 :
2494 : !=======================================================================
2495 :
2496 : ! Compute global dot product of two grid vectors, each split into X and Y components
2497 :
2498 0 : function global_dot_product (nx_block , ny_block , &
2499 : icellU , & ! LCOV_EXCL_LINE
2500 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2501 : vector1_x , vector1_y, & ! LCOV_EXCL_LINE
2502 : vector2_x , vector2_y) & ! LCOV_EXCL_LINE
2503 0 : result(dot_product)
2504 :
2505 : use ice_domain, only: distrb_info, ns_boundary_type
2506 : use ice_domain_size, only: max_blocks
2507 : use ice_fileunits, only: bfbflag
2508 :
2509 : integer (kind=int_kind), intent(in) :: &
2510 : nx_block, ny_block ! block dimensions
2511 :
2512 : integer (kind=int_kind), dimension (max_blocks), intent(in) :: &
2513 : icellU ! total count when iceumask is true
2514 :
2515 : integer (kind=int_kind), dimension (nx_block*ny_block,max_blocks), intent(in) :: &
2516 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2517 : indxUj ! compressed index in j-direction
2518 :
2519 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
2520 : vector1_x , & ! x-component of first vector ! LCOV_EXCL_LINE
2521 : vector1_y , & ! y-component of first vector ! LCOV_EXCL_LINE
2522 : vector2_x , & ! x-component of second vector ! LCOV_EXCL_LINE
2523 : vector2_y ! y-component of second vector
2524 :
2525 : real (kind=dbl_kind) :: &
2526 : dot_product ! l^2 norm of vector field
2527 :
2528 : ! local variables
2529 :
2530 : integer (kind=int_kind) :: &
2531 : i, j, ij, iblk
2532 :
2533 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
2534 0 : prod ! temporary array
2535 :
2536 : real (kind=dbl_kind), dimension(max_blocks) :: &
2537 0 : dot ! temporary scalar for accumulating the result
2538 :
2539 : character(len=*), parameter :: subname = '(global_dot_product)'
2540 :
2541 0 : prod = c0
2542 0 : dot = c0
2543 :
2544 0 : !$OMP PARALLEL DO PRIVATE(i, j, ij, iblk)
2545 0 : do iblk = 1, nblocks
2546 0 : do ij = 1, icellU(iblk)
2547 0 : i = indxUi(ij, iblk)
2548 0 : j = indxUj(ij, iblk)
2549 0 : prod(i,j,iblk) = vector1_x(i,j,iblk)*vector2_x(i,j,iblk) + vector1_y(i,j,iblk)*vector2_y(i,j,iblk)
2550 0 : dot(iblk) = dot(iblk) + prod(i,j,iblk)
2551 : enddo ! ij
2552 : enddo
2553 : !$OMP END PARALLEL DO
2554 :
2555 : ! Use faster local summation result for several bfbflag settings.
2556 : ! The local implementation sums over each block, sums over local
2557 : ! blocks, and calls global_sum on a scalar and should be just as accurate as
2558 : ! bfbflag = 'off', 'lsum8', and 'lsum4' without the extra copies and overhead
2559 : ! in the more general array global_sum. But use the array global_sum
2560 : ! if bfbflag is more strict or for tripole grids (requires special masking)
2561 0 : if (ns_boundary_type /= 'tripole' .and. ns_boundary_type /= 'tripoleT' .and. &
2562 : (bfbflag == 'off' .or. bfbflag == 'lsum8' .or. bfbflag == 'lsum4')) then
2563 0 : dot_product = global_sum(sum(dot), distrb_info)
2564 : else
2565 0 : dot_product = global_sum(prod, distrb_info, field_loc_NEcorner)
2566 : endif
2567 :
2568 0 : end function global_dot_product
2569 :
2570 : !=======================================================================
2571 :
2572 : ! Convert a grid function (tpu,tpv) to a one dimensional vector
2573 :
2574 0 : subroutine arrays_to_vec (nx_block, ny_block , &
2575 : nblocks , max_blocks, & ! LCOV_EXCL_LINE
2576 : icellU , ntot , & ! LCOV_EXCL_LINE
2577 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2578 : tpu , tpv , & ! LCOV_EXCL_LINE
2579 0 : outvec)
2580 :
2581 : integer (kind=int_kind), intent(in) :: &
2582 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2583 : nblocks, & ! nb of blocks ! LCOV_EXCL_LINE
2584 : max_blocks, & ! max nb of blocks ! LCOV_EXCL_LINE
2585 : ntot ! size of problem for Anderson
2586 :
2587 : integer (kind=int_kind), dimension (max_blocks), intent(in) :: &
2588 : icellU
2589 :
2590 : integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: &
2591 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2592 : indxUj ! compressed index in j-direction
2593 :
2594 : real (kind=dbl_kind), dimension (nx_block,ny_block, max_blocks), intent(in) :: &
2595 : tpu , & ! x-component of vector ! LCOV_EXCL_LINE
2596 : tpv ! y-component of vector
2597 :
2598 : real (kind=dbl_kind), dimension (ntot), intent(out) :: &
2599 : outvec ! output 1D vector
2600 :
2601 : ! local variables
2602 :
2603 : integer (kind=int_kind) :: &
2604 : i, j, iblk, tot, ij
2605 :
2606 : character(len=*), parameter :: subname = '(arrays_to_vec)'
2607 :
2608 : !-----------------------------------------------------------------
2609 : ! form vector (converts from max_blocks arrays to single vector)
2610 : !-----------------------------------------------------------------
2611 :
2612 0 : outvec(:) = c0
2613 0 : tot = 0
2614 :
2615 0 : do iblk = 1, nblocks
2616 0 : do ij = 1, icellU(iblk)
2617 0 : i = indxUi(ij, iblk)
2618 0 : j = indxUj(ij, iblk)
2619 0 : tot = tot + 1
2620 0 : outvec(tot) = tpu(i, j, iblk)
2621 0 : tot = tot + 1
2622 0 : outvec(tot) = tpv(i, j, iblk)
2623 : enddo
2624 : enddo ! ij
2625 :
2626 0 : end subroutine arrays_to_vec
2627 :
2628 : !=======================================================================
2629 :
2630 : ! Convert one dimensional vector to a grid function (tpu,tpv)
2631 :
2632 0 : subroutine vec_to_arrays (nx_block, ny_block , &
2633 : nblocks , max_blocks, & ! LCOV_EXCL_LINE
2634 : icellU , ntot , & ! LCOV_EXCL_LINE
2635 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2636 : invec , & ! LCOV_EXCL_LINE
2637 0 : tpu , tpv)
2638 :
2639 : integer (kind=int_kind), intent(in) :: &
2640 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2641 : nblocks, & ! nb of blocks ! LCOV_EXCL_LINE
2642 : max_blocks, & ! max nb of blocks ! LCOV_EXCL_LINE
2643 : ntot ! size of problem for Anderson
2644 :
2645 : integer (kind=int_kind), dimension (max_blocks), intent(in) :: &
2646 : icellU
2647 :
2648 : integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: &
2649 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2650 : indxUj ! compressed index in j-direction
2651 :
2652 : real (kind=dbl_kind), dimension (ntot), intent(in) :: &
2653 : invec ! input 1D vector
2654 :
2655 : real (kind=dbl_kind), dimension (nx_block,ny_block, max_blocks), intent(out) :: &
2656 : tpu , & ! x-component of vector ! LCOV_EXCL_LINE
2657 : tpv ! y-component of vector
2658 :
2659 : ! local variables
2660 :
2661 : integer (kind=int_kind) :: &
2662 : i, j, iblk, tot, ij
2663 :
2664 : character(len=*), parameter :: subname = '(vec_to_arrays)'
2665 :
2666 : !-----------------------------------------------------------------
2667 : ! form arrays (converts from vector to the max_blocks arrays)
2668 : !-----------------------------------------------------------------
2669 :
2670 0 : tpu(:,:,:) = c0
2671 0 : tpv(:,:,:) = c0
2672 0 : tot = 0
2673 :
2674 0 : do iblk = 1, nblocks
2675 0 : do ij = 1, icellU(iblk)
2676 0 : i = indxUi(ij, iblk)
2677 0 : j = indxUj(ij, iblk)
2678 0 : tot = tot + 1
2679 0 : tpu(i, j, iblk) = invec(tot)
2680 0 : tot = tot + 1
2681 0 : tpv(i, j, iblk) = invec(tot)
2682 : enddo
2683 : enddo! ij
2684 :
2685 0 : end subroutine vec_to_arrays
2686 :
2687 : !=======================================================================
2688 :
2689 : ! Update Q and R factors after deletion of the 1st column of G_diff
2690 : !
2691 : ! author: P. Blain ECCC
2692 : !
2693 : ! adapted from :
2694 : ! H. F. Walker, “Anderson Acceleration: Algorithms and Implementations”
2695 : ! [Online]. Available: https://users.wpi.edu/~walker/Papers/anderson_accn_algs_imps.pdf
2696 :
2697 0 : subroutine qr_delete(Q, R)
2698 :
2699 : real (kind=dbl_kind), intent(inout) :: &
2700 : Q(:,:), & ! Q factor ! LCOV_EXCL_LINE
2701 : R(:,:) ! R factor
2702 :
2703 : ! local variables
2704 :
2705 : integer (kind=int_kind) :: &
2706 : i, j, k, & ! loop indices ! LCOV_EXCL_LINE
2707 : m, n ! size of Q matrix
2708 :
2709 : real (kind=dbl_kind) :: &
2710 0 : temp, c, s
2711 :
2712 : character(len=*), parameter :: subname = '(qr_delete)'
2713 :
2714 0 : n = size(Q, 1)
2715 0 : m = size(Q, 2)
2716 0 : do i = 1, m-1
2717 0 : temp = sqrt(R(i, i+1)**2 + R(i+1, i+1)**2)
2718 0 : c = R(i , i+1) / temp
2719 0 : s = R(i+1, i+1) / temp
2720 0 : R(i , i+1) = temp
2721 0 : R(i+1, i+1) = 0
2722 0 : if (i < m-1) then
2723 0 : do j = i+2, m
2724 0 : temp = c*R(i, j) + s*R(i+1, j)
2725 0 : R(i+1, j) = -s*R(i, j) + c*R(i+1, j)
2726 0 : R(i , j) = temp
2727 : enddo
2728 : endif
2729 0 : do k = 1, n
2730 0 : temp = c*Q(k, i) + s*Q(k, i+1);
2731 0 : Q(k, i+1) = -s*Q(k, i) + c*Q(k, i+1);
2732 0 : Q(k, i) = temp
2733 : enddo
2734 : enddo
2735 0 : R(:, 1:m-1) = R(:, 2:m)
2736 :
2737 0 : end subroutine qr_delete
2738 :
2739 : !=======================================================================
2740 :
2741 : ! FGMRES: Flexible generalized minimum residual method (with restarts).
2742 : ! Solves the linear system A x = b using GMRES with a varying (right) preconditioner
2743 : !
2744 : ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC
2745 :
2746 0 : subroutine fgmres (zetax2 , etax2 , &
2747 : Cb , vrel , & ! LCOV_EXCL_LINE
2748 : umassdti , & ! LCOV_EXCL_LINE
2749 : halo_info_mask , & ! LCOV_EXCL_LINE
2750 : bx , by , & ! LCOV_EXCL_LINE
2751 : diagx , diagy , & ! LCOV_EXCL_LINE
2752 : tolerance, maxinner, & ! LCOV_EXCL_LINE
2753 : maxouter , & ! LCOV_EXCL_LINE
2754 : solx , soly , & ! LCOV_EXCL_LINE
2755 : nbiter)
2756 :
2757 : use ice_boundary, only: ice_HaloUpdate
2758 : use ice_domain, only: maskhalo_dyn, halo_info
2759 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2760 :
2761 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
2762 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
2763 : etax2 ! etax2 = 2*eta (shear viscosity)
2764 :
2765 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
2766 : vrel , & ! coefficient for tauw ! LCOV_EXCL_LINE
2767 : Cb , & ! seabed stress coefficient ! LCOV_EXCL_LINE
2768 : umassdti ! mass of U-cell/dte (kg/m^2 s)
2769 :
2770 : type (ice_halo), intent(in) :: &
2771 : halo_info_mask ! ghost cell update info for masked halo
2772 :
2773 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
2774 : bx , & ! Right hand side of the linear system (x components) ! LCOV_EXCL_LINE
2775 : by , & ! Right hand side of the linear system (y components) ! LCOV_EXCL_LINE
2776 : diagx , & ! Diagonal of the system matrix (x components) ! LCOV_EXCL_LINE
2777 : diagy ! Diagonal of the system matrix (y components)
2778 :
2779 : real (kind=dbl_kind), intent(in) :: &
2780 : tolerance ! Tolerance to achieve. The algorithm terminates when the relative
2781 : ! residual is below tolerance
2782 :
2783 : integer (kind=int_kind), intent(in) :: &
2784 : maxinner, & ! Restart the method every maxinner inner (Arnoldi) iterations ! LCOV_EXCL_LINE
2785 : maxouter ! Maximum number of outer (restarts) iterations
2786 : ! Iteration will stop after maxinner*maxouter Arnoldi steps
2787 : ! even if the specified tolerance has not been achieved
2788 :
2789 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(inout) :: &
2790 : solx , & ! Initial guess on input, approximate solution on output (x components) ! LCOV_EXCL_LINE
2791 : soly ! Initial guess on input, approximate solution on output (y components)
2792 :
2793 : integer (kind=int_kind), intent(out) :: &
2794 : nbiter ! Total number of Arnoldi iterations performed
2795 :
2796 : ! local variables
2797 :
2798 : integer (kind=int_kind) :: &
2799 : iblk , & ! block index ! LCOV_EXCL_LINE
2800 : ij , & ! index for indx[t|u][i|j] ! LCOV_EXCL_LINE
2801 : i, j ! grid indices
2802 :
2803 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
2804 : workspace_x , & ! work vector (x components) ! LCOV_EXCL_LINE
2805 0 : workspace_y ! work vector (y components)
2806 :
2807 : real (kind=dbl_kind), dimension (max_blocks) :: &
2808 0 : norm_squared ! array to accumulate squared norm of grid function over blocks
2809 :
2810 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks, maxinner+1) :: &
2811 : arnoldi_basis_x , & ! Arnoldi basis (x components) ! LCOV_EXCL_LINE
2812 0 : arnoldi_basis_y ! Arnoldi basis (y components)
2813 :
2814 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks, maxinner) :: &
2815 : orig_basis_x , & ! original basis (x components) ! LCOV_EXCL_LINE
2816 0 : orig_basis_y ! original basis (y components)
2817 :
2818 : real (kind=dbl_kind) :: &
2819 : norm_residual , & ! current L^2 norm of residual vector ! LCOV_EXCL_LINE
2820 : inverse_norm , & ! inverse of the norm of a vector ! LCOV_EXCL_LINE
2821 : norm_rhs , & ! L^2 norm of right-hand-side vector ! LCOV_EXCL_LINE
2822 0 : nu, t ! local temporary values
2823 :
2824 : integer (kind=int_kind) :: &
2825 : initer , & ! inner (Arnoldi) loop counter ! LCOV_EXCL_LINE
2826 : outiter , & ! outer (restarts) loop counter ! LCOV_EXCL_LINE
2827 : nextit , & ! nextit == initer+1 ! LCOV_EXCL_LINE
2828 : it, k, ii, jj ! reusable loop counters
2829 :
2830 : real (kind=dbl_kind), dimension(maxinner+1) :: &
2831 : rot_cos , & ! cosine elements of Givens rotations ! LCOV_EXCL_LINE
2832 : rot_sin , & ! sine elements of Givens rotations ! LCOV_EXCL_LINE
2833 0 : rhs_hess ! right hand side vector of the Hessenberg (least squares) system
2834 :
2835 : real (kind=dbl_kind), dimension(maxinner+1, maxinner) :: &
2836 0 : hessenberg ! system matrix of the Hessenberg (least squares) system
2837 :
2838 : character (len=char_len) :: &
2839 : precond_type ! type of preconditioner
2840 :
2841 : real (kind=dbl_kind) :: &
2842 0 : relative_tolerance ! relative_tolerance, i.e. tolerance*norm(initial residual)
2843 :
2844 : character(len=*), parameter :: subname = '(fgmres)'
2845 :
2846 : ! Here we go !
2847 :
2848 : ! Initialize
2849 0 : outiter = 0
2850 0 : nbiter = 0
2851 :
2852 0 : norm_squared = c0
2853 0 : precond_type = precond
2854 :
2855 : ! Cells with no ice should be zero-initialized
2856 0 : workspace_x = c0
2857 0 : workspace_y = c0
2858 0 : arnoldi_basis_x = c0
2859 0 : arnoldi_basis_y = c0
2860 :
2861 : ! solution is zero if RHS is zero
2862 : norm_rhs = global_norm(nx_block, ny_block, &
2863 : icellU , & ! LCOV_EXCL_LINE
2864 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2865 0 : bx , by )
2866 0 : if (norm_rhs == c0) then
2867 0 : solx = bx
2868 0 : soly = by
2869 0 : return
2870 : endif
2871 :
2872 : ! Residual of the initial iterate
2873 :
2874 0 : !$OMP PARALLEL DO PRIVATE(iblk)
2875 0 : do iblk = 1, nblocks
2876 : call matvec (nx_block , ny_block , &
2877 : icellU (iblk) , icellT (iblk), & ! LCOV_EXCL_LINE
2878 : indxUi (:,iblk) , indxUj (:,iblk), & ! LCOV_EXCL_LINE
2879 : indxTi (:,iblk) , indxTj (:,iblk), & ! LCOV_EXCL_LINE
2880 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
2881 : dxhy (:,:,iblk) , dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
2882 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
2883 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
2884 : solx (:,:,iblk) , soly (:,:,iblk), & ! LCOV_EXCL_LINE
2885 : vrel (:,:,iblk) , Cb (:,:,iblk), & ! LCOV_EXCL_LINE
2886 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
2887 : umassdti (:,:,iblk) , fmU (:,:,iblk), & ! LCOV_EXCL_LINE
2888 : uarear (:,:,iblk) , & ! LCOV_EXCL_LINE
2889 0 : workspace_x(:,:,iblk) , workspace_y(:,:,iblk))
2890 : call residual_vec (nx_block , ny_block , &
2891 : icellU (iblk), & ! LCOV_EXCL_LINE
2892 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
2893 : bx (:,:,iblk), by (:,:,iblk), & ! LCOV_EXCL_LINE
2894 : workspace_x(:,:,iblk), workspace_y(:,:,iblk), & ! LCOV_EXCL_LINE
2895 : arnoldi_basis_x (:,:,iblk, 1), & ! LCOV_EXCL_LINE
2896 0 : arnoldi_basis_y (:,:,iblk, 1))
2897 : enddo
2898 : !$OMP END PARALLEL DO
2899 :
2900 : ! Start outer (restarts) loop
2901 0 : do
2902 : ! Compute norm of initial residual
2903 : norm_residual = global_norm(nx_block, ny_block, &
2904 : icellU , & ! LCOV_EXCL_LINE
2905 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2906 : arnoldi_basis_x(:,:,:, 1), & ! LCOV_EXCL_LINE
2907 0 : arnoldi_basis_y(:,:,:, 1))
2908 :
2909 0 : if (my_task == master_task .and. monitor_fgmres) then
2910 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_fgmres: iter_fgmres= ", nbiter, &
2911 0 : " fgmres_L2norm= ", norm_residual
2912 : endif
2913 :
2914 : ! Current guess is a good enough solution TODO: reactivate and test this
2915 : ! if (norm_residual < tolerance) then
2916 : ! return
2917 : ! end if
2918 :
2919 : ! Normalize the first Arnoldi vector
2920 0 : inverse_norm = c1 / norm_residual
2921 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
2922 0 : do iblk = 1, nblocks
2923 0 : do ij = 1, icellU(iblk)
2924 0 : i = indxUi(ij, iblk)
2925 0 : j = indxUj(ij, iblk)
2926 :
2927 0 : arnoldi_basis_x(i, j, iblk, 1) = arnoldi_basis_x(i, j, iblk, 1) * inverse_norm
2928 0 : arnoldi_basis_y(i, j, iblk, 1) = arnoldi_basis_y(i, j, iblk, 1) * inverse_norm
2929 : enddo ! ij
2930 : enddo
2931 : !$OMP END PARALLEL DO
2932 :
2933 0 : if (outiter == 0) then
2934 0 : relative_tolerance = tolerance * norm_residual
2935 : end if
2936 :
2937 : ! Initialize 1-st term of RHS of Hessenberg system
2938 0 : rhs_hess(1) = norm_residual
2939 0 : rhs_hess(2:) = c0
2940 :
2941 0 : initer = 0
2942 :
2943 : ! Start of inner (Arnoldi) loop
2944 0 : do
2945 :
2946 0 : nbiter = nbiter + 1
2947 0 : initer = initer + 1
2948 0 : nextit = initer + 1
2949 : ! precondition the current Arnoldi vector
2950 : call precondition(zetax2 , etax2 , &
2951 : Cb , vrel , & ! LCOV_EXCL_LINE
2952 : umassdti , & ! LCOV_EXCL_LINE
2953 : halo_info_mask, & ! LCOV_EXCL_LINE
2954 : arnoldi_basis_x(:,:,:,initer), & ! LCOV_EXCL_LINE
2955 : arnoldi_basis_y(:,:,:,initer), & ! LCOV_EXCL_LINE
2956 : diagx , diagy , & ! LCOV_EXCL_LINE
2957 : precond_type, & ! LCOV_EXCL_LINE
2958 0 : workspace_x , workspace_y)
2959 0 : orig_basis_x(:,:,:,initer) = workspace_x
2960 0 : orig_basis_y(:,:,:,initer) = workspace_y
2961 :
2962 : ! Update workspace with boundary values
2963 0 : call stack_fields(workspace_x, workspace_y, fld2)
2964 0 : call ice_timer_start(timer_bound)
2965 0 : if (maskhalo_dyn) then
2966 : call ice_HaloUpdate (fld2, halo_info_mask, &
2967 0 : field_loc_NEcorner, field_type_vector)
2968 : else
2969 : call ice_HaloUpdate (fld2, halo_info, &
2970 0 : field_loc_NEcorner, field_type_vector)
2971 : endif
2972 0 : call ice_timer_stop(timer_bound)
2973 0 : call unstack_fields(fld2, workspace_x, workspace_y)
2974 :
2975 0 : !$OMP PARALLEL DO PRIVATE(iblk)
2976 0 : do iblk = 1, nblocks
2977 : call matvec (nx_block , ny_block , &
2978 : icellU (iblk) , icellT (iblk), & ! LCOV_EXCL_LINE
2979 : indxUi (:,iblk) , indxUj (:,iblk), & ! LCOV_EXCL_LINE
2980 : indxTi (:,iblk) , indxTj (:,iblk), & ! LCOV_EXCL_LINE
2981 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
2982 : dxhy (:,:,iblk) , dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
2983 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
2984 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
2985 : workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & ! LCOV_EXCL_LINE
2986 : vrel (:,:,iblk) , Cb (:,:,iblk), & ! LCOV_EXCL_LINE
2987 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
2988 : umassdti (:,:,iblk) , fmU (:,:,iblk), & ! LCOV_EXCL_LINE
2989 : uarear (:,:,iblk) , & ! LCOV_EXCL_LINE
2990 : arnoldi_basis_x(:,:,iblk,nextit), & ! LCOV_EXCL_LINE
2991 0 : arnoldi_basis_y(:,:,iblk,nextit))
2992 : enddo
2993 : !$OMP END PARALLEL DO
2994 :
2995 : ! Orthogonalize the new vector
2996 : call orthogonalize(ortho_type , initer , &
2997 : nextit , maxinner , & ! LCOV_EXCL_LINE
2998 : arnoldi_basis_x, arnoldi_basis_y, & ! LCOV_EXCL_LINE
2999 0 : hessenberg)
3000 :
3001 : ! Compute norm of new Arnoldi vector and update Hessenberg matrix
3002 0 : hessenberg(nextit,initer) = global_norm(nx_block, ny_block, &
3003 : icellU , & ! LCOV_EXCL_LINE
3004 : indxUi , indxUj , & ! LCOV_EXCL_LINE
3005 : arnoldi_basis_x(:,:,:, nextit), & ! LCOV_EXCL_LINE
3006 0 : arnoldi_basis_y(:,:,:, nextit))
3007 :
3008 : ! Watch out for happy breakdown
3009 0 : if (.not. almost_zero( hessenberg(nextit,initer) ) ) then
3010 : ! Normalize next Arnoldi vector
3011 0 : inverse_norm = c1 / hessenberg(nextit,initer)
3012 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3013 0 : do iblk = 1, nblocks
3014 0 : do ij = 1, icellU(iblk)
3015 0 : i = indxUi(ij, iblk)
3016 0 : j = indxUj(ij, iblk)
3017 :
3018 0 : arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit)*inverse_norm
3019 0 : arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit)*inverse_norm
3020 : enddo ! ij
3021 : enddo
3022 : !$OMP END PARALLEL DO
3023 : end if
3024 :
3025 : ! Apply previous Givens rotation to the last column of the Hessenberg matrix
3026 0 : if (initer > 1) then
3027 0 : do k = 2, initer
3028 0 : t = hessenberg(k-1, initer)
3029 0 : hessenberg(k-1, initer) = rot_cos(k-1)*t + rot_sin(k-1)*hessenberg(k, initer)
3030 0 : hessenberg(k, initer) = -rot_sin(k-1)*t + rot_cos(k-1)*hessenberg(k, initer)
3031 : end do
3032 : end if
3033 :
3034 : ! Compute and apply new Givens rotation
3035 0 : nu = sqrt(hessenberg(initer,initer)**2 + hessenberg(nextit,initer)**2)
3036 0 : if (.not. almost_zero(nu)) then
3037 0 : rot_cos(initer) = hessenberg(initer,initer) / nu
3038 0 : rot_sin(initer) = hessenberg(nextit,initer) / nu
3039 :
3040 0 : rhs_hess(nextit) = -rot_sin(initer) * rhs_hess(initer)
3041 0 : rhs_hess(initer) = rot_cos(initer) * rhs_hess(initer)
3042 :
3043 0 : hessenberg(initer,initer) = rot_cos(initer) * hessenberg(initer,initer) + rot_sin(initer) * hessenberg(nextit,initer)
3044 : end if
3045 :
3046 : ! Check for convergence
3047 0 : norm_residual = abs(rhs_hess(nextit))
3048 :
3049 0 : if (my_task == master_task .and. monitor_fgmres) then
3050 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_fgmres: iter_fgmres= ", nbiter, &
3051 0 : " fgmres_L2norm= ", norm_residual
3052 : endif
3053 :
3054 0 : if ((initer >= maxinner) .or. (norm_residual <= relative_tolerance)) then
3055 0 : exit
3056 : endif
3057 :
3058 : end do ! end of inner (Arnoldi) loop
3059 :
3060 : ! At this point either the maximum number of inner iterations
3061 : ! was reached or the absolute residual is below the scaled tolerance.
3062 :
3063 : ! Solve the (now upper triangular) system "hessenberg * sol_hess = rhs_hess"
3064 : ! (sol_hess is stored in rhs_hess)
3065 0 : rhs_hess(initer) = rhs_hess(initer) / hessenberg(initer,initer)
3066 0 : do ii = 2, initer
3067 0 : k = initer - ii + 1
3068 0 : t = rhs_hess(k)
3069 0 : do j = k + 1, initer
3070 0 : t = t - hessenberg(k,j) * rhs_hess(j)
3071 : end do
3072 0 : rhs_hess(k) = t / hessenberg(k,k)
3073 : end do
3074 :
3075 : ! Form linear combination to get new solution iterate
3076 0 : do it = 1, initer
3077 0 : t = rhs_hess(it)
3078 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3079 0 : do iblk = 1, nblocks
3080 0 : do ij = 1, icellU(iblk)
3081 0 : i = indxUi(ij, iblk)
3082 0 : j = indxUj(ij, iblk)
3083 :
3084 0 : solx(i, j, iblk) = solx(i, j, iblk) + t * orig_basis_x(i, j, iblk, it)
3085 0 : soly(i, j, iblk) = soly(i, j, iblk) + t * orig_basis_y(i, j, iblk, it)
3086 : enddo ! ij
3087 : enddo
3088 : !$OMP END PARALLEL DO
3089 : end do
3090 :
3091 : ! Increment outer loop counter and check for convergence
3092 0 : outiter = outiter + 1
3093 0 : if (norm_residual <= relative_tolerance .or. outiter >= maxouter) then
3094 0 : return
3095 : end if
3096 :
3097 : ! Solution is not convergent : compute residual vector and continue.
3098 :
3099 : ! The residual vector is computed here using (see Saad p. 177) :
3100 : ! \begin{equation}
3101 : ! r = V_{m+1} * Q_m^T * (\gamma_{m+1} * e_{m+1})
3102 : ! \end{equation}
3103 : ! where :
3104 : ! $r$ is the residual
3105 : ! $V_{m+1}$ is a matrix whose columns are the Arnoldi vectors from 1 to nextit (m+1)
3106 : ! $Q_m$ is the product of the Givens rotation : Q_m = G_m G_{m-1} ... G_1
3107 : ! $gamma_{m+1}$ is the last element of rhs_hess
3108 : ! $e_{m+1})$ is the unit vector (0, 0, ..., 1)^T \in \reals^{m+1}
3109 :
3110 : ! Apply the Givens rotation in reverse order to g := \gamma_{m+1} * e_{m+1},
3111 : ! store the result in rhs_hess
3112 0 : do it = 1, initer
3113 0 : jj = nextit - it + 1
3114 0 : rhs_hess(jj-1) = -rot_sin(jj-1) * rhs_hess(jj) ! + rot_cos(jj-1) * g(jj-1) (== 0)
3115 0 : rhs_hess(jj) = rot_cos(jj-1) * rhs_hess(jj) ! + rot_sin(jj-1) * g(jj-1) (== 0)
3116 : end do
3117 :
3118 : ! Compute the residual by multiplying V_{m+1} and rhs_hess
3119 0 : workspace_x = c0
3120 0 : workspace_y = c0
3121 0 : do it = 1, nextit
3122 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3123 0 : do iblk = 1, nblocks
3124 0 : do ij = 1, icellU(iblk)
3125 0 : i = indxUi(ij, iblk)
3126 0 : j = indxUj(ij, iblk)
3127 :
3128 0 : workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_x(i, j, iblk, it)
3129 0 : workspace_y(i, j, iblk) = workspace_y(i, j, iblk) + rhs_hess(it) * arnoldi_basis_y(i, j, iblk, it)
3130 : enddo ! ij
3131 : enddo
3132 : !$OMP END PARALLEL DO
3133 0 : arnoldi_basis_x(:,:,:,1) = workspace_x
3134 0 : arnoldi_basis_y(:,:,:,1) = workspace_y
3135 : end do
3136 : end do ! end of outer (restarts) loop
3137 :
3138 : end subroutine fgmres
3139 :
3140 : !=======================================================================
3141 :
3142 : ! PGMRES: Right-preconditioned generalized minimum residual method (with restarts).
3143 : ! Solves the linear A x = b using GMRES with a right preconditioner
3144 : !
3145 : ! authors: Stéphane Gaudreault, Abdessamad Qaddouri, Philippe Blain, ECCC
3146 :
3147 0 : subroutine pgmres (zetax2 , etax2 , &
3148 : Cb , vrel , & ! LCOV_EXCL_LINE
3149 : umassdti , & ! LCOV_EXCL_LINE
3150 : halo_info_mask, & ! LCOV_EXCL_LINE
3151 : bx , by , & ! LCOV_EXCL_LINE
3152 : diagx , diagy , & ! LCOV_EXCL_LINE
3153 : tolerance, maxinner, & ! LCOV_EXCL_LINE
3154 : maxouter , & ! LCOV_EXCL_LINE
3155 : solx , soly , & ! LCOV_EXCL_LINE
3156 : nbiter)
3157 :
3158 : use ice_boundary, only: ice_HaloUpdate
3159 : use ice_domain, only: maskhalo_dyn, halo_info
3160 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
3161 :
3162 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
3163 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
3164 : etax2 ! etax2 = 2*eta (shear viscosity)
3165 :
3166 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
3167 : vrel , & ! coefficient for tauw ! LCOV_EXCL_LINE
3168 : Cb , & ! seabed stress coefficient ! LCOV_EXCL_LINE
3169 : umassdti ! mass of U-cell/dte (kg/m^2 s)
3170 :
3171 : type (ice_halo), intent(in) :: &
3172 : halo_info_mask ! ghost cell update info for masked halo
3173 :
3174 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
3175 : bx , & ! Right hand side of the linear system (x components) ! LCOV_EXCL_LINE
3176 : by ! Right hand side of the linear system (y components)
3177 :
3178 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
3179 : diagx , & ! Diagonal of the system matrix (x components) ! LCOV_EXCL_LINE
3180 : diagy ! Diagonal of the system matrix (y components)
3181 :
3182 : real (kind=dbl_kind), intent(in) :: &
3183 : tolerance ! Tolerance to achieve. The algorithm terminates when the relative
3184 : ! residual is below tolerance
3185 :
3186 : integer (kind=int_kind), intent(in) :: &
3187 : maxinner, & ! Restart the method every maxinner inner (Arnoldi) iterations ! LCOV_EXCL_LINE
3188 : maxouter ! Maximum number of outer (restarts) iterations
3189 : ! Iteration will stop after maxinner*maxouter Arnoldi steps
3190 : ! even if the specified tolerance has not been achieved
3191 :
3192 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(inout) :: &
3193 : solx , & ! Initial guess on input, approximate solution on output (x components) ! LCOV_EXCL_LINE
3194 : soly ! Initial guess on input, approximate solution on output (y components)
3195 :
3196 : integer (kind=int_kind), intent(out) :: &
3197 : nbiter ! Total number of Arnoldi iterations performed
3198 :
3199 : ! local variables
3200 :
3201 : integer (kind=int_kind) :: &
3202 : iblk , & ! block index ! LCOV_EXCL_LINE
3203 : ij , & ! index for indx[t|u][i|j] ! LCOV_EXCL_LINE
3204 : i, j ! grid indices
3205 :
3206 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
3207 : workspace_x , & ! work vector (x components) ! LCOV_EXCL_LINE
3208 0 : workspace_y ! work vector (y components)
3209 :
3210 : real (kind=dbl_kind), dimension (max_blocks) :: &
3211 0 : norm_squared ! array to accumulate squared norm of grid function over blocks
3212 :
3213 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks, maxinner+1) :: &
3214 : arnoldi_basis_x , & ! Arnoldi basis (x components) ! LCOV_EXCL_LINE
3215 0 : arnoldi_basis_y ! Arnoldi basis (y components)
3216 :
3217 : real (kind=dbl_kind) :: &
3218 : norm_residual , & ! current L^2 norm of residual vector ! LCOV_EXCL_LINE
3219 : inverse_norm , & ! inverse of the norm of a vector ! LCOV_EXCL_LINE
3220 0 : nu, t ! local temporary values
3221 :
3222 : integer (kind=int_kind) :: &
3223 : initer , & ! inner (Arnoldi) loop counter ! LCOV_EXCL_LINE
3224 : outiter , & ! outer (restarts) loop counter ! LCOV_EXCL_LINE
3225 : nextit , & ! nextit == initer+1 ! LCOV_EXCL_LINE
3226 : it, k, ii, jj ! reusable loop counters
3227 :
3228 : real (kind=dbl_kind), dimension(maxinner+1) :: &
3229 : rot_cos , & ! cosine elements of Givens rotations ! LCOV_EXCL_LINE
3230 : rot_sin , & ! sine elements of Givens rotations ! LCOV_EXCL_LINE
3231 0 : rhs_hess ! right hand side vector of the Hessenberg (least squares) system
3232 :
3233 : real (kind=dbl_kind), dimension(maxinner+1, maxinner) :: &
3234 0 : hessenberg ! system matrix of the Hessenberg (least squares) system
3235 :
3236 : character(len=char_len) :: &
3237 : precond_type , & ! type of preconditioner ! LCOV_EXCL_LINE
3238 : ortho_type ! type of orthogonalization
3239 :
3240 : real (kind=dbl_kind) :: &
3241 0 : relative_tolerance ! relative_tolerance, i.e. tolerance*norm(initial residual)
3242 :
3243 : character(len=*), parameter :: subname = '(pgmres)'
3244 :
3245 : ! Here we go !
3246 :
3247 : ! Initialize
3248 0 : outiter = 0
3249 0 : nbiter = 0
3250 :
3251 0 : norm_squared = c0
3252 0 : precond_type = 'diag' ! Jacobi preconditioner
3253 0 : ortho_type = 'cgs' ! classical gram-schmidt TODO: try with MGS
3254 :
3255 : ! Cells with no ice should be zero-initialized
3256 0 : workspace_x = c0
3257 0 : workspace_y = c0
3258 0 : arnoldi_basis_x = c0
3259 0 : arnoldi_basis_y = c0
3260 :
3261 : ! Residual of the initial iterate
3262 :
3263 0 : !$OMP PARALLEL DO PRIVATE(iblk)
3264 0 : do iblk = 1, nblocks
3265 : call matvec (nx_block , ny_block , &
3266 : icellU (iblk) , icellT (iblk), & ! LCOV_EXCL_LINE
3267 : indxUi (:,iblk) , indxUj (:,iblk), & ! LCOV_EXCL_LINE
3268 : indxTi (:,iblk) , indxTj (:,iblk), & ! LCOV_EXCL_LINE
3269 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
3270 : dxhy (:,:,iblk) , dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
3271 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
3272 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
3273 : solx (:,:,iblk) , soly (:,:,iblk), & ! LCOV_EXCL_LINE
3274 : vrel (:,:,iblk) , Cb (:,:,iblk), & ! LCOV_EXCL_LINE
3275 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
3276 : umassdti (:,:,iblk) , fmU (:,:,iblk), & ! LCOV_EXCL_LINE
3277 : uarear (:,:,iblk) , & ! LCOV_EXCL_LINE
3278 0 : workspace_x(:,:,iblk) , workspace_y(:,:,iblk))
3279 : call residual_vec (nx_block , ny_block , &
3280 : icellU (iblk), & ! LCOV_EXCL_LINE
3281 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
3282 : bx (:,:,iblk), by (:,:,iblk), & ! LCOV_EXCL_LINE
3283 : workspace_x(:,:,iblk), workspace_y(:,:,iblk), & ! LCOV_EXCL_LINE
3284 : arnoldi_basis_x (:,:,iblk, 1), & ! LCOV_EXCL_LINE
3285 0 : arnoldi_basis_y (:,:,iblk, 1))
3286 : enddo
3287 : !$OMP END PARALLEL DO
3288 :
3289 : ! Start outer (restarts) loop
3290 0 : do
3291 : ! Compute norm of initial residual
3292 : norm_residual = global_norm(nx_block, ny_block, &
3293 : icellU , & ! LCOV_EXCL_LINE
3294 : indxUi , indxUj , & ! LCOV_EXCL_LINE
3295 : arnoldi_basis_x(:,:,:, 1), & ! LCOV_EXCL_LINE
3296 0 : arnoldi_basis_y(:,:,:, 1))
3297 :
3298 0 : if (my_task == master_task .and. monitor_pgmres) then
3299 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_pgmres: iter_pgmres= ", nbiter, &
3300 0 : " pgmres_L2norm= ", norm_residual
3301 : endif
3302 :
3303 : ! Current guess is a good enough solution
3304 : ! if (norm_residual < tolerance) then
3305 : ! return
3306 : ! end if
3307 :
3308 : ! Normalize the first Arnoldi vector
3309 0 : inverse_norm = c1 / norm_residual
3310 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3311 0 : do iblk = 1, nblocks
3312 0 : do ij = 1, icellU(iblk)
3313 0 : i = indxUi(ij, iblk)
3314 0 : j = indxUj(ij, iblk)
3315 :
3316 0 : arnoldi_basis_x(i, j, iblk, 1) = arnoldi_basis_x(i, j, iblk, 1) * inverse_norm
3317 0 : arnoldi_basis_y(i, j, iblk, 1) = arnoldi_basis_y(i, j, iblk, 1) * inverse_norm
3318 : enddo ! ij
3319 : enddo
3320 : !$OMP END PARALLEL DO
3321 :
3322 0 : if (outiter == 0) then
3323 0 : relative_tolerance = tolerance * norm_residual
3324 : end if
3325 :
3326 : ! Initialize 1-st term of RHS of Hessenberg system
3327 0 : rhs_hess(1) = norm_residual
3328 0 : rhs_hess(2:) = c0
3329 :
3330 0 : initer = 0
3331 :
3332 : ! Start of inner (Arnoldi) loop
3333 0 : do
3334 :
3335 0 : nbiter = nbiter + 1
3336 0 : initer = initer + 1
3337 0 : nextit = initer + 1
3338 :
3339 : ! precondition the current Arnoldi vector
3340 : call precondition(zetax2 , etax2 , &
3341 : Cb , vrel , & ! LCOV_EXCL_LINE
3342 : umassdti , & ! LCOV_EXCL_LINE
3343 : halo_info_mask, & ! LCOV_EXCL_LINE
3344 : arnoldi_basis_x(:,:,:,initer), & ! LCOV_EXCL_LINE
3345 : arnoldi_basis_y(:,:,:,initer), & ! LCOV_EXCL_LINE
3346 : diagx , diagy , & ! LCOV_EXCL_LINE
3347 : precond_type, & ! LCOV_EXCL_LINE
3348 0 : workspace_x , workspace_y)
3349 :
3350 : ! Update workspace with boundary values
3351 0 : call stack_fields(workspace_x, workspace_y, fld2)
3352 0 : call ice_timer_start(timer_bound)
3353 0 : if (maskhalo_dyn) then
3354 : call ice_HaloUpdate (fld2, halo_info_mask, &
3355 0 : field_loc_NEcorner, field_type_vector)
3356 : else
3357 : call ice_HaloUpdate (fld2, halo_info, &
3358 0 : field_loc_NEcorner, field_type_vector)
3359 : endif
3360 0 : call ice_timer_stop(timer_bound)
3361 0 : call unstack_fields(fld2, workspace_x, workspace_y)
3362 :
3363 0 : !$OMP PARALLEL DO PRIVATE(iblk)
3364 0 : do iblk = 1, nblocks
3365 : call matvec (nx_block , ny_block , &
3366 : icellU (iblk) , icellT (iblk), & ! LCOV_EXCL_LINE
3367 : indxUi (:,iblk) , indxUj (:,iblk), & ! LCOV_EXCL_LINE
3368 : indxTi (:,iblk) , indxTj (:,iblk), & ! LCOV_EXCL_LINE
3369 : dxT (:,:,iblk) , dyT (:,:,iblk), & ! LCOV_EXCL_LINE
3370 : dxhy (:,:,iblk) , dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
3371 : cxp (:,:,iblk) , cyp (:,:,iblk), & ! LCOV_EXCL_LINE
3372 : cxm (:,:,iblk) , cym (:,:,iblk), & ! LCOV_EXCL_LINE
3373 : workspace_x(:,:,iblk) , workspace_y(:,:,iblk), & ! LCOV_EXCL_LINE
3374 : vrel (:,:,iblk) , Cb (:,:,iblk), & ! LCOV_EXCL_LINE
3375 : zetax2 (:,:,iblk,:), etax2 (:,:,iblk,:), & ! LCOV_EXCL_LINE
3376 : umassdti (:,:,iblk) , fmU (:,:,iblk), & ! LCOV_EXCL_LINE
3377 : uarear (:,:,iblk) , & ! LCOV_EXCL_LINE
3378 : arnoldi_basis_x(:,:,iblk,nextit), & ! LCOV_EXCL_LINE
3379 0 : arnoldi_basis_y(:,:,iblk,nextit))
3380 : enddo
3381 : !$OMP END PARALLEL DO
3382 :
3383 : ! Orthogonalize the new vector
3384 : call orthogonalize(ortho_type , initer , &
3385 : nextit , maxinner , & ! LCOV_EXCL_LINE
3386 : arnoldi_basis_x, arnoldi_basis_y, & ! LCOV_EXCL_LINE
3387 0 : hessenberg)
3388 :
3389 : ! Compute norm of new Arnoldi vector and update Hessenberg matrix
3390 0 : hessenberg(nextit,initer) = global_norm(nx_block, ny_block, &
3391 : icellU , & ! LCOV_EXCL_LINE
3392 : indxUi , indxUj , & ! LCOV_EXCL_LINE
3393 : arnoldi_basis_x(:,:,:, nextit), & ! LCOV_EXCL_LINE
3394 0 : arnoldi_basis_y(:,:,:, nextit))
3395 :
3396 : ! Watch out for happy breakdown
3397 0 : if (.not. almost_zero( hessenberg(nextit,initer) ) ) then
3398 : ! Normalize next Arnoldi vector
3399 0 : inverse_norm = c1 / hessenberg(nextit,initer)
3400 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3401 0 : do iblk = 1, nblocks
3402 0 : do ij = 1, icellU(iblk)
3403 0 : i = indxUi(ij, iblk)
3404 0 : j = indxUj(ij, iblk)
3405 :
3406 0 : arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit)*inverse_norm
3407 0 : arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit)*inverse_norm
3408 : enddo ! ij
3409 : enddo
3410 : !$OMP END PARALLEL DO
3411 : end if
3412 :
3413 : ! Apply previous Givens rotation to the last column of the Hessenberg matrix
3414 0 : if (initer > 1) then
3415 0 : do k = 2, initer
3416 0 : t = hessenberg(k-1, initer)
3417 0 : hessenberg(k-1, initer) = rot_cos(k-1)*t + rot_sin(k-1)*hessenberg(k, initer)
3418 0 : hessenberg(k, initer) = -rot_sin(k-1)*t + rot_cos(k-1)*hessenberg(k, initer)
3419 : end do
3420 : end if
3421 :
3422 : ! Compute and apply new Givens rotation
3423 0 : nu = sqrt(hessenberg(initer,initer)**2 + hessenberg(nextit,initer)**2)
3424 0 : if (.not. almost_zero(nu)) then
3425 0 : rot_cos(initer) = hessenberg(initer,initer) / nu
3426 0 : rot_sin(initer) = hessenberg(nextit,initer) / nu
3427 :
3428 0 : rhs_hess(nextit) = -rot_sin(initer) * rhs_hess(initer)
3429 0 : rhs_hess(initer) = rot_cos(initer) * rhs_hess(initer)
3430 :
3431 0 : hessenberg(initer,initer) = rot_cos(initer) * hessenberg(initer,initer) + rot_sin(initer) * hessenberg(nextit,initer)
3432 : end if
3433 :
3434 : ! Check for convergence
3435 0 : norm_residual = abs(rhs_hess(nextit))
3436 :
3437 0 : if (my_task == master_task .and. monitor_pgmres) then
3438 0 : write(nu_diag, '(a,i4,a,d26.16)') "monitor_pgmres: iter_pgmres= ", nbiter, &
3439 0 : " pgmres_L2norm= ", norm_residual
3440 : endif
3441 :
3442 0 : if ((initer >= maxinner) .or. (norm_residual <= relative_tolerance)) then
3443 0 : exit
3444 : endif
3445 :
3446 : end do ! end of inner (Arnoldi) loop
3447 :
3448 : ! At this point either the maximum number of inner iterations
3449 : ! was reached or the absolute residual is below the scaled tolerance.
3450 :
3451 : ! Solve the (now upper triangular) system "hessenberg * sol_hess = rhs_hess"
3452 : ! (sol_hess is stored in rhs_hess)
3453 0 : rhs_hess(initer) = rhs_hess(initer) / hessenberg(initer,initer)
3454 0 : do ii = 2, initer
3455 0 : k = initer - ii + 1
3456 0 : t = rhs_hess(k)
3457 0 : do j = k + 1, initer
3458 0 : t = t - hessenberg(k,j) * rhs_hess(j)
3459 : end do
3460 0 : rhs_hess(k) = t / hessenberg(k,k)
3461 : end do
3462 :
3463 : ! Form linear combination to get new solution iterate
3464 0 : workspace_x = c0
3465 0 : workspace_y = c0
3466 0 : do it = 1, initer
3467 0 : t = rhs_hess(it)
3468 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3469 0 : do iblk = 1, nblocks
3470 0 : do ij = 1, icellU(iblk)
3471 0 : i = indxUi(ij, iblk)
3472 0 : j = indxUj(ij, iblk)
3473 :
3474 0 : workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + t * arnoldi_basis_x(i, j, iblk, it)
3475 0 : workspace_y(i, j, iblk) = workspace_y(i, j, iblk) + t * arnoldi_basis_y(i, j, iblk, it)
3476 : enddo ! ij
3477 : enddo
3478 : !$OMP END PARALLEL DO
3479 : end do
3480 :
3481 : ! Call preconditioner
3482 : call precondition(zetax2 , etax2 , &
3483 : Cb , vrel , & ! LCOV_EXCL_LINE
3484 : umassdti , & ! LCOV_EXCL_LINE
3485 : halo_info_mask, & ! LCOV_EXCL_LINE
3486 : workspace_x , workspace_y, & ! LCOV_EXCL_LINE
3487 : diagx , diagy , & ! LCOV_EXCL_LINE
3488 : precond_type, & ! LCOV_EXCL_LINE
3489 0 : workspace_x , workspace_y)
3490 :
3491 0 : solx = solx + workspace_x
3492 0 : soly = soly + workspace_y
3493 :
3494 : ! Increment outer loop counter and check for convergence
3495 0 : outiter = outiter + 1
3496 0 : if (norm_residual <= relative_tolerance .or. outiter >= maxouter) then
3497 0 : return
3498 : end if
3499 :
3500 : ! Solution is not convergent : compute residual vector and continue.
3501 :
3502 : ! The residual vector is computed here using (see Saad p. 177) :
3503 : ! \begin{equation}
3504 : ! r = V_{m+1} * Q_m^T * (\gamma_{m+1} * e_{m+1})
3505 : ! \end{equation}
3506 : ! where :
3507 : ! $r$ is the residual
3508 : ! $V_{m+1}$ is a matrix whose columns are the Arnoldi vectors from 1 to nextit (m+1)
3509 : ! $Q_m$ is the product of the Givens rotation : Q_m = G_m G_{m-1} ... G_1
3510 : ! $gamma_{m+1}$ is the last element of rhs_hess
3511 : ! $e_{m+1})$ is the unit vector (0, 0, ..., 1)^T \in \reals^{m+1}
3512 :
3513 : ! Apply the Givens rotation in reverse order to g := \gamma_{m+1} * e_{m+1},
3514 : ! store the result in rhs_hess
3515 0 : do it = 1, initer
3516 0 : jj = nextit - it + 1
3517 0 : rhs_hess(jj-1) = -rot_sin(jj-1) * rhs_hess(jj) ! + rot_cos(jj-1) * g(jj-1) (== 0)
3518 0 : rhs_hess(jj) = rot_cos(jj-1) * rhs_hess(jj) ! + rot_sin(jj-1) * g(jj-1) (== 0)
3519 : end do
3520 :
3521 : ! Compute the residual by multiplying V_{m+1} and rhs_hess
3522 0 : workspace_x = c0
3523 0 : workspace_y = c0
3524 0 : do it = 1, nextit
3525 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3526 0 : do iblk = 1, nblocks
3527 0 : do ij = 1, icellU(iblk)
3528 0 : i = indxUi(ij, iblk)
3529 0 : j = indxUj(ij, iblk)
3530 :
3531 0 : workspace_x(i, j, iblk) = workspace_x(i, j, iblk) + rhs_hess(it) * arnoldi_basis_x(i, j, iblk, it)
3532 0 : workspace_y(i, j, iblk) = workspace_y(i, j, iblk) + rhs_hess(it) * arnoldi_basis_y(i, j, iblk, it)
3533 : enddo ! ij
3534 : enddo
3535 : !$OMP END PARALLEL DO
3536 0 : arnoldi_basis_x(:,:,:,1) = workspace_x
3537 0 : arnoldi_basis_y(:,:,:,1) = workspace_y
3538 : end do
3539 : end do ! end of outer (restarts) loop
3540 :
3541 : end subroutine pgmres
3542 :
3543 : !=======================================================================
3544 :
3545 : ! Generic routine to precondition a vector
3546 : !
3547 : ! authors: Philippe Blain, ECCC
3548 :
3549 0 : subroutine precondition(zetax2 , etax2, &
3550 : Cb , vrel , & ! LCOV_EXCL_LINE
3551 : umassdti , & ! LCOV_EXCL_LINE
3552 : halo_info_mask, & ! LCOV_EXCL_LINE
3553 : vx , vy , & ! LCOV_EXCL_LINE
3554 : diagx , diagy, & ! LCOV_EXCL_LINE
3555 : precond_type, & ! LCOV_EXCL_LINE
3556 0 : wx , wy)
3557 :
3558 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: &
3559 : zetax2 , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
3560 : etax2 ! etax2 = 2*eta (shear viscosity)
3561 :
3562 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks), intent(in) :: &
3563 : vrel , & ! coefficient for tauw ! LCOV_EXCL_LINE
3564 : Cb , & ! seabed stress coefficient ! LCOV_EXCL_LINE
3565 : umassdti ! mass of U-cell/dte (kg/m^2 s)
3566 :
3567 : type (ice_halo), intent(in) :: &
3568 : halo_info_mask ! ghost cell update info for masked halo
3569 :
3570 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
3571 : vx , & ! input vector (x components) ! LCOV_EXCL_LINE
3572 : vy ! input vector (y components)
3573 :
3574 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: &
3575 : diagx , & ! diagonal of the system matrix (x components) ! LCOV_EXCL_LINE
3576 : diagy ! diagonal of the system matrix (y components)
3577 :
3578 : character (len=char_len), intent(in) :: &
3579 : precond_type ! type of preconditioner
3580 :
3581 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
3582 : wx , & ! preconditionned vector (x components) ! LCOV_EXCL_LINE
3583 : wy ! preconditionned vector (y components)
3584 :
3585 : ! local variables
3586 :
3587 : integer (kind=int_kind) :: &
3588 : iblk , & ! block index ! LCOV_EXCL_LINE
3589 : ij , & ! compressed index ! LCOV_EXCL_LINE
3590 : i, j ! grid indices
3591 :
3592 : real (kind=dbl_kind) :: &
3593 0 : tolerance ! Tolerance for PGMRES
3594 :
3595 : integer (kind=int_kind) :: &
3596 : maxinner ! Restart parameter for PGMRES
3597 :
3598 : integer (kind=int_kind) :: &
3599 : maxouter ! Maximum number of outer iterations for PGMRES
3600 :
3601 : integer (kind=int_kind) :: &
3602 : nbiter ! Total number of iteration PGMRES performed
3603 :
3604 : character(len=*), parameter :: subname = '(precondition)'
3605 :
3606 0 : if (precond_type == 'ident') then ! identity (no preconditioner)
3607 0 : wx = vx
3608 0 : wy = vy
3609 0 : elseif (precond_type == 'diag') then ! Jacobi preconditioner (diagonal)
3610 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3611 0 : do iblk = 1, nblocks
3612 0 : do ij = 1, icellU(iblk)
3613 0 : i = indxUi(ij, iblk)
3614 0 : j = indxUj(ij, iblk)
3615 :
3616 0 : wx(i,j,iblk) = vx(i,j,iblk) / diagx(i,j,iblk)
3617 0 : wy(i,j,iblk) = vy(i,j,iblk) / diagy(i,j,iblk)
3618 : enddo ! ij
3619 : enddo
3620 : !$OMP END PARALLEL DO
3621 0 : elseif (precond_type == 'pgmres') then ! PGMRES (Jacobi-preconditioned GMRES)
3622 : ! Initialize preconditioned vector to 0 ! TODO: try with wx = vx or vx/diagx
3623 0 : wx = c0
3624 0 : wy = c0
3625 0 : tolerance = reltol_pgmres
3626 0 : maxinner = dim_pgmres
3627 0 : maxouter = maxits_pgmres
3628 : call pgmres (zetax2, etax2 , &
3629 : Cb , vrel , & ! LCOV_EXCL_LINE
3630 : umassdti , & ! LCOV_EXCL_LINE
3631 : halo_info_mask , & ! LCOV_EXCL_LINE
3632 : vx , vy , & ! LCOV_EXCL_LINE
3633 : diagx , diagy , & ! LCOV_EXCL_LINE
3634 : tolerance, maxinner, & ! LCOV_EXCL_LINE
3635 : maxouter , & ! LCOV_EXCL_LINE
3636 : wx , wy , & ! LCOV_EXCL_LINE
3637 0 : nbiter)
3638 : else
3639 : call abort_ice(error_message='wrong preconditioner in ' // subname, &
3640 0 : file=__FILE__, line=__LINE__)
3641 : endif
3642 0 : end subroutine precondition
3643 :
3644 : !=======================================================================
3645 :
3646 : ! Generic routine to orthogonalize a vector (arnoldi_basis_[xy](:, :, :, nextit))
3647 : ! against a set of vectors (arnoldi_basis_[xy](:, :, :, 1:initer))
3648 : !
3649 : ! authors: Philippe Blain, ECCC
3650 :
3651 0 : subroutine orthogonalize(ortho_type , initer , &
3652 : nextit , maxinner , & ! LCOV_EXCL_LINE
3653 : arnoldi_basis_x, arnoldi_basis_y, & ! LCOV_EXCL_LINE
3654 0 : hessenberg)
3655 :
3656 : character(len=*), intent(in) :: &
3657 : ortho_type ! type of orthogonalization
3658 :
3659 : integer (kind=int_kind), intent(in) :: &
3660 : initer , & ! inner (Arnoldi) loop counter ! LCOV_EXCL_LINE
3661 : nextit , & ! nextit == initer+1 ! LCOV_EXCL_LINE
3662 : maxinner ! Restart the method every maxinner inner iterations
3663 :
3664 : real (kind=dbl_kind), dimension(nx_block, ny_block, max_blocks, maxinner+1), intent(inout) :: &
3665 : arnoldi_basis_x , & ! arnoldi basis (x components) ! LCOV_EXCL_LINE
3666 : arnoldi_basis_y ! arnoldi basis (y components)
3667 :
3668 : real (kind=dbl_kind), dimension(maxinner+1, maxinner), intent(inout) :: &
3669 : hessenberg ! system matrix of the Hessenberg (least squares) system
3670 :
3671 : ! local variables
3672 :
3673 : integer (kind=int_kind) :: &
3674 : it , & ! reusable loop counter ! LCOV_EXCL_LINE
3675 : iblk , & ! block index ! LCOV_EXCL_LINE
3676 : ij , & ! compressed index ! LCOV_EXCL_LINE
3677 : i, j ! grid indices
3678 :
3679 : character(len=*), parameter :: subname = '(orthogonalize)'
3680 :
3681 0 : if (trim(ortho_type) == 'cgs') then ! Classical Gram-Schmidt
3682 : ! Classical Gram-Schmidt orthogonalisation process
3683 : ! First loop of Gram-Schmidt (compute coefficients)
3684 0 : do it = 1, initer
3685 0 : hessenberg(it, initer) = global_dot_product(nx_block, ny_block, &
3686 : icellU , & ! LCOV_EXCL_LINE
3687 : indxUi , indxUj , & ! LCOV_EXCL_LINE
3688 : arnoldi_basis_x(:,:,:, it) , & ! LCOV_EXCL_LINE
3689 : arnoldi_basis_y(:,:,:, it) , & ! LCOV_EXCL_LINE
3690 : arnoldi_basis_x(:,:,:, nextit), & ! LCOV_EXCL_LINE
3691 0 : arnoldi_basis_y(:,:,:, nextit))
3692 : end do
3693 : ! Second loop of Gram-Schmidt (orthonormalize)
3694 0 : do it = 1, initer
3695 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3696 0 : do iblk = 1, nblocks
3697 0 : do ij = 1, icellU(iblk)
3698 0 : i = indxUi(ij, iblk)
3699 0 : j = indxUj(ij, iblk)
3700 :
3701 : arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit) &
3702 0 : - hessenberg(it,initer) * arnoldi_basis_x(i, j, iblk, it)
3703 : arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit) &
3704 0 : - hessenberg(it,initer) * arnoldi_basis_y(i, j, iblk, it)
3705 : enddo ! ij
3706 : enddo
3707 : !$OMP END PARALLEL DO
3708 : end do
3709 0 : elseif (trim(ortho_type) == 'mgs') then ! Modified Gram-Schmidt
3710 : ! Modified Gram-Schmidt orthogonalisation process
3711 0 : do it = 1, initer
3712 0 : hessenberg(it, initer) = global_dot_product(nx_block, ny_block, &
3713 : icellU , & ! LCOV_EXCL_LINE
3714 : indxUi , indxUj , & ! LCOV_EXCL_LINE
3715 : arnoldi_basis_x(:,:,:, it) , & ! LCOV_EXCL_LINE
3716 : arnoldi_basis_y(:,:,:, it) , & ! LCOV_EXCL_LINE
3717 : arnoldi_basis_x(:,:,:, nextit), & ! LCOV_EXCL_LINE
3718 0 : arnoldi_basis_y(:,:,:, nextit))
3719 :
3720 0 : !$OMP PARALLEL DO PRIVATE(iblk,ij,i,j)
3721 0 : do iblk = 1, nblocks
3722 0 : do ij = 1, icellU(iblk)
3723 0 : i = indxUi(ij, iblk)
3724 0 : j = indxUj(ij, iblk)
3725 :
3726 : arnoldi_basis_x(i, j, iblk, nextit) = arnoldi_basis_x(i, j, iblk, nextit) &
3727 0 : - hessenberg(it,initer) * arnoldi_basis_x(i, j, iblk, it)
3728 : arnoldi_basis_y(i, j, iblk, nextit) = arnoldi_basis_y(i, j, iblk, nextit) &
3729 0 : - hessenberg(it,initer) * arnoldi_basis_y(i, j, iblk, it)
3730 : enddo ! ij
3731 : enddo
3732 : !$OMP END PARALLEL DO
3733 : end do
3734 : else
3735 : call abort_ice(error_message='wrong orthonalization in ' // subname, &
3736 0 : file=__FILE__, line=__LINE__)
3737 : endif
3738 :
3739 0 : end subroutine orthogonalize
3740 :
3741 : !=======================================================================
3742 :
3743 : ! Check if value A is close to zero, up to machine precision
3744 : !
3745 : !author
3746 : ! Stéphane Gaudreault, ECCC -- June 2014
3747 : !
3748 : !revision
3749 : ! v4-80 - Gaudreault S. - gfortran compatibility
3750 : ! 2019 - Philippe Blain, ECCC - converted to CICE standards
3751 :
3752 0 : logical function almost_zero(A) result(retval)
3753 :
3754 : real (kind=dbl_kind), intent(in) :: A
3755 :
3756 : ! local variables
3757 :
3758 : character(len=*), parameter :: subname = '(almost_zero)'
3759 :
3760 : integer (kind=int8_kind) :: aBit
3761 : integer (kind=int8_kind), parameter :: two_complement = int(Z'80000000', kind=int8_kind)
3762 0 : aBit = 0
3763 0 : aBit = transfer(A, aBit)
3764 0 : if (aBit < 0) then
3765 0 : aBit = two_complement - aBit
3766 : end if
3767 : ! lexicographic order test with a tolerance of 1 adjacent float
3768 0 : retval = (abs(aBit) <= 1)
3769 :
3770 0 : end function almost_zero
3771 :
3772 : !=======================================================================
3773 :
3774 : end module ice_dyn_vp
3775 :
3776 : !=======================================================================
|
