Line data Source code
1 : !=======================================================================
2 : !
3 : ! Elastic-anisotropic sea ice dynamics model
4 : ! Computes ice velocity and deformation
5 : !
6 : ! See:
7 : !
8 : ! Wilchinsky, A.V. and D.L. Feltham (2006). Modelling the rheology of
9 : ! sea ice as a collection of diamond-shaped floes.
10 : ! Journal of Non-Newtonian Fluid Mechanics, 138(1), 22-32.
11 : !
12 : ! Tsamados, M., D.L. Feltham, and A.V. Wilchinsky (2013). Impact on new
13 : ! anisotropic rheology on simulations of Arctic sea ice. JGR, 118, 91-107.
14 : !
15 : ! authors: Michel Tsamados, CPOM
16 : ! David Schroeder, CPOM
17 :
18 : module ice_dyn_eap
19 :
20 : use ice_kinds_mod
21 : use ice_blocks, only: nx_block, ny_block
22 : use ice_communicate, only: my_task, master_task
23 : use ice_domain_size, only: max_blocks, ncat
24 : use ice_constants, only: c0, c1, c2, c3, c4, c12, p1, p2, p5, &
25 : p001, p027, p055, p111, p166, p222, p25, p333
26 : use ice_fileunits, only: nu_diag, nu_dump_eap, nu_restart_eap
27 : use ice_exit, only: abort_ice
28 : use ice_flux, only: rdg_shear
29 : ! use ice_timers, only: &
30 : ! ice_timer_start, ice_timer_stop, & ! LCOV_EXCL_LINE
31 : ! timer_tmp1, timer_tmp2, timer_tmp3, timer_tmp4, & ! LCOV_EXCL_LINE
32 : ! timer_tmp5, timer_tmp6, timer_tmp7, timer_tmp8, timer_tmp9
33 :
34 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
35 : use icepack_intfc, only: icepack_query_parameters
36 : use icepack_intfc, only: icepack_ice_strength
37 :
38 : implicit none
39 : private
40 : public :: eap, init_eap, write_restart_eap, read_restart_eap
41 :
42 : ! Look-up table needed for calculating structure tensor
43 : integer (int_kind), parameter :: &
44 : nx_yield = 41, & ! LCOV_EXCL_LINE
45 : ny_yield = 41, & ! LCOV_EXCL_LINE
46 : na_yield = 21
47 :
48 : real (kind=dbl_kind), dimension (nx_yield,ny_yield,na_yield) :: &
49 : s11r, s12r, s22r, s11s, s12s, s22s
50 :
51 : real (kind=dbl_kind), dimension (:,:,:), allocatable :: &
52 : a11_1, a11_2, a11_3, a11_4, & ! components of ! LCOV_EXCL_LINE
53 : a12_1, a12_2, a12_3, a12_4 ! structure tensor
54 :
55 : ! history
56 : real (kind=dbl_kind), dimension(:,:,:), allocatable, public :: &
57 : e11 , & ! components of strain rate tensor (1/s) ! LCOV_EXCL_LINE
58 : e12 , & ! LCOV_EXCL_LINE
59 : e22 , & ! LCOV_EXCL_LINE
60 : yieldstress11, & ! components of yield stress tensor (kg/s^2) ! LCOV_EXCL_LINE
61 : yieldstress12, & ! LCOV_EXCL_LINE
62 : yieldstress22, & ! LCOV_EXCL_LINE
63 : s11 , & ! components of stress tensor (kg/s^2) ! LCOV_EXCL_LINE
64 : s12 , & ! LCOV_EXCL_LINE
65 : s22 , & ! LCOV_EXCL_LINE
66 : a11 , & ! components of structure tensor () ! LCOV_EXCL_LINE
67 : a12
68 :
69 : ! private for reuse, set in init_eap
70 :
71 : real (kind=dbl_kind) :: &
72 : puny, pi, pi2, piq, pih
73 :
74 : real (kind=dbl_kind), parameter :: &
75 : kfriction = 0.45_dbl_kind
76 :
77 : real (kind=dbl_kind), save :: &
78 : invdx, invdy, invda, invsin
79 :
80 :
81 : !=======================================================================
82 :
83 : contains
84 :
85 : !=======================================================================
86 : ! Elastic-anisotropic-plastic dynamics driver
87 : ! based on subroutine evp
88 :
89 0 : subroutine eap (dt)
90 :
91 : #ifdef CICE_IN_NEMO
92 : ! Wind stress is set during this routine from the values supplied
93 : ! via NEMO (unless calc_strair is true). These values are supplied
94 : ! rotated on u grid and multiplied by aice. strairxT = 0 in this
95 : ! case so operations in dyn_prep1 are pointless but carried out to
96 : ! minimise code changes.
97 : #endif
98 :
99 : use ice_arrays_column, only: Cdn_ocn
100 : use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
101 : ice_HaloDestroy
102 : use ice_blocks, only: block, get_block
103 : use ice_constants, only: field_loc_center, field_loc_NEcorner, &
104 : field_type_scalar, field_type_vector
105 : use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn
106 : use ice_dyn_shared, only: fcor_blk, ndte, dtei, &
107 : denom1, uvel_init, vvel_init, arlx1i, & ! LCOV_EXCL_LINE
108 : dyn_prep1, dyn_prep2, stepu, dyn_finish, & ! LCOV_EXCL_LINE
109 : seabed_stress_factor_LKD, seabed_stress_factor_prob, & ! LCOV_EXCL_LINE
110 : seabed_stress_method, seabed_stress, & ! LCOV_EXCL_LINE
111 : stack_fields, unstack_fields, iceTmask, iceUmask, & ! LCOV_EXCL_LINE
112 : fld2, fld3, fld4
113 : use ice_flux, only: rdg_conv, strairxT, strairyT, &
114 : strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, & ! LCOV_EXCL_LINE
115 : strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, & ! LCOV_EXCL_LINE
116 : strax, stray, & ! LCOV_EXCL_LINE
117 : TbU, hwater, & ! LCOV_EXCL_LINE
118 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
119 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
120 : stress12_1, stress12_2, stress12_3, stress12_4
121 : use ice_grid, only: tmask, umask, dxT, dyT, dxhy, dyhx, cxp, cyp, cxm, cym, &
122 : tarear, uarear, grid_average_X2Y, & ! LCOV_EXCL_LINE
123 : grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
124 : use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, &
125 : aice_init, aice0, aicen, vicen, strength
126 : use ice_timers, only: timer_dynamics, timer_bound, &
127 : ice_timer_start, ice_timer_stop
128 :
129 : real (kind=dbl_kind), intent(in) :: &
130 : dt ! time step
131 :
132 : ! local variables
133 :
134 : integer (kind=int_kind) :: &
135 : ksub , & ! subcycle step ! LCOV_EXCL_LINE
136 : iblk , & ! block index ! LCOV_EXCL_LINE
137 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
138 : i, j, ij
139 :
140 : integer (kind=int_kind), dimension(max_blocks) :: &
141 : icellT , & ! no. of cells where iceTmask = .true. ! LCOV_EXCL_LINE
142 0 : icellU ! no. of cells where iceUmask = .true.
143 :
144 : integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: &
145 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
146 : indxTj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
147 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
148 0 : indxUj ! compressed index in j-direction
149 :
150 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
151 : uocnU , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
152 : vocnU , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
153 : ss_tltxU , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
154 : ss_tltyU , & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
155 : cdn_ocnU , & ! ocn drag coefficient ! LCOV_EXCL_LINE
156 : tmass , & ! total mass of ice and snow (kg/m^2) ! LCOV_EXCL_LINE
157 : waterxU , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
158 : wateryU , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
159 : forcexU , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
160 : forceyU , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
161 : umass , & ! total mass of ice and snow (u grid) ! LCOV_EXCL_LINE
162 0 : umassdti ! mass of U-cell/dte (kg/m^2 s)
163 :
164 : real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
165 0 : strtmp ! stress combinations for momentum equation
166 :
167 : logical (kind=log_kind) :: &
168 : calc_strair
169 :
170 : integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: &
171 0 : halomask ! ice mask for halo update
172 :
173 : type (ice_halo) :: &
174 : halo_info_mask ! ghost cell update info for masked halo
175 :
176 : type (block) :: &
177 : this_block ! block information for current block
178 :
179 : character(len=*), parameter :: subname = '(eap)'
180 :
181 0 : call ice_timer_start(timer_dynamics) ! dynamics
182 :
183 : !-----------------------------------------------------------------
184 : ! Initialize
185 : !-----------------------------------------------------------------
186 :
187 : ! This call is needed only if dt changes during runtime.
188 : ! call set_evp_parameters (dt)
189 :
190 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) SCHEDULE(runtime)
191 0 : do iblk = 1, nblocks
192 0 : do j = 1, ny_block
193 0 : do i = 1, nx_block
194 0 : rdg_conv (i,j,iblk) = c0
195 0 : rdg_shear(i,j,iblk) = c0 ! always zero. Could be moved
196 0 : divu (i,j,iblk) = c0
197 0 : shear(i,j,iblk) = c0
198 0 : e11(i,j,iblk) = c0
199 0 : e12(i,j,iblk) = c0
200 0 : e22(i,j,iblk) = c0
201 0 : s11(i,j,iblk) = c0
202 0 : s12(i,j,iblk) = c0
203 0 : s22(i,j,iblk) = c0
204 0 : yieldstress11(i,j,iblk) = c0
205 0 : yieldstress12(i,j,iblk) = c0
206 0 : yieldstress22(i,j,iblk) = c0
207 : enddo
208 : enddo
209 :
210 : !-----------------------------------------------------------------
211 : ! preparation for dynamics
212 : !-----------------------------------------------------------------
213 :
214 0 : this_block = get_block(blocks_ice(iblk),iblk)
215 0 : ilo = this_block%ilo
216 0 : ihi = this_block%ihi
217 0 : jlo = this_block%jlo
218 0 : jhi = this_block%jhi
219 :
220 : call dyn_prep1 (nx_block, ny_block, &
221 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
222 : aice (:,:,iblk), vice (:,:,iblk), & ! LCOV_EXCL_LINE
223 : vsno (:,:,iblk), tmask (:,:,iblk), & ! LCOV_EXCL_LINE
224 0 : tmass (:,:,iblk), iceTmask(:,:,iblk))
225 :
226 : enddo ! iblk
227 : !$OMP END PARALLEL DO
228 :
229 0 : call ice_timer_start(timer_bound)
230 : call ice_HaloUpdate (iceTmask, halo_info, &
231 0 : field_loc_center, field_type_scalar)
232 0 : call ice_timer_stop(timer_bound)
233 :
234 : !-----------------------------------------------------------------
235 : ! convert fields from T to U grid
236 : !-----------------------------------------------------------------
237 :
238 0 : call stack_fields(tmass, aice_init, cdn_ocn, fld3)
239 : call ice_HaloUpdate (fld3, halo_info, &
240 0 : field_loc_center, field_type_scalar)
241 0 : call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
242 : call ice_HaloUpdate (fld4, halo_info, &
243 0 : field_loc_center, field_type_vector)
244 0 : call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
245 0 : call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)
246 :
247 0 : call grid_average_X2Y('S', tmass , 'T' , umass , 'U')
248 0 : call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
249 0 : call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnU, 'U')
250 0 : call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
251 0 : call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
252 0 : call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
253 0 : call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')
254 :
255 : !----------------------------------------------------------------
256 : ! Set wind stress to values supplied via NEMO or other forcing
257 : ! This wind stress is rotated on u grid and multiplied by aice
258 : !----------------------------------------------------------------
259 :
260 0 : call icepack_query_parameters(calc_strair_out=calc_strair)
261 0 : call icepack_warnings_flush(nu_diag)
262 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
263 0 : file=__FILE__, line=__LINE__)
264 :
265 0 : if (.not. calc_strair) then
266 0 : call grid_average_X2Y('F', strax, grid_atm_dynu, strairxU, 'U')
267 0 : call grid_average_X2Y('F', stray, grid_atm_dynv, strairyU, 'U')
268 : else
269 : call ice_HaloUpdate (strairxT, halo_info, &
270 0 : field_loc_center, field_type_vector)
271 : call ice_HaloUpdate (strairyT, halo_info, &
272 0 : field_loc_center, field_type_vector)
273 0 : call grid_average_X2Y('F', strairxT, 'T', strairxU, 'U')
274 0 : call grid_average_X2Y('F', strairyT, 'T', strairyU, 'U')
275 : endif
276 :
277 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,ij,i,j) SCHEDULE(runtime)
278 0 : do iblk = 1, nblocks
279 :
280 : !-----------------------------------------------------------------
281 : ! more preparation for dynamics
282 : !-----------------------------------------------------------------
283 :
284 0 : this_block = get_block(blocks_ice(iblk),iblk)
285 0 : ilo = this_block%ilo
286 0 : ihi = this_block%ihi
287 0 : jlo = this_block%jlo
288 0 : jhi = this_block%jhi
289 :
290 : call dyn_prep2 (nx_block, ny_block, &
291 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
292 : icellT (iblk), icellU (iblk), & ! LCOV_EXCL_LINE
293 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
294 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
295 : aiU (:,:,iblk), umass (:,:,iblk), & ! LCOV_EXCL_LINE
296 : umassdti (:,:,iblk), fcor_blk (:,:,iblk), & ! LCOV_EXCL_LINE
297 : umask (:,:,iblk), & ! LCOV_EXCL_LINE
298 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
299 : strairxU (:,:,iblk), strairyU (:,:,iblk), & ! LCOV_EXCL_LINE
300 : ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & ! LCOV_EXCL_LINE
301 : iceTmask (:,:,iblk), iceUmask (:,:,iblk), & ! LCOV_EXCL_LINE
302 : fmU (:,:,iblk), dt, & ! LCOV_EXCL_LINE
303 : strtltxU (:,:,iblk), strtltyU (:,:,iblk), & ! LCOV_EXCL_LINE
304 : strocnxU (:,:,iblk), strocnyU (:,:,iblk), & ! LCOV_EXCL_LINE
305 : strintxU (:,:,iblk), strintyU (:,:,iblk), & ! LCOV_EXCL_LINE
306 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
307 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
308 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
309 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
310 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
311 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
312 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
313 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
314 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
315 : uvel_init (:,:,iblk), vvel_init (:,:,iblk), & ! LCOV_EXCL_LINE
316 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
317 0 : TbU (:,:,iblk))
318 :
319 : !-----------------------------------------------------------------
320 : ! Initialize structure tensor
321 : !-----------------------------------------------------------------
322 :
323 0 : do j = 1, ny_block
324 0 : do i = 1, nx_block
325 0 : if (.not.iceTmask(i,j,iblk)) then
326 0 : if (tmask(i,j,iblk)) then
327 : ! structure tensor
328 0 : a11_1(i,j,iblk) = p5
329 0 : a11_2(i,j,iblk) = p5
330 0 : a11_3(i,j,iblk) = p5
331 0 : a11_4(i,j,iblk) = p5
332 : else
333 0 : a11_1(i,j,iblk) = c0
334 0 : a11_2(i,j,iblk) = c0
335 0 : a11_3(i,j,iblk) = c0
336 0 : a11_4(i,j,iblk) = c0
337 : endif
338 0 : a12_1(i,j,iblk) = c0
339 0 : a12_2(i,j,iblk) = c0
340 0 : a12_3(i,j,iblk) = c0
341 0 : a12_4(i,j,iblk) = c0
342 : endif ! iceTmask
343 : enddo ! i
344 : enddo ! j
345 :
346 : !-----------------------------------------------------------------
347 : ! ice strength
348 : ! New strength used in Ukita Moritz rheology
349 : !-----------------------------------------------------------------
350 :
351 0 : strength(:,:,iblk) = c0 ! initialize
352 0 : do ij = 1, icellT(iblk)
353 0 : i = indxTi(ij, iblk)
354 0 : j = indxTj(ij, iblk)
355 : call icepack_ice_strength(ncat=ncat, &
356 : aice = aice (i,j, iblk), & ! LCOV_EXCL_LINE
357 : vice = vice (i,j, iblk), & ! LCOV_EXCL_LINE
358 : aice0 = aice0 (i,j, iblk), & ! LCOV_EXCL_LINE
359 : aicen = aicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
360 : vicen = vicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
361 0 : strength = strength(i,j, iblk) )
362 : enddo ! ij
363 : enddo ! iblk
364 : !$OMP END PARALLEL DO
365 :
366 0 : call icepack_warnings_flush(nu_diag)
367 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
368 0 : file=__FILE__, line=__LINE__)
369 :
370 0 : call ice_timer_start(timer_bound)
371 : call ice_HaloUpdate (strength, halo_info, &
372 0 : field_loc_center, field_type_scalar)
373 : ! velocities may have changed in dyn_prep2
374 0 : call stack_fields(uvel, vvel, fld2)
375 : call ice_HaloUpdate (fld2, halo_info, &
376 0 : field_loc_NEcorner, field_type_vector)
377 0 : call unstack_fields(fld2, uvel, vvel)
378 0 : call ice_timer_stop(timer_bound)
379 :
380 0 : if (maskhalo_dyn) then
381 0 : call ice_timer_start(timer_bound)
382 0 : halomask = 0
383 0 : where (iceUmask) halomask = 1
384 0 : call ice_HaloUpdate (halomask, halo_info, &
385 0 : field_loc_center, field_type_scalar)
386 0 : call ice_timer_stop(timer_bound)
387 0 : call ice_HaloMask(halo_info_mask, halo_info, halomask)
388 : endif
389 :
390 : !-----------------------------------------------------------------
391 : ! seabed stress factor TbU (TbU is part of Cb coefficient)
392 : !-----------------------------------------------------------------
393 :
394 0 : if (seabed_stress) then
395 0 : if ( seabed_stress_method == 'LKD' ) then
396 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
397 0 : do iblk = 1, nblocks
398 : call seabed_stress_factor_LKD (nx_block , ny_block , &
399 : icellU (iblk), & ! LCOV_EXCL_LINE
400 : indxUi (:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
401 : vice (:,:,iblk), aice(:,:,iblk), & ! LCOV_EXCL_LINE
402 0 : hwater(:,:,iblk), TbU (:,:,iblk))
403 : enddo
404 : !$OMP END PARALLEL DO
405 :
406 0 : elseif ( seabed_stress_method == 'probabilistic' ) then
407 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
408 0 : do iblk = 1, nblocks
409 : call seabed_stress_factor_prob (nx_block , ny_block , &
410 : icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & ! LCOV_EXCL_LINE
411 : icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
412 : aicen(:,:,:,iblk), vicen(:,:,:,iblk), & ! LCOV_EXCL_LINE
413 0 : hwater (:,:,iblk), TbU (:,:,iblk))
414 : enddo
415 : !$OMP END PARALLEL DO
416 : endif
417 : endif
418 :
419 0 : do ksub = 1,ndte ! subcycling
420 :
421 : !-----------------------------------------------------------------
422 : ! stress tensor equation, total surface stress
423 : !-----------------------------------------------------------------
424 :
425 0 : !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime)
426 0 : do iblk = 1, nblocks
427 :
428 : ! call ice_timer_start(timer_tmp1,iblk)
429 : call stress_eap (nx_block, ny_block, &
430 : ksub, ndte, & ! LCOV_EXCL_LINE
431 : icellT (iblk), & ! LCOV_EXCL_LINE
432 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
433 : arlx1i, denom1, & ! LCOV_EXCL_LINE
434 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
435 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
436 : dxhy (:,:,iblk), dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
437 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
438 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
439 : tarear (:,:,iblk), strength (:,:,iblk), & ! LCOV_EXCL_LINE
440 : a11_1 (:,:,iblk), a11_2 (:,:,iblk), & ! LCOV_EXCL_LINE
441 : a11_3 (:,:,iblk), a11_4 (:,:,iblk), & ! LCOV_EXCL_LINE
442 : a12_1 (:,:,iblk), a12_2 (:,:,iblk), & ! LCOV_EXCL_LINE
443 : a12_3 (:,:,iblk), a12_4 (:,:,iblk), & ! LCOV_EXCL_LINE
444 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
445 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
446 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
447 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
448 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
449 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
450 : shear (:,:,iblk), divu (:,:,iblk), & ! LCOV_EXCL_LINE
451 : e11 (:,:,iblk), e12 (:,:,iblk), & ! LCOV_EXCL_LINE
452 : e22 (:,:,iblk), & ! LCOV_EXCL_LINE
453 : s11 (:,:,iblk), s12 (:,:,iblk), & ! LCOV_EXCL_LINE
454 : s22 (:,:,iblk), & ! LCOV_EXCL_LINE
455 : yieldstress11 (:,:,iblk), & ! LCOV_EXCL_LINE
456 : yieldstress12 (:,:,iblk), & ! LCOV_EXCL_LINE
457 : yieldstress22 (:,:,iblk), & ! LCOV_EXCL_LINE
458 : ! rdg_conv (:,:,iblk), rdg_shear (:,:,iblk), & ! LCOV_EXCL_LINE
459 : rdg_conv (:,:,iblk), & ! LCOV_EXCL_LINE
460 0 : strtmp (:,:,:))
461 : ! call ice_timer_stop(timer_tmp1,iblk)
462 :
463 : !-----------------------------------------------------------------
464 : ! momentum equation
465 : !-----------------------------------------------------------------
466 :
467 : ! call ice_timer_start(timer_tmp2,iblk)
468 : call stepu (nx_block, ny_block, &
469 : icellU (iblk), Cdn_ocnU (:,:,iblk), & ! LCOV_EXCL_LINE
470 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
471 : aiU (:,:,iblk), strtmp (:,:,:), & ! LCOV_EXCL_LINE
472 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
473 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
474 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
475 : umassdti (:,:,iblk), fmU (:,:,iblk), & ! LCOV_EXCL_LINE
476 : uarear (:,:,iblk), & ! LCOV_EXCL_LINE
477 : strintxU (:,:,iblk), strintyU (:,:,iblk), & ! LCOV_EXCL_LINE
478 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
479 : uvel_init(:,:,iblk), vvel_init(:,:,iblk), & ! LCOV_EXCL_LINE
480 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
481 0 : TbU (:,:,iblk))
482 : ! call ice_timer_stop(timer_tmp2,iblk)
483 :
484 : !-----------------------------------------------------------------
485 : ! evolution of structure tensor A
486 : !-----------------------------------------------------------------
487 :
488 : ! call ice_timer_start(timer_tmp3,iblk)
489 0 : if (mod(ksub,10) == 1) then ! only called every 10th timestep
490 : call stepa (nx_block , ny_block , &
491 : dtei , icellT (iblk), & ! LCOV_EXCL_LINE
492 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
493 : a11 (:,:,iblk), a12 (:,:,iblk), & ! LCOV_EXCL_LINE
494 : a11_1 (:,:,iblk), a11_2 (:,:,iblk), & ! LCOV_EXCL_LINE
495 : a11_3 (:,:,iblk), a11_4 (:,:,iblk), & ! LCOV_EXCL_LINE
496 : a12_1 (:,:,iblk), a12_2 (:,:,iblk), & ! LCOV_EXCL_LINE
497 : a12_3 (:,:,iblk), a12_4 (:,:,iblk), & ! LCOV_EXCL_LINE
498 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
499 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
500 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
501 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
502 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
503 0 : stress12_3(:,:,iblk), stress12_4(:,:,iblk))
504 : endif
505 : ! call ice_timer_stop(timer_tmp3,iblk)
506 : enddo
507 : !$OMP END PARALLEL DO
508 :
509 0 : call stack_fields(uvel, vvel, fld2)
510 0 : call ice_timer_start(timer_bound)
511 0 : if (maskhalo_dyn) then
512 : call ice_HaloUpdate (fld2, halo_info_mask, &
513 0 : field_loc_NEcorner, field_type_vector)
514 : else
515 : call ice_HaloUpdate (fld2, halo_info, &
516 0 : field_loc_NEcorner, field_type_vector)
517 : endif
518 0 : call ice_timer_stop(timer_bound)
519 0 : call unstack_fields(fld2, uvel, vvel)
520 :
521 : enddo ! subcycling
522 :
523 0 : if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)
524 :
525 : !-----------------------------------------------------------------
526 : ! ice-ocean stress
527 : !-----------------------------------------------------------------
528 :
529 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
530 0 : do iblk = 1, nblocks
531 :
532 : call dyn_finish &
533 : (nx_block, ny_block, & ! LCOV_EXCL_LINE
534 : icellU (iblk), Cdn_ocnU(:,:,iblk), & ! LCOV_EXCL_LINE
535 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
536 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
537 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
538 : aiU (:,:,iblk), fmU (:,:,iblk), & ! LCOV_EXCL_LINE
539 0 : strocnxU(:,:,iblk), strocnyU(:,:,iblk))
540 :
541 : enddo
542 : !$OMP END PARALLEL DO
543 :
544 0 : call ice_timer_stop(timer_dynamics) ! dynamics
545 :
546 0 : end subroutine eap
547 :
548 : !=======================================================================
549 : ! Initialize parameters and variables needed for the eap dynamics
550 : ! (based on init_dyn)
551 :
552 0 : subroutine init_eap
553 :
554 : use ice_blocks, only: nx_block, ny_block
555 : use ice_domain, only: nblocks
556 : use ice_calendar, only: dt_dyn
557 : use ice_dyn_shared, only: init_dyn_shared
558 :
559 : ! local variables
560 :
561 : integer (kind=int_kind) :: &
562 : i, j, & ! LCOV_EXCL_LINE
563 : iblk ! block index
564 :
565 : real (kind=dbl_kind), parameter :: &
566 : eps6 = 1.0e-6_dbl_kind
567 :
568 : integer (kind=int_kind) :: &
569 : ix, iy, iz, ia, ierr
570 :
571 : integer (kind=int_kind), parameter :: &
572 : nz = 100
573 :
574 : real (kind=dbl_kind) :: &
575 : ainit, xinit, yinit, zinit, & ! LCOV_EXCL_LINE
576 : da, dx, dy, dz, & ! LCOV_EXCL_LINE
577 0 : phi
578 :
579 0 : real (kind=dbl_kind) :: invstressconviso
580 :
581 : character(len=*), parameter :: subname = '(init_eap)'
582 :
583 : call icepack_query_parameters(puny_out=puny, &
584 0 : pi_out=pi, pi2_out=pi2, piq_out=piq, pih_out=pih)
585 0 : call icepack_warnings_flush(nu_diag)
586 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
587 0 : file=__FILE__, line=__LINE__)
588 :
589 0 : phi = pi/c12 ! diamond shaped floe smaller angle (default phi = 30 deg)
590 :
591 0 : call init_dyn_shared(dt_dyn)
592 :
593 : allocate( a11_1 (nx_block,ny_block,max_blocks), &
594 : a11_2 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
595 : a11_3 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
596 : a11_4 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
597 : a12_1 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
598 : a12_2 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
599 : a12_3 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
600 : a12_4 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
601 : e11 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
602 : e12 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
603 : e22 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
604 : yieldstress11(nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
605 : yieldstress12(nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
606 : yieldstress22(nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
607 : s11 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
608 : s12 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
609 : s22 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
610 : a11 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
611 : a12 (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
612 0 : stat=ierr)
613 0 : if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory')
614 :
615 :
616 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
617 0 : do iblk = 1, nblocks
618 0 : do j = 1, ny_block
619 0 : do i = 1, nx_block
620 0 : e11 (i,j,iblk) = c0
621 0 : e12 (i,j,iblk) = c0
622 0 : e22 (i,j,iblk) = c0
623 0 : s11 (i,j,iblk) = c0
624 0 : s12 (i,j,iblk) = c0
625 0 : s22 (i,j,iblk) = c0
626 0 : yieldstress11(i,j,iblk) = c0
627 0 : yieldstress12(i,j,iblk) = c0
628 0 : yieldstress22(i,j,iblk) = c0
629 0 : a11_1 (i,j,iblk) = p5
630 0 : a11_2 (i,j,iblk) = p5
631 0 : a11_3 (i,j,iblk) = p5
632 0 : a11_4 (i,j,iblk) = p5
633 0 : a12_1 (i,j,iblk) = c0
634 0 : a12_2 (i,j,iblk) = c0
635 0 : a12_3 (i,j,iblk) = c0
636 0 : a12_4 (i,j,iblk) = c0
637 0 : rdg_shear (i,j,iblk) = c0
638 : enddo ! i
639 : enddo ! j
640 : enddo ! iblk
641 : !$OMP END PARALLEL DO
642 :
643 : !-----------------------------------------------------------------
644 : ! create lookup table for eap dynamics (see Appendix A1)
645 : !-----------------------------------------------------------------
646 :
647 0 : da = p5/real(na_yield-1,kind=dbl_kind)
648 0 : ainit = p5 - da
649 0 : dx = pi/real(nx_yield-1,kind=dbl_kind)
650 0 : xinit = pi + piq - dx
651 0 : dz = pi/real(nz,kind=dbl_kind)
652 0 : zinit = -pih
653 0 : dy = pi/real(ny_yield-1,kind=dbl_kind)
654 0 : yinit = -dy
655 0 : invdx = c1/dx
656 0 : invdy = c1/dy
657 0 : invda = c1/da
658 :
659 0 : do ia=1,na_yield
660 0 : do ix=1,nx_yield
661 0 : do iy=1,ny_yield
662 0 : s11r(ix,iy,ia) = c0
663 0 : s12r(ix,iy,ia) = c0
664 0 : s22r(ix,iy,ia) = c0
665 0 : s11s(ix,iy,ia) = c0
666 0 : s12s(ix,iy,ia) = c0
667 0 : s22s(ix,iy,ia) = c0
668 0 : if (ia <= na_yield-1) then
669 0 : do iz=1,nz
670 0 : s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
671 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
672 0 : s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
673 0 : s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
674 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
675 0 : s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
676 0 : s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
677 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
678 0 : s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
679 0 : s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
680 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
681 0 : s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
682 0 : s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
683 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
684 0 : s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
685 0 : s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
686 : exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & ! LCOV_EXCL_LINE
687 0 : s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(c2*phi)
688 : enddo
689 0 : if (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = c0
690 0 : if (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = c0
691 0 : if (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = c0
692 0 : if (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = c0
693 0 : if (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = c0
694 0 : if (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = c0
695 : else
696 0 : s11r(ix,iy,ia) = p5*s11kr(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
697 0 : s12r(ix,iy,ia) = p5*s12kr(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
698 0 : s22r(ix,iy,ia) = p5*s22kr(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
699 0 : s11s(ix,iy,ia) = p5*s11ks(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
700 0 : s12s(ix,iy,ia) = p5*s12ks(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
701 0 : s22s(ix,iy,ia) = p5*s22ks(xinit+ix*dx,yinit+iy*dy,c0,phi)/sin(c2*phi)
702 0 : if (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = c0
703 0 : if (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = c0
704 0 : if (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = c0
705 0 : if (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = c0
706 0 : if (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = c0
707 0 : if (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = c0
708 : endif
709 : enddo
710 : enddo
711 : enddo
712 :
713 : ! Factor to maintain the same stress as in EVP (see Section 3)
714 : ! Can be set to 1 otherwise
715 :
716 0 : invstressconviso = c1/(c1+kfriction*kfriction)
717 0 : invsin = c1/sin(pi2/c12) * invstressconviso
718 :
719 0 : end subroutine init_eap
720 :
721 : !=======================================================================
722 : ! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
723 :
724 0 : FUNCTION w1(a)
725 : double precision, intent(in) :: a
726 :
727 : real (kind=dbl_kind) :: w1
728 : character(len=*), parameter :: subname = '(w1)'
729 :
730 : w1 = - 223.87569446_dbl_kind &
731 : + 2361.2198663_dbl_kind*a & ! LCOV_EXCL_LINE
732 : - 10606.56079975_dbl_kind*a*a & ! LCOV_EXCL_LINE
733 : + 26315.50025642_dbl_kind*a*a*a & ! LCOV_EXCL_LINE
734 : - 38948.30444297_dbl_kind*a*a*a*a & ! LCOV_EXCL_LINE
735 : + 34397.72407466_dbl_kind*a*a*a*a*a & ! LCOV_EXCL_LINE
736 : - 16789.98003081_dbl_kind*a*a*a*a*a*a & ! LCOV_EXCL_LINE
737 0 : + 3495.82839237_dbl_kind*a*a*a*a*a*a*a
738 :
739 0 : end FUNCTION w1
740 :
741 : !=======================================================================
742 : ! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
743 :
744 0 : FUNCTION w2(a)
745 : double precision, intent(in) :: a
746 :
747 : real (kind=dbl_kind) :: w2
748 : character(len=*), parameter :: subname = '(w2)'
749 :
750 : w2 = - 6670.68911883_dbl_kind &
751 : + 70222.33061536_dbl_kind*a & ! LCOV_EXCL_LINE
752 : - 314871.71525448_dbl_kind*a*a & ! LCOV_EXCL_LINE
753 : + 779570.02793492_dbl_kind*a*a*a & ! LCOV_EXCL_LINE
754 : - 1151098.82436864_dbl_kind*a*a*a*a & ! LCOV_EXCL_LINE
755 : + 1013896.59464498_dbl_kind*a*a*a*a*a & ! LCOV_EXCL_LINE
756 : - 493379.44906738_dbl_kind*a*a*a*a*a*a & ! LCOV_EXCL_LINE
757 0 : + 102356.551518_dbl_kind*a*a*a*a*a*a*a
758 :
759 0 : end FUNCTION w2
760 :
761 : !=======================================================================
762 : ! Function : s11kr
763 :
764 0 : FUNCTION s11kr(x,y,z,phi)
765 :
766 : real (kind=dbl_kind), intent(in) :: &
767 : x,y,z,phi
768 :
769 : real (kind=dbl_kind) :: &
770 0 : s11kr, p
771 :
772 : real (kind=dbl_kind) :: &
773 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
774 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
775 : ! t1t2i11, t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
776 : ! t2t1i11, t2t1i12, t2t1i21, t2t1i22, & ! LCOV_EXCL_LINE
777 : d11, d12, d22, & ! LCOV_EXCL_LINE
778 : IIn1t2, IIn2t1, & ! LCOV_EXCL_LINE
779 : ! IIt1t2, & ! LCOV_EXCL_LINE
780 0 : Hen1t2, Hen2t1
781 :
782 : character(len=*), parameter :: subname = '(s11kr)'
783 :
784 0 : p = phi
785 :
786 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
787 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
788 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
789 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
790 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
791 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
792 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
793 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
794 : ! t1t2i11 = cos(z-p) * cos(z+p)
795 : ! t1t2i12 = cos(z-p) * sin(z+p)
796 : ! t1t2i21 = sin(z-p) * cos(z+p)
797 : ! t1t2i22 = sin(z-p) * sin(z+p)
798 : ! t2t1i11 = cos(z+p) * cos(z-p)
799 : ! t2t1i12 = cos(z+p) * sin(z-p)
800 : ! t2t1i21 = sin(z+p) * cos(z-p)
801 : ! t2t1i22 = sin(z+p) * sin(z-p)
802 : ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
803 : ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else x=x-pi/2
804 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
805 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
806 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
807 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
808 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
809 : ! IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
810 :
811 0 : if (-IIn1t2>=puny) then
812 0 : Hen1t2 = c1
813 : else
814 0 : Hen1t2 = c0
815 : endif
816 :
817 0 : if (-IIn2t1>=puny) then
818 0 : Hen2t1 = c1
819 : else
820 0 : Hen2t1 = c0
821 : endif
822 :
823 0 : s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11)
824 :
825 0 : end FUNCTION s11kr
826 :
827 : !=======================================================================
828 : ! Function : s12kr
829 :
830 0 : FUNCTION s12kr(x,y,z,phi)
831 :
832 : real (kind=dbl_kind), intent(in) :: &
833 : x,y,z,phi
834 :
835 : real (kind=dbl_kind) :: &
836 0 : s12kr, s12r0, s21r0, p
837 :
838 : real (kind=dbl_kind) :: &
839 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
840 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
841 : ! t1t2i11, t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
842 : ! t2t1i11, t2t1i12, t2t1i21, t2t1i22, & ! LCOV_EXCL_LINE
843 : d11, d12, d22, & ! LCOV_EXCL_LINE
844 : IIn1t2, IIn2t1, & ! LCOV_EXCL_LINE
845 : ! IIt1t2, & ! LCOV_EXCL_LINE
846 0 : Hen1t2, Hen2t1
847 :
848 : character(len=*), parameter :: subname = '(s12kr)'
849 :
850 0 : p = phi
851 :
852 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
853 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
854 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
855 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
856 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
857 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
858 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
859 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
860 : ! t1t2i11 = cos(z-p) * cos(z+p)
861 : ! t1t2i12 = cos(z-p) * sin(z+p)
862 : ! t1t2i21 = sin(z-p) * cos(z+p)
863 : ! t1t2i22 = sin(z-p) * sin(z+p)
864 : ! t2t1i11 = cos(z+p) * cos(z-p)
865 : ! t2t1i12 = cos(z+p) * sin(z-p)
866 : ! t2t1i21 = sin(z+p) * cos(z-p)
867 : ! t2t1i22 = sin(z+p) * sin(z-p)
868 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
869 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
870 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
871 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
872 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
873 : ! IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
874 :
875 0 : if (-IIn1t2>=puny) then
876 0 : Hen1t2 = c1
877 : else
878 0 : Hen1t2 = c0
879 : endif
880 :
881 0 : if (-IIn2t1>=puny) then
882 0 : Hen2t1 = c1
883 : else
884 0 : Hen2t1 = c0
885 : endif
886 :
887 0 : s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12)
888 0 : s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21)
889 0 : s12kr=p5*(s12r0+s21r0)
890 :
891 0 : end FUNCTION s12kr
892 :
893 : !=======================================================================
894 : ! Function : s22r
895 :
896 0 : FUNCTION s22kr(x,y,z,phi)
897 :
898 : real (kind=dbl_kind), intent(in) :: &
899 : x,y,z,phi
900 :
901 : real (kind=dbl_kind) :: &
902 0 : s22kr, p
903 :
904 : real (kind=dbl_kind) :: &
905 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
906 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
907 : ! t1t2i11, t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
908 : ! t2t1i11, t2t1i12, t2t1i21, t2t1i22, & ! LCOV_EXCL_LINE
909 : d11, d12, d22, & ! LCOV_EXCL_LINE
910 : IIn1t2, IIn2t1, & ! LCOV_EXCL_LINE
911 : ! IIt1t2, & ! LCOV_EXCL_LINE
912 0 : Hen1t2, Hen2t1
913 :
914 : character(len=*), parameter :: subname = '(s22kr)'
915 :
916 0 : p = phi
917 :
918 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
919 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
920 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
921 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
922 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
923 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
924 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
925 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
926 : ! t1t2i11 = cos(z-p) * cos(z+p)
927 : ! t1t2i12 = cos(z-p) * sin(z+p)
928 : ! t1t2i21 = sin(z-p) * cos(z+p)
929 : ! t1t2i22 = sin(z-p) * sin(z+p)
930 : ! t2t1i11 = cos(z+p) * cos(z-p)
931 : ! t2t1i12 = cos(z+p) * sin(z-p)
932 : ! t2t1i21 = sin(z+p) * cos(z-p)
933 : ! t2t1i22 = sin(z+p) * sin(z-p)
934 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
935 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
936 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
937 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
938 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
939 : ! IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
940 :
941 0 : if (-IIn1t2>=puny) then
942 0 : Hen1t2 = c1
943 : else
944 0 : Hen1t2 = c0
945 : endif
946 :
947 0 : if (-IIn2t1>=puny) then
948 0 : Hen2t1 = c1
949 : else
950 0 : Hen2t1 = c0
951 : endif
952 :
953 0 : s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22)
954 :
955 0 : end FUNCTION s22kr
956 :
957 : !=======================================================================
958 : ! Function : s11ks
959 :
960 0 : FUNCTION s11ks(x,y,z,phi)
961 :
962 : real (kind=dbl_kind), intent(in):: &
963 : x,y,z,phi
964 :
965 : real (kind=dbl_kind) :: &
966 0 : s11ks, p
967 :
968 : real (kind=dbl_kind) :: &
969 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
970 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
971 : t1t2i11, & ! LCOV_EXCL_LINE
972 : t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
973 : t2t1i11, & ! LCOV_EXCL_LINE
974 : ! t2t1i12, t2t1i21, t2t1i22, & ! LCOV_EXCL_LINE
975 : d11, d12, d22, & ! LCOV_EXCL_LINE
976 : IIn1t2, IIn2t1, IIt1t2, & ! LCOV_EXCL_LINE
977 0 : Hen1t2, Hen2t1
978 :
979 : character(len=*), parameter :: subname = '(s11ks)'
980 :
981 0 : p = phi
982 :
983 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
984 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
985 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
986 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
987 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
988 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
989 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
990 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
991 0 : t1t2i11 = cos(z-p) * cos(z+p)
992 0 : t1t2i12 = cos(z-p) * sin(z+p)
993 0 : t1t2i21 = sin(z-p) * cos(z+p)
994 0 : t1t2i22 = sin(z-p) * sin(z+p)
995 0 : t2t1i11 = cos(z+p) * cos(z-p)
996 : ! t2t1i12 = cos(z+p) * sin(z-p)
997 : ! t2t1i21 = sin(z+p) * cos(z-p)
998 : ! t2t1i22 = sin(z+p) * sin(z-p)
999 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
1000 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
1001 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
1002 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
1003 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
1004 0 : IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
1005 :
1006 0 : if (-IIn1t2>=puny) then
1007 0 : Hen1t2 = c1
1008 : else
1009 0 : Hen1t2 = c0
1010 : endif
1011 :
1012 0 : if (-IIn2t1>=puny) then
1013 0 : Hen2t1 = c1
1014 : else
1015 0 : Hen2t1 = c0
1016 : endif
1017 :
1018 0 : s11ks = sign(c1,IIt1t2+puny)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11)
1019 :
1020 0 : end FUNCTION s11ks
1021 :
1022 : !=======================================================================
1023 : ! Function : s12ks
1024 :
1025 0 : FUNCTION s12ks(x,y,z,phi)
1026 :
1027 : real (kind=dbl_kind), intent(in) :: &
1028 : x,y,z,phi
1029 :
1030 : real (kind=dbl_kind) :: &
1031 0 : s12ks,s12s0,s21s0,p
1032 :
1033 : real (kind=dbl_kind) :: &
1034 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
1035 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
1036 : t1t2i11, t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
1037 : ! t2t1i11, t2t1i22, & ! LCOV_EXCL_LINE
1038 : t2t1i12, t2t1i21, & ! LCOV_EXCL_LINE
1039 : d11, d12, d22, & ! LCOV_EXCL_LINE
1040 : IIn1t2, IIn2t1, IIt1t2, & ! LCOV_EXCL_LINE
1041 0 : Hen1t2, Hen2t1
1042 :
1043 : character(len=*), parameter :: subname = '(s12ks)'
1044 :
1045 0 : p =phi
1046 :
1047 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
1048 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
1049 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
1050 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
1051 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
1052 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
1053 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
1054 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
1055 0 : t1t2i11 = cos(z-p) * cos(z+p)
1056 0 : t1t2i12 = cos(z-p) * sin(z+p)
1057 0 : t1t2i21 = sin(z-p) * cos(z+p)
1058 0 : t1t2i22 = sin(z-p) * sin(z+p)
1059 : ! t2t1i11 = cos(z+p) * cos(z-p)
1060 0 : t2t1i12 = cos(z+p) * sin(z-p)
1061 0 : t2t1i21 = sin(z+p) * cos(z-p)
1062 : ! t2t1i22 = sin(z+p) * sin(z-p)
1063 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
1064 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
1065 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
1066 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
1067 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
1068 0 : IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
1069 :
1070 0 : if (-IIn1t2>=puny) then
1071 0 : Hen1t2 = c1
1072 : else
1073 0 : Hen1t2 = c0
1074 : endif
1075 :
1076 0 : if (-IIn2t1>=puny) then
1077 0 : Hen2t1 = c1
1078 : else
1079 0 : Hen2t1 = c0
1080 : endif
1081 :
1082 0 : s12s0 = sign(c1,IIt1t2+puny)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12)
1083 0 : s21s0 = sign(c1,IIt1t2+puny)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21)
1084 0 : s12ks=p5*(s12s0+s21s0)
1085 :
1086 0 : end FUNCTION s12ks
1087 :
1088 : !=======================================================================
1089 : ! Function : s22ks
1090 :
1091 0 : FUNCTION s22ks(x,y,z,phi)
1092 :
1093 : real (kind=dbl_kind), intent(in) :: &
1094 : x,y,z,phi
1095 :
1096 : real (kind=dbl_kind) :: &
1097 0 : s22ks,p
1098 :
1099 : real (kind=dbl_kind) :: &
1100 : n1t2i11, n1t2i12, n1t2i21, n1t2i22, & ! LCOV_EXCL_LINE
1101 : n2t1i11, n2t1i12, n2t1i21, n2t1i22, & ! LCOV_EXCL_LINE
1102 : t1t2i11, t1t2i12, t1t2i21, t1t2i22, & ! LCOV_EXCL_LINE
1103 : ! t2t1i11, t2t1i12, t2t1i21, & ! LCOV_EXCL_LINE
1104 : t2t1i22, & ! LCOV_EXCL_LINE
1105 : d11, d12, d22, & ! LCOV_EXCL_LINE
1106 : IIn1t2, IIn2t1, IIt1t2, & ! LCOV_EXCL_LINE
1107 0 : Hen1t2, Hen2t1
1108 :
1109 : character(len=*), parameter :: subname = '(s22ks)'
1110 :
1111 0 : p = phi
1112 :
1113 0 : n1t2i11 = cos(z+pih-p) * cos(z+p)
1114 0 : n1t2i12 = cos(z+pih-p) * sin(z+p)
1115 0 : n1t2i21 = sin(z+pih-p) * cos(z+p)
1116 0 : n1t2i22 = sin(z+pih-p) * sin(z+p)
1117 0 : n2t1i11 = cos(z-pih+p) * cos(z-p)
1118 0 : n2t1i12 = cos(z-pih+p) * sin(z-p)
1119 0 : n2t1i21 = sin(z-pih+p) * cos(z-p)
1120 0 : n2t1i22 = sin(z-pih+p) * sin(z-p)
1121 0 : t1t2i11 = cos(z-p) * cos(z+p)
1122 0 : t1t2i12 = cos(z-p) * sin(z+p)
1123 0 : t1t2i21 = sin(z-p) * cos(z+p)
1124 0 : t1t2i22 = sin(z-p) * sin(z+p)
1125 : ! t2t1i11 = cos(z+p) * cos(z-p)
1126 : ! t2t1i12 = cos(z+p) * sin(z-p)
1127 : ! t2t1i21 = sin(z+p) * cos(z-p)
1128 0 : t2t1i22 = sin(z+p) * sin(z-p)
1129 0 : d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
1130 0 : d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
1131 0 : d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
1132 0 : IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
1133 0 : IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
1134 0 : IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
1135 :
1136 0 : if (-IIn1t2>=puny) then
1137 0 : Hen1t2 = c1
1138 : else
1139 0 : Hen1t2 = c0
1140 : endif
1141 :
1142 0 : if (-IIn2t1>=puny) then
1143 0 : Hen2t1 = c1
1144 : else
1145 0 : Hen2t1 = c0
1146 : endif
1147 :
1148 0 : s22ks = sign(c1,IIt1t2+puny)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22)
1149 :
1150 0 : end FUNCTION s22ks
1151 :
1152 : !=======================================================================
1153 : ! Computes the rates of strain and internal stress components for
1154 : ! each of the four corners on each T-grid cell.
1155 : ! Computes stress terms for the momentum equation
1156 : ! (based on subroutine stress)
1157 :
1158 0 : subroutine stress_eap (nx_block, ny_block, &
1159 : ksub, ndte, & ! LCOV_EXCL_LINE
1160 : icellT, & ! LCOV_EXCL_LINE
1161 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1162 : arlx1i, denom1, & ! LCOV_EXCL_LINE
1163 : uvel, vvel, & ! LCOV_EXCL_LINE
1164 : dxT, dyT, & ! LCOV_EXCL_LINE
1165 : dxhy, dyhx, & ! LCOV_EXCL_LINE
1166 : cxp, cyp, & ! LCOV_EXCL_LINE
1167 : cxm, cym, & ! LCOV_EXCL_LINE
1168 : tarear, strength, & ! LCOV_EXCL_LINE
1169 : a11_1, a11_2, a11_3, a11_4, & ! LCOV_EXCL_LINE
1170 : a12_1, a12_2, a12_3, a12_4, & ! LCOV_EXCL_LINE
1171 : stressp_1, stressp_2, & ! LCOV_EXCL_LINE
1172 : stressp_3, stressp_4, & ! LCOV_EXCL_LINE
1173 : stressm_1, stressm_2, & ! LCOV_EXCL_LINE
1174 : stressm_3, stressm_4, & ! LCOV_EXCL_LINE
1175 : stress12_1, stress12_2, & ! LCOV_EXCL_LINE
1176 : stress12_3, stress12_4, & ! LCOV_EXCL_LINE
1177 : shear, divu, & ! LCOV_EXCL_LINE
1178 : e11, e12, & ! LCOV_EXCL_LINE
1179 : e22, & ! LCOV_EXCL_LINE
1180 : s11, s12, & ! LCOV_EXCL_LINE
1181 : s22, & ! LCOV_EXCL_LINE
1182 : yieldstress11, & ! LCOV_EXCL_LINE
1183 : yieldstress12, & ! LCOV_EXCL_LINE
1184 : yieldstress22, & ! LCOV_EXCL_LINE
1185 : ! rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
1186 : rdg_conv, & ! LCOV_EXCL_LINE
1187 0 : strtmp)
1188 :
1189 : integer (kind=int_kind), intent(in) :: &
1190 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1191 : ksub , & ! subcycling step ! LCOV_EXCL_LINE
1192 : ndte , & ! number of subcycles ! LCOV_EXCL_LINE
1193 : icellT ! no. of cells where iceTmask = .true.
1194 :
1195 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1196 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1197 : indxTj ! compressed index in j-direction
1198 :
1199 : real (kind=dbl_kind), intent(in) :: &
1200 : arlx1i , & ! dte/2T (original) or 1/alpha1 (revised) ! LCOV_EXCL_LINE
1201 : denom1 ! constant for stress equation
1202 :
1203 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1204 : strength , & ! ice strength (N/m) ! LCOV_EXCL_LINE
1205 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1206 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1207 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1208 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1209 : dxhy , & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
1210 : dyhx , & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
1211 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1212 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1213 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1214 : cxm , & ! 0.5*HTN - 1.5*HTS ! LCOV_EXCL_LINE
1215 : tarear ! 1/tarea
1216 :
1217 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1218 : stressp_1, stressp_2, stressp_3, stressp_4, & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1219 : stressm_1, stressm_2, stressm_3, stressm_4, & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1220 : stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12
1221 :
1222 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1223 : a11_1, a11_2, a11_3, a11_4, & ! structure tensor ! LCOV_EXCL_LINE
1224 : a12_1, a12_2, a12_3, a12_4 ! structure tensor
1225 :
1226 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1227 : shear , & ! strain rate II component (1/s) ! LCOV_EXCL_LINE
1228 : divu , & ! strain rate I component, velocity divergence (1/s) ! LCOV_EXCL_LINE
1229 : e11 , & ! components of strain rate tensor (1/s) ! LCOV_EXCL_LINE
1230 : e12 , & ! ! LCOV_EXCL_LINE
1231 : e22 , & ! ! LCOV_EXCL_LINE
1232 : s11 , & ! components of stress tensor (kg/s^2) ! LCOV_EXCL_LINE
1233 : s12 , & ! ! LCOV_EXCL_LINE
1234 : s22 , & ! ! LCOV_EXCL_LINE
1235 : yieldstress11, & ! components of yield stress tensor (kg/s^2) ! LCOV_EXCL_LINE
1236 : yieldstress12, & ! LCOV_EXCL_LINE
1237 : yieldstress22, & ! LCOV_EXCL_LINE
1238 : rdg_conv ! convergence term for ridging (1/s)
1239 : ! rdg_shear ! shear term for ridging (1/s)
1240 :
1241 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: &
1242 : strtmp ! stress combinations
1243 :
1244 : ! local variables
1245 :
1246 : integer (kind=int_kind) :: &
1247 : i, j, ij
1248 :
1249 : real (kind=dbl_kind) :: &
1250 : stressptmp_1, stressptmp_2, stressptmp_3, stressptmp_4, & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1251 : stressmtmp_1, stressmtmp_2, stressmtmp_3, stressmtmp_4, & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1252 0 : stress12tmp_1,stress12tmp_2,stress12tmp_3,stress12tmp_4 ! sigma12
1253 :
1254 : real (kind=dbl_kind) :: &
1255 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1256 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1257 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1258 : ssigpn, ssigps, ssigpe, ssigpw , & ! LCOV_EXCL_LINE
1259 : ssigmn, ssigms, ssigme, ssigmw , & ! LCOV_EXCL_LINE
1260 : ssig12n, ssig12s, ssig12e, ssig12w , & ! LCOV_EXCL_LINE
1261 : ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122, & ! LCOV_EXCL_LINE
1262 : csigpne, csigpnw, csigpse, csigpsw , & ! LCOV_EXCL_LINE
1263 : csigmne, csigmnw, csigmse, csigmsw , & ! LCOV_EXCL_LINE
1264 : csig12ne, csig12nw, csig12se, csig12sw , & ! LCOV_EXCL_LINE
1265 : str12ew, str12we, str12ns, str12sn , & ! LCOV_EXCL_LINE
1266 0 : strp_tmp, strm_tmp
1267 :
1268 : real (kind=dbl_kind) :: &
1269 : alpharne, alpharnw, alpharsw, alpharse, & ! LCOV_EXCL_LINE
1270 0 : alphasne, alphasnw, alphassw, alphasse
1271 :
1272 : character(len=*), parameter :: subname = '(stress_eap)'
1273 :
1274 : !-----------------------------------------------------------------
1275 : ! Initialize
1276 : !-----------------------------------------------------------------
1277 :
1278 0 : strtmp(:,:,:) = c0
1279 :
1280 0 : do ij = 1, icellT
1281 0 : i = indxTi(ij)
1282 0 : j = indxTj(ij)
1283 :
1284 : !-----------------------------------------------------------------
1285 : ! strain rates
1286 : ! NOTE these are actually strain rates * area (m^2/s)
1287 : !-----------------------------------------------------------------
1288 :
1289 : ! divergence = e_11 + e_22
1290 0 : divune = cyp(i,j)*uvel(i ,j ) - dyT(i,j)*uvel(i-1,j ) &
1291 0 : + cxp(i,j)*vvel(i ,j ) - dxT(i,j)*vvel(i ,j-1)
1292 0 : divunw = cym(i,j)*uvel(i-1,j ) + dyT(i,j)*uvel(i ,j ) &
1293 0 : + cxp(i,j)*vvel(i-1,j ) - dxT(i,j)*vvel(i-1,j-1)
1294 0 : divusw = cym(i,j)*uvel(i-1,j-1) + dyT(i,j)*uvel(i ,j-1) &
1295 0 : + cxm(i,j)*vvel(i-1,j-1) + dxT(i,j)*vvel(i-1,j )
1296 0 : divuse = cyp(i,j)*uvel(i ,j-1) - dyT(i,j)*uvel(i-1,j-1) &
1297 0 : + cxm(i,j)*vvel(i ,j-1) + dxT(i,j)*vvel(i ,j )
1298 :
1299 : ! tension strain rate = e_11 - e_22
1300 0 : tensionne = -cym(i,j)*uvel(i ,j ) - dyT(i,j)*uvel(i-1,j ) &
1301 0 : + cxm(i,j)*vvel(i ,j ) + dxT(i,j)*vvel(i ,j-1)
1302 0 : tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyT(i,j)*uvel(i ,j ) &
1303 0 : + cxm(i,j)*vvel(i-1,j ) + dxT(i,j)*vvel(i-1,j-1)
1304 0 : tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyT(i,j)*uvel(i ,j-1) &
1305 0 : + cxp(i,j)*vvel(i-1,j-1) - dxT(i,j)*vvel(i-1,j )
1306 0 : tensionse = -cym(i,j)*uvel(i ,j-1) - dyT(i,j)*uvel(i-1,j-1) &
1307 0 : + cxp(i,j)*vvel(i ,j-1) - dxT(i,j)*vvel(i ,j )
1308 :
1309 : ! shearing strain rate = 2*e_12
1310 0 : shearne = -cym(i,j)*vvel(i ,j ) - dyT(i,j)*vvel(i-1,j ) &
1311 0 : - cxm(i,j)*uvel(i ,j ) - dxT(i,j)*uvel(i ,j-1)
1312 0 : shearnw = -cyp(i,j)*vvel(i-1,j ) + dyT(i,j)*vvel(i ,j ) &
1313 0 : - cxm(i,j)*uvel(i-1,j ) - dxT(i,j)*uvel(i-1,j-1)
1314 0 : shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyT(i,j)*vvel(i ,j-1) &
1315 0 : - cxp(i,j)*uvel(i-1,j-1) + dxT(i,j)*uvel(i-1,j )
1316 0 : shearse = -cym(i,j)*vvel(i ,j-1) - dyT(i,j)*vvel(i-1,j-1) &
1317 0 : - cxp(i,j)*uvel(i ,j-1) + dxT(i,j)*uvel(i ,j )
1318 :
1319 : !-----------------------------------------------------------------
1320 : ! Stress updated depending on strain rate and structure tensor
1321 : !-----------------------------------------------------------------
1322 :
1323 : ! ne
1324 : call update_stress_rdg (ksub, ndte, divune, tensionne, &
1325 : shearne, a11_1(i,j), a12_1(i,j), & ! LCOV_EXCL_LINE
1326 : stressptmp_1, stressmtmp_1, & ! LCOV_EXCL_LINE
1327 : stress12tmp_1, strength(i,j), & ! LCOV_EXCL_LINE
1328 0 : alpharne, alphasne)
1329 : ! nw
1330 : call update_stress_rdg (ksub, ndte, divunw, tensionnw, &
1331 : shearnw, a11_2(i,j), a12_2(i,j), & ! LCOV_EXCL_LINE
1332 : stressptmp_2, stressmtmp_2, & ! LCOV_EXCL_LINE
1333 : stress12tmp_2, strength(i,j), & ! LCOV_EXCL_LINE
1334 0 : alpharnw, alphasnw)
1335 : ! sw
1336 : call update_stress_rdg (ksub, ndte, divusw, tensionsw, &
1337 : shearsw, a11_3(i,j), a12_3(i,j), & ! LCOV_EXCL_LINE
1338 : stressptmp_3, stressmtmp_3, & ! LCOV_EXCL_LINE
1339 : stress12tmp_3, strength(i,j), & ! LCOV_EXCL_LINE
1340 0 : alpharsw, alphassw)
1341 : ! se
1342 : call update_stress_rdg (ksub, ndte, divuse, tensionse, &
1343 : shearse, a11_4(i,j), a12_4(i,j), & ! LCOV_EXCL_LINE
1344 : stressptmp_4, stressmtmp_4, & ! LCOV_EXCL_LINE
1345 : stress12tmp_4, strength(i,j), & ! LCOV_EXCL_LINE
1346 0 : alpharse, alphasse)
1347 :
1348 : !-----------------------------------------------------------------
1349 : ! on last subcycle, save quantities for mechanical redistribution
1350 : !-----------------------------------------------------------------
1351 :
1352 0 : if (ksub == ndte) then
1353 :
1354 : ! diagnostic only
1355 : ! shear = sqrt(tension**2 + shearing**2)
1356 0 : shear(i,j) = p25*tarear(i,j)*sqrt( &
1357 : (tensionne + tensionnw + tensionse + tensionsw)**2 & ! LCOV_EXCL_LINE
1358 0 : + (shearne + shearnw + shearse + shearsw)**2)
1359 :
1360 0 : divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
1361 0 : rdg_conv(i,j) = -min(p25*(alpharne + alpharnw &
1362 0 : + alpharsw + alpharse),c0) * tarear(i,j)
1363 : !rdg_shear=0 for computing closing_net in ridge_prep
1364 : !rdg_shear(i,j) = p25*(alphasne + alphasnw &
1365 : ! + alphassw + alphasse) * tarear(i,j)
1366 : endif
1367 :
1368 0 : e11(i,j) = p5*p25*(divune + divunw + divuse + divusw + &
1369 0 : tensionne + tensionnw + tensionse + tensionsw) * tarear(i,j)
1370 :
1371 0 : e12(i,j) = p5*p25*(shearne + shearnw + shearse + shearsw) * tarear(i,j)
1372 :
1373 0 : e22(i,j) = p5*p25*(divune + divunw + divuse + divusw - &
1374 0 : tensionne - tensionnw - tensionse - tensionsw) * tarear(i,j)
1375 :
1376 : !-----------------------------------------------------------------
1377 : ! elastic relaxation, see Eq. A12-A14
1378 : !-----------------------------------------------------------------
1379 :
1380 0 : stressp_1(i,j) = (stressp_1(i,j) + stressptmp_1*arlx1i) &
1381 0 : * denom1
1382 0 : stressp_2(i,j) = (stressp_2(i,j) + stressptmp_2*arlx1i) &
1383 0 : * denom1
1384 0 : stressp_3(i,j) = (stressp_3(i,j) + stressptmp_3*arlx1i) &
1385 0 : * denom1
1386 0 : stressp_4(i,j) = (stressp_4(i,j) + stressptmp_4*arlx1i) &
1387 0 : * denom1
1388 :
1389 0 : stressm_1(i,j) = (stressm_1(i,j) + stressmtmp_1*arlx1i) &
1390 0 : * denom1
1391 0 : stressm_2(i,j) = (stressm_2(i,j) + stressmtmp_2*arlx1i) &
1392 0 : * denom1
1393 0 : stressm_3(i,j) = (stressm_3(i,j) + stressmtmp_3*arlx1i) &
1394 0 : * denom1
1395 0 : stressm_4(i,j) = (stressm_4(i,j) + stressmtmp_4*arlx1i) &
1396 0 : * denom1
1397 :
1398 0 : stress12_1(i,j) = (stress12_1(i,j) + stress12tmp_1*arlx1i) &
1399 0 : * denom1
1400 0 : stress12_2(i,j) = (stress12_2(i,j) + stress12tmp_2*arlx1i) &
1401 0 : * denom1
1402 0 : stress12_3(i,j) = (stress12_3(i,j) + stress12tmp_3*arlx1i) &
1403 0 : * denom1
1404 0 : stress12_4(i,j) = (stress12_4(i,j) + stress12tmp_4*arlx1i) &
1405 0 : * denom1
1406 :
1407 0 : s11(i,j) = p5 * p25 * (stressp_1 (i,j) + stressp_2 (i,j) &
1408 : + stressp_3 (i,j) + stressp_4 (i,j) & ! LCOV_EXCL_LINE
1409 : + stressm_1 (i,j) + stressm_2 (i,j) & ! LCOV_EXCL_LINE
1410 0 : + stressm_3 (i,j) + stressm_4 (i,j))
1411 0 : s22(i,j) = p5 * p25 * (stressp_1 (i,j) + stressp_2 (i,j) &
1412 : + stressp_3 (i,j) + stressp_4 (i,j) & ! LCOV_EXCL_LINE
1413 : - stressm_1 (i,j) - stressm_2 (i,j) & ! LCOV_EXCL_LINE
1414 0 : - stressm_3 (i,j) - stressm_4 (i,j))
1415 0 : s12(i,j) = p25 * (stress12_1(i,j) + stress12_2(i,j) &
1416 0 : + stress12_3(i,j) + stress12_4(i,j))
1417 :
1418 0 : yieldstress11(i,j) = p5 * p25 * (stressptmp_1 + stressptmp_2 &
1419 : + stressptmp_3 + stressptmp_4 & ! LCOV_EXCL_LINE
1420 : + stressmtmp_1 + stressmtmp_2 & ! LCOV_EXCL_LINE
1421 0 : + stressmtmp_3 + stressmtmp_4)
1422 0 : yieldstress22(i,j) = p5 * p25 * (stressptmp_1 + stressptmp_2 &
1423 : + stressptmp_3 + stressptmp_4 & ! LCOV_EXCL_LINE
1424 : - stressmtmp_1 - stressmtmp_2 & ! LCOV_EXCL_LINE
1425 0 : - stressmtmp_3 - stressmtmp_4)
1426 0 : yieldstress12(i,j) = p25 * (stress12tmp_1 + stress12tmp_2 &
1427 0 : + stress12tmp_3 + stress12tmp_4)
1428 :
1429 : !-----------------------------------------------------------------
1430 : ! Eliminate underflows.
1431 : ! The following code is commented out because it is relatively
1432 : ! expensive and most compilers include a flag that accomplishes
1433 : ! the same thing more efficiently. This code is cheaper than
1434 : ! handling underflows if the compiler lacks a flag; uncomment
1435 : ! it in that case. The compiler flag is often described with the
1436 : ! phrase "flush to zero".
1437 : !-----------------------------------------------------------------
1438 :
1439 : ! stressp_1(i,j) = sign(max(abs(stressp_1(i,j)),puny),stressp_1(i,j))
1440 : ! stressp_2(i,j) = sign(max(abs(stressp_2(i,j)),puny),stressp_2(i,j))
1441 : ! stressp_3(i,j) = sign(max(abs(stressp_3(i,j)),puny),stressp_3(i,j))
1442 : ! stressp_4(i,j) = sign(max(abs(stressp_4(i,j)),puny),stressp_4(i,j))
1443 :
1444 : ! stressm_1(i,j) = sign(max(abs(stressm_1(i,j)),puny),stressm_1(i,j))
1445 : ! stressm_2(i,j) = sign(max(abs(stressm_2(i,j)),puny),stressm_2(i,j))
1446 : ! stressm_3(i,j) = sign(max(abs(stressm_3(i,j)),puny),stressm_3(i,j))
1447 : ! stressm_4(i,j) = sign(max(abs(stressm_4(i,j)),puny),stressm_4(i,j))
1448 :
1449 : ! stress12_1(i,j) = sign(max(abs(stress12_1(i,j)),puny),stress12_1(i,j))
1450 : ! stress12_2(i,j) = sign(max(abs(stress12_2(i,j)),puny),stress12_2(i,j))
1451 : ! stress12_3(i,j) = sign(max(abs(stress12_3(i,j)),puny),stress12_3(i,j))
1452 : ! stress12_4(i,j) = sign(max(abs(stress12_4(i,j)),puny),stress12_4(i,j))
1453 :
1454 : !-----------------------------------------------------------------
1455 : ! combinations of the stresses for the momentum equation ! kg/s^2
1456 : !-----------------------------------------------------------------
1457 :
1458 0 : ssigpn = stressp_1(i,j) + stressp_2(i,j)
1459 0 : ssigps = stressp_3(i,j) + stressp_4(i,j)
1460 0 : ssigpe = stressp_1(i,j) + stressp_4(i,j)
1461 0 : ssigpw = stressp_2(i,j) + stressp_3(i,j)
1462 0 : ssigp1 =(stressp_1(i,j) + stressp_3(i,j))*p055
1463 0 : ssigp2 =(stressp_2(i,j) + stressp_4(i,j))*p055
1464 :
1465 0 : ssigmn = stressm_1(i,j) + stressm_2(i,j)
1466 0 : ssigms = stressm_3(i,j) + stressm_4(i,j)
1467 0 : ssigme = stressm_1(i,j) + stressm_4(i,j)
1468 0 : ssigmw = stressm_2(i,j) + stressm_3(i,j)
1469 0 : ssigm1 =(stressm_1(i,j) + stressm_3(i,j))*p055
1470 0 : ssigm2 =(stressm_2(i,j) + stressm_4(i,j))*p055
1471 :
1472 0 : ssig12n = stress12_1(i,j) + stress12_2(i,j)
1473 0 : ssig12s = stress12_3(i,j) + stress12_4(i,j)
1474 0 : ssig12e = stress12_1(i,j) + stress12_4(i,j)
1475 0 : ssig12w = stress12_2(i,j) + stress12_3(i,j)
1476 0 : ssig121 =(stress12_1(i,j) + stress12_3(i,j))*p111
1477 0 : ssig122 =(stress12_2(i,j) + stress12_4(i,j))*p111
1478 :
1479 0 : csigpne = p111*stressp_1(i,j) + ssigp2 + p027*stressp_3(i,j)
1480 0 : csigpnw = p111*stressp_2(i,j) + ssigp1 + p027*stressp_4(i,j)
1481 0 : csigpsw = p111*stressp_3(i,j) + ssigp2 + p027*stressp_1(i,j)
1482 0 : csigpse = p111*stressp_4(i,j) + ssigp1 + p027*stressp_2(i,j)
1483 :
1484 0 : csigmne = p111*stressm_1(i,j) + ssigm2 + p027*stressm_3(i,j)
1485 0 : csigmnw = p111*stressm_2(i,j) + ssigm1 + p027*stressm_4(i,j)
1486 0 : csigmsw = p111*stressm_3(i,j) + ssigm2 + p027*stressm_1(i,j)
1487 0 : csigmse = p111*stressm_4(i,j) + ssigm1 + p027*stressm_2(i,j)
1488 :
1489 0 : csig12ne = p222*stress12_1(i,j) + ssig122 &
1490 0 : + p055*stress12_3(i,j)
1491 0 : csig12nw = p222*stress12_2(i,j) + ssig121 &
1492 0 : + p055*stress12_4(i,j)
1493 0 : csig12sw = p222*stress12_3(i,j) + ssig122 &
1494 0 : + p055*stress12_1(i,j)
1495 0 : csig12se = p222*stress12_4(i,j) + ssig121 &
1496 0 : + p055*stress12_2(i,j)
1497 :
1498 0 : str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
1499 0 : str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
1500 0 : str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
1501 0 : str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)
1502 :
1503 : !-----------------------------------------------------------------
1504 : ! for dF/dx (u momentum)
1505 : !-----------------------------------------------------------------
1506 :
1507 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
1508 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)
1509 :
1510 : ! northeast (i,j)
1511 0 : strtmp(i,j,1) = -strp_tmp - strm_tmp - str12ew &
1512 0 : + dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne
1513 :
1514 : ! northwest (i+1,j)
1515 0 : strtmp(i,j,2) = strp_tmp + strm_tmp - str12we &
1516 0 : + dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw
1517 :
1518 0 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
1519 0 : strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)
1520 :
1521 : ! southeast (i,j+1)
1522 0 : strtmp(i,j,3) = -strp_tmp - strm_tmp + str12ew &
1523 0 : + dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se
1524 :
1525 : ! southwest (i+1,j+1)
1526 0 : strtmp(i,j,4) = strp_tmp + strm_tmp + str12we &
1527 0 : + dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw
1528 :
1529 : !-----------------------------------------------------------------
1530 : ! for dF/dy (v momentum)
1531 : !-----------------------------------------------------------------
1532 :
1533 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
1534 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)
1535 :
1536 : ! northeast (i,j)
1537 0 : strtmp(i,j,5) = -strp_tmp + strm_tmp - str12ns &
1538 0 : - dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne
1539 :
1540 : ! southeast (i,j+1)
1541 0 : strtmp(i,j,6) = strp_tmp - strm_tmp - str12sn &
1542 0 : - dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se
1543 :
1544 0 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
1545 0 : strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)
1546 :
1547 : ! northwest (i+1,j)
1548 0 : strtmp(i,j,7) = -strp_tmp + strm_tmp + str12ns &
1549 0 : - dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw
1550 :
1551 : ! southwest (i+1,j+1)
1552 0 : strtmp(i,j,8) = strp_tmp - strm_tmp + str12sn &
1553 0 : - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw
1554 :
1555 : enddo ! ij
1556 :
1557 0 : end subroutine stress_eap
1558 :
1559 : !=======================================================================
1560 : ! Updates the stress depending on values of strain rate and structure
1561 : ! tensor and for ksub=ndte it computes closing and sliding rate
1562 :
1563 0 : subroutine update_stress_rdg (ksub, ndte, divu, tension, &
1564 : shear, a11, a12, & ! LCOV_EXCL_LINE
1565 : stressp, stressm, & ! LCOV_EXCL_LINE
1566 : stress12, strength, & ! LCOV_EXCL_LINE
1567 : alphar, alphas)
1568 :
1569 : integer (kind=int_kind), intent(in) :: &
1570 : ksub, & ! LCOV_EXCL_LINE
1571 : ndte
1572 :
1573 : real (kind=dbl_kind), intent(in) :: &
1574 : a11, a12, & ! LCOV_EXCL_LINE
1575 : divu, tension, shear, & ! LCOV_EXCL_LINE
1576 : strength
1577 :
1578 : real (kind=dbl_kind), intent(out) :: &
1579 : stressp, stressm, stress12, & ! LCOV_EXCL_LINE
1580 : alphar, alphas
1581 :
1582 : ! local variables
1583 :
1584 : integer (kind=int_kind) :: &
1585 : kx ,ky, ka
1586 :
1587 : real (kind=dbl_kind) :: &
1588 : stemp11r, stemp12r, stemp22r, & ! LCOV_EXCL_LINE
1589 : stemp11s, stemp12s, stemp22s, & ! LCOV_EXCL_LINE
1590 : a22, Qd11Qd11, Qd11Qd12, Qd12Qd12, & ! LCOV_EXCL_LINE
1591 : Q11Q11, Q11Q12, Q12Q12, & ! LCOV_EXCL_LINE
1592 : dtemp11, dtemp12, dtemp22, & ! LCOV_EXCL_LINE
1593 : rotstemp11r, rotstemp12r, rotstemp22r, & ! LCOV_EXCL_LINE
1594 : rotstemp11s, rotstemp12s, rotstemp22s, & ! LCOV_EXCL_LINE
1595 : sig11, sig12, sig22, & ! LCOV_EXCL_LINE
1596 : sgprm11, sgprm12, sgprm22, & ! LCOV_EXCL_LINE
1597 : Angle_denom_gamma, Angle_denom_alpha, & ! LCOV_EXCL_LINE
1598 : Tany_1, Tany_2, & ! LCOV_EXCL_LINE
1599 : x, y, dx, dy, da, & ! LCOV_EXCL_LINE
1600 : dtemp1, dtemp2, atempprime, & ! LCOV_EXCL_LINE
1601 0 : kxw, kyw, kaw
1602 :
1603 : ! tcraig, temporary, should be moved to namelist
1604 : ! turns on interpolation in stress_rdg
1605 : logical(kind=log_kind), parameter :: &
1606 : interpolate_stress_rdg = .false.
1607 :
1608 : character(len=*), parameter :: subname = '(update_stress_rdg)'
1609 :
1610 : ! compute eigenvalues, eigenvectors and angles for structure tensor, strain rates
1611 :
1612 : ! 1) structure tensor
1613 :
1614 0 : a22 = c1-a11
1615 :
1616 : ! gamma: angle between general coordinates and principal axis of A
1617 : ! here Tan2gamma = 2 a12 / (a11 - a22)
1618 :
1619 0 : Q11Q11 = c1
1620 0 : Q12Q12 = puny
1621 0 : Q11Q12 = puny
1622 :
1623 0 : if ((ABS(a11 - a22) > puny).or.(ABS(a12) > puny)) then
1624 : Angle_denom_gamma = sqrt( ( a11 - a22 )*( a11 - a22) + &
1625 0 : c4*a12*a12 )
1626 :
1627 0 : Q11Q11 = p5 + ( a11 - a22 )*p5/Angle_denom_gamma !Cos^2
1628 0 : Q12Q12 = p5 - ( a11 - a22 )*p5/Angle_denom_gamma !Sin^2
1629 0 : Q11Q12 = a12/Angle_denom_gamma !CosSin
1630 : endif
1631 :
1632 : ! rotation Q*atemp*Q^T
1633 0 : atempprime = Q11Q11*a11 + c2*Q11Q12*a12 + Q12Q12*a22
1634 :
1635 : ! make first principal value the largest
1636 0 : atempprime = max(atempprime, c1 - atempprime)
1637 :
1638 : ! 2) strain rate
1639 :
1640 0 : dtemp11 = p5*(divu + tension)
1641 0 : dtemp12 = shear*p5
1642 0 : dtemp22 = p5*(divu - tension)
1643 :
1644 : ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)
1645 :
1646 0 : Qd11Qd11 = c1
1647 0 : Qd12Qd12 = puny
1648 0 : Qd11Qd12 = puny
1649 :
1650 0 : if ((ABS( dtemp11 - dtemp22) > puny).or.(ABS(dtemp12) > puny)) then
1651 : Angle_denom_alpha = sqrt( ( dtemp11 - dtemp22 )* &
1652 0 : ( dtemp11 - dtemp22 ) + c4*dtemp12*dtemp12)
1653 :
1654 0 : Qd11Qd11 = p5 + ( dtemp11 - dtemp22 )*p5/Angle_denom_alpha !Cos^2
1655 0 : Qd12Qd12 = p5 - ( dtemp11 - dtemp22 )*p5/Angle_denom_alpha !Sin^2
1656 0 : Qd11Qd12 = dtemp12/Angle_denom_alpha !CosSin
1657 : endif
1658 :
1659 0 : dtemp1 = Qd11Qd11*dtemp11 + c2*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22
1660 0 : dtemp2 = Qd12Qd12*dtemp11 - c2*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22
1661 :
1662 : ! In cos and sin values
1663 0 : x = c0
1664 :
1665 0 : if ((ABS(dtemp1) > puny).or.(ABS(dtemp2) > puny)) then
1666 : ! invleng = c1/sqrt(dtemp1*dtemp1 + dtemp2*dtemp2) ! not sure if this is neccessary
1667 : ! dtemp1 = dtemp1*invleng
1668 : ! dtemp2 = dtemp2*invleng
1669 0 : if (dtemp1 == c0) then
1670 0 : x = pih
1671 : else
1672 0 : x = atan2(dtemp2,dtemp1)
1673 : endif
1674 : endif
1675 :
1676 : !echmod to ensure the angle lies between pi/4 and 9 pi/4
1677 0 : if (x < piq) x = x + pi2
1678 : !echmod require 0 <= x < (nx_yield-1)*dx = 2 pi
1679 : ! x = mod(x+pi2, pi2)
1680 : ! y: angle between major principal axis of strain rate and structure tensor
1681 : ! y = gamma - alpha
1682 : ! Expressesed componently with
1683 : ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
1684 :
1685 0 : Tany_1 = Q11Q12 - Qd11Qd12
1686 0 : Tany_2 = Q11Q11 - Qd12Qd12
1687 :
1688 0 : y = c0
1689 :
1690 0 : if ((ABS(Tany_1) > puny).or.(ABS(Tany_2) > puny)) then
1691 : ! invleng = c1/sqrt(Tany_1*Tany_1 + Tany_2*Tany_2) ! not sure if this is neccessary
1692 : ! Tany_1 = Tany_1*invleng
1693 : ! Tany_2 = Tany_2*invleng
1694 0 : if (Tany_2 == c0) then
1695 0 : y = pih
1696 : else
1697 0 : y = atan2(Tany_1,Tany_2)
1698 : endif
1699 : endif
1700 :
1701 : ! to make sure y is between 0 and pi
1702 :
1703 0 : if (y > pi) y = y - pi
1704 0 : if (y < 0) y = y + pi
1705 :
1706 : if (interpolate_stress_rdg) then
1707 :
1708 : ! Interpolated lookup
1709 :
1710 : ! if (x>=9*pi/4) x=9*pi/4-puny; end
1711 : ! if (y>=pi/2) y=pi/2-puny; end
1712 : ! if (atempprime>=1.0), atempprime=1.0-puny; end
1713 :
1714 : ! % need 8 coords and 8 weights
1715 : ! % range in kx
1716 :
1717 : kx = int((x-piq-pi)*invdx) + 1
1718 : kxw = c1 - ((x-piq-pi)*invdx - (kx-1))
1719 :
1720 : ky = int(y*invdy) + 1
1721 : kyw = c1 - (y*invdy - (ky-1))
1722 :
1723 : ka = int((atempprime-p5)*invda) + 1
1724 : kaw = c1 - ((atempprime-p5)*invda - (ka-1))
1725 :
1726 : ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1)
1727 :
1728 : stemp11r = kxw* kyw * kaw * s11r(kx ,ky ,ka ) &
1729 : + (c1-kxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1730 : + kxw * (c1-kyw) * kaw * s11r(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1731 : + kxw * kyw * (c1-kaw) * s11r(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1732 : + (c1-kxw) * (c1-kyw) * kaw * s11r(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1733 : + (c1-kxw) * kyw * (c1-kaw) * s11r(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1734 : + kxw * (c1-kyw) * (c1-kaw) * s11r(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1735 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s11r(kx+1,ky+1,ka+1)
1736 :
1737 : stemp12r = kxw* kyw * kaw * s12r(kx ,ky ,ka ) &
1738 : + (c1-kxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1739 : + kxw * (c1-kyw) * kaw * s12r(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1740 : + kxw * kyw * (c1-kaw) * s12r(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1741 : + (c1-kxw) * (c1-kyw) * kaw * s12r(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1742 : + (c1-kxw) * kyw * (c1-kaw) * s12r(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1743 : + kxw * (c1-kyw) * (c1-kaw) * s12r(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1744 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s12r(kx+1,ky+1,ka+1)
1745 :
1746 : stemp22r = kxw * kyw * kaw * s22r(kx ,ky ,ka ) &
1747 : + (c1-kxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1748 : + kxw * (c1-kyw) * kaw * s22r(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1749 : + kxw * kyw * (c1-kaw) * s22r(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1750 : + (c1-kxw) * (c1-kyw) * kaw * s22r(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1751 : + (c1-kxw) * kyw * (c1-kaw) * s22r(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1752 : + kxw * (c1-kyw) * (c1-kaw) * s22r(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1753 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s22r(kx+1,ky+1,ka+1)
1754 :
1755 : stemp11s = kxw* kyw * kaw * s11s(kx ,ky ,ka ) &
1756 : + (c1-kxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1757 : + kxw * (c1-kyw) * kaw * s11s(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1758 : + kxw * kyw * (c1-kaw) * s11s(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1759 : + (c1-kxw) * (c1-kyw) * kaw * s11s(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1760 : + (c1-kxw) * kyw * (c1-kaw) * s11s(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1761 : + kxw * (c1-kyw) * (c1-kaw) * s11s(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1762 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s11s(kx+1,ky+1,ka+1)
1763 :
1764 : stemp12s = kxw* kyw * kaw * s12s(kx ,ky ,ka ) &
1765 : + (c1-kxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1766 : + kxw * (c1-kyw) * kaw * s12s(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1767 : + kxw * kyw * (c1-kaw) * s12s(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1768 : + (c1-kxw) * (c1-kyw) * kaw * s12s(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1769 : + (c1-kxw) * kyw * (c1-kaw) * s12s(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1770 : + kxw * (c1-kyw) * (c1-kaw) * s12s(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1771 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s12s(kx+1,ky+1,ka+1)
1772 :
1773 : stemp22s = kxw* kyw * kaw * s22s(kx ,ky ,ka ) &
1774 : + (c1-kxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & ! LCOV_EXCL_LINE
1775 : + kxw * (c1-kyw) * kaw * s22s(kx ,ky+1,ka ) & ! LCOV_EXCL_LINE
1776 : + kxw * kyw * (c1-kaw) * s22s(kx ,ky ,ka+1) & ! LCOV_EXCL_LINE
1777 : + (c1-kxw) * (c1-kyw) * kaw * s22s(kx+1,ky+1,ka ) & ! LCOV_EXCL_LINE
1778 : + (c1-kxw) * kyw * (c1-kaw) * s22s(kx+1,ky ,ka+1) & ! LCOV_EXCL_LINE
1779 : + kxw * (c1-kyw) * (c1-kaw) * s22s(kx ,ky+1,ka+1) & ! LCOV_EXCL_LINE
1780 : + (c1-kxw) * (c1-kyw) * (c1-kaw) * s22s(kx+1,ky+1,ka+1)
1781 :
1782 : else
1783 :
1784 0 : kx = int((x-piq-pi)*invdx) + 1
1785 0 : ky = int(y*invdy) + 1
1786 0 : ka = int((atempprime-p5)*invda) + 1
1787 :
1788 : ! Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1)
1789 :
1790 0 : stemp11r = s11r(kx,ky,ka)
1791 0 : stemp12r = s12r(kx,ky,ka)
1792 0 : stemp22r = s22r(kx,ky,ka)
1793 :
1794 0 : stemp11s = s11s(kx,ky,ka)
1795 0 : stemp12s = s12s(kx,ky,ka)
1796 0 : stemp22s = s22s(kx,ky,ka)
1797 :
1798 : endif
1799 :
1800 : ! Calculate mean ice stress over a collection of floes (Equation 3)
1801 :
1802 : stressp = strength*(stemp11r + kfriction*stemp11s &
1803 0 : + stemp22r + kfriction*stemp22s) * invsin
1804 0 : stress12 = strength*(stemp12r + kfriction*stemp12s) * invsin
1805 : stressm = strength*(stemp11r + kfriction*stemp11s &
1806 0 : - stemp22r - kfriction*stemp22s) * invsin
1807 :
1808 : ! Back - rotation of the stress from principal axes into general coordinates
1809 :
1810 : ! Update stress
1811 :
1812 0 : sig11 = p5*(stressp + stressm)
1813 0 : sig12 = stress12
1814 0 : sig22 = p5*(stressp - stressm)
1815 :
1816 0 : sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 - c2*Q11Q12 *sig12
1817 0 : sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12
1818 0 : sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 + c2*Q11Q12 *sig12
1819 :
1820 0 : stressp = sgprm11 + sgprm22
1821 0 : stress12 = sgprm12
1822 0 : stressm = sgprm11 - sgprm22
1823 :
1824 : ! Compute ridging and sliding functions in general coordinates (Equation 11)
1825 :
1826 0 : if (ksub == ndte) then
1827 : rotstemp11r = Q11Q11*stemp11r - c2*Q11Q12* stemp12r &
1828 0 : + Q12Q12*stemp22r
1829 : rotstemp12r = Q11Q11*stemp12r + Q11Q12*(stemp11r-stemp22r) &
1830 0 : - Q12Q12*stemp12r
1831 : rotstemp22r = Q12Q12*stemp11r + c2*Q11Q12* stemp12r &
1832 0 : + Q11Q11*stemp22r
1833 :
1834 : rotstemp11s = Q11Q11*stemp11s - c2*Q11Q12* stemp12s &
1835 0 : + Q12Q12*stemp22s
1836 : rotstemp12s = Q11Q11*stemp12s + Q11Q12*(stemp11s-stemp22s) &
1837 0 : - Q12Q12*stemp12s
1838 : rotstemp22s = Q12Q12*stemp11s + c2*Q11Q12* stemp12s &
1839 0 : + Q11Q11*stemp22s
1840 :
1841 : alphar = rotstemp11r*dtemp11 + c2*rotstemp12r*dtemp12 &
1842 0 : + rotstemp22r*dtemp22
1843 : alphas = rotstemp11s*dtemp11 + c2*rotstemp12s*dtemp12 &
1844 0 : + rotstemp22s*dtemp22
1845 : endif
1846 :
1847 0 : end subroutine update_stress_rdg
1848 :
1849 : !=======================================================================
1850 : ! Solves evolution equation for structure tensor (A19, A20)
1851 :
1852 0 : subroutine stepa (nx_block, ny_block, &
1853 : dtei, icellT, & ! LCOV_EXCL_LINE
1854 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1855 : a11, a12, & ! LCOV_EXCL_LINE
1856 : a11_1, a11_2, a11_3, a11_4, & ! LCOV_EXCL_LINE
1857 : a12_1, a12_2, a12_3, a12_4, & ! LCOV_EXCL_LINE
1858 : stressp_1, stressp_2, & ! LCOV_EXCL_LINE
1859 : stressp_3, stressp_4, & ! LCOV_EXCL_LINE
1860 : stressm_1, stressm_2, & ! LCOV_EXCL_LINE
1861 : stressm_3, stressm_4, & ! LCOV_EXCL_LINE
1862 : stress12_1, stress12_2, & ! LCOV_EXCL_LINE
1863 0 : stress12_3, stress12_4)
1864 :
1865 : integer (kind=int_kind), intent(in) :: &
1866 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1867 : icellT ! no. of cells where iceTmask = .true.
1868 :
1869 : real (kind=dbl_kind), intent(in) :: &
1870 : dtei ! 1/dte, where dte is subcycling timestep (1/s)
1871 :
1872 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1873 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1874 : indxTj ! compressed index in j-direction
1875 :
1876 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1877 : ! ice stress tensor (kg/s^2) in each corner of T cell
1878 : stressp_1, stressp_2, stressp_3, stressp_4, & ! sigma11+sigma22
1879 : stressm_1, stressm_2, stressm_3, stressm_4, & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1880 : stress12_1, stress12_2, stress12_3, stress12_4 ! sigma12
1881 :
1882 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1883 : ! structure tensor () in each corner of T cell
1884 : a11, a12, a11_1, a11_2, a11_3, a11_4, & ! components of
1885 : a12_1, a12_2, a12_3, a12_4 ! structure tensor ()
1886 :
1887 : ! local variables
1888 :
1889 : integer (kind=int_kind) :: &
1890 : i, j, ij
1891 :
1892 : real (kind=dbl_kind) :: &
1893 : mresult11, mresult12, & ! LCOV_EXCL_LINE
1894 0 : dteikth, p5kth
1895 :
1896 : real (kind=dbl_kind), parameter :: &
1897 : kth = p2*p001
1898 :
1899 : character(len=*), parameter :: subname = '(stepa)'
1900 :
1901 0 : dteikth = c1 / (dtei + kth)
1902 0 : p5kth = p5 * kth
1903 :
1904 0 : do ij = 1, icellT
1905 0 : i = indxTi(ij)
1906 0 : j = indxTj(ij)
1907 :
1908 : ! ne
1909 0 : call calc_ffrac(stressp_1(i,j), stressm_1(i,j), &
1910 : stress12_1(i,j), & ! LCOV_EXCL_LINE
1911 : a11_1(i,j), a12_1(i,j), & ! LCOV_EXCL_LINE
1912 0 : mresult11, mresult12)
1913 :
1914 0 : a11_1(i,j) = (a11_1(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
1915 0 : a12_1(i,j) = (a12_1(i,j)*dtei - mresult12) * dteikth ! implicit
1916 :
1917 :
1918 : ! nw
1919 0 : call calc_ffrac(stressp_2(i,j), stressm_2(i,j), &
1920 : stress12_2(i,j), & ! LCOV_EXCL_LINE
1921 : a11_2(i,j), a12_2(i,j), & ! LCOV_EXCL_LINE
1922 0 : mresult11, mresult12)
1923 :
1924 0 : a11_2(i,j) = (a11_2(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
1925 0 : a12_2(i,j) = (a12_2(i,j)*dtei - mresult12) * dteikth ! implicit
1926 :
1927 : ! sw
1928 0 : call calc_ffrac(stressp_3(i,j), stressm_3(i,j), &
1929 : stress12_3(i,j), & ! LCOV_EXCL_LINE
1930 : a11_3(i,j), a12_3(i,j), & ! LCOV_EXCL_LINE
1931 0 : mresult11, mresult12)
1932 :
1933 0 : a11_3(i,j) = (a11_3(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
1934 0 : a12_3(i,j) = (a12_3(i,j)*dtei - mresult12) * dteikth ! implicit
1935 :
1936 : ! se
1937 0 : call calc_ffrac(stressp_4(i,j), stressm_4(i,j), &
1938 : stress12_4(i,j), & ! LCOV_EXCL_LINE
1939 : a11_4(i,j), a12_4(i,j), & ! LCOV_EXCL_LINE
1940 0 : mresult11, mresult12)
1941 :
1942 0 : a11_4(i,j) = (a11_4(i,j)*dtei + p5kth - mresult11) * dteikth ! implicit
1943 0 : a12_4(i,j) = (a12_4(i,j)*dtei - mresult12) * dteikth ! implicit
1944 :
1945 : ! average structure tensor
1946 :
1947 0 : a11(i,j) = p25*(a11_1(i,j) + a11_2(i,j) + a11_3(i,j) + a11_4(i,j))
1948 0 : a12(i,j) = p25*(a12_1(i,j) + a12_2(i,j) + a12_3(i,j) + a12_4(i,j))
1949 :
1950 : enddo ! ij
1951 :
1952 0 : end subroutine stepa
1953 :
1954 : !=======================================================================
1955 : ! computes term in evolution equation for structure tensor which determines
1956 : ! the ice floe re-orientation due to fracture
1957 : ! Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2
1958 :
1959 :
1960 0 : subroutine calc_ffrac (stressp, stressm, &
1961 : stress12, & ! LCOV_EXCL_LINE
1962 : a1x, a2x, & ! LCOV_EXCL_LINE
1963 : mresult1, mresult2)
1964 :
1965 : real (kind=dbl_kind), intent(in) :: &
1966 : stressp, stressm, stress12, a1x, a2x
1967 :
1968 : real (kind=dbl_kind), intent(out) :: &
1969 : mresult1, mresult2
1970 :
1971 : ! local variables
1972 :
1973 : real (kind=dbl_kind) :: &
1974 : sigma11, sigma12, sigma22, & ! LCOV_EXCL_LINE
1975 : gamma, sigma_1, sigma_2, & ! LCOV_EXCL_LINE
1976 0 : Q11, Q12, Q11Q11, Q11Q12, Q12Q12
1977 :
1978 : real (kind=dbl_kind), parameter :: &
1979 : kfrac = p001, & ! LCOV_EXCL_LINE
1980 : threshold = c3*p1
1981 :
1982 : character(len=*), parameter :: subname = '(calc_ffrac)'
1983 :
1984 0 : sigma11 = p5*(stressp+stressm)
1985 0 : sigma12 = stress12
1986 0 : sigma22 = p5*(stressp-stressm)
1987 :
1988 : ! if ((sigma11-sigma22) == c0) then sigma11-sigma22 == 0 => stressn ==0
1989 0 : if (stressm == c0) then
1990 0 : gamma = p5*(pih)
1991 : else
1992 0 : gamma = p5*atan2((c2*sigma12),(sigma11-sigma22))
1993 : endif
1994 :
1995 : ! rotate tensor to get into sigma principal axis
1996 :
1997 0 : Q11 = cos(gamma)
1998 0 : Q12 = sin(gamma)
1999 :
2000 0 : Q11Q11 = Q11*Q11
2001 0 : Q11Q12 = Q11*Q12
2002 0 : Q12Q12 = Q12*Q12
2003 :
2004 : sigma_1 = Q11Q11*sigma11 + c2*Q11Q12*sigma12 &
2005 0 : + Q12Q12*sigma22 ! S(1,1)
2006 : sigma_2 = Q12Q12*sigma11 - c2*Q11Q12*sigma12 &
2007 0 : + Q11Q11*sigma22 ! S(2,2)
2008 :
2009 : ! Pure divergence
2010 0 : if ((sigma_1 >= c0).and.(sigma_2 >= c0)) then
2011 0 : mresult1 = c0
2012 0 : mresult2 = c0
2013 :
2014 : ! Unconfined compression: cracking of blocks not along the axial splitting direction
2015 : ! which leads to the loss of their shape, so we again model it through diffusion
2016 0 : elseif ((sigma_1 >= c0).and.(sigma_2 < c0)) then
2017 0 : mresult1 = kfrac * (a1x - Q12Q12)
2018 0 : mresult2 = kfrac * (a2x + Q11Q12)
2019 :
2020 : ! Shear faulting
2021 0 : elseif (sigma_2 == c0) then
2022 0 : mresult1 = c0
2023 0 : mresult2 = c0
2024 0 : elseif ((sigma_1 <= c0).and.(sigma_1/sigma_2 <= threshold)) then
2025 0 : mresult1 = kfrac * (a1x - Q12Q12)
2026 0 : mresult2 = kfrac * (a2x + Q11Q12)
2027 :
2028 : ! Horizontal spalling
2029 : else
2030 0 : mresult1 = c0
2031 0 : mresult2 = c0
2032 : endif
2033 :
2034 0 : end subroutine calc_ffrac
2035 :
2036 : !=======================================================================
2037 : !---! these subroutines write/read Fortran unformatted data files ..
2038 : !=======================================================================
2039 : ! Dumps all values needed for a restart
2040 :
2041 0 : subroutine write_restart_eap ()
2042 :
2043 : use ice_restart, only: write_restart_field
2044 :
2045 : ! local variables
2046 :
2047 : logical (kind=log_kind) :: diag
2048 :
2049 : character(len=*), parameter :: subname = '(write_restart_eap)'
2050 :
2051 0 : diag = .true.
2052 :
2053 : !-----------------------------------------------------------------
2054 : ! structure tensor
2055 : !-----------------------------------------------------------------
2056 :
2057 0 : call write_restart_field(nu_dump_eap,0,a11_1,'ruf8','a11_1',1,diag)
2058 0 : call write_restart_field(nu_dump_eap,0,a11_3,'ruf8','a11_3',1,diag)
2059 0 : call write_restart_field(nu_dump_eap,0,a11_2,'ruf8','a11_2',1,diag)
2060 0 : call write_restart_field(nu_dump_eap,0,a11_4,'ruf8','a11_4',1,diag)
2061 :
2062 0 : call write_restart_field(nu_dump_eap,0,a12_1,'ruf8','a12_1',1,diag)
2063 0 : call write_restart_field(nu_dump_eap,0,a12_3,'ruf8','a12_3',1,diag)
2064 0 : call write_restart_field(nu_dump_eap,0,a12_2,'ruf8','a12_2',1,diag)
2065 0 : call write_restart_field(nu_dump_eap,0,a12_4,'ruf8','a12_4',1,diag)
2066 :
2067 0 : end subroutine write_restart_eap
2068 :
2069 : !=======================================================================
2070 : ! Reads all values needed for elastic anisotropic plastic dynamics restart
2071 :
2072 0 : subroutine read_restart_eap()
2073 :
2074 : use ice_blocks, only: nghost
2075 : use ice_boundary, only: ice_HaloUpdate_stress
2076 : use ice_constants, only: &
2077 : field_loc_center, field_type_scalar
2078 : use ice_domain, only: nblocks, halo_info
2079 : use ice_grid, only: grid_type
2080 : use ice_restart, only: read_restart_field
2081 :
2082 : ! local variables
2083 :
2084 : integer (kind=int_kind) :: &
2085 : i, j, iblk
2086 :
2087 : logical (kind=log_kind) :: &
2088 : diag
2089 :
2090 : character(len=*), parameter :: subname = '(read_restart_eap)'
2091 :
2092 0 : diag = .true.
2093 :
2094 : !-----------------------------------------------------------------
2095 : ! Structure tensor must be read and scattered in pairs in order
2096 : ! to properly match corner values across a tripole grid cut.
2097 : !-----------------------------------------------------------------
2098 0 : if (my_task == master_task) write(nu_diag,*) &
2099 0 : 'structure tensor restart data'
2100 :
2101 : call read_restart_field(nu_restart_eap,0,a11_1,'ruf8', &
2102 0 : 'a11_1',1,diag,field_loc_center,field_type_scalar) ! a11_1
2103 : call read_restart_field(nu_restart_eap,0,a11_3,'ruf8', &
2104 0 : 'a11_3',1,diag,field_loc_center,field_type_scalar) ! a11_3
2105 : call read_restart_field(nu_restart_eap,0,a11_2,'ruf8', &
2106 0 : 'a11_2',1,diag,field_loc_center,field_type_scalar) ! a11_2
2107 : call read_restart_field(nu_restart_eap,0,a11_4,'ruf8', &
2108 0 : 'a11_4',1,diag,field_loc_center,field_type_scalar) ! a11_4
2109 :
2110 : call read_restart_field(nu_restart_eap,0,a12_1,'ruf8', &
2111 0 : 'a12_1',1,diag,field_loc_center,field_type_scalar) ! a12_1
2112 : call read_restart_field(nu_restart_eap,0,a12_3,'ruf8', &
2113 0 : 'a12_3',1,diag,field_loc_center,field_type_scalar) ! a12_3
2114 : call read_restart_field(nu_restart_eap,0,a12_2,'ruf8', &
2115 0 : 'a12_2',1,diag,field_loc_center,field_type_scalar) ! a12_2
2116 : call read_restart_field(nu_restart_eap,0,a12_4,'ruf8', &
2117 0 : 'a12_4',1,diag,field_loc_center,field_type_scalar) ! a12_4
2118 :
2119 0 : if (trim(grid_type) == 'tripole') then
2120 :
2121 : call ice_HaloUpdate_stress(a11_1, a11_3, halo_info, &
2122 0 : field_loc_center, field_type_scalar)
2123 : call ice_HaloUpdate_stress(a11_3, a11_1, halo_info, &
2124 0 : field_loc_center, field_type_scalar)
2125 : call ice_HaloUpdate_stress(a11_2, a11_4, halo_info, &
2126 0 : field_loc_center, field_type_scalar)
2127 : call ice_HaloUpdate_stress(a11_4, a11_2, halo_info, &
2128 0 : field_loc_center, field_type_scalar)
2129 : call ice_HaloUpdate_stress(a12_1, a12_3, halo_info, &
2130 0 : field_loc_center, field_type_scalar)
2131 : call ice_HaloUpdate_stress(a12_3, a12_1, halo_info, &
2132 0 : field_loc_center, field_type_scalar)
2133 : call ice_HaloUpdate_stress(a12_2, a12_4, halo_info, &
2134 0 : field_loc_center, field_type_scalar)
2135 : call ice_HaloUpdate_stress(a12_4, a12_2, halo_info, &
2136 0 : field_loc_center, field_type_scalar)
2137 :
2138 : endif
2139 :
2140 : !-----------------------------------------------------------------
2141 : ! Ensure unused values in west and south ghost cells are 0
2142 : !-----------------------------------------------------------------
2143 :
2144 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
2145 0 : do iblk = 1, nblocks
2146 0 : do j = 1, nghost
2147 0 : do i = 1, nx_block
2148 0 : a11_1 (i,j,iblk) = c0
2149 0 : a11_2 (i,j,iblk) = c0
2150 0 : a11_3 (i,j,iblk) = c0
2151 0 : a11_4 (i,j,iblk) = c0
2152 0 : a12_1 (i,j,iblk) = c0
2153 0 : a12_2 (i,j,iblk) = c0
2154 0 : a12_3 (i,j,iblk) = c0
2155 0 : a12_4 (i,j,iblk) = c0
2156 : enddo
2157 : enddo
2158 0 : do j = 1, ny_block
2159 0 : do i = 1, nghost
2160 0 : a11_1 (i,j,iblk) = c0
2161 0 : a11_2 (i,j,iblk) = c0
2162 0 : a11_3 (i,j,iblk) = c0
2163 0 : a11_4 (i,j,iblk) = c0
2164 0 : a12_1 (i,j,iblk) = c0
2165 0 : a12_2 (i,j,iblk) = c0
2166 0 : a12_3 (i,j,iblk) = c0
2167 0 : a12_4 (i,j,iblk) = c0
2168 : enddo
2169 : enddo
2170 : enddo
2171 : !$OMP END PARALLEL DO
2172 :
2173 0 : end subroutine read_restart_eap
2174 :
2175 : !=======================================================================
2176 :
2177 : end module ice_dyn_eap
2178 :
2179 : !=======================================================================
|
