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