Line data Source code
1 : !=======================================================================
2 : !
3 : ! Elastic-viscous-plastic sea ice dynamics model
4 : ! Computes ice velocity and deformation
5 : !
6 : ! See:
7 : !
8 : ! Hunke, E. C., and J. K. Dukowicz (1997). An elastic-viscous-plastic model
9 : ! for sea ice dynamics. J. Phys. Oceanogr., 27, 1849-1867.
10 : !
11 : ! Hunke, E. C. (2001). Viscous-Plastic Sea Ice Dynamics with the EVP Model:
12 : ! Linearization Issues. J. Comput. Phys., 170, 18-38.
13 : !
14 : ! Hunke, E. C., and J. K. Dukowicz (2002). The Elastic-Viscous-Plastic
15 : ! Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates
16 : ! on a Sphere - Incorporation of Metric Terms. Mon. Weather Rev.,
17 : ! 130, 1848-1865.
18 : !
19 : ! Hunke, E. C., and J. K. Dukowicz (2003). The sea ice momentum
20 : ! equation in the free drift regime. Los Alamos Tech. Rep. LA-UR-03-2219.
21 : !
22 : ! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys.
23 : ! Oceanogr., 9, 817-846.
24 : !
25 : ! Bouillon, S., T. Fichefet, V. Legat and G. Madec (2013). The
26 : ! elastic-viscous-plastic method revisited. Ocean Model., 71, 2-12.
27 : !
28 : ! author: Elizabeth C. Hunke, LANL
29 : !
30 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb (LANL)
31 : ! 2004: Block structure added by William Lipscomb
32 : ! 2005: Removed boundary calls for stress arrays (WHL)
33 : ! 2006: Streamlined for efficiency by Elizabeth Hunke
34 : ! Converted to free source form (F90)
35 :
36 : module ice_dyn_evp
37 :
38 : use ice_kinds_mod
39 : use ice_communicate, only: my_task, master_task
40 : use ice_constants, only: field_loc_center, field_loc_NEcorner, &
41 : field_loc_Nface, field_loc_Eface, & ! LCOV_EXCL_LINE
42 : field_type_scalar, field_type_vector
43 : use ice_constants, only: c0, p027, p055, p111, p166, &
44 : p222, p25, p333, p5, c1
45 : use ice_dyn_shared, only: stepu, stepuv_CD, stepu_C, stepv_C, &
46 : dyn_prep1, dyn_prep2, dyn_finish, & ! LCOV_EXCL_LINE
47 : ndte, yield_curve, ecci, denom1, arlx1i, fcor_blk, fcorE_blk, fcorN_blk, & ! LCOV_EXCL_LINE
48 : uvel_init, vvel_init, uvelE_init, vvelE_init, uvelN_init, vvelN_init, & ! LCOV_EXCL_LINE
49 : seabed_stress_factor_LKD, seabed_stress_factor_prob, seabed_stress_method, & ! LCOV_EXCL_LINE
50 : seabed_stress, Ktens, revp
51 : use ice_fileunits, only: nu_diag
52 : use ice_exit, only: abort_ice
53 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
54 : use icepack_intfc, only: icepack_ice_strength, icepack_query_parameters
55 :
56 : implicit none
57 : private
58 : ! all c or cd
59 : real (kind=dbl_kind), allocatable :: &
60 : uocnN (:,:,:) , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
61 : vocnN (:,:,:) , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
62 : ss_tltxN (:,:,:) , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
63 : ss_tltyN (:,:,:) , & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
64 : cdn_ocnN (:,:,:) , & ! ocn drag coefficient ! LCOV_EXCL_LINE
65 : waterxN (:,:,:) , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
66 : wateryN (:,:,:) , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
67 : forcexN (:,:,:) , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
68 : forceyN (:,:,:) , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
69 : aiN (:,:,:) , & ! ice fraction on N-grid ! LCOV_EXCL_LINE
70 : nmass (:,:,:) , & ! total mass of ice and snow (N grid) ! LCOV_EXCL_LINE
71 : nmassdti (:,:,:) ! mass of N-cell/dte (kg/m^2 s)
72 : ! all c or d
73 : real (kind=dbl_kind), allocatable :: &
74 : uocnE (:,:,:) , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
75 : vocnE (:,:,:) , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
76 : ss_tltxE (:,:,:) , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
77 : ss_tltyE (:,:,:) , & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
78 : cdn_ocnE (:,:,:) , & ! ocn drag coefficient ! LCOV_EXCL_LINE
79 : waterxE (:,:,:) , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
80 : wateryE (:,:,:) , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
81 : forcexE (:,:,:) , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
82 : forceyE (:,:,:) , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
83 : aiE (:,:,:) , & ! ice fraction on E-grid ! LCOV_EXCL_LINE
84 : emass (:,:,:) , & ! total mass of ice and snow (E grid) ! LCOV_EXCL_LINE
85 : emassdti (:,:,:) ! mass of E-cell/dte (kg/m^2 s)
86 :
87 : real (kind=dbl_kind), allocatable :: &
88 : strengthU(:,:,:) , & ! strength averaged to U points ! LCOV_EXCL_LINE
89 : divergU (:,:,:) , & ! div array on U points, differentiate from divu ! LCOV_EXCL_LINE
90 : tensionU (:,:,:) , & ! tension array on U points ! LCOV_EXCL_LINE
91 : shearU (:,:,:) , & ! shear array on U points ! LCOV_EXCL_LINE
92 : deltaU (:,:,:) , & ! delta array on U points ! LCOV_EXCL_LINE
93 : zetax2T (:,:,:) , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
94 : zetax2U (:,:,:) , & ! zetax2T averaged to U points ! LCOV_EXCL_LINE
95 : etax2T (:,:,:) , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
96 : etax2U (:,:,:) ! etax2T averaged to U points
97 :
98 : real (kind=dbl_kind), allocatable :: &
99 : uocnU (:,:,:) , & ! i ocean current (m/s) ! LCOV_EXCL_LINE
100 : vocnU (:,:,:) , & ! j ocean current (m/s) ! LCOV_EXCL_LINE
101 : ss_tltxU (:,:,:) , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
102 : ss_tltyU (:,:,:) , & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
103 : cdn_ocnU (:,:,:) , & ! ocn drag coefficient ! LCOV_EXCL_LINE
104 : tmass (:,:,:) , & ! total mass of ice and snow (kg/m^2) ! LCOV_EXCL_LINE
105 : waterxU (:,:,:) , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
106 : wateryU (:,:,:) , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
107 : forcexU (:,:,:) , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
108 : forceyU (:,:,:) , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
109 : umass (:,:,:) , & ! total mass of ice and snow (u grid) ! LCOV_EXCL_LINE
110 : umassdti (:,:,:) ! mass of U-cell/dte (kg/m^2 s)
111 :
112 : public :: evp, init_evp
113 :
114 : !=======================================================================
115 :
116 : contains
117 :
118 : !=======================================================================
119 : ! Elastic-viscous-plastic dynamics driver
120 : !
121 37 : subroutine init_evp
122 : use ice_blocks, only: nx_block, ny_block
123 : use ice_domain_size, only: max_blocks
124 : use ice_grid, only: grid_ice
125 : use ice_calendar, only: dt_dyn
126 : use ice_dyn_shared, only: init_dyn_shared
127 :
128 : !allocate c and cd grid var. Follow structucre of eap
129 : integer (int_kind) :: ierr
130 :
131 : character(len=*), parameter :: subname = '(alloc_dyn_evp)'
132 :
133 37 : call init_dyn_shared(dt_dyn)
134 :
135 : allocate( uocnU (nx_block,ny_block,max_blocks), & ! i ocean current (m/s)
136 : vocnU (nx_block,ny_block,max_blocks), & ! j ocean current (m/s) ! LCOV_EXCL_LINE
137 : ss_tltxU (nx_block,ny_block,max_blocks), & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
138 : ss_tltyU (nx_block,ny_block,max_blocks), & ! sea surface slope, y-direction (m/m) ! LCOV_EXCL_LINE
139 : cdn_ocnU (nx_block,ny_block,max_blocks), & ! ocn drag coefficient ! LCOV_EXCL_LINE
140 : tmass (nx_block,ny_block,max_blocks), & ! total mass of ice and snow (kg/m^2) ! LCOV_EXCL_LINE
141 : waterxU (nx_block,ny_block,max_blocks), & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
142 : wateryU (nx_block,ny_block,max_blocks), & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
143 : forcexU (nx_block,ny_block,max_blocks), & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
144 : forceyU (nx_block,ny_block,max_blocks), & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
145 : umass (nx_block,ny_block,max_blocks), & ! total mass of ice and snow (u grid) ! LCOV_EXCL_LINE
146 : umassdti (nx_block,ny_block,max_blocks), & ! mass of U-cell/dte (kg/m^2 s) ! LCOV_EXCL_LINE
147 37 : stat=ierr)
148 37 : if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory B-Grid evp')
149 :
150 :
151 37 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
152 :
153 : allocate( strengthU(nx_block,ny_block,max_blocks), &
154 : divergU (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
155 : tensionU (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
156 : shearU (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
157 : deltaU (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
158 : zetax2T (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
159 : zetax2U (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
160 : etax2T (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
161 : etax2U (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
162 0 : stat=ierr)
163 0 : if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory U evp')
164 :
165 : allocate( uocnN (nx_block,ny_block,max_blocks), &
166 : vocnN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
167 : ss_tltxN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
168 : ss_tltyN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
169 : cdn_ocnN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
170 : waterxN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
171 : wateryN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
172 : forcexN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
173 : forceyN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
174 : aiN (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
175 : nmass (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
176 : nmassdti (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
177 0 : stat=ierr)
178 0 : if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory N evp')
179 :
180 : allocate( uocnE (nx_block,ny_block,max_blocks), &
181 : vocnE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
182 : ss_tltxE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
183 : ss_tltyE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
184 : cdn_ocnE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
185 : waterxE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
186 : wateryE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
187 : forcexE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
188 : forceyE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
189 : aiE (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
190 : emass (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
191 : emassdti (nx_block,ny_block,max_blocks), & ! LCOV_EXCL_LINE
192 0 : stat=ierr)
193 0 : if (ierr/=0) call abort_ice(subname//' ERROR: Out of memory E evp')
194 :
195 : endif
196 :
197 37 : end subroutine init_evp
198 :
199 : #ifdef CICE_IN_NEMO
200 : ! Wind stress is set during this routine from the values supplied
201 : ! via NEMO (unless calc_strair is true). These values are supplied
202 : ! rotated on u grid and multiplied by aice. strairxT = 0 in this
203 : ! case so operations in dyn_prep1 are pointless but carried out to
204 : ! minimise code changes.
205 : #endif
206 : !
207 : ! author: Elizabeth C. Hunke, LANL
208 :
209 5784 : subroutine evp (dt)
210 :
211 : use ice_arrays_column, only: Cdn_ocn
212 : use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
213 : ice_HaloDestroy, ice_HaloUpdate_stress
214 : use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
215 : use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn
216 : use ice_domain_size, only: max_blocks, ncat, nx_global, ny_global
217 : use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
218 : strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, & ! LCOV_EXCL_LINE
219 : strtltxU, strtltyU, strocnxU, strocnyU, strintxU, strintyU, taubxU, taubyU, & ! LCOV_EXCL_LINE
220 : strax, stray, & ! LCOV_EXCL_LINE
221 : TbU, hwater, & ! LCOV_EXCL_LINE
222 : strairxN, strairyN, fmN, & ! LCOV_EXCL_LINE
223 : strtltxN, strtltyN, strocnxN, strocnyN, strintxN, strintyN, taubxN, taubyN, & ! LCOV_EXCL_LINE
224 : TbN, & ! LCOV_EXCL_LINE
225 : strairxE, strairyE, fmE, & ! LCOV_EXCL_LINE
226 : strtltxE, strtltyE, strocnxE, strocnyE, strintxE, strintyE, taubxE, taubyE, & ! LCOV_EXCL_LINE
227 : TbE, & ! LCOV_EXCL_LINE
228 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
229 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
230 : stress12_1, stress12_2, stress12_3, stress12_4, & ! LCOV_EXCL_LINE
231 : stresspT, stressmT, stress12T, & ! LCOV_EXCL_LINE
232 : stresspU, stressmU, stress12U
233 : use ice_grid, only: tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, &
234 : dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, & ! LCOV_EXCL_LINE
235 : ratiodxN, ratiodxNr, ratiodyE, ratiodyEr, & ! LCOV_EXCL_LINE
236 : dxhy, dyhx, cxp, cyp, cxm, cym, & ! LCOV_EXCL_LINE
237 : tarear, uarear, earear, narear, grid_average_X2Y, uarea, & ! LCOV_EXCL_LINE
238 : grid_type, grid_ice, & ! LCOV_EXCL_LINE
239 : grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
240 : use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, uvelN, vvelN, &
241 : uvelE, vvelE, divu, shear, & ! LCOV_EXCL_LINE
242 : aice_init, aice0, aicen, vicen, strength
243 : use ice_timers, only: timer_dynamics, timer_bound, &
244 : ice_timer_start, ice_timer_stop, timer_evp_1d, timer_evp_2d
245 : use ice_dyn_evp_1d, only: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_kernel, &
246 : ice_dyn_evp_1d_copyout
247 : use ice_dyn_shared, only: evp_algorithm, stack_fields, unstack_fields, &
248 : DminTarea, visc_method, deformations, deformationsC_T, deformationsCD_T, & ! LCOV_EXCL_LINE
249 : strain_rates_U, & ! LCOV_EXCL_LINE
250 : iceTmask, iceUmask, iceEmask, iceNmask, & ! LCOV_EXCL_LINE
251 : dyn_haloUpdate, fld2, fld3, fld4
252 :
253 : real (kind=dbl_kind), intent(in) :: &
254 : dt ! time step
255 :
256 : ! local variables
257 :
258 : integer (kind=int_kind) :: &
259 : ksub , & ! subcycle step ! LCOV_EXCL_LINE
260 : iblk , & ! block index ! LCOV_EXCL_LINE
261 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
262 : i, j, ij ! local indices
263 :
264 : integer (kind=int_kind), dimension(max_blocks) :: &
265 : icellT , & ! no. of cells where iceTmask = .true. ! LCOV_EXCL_LINE
266 : icellN , & ! no. of cells where iceNmask = .true. ! LCOV_EXCL_LINE
267 : icellE , & ! no. of cells where iceEmask = .true. ! LCOV_EXCL_LINE
268 11568 : icellU ! no. of cells where iceUmask = .true.
269 :
270 : integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks) :: &
271 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
272 : indxTj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
273 : indxEi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
274 : indxEj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
275 : indxNi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
276 : indxNj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
277 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
278 11568 : indxUj ! compressed index in j-direction
279 :
280 : real (kind=dbl_kind), dimension(nx_block,ny_block,8):: &
281 10021008 : strtmp ! stress combinations for momentum equation
282 :
283 : logical (kind=log_kind) :: &
284 : calc_strair ! calculate air/ice stress
285 :
286 : integer (kind=int_kind), dimension (nx_block,ny_block,max_blocks) :: &
287 10128 : halomask ! generic halo mask
288 :
289 : type (ice_halo) :: &
290 : halo_info_mask ! ghost cell update info for masked halo
291 :
292 : type (block) :: &
293 : this_block ! block information for current block
294 :
295 : character(len=*), parameter :: subname = '(evp)'
296 :
297 5784 : call ice_timer_start(timer_dynamics) ! dynamics
298 :
299 : !-----------------------------------------------------------------
300 : ! Initialize
301 : !-----------------------------------------------------------------
302 :
303 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
304 :
305 0 : strengthU(:,:,:) = c0
306 0 : divergU (:,:,:) = c0
307 0 : tensionU (:,:,:) = c0
308 0 : shearU (:,:,:) = c0
309 0 : deltaU (:,:,:) = c0
310 0 : zetax2T (:,:,:) = c0
311 0 : zetax2U (:,:,:) = c0
312 0 : etax2T (:,:,:) = c0
313 0 : etax2U (:,:,:) = c0
314 :
315 : endif
316 :
317 : ! This call is needed only if dt changes during runtime.
318 : ! call set_evp_parameters (dt)
319 :
320 : !-----------------------------------------------------------------
321 : ! boundary updates
322 : ! commented out because the ghost cells are freshly
323 : ! updated after cleanup_itd
324 : !-----------------------------------------------------------------
325 :
326 : ! call ice_timer_start(timer_bound)
327 : ! call ice_HaloUpdate (aice, halo_info, &
328 : ! field_loc_center, field_type_scalar)
329 : ! call ice_HaloUpdate (vice, halo_info, &
330 : ! field_loc_center, field_type_scalar)
331 : ! call ice_HaloUpdate (vsno, halo_info, &
332 : ! field_loc_center, field_type_scalar)
333 : ! call ice_timer_stop(timer_bound)
334 :
335 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) SCHEDULE(runtime)
336 8688 : do iblk = 1, nblocks
337 :
338 204456 : do j = 1, ny_block
339 7151880 : do i = 1, nx_block
340 6947424 : rdg_conv (i,j,iblk) = c0
341 6947424 : rdg_shear(i,j,iblk) = c0
342 6947424 : divu (i,j,iblk) = c0
343 7146096 : shear(i,j,iblk) = c0
344 : enddo
345 : enddo
346 :
347 : !-----------------------------------------------------------------
348 : ! preparation for dynamics
349 : !-----------------------------------------------------------------
350 :
351 5784 : this_block = get_block(blocks_ice(iblk),iblk)
352 5784 : ilo = this_block%ilo
353 5784 : ihi = this_block%ihi
354 5784 : jlo = this_block%jlo
355 5784 : jhi = this_block%jhi
356 :
357 : call dyn_prep1 (nx_block, ny_block, &
358 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
359 : aice (:,:,iblk), vice (:,:,iblk), & ! LCOV_EXCL_LINE
360 : vsno (:,:,iblk), tmask (:,:,iblk), & ! LCOV_EXCL_LINE
361 11568 : tmass (:,:,iblk), iceTmask(:,:,iblk))
362 :
363 : enddo ! iblk
364 : !$OMP END PARALLEL DO
365 :
366 5784 : call ice_timer_start(timer_bound)
367 : call ice_HaloUpdate (iceTmask, halo_info, &
368 5784 : field_loc_center, field_type_scalar)
369 5784 : call ice_timer_stop(timer_bound)
370 :
371 : !-----------------------------------------------------------------
372 : ! convert fields from T to U grid
373 : !-----------------------------------------------------------------
374 :
375 5784 : call stack_fields(tmass, aice_init, cdn_ocn, fld3)
376 : call ice_HaloUpdate (fld3, halo_info, &
377 5784 : field_loc_center, field_type_scalar)
378 5784 : call stack_fields(uocn, vocn, ss_tltx, ss_tlty, fld4)
379 : call ice_HaloUpdate (fld4, halo_info, &
380 5784 : field_loc_center, field_type_vector)
381 5784 : call unstack_fields(fld3, tmass, aice_init, cdn_ocn)
382 5784 : call unstack_fields(fld4, uocn, vocn, ss_tltx, ss_tlty)
383 :
384 5784 : call grid_average_X2Y('S', tmass , 'T' , umass , 'U')
385 5784 : call grid_average_X2Y('S', aice_init, 'T' , aiU , 'U')
386 5784 : call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnU, 'U')
387 5784 : call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnU , 'U')
388 5784 : call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnU , 'U')
389 5784 : call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxU, 'U')
390 5784 : call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyU, 'U')
391 :
392 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
393 0 : call grid_average_X2Y('S', tmass , 'T' , emass , 'E')
394 0 : call grid_average_X2Y('S', aice_init, 'T' , aie , 'E')
395 0 : call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnE, 'E')
396 0 : call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnE , 'E')
397 0 : call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnE , 'E')
398 0 : call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxE, 'E')
399 0 : call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyE, 'E')
400 0 : call grid_average_X2Y('S', tmass , 'T' , nmass , 'N')
401 0 : call grid_average_X2Y('S', aice_init, 'T' , ain , 'N')
402 0 : call grid_average_X2Y('S', cdn_ocn , 'T' , cdn_ocnN, 'N')
403 0 : call grid_average_X2Y('S', uocn , grid_ocn_dynu, uocnN , 'N')
404 0 : call grid_average_X2Y('S', vocn , grid_ocn_dynv, vocnN , 'N')
405 0 : call grid_average_X2Y('S', ss_tltx , grid_ocn_dynu, ss_tltxN, 'N')
406 0 : call grid_average_X2Y('S', ss_tlty , grid_ocn_dynv, ss_tltyN, 'N')
407 : endif
408 : !----------------------------------------------------------------
409 : ! Set wind stress to values supplied via NEMO or other forcing
410 : ! Map T to U, N, E as needed
411 : ! This wind stress is rotated on u grid and multiplied by aice in coupler
412 : !----------------------------------------------------------------
413 5784 : call icepack_query_parameters(calc_strair_out=calc_strair)
414 5784 : call icepack_warnings_flush(nu_diag)
415 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
416 0 : file=__FILE__, line=__LINE__)
417 :
418 5784 : if (.not. calc_strair) then
419 0 : call grid_average_X2Y('F', strax, grid_atm_dynu, strairxU, 'U')
420 0 : call grid_average_X2Y('F', stray, grid_atm_dynv, strairyU, 'U')
421 : else
422 : call ice_HaloUpdate (strairxT, halo_info, &
423 5784 : field_loc_center, field_type_vector)
424 : call ice_HaloUpdate (strairyT, halo_info, &
425 5784 : field_loc_center, field_type_vector)
426 5784 : call grid_average_X2Y('F', strairxT, 'T', strairxU, 'U')
427 5784 : call grid_average_X2Y('F', strairyT, 'T', strairyU, 'U')
428 : endif
429 :
430 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
431 0 : if (.not. calc_strair) then
432 0 : call grid_average_X2Y('F', strax , grid_atm_dynu, strairxN, 'N')
433 0 : call grid_average_X2Y('F', stray , grid_atm_dynv, strairyN, 'N')
434 0 : call grid_average_X2Y('F', strax , grid_atm_dynu, strairxE, 'E')
435 0 : call grid_average_X2Y('F', stray , grid_atm_dynv, strairyE, 'E')
436 : else
437 0 : call grid_average_X2Y('F', strairxT, 'T' , strairxN, 'N')
438 0 : call grid_average_X2Y('F', strairyT, 'T' , strairyN, 'N')
439 0 : call grid_average_X2Y('F', strairxT, 'T' , strairxE, 'E')
440 0 : call grid_average_X2Y('F', strairyT, 'T' , strairyE, 'E')
441 : endif
442 : endif
443 :
444 5784 : if (trim(grid_ice) == 'B') then
445 2880 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,ij,i,j) SCHEDULE(runtime)
446 8688 : do iblk = 1, nblocks
447 :
448 : !-----------------------------------------------------------------
449 : ! more preparation for dynamics
450 : !-----------------------------------------------------------------
451 :
452 5784 : this_block = get_block(blocks_ice(iblk),iblk)
453 5784 : ilo = this_block%ilo
454 5784 : ihi = this_block%ihi
455 5784 : jlo = this_block%jlo
456 5784 : jhi = this_block%jhi
457 :
458 : call dyn_prep2 (nx_block, ny_block, &
459 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
460 : icellT (iblk), icellU (iblk), & ! LCOV_EXCL_LINE
461 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
462 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
463 : aiU (:,:,iblk), umass (:,:,iblk), & ! LCOV_EXCL_LINE
464 : umassdti (:,:,iblk), fcor_blk (:,:,iblk), & ! LCOV_EXCL_LINE
465 : umask (:,:,iblk), & ! LCOV_EXCL_LINE
466 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
467 : strairxU (:,:,iblk), strairyU (:,:,iblk), & ! LCOV_EXCL_LINE
468 : ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & ! LCOV_EXCL_LINE
469 : iceTmask (:,:,iblk), iceUmask (:,:,iblk), & ! LCOV_EXCL_LINE
470 : fmU (:,:,iblk), dt , & ! LCOV_EXCL_LINE
471 : strtltxU (:,:,iblk), strtltyU (:,:,iblk), & ! LCOV_EXCL_LINE
472 : strocnxU (:,:,iblk), strocnyU (:,:,iblk), & ! LCOV_EXCL_LINE
473 : strintxU (:,:,iblk), strintyU (:,:,iblk), & ! LCOV_EXCL_LINE
474 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
475 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
476 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
477 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
478 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
479 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
480 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
481 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
482 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
483 : uvel_init (:,:,iblk), vvel_init (:,:,iblk), & ! LCOV_EXCL_LINE
484 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
485 5784 : TbU (:,:,iblk))
486 :
487 : !-----------------------------------------------------------------
488 : ! ice strength
489 : !-----------------------------------------------------------------
490 :
491 7151880 : strength(:,:,iblk) = c0 ! initialize
492 6077205 : do ij = 1, icellT(iblk)
493 6065637 : i = indxTi(ij, iblk)
494 6065637 : j = indxTj(ij, iblk)
495 : call icepack_ice_strength(ncat = ncat, &
496 : aice = aice (i,j, iblk), & ! LCOV_EXCL_LINE
497 : vice = vice (i,j, iblk), & ! LCOV_EXCL_LINE
498 : aice0 = aice0 (i,j, iblk), & ! LCOV_EXCL_LINE
499 : aicen = aicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
500 : vicen = vicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
501 6071421 : strength = strength(i,j, iblk))
502 : enddo ! ij
503 :
504 : enddo ! iblk
505 : !$OMP END PARALLEL DO
506 0 : elseif (trim(grid_ice) == 'CD' .or. grid_ice == 'C') then
507 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,ij,i,j) SCHEDULE(runtime)
508 0 : do iblk = 1, nblocks
509 :
510 0 : this_block = get_block(blocks_ice(iblk),iblk)
511 0 : ilo = this_block%ilo
512 0 : ihi = this_block%ihi
513 0 : jlo = this_block%jlo
514 0 : jhi = this_block%jhi
515 :
516 : call dyn_prep2 (nx_block, ny_block, &
517 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
518 : icellT (iblk), icellU (iblk), & ! LCOV_EXCL_LINE
519 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
520 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
521 : aiU (:,:,iblk), umass (:,:,iblk), & ! LCOV_EXCL_LINE
522 : umassdti (:,:,iblk), fcor_blk (:,:,iblk), & ! LCOV_EXCL_LINE
523 : umaskCD (:,:,iblk), & ! LCOV_EXCL_LINE
524 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
525 : strairxU (:,:,iblk), strairyU (:,:,iblk), & ! LCOV_EXCL_LINE
526 : ss_tltxU (:,:,iblk), ss_tltyU (:,:,iblk), & ! LCOV_EXCL_LINE
527 : iceTmask (:,:,iblk), iceUmask (:,:,iblk), & ! LCOV_EXCL_LINE
528 : fmU (:,:,iblk), dt, & ! LCOV_EXCL_LINE
529 : strtltxU (:,:,iblk), strtltyU (:,:,iblk), & ! LCOV_EXCL_LINE
530 : strocnxU (:,:,iblk), strocnyU (:,:,iblk), & ! LCOV_EXCL_LINE
531 : strintxU (:,:,iblk), strintyU (:,:,iblk), & ! LCOV_EXCL_LINE
532 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
533 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
534 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
535 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
536 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
537 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
538 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
539 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
540 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
541 : uvel_init (:,:,iblk), vvel_init (:,:,iblk), & ! LCOV_EXCL_LINE
542 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
543 0 : TbU (:,:,iblk))
544 :
545 : !-----------------------------------------------------------------
546 : ! ice strength
547 : !-----------------------------------------------------------------
548 :
549 0 : strength(:,:,iblk) = c0 ! initialize
550 0 : do ij = 1, icellT(iblk)
551 0 : i = indxTi(ij, iblk)
552 0 : j = indxTj(ij, iblk)
553 : call icepack_ice_strength(ncat = ncat, &
554 : aice = aice (i,j, iblk), & ! LCOV_EXCL_LINE
555 : vice = vice (i,j, iblk), & ! LCOV_EXCL_LINE
556 : aice0 = aice0 (i,j, iblk), & ! LCOV_EXCL_LINE
557 : aicen = aicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
558 : vicen = vicen (i,j,:,iblk), & ! LCOV_EXCL_LINE
559 0 : strength = strength(i,j, iblk) )
560 : enddo ! ij
561 :
562 :
563 : !-----------------------------------------------------------------
564 : ! more preparation for dynamics on N grid
565 : !-----------------------------------------------------------------
566 :
567 : call dyn_prep2 (nx_block, ny_block, &
568 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
569 : icellT (iblk), icellN (iblk), & ! LCOV_EXCL_LINE
570 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
571 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
572 : aiN (:,:,iblk), nmass (:,:,iblk), & ! LCOV_EXCL_LINE
573 : nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & ! LCOV_EXCL_LINE
574 : nmask (:,:,iblk), & ! LCOV_EXCL_LINE
575 : uocnN (:,:,iblk), vocnN (:,:,iblk), & ! LCOV_EXCL_LINE
576 : strairxN (:,:,iblk), strairyN (:,:,iblk), & ! LCOV_EXCL_LINE
577 : ss_tltxN (:,:,iblk), ss_tltyN (:,:,iblk), & ! LCOV_EXCL_LINE
578 : iceTmask (:,:,iblk), iceNmask (:,:,iblk), & ! LCOV_EXCL_LINE
579 : fmN (:,:,iblk), dt , & ! LCOV_EXCL_LINE
580 : strtltxN (:,:,iblk), strtltyN (:,:,iblk), & ! LCOV_EXCL_LINE
581 : strocnxN (:,:,iblk), strocnyN (:,:,iblk), & ! LCOV_EXCL_LINE
582 : strintxN (:,:,iblk), strintyN (:,:,iblk), & ! LCOV_EXCL_LINE
583 : taubxN (:,:,iblk), taubyN (:,:,iblk), & ! LCOV_EXCL_LINE
584 : waterxN (:,:,iblk), wateryN (:,:,iblk), & ! LCOV_EXCL_LINE
585 : forcexN (:,:,iblk), forceyN (:,:,iblk), & ! LCOV_EXCL_LINE
586 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
587 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
588 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
589 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
590 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
591 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
592 : uvelN_init(:,:,iblk), vvelN_init(:,:,iblk), & ! LCOV_EXCL_LINE
593 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
594 0 : TbN (:,:,iblk))
595 :
596 : !-----------------------------------------------------------------
597 : ! more preparation for dynamics on E grid
598 : !-----------------------------------------------------------------
599 :
600 : call dyn_prep2 (nx_block, ny_block, &
601 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
602 : icellT (iblk), icellE (iblk), & ! LCOV_EXCL_LINE
603 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
604 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
605 : aiE (:,:,iblk), emass (:,:,iblk), & ! LCOV_EXCL_LINE
606 : emassdti (:,:,iblk), fcorE_blk (:,:,iblk), & ! LCOV_EXCL_LINE
607 : emask (:,:,iblk), & ! LCOV_EXCL_LINE
608 : uocnE (:,:,iblk), vocnE (:,:,iblk), & ! LCOV_EXCL_LINE
609 : strairxE (:,:,iblk), strairyE (:,:,iblk), & ! LCOV_EXCL_LINE
610 : ss_tltxE (:,:,iblk), ss_tltyE (:,:,iblk), & ! LCOV_EXCL_LINE
611 : iceTmask (:,:,iblk), iceEmask (:,:,iblk), & ! LCOV_EXCL_LINE
612 : fmE (:,:,iblk), dt , & ! LCOV_EXCL_LINE
613 : strtltxE (:,:,iblk), strtltyE (:,:,iblk), & ! LCOV_EXCL_LINE
614 : strocnxE (:,:,iblk), strocnyE (:,:,iblk), & ! LCOV_EXCL_LINE
615 : strintxE (:,:,iblk), strintyE (:,:,iblk), & ! LCOV_EXCL_LINE
616 : taubxE (:,:,iblk), taubyE (:,:,iblk), & ! LCOV_EXCL_LINE
617 : waterxE (:,:,iblk), wateryE (:,:,iblk), & ! LCOV_EXCL_LINE
618 : forcexE (:,:,iblk), forceyE (:,:,iblk), & ! LCOV_EXCL_LINE
619 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
620 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
621 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
622 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
623 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
624 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
625 : uvelE_init(:,:,iblk), vvelE_init(:,:,iblk), & ! LCOV_EXCL_LINE
626 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
627 0 : TbE (:,:,iblk))
628 :
629 :
630 0 : do i=1,nx_block
631 0 : do j=1,ny_block
632 0 : if (.not.iceUmask(i,j,iblk)) then
633 0 : stresspU (i,j,iblk) = c0
634 0 : stressmU (i,j,iblk) = c0
635 0 : stress12U(i,j,iblk) = c0
636 : endif
637 0 : if (.not.iceTmask(i,j,iblk)) then
638 0 : stresspT (i,j,iblk) = c0
639 0 : stressmT (i,j,iblk) = c0
640 0 : stress12T(i,j,iblk) = c0
641 : endif
642 : enddo
643 : enddo
644 : enddo ! iblk
645 : !$OMP END PARALLEL DO
646 :
647 : endif ! grid_ice
648 :
649 5784 : call icepack_warnings_flush(nu_diag)
650 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
651 0 : file=__FILE__, line=__LINE__)
652 :
653 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
654 :
655 0 : call ice_timer_start(timer_bound)
656 : call ice_HaloUpdate (uvelE, halo_info, &
657 0 : field_loc_Eface, field_type_vector)
658 : call ice_HaloUpdate (vvelN, halo_info, &
659 0 : field_loc_Nface, field_type_vector)
660 0 : call ice_timer_stop(timer_bound)
661 :
662 0 : if (grid_ice == 'C') then
663 0 : call grid_average_X2Y('A', uvelE, 'E', uvelN, 'N')
664 0 : call grid_average_X2Y('A', vvelN, 'N', vvelE, 'E')
665 0 : uvelN(:,:,:) = uvelN(:,:,:)*npm(:,:,:)
666 0 : vvelE(:,:,:) = vvelE(:,:,:)*epm(:,:,:)
667 : endif
668 :
669 0 : call ice_timer_start(timer_bound)
670 : call ice_HaloUpdate (uvelN, halo_info, &
671 0 : field_loc_Nface, field_type_vector)
672 : call ice_HaloUpdate (vvelE, halo_info, &
673 0 : field_loc_Eface, field_type_vector)
674 0 : call ice_timer_stop(timer_bound)
675 :
676 0 : call grid_average_X2Y('S', uvelE, 'E', uvel, 'U')
677 0 : call grid_average_X2Y('S', vvelN, 'N', vvel, 'U')
678 0 : uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:)
679 0 : vvel(:,:,:) = vvel(:,:,:)*uvm(:,:,:)
680 : endif
681 :
682 5784 : call ice_timer_start(timer_bound)
683 : call ice_HaloUpdate (strength, halo_info, &
684 5784 : field_loc_center, field_type_scalar)
685 :
686 : ! velocities may have changed in dyn_prep2
687 5784 : call stack_fields(uvel, vvel, fld2)
688 : call ice_HaloUpdate (fld2, halo_info, &
689 5784 : field_loc_NEcorner, field_type_vector)
690 5784 : call unstack_fields(fld2, uvel, vvel)
691 5784 : call ice_timer_stop(timer_bound)
692 :
693 5784 : if (maskhalo_dyn) then
694 0 : halomask = 0
695 0 : if (grid_ice == 'B') then
696 0 : where (iceUmask) halomask = 1
697 0 : elseif (grid_ice == 'C' .or. grid_ice == 'CD') then
698 0 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block,i,j) SCHEDULE(runtime)
699 0 : do iblk = 1, nblocks
700 0 : this_block = get_block(blocks_ice(iblk),iblk)
701 0 : ilo = this_block%ilo
702 0 : ihi = this_block%ihi
703 0 : jlo = this_block%jlo
704 0 : jhi = this_block%jhi
705 0 : do j = jlo,jhi
706 0 : do i = ilo,ihi
707 : if (iceTmask(i ,j ,iblk) .or. &
708 : iceTmask(i-1,j ,iblk) .or. & ! LCOV_EXCL_LINE
709 : iceTmask(i+1,j ,iblk) .or. & ! LCOV_EXCL_LINE
710 : iceTmask(i ,j-1,iblk) .or. & ! LCOV_EXCL_LINE
711 0 : iceTmask(i ,j+1,iblk)) then
712 0 : halomask(i,j,iblk) = 1
713 : endif
714 : enddo
715 : enddo
716 : enddo
717 : !$OMP END PARALLEL DO
718 : endif
719 0 : call ice_timer_start(timer_bound)
720 0 : call ice_HaloUpdate (halomask, halo_info, &
721 0 : field_loc_center, field_type_scalar)
722 0 : call ice_timer_stop(timer_bound)
723 0 : call ice_HaloMask(halo_info_mask, halo_info, halomask)
724 : endif
725 :
726 : !-----------------------------------------------------------------
727 : ! seabed stress factor TbU (TbU is part of Cb coefficient)
728 : !-----------------------------------------------------------------
729 :
730 5784 : if (seabed_stress) then
731 :
732 0 : if (grid_ice == "B") then
733 :
734 0 : if ( seabed_stress_method == 'LKD' ) then
735 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
736 0 : do iblk = 1, nblocks
737 : call seabed_stress_factor_LKD (nx_block , ny_block , &
738 : icellU (iblk), & ! LCOV_EXCL_LINE
739 : indxUi (:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
740 : vice (:,:,iblk), aice(:,:,iblk), & ! LCOV_EXCL_LINE
741 0 : hwater(:,:,iblk), TbU (:,:,iblk))
742 : enddo
743 : !$OMP END PARALLEL DO
744 :
745 0 : elseif ( seabed_stress_method == 'probabilistic' ) then
746 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
747 0 : do iblk = 1, nblocks
748 : call seabed_stress_factor_prob (nx_block , ny_block , &
749 : icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & ! LCOV_EXCL_LINE
750 : icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
751 : aicen(:,:,:,iblk), vicen(:,:,:,iblk) , & ! LCOV_EXCL_LINE
752 0 : hwater (:,:,iblk), TbU (:,:,iblk))
753 : enddo
754 : !$OMP END PARALLEL DO
755 : endif
756 :
757 0 : elseif (grid_ice == "C" .or. grid_ice == "CD") then
758 :
759 0 : if ( seabed_stress_method == 'LKD' ) then
760 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
761 0 : do iblk = 1, nblocks
762 : call seabed_stress_factor_LKD (nx_block , ny_block, &
763 : icellE (iblk), & ! LCOV_EXCL_LINE
764 : indxEi (:,iblk), indxEj(:,iblk), & ! LCOV_EXCL_LINE
765 : vice (:,:,iblk), aice(:,:,iblk), & ! LCOV_EXCL_LINE
766 0 : hwater(:,:,iblk), TbE (:,:,iblk))
767 : call seabed_stress_factor_LKD (nx_block , ny_block, &
768 : icellN (iblk), & ! LCOV_EXCL_LINE
769 : indxNi (:,iblk), indxNj(:,iblk), & ! LCOV_EXCL_LINE
770 : vice (:,:,iblk), aice(:,:,iblk), & ! LCOV_EXCL_LINE
771 0 : hwater(:,:,iblk), TbN (:,:,iblk))
772 : enddo
773 : !$OMP END PARALLEL DO
774 :
775 0 : elseif ( seabed_stress_method == 'probabilistic' ) then
776 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
777 0 : do iblk = 1, nblocks
778 : call seabed_stress_factor_prob (nx_block , ny_block , &
779 : icellT(iblk), indxTi(:,iblk), indxTj(:,iblk), & ! LCOV_EXCL_LINE
780 : icellU(iblk), indxUi(:,iblk), indxUj(:,iblk), & ! LCOV_EXCL_LINE
781 : aicen(:,:,:,iblk), vicen(:,:,:,iblk) , & ! LCOV_EXCL_LINE
782 : hwater (:,:,iblk), TbU (:,:,iblk) , & ! LCOV_EXCL_LINE
783 : TbE (:,:,iblk), TbN (:,:,iblk) , & ! LCOV_EXCL_LINE
784 : icellE(iblk), indxEi(:,iblk), indxEj(:,iblk), & ! LCOV_EXCL_LINE
785 0 : icellN(iblk), indxNi(:,iblk), indxNj(:,iblk) )
786 : enddo
787 : !$OMP END PARALLEL DO
788 : endif
789 :
790 : endif
791 :
792 : endif
793 :
794 5784 : if (evp_algorithm == "shared_mem_1d" ) then
795 :
796 0 : if (trim(grid_type) == 'tripole') then
797 : call abort_ice(trim(subname)//' &
798 : & Kernel not tested on tripole grid. Set evp_algorithm=standard_2d') ! LCOV_EXCL_LINE
799 : endif
800 :
801 0 : call ice_timer_start(timer_evp_1d)
802 : call ice_dyn_evp_1d_copyin( &
803 : nx_block,ny_block,nblocks,nx_global+2*nghost,ny_global+2*nghost, & ! LCOV_EXCL_LINE
804 : iceTmask, iceUmask, & ! LCOV_EXCL_LINE
805 : cdn_ocnU,aiU,uocnU,vocnU,forcexU,forceyU,TbU, & ! LCOV_EXCL_LINE
806 : umassdti,fmU,uarear,tarear,strintxU,strintyU,uvel_init,vvel_init,& ! LCOV_EXCL_LINE
807 : strength,uvel,vvel,dxT,dyT, & ! LCOV_EXCL_LINE
808 : stressp_1 ,stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
809 : stressm_1 ,stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
810 0 : stress12_1,stress12_2,stress12_3,stress12_4 )
811 0 : call ice_dyn_evp_1d_kernel()
812 : call ice_dyn_evp_1d_copyout( &
813 : nx_block,ny_block,nblocks,nx_global+2*nghost,ny_global+2*nghost, & ! LCOV_EXCL_LINE
814 : !strocn uvel,vvel, strocnxU,strocnyU, strintxU,strintyU, & ! LCOV_EXCL_LINE
815 : uvel,vvel, strintxU,strintyU, & ! LCOV_EXCL_LINE
816 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
817 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
818 : stress12_1,stress12_2,stress12_3,stress12_4, & ! LCOV_EXCL_LINE
819 0 : divu,rdg_conv,rdg_shear,shear,taubxU,taubyU )
820 0 : call ice_timer_stop(timer_evp_1d)
821 :
822 : else ! evp_algorithm == standard_2d (Standard CICE)
823 :
824 5784 : call ice_timer_start(timer_evp_2d)
825 :
826 5784 : if (grid_ice == "B") then
827 :
828 1393944 : do ksub = 1,ndte ! subcycling
829 :
830 691200 : !$OMP PARALLEL DO PRIVATE(iblk,strtmp) SCHEDULE(runtime)
831 2085120 : do iblk = 1, nblocks
832 :
833 : !-----------------------------------------------------------------
834 : ! stress tensor equation, total surface stress
835 : !-----------------------------------------------------------------
836 : call stress (nx_block , ny_block , &
837 : icellT (iblk), & ! LCOV_EXCL_LINE
838 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
839 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
840 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
841 : dxhy (:,:,iblk), dyhx (:,:,iblk), & ! LCOV_EXCL_LINE
842 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
843 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
844 : DminTarea (:,:,iblk), & ! LCOV_EXCL_LINE
845 : strength (:,:,iblk), & ! LCOV_EXCL_LINE
846 : stressp_1 (:,:,iblk), stressp_2 (:,:,iblk), & ! LCOV_EXCL_LINE
847 : stressp_3 (:,:,iblk), stressp_4 (:,:,iblk), & ! LCOV_EXCL_LINE
848 : stressm_1 (:,:,iblk), stressm_2 (:,:,iblk), & ! LCOV_EXCL_LINE
849 : stressm_3 (:,:,iblk), stressm_4 (:,:,iblk), & ! LCOV_EXCL_LINE
850 : stress12_1(:,:,iblk), stress12_2(:,:,iblk), & ! LCOV_EXCL_LINE
851 : stress12_3(:,:,iblk), stress12_4(:,:,iblk), & ! LCOV_EXCL_LINE
852 1388160 : strtmp (:,:,:) )
853 :
854 : !-----------------------------------------------------------------
855 : ! momentum equation
856 : !-----------------------------------------------------------------
857 : call stepu (nx_block , ny_block , &
858 : icellU (iblk), Cdn_ocnU(:,:,iblk), & ! LCOV_EXCL_LINE
859 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
860 : aiU (:,:,iblk), strtmp (:,:,:), & ! LCOV_EXCL_LINE
861 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
862 : waterxU (:,:,iblk), wateryU (:,:,iblk), & ! LCOV_EXCL_LINE
863 : forcexU (:,:,iblk), forceyU (:,:,iblk), & ! LCOV_EXCL_LINE
864 : umassdti (:,:,iblk), fmU (:,:,iblk), & ! LCOV_EXCL_LINE
865 : uarear (:,:,iblk), & ! LCOV_EXCL_LINE
866 : strintxU (:,:,iblk), strintyU(:,:,iblk), & ! LCOV_EXCL_LINE
867 : taubxU (:,:,iblk), taubyU (:,:,iblk), & ! LCOV_EXCL_LINE
868 : uvel_init(:,:,iblk), vvel_init(:,:,iblk),& ! LCOV_EXCL_LINE
869 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
870 2776320 : TbU (:,:,iblk))
871 :
872 : enddo ! iblk
873 : !$OMP END PARALLEL DO
874 :
875 : ! U fields at NE corner
876 : ! calls ice_haloUpdate, controls bundles and masks
877 : call dyn_haloUpdate (halo_info, halo_info_mask, &
878 : field_loc_NEcorner, field_type_vector, & ! LCOV_EXCL_LINE
879 1393944 : uvel, vvel)
880 :
881 : enddo ! sub cycling
882 :
883 : !-----------------------------------------------------------------
884 : ! save quantities for mechanical redistribution
885 : !-----------------------------------------------------------------
886 :
887 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
888 8688 : do iblk = 1, nblocks
889 : call deformations (nx_block , ny_block , &
890 : icellT (iblk), & ! LCOV_EXCL_LINE
891 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
892 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
893 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
894 : cxp (:,:,iblk), cyp (:,:,iblk), & ! LCOV_EXCL_LINE
895 : cxm (:,:,iblk), cym (:,:,iblk), & ! LCOV_EXCL_LINE
896 : tarear (:,:,iblk), & ! LCOV_EXCL_LINE
897 : shear (:,:,iblk), divu (:,:,iblk), & ! LCOV_EXCL_LINE
898 11568 : rdg_conv(:,:,iblk), rdg_shear(:,:,iblk) )
899 : enddo
900 : !$OMP END PARALLEL DO
901 :
902 :
903 0 : elseif (grid_ice == "C") then
904 :
905 0 : do ksub = 1,ndte ! subcycling
906 :
907 0 : !$OMP PARALLEL DO PRIVATE(iblk)
908 0 : do iblk = 1, nblocks
909 :
910 : !-----------------------------------------------------------------
911 : ! strain rates at U point
912 : ! NOTE these are actually strain rates * area (m^2/s)
913 : !-----------------------------------------------------------------
914 : call strain_rates_U (nx_block , ny_block , &
915 : icellU (iblk), & ! LCOV_EXCL_LINE
916 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
917 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
918 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
919 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
920 : dxE (:,:,iblk), dyN (:,:,iblk), & ! LCOV_EXCL_LINE
921 : dxU (:,:,iblk), dyU (:,:,iblk), & ! LCOV_EXCL_LINE
922 : ratiodxN(:,:,iblk), ratiodxNr(:,:,iblk), & ! LCOV_EXCL_LINE
923 : ratiodyE(:,:,iblk), ratiodyEr(:,:,iblk), & ! LCOV_EXCL_LINE
924 : epm (:,:,iblk), npm (:,:,iblk), & ! LCOV_EXCL_LINE
925 : divergU (:,:,iblk), tensionU (:,:,iblk), & ! LCOV_EXCL_LINE
926 0 : shearU (:,:,iblk), deltaU (:,:,iblk) )
927 :
928 : enddo ! iblk
929 : !$OMP END PARALLEL DO
930 :
931 : ! calls ice_haloUpdate, controls bundles and masks
932 : call dyn_haloUpdate (halo_info, halo_info_mask, &
933 : field_loc_NEcorner, field_type_scalar, & ! LCOV_EXCL_LINE
934 0 : shearU)
935 :
936 0 : !$OMP PARALLEL DO PRIVATE(iblk)
937 0 : do iblk = 1, nblocks
938 : call stressC_T (nx_block , ny_block , &
939 : icellT (iblk), & ! LCOV_EXCL_LINE
940 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
941 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
942 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
943 : dxN (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
944 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
945 : uarea (:,:,iblk), DminTarea (:,:,iblk), & ! LCOV_EXCL_LINE
946 : strength (:,:,iblk), shearU (:,:,iblk), & ! LCOV_EXCL_LINE
947 : zetax2T (:,:,iblk), etax2T (:,:,iblk), & ! LCOV_EXCL_LINE
948 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
949 0 : stress12T (:,:,iblk))
950 :
951 : enddo
952 : !$OMP END PARALLEL DO
953 :
954 : ! calls ice_haloUpdate, controls bundles and masks
955 : call dyn_haloUpdate (halo_info, halo_info_mask, &
956 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
957 0 : zetax2T, etax2T, stresspT, stressmT)
958 :
959 0 : if (visc_method == 'avg_strength') then
960 0 : call grid_average_X2Y('S', strength, 'T', strengthU, 'U')
961 0 : elseif (visc_method == 'avg_zeta') then
962 0 : call grid_average_X2Y('S', etax2T , 'T', etax2U , 'U')
963 : endif
964 :
965 0 : !$OMP PARALLEL DO PRIVATE(iblk)
966 0 : do iblk = 1, nblocks
967 : call stressC_U (nx_block , ny_block , &
968 : icellU (iblk), & ! LCOV_EXCL_LINE
969 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
970 : uarea (:,:,iblk), & ! LCOV_EXCL_LINE
971 : etax2U (:,:,iblk), deltaU (:,:,iblk), & ! LCOV_EXCL_LINE
972 : strengthU (:,:,iblk), shearU (:,:,iblk), & ! LCOV_EXCL_LINE
973 0 : stress12U (:,:,iblk))
974 : enddo
975 : !$OMP END PARALLEL DO
976 :
977 : ! calls ice_haloUpdate, controls bundles and masks
978 : call dyn_haloUpdate (halo_info , halo_info_mask, &
979 : field_loc_NEcorner, field_type_scalar, & ! LCOV_EXCL_LINE
980 0 : stress12U)
981 :
982 0 : !$OMP PARALLEL DO PRIVATE(iblk)
983 0 : do iblk = 1, nblocks
984 :
985 : call div_stress_Ex (nx_block , ny_block , &
986 : icellE (iblk), & ! LCOV_EXCL_LINE
987 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
988 : dxE (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
989 : dxU (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
990 : earear (:,:,iblk) , & ! LCOV_EXCL_LINE
991 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
992 0 : stress12U (:,:,iblk), strintxE (:,:,iblk) )
993 :
994 : call div_stress_Ny (nx_block , ny_block , &
995 : icellN (iblk), & ! LCOV_EXCL_LINE
996 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
997 : dxN (:,:,iblk), dyN (:,:,iblk), & ! LCOV_EXCL_LINE
998 : dxT (:,:,iblk), dyU (:,:,iblk), & ! LCOV_EXCL_LINE
999 : narear (:,:,iblk) , & ! LCOV_EXCL_LINE
1000 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
1001 0 : stress12U (:,:,iblk), strintyN (:,:,iblk) )
1002 :
1003 : enddo
1004 : !$OMP END PARALLEL DO
1005 :
1006 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1007 0 : do iblk = 1, nblocks
1008 :
1009 : call stepu_C (nx_block , ny_block , & ! u, E point
1010 : icellE (iblk), Cdn_ocnE (:,:,iblk), & ! LCOV_EXCL_LINE
1011 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
1012 : aiE (:,:,iblk), & ! LCOV_EXCL_LINE
1013 : uocnE (:,:,iblk), vocnE (:,:,iblk), & ! LCOV_EXCL_LINE
1014 : waterxE (:,:,iblk), forcexE (:,:,iblk), & ! LCOV_EXCL_LINE
1015 : emassdti (:,:,iblk), fmE (:,:,iblk), & ! LCOV_EXCL_LINE
1016 : strintxE (:,:,iblk), taubxE (:,:,iblk), & ! LCOV_EXCL_LINE
1017 : uvelE_init(:,:,iblk), & ! LCOV_EXCL_LINE
1018 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1019 0 : TbE (:,:,iblk))
1020 :
1021 : call stepv_C (nx_block, ny_block, & ! v, N point
1022 : icellN (iblk), Cdn_ocnN (:,:,iblk), & ! LCOV_EXCL_LINE
1023 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
1024 : aiN (:,:,iblk), & ! LCOV_EXCL_LINE
1025 : uocnN (:,:,iblk), vocnN (:,:,iblk), & ! LCOV_EXCL_LINE
1026 : wateryN (:,:,iblk), forceyN (:,:,iblk), & ! LCOV_EXCL_LINE
1027 : nmassdti (:,:,iblk), fmN (:,:,iblk), & ! LCOV_EXCL_LINE
1028 : strintyN (:,:,iblk), taubyN (:,:,iblk), & ! LCOV_EXCL_LINE
1029 : vvelN_init(:,:,iblk), & ! LCOV_EXCL_LINE
1030 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1031 0 : TbN (:,:,iblk))
1032 : enddo
1033 : !$OMP END PARALLEL DO
1034 :
1035 : ! calls ice_haloUpdate, controls bundles and masks
1036 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1037 : field_loc_Eface, field_type_vector, & ! LCOV_EXCL_LINE
1038 0 : uvelE)
1039 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1040 : field_loc_Nface, field_type_vector, & ! LCOV_EXCL_LINE
1041 0 : vvelN)
1042 :
1043 0 : call grid_average_X2Y('A', uvelE, 'E', uvelN, 'N')
1044 0 : call grid_average_X2Y('A', vvelN, 'N', vvelE, 'E')
1045 0 : uvelN(:,:,:) = uvelN(:,:,:)*npm(:,:,:)
1046 0 : vvelE(:,:,:) = vvelE(:,:,:)*epm(:,:,:)
1047 :
1048 : ! calls ice_haloUpdate, controls bundles and masks
1049 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1050 : field_loc_Nface, field_type_vector, & ! LCOV_EXCL_LINE
1051 0 : uvelN)
1052 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1053 : field_loc_Eface, field_type_vector, & ! LCOV_EXCL_LINE
1054 0 : vvelE)
1055 :
1056 0 : call grid_average_X2Y('S', uvelE, 'E', uvel, 'U')
1057 0 : call grid_average_X2Y('S', vvelN, 'N', vvel, 'U')
1058 :
1059 0 : uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:)
1060 0 : vvel(:,:,:) = vvel(:,:,:)*uvm(:,:,:)
1061 : ! U fields at NE corner
1062 : ! calls ice_haloUpdate, controls bundles and masks
1063 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1064 : field_loc_NEcorner, field_type_vector, & ! LCOV_EXCL_LINE
1065 0 : uvel, vvel)
1066 :
1067 : enddo ! subcycling
1068 :
1069 : !-----------------------------------------------------------------
1070 : ! save quantities for mechanical redistribution
1071 : !-----------------------------------------------------------------
1072 :
1073 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
1074 0 : do iblk = 1, nblocks
1075 : call deformationsC_T (nx_block , ny_block , &
1076 : icellT (iblk), & ! LCOV_EXCL_LINE
1077 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
1078 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1079 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1080 : dxN (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
1081 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
1082 : tarear (:,:,iblk), uarea (:,:,iblk), & ! LCOV_EXCL_LINE
1083 : shearU (:,:,iblk), & ! LCOV_EXCL_LINE
1084 : shear (:,:,iblk), divu (:,:,iblk), & ! LCOV_EXCL_LINE
1085 0 : rdg_conv(:,:,iblk), rdg_shear(:,:,iblk))
1086 : enddo
1087 : !$OMP END PARALLEL DO
1088 :
1089 0 : elseif (grid_ice == "CD") then
1090 :
1091 0 : do ksub = 1,ndte ! subcycling
1092 :
1093 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1094 0 : do iblk = 1, nblocks
1095 : call stressCD_T (nx_block , ny_block , &
1096 : icellT (iblk), & ! LCOV_EXCL_LINE
1097 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
1098 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1099 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1100 : dxN (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
1101 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
1102 : DminTarea (:,:,iblk), & ! LCOV_EXCL_LINE
1103 : strength (:,:,iblk), & ! LCOV_EXCL_LINE
1104 : zetax2T (:,:,iblk), etax2T (:,:,iblk), & ! LCOV_EXCL_LINE
1105 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
1106 0 : stress12T (:,:,iblk) )
1107 :
1108 : enddo
1109 : !$OMP END PARALLEL DO
1110 :
1111 : ! calls ice_haloUpdate, controls bundles and masks
1112 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1113 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
1114 0 : zetax2T, etax2T)
1115 :
1116 0 : if (visc_method == 'avg_strength') then
1117 0 : call grid_average_X2Y('S', strength, 'T', strengthU, 'U')
1118 0 : elseif (visc_method == 'avg_zeta') then
1119 0 : call grid_average_X2Y('S', zetax2T , 'T', zetax2U , 'U')
1120 0 : call grid_average_X2Y('S', etax2T , 'T', etax2U , 'U')
1121 : endif
1122 :
1123 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1124 0 : do iblk = 1, nblocks
1125 : !-----------------------------------------------------------------
1126 : ! strain rates at U point
1127 : ! NOTE these are actually strain rates * area (m^2/s)
1128 : !-----------------------------------------------------------------
1129 : call strain_rates_U (nx_block , ny_block , &
1130 : icellU (iblk), & ! LCOV_EXCL_LINE
1131 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
1132 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1133 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1134 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
1135 : dxE (:,:,iblk), dyN (:,:,iblk), & ! LCOV_EXCL_LINE
1136 : dxU (:,:,iblk), dyU (:,:,iblk), & ! LCOV_EXCL_LINE
1137 : ratiodxN (:,:,iblk), ratiodxNr(:,:,iblk), & ! LCOV_EXCL_LINE
1138 : ratiodyE (:,:,iblk), ratiodyEr(:,:,iblk), & ! LCOV_EXCL_LINE
1139 : epm (:,:,iblk), npm (:,:,iblk), & ! LCOV_EXCL_LINE
1140 : divergU (:,:,iblk), tensionU (:,:,iblk), & ! LCOV_EXCL_LINE
1141 0 : shearU (:,:,iblk), DeltaU (:,:,iblk) )
1142 :
1143 : call stressCD_U (nx_block , ny_block , &
1144 : icellU (iblk), & ! LCOV_EXCL_LINE
1145 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
1146 : uarea (:,:,iblk), & ! LCOV_EXCL_LINE
1147 : zetax2U (:,:,iblk), etax2U (:,:,iblk), & ! LCOV_EXCL_LINE
1148 : strengthU(:,:,iblk), & ! LCOV_EXCL_LINE
1149 : divergU (:,:,iblk), tensionU (:,:,iblk), & ! LCOV_EXCL_LINE
1150 : shearU (:,:,iblk), DeltaU (:,:,iblk), & ! LCOV_EXCL_LINE
1151 : stresspU (:,:,iblk), stressmU (:,:,iblk), & ! LCOV_EXCL_LINE
1152 0 : stress12U(:,:,iblk))
1153 : enddo
1154 : !$OMP END PARALLEL DO
1155 :
1156 : ! calls ice_haloUpdate, controls bundles and masks
1157 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1158 : field_loc_center, field_type_scalar, & ! LCOV_EXCL_LINE
1159 0 : stresspT, stressmT, stress12T)
1160 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1161 : field_loc_NEcorner,field_type_scalar, & ! LCOV_EXCL_LINE
1162 0 : stresspU, stressmU, stress12U)
1163 :
1164 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1165 0 : do iblk = 1, nblocks
1166 :
1167 : call div_stress_Ex (nx_block , ny_block , &
1168 : icellE (iblk), & ! LCOV_EXCL_LINE
1169 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
1170 : dxE (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
1171 : dxU (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
1172 : earear (:,:,iblk) , & ! LCOV_EXCL_LINE
1173 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
1174 0 : stress12U (:,:,iblk), strintxE (:,:,iblk) )
1175 :
1176 : call div_stress_Ey (nx_block , ny_block , &
1177 : icellE (iblk), & ! LCOV_EXCL_LINE
1178 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
1179 : dxE (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
1180 : dxU (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
1181 : earear (:,:,iblk) , & ! LCOV_EXCL_LINE
1182 : stresspU (:,:,iblk), stressmU (:,:,iblk), & ! LCOV_EXCL_LINE
1183 0 : stress12T (:,:,iblk), strintyE (:,:,iblk) )
1184 :
1185 : call div_stress_Nx (nx_block , ny_block , &
1186 : icellN (iblk), & ! LCOV_EXCL_LINE
1187 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
1188 : dxN (:,:,iblk), dyN (:,:,iblk), & ! LCOV_EXCL_LINE
1189 : dxT (:,:,iblk), dyU (:,:,iblk), & ! LCOV_EXCL_LINE
1190 : narear (:,:,iblk) , & ! LCOV_EXCL_LINE
1191 : stresspU (:,:,iblk), stressmU (:,:,iblk), & ! LCOV_EXCL_LINE
1192 0 : stress12T (:,:,iblk), strintxN (:,:,iblk) )
1193 :
1194 : call div_stress_Ny (nx_block , ny_block , &
1195 : icellN (iblk), & ! LCOV_EXCL_LINE
1196 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
1197 : dxN (:,:,iblk), dyN (:,:,iblk), & ! LCOV_EXCL_LINE
1198 : dxT (:,:,iblk), dyU (:,:,iblk), & ! LCOV_EXCL_LINE
1199 : narear (:,:,iblk) , & ! LCOV_EXCL_LINE
1200 : stresspT (:,:,iblk), stressmT (:,:,iblk), & ! LCOV_EXCL_LINE
1201 0 : stress12U (:,:,iblk), strintyN (:,:,iblk) )
1202 :
1203 : enddo
1204 : !$OMP END PARALLEL DO
1205 :
1206 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1207 0 : do iblk = 1, nblocks
1208 :
1209 : call stepuv_CD (nx_block , ny_block , & ! E point
1210 : icellE (iblk), Cdn_ocnE (:,:,iblk), & ! LCOV_EXCL_LINE
1211 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
1212 : aiE (:,:,iblk), & ! LCOV_EXCL_LINE
1213 : uocnE (:,:,iblk), vocnE (:,:,iblk), & ! LCOV_EXCL_LINE
1214 : waterxE (:,:,iblk), wateryE (:,:,iblk), & ! LCOV_EXCL_LINE
1215 : forcexE (:,:,iblk), forceyE (:,:,iblk), & ! LCOV_EXCL_LINE
1216 : emassdti (:,:,iblk), fmE (:,:,iblk), & ! LCOV_EXCL_LINE
1217 : strintxE (:,:,iblk), strintyE (:,:,iblk), & ! LCOV_EXCL_LINE
1218 : taubxE (:,:,iblk), taubyE (:,:,iblk), & ! LCOV_EXCL_LINE
1219 : uvelE_init(:,:,iblk), vvelE_init(:,:,iblk), & ! LCOV_EXCL_LINE
1220 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1221 0 : TbE (:,:,iblk))
1222 :
1223 : call stepuv_CD (nx_block , ny_block , & ! N point
1224 : icellN (iblk), Cdn_ocnN (:,:,iblk), & ! LCOV_EXCL_LINE
1225 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
1226 : aiN (:,:,iblk), & ! LCOV_EXCL_LINE
1227 : uocnN (:,:,iblk), vocnN (:,:,iblk), & ! LCOV_EXCL_LINE
1228 : waterxN (:,:,iblk), wateryN (:,:,iblk), & ! LCOV_EXCL_LINE
1229 : forcexN (:,:,iblk), forceyN (:,:,iblk), & ! LCOV_EXCL_LINE
1230 : nmassdti (:,:,iblk), fmN (:,:,iblk), & ! LCOV_EXCL_LINE
1231 : strintxN (:,:,iblk), strintyN (:,:,iblk), & ! LCOV_EXCL_LINE
1232 : taubxN (:,:,iblk), taubyN (:,:,iblk), & ! LCOV_EXCL_LINE
1233 : uvelN_init(:,:,iblk), vvelN_init(:,:,iblk), & ! LCOV_EXCL_LINE
1234 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1235 0 : TbN (:,:,iblk))
1236 : enddo
1237 : !$OMP END PARALLEL DO
1238 :
1239 : ! calls ice_haloUpdate, controls bundles and masks
1240 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1241 : field_loc_Eface, field_type_vector, & ! LCOV_EXCL_LINE
1242 0 : uvelE, vvelE)
1243 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1244 : field_loc_Nface, field_type_vector, & ! LCOV_EXCL_LINE
1245 0 : uvelN, vvelN)
1246 :
1247 0 : call grid_average_X2Y('S', uvelE, 'E', uvel, 'U')
1248 0 : call grid_average_X2Y('S', vvelN, 'N', vvel, 'U')
1249 :
1250 0 : uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:)
1251 0 : vvel(:,:,:) = vvel(:,:,:)*uvm(:,:,:)
1252 : ! U fields at NE corner
1253 : ! calls ice_haloUpdate, controls bundles and masks
1254 : call dyn_haloUpdate (halo_info, halo_info_mask, &
1255 : field_loc_NEcorner, field_type_vector, & ! LCOV_EXCL_LINE
1256 0 : uvel, vvel)
1257 :
1258 : enddo ! subcycling
1259 :
1260 : !-----------------------------------------------------------------
1261 : ! save quantities for mechanical redistribution
1262 : !-----------------------------------------------------------------
1263 :
1264 0 : !$OMP PARALLEL DO PRIVATE(iblk)
1265 0 : do iblk = 1, nblocks
1266 : call deformationsCD_T (nx_block , ny_block , &
1267 : icellT (iblk), & ! LCOV_EXCL_LINE
1268 : indxTi (:,iblk), indxTj (:,iblk), & ! LCOV_EXCL_LINE
1269 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1270 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1271 : dxN (:,:,iblk), dyE (:,:,iblk), & ! LCOV_EXCL_LINE
1272 : dxT (:,:,iblk), dyT (:,:,iblk), & ! LCOV_EXCL_LINE
1273 : tarear (:,:,iblk), & ! LCOV_EXCL_LINE
1274 : shear (:,:,iblk), divu (:,:,iblk), & ! LCOV_EXCL_LINE
1275 0 : rdg_conv(:,:,iblk), rdg_shear(:,:,iblk))
1276 : enddo
1277 : !$OMP END PARALLEL DO
1278 : endif ! grid_ice
1279 :
1280 5784 : call ice_timer_stop(timer_evp_2d)
1281 : endif ! evp_algorithm
1282 :
1283 5784 : if (maskhalo_dyn) then
1284 0 : call ice_HaloDestroy(halo_info_mask)
1285 : endif
1286 :
1287 : ! Force symmetry across the tripole seam
1288 5784 : if (trim(grid_type) == 'tripole') then
1289 : ! TODO: C/CD-grid
1290 0 : if (maskhalo_dyn) then
1291 : !-------------------------------------------------------
1292 : ! set halomask to zero because ice_HaloMask always keeps
1293 : ! local copies AND tripole zipper communication
1294 : !-------------------------------------------------------
1295 0 : halomask = 0
1296 0 : call ice_HaloMask(halo_info_mask, halo_info, halomask)
1297 :
1298 : call ice_HaloUpdate_stress(stressp_1 , stressp_3 , halo_info_mask, &
1299 0 : field_loc_center, field_type_scalar)
1300 : call ice_HaloUpdate_stress(stressp_3 , stressp_1 , halo_info_mask, &
1301 0 : field_loc_center, field_type_scalar)
1302 : call ice_HaloUpdate_stress(stressp_2 , stressp_4 , halo_info_mask, &
1303 0 : field_loc_center, field_type_scalar)
1304 : call ice_HaloUpdate_stress(stressp_4 , stressp_2 , halo_info_mask, &
1305 0 : field_loc_center, field_type_scalar)
1306 :
1307 : call ice_HaloUpdate_stress(stressm_1 , stressm_3 , halo_info_mask, &
1308 0 : field_loc_center, field_type_scalar)
1309 : call ice_HaloUpdate_stress(stressm_3 , stressm_1 , halo_info_mask, &
1310 0 : field_loc_center, field_type_scalar)
1311 : call ice_HaloUpdate_stress(stressm_2 , stressm_4 , halo_info_mask, &
1312 0 : field_loc_center, field_type_scalar)
1313 : call ice_HaloUpdate_stress(stressm_4 , stressm_2 , halo_info_mask, &
1314 0 : field_loc_center, field_type_scalar)
1315 :
1316 : call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info_mask, &
1317 0 : field_loc_center, field_type_scalar)
1318 : call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info_mask, &
1319 0 : field_loc_center, field_type_scalar)
1320 : call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info_mask, &
1321 0 : field_loc_center, field_type_scalar)
1322 : call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info_mask, &
1323 0 : field_loc_center, field_type_scalar)
1324 :
1325 0 : call ice_HaloDestroy(halo_info_mask)
1326 : else
1327 : call ice_HaloUpdate_stress(stressp_1 , stressp_3 , halo_info, &
1328 0 : field_loc_center, field_type_scalar)
1329 : call ice_HaloUpdate_stress(stressp_3 , stressp_1 , halo_info, &
1330 0 : field_loc_center, field_type_scalar)
1331 : call ice_HaloUpdate_stress(stressp_2 , stressp_4 , halo_info, &
1332 0 : field_loc_center, field_type_scalar)
1333 : call ice_HaloUpdate_stress(stressp_4 , stressp_2 , halo_info, &
1334 0 : field_loc_center, field_type_scalar)
1335 :
1336 : call ice_HaloUpdate_stress(stressm_1 , stressm_3 , halo_info, &
1337 0 : field_loc_center, field_type_scalar)
1338 : call ice_HaloUpdate_stress(stressm_3 , stressm_1 , halo_info, &
1339 0 : field_loc_center, field_type_scalar)
1340 : call ice_HaloUpdate_stress(stressm_2 , stressm_4 , halo_info, &
1341 0 : field_loc_center, field_type_scalar)
1342 : call ice_HaloUpdate_stress(stressm_4 , stressm_2 , halo_info, &
1343 0 : field_loc_center, field_type_scalar)
1344 :
1345 : call ice_HaloUpdate_stress(stress12_1, stress12_3, halo_info, &
1346 0 : field_loc_center, field_type_scalar)
1347 : call ice_HaloUpdate_stress(stress12_3, stress12_1, halo_info, &
1348 0 : field_loc_center, field_type_scalar)
1349 : call ice_HaloUpdate_stress(stress12_2, stress12_4, halo_info, &
1350 0 : field_loc_center, field_type_scalar)
1351 : call ice_HaloUpdate_stress(stress12_4, stress12_2, halo_info, &
1352 0 : field_loc_center, field_type_scalar)
1353 : endif ! maskhalo
1354 : endif ! tripole
1355 :
1356 : !-----------------------------------------------------------------
1357 : ! ice-ocean stress
1358 : !-----------------------------------------------------------------
1359 :
1360 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
1361 8688 : do iblk = 1, nblocks
1362 : call dyn_finish &
1363 : (nx_block , ny_block , & ! LCOV_EXCL_LINE
1364 : icellU (iblk), Cdn_ocnU(:,:,iblk), & ! LCOV_EXCL_LINE
1365 : indxUi (:,iblk), indxUj (:,iblk), & ! LCOV_EXCL_LINE
1366 : uvel (:,:,iblk), vvel (:,:,iblk), & ! LCOV_EXCL_LINE
1367 : uocnU (:,:,iblk), vocnU (:,:,iblk), & ! LCOV_EXCL_LINE
1368 : aiU (:,:,iblk), fmU (:,:,iblk), & ! LCOV_EXCL_LINE
1369 11568 : strocnxU(:,:,iblk), strocnyU(:,:,iblk))
1370 : enddo
1371 : !$OMP END PARALLEL DO
1372 :
1373 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
1374 :
1375 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
1376 0 : do iblk = 1, nblocks
1377 :
1378 : call dyn_finish &
1379 : (nx_block , ny_block , & ! LCOV_EXCL_LINE
1380 : icellN (iblk), Cdn_ocnN(:,:,iblk), & ! LCOV_EXCL_LINE
1381 : indxNi (:,iblk), indxNj (:,iblk), & ! LCOV_EXCL_LINE
1382 : uvelN (:,:,iblk), vvelN (:,:,iblk), & ! LCOV_EXCL_LINE
1383 : uocnN (:,:,iblk), vocnN (:,:,iblk), & ! LCOV_EXCL_LINE
1384 : aiN (:,:,iblk), fmN (:,:,iblk), & ! LCOV_EXCL_LINE
1385 0 : strocnxN(:,:,iblk), strocnyN(:,:,iblk))
1386 :
1387 : call dyn_finish &
1388 : (nx_block , ny_block , & ! LCOV_EXCL_LINE
1389 : icellE (iblk), Cdn_ocnE(:,:,iblk), & ! LCOV_EXCL_LINE
1390 : indxEi (:,iblk), indxEj (:,iblk), & ! LCOV_EXCL_LINE
1391 : uvelE (:,:,iblk), vvelE (:,:,iblk), & ! LCOV_EXCL_LINE
1392 : uocnE (:,:,iblk), vocnE (:,:,iblk), & ! LCOV_EXCL_LINE
1393 : aiE (:,:,iblk), fmE (:,:,iblk), & ! LCOV_EXCL_LINE
1394 0 : strocnxE(:,:,iblk), strocnyE(:,:,iblk))
1395 :
1396 : enddo
1397 : !$OMP END PARALLEL DO
1398 :
1399 : endif
1400 :
1401 5784 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
1402 : call ice_HaloUpdate (strintxE, halo_info, &
1403 0 : field_loc_Eface, field_type_vector)
1404 : call ice_HaloUpdate (strintyN, halo_info, &
1405 0 : field_loc_Nface, field_type_vector)
1406 0 : call ice_timer_stop(timer_bound)
1407 0 : call grid_average_X2Y('S', strintxE, 'E', strintxU, 'U') ! diagnostic
1408 0 : call grid_average_X2Y('S', strintyN, 'N', strintyU, 'U') ! diagnostic
1409 : endif
1410 :
1411 5784 : call ice_timer_stop(timer_dynamics) ! dynamics
1412 :
1413 5784 : end subroutine evp
1414 :
1415 : !=======================================================================
1416 : ! Computes the rates of strain and internal stress components for
1417 : ! each of the four corners on each T-grid cell.
1418 : ! Computes stress terms for the momentum equation
1419 : !
1420 : ! author: Elizabeth C. Hunke, LANL
1421 :
1422 5506560 : subroutine stress (nx_block, ny_block, &
1423 : icellT, & ! LCOV_EXCL_LINE
1424 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1425 : uvel, vvel, & ! LCOV_EXCL_LINE
1426 : dxT, dyT, & ! LCOV_EXCL_LINE
1427 : dxhy, dyhx, & ! LCOV_EXCL_LINE
1428 : cxp, cyp, & ! LCOV_EXCL_LINE
1429 : cxm, cym, & ! LCOV_EXCL_LINE
1430 : DminTarea, & ! LCOV_EXCL_LINE
1431 : strength, & ! LCOV_EXCL_LINE
1432 : stressp_1, stressp_2, & ! LCOV_EXCL_LINE
1433 : stressp_3, stressp_4, & ! LCOV_EXCL_LINE
1434 : stressm_1, stressm_2, & ! LCOV_EXCL_LINE
1435 : stressm_3, stressm_4, & ! LCOV_EXCL_LINE
1436 : stress12_1, stress12_2, & ! LCOV_EXCL_LINE
1437 : stress12_3, stress12_4, & ! LCOV_EXCL_LINE
1438 5506560 : str )
1439 :
1440 : use ice_dyn_shared, only: strain_rates, visc_replpress, capping
1441 :
1442 : integer (kind=int_kind), intent(in) :: &
1443 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1444 : icellT ! no. of cells where iceTmask = .true.
1445 :
1446 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1447 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1448 : indxTj ! compressed index in j-direction
1449 :
1450 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1451 : strength , & ! ice strength (N/m) ! LCOV_EXCL_LINE
1452 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1453 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1454 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1455 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1456 : dxhy , & ! 0.5*(HTE - HTW) ! LCOV_EXCL_LINE
1457 : dyhx , & ! 0.5*(HTN - HTS) ! LCOV_EXCL_LINE
1458 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1459 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1460 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1461 : cxm , & ! 0.5*HTN - 1.5*HTS ! LCOV_EXCL_LINE
1462 : DminTarea ! deltaminEVP*tarea
1463 :
1464 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1465 : stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1466 : stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1467 : stress12_1,stress12_2,stress12_3,stress12_4 ! sigma12
1468 :
1469 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(out) :: &
1470 : str ! stress combinations
1471 :
1472 : ! local variables
1473 :
1474 : integer (kind=int_kind) :: &
1475 : i, j, ij
1476 :
1477 : real (kind=dbl_kind) :: &
1478 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1479 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1480 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1481 : Deltane, Deltanw, Deltase, Deltasw , & ! Delt ! LCOV_EXCL_LINE
1482 : zetax2ne, zetax2nw, zetax2se, zetax2sw , & ! 2 x zeta (bulk visc) ! LCOV_EXCL_LINE
1483 : etax2ne, etax2nw, etax2se, etax2sw , & ! 2 x eta (shear visc) ! LCOV_EXCL_LINE
1484 : rep_prsne, rep_prsnw, rep_prsse, rep_prssw, & ! replacement pressure ! LCOV_EXCL_LINE
1485 : ! puny , & ! puny ! LCOV_EXCL_LINE
1486 : ssigpn, ssigps, ssigpe, ssigpw , & ! LCOV_EXCL_LINE
1487 : ssigmn, ssigms, ssigme, ssigmw , & ! LCOV_EXCL_LINE
1488 : ssig12n, ssig12s, ssig12e, ssig12w , & ! LCOV_EXCL_LINE
1489 : ssigp1, ssigp2, ssigm1, ssigm2, ssig121, ssig122, & ! LCOV_EXCL_LINE
1490 : csigpne, csigpnw, csigpse, csigpsw , & ! LCOV_EXCL_LINE
1491 : csigmne, csigmnw, csigmse, csigmsw , & ! LCOV_EXCL_LINE
1492 : csig12ne, csig12nw, csig12se, csig12sw , & ! LCOV_EXCL_LINE
1493 : str12ew, str12we, str12ns, str12sn , & ! LCOV_EXCL_LINE
1494 1382400 : strp_tmp, strm_tmp
1495 :
1496 : character(len=*), parameter :: subname = '(stress)'
1497 :
1498 : !-----------------------------------------------------------------
1499 : ! Initialize
1500 : !-----------------------------------------------------------------
1501 :
1502 31083240960 : str(:,:,:) = c0
1503 :
1504 1718695440 : do ij = 1, icellT
1505 1713188880 : i = indxTi(ij)
1506 1713188880 : j = indxTj(ij)
1507 :
1508 : !-----------------------------------------------------------------
1509 : ! strain rates
1510 : ! NOTE these are actually strain rates * area (m^2/s)
1511 : !-----------------------------------------------------------------
1512 :
1513 : call strain_rates (nx_block, ny_block, &
1514 : i, j, & ! LCOV_EXCL_LINE
1515 : uvel, vvel, & ! LCOV_EXCL_LINE
1516 : dxT, dyT, & ! LCOV_EXCL_LINE
1517 : cxp, cyp, & ! LCOV_EXCL_LINE
1518 : cxm, cym, & ! LCOV_EXCL_LINE
1519 : divune, divunw, & ! LCOV_EXCL_LINE
1520 : divuse, divusw, & ! LCOV_EXCL_LINE
1521 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1522 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1523 : shearne, shearnw, & ! LCOV_EXCL_LINE
1524 : shearse, shearsw, & ! LCOV_EXCL_LINE
1525 : Deltane, Deltanw, & ! LCOV_EXCL_LINE
1526 1713188880 : Deltase, Deltasw )
1527 :
1528 : !-----------------------------------------------------------------
1529 : ! viscosities and replacement pressure
1530 : !-----------------------------------------------------------------
1531 :
1532 149956080 : call visc_replpress (strength(i,j), DminTarea(i,j), Deltane, &
1533 1713188880 : zetax2ne, etax2ne, rep_prsne, capping)
1534 :
1535 149956080 : call visc_replpress (strength(i,j), DminTarea(i,j), Deltanw, &
1536 1713188880 : zetax2nw, etax2nw, rep_prsnw, capping)
1537 :
1538 149956080 : call visc_replpress (strength(i,j), DminTarea(i,j), Deltasw, &
1539 1713188880 : zetax2sw, etax2sw, rep_prssw, capping)
1540 :
1541 149956080 : call visc_replpress (strength(i,j), DminTarea(i,j), Deltase, &
1542 1713188880 : zetax2se, etax2se, rep_prsse, capping)
1543 :
1544 : !-----------------------------------------------------------------
1545 : ! the stresses ! kg/s^2
1546 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
1547 : !-----------------------------------------------------------------
1548 :
1549 : ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code
1550 :
1551 149956080 : stressp_1 (i,j) = (stressp_1 (i,j)*(c1-arlx1i*revp) &
1552 1713188880 : + arlx1i*(zetax2ne*divune - rep_prsne)) * denom1
1553 149956080 : stressp_2 (i,j) = (stressp_2 (i,j)*(c1-arlx1i*revp) &
1554 1713188880 : + arlx1i*(zetax2nw*divunw - rep_prsnw)) * denom1
1555 149956080 : stressp_3 (i,j) = (stressp_3 (i,j)*(c1-arlx1i*revp) &
1556 1713188880 : + arlx1i*(zetax2sw*divusw - rep_prssw)) * denom1
1557 149956080 : stressp_4 (i,j) = (stressp_4 (i,j)*(c1-arlx1i*revp) &
1558 1713188880 : + arlx1i*(zetax2se*divuse - rep_prsse)) * denom1
1559 :
1560 149956080 : stressm_1 (i,j) = (stressm_1 (i,j)*(c1-arlx1i*revp) &
1561 1713188880 : + arlx1i*etax2ne*tensionne) * denom1
1562 149956080 : stressm_2 (i,j) = (stressm_2 (i,j)*(c1-arlx1i*revp) &
1563 1713188880 : + arlx1i*etax2nw*tensionnw) * denom1
1564 149956080 : stressm_3 (i,j) = (stressm_3 (i,j)*(c1-arlx1i*revp) &
1565 1713188880 : + arlx1i*etax2sw*tensionsw) * denom1
1566 149956080 : stressm_4 (i,j) = (stressm_4 (i,j)*(c1-arlx1i*revp) &
1567 1713188880 : + arlx1i*etax2se*tensionse) * denom1
1568 :
1569 149956080 : stress12_1(i,j) = (stress12_1(i,j)*(c1-arlx1i*revp) &
1570 1713188880 : + arlx1i*p5*etax2ne*shearne) * denom1
1571 149956080 : stress12_2(i,j) = (stress12_2(i,j)*(c1-arlx1i*revp) &
1572 1713188880 : + arlx1i*p5*etax2nw*shearnw) * denom1
1573 149956080 : stress12_3(i,j) = (stress12_3(i,j)*(c1-arlx1i*revp) &
1574 1713188880 : + arlx1i*p5*etax2sw*shearsw) * denom1
1575 149956080 : stress12_4(i,j) = (stress12_4(i,j)*(c1-arlx1i*revp) &
1576 1713188880 : + arlx1i*p5*etax2se*shearse) * denom1
1577 :
1578 : !-----------------------------------------------------------------
1579 : ! Eliminate underflows.
1580 : ! The following code is commented out because it is relatively
1581 : ! expensive and most compilers include a flag that accomplishes
1582 : ! the same thing more efficiently. This code is cheaper than
1583 : ! handling underflows if the compiler lacks a flag; uncomment
1584 : ! it in that case. The compiler flag is often described with the
1585 : ! phrase "flush to zero".
1586 : !-----------------------------------------------------------------
1587 :
1588 : ! call icepack_query_parameters(puny_out=puny)
1589 : ! call icepack_warnings_flush(nu_diag)
1590 : ! if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1591 : ! file=__FILE__, line=__LINE__)
1592 :
1593 : ! stressp_1(i,j) = sign(max(abs(stressp_1(i,j)),puny),stressp_1(i,j))
1594 : ! stressp_2(i,j) = sign(max(abs(stressp_2(i,j)),puny),stressp_2(i,j))
1595 : ! stressp_3(i,j) = sign(max(abs(stressp_3(i,j)),puny),stressp_3(i,j))
1596 : ! stressp_4(i,j) = sign(max(abs(stressp_4(i,j)),puny),stressp_4(i,j))
1597 :
1598 : ! stressm_1(i,j) = sign(max(abs(stressm_1(i,j)),puny),stressm_1(i,j))
1599 : ! stressm_2(i,j) = sign(max(abs(stressm_2(i,j)),puny),stressm_2(i,j))
1600 : ! stressm_3(i,j) = sign(max(abs(stressm_3(i,j)),puny),stressm_3(i,j))
1601 : ! stressm_4(i,j) = sign(max(abs(stressm_4(i,j)),puny),stressm_4(i,j))
1602 :
1603 : ! stress12_1(i,j) = sign(max(abs(stress12_1(i,j)),puny),stress12_1(i,j))
1604 : ! stress12_2(i,j) = sign(max(abs(stress12_2(i,j)),puny),stress12_2(i,j))
1605 : ! stress12_3(i,j) = sign(max(abs(stress12_3(i,j)),puny),stress12_3(i,j))
1606 : ! stress12_4(i,j) = sign(max(abs(stress12_4(i,j)),puny),stress12_4(i,j))
1607 :
1608 : !-----------------------------------------------------------------
1609 : ! combinations of the stresses for the momentum equation ! kg/s^2
1610 : !-----------------------------------------------------------------
1611 :
1612 1713188880 : ssigpn = stressp_1(i,j) + stressp_2(i,j)
1613 1713188880 : ssigps = stressp_3(i,j) + stressp_4(i,j)
1614 1713188880 : ssigpe = stressp_1(i,j) + stressp_4(i,j)
1615 1713188880 : ssigpw = stressp_2(i,j) + stressp_3(i,j)
1616 1713188880 : ssigp1 =(stressp_1(i,j) + stressp_3(i,j))*p055
1617 1713188880 : ssigp2 =(stressp_2(i,j) + stressp_4(i,j))*p055
1618 :
1619 1713188880 : ssigmn = stressm_1(i,j) + stressm_2(i,j)
1620 1713188880 : ssigms = stressm_3(i,j) + stressm_4(i,j)
1621 1713188880 : ssigme = stressm_1(i,j) + stressm_4(i,j)
1622 1713188880 : ssigmw = stressm_2(i,j) + stressm_3(i,j)
1623 1713188880 : ssigm1 =(stressm_1(i,j) + stressm_3(i,j))*p055
1624 1713188880 : ssigm2 =(stressm_2(i,j) + stressm_4(i,j))*p055
1625 :
1626 1713188880 : ssig12n = stress12_1(i,j) + stress12_2(i,j)
1627 1713188880 : ssig12s = stress12_3(i,j) + stress12_4(i,j)
1628 1713188880 : ssig12e = stress12_1(i,j) + stress12_4(i,j)
1629 1713188880 : ssig12w = stress12_2(i,j) + stress12_3(i,j)
1630 1713188880 : ssig121 =(stress12_1(i,j) + stress12_3(i,j))*p111
1631 1713188880 : ssig122 =(stress12_2(i,j) + stress12_4(i,j))*p111
1632 :
1633 1713188880 : csigpne = p111*stressp_1(i,j) + ssigp2 + p027*stressp_3(i,j)
1634 1713188880 : csigpnw = p111*stressp_2(i,j) + ssigp1 + p027*stressp_4(i,j)
1635 1713188880 : csigpsw = p111*stressp_3(i,j) + ssigp2 + p027*stressp_1(i,j)
1636 1713188880 : csigpse = p111*stressp_4(i,j) + ssigp1 + p027*stressp_2(i,j)
1637 :
1638 1713188880 : csigmne = p111*stressm_1(i,j) + ssigm2 + p027*stressm_3(i,j)
1639 1713188880 : csigmnw = p111*stressm_2(i,j) + ssigm1 + p027*stressm_4(i,j)
1640 1713188880 : csigmsw = p111*stressm_3(i,j) + ssigm2 + p027*stressm_1(i,j)
1641 1713188880 : csigmse = p111*stressm_4(i,j) + ssigm1 + p027*stressm_2(i,j)
1642 :
1643 149956080 : csig12ne = p222*stress12_1(i,j) + ssig122 &
1644 1713188880 : + p055*stress12_3(i,j)
1645 149956080 : csig12nw = p222*stress12_2(i,j) + ssig121 &
1646 1713188880 : + p055*stress12_4(i,j)
1647 149956080 : csig12sw = p222*stress12_3(i,j) + ssig122 &
1648 1713188880 : + p055*stress12_1(i,j)
1649 149956080 : csig12se = p222*stress12_4(i,j) + ssig121 &
1650 1713188880 : + p055*stress12_2(i,j)
1651 :
1652 1713188880 : str12ew = p5*dxT(i,j)*(p333*ssig12e + p166*ssig12w)
1653 1713188880 : str12we = p5*dxT(i,j)*(p333*ssig12w + p166*ssig12e)
1654 1713188880 : str12ns = p5*dyT(i,j)*(p333*ssig12n + p166*ssig12s)
1655 1713188880 : str12sn = p5*dyT(i,j)*(p333*ssig12s + p166*ssig12n)
1656 :
1657 : !-----------------------------------------------------------------
1658 : ! for dF/dx (u momentum)
1659 : !-----------------------------------------------------------------
1660 1713188880 : strp_tmp = p25*dyT(i,j)*(p333*ssigpn + p166*ssigps)
1661 1713188880 : strm_tmp = p25*dyT(i,j)*(p333*ssigmn + p166*ssigms)
1662 :
1663 : ! northeast (i,j)
1664 149956080 : str(i,j,1) = -strp_tmp - strm_tmp - str12ew &
1665 1713188880 : + dxhy(i,j)*(-csigpne + csigmne) + dyhx(i,j)*csig12ne
1666 :
1667 : ! northwest (i+1,j)
1668 149956080 : str(i,j,2) = strp_tmp + strm_tmp - str12we &
1669 1713188880 : + dxhy(i,j)*(-csigpnw + csigmnw) + dyhx(i,j)*csig12nw
1670 :
1671 1713188880 : strp_tmp = p25*dyT(i,j)*(p333*ssigps + p166*ssigpn)
1672 1713188880 : strm_tmp = p25*dyT(i,j)*(p333*ssigms + p166*ssigmn)
1673 :
1674 : ! southeast (i,j+1)
1675 149956080 : str(i,j,3) = -strp_tmp - strm_tmp + str12ew &
1676 1713188880 : + dxhy(i,j)*(-csigpse + csigmse) + dyhx(i,j)*csig12se
1677 :
1678 : ! southwest (i+1,j+1)
1679 149956080 : str(i,j,4) = strp_tmp + strm_tmp + str12we &
1680 1713188880 : + dxhy(i,j)*(-csigpsw + csigmsw) + dyhx(i,j)*csig12sw
1681 :
1682 : !-----------------------------------------------------------------
1683 : ! for dF/dy (v momentum)
1684 : !-----------------------------------------------------------------
1685 1713188880 : strp_tmp = p25*dxT(i,j)*(p333*ssigpe + p166*ssigpw)
1686 1713188880 : strm_tmp = p25*dxT(i,j)*(p333*ssigme + p166*ssigmw)
1687 :
1688 : ! northeast (i,j)
1689 149956080 : str(i,j,5) = -strp_tmp + strm_tmp - str12ns &
1690 1713188880 : - dyhx(i,j)*(csigpne + csigmne) + dxhy(i,j)*csig12ne
1691 :
1692 : ! southeast (i,j+1)
1693 149956080 : str(i,j,6) = strp_tmp - strm_tmp - str12sn &
1694 1713188880 : - dyhx(i,j)*(csigpse + csigmse) + dxhy(i,j)*csig12se
1695 :
1696 1713188880 : strp_tmp = p25*dxT(i,j)*(p333*ssigpw + p166*ssigpe)
1697 1713188880 : strm_tmp = p25*dxT(i,j)*(p333*ssigmw + p166*ssigme)
1698 :
1699 : ! northwest (i+1,j)
1700 149956080 : str(i,j,7) = -strp_tmp + strm_tmp + str12ns &
1701 1713188880 : - dyhx(i,j)*(csigpnw + csigmnw) + dxhy(i,j)*csig12nw
1702 :
1703 : ! southwest (i+1,j+1)
1704 149956080 : str(i,j,8) = strp_tmp - strm_tmp + str12sn &
1705 2468475840 : - dyhx(i,j)*(csigpsw + csigmsw) + dxhy(i,j)*csig12sw
1706 :
1707 : enddo ! ij
1708 :
1709 5506560 : end subroutine stress
1710 :
1711 : !=======================================================================
1712 : ! Computes the strain rates and internal stress components for C grid
1713 : !
1714 : ! author: JF Lemieux, ECCC
1715 : ! updated: D. Bailey, NCAR
1716 : ! Nov 2021
1717 : !
1718 : ! Bouillon, S., T. Fichefet, V. Legat and G. Madec (2013). The
1719 : ! elastic-viscous-plastic method revisited. Ocean Model., 71, 2-12.
1720 : !
1721 : ! Kimmritz, M., S. Danilov and M. Losch (2016). The adaptive EVP method
1722 : ! for solving the sea ice momentum equation. Ocean Model., 101, 59-67.
1723 :
1724 0 : subroutine stressC_T (nx_block, ny_block , &
1725 : icellT , & ! LCOV_EXCL_LINE
1726 : indxTi , indxTj , & ! LCOV_EXCL_LINE
1727 : uvelE , vvelE , & ! LCOV_EXCL_LINE
1728 : uvelN , vvelN , & ! LCOV_EXCL_LINE
1729 : dxN , dyE , & ! LCOV_EXCL_LINE
1730 : dxT , dyT , & ! LCOV_EXCL_LINE
1731 : uarea , DminTarea , & ! LCOV_EXCL_LINE
1732 : strength , shearU , & ! LCOV_EXCL_LINE
1733 : zetax2T , etax2T , & ! LCOV_EXCL_LINE
1734 : stresspT , stressmT , & ! LCOV_EXCL_LINE
1735 0 : stress12T)
1736 :
1737 : use ice_dyn_shared, only: strain_rates_T, capping, &
1738 : visc_replpress, e_factor
1739 :
1740 : integer (kind=int_kind), intent(in) :: &
1741 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1742 : icellT ! no. of cells where iceTmask = .true.
1743 :
1744 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1745 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1746 : indxTj ! compressed index in j-direction
1747 :
1748 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1749 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1750 : vvelE , & ! y-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1751 : uvelN , & ! x-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1752 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1753 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
1754 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
1755 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1756 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1757 : strength , & ! ice strength (N/m) ! LCOV_EXCL_LINE
1758 : shearU , & ! shearU local for this routine ! LCOV_EXCL_LINE
1759 : uarea , & ! area of u cell ! LCOV_EXCL_LINE
1760 : DminTarea ! deltaminEVP*tarea
1761 :
1762 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1763 : zetax2T , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
1764 : etax2T , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
1765 : stresspT , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1766 : stressmT , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1767 : stress12T ! sigma12
1768 :
1769 : ! local variables
1770 :
1771 : integer (kind=int_kind) :: &
1772 : i, j, ij
1773 :
1774 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1775 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
1776 0 : tensionT ! tension at T point
1777 :
1778 : real (kind=dbl_kind) :: &
1779 : shearTsqr , & ! strain rates squared at T point ! LCOV_EXCL_LINE
1780 : shearT , & ! strain rate at T point ! LCOV_EXCL_LINE
1781 : DeltaT , & ! delt at T point ! LCOV_EXCL_LINE
1782 : uareaavgr , & ! 1 / uarea avg ! LCOV_EXCL_LINE
1783 0 : rep_prsT ! replacement pressure at T point
1784 :
1785 : character(len=*), parameter :: subname = '(stressC_T)'
1786 :
1787 : !-----------------------------------------------------------------
1788 : ! Initialize
1789 : !-----------------------------------------------------------------
1790 :
1791 : call strain_rates_T (nx_block , ny_block , &
1792 : icellT , & ! LCOV_EXCL_LINE
1793 : indxTi(:) , indxTj (:) , & ! LCOV_EXCL_LINE
1794 : uvelE (:,:), vvelE (:,:), & ! LCOV_EXCL_LINE
1795 : uvelN (:,:), vvelN (:,:), & ! LCOV_EXCL_LINE
1796 : dxN (:,:), dyE (:,:), & ! LCOV_EXCL_LINE
1797 : dxT (:,:), dyT (:,:), & ! LCOV_EXCL_LINE
1798 0 : divT (:,:), tensionT(:,:))
1799 :
1800 0 : do ij = 1, icellT
1801 0 : i = indxTi(ij)
1802 0 : j = indxTj(ij)
1803 :
1804 : !-----------------------------------------------------------------
1805 : ! Square of shear strain rate at T obtained from interpolation of
1806 : ! U point values (Bouillon et al., 2013, Kimmritz et al., 2016
1807 : !-----------------------------------------------------------------
1808 :
1809 0 : uareaavgr = c1/(uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))
1810 :
1811 0 : shearTsqr = (shearU(i ,j )**2 * uarea(i ,j ) &
1812 : + shearU(i ,j-1)**2 * uarea(i ,j-1) & ! LCOV_EXCL_LINE
1813 : + shearU(i-1,j-1)**2 * uarea(i-1,j-1) & ! LCOV_EXCL_LINE
1814 : + shearU(i-1,j )**2 * uarea(i-1,j )) & ! LCOV_EXCL_LINE
1815 0 : * uareaavgr
1816 :
1817 0 : shearT = (shearU(i ,j ) * uarea(i ,j ) &
1818 : + shearU(i ,j-1) * uarea(i ,j-1) & ! LCOV_EXCL_LINE
1819 : + shearU(i-1,j-1) * uarea(i-1,j-1) & ! LCOV_EXCL_LINE
1820 : + shearU(i-1,j ) * uarea(i-1,j )) & ! LCOV_EXCL_LINE
1821 0 : * uareaavgr
1822 :
1823 0 : DeltaT = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr))
1824 :
1825 : !-----------------------------------------------------------------
1826 : ! viscosities and replacement pressure at T point
1827 : !-----------------------------------------------------------------
1828 :
1829 0 : call visc_replpress (strength(i,j), DminTarea(i,j), DeltaT, &
1830 0 : zetax2T (i,j), etax2T(i,j), rep_prsT, capping)
1831 :
1832 : !-----------------------------------------------------------------
1833 : ! the stresses ! kg/s^2
1834 : !-----------------------------------------------------------------
1835 :
1836 : ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code
1837 :
1838 0 : stresspT(i,j) = (stresspT (i,j)*(c1-arlx1i*revp) &
1839 0 : + arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1
1840 :
1841 0 : stressmT(i,j) = (stressmT (i,j)*(c1-arlx1i*revp) &
1842 0 : + arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1
1843 :
1844 0 : stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) &
1845 0 : + arlx1i*p5*etax2T(i,j)*shearT ) * denom1
1846 :
1847 : enddo ! ij
1848 :
1849 0 : end subroutine stressC_T
1850 :
1851 : !=======================================================================
1852 : !
1853 : ! Computes the strain rates and internal stress components for U points
1854 : !
1855 : ! author: JF Lemieux, ECCC
1856 : ! Nov 2021
1857 : !
1858 : ! Bouillon, S., T. Fichefet, V. Legat and G. Madec (2013). The
1859 : ! elastic-viscous-plastic method revisited. Ocean Model., 71, 2-12.
1860 : !
1861 : ! Kimmritz, M., S. Danilov and M. Losch (2016). The adaptive EVP method
1862 : ! for solving the sea ice momentum equation. Ocean Model., 101, 59-67.
1863 :
1864 0 : subroutine stressC_U (nx_block , ny_block ,&
1865 : icellU ,& ! LCOV_EXCL_LINE
1866 : indxUi , indxUj ,& ! LCOV_EXCL_LINE
1867 : uarea , & ! LCOV_EXCL_LINE
1868 : etax2U , deltaU ,& ! LCOV_EXCL_LINE
1869 : strengthU, shearU ,& ! LCOV_EXCL_LINE
1870 0 : stress12U)
1871 :
1872 : use ice_dyn_shared, only: visc_replpress, &
1873 : visc_method, deltaminEVP, capping
1874 :
1875 : integer (kind=int_kind), intent(in) :: &
1876 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1877 : icellU ! no. of cells where iceUmask = 1
1878 :
1879 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1880 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1881 : indxUj ! compressed index in j-direction
1882 :
1883 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1884 : uarea , & ! area of U point ! LCOV_EXCL_LINE
1885 : etax2U , & ! 2*eta at the U point ! LCOV_EXCL_LINE
1886 : shearU , & ! shearU array ! LCOV_EXCL_LINE
1887 : deltaU , & ! deltaU array ! LCOV_EXCL_LINE
1888 : strengthU ! ice strength at the U point
1889 :
1890 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1891 : stress12U ! sigma12
1892 :
1893 : ! local variables
1894 :
1895 : integer (kind=int_kind) :: &
1896 : i, j, ij
1897 :
1898 : real (kind=dbl_kind) :: &
1899 : lzetax2U , & ! bulk viscosity at U point ! LCOV_EXCL_LINE
1900 : letax2U , & ! shear viscosity at U point ! LCOV_EXCL_LINE
1901 : lrep_prsU, & ! replacement pressure at U point ! LCOV_EXCL_LINE
1902 0 : DminUarea ! Dmin on U
1903 :
1904 : character(len=*), parameter :: subname = '(stressC_U)'
1905 :
1906 : !-----------------------------------------------------------------
1907 : ! viscosities and replacement pressure at U point
1908 : ! avg_zeta: Bouillon et al. 2013, C1 method of Kimmritz et al. 2016
1909 : ! avg_strength: C2 method of Kimmritz et al. 2016
1910 : ! if outside do and stress12U equation repeated in each loop for performance
1911 : !-----------------------------------------------------------------
1912 :
1913 0 : if (visc_method == 'avg_zeta') then
1914 0 : do ij = 1, icellU
1915 0 : i = indxUi(ij)
1916 0 : j = indxUj(ij)
1917 0 : stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) &
1918 0 : + arlx1i*p5*etax2U(i,j)*shearU(i,j)) * denom1
1919 : enddo
1920 :
1921 0 : elseif (visc_method == 'avg_strength') then
1922 0 : do ij = 1, icellU
1923 0 : i = indxUi(ij)
1924 0 : j = indxUj(ij)
1925 0 : DminUarea = deltaminEVP*uarea(i,j)
1926 : ! only need etax2U here, but other terms are calculated with etax2U
1927 : ! minimal extra calculations here even though it seems like there is
1928 0 : call visc_replpress (strengthU(i,j), DminUarea, deltaU(i,j), &
1929 0 : lzetax2U , letax2U , lrep_prsU , capping)
1930 0 : stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) &
1931 0 : + arlx1i*p5*letax2U*shearU(i,j)) * denom1
1932 : enddo
1933 :
1934 : endif
1935 :
1936 0 : end subroutine stressC_U
1937 :
1938 : !=======================================================================
1939 : ! Computes the strain rates and internal stress components for T points
1940 : !
1941 : ! author: JF Lemieux, ECCC
1942 : ! Nov 2021
1943 :
1944 0 : subroutine stressCD_T (nx_block, ny_block , &
1945 : icellT , & ! LCOV_EXCL_LINE
1946 : indxTi , indxTj , & ! LCOV_EXCL_LINE
1947 : uvelE , vvelE , & ! LCOV_EXCL_LINE
1948 : uvelN , vvelN , & ! LCOV_EXCL_LINE
1949 : dxN , dyE , & ! LCOV_EXCL_LINE
1950 : dxT , dyT , & ! LCOV_EXCL_LINE
1951 : DminTarea, & ! LCOV_EXCL_LINE
1952 : strength, & ! LCOV_EXCL_LINE
1953 : zetax2T , etax2T , & ! LCOV_EXCL_LINE
1954 : stresspT, stressmT , & ! LCOV_EXCL_LINE
1955 0 : stress12T)
1956 :
1957 : use ice_dyn_shared, only: strain_rates_T, capping, &
1958 : visc_replpress
1959 :
1960 : integer (kind=int_kind), intent(in) :: &
1961 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1962 : icellT ! no. of cells where iceTmask = .true.
1963 :
1964 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1965 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1966 : indxTj ! compressed index in j-direction
1967 :
1968 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1969 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1970 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1971 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1972 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1973 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
1974 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
1975 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1976 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1977 : strength , & ! ice strength (N/m) ! LCOV_EXCL_LINE
1978 : DminTarea ! deltaminEVP*tarea
1979 :
1980 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1981 : zetax2T , & ! zetax2 = 2*zeta (bulk viscosity) ! LCOV_EXCL_LINE
1982 : etax2T , & ! etax2 = 2*eta (shear viscosity) ! LCOV_EXCL_LINE
1983 : stresspT , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
1984 : stressmT , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
1985 : stress12T ! sigma12
1986 :
1987 : ! local variables
1988 :
1989 : integer (kind=int_kind) :: &
1990 : i, j, ij
1991 :
1992 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1993 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
1994 : tensionT , & ! tension at T point ! LCOV_EXCL_LINE
1995 : shearT , & ! shear at T point ! LCOV_EXCL_LINE
1996 0 : DeltaT ! delt at T point
1997 :
1998 : real (kind=dbl_kind) :: &
1999 0 : rep_prsT ! replacement pressure at T point
2000 :
2001 : character(len=*), parameter :: subname = '(stressCD_T)'
2002 :
2003 : !-----------------------------------------------------------------
2004 : ! strain rates at T point
2005 : ! NOTE these are actually strain rates * area (m^2/s)
2006 : !-----------------------------------------------------------------
2007 :
2008 : call strain_rates_T (nx_block , ny_block , &
2009 : icellT , & ! LCOV_EXCL_LINE
2010 : indxTi(:) , indxTj (:) , & ! LCOV_EXCL_LINE
2011 : uvelE (:,:), vvelE (:,:), & ! LCOV_EXCL_LINE
2012 : uvelN (:,:), vvelN (:,:), & ! LCOV_EXCL_LINE
2013 : dxN (:,:), dyE (:,:), & ! LCOV_EXCL_LINE
2014 : dxT (:,:), dyT (:,:), & ! LCOV_EXCL_LINE
2015 : divT (:,:), tensionT(:,:), & ! LCOV_EXCL_LINE
2016 0 : shearT(:,:), DeltaT (:,:))
2017 :
2018 0 : do ij = 1, icellT
2019 0 : i = indxTi(ij)
2020 0 : j = indxTj(ij)
2021 :
2022 : !-----------------------------------------------------------------
2023 : ! viscosities and replacement pressure at T point
2024 : !-----------------------------------------------------------------
2025 :
2026 0 : call visc_replpress (strength(i,j), DminTarea(i,j), DeltaT(i,j), &
2027 0 : zetax2T (i,j), etax2T(i,j), rep_prsT , capping)
2028 :
2029 : !-----------------------------------------------------------------
2030 : ! the stresses ! kg/s^2
2031 : !-----------------------------------------------------------------
2032 :
2033 : ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code
2034 :
2035 0 : stresspT(i,j) = (stresspT (i,j)*(c1-arlx1i*revp) &
2036 0 : + arlx1i*(zetax2T(i,j)*divT(i,j) - rep_prsT)) * denom1
2037 :
2038 0 : stressmT(i,j) = (stressmT (i,j)*(c1-arlx1i*revp) &
2039 0 : + arlx1i*etax2T(i,j)*tensionT(i,j)) * denom1
2040 :
2041 0 : stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) &
2042 0 : + arlx1i*p5*etax2T(i,j)*shearT(i,j)) * denom1
2043 :
2044 : enddo ! ij
2045 :
2046 0 : end subroutine stressCD_T
2047 :
2048 : !=======================================================================
2049 : ! Computes the strain rates and internal stress components for U points
2050 : !
2051 : ! author: JF Lemieux, ECCC
2052 : ! Nov 2021
2053 :
2054 0 : subroutine stressCD_U (nx_block, ny_block, &
2055 : icellU , & ! LCOV_EXCL_LINE
2056 : indxUi , indxUj , & ! LCOV_EXCL_LINE
2057 : uarea , & ! LCOV_EXCL_LINE
2058 : zetax2U , etax2U , & ! LCOV_EXCL_LINE
2059 : strengthU , & ! LCOV_EXCL_LINE
2060 : divergU , tensionU, & ! LCOV_EXCL_LINE
2061 : shearU , deltaU , & ! LCOV_EXCL_LINE
2062 : stresspU , stressmU, & ! LCOV_EXCL_LINE
2063 0 : stress12U)
2064 :
2065 : use ice_dyn_shared, only: visc_replpress, &
2066 : visc_method, deltaminEVP, capping
2067 :
2068 : integer (kind=int_kind), intent(in) :: &
2069 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2070 : icellU ! no. of cells where iceUmask = 1
2071 :
2072 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2073 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2074 : indxUj ! compressed index in j-direction
2075 :
2076 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2077 : uarea , & ! area of U-cell (m^2) ! LCOV_EXCL_LINE
2078 : zetax2U , & ! 2*zeta at U point ! LCOV_EXCL_LINE
2079 : etax2U , & ! 2*eta at U point ! LCOV_EXCL_LINE
2080 : strengthU, & ! ice strength at U point ! LCOV_EXCL_LINE
2081 : divergU , & ! div at U point ! LCOV_EXCL_LINE
2082 : tensionU , & ! tension at U point ! LCOV_EXCL_LINE
2083 : shearU , & ! shear at U point ! LCOV_EXCL_LINE
2084 : deltaU ! delt at U point
2085 :
2086 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
2087 : stresspU , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
2088 : stressmU , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
2089 : stress12U ! sigma12
2090 :
2091 : ! local variables
2092 :
2093 : integer (kind=int_kind) :: &
2094 : i, j, ij
2095 :
2096 : real (kind=dbl_kind) :: &
2097 : lzetax2U , & ! bulk viscosity at U point ! LCOV_EXCL_LINE
2098 : letax2U , & ! shear viscosity at U point ! LCOV_EXCL_LINE
2099 : lrep_prsU , & ! replacement pressure at U point ! LCOV_EXCL_LINE
2100 0 : DminUarea ! Dmin on U
2101 :
2102 : character(len=*), parameter :: subname = '(stressCD_U)'
2103 :
2104 0 : do ij = 1, icellU
2105 0 : i = indxUi(ij)
2106 0 : j = indxUj(ij)
2107 :
2108 : !-----------------------------------------------------------------
2109 : ! viscosities and replacement pressure at U point
2110 : ! avg_zeta: Bouillon et al. 2013, C1 method of Kimmritz et al. 2016
2111 : ! avg_strength: C2 method of Kimmritz et al. 2016
2112 : !-----------------------------------------------------------------
2113 :
2114 0 : if (visc_method == 'avg_zeta') then
2115 0 : lzetax2U = zetax2U(i,j)
2116 0 : letax2U = etax2U(i,j)
2117 0 : lrep_prsU = (c1-Ktens)/(c1+Ktens)*lzetax2U*deltaU(i,j)
2118 :
2119 0 : elseif (visc_method == 'avg_strength') then
2120 0 : DminUarea = deltaminEVP*uarea(i,j)
2121 : ! only need etax2U here, but other terms are calculated with etax2U
2122 : ! minimal extra calculations here even though it seems like there is
2123 0 : call visc_replpress (strengthU(i,j), DminUarea, deltaU(i,j), &
2124 0 : lzetax2U , letax2U , lrep_prsU , capping)
2125 : endif
2126 :
2127 : !-----------------------------------------------------------------
2128 : ! the stresses ! kg/s^2
2129 : !-----------------------------------------------------------------
2130 :
2131 : ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code
2132 :
2133 0 : stresspU(i,j) = (stresspU (i,j)*(c1-arlx1i*revp) &
2134 0 : + arlx1i*(lzetax2U*divergU(i,j) - lrep_prsU)) * denom1
2135 :
2136 0 : stressmU(i,j) = (stressmU (i,j)*(c1-arlx1i*revp) &
2137 0 : + arlx1i*letax2U*tensionU(i,j)) * denom1
2138 :
2139 0 : stress12U(i,j) = (stress12U(i,j)*(c1-arlx1i*revp) &
2140 0 : + arlx1i*p5*letax2U*shearU(i,j)) * denom1
2141 :
2142 : enddo ! ij
2143 :
2144 0 : end subroutine stressCD_U
2145 :
2146 : !=======================================================================
2147 : ! Computes divergence of stress tensor at the E or N point for the mom equation
2148 : !
2149 : ! author: JF Lemieux, ECCC
2150 : ! Nov 2021
2151 : !
2152 : ! Hunke, E. C., and J. K. Dukowicz (2002). The Elastic-Viscous-Plastic
2153 : ! Sea Ice Dynamics Model in General Orthogonal Curvilinear Coordinates
2154 : ! on a Sphere - Incorporation of Metric Terms. Mon. Weather Rev.,
2155 : ! 130, 1848-1865.
2156 : !
2157 : ! Bouillon, S., M. Morales Maqueda, V. Legat and T. Fichefet (2009). An
2158 : ! elastic-viscous-plastic sea ice model formulated on Arakawa B and C grids.
2159 : ! Ocean Model., 27, 174-184.
2160 :
2161 0 : subroutine div_stress_Ex(nx_block, ny_block, &
2162 : icell , & ! LCOV_EXCL_LINE
2163 : indxi , indxj , & ! LCOV_EXCL_LINE
2164 : dxE , dyE , & ! LCOV_EXCL_LINE
2165 : dxU , dyT , & ! LCOV_EXCL_LINE
2166 : arear , & ! LCOV_EXCL_LINE
2167 : stressp , stressm , & ! LCOV_EXCL_LINE
2168 : stress12, & ! LCOV_EXCL_LINE
2169 0 : strintx )
2170 :
2171 :
2172 : integer (kind=int_kind), intent(in) :: &
2173 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2174 : icell ! no. of cells where epm (or npm) = 1
2175 :
2176 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2177 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2178 : indxj ! compressed index in j-direction
2179 :
2180 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2181 : dxE , & ! width of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2182 : dyE , & ! height of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2183 : dxU , & ! width of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2184 : dyT , & ! height of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2185 : arear ! earear or narear
2186 :
2187 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(in) :: &
2188 : stressp , & ! stressp (U or T) used for strintx calculation ! LCOV_EXCL_LINE
2189 : stressm , & ! stressm (U or T) used for strintx calculation ! LCOV_EXCL_LINE
2190 : stress12 ! stress12 (U or T) used for strintx calculation
2191 :
2192 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(out) :: &
2193 : strintx ! div of stress tensor for u component
2194 :
2195 : ! local variables
2196 :
2197 : integer (kind=int_kind) :: &
2198 : i, j, ij
2199 :
2200 : character(len=*), parameter :: subname = '(div_stress_Ex)'
2201 :
2202 0 : do ij = 1, icell
2203 0 : i = indxi(ij)
2204 0 : j = indxj(ij)
2205 0 : strintx(i,j) = arear(i,j) * &
2206 : ( p5 * dyE(i,j) * ( stressp(i+1,j ) - stressp (i ,j ) ) & ! LCOV_EXCL_LINE
2207 : + (p5/ dyE(i,j)) * ( (dyT(i+1,j )**2) * stressm (i+1,j ) & ! LCOV_EXCL_LINE
2208 : -(dyT(i ,j )**2) * stressm (i ,j ) ) & ! LCOV_EXCL_LINE
2209 : + (c1/ dxE(i,j)) * ( (dxU(i ,j )**2) * stress12(i ,j ) & ! LCOV_EXCL_LINE
2210 0 : -(dxU(i ,j-1)**2) * stress12(i ,j-1) ) )
2211 : enddo
2212 :
2213 0 : end subroutine div_stress_Ex
2214 :
2215 : !=======================================================================
2216 0 : subroutine div_stress_Ey(nx_block, ny_block, &
2217 : icell , & ! LCOV_EXCL_LINE
2218 : indxi , indxj , & ! LCOV_EXCL_LINE
2219 : dxE , dyE , & ! LCOV_EXCL_LINE
2220 : dxU , dyT , & ! LCOV_EXCL_LINE
2221 : arear , & ! LCOV_EXCL_LINE
2222 : stressp , stressm , & ! LCOV_EXCL_LINE
2223 : stress12, & ! LCOV_EXCL_LINE
2224 0 : strinty )
2225 :
2226 : integer (kind=int_kind), intent(in) :: &
2227 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2228 : icell ! no. of cells where epm (or npm) = 1
2229 :
2230 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2231 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2232 : indxj ! compressed index in j-direction
2233 :
2234 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2235 : dxE , & ! width of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2236 : dyE , & ! height of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2237 : dxU , & ! width of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2238 : dyT , & ! height of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2239 : arear ! earear or narear
2240 :
2241 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(in) :: &
2242 : stressp , & ! stressp (U or T) used for strinty calculation ! LCOV_EXCL_LINE
2243 : stressm , & ! stressm (U or T) used for strinty calculation ! LCOV_EXCL_LINE
2244 : stress12 ! stress12 (U or T) used for strinty calculation
2245 :
2246 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(out) :: &
2247 : strinty ! div of stress tensor for v component
2248 :
2249 : ! local variables
2250 :
2251 : integer (kind=int_kind) :: &
2252 : i, j, ij
2253 :
2254 : character(len=*), parameter :: subname = '(div_stress_Ey)'
2255 :
2256 0 : do ij = 1, icell
2257 0 : i = indxi(ij)
2258 0 : j = indxj(ij)
2259 0 : strinty(i,j) = arear(i,j) * &
2260 : ( p5 * dxE(i,j) * ( stressp(i ,j ) - stressp (i ,j-1) ) & ! LCOV_EXCL_LINE
2261 : - (p5/ dxE(i,j)) * ( (dxU(i ,j )**2) * stressm (i ,j ) & ! LCOV_EXCL_LINE
2262 : -(dxU(i ,j-1)**2) * stressm (i ,j-1) ) & ! LCOV_EXCL_LINE
2263 : + (c1/ dyE(i,j)) * ( (dyT(i+1,j )**2) * stress12(i+1,j ) & ! LCOV_EXCL_LINE
2264 0 : -(dyT(i ,j )**2) * stress12(i ,j ) ) )
2265 : enddo
2266 :
2267 0 : end subroutine div_stress_Ey
2268 :
2269 : !=======================================================================
2270 0 : subroutine div_stress_Nx(nx_block, ny_block, &
2271 : icell , & ! LCOV_EXCL_LINE
2272 : indxi , indxj , & ! LCOV_EXCL_LINE
2273 : dxN , dyN , & ! LCOV_EXCL_LINE
2274 : dxT , dyU , & ! LCOV_EXCL_LINE
2275 : arear , & ! LCOV_EXCL_LINE
2276 : stressp , stressm , & ! LCOV_EXCL_LINE
2277 : stress12, & ! LCOV_EXCL_LINE
2278 0 : strintx )
2279 :
2280 : integer (kind=int_kind), intent(in) :: &
2281 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2282 : icell ! no. of cells where epm (or npm) = 1
2283 :
2284 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2285 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2286 : indxj ! compressed index in j-direction
2287 :
2288 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2289 : dxN , & ! width of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2290 : dyN , & ! height of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2291 : dxT , & ! width of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2292 : dyU , & ! height of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2293 : arear ! earear or narear
2294 :
2295 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(in) :: &
2296 : stressp , & ! stressp (U or T) used for strintx calculation ! LCOV_EXCL_LINE
2297 : stressm , & ! stressm (U or T) used for strintx calculation ! LCOV_EXCL_LINE
2298 : stress12 ! stress12 (U or T) used for strintx calculation
2299 :
2300 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(out) :: &
2301 : strintx ! div of stress tensor for u component
2302 :
2303 : ! local variables
2304 :
2305 : integer (kind=int_kind) :: &
2306 : i, j, ij
2307 :
2308 : character(len=*), parameter :: subname = '(div_stress_Nx)'
2309 :
2310 0 : do ij = 1, icell
2311 0 : i = indxi(ij)
2312 0 : j = indxj(ij)
2313 0 : strintx(i,j) = arear(i,j) * &
2314 : ( p5 * dyN(i,j) * ( stressp(i ,j ) - stressp (i-1,j ) ) & ! LCOV_EXCL_LINE
2315 : + (p5/ dyN(i,j)) * ( (dyU(i ,j )**2) * stressm (i ,j ) & ! LCOV_EXCL_LINE
2316 : -(dyU(i-1,j )**2) * stressm (i-1,j ) ) & ! LCOV_EXCL_LINE
2317 : + (c1/ dxN(i,j)) * ( (dxT(i ,j+1)**2) * stress12(i ,j+1) & ! LCOV_EXCL_LINE
2318 0 : -(dxT(i ,j )**2) * stress12(i ,j ) ) )
2319 : enddo
2320 :
2321 0 : end subroutine div_stress_Nx
2322 :
2323 : !=======================================================================
2324 0 : subroutine div_stress_Ny(nx_block, ny_block, &
2325 : icell , & ! LCOV_EXCL_LINE
2326 : indxi , indxj , & ! LCOV_EXCL_LINE
2327 : dxN , dyN , & ! LCOV_EXCL_LINE
2328 : dxT , dyU , & ! LCOV_EXCL_LINE
2329 : arear , & ! LCOV_EXCL_LINE
2330 : stressp , stressm , & ! LCOV_EXCL_LINE
2331 : stress12, & ! LCOV_EXCL_LINE
2332 0 : strinty )
2333 :
2334 : integer (kind=int_kind), intent(in) :: &
2335 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2336 : icell ! no. of cells where epm (or npm) = 1
2337 :
2338 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2339 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2340 : indxj ! compressed index in j-direction
2341 :
2342 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2343 : dxN , & ! width of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2344 : dyN , & ! height of E or N-cell through the middle (m) ! LCOV_EXCL_LINE
2345 : dxT , & ! width of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2346 : dyU , & ! height of T or U-cell through the middle (m) ! LCOV_EXCL_LINE
2347 : arear ! earear or narear
2348 :
2349 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(in) :: &
2350 : stressp , & ! stressp (U or T) used for strinty calculation ! LCOV_EXCL_LINE
2351 : stressm , & ! stressm (U or T) used for strinty calculation ! LCOV_EXCL_LINE
2352 : stress12 ! stress12 (U or T) used for strinty calculation
2353 :
2354 : real (kind=dbl_kind), optional, dimension (nx_block,ny_block), intent(out) :: &
2355 : strinty ! div of stress tensor for v component
2356 :
2357 : ! local variables
2358 :
2359 : integer (kind=int_kind) :: &
2360 : i, j, ij
2361 :
2362 : character(len=*), parameter :: subname = '(div_stress_Ny)'
2363 :
2364 0 : do ij = 1, icell
2365 0 : i = indxi(ij)
2366 0 : j = indxj(ij)
2367 0 : strinty(i,j) = arear(i,j) * &
2368 : ( p5 * dxN(i,j) * ( stressp(i ,j+1) - stressp (i ,j ) ) & ! LCOV_EXCL_LINE
2369 : - (p5/ dxN(i,j)) * ( (dxT(i ,j+1)**2) * stressm (i ,j+1) & ! LCOV_EXCL_LINE
2370 : -(dxT(i ,j )**2) * stressm (i ,j ) ) & ! LCOV_EXCL_LINE
2371 : + (c1/ dyN(i,j)) * ( (dyU(i ,j )**2) * stress12(i ,j ) & ! LCOV_EXCL_LINE
2372 0 : -(dyU(i-1,j )**2) * stress12(i-1,j ) ) )
2373 : enddo
2374 :
2375 0 : end subroutine div_stress_Ny
2376 :
2377 : !=======================================================================
2378 :
2379 : end module ice_dyn_evp
2380 :
2381 : !=======================================================================
|
