Line data Source code
1 : !=======================================================================
2 :
3 : ! Driver for ice mechanical redistribution (ridging)
4 : !
5 : ! See these references:
6 : !
7 : ! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength
8 : ! in modeling the thickness distribution of Arctic sea ice,
9 : ! J. Geophys. Res., 100, 18,611-18,626.
10 : !
11 : ! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice
12 : ! cover, Mon. Wea. Rev., 108, 1943-1973, 1980.
13 : !
14 : ! Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki, 2007:
15 : ! Improving ridging schemes for high-resolution sea ice models.
16 : ! J. Geophys. Res. 112, C03S91, doi:10.1029/2005JC003355.
17 : !
18 : ! Rothrock, D. A., 1975: The energetics of the plastic deformation of
19 : ! pack ice by ridging, J. Geophys. Res., 80, 4514-4519.
20 : !
21 : ! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,
22 : ! 1975: The thickness distribution of sea ice, J. Geophys. Res.,
23 : ! 80, 4501-4513.
24 : !
25 : ! authors: William H. Lipscomb, LANL
26 : ! Elizabeth C. Hunke, LANL
27 : !
28 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
29 : ! 2004: Block structure added by William Lipscomb
30 : ! 2006: New options for participation and redistribution (WHL)
31 : ! 2006: Streamlined for efficiency by Elizabeth Hunke
32 : ! Converted to free source form (F90)
33 :
34 : module icepack_mechred
35 :
36 : use icepack_kinds
37 : use icepack_parameters, only: c0, c1, c2, c10, c25, Cf, Cp, Pstar, Cstar
38 : use icepack_parameters, only: p05, p15, p25, p333, p5
39 : use icepack_parameters, only: puny, Lfresh, rhoi, rhos
40 : use icepack_parameters, only: icepack_chkoptargflag
41 :
42 : use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
43 : use icepack_parameters, only: conserv_check, z_tracers
44 : use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, n_aero
45 : use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
46 : use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
47 : use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice
48 : use icepack_tracers, only: nt_apnd, nt_hpnd
49 : use icepack_tracers, only: n_iso, bio_index
50 : use icepack_tracers, only: icepack_compute_tracers
51 :
52 : use icepack_warnings, only: warnstr, icepack_warnings_add
53 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
54 :
55 : use icepack_itd, only: column_sum
56 : use icepack_itd, only: column_conservation_check
57 : use icepack_itd, only: cleanup_itd
58 :
59 : implicit none
60 :
61 : private
62 : public :: icepack_ice_strength, &
63 : icepack_step_ridge
64 :
65 : real (kind=dbl_kind), parameter :: &
66 : exp_argmax = 100.0_dbl_kind, & ! maximum argument of exponential for underflow ! LCOV_EXCL_LINE
67 : Cs = p25 , & ! fraction of shear energy contrbtng to ridging ! LCOV_EXCL_LINE
68 : fsnowrdg = p5 , & ! snow fraction that survives in ridging ! LCOV_EXCL_LINE
69 : Gstar = p15 , & ! max value of G(h) that participates ! LCOV_EXCL_LINE
70 : ! (krdg_partic = 0)
71 : astar = p05 , & ! e-folding scale for G(h) participation
72 : ! (krdg_partic = 1)
73 : maxraft= c1 , & ! max value of hrmin - hi = max thickness
74 : ! of ice that rafts (m)
75 : Hstar = c25 ! determines mean thickness of ridged ice (m)
76 : ! (krdg_redist = 0)
77 : ! Flato & Hibler (1995) have Hstar = 100
78 :
79 : !=======================================================================
80 :
81 : contains
82 :
83 : !=======================================================================
84 :
85 : ! Compute changes in the ice thickness distribution due to divergence
86 : ! and shear.
87 : !
88 : ! author: William H. Lipscomb, LANL
89 :
90 1498353240 : subroutine ridge_ice (dt, ndtd, &
91 1498353240 : hin_max, & ! LCOV_EXCL_LINE
92 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
93 1498353240 : aicen, trcrn, & ! LCOV_EXCL_LINE
94 1498353240 : vicen, vsnon, & ! LCOV_EXCL_LINE
95 : aice0, & ! LCOV_EXCL_LINE
96 1498353240 : trcr_depend, trcr_base, & ! LCOV_EXCL_LINE
97 1498353240 : n_trcr_strata, & ! LCOV_EXCL_LINE
98 1498353240 : nt_strata, & ! LCOV_EXCL_LINE
99 : krdg_partic, krdg_redist,& ! LCOV_EXCL_LINE
100 : mu_rdg, tr_brine, & ! LCOV_EXCL_LINE
101 : dardg1dt, dardg2dt, & ! LCOV_EXCL_LINE
102 : dvirdgdt, opening, & ! LCOV_EXCL_LINE
103 1498353240 : fpond, flux_bio, & ! LCOV_EXCL_LINE
104 : fresh, fhocn, & ! LCOV_EXCL_LINE
105 2996706480 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
106 1498353240 : aparticn, krdgn, & ! LCOV_EXCL_LINE
107 1498353240 : aredistn, vredistn, & ! LCOV_EXCL_LINE
108 1498353240 : dardg1ndt, dardg2ndt, & ! LCOV_EXCL_LINE
109 1498353240 : dvirdgndt, Tf, & ! LCOV_EXCL_LINE
110 1498353240 : araftn, vraftn, & ! LCOV_EXCL_LINE
111 : closing )
112 :
113 : integer (kind=int_kind), intent(in) :: &
114 : ndtd ! number of dynamics subcycles
115 :
116 : real (kind=dbl_kind), intent(in) :: &
117 : mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) ! LCOV_EXCL_LINE
118 : dt ! time step
119 :
120 : real (kind=dbl_kind), intent(in) :: &
121 : Tf ! freezing temperature
122 :
123 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
124 : hin_max ! category limits (m)
125 :
126 : real (kind=dbl_kind), intent(in) :: &
127 : rdg_conv , & ! normalized energy dissipation due to convergence (1/s) ! LCOV_EXCL_LINE
128 : rdg_shear ! normalized energy dissipation due to shear (1/s)
129 :
130 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
131 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
132 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
133 : vsnon ! volume per unit area of snow (m)
134 :
135 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
136 : trcrn ! ice tracers
137 :
138 : real (kind=dbl_kind), intent(inout) :: &
139 : aice0 ! concentration of open water
140 :
141 : integer (kind=int_kind), dimension (:), intent(in) :: &
142 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
143 : n_trcr_strata ! number of underlying tracer layers
144 :
145 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
146 : trcr_base ! = 0 or 1 depending on tracer dependency
147 : ! argument 2: (1) aice, (2) vice, (3) vsno
148 :
149 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
150 : nt_strata ! indices of underlying tracer layers
151 :
152 : integer (kind=int_kind), intent(in) :: &
153 : krdg_partic, & ! selects participation function ! LCOV_EXCL_LINE
154 : krdg_redist ! selects redistribution function
155 :
156 : logical (kind=log_kind), intent(in) :: &
157 : tr_brine ! if .true., brine height differs from ice thickness
158 :
159 : ! optional history fields
160 : real (kind=dbl_kind), intent(inout), optional :: &
161 : dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
162 : dardg2dt , & ! rate of fractional area gain by new ridges (1/s) ! LCOV_EXCL_LINE
163 : dvirdgdt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
164 : opening , & ! rate of opening due to divergence/shear (1/s) ! LCOV_EXCL_LINE
165 : closing , & ! rate of closing due to divergence/shear (1/s) ! LCOV_EXCL_LINE
166 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
167 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
168 : fhocn ! net heat flux to ocean (W/m^2)
169 :
170 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
171 : dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
172 : dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) ! LCOV_EXCL_LINE
173 : dvirdgndt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
174 : aparticn , & ! participation function ! LCOV_EXCL_LINE
175 : krdgn , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
176 : araftn , & ! rafting ice area ! LCOV_EXCL_LINE
177 : vraftn , & ! rafting ice volume ! LCOV_EXCL_LINE
178 : aredistn , & ! redistribution function: fraction of new ridge area ! LCOV_EXCL_LINE
179 : vredistn ! redistribution function: fraction of new ridge volume
180 :
181 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
182 : faero_ocn ! aerosol flux to ocean (kg/m^2/s)
183 :
184 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
185 : flux_bio ! biological and zaerosol flux to ocean (kg/m^2/s)
186 :
187 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
188 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
189 :
190 : ! local variables
191 :
192 : real (kind=dbl_kind), dimension (ncat) :: &
193 2996706480 : eicen ! energy of melting for each ice layer (J/m^2)
194 :
195 : real (kind=dbl_kind), dimension (ncat) :: &
196 2996706480 : esnon, & ! energy of melting for each snow layer (J/m^2) ! LCOV_EXCL_LINE
197 2996706480 : vbrin, & ! ice volume with defined by brine height (m) ! LCOV_EXCL_LINE
198 2996706480 : sicen ! Bulk salt in h ice (ppt*m)
199 :
200 : real (kind=dbl_kind) :: &
201 : asum , & ! sum of ice and open water area ! LCOV_EXCL_LINE
202 : aksum , & ! ratio of area removed to area ridged ! LCOV_EXCL_LINE
203 : msnow_mlt , & ! mass of snow added to ocean (kg m-2) ! LCOV_EXCL_LINE
204 : esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) ! LCOV_EXCL_LINE
205 : mpond , & ! mass of pond added to ocean (kg m-2) ! LCOV_EXCL_LINE
206 : closing_net, & ! net rate at which area is removed (1/s) ! LCOV_EXCL_LINE
207 : ! (ridging ice area - area of new ridges) / dt
208 : divu_adv , & ! divu as implied by transport scheme (1/s)
209 : opning , & ! rate of opening due to divergence/shear ! LCOV_EXCL_LINE
210 : ! opning is a local variable;
211 : ! opening is the history diagnostic variable
212 : ardg1 , & ! fractional area loss by ridging ice
213 : ardg2 , & ! fractional area gain by new ridges ! LCOV_EXCL_LINE
214 : virdg , & ! ice volume ridged ! LCOV_EXCL_LINE
215 : aopen ! area opening due to divergence/shear
216 :
217 : real (kind=dbl_kind), dimension (n_aero) :: &
218 2996706480 : maero ! aerosol mass added to ocean (kg m-2)
219 :
220 : real (kind=dbl_kind), dimension (nbtrcr) :: &
221 2996706480 : mbio ! bio mass added to ocean (mmol or kg m-2)
222 :
223 : real (kind=dbl_kind), dimension (n_iso) :: &
224 2996706480 : miso ! isotope mass added to ocean (kg m-2)
225 :
226 : real (kind=dbl_kind), dimension (0:ncat) :: &
227 1498353240 : apartic ! participation function; fraction of ridging
228 : ! and closing associated w/ category n
229 :
230 : real (kind=dbl_kind), dimension (ncat) :: &
231 2996706480 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
232 2996706480 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
233 2996706480 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
234 2996706480 : krdg , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
235 2996706480 : ardg1n , & ! area of ice ridged ! LCOV_EXCL_LINE
236 2996706480 : ardg2n , & ! area of new ridges ! LCOV_EXCL_LINE
237 2996706480 : virdgn , & ! ridging ice volume ! LCOV_EXCL_LINE
238 1498353240 : mraftn ! rafting ice mask
239 :
240 : real (kind=dbl_kind) :: &
241 : vice_init, vice_final, & ! ice volume summed over categories ! LCOV_EXCL_LINE
242 : vsno_init, vsno_final, & ! snow volume summed over categories ! LCOV_EXCL_LINE
243 : eice_init, eice_final, & ! ice energy summed over layers ! LCOV_EXCL_LINE
244 : vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories ! LCOV_EXCL_LINE
245 : sice_init ,sice_final, & ! ice bulk salinity summed over categories ! LCOV_EXCL_LINE
246 : esno_init, esno_final ! snow energy summed over layers
247 :
248 : integer (kind=int_kind), parameter :: &
249 : nitermax = 20 ! max number of ridging iterations
250 :
251 : integer (kind=int_kind) :: &
252 : n , & ! thickness category index ! LCOV_EXCL_LINE
253 : niter , & ! iteration counter ! LCOV_EXCL_LINE
254 : k , & ! vertical index ! LCOV_EXCL_LINE
255 : it ! tracer index
256 :
257 : real (kind=dbl_kind) :: &
258 : dti ! 1 / dt
259 :
260 : logical (kind=log_kind) :: &
261 : iterate_ridging ! if true, repeat the ridging
262 :
263 : character (len=char_len) :: &
264 : fieldid ! field identifier
265 :
266 : character(len=*),parameter :: subname='(ridge_ice)'
267 :
268 : !-----------------------------------------------------------------
269 : ! Initialize
270 : !-----------------------------------------------------------------
271 :
272 1498353240 : msnow_mlt = c0
273 1498353240 : esnow_mlt = c0
274 2943906768 : maero (:) = c0
275 2453087592 : mbio (:) = c0
276 1542162072 : miso (:) = c0
277 1498353240 : mpond = c0
278 1498353240 : ardg1 = c0
279 1498353240 : ardg2 = c0
280 1498353240 : virdg = c0
281 8967638592 : ardg1n(:) = c0
282 8967638592 : ardg2n(:) = c0
283 8967638592 : virdgn(:) = c0
284 8967638592 : mraftn(:) = c0
285 1498353240 : aopen = c0
286 :
287 : !-----------------------------------------------------------------
288 : ! Compute area of ice plus open water before ridging.
289 : !-----------------------------------------------------------------
290 :
291 1498353240 : call asum_ridging (aicen, aice0, asum)
292 1498353240 : if (icepack_warnings_aborted(subname)) return
293 :
294 : !-----------------------------------------------------------------
295 : ! Compute the area opening and closing.
296 : !-----------------------------------------------------------------
297 :
298 1498353240 : if (present(opening) .and. present(closing)) then
299 :
300 0 : opning = opening
301 0 : closing_net = closing
302 0 : divu_adv = opening - closing
303 :
304 : else
305 :
306 : call ridge_prep (dt, hin_max, &
307 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
308 : asum, closing_net, & ! LCOV_EXCL_LINE
309 1498353240 : divu_adv, opning)
310 :
311 : endif
312 :
313 1498353240 : if (icepack_warnings_aborted(subname)) return
314 :
315 : !-----------------------------------------------------------------
316 : ! Compute initial values of conserved quantities.
317 : !-----------------------------------------------------------------
318 :
319 1498353240 : if (conserv_check) then
320 :
321 122395728 : do n = 1, ncat
322 104910624 : eicen(n) = c0
323 104910624 : esnon(n) = c0
324 104910624 : sicen(n) = c0
325 104910624 : vbrin(n) = c0
326 :
327 839284992 : do k = 1, nilyr
328 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
329 734374368 : * vicen(n)/real(nilyr,kind=dbl_kind)
330 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
331 839284992 : * vicen(n)/real(nilyr,kind=dbl_kind)
332 : enddo
333 209821248 : do k = 1, nslyr
334 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
335 209821248 : * vsnon(n)/real(nslyr,kind=dbl_kind)
336 : enddo
337 104910624 : vbrin(n) = vicen(n)
338 122395728 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vicen(n)
339 : enddo ! n
340 :
341 : call column_sum (ncat, &
342 17485104 : vicen, vice_init)
343 17485104 : if (icepack_warnings_aborted(subname)) return
344 : call column_sum (ncat, &
345 17485104 : vsnon, vsno_init)
346 17485104 : if (icepack_warnings_aborted(subname)) return
347 : call column_sum (ncat, &
348 17485104 : eicen, eice_init)
349 17485104 : if (icepack_warnings_aborted(subname)) return
350 : call column_sum (ncat, &
351 17485104 : esnon, esno_init)
352 17485104 : if (icepack_warnings_aborted(subname)) return
353 : call column_sum (ncat, &
354 17485104 : sicen, sice_init)
355 17485104 : if (icepack_warnings_aborted(subname)) return
356 : call column_sum (ncat, &
357 17485104 : vbrin, vbri_init)
358 17485104 : if (icepack_warnings_aborted(subname)) return
359 :
360 : endif
361 :
362 1498353257 : rdg_iteration: do niter = 1, nitermax
363 :
364 : !-----------------------------------------------------------------
365 : ! Compute the thickness distribution of ridging ice
366 : ! and various quantities associated with the new ridged ice.
367 : !-----------------------------------------------------------------
368 :
369 : call ridge_itd (aice0, &
370 : aicen, vicen, & ! LCOV_EXCL_LINE
371 : krdg_partic, krdg_redist,& ! LCOV_EXCL_LINE
372 : mu_rdg, & ! LCOV_EXCL_LINE
373 : aksum, apartic, & ! LCOV_EXCL_LINE
374 : hrmin, hrmax, & ! LCOV_EXCL_LINE
375 : hrexp, krdg, & ! LCOV_EXCL_LINE
376 : aparticn, krdgn, & ! LCOV_EXCL_LINE
377 1498353257 : mraftn)
378 1498353257 : if (icepack_warnings_aborted(subname)) return
379 :
380 : !-----------------------------------------------------------------
381 : ! Redistribute area, volume, and energy.
382 : !-----------------------------------------------------------------
383 :
384 : call ridge_shift (dt, hin_max, &
385 : aicen, trcrn, & ! LCOV_EXCL_LINE
386 : vicen, vsnon, & ! LCOV_EXCL_LINE
387 : aice0, trcr_depend, & ! LCOV_EXCL_LINE
388 : trcr_base, n_trcr_strata,& ! LCOV_EXCL_LINE
389 : nt_strata, krdg_redist, & ! LCOV_EXCL_LINE
390 : aksum, apartic, & ! LCOV_EXCL_LINE
391 : hrmin, hrmax, & ! LCOV_EXCL_LINE
392 : hrexp, krdg, & ! LCOV_EXCL_LINE
393 : closing_net, opning, & ! LCOV_EXCL_LINE
394 : ardg1, ardg2, & ! LCOV_EXCL_LINE
395 : virdg, aopen, & ! LCOV_EXCL_LINE
396 : ardg1n, ardg2n, & ! LCOV_EXCL_LINE
397 : virdgn, & ! LCOV_EXCL_LINE
398 : msnow_mlt, esnow_mlt, & ! LCOV_EXCL_LINE
399 : maero, miso, & ! LCOV_EXCL_LINE
400 : mpond, Tf, & ! LCOV_EXCL_LINE
401 : aredistn, vredistn, & ! LCOV_EXCL_LINE
402 1498353257 : mbio)
403 1498353257 : if (icepack_warnings_aborted(subname)) return
404 :
405 : !-----------------------------------------------------------------
406 : ! Make sure the new area = 1. If not (because the closing
407 : ! and opening rates were reduced above), prepare to ridge again
408 : ! with new rates.
409 : !-----------------------------------------------------------------
410 :
411 1498353257 : call asum_ridging (aicen, aice0, asum)
412 1498353257 : if (icepack_warnings_aborted(subname)) return
413 :
414 1498353257 : if (abs(asum - c1) < puny) then
415 1498353240 : iterate_ridging = .false.
416 1498353240 : closing_net = c0
417 1498353240 : opning = c0
418 : else
419 17 : iterate_ridging = .true.
420 17 : divu_adv = (c1 - asum) / dt
421 17 : closing_net = max(c0, -divu_adv)
422 17 : opning = max(c0, divu_adv)
423 : endif
424 :
425 : !-----------------------------------------------------------------
426 : ! If done, exit. If not, prepare to ridge again.
427 : !-----------------------------------------------------------------
428 :
429 1498353257 : if (.not.iterate_ridging) then
430 1498353240 : exit rdg_iteration
431 : endif
432 :
433 17 : if (niter == nitermax) then
434 0 : write(warnstr,*) ' '
435 0 : call icepack_warnings_add(warnstr)
436 0 : write(warnstr,*) subname, 'Exceeded max number of ridging iterations'
437 0 : call icepack_warnings_add(warnstr)
438 0 : write(warnstr,*) subname, 'max =',nitermax
439 0 : call icepack_warnings_add(warnstr)
440 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
441 0 : call icepack_warnings_add(subname//" ridge_ice: Exceeded max number of ridging iterations" )
442 0 : return
443 : endif
444 :
445 : enddo rdg_iteration ! niter
446 :
447 : !-----------------------------------------------------------------
448 : ! Compute final values of conserved quantities.
449 : ! Check for conservation (allowing for snow thrown into ocean).
450 : !-----------------------------------------------------------------
451 :
452 1498353240 : if (conserv_check) then
453 :
454 122395728 : do n = 1, ncat
455 104910624 : eicen(n) = c0
456 104910624 : esnon(n) = c0
457 104910624 : sicen(n) = c0
458 104910624 : vbrin(n) = c0
459 :
460 839284992 : do k = 1, nilyr
461 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
462 734374368 : * vicen(n)/real(nilyr,kind=dbl_kind)
463 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
464 839284992 : * vicen(n)/real(nilyr,kind=dbl_kind)
465 : enddo
466 209821248 : do k = 1, nslyr
467 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
468 209821248 : * vsnon(n)/real(nslyr,kind=dbl_kind)
469 : enddo
470 104910624 : vbrin(n) = vicen(n)
471 122395728 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vbrin(n)
472 : enddo ! n
473 :
474 : call column_sum (ncat, &
475 17485104 : vicen, vice_final)
476 17485104 : if (icepack_warnings_aborted(subname)) return
477 : call column_sum (ncat, &
478 17485104 : vsnon, vsno_final)
479 17485104 : if (icepack_warnings_aborted(subname)) return
480 : call column_sum (ncat, &
481 17485104 : eicen, eice_final)
482 17485104 : if (icepack_warnings_aborted(subname)) return
483 : call column_sum (ncat, &
484 17485104 : esnon, esno_final)
485 17485104 : if (icepack_warnings_aborted(subname)) return
486 : call column_sum (ncat, &
487 17485104 : sicen, sice_final)
488 17485104 : if (icepack_warnings_aborted(subname)) return
489 : call column_sum (ncat, &
490 17485104 : vbrin, vbri_final)
491 17485104 : if (icepack_warnings_aborted(subname)) return
492 :
493 17485104 : vsno_final = vsno_final + msnow_mlt/rhos
494 17485104 : esno_final = esno_final + esnow_mlt
495 :
496 17485104 : fieldid = subname//':vice'
497 : call column_conservation_check (fieldid, &
498 : vice_init, vice_final, & ! LCOV_EXCL_LINE
499 17485104 : puny)
500 17485104 : if (icepack_warnings_aborted(subname)) return
501 17485104 : fieldid = subname//':vsno'
502 : call column_conservation_check (fieldid, &
503 : vsno_init, vsno_final, & ! LCOV_EXCL_LINE
504 17485104 : puny)
505 17485104 : if (icepack_warnings_aborted(subname)) return
506 17485104 : fieldid = subname//':eice'
507 : call column_conservation_check (fieldid, &
508 : eice_init, eice_final, & ! LCOV_EXCL_LINE
509 17485104 : puny*Lfresh*rhoi)
510 17485104 : if (icepack_warnings_aborted(subname)) return
511 17485104 : fieldid = subname//':esno'
512 : call column_conservation_check (fieldid, &
513 : esno_init, esno_final, & ! LCOV_EXCL_LINE
514 17485104 : puny*Lfresh*rhos)
515 17485104 : if (icepack_warnings_aborted(subname)) return
516 17485104 : fieldid = subname//':sice'
517 : call column_conservation_check (fieldid, &
518 : sice_init, sice_final, & ! LCOV_EXCL_LINE
519 17485104 : puny)
520 17485104 : if (icepack_warnings_aborted(subname)) return
521 17485104 : fieldid = subname//':vbrin'
522 : call column_conservation_check (fieldid, &
523 : vbri_init, vbri_final, & ! LCOV_EXCL_LINE
524 17485104 : puny*c10)
525 17485104 : if (icepack_warnings_aborted(subname)) return
526 :
527 : endif ! conserv_check
528 :
529 : !-----------------------------------------------------------------
530 : ! Compute ridging diagnostics.
531 : !-----------------------------------------------------------------
532 :
533 1498353240 : dti = c1/dt
534 :
535 1498353240 : if (present(dardg1dt)) then
536 1498353240 : dardg1dt = ardg1*dti
537 : endif
538 1498353240 : if (present(dardg2dt)) then
539 1498353240 : dardg2dt = ardg2*dti
540 : endif
541 1498353240 : if (present(dvirdgdt)) then
542 1498353240 : dvirdgdt = virdg*dti
543 : endif
544 1498353240 : if (present(opening)) then
545 1498353240 : opening = aopen*dti
546 : endif
547 1498353240 : if (present(dardg1ndt)) then
548 8967638592 : do n = 1, ncat
549 8967638592 : dardg1ndt(n) = ardg1n(n)*dti
550 : enddo
551 : endif
552 1498353240 : if (present(dardg2ndt)) then
553 8967638592 : do n = 1, ncat
554 8967638592 : dardg2ndt(n) = ardg2n(n)*dti
555 : enddo
556 : endif
557 1498353240 : if (present(dvirdgndt)) then
558 8967638592 : do n = 1, ncat
559 8967638592 : dvirdgndt(n) = virdgn(n)*dti
560 : enddo
561 : endif
562 1498353240 : if (present(araftn)) then
563 8967638592 : do n = 1, ncat
564 8967638592 : araftn(n) = mraftn(n)*ardg2n(n)
565 : ! araftn(n) = mraftn(n)*ardg1n(n)*p5
566 : enddo
567 : endif
568 1498353240 : if (present(vraftn)) then
569 8967638592 : do n = 1, ncat
570 8967638592 : vraftn(n) = mraftn(n)*virdgn(n)
571 : enddo
572 : endif
573 :
574 : !-----------------------------------------------------------------
575 : ! Update fresh water and heat fluxes due to snow melt.
576 : !-----------------------------------------------------------------
577 :
578 : ! use thermodynamic time step (ndtd*dt here) to average properly
579 1498353240 : dti = c1/(ndtd*dt)
580 :
581 1498353240 : if (present(fresh)) then
582 1498353240 : fresh = fresh + msnow_mlt*dti
583 : endif
584 1498353240 : if (present(fhocn)) then
585 1498353240 : fhocn = fhocn + esnow_mlt*dti
586 : endif
587 1498353240 : if (present(faero_ocn)) then
588 2943906768 : do it = 1, n_aero
589 2891107056 : faero_ocn(it) = faero_ocn(it) + maero(it)*dti
590 : enddo
591 : endif
592 1498353240 : if (present(flux_bio)) then
593 2453087592 : do it = 1, nbtrcr
594 1007534064 : flux_bio(it) = flux_bio(it) + mbio(it)*dti
595 : enddo
596 : endif
597 1498353240 : if (present(fiso_ocn)) then
598 1496431800 : if (tr_iso) then
599 : ! check size fiso_ocn vs n_iso ???
600 58411776 : do it = 1, n_iso
601 58411776 : fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti
602 : enddo
603 : endif
604 : endif
605 1498353240 : if (present(fpond)) then
606 1498353240 : fpond = fpond - mpond ! units change later
607 : endif
608 :
609 : !-----------------------------------------------------------------
610 : ! Check for fractional ice area > 1.
611 : !-----------------------------------------------------------------
612 :
613 1498353240 : if (abs(asum - c1) > puny) then
614 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
615 0 : call icepack_warnings_add(subname//" total area > 1" )
616 :
617 0 : write(warnstr,*) ' '
618 0 : call icepack_warnings_add(warnstr)
619 0 : write(warnstr,*) subname, 'Ridging error: total area > 1'
620 0 : call icepack_warnings_add(warnstr)
621 0 : write(warnstr,*) subname, 'area:', asum
622 0 : call icepack_warnings_add(warnstr)
623 0 : write(warnstr,*) subname, 'n, aicen:'
624 0 : call icepack_warnings_add(warnstr)
625 0 : write(warnstr,*) subname, 0, aice0
626 0 : call icepack_warnings_add(warnstr)
627 0 : do n = 1, ncat
628 0 : write(warnstr,*) subname, n, aicen(n)
629 0 : call icepack_warnings_add(warnstr)
630 : enddo
631 0 : return
632 : endif
633 :
634 : end subroutine ridge_ice
635 :
636 : !=======================================================================
637 :
638 : ! Find the total area of ice plus open water in each grid cell.
639 : !
640 : ! This is similar to the aggregate_area subroutine except that the
641 : ! total area can be greater than 1, so the open water area is
642 : ! included in the sum instead of being computed as a residual.
643 : !
644 : ! author: William H. Lipscomb, LANL
645 :
646 3339652841 : subroutine asum_ridging (aicen, aice0, asum)
647 :
648 : real (kind=dbl_kind), dimension (:), intent(in) :: &
649 : aicen ! concentration of ice in each category
650 :
651 : real (kind=dbl_kind), intent(in) :: &
652 : aice0 ! concentration of open water
653 :
654 : real (kind=dbl_kind), intent(out):: &
655 : asum ! sum of ice and open water area
656 :
657 : ! local variables
658 :
659 : integer (kind=int_kind) :: n
660 :
661 : character(len=*),parameter :: subname='(asum_ridging)'
662 :
663 3339652841 : asum = aice0
664 20006559342 : do n = 1, ncat
665 20006559342 : asum = asum + aicen(n)
666 : enddo
667 :
668 3339652841 : end subroutine asum_ridging
669 :
670 : !=======================================================================
671 :
672 : ! Initialize arrays, compute area of closing and opening
673 : !
674 : ! author: William H. Lipscomb, LANL
675 :
676 1498353240 : subroutine ridge_prep (dt, hin_max, &
677 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
678 : asum, closing_net, & ! LCOV_EXCL_LINE
679 : divu_adv, opning)
680 :
681 : real (kind=dbl_kind), intent(in) :: &
682 : dt ! time step (s)
683 :
684 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
685 : hin_max ! category limits (m)
686 :
687 : real (kind=dbl_kind), intent(in) :: &
688 : rdg_conv , & ! normalized energy dissipation due to convergence (1/s) ! LCOV_EXCL_LINE
689 : rdg_shear ! normalized energy dissipation due to shear (1/s)
690 :
691 : real (kind=dbl_kind), intent(inout):: &
692 : asum ! sum of ice and open water area
693 :
694 : real (kind=dbl_kind), intent(out):: &
695 : closing_net, & ! net rate at which area is removed (1/s) ! LCOV_EXCL_LINE
696 : divu_adv , & ! divu as implied by transport scheme (1/s) ! LCOV_EXCL_LINE
697 : opning ! rate of opening due to divergence/shear
698 :
699 : ! local variables
700 :
701 : real (kind=dbl_kind), parameter :: &
702 : big = 1.0e+8_dbl_kind
703 :
704 : character(len=*),parameter :: subname='(ridge_prep)'
705 :
706 : ! Set hin_max(ncat) to a big value to ensure that all ridged ice
707 : ! is thinner than hin_max(ncat).
708 1498353240 : hin_max(ncat) = big
709 :
710 : !-----------------------------------------------------------------
711 : ! Compute the net rate of closing due to convergence
712 : ! and shear, based on Flato and Hibler (1995).
713 : !
714 : ! For the elliptical yield curve:
715 : ! rdg_conv = -min (divu, 0)
716 : ! rdg_shear = (1/2) * (Delta - abs(divu))
717 : ! Note that the shear term also accounts for divergence.
718 : !
719 : ! The energy dissipation rate is equal to the net closing rate
720 : ! times the ice strength.
721 : !
722 : ! NOTE: The NET closing rate is equal to the rate that open water
723 : ! area is removed, plus the rate at which ice area is removed by
724 : ! ridging, minus the rate at which area is added in new ridges.
725 : ! The GROSS closing rate is equal to the first two terms (open
726 : ! water closing and thin ice ridging) without the third term
727 : ! (thick, newly ridged ice).
728 : !
729 : ! rdg_conv is calculated differently in EAP (update_ice_rdg) and
730 : ! represents closing_net directly. In that case, rdg_shear=0.
731 : !-----------------------------------------------------------------
732 :
733 1498353240 : closing_net = Cs*rdg_shear + rdg_conv
734 :
735 : !-----------------------------------------------------------------
736 : ! Compute divu_adv, the divergence rate given by the transport/
737 : ! advection scheme, which may not be equal to divu as computed
738 : ! from the velocity field.
739 : !
740 : ! If divu_adv < 0, make sure the closing rate is large enough
741 : ! to give asum = 1.0 after ridging.
742 : !-----------------------------------------------------------------
743 :
744 1498353240 : divu_adv = (c1-asum) / dt
745 :
746 1498353240 : if (divu_adv < c0) closing_net = max(closing_net, -divu_adv)
747 :
748 : !-----------------------------------------------------------------
749 : ! Compute the (non-negative) opening rate that will give
750 : ! asum = 1.0 after ridging.
751 : !-----------------------------------------------------------------
752 :
753 1498353240 : opning = closing_net + divu_adv
754 :
755 1498353240 : end subroutine ridge_prep
756 :
757 : !=======================================================================
758 :
759 : ! Compute the thickness distribution of the ice and open water
760 : ! participating in ridging and of the resulting ridges.
761 : !
762 : ! This version includes new options for ridging participation and
763 : ! redistribution.
764 : ! The new participation scheme (krdg_partic = 1) improves stability
765 : ! by increasing the time scale for large changes in ice strength.
766 : ! The new exponential redistribution function (krdg_redist = 1) improves
767 : ! agreement between ITDs of modeled and observed ridges.
768 : !
769 : ! author: William H. Lipscomb, LANL
770 : !
771 : ! 2006: Changed subroutine name to ridge_itd
772 : ! Added new options for ridging participation and redistribution.
773 :
774 1841299601 : subroutine ridge_itd (aice0, &
775 1841299601 : aicen, vicen, & ! LCOV_EXCL_LINE
776 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
777 : mu_rdg, & ! LCOV_EXCL_LINE
778 1841299601 : aksum, apartic, & ! LCOV_EXCL_LINE
779 1841299601 : hrmin, hrmax, & ! LCOV_EXCL_LINE
780 1841299601 : hrexp, krdg, & ! LCOV_EXCL_LINE
781 1841299601 : aparticn, krdgn, & ! LCOV_EXCL_LINE
782 1841299601 : mraft)
783 :
784 : real (kind=dbl_kind), intent(in) :: &
785 : mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) ! LCOV_EXCL_LINE
786 : aice0 ! concentration of open water
787 :
788 : real (kind=dbl_kind), dimension (:), intent(in) :: &
789 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
790 : vicen ! volume per unit area of ice (m)
791 :
792 : integer (kind=int_kind), intent(in) :: &
793 : krdg_partic , & ! selects participation function ! LCOV_EXCL_LINE
794 : krdg_redist ! selects redistribution function
795 :
796 : real (kind=dbl_kind), intent(out):: &
797 : aksum ! ratio of area removed to area ridged
798 :
799 : real (kind=dbl_kind), dimension (0:ncat), intent(out) :: &
800 : apartic ! participation function; fraction of ridging
801 : ! and closing associated w/ category n
802 :
803 : real (kind=dbl_kind), dimension (:), intent(out) :: &
804 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
805 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
806 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
807 : krdg ! mean ridge thickness/thickness of ridging ice
808 :
809 : ! diagnostic, category values
810 : real (kind=dbl_kind), dimension(:), intent(out), optional :: &
811 : aparticn, & ! participation function ! LCOV_EXCL_LINE
812 : krdgn ! mean ridge thickness/thickness of ridging ice
813 :
814 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
815 : mraft ! rafting ice mask
816 :
817 : ! local variables
818 :
819 : integer (kind=int_kind) :: &
820 : n ! thickness category index
821 :
822 : real (kind=dbl_kind), parameter :: &
823 : Gstari = c1/Gstar, & ! LCOV_EXCL_LINE
824 : astari = c1/astar
825 :
826 : real (kind=dbl_kind), dimension(-1:ncat) :: &
827 3682599202 : Gsum ! Gsum(n) = sum of areas in categories 0 to n
828 :
829 : real (kind=dbl_kind) :: &
830 : work ! temporary work variable
831 :
832 : real (kind=dbl_kind) :: &
833 : hi , & ! ice thickness for each cat (m) ! LCOV_EXCL_LINE
834 : hrmean , & ! mean ridge thickness (m) ! LCOV_EXCL_LINE
835 : xtmp ! temporary variable
836 :
837 : character(len=*),parameter :: subname='(ridge_itd)'
838 :
839 : !-----------------------------------------------------------------
840 : ! Initialize
841 : !-----------------------------------------------------------------
842 :
843 14721519952 : Gsum (-1:ncat) = c0 ! initialize
844 :
845 1841299601 : Gsum (-1) = c0 ! by definition
846 1841299601 : if (aice0 > puny) then
847 1588184065 : Gsum(0) = aice0
848 : else
849 253115536 : Gsum(0) = Gsum(-1)
850 : endif
851 1841299601 : apartic(0) = c0
852 :
853 11038920750 : do n = 1, ncat
854 9197621149 : apartic(n) = c0
855 9197621149 : hrmin (n) = c0
856 9197621149 : hrmax (n) = c0
857 9197621149 : hrexp (n) = c0
858 9197621149 : krdg (n) = c1
859 :
860 : !-----------------------------------------------------------------
861 : ! Compute the thickness distribution of ice participating in ridging.
862 : !-----------------------------------------------------------------
863 :
864 : !-----------------------------------------------------------------
865 : ! First compute the cumulative thickness distribution function Gsum,
866 : ! where Gsum(n) is the fractional area in categories 0 to n.
867 : ! Ignore categories with very small areas.
868 : !-----------------------------------------------------------------
869 :
870 11038920750 : if (aicen(n) > puny) then
871 3288082833 : Gsum(n) = Gsum(n-1) + aicen(n)
872 : else
873 5909538316 : Gsum(n) = Gsum(n-1)
874 : endif
875 : enddo
876 :
877 : ! normalize
878 :
879 1841299601 : if (Gsum(ncat) > c0) then
880 1841299601 : work = c1 / Gsum(ncat)
881 12880220351 : do n = 0, ncat
882 12880220351 : Gsum(n) = Gsum(n) * work
883 : enddo
884 : endif
885 :
886 : !-----------------------------------------------------------------
887 : ! Compute the participation function apartic; this is analogous to
888 : ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
889 : !
890 : ! area lost from category n due to ridging/closing
891 : ! apartic(n) = --------------------------------------------------
892 : ! total area lost due to ridging/closing
893 : !
894 : !-----------------------------------------------------------------
895 :
896 1841299601 : if (krdg_partic == 0) then ! Thornike et al. 1975 formulation
897 :
898 : !-----------------------------------------------------------------
899 : ! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
900 : ! The expressions for apartic are found by integrating b(h)g(h) between
901 : ! the category boundaries.
902 : !-----------------------------------------------------------------
903 :
904 282317088 : do n = 0, ncat
905 282317088 : if (Gsum(n) < Gstar) then
906 : apartic(n) = Gstari*(Gsum(n ) - Gsum(n-1)) * &
907 34906415 : (c2 - Gstari*(Gsum(n-1) + Gsum(n )))
908 197197969 : elseif (Gsum(n-1) < Gstar) then
909 : apartic(n) = Gstari*(Gstar - Gsum(n-1)) * &
910 50212704 : (c2 - Gstari*(Gstar + Gsum(n-1)))
911 : endif
912 : enddo ! n
913 :
914 1791086897 : elseif (krdg_partic==1) then ! exponential dependence on G(h)
915 :
916 : !-----------------------------------------------------------------
917 : ! b(h) = exp(-G(h)/astar)
918 : ! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)].
919 : ! The expression for apartic is found by integrating b(h)g(h)
920 : ! between the category boundaries.
921 : !-----------------------------------------------------------------
922 :
923 : ! precompute exponential terms using Gsum as work array
924 1791086897 : xtmp = c1 / (c1 - exp(-astari))
925 1791086897 : Gsum(-1) = exp(-Gsum(-1)*astari) * xtmp
926 12597903263 : do n = 0, ncat
927 10806816366 : Gsum(n) = exp(-Gsum(n)*astari) * xtmp
928 12597903263 : apartic(n) = Gsum(n-1) - Gsum(n)
929 : enddo ! n
930 :
931 : endif ! krdg_partic
932 :
933 : !-----------------------------------------------------------------
934 : ! Compute variables related to ITD of ridged ice:
935 : !
936 : ! krdg = mean ridge thickness / thickness of ridging ice
937 : ! hrmin = min ridge thickness
938 : ! hrmax = max ridge thickness (krdg_redist = 0)
939 : ! hrexp = ridge e-folding scale (krdg_redist = 1)
940 : !----------------------------------------------------------------
941 :
942 1841299601 : if (krdg_redist == 0) then ! Hibler 1980 formulation
943 :
944 : !-----------------------------------------------------------------
945 : ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
946 : !
947 : ! This parameterization is a modified version of Hibler (1980).
948 : ! In the original paper the min ridging thickness is hrmin = 2*hi,
949 : ! and the max thickness is hrmax = 2*sqrt(hi*Hstar).
950 : !
951 : ! Here the min thickness is hrmin = min(2*hi, hi+maxraft),
952 : ! so thick ridging ice is not required to raft.
953 : !
954 : !-----------------------------------------------------------------
955 :
956 232104384 : do n = 1, ncat
957 232104384 : if (aicen(n) > puny) then
958 101512155 : hi = vicen(n) / aicen(n)
959 101512155 : hrmin(n) = min(c2*hi, hi + maxraft)
960 101512155 : hrmax(n) = c2*sqrt(Hstar*hi)
961 101512155 : hrmax(n) = max(hrmax(n), hrmin(n)+puny)
962 101512155 : hrmean = p5 * (hrmin(n) + hrmax(n))
963 101512155 : krdg(n) = hrmean / hi
964 :
965 : ! diagnostic rafting mask not implemented
966 : endif
967 : enddo ! n
968 :
969 : else ! krdg_redist = 1; exponential redistribution
970 :
971 : !-----------------------------------------------------------------
972 : ! The ridge ITD is a negative exponential:
973 : !
974 : ! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin
975 : !
976 : ! where hrmin is the minimum thickness of ridging ice and
977 : ! hrexp is the e-folding thickness.
978 : !
979 : ! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
980 : ! That is, the minimum ridge thickness results from rafting,
981 : ! unless the ice is thicker than maxraft.
982 : !
983 : ! Also, assume that hrexp = mu_rdg*sqrt(hi).
984 : ! The parameter mu_rdg is tuned to give e-folding scales mostly
985 : ! in the range 2-4 m as observed by upward-looking sonar.
986 : !
987 : ! Values of mu_rdg in the right column give ice strengths
988 : ! roughly equal to values of Hstar in the left column
989 : ! (within ~10 kN/m for typical ITDs):
990 : !
991 : ! Hstar mu_rdg
992 : !
993 : ! 25 3.0
994 : ! 50 4.0
995 : ! 75 5.0
996 : ! 100 6.0
997 : !-----------------------------------------------------------------
998 :
999 10806816366 : do n = 1, ncat
1000 10806816366 : if (aicen(n) > puny) then
1001 3186570678 : hi = vicen(n) / aicen(n)
1002 3186570678 : hi = max(hi,puny)
1003 3186570678 : hrmin(n) = min(c2*hi, hi + maxraft)
1004 3186570678 : hrexp(n) = mu_rdg * sqrt(hi)
1005 3186570678 : krdg(n) = (hrmin(n) + hrexp(n)) / hi
1006 :
1007 : !echmod: check computational efficiency
1008 : ! diagnostic rafting mask
1009 3186570678 : if (present(mraft)) then
1010 1524526828 : mraft(n) = max(c0, sign(c1, hi+maxraft-hrmin(n)))
1011 1524526828 : xtmp = mraft(n)*((c2*hi+hrexp(n))/hi - krdg(n))
1012 1524526828 : mraft(n) = max(c0, sign(c1, puny-abs(xtmp)))
1013 : endif
1014 : endif
1015 : enddo
1016 :
1017 : endif ! krdg_redist
1018 :
1019 : !----------------------------------------------------------------
1020 : ! Compute aksum = net ice area removed / total area participating.
1021 : ! For instance, if a unit area of ice with h = 1 participates in
1022 : ! ridging to form a ridge with a = 1/3 and h = 3, then
1023 : ! aksum = 1 - 1/3 = 2/3.
1024 : !----------------------------------------------------------------
1025 :
1026 1841299601 : aksum = apartic(0) ! area participating = area removed
1027 :
1028 11038920750 : do n = 1, ncat
1029 : ! area participating > area removed
1030 11038920750 : aksum = aksum + apartic(n) * (c1 - c1/krdg(n))
1031 : enddo
1032 :
1033 : ! diagnostics
1034 1841299601 : if (present(aparticn)) then
1035 8967638694 : do n = 1, ncat
1036 8967638694 : aparticn(n) = apartic(n)
1037 : enddo
1038 : endif
1039 1841299601 : if (present(krdgn)) then
1040 8967638694 : do n = 1, ncat
1041 8967638694 : krdgn(n) = krdg(n)
1042 : enddo
1043 : endif
1044 :
1045 1841299601 : end subroutine ridge_itd
1046 :
1047 : !=======================================================================
1048 :
1049 : ! Remove area, volume, and energy from each ridging category
1050 : ! and add to thicker ice categories.
1051 : !
1052 : ! Tracers: Ridging conserves ice volume and therefore conserves volume
1053 : ! tracers. It does not conserve ice area, and therefore a portion of area
1054 : ! tracers are lost (corresponding to the net closing). Area tracers on
1055 : ! ice that participates in ridging are carried onto the resulting ridged
1056 : ! ice (except the portion that are lost due to closing). Therefore,
1057 : ! tracers must be decremented if they are lost to the ocean during ridging
1058 : ! (e.g. snow, ponds) or if they are being carried only on the level ice
1059 : ! area.
1060 : !
1061 : ! author: William H. Lipscomb, LANL
1062 :
1063 2996706514 : subroutine ridge_shift (dt, hin_max, &
1064 2996706514 : aicen, trcrn, & ! LCOV_EXCL_LINE
1065 1498353257 : vicen, vsnon, & ! LCOV_EXCL_LINE
1066 1498353257 : aice0, trcr_depend, & ! LCOV_EXCL_LINE
1067 1498353257 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
1068 1498353257 : nt_strata, krdg_redist, & ! LCOV_EXCL_LINE
1069 1498353257 : aksum, apartic, & ! LCOV_EXCL_LINE
1070 1498353257 : hrmin, hrmax, & ! LCOV_EXCL_LINE
1071 2996706514 : hrexp, krdg, & ! LCOV_EXCL_LINE
1072 : closing_net, opning, & ! LCOV_EXCL_LINE
1073 : ardg1, ardg2, & ! LCOV_EXCL_LINE
1074 : virdg, aopen, & ! LCOV_EXCL_LINE
1075 1498353257 : ardg1nn, ardg2nn, & ! LCOV_EXCL_LINE
1076 1498353257 : virdgnn, & ! LCOV_EXCL_LINE
1077 : msnow_mlt, esnow_mlt, & ! LCOV_EXCL_LINE
1078 1498353257 : maero, miso, & ! LCOV_EXCL_LINE
1079 : mpond, Tf, & ! LCOV_EXCL_LINE
1080 1498353257 : aredistn, vredistn, & ! LCOV_EXCL_LINE
1081 1498353257 : mbio)
1082 :
1083 : integer (kind=int_kind), intent(in) :: &
1084 : krdg_redist ! selects redistribution function
1085 :
1086 : real (kind=dbl_kind), intent(in) :: &
1087 : dt ! time step (s)
1088 :
1089 : real (kind=dbl_kind), intent(in) :: &
1090 : Tf ! freezing temperature
1091 :
1092 : integer (kind=int_kind), dimension (:), intent(in) :: &
1093 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
1094 : n_trcr_strata ! number of underlying tracer layers
1095 :
1096 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1097 : trcr_base ! = 0 or 1 depending on tracer dependency
1098 : ! argument 2: (1) aice, (2) vice, (3) vsno
1099 :
1100 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
1101 : nt_strata ! indices of underlying tracer layers
1102 :
1103 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
1104 : hin_max ! category limits (m)
1105 :
1106 : real (kind=dbl_kind), intent(inout) :: &
1107 : aice0 ! concentration of open water
1108 :
1109 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1110 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1111 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1112 : vsnon ! volume per unit area of snow (m)
1113 :
1114 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1115 : trcrn ! ice tracers
1116 :
1117 : real (kind=dbl_kind), intent(in) :: &
1118 : aksum ! ratio of area removed to area ridged
1119 :
1120 : real (kind=dbl_kind), dimension (0:ncat), intent(in) :: &
1121 : apartic ! participation function; fraction of ridging
1122 : ! and closing associated w/ category n
1123 :
1124 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1125 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
1126 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
1127 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
1128 : krdg ! mean ridge thickness/thickness of ridging ice
1129 :
1130 : real (kind=dbl_kind), intent(inout) :: &
1131 : closing_net, & ! net rate at which area is removed (1/s) ! LCOV_EXCL_LINE
1132 : opning , & ! rate of opening due to divergence/shear (1/s) ! LCOV_EXCL_LINE
1133 : ardg1 , & ! fractional area loss by ridging ice ! LCOV_EXCL_LINE
1134 : ardg2 , & ! fractional area gain by new ridges ! LCOV_EXCL_LINE
1135 : virdg , & ! ice volume ridged (m) ! LCOV_EXCL_LINE
1136 : aopen ! area opened due to divergence/shear
1137 :
1138 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1139 : ardg1nn , & ! area of ice ridged ! LCOV_EXCL_LINE
1140 : ardg2nn , & ! area of new ridges ! LCOV_EXCL_LINE
1141 : virdgnn ! ridging ice volume
1142 :
1143 : real (kind=dbl_kind), intent(inout) :: &
1144 : msnow_mlt , & ! mass of snow added to ocean (kg m-2) ! LCOV_EXCL_LINE
1145 : esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) ! LCOV_EXCL_LINE
1146 : mpond ! mass of pond added to ocean (kg m-2)
1147 :
1148 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1149 : maero ! aerosol mass added to ocean (kg m-2)
1150 :
1151 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1152 : miso ! isotope mass added to ocean (kg m-2)
1153 :
1154 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1155 : mbio ! biology and zaerosol mass added to ocean (kg m-2)
1156 :
1157 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
1158 : aredistn , & ! redistribution function: fraction of new ridge area ! LCOV_EXCL_LINE
1159 : vredistn ! redistribution function: fraction of new ridge volume
1160 :
1161 : ! local variables
1162 :
1163 : integer (kind=int_kind) :: &
1164 : n, nr , & ! thickness category indices ! LCOV_EXCL_LINE
1165 : k , & ! ice layer index ! LCOV_EXCL_LINE
1166 : it , & ! tracer index ! LCOV_EXCL_LINE
1167 : ntr , & ! tracer index ! LCOV_EXCL_LINE
1168 : itl ! loop index
1169 :
1170 : real (kind=dbl_kind), dimension (ncat) :: &
1171 1498353257 : aicen_init , & ! ice area before ridging ! LCOV_EXCL_LINE
1172 2996706514 : vicen_init , & ! ice volume before ridging ! LCOV_EXCL_LINE
1173 2996706514 : vsnon_init ! snow volume before ridging
1174 :
1175 : real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
1176 1498353257 : atrcrn ! aicen*trcrn
1177 :
1178 : real (kind=dbl_kind), dimension(3) :: &
1179 : trfactor ! base quantity on which tracers are carried
1180 :
1181 : real (kind=dbl_kind) :: &
1182 : work , & ! temporary variable ! LCOV_EXCL_LINE
1183 : expL_arg , & ! temporary exp arg values ! LCOV_EXCL_LINE
1184 : expR_arg , & ! temporary exp arg values ! LCOV_EXCL_LINE
1185 : closing_gross ! rate at which area removed, not counting
1186 : ! area of new ridges
1187 :
1188 : ! ECH note: the following arrays only need be defined on iridge cells
1189 : real (kind=dbl_kind) :: &
1190 : afrac , & ! fraction of category area ridged ! LCOV_EXCL_LINE
1191 : ardg1n , & ! area of ice ridged ! LCOV_EXCL_LINE
1192 : ardg2n , & ! area of new ridges ! LCOV_EXCL_LINE
1193 : virdgn , & ! ridging ice volume ! LCOV_EXCL_LINE
1194 : vsrdgn , & ! ridging snow volume ! LCOV_EXCL_LINE
1195 : dhr , & ! hrmax - hrmin ! LCOV_EXCL_LINE
1196 : dhr2 , & ! hrmax^2 - hrmin^2 ! LCOV_EXCL_LINE
1197 : farea , & ! fraction of new ridge area going to nr ! LCOV_EXCL_LINE
1198 : fvol ! fraction of new ridge volume going to nr
1199 :
1200 : real (kind=dbl_kind) :: &
1201 : esrdgn ! ridging snow energy
1202 :
1203 : real (kind=dbl_kind) :: &
1204 : hi1 , & ! thickness of ridging ice ! LCOV_EXCL_LINE
1205 : hexp , & ! ridge e-folding thickness ! LCOV_EXCL_LINE
1206 : hL, hR , & ! left and right limits of integration ! LCOV_EXCL_LINE
1207 : expL, expR , & ! exponentials involving hL, hR ! LCOV_EXCL_LINE
1208 : tmpfac , & ! factor by which opening/closing rates are cut ! LCOV_EXCL_LINE
1209 : wk1 , & ! work variable ! LCOV_EXCL_LINE
1210 : dzssl , & ! fraction of snow surface biotracers ! LCOV_EXCL_LINE
1211 : dzint ! fraction of interior snow biotracers
1212 :
1213 : character(len=*),parameter :: subname='(ridge_shift)'
1214 :
1215 8967638694 : do n = 1, ncat
1216 :
1217 : !-----------------------------------------------------------------
1218 : ! Save initial state variables
1219 : !-----------------------------------------------------------------
1220 :
1221 7469285437 : aicen_init(n) = aicen(n)
1222 7469285437 : vicen_init(n) = vicen(n)
1223 7469285437 : vsnon_init(n) = vsnon(n)
1224 :
1225 : !-----------------------------------------------------------------
1226 : ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
1227 : !-----------------------------------------------------------------
1228 :
1229 >24505*10^7 : do it = 1, ntrcr
1230 : atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
1231 : + trcr_base(it,2) * vicen(n) & ! LCOV_EXCL_LINE
1232 >23608*10^7 : + trcr_base(it,3) * vsnon(n))
1233 >24355*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1234 56790078440 : do itl = 1, n_trcr_strata(it)
1235 35441343713 : ntr = nt_strata(it,itl)
1236 56790078440 : atrcrn(it,n) = atrcrn(it,n) * trcrn(ntr,n)
1237 : enddo
1238 : endif
1239 : enddo
1240 :
1241 : enddo ! ncat
1242 :
1243 : !-----------------------------------------------------------------
1244 : ! Based on the ITD of ridging and ridged ice, convert the net
1245 : ! closing rate to a gross closing rate.
1246 : ! NOTE: 0 < aksum <= 1
1247 : !-----------------------------------------------------------------
1248 :
1249 1498353257 : closing_gross = closing_net / aksum
1250 :
1251 : !-----------------------------------------------------------------
1252 : ! Reduce the closing rate if more than 100% of the open water
1253 : ! would be removed. Reduce the opening rate proportionately.
1254 : !-----------------------------------------------------------------
1255 :
1256 1498353257 : if (apartic(0) > c0) then
1257 1383018030 : wk1 = apartic(0) * closing_gross * dt
1258 1383018030 : if (wk1 > aice0) then
1259 17 : tmpfac = aice0 / wk1
1260 17 : closing_gross = closing_gross * tmpfac
1261 17 : opning = opning * tmpfac
1262 : endif
1263 : endif
1264 :
1265 : !-----------------------------------------------------------------
1266 : ! Reduce the closing rate if more than 100% of any ice category
1267 : ! would be removed. Reduce the opening rate proportionately.
1268 : !-----------------------------------------------------------------
1269 8967638694 : do n = 1, ncat
1270 8967638694 : if (aicen(n) > puny .and. apartic(n) > c0) then
1271 1559214998 : wk1 = apartic(n) * closing_gross * dt
1272 1559214998 : if (wk1 > aicen(n)) then
1273 0 : tmpfac = aicen(n) / wk1
1274 0 : closing_gross = closing_gross * tmpfac
1275 0 : opning = opning * tmpfac
1276 : endif
1277 : endif
1278 : enddo ! n
1279 :
1280 : !-----------------------------------------------------------------
1281 : ! Compute change in open water area due to closing and opening.
1282 : !-----------------------------------------------------------------
1283 :
1284 1498353257 : aice0 = aice0 - apartic(0)*closing_gross*dt + opning*dt
1285 :
1286 1498353257 : if (aice0 < -puny) then
1287 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1288 0 : call icepack_warnings_add(subname//' Ridging error: aice0 < 0')
1289 0 : write(warnstr,*) subname, 'aice0:', aice0
1290 0 : call icepack_warnings_add(warnstr)
1291 0 : return
1292 :
1293 1498353257 : elseif (aice0 < c0) then ! roundoff error
1294 1 : aice0 = c0
1295 : endif
1296 :
1297 1498353257 : aopen = opning*dt ! optional diagnostic
1298 :
1299 : !-----------------------------------------------------------------
1300 : ! Compute the area, volume, and energy of ice ridging in each
1301 : ! category, along with the area of the resulting ridge.
1302 : !-----------------------------------------------------------------
1303 :
1304 8967638694 : do n = 1, ncat
1305 :
1306 : !-----------------------------------------------------------------
1307 : ! Identify grid cells with nonzero ridging
1308 : !-----------------------------------------------------------------
1309 :
1310 : if (aicen_init(n) > puny .and. apartic(n) > c0 &
1311 8967638694 : .and. closing_gross > c0) then
1312 :
1313 : !-----------------------------------------------------------------
1314 : ! Compute area of ridging ice (ardg1n) and of new ridge (ardg2n).
1315 : ! Make sure ridging fraction <=1. (Roundoff errors can give
1316 : ! ardg1 slightly greater than aicen.)
1317 : !-----------------------------------------------------------------
1318 :
1319 1426636056 : ardg1n = apartic(n)*closing_gross*dt
1320 :
1321 1426636056 : if (ardg1n > aicen_init(n) + puny) then
1322 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1323 0 : call icepack_warnings_add(subname//' Ridging error: ardg > aicen')
1324 0 : write(warnstr,*) subname, 'n, ardg, aicen:', &
1325 0 : n, ardg1n, aicen_init(n)
1326 0 : call icepack_warnings_add(warnstr)
1327 0 : return
1328 : else
1329 1426636056 : ardg1n = min(aicen_init(n), ardg1n)
1330 : endif
1331 :
1332 1426636056 : ardg2n = ardg1n / krdg(n)
1333 1426636056 : afrac = ardg1n / aicen_init(n)
1334 :
1335 : !-----------------------------------------------------------------
1336 : ! Subtract area, volume, and energy from ridging category n.
1337 : ! Note: Tracer values are unchanged.
1338 : !-----------------------------------------------------------------
1339 :
1340 1426636056 : virdgn = vicen_init(n) * afrac
1341 1426636056 : vsrdgn = vsnon_init(n) * afrac
1342 :
1343 1426636056 : aicen(n) = aicen(n) - ardg1n
1344 1426636056 : vicen(n) = vicen(n) - virdgn
1345 1426636056 : vsnon(n) = vsnon(n) - vsrdgn
1346 :
1347 : !-----------------------------------------------------------------
1348 : ! Increment ridging diagnostics
1349 : !-----------------------------------------------------------------
1350 :
1351 1426636056 : ardg1 = ardg1 + ardg1n
1352 1426636056 : ardg2 = ardg2 + ardg2n
1353 1426636056 : virdg = virdg + virdgn
1354 :
1355 1426636056 : ardg1nn(n) = ardg1n
1356 1426636056 : ardg2nn(n) = ardg2n
1357 1426636056 : virdgnn(n) = virdgn
1358 :
1359 : !-----------------------------------------------------------------
1360 : ! Place part of the snow and tracer lost by ridging into the ocean.
1361 : !-----------------------------------------------------------------
1362 :
1363 1426636056 : msnow_mlt = msnow_mlt + rhos*vsrdgn*(c1-fsnowrdg)
1364 :
1365 1426636056 : if (tr_aero) then
1366 38121522 : do it = 1, n_aero
1367 : maero(it) = maero(it) &
1368 : + vsrdgn*(c1-fsnowrdg) & ! LCOV_EXCL_LINE
1369 : *(trcrn(nt_aero +4*(it-1),n) & ! LCOV_EXCL_LINE
1370 38121522 : + trcrn(nt_aero+1+4*(it-1),n))
1371 : enddo
1372 : endif
1373 :
1374 1426636056 : if (tr_iso) then
1375 60818308 : do it = 1, n_iso
1376 : miso(it) = miso(it) + vsrdgn*(c1-fsnowrdg) &
1377 : * (trcrn(nt_isosno+it-1,n) & ! LCOV_EXCL_LINE
1378 60818308 : + trcrn(nt_isoice+it-1,n))
1379 : enddo
1380 : endif
1381 :
1382 1426636056 : if (z_tracers .and. nbtrcr > 0) then
1383 32929030 : dzssl = p5/real(nslyr,kind=dbl_kind)
1384 32929030 : dzint = c1-dzssl
1385 577485479 : do it = 1, nbtrcr
1386 : mbio(it) = mbio(it) + vsrdgn*(c1-fsnowrdg) &
1387 : * (trcrn(bio_index(it) + nblyr + 1,n) * dzssl & ! LCOV_EXCL_LINE
1388 577485479 : + trcrn(bio_index(it) + nblyr + 2,n) * dzint)
1389 : enddo
1390 : endif
1391 :
1392 1426636056 : if (tr_pond_topo) then
1393 : mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
1394 11322585 : * trcrn(nt_hpnd,n)
1395 : endif
1396 :
1397 : !-----------------------------------------------------------------
1398 : ! Compute quantities used to apportion ice among categories
1399 : ! in the nr loop below
1400 : !-----------------------------------------------------------------
1401 :
1402 1426636056 : dhr = hrmax(n) - hrmin(n)
1403 1426636056 : dhr2 = hrmax(n) * hrmax(n) - hrmin(n) * hrmin(n)
1404 :
1405 : !-----------------------------------------------------------------
1406 : ! Increment energy needed to melt snow in ocean.
1407 : ! Note that esnow_mlt < 0; the ocean must cool to melt snow.
1408 : !-----------------------------------------------------------------
1409 :
1410 3032687236 : do k = 1, nslyr
1411 : esrdgn = vsrdgn * trcrn(nt_qsno+k-1,n) &
1412 1606051180 : / real(nslyr,kind=dbl_kind)
1413 3032687236 : esnow_mlt = esnow_mlt + esrdgn*(c1-fsnowrdg)
1414 : enddo
1415 :
1416 : !-----------------------------------------------------------------
1417 : ! Subtract area- and volume-weighted tracers from category n.
1418 : !-----------------------------------------------------------------
1419 :
1420 41206954080 : do it = 1, ntrcr
1421 :
1422 39780318024 : trfactor(1) = trcr_base(it,1)*ardg1n
1423 39780318024 : trfactor(2) = trcr_base(it,2)*virdgn
1424 39780318024 : trfactor(3) = trcr_base(it,3)*vsrdgn
1425 :
1426 39780318024 : work = c0
1427 >15912*10^7 : do k = 1, 3
1428 >15912*10^7 : work = work + trfactor(k)*trcrn(it,n)
1429 : enddo
1430 39780318024 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1431 10465694476 : do itl = 1, n_trcr_strata(it)
1432 6535397755 : ntr = nt_strata(it,itl)
1433 10465694476 : work = work * trcrn(ntr,n)
1434 : enddo
1435 : endif
1436 41206954080 : atrcrn(it,n) = atrcrn(it,n) - work
1437 :
1438 : enddo ! ntrcr
1439 :
1440 : !-----------------------------------------------------------------
1441 : ! Add area, volume, and energy of new ridge to each category nr.
1442 : !-----------------------------------------------------------------
1443 :
1444 8605053881 : do nr = 1, ncat
1445 :
1446 7178417825 : if (krdg_redist == 0) then ! Hibler 1980 formulation
1447 :
1448 : !-----------------------------------------------------------------
1449 : ! Compute the fraction of ridged ice area and volume going to
1450 : ! thickness category nr.
1451 : !-----------------------------------------------------------------
1452 :
1453 139095865 : if (hrmin(n) >= hin_max(nr) .or. &
1454 : hrmax(n) <= hin_max(nr-1)) then
1455 37067090 : hL = c0
1456 37067090 : hR = c0
1457 : else
1458 102028775 : hL = max (hrmin(n), hin_max(nr-1))
1459 102028775 : hR = min (hrmax(n), hin_max(nr))
1460 : endif
1461 :
1462 139095865 : farea = (hR-hL) / dhr
1463 139095865 : fvol = (hR*hR - hL*hL) / dhr2
1464 :
1465 : else ! krdg_redist = 1; 2005 exponential formulation
1466 :
1467 : !-----------------------------------------------------------------
1468 : ! Compute the fraction of ridged ice area and volume going to
1469 : ! thickness category nr.
1470 : !-----------------------------------------------------------------
1471 :
1472 7039321960 : if (nr < ncat) then
1473 :
1474 5643027905 : hi1 = hrmin(n)
1475 5643027905 : hexp = hrexp(n)
1476 :
1477 5643027905 : if (hi1 >= hin_max(nr)) then
1478 3537364797 : farea = c0
1479 3537364797 : fvol = c0
1480 : else
1481 2105663108 : hL = max (hi1, hin_max(nr-1))
1482 2105663108 : hR = hin_max(nr)
1483 2105663108 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1484 2105663108 : expR_arg = min(((hR-hi1)/hexp),exp_argmax)
1485 2105663108 : expL = exp(-(expL_arg))
1486 2105663108 : expR = exp(-(expR_arg))
1487 2105663108 : farea = expL - expR
1488 : fvol = ((hL + hexp)*expL &
1489 2105663108 : - (hR + hexp)*expR) / (hi1 + hexp)
1490 : endif
1491 :
1492 : else ! nr = ncat
1493 :
1494 1396294055 : hi1 = hrmin(n)
1495 1396294055 : hexp = hrexp(n)
1496 :
1497 1396294055 : hL = max (hi1, hin_max(nr-1))
1498 1396294055 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1499 1396294055 : expL = exp(-(expL_arg))
1500 1396294055 : farea = expL
1501 1396294055 : fvol = (hL + hexp)*expL / (hi1 + hexp)
1502 :
1503 : endif ! nr < ncat
1504 :
1505 : ! diagnostics
1506 7039321960 : if (n ==1) then ! only for thinnest ridging ice
1507 1403422502 : if (present(aredistn)) then
1508 1403422502 : aredistn(nr) = farea*ardg2n
1509 : endif
1510 1403422502 : if (present(vredistn)) then
1511 1403422502 : vredistn(nr) = fvol*virdgn
1512 : endif
1513 : endif
1514 :
1515 : endif ! krdg_redist
1516 :
1517 : !-----------------------------------------------------------------
1518 : ! Transfer ice area, ice volume, and snow volume to category nr.
1519 : !-----------------------------------------------------------------
1520 :
1521 7178417825 : aicen(nr) = aicen(nr) + farea*ardg2n
1522 7178417825 : vicen(nr) = vicen(nr) + fvol *virdgn
1523 7178417825 : vsnon(nr) = vsnon(nr) + fvol *vsrdgn*fsnowrdg
1524 :
1525 : !-----------------------------------------------------------------
1526 : ! Transfer area-weighted and volume-weighted tracers to category nr.
1527 : ! Note: The global sum aicen*trcrn of ice area tracers
1528 : ! (trcr_depend = 0) is not conserved by ridging.
1529 : ! However, ridging conserves the global sum of volume
1530 : ! tracers (trcr_depend = 1 or 2).
1531 : ! Tracers associated with level ice, or that are otherwise lost
1532 : ! from ridging ice, are not transferred.
1533 : ! We assume that all pond water is lost from ridging ice.
1534 : !-----------------------------------------------------------------
1535 :
1536 >20872*10^7 : do it = 1, ntrcr
1537 :
1538 >20011*10^7 : if (it /= nt_alvl .and. it /= nt_vlvl) then
1539 >18591*10^7 : trfactor(1) = trcr_base(it,1)*ardg2n*farea
1540 >18591*10^7 : trfactor(2) = trcr_base(it,2)*virdgn*fvol
1541 >18591*10^7 : trfactor(3) = trcr_base(it,3)*vsrdgn*fvol*fsnowrdg
1542 : else
1543 14205569465 : trfactor(1) = c0
1544 14205569465 : trfactor(2) = c0
1545 14205569465 : trfactor(3) = c0
1546 : endif
1547 :
1548 >20011*10^7 : work = c0
1549 >80047*10^7 : do k = 1, 3
1550 >80047*10^7 : work = work + trfactor(k)*trcrn(it,n)
1551 : enddo
1552 >20011*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1553 52745995520 : do itl = 1, n_trcr_strata(it)
1554 32932279445 : ntr = nt_strata(it,itl)
1555 52745995520 : if (ntr == nt_fbri) then ! brine fraction only
1556 0 : work = work * trcrn(ntr,n)
1557 : else
1558 32932279445 : work = c0
1559 : endif
1560 : enddo
1561 : endif
1562 >20729*10^7 : atrcrn(it,nr) = atrcrn(it,nr) + work
1563 :
1564 : enddo ! ntrcr
1565 :
1566 : enddo ! nr (new ridges)
1567 :
1568 : endif ! nonzero ridging
1569 :
1570 : enddo ! n (ridging categories)
1571 :
1572 : !-----------------------------------------------------------------
1573 : ! Compute new tracers
1574 : !-----------------------------------------------------------------
1575 :
1576 8967638694 : do n = 1, ncat
1577 : call icepack_compute_tracers (trcr_depend, &
1578 : atrcrn(:,n), aicen(n), & ! LCOV_EXCL_LINE
1579 : vicen(n), vsnon(n), & ! LCOV_EXCL_LINE
1580 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
1581 7469285437 : nt_strata, trcrn(:,n), Tf)
1582 8967638694 : if (icepack_warnings_aborted(subname)) return
1583 : enddo
1584 :
1585 : end subroutine ridge_shift
1586 :
1587 : !=======================================================================
1588 : !autodocument_start icepack_ice_strength
1589 : ! Compute the strength of the ice pack, defined as the energy (J m-2)
1590 : ! dissipated per unit area removed from the ice pack under compression,
1591 : ! and assumed proportional to the change in potential energy caused
1592 : ! by ridging.
1593 : !
1594 : ! See Rothrock (1975) and Hibler (1980).
1595 : !
1596 : ! For simpler strength parameterization, see this reference:
1597 : ! Hibler, W. D. III, 1979: A dynamic-thermodynamic sea ice model,
1598 : ! J. Phys. Oceanog., 9, 817-846.
1599 : !
1600 : ! authors: William H. Lipscomb, LANL
1601 : ! Elizabeth C. Hunke, LANL
1602 :
1603 381893651 : subroutine icepack_ice_strength(aice, vice, &
1604 381893651 : aice0, aicen, & ! LCOV_EXCL_LINE
1605 381893651 : vicen, & ! LCOV_EXCL_LINE
1606 : strength)
1607 :
1608 : real (kind=dbl_kind), intent(in) :: &
1609 : aice , & ! concentration of ice ! LCOV_EXCL_LINE
1610 : vice , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1611 : aice0 ! concentration of open water
1612 :
1613 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1614 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1615 : vicen ! volume per unit area of ice (m)
1616 :
1617 : real (kind=dbl_kind), intent(inout) :: &
1618 : strength ! ice strength (N/m)
1619 :
1620 : !autodocument_end
1621 :
1622 : ! local variables
1623 :
1624 : real (kind=dbl_kind) :: &
1625 : asum , & ! sum of ice and open water area ! LCOV_EXCL_LINE
1626 : aksum ! ratio of area removed to area ridged
1627 :
1628 : real (kind=dbl_kind), dimension (0:ncat) :: &
1629 381893651 : apartic ! participation function; fraction of ridging
1630 : ! and closing associated w/ category n
1631 :
1632 : real (kind=dbl_kind), dimension (ncat) :: &
1633 763787302 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
1634 763787302 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
1635 763787302 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
1636 381893651 : krdg ! mean ridge thickness/thickness of ridging ice
1637 :
1638 : integer (kind=int_kind) :: &
1639 : n ! thickness category index
1640 :
1641 : real (kind=dbl_kind) :: &
1642 : hi , & ! ice thickness (m) ! LCOV_EXCL_LINE
1643 : h2rdg , & ! mean value of h^2 for new ridge ! LCOV_EXCL_LINE
1644 : dh2rdg ! change in mean value of h^2 per unit area
1645 : ! consumed by ridging
1646 :
1647 : character(len=*),parameter :: subname='(icepack_ice_strength)'
1648 :
1649 381893651 : if (kstrength == 1) then ! Rothrock 1975 formulation
1650 :
1651 : !-----------------------------------------------------------------
1652 : ! Compute thickness distribution of ridging and ridged ice.
1653 : !-----------------------------------------------------------------
1654 :
1655 342946344 : call asum_ridging (aicen, aice0, asum)
1656 342946344 : if (icepack_warnings_aborted(subname)) return
1657 :
1658 : call ridge_itd (aice0, &
1659 : aicen, vicen, & ! LCOV_EXCL_LINE
1660 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
1661 : mu_rdg, & ! LCOV_EXCL_LINE
1662 : aksum, apartic, & ! LCOV_EXCL_LINE
1663 : hrmin, hrmax, & ! LCOV_EXCL_LINE
1664 342946344 : hrexp, krdg)
1665 342946344 : if (icepack_warnings_aborted(subname)) return
1666 :
1667 : !-----------------------------------------------------------------
1668 : ! Compute ice strength based on change in potential energy,
1669 : ! as in Rothrock (1975)
1670 : !-----------------------------------------------------------------
1671 :
1672 342946344 : if (krdg_redist==0) then ! Hibler 1980 formulation
1673 :
1674 0 : do n = 1, ncat
1675 0 : if (aicen(n) > puny .and. apartic(n) > c0)then
1676 0 : hi = vicen(n) / aicen(n)
1677 : h2rdg = p333 * (hrmax(n)**3 - hrmin(n)**3) &
1678 0 : / (hrmax(n) - hrmin(n))
1679 0 : dh2rdg = -hi*hi + h2rdg/krdg(n)
1680 0 : strength = strength + apartic(n) * dh2rdg
1681 : endif ! aicen > puny
1682 : enddo ! n
1683 :
1684 342946344 : elseif (krdg_redist==1) then ! exponential formulation
1685 :
1686 2071282056 : do n = 1, ncat
1687 2071282056 : if (aicen(n) > puny .and. apartic(n) > c0) then
1688 1662043850 : hi = vicen(n) / aicen(n)
1689 : h2rdg = hrmin(n)*hrmin(n) &
1690 : + c2*hrmin(n)*hrexp(n) & ! LCOV_EXCL_LINE
1691 1662043850 : + c2*hrexp(n)*hrexp(n)
1692 1662043850 : dh2rdg = -hi*hi + h2rdg/krdg(n)
1693 1662043850 : strength = strength + apartic(n) * dh2rdg
1694 : endif
1695 : enddo ! n
1696 :
1697 : endif ! krdg_redist
1698 :
1699 342946344 : strength = Cf * Cp * strength / aksum
1700 : ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow)
1701 : ! Cf accounts for frictional dissipation
1702 :
1703 : else ! kstrength /= 1: Hibler (1979) form
1704 :
1705 : !-----------------------------------------------------------------
1706 : ! Compute ice strength as in Hibler (1979)
1707 : !-----------------------------------------------------------------
1708 :
1709 38947307 : strength = Pstar*vice*exp(-Cstar*(c1-aice))
1710 :
1711 : endif ! kstrength
1712 :
1713 : end subroutine icepack_ice_strength
1714 :
1715 : !=======================================================================
1716 : !autodocument_start icepack_step_ridge
1717 : ! Computes sea ice mechanical deformation
1718 : !
1719 : ! authors: William H. Lipscomb, LANL
1720 : ! Elizabeth C. Hunke, LANL
1721 :
1722 1498353240 : subroutine icepack_step_ridge(dt, ndtd, &
1723 1498353240 : hin_max, & ! LCOV_EXCL_LINE
1724 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
1725 1498353240 : aicen, & ! LCOV_EXCL_LINE
1726 1498353240 : trcrn, & ! LCOV_EXCL_LINE
1727 1498353240 : vicen, vsnon, & ! LCOV_EXCL_LINE
1728 1498353240 : aice0, trcr_depend, & ! LCOV_EXCL_LINE
1729 1498353240 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
1730 1498353240 : nt_strata, & ! LCOV_EXCL_LINE
1731 : dardg1dt, dardg2dt, & ! LCOV_EXCL_LINE
1732 : dvirdgdt, opening, & ! LCOV_EXCL_LINE
1733 : fpond, & ! LCOV_EXCL_LINE
1734 : fresh, fhocn, & ! LCOV_EXCL_LINE
1735 2996706480 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
1736 1498353240 : aparticn, krdgn, & ! LCOV_EXCL_LINE
1737 1498353240 : aredistn, vredistn, & ! LCOV_EXCL_LINE
1738 1498353240 : dardg1ndt, dardg2ndt, & ! LCOV_EXCL_LINE
1739 1498353240 : dvirdgndt, & ! LCOV_EXCL_LINE
1740 1498353240 : araftn, vraftn, & ! LCOV_EXCL_LINE
1741 : aice, fsalt, & ! LCOV_EXCL_LINE
1742 1498353240 : first_ice, fzsal, & ! LCOV_EXCL_LINE
1743 1498353240 : flux_bio, closing, & ! LCOV_EXCL_LINE
1744 : Tf, & ! LCOV_EXCL_LINE
1745 : docleanup, dorebin)
1746 :
1747 : real (kind=dbl_kind), intent(in) :: &
1748 : dt ! time step
1749 :
1750 : real (kind=dbl_kind), intent(in) :: &
1751 : Tf ! freezing temperature
1752 :
1753 : integer (kind=int_kind), intent(in) :: &
1754 : ndtd ! number of dynamics supercycles
1755 :
1756 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
1757 : hin_max ! category limits (m)
1758 :
1759 : integer (kind=int_kind), dimension (:), intent(in) :: &
1760 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
1761 : n_trcr_strata ! number of underlying tracer layers
1762 :
1763 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1764 : trcr_base ! = 0 or 1 depending on tracer dependency
1765 : ! argument 2: (1) aice, (2) vice, (3) vsno
1766 :
1767 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
1768 : nt_strata ! indices of underlying tracer layers
1769 :
1770 : real (kind=dbl_kind), intent(inout) :: &
1771 : aice , & ! sea ice concentration ! LCOV_EXCL_LINE
1772 : aice0 , & ! concentration of open water ! LCOV_EXCL_LINE
1773 : rdg_conv , & ! convergence term for ridging (1/s) ! LCOV_EXCL_LINE
1774 : rdg_shear, & ! shear term for ridging (1/s) ! LCOV_EXCL_LINE
1775 : dardg1dt , & ! rate of area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
1776 : dardg2dt , & ! rate of area gain by new ridges (1/s) ! LCOV_EXCL_LINE
1777 : dvirdgdt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
1778 : opening , & ! rate of opening due to divergence/shear (1/s) ! LCOV_EXCL_LINE
1779 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
1780 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1781 : fsalt , & ! salt flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1782 : fhocn ! net heat flux to ocean (W/m^2)
1783 :
1784 : real (kind=dbl_kind), intent(inout), optional :: &
1785 : fzsal ! zsalinity flux to ocean(kg/m^2/s) (deprecated)
1786 :
1787 : real (kind=dbl_kind), intent(inout), optional :: &
1788 : closing ! rate of closing due to divergence/shear (1/s)
1789 :
1790 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1791 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1792 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1793 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
1794 : dardg1ndt, & ! rate of area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
1795 : dardg2ndt, & ! rate of area gain by new ridges (1/s) ! LCOV_EXCL_LINE
1796 : dvirdgndt, & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
1797 : aparticn , & ! participation function ! LCOV_EXCL_LINE
1798 : krdgn , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
1799 : araftn , & ! rafting ice area ! LCOV_EXCL_LINE
1800 : vraftn , & ! rafting ice volume ! LCOV_EXCL_LINE
1801 : aredistn , & ! redistribution function: fraction of new ridge area ! LCOV_EXCL_LINE
1802 : vredistn , & ! redistribution function: fraction of new ridge volume ! LCOV_EXCL_LINE
1803 : faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1804 : flux_bio ! all bio fluxes to ocean
1805 :
1806 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
1807 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
1808 :
1809 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
1810 : trcrn ! tracers
1811 :
1812 : !logical (kind=log_kind), intent(in) :: &
1813 : !tr_pond_topo,& ! if .true., use explicit topography-based ponds ! LCOV_EXCL_LINE
1814 : !tr_aero ,& ! if .true., use aerosol tracers ! LCOV_EXCL_LINE
1815 : !tr_brine !,& ! if .true., brine height differs from ice thickness ! LCOV_EXCL_LINE
1816 :
1817 : logical (kind=log_kind), dimension(:), intent(inout) :: &
1818 : first_ice ! true until ice forms
1819 :
1820 : logical (kind=log_kind), intent(in), optional :: &
1821 : docleanup, & ! if false, do not call cleanup_itd (default true) ! LCOV_EXCL_LINE
1822 : dorebin ! if false, do not call rebin in cleanup_itd (default true)
1823 :
1824 : !autodocument_end
1825 :
1826 : ! local variables
1827 :
1828 : real (kind=dbl_kind) :: &
1829 : dtt ! thermo time step
1830 :
1831 : logical (kind=log_kind) :: &
1832 : ldocleanup, &! if true, call cleanup_itd ! LCOV_EXCL_LINE
1833 : ldorebin ! if true, call rebin in cleanup_itd
1834 :
1835 : logical (kind=log_kind), save :: &
1836 : first_call = .true. ! first call flag
1837 :
1838 : character(len=*),parameter :: subname='(icepack_step_ridge)'
1839 :
1840 : !-----------------------------------------------------------------
1841 : ! Check optional arguments
1842 : !-----------------------------------------------------------------
1843 :
1844 1498353240 : if (icepack_chkoptargflag(first_call)) then
1845 7445 : if (tr_iso) then
1846 282 : if (.not.(present(fiso_ocn))) then
1847 0 : call icepack_warnings_add(subname//' error in fiso_ocn argument, tr_iso=T')
1848 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1849 0 : return
1850 : endif
1851 : endif
1852 : endif
1853 :
1854 1498353240 : if (present(docleanup)) then
1855 0 : ldocleanup = docleanup
1856 : else
1857 1498353240 : ldocleanup = .true.
1858 : endif
1859 :
1860 1498353240 : if (present(dorebin)) then
1861 0 : ldorebin = dorebin
1862 : else
1863 1498353240 : ldorebin = .true.
1864 : endif
1865 :
1866 : !-----------------------------------------------------------------
1867 : ! Identify ice-ocean cells.
1868 : ! Note: We can not limit the loop here using aice>puny because
1869 : ! aice has not yet been updated since the transport (and
1870 : ! it may be out of whack, which the ridging helps fix).-ECH
1871 : !-----------------------------------------------------------------
1872 :
1873 : call ridge_ice (dt, ndtd, &
1874 : hin_max, & ! LCOV_EXCL_LINE
1875 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
1876 : aicen, & ! LCOV_EXCL_LINE
1877 : trcrn, & ! LCOV_EXCL_LINE
1878 : vicen, vsnon, & ! LCOV_EXCL_LINE
1879 : aice0, & ! LCOV_EXCL_LINE
1880 : trcr_depend, & ! LCOV_EXCL_LINE
1881 : trcr_base, & ! LCOV_EXCL_LINE
1882 : n_trcr_strata, & ! LCOV_EXCL_LINE
1883 : nt_strata, & ! LCOV_EXCL_LINE
1884 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
1885 : mu_rdg, tr_brine, & ! LCOV_EXCL_LINE
1886 : dardg1dt, dardg2dt, & ! LCOV_EXCL_LINE
1887 : dvirdgdt, opening, & ! LCOV_EXCL_LINE
1888 : fpond, flux_bio, & ! LCOV_EXCL_LINE
1889 : fresh, fhocn, & ! LCOV_EXCL_LINE
1890 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
1891 : aparticn, krdgn, & ! LCOV_EXCL_LINE
1892 : aredistn, vredistn, & ! LCOV_EXCL_LINE
1893 : dardg1ndt, dardg2ndt, & ! LCOV_EXCL_LINE
1894 : dvirdgndt, Tf, & ! LCOV_EXCL_LINE
1895 : araftn, vraftn, & ! LCOV_EXCL_LINE
1896 1498353240 : closing )
1897 1498353240 : if (icepack_warnings_aborted(subname)) return
1898 :
1899 : !-----------------------------------------------------------------
1900 : ! ITD cleanup: Rebin thickness categories if necessary, and remove
1901 : ! categories with very small areas.
1902 : !-----------------------------------------------------------------
1903 :
1904 1498353240 : if (ldocleanup) then
1905 1498353240 : dtt = dt * ndtd ! for proper averaging over thermo timestep
1906 : call cleanup_itd(dtt, hin_max, &
1907 : aicen, trcrn, & ! LCOV_EXCL_LINE
1908 : vicen, vsnon, & ! LCOV_EXCL_LINE
1909 : aice0, aice, & ! LCOV_EXCL_LINE
1910 : tr_aero, & ! LCOV_EXCL_LINE
1911 : tr_pond_topo, & ! LCOV_EXCL_LINE
1912 : first_ice, & ! LCOV_EXCL_LINE
1913 : trcr_depend, trcr_base, & ! LCOV_EXCL_LINE
1914 : n_trcr_strata, nt_strata, & ! LCOV_EXCL_LINE
1915 : fpond, fresh, & ! LCOV_EXCL_LINE
1916 : fsalt, fhocn, & ! LCOV_EXCL_LINE
1917 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
1918 : flux_bio, Tf, & ! LCOV_EXCL_LINE
1919 1498353240 : dorebin = ldorebin)
1920 1498353240 : if (icepack_warnings_aborted(subname)) return
1921 : endif
1922 :
1923 1498353240 : first_call = .false.
1924 :
1925 : end subroutine icepack_step_ridge
1926 :
1927 : !=======================================================================
1928 :
1929 : end module icepack_mechred
1930 :
1931 : !=======================================================================
|
