Line data Source code
1 : !=======================================================================
2 : !
3 : ! Elastic-viscous-plastic sea ice dynamics model (1D implementations)
4 : ! Computes ice velocity and deformation
5 : !
6 : ! authors: Jacob Weismann Poulsen, DMI
7 : ! Mads Hvid Ribergaard, DMI
8 :
9 : module ice_dyn_evp_1d
10 :
11 : use ice_kinds_mod
12 : use ice_fileunits, only : nu_diag
13 : use ice_exit, only : abort_ice
14 : use icepack_intfc, only : icepack_query_parameters
15 : use icepack_intfc, only : icepack_warnings_flush, &
16 : icepack_warnings_aborted
17 :
18 : implicit none
19 : private
20 : public :: ice_dyn_evp_1d_copyin, ice_dyn_evp_1d_copyout, &
21 : ice_dyn_evp_1d_kernel
22 :
23 : integer(kind=int_kind) :: NA_len, NAVEL_len, domp_iam, domp_nt
24 : #if defined (_OPENMP)
25 : real(kind=dbl_kind) :: rdomp_iam, rdomp_nt
26 : !$OMP THREADPRIVATE(domp_iam, domp_nt, rdomp_iam, rdomp_nt)
27 : #endif
28 : logical(kind=log_kind), dimension(:), allocatable :: skiptcell, skipucell
29 : integer(kind=int_kind), dimension(:), allocatable :: ee, ne, se, &
30 : nw, sw, sse, indi, indj, indij, halo_parent
31 : real(kind=dbl_kind), dimension(:), allocatable :: cdn_ocn, aiu, &
32 : uocn, vocn, forcex, forcey, Tbu, tarear, umassdti, fm, uarear, & ! LCOV_EXCL_LINE
33 : strintx, strinty, uvel_init, vvel_init, strength, uvel, vvel, & ! LCOV_EXCL_LINE
34 : dxT, dyT, stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, & ! LCOV_EXCL_LINE
35 : stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, & ! LCOV_EXCL_LINE
36 : stress12_3, stress12_4, divu, rdg_conv, rdg_shear, shear, taubx, & ! LCOV_EXCL_LINE
37 : tauby, str1, str2, str3, str4, str5, str6, str7, str8, HTE, HTN, & ! LCOV_EXCL_LINE
38 : HTEm1, HTNm1
39 : integer, parameter :: JPIM = selected_int_kind(9)
40 :
41 : interface evp1d_stress
42 : module procedure stress_iter
43 : module procedure stress_last
44 : end interface
45 :
46 : interface evp1d_stepu
47 : module procedure stepu_iter
48 : module procedure stepu_last
49 : end interface
50 :
51 : !=======================================================================
52 :
53 : contains
54 :
55 : !=======================================================================
56 :
57 0 : subroutine domp_init
58 : #if defined (_OPENMP)
59 :
60 : use omp_lib, only : omp_get_thread_num, omp_get_num_threads
61 : #endif
62 :
63 : implicit none
64 :
65 : character(len=*), parameter :: subname = '(domp_init)'
66 :
67 : !$OMP PARALLEL DEFAULT(none)
68 : #if defined (_OPENMP)
69 : domp_iam = omp_get_thread_num()
70 : rdomp_iam = real(domp_iam, dbl_kind)
71 : domp_nt = omp_get_num_threads()
72 : rdomp_nt = real(domp_nt, dbl_kind)
73 : #else
74 0 : domp_iam = 0
75 0 : domp_nt = 1
76 : #endif
77 : !$OMP END PARALLEL
78 :
79 0 : end subroutine domp_init
80 :
81 : !=======================================================================
82 :
83 0 : subroutine domp_get_domain(lower, upper, d_lower, d_upper)
84 : #if defined (_OPENMP)
85 :
86 : use omp_lib, only : omp_in_parallel
87 : use ice_constants, only : p5
88 : #endif
89 :
90 : implicit none
91 :
92 : integer(kind=JPIM), intent(in) :: lower, upper
93 : integer(kind=JPIM), intent(out) :: d_lower, d_upper
94 :
95 : ! local variables
96 : #if defined (_OPENMP)
97 :
98 0 : real(kind=dbl_kind) :: dlen
99 : #endif
100 :
101 : character(len=*), parameter :: subname = '(domp_get_domain)'
102 :
103 : ! proper action in "null" case
104 0 : if (upper <= 0 .or. upper < lower) then
105 0 : d_lower = 0
106 0 : d_upper = -1
107 0 : return
108 : end if
109 :
110 : ! proper action in serial case
111 0 : d_lower = lower
112 0 : d_upper = upper
113 : #if defined (_OPENMP)
114 :
115 0 : if (omp_in_parallel()) then
116 0 : dlen = real((upper - lower + 1), dbl_kind)
117 0 : d_lower = lower + floor(((rdomp_iam * dlen + p5) / rdomp_nt), JPIM)
118 0 : d_upper = lower - 1 + floor(((rdomp_iam * dlen + dlen + p5) / rdomp_nt), JPIM)
119 : end if
120 : #endif
121 :
122 : end subroutine domp_get_domain
123 :
124 : !=======================================================================
125 :
126 0 : subroutine stress_iter(NA_len, ee, ne, se, lb, ub, uvel, vvel, dxT, &
127 : dyT, hte, htn, htem1, htnm1, strength, stressp_1, stressp_2, & ! LCOV_EXCL_LINE
128 : stressp_3, stressp_4, stressm_1, stressm_2, stressm_3, & ! LCOV_EXCL_LINE
129 : stressm_4, stress12_1, stress12_2, stress12_3, stress12_4, str1, & ! LCOV_EXCL_LINE
130 0 : str2, str3, str4, str5, str6, str7, str8, skiptcell)
131 :
132 : use ice_kinds_mod
133 : use ice_constants, only : p027, p055, p111, p166, p222, p25, &
134 : p333, p5, c1p5, c1
135 : use ice_dyn_shared, only : ecci, denom1, arlx1i, Ktens, revp, &
136 : deltaminEVP
137 :
138 : implicit none
139 :
140 : integer(kind=int_kind), intent(in) :: NA_len, lb, ub
141 : integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
142 : ee, ne, se
143 : real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
144 : strength, uvel, vvel, dxT, dyT, hte, htn, htem1, htnm1
145 : logical(kind=log_kind), intent(in), dimension(:) :: skiptcell
146 : real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
147 : stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, & ! LCOV_EXCL_LINE
148 : stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, & ! LCOV_EXCL_LINE
149 : stress12_3, stress12_4
150 : real(kind=dbl_kind), dimension(:), intent(out), contiguous :: &
151 : str1, str2, str3, str4, str5, str6, str7, str8
152 :
153 : ! local variables
154 :
155 : integer(kind=int_kind) :: iw, il, iu
156 0 : real(kind=dbl_kind) :: divune, divunw, divuse, divusw, &
157 : tensionne, tensionnw, tensionse, tensionsw, shearne, shearnw, & ! LCOV_EXCL_LINE
158 : shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw, c0ne, & ! LCOV_EXCL_LINE
159 : c0nw, c0se, c0sw, c1ne, c1nw, c1se, c1sw, ssigpn, ssigps, & ! LCOV_EXCL_LINE
160 : ssigpe, ssigpw, ssigmn, ssigms, ssigme, ssigmw, ssig12n, & ! LCOV_EXCL_LINE
161 : ssig12s, ssig12e, ssig12w, ssigp1, ssigp2, ssigm1, ssigm2, & ! LCOV_EXCL_LINE
162 : ssig121, ssig122, csigpne, csigpnw, csigpse, csigpsw, & ! LCOV_EXCL_LINE
163 : csigmne, csigmnw, csigmse, csigmsw, csig12ne, csig12nw, & ! LCOV_EXCL_LINE
164 : csig12se, csig12sw, str12ew, str12we, str12ns, str12sn, & ! LCOV_EXCL_LINE
165 : strp_tmp, strm_tmp, tmp_uvel_ee, tmp_vvel_se, tmp_vvel_ee, & ! LCOV_EXCL_LINE
166 : tmp_vvel_ne, tmp_uvel_ne, tmp_uvel_se, dxhy, dyhx, cxp, cyp, & ! LCOV_EXCL_LINE
167 0 : cxm, cym, tmparea, DminTarea
168 :
169 : character(len=*), parameter :: subname = '(stress_iter)'
170 :
171 : #ifdef _OPENACC
172 : !$acc parallel &
173 : !$acc present(ee, ne, se, strength, uvel, vvel, dxT, dyT, hte, & ! LCOV_EXCL_LINE
174 : !$acc htn, htem1, htnm1, str1, str2, str3, str4, str5, str6, & ! LCOV_EXCL_LINE
175 : !$acc str7, str8, stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
176 : !$acc stressm_1, stressm_2, stressm_3, stressm_4, stress12_1, & ! LCOV_EXCL_LINE
177 : !$acc stress12_2, stress12_3, stress12_4, skiptcell)
178 : !$acc loop
179 : do iw = 1, NA_len
180 : #else
181 0 : call domp_get_domain(lb, ub, il, iu)
182 0 : do iw = il, iu
183 : #endif
184 :
185 0 : if (skiptcell(iw)) cycle
186 :
187 0 : tmparea = dxT(iw) * dyT(iw) ! necessary to split calc of DminTarea. Otherwize not binary identical
188 0 : DminTarea = deltaminEVP * tmparea
189 0 : dxhy = p5 * (hte(iw) - htem1(iw))
190 0 : dyhx = p5 * (htn(iw) - htnm1(iw))
191 0 : cxp = c1p5 * htn(iw) - p5 * htnm1(iw)
192 0 : cyp = c1p5 * hte(iw) - p5 * htem1(iw)
193 0 : cxm = -(c1p5 * htnm1(iw) - p5 * htn(iw))
194 0 : cym = -(c1p5 * htem1(iw) - p5 * hte(iw))
195 :
196 : !--------------------------------------------------------------
197 : ! strain rates
198 : ! NOTE: these are actually strain rates * area (m^2/s)
199 : !--------------------------------------------------------------
200 :
201 0 : tmp_uvel_ne = uvel(ne(iw))
202 0 : tmp_uvel_se = uvel(se(iw))
203 0 : tmp_uvel_ee = uvel(ee(iw))
204 :
205 0 : tmp_vvel_ee = vvel(ee(iw))
206 0 : tmp_vvel_se = vvel(se(iw))
207 0 : tmp_vvel_ne = vvel(ne(iw))
208 : ! divergence = e_11 + e_22
209 0 : divune = cyp * uvel(iw) - dyT(iw) * tmp_uvel_ee &
210 0 : + cxp * vvel(iw) - dxT(iw) * tmp_vvel_se
211 0 : divunw = cym * tmp_uvel_ee + dyT(iw) * uvel(iw) &
212 0 : + cxp * tmp_vvel_ee - dxT(iw) * tmp_vvel_ne
213 0 : divusw = cym * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
214 0 : + cxm * tmp_vvel_ne + dxT(iw) * tmp_vvel_ee
215 0 : divuse = cyp * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
216 0 : + cxm * tmp_vvel_se + dxT(iw) * vvel(iw)
217 :
218 : ! tension strain rate = e_11 - e_22
219 0 : tensionne = -cym * uvel(iw) - dyT(iw) * tmp_uvel_ee &
220 0 : + cxm * vvel(iw) + dxT(iw) * tmp_vvel_se
221 0 : tensionnw = -cyp * tmp_uvel_ee + dyT(iw) * uvel(iw) &
222 0 : + cxm * tmp_vvel_ee + dxT(iw) * tmp_vvel_ne
223 0 : tensionsw = -cyp * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
224 0 : + cxp * tmp_vvel_ne - dxT(iw) * tmp_vvel_ee
225 0 : tensionse = -cym * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
226 0 : + cxp * tmp_vvel_se - dxT(iw) * vvel(iw)
227 :
228 : ! shearing strain rate = 2 * e_12
229 0 : shearne = -cym * vvel(iw) - dyT(iw) * tmp_vvel_ee &
230 0 : - cxm * uvel(iw) - dxT(iw) * tmp_uvel_se
231 0 : shearnw = -cyp * tmp_vvel_ee + dyT(iw) * vvel(iw) &
232 0 : - cxm * tmp_uvel_ee - dxT(iw) * tmp_uvel_ne
233 0 : shearsw = -cyp * tmp_vvel_ne + dyT(iw) * tmp_vvel_se &
234 0 : - cxp * tmp_uvel_ne + dxT(iw) * tmp_uvel_ee
235 0 : shearse = -cym * tmp_vvel_se - dyT(iw) * tmp_vvel_ne &
236 0 : - cxp * tmp_uvel_se + dxT(iw) * uvel(iw)
237 :
238 : ! Delta (in the denominator of zeta and eta)
239 0 : Deltane = sqrt(divune**2 + ecci * (tensionne**2 + shearne**2))
240 0 : Deltanw = sqrt(divunw**2 + ecci * (tensionnw**2 + shearnw**2))
241 0 : Deltasw = sqrt(divusw**2 + ecci * (tensionsw**2 + shearsw**2))
242 0 : Deltase = sqrt(divuse**2 + ecci * (tensionse**2 + shearse**2))
243 :
244 : !--------------------------------------------------------------
245 : ! replacement pressure/Delta (kg/s)
246 : ! save replacement pressure for principal stress calculation
247 : !--------------------------------------------------------------
248 :
249 0 : c0ne = strength(iw) / max(Deltane, DminTarea)
250 0 : c0nw = strength(iw) / max(Deltanw, DminTarea)
251 0 : c0sw = strength(iw) / max(Deltasw, DminTarea)
252 0 : c0se = strength(iw) / max(Deltase, DminTarea)
253 :
254 0 : c1ne = c0ne * arlx1i
255 0 : c1nw = c0nw * arlx1i
256 0 : c1sw = c0sw * arlx1i
257 0 : c1se = c0se * arlx1i
258 :
259 0 : c0ne = c1ne * ecci
260 0 : c0nw = c1nw * ecci
261 0 : c0sw = c1sw * ecci
262 0 : c0se = c1se * ecci
263 :
264 : !--------------------------------------------------------------
265 : ! the stresses (kg/s^2)
266 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
267 : !--------------------------------------------------------------
268 :
269 0 : stressp_1(iw) = (stressp_1(iw) * (c1 - arlx1i * revp) &
270 0 : + c1ne * (divune * (c1 + Ktens) - Deltane * (c1 - Ktens))) * denom1
271 0 : stressp_2(iw) = (stressp_2(iw) * (c1 - arlx1i * revp) &
272 0 : + c1nw * (divunw * (c1 + Ktens) - Deltanw * (c1 - Ktens))) * denom1
273 0 : stressp_3(iw) = (stressp_3(iw) * (c1 - arlx1i * revp) &
274 0 : + c1sw * (divusw * (c1 + Ktens) - Deltasw * (c1 - Ktens))) * denom1
275 0 : stressp_4(iw) = (stressp_4(iw) * (c1 - arlx1i * revp) &
276 0 : + c1se * (divuse * (c1 + Ktens) - Deltase * (c1 - Ktens))) * denom1
277 :
278 0 : stressm_1(iw) = (stressm_1(iw) * (c1 - arlx1i * revp) + c0ne * tensionne * (c1 + Ktens)) * denom1
279 0 : stressm_2(iw) = (stressm_2(iw) * (c1 - arlx1i * revp) + c0nw * tensionnw * (c1 + Ktens)) * denom1
280 0 : stressm_3(iw) = (stressm_3(iw) * (c1 - arlx1i * revp) + c0sw * tensionsw * (c1 + Ktens)) * denom1
281 0 : stressm_4(iw) = (stressm_4(iw) * (c1 - arlx1i * revp) + c0se * tensionse * (c1 + Ktens)) * denom1
282 :
283 0 : stress12_1(iw) = (stress12_1(iw) * (c1 - arlx1i * revp) + c0ne * shearne * p5 * (c1 + Ktens)) * denom1
284 0 : stress12_2(iw) = (stress12_2(iw) * (c1 - arlx1i * revp) + c0nw * shearnw * p5 * (c1 + Ktens)) * denom1
285 0 : stress12_3(iw) = (stress12_3(iw) * (c1 - arlx1i * revp) + c0sw * shearsw * p5 * (c1 + Ktens)) * denom1
286 0 : stress12_4(iw) = (stress12_4(iw) * (c1 - arlx1i * revp) + c0se * shearse * p5 * (c1 + Ktens)) * denom1
287 :
288 : !--------------------------------------------------------------
289 : ! combinations of the stresses for the momentum equation
290 : ! (kg/s^2)
291 : !--------------------------------------------------------------
292 :
293 0 : ssigpn = stressp_1(iw) + stressp_2(iw)
294 0 : ssigps = stressp_3(iw) + stressp_4(iw)
295 0 : ssigpe = stressp_1(iw) + stressp_4(iw)
296 0 : ssigpw = stressp_2(iw) + stressp_3(iw)
297 0 : ssigp1 = (stressp_1(iw) + stressp_3(iw)) * p055
298 0 : ssigp2 = (stressp_2(iw) + stressp_4(iw)) * p055
299 :
300 0 : ssigmn = stressm_1(iw) + stressm_2(iw)
301 0 : ssigms = stressm_3(iw) + stressm_4(iw)
302 0 : ssigme = stressm_1(iw) + stressm_4(iw)
303 0 : ssigmw = stressm_2(iw) + stressm_3(iw)
304 0 : ssigm1 = (stressm_1(iw) + stressm_3(iw)) * p055
305 0 : ssigm2 = (stressm_2(iw) + stressm_4(iw)) * p055
306 :
307 0 : ssig12n = stress12_1(iw) + stress12_2(iw)
308 0 : ssig12s = stress12_3(iw) + stress12_4(iw)
309 0 : ssig12e = stress12_1(iw) + stress12_4(iw)
310 0 : ssig12w = stress12_2(iw) + stress12_3(iw)
311 0 : ssig121 = (stress12_1(iw) + stress12_3(iw)) * p111
312 0 : ssig122 = (stress12_2(iw) + stress12_4(iw)) * p111
313 :
314 0 : csigpne = p111 * stressp_1(iw) + ssigp2 + p027 * stressp_3(iw)
315 0 : csigpnw = p111 * stressp_2(iw) + ssigp1 + p027 * stressp_4(iw)
316 0 : csigpsw = p111 * stressp_3(iw) + ssigp2 + p027 * stressp_1(iw)
317 0 : csigpse = p111 * stressp_4(iw) + ssigp1 + p027 * stressp_2(iw)
318 :
319 0 : csigmne = p111 * stressm_1(iw) + ssigm2 + p027 * stressm_3(iw)
320 0 : csigmnw = p111 * stressm_2(iw) + ssigm1 + p027 * stressm_4(iw)
321 0 : csigmsw = p111 * stressm_3(iw) + ssigm2 + p027 * stressm_1(iw)
322 0 : csigmse = p111 * stressm_4(iw) + ssigm1 + p027 * stressm_2(iw)
323 :
324 0 : csig12ne = p222 * stress12_1(iw) + ssig122 + p055 * stress12_3(iw)
325 0 : csig12nw = p222 * stress12_2(iw) + ssig121 + p055 * stress12_4(iw)
326 0 : csig12sw = p222 * stress12_3(iw) + ssig122 + p055 * stress12_1(iw)
327 0 : csig12se = p222 * stress12_4(iw) + ssig121 + p055 * stress12_2(iw)
328 :
329 0 : str12ew = p5 * dxT(iw) * (p333 * ssig12e + p166 * ssig12w)
330 0 : str12we = p5 * dxT(iw) * (p333 * ssig12w + p166 * ssig12e)
331 0 : str12ns = p5 * dyT(iw) * (p333 * ssig12n + p166 * ssig12s)
332 0 : str12sn = p5 * dyT(iw) * (p333 * ssig12s + p166 * ssig12n)
333 :
334 : !--------------------------------------------------------------
335 : ! for dF/dx (u momentum)
336 : !--------------------------------------------------------------
337 :
338 0 : strp_tmp = p25 * dyT(iw) * (p333 * ssigpn + p166 * ssigps)
339 0 : strm_tmp = p25 * dyT(iw) * (p333 * ssigmn + p166 * ssigms)
340 :
341 : ! northeast (i,j)
342 0 : str1(iw) = -strp_tmp - strm_tmp - str12ew &
343 0 : + dxhy * (-csigpne + csigmne) + dyhx * csig12ne
344 :
345 : ! northwest (i+1,j)
346 0 : str2(iw) = strp_tmp + strm_tmp - str12we &
347 0 : + dxhy * (-csigpnw + csigmnw) + dyhx * csig12nw
348 :
349 0 : strp_tmp = p25 * dyT(iw) * (p333 * ssigps + p166 * ssigpn)
350 0 : strm_tmp = p25 * dyT(iw) * (p333 * ssigms + p166 * ssigmn)
351 :
352 : ! southeast (i,j+1)
353 0 : str3(iw) = -strp_tmp - strm_tmp + str12ew &
354 0 : + dxhy * (-csigpse + csigmse) + dyhx * csig12se
355 :
356 : ! southwest (i+1,j+1)
357 0 : str4(iw) = strp_tmp + strm_tmp + str12we &
358 0 : + dxhy * (-csigpsw + csigmsw) + dyhx * csig12sw
359 :
360 : !--------------------------------------------------------------
361 : ! for dF/dy (v momentum)
362 : !--------------------------------------------------------------
363 :
364 0 : strp_tmp = p25 * dxT(iw) * (p333 * ssigpe + p166 * ssigpw)
365 0 : strm_tmp = p25 * dxT(iw) * (p333 * ssigme + p166 * ssigmw)
366 :
367 : ! northeast (i,j)
368 0 : str5(iw) = -strp_tmp + strm_tmp - str12ns &
369 0 : - dyhx * (csigpne + csigmne) + dxhy * csig12ne
370 :
371 : ! southeast (i,j+1)
372 0 : str6(iw) = strp_tmp - strm_tmp - str12sn &
373 0 : - dyhx * (csigpse + csigmse) + dxhy * csig12se
374 :
375 0 : strp_tmp = p25 * dxT(iw) * (p333 * ssigpw + p166 * ssigpe)
376 0 : strm_tmp = p25 * dxT(iw) * (p333 * ssigmw + p166 * ssigme)
377 :
378 : ! northwest (i+1,j)
379 0 : str7(iw) = -strp_tmp + strm_tmp + str12ns &
380 0 : - dyhx * (csigpnw + csigmnw) + dxhy * csig12nw
381 :
382 : ! southwest (i+1,j+1)
383 0 : str8(iw) = strp_tmp - strm_tmp + str12sn &
384 0 : - dyhx * (csigpsw + csigmsw) + dxhy * csig12sw
385 :
386 : end do
387 : #ifdef _OPENACC
388 : !$acc end parallel
389 : #endif
390 :
391 0 : end subroutine stress_iter
392 :
393 : !=======================================================================
394 :
395 0 : subroutine stress_last(NA_len, ee, ne, se, lb, ub, uvel, vvel, dxT, &
396 : dyT, hte, htn, htem1, htnm1, strength, stressp_1, stressp_2, & ! LCOV_EXCL_LINE
397 : stressp_3, stressp_4, stressm_1, stressm_2, stressm_3, & ! LCOV_EXCL_LINE
398 : stressm_4, stress12_1, stress12_2, stress12_3, stress12_4, str1, & ! LCOV_EXCL_LINE
399 : str2, str3, str4, str5, str6, str7, str8, skiptcell, tarear, divu, & ! LCOV_EXCL_LINE
400 0 : rdg_conv, rdg_shear, shear)
401 :
402 : use ice_kinds_mod
403 : use ice_constants, only : p027, p055, p111, p166, p222, p25, &
404 : p333, p5, c1p5, c1, c0
405 : use ice_dyn_shared, only : ecci, denom1, arlx1i, Ktens, revp,&
406 : deltaminEVP
407 :
408 : implicit none
409 :
410 : integer(kind=int_kind), intent(in) :: NA_len, lb, ub
411 : integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
412 : ee, ne, se
413 : real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
414 : strength, uvel, vvel, dxT, dyT, hte, htn, htem1, htnm1, tarear
415 : logical(kind=log_kind), intent(in), dimension(:) :: skiptcell
416 : real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
417 : stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, & ! LCOV_EXCL_LINE
418 : stressm_2, stressm_3, stressm_4, stress12_1, stress12_2, & ! LCOV_EXCL_LINE
419 : stress12_3, stress12_4
420 : real(kind=dbl_kind), dimension(:), intent(out), contiguous :: &
421 : str1, str2, str3, str4, str5, str6, str7, str8, divu, & ! LCOV_EXCL_LINE
422 : rdg_conv, rdg_shear, shear
423 :
424 : ! local variables
425 :
426 : integer(kind=int_kind) :: iw, il, iu
427 0 : real(kind=dbl_kind) :: divune, divunw, divuse, divusw, &
428 : tensionne, tensionnw, tensionse, tensionsw, shearne, shearnw, & ! LCOV_EXCL_LINE
429 : shearse, shearsw, Deltane, Deltanw, Deltase, Deltasw, c0ne, & ! LCOV_EXCL_LINE
430 : c0nw, c0se, c0sw, c1ne, c1nw, c1se, c1sw, ssigpn, ssigps, & ! LCOV_EXCL_LINE
431 : ssigpe, ssigpw, ssigmn, ssigms, ssigme, ssigmw, ssig12n, & ! LCOV_EXCL_LINE
432 : ssig12s, ssig12e, ssig12w, ssigp1, ssigp2, ssigm1, ssigm2, & ! LCOV_EXCL_LINE
433 : ssig121, ssig122, csigpne, csigpnw, csigpse, csigpsw, & ! LCOV_EXCL_LINE
434 : csigmne, csigmnw, csigmse, csigmsw, csig12ne, csig12nw, & ! LCOV_EXCL_LINE
435 : csig12se, csig12sw, str12ew, str12we, str12ns, str12sn, & ! LCOV_EXCL_LINE
436 : strp_tmp, strm_tmp, tmp_uvel_ee, tmp_vvel_se, tmp_vvel_ee, & ! LCOV_EXCL_LINE
437 : tmp_vvel_ne, tmp_uvel_ne, tmp_uvel_se, dxhy, dyhx, cxp, cyp, & ! LCOV_EXCL_LINE
438 0 : cxm, cym, tmparea, DminTarea
439 :
440 : character(len=*), parameter :: subname = '(stress_last)'
441 :
442 : #ifdef _OPENACC
443 : !$acc parallel &
444 : !$acc present(ee, ne, se, strength, uvel, vvel, dxT, dyT, hte, & ! LCOV_EXCL_LINE
445 : !$acc htn, htem1, htnm1, str1, str2, str3, str4, str5, str6, & ! LCOV_EXCL_LINE
446 : !$acc str7, str8, stressp_1, stressp_2, stressp_3, stressp_4, & ! LCOV_EXCL_LINE
447 : !$acc stressm_1, stressm_2, stressm_3, stressm_4, stress12_1, & ! LCOV_EXCL_LINE
448 : !$acc stress12_2, stress12_3, stress12_4, tarear, divu, & ! LCOV_EXCL_LINE
449 : !$acc rdg_conv, rdg_shear, shear, skiptcell)
450 : !$acc loop
451 : do iw = 1, NA_len
452 : #else
453 0 : call domp_get_domain(lb, ub, il, iu)
454 0 : do iw = il, iu
455 : #endif
456 :
457 0 : if (skiptcell(iw)) cycle
458 :
459 0 : tmparea = dxT(iw) * dyT(iw) ! necessary to split calc of DminTarea. Otherwize not binary identical
460 0 : DminTarea = deltaminEVP * tmparea
461 0 : dxhy = p5 * (hte(iw) - htem1(iw))
462 0 : dyhx = p5 * (htn(iw) - htnm1(iw))
463 0 : cxp = c1p5 * htn(iw) - p5 * htnm1(iw)
464 0 : cyp = c1p5 * hte(iw) - p5 * htem1(iw)
465 0 : cxm = -(c1p5 * htnm1(iw) - p5 * htn(iw))
466 0 : cym = -(c1p5 * htem1(iw) - p5 * hte(iw))
467 :
468 : !--------------------------------------------------------------
469 : ! strain rates
470 : ! NOTE: these are actually strain rates * area (m^2/s)
471 : !--------------------------------------------------------------
472 :
473 0 : tmp_uvel_ne = uvel(ne(iw))
474 0 : tmp_uvel_se = uvel(se(iw))
475 0 : tmp_uvel_ee = uvel(ee(iw))
476 :
477 0 : tmp_vvel_ee = vvel(ee(iw))
478 0 : tmp_vvel_se = vvel(se(iw))
479 0 : tmp_vvel_ne = vvel(ne(iw))
480 :
481 : ! divergence = e_11 + e_22
482 0 : divune = cyp * uvel(iw) - dyT(iw) * tmp_uvel_ee &
483 0 : + cxp * vvel(iw) - dxT(iw) * tmp_vvel_se
484 0 : divunw = cym * tmp_uvel_ee + dyT(iw) * uvel(iw) &
485 0 : + cxp * tmp_vvel_ee - dxT(iw) * tmp_vvel_ne
486 0 : divusw = cym * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
487 0 : + cxm * tmp_vvel_ne + dxT(iw) * tmp_vvel_ee
488 0 : divuse = cyp * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
489 0 : + cxm * tmp_vvel_se + dxT(iw) * vvel(iw)
490 :
491 : ! tension strain rate = e_11 - e_22
492 0 : tensionne = -cym * uvel(iw) - dyT(iw) * tmp_uvel_ee &
493 0 : + cxm * vvel(iw) + dxT(iw) * tmp_vvel_se
494 0 : tensionnw = -cyp * tmp_uvel_ee + dyT(iw) * uvel(iw) &
495 0 : + cxm * tmp_vvel_ee + dxT(iw) * tmp_vvel_ne
496 0 : tensionsw = -cyp * tmp_uvel_ne + dyT(iw) * tmp_uvel_se &
497 0 : + cxp * tmp_vvel_ne - dxT(iw) * tmp_vvel_ee
498 0 : tensionse = -cym * tmp_uvel_se - dyT(iw) * tmp_uvel_ne &
499 0 : + cxp * tmp_vvel_se - dxT(iw) * vvel(iw)
500 :
501 : ! shearing strain rate = 2 * e_12
502 0 : shearne = -cym * vvel(iw) - dyT(iw) * tmp_vvel_ee &
503 0 : - cxm * uvel(iw) - dxT(iw) * tmp_uvel_se
504 0 : shearnw = -cyp * tmp_vvel_ee + dyT(iw) * vvel(iw) &
505 0 : - cxm * tmp_uvel_ee - dxT(iw) * tmp_uvel_ne
506 0 : shearsw = -cyp * tmp_vvel_ne + dyT(iw) * tmp_vvel_se &
507 0 : - cxp * tmp_uvel_ne + dxT(iw) * tmp_uvel_ee
508 0 : shearse = -cym * tmp_vvel_se - dyT(iw) * tmp_vvel_ne &
509 0 : - cxp * tmp_uvel_se + dxT(iw) * uvel(iw)
510 :
511 : ! Delta (in the denominator of zeta and eta)
512 0 : Deltane = sqrt(divune**2 + ecci * (tensionne**2 + shearne**2))
513 0 : Deltanw = sqrt(divunw**2 + ecci * (tensionnw**2 + shearnw**2))
514 0 : Deltasw = sqrt(divusw**2 + ecci * (tensionsw**2 + shearsw**2))
515 0 : Deltase = sqrt(divuse**2 + ecci * (tensionse**2 + shearse**2))
516 :
517 : !--------------------------------------------------------------
518 : ! on last subcycle, save quantities for mechanical
519 : ! redistribution
520 : !--------------------------------------------------------------
521 :
522 0 : divu(iw) = p25 * (divune + divunw + divuse + divusw) * tarear(iw)
523 0 : rdg_conv(iw) = -min(divu(iw), c0) ! TODO: Could move outside the entire kernel
524 0 : rdg_shear(iw) = p5 * (p25 * (Deltane + Deltanw + Deltase + Deltasw) * tarear(iw) - abs(divu(iw)))
525 :
526 : ! diagnostic only
527 : ! shear = sqrt(tension**2 + shearing**2)
528 0 : shear(iw) = p25 * tarear(iw) * sqrt((tensionne + tensionnw + tensionse + tensionsw)**2 &
529 0 : + (shearne + shearnw + shearse + shearsw)**2)
530 :
531 : !--------------------------------------------------------------
532 : ! replacement pressure/Delta (kg/s)
533 : ! save replacement pressure for principal stress calculation
534 : !--------------------------------------------------------------
535 :
536 0 : c0ne = strength(iw) / max(Deltane, DminTarea)
537 0 : c0nw = strength(iw) / max(Deltanw, DminTarea)
538 0 : c0sw = strength(iw) / max(Deltasw, DminTarea)
539 0 : c0se = strength(iw) / max(Deltase, DminTarea)
540 :
541 0 : c1ne = c0ne * arlx1i
542 0 : c1nw = c0nw * arlx1i
543 0 : c1sw = c0sw * arlx1i
544 0 : c1se = c0se * arlx1i
545 :
546 0 : c0ne = c1ne * ecci
547 0 : c0nw = c1nw * ecci
548 0 : c0sw = c1sw * ecci
549 0 : c0se = c1se * ecci
550 :
551 : !--------------------------------------------------------------
552 : ! the stresses (kg/s^2)
553 : ! (1) northeast, (2) northwest, (3) southwest, (4) southeast
554 : !--------------------------------------------------------------
555 :
556 0 : stressp_1(iw) = (stressp_1(iw) * (c1 - arlx1i * revp) &
557 0 : + c1ne * (divune * (c1 + Ktens) - Deltane * (c1 - Ktens))) * denom1
558 0 : stressp_2(iw) = (stressp_2(iw) * (c1 - arlx1i * revp) &
559 0 : + c1nw * (divunw * (c1 + Ktens) - Deltanw * (c1 - Ktens))) * denom1
560 0 : stressp_3(iw) = (stressp_3(iw) * (c1 - arlx1i * revp) &
561 0 : + c1sw * (divusw * (c1 + Ktens) - Deltasw * (c1 - Ktens))) * denom1
562 0 : stressp_4(iw) = (stressp_4(iw) * (c1 - arlx1i * revp) &
563 0 : + c1se * (divuse * (c1 + Ktens) - Deltase * (c1 - Ktens))) * denom1
564 :
565 0 : stressm_1(iw) = (stressm_1(iw) * (c1 - arlx1i * revp) + c0ne * tensionne * (c1 + Ktens)) * denom1
566 0 : stressm_2(iw) = (stressm_2(iw) * (c1 - arlx1i * revp) + c0nw * tensionnw * (c1 + Ktens)) * denom1
567 0 : stressm_3(iw) = (stressm_3(iw) * (c1 - arlx1i * revp) + c0sw * tensionsw * (c1 + Ktens)) * denom1
568 0 : stressm_4(iw) = (stressm_4(iw) * (c1 - arlx1i * revp) + c0se * tensionse * (c1 + Ktens)) * denom1
569 :
570 0 : stress12_1(iw) = (stress12_1(iw) * (c1 - arlx1i * revp) + c0ne * shearne * p5 * (c1 + Ktens)) * denom1
571 0 : stress12_2(iw) = (stress12_2(iw) * (c1 - arlx1i * revp) + c0nw * shearnw * p5 * (c1 + Ktens)) * denom1
572 0 : stress12_3(iw) = (stress12_3(iw) * (c1 - arlx1i * revp) + c0sw * shearsw * p5 * (c1 + Ktens)) * denom1
573 0 : stress12_4(iw) = (stress12_4(iw) * (c1 - arlx1i * revp) + c0se * shearse * p5 * (c1 + Ktens)) * denom1
574 :
575 : !--------------------------------------------------------------
576 : ! combinations of the stresses for the momentum equation
577 : ! (kg/s^2)
578 : !--------------------------------------------------------------
579 :
580 0 : ssigpn = stressp_1(iw) + stressp_2(iw)
581 0 : ssigps = stressp_3(iw) + stressp_4(iw)
582 0 : ssigpe = stressp_1(iw) + stressp_4(iw)
583 0 : ssigpw = stressp_2(iw) + stressp_3(iw)
584 0 : ssigp1 = (stressp_1(iw) + stressp_3(iw)) * p055
585 0 : ssigp2 = (stressp_2(iw) + stressp_4(iw)) * p055
586 :
587 0 : ssigmn = stressm_1(iw) + stressm_2(iw)
588 0 : ssigms = stressm_3(iw) + stressm_4(iw)
589 0 : ssigme = stressm_1(iw) + stressm_4(iw)
590 0 : ssigmw = stressm_2(iw) + stressm_3(iw)
591 0 : ssigm1 = (stressm_1(iw) + stressm_3(iw)) * p055
592 0 : ssigm2 = (stressm_2(iw) + stressm_4(iw)) * p055
593 :
594 0 : ssig12n = stress12_1(iw) + stress12_2(iw)
595 0 : ssig12s = stress12_3(iw) + stress12_4(iw)
596 0 : ssig12e = stress12_1(iw) + stress12_4(iw)
597 0 : ssig12w = stress12_2(iw) + stress12_3(iw)
598 0 : ssig121 = (stress12_1(iw) + stress12_3(iw)) * p111
599 0 : ssig122 = (stress12_2(iw) + stress12_4(iw)) * p111
600 :
601 0 : csigpne = p111 * stressp_1(iw) + ssigp2 + p027 * stressp_3(iw)
602 0 : csigpnw = p111 * stressp_2(iw) + ssigp1 + p027 * stressp_4(iw)
603 0 : csigpsw = p111 * stressp_3(iw) + ssigp2 + p027 * stressp_1(iw)
604 0 : csigpse = p111 * stressp_4(iw) + ssigp1 + p027 * stressp_2(iw)
605 :
606 0 : csigmne = p111 * stressm_1(iw) + ssigm2 + p027 * stressm_3(iw)
607 0 : csigmnw = p111 * stressm_2(iw) + ssigm1 + p027 * stressm_4(iw)
608 0 : csigmsw = p111 * stressm_3(iw) + ssigm2 + p027 * stressm_1(iw)
609 0 : csigmse = p111 * stressm_4(iw) + ssigm1 + p027 * stressm_2(iw)
610 :
611 0 : csig12ne = p222 * stress12_1(iw) + ssig122 + p055 * stress12_3(iw)
612 0 : csig12nw = p222 * stress12_2(iw) + ssig121 + p055 * stress12_4(iw)
613 0 : csig12sw = p222 * stress12_3(iw) + ssig122 + p055 * stress12_1(iw)
614 0 : csig12se = p222 * stress12_4(iw) + ssig121 + p055 * stress12_2(iw)
615 :
616 0 : str12ew = p5 * dxT(iw) * (p333 * ssig12e + p166 * ssig12w)
617 0 : str12we = p5 * dxT(iw) * (p333 * ssig12w + p166 * ssig12e)
618 0 : str12ns = p5 * dyT(iw) * (p333 * ssig12n + p166 * ssig12s)
619 0 : str12sn = p5 * dyT(iw) * (p333 * ssig12s + p166 * ssig12n)
620 :
621 : !--------------------------------------------------------------
622 : ! for dF/dx (u momentum)
623 : !--------------------------------------------------------------
624 :
625 0 : strp_tmp = p25 * dyT(iw) * (p333 * ssigpn + p166 * ssigps)
626 0 : strm_tmp = p25 * dyT(iw) * (p333 * ssigmn + p166 * ssigms)
627 :
628 : ! northeast (i,j)
629 0 : str1(iw) = -strp_tmp - strm_tmp - str12ew &
630 0 : + dxhy * (-csigpne + csigmne) + dyhx * csig12ne
631 :
632 : ! northwest (i+1,j)
633 0 : str2(iw) = strp_tmp + strm_tmp - str12we &
634 0 : + dxhy * (-csigpnw + csigmnw) + dyhx * csig12nw
635 :
636 0 : strp_tmp = p25 * dyT(iw) * (p333 * ssigps + p166 * ssigpn)
637 0 : strm_tmp = p25 * dyT(iw) * (p333 * ssigms + p166 * ssigmn)
638 :
639 : ! southeast (i,j+1)
640 0 : str3(iw) = -strp_tmp - strm_tmp + str12ew &
641 0 : + dxhy * (-csigpse + csigmse) + dyhx * csig12se
642 :
643 : ! southwest (i+1,j+1)
644 0 : str4(iw) = strp_tmp + strm_tmp + str12we &
645 0 : + dxhy * (-csigpsw + csigmsw) + dyhx * csig12sw
646 :
647 : !--------------------------------------------------------------
648 : ! for dF/dy (v momentum)
649 : !--------------------------------------------------------------
650 :
651 0 : strp_tmp = p25 * dxT(iw) * (p333 * ssigpe + p166 * ssigpw)
652 0 : strm_tmp = p25 * dxT(iw) * (p333 * ssigme + p166 * ssigmw)
653 :
654 : ! northeast (i,j)
655 0 : str5(iw) = -strp_tmp + strm_tmp - str12ns &
656 0 : - dyhx * (csigpne + csigmne) + dxhy * csig12ne
657 :
658 : ! southeast (i,j+1)
659 0 : str6(iw) = strp_tmp - strm_tmp - str12sn &
660 0 : - dyhx * (csigpse + csigmse) + dxhy * csig12se
661 :
662 0 : strp_tmp = p25 * dxT(iw) * (p333 * ssigpw + p166 * ssigpe)
663 0 : strm_tmp = p25 * dxT(iw) * (p333 * ssigmw + p166 * ssigme)
664 :
665 : ! northwest (i+1,j)
666 0 : str7(iw) = -strp_tmp + strm_tmp + str12ns &
667 0 : - dyhx * (csigpnw + csigmnw) + dxhy * csig12nw
668 :
669 : ! southwest (i+1,j+1)
670 0 : str8(iw) = strp_tmp - strm_tmp + str12sn &
671 0 : - dyhx * (csigpsw + csigmsw) + dxhy * csig12sw
672 :
673 : end do
674 : #ifdef _OPENACC
675 : !$acc end parallel
676 : #endif
677 :
678 0 : end subroutine stress_last
679 :
680 : !=======================================================================
681 :
682 0 : subroutine stepu_iter(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, &
683 : forcex, forcey, umassdti, fm, uarear, Tbu, uvel_init, vvel_init, & ! LCOV_EXCL_LINE
684 : uvel, vvel, str1, str2, str3, str4, str5, str6, str7, str8, nw, & ! LCOV_EXCL_LINE
685 0 : sw, sse, skipucell)
686 :
687 : use ice_kinds_mod
688 : use ice_constants, only : c0, c1
689 : use ice_dyn_shared, only : brlx, revp, u0, cosw, sinw
690 :
691 : implicit none
692 :
693 : integer(kind=int_kind), intent(in) :: NA_len, lb, ub
694 : real(kind=dbl_kind), intent(in) :: rhow
695 : logical(kind=log_kind), intent(in), dimension(:) :: skipucell
696 : integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
697 : nw, sw, sse
698 : real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
699 : uvel_init, vvel_init, aiu, forcex, forcey, umassdti, Tbu, & ! LCOV_EXCL_LINE
700 : uocn, vocn, fm, uarear, Cw, str1, str2, str3, str4, str5, & ! LCOV_EXCL_LINE
701 : str6, str7, str8
702 : real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
703 : uvel, vvel
704 :
705 : ! local variables
706 :
707 : integer(kind=int_kind) :: iw, il, iu
708 0 : real(kind=dbl_kind) :: uold, vold, vrel, cca, ccb, ab2, cc1, &
709 : cc2, taux, tauy, Cb, tmp_str2_nw, tmp_str3_sse, tmp_str4_sw, & ! LCOV_EXCL_LINE
710 : tmp_str6_sse, tmp_str7_nw, tmp_str8_sw, waterx, watery, & ! LCOV_EXCL_LINE
711 0 : tmp_strintx, tmp_strinty
712 :
713 : character(len=*), parameter :: subname = '(stepu_iter)'
714 :
715 : #ifdef _OPENACC
716 : !$acc parallel &
717 : !$acc present(Cw, aiu, uocn, vocn, forcex, forcey, umassdti, fm, & ! LCOV_EXCL_LINE
718 : !$acc uarear, Tbu, uvel_init, vvel_init, nw, sw, sse, skipucell, & ! LCOV_EXCL_LINE
719 : !$acc str1, str2, str3, str4, str5, str6, str7, str8, uvel, & ! LCOV_EXCL_LINE
720 : !$acc vvel)
721 : !$acc loop
722 : do iw = 1, NA_len
723 : #else
724 0 : call domp_get_domain(lb, ub, il, iu)
725 0 : do iw = il, iu
726 : #endif
727 :
728 0 : if (skipucell(iw)) cycle
729 :
730 0 : uold = uvel(iw)
731 0 : vold = vvel(iw)
732 :
733 0 : vrel = aiu(iw) * rhow * Cw(iw) * sqrt((uocn(iw) - uold)**2 + (vocn(iw) - vold)**2)
734 :
735 0 : waterx = uocn(iw) * cosw - vocn(iw) * sinw * sign(c1, fm(iw))
736 0 : watery = vocn(iw) * cosw + uocn(iw) * sinw * sign(c1, fm(iw))
737 :
738 0 : taux = vrel * waterx
739 0 : tauy = vrel * watery
740 :
741 0 : Cb = Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
742 :
743 0 : cca = (brlx + revp) * umassdti(iw) + vrel * cosw + Cb
744 0 : ccb = fm(iw) + sign(c1, fm(iw)) * vrel * sinw
745 :
746 0 : ab2 = cca**2 + ccb**2
747 :
748 0 : tmp_str2_nw = str2(nw(iw))
749 0 : tmp_str3_sse = str3(sse(iw))
750 0 : tmp_str4_sw = str4(sw(iw))
751 0 : tmp_str6_sse = str6(sse(iw))
752 0 : tmp_str7_nw = str7(nw(iw))
753 0 : tmp_str8_sw = str8(sw(iw))
754 :
755 0 : tmp_strintx = uarear(iw) * (str1(iw) + tmp_str2_nw + tmp_str3_sse + tmp_str4_sw)
756 0 : tmp_strinty = uarear(iw) * (str5(iw) + tmp_str6_sse + tmp_str7_nw + tmp_str8_sw)
757 :
758 0 : cc1 = tmp_strintx + forcex(iw) + taux &
759 0 : + umassdti(iw) * (brlx * uold + revp * uvel_init(iw))
760 0 : cc2 = tmp_strinty + forcey(iw) + tauy &
761 0 : + umassdti(iw) * (brlx * vold + revp * vvel_init(iw))
762 :
763 0 : uvel(iw) = (cca * cc1 + ccb * cc2) / ab2
764 0 : vvel(iw) = (cca * cc2 - ccb * cc1) / ab2
765 :
766 : end do
767 : #ifdef _OPENACC
768 : !$acc end parallel
769 : #endif
770 :
771 0 : end subroutine stepu_iter
772 :
773 : !=======================================================================
774 :
775 0 : subroutine stepu_last(NA_len, rhow, lb, ub, Cw, aiu, uocn, vocn, &
776 : forcex, forcey, umassdti, fm, uarear, Tbu, uvel_init, vvel_init, & ! LCOV_EXCL_LINE
777 : uvel, vvel, str1, str2, str3, str4, str5, str6, str7, str8, nw, & ! LCOV_EXCL_LINE
778 0 : sw, sse, skipucell, strintx, strinty, taubx, tauby)
779 :
780 : use ice_kinds_mod
781 : use ice_constants, only : c0, c1
782 : use ice_dyn_shared, only : brlx, revp, u0, cosw, sinw
783 :
784 : implicit none
785 :
786 : integer(kind=int_kind), intent(in) :: NA_len, lb, ub
787 : real(kind=dbl_kind), intent(in) :: rhow
788 : logical(kind=log_kind), intent(in), dimension(:) :: skipucell
789 : integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
790 : nw, sw, sse
791 : real(kind=dbl_kind), dimension(:), intent(in), contiguous :: &
792 : uvel_init, vvel_init, aiu, forcex, forcey, umassdti, Tbu, & ! LCOV_EXCL_LINE
793 : uocn, vocn, fm, uarear, Cw, str1, str2, str3, str4, str5, & ! LCOV_EXCL_LINE
794 : str6, str7, str8
795 : real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
796 : uvel, vvel, strintx, strinty, taubx, tauby
797 :
798 : ! local variables
799 :
800 : integer(kind=int_kind) :: iw, il, iu
801 0 : real(kind=dbl_kind) :: uold, vold, vrel, cca, ccb, ab2, cc1, &
802 : cc2, taux, tauy, Cb, tmp_str2_nw, tmp_str3_sse, tmp_str4_sw, & ! LCOV_EXCL_LINE
803 0 : tmp_str6_sse, tmp_str7_nw, tmp_str8_sw, waterx, watery
804 :
805 : character(len=*), parameter :: subname = '(stepu_last)'
806 :
807 : #ifdef _OPENACC
808 : !$acc parallel &
809 : !$acc present(Cw, aiu, uocn, vocn, forcex, forcey, umassdti, fm, & ! LCOV_EXCL_LINE
810 : !$acc uarear, Tbu, uvel_init, vvel_init, nw, sw, sse, skipucell, & ! LCOV_EXCL_LINE
811 : !$acc str1, str2, str3, str4, str5, str6, str7, str8, uvel, & ! LCOV_EXCL_LINE
812 : !$acc vvel, strintx, strinty, taubx, tauby)
813 : !$acc loop
814 : do iw = 1, NA_len
815 : #else
816 0 : call domp_get_domain(lb, ub, il, iu)
817 0 : do iw = il, iu
818 : #endif
819 :
820 0 : if (skipucell(iw)) cycle
821 :
822 0 : uold = uvel(iw)
823 0 : vold = vvel(iw)
824 :
825 0 : vrel = aiu(iw) * rhow * Cw(iw) * sqrt((uocn(iw) - uold)**2 + (vocn(iw) - vold)**2)
826 :
827 0 : waterx = uocn(iw) * cosw - vocn(iw) * sinw * sign(c1, fm(iw))
828 0 : watery = vocn(iw) * cosw + uocn(iw) * sinw * sign(c1, fm(iw))
829 :
830 0 : taux = vrel * waterx
831 0 : tauy = vrel * watery
832 :
833 0 : Cb = Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
834 :
835 0 : cca = (brlx + revp) * umassdti(iw) + vrel * cosw + Cb
836 0 : ccb = fm(iw) + sign(c1, fm(iw)) * vrel * sinw
837 :
838 0 : ab2 = cca**2 + ccb**2
839 :
840 0 : tmp_str2_nw = str2(nw(iw))
841 0 : tmp_str3_sse = str3(sse(iw))
842 0 : tmp_str4_sw = str4(sw(iw))
843 0 : tmp_str6_sse = str6(sse(iw))
844 0 : tmp_str7_nw = str7(nw(iw))
845 0 : tmp_str8_sw = str8(sw(iw))
846 :
847 0 : strintx(iw) = uarear(iw) * (str1(iw) + tmp_str2_nw + tmp_str3_sse + tmp_str4_sw)
848 0 : strinty(iw) = uarear(iw) * (str5(iw) + tmp_str6_sse + tmp_str7_nw + tmp_str8_sw)
849 :
850 0 : cc1 = strintx(iw) + forcex(iw) + taux &
851 0 : + umassdti(iw) * (brlx * uold + revp * uvel_init(iw))
852 0 : cc2 = strinty(iw) + forcey(iw) + tauy &
853 0 : + umassdti(iw) * (brlx * vold + revp * vvel_init(iw))
854 :
855 0 : uvel(iw) = (cca * cc1 + ccb * cc2) / ab2
856 0 : vvel(iw) = (cca * cc2 - ccb * cc1) / ab2
857 :
858 : ! calculate seabed stress component for outputs
859 0 : taubx(iw) = -uvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
860 0 : tauby(iw) = -vvel(iw) * Tbu(iw) / (sqrt(uold**2 + vold**2) + u0)
861 :
862 : end do
863 : #ifdef _OPENACC
864 : !$acc end parallel
865 : #endif
866 :
867 0 : end subroutine stepu_last
868 :
869 : !=======================================================================
870 :
871 0 : subroutine evp1d_halo_update(NAVEL_len, lb, ub, uvel, vvel, &
872 0 : halo_parent)
873 :
874 : use ice_kinds_mod
875 :
876 : implicit none
877 :
878 : integer(kind=int_kind), intent(in) :: NAVEL_len, lb, ub
879 : integer(kind=int_kind), dimension(:), intent(in), contiguous :: &
880 : halo_parent
881 : real(kind=dbl_kind), dimension(:), intent(inout), contiguous :: &
882 : uvel, vvel
883 :
884 : ! local variables
885 :
886 : integer (kind=int_kind) :: iw, il, iu
887 :
888 : character(len=*), parameter :: subname = '(evp1d_halo_update)'
889 :
890 : #ifdef _OPENACC
891 : !$acc parallel &
892 : !$acc present(uvel, vvel)
893 : !$acc loop
894 : do iw = 1, NAVEL_len
895 : if (halo_parent(iw) == 0) cycle
896 : uvel(iw) = uvel(halo_parent(iw))
897 : vvel(iw) = vvel(halo_parent(iw))
898 : end do
899 : !$acc end parallel
900 : #else
901 0 : call domp_get_domain(lb, ub, il, iu)
902 0 : do iw = il, iu
903 0 : if (halo_parent(iw) == 0) cycle
904 0 : uvel(iw) = uvel(halo_parent(iw))
905 0 : vvel(iw) = vvel(halo_parent(iw))
906 : end do
907 0 : call domp_get_domain(ub + 1, NAVEL_len, il, iu)
908 0 : do iw = il, iu
909 0 : if (halo_parent(iw) == 0) cycle
910 0 : uvel(iw) = uvel(halo_parent(iw))
911 0 : vvel(iw) = vvel(halo_parent(iw))
912 : end do
913 : #endif
914 :
915 0 : end subroutine evp1d_halo_update
916 :
917 : !=======================================================================
918 :
919 0 : subroutine alloc1d(na)
920 :
921 : implicit none
922 :
923 : integer(kind=int_kind), intent(in) :: na
924 :
925 : ! local variables
926 :
927 : integer(kind=int_kind) :: ierr
928 :
929 : character(len=*), parameter :: subname = '(alloc1d)'
930 :
931 : allocate( &
932 : ! helper indices for neighbours
933 : indj(1:na), indi(1:na), ee(1:na), ne(1:na), se(1:na), &
934 : nw(1:na), sw(1:na), sse(1:na), skipucell(1:na), & ! LCOV_EXCL_LINE
935 : skiptcell(1:na), & ! LCOV_EXCL_LINE
936 : ! grid distances and their "-1 neighbours"
937 : HTE(1:na), HTN(1:na), HTEm1(1:na), HTNm1(1:na), &
938 : ! T cells
939 : strength(1:na), dxT(1:na), dyT(1:na), tarear(1:na), &
940 : stressp_1(1:na), stressp_2(1:na), stressp_3(1:na), & ! LCOV_EXCL_LINE
941 : stressp_4(1:na), stressm_1(1:na), stressm_2(1:na), & ! LCOV_EXCL_LINE
942 : stressm_3(1:na), stressm_4(1:na), stress12_1(1:na), & ! LCOV_EXCL_LINE
943 : stress12_2(1:na), stress12_3(1:na), stress12_4(1:na), & ! LCOV_EXCL_LINE
944 : divu(1:na), rdg_conv(1:na), rdg_shear(1:na), shear(1:na), & ! LCOV_EXCL_LINE
945 : ! U cells
946 : cdn_ocn(1:na), aiu(1:na), uocn(1:na), vocn(1:na), &
947 : forcex(1:na), forcey(1:na), Tbu(1:na), umassdti(1:na), & ! LCOV_EXCL_LINE
948 : fm(1:na), uarear(1:na), strintx(1:na), strinty(1:na), & ! LCOV_EXCL_LINE
949 : uvel_init(1:na), vvel_init(1:na), taubx(1:na), tauby(1:na), & ! LCOV_EXCL_LINE
950 : ! error handling
951 0 : stat=ierr)
952 :
953 0 : if (ierr /= 0) call abort_ice(subname &
954 0 : // ' ERROR: could not allocate 1D arrays')
955 :
956 0 : end subroutine alloc1d
957 :
958 : !=======================================================================
959 :
960 0 : subroutine alloc1d_navel(navel)
961 :
962 : implicit none
963 :
964 : integer(kind=int_kind), intent(in) :: navel
965 :
966 : ! local variables
967 :
968 : integer(kind=int_kind) :: ierr
969 :
970 : character(len=*), parameter :: subname = '(alloc1d_navel)'
971 :
972 : allocate(uvel(1:navel), vvel(1:navel), indij(1:navel), &
973 : halo_parent(1:navel), str1(1:navel), str2(1:navel), & ! LCOV_EXCL_LINE
974 : str3(1:navel), str4(1:navel), str5(1:navel), str6(1:navel), & ! LCOV_EXCL_LINE
975 0 : str7(1:navel), str8(1:navel), stat=ierr)
976 :
977 0 : if (ierr /= 0) call abort_ice(subname &
978 0 : // ' ERROR: could not allocate 1D arrays')
979 :
980 0 : end subroutine alloc1d_navel
981 :
982 : !=======================================================================
983 :
984 0 : subroutine dealloc1d
985 :
986 : implicit none
987 :
988 : ! local variables
989 :
990 : integer(kind=int_kind) :: ierr
991 :
992 : character(len=*), parameter :: subname = '(dealloc1d)'
993 :
994 : deallocate( &
995 : ! helper indices for neighbours
996 : indj, indi, ee, ne, se, nw, sw, sse, skipucell, skiptcell, &
997 : ! grid distances and their "-1 neighbours"
998 : HTE, HTN, HTEm1, HTNm1, &
999 : ! T cells
1000 : strength, dxT, dyT, tarear, stressp_1, stressp_2, stressp_3, &
1001 : stressp_4, stressm_1, stressm_2, stressm_3, stressm_4, & ! LCOV_EXCL_LINE
1002 : stress12_1, stress12_2, stress12_3, stress12_4, str1, str2, & ! LCOV_EXCL_LINE
1003 : str3, str4, str5, str6, str7, str8, divu, rdg_conv, & ! LCOV_EXCL_LINE
1004 : rdg_shear, shear, & ! LCOV_EXCL_LINE
1005 : ! U cells
1006 : cdn_ocn, aiu, uocn, vocn, forcex, forcey, Tbu, umassdti, fm, &
1007 : uarear, strintx, strinty, uvel_init, vvel_init, taubx, tauby, & ! LCOV_EXCL_LINE
1008 : uvel, vvel, indij, halo_parent, & ! LCOV_EXCL_LINE
1009 : ! error handling
1010 0 : stat=ierr)
1011 :
1012 0 : if (ierr /= 0) call abort_ice(subname &
1013 0 : // ' ERROR: could not deallocate 1D arrays')
1014 :
1015 0 : end subroutine dealloc1d
1016 :
1017 : !=======================================================================
1018 :
1019 0 : subroutine ice_dyn_evp_1d_copyin(nx, ny, nblk, nx_glob, ny_glob, &
1020 : I_iceTmask, I_iceUmask, I_cdn_ocn, I_aiu, I_uocn, I_vocn, & ! LCOV_EXCL_LINE
1021 : I_forcex, I_forcey, I_Tbu, I_umassdti, I_fm, I_uarear, I_tarear, & ! LCOV_EXCL_LINE
1022 : I_strintx, I_strinty, I_uvel_init, I_vvel_init, I_strength, & ! LCOV_EXCL_LINE
1023 : I_uvel, I_vvel, I_dxT, I_dyT, I_stressp_1, I_stressp_2, & ! LCOV_EXCL_LINE
1024 : I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, & ! LCOV_EXCL_LINE
1025 : I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, & ! LCOV_EXCL_LINE
1026 0 : I_stress12_4)
1027 :
1028 : use ice_gather_scatter, only : gather_global_ext
1029 : use ice_domain, only : distrb_info
1030 : use ice_communicate, only : my_task, master_task
1031 : use ice_grid, only : G_HTE, G_HTN
1032 : use ice_constants, only : c0
1033 :
1034 : implicit none
1035 :
1036 : integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob
1037 : logical(kind=log_kind), dimension(nx, ny, nblk), intent(in) :: &
1038 : I_iceTmask, I_iceUmask
1039 : real(kind=dbl_kind), dimension(nx, ny, nblk), intent(in) :: &
1040 : I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, I_Tbu, & ! LCOV_EXCL_LINE
1041 : I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, I_strinty, & ! LCOV_EXCL_LINE
1042 : I_uvel_init, I_vvel_init, I_strength, I_uvel, I_vvel, I_dxT, & ! LCOV_EXCL_LINE
1043 : I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, I_stressp_4, & ! LCOV_EXCL_LINE
1044 : I_stressm_1, I_stressm_2, I_stressm_3, I_stressm_4, & ! LCOV_EXCL_LINE
1045 : I_stress12_1, I_stress12_2, I_stress12_3, I_stress12_4
1046 :
1047 : ! local variables
1048 :
1049 : logical(kind=log_kind), dimension(nx_glob, ny_glob) :: &
1050 0 : G_iceTmask, G_iceUmask
1051 : real(kind=dbl_kind), dimension(nx_glob, ny_glob) :: &
1052 : G_cdn_ocn, G_aiu, G_uocn, G_vocn, G_forcex, G_forcey, G_Tbu, & ! LCOV_EXCL_LINE
1053 : G_umassdti, G_fm, G_uarear, G_tarear, G_strintx, G_strinty, & ! LCOV_EXCL_LINE
1054 : G_uvel_init, G_vvel_init, G_strength, G_uvel, G_vvel, G_dxT, & ! LCOV_EXCL_LINE
1055 : G_dyT, G_stressp_1, G_stressp_2, G_stressp_3, G_stressp_4, & ! LCOV_EXCL_LINE
1056 : G_stressm_1, G_stressm_2, G_stressm_3, G_stressm_4, & ! LCOV_EXCL_LINE
1057 0 : G_stress12_1, G_stress12_2, G_stress12_3, G_stress12_4
1058 :
1059 : character(len=*), parameter :: &
1060 : subname = '(ice_dyn_evp_1d_copyin)'
1061 :
1062 0 : call gather_global_ext(G_iceTmask, I_iceTmask, master_task, distrb_info )
1063 0 : call gather_global_ext(G_iceUmask, I_iceUmask, master_task, distrb_info )
1064 0 : call gather_global_ext(G_cdn_ocn, I_cdn_ocn, master_task, distrb_info )
1065 0 : call gather_global_ext(G_aiu, I_aiu, master_task, distrb_info )
1066 0 : call gather_global_ext(G_uocn, I_uocn, master_task, distrb_info )
1067 0 : call gather_global_ext(G_vocn, I_vocn, master_task, distrb_info )
1068 0 : call gather_global_ext(G_forcex, I_forcex, master_task, distrb_info )
1069 0 : call gather_global_ext(G_forcey, I_forcey, master_task, distrb_info )
1070 0 : call gather_global_ext(G_Tbu, I_Tbu, master_task, distrb_info )
1071 0 : call gather_global_ext(G_umassdti, I_umassdti, master_task, distrb_info )
1072 0 : call gather_global_ext(G_fm, I_fm, master_task, distrb_info )
1073 0 : call gather_global_ext(G_uarear, I_uarear, master_task, distrb_info )
1074 0 : call gather_global_ext(G_tarear, I_tarear, master_task, distrb_info )
1075 0 : call gather_global_ext(G_strintx, I_strintx, master_task, distrb_info )
1076 0 : call gather_global_ext(G_strinty, I_strinty, master_task, distrb_info )
1077 0 : call gather_global_ext(G_uvel_init, I_uvel_init, master_task, distrb_info )
1078 0 : call gather_global_ext(G_vvel_init, I_vvel_init, master_task, distrb_info )
1079 0 : call gather_global_ext(G_strength, I_strength, master_task, distrb_info )
1080 0 : call gather_global_ext(G_uvel, I_uvel, master_task, distrb_info, c0)
1081 0 : call gather_global_ext(G_vvel, I_vvel, master_task, distrb_info, c0)
1082 0 : call gather_global_ext(G_dxT, I_dxT, master_task, distrb_info )
1083 0 : call gather_global_ext(G_dyT, I_dyT, master_task, distrb_info )
1084 0 : call gather_global_ext(G_stressp_1, I_stressp_1, master_task, distrb_info )
1085 0 : call gather_global_ext(G_stressp_2, I_stressp_2, master_task, distrb_info )
1086 0 : call gather_global_ext(G_stressp_3, I_stressp_3, master_task, distrb_info )
1087 0 : call gather_global_ext(G_stressp_4, I_stressp_4, master_task, distrb_info )
1088 0 : call gather_global_ext(G_stressm_1, I_stressm_1, master_task, distrb_info )
1089 0 : call gather_global_ext(G_stressm_2, I_stressm_2, master_task, distrb_info )
1090 0 : call gather_global_ext(G_stressm_3, I_stressm_3, master_task, distrb_info )
1091 0 : call gather_global_ext(G_stressm_4, I_stressm_4, master_task, distrb_info )
1092 0 : call gather_global_ext(G_stress12_1, I_stress12_1, master_task, distrb_info )
1093 0 : call gather_global_ext(G_stress12_2, I_stress12_2, master_task, distrb_info )
1094 0 : call gather_global_ext(G_stress12_3, I_stress12_3, master_task, distrb_info )
1095 0 : call gather_global_ext(G_stress12_4, I_stress12_4, master_task, distrb_info )
1096 :
1097 : ! all calculations id done on master task
1098 0 : if (my_task == master_task) then
1099 : ! find number of active points and allocate 1D vectors
1100 0 : call calc_na(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask)
1101 0 : call alloc1d(NA_len)
1102 0 : call calc_2d_indices(nx_glob, ny_glob, NA_len, G_iceTmask, G_iceUmask)
1103 0 : call calc_navel(nx_glob, ny_glob, NA_len, NAVEL_len)
1104 0 : call alloc1d_navel(NAVEL_len)
1105 : ! initialize OpenMP. FIXME: ought to be called from main
1106 0 : call domp_init()
1107 : !$OMP PARALLEL DEFAULT(shared)
1108 0 : call numainit(1, NA_len, NAVEL_len)
1109 : !$OMP END PARALLEL
1110 : ! map 2D arrays to 1D arrays
1111 : call convert_2d_1d(nx_glob, ny_glob, NA_len, NAVEL_len, &
1112 : G_HTE, G_HTN, G_cdn_ocn, G_aiu, G_uocn, G_vocn, G_forcex, & ! LCOV_EXCL_LINE
1113 : G_forcey, G_Tbu, G_umassdti, G_fm, G_uarear, G_tarear, & ! LCOV_EXCL_LINE
1114 : G_strintx, G_strinty, G_uvel_init, G_vvel_init, & ! LCOV_EXCL_LINE
1115 : G_strength, G_uvel, G_vvel, G_dxT, G_dyT, G_stressp_1, & ! LCOV_EXCL_LINE
1116 : G_stressp_2, G_stressp_3, G_stressp_4, G_stressm_1, & ! LCOV_EXCL_LINE
1117 : G_stressm_2, G_stressm_3, G_stressm_4, G_stress12_1, & ! LCOV_EXCL_LINE
1118 0 : G_stress12_2, G_stress12_3, G_stress12_4)
1119 0 : call calc_halo_parent(nx_glob, ny_glob, NA_len, NAVEL_len, G_iceTmask)
1120 : end if
1121 :
1122 0 : end subroutine ice_dyn_evp_1d_copyin
1123 :
1124 : !=======================================================================
1125 :
1126 0 : subroutine ice_dyn_evp_1d_copyout(nx, ny, nblk, nx_glob, ny_glob, &
1127 : I_uvel, I_vvel, I_strintx, I_strinty, I_stressp_1, I_stressp_2, & ! LCOV_EXCL_LINE
1128 : I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, & ! LCOV_EXCL_LINE
1129 : I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, & ! LCOV_EXCL_LINE
1130 : I_stress12_4, I_divu, I_rdg_conv, I_rdg_shear, I_shear, I_taubx, & ! LCOV_EXCL_LINE
1131 0 : I_tauby)
1132 :
1133 : use ice_constants, only : c0
1134 : use ice_gather_scatter, only : scatter_global_ext
1135 : use ice_domain, only : distrb_info
1136 : use ice_communicate, only : my_task, master_task
1137 :
1138 : implicit none
1139 :
1140 : integer(int_kind), intent(in) :: nx, ny, nblk, nx_glob, ny_glob
1141 : real(dbl_kind), dimension(nx, ny, nblk), intent(out) :: I_uvel, &
1142 : I_vvel, I_strintx, I_strinty, I_stressp_1, I_stressp_2, & ! LCOV_EXCL_LINE
1143 : I_stressp_3, I_stressp_4, I_stressm_1, I_stressm_2, & ! LCOV_EXCL_LINE
1144 : I_stressm_3, I_stressm_4, I_stress12_1, I_stress12_2, & ! LCOV_EXCL_LINE
1145 : I_stress12_3, I_stress12_4, I_divu, I_rdg_conv, I_rdg_shear, & ! LCOV_EXCL_LINE
1146 : I_shear, I_taubx, I_tauby
1147 :
1148 : ! local variables
1149 :
1150 : integer(int_kind) :: iw, lo, up, j, i
1151 0 : real(dbl_kind), dimension(nx_glob, ny_glob) :: G_uvel, G_vvel, &
1152 : G_strintx, G_strinty, G_stressp_1, G_stressp_2, G_stressp_3, & ! LCOV_EXCL_LINE
1153 : G_stressp_4, G_stressm_1, G_stressm_2, G_stressm_3, & ! LCOV_EXCL_LINE
1154 : G_stressm_4, G_stress12_1, G_stress12_2, G_stress12_3, & ! LCOV_EXCL_LINE
1155 : G_stress12_4, G_divu, G_rdg_conv, G_rdg_shear, G_shear, & ! LCOV_EXCL_LINE
1156 0 : G_taubx, G_tauby
1157 :
1158 : character(len=*), parameter :: &
1159 : subname = '(ice_dyn_evp_1d_copyout)'
1160 :
1161 : ! remap 1D arrays into 2D arrays
1162 0 : if (my_task == master_task) then
1163 :
1164 0 : G_uvel = c0
1165 0 : G_vvel = c0
1166 0 : G_strintx = c0
1167 0 : G_strinty = c0
1168 0 : G_stressp_1 = c0
1169 0 : G_stressp_2 = c0
1170 0 : G_stressp_3 = c0
1171 0 : G_stressp_4 = c0
1172 0 : G_stressm_1 = c0
1173 0 : G_stressm_2 = c0
1174 0 : G_stressm_3 = c0
1175 0 : G_stressm_4 = c0
1176 0 : G_stress12_1 = c0
1177 0 : G_stress12_2 = c0
1178 0 : G_stress12_3 = c0
1179 0 : G_stress12_4 = c0
1180 0 : G_divu = c0
1181 0 : G_rdg_conv = c0
1182 0 : G_rdg_shear = c0
1183 0 : G_shear = c0
1184 0 : G_taubx = c0
1185 0 : G_tauby = c0
1186 :
1187 0 : !$OMP PARALLEL PRIVATE(iw, lo, up, j, i)
1188 0 : call domp_get_domain(1, NA_len, lo, up)
1189 0 : do iw = lo, up
1190 : ! get 2D indices
1191 0 : i = indi(iw)
1192 0 : j = indj(iw)
1193 : ! remap
1194 0 : G_strintx(i, j) = strintx(iw)
1195 0 : G_strinty(i, j) = strinty(iw)
1196 0 : G_stressp_1(i, j) = stressp_1(iw)
1197 0 : G_stressp_2(i, j) = stressp_2(iw)
1198 0 : G_stressp_3(i, j) = stressp_3(iw)
1199 0 : G_stressp_4(i, j) = stressp_4(iw)
1200 0 : G_stressm_1(i, j) = stressm_1(iw)
1201 0 : G_stressm_2(i, j) = stressm_2(iw)
1202 0 : G_stressm_3(i, j) = stressm_3(iw)
1203 0 : G_stressm_4(i, j) = stressm_4(iw)
1204 0 : G_stress12_1(i, j) = stress12_1(iw)
1205 0 : G_stress12_2(i, j) = stress12_2(iw)
1206 0 : G_stress12_3(i, j) = stress12_3(iw)
1207 0 : G_stress12_4(i, j) = stress12_4(iw)
1208 0 : G_divu(i, j) = divu(iw)
1209 0 : G_rdg_conv(i, j) = rdg_conv(iw)
1210 0 : G_rdg_shear(i, j) = rdg_shear(iw)
1211 0 : G_shear(i, j) = shear(iw)
1212 0 : G_taubx(i, j) = taubx(iw)
1213 0 : G_tauby(i, j) = tauby(iw)
1214 0 : G_uvel(i, j) = uvel(iw)
1215 0 : G_vvel(i, j) = vvel(iw)
1216 : end do
1217 0 : call domp_get_domain(NA_len + 1, NAVEL_len, lo, up)
1218 0 : do iw = lo, up
1219 : ! get 2D indices
1220 0 : j = int((indij(iw) - 1) / (nx_glob)) + 1
1221 0 : i = indij(iw) - (j - 1) * nx_glob
1222 : ! remap
1223 0 : G_uvel(i, j) = uvel(iw)
1224 0 : G_vvel(i, j) = vvel(iw)
1225 : end do
1226 : !$OMP END PARALLEL
1227 :
1228 0 : call dealloc1d()
1229 :
1230 : end if
1231 :
1232 : ! scatter data on all tasks
1233 0 : call scatter_global_ext(I_uvel, G_uvel, master_task, distrb_info)
1234 0 : call scatter_global_ext(I_vvel, G_vvel, master_task, distrb_info)
1235 0 : call scatter_global_ext(I_strintx, G_strintx, master_task, distrb_info)
1236 0 : call scatter_global_ext(I_strinty, G_strinty, master_task, distrb_info)
1237 0 : call scatter_global_ext(I_stressp_1, G_stressp_1, master_task, distrb_info)
1238 0 : call scatter_global_ext(I_stressp_2, G_stressp_2, master_task, distrb_info)
1239 0 : call scatter_global_ext(I_stressp_3, G_stressp_3, master_task, distrb_info)
1240 0 : call scatter_global_ext(I_stressp_4, G_stressp_4, master_task, distrb_info)
1241 0 : call scatter_global_ext(I_stressm_1, G_stressm_1, master_task, distrb_info)
1242 0 : call scatter_global_ext(I_stressm_2, G_stressm_2, master_task, distrb_info)
1243 0 : call scatter_global_ext(I_stressm_3, G_stressm_3, master_task, distrb_info)
1244 0 : call scatter_global_ext(I_stressm_4, G_stressm_4, master_task, distrb_info)
1245 0 : call scatter_global_ext(I_stress12_1, G_stress12_1, master_task, distrb_info)
1246 0 : call scatter_global_ext(I_stress12_2, G_stress12_2, master_task, distrb_info)
1247 0 : call scatter_global_ext(I_stress12_3, G_stress12_3, master_task, distrb_info)
1248 0 : call scatter_global_ext(I_stress12_4, G_stress12_4, master_task, distrb_info)
1249 0 : call scatter_global_ext(I_divu, G_divu, master_task, distrb_info)
1250 0 : call scatter_global_ext(I_rdg_conv, G_rdg_conv, master_task, distrb_info)
1251 0 : call scatter_global_ext(I_rdg_shear, G_rdg_shear, master_task, distrb_info)
1252 0 : call scatter_global_ext(I_shear, G_shear, master_task, distrb_info)
1253 0 : call scatter_global_ext(I_taubx, G_taubx, master_task, distrb_info)
1254 0 : call scatter_global_ext(I_tauby, G_tauby, master_task, distrb_info)
1255 :
1256 0 : end subroutine ice_dyn_evp_1d_copyout
1257 :
1258 : !=======================================================================
1259 :
1260 0 : subroutine ice_dyn_evp_1d_kernel
1261 :
1262 : use ice_constants, only : c0
1263 : use ice_dyn_shared, only : ndte
1264 : use ice_communicate, only : my_task, master_task
1265 :
1266 : implicit none
1267 :
1268 : ! local variables
1269 :
1270 0 : real(kind=dbl_kind) :: rhow
1271 : integer(kind=int_kind) :: ksub
1272 :
1273 : character(len=*), parameter :: &
1274 : subname = '(ice_dyn_evp_1d_kernel)'
1275 :
1276 : ! all calculations is done on master task
1277 0 : if (my_task == master_task) then
1278 :
1279 : ! read constants
1280 0 : call icepack_query_parameters(rhow_out = rhow)
1281 0 : call icepack_warnings_flush(nu_diag)
1282 0 : if (icepack_warnings_aborted()) then
1283 : call abort_ice(error_message=subname, file=__FILE__, &
1284 0 : line=__LINE__)
1285 : end if
1286 :
1287 0 : if (ndte < 2) call abort_ice(subname &
1288 0 : // ' ERROR: ndte must be 2 or higher for this kernel')
1289 :
1290 : ! tcraig, turn off the OMP directives here, Jan, 2022
1291 : ! This produces non bit-for-bit results with different thread counts.
1292 : ! Seems like there isn't an opportunity for safe threading here ???
1293 : !$XXXOMP PARALLEL PRIVATE(ksub)
1294 0 : do ksub = 1, ndte - 1
1295 : call evp1d_stress(NA_len, ee, ne, se, 1, NA_len, uvel, &
1296 : vvel, dxT, dyT, hte, htn, htem1, htnm1, strength, & ! LCOV_EXCL_LINE
1297 : stressp_1, stressp_2, stressp_3, stressp_4, stressm_1, & ! LCOV_EXCL_LINE
1298 : stressm_2, stressm_3, stressm_4, stress12_1, & ! LCOV_EXCL_LINE
1299 : stress12_2, stress12_3, stress12_4, str1, str2, str3, & ! LCOV_EXCL_LINE
1300 0 : str4, str5, str6, str7, str8, skiptcell)
1301 : !$XXXOMP BARRIER
1302 : call evp1d_stepu(NA_len, rhow, 1, NA_len, cdn_ocn, aiu, &
1303 : uocn, vocn, forcex, forcey, umassdti, fm, uarear, Tbu, & ! LCOV_EXCL_LINE
1304 : uvel_init, vvel_init, uvel, vvel, str1, str2, str3, & ! LCOV_EXCL_LINE
1305 0 : str4, str5, str6, str7, str8, nw, sw, sse, skipucell)
1306 : !$XXXOMP BARRIER
1307 : call evp1d_halo_update(NAVEL_len, 1, NA_len, uvel, vvel, &
1308 0 : halo_parent)
1309 : !$XXXOMP BARRIER
1310 : end do
1311 :
1312 : call evp1d_stress(NA_len, ee, ne, se, 1, NA_len, uvel, vvel, &
1313 : dxT, dyT, hte, htn, htem1, htnm1, strength, stressp_1, & ! LCOV_EXCL_LINE
1314 : stressp_2, stressp_3, stressp_4, stressm_1, stressm_2, & ! LCOV_EXCL_LINE
1315 : stressm_3, stressm_4, stress12_1, stress12_2, stress12_3, & ! LCOV_EXCL_LINE
1316 : stress12_4, str1, str2, str3, str4, str5, str6, str7, & ! LCOV_EXCL_LINE
1317 0 : str8, skiptcell, tarear, divu, rdg_conv, rdg_shear, shear)
1318 : !$XXXOMP BARRIER
1319 : call evp1d_stepu(NA_len, rhow, 1, NA_len, cdn_ocn, aiu, uocn, &
1320 : vocn, forcex, forcey, umassdti, fm, uarear, Tbu, & ! LCOV_EXCL_LINE
1321 : uvel_init, vvel_init, uvel, vvel, str1, str2, str3, str4, & ! LCOV_EXCL_LINE
1322 : str5, str6, str7, str8, nw, sw, sse, skipucell, strintx, & ! LCOV_EXCL_LINE
1323 0 : strinty, taubx, tauby)
1324 : !$XXXOMP BARRIER
1325 : call evp1d_halo_update(NAVEL_len, 1, NA_len, uvel, vvel, &
1326 0 : halo_parent)
1327 : !$XXXOMP END PARALLEL
1328 :
1329 : end if ! master task
1330 :
1331 0 : end subroutine ice_dyn_evp_1d_kernel
1332 :
1333 : !=======================================================================
1334 :
1335 0 : subroutine calc_na(nx, ny, na, iceTmask, iceUmask)
1336 : ! Calculate number of active points
1337 :
1338 : use ice_blocks, only : nghost
1339 :
1340 : implicit none
1341 :
1342 : integer(kind=int_kind), intent(in) :: nx, ny
1343 : logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
1344 : iceTmask, iceUmask
1345 : integer(kind=int_kind), intent(out) :: na
1346 :
1347 : ! local variables
1348 :
1349 : integer(kind=int_kind) :: i, j
1350 :
1351 : character(len=*), parameter :: subname = '(calc_na)'
1352 :
1353 0 : na = 0
1354 : ! NOTE: T mask includes northern and eastern ghost cells
1355 0 : do j = 1 + nghost, ny
1356 0 : do i = 1 + nghost, nx
1357 0 : if (iceTmask(i,j) .or. iceUmask(i,j)) na = na + 1
1358 : end do
1359 : end do
1360 :
1361 0 : end subroutine calc_na
1362 :
1363 : !=======================================================================
1364 :
1365 0 : subroutine calc_2d_indices(nx, ny, na, iceTmask, iceUmask)
1366 :
1367 : use ice_blocks, only : nghost
1368 :
1369 : implicit none
1370 :
1371 : integer(kind=int_kind), intent(in) :: nx, ny, na
1372 : logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
1373 : iceTmask, iceUmask
1374 :
1375 : ! local variables
1376 :
1377 : integer(kind=int_kind) :: i, j, Nmaskt
1378 :
1379 : character(len=*), parameter :: subname = '(calc_2d_indices)'
1380 :
1381 0 : skipucell(:) = .false.
1382 0 : skiptcell(:) = .false.
1383 0 : indi = 0
1384 0 : indj = 0
1385 0 : Nmaskt = 0
1386 : ! NOTE: T mask includes northern and eastern ghost cells
1387 0 : do j = 1 + nghost, ny
1388 0 : do i = 1 + nghost, nx
1389 0 : if (iceTmask(i,j) .or. iceUmask(i,j)) then
1390 0 : Nmaskt = Nmaskt + 1
1391 0 : indi(Nmaskt) = i
1392 0 : indj(Nmaskt) = j
1393 0 : if (.not. iceTmask(i,j)) skiptcell(Nmaskt) = .true.
1394 0 : if (.not. iceUmask(i,j)) skipucell(Nmaskt) = .true.
1395 : ! NOTE: U mask does not include northern and eastern
1396 : ! ghost cells. Skip northern and eastern ghost cells
1397 0 : if (i == nx) skipucell(Nmaskt) = .true.
1398 0 : if (j == ny) skipucell(Nmaskt) = .true.
1399 : end if
1400 : end do
1401 : end do
1402 :
1403 0 : end subroutine calc_2d_indices
1404 :
1405 : !=======================================================================
1406 :
1407 0 : subroutine calc_navel(nx_block, ny_block, na, navel)
1408 : ! Calculate number of active points, including halo points
1409 :
1410 : implicit none
1411 :
1412 : integer(kind=int_kind), intent(in) :: nx_block, ny_block, na
1413 : integer(kind=int_kind), intent(out) :: navel
1414 :
1415 : ! local variables
1416 :
1417 : integer(kind=int_kind) :: iw, i, j
1418 0 : integer(kind=int_kind), dimension(1:na) :: Iin, Iee, Ine, Ise, &
1419 0 : Inw, Isw, Isse
1420 0 : integer(kind=int_kind), dimension(1:7 * na) :: util1, util2
1421 :
1422 : character(len=*), parameter :: subname = '(calc_navel)'
1423 :
1424 : ! calculate additional 1D indices used for finite differences
1425 0 : do iw = 1, na
1426 : ! get 2D indices
1427 0 : i = indi(iw)
1428 0 : j = indj(iw)
1429 : ! calculate 1D indices
1430 0 : Iin(iw) = i + (j - 1) * nx_block ! ( 0, 0) target point
1431 0 : Iee(iw) = i - 1 + (j - 1) * nx_block ! (-1, 0)
1432 0 : Ine(iw) = i - 1 + (j - 2) * nx_block ! (-1, -1)
1433 0 : Ise(iw) = i + (j - 2) * nx_block ! ( 0, -1)
1434 0 : Inw(iw) = i + 1 + (j - 1) * nx_block ! (+1, 0)
1435 0 : Isw(iw) = i + 1 + (j - 0) * nx_block ! (+1, +1)
1436 0 : Isse(iw) = i + (j - 0) * nx_block ! ( 0, +1)
1437 : end do
1438 :
1439 : ! find number of points needed for finite difference calculations
1440 0 : call union(Iin, Iee, na, na, util1, i )
1441 0 : call union(util1, Ine, i, na, util2, j )
1442 0 : call union(util2, Ise, j, na, util1, i )
1443 0 : call union(util1, Inw, i, na, util2, j )
1444 0 : call union(util2, Isw, j, na, util1, i )
1445 0 : call union(util1, Isse, i, na, util2, navel)
1446 :
1447 0 : end subroutine calc_navel
1448 :
1449 : !=======================================================================
1450 :
1451 0 : subroutine convert_2d_1d(nx, ny, na, navel, I_HTE, I_HTN, &
1452 : I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, I_Tbu, & ! LCOV_EXCL_LINE
1453 : I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, I_strinty, & ! LCOV_EXCL_LINE
1454 : I_uvel_init, I_vvel_init, I_strength, I_uvel, I_vvel, I_dxT, & ! LCOV_EXCL_LINE
1455 : I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, I_stressp_4, & ! LCOV_EXCL_LINE
1456 : I_stressm_1, I_stressm_2, I_stressm_3, I_stressm_4, & ! LCOV_EXCL_LINE
1457 0 : I_stress12_1, I_stress12_2, I_stress12_3, I_stress12_4)
1458 :
1459 : implicit none
1460 :
1461 : integer(kind=int_kind), intent(in) :: nx, ny, na, navel
1462 : real (kind=dbl_kind), dimension(nx, ny), intent(in) :: I_HTE, &
1463 : I_HTN, I_cdn_ocn, I_aiu, I_uocn, I_vocn, I_forcex, I_forcey, & ! LCOV_EXCL_LINE
1464 : I_Tbu, I_umassdti, I_fm, I_uarear, I_tarear, I_strintx, & ! LCOV_EXCL_LINE
1465 : I_strinty, I_uvel_init, I_vvel_init, I_strength, I_uvel, & ! LCOV_EXCL_LINE
1466 : I_vvel, I_dxT, I_dyT, I_stressp_1, I_stressp_2, I_stressp_3, & ! LCOV_EXCL_LINE
1467 : I_stressp_4, I_stressm_1, I_stressm_2, I_stressm_3, & ! LCOV_EXCL_LINE
1468 : I_stressm_4, I_stress12_1, I_stress12_2, I_stress12_3, & ! LCOV_EXCL_LINE
1469 : I_stress12_4
1470 :
1471 : ! local variables
1472 :
1473 : integer(kind=int_kind) :: iw, lo, up, j, i, nachk
1474 0 : integer(kind=int_kind), dimension(1:na) :: Iin, Iee, Ine, Ise, &
1475 0 : Inw, Isw, Isse
1476 0 : integer(kind=int_kind), dimension(1:7 * na) :: util1, util2
1477 :
1478 : character(len=*), parameter :: subname = '(convert_2d_1d)'
1479 :
1480 : ! calculate additional 1D indices used for finite differences
1481 0 : do iw = 1, na
1482 : ! get 2D indices
1483 0 : i = indi(iw)
1484 0 : j = indj(iw)
1485 : ! calculate 1D indices
1486 0 : Iin(iw) = i + (j - 1) * nx ! ( 0, 0) target point
1487 0 : Iee(iw) = i - 1 + (j - 1) * nx ! (-1, 0)
1488 0 : Ine(iw) = i - 1 + (j - 2) * nx ! (-1,-1)
1489 0 : Ise(iw) = i + (j - 2) * nx ! ( 0,-1)
1490 0 : Inw(iw) = i + 1 + (j - 1) * nx ! (+1, 0)
1491 0 : Isw(iw) = i + 1 + (j - 0) * nx ! (+1,+1)
1492 0 : Isse(iw) = i + (j - 0) * nx ! ( 0,+1)
1493 : end do
1494 :
1495 : ! find number of points needed for finite difference calculations
1496 0 : call union(Iin, Iee, na, na, util1, i )
1497 0 : call union(util1, Ine, i, na, util2, j )
1498 0 : call union(util2, Ise, j, na, util1, i )
1499 0 : call union(util1, Inw, i, na, util2, j )
1500 0 : call union(util2, Isw, j, na, util1, i )
1501 0 : call union(util1, Isse, i, na, util2, nachk)
1502 :
1503 : ! index vector with sorted target points
1504 0 : do iw = 1, na
1505 0 : indij(iw) = Iin(iw)
1506 : end do
1507 :
1508 : ! sorted additional points
1509 0 : call setdiff(util2, Iin, navel, na, util1, j)
1510 0 : do iw = na + 1, navel
1511 0 : indij(iw) = util1(iw - na)
1512 : end do
1513 :
1514 : ! indices for additional points needed for uvel and vvel
1515 0 : call findXinY(Iee, indij, na, navel, ee)
1516 0 : call findXinY(Ine, indij, na, navel, ne)
1517 0 : call findXinY(Ise, indij, na, navel, se)
1518 0 : call findXinY(Inw, indij, na, navel, nw)
1519 0 : call findXinY(Isw, indij, na, navel, sw)
1520 0 : call findXinY(Isse, indij, na, navel, sse)
1521 :
1522 0 : !$OMP PARALLEL PRIVATE(iw, lo, up, j, i)
1523 : ! write 1D arrays from 2D arrays (target points)
1524 0 : call domp_get_domain(1, na, lo, up)
1525 0 : do iw = lo, up
1526 : ! get 2D indices
1527 0 : i = indi(iw)
1528 0 : j = indj(iw)
1529 : ! map
1530 0 : uvel(iw) = I_uvel(i, j)
1531 0 : vvel(iw) = I_vvel(i, j)
1532 0 : cdn_ocn(iw) = I_cdn_ocn(i, j)
1533 0 : aiu(iw) = I_aiu(i, j)
1534 0 : uocn(iw) = I_uocn(i, j)
1535 0 : vocn(iw) = I_vocn(i, j)
1536 0 : forcex(iw) = I_forcex(i, j)
1537 0 : forcey(iw) = I_forcey(i, j)
1538 0 : Tbu(iw) = I_Tbu(i, j)
1539 0 : umassdti(iw) = I_umassdti(i, j)
1540 0 : fm(iw) = I_fm(i, j)
1541 0 : tarear(iw) = I_tarear(i, j)
1542 0 : uarear(iw) = I_uarear(i, j)
1543 0 : strintx(iw) = I_strintx(i, j)
1544 0 : strinty(iw) = I_strinty(i, j)
1545 0 : uvel_init(iw) = I_uvel_init(i, j)
1546 0 : vvel_init(iw) = I_vvel_init(i, j)
1547 0 : strength(iw) = I_strength(i, j)
1548 0 : dxT(iw) = I_dxT(i, j)
1549 0 : dyT(iw) = I_dyT(i, j)
1550 0 : stressp_1(iw) = I_stressp_1(i, j)
1551 0 : stressp_2(iw) = I_stressp_2(i, j)
1552 0 : stressp_3(iw) = I_stressp_3(i, j)
1553 0 : stressp_4(iw) = I_stressp_4(i, j)
1554 0 : stressm_1(iw) = I_stressm_1(i, j)
1555 0 : stressm_2(iw) = I_stressm_2(i, j)
1556 0 : stressm_3(iw) = I_stressm_3(i, j)
1557 0 : stressm_4(iw) = I_stressm_4(i, j)
1558 0 : stress12_1(iw) = I_stress12_1(i, j)
1559 0 : stress12_2(iw) = I_stress12_2(i, j)
1560 0 : stress12_3(iw) = I_stress12_3(i, j)
1561 0 : stress12_4(iw) = I_stress12_4(i, j)
1562 0 : HTE(iw) = I_HTE(i, j)
1563 0 : HTN(iw) = I_HTN(i, j)
1564 0 : HTEm1(iw) = I_HTE(i - 1, j)
1565 0 : HTNm1(iw) = I_HTN(i, j - 1)
1566 : end do
1567 : ! write 1D arrays from 2D arrays (additional points)
1568 0 : call domp_get_domain(na + 1, navel, lo, up)
1569 0 : do iw = lo, up
1570 : ! get 2D indices
1571 0 : j = int((indij(iw) - 1) / (nx)) + 1
1572 0 : i = indij(iw) - (j - 1) * nx
1573 : ! map
1574 0 : uvel(iw) = I_uvel(i, j)
1575 0 : vvel(iw) = I_vvel(i, j)
1576 : end do
1577 : !$OMP END PARALLEL
1578 :
1579 0 : end subroutine convert_2d_1d
1580 :
1581 : !=======================================================================
1582 :
1583 0 : subroutine calc_halo_parent(nx, ny, na, navel, I_iceTmask)
1584 :
1585 : implicit none
1586 :
1587 : integer(kind=int_kind), intent(in) :: nx, ny, na, navel
1588 : logical(kind=log_kind), dimension(nx, ny), intent(in) :: &
1589 : I_iceTmask
1590 :
1591 : ! local variables
1592 :
1593 : integer(kind=int_kind) :: iw, i, j
1594 0 : integer(kind=int_kind), dimension(1:navel) :: Ihalo
1595 :
1596 : character(len=*), parameter :: subname = '(calc_halo_parent)'
1597 :
1598 : !-----------------------------------------------------------------
1599 : ! Indices for halo update:
1600 : ! 0: no halo point
1601 : ! >0: index for halo point parent, related to indij vector
1602 : !
1603 : ! TODO: Implement for nghost > 1
1604 : ! TODO: Implement for tripole grids
1605 : !-----------------------------------------------------------------
1606 :
1607 0 : Ihalo(:) = 0
1608 0 : halo_parent(:) = 0
1609 :
1610 0 : do iw = 1, navel
1611 0 : j = int((indij(iw) - 1) / (nx)) + 1
1612 0 : i = indij(iw) - (j - 1) * nx
1613 : ! if within ghost zone
1614 0 : if (i == nx .and. I_iceTmask(2, j) ) Ihalo(iw) = 2 + (j - 1) * nx
1615 0 : if (i == 1 .and. I_iceTmask(nx - 1, j) ) Ihalo(iw) = (nx - 1) + (j - 1) * nx
1616 0 : if (j == ny .and. I_iceTmask(i, 2) ) Ihalo(iw) = i + nx
1617 0 : if (j == 1 .and. I_iceTmask(i, ny - 1) ) Ihalo(iw) = i + (ny - 2) * nx
1618 : end do
1619 :
1620 : ! relate halo indices to indij vector
1621 0 : call findXinY_halo(Ihalo, indij, navel, navel, halo_parent)
1622 :
1623 0 : end subroutine calc_halo_parent
1624 :
1625 : !=======================================================================
1626 :
1627 0 : subroutine union(x, y, nx, ny, xy, nxy)
1628 : ! Find union (xy) of two sorted integer vectors (x and y), i.e.
1629 : ! combined values of the two vectors with no repetitions
1630 :
1631 : implicit none
1632 :
1633 : integer(int_kind), intent(in) :: nx, ny
1634 : integer(int_kind), intent(in) :: x(1:nx), y(1:ny)
1635 : integer(int_kind), intent(out) :: xy(1:nx + ny)
1636 : integer(int_kind), intent(out) :: nxy
1637 :
1638 : ! local variables
1639 :
1640 : integer(int_kind) :: i, j, k
1641 :
1642 : character(len=*), parameter :: subname = '(union)'
1643 :
1644 0 : i = 1
1645 0 : j = 1
1646 0 : k = 1
1647 0 : do while (i <= nx .and. j <= ny)
1648 0 : if (x(i) < y(j)) then
1649 0 : xy(k) = x(i)
1650 0 : i = i + 1
1651 0 : else if (x(i) > y(j)) then
1652 0 : xy(k) = y(j)
1653 0 : j = j + 1
1654 : else
1655 0 : xy(k) = x(i)
1656 0 : i = i + 1
1657 0 : j = j + 1
1658 : end if
1659 0 : k = k + 1
1660 : end do
1661 :
1662 : ! the rest
1663 0 : do while (i <= nx)
1664 0 : xy(k) = x(i)
1665 0 : i = i + 1
1666 0 : k = k + 1
1667 : end do
1668 0 : do while (j <= ny)
1669 0 : xy(k) = y(j)
1670 0 : j = j + 1
1671 0 : k = k + 1
1672 : end do
1673 0 : nxy = k - 1
1674 :
1675 0 : end subroutine union
1676 :
1677 : !=======================================================================
1678 :
1679 0 : subroutine setdiff(x, y, nx, ny, xy, nxy)
1680 : ! Find element (xy) of two sorted integer vectors (x and y) that
1681 : ! are in x, but not in y, or in y, but not in x
1682 :
1683 : implicit none
1684 :
1685 : integer(kind=int_kind), intent(in) :: nx, ny
1686 : integer(kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
1687 : integer(kind=int_kind), intent(out) :: xy(1:nx + ny)
1688 : integer(kind=int_kind), intent(out) :: nxy
1689 :
1690 : ! local variables
1691 :
1692 : integer(kind=int_kind) :: i, j, k
1693 :
1694 : character(len=*), parameter :: subname = '(setdiff)'
1695 :
1696 0 : i = 1
1697 0 : j = 1
1698 0 : k = 1
1699 0 : do while (i <= nx .and. j <= ny)
1700 0 : if (x(i) < y(j)) then
1701 0 : xy(k) = x(i)
1702 0 : i = i + 1
1703 0 : k = k + 1
1704 0 : else if (x(i) > y(j)) then
1705 0 : xy(k) = y(j)
1706 0 : j = j + 1
1707 0 : k = k + 1
1708 : else
1709 0 : i = i + 1
1710 0 : j = j + 1
1711 : end if
1712 : end do
1713 :
1714 : ! the rest
1715 0 : do while (i <= nx)
1716 0 : xy(k) = x(i)
1717 0 : i = i + 1
1718 0 : k = k + 1
1719 : end do
1720 0 : do while (j <= ny)
1721 0 : xy(k) = y(j)
1722 0 : j = j + 1
1723 0 : k = k + 1
1724 : end do
1725 0 : nxy = k - 1
1726 :
1727 0 : end subroutine setdiff
1728 :
1729 : !========================================================================
1730 :
1731 0 : subroutine findXinY(x, y, nx, ny, indx)
1732 : ! Find indx vector so that x(1:na) = y(indx(1:na))
1733 : !
1734 : ! Conditions:
1735 : ! * EVERY item in x is found in y
1736 : ! * x(1:nx) is a sorted integer vector
1737 : ! * y(1:ny) consists of two sorted integer vectors:
1738 : ! [y(1:nx); y(nx + 1:ny)]
1739 : ! * ny >= nx
1740 : !
1741 : ! Return: indx(1:na)
1742 :
1743 : implicit none
1744 :
1745 : integer (kind=int_kind), intent(in) :: nx, ny
1746 : integer (kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
1747 : integer (kind=int_kind), intent(out) :: indx(1:nx)
1748 :
1749 : ! local variables
1750 :
1751 : integer (kind=int_kind) :: i, j1, j2
1752 :
1753 : character(len=*), parameter :: subname = '(findXinY)'
1754 :
1755 0 : i = 1
1756 0 : j1 = 1
1757 0 : j2 = nx + 1
1758 0 : do while (i <= nx)
1759 0 : if (x(i) == y(j1)) then
1760 0 : indx(i) = j1
1761 0 : i = i + 1
1762 0 : j1 = j1 + 1
1763 0 : else if (x(i) == y(j2)) then
1764 0 : indx(i) = j2
1765 0 : i = i + 1
1766 0 : j2 = j2 + 1
1767 0 : else if (x(i) > y(j1)) then
1768 0 : j1 = j1 + 1
1769 0 : else if (x(i) > y(j2)) then
1770 0 : j2 = j2 + 1
1771 : else
1772 : call abort_ice(subname &
1773 0 : // ': ERROR: conditions not met')
1774 : end if
1775 : end do
1776 :
1777 0 : end subroutine findXinY
1778 :
1779 : !=======================================================================
1780 :
1781 0 : subroutine findXinY_halo(x, y, nx, ny, indx)
1782 : ! Find indx vector so that x(1:na) = y(indx(1:na))
1783 : !
1784 : ! Conditions:
1785 : ! * EVERY item in x is found in y,
1786 : ! except for x == 0, where indx = 0 is returned
1787 : ! * x(1:nx) is a non-sorted integer vector
1788 : ! * y(1:ny) is a sorted integer vector
1789 : ! * ny >= nx
1790 : !
1791 : ! Return: indx(1:na)
1792 :
1793 : implicit none
1794 :
1795 : integer (kind=int_kind), intent(in) :: nx, ny
1796 : integer (kind=int_kind), intent(in) :: x(1:nx), y(1:ny)
1797 : integer (kind=int_kind), intent(out) :: indx(1:nx)
1798 :
1799 : ! local variables
1800 :
1801 : integer (kind=int_kind) :: i, j1, nloop
1802 :
1803 : character(len=*), parameter :: subname = '(findXinY_halo)'
1804 :
1805 0 : nloop = 1
1806 0 : i = 1
1807 0 : j1 = int((ny + 1) / 2) ! initial guess in the middle
1808 0 : do while (i <= nx)
1809 0 : if (x(i) == 0) then
1810 0 : indx(i) = 0
1811 0 : i = i + 1
1812 0 : nloop = 1
1813 0 : else if (x(i) == y(j1)) then
1814 0 : indx(i) = j1
1815 0 : i = i + 1
1816 0 : j1 = j1 + 1
1817 : ! initial guess in the middle
1818 0 : if (j1 > ny) j1 = int((ny + 1) / 2)
1819 0 : nloop = 1
1820 0 : else if (x(i) < y(j1)) then
1821 0 : j1 = 1
1822 0 : else if (x(i) > y(j1)) then
1823 0 : j1 = j1 + 1
1824 0 : if (j1 > ny) then
1825 0 : j1 = 1
1826 0 : nloop = nloop + 1
1827 0 : if (nloop > 2) then
1828 : ! stop for infinite loop. This check should not be
1829 : ! necessary for halo
1830 0 : call abort_ice(subname // ' ERROR: too many loops')
1831 : end if
1832 : end if
1833 : end if
1834 : end do
1835 :
1836 0 : end subroutine findXinY_halo
1837 :
1838 : !=======================================================================
1839 :
1840 0 : subroutine numainit(l, u, uu)
1841 :
1842 : use ice_constants, only : c0
1843 :
1844 : implicit none
1845 :
1846 : integer(kind=int_kind), intent(in) :: l, u, uu
1847 :
1848 : ! local variables
1849 :
1850 : integer(kind=int_kind) :: lo, up
1851 :
1852 : character(len=*), parameter :: subname = '(numainit)'
1853 :
1854 0 : call domp_get_domain(l, u, lo, up)
1855 0 : ee(lo:up) = 0
1856 0 : ne(lo:up) = 0
1857 0 : se(lo:up) = 0
1858 0 : sse(lo:up) = 0
1859 0 : nw(lo:up) = 0
1860 0 : sw(lo:up) = 0
1861 0 : halo_parent(lo:up) = 0
1862 0 : strength(lo:up) = c0
1863 0 : uvel(lo:up) = c0
1864 0 : vvel(lo:up) = c0
1865 0 : uvel_init(lo:up) = c0
1866 0 : vvel_init(lo:up) = c0
1867 0 : uocn(lo:up) = c0
1868 0 : vocn(lo:up) = c0
1869 0 : dxT(lo:up) = c0
1870 0 : dyT(lo:up) = c0
1871 0 : HTE(lo:up) = c0
1872 0 : HTN(lo:up) = c0
1873 0 : HTEm1(lo:up) = c0
1874 0 : HTNm1(lo:up) = c0
1875 0 : stressp_1(lo:up) = c0
1876 0 : stressp_2(lo:up) = c0
1877 0 : stressp_3(lo:up) = c0
1878 0 : stressp_4(lo:up) = c0
1879 0 : stressm_1(lo:up) = c0
1880 0 : stressm_2(lo:up) = c0
1881 0 : stressm_3(lo:up) = c0
1882 0 : stressm_4(lo:up) = c0
1883 0 : stress12_1(lo:up) = c0
1884 0 : stress12_2(lo:up) = c0
1885 0 : stress12_3(lo:up) = c0
1886 0 : stress12_4(lo:up) = c0
1887 0 : tarear(lo:up) = c0
1888 0 : Tbu(lo:up) = c0
1889 0 : taubx(lo:up) = c0
1890 0 : tauby(lo:up) = c0
1891 0 : divu(lo:up) = c0
1892 0 : rdg_conv(lo:up) = c0
1893 0 : rdg_shear(lo:up) = c0
1894 0 : shear(lo:up) = c0
1895 0 : str1(lo:up) = c0
1896 0 : str2(lo:up) = c0
1897 0 : str3(lo:up) = c0
1898 0 : str4(lo:up) = c0
1899 0 : str5(lo:up) = c0
1900 0 : str6(lo:up) = c0
1901 0 : str7(lo:up) = c0
1902 0 : str8(lo:up) = c0
1903 :
1904 0 : call domp_get_domain(u + 1, uu, lo, up)
1905 0 : halo_parent(lo:up) = 0
1906 0 : uvel(lo:up) = c0
1907 0 : vvel(lo:up) = c0
1908 0 : str1(lo:up) = c0
1909 0 : str2(lo:up) = c0
1910 0 : str3(lo:up) = c0
1911 0 : str4(lo:up) = c0
1912 0 : str5(lo:up) = c0
1913 0 : str6(lo:up) = c0
1914 0 : str7(lo:up) = c0
1915 0 : str8(lo:up) = c0
1916 :
1917 0 : end subroutine numainit
1918 :
1919 : !=======================================================================
1920 :
1921 : end module ice_dyn_evp_1d
|
