Line data Source code
1 : !=======================================================================
2 :
3 : ! Elastic-viscous-plastic sea ice dynamics model code shared with other
4 : ! approaches
5 : !
6 : ! author: Elizabeth C. Hunke, LANL
7 : !
8 : ! 2013: Split from ice_dyn_evp.F90 by Elizabeth Hunke
9 :
10 : module ice_dyn_shared
11 :
12 : use ice_kinds_mod
13 : use ice_communicate, only: my_task, master_task, get_num_procs
14 : use ice_constants, only: c0, c1, c2, c3, c4, c6
15 : use ice_constants, only: omega, spval_dbl, p01, p001, p5
16 : use ice_blocks, only: nx_block, ny_block
17 : use ice_domain_size, only: max_blocks
18 : use ice_fileunits, only: nu_diag
19 : use ice_exit, only: abort_ice
20 : use ice_grid, only: grid_ice
21 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
22 : use icepack_intfc, only: icepack_query_parameters
23 :
24 : implicit none
25 : private
26 : public :: set_evp_parameters, stepu, stepuv_CD, stepu_C, stepv_C, &
27 : principal_stress, init_dyn_shared, dyn_prep1, dyn_prep2, dyn_finish, & ! LCOV_EXCL_LINE
28 : seabed_stress_factor_LKD, seabed_stress_factor_prob, & ! LCOV_EXCL_LINE
29 : deformations, deformationsC_T, deformationsCD_T, & ! LCOV_EXCL_LINE
30 : strain_rates, strain_rates_T, strain_rates_U, & ! LCOV_EXCL_LINE
31 : visc_replpress, & ! LCOV_EXCL_LINE
32 : dyn_haloUpdate, & ! LCOV_EXCL_LINE
33 : stack_fields, unstack_fields
34 :
35 : ! namelist parameters
36 :
37 : integer (kind=int_kind), public :: &
38 : kdyn , & ! type of dynamics ( -1, 0 = off, 1 = evp, 2 = eap ) ! LCOV_EXCL_LINE
39 : kridge , & ! set to "-1" to turn off ridging ! LCOV_EXCL_LINE
40 : ndte ! number of subcycles
41 :
42 : character (len=char_len), public :: &
43 : coriolis , & ! 'constant', 'zero', or 'latitude' ! LCOV_EXCL_LINE
44 : ssh_stress ! 'geostrophic' or 'coupled'
45 :
46 : logical (kind=log_kind), public :: &
47 : revised_evp ! if true, use revised evp procedure
48 :
49 : character (len=char_len), public :: &
50 : evp_algorithm ! standard_2d = 2D org version (standard)
51 : ! shared_mem_1d = 1d without mpi call and refactorization to 1d
52 :
53 : real (kind=dbl_kind), public :: &
54 : elasticDamp ! coefficient for calculating the parameter E, elastic damping parameter
55 :
56 : ! other EVP parameters
57 :
58 : character (len=char_len), public :: &
59 : yield_curve , & ! 'ellipse' ('teardrop' needs further testing) ! LCOV_EXCL_LINE
60 : visc_method , & ! method for viscosity calc at U points (C, CD grids) ! LCOV_EXCL_LINE
61 : seabed_stress_method ! method for seabed stress calculation
62 : ! LKD: Lemieux et al. 2015, probabilistic: Dupont et al. 2022
63 :
64 : real (kind=dbl_kind), parameter, public :: &
65 : u0 = 5e-5_dbl_kind, & ! residual velocity for seabed stress (m/s) ! LCOV_EXCL_LINE
66 : cosw = c1 , & ! cos(ocean turning angle) ! turning angle = 0 ! LCOV_EXCL_LINE
67 : sinw = c0 , & ! sin(ocean turning angle) ! turning angle = 0 ! LCOV_EXCL_LINE
68 : a_min = p001 , & ! minimum ice area ! LCOV_EXCL_LINE
69 : m_min = p01 ! minimum ice mass (kg/m^2)
70 :
71 : real (kind=dbl_kind), public :: &
72 : revp , & ! 0 for classic EVP, 1 for revised EVP ! LCOV_EXCL_LINE
73 : e_yieldcurve, & ! VP aspect ratio of elliptical yield curve ! LCOV_EXCL_LINE
74 : e_plasticpot, & ! VP aspect ratio of elliptical plastic potential ! LCOV_EXCL_LINE
75 : epp2i , & ! 1/(e_plasticpot)^2 ! LCOV_EXCL_LINE
76 : e_factor , & ! (e_yieldcurve)^2/(e_plasticpot)^4 ! LCOV_EXCL_LINE
77 : ecci , & ! temporary for 1d evp ! LCOV_EXCL_LINE
78 : deltaminEVP , & ! minimum delta for viscosities (EVP) ! LCOV_EXCL_LINE
79 : deltaminVP , & ! minimum delta for viscosities (VP) ! LCOV_EXCL_LINE
80 : capping , & ! capping of viscosities (1=Hibler79, 0=Kreyscher2000) ! LCOV_EXCL_LINE
81 : dtei , & ! ndte/dt, where dt/ndte is subcycling timestep (1/s) ! LCOV_EXCL_LINE
82 : denom1 ! constants for stress equation
83 :
84 : real (kind=dbl_kind), public :: & ! Bouillon et al relaxation constants
85 : arlx , & ! alpha for stressp ! LCOV_EXCL_LINE
86 : arlx1i , & ! (inverse of alpha) for stressp ! LCOV_EXCL_LINE
87 : brlx ! beta for momentum
88 :
89 : real (kind=dbl_kind), allocatable, public :: &
90 : fcor_blk(:,:,:) ! Coriolis parameter (1/s)
91 :
92 : real (kind=dbl_kind), allocatable, public :: &
93 : fcorE_blk(:,:,:), & ! Coriolis parameter at E points (1/s) ! LCOV_EXCL_LINE
94 : fcorN_blk(:,:,:) ! Coriolis parameter at N points (1/s)
95 :
96 : real (kind=dbl_kind), allocatable, public :: &
97 : fld2(:,:,:,:), & ! 2 bundled fields ! LCOV_EXCL_LINE
98 : fld3(:,:,:,:), & ! 3 bundled fields ! LCOV_EXCL_LINE
99 : fld4(:,:,:,:) ! 4 bundled fields
100 :
101 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
102 : uvel_init , & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
103 : vvel_init ! y-component of velocity (m/s), beginning of timestep
104 :
105 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
106 : uvelN_init , & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
107 : vvelN_init ! y-component of velocity (m/s), beginning of timestep
108 :
109 : real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
110 : uvelE_init , & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
111 : vvelE_init ! y-component of velocity (m/s), beginning of timestep
112 :
113 : logical (kind=log_kind), dimension (:,:,:), allocatable, public :: &
114 : iceTmask, & ! ice extent mask (T-cell) ! LCOV_EXCL_LINE
115 : iceUmask, & ! ice extent mask (U-cell) ! LCOV_EXCL_LINE
116 : iceNmask, & ! ice extent mask (N-cell) ! LCOV_EXCL_LINE
117 : iceEmask ! ice extent mask (E-cell)
118 :
119 : real (kind=dbl_kind), allocatable, public :: &
120 : DminTarea(:,:,:) ! deltamin * tarea (m^2/s)
121 :
122 : ! ice isotropic tensile strength parameter
123 : real (kind=dbl_kind), public :: &
124 : Ktens ! T=Ktens*P (tensile strength: see Konig and Holland, 2010)
125 :
126 : ! seabed (basal) stress parameters and settings
127 : logical (kind=log_kind), public :: &
128 : seabed_stress ! if true, seabed stress for landfast on
129 :
130 : real (kind=dbl_kind), public :: &
131 : k1 , & ! 1st free parameter for seabed1 grounding parameterization ! LCOV_EXCL_LINE
132 : k2 , & ! second free parameter (N/m^3) for seabed1 grounding parametrization ! LCOV_EXCL_LINE
133 : alphab , & ! alphab=Cb factor in Lemieux et al 2015 ! LCOV_EXCL_LINE
134 : threshold_hw ! max water depth for grounding
135 : ! see keel data from Amundrud et al. 2004 (JGR)
136 :
137 : interface strain_rates_T
138 : module procedure strain_rates_Tdt
139 : module procedure strain_rates_Tdtsd
140 : end interface
141 :
142 : interface dyn_haloUpdate
143 : module procedure dyn_haloUpdate1
144 : module procedure dyn_haloUpdate2
145 : module procedure dyn_haloUpdate3
146 : module procedure dyn_haloUpdate4
147 : module procedure dyn_haloUpdate5
148 : end interface
149 :
150 : interface stack_fields
151 : module procedure stack_fields2
152 : module procedure stack_fields3
153 : module procedure stack_fields4
154 : module procedure stack_fields5
155 : end interface
156 :
157 : interface unstack_fields
158 : module procedure unstack_fields2
159 : module procedure unstack_fields3
160 : module procedure unstack_fields4
161 : module procedure unstack_fields5
162 : end interface
163 :
164 : !=======================================================================
165 :
166 : contains
167 :
168 : !=======================================================================
169 : !
170 : ! Allocate space for all variables
171 : !
172 37 : subroutine alloc_dyn_shared
173 :
174 : integer (int_kind) :: ierr
175 :
176 : character(len=*), parameter :: subname = '(alloc_dyn_shared)'
177 :
178 : allocate( &
179 : uvel_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
180 : vvel_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
181 : iceTmask (nx_block,ny_block,max_blocks), & ! T mask for dynamics ! LCOV_EXCL_LINE
182 : iceUmask (nx_block,ny_block,max_blocks), & ! U mask for dynamics ! LCOV_EXCL_LINE
183 : fcor_blk (nx_block,ny_block,max_blocks), & ! Coriolis ! LCOV_EXCL_LINE
184 : DminTarea (nx_block,ny_block,max_blocks), & ! ! LCOV_EXCL_LINE
185 37 : stat=ierr)
186 37 : if (ierr/=0) call abort_ice(subname//': Out of memory')
187 :
188 : allocate( &
189 : fld2(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
190 : fld3(nx_block,ny_block,3,max_blocks), & ! LCOV_EXCL_LINE
191 : fld4(nx_block,ny_block,4,max_blocks), & ! LCOV_EXCL_LINE
192 37 : stat=ierr)
193 37 : if (ierr/=0) call abort_ice(subname//': Out of memory')
194 :
195 37 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
196 : allocate( &
197 : uvelE_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
198 : vvelE_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
199 : uvelN_init (nx_block,ny_block,max_blocks), & ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
200 : vvelN_init (nx_block,ny_block,max_blocks), & ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
201 : iceEmask (nx_block,ny_block,max_blocks), & ! T mask for dynamics ! LCOV_EXCL_LINE
202 : iceNmask (nx_block,ny_block,max_blocks), & ! U mask for dynamics ! LCOV_EXCL_LINE
203 : fcorE_blk (nx_block,ny_block,max_blocks), & ! Coriolis ! LCOV_EXCL_LINE
204 : fcorN_blk (nx_block,ny_block,max_blocks), & ! Coriolis ! LCOV_EXCL_LINE
205 0 : stat=ierr)
206 0 : if (ierr/=0) call abort_ice(subname//': Out of memory')
207 : endif
208 :
209 37 : end subroutine alloc_dyn_shared
210 :
211 : !=======================================================================
212 : ! Initialize parameters and variables needed for the dynamics
213 : ! author: Elizabeth C. Hunke, LANL
214 :
215 37 : subroutine init_dyn_shared (dt)
216 :
217 : use ice_blocks, only: nx_block, ny_block
218 : use ice_domain, only: nblocks, halo_dynbundle
219 : use ice_domain_size, only: max_blocks
220 : use ice_flux, only: &
221 : stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
222 : stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
223 : stress12_1, stress12_2, stress12_3, stress12_4, & ! LCOV_EXCL_LINE
224 : stresspT, stressmT, stress12T, & ! LCOV_EXCL_LINE
225 : stresspU, stressmU, stress12U
226 : use ice_state, only: uvel, vvel, uvelE, vvelE, uvelN, vvelN
227 : use ice_grid, only: ULAT, NLAT, ELAT, tarea
228 :
229 : real (kind=dbl_kind), intent(in) :: &
230 : dt ! time step
231 :
232 : ! local variables
233 :
234 : integer (kind=int_kind) :: &
235 : i, j , & ! indices ! LCOV_EXCL_LINE
236 : nprocs, & ! number of processors ! LCOV_EXCL_LINE
237 : iblk ! block index
238 :
239 : character(len=*), parameter :: subname = '(init_dyn_shared)'
240 :
241 37 : call set_evp_parameters (dt)
242 : ! allocate dyn shared (init_uvel,init_vvel)
243 37 : call alloc_dyn_shared
244 : ! Set halo_dynbundle, this is empirical at this point, could become namelist
245 37 : halo_dynbundle = .true.
246 37 : nprocs = get_num_procs()
247 37 : if (nx_block*ny_block/nprocs > 100) halo_dynbundle = .false.
248 :
249 37 : if (my_task == master_task) then
250 7 : write(nu_diag,*) 'dt = ',dt
251 7 : write(nu_diag,*) 'dt_subcyle = ',dt/real(ndte,kind=dbl_kind)
252 7 : write(nu_diag,*) 'tdamp =', elasticDamp * dt
253 7 : write(nu_diag,*) 'halo_dynbundle =', halo_dynbundle
254 : endif
255 :
256 20 : !$OMP PARALLEL DO PRIVATE(iblk,i,j) SCHEDULE(runtime)
257 50 : do iblk = 1, nblocks
258 1276 : do j = 1, ny_block
259 50267 : do i = 1, nx_block
260 :
261 : ! velocity
262 49028 : uvel(i,j,iblk) = c0 ! m/s
263 49028 : vvel(i,j,iblk) = c0 ! m/s
264 49028 : if (grid_ice == 'CD' .or. grid_ice == 'C') then ! extra velocity variables
265 0 : uvelE(i,j,iblk) = c0
266 0 : vvelE(i,j,iblk) = c0
267 0 : uvelN(i,j,iblk) = c0
268 0 : vvelN(i,j,iblk) = c0
269 : endif
270 :
271 :
272 : ! Coriolis parameter
273 49028 : if (trim(coriolis) == 'constant') then
274 0 : fcor_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
275 49028 : else if (trim(coriolis) == 'zero') then
276 0 : fcor_blk(i,j,iblk) = c0
277 : else
278 49028 : fcor_blk(i,j,iblk) = c2*omega*sin(ULAT(i,j,iblk)) ! 1/s
279 : endif
280 :
281 49028 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
282 :
283 0 : if (trim(coriolis) == 'constant') then
284 0 : fcorE_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
285 0 : fcorN_blk(i,j,iblk) = 1.46e-4_dbl_kind ! Hibler 1979, N. Hem; 1/s
286 0 : else if (trim(coriolis) == 'zero') then
287 0 : fcorE_blk(i,j,iblk) = c0
288 0 : fcorN_blk(i,j,iblk) = c0
289 : else
290 0 : fcorE_blk(i,j,iblk) = c2*omega*sin(ELAT(i,j,iblk)) ! 1/s
291 0 : fcorN_blk(i,j,iblk) = c2*omega*sin(NLAT(i,j,iblk)) ! 1/s
292 : endif
293 :
294 : endif
295 :
296 : ! stress tensor, kg/s^2
297 49028 : stressp_1 (i,j,iblk) = c0
298 49028 : stressp_2 (i,j,iblk) = c0
299 49028 : stressp_3 (i,j,iblk) = c0
300 49028 : stressp_4 (i,j,iblk) = c0
301 49028 : stressm_1 (i,j,iblk) = c0
302 49028 : stressm_2 (i,j,iblk) = c0
303 49028 : stressm_3 (i,j,iblk) = c0
304 49028 : stressm_4 (i,j,iblk) = c0
305 49028 : stress12_1(i,j,iblk) = c0
306 49028 : stress12_2(i,j,iblk) = c0
307 49028 : stress12_3(i,j,iblk) = c0
308 49028 : stress12_4(i,j,iblk) = c0
309 :
310 49028 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
311 0 : stresspT (i,j,iblk) = c0
312 0 : stressmT (i,j,iblk) = c0
313 0 : stress12T (i,j,iblk) = c0
314 0 : stresspU (i,j,iblk) = c0
315 0 : stressmU (i,j,iblk) = c0
316 0 : stress12U (i,j,iblk) = c0
317 : endif
318 :
319 49028 : if (kdyn == 1) then
320 49028 : DminTarea(i,j,iblk) = deltaminEVP*tarea(i,j,iblk)
321 0 : elseif (kdyn == 3) then
322 0 : DminTarea(i,j,iblk) = deltaminVP*tarea(i,j,iblk)
323 : endif
324 :
325 : ! ice extent mask on velocity points
326 49028 : iceUmask(i,j,iblk) = .false.
327 50234 : if (grid_ice == 'CD' .or. grid_ice == 'C') then
328 0 : iceEmask(i,j,iblk) = .false.
329 0 : iceNmask(i,j,iblk) = .false.
330 : end if
331 : enddo ! i
332 : enddo ! j
333 : enddo ! iblk
334 : !$OMP END PARALLEL DO
335 :
336 37 : end subroutine init_dyn_shared
337 :
338 : !=======================================================================
339 : ! Set parameters needed for the evp dynamics.
340 : ! Note: This subroutine is currently called only during initialization.
341 : ! If the dynamics time step can vary during runtime, it should
342 : ! be called whenever the time step changes.
343 : !
344 : ! author: Elizabeth C. Hunke, LANL
345 :
346 37 : subroutine set_evp_parameters (dt)
347 :
348 : real (kind=dbl_kind), intent(in) :: &
349 : dt ! time step
350 :
351 : ! local variables
352 :
353 : character(len=*), parameter :: subname = '(set_evp_parameters)'
354 :
355 : ! elastic time step
356 37 : dtei = real(ndte,kind=dbl_kind)/dt
357 :
358 : ! variables for elliptical yield curve and plastic potential
359 37 : epp2i = c1/e_plasticpot**2
360 37 : e_factor = e_yieldcurve**2 / e_plasticpot**4
361 37 : ecci = c1/e_yieldcurve**2 ! temporary for 1d evp
362 :
363 37 : if (revised_evp) then ! Bouillon et al, Ocean Mod 2013
364 0 : revp = c1
365 0 : denom1 = c1
366 0 : arlx1i = c1/arlx
367 : else ! Hunke, JCP 2013 with modified stress eq
368 37 : revp = c0
369 37 : arlx = c2 * elasticDamp * real(ndte,kind=dbl_kind)
370 37 : arlx1i = c1/arlx
371 37 : brlx = real(ndte,kind=dbl_kind)
372 37 : denom1 = c1/(c1+arlx1i)
373 : endif
374 37 : if (my_task == master_task) then
375 7 : write (nu_diag,*) 'arlx, arlxi, brlx, denom1', &
376 14 : arlx, arlx1i, brlx, denom1
377 : endif
378 :
379 37 : end subroutine set_evp_parameters
380 :
381 : !=======================================================================
382 : ! Computes quantities needed in the stress tensor (sigma)
383 : ! and momentum (u) equations, but which do not change during
384 : ! the thermodynamics/transport time step:
385 : ! ice mass and ice extent masks
386 : !
387 : ! author: Elizabeth C. Hunke, LANL
388 :
389 22944 : subroutine dyn_prep1 (nx_block, ny_block, &
390 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
391 : aice, vice, & ! LCOV_EXCL_LINE
392 : vsno, Tmask, & ! LCOV_EXCL_LINE
393 22944 : Tmass, iceTmask)
394 :
395 : integer (kind=int_kind), intent(in) :: &
396 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
397 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
398 :
399 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
400 : aice , & ! concentration of ice ! LCOV_EXCL_LINE
401 : vice , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
402 : vsno ! volume per unit area of snow (m)
403 :
404 : logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: &
405 : Tmask ! land/boundary mask, thickness (T-cell)
406 :
407 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
408 : Tmass ! total mass of ice and snow (kg/m^2)
409 :
410 : logical (kind=log_kind), dimension (nx_block,ny_block), intent(out) :: &
411 : iceTmask ! ice extent mask (T-cell)
412 :
413 : ! local variables
414 :
415 : integer (kind=int_kind) :: &
416 : i, j
417 :
418 : real (kind=dbl_kind) :: &
419 5760 : rhoi, rhos
420 :
421 : logical (kind=log_kind), dimension(nx_block,ny_block) :: &
422 45888 : tmphm ! temporary mask
423 :
424 : character(len=*), parameter :: subname = '(dyn_prep1)'
425 :
426 22944 : call icepack_query_parameters(rhos_out=rhos, rhoi_out=rhoi)
427 22944 : call icepack_warnings_flush(nu_diag)
428 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
429 0 : file=__FILE__, line=__LINE__)
430 :
431 753576 : do j = 1, ny_block
432 16186320 : do i = 1, nx_block
433 :
434 : !-----------------------------------------------------------------
435 : ! total mass of ice and snow, centered in T-cell
436 : ! NOTE: vice and vsno must be up to date in all grid cells,
437 : ! including ghost cells
438 : !-----------------------------------------------------------------
439 15432744 : if (Tmask(i,j)) then
440 12324840 : Tmass(i,j) = (rhoi*vice(i,j) + rhos*vsno(i,j)) ! kg/m^2
441 : else
442 3107904 : Tmass(i,j) = c0
443 : endif
444 :
445 : !-----------------------------------------------------------------
446 : ! ice extent mask (T-cells)
447 : !-----------------------------------------------------------------
448 4821120 : tmphm(i,j) = Tmask(i,j) .and. (aice (i,j) > a_min) &
449 15432744 : .and. (Tmass(i,j) > m_min)
450 :
451 : !-----------------------------------------------------------------
452 : ! augmented mask (land + open ocean)
453 : !-----------------------------------------------------------------
454 16163376 : iceTmask (i,j) = .false.
455 :
456 : enddo
457 : enddo
458 :
459 707688 : do j = jlo, jhi
460 13826928 : do i = ilo, ihi
461 :
462 : ! extend ice extent mask (T-cells) to points around pack
463 25056000 : if (tmphm(i-1,j+1) .or. tmphm(i,j+1) .or. tmphm(i+1,j+1) .or. &
464 : tmphm(i-1,j) .or. tmphm(i,j) .or. tmphm(i+1,j) .or. & ! LCOV_EXCL_LINE
465 50703240 : tmphm(i-1,j-1) .or. tmphm(i,j-1) .or. tmphm(i+1,j-1) ) then
466 7176702 : iceTmask(i,j) = .true.
467 : endif
468 :
469 13803984 : if (.not.Tmask(i,j)) iceTmask(i,j) = .false.
470 :
471 : enddo
472 : enddo
473 :
474 22944 : end subroutine dyn_prep1
475 :
476 : !=======================================================================
477 : ! Computes quantities needed in the stress tensor (sigma)
478 : ! and momentum (u) equations, but which do not change during
479 : ! the thermodynamics/transport time step:
480 : ! --wind stress shift to U grid,
481 : ! --ice mass and ice extent masks,
482 : ! initializes ice velocity for new points to ocean sfc current
483 : !
484 : ! author: Elizabeth C. Hunke, LANL
485 :
486 22944 : subroutine dyn_prep2 (nx_block, ny_block, &
487 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
488 : icellT, icellX, & ! LCOV_EXCL_LINE
489 : indxTi, indxTj, & ! LCOV_EXCL_LINE
490 : indxXi, indxXj, & ! LCOV_EXCL_LINE
491 : aiX, Xmass, & ! LCOV_EXCL_LINE
492 : Xmassdti, fcor, & ! LCOV_EXCL_LINE
493 : Xmask, & ! LCOV_EXCL_LINE
494 : uocn, vocn, & ! LCOV_EXCL_LINE
495 : strairx, strairy, & ! LCOV_EXCL_LINE
496 : ss_tltx, ss_tlty, & ! LCOV_EXCL_LINE
497 : iceTmask, iceXmask, & ! LCOV_EXCL_LINE
498 : fm, dt, & ! LCOV_EXCL_LINE
499 : strtltx, strtlty, & ! LCOV_EXCL_LINE
500 : strocnx, strocny, & ! LCOV_EXCL_LINE
501 : strintx, strinty, & ! LCOV_EXCL_LINE
502 : taubx, tauby, & ! LCOV_EXCL_LINE
503 : waterx, watery, & ! LCOV_EXCL_LINE
504 : forcex, forcey, & ! LCOV_EXCL_LINE
505 : stressp_1, stressp_2, & ! LCOV_EXCL_LINE
506 : stressp_3, stressp_4, & ! LCOV_EXCL_LINE
507 : stressm_1, stressm_2, & ! LCOV_EXCL_LINE
508 : stressm_3, stressm_4, & ! LCOV_EXCL_LINE
509 : stress12_1, stress12_2, & ! LCOV_EXCL_LINE
510 : stress12_3, stress12_4, & ! LCOV_EXCL_LINE
511 : uvel_init, vvel_init, & ! LCOV_EXCL_LINE
512 : uvel, vvel, & ! LCOV_EXCL_LINE
513 22944 : TbU)
514 :
515 : integer (kind=int_kind), intent(in) :: &
516 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
517 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
518 :
519 : integer (kind=int_kind), intent(out) :: &
520 : icellT , & ! no. of cells where iceTmask = .true. ! LCOV_EXCL_LINE
521 : icellX ! no. of cells where iceXmask = .true.
522 :
523 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(out) :: &
524 : indxTi , & ! compressed index in i-direction on T grid ! LCOV_EXCL_LINE
525 : indxTj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
526 : indxXi , & ! compressed index in i-direction on X grid, grid depends on call ! LCOV_EXCL_LINE
527 : indxXj ! compressed index in j-direction
528 :
529 : logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: &
530 : Xmask ! land/boundary mask, thickness (X-grid-cell)
531 :
532 : logical (kind=log_kind), dimension (nx_block,ny_block), intent(in) :: &
533 : iceTmask ! ice extent mask (T-cell)
534 :
535 : logical (kind=log_kind), dimension (nx_block,ny_block), intent(inout) :: &
536 : iceXmask ! ice extent mask (X-grid-cell)
537 :
538 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
539 : aiX , & ! ice fraction on u-grid (X grid) ! LCOV_EXCL_LINE
540 : Xmass , & ! total mass of ice and snow (X grid) ! LCOV_EXCL_LINE
541 : fcor , & ! Coriolis parameter (1/s) ! LCOV_EXCL_LINE
542 : strairx , & ! stress on ice by air, x-direction ! LCOV_EXCL_LINE
543 : strairy , & ! stress on ice by air, y-direction ! LCOV_EXCL_LINE
544 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
545 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
546 : ss_tltx , & ! sea surface slope, x-direction (m/m) ! LCOV_EXCL_LINE
547 : ss_tlty ! sea surface slope, y-direction
548 :
549 : real (kind=dbl_kind), intent(in) :: &
550 : dt ! time step
551 :
552 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out) :: &
553 : TbU, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
554 : uvel_init,& ! x-component of velocity (m/s), beginning of time step ! LCOV_EXCL_LINE
555 : vvel_init,& ! y-component of velocity (m/s), beginning of time step ! LCOV_EXCL_LINE
556 : Xmassdti, & ! mass of X-grid-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
557 : waterx , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
558 : watery , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
559 : forcex , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
560 : forcey ! work array: combined atm stress and ocn tilt, y
561 :
562 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
563 : fm , & ! Coriolis param. * mass in U-cell (kg/s) ! LCOV_EXCL_LINE
564 : stressp_1, stressp_2, stressp_3, stressp_4 , & ! sigma11+sigma22 ! LCOV_EXCL_LINE
565 : stressm_1, stressm_2, stressm_3, stressm_4 , & ! sigma11-sigma22 ! LCOV_EXCL_LINE
566 : stress12_1,stress12_2,stress12_3,stress12_4, & ! sigma12 ! LCOV_EXCL_LINE
567 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
568 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
569 : strtltx , & ! stress due to sea surface slope, x-direction ! LCOV_EXCL_LINE
570 : strtlty , & ! stress due to sea surface slope, y-direction ! LCOV_EXCL_LINE
571 : strocnx , & ! ice-ocean stress, x-direction ! LCOV_EXCL_LINE
572 : strocny , & ! ice-ocean stress, y-direction ! LCOV_EXCL_LINE
573 : strintx , & ! divergence of internal ice stress, x (N/m^2) ! LCOV_EXCL_LINE
574 : strinty , & ! divergence of internal ice stress, y (N/m^2) ! LCOV_EXCL_LINE
575 : taubx , & ! seabed stress, x-direction (N/m^2) ! LCOV_EXCL_LINE
576 : tauby ! seabed stress, y-direction (N/m^2)
577 :
578 : ! local variables
579 :
580 : integer (kind=int_kind) :: &
581 : i, j, ij
582 :
583 : real (kind=dbl_kind) :: &
584 5760 : gravit
585 :
586 : logical (kind=log_kind), dimension(nx_block,ny_block) :: &
587 45888 : iceXmask_old ! old-time iceXmask
588 :
589 : character(len=*), parameter :: subname = '(dyn_prep2)'
590 :
591 : !-----------------------------------------------------------------
592 : ! Initialize
593 : !-----------------------------------------------------------------
594 :
595 753576 : do j = 1, ny_block
596 16186320 : do i = 1, nx_block
597 15432744 : waterx (i,j) = c0
598 15432744 : watery (i,j) = c0
599 15432744 : forcex (i,j) = c0
600 15432744 : forcey (i,j) = c0
601 15432744 : Xmassdti (i,j) = c0
602 15432744 : TbU (i,j) = c0
603 15432744 : taubx (i,j) = c0
604 15432744 : tauby (i,j) = c0
605 :
606 16163376 : if (.not.iceTmask(i,j)) then
607 7916552 : stressp_1 (i,j) = c0
608 7916552 : stressp_2 (i,j) = c0
609 7916552 : stressp_3 (i,j) = c0
610 7916552 : stressp_4 (i,j) = c0
611 7916552 : stressm_1 (i,j) = c0
612 7916552 : stressm_2 (i,j) = c0
613 7916552 : stressm_3 (i,j) = c0
614 7916552 : stressm_4 (i,j) = c0
615 7916552 : stress12_1(i,j) = c0
616 7916552 : stress12_2(i,j) = c0
617 7916552 : stress12_3(i,j) = c0
618 7916552 : stress12_4(i,j) = c0
619 : endif
620 : enddo ! i
621 : enddo ! j
622 :
623 : !-----------------------------------------------------------------
624 : ! Identify cells where iceTmask = .true.
625 : ! Note: The icellT mask includes north and east ghost cells
626 : ! where stresses are needed.
627 : !-----------------------------------------------------------------
628 :
629 22944 : icellT = 0
630 730632 : do j = jlo, jhi+1
631 14983680 : do i = ilo, ihi+1
632 14960736 : if (iceTmask(i,j)) then
633 7138287 : icellT = icellT + 1
634 7138287 : indxTi(icellT) = i
635 7138287 : indxTj(icellT) = j
636 : endif
637 : enddo
638 : enddo
639 :
640 : !-----------------------------------------------------------------
641 : ! Define iceXmask
642 : ! Identify cells where iceXmask is true
643 : ! Initialize velocity where needed
644 : !-----------------------------------------------------------------
645 :
646 22944 : icellX = 0
647 :
648 707688 : do j = jlo, jhi
649 13826928 : do i = ilo, ihi
650 13119240 : iceXmask_old(i,j) = iceXmask(i,j) ! save
651 : ! ice extent mask (U-cells)
652 4176000 : iceXmask(i,j) = (Xmask(i,j)) .and. (aiX (i,j) > a_min) &
653 13119240 : .and. (Xmass(i,j) > m_min)
654 :
655 13119240 : if (iceXmask(i,j)) then
656 6452519 : icellX = icellX + 1
657 6452519 : indxXi(icellX) = i
658 6452519 : indxXj(icellX) = j
659 :
660 : ! initialize velocity for new ice points to ocean sfc current
661 6452519 : if (.not. iceXmask_old(i,j)) then
662 15917 : uvel(i,j) = uocn(i,j)
663 15917 : vvel(i,j) = vocn(i,j)
664 : endif
665 : else
666 : ! set velocity and stresses to zero for masked-out points
667 6666721 : uvel(i,j) = c0
668 6666721 : vvel(i,j) = c0
669 6666721 : strintx(i,j) = c0
670 6666721 : strinty(i,j) = c0
671 6666721 : strocnx(i,j) = c0
672 6666721 : strocny(i,j) = c0
673 : endif
674 :
675 13119240 : uvel_init(i,j) = uvel(i,j)
676 13803984 : vvel_init(i,j) = vvel(i,j)
677 : enddo
678 : enddo
679 :
680 : !-----------------------------------------------------------------
681 : ! Define variables for momentum equation
682 : !-----------------------------------------------------------------
683 :
684 22944 : if (trim(ssh_stress) == 'coupled') then
685 0 : call icepack_query_parameters(gravit_out=gravit)
686 0 : call icepack_warnings_flush(nu_diag)
687 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
688 0 : file=__FILE__, line=__LINE__)
689 : endif
690 :
691 6475463 : do ij = 1, icellX
692 6452519 : i = indxXi(ij)
693 6452519 : j = indxXj(ij)
694 :
695 6452519 : Xmassdti(i,j) = Xmass(i,j)/dt ! kg/m^2 s
696 :
697 6452519 : fm(i,j) = fcor(i,j)*Xmass(i,j) ! Coriolis * mass
698 :
699 : ! for ocean stress
700 6452519 : waterx(i,j) = uocn(i,j)*cosw - vocn(i,j)*sinw*sign(c1,fm(i,j))
701 6452519 : watery(i,j) = vocn(i,j)*cosw + uocn(i,j)*sinw*sign(c1,fm(i,j))
702 :
703 : ! combine tilt with wind stress
704 6452519 : if (trim(ssh_stress) == 'geostrophic') then
705 : ! calculate tilt from geostrophic currents if needed
706 6452519 : strtltx(i,j) = -fm(i,j)*vocn(i,j)
707 6452519 : strtlty(i,j) = fm(i,j)*uocn(i,j)
708 0 : elseif (trim(ssh_stress) == 'coupled') then
709 0 : strtltx(i,j) = -gravit*Xmass(i,j)*ss_tltx(i,j)
710 0 : strtlty(i,j) = -gravit*Xmass(i,j)*ss_tlty(i,j)
711 : else
712 : call abort_ice(subname//' ERROR: unknown ssh_stress='//trim(ssh_stress), &
713 0 : file=__FILE__, line=__LINE__)
714 : endif
715 :
716 6452519 : forcex(i,j) = strairx(i,j) + strtltx(i,j)
717 6475463 : forcey(i,j) = strairy(i,j) + strtlty(i,j)
718 : enddo
719 :
720 22944 : end subroutine dyn_prep2
721 :
722 : !=======================================================================
723 : ! Calculation of the surface stresses
724 : ! Integration of the momentum equation to find velocity (u,v)
725 : !
726 : ! author: Elizabeth C. Hunke, LANL
727 :
728 5506560 : subroutine stepu (nx_block, ny_block, &
729 : icellU, Cw, & ! LCOV_EXCL_LINE
730 : indxUi, indxUj, & ! LCOV_EXCL_LINE
731 : aiX, str, & ! LCOV_EXCL_LINE
732 : uocn, vocn, & ! LCOV_EXCL_LINE
733 : waterx, watery, & ! LCOV_EXCL_LINE
734 : forcex, forcey, & ! LCOV_EXCL_LINE
735 : Umassdti, fm, & ! LCOV_EXCL_LINE
736 : uarear, & ! LCOV_EXCL_LINE
737 : strintx, strinty, & ! LCOV_EXCL_LINE
738 : taubx, tauby, & ! LCOV_EXCL_LINE
739 : uvel_init, vvel_init,& ! LCOV_EXCL_LINE
740 : uvel, vvel, & ! LCOV_EXCL_LINE
741 5506560 : TbU)
742 :
743 : integer (kind=int_kind), intent(in) :: &
744 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
745 : icellU ! total count when iceUmask is true
746 :
747 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
748 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
749 : indxUj ! compressed index in j-direction
750 :
751 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
752 : TbU, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
753 : uvel_init,& ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
754 : vvel_init,& ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
755 : aiX , & ! ice fraction on X-grid ! LCOV_EXCL_LINE
756 : waterx , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
757 : watery , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
758 : forcex , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
759 : forcey , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
760 : Umassdti, & ! mass of U-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
761 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
762 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
763 : fm , & ! Coriolis param. * mass in U-cell (kg/s) ! LCOV_EXCL_LINE
764 : uarear ! 1/uarea
765 :
766 : real (kind=dbl_kind), dimension(nx_block,ny_block,8), intent(in) :: &
767 : str ! temporary
768 :
769 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
770 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
771 : vvel ! y-component of velocity (m/s)
772 :
773 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
774 : strintx , & ! divergence of internal ice stress, x (N/m^2) ! LCOV_EXCL_LINE
775 : strinty , & ! divergence of internal ice stress, y (N/m^2) ! LCOV_EXCL_LINE
776 : taubx , & ! seabed stress, x-direction (N/m^2) ! LCOV_EXCL_LINE
777 : tauby ! seabed stress, y-direction (N/m^2)
778 :
779 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
780 : Cw ! ocean-ice neutral drag coefficient
781 :
782 : ! local variables
783 :
784 : integer (kind=int_kind) :: &
785 : i, j, ij
786 :
787 : real (kind=dbl_kind) :: &
788 : uold, vold , & ! old-time uvel, vvel ! LCOV_EXCL_LINE
789 : vrel , & ! relative ice-ocean velocity ! LCOV_EXCL_LINE
790 : cca,ccb,ab2,cc1,cc2,& ! intermediate variables ! LCOV_EXCL_LINE
791 : taux, tauy , & ! part of ocean stress term ! LCOV_EXCL_LINE
792 : Cb , & ! complete seabed (basal) stress coeff ! LCOV_EXCL_LINE
793 1382400 : rhow !
794 :
795 : character(len=*), parameter :: subname = '(stepu)'
796 :
797 : !-----------------------------------------------------------------
798 : ! integrate the momentum equation
799 : !-----------------------------------------------------------------
800 :
801 5506560 : call icepack_query_parameters(rhow_out=rhow)
802 5506560 : call icepack_warnings_flush(nu_diag)
803 5506560 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
804 0 : file=__FILE__, line=__LINE__)
805 :
806 1554111120 : do ij =1, icellU
807 1548604560 : i = indxUi(ij)
808 1548604560 : j = indxUj(ij)
809 :
810 1548604560 : uold = uvel(i,j)
811 1548604560 : vold = vvel(i,j)
812 :
813 : ! (magnitude of relative ocean current)*rhow*drag*aice
814 113370480 : vrel = aiX(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uold)**2 + &
815 1548604560 : (vocn(i,j) - vold)**2) ! m/s
816 : ! ice/ocean stress
817 1548604560 : taux = vrel*waterx(i,j) ! NOTE this is not the entire
818 1548604560 : tauy = vrel*watery(i,j) ! ocn stress term
819 :
820 1548604560 : Cb = TbU(i,j) / (sqrt(uold**2 + vold**2) + u0) ! for seabed stress
821 : ! revp = 0 for classic evp, 1 for revised evp
822 1548604560 : cca = (brlx + revp)*Umassdti(i,j) + vrel * cosw + Cb ! kg/m^2 s
823 :
824 1548604560 : ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s
825 :
826 1548604560 : ab2 = cca**2 + ccb**2
827 :
828 : ! divergence of the internal stress tensor
829 113370480 : strintx(i,j) = uarear(i,j)* &
830 1548604560 : (str(i,j,1) + str(i+1,j,2) + str(i,j+1,3) + str(i+1,j+1,4))
831 113370480 : strinty(i,j) = uarear(i,j)* &
832 1548604560 : (str(i,j,5) + str(i,j+1,6) + str(i+1,j,7) + str(i+1,j+1,8))
833 :
834 : ! finally, the velocity components
835 113370480 : cc1 = strintx(i,j) + forcex(i,j) + taux &
836 1548604560 : + Umassdti(i,j)*(brlx*uold + revp*uvel_init(i,j))
837 113370480 : cc2 = strinty(i,j) + forcey(i,j) + tauy &
838 1548604560 : + Umassdti(i,j)*(brlx*vold + revp*vvel_init(i,j))
839 :
840 1548604560 : uvel(i,j) = (cca*cc1 + ccb*cc2) / ab2 ! m/s
841 1548604560 : vvel(i,j) = (cca*cc2 - ccb*cc1) / ab2
842 :
843 : ! calculate seabed stress component for outputs
844 : ! only needed on last iteration.
845 1548604560 : taubx(i,j) = -uvel(i,j)*Cb
846 1554111120 : tauby(i,j) = -vvel(i,j)*Cb
847 : enddo ! ij
848 :
849 5506560 : end subroutine stepu
850 :
851 : !=======================================================================
852 : ! Integration of the momentum equation to find velocity (u,v) at E and N locations
853 :
854 0 : subroutine stepuv_CD (nx_block, ny_block, &
855 : icell, Cw, & ! LCOV_EXCL_LINE
856 : indxi, indxj, & ! LCOV_EXCL_LINE
857 : aiX, & ! LCOV_EXCL_LINE
858 : uocn, vocn, & ! LCOV_EXCL_LINE
859 : waterx, watery, & ! LCOV_EXCL_LINE
860 : forcex, forcey, & ! LCOV_EXCL_LINE
861 : massdti, fm, & ! LCOV_EXCL_LINE
862 : strintx, strinty, & ! LCOV_EXCL_LINE
863 : taubx, tauby, & ! LCOV_EXCL_LINE
864 : uvel_init, vvel_init,& ! LCOV_EXCL_LINE
865 : uvel, vvel, & ! LCOV_EXCL_LINE
866 0 : Tb)
867 :
868 : integer (kind=int_kind), intent(in) :: &
869 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
870 : icell ! total count when ice[en]mask is true
871 :
872 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
873 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
874 : indxj ! compressed index in j-direction
875 :
876 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
877 : Tb, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
878 : uvel_init,& ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
879 : vvel_init,& ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
880 : aiX , & ! ice fraction on X-grid ! LCOV_EXCL_LINE
881 : waterx , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
882 : watery , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
883 : forcex , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
884 : forcey , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
885 : massdti , & ! mass of [EN]-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
886 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
887 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
888 : fm , & ! Coriolis param. * mass in [EN]-cell (kg/s) ! LCOV_EXCL_LINE
889 : strintx , & ! divergence of internal ice stress, x (N/m^2) ! LCOV_EXCL_LINE
890 : strinty ! divergence of internal ice stress, y (N/m^2)
891 :
892 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
893 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
894 : vvel ! y-component of velocity (m/s)
895 :
896 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
897 : taubx , & ! seabed stress, x-direction (N/m^2) ! LCOV_EXCL_LINE
898 : tauby ! seabed stress, y-direction (N/m^2)
899 :
900 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
901 : Cw ! ocean-ice neutral drag coefficient
902 :
903 : ! local variables
904 :
905 : integer (kind=int_kind) :: &
906 : i, j, ij
907 :
908 : real (kind=dbl_kind) :: &
909 : uold, vold , & ! old-time uvel, vvel ! LCOV_EXCL_LINE
910 : vrel , & ! relative ice-ocean velocity ! LCOV_EXCL_LINE
911 : cca,ccb,ccc,ab2 , & ! intermediate variables ! LCOV_EXCL_LINE
912 : cc1,cc2 , & ! " ! LCOV_EXCL_LINE
913 : taux, tauy , & ! part of ocean stress term ! LCOV_EXCL_LINE
914 : Cb , & ! complete seabed (basal) stress coeff ! LCOV_EXCL_LINE
915 0 : rhow !
916 :
917 : character(len=*), parameter :: subname = '(stepuv_CD)'
918 :
919 : !-----------------------------------------------------------------
920 : ! integrate the momentum equation
921 : !-----------------------------------------------------------------
922 :
923 0 : call icepack_query_parameters(rhow_out=rhow)
924 0 : call icepack_warnings_flush(nu_diag)
925 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
926 0 : file=__FILE__, line=__LINE__)
927 :
928 0 : do ij =1, icell
929 0 : i = indxi(ij)
930 0 : j = indxj(ij)
931 :
932 0 : uold = uvel(i,j)
933 0 : vold = vvel(i,j)
934 :
935 : ! (magnitude of relative ocean current)*rhow*drag*aice
936 0 : vrel = aiX(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uold)**2 + &
937 0 : (vocn(i,j) - vold)**2) ! m/s
938 : ! ice/ocean stress
939 0 : taux = vrel*waterx(i,j) ! NOTE this is not the entire
940 0 : tauy = vrel*watery(i,j) ! ocn stress term
941 :
942 0 : ccc = sqrt(uold**2 + vold**2) + u0
943 0 : Cb = Tb(i,j) / ccc ! for seabed stress
944 : ! revp = 0 for classic evp, 1 for revised evp
945 0 : cca = (brlx + revp)*massdti(i,j) + vrel * cosw + Cb ! kg/m^2 s
946 :
947 0 : ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s
948 :
949 0 : ab2 = cca**2 + ccb**2
950 :
951 : ! compute the velocity components
952 0 : cc1 = strintx(i,j) + forcex(i,j) + taux &
953 0 : + massdti(i,j)*(brlx*uold + revp*uvel_init(i,j))
954 0 : cc2 = strinty(i,j) + forcey(i,j) + tauy &
955 0 : + massdti(i,j)*(brlx*vold + revp*vvel_init(i,j))
956 0 : uvel(i,j) = (cca*cc1 + ccb*cc2) / ab2 ! m/s
957 0 : vvel(i,j) = (cca*cc2 - ccb*cc1) / ab2
958 :
959 : ! calculate seabed stress component for outputs
960 : ! only needed on last iteration.
961 0 : taubx(i,j) = -uvel(i,j)*Cb
962 0 : tauby(i,j) = -vvel(i,j)*Cb
963 :
964 : enddo ! ij
965 :
966 0 : end subroutine stepuv_CD
967 :
968 : !=======================================================================
969 : ! Integration of the momentum equation to find velocity u at E location on C grid
970 :
971 0 : subroutine stepu_C (nx_block, ny_block, &
972 : icell, Cw, & ! LCOV_EXCL_LINE
973 : indxi, indxj, & ! LCOV_EXCL_LINE
974 : aiX, & ! LCOV_EXCL_LINE
975 : uocn, vocn, & ! LCOV_EXCL_LINE
976 : waterx, forcex, & ! LCOV_EXCL_LINE
977 : massdti, fm, & ! LCOV_EXCL_LINE
978 : strintx, taubx, & ! LCOV_EXCL_LINE
979 : uvel_init, & ! LCOV_EXCL_LINE
980 : uvel, vvel, & ! LCOV_EXCL_LINE
981 0 : Tb)
982 :
983 : integer (kind=int_kind), intent(in) :: &
984 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
985 : icell ! total count when ice[en]mask is true
986 :
987 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
988 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
989 : indxj ! compressed index in j-direction
990 :
991 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
992 : Tb, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
993 : uvel_init,& ! x-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
994 : aiX , & ! ice fraction on X-grid ! LCOV_EXCL_LINE
995 : waterx , & ! for ocean stress calculation, x (m/s) ! LCOV_EXCL_LINE
996 : forcex , & ! work array: combined atm stress and ocn tilt, x ! LCOV_EXCL_LINE
997 : massdti , & ! mass of e-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
998 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
999 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
1000 : fm , & ! Coriolis param. * mass in e-cell (kg/s) ! LCOV_EXCL_LINE
1001 : strintx , & ! divergence of internal ice stress, x (N/m^2) ! LCOV_EXCL_LINE
1002 : Cw , & ! ocean-ice neutral drag coefficient ! LCOV_EXCL_LINE
1003 : vvel ! y-component of velocity (m/s) interpolated to E location
1004 :
1005 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1006 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1007 : taubx ! seabed stress, x-direction (N/m^2)
1008 :
1009 : ! local variables
1010 :
1011 : integer (kind=int_kind) :: &
1012 : i, j, ij
1013 :
1014 : real (kind=dbl_kind) :: &
1015 : uold, vold , & ! old-time uvel, vvel ! LCOV_EXCL_LINE
1016 : vrel , & ! relative ice-ocean velocity ! LCOV_EXCL_LINE
1017 : cca,ccb,ccc,cc1 , & ! intermediate variables ! LCOV_EXCL_LINE
1018 : taux , & ! part of ocean stress term ! LCOV_EXCL_LINE
1019 : Cb , & ! complete seabed (basal) stress coeff ! LCOV_EXCL_LINE
1020 0 : rhow !
1021 :
1022 : character(len=*), parameter :: subname = '(stepu_C)'
1023 :
1024 : !-----------------------------------------------------------------
1025 : ! integrate the momentum equation
1026 : !-----------------------------------------------------------------
1027 :
1028 0 : call icepack_query_parameters(rhow_out=rhow)
1029 0 : call icepack_warnings_flush(nu_diag)
1030 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1031 0 : file=__FILE__, line=__LINE__)
1032 :
1033 0 : do ij =1, icell
1034 0 : i = indxi(ij)
1035 0 : j = indxj(ij)
1036 :
1037 0 : uold = uvel(i,j)
1038 0 : vold = vvel(i,j)
1039 :
1040 : ! (magnitude of relative ocean current)*rhow*drag*aice
1041 0 : vrel = aiX(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uold)**2 + &
1042 0 : (vocn(i,j) - vold)**2) ! m/s
1043 : ! ice/ocean stress
1044 0 : taux = vrel*waterx(i,j) ! NOTE this is not the entire
1045 :
1046 0 : ccc = sqrt(uold**2 + vold**2) + u0
1047 0 : Cb = Tb(i,j) / ccc ! for seabed stress
1048 : ! revp = 0 for classic evp, 1 for revised evp
1049 0 : cca = (brlx + revp)*massdti(i,j) + vrel * cosw + Cb ! kg/m^2 s
1050 :
1051 0 : ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s
1052 :
1053 : ! compute the velocity components
1054 0 : cc1 = strintx(i,j) + forcex(i,j) + taux &
1055 0 : + massdti(i,j)*(brlx*uold + revp*uvel_init(i,j))
1056 :
1057 0 : uvel(i,j) = (ccb*vold + cc1) / cca ! m/s
1058 :
1059 : ! calculate seabed stress component for outputs
1060 : ! only needed on last iteration.
1061 0 : taubx(i,j) = -uvel(i,j)*Cb
1062 :
1063 : enddo ! ij
1064 :
1065 0 : end subroutine stepu_C
1066 :
1067 : !=======================================================================
1068 : ! Integration of the momentum equation to find velocity v at N location on C grid
1069 :
1070 0 : subroutine stepv_C (nx_block, ny_block, &
1071 : icell, Cw, & ! LCOV_EXCL_LINE
1072 : indxi, indxj, & ! LCOV_EXCL_LINE
1073 : aiX, & ! LCOV_EXCL_LINE
1074 : uocn, vocn, & ! LCOV_EXCL_LINE
1075 : watery, forcey, & ! LCOV_EXCL_LINE
1076 : massdti, fm, & ! LCOV_EXCL_LINE
1077 : strinty, tauby, & ! LCOV_EXCL_LINE
1078 : vvel_init, & ! LCOV_EXCL_LINE
1079 : uvel, vvel, & ! LCOV_EXCL_LINE
1080 0 : Tb)
1081 :
1082 : integer (kind=int_kind), intent(in) :: &
1083 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1084 : icell ! total count when ice[en]mask is true
1085 :
1086 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1087 : indxi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1088 : indxj ! compressed index in j-direction
1089 :
1090 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1091 : Tb, & ! seabed stress factor (N/m^2) ! LCOV_EXCL_LINE
1092 : vvel_init,& ! y-component of velocity (m/s), beginning of timestep ! LCOV_EXCL_LINE
1093 : aiX , & ! ice fraction on X-grid ! LCOV_EXCL_LINE
1094 : watery , & ! for ocean stress calculation, y (m/s) ! LCOV_EXCL_LINE
1095 : forcey , & ! work array: combined atm stress and ocn tilt, y ! LCOV_EXCL_LINE
1096 : massdti , & ! mass of n-cell/dt (kg/m^2 s) ! LCOV_EXCL_LINE
1097 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
1098 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
1099 : fm , & ! Coriolis param. * mass in n-cell (kg/s) ! LCOV_EXCL_LINE
1100 : strinty , & ! divergence of internal ice stress, y (N/m^2) ! LCOV_EXCL_LINE
1101 : Cw , & ! ocean-ice neutral drag coefficient ! LCOV_EXCL_LINE
1102 : uvel ! x-component of velocity (m/s) interpolated to N location
1103 :
1104 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1105 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1106 : tauby ! seabed stress, y-direction (N/m^2)
1107 :
1108 : ! local variables
1109 :
1110 : integer (kind=int_kind) :: &
1111 : i, j, ij
1112 :
1113 : real (kind=dbl_kind) :: &
1114 : uold, vold , & ! old-time uvel, vvel ! LCOV_EXCL_LINE
1115 : vrel , & ! relative ice-ocean velocity ! LCOV_EXCL_LINE
1116 : cca,ccb,ccc,cc2 , & ! intermediate variables ! LCOV_EXCL_LINE
1117 : tauy , & ! part of ocean stress term ! LCOV_EXCL_LINE
1118 : Cb , & ! complete seabed (basal) stress coeff ! LCOV_EXCL_LINE
1119 0 : rhow !
1120 :
1121 : character(len=*), parameter :: subname = '(stepv_C)'
1122 :
1123 : !-----------------------------------------------------------------
1124 : ! integrate the momentum equation
1125 : !-----------------------------------------------------------------
1126 :
1127 0 : call icepack_query_parameters(rhow_out=rhow)
1128 0 : call icepack_warnings_flush(nu_diag)
1129 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1130 0 : file=__FILE__, line=__LINE__)
1131 :
1132 0 : do ij =1, icell
1133 0 : i = indxi(ij)
1134 0 : j = indxj(ij)
1135 :
1136 0 : uold = uvel(i,j)
1137 0 : vold = vvel(i,j)
1138 :
1139 : ! (magnitude of relative ocean current)*rhow*drag*aice
1140 0 : vrel = aiX(i,j)*rhow*Cw(i,j)*sqrt((uocn(i,j) - uold)**2 + &
1141 0 : (vocn(i,j) - vold)**2) ! m/s
1142 : ! ice/ocean stress
1143 0 : tauy = vrel*watery(i,j) ! NOTE this is not the entire ocn stress
1144 :
1145 0 : ccc = sqrt(uold**2 + vold**2) + u0
1146 0 : Cb = Tb(i,j) / ccc ! for seabed stress
1147 : ! revp = 0 for classic evp, 1 for revised evp
1148 0 : cca = (brlx + revp)*massdti(i,j) + vrel * cosw + Cb ! kg/m^2 s
1149 :
1150 0 : ccb = fm(i,j) + sign(c1,fm(i,j)) * vrel * sinw ! kg/m^2 s
1151 :
1152 : ! compute the velocity components
1153 0 : cc2 = strinty(i,j) + forcey(i,j) + tauy &
1154 0 : + massdti(i,j)*(brlx*vold + revp*vvel_init(i,j))
1155 :
1156 0 : vvel(i,j) = (-ccb*uold + cc2) / cca
1157 :
1158 : ! calculate seabed stress component for outputs
1159 : ! only needed on last iteration.
1160 0 : tauby(i,j) = -vvel(i,j)*Cb
1161 :
1162 : enddo ! ij
1163 :
1164 0 : end subroutine stepv_C
1165 :
1166 : !=======================================================================
1167 : ! Calculation of the ice-ocean stress.
1168 : ! ...the sign will be reversed later...
1169 : !
1170 : ! author: Elizabeth C. Hunke, LANL
1171 :
1172 22944 : subroutine dyn_finish (nx_block, ny_block, &
1173 : icellU, Cw, & ! LCOV_EXCL_LINE
1174 : indxUi, indxUj, & ! LCOV_EXCL_LINE
1175 : uvel, vvel, & ! LCOV_EXCL_LINE
1176 : uocn, vocn, & ! LCOV_EXCL_LINE
1177 : aiX, fm, & ! LCOV_EXCL_LINE
1178 22944 : strocnx, strocny)
1179 :
1180 : integer (kind=int_kind), intent(in) :: &
1181 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1182 : icellU ! total count when iceUmask is true
1183 :
1184 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1185 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1186 : indxUj ! compressed index in j-direction
1187 :
1188 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1189 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1190 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1191 : uocn , & ! ocean current, x-direction (m/s) ! LCOV_EXCL_LINE
1192 : vocn , & ! ocean current, y-direction (m/s) ! LCOV_EXCL_LINE
1193 : aiX , & ! ice fraction on X-grid ! LCOV_EXCL_LINE
1194 : fm ! Coriolis param. * mass in U-cell (kg/s)
1195 :
1196 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1197 : strocnx , & ! ice-ocean stress, x-direction ! LCOV_EXCL_LINE
1198 : strocny ! ice-ocean stress, y-direction
1199 :
1200 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1201 : Cw ! ocean-ice neutral drag coefficient
1202 :
1203 : ! local variables
1204 :
1205 : integer (kind=int_kind) :: &
1206 : i, j, ij
1207 :
1208 : real (kind=dbl_kind) :: &
1209 : vrel , & ! ! LCOV_EXCL_LINE
1210 5760 : rhow !
1211 :
1212 : character(len=*), parameter :: subname = '(dyn_finish)'
1213 :
1214 22944 : call icepack_query_parameters(rhow_out=rhow)
1215 22944 : call icepack_warnings_flush(nu_diag)
1216 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1217 0 : file=__FILE__, line=__LINE__)
1218 :
1219 : ! ocean-ice stress for coupling
1220 6475463 : do ij =1, icellU
1221 6452519 : i = indxUi(ij)
1222 6452519 : j = indxUj(ij)
1223 :
1224 472377 : vrel = rhow*Cw(i,j)*sqrt((uocn(i,j) - uvel(i,j))**2 + &
1225 6452519 : (vocn(i,j) - vvel(i,j))**2) ! m/s
1226 :
1227 : ! strocnx(i,j) = strocnx(i,j) &
1228 : ! - vrel*(uvel(i,j)*cosw - vvel(i,j)*sinw) * aiX(i,j)
1229 : ! strocny(i,j) = strocny(i,j) &
1230 : ! - vrel*(vvel(i,j)*cosw + uvel(i,j)*sinw) * aiX(i,j)
1231 :
1232 : ! update strocnx to most recent iterate and complete the term
1233 6452519 : vrel = vrel * aiX(i,j)
1234 472377 : strocnx(i,j) = vrel*((uocn(i,j) - uvel(i,j))*cosw &
1235 6452519 : - (vocn(i,j) - vvel(i,j))*sinw*sign(c1,fm(i,j)))
1236 472377 : strocny(i,j) = vrel*((vocn(i,j) - vvel(i,j))*cosw &
1237 6475463 : + (uocn(i,j) - uvel(i,j))*sinw*sign(c1,fm(i,j)))
1238 :
1239 : ! Hibler/Bryan stress
1240 : ! the sign is reversed later, therefore negative here
1241 : ! strocnx(i,j) = -(strairx(i,j) + strintx(i,j))
1242 : ! strocny(i,j) = -(strairy(i,j) + strinty(i,j))
1243 :
1244 : enddo
1245 :
1246 22944 : end subroutine dyn_finish
1247 :
1248 : !=======================================================================
1249 : ! Computes seabed (basal) stress factor TbU (landfast ice) based on mean
1250 : ! thickness and bathymetry data. LKD refers to linear keel draft. This
1251 : ! parameterization assumes that the largest keel draft varies linearly
1252 : ! with the mean thickness.
1253 : !
1254 : ! Lemieux, J. F., B. Tremblay, F. Dupont, M. Plante, G.C. Smith, D. Dumont (2015).
1255 : ! A basal stress parameterization form modeling landfast ice, J. Geophys. Res.
1256 : ! Oceans, 120, 3157-3173.
1257 : !
1258 : ! Lemieux, J. F., F. Dupont, P. Blain, F. Roy, G.C. Smith, G.M. Flato (2016).
1259 : ! Improving the simulation of landfast ice by combining tensile strength and a
1260 : ! parameterization for grounded ridges, J. Geophys. Res. Oceans, 121, 7354-7368.
1261 : !
1262 : ! author: JF Lemieux, Philippe Blain (ECCC)
1263 : !
1264 : ! note1: TbU is a part of the Cb as defined in Lemieux et al. 2015 and 2016.
1265 : ! note2: Seabed stress (better name) was called basal stress in Lemieux et al. 2015
1266 :
1267 0 : subroutine seabed_stress_factor_LKD (nx_block, ny_block, &
1268 : icellU, & ! LCOV_EXCL_LINE
1269 : indxUi, indxUj, & ! LCOV_EXCL_LINE
1270 : vice, aice, & ! LCOV_EXCL_LINE
1271 : hwater, TbU, & ! LCOV_EXCL_LINE
1272 : grid_location)
1273 :
1274 : use ice_grid, only: grid_neighbor_min, grid_neighbor_max
1275 :
1276 : integer (kind=int_kind), intent(in) :: &
1277 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1278 : icellU ! no. of cells where ice[uen]mask = 1
1279 :
1280 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1281 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1282 : indxUj ! compressed index in j-direction
1283 :
1284 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1285 : aice , & ! concentration of ice at tracer location ! LCOV_EXCL_LINE
1286 : vice , & ! volume per unit area of ice at tracer location (m) ! LCOV_EXCL_LINE
1287 : hwater ! water depth at tracer location (m)
1288 :
1289 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1290 : TbU ! seabed stress factor at 'grid_location' (N/m^2)
1291 :
1292 : character(len=*), optional, intent(inout) :: &
1293 : grid_location ! grid location (U, E, N), U assumed if not present
1294 :
1295 : real (kind=dbl_kind) :: &
1296 : au , & ! concentration of ice at u location ! LCOV_EXCL_LINE
1297 : hu , & ! volume per unit area of ice at u location (mean thickness, m) ! LCOV_EXCL_LINE
1298 : hwu , & ! water depth at u location (m) ! LCOV_EXCL_LINE
1299 : docalc_tbu, & ! logical as real (C0,C1) decides whether c0 is 0 or ! LCOV_EXCL_LINE
1300 0 : hcu ! critical thickness at u location (m)
1301 :
1302 : integer (kind=int_kind) :: &
1303 : i, j, ij
1304 :
1305 : character(len=char_len) :: &
1306 : l_grid_location ! local version of 'grid_location'
1307 :
1308 : character(len=*), parameter :: subname = '(seabed_stress_factor_LKD)'
1309 :
1310 : ! Assume U location (NE corner) if grid_location not present
1311 0 : if (.not. (present(grid_location))) then
1312 0 : l_grid_location = 'U'
1313 : else
1314 0 : l_grid_location = grid_location
1315 : endif
1316 :
1317 0 : do ij = 1, icellU
1318 0 : i = indxUi(ij)
1319 0 : j = indxUj(ij)
1320 :
1321 : ! convert quantities to grid_location
1322 :
1323 0 : hwu = grid_neighbor_min(hwater, i, j, l_grid_location)
1324 :
1325 0 : docalc_tbu = merge(c1,c0,hwu < threshold_hw)
1326 :
1327 :
1328 0 : au = grid_neighbor_max(aice, i, j, l_grid_location)
1329 0 : hu = grid_neighbor_max(vice, i, j, l_grid_location)
1330 :
1331 : ! 1- calculate critical thickness
1332 0 : hcu = au * hwu / k1
1333 :
1334 : ! 2- calculate seabed stress factor
1335 0 : TbU(i,j) = docalc_tbu*k2 * max(c0,(hu - hcu)) * exp(-alphab * (c1 - au))
1336 :
1337 : enddo ! ij
1338 :
1339 0 : end subroutine seabed_stress_factor_LKD
1340 :
1341 : !=======================================================================
1342 : ! Computes seabed (basal) stress factor TbU (landfast ice) based on
1343 : ! probability of contact between the ITD and the seabed. The water depth
1344 : ! could take into account variations of the SSH. In the simplest
1345 : ! formulation, hwater is simply the value of the bathymetry. To calculate
1346 : ! the probability of contact, it is assumed that the bathymetry follows
1347 : ! a normal distribution with sigma_b = 2.5d0. An improvement would
1348 : ! be to provide the distribution based on high resolution data.
1349 : !
1350 : ! Dupont, F., D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, A. Caya (2022).
1351 : ! A probabilistic seabed-ice keel interaction model, The Cryosphere, 16,
1352 : ! 1963-1977.
1353 : !
1354 : ! authors: D. Dumont, J.F. Lemieux, E. Dumas-Lefebvre, F. Dupont
1355 : !
1356 0 : subroutine seabed_stress_factor_prob (nx_block, ny_block, &
1357 : icellT, indxTi, indxTj, & ! LCOV_EXCL_LINE
1358 : icellU, indxUi, indxUj, & ! LCOV_EXCL_LINE
1359 : aicen, vicen, & ! LCOV_EXCL_LINE
1360 : hwater, TbU, & ! LCOV_EXCL_LINE
1361 : TbE, TbN, & ! LCOV_EXCL_LINE
1362 : icellE, indxEi, indxEj, & ! LCOV_EXCL_LINE
1363 0 : icellN, indxNi, indxNj)
1364 : ! use modules
1365 :
1366 : use ice_arrays_column, only: hin_max
1367 : use ice_domain_size, only: ncat
1368 : use ice_grid, only: grid_neighbor_min, grid_neighbor_max
1369 :
1370 : integer (kind=int_kind), intent(in) :: &
1371 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1372 : icellT, icellU ! no. of cells where ice[tu]mask = 1
1373 :
1374 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1375 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1376 : indxTj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
1377 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1378 : indxUj ! compressed index in j-direction
1379 :
1380 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1381 : hwater ! water depth at tracer location (m)
1382 :
1383 : real (kind=dbl_kind), dimension (nx_block,ny_block,ncat), intent(in) :: &
1384 : aicen, & ! partial concentration for last thickness category in ITD ! LCOV_EXCL_LINE
1385 : vicen ! partial volume for last thickness category in ITD (m)
1386 :
1387 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1388 : TbU ! seabed stress factor at U location (N/m^2)
1389 :
1390 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout), optional :: &
1391 : TbE, & ! seabed stress factor at E location (N/m^2) ! LCOV_EXCL_LINE
1392 : TbN ! seabed stress factor at N location (N/m^2)
1393 :
1394 : integer (kind=int_kind), intent(in), optional :: &
1395 : icellE, icellN ! no. of cells where ice[en]mask = 1
1396 :
1397 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in), optional :: &
1398 : indxEi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1399 : indxEj , & ! compressed index in j-direction ! LCOV_EXCL_LINE
1400 : indxNi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1401 : indxNj ! compressed index in j-direction
1402 :
1403 : ! local variables
1404 :
1405 : integer (kind=int_kind) :: &
1406 : i, j, ij, ii, n
1407 :
1408 : integer, parameter :: &
1409 : ncat_b = 100, & ! number of bathymetry categories ! LCOV_EXCL_LINE
1410 : ncat_i = 100 ! number of ice thickness categories (log-normal)
1411 :
1412 : real (kind=dbl_kind), parameter :: &
1413 : max_depth = 50.0_dbl_kind, & ! initial range of log-normal distribution ! LCOV_EXCL_LINE
1414 : mu_s = 0.1_dbl_kind, & ! friction coefficient ! LCOV_EXCL_LINE
1415 : sigma_b = 2.5d0 ! Standard deviation of bathymetry
1416 :
1417 : real (kind=dbl_kind), dimension(ncat_i) :: & ! log-normal for ice thickness
1418 : x_k, & ! center of thickness categories (m) ! LCOV_EXCL_LINE
1419 : g_k, & ! probability density function (thickness, 1/m) ! LCOV_EXCL_LINE
1420 0 : P_x ! probability for each thickness category
1421 :
1422 : real (kind=dbl_kind), dimension(ncat_b) :: & ! normal dist for bathymetry
1423 : y_n, & ! center of bathymetry categories (m) ! LCOV_EXCL_LINE
1424 : b_n, & ! probability density function (bathymetry, 1/m) ! LCOV_EXCL_LINE
1425 0 : P_y ! probability for each bathymetry category
1426 :
1427 : real (kind=dbl_kind), dimension(ncat) :: &
1428 0 : vcat, acat ! vice, aice temporary arrays
1429 :
1430 : integer, dimension(ncat_b) :: &
1431 : tmp ! Temporary vector tmp = merge(1,0,gt)
1432 :
1433 : logical, dimension (ncat_b) :: &
1434 : gt !
1435 :
1436 : real (kind=dbl_kind) :: &
1437 : wid_i, wid_b , & ! parameters for PDFs ! LCOV_EXCL_LINE
1438 : mu_i, sigma_i , & ! ! LCOV_EXCL_LINE
1439 : mu_b, m_i, v_i, & ! ! LCOV_EXCL_LINE
1440 : atot, x_kmax , & ! ! LCOV_EXCL_LINE
1441 : cut , & ! ! LCOV_EXCL_LINE
1442 : rhoi, rhow , & ! ! LCOV_EXCL_LINE
1443 : gravit , & ! ! LCOV_EXCL_LINE
1444 0 : pi, puny !
1445 :
1446 : real (kind=dbl_kind), dimension(ncat_i) :: &
1447 0 : tb_tmp
1448 :
1449 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1450 0 : Tbt ! seabed stress factor at t point (N/m^2)
1451 :
1452 : character(len=*), parameter :: subname = '(seabed_stress_factor_prob)'
1453 :
1454 0 : call icepack_query_parameters(rhow_out=rhow, rhoi_out=rhoi)
1455 0 : call icepack_query_parameters(gravit_out=gravit)
1456 0 : call icepack_query_parameters(pi_out=pi)
1457 0 : call icepack_query_parameters(puny_out=puny)
1458 :
1459 0 : Tbt=c0
1460 :
1461 0 : do ij = 1, icellT
1462 0 : i = indxTi(ij)
1463 0 : j = indxTj(ij)
1464 :
1465 0 : atot = sum(aicen(i,j,1:ncat))
1466 :
1467 0 : if (atot > 0.05_dbl_kind .and. hwater(i,j) < max_depth) then
1468 :
1469 0 : mu_b = hwater(i,j) ! mean of PDF (normal dist) bathymetry
1470 0 : wid_i = max_depth/ncat_i ! width of ice categories
1471 0 : wid_b = c6*sigma_b/ncat_b ! width of bathymetry categories (6 sigma_b = 2x3 sigma_b)
1472 :
1473 0 : x_k = (/( wid_i*( real(i,kind=dbl_kind) - p5 ), i=1, ncat_i )/)
1474 0 : y_n = (/( ( mu_b-c3*sigma_b )+( real(i,kind=dbl_kind) - p5 )*( c6*sigma_b/ncat_b ), i=1, ncat_b )/)
1475 :
1476 0 : vcat(1:ncat) = vicen(i,j,1:ncat)
1477 0 : acat(1:ncat) = aicen(i,j,1:ncat)
1478 :
1479 0 : m_i = sum(vcat)
1480 :
1481 0 : v_i=c0
1482 0 : do n =1, ncat
1483 0 : v_i = v_i + vcat(n)**2 / (max(acat(n), puny))
1484 : enddo
1485 0 : v_i = max((v_i - m_i**2), puny)
1486 :
1487 0 : mu_i = log(m_i/sqrt(c1 + v_i/m_i**2)) ! parameters for the log-normal
1488 0 : sigma_i = sqrt(log(c1 + v_i/m_i**2))
1489 :
1490 : ! max thickness associated with percentile of log-normal PDF
1491 : ! x_kmax=x997 was obtained from an optimization procedure (Dupont et al. 2022)
1492 :
1493 0 : x_kmax = exp(mu_i + sqrt(c2*sigma_i)*1.9430d0)
1494 :
1495 : ! Set x_kmax to hlev of the last category where there is ice
1496 : ! when there is no ice in the last category
1497 0 : cut = x_k(ncat_i)
1498 0 : do n = ncat,-1,1
1499 0 : if (acat(n) < puny) then
1500 0 : cut = hin_max(n-1)
1501 : else
1502 0 : exit
1503 : endif
1504 : enddo
1505 0 : x_kmax = min(cut, x_kmax)
1506 :
1507 0 : g_k = exp(-(log(x_k) - mu_i) ** 2 / (c2 * sigma_i ** 2)) / (x_k * sigma_i * sqrt(c2 * pi))
1508 :
1509 0 : b_n = exp(-(y_n - mu_b) ** 2 / (c2 * sigma_b ** 2)) / (sigma_b * sqrt(c2 * pi))
1510 :
1511 0 : P_x = g_k*wid_i
1512 0 : P_y = b_n*wid_b
1513 :
1514 0 : do n =1, ncat_i
1515 0 : if (x_k(n) > x_kmax) P_x(n)=c0
1516 : enddo
1517 :
1518 : ! calculate Tb factor at t-location
1519 0 : do n=1, ncat_i
1520 0 : gt = (y_n <= rhoi*x_k(n)/rhow)
1521 0 : tmp = merge(1,0,gt)
1522 0 : ii = sum(tmp)
1523 0 : if (ii == 0) then
1524 0 : tb_tmp(n) = c0
1525 : else
1526 0 : tb_tmp(n) = max(mu_s*gravit*P_x(n)*sum(P_y(1:ii)*(rhoi*x_k(n) - rhow*y_n(1:ii))),c0)
1527 : endif
1528 : enddo
1529 0 : Tbt(i,j) = sum(tb_tmp)*exp(-alphab * (c1 - atot))
1530 : endif
1531 : enddo
1532 :
1533 0 : if (grid_ice == "B") then
1534 0 : do ij = 1, icellU
1535 0 : i = indxUi(ij)
1536 0 : j = indxUj(ij)
1537 : ! convert quantities to U-location
1538 0 : TbU(i,j) = grid_neighbor_max(Tbt, i, j, 'U')
1539 : enddo ! ij
1540 0 : elseif (grid_ice == "C" .or. grid_ice == "CD") then
1541 : if (present(Tbe) .and. present(TbN) .and. &
1542 : present(icellE) .and. present(icellN) .and. & ! LCOV_EXCL_LINE
1543 : present(indxEi) .and. present(indxEj) .and. & ! LCOV_EXCL_LINE
1544 0 : present(indxNi) .and. present(indxNj)) then
1545 :
1546 0 : do ij = 1, icellE
1547 0 : i = indxEi(ij)
1548 0 : j = indxEj(ij)
1549 : ! convert quantities to E-location
1550 0 : TbE(i,j) = grid_neighbor_max(Tbt, i, j, 'E')
1551 : enddo
1552 0 : do ij = 1, icellN
1553 0 : i = indxNi(ij)
1554 0 : j = indxNj(ij)
1555 : ! convert quantities to N-location
1556 0 : TbN(i,j) = grid_neighbor_max(Tbt, i, j, 'N')
1557 : enddo
1558 :
1559 : else
1560 0 : call abort_ice(subname // ' insufficient number of arguments for grid_ice:' // grid_ice)
1561 : endif
1562 : endif
1563 :
1564 0 : end subroutine seabed_stress_factor_prob
1565 :
1566 : !=======================================================================
1567 : ! Computes principal stresses for comparison with the theoretical
1568 : ! yield curve
1569 : !
1570 : ! author: Elizabeth C. Hunke, LANL
1571 :
1572 636 : subroutine principal_stress(nx_block, ny_block, &
1573 : stressp, stressm, & ! LCOV_EXCL_LINE
1574 : stress12, strength, & ! LCOV_EXCL_LINE
1575 : sig1, sig2, & ! LCOV_EXCL_LINE
1576 636 : sigP)
1577 :
1578 : integer (kind=int_kind), intent(in) :: &
1579 : nx_block, ny_block ! block dimensions
1580 :
1581 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1582 : stressp , & ! sigma11 + sigma22 ! LCOV_EXCL_LINE
1583 : stressm , & ! sigma11 - sigma22 ! LCOV_EXCL_LINE
1584 : stress12 , & ! sigma12 ! LCOV_EXCL_LINE
1585 : strength ! for normalization of sig1 and sig2
1586 :
1587 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
1588 : sig1 , & ! normalized principal stress component ! LCOV_EXCL_LINE
1589 : sig2 , & ! normalized principal stress component ! LCOV_EXCL_LINE
1590 : sigP ! internal ice pressure (N/m)
1591 :
1592 : ! local variables
1593 :
1594 : integer (kind=int_kind) :: &
1595 : i, j
1596 :
1597 : real (kind=dbl_kind) :: &
1598 32 : puny
1599 :
1600 : character(len=*), parameter :: subname = '(principal_stress)'
1601 :
1602 636 : call icepack_query_parameters(puny_out=puny)
1603 636 : call icepack_warnings_flush(nu_diag)
1604 636 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1605 0 : file=__FILE__, line=__LINE__)
1606 :
1607 20622 : do j = 1, ny_block
1608 291680 : do i = 1, nx_block
1609 291044 : if (strength(i,j) > puny) then
1610 : ! ice internal pressure
1611 22166 : sigP(i,j) = -p5*stressp(i,j)
1612 :
1613 : ! normalized principal stresses
1614 0 : sig1(i,j) = (p5*(stressp(i,j) &
1615 : + sqrt(stressm(i,j)**2+c4*stress12(i,j)**2))) & ! LCOV_EXCL_LINE
1616 22166 : / strength(i,j)
1617 0 : sig2(i,j) = (p5*(stressp(i,j) &
1618 : - sqrt(stressm(i,j)**2+c4*stress12(i,j)**2))) & ! LCOV_EXCL_LINE
1619 22166 : / strength(i,j)
1620 : else
1621 248892 : sig1(i,j) = spval_dbl
1622 248892 : sig2(i,j) = spval_dbl
1623 248892 : sigP(i,j) = spval_dbl
1624 : endif
1625 : enddo
1626 : enddo
1627 :
1628 636 : end subroutine principal_stress
1629 :
1630 : !=======================================================================
1631 : ! Compute deformations for mechanical redistribution
1632 : !
1633 : ! author: Elizabeth C. Hunke, LANL
1634 : !
1635 : ! 2019: subroutine created by Philippe Blain, ECCC
1636 :
1637 22944 : subroutine deformations (nx_block, ny_block, &
1638 : icellT, & ! LCOV_EXCL_LINE
1639 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1640 : uvel, vvel, & ! LCOV_EXCL_LINE
1641 : dxT, dyT, & ! LCOV_EXCL_LINE
1642 : cxp, cyp, & ! LCOV_EXCL_LINE
1643 : cxm, cym, & ! LCOV_EXCL_LINE
1644 : tarear, & ! LCOV_EXCL_LINE
1645 : shear, divu, & ! LCOV_EXCL_LINE
1646 22944 : rdg_conv, rdg_shear )
1647 :
1648 : use ice_constants, only: p25, p5
1649 :
1650 : integer (kind=int_kind), intent(in) :: &
1651 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1652 : icellT ! no. of cells where iceTmask = .true.
1653 :
1654 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1655 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1656 : indxTj ! compressed index in j-direction
1657 :
1658 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1659 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1660 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1661 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1662 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1663 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1664 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1665 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1666 : cxm , & ! 0.5*HTN - 1.5*HTS ! LCOV_EXCL_LINE
1667 : tarear ! 1/tarea
1668 :
1669 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1670 : shear , & ! strain rate II component (1/s) ! LCOV_EXCL_LINE
1671 : divu , & ! strain rate I component, velocity divergence (1/s) ! LCOV_EXCL_LINE
1672 : rdg_conv , & ! convergence term for ridging (1/s) ! LCOV_EXCL_LINE
1673 : rdg_shear ! shear term for ridging (1/s)
1674 :
1675 : ! local variables
1676 :
1677 : integer (kind=int_kind) :: &
1678 : i, j, ij
1679 :
1680 : real (kind=dbl_kind) :: & ! at each corner :
1681 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1682 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1683 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1684 : Deltane, Deltanw, Deltase, Deltasw , & ! Delta ! LCOV_EXCL_LINE
1685 5760 : tmp ! useful combination
1686 :
1687 : character(len=*), parameter :: subname = '(deformations)'
1688 :
1689 7161231 : do ij = 1, icellT
1690 7138287 : i = indxTi(ij)
1691 7138287 : j = indxTj(ij)
1692 :
1693 : !-----------------------------------------------------------------
1694 : ! strain rates
1695 : ! NOTE these are actually strain rates * area (m^2/s)
1696 : !-----------------------------------------------------------------
1697 : call strain_rates (nx_block, ny_block, &
1698 : i, j, & ! LCOV_EXCL_LINE
1699 : uvel, vvel, & ! LCOV_EXCL_LINE
1700 : dxT, dyT, & ! LCOV_EXCL_LINE
1701 : cxp, cyp, & ! LCOV_EXCL_LINE
1702 : cxm, cym, & ! LCOV_EXCL_LINE
1703 : divune, divunw, & ! LCOV_EXCL_LINE
1704 : divuse, divusw, & ! LCOV_EXCL_LINE
1705 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1706 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1707 : shearne, shearnw, & ! LCOV_EXCL_LINE
1708 : shearse, shearsw, & ! LCOV_EXCL_LINE
1709 : Deltane, Deltanw, & ! LCOV_EXCL_LINE
1710 7138287 : Deltase, Deltasw )
1711 : !-----------------------------------------------------------------
1712 : ! deformations for mechanical redistribution
1713 : !-----------------------------------------------------------------
1714 7138287 : divu(i,j) = p25*(divune + divunw + divuse + divusw) * tarear(i,j)
1715 7138287 : tmp = p25*(Deltane + Deltanw + Deltase + Deltasw) * tarear(i,j)
1716 7138287 : rdg_conv(i,j) = -min(divu(i,j),c0)
1717 7138287 : rdg_shear(i,j) = p5*(tmp-abs(divu(i,j)))
1718 :
1719 : ! diagnostic only
1720 : ! shear = sqrt(tension**2 + shearing**2)
1721 624817 : shear(i,j) = p25*tarear(i,j)*sqrt( &
1722 : (tensionne + tensionnw + tensionse + tensionsw)**2 + & ! LCOV_EXCL_LINE
1723 7786048 : (shearne + shearnw + shearse + shearsw )**2)
1724 :
1725 : enddo ! ij
1726 :
1727 22944 : end subroutine deformations
1728 :
1729 : !=======================================================================
1730 : ! Compute deformations for mechanical redistribution at T point
1731 : !
1732 : ! author: JF Lemieux, ECCC
1733 : ! Nov 2021
1734 :
1735 0 : subroutine deformationsCD_T (nx_block, ny_block, &
1736 : icellT, & ! LCOV_EXCL_LINE
1737 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1738 : uvelE, vvelE, & ! LCOV_EXCL_LINE
1739 : uvelN, vvelN, & ! LCOV_EXCL_LINE
1740 : dxN, dyE, & ! LCOV_EXCL_LINE
1741 : dxT, dyT, & ! LCOV_EXCL_LINE
1742 : tarear, & ! LCOV_EXCL_LINE
1743 : shear, divu, & ! LCOV_EXCL_LINE
1744 0 : rdg_conv, rdg_shear )
1745 :
1746 : use ice_constants, only: p5
1747 :
1748 : integer (kind=int_kind), intent(in) :: &
1749 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1750 : icellT ! no. of cells where iceTmask = .true.
1751 :
1752 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1753 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1754 : indxTj ! compressed index in j-direction
1755 :
1756 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1757 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1758 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1759 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1760 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1761 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
1762 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
1763 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1764 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1765 : tarear ! 1/tarea
1766 :
1767 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1768 : shear , & ! strain rate II component (1/s) ! LCOV_EXCL_LINE
1769 : divu , & ! strain rate I component, velocity divergence (1/s) ! LCOV_EXCL_LINE
1770 : rdg_conv , & ! convergence term for ridging (1/s) ! LCOV_EXCL_LINE
1771 : rdg_shear ! shear term for ridging (1/s)
1772 :
1773 : ! local variables
1774 :
1775 : integer (kind=int_kind) :: &
1776 : i, j, ij
1777 :
1778 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1779 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
1780 : tensionT , & ! tension at T point ! LCOV_EXCL_LINE
1781 : shearT , & ! shear at T point ! LCOV_EXCL_LINE
1782 0 : DeltaT ! delt at T point
1783 :
1784 : real (kind=dbl_kind) :: &
1785 0 : tmp ! useful combination
1786 :
1787 : character(len=*), parameter :: subname = '(deformations_T)'
1788 :
1789 : !-----------------------------------------------------------------
1790 : ! strain rates
1791 : ! NOTE these are actually strain rates * area (m^2/s)
1792 : !-----------------------------------------------------------------
1793 :
1794 : call strain_rates_T (nx_block , ny_block , &
1795 : icellT , & ! LCOV_EXCL_LINE
1796 : indxTi(:) , indxTj (:) , & ! LCOV_EXCL_LINE
1797 : uvelE (:,:), vvelE (:,:), & ! LCOV_EXCL_LINE
1798 : uvelN (:,:), vvelN (:,:), & ! LCOV_EXCL_LINE
1799 : dxN (:,:), dyE (:,:), & ! LCOV_EXCL_LINE
1800 : dxT (:,:), dyT (:,:), & ! LCOV_EXCL_LINE
1801 : divT (:,:), tensionT(:,:), & ! LCOV_EXCL_LINE
1802 0 : shearT(:,:), DeltaT (:,:) )
1803 :
1804 0 : do ij = 1, icellT
1805 0 : i = indxTi(ij)
1806 0 : j = indxTj(ij)
1807 :
1808 : !-----------------------------------------------------------------
1809 : ! deformations for mechanical redistribution
1810 : !-----------------------------------------------------------------
1811 0 : divu(i,j) = divT(i,j) * tarear(i,j)
1812 0 : tmp = Deltat(i,j) * tarear(i,j)
1813 0 : rdg_conv(i,j) = -min(divu(i,j),c0)
1814 0 : rdg_shear(i,j) = p5*(tmp-abs(divu(i,j)))
1815 :
1816 : ! diagnostic only
1817 : ! shear = sqrt(tension**2 + shearing**2)
1818 0 : shear(i,j) = tarear(i,j)*sqrt( tensionT(i,j)**2 + shearT(i,j)**2 )
1819 :
1820 : enddo ! ij
1821 :
1822 0 : end subroutine deformationsCD_T
1823 :
1824 :
1825 : !=======================================================================
1826 : ! Compute deformations for mechanical redistribution at T point
1827 : !
1828 : ! author: JF Lemieux, ECCC
1829 : ! Nov 2021
1830 :
1831 0 : subroutine deformationsC_T (nx_block, ny_block, &
1832 : icellT, & ! LCOV_EXCL_LINE
1833 : indxTi, indxTj, & ! LCOV_EXCL_LINE
1834 : uvelE, vvelE, & ! LCOV_EXCL_LINE
1835 : uvelN, vvelN, & ! LCOV_EXCL_LINE
1836 : dxN, dyE, & ! LCOV_EXCL_LINE
1837 : dxT, dyT, & ! LCOV_EXCL_LINE
1838 : tarear, uarea, & ! LCOV_EXCL_LINE
1839 : shearU, & ! LCOV_EXCL_LINE
1840 : shear, divu, & ! LCOV_EXCL_LINE
1841 0 : rdg_conv, rdg_shear )
1842 :
1843 : use ice_constants, only: p5
1844 :
1845 : integer (kind=int_kind), intent(in) :: &
1846 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1847 : icellT ! no. of cells where iceTmask = .true.
1848 :
1849 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
1850 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
1851 : indxTj ! compressed index in j-direction
1852 :
1853 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1854 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1855 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1856 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
1857 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
1858 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
1859 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
1860 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1861 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1862 : tarear , & ! 1/tarea ! LCOV_EXCL_LINE
1863 : uarea , & ! area of u cell ! LCOV_EXCL_LINE
1864 : shearU ! shearU
1865 :
1866 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(inout) :: &
1867 : shear , & ! strain rate II component (1/s) ! LCOV_EXCL_LINE
1868 : divu , & ! strain rate I component, velocity divergence (1/s) ! LCOV_EXCL_LINE
1869 : rdg_conv , & ! convergence term for ridging (1/s) ! LCOV_EXCL_LINE
1870 : rdg_shear ! shear term for ridging (1/s)
1871 :
1872 : ! local variables
1873 :
1874 : integer (kind=int_kind) :: &
1875 : i, j, ij
1876 :
1877 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
1878 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
1879 : tensionT , & ! tension at T point ! LCOV_EXCL_LINE
1880 : shearT , & ! shear at T point ! LCOV_EXCL_LINE
1881 0 : DeltaT ! delt at T point
1882 :
1883 : real (kind=dbl_kind) :: &
1884 : tmp , & ! useful combination ! LCOV_EXCL_LINE
1885 0 : shearTsqr ! strain rates squared at T point
1886 :
1887 : character(len=*), parameter :: subname = '(deformations_T2)'
1888 :
1889 : !-----------------------------------------------------------------
1890 : ! strain rates
1891 : ! NOTE these are actually strain rates * area (m^2/s)
1892 : !-----------------------------------------------------------------
1893 :
1894 : call strain_rates_T (nx_block , ny_block , &
1895 : icellT , & ! LCOV_EXCL_LINE
1896 : indxTi(:) , indxTj (:) , & ! LCOV_EXCL_LINE
1897 : uvelE (:,:), vvelE (:,:), & ! LCOV_EXCL_LINE
1898 : uvelN (:,:), vvelN (:,:), & ! LCOV_EXCL_LINE
1899 : dxN (:,:), dyE (:,:), & ! LCOV_EXCL_LINE
1900 : dxT (:,:), dyT (:,:), & ! LCOV_EXCL_LINE
1901 : divT (:,:), tensionT(:,:), & ! LCOV_EXCL_LINE
1902 0 : shearT(:,:), DeltaT (:,:) )
1903 :
1904 : ! DeltaT is calc by strain_rates_T but replaced by calculation below.
1905 :
1906 0 : do ij = 1, icellT
1907 0 : i = indxTi(ij)
1908 0 : j = indxTj(ij)
1909 :
1910 : !-----------------------------------------------------------------
1911 : ! deformations for mechanical redistribution
1912 : !-----------------------------------------------------------------
1913 :
1914 0 : shearTsqr = (shearU(i ,j )**2 * uarea(i ,j ) &
1915 : + shearU(i ,j-1)**2 * uarea(i ,j-1) & ! LCOV_EXCL_LINE
1916 : + shearU(i-1,j-1)**2 * uarea(i-1,j-1) & ! LCOV_EXCL_LINE
1917 : + shearU(i-1,j )**2 * uarea(i-1,j )) & ! LCOV_EXCL_LINE
1918 0 : / (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))
1919 :
1920 0 : DeltaT(i,j) = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr))
1921 :
1922 0 : divu(i,j) = divT(i,j) * tarear(i,j)
1923 0 : tmp = DeltaT(i,j) * tarear(i,j)
1924 0 : rdg_conv(i,j) = -min(divu(i,j),c0)
1925 0 : rdg_shear(i,j) = p5*(tmp-abs(divu(i,j)))
1926 :
1927 : ! diagnostic only...maybe we dont want to use shearTsqr here????
1928 : ! shear = sqrt(tension**2 + shearing**2)
1929 0 : shear(i,j) = tarear(i,j)*sqrt( tensionT(i,j)**2 + shearT(i,j)**2 )
1930 :
1931 : enddo ! ij
1932 :
1933 0 : end subroutine deformationsC_T
1934 :
1935 : !=======================================================================
1936 : ! Compute strain rates
1937 : !
1938 : ! author: Elizabeth C. Hunke, LANL
1939 : !
1940 : ! 2019: subroutine created by Philippe Blain, ECCC
1941 :
1942 1720327167 : subroutine strain_rates (nx_block, ny_block, &
1943 : i, j, & ! LCOV_EXCL_LINE
1944 : uvel, vvel, & ! LCOV_EXCL_LINE
1945 : dxT, dyT, & ! LCOV_EXCL_LINE
1946 : cxp, cyp, & ! LCOV_EXCL_LINE
1947 : cxm, cym, & ! LCOV_EXCL_LINE
1948 : divune, divunw, & ! LCOV_EXCL_LINE
1949 : divuse, divusw, & ! LCOV_EXCL_LINE
1950 : tensionne, tensionnw, & ! LCOV_EXCL_LINE
1951 : tensionse, tensionsw, & ! LCOV_EXCL_LINE
1952 : shearne, shearnw, & ! LCOV_EXCL_LINE
1953 : shearse, shearsw, & ! LCOV_EXCL_LINE
1954 : Deltane, Deltanw, & ! LCOV_EXCL_LINE
1955 : Deltase, Deltasw )
1956 :
1957 : integer (kind=int_kind), intent(in) :: &
1958 : nx_block, ny_block ! block dimensions
1959 :
1960 : integer (kind=int_kind), intent(in) :: &
1961 : i, j ! indices
1962 :
1963 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
1964 : uvel , & ! x-component of velocity (m/s) ! LCOV_EXCL_LINE
1965 : vvel , & ! y-component of velocity (m/s) ! LCOV_EXCL_LINE
1966 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
1967 : dyT , & ! height of T-cell through the middle (m) ! LCOV_EXCL_LINE
1968 : cyp , & ! 1.5*HTE - 0.5*HTW ! LCOV_EXCL_LINE
1969 : cxp , & ! 1.5*HTN - 0.5*HTS ! LCOV_EXCL_LINE
1970 : cym , & ! 0.5*HTE - 1.5*HTW ! LCOV_EXCL_LINE
1971 : cxm ! 0.5*HTN - 1.5*HTS
1972 :
1973 : real (kind=dbl_kind), intent(out):: & ! at each corner :
1974 : divune, divunw, divuse, divusw , & ! divergence ! LCOV_EXCL_LINE
1975 : tensionne, tensionnw, tensionse, tensionsw, & ! tension ! LCOV_EXCL_LINE
1976 : shearne, shearnw, shearse, shearsw , & ! shearing ! LCOV_EXCL_LINE
1977 : Deltane, Deltanw, Deltase, Deltasw ! Delta
1978 :
1979 : character(len=*), parameter :: subname = '(strain_rates)'
1980 :
1981 : !-----------------------------------------------------------------
1982 : ! strain rates
1983 : ! NOTE these are actually strain rates * area (m^2/s)
1984 : !-----------------------------------------------------------------
1985 :
1986 : ! divergence = e_11 + e_22
1987 1204647176 : divune = cyp(i,j)*uvel(i ,j ) - dyT(i,j)*uvel(i-1,j ) &
1988 2924974343 : + cxp(i,j)*vvel(i ,j ) - dxT(i,j)*vvel(i ,j-1)
1989 1204647176 : divunw = cym(i,j)*uvel(i-1,j ) + dyT(i,j)*uvel(i ,j ) &
1990 2924974343 : + cxp(i,j)*vvel(i-1,j ) - dxT(i,j)*vvel(i-1,j-1)
1991 1204647176 : divusw = cym(i,j)*uvel(i-1,j-1) + dyT(i,j)*uvel(i ,j-1) &
1992 2924974343 : + cxm(i,j)*vvel(i-1,j-1) + dxT(i,j)*vvel(i-1,j )
1993 1204647176 : divuse = cyp(i,j)*uvel(i ,j-1) - dyT(i,j)*uvel(i-1,j-1) &
1994 2924974343 : + cxm(i,j)*vvel(i ,j-1) + dxT(i,j)*vvel(i ,j )
1995 :
1996 : ! tension strain rate = e_11 - e_22
1997 1204647176 : tensionne = -cym(i,j)*uvel(i ,j ) - dyT(i,j)*uvel(i-1,j ) &
1998 2924974343 : + cxm(i,j)*vvel(i ,j ) + dxT(i,j)*vvel(i ,j-1)
1999 1204647176 : tensionnw = -cyp(i,j)*uvel(i-1,j ) + dyT(i,j)*uvel(i ,j ) &
2000 2924974343 : + cxm(i,j)*vvel(i-1,j ) + dxT(i,j)*vvel(i-1,j-1)
2001 1204647176 : tensionsw = -cyp(i,j)*uvel(i-1,j-1) + dyT(i,j)*uvel(i ,j-1) &
2002 2924974343 : + cxp(i,j)*vvel(i-1,j-1) - dxT(i,j)*vvel(i-1,j )
2003 1204647176 : tensionse = -cym(i,j)*uvel(i ,j-1) - dyT(i,j)*uvel(i-1,j-1) &
2004 2924974343 : + cxp(i,j)*vvel(i ,j-1) - dxT(i,j)*vvel(i ,j )
2005 :
2006 : ! shearing strain rate = 2*e_12
2007 1204647176 : shearne = -cym(i,j)*vvel(i ,j ) - dyT(i,j)*vvel(i-1,j ) &
2008 2924974343 : - cxm(i,j)*uvel(i ,j ) - dxT(i,j)*uvel(i ,j-1)
2009 1204647176 : shearnw = -cyp(i,j)*vvel(i-1,j ) + dyT(i,j)*vvel(i ,j ) &
2010 2924974343 : - cxm(i,j)*uvel(i-1,j ) - dxT(i,j)*uvel(i-1,j-1)
2011 1204647176 : shearsw = -cyp(i,j)*vvel(i-1,j-1) + dyT(i,j)*vvel(i ,j-1) &
2012 2924974343 : - cxp(i,j)*uvel(i-1,j-1) + dxT(i,j)*uvel(i-1,j )
2013 1204647176 : shearse = -cym(i,j)*vvel(i ,j-1) - dyT(i,j)*vvel(i-1,j-1) &
2014 2924974343 : - cxp(i,j)*uvel(i ,j-1) + dxT(i,j)*uvel(i ,j )
2015 :
2016 : ! Delta (in the denominator of zeta, eta)
2017 1720327167 : Deltane = sqrt(divune**2 + e_factor*(tensionne**2 + shearne**2))
2018 1720327167 : Deltanw = sqrt(divunw**2 + e_factor*(tensionnw**2 + shearnw**2))
2019 1720327167 : Deltasw = sqrt(divusw**2 + e_factor*(tensionsw**2 + shearsw**2))
2020 1720327167 : Deltase = sqrt(divuse**2 + e_factor*(tensionse**2 + shearse**2))
2021 :
2022 1720327167 : end subroutine strain_rates
2023 :
2024 : !=======================================================================
2025 : ! Compute dtsd (div, tension, shear, delta) strain rates at the T point
2026 : !
2027 : ! author: JF Lemieux, ECCC
2028 : ! Nov 2021
2029 :
2030 0 : subroutine strain_rates_Tdtsd (nx_block, ny_block, &
2031 : icellT, & ! LCOV_EXCL_LINE
2032 : indxTi, indxTj, & ! LCOV_EXCL_LINE
2033 : uvelE, vvelE, & ! LCOV_EXCL_LINE
2034 : uvelN, vvelN, & ! LCOV_EXCL_LINE
2035 : dxN, dyE, & ! LCOV_EXCL_LINE
2036 : dxT, dyT, & ! LCOV_EXCL_LINE
2037 : divT, tensionT, & ! LCOV_EXCL_LINE
2038 0 : shearT, DeltaT )
2039 :
2040 : integer (kind=int_kind), intent(in) :: &
2041 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2042 : icellT
2043 :
2044 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2045 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2046 : indxTj ! compressed index in j-direction
2047 :
2048 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2049 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2050 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2051 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2052 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2053 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
2054 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
2055 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
2056 : dyT ! height of T-cell through the middle (m)
2057 :
2058 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
2059 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
2060 : tensionT , & ! tension at T point ! LCOV_EXCL_LINE
2061 : shearT , & ! shear at T point ! LCOV_EXCL_LINE
2062 : DeltaT ! strain rates at the T point
2063 :
2064 : ! local variables
2065 :
2066 : integer (kind=int_kind) :: &
2067 : ij, i, j ! indices
2068 :
2069 : character(len=*), parameter :: subname = '(strain_rates_Tdtsd)'
2070 :
2071 : !-----------------------------------------------------------------
2072 : ! strain rates
2073 : ! NOTE these are actually strain rates * area (m^2/s)
2074 : !-----------------------------------------------------------------
2075 :
2076 : ! compute divT, tensionT
2077 : call strain_rates_Tdt (nx_block, ny_block, &
2078 : icellT, & ! LCOV_EXCL_LINE
2079 : indxTi, indxTj, & ! LCOV_EXCL_LINE
2080 : uvelE, vvelE, & ! LCOV_EXCL_LINE
2081 : uvelN, vvelN, & ! LCOV_EXCL_LINE
2082 : dxN, dyE, & ! LCOV_EXCL_LINE
2083 : dxT, dyT, & ! LCOV_EXCL_LINE
2084 0 : divT, tensionT )
2085 :
2086 0 : shearT (:,:) = c0
2087 0 : deltaT (:,:) = c0
2088 :
2089 0 : do ij = 1, icellT
2090 0 : i = indxTi(ij)
2091 0 : j = indxTj(ij)
2092 :
2093 : ! shearing strain rate = 2*e_12
2094 0 : shearT(i,j) = (dxT(i,j)**2)*(uvelN(i,j)/dxN(i,j) - uvelN(i,j-1)/dxN(i,j-1)) &
2095 0 : + (dyT(i,j)**2)*(vvelE(i,j)/dyE(i,j) - vvelE(i-1,j)/dyE(i-1,j))
2096 :
2097 : ! Delta (in the denominator of zeta, eta)
2098 0 : DeltaT(i,j) = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearT(i,j)**2))
2099 :
2100 : enddo
2101 :
2102 0 : end subroutine strain_rates_Tdtsd
2103 :
2104 : !=======================================================================
2105 : ! Compute the dt (div, tension) strain rates at the T point
2106 : !
2107 : ! author: JF Lemieux, ECCC
2108 : ! Nov 2021
2109 :
2110 0 : subroutine strain_rates_Tdt (nx_block, ny_block, &
2111 : icellT, & ! LCOV_EXCL_LINE
2112 : indxTi, indxTj, & ! LCOV_EXCL_LINE
2113 : uvelE, vvelE, & ! LCOV_EXCL_LINE
2114 : uvelN, vvelN, & ! LCOV_EXCL_LINE
2115 : dxN, dyE, & ! LCOV_EXCL_LINE
2116 : dxT, dyT, & ! LCOV_EXCL_LINE
2117 0 : divT, tensionT )
2118 :
2119 : integer (kind=int_kind), intent(in) :: &
2120 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2121 : icellT
2122 :
2123 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2124 : indxTi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2125 : indxTj ! compressed index in j-direction
2126 :
2127 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2128 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2129 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2130 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2131 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2132 : dxN , & ! width of N-cell through the middle (m) ! LCOV_EXCL_LINE
2133 : dyE , & ! height of E-cell through the middle (m) ! LCOV_EXCL_LINE
2134 : dxT , & ! width of T-cell through the middle (m) ! LCOV_EXCL_LINE
2135 : dyT ! height of T-cell through the middle (m)
2136 :
2137 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
2138 : divT , & ! divergence at T point ! LCOV_EXCL_LINE
2139 : tensionT ! tension at T point
2140 :
2141 : ! local variables
2142 :
2143 : integer (kind=int_kind) :: &
2144 : ij, i, j ! indices
2145 :
2146 : character(len=*), parameter :: subname = '(strain_rates_Tdt)'
2147 :
2148 : !-----------------------------------------------------------------
2149 : ! strain rates
2150 : ! NOTE these are actually strain rates * area (m^2/s)
2151 : !-----------------------------------------------------------------
2152 :
2153 0 : divT (:,:) = c0
2154 0 : tensionT(:,:) = c0
2155 :
2156 0 : do ij = 1, icellT
2157 0 : i = indxTi(ij)
2158 0 : j = indxTj(ij)
2159 :
2160 : ! divergence = e_11 + e_22
2161 0 : divT (i,j)= dyE(i,j)*uvelE(i ,j ) - dyE(i-1,j)*uvelE(i-1,j ) &
2162 0 : + dxN(i,j)*vvelN(i ,j ) - dxN(i,j-1)*vvelN(i ,j-1)
2163 :
2164 : ! tension strain rate = e_11 - e_22
2165 0 : tensionT(i,j) = (dyT(i,j)**2)*(uvelE(i,j)/dyE(i,j) - uvelE(i-1,j)/dyE(i-1,j)) &
2166 0 : - (dxT(i,j)**2)*(vvelN(i,j)/dxN(i,j) - vvelN(i,j-1)/dxN(i,j-1))
2167 :
2168 : enddo
2169 :
2170 0 : end subroutine strain_rates_Tdt
2171 :
2172 : !=======================================================================
2173 : ! Compute strain rates at the U point including boundary conditions
2174 : !
2175 : ! author: JF Lemieux, ECCC
2176 : ! Nov 2021
2177 :
2178 0 : subroutine strain_rates_U (nx_block, ny_block, &
2179 : icellU, & ! LCOV_EXCL_LINE
2180 : indxUi, indxUj, & ! LCOV_EXCL_LINE
2181 : uvelE, vvelE, & ! LCOV_EXCL_LINE
2182 : uvelN, vvelN, & ! LCOV_EXCL_LINE
2183 : uvelU, vvelU, & ! LCOV_EXCL_LINE
2184 : dxE, dyN, & ! LCOV_EXCL_LINE
2185 : dxU, dyU, & ! LCOV_EXCL_LINE
2186 : ratiodxN, ratiodxNr, & ! LCOV_EXCL_LINE
2187 : ratiodyE, ratiodyEr, & ! LCOV_EXCL_LINE
2188 : epm, npm, & ! LCOV_EXCL_LINE
2189 : divergU, tensionU, & ! LCOV_EXCL_LINE
2190 0 : shearU, DeltaU )
2191 :
2192 : integer (kind=int_kind), intent(in) :: &
2193 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2194 : icellU
2195 :
2196 : integer (kind=int_kind), dimension (nx_block*ny_block), intent(in) :: &
2197 : indxUi , & ! compressed index in i-direction ! LCOV_EXCL_LINE
2198 : indxUj ! compressed index in j-direction
2199 :
2200 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(in) :: &
2201 : uvelE , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2202 : vvelE , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2203 : uvelN , & ! x-component of velocity (m/s) at the E point ! LCOV_EXCL_LINE
2204 : vvelN , & ! y-component of velocity (m/s) at the N point ! LCOV_EXCL_LINE
2205 : uvelU , & ! x-component of velocity (m/s) interp. at U point ! LCOV_EXCL_LINE
2206 : vvelU , & ! y-component of velocity (m/s) interp. at U point ! LCOV_EXCL_LINE
2207 : dxE , & ! width of E-cell through the middle (m) ! LCOV_EXCL_LINE
2208 : dyN , & ! height of N-cell through the middle (m) ! LCOV_EXCL_LINE
2209 : dxU , & ! width of U-cell through the middle (m) ! LCOV_EXCL_LINE
2210 : dyU , & ! height of U-cell through the middle (m) ! LCOV_EXCL_LINE
2211 : ratiodxN , & ! -dxN(i+1,j)/dxN(i,j) for BCs ! LCOV_EXCL_LINE
2212 : ratiodxNr, & ! -dxN(i,j)/dxN(i+1,j) for BCs ! LCOV_EXCL_LINE
2213 : ratiodyE , & ! -dyE(i,j+1)/dyE(i,j) for BCs ! LCOV_EXCL_LINE
2214 : ratiodyEr, & ! -dyE(i,j)/dyE(i,j+1) for BCs ! LCOV_EXCL_LINE
2215 : epm , & ! E-cell mask ! LCOV_EXCL_LINE
2216 : npm ! N-cell mask
2217 :
2218 : real (kind=dbl_kind), dimension (nx_block,ny_block), intent(out):: &
2219 : divergU , & ! divergence at U point ! LCOV_EXCL_LINE
2220 : tensionU , & ! tension at U point ! LCOV_EXCL_LINE
2221 : shearU , & ! shear at U point ! LCOV_EXCL_LINE
2222 : DeltaU ! delt at the U point
2223 :
2224 : ! local variables
2225 :
2226 : integer (kind=int_kind) :: &
2227 : ij, i, j ! indices
2228 :
2229 : real (kind=dbl_kind) :: &
2230 0 : uNip1j, uNij, vEijp1, vEij, uEijp1, uEij, vNip1j, vNij
2231 :
2232 : character(len=*), parameter :: subname = '(strain_rates_U)'
2233 :
2234 : !-----------------------------------------------------------------
2235 : ! strain rates
2236 : ! NOTE these are actually strain rates * area (m^2/s)
2237 : !-----------------------------------------------------------------
2238 :
2239 0 : divergU (:,:) = c0
2240 0 : tensionU(:,:) = c0
2241 0 : shearU (:,:) = c0
2242 0 : deltaU (:,:) = c0
2243 :
2244 0 : do ij = 1, icellU
2245 0 : i = indxUi(ij)
2246 0 : j = indxUj(ij)
2247 :
2248 0 : uNip1j = uvelN(i+1,j) * npm(i+1,j) &
2249 0 : +(npm(i,j)-npm(i+1,j)) * npm(i,j) * ratiodxN(i,j) * uvelN(i,j)
2250 0 : uNij = uvelN(i,j) * npm(i,j) &
2251 0 : +(npm(i+1,j)-npm(i,j)) * npm(i+1,j) * ratiodxNr(i,j) * uvelN(i+1,j)
2252 0 : vEijp1 = vvelE(i,j+1) * epm(i,j+1) &
2253 0 : +(epm(i,j)-epm(i,j+1)) * epm(i,j) * ratiodyE(i,j) * vvelE(i,j)
2254 0 : vEij = vvelE(i,j) * epm(i,j) &
2255 0 : +(epm(i,j+1)-epm(i,j)) * epm(i,j+1) * ratiodyEr(i,j) * vvelE(i,j+1)
2256 :
2257 : ! divergence = e_11 + e_22
2258 0 : divergU (i,j) = dyU(i,j) * ( uNip1j - uNij ) &
2259 : + uvelU(i,j) * ( dyN(i+1,j) - dyN(i,j) ) & ! LCOV_EXCL_LINE
2260 : + dxU(i,j) * ( vEijp1 - vEij ) & ! LCOV_EXCL_LINE
2261 0 : + vvelU(i,j) * ( dxE(i,j+1) - dxE(i,j) )
2262 :
2263 : ! tension strain rate = e_11 - e_22
2264 0 : tensionU(i,j) = dyU(i,j) * ( uNip1j - uNij ) &
2265 : - uvelU(i,j) * ( dyN(i+1,j) - dyN(i,j) ) & ! LCOV_EXCL_LINE
2266 : - dxU(i,j) * ( vEijp1 - vEij ) & ! LCOV_EXCL_LINE
2267 0 : + vvelU(i,j) * ( dxE(i,j+1) - dxE(i,j) )
2268 :
2269 0 : uEijp1 = uvelE(i,j+1) * epm(i,j+1) &
2270 0 : +(epm(i,j)-epm(i,j+1)) * epm(i,j) * ratiodyE(i,j) * uvelE(i,j)
2271 0 : uEij = uvelE(i,j) * epm(i,j) &
2272 0 : +(epm(i,j+1)-epm(i,j)) * epm(i,j+1) * ratiodyEr(i,j) * uvelE(i,j+1)
2273 0 : vNip1j = vvelN(i+1,j) * npm(i+1,j) &
2274 0 : +(npm(i,j)-npm(i+1,j)) * npm(i,j) * ratiodxN(i,j) * vvelN(i,j)
2275 0 : vNij = vvelN(i,j) * npm(i,j) &
2276 0 : +(npm(i+1,j)-npm(i,j)) * npm(i+1,j) * ratiodxNr(i,j) * vvelN(i+1,j)
2277 :
2278 : ! shearing strain rate = 2*e_12
2279 0 : shearU(i,j) = dxU(i,j) * ( uEijp1 - uEij ) &
2280 : - uvelU(i,j) * ( dxE(i,j+1) - dxE(i,j) ) & ! LCOV_EXCL_LINE
2281 : + dyU(i,j) * ( vNip1j - vNij ) & ! LCOV_EXCL_LINE
2282 0 : - vvelU(i,j) * ( dyN(i+1,j) - dyN(i,j) )
2283 :
2284 : ! Delta (in the denominator of zeta, eta)
2285 0 : DeltaU(i,j) = sqrt(divergU(i,j)**2 + e_factor*(tensionU(i,j)**2 + shearU(i,j)**2))
2286 :
2287 : enddo
2288 :
2289 0 : end subroutine strain_rates_U
2290 :
2291 : !=======================================================================
2292 : ! Computes viscosities and replacement pressure for stress
2293 : ! calculations. Note that tensile strength is included here.
2294 : !
2295 : ! Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. J. Phys.
2296 : ! Oceanogr., 9, 817-846.
2297 : !
2298 : ! Konig Beatty, C. and Holland, D. M. (2010). Modeling landfast ice by
2299 : ! adding tensile strength. J. Phys. Oceanogr. 40, 185-198.
2300 : !
2301 : ! Lemieux, J. F. et al. (2016). Improving the simulation of landfast ice
2302 : ! by combining tensile strength and a parameterization for grounded ridges.
2303 : ! J. Geophys. Res. Oceans, 121, 7354-7368.
2304 :
2305 6852755520 : subroutine visc_replpress(strength, DminArea, Delta, &
2306 : zetax2, etax2, rep_prs, capping)
2307 :
2308 : real (kind=dbl_kind), intent(in):: &
2309 : strength, & ! ! LCOV_EXCL_LINE
2310 : DminArea !
2311 :
2312 : real (kind=dbl_kind), intent(in):: &
2313 : Delta , & ! ! LCOV_EXCL_LINE
2314 : capping !
2315 :
2316 : real (kind=dbl_kind), intent(out):: &
2317 : zetax2 , & ! bulk viscosity ! LCOV_EXCL_LINE
2318 : etax2 , & ! shear viscosity ! LCOV_EXCL_LINE
2319 : rep_prs ! replacement pressure
2320 :
2321 : ! local variables
2322 : real (kind=dbl_kind) :: &
2323 599824320 : tmpcalc ! temporary
2324 :
2325 : character(len=*), parameter :: subname = '(visc_replpress)'
2326 :
2327 : ! NOTE: for comp. efficiency 2 x zeta and 2 x eta are used in the code
2328 :
2329 : tmpcalc = capping *(strength/max(Delta,DminArea))+ &
2330 6852755520 : (c1-capping)*(strength/(Delta + DminArea))
2331 6852755520 : zetax2 = (c1+Ktens)*tmpcalc
2332 6852755520 : rep_prs = (c1-Ktens)*tmpcalc*Delta
2333 6852755520 : etax2 = epp2i*zetax2
2334 :
2335 6852755520 : end subroutine visc_replpress
2336 :
2337 : !=======================================================================
2338 : ! Do a halo update on 1 field
2339 :
2340 0 : subroutine dyn_haloUpdate1(halo_info, halo_info_mask, field_loc, field_type, fld1)
2341 :
2342 : use ice_boundary, only: ice_halo, ice_HaloUpdate
2343 : use ice_domain, only: maskhalo_dyn, halo_dynbundle
2344 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2345 :
2346 : type (ice_halo), intent(in) :: &
2347 : halo_info , & ! standard unmasked halo ! LCOV_EXCL_LINE
2348 : halo_info_mask ! masked halo
2349 :
2350 : integer (kind=int_kind), intent(in) :: &
2351 : field_loc , & ! field loc ! LCOV_EXCL_LINE
2352 : field_type ! field_type
2353 :
2354 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
2355 : fld1 ! fields to halo
2356 :
2357 : ! local variables
2358 :
2359 : character(len=*), parameter :: subname = '(dyn_haloUpdate1)'
2360 :
2361 0 : call ice_timer_start(timer_bound)
2362 :
2363 0 : if (maskhalo_dyn) then
2364 0 : call ice_HaloUpdate (fld1 , halo_info_mask, &
2365 0 : field_loc, field_type)
2366 : else
2367 0 : call ice_HaloUpdate (fld1 , halo_info , &
2368 0 : field_loc, field_type)
2369 : endif
2370 :
2371 0 : call ice_timer_stop(timer_bound)
2372 :
2373 0 : end subroutine dyn_haloUpdate1
2374 :
2375 : !=======================================================================
2376 : ! Do a halo update on 2 fields
2377 :
2378 1388160 : subroutine dyn_haloUpdate2(halo_info, halo_info_mask, field_loc, field_type, fld1, fld2)
2379 :
2380 : use ice_boundary, only: ice_halo, ice_HaloUpdate
2381 : use ice_domain, only: maskhalo_dyn, halo_dynbundle
2382 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2383 :
2384 : type (ice_halo), intent(in) :: &
2385 : halo_info , & ! standard unmasked halo ! LCOV_EXCL_LINE
2386 : halo_info_mask ! masked halo
2387 :
2388 : integer (kind=int_kind), intent(in) :: &
2389 : field_loc , & ! field loc ! LCOV_EXCL_LINE
2390 : field_type ! field_type
2391 :
2392 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
2393 : fld1 , & ! fields to halo ! LCOV_EXCL_LINE
2394 : fld2 !
2395 :
2396 : ! local variables
2397 :
2398 : real (kind=dbl_kind), dimension (nx_block,ny_block,2,max_blocks) :: &
2399 2406424320 : fldbundle ! work array for boundary updates
2400 :
2401 : character(len=*), parameter :: subname = '(dyn_haloUpdate2)'
2402 :
2403 1388160 : call ice_timer_start(timer_bound)
2404 : ! single process performs better without bundling fields
2405 1388160 : if (halo_dynbundle) then
2406 :
2407 230400 : call stack_fields(fld1, fld2, fldbundle)
2408 230400 : if (maskhalo_dyn) then
2409 0 : call ice_HaloUpdate (fldbundle, halo_info_mask, &
2410 0 : field_loc, field_type)
2411 : else
2412 0 : call ice_HaloUpdate (fldbundle, halo_info , &
2413 230400 : field_loc, field_type)
2414 : endif
2415 230400 : call unstack_fields(fldbundle, fld1, fld2)
2416 :
2417 : else
2418 :
2419 1157760 : if (maskhalo_dyn) then
2420 0 : call ice_HaloUpdate (fld1 , halo_info_mask, &
2421 0 : field_loc, field_type)
2422 0 : call ice_HaloUpdate (fld2 , halo_info_mask, &
2423 0 : field_loc, field_type)
2424 : else
2425 0 : call ice_HaloUpdate (fld1 , halo_info , &
2426 1157760 : field_loc, field_type)
2427 0 : call ice_HaloUpdate (fld2 , halo_info , &
2428 1157760 : field_loc, field_type)
2429 : endif
2430 :
2431 : endif
2432 1388160 : call ice_timer_stop(timer_bound)
2433 :
2434 1388160 : end subroutine dyn_haloUpdate2
2435 :
2436 : !=======================================================================
2437 : ! Do a halo update on 3 fields
2438 :
2439 0 : subroutine dyn_haloUpdate3(halo_info, halo_info_mask, field_loc, field_type, fld1, fld2, fld3)
2440 :
2441 : use ice_boundary, only: ice_halo, ice_HaloUpdate
2442 : use ice_domain, only: maskhalo_dyn, halo_dynbundle
2443 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2444 :
2445 : type (ice_halo), intent(in) :: &
2446 : halo_info , & ! standard unmasked halo ! LCOV_EXCL_LINE
2447 : halo_info_mask ! masked halo
2448 :
2449 : integer (kind=int_kind), intent(in) :: &
2450 : field_loc , & ! field loc ! LCOV_EXCL_LINE
2451 : field_type ! field_type
2452 :
2453 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
2454 : fld1 , & ! fields to halo ! LCOV_EXCL_LINE
2455 : fld2 , & ! ! LCOV_EXCL_LINE
2456 : fld3 !
2457 :
2458 : ! local variables
2459 :
2460 : real (kind=dbl_kind), dimension (nx_block,ny_block,3,max_blocks) :: &
2461 0 : fldbundle ! work array for boundary updates
2462 :
2463 : character(len=*), parameter :: subname = '(dyn_haloUpdate3)'
2464 :
2465 0 : call ice_timer_start(timer_bound)
2466 : ! single process performs better without bundling fields
2467 0 : if (halo_dynbundle) then
2468 :
2469 0 : call stack_fields(fld1, fld2, fld3, fldbundle)
2470 0 : if (maskhalo_dyn) then
2471 0 : call ice_HaloUpdate (fldbundle, halo_info_mask, &
2472 0 : field_loc, field_type)
2473 : else
2474 0 : call ice_HaloUpdate (fldbundle, halo_info , &
2475 0 : field_loc, field_type)
2476 : endif
2477 0 : call unstack_fields(fldbundle, fld1, fld2, fld3)
2478 :
2479 : else
2480 :
2481 0 : if (maskhalo_dyn) then
2482 0 : call ice_HaloUpdate (fld1 , halo_info_mask, &
2483 0 : field_loc, field_type)
2484 0 : call ice_HaloUpdate (fld2 , halo_info_mask, &
2485 0 : field_loc, field_type)
2486 0 : call ice_HaloUpdate (fld3 , halo_info_mask, &
2487 0 : field_loc, field_type)
2488 : else
2489 0 : call ice_HaloUpdate (fld1 , halo_info , &
2490 0 : field_loc, field_type)
2491 0 : call ice_HaloUpdate (fld2 , halo_info , &
2492 0 : field_loc, field_type)
2493 0 : call ice_HaloUpdate (fld3 , halo_info , &
2494 0 : field_loc, field_type)
2495 : endif
2496 :
2497 : endif
2498 0 : call ice_timer_stop(timer_bound)
2499 :
2500 0 : end subroutine dyn_haloUpdate3
2501 :
2502 : !=======================================================================
2503 : ! Do a halo update on 4 fields
2504 :
2505 0 : subroutine dyn_haloUpdate4(halo_info, halo_info_mask, field_loc, field_type, fld1, fld2, fld3, fld4)
2506 :
2507 : use ice_boundary, only: ice_halo, ice_HaloUpdate
2508 : use ice_domain, only: maskhalo_dyn, halo_dynbundle
2509 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2510 :
2511 : type (ice_halo), intent(in) :: &
2512 : halo_info , & ! standard unmasked halo ! LCOV_EXCL_LINE
2513 : halo_info_mask ! masked halo
2514 :
2515 : integer (kind=int_kind), intent(in) :: &
2516 : field_loc, & ! field loc ! LCOV_EXCL_LINE
2517 : field_type ! field_type
2518 :
2519 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
2520 : fld1 , & ! fields to halo ! LCOV_EXCL_LINE
2521 : fld2 , & ! ! LCOV_EXCL_LINE
2522 : fld3 , & ! ! LCOV_EXCL_LINE
2523 : fld4 !
2524 :
2525 : ! local variables
2526 :
2527 : real (kind=dbl_kind), dimension (nx_block,ny_block,4,max_blocks) :: &
2528 0 : fldbundle ! work array for boundary updates
2529 :
2530 : character(len=*), parameter :: subname = '(dyn_haloUpdate4)'
2531 :
2532 0 : call ice_timer_start(timer_bound)
2533 : ! single process performs better without bundling fields
2534 0 : if (halo_dynbundle) then
2535 :
2536 0 : call stack_fields(fld1, fld2, fld3, fld4, fldbundle)
2537 0 : if (maskhalo_dyn) then
2538 0 : call ice_HaloUpdate (fldbundle, halo_info_mask, &
2539 0 : field_loc, field_type)
2540 : else
2541 0 : call ice_HaloUpdate (fldbundle, halo_info , &
2542 0 : field_loc, field_type)
2543 : endif
2544 0 : call unstack_fields(fldbundle, fld1, fld2, fld3, fld4)
2545 :
2546 : else
2547 :
2548 0 : if (maskhalo_dyn) then
2549 0 : call ice_HaloUpdate (fld1 , halo_info_mask, &
2550 0 : field_loc, field_type)
2551 0 : call ice_HaloUpdate (fld2 , halo_info_mask, &
2552 0 : field_loc, field_type)
2553 0 : call ice_HaloUpdate (fld3 , halo_info_mask, &
2554 0 : field_loc, field_type)
2555 0 : call ice_HaloUpdate (fld4 , halo_info_mask, &
2556 0 : field_loc, field_type)
2557 : else
2558 0 : call ice_HaloUpdate (fld1 , halo_info , &
2559 0 : field_loc, field_type)
2560 0 : call ice_HaloUpdate (fld2 , halo_info , &
2561 0 : field_loc, field_type)
2562 0 : call ice_HaloUpdate (fld3 , halo_info , &
2563 0 : field_loc, field_type)
2564 0 : call ice_HaloUpdate (fld4 , halo_info , &
2565 0 : field_loc, field_type)
2566 : endif
2567 :
2568 : endif
2569 0 : call ice_timer_stop(timer_bound)
2570 :
2571 0 : end subroutine dyn_haloUpdate4
2572 :
2573 : !=======================================================================
2574 : ! Do a halo update on 5 fields
2575 :
2576 0 : subroutine dyn_haloUpdate5(halo_info, halo_info_mask, field_loc, field_type, fld1, fld2, fld3, fld4, fld5)
2577 :
2578 : use ice_boundary, only: ice_halo, ice_HaloUpdate
2579 : use ice_domain, only: maskhalo_dyn, halo_dynbundle
2580 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
2581 :
2582 : type (ice_halo), intent(in) :: &
2583 : halo_info , & ! standard unmasked halo ! LCOV_EXCL_LINE
2584 : halo_info_mask ! masked halo
2585 :
2586 : integer (kind=int_kind), intent(in) :: &
2587 : field_loc , & ! field loc ! LCOV_EXCL_LINE
2588 : field_type ! field_type
2589 :
2590 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(inout) :: &
2591 : fld1 , & ! fields to halo ! LCOV_EXCL_LINE
2592 : fld2 , & ! ! LCOV_EXCL_LINE
2593 : fld3 , & ! ! LCOV_EXCL_LINE
2594 : fld4 , & ! ! LCOV_EXCL_LINE
2595 : fld5 !
2596 :
2597 : ! local variables
2598 :
2599 : real (kind=dbl_kind), dimension (nx_block,ny_block,5,max_blocks) :: &
2600 0 : fldbundle ! work array for boundary updates
2601 :
2602 : character(len=*), parameter :: subname = '(dyn_haloUpdate5)'
2603 :
2604 0 : call ice_timer_start(timer_bound)
2605 : ! single process performs better without bundling fields
2606 0 : if (halo_dynbundle) then
2607 :
2608 0 : call stack_fields(fld1, fld2, fld3, fld4, fld5, fldbundle)
2609 0 : if (maskhalo_dyn) then
2610 0 : call ice_HaloUpdate (fldbundle, halo_info_mask, &
2611 0 : field_loc, field_type)
2612 : else
2613 0 : call ice_HaloUpdate (fldbundle, halo_info , &
2614 0 : field_loc, field_type)
2615 : endif
2616 0 : call unstack_fields(fldbundle, fld1, fld2, fld3, fld4, fld5)
2617 :
2618 : else
2619 :
2620 0 : if (maskhalo_dyn) then
2621 0 : call ice_HaloUpdate (fld1 , halo_info_mask, &
2622 0 : field_loc, field_type)
2623 0 : call ice_HaloUpdate (fld2 , halo_info_mask, &
2624 0 : field_loc, field_type)
2625 0 : call ice_HaloUpdate (fld3 , halo_info_mask, &
2626 0 : field_loc, field_type)
2627 0 : call ice_HaloUpdate (fld4 , halo_info_mask, &
2628 0 : field_loc, field_type)
2629 0 : call ice_HaloUpdate (fld5 , halo_info_mask, &
2630 0 : field_loc, field_type)
2631 : else
2632 0 : call ice_HaloUpdate (fld1 , halo_info , &
2633 0 : field_loc, field_type)
2634 0 : call ice_HaloUpdate (fld2 , halo_info , &
2635 0 : field_loc, field_type)
2636 0 : call ice_HaloUpdate (fld3 , halo_info , &
2637 0 : field_loc, field_type)
2638 0 : call ice_HaloUpdate (fld4 , halo_info , &
2639 0 : field_loc, field_type)
2640 0 : call ice_HaloUpdate (fld5 , halo_info , &
2641 0 : field_loc, field_type)
2642 : endif
2643 :
2644 : endif
2645 0 : call ice_timer_stop(timer_bound)
2646 :
2647 0 : end subroutine dyn_haloUpdate5
2648 :
2649 : !=======================================================================
2650 : ! Load fields into array for boundary updates
2651 :
2652 236184 : subroutine stack_fields2(fld1, fld2, fldbundle)
2653 :
2654 : use ice_domain, only: nblocks
2655 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2656 :
2657 : real (kind=dbl_kind), dimension (:,:,:), intent(in) :: &
2658 : fld1 , & ! fields to stack ! LCOV_EXCL_LINE
2659 : fld2 !
2660 :
2661 : real (kind=dbl_kind), dimension (:,:,:,:), intent(out) :: &
2662 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2663 :
2664 : ! local variables
2665 :
2666 : integer (kind=int_kind) :: &
2667 : iblk ! block index
2668 :
2669 : character(len=*), parameter :: subname = '(stack_fields2)'
2670 :
2671 236184 : call ice_timer_start(timer_bundbound)
2672 233280 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2673 8688 : do iblk = 1, nblocks
2674 7151880 : fldbundle(:,:,1,iblk) = fld1(:,:,iblk)
2675 7388064 : fldbundle(:,:,2,iblk) = fld2(:,:,iblk)
2676 : enddo
2677 : !$OMP END PARALLEL DO
2678 236184 : call ice_timer_stop(timer_bundbound)
2679 :
2680 236184 : end subroutine stack_fields2
2681 :
2682 : !=======================================================================
2683 : ! Load fields into array for boundary updates
2684 :
2685 5784 : subroutine stack_fields3(fld1, fld2, fld3, fldbundle)
2686 :
2687 : use ice_domain, only: nblocks
2688 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2689 :
2690 : real (kind=dbl_kind), dimension (:,:,:), intent(in) :: &
2691 : fld1 , & ! fields to stack ! LCOV_EXCL_LINE
2692 : fld2 , & ! ! LCOV_EXCL_LINE
2693 : fld3 !
2694 :
2695 : real (kind=dbl_kind), dimension (:,:,:,:), intent(out) :: &
2696 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2697 :
2698 : ! local variables
2699 :
2700 : integer (kind=int_kind) :: &
2701 : iblk ! block index
2702 :
2703 : character(len=*), parameter :: subname = '(stack_fields3)'
2704 :
2705 5784 : call ice_timer_start(timer_bundbound)
2706 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2707 8688 : do iblk = 1, nblocks
2708 7151880 : fldbundle(:,:,1,iblk) = fld1(:,:,iblk)
2709 7151880 : fldbundle(:,:,2,iblk) = fld2(:,:,iblk)
2710 7157664 : fldbundle(:,:,3,iblk) = fld3(:,:,iblk)
2711 : enddo
2712 : !$OMP END PARALLEL DO
2713 5784 : call ice_timer_stop(timer_bundbound)
2714 :
2715 5784 : end subroutine stack_fields3
2716 :
2717 : !=======================================================================
2718 : ! Load fields into array for boundary updates
2719 :
2720 5784 : subroutine stack_fields4(fld1, fld2, fld3, fld4, fldbundle)
2721 :
2722 : use ice_domain, only: nblocks
2723 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2724 :
2725 : real (kind=dbl_kind), dimension (:,:,:), intent(in) :: &
2726 : fld1 , & ! fields to stack ! LCOV_EXCL_LINE
2727 : fld2 , & ! ! LCOV_EXCL_LINE
2728 : fld3 , & ! ! LCOV_EXCL_LINE
2729 : fld4 !
2730 :
2731 : real (kind=dbl_kind), dimension (:,:,:,:), intent(out) :: &
2732 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2733 :
2734 : ! local variables
2735 :
2736 : integer (kind=int_kind) :: &
2737 : iblk ! block index
2738 :
2739 : character(len=*), parameter :: subname = '(stack_fields4)'
2740 :
2741 5784 : call ice_timer_start(timer_bundbound)
2742 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2743 8688 : do iblk = 1, nblocks
2744 7151880 : fldbundle(:,:,1,iblk) = fld1(:,:,iblk)
2745 7151880 : fldbundle(:,:,2,iblk) = fld2(:,:,iblk)
2746 7151880 : fldbundle(:,:,3,iblk) = fld3(:,:,iblk)
2747 7157664 : fldbundle(:,:,4,iblk) = fld4(:,:,iblk)
2748 : enddo
2749 : !$OMP END PARALLEL DO
2750 5784 : call ice_timer_stop(timer_bundbound)
2751 :
2752 5784 : end subroutine stack_fields4
2753 :
2754 : !=======================================================================
2755 : ! Load fields into array for boundary updates
2756 :
2757 0 : subroutine stack_fields5(fld1, fld2, fld3, fld4, fld5, fldbundle)
2758 :
2759 : use ice_domain, only: nblocks
2760 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2761 :
2762 : real (kind=dbl_kind), dimension (:,:,:), intent(in) :: &
2763 : fld1 , & ! fields to stack ! LCOV_EXCL_LINE
2764 : fld2 , & ! ! LCOV_EXCL_LINE
2765 : fld3 , & ! ! LCOV_EXCL_LINE
2766 : fld4 , & ! ! LCOV_EXCL_LINE
2767 : fld5 !
2768 :
2769 : real (kind=dbl_kind), dimension (:,:,:,:), intent(out) :: &
2770 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2771 :
2772 : ! local variables
2773 :
2774 : integer (kind=int_kind) :: &
2775 : iblk ! block index
2776 :
2777 : character(len=*), parameter :: subname = '(stack_fields5)'
2778 :
2779 0 : call ice_timer_start(timer_bundbound)
2780 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2781 0 : do iblk = 1, nblocks
2782 0 : fldbundle(:,:,1,iblk) = fld1(:,:,iblk)
2783 0 : fldbundle(:,:,2,iblk) = fld2(:,:,iblk)
2784 0 : fldbundle(:,:,3,iblk) = fld3(:,:,iblk)
2785 0 : fldbundle(:,:,4,iblk) = fld4(:,:,iblk)
2786 0 : fldbundle(:,:,5,iblk) = fld5(:,:,iblk)
2787 : enddo
2788 : !$OMP END PARALLEL DO
2789 0 : call ice_timer_stop(timer_bundbound)
2790 :
2791 0 : end subroutine stack_fields5
2792 :
2793 : !=======================================================================
2794 : ! Unload fields from array after boundary updates
2795 :
2796 236184 : subroutine unstack_fields2(fldbundle, fld1, fld2)
2797 :
2798 : use ice_domain, only: nblocks
2799 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2800 :
2801 : real (kind=dbl_kind), dimension (:,:,:,:), intent(in) :: &
2802 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2803 :
2804 : real (kind=dbl_kind), dimension (:,:,:), intent(out) :: &
2805 : fld1 , & ! fields to unstack ! LCOV_EXCL_LINE
2806 : fld2 !
2807 :
2808 : ! local variables
2809 :
2810 : integer (kind=int_kind) :: &
2811 : iblk ! block index
2812 :
2813 : character(len=*), parameter :: subname = '(unstack_fields2)'
2814 :
2815 236184 : call ice_timer_start(timer_bundbound)
2816 233280 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2817 8688 : do iblk = 1, nblocks
2818 7151880 : fld1(:,:,iblk) = fldbundle(:,:,1,iblk)
2819 7388064 : fld2(:,:,iblk) = fldbundle(:,:,2,iblk)
2820 : enddo
2821 : !$OMP END PARALLEL DO
2822 236184 : call ice_timer_stop(timer_bundbound)
2823 :
2824 236184 : end subroutine unstack_fields2
2825 :
2826 : !=======================================================================
2827 : ! Unload fields from array after boundary updates
2828 :
2829 5784 : subroutine unstack_fields3(fldbundle, fld1, fld2, fld3)
2830 :
2831 : use ice_domain, only: nblocks
2832 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2833 :
2834 : real (kind=dbl_kind), dimension (:,:,:,:), intent(in) :: &
2835 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2836 :
2837 : real (kind=dbl_kind), dimension (:,:,:), intent(out) :: &
2838 : fld1 , & ! fields to unstack ! LCOV_EXCL_LINE
2839 : fld2 , & ! ! LCOV_EXCL_LINE
2840 : fld3 !
2841 :
2842 : ! local variables
2843 :
2844 : integer (kind=int_kind) :: &
2845 : iblk ! block index
2846 :
2847 : character(len=*), parameter :: subname = '(unstack_fields3)'
2848 :
2849 5784 : call ice_timer_start(timer_bundbound)
2850 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2851 8688 : do iblk = 1, nblocks
2852 7151880 : fld1(:,:,iblk) = fldbundle(:,:,1,iblk)
2853 7151880 : fld2(:,:,iblk) = fldbundle(:,:,2,iblk)
2854 7157664 : fld3(:,:,iblk) = fldbundle(:,:,3,iblk)
2855 : enddo
2856 : !$OMP END PARALLEL DO
2857 5784 : call ice_timer_stop(timer_bundbound)
2858 :
2859 5784 : end subroutine unstack_fields3
2860 :
2861 : !=======================================================================
2862 : ! Unload fields from array after boundary updates
2863 :
2864 5784 : subroutine unstack_fields4(fldbundle, fld1, fld2, fld3, fld4)
2865 :
2866 : use ice_domain, only: nblocks
2867 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2868 :
2869 : real (kind=dbl_kind), dimension (:,:,:,:), intent(in) :: &
2870 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2871 :
2872 : real (kind=dbl_kind), dimension (:,:,:), intent(out) :: &
2873 : fld1 , & ! fields to unstack ! LCOV_EXCL_LINE
2874 : fld2 , & ! ! LCOV_EXCL_LINE
2875 : fld3 , & ! ! LCOV_EXCL_LINE
2876 : fld4 !
2877 :
2878 : ! local variables
2879 :
2880 : integer (kind=int_kind) :: &
2881 : iblk ! block index
2882 :
2883 : character(len=*), parameter :: subname = '(unstack_fields4)'
2884 :
2885 5784 : call ice_timer_start(timer_bundbound)
2886 2880 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2887 8688 : do iblk = 1, nblocks
2888 7151880 : fld1(:,:,iblk) = fldbundle(:,:,1,iblk)
2889 7151880 : fld2(:,:,iblk) = fldbundle(:,:,2,iblk)
2890 7151880 : fld3(:,:,iblk) = fldbundle(:,:,3,iblk)
2891 7157664 : fld4(:,:,iblk) = fldbundle(:,:,4,iblk)
2892 : enddo
2893 : !$OMP END PARALLEL DO
2894 5784 : call ice_timer_stop(timer_bundbound)
2895 :
2896 5784 : end subroutine unstack_fields4
2897 :
2898 : !=======================================================================
2899 : ! Unload fields from array after boundary updates
2900 :
2901 0 : subroutine unstack_fields5(fldbundle, fld1, fld2, fld3, fld4, fld5)
2902 :
2903 : use ice_domain, only: nblocks
2904 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bundbound
2905 :
2906 : real (kind=dbl_kind), dimension (:,:,:,:), intent(in) :: &
2907 : fldbundle ! work array for boundary updates (i,j,n,iblk)
2908 :
2909 : real (kind=dbl_kind), dimension (:,:,:), intent(out) :: &
2910 : fld1 , & ! fields to unstack ! LCOV_EXCL_LINE
2911 : fld2 , & ! ! LCOV_EXCL_LINE
2912 : fld3 , & ! ! LCOV_EXCL_LINE
2913 : fld4 , & ! ! LCOV_EXCL_LINE
2914 : fld5 !
2915 :
2916 : ! local variables
2917 :
2918 : integer (kind=int_kind) :: &
2919 : iblk ! block index
2920 :
2921 : character(len=*), parameter :: subname = '(unstack_fields5)'
2922 :
2923 0 : call ice_timer_start(timer_bundbound)
2924 0 : !$OMP PARALLEL DO PRIVATE(iblk) SCHEDULE(runtime)
2925 0 : do iblk = 1, nblocks
2926 0 : fld1(:,:,iblk) = fldbundle(:,:,1,iblk)
2927 0 : fld2(:,:,iblk) = fldbundle(:,:,2,iblk)
2928 0 : fld3(:,:,iblk) = fldbundle(:,:,3,iblk)
2929 0 : fld4(:,:,iblk) = fldbundle(:,:,4,iblk)
2930 0 : fld5(:,:,iblk) = fldbundle(:,:,5,iblk)
2931 : enddo
2932 : !$OMP END PARALLEL DO
2933 0 : call ice_timer_stop(timer_bundbound)
2934 :
2935 0 : end subroutine unstack_fields5
2936 :
2937 : !=======================================================================
2938 :
2939 : end module ice_dyn_shared
2940 :
2941 : !=======================================================================
|
