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 :
41 : use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
42 : use icepack_parameters, only: heat_capacity, conserv_check
43 :
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 :: ridge_ice, &
62 : asum_ridging, & ! LCOV_EXCL_LINE
63 : ridge_itd, & ! LCOV_EXCL_LINE
64 : icepack_ice_strength, & ! LCOV_EXCL_LINE
65 : icepack_step_ridge
66 :
67 : real (kind=dbl_kind), parameter :: &
68 : exp_argmax = 100.0_dbl_kind, & ! maximum argument of exponential for underflow ! LCOV_EXCL_LINE
69 : Cs = p25 , & ! fraction of shear energy contrbtng to ridging ! LCOV_EXCL_LINE
70 : fsnowrdg = p5 , & ! snow fraction that survives in ridging ! LCOV_EXCL_LINE
71 : Gstar = p15 , & ! max value of G(h) that participates ! LCOV_EXCL_LINE
72 : ! (krdg_partic = 0)
73 : astar = p05 , & ! e-folding scale for G(h) participation
74 : !echmod astar = p1 , & ! e-folding scale for G(h) participation ! LCOV_EXCL_LINE
75 : ! (krdg_partic = 1)
76 : maxraft= c1 , & ! max value of hrmin - hi = max thickness
77 : ! of ice that rafts (m)
78 : Hstar = c25 ! determines mean thickness of ridged ice (m)
79 : ! (krdg_redist = 0)
80 : ! Flato & Hibler (1995) have Hstar = 100
81 :
82 : !=======================================================================
83 :
84 : contains
85 :
86 : !=======================================================================
87 :
88 : ! Compute changes in the ice thickness distribution due to divergence
89 : ! and shear.
90 : !
91 : ! author: William H. Lipscomb, LANL
92 :
93 985712928 : subroutine ridge_ice (dt, ndtd, &
94 : ncat, n_aero, & ! LCOV_EXCL_LINE
95 : nilyr, nslyr, & ! LCOV_EXCL_LINE
96 : ntrcr, hin_max, & ! LCOV_EXCL_LINE
97 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
98 : aicen, trcrn, & ! LCOV_EXCL_LINE
99 : vicen, vsnon, & ! LCOV_EXCL_LINE
100 : aice0, & ! LCOV_EXCL_LINE
101 : trcr_depend, trcr_base, & ! LCOV_EXCL_LINE
102 : n_trcr_strata, & ! LCOV_EXCL_LINE
103 : nt_strata, & ! LCOV_EXCL_LINE
104 : krdg_partic, krdg_redist,& ! LCOV_EXCL_LINE
105 : mu_rdg, tr_brine, & ! LCOV_EXCL_LINE
106 : dardg1dt, dardg2dt, & ! LCOV_EXCL_LINE
107 : dvirdgdt, opening, & ! LCOV_EXCL_LINE
108 : fpond, & ! LCOV_EXCL_LINE
109 : fresh, fhocn, & ! LCOV_EXCL_LINE
110 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
111 : aparticn, krdgn, & ! LCOV_EXCL_LINE
112 : aredistn, vredistn, & ! LCOV_EXCL_LINE
113 : dardg1ndt, dardg2ndt, & ! LCOV_EXCL_LINE
114 : dvirdgndt, & ! LCOV_EXCL_LINE
115 : araftn, vraftn, & ! LCOV_EXCL_LINE
116 : closing_flag,closing )
117 :
118 : integer (kind=int_kind), intent(in) :: &
119 : ndtd , & ! number of dynamics subcycles ! LCOV_EXCL_LINE
120 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
121 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
122 : nslyr , & ! number of snow layers ! LCOV_EXCL_LINE
123 : n_aero, & ! number of aerosol tracers ! LCOV_EXCL_LINE
124 : ntrcr ! number of tracers in use
125 :
126 : real (kind=dbl_kind), intent(in) :: &
127 : mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) ! LCOV_EXCL_LINE
128 : dt ! time step
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 : closing_flag, &! flag if closing is valid ! LCOV_EXCL_LINE
165 : tr_brine ! if .true., brine height differs from ice thickness
166 :
167 : ! optional history fields
168 : real (kind=dbl_kind), intent(inout), optional :: &
169 : dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
170 : dardg2dt , & ! rate of fractional area gain by new ridges (1/s) ! LCOV_EXCL_LINE
171 : dvirdgdt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
172 : opening , & ! rate of opening due to divergence/shear (1/s) ! LCOV_EXCL_LINE
173 : closing , & ! rate of closing due to divergence/shear (1/s) ! LCOV_EXCL_LINE
174 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
175 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
176 : fhocn ! net heat flux to ocean (W/m^2)
177 :
178 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
179 : dardg1ndt , & ! rate of fractional area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
180 : dardg2ndt , & ! rate of fractional area gain by new ridges (1/s) ! LCOV_EXCL_LINE
181 : dvirdgndt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
182 : aparticn , & ! participation function ! LCOV_EXCL_LINE
183 : krdgn , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
184 : araftn , & ! rafting ice area ! LCOV_EXCL_LINE
185 : vraftn , & ! rafting ice volume ! LCOV_EXCL_LINE
186 : aredistn , & ! redistribution function: fraction of new ridge area ! LCOV_EXCL_LINE
187 : vredistn ! redistribution function: fraction of new ridge volume
188 :
189 : real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
190 : faero_ocn ! aerosol flux to ocean (kg/m^2/s)
191 :
192 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
193 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
194 :
195 : ! local variables
196 :
197 : real (kind=dbl_kind), dimension (ncat) :: &
198 2356095840 : eicen ! energy of melting for each ice layer (J/m^2)
199 :
200 : real (kind=dbl_kind), dimension (ncat) :: &
201 : esnon, & ! energy of melting for each snow layer (J/m^2) ! LCOV_EXCL_LINE
202 : vbrin, & ! ice volume with defined by brine height (m) ! LCOV_EXCL_LINE
203 2356095840 : sicen ! Bulk salt in h ice (ppt*m)
204 :
205 : real (kind=dbl_kind) :: &
206 : asum , & ! sum of ice and open water area ! LCOV_EXCL_LINE
207 : aksum , & ! ratio of area removed to area ridged ! LCOV_EXCL_LINE
208 : msnow_mlt , & ! mass of snow added to ocean (kg m-2) ! LCOV_EXCL_LINE
209 : esnow_mlt , & ! energy needed to melt snow in ocean (J m-2) ! LCOV_EXCL_LINE
210 : mpond , & ! mass of pond added to ocean (kg m-2) ! LCOV_EXCL_LINE
211 : closing_net, & ! net rate at which area is removed (1/s) ! LCOV_EXCL_LINE
212 : ! (ridging ice area - area of new ridges) / dt
213 98185008 : divu_adv , & ! divu as implied by transport scheme (1/s)
214 : opning , & ! rate of opening due to divergence/shear ! LCOV_EXCL_LINE
215 : ! opning is a local variable;
216 : ! opening is the history diagnostic variable
217 98185008 : ardg1 , & ! fractional area loss by ridging ice
218 : ardg2 , & ! fractional area gain by new ridges ! LCOV_EXCL_LINE
219 : virdg , & ! ice volume ridged ! LCOV_EXCL_LINE
220 98185008 : aopen ! area opening due to divergence/shear
221 :
222 : real (kind=dbl_kind), dimension (n_aero) :: &
223 1876507296 : maero ! aerosol mass added to ocean (kg m-2)
224 :
225 : real (kind=dbl_kind), dimension (n_iso) :: &
226 1873817280 : miso ! isotope mass added to ocean (kg m-2)
227 :
228 : real (kind=dbl_kind), dimension (0:ncat) :: &
229 2650650864 : apartic ! participation function; fraction of ridging
230 : ! and closing associated w/ category n
231 :
232 : real (kind=dbl_kind), dimension (ncat) :: &
233 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
234 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
235 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
236 : krdg , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
237 : ardg1n , & ! area of ice ridged ! LCOV_EXCL_LINE
238 : ardg2n , & ! area of new ridges ! LCOV_EXCL_LINE
239 : virdgn , & ! ridging ice volume ! LCOV_EXCL_LINE
240 1566752928 : mraftn ! rafting ice mask
241 :
242 : real (kind=dbl_kind) :: &
243 : vice_init, vice_final, & ! ice volume summed over categories ! LCOV_EXCL_LINE
244 : vsno_init, vsno_final, & ! snow volume summed over categories ! LCOV_EXCL_LINE
245 : eice_init, eice_final, & ! ice energy summed over layers ! LCOV_EXCL_LINE
246 : vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories ! LCOV_EXCL_LINE
247 : sice_init ,sice_final, & ! ice bulk salinity summed over categories ! LCOV_EXCL_LINE
248 98185008 : esno_init, esno_final ! snow energy summed over layers
249 :
250 : integer (kind=int_kind), parameter :: &
251 : nitermax = 20 ! max number of ridging iterations
252 :
253 : integer (kind=int_kind) :: &
254 : n , & ! thickness category index ! LCOV_EXCL_LINE
255 : niter , & ! iteration counter ! LCOV_EXCL_LINE
256 : k , & ! vertical index ! LCOV_EXCL_LINE
257 : it ! tracer index
258 :
259 : real (kind=dbl_kind) :: &
260 98185008 : dti ! 1 / dt
261 :
262 : logical (kind=log_kind) :: &
263 : iterate_ridging ! if true, repeat the ridging
264 :
265 : character (len=char_len) :: &
266 : fieldid ! field identifier
267 :
268 : character(len=*),parameter :: subname='(ridge_ice)'
269 :
270 : !-----------------------------------------------------------------
271 : ! Initialize
272 : !-----------------------------------------------------------------
273 :
274 985712928 : msnow_mlt = c0
275 985712928 : esnow_mlt = c0
276 1029329616 : maero (:) = c0
277 1046814720 : miso (:) = c0
278 985712928 : mpond = c0
279 985712928 : ardg1 = c0
280 985712928 : ardg2 = c0
281 985712928 : virdg = c0
282 5886032400 : ardg1n(:) = c0
283 5886032400 : ardg2n(:) = c0
284 5886032400 : virdgn(:) = c0
285 5886032400 : mraftn(:) = c0
286 985712928 : aopen = c0
287 :
288 : !-----------------------------------------------------------------
289 : ! Compute area of ice plus open water before ridging.
290 : !-----------------------------------------------------------------
291 :
292 985712928 : call asum_ridging (ncat, aicen, aice0, asum)
293 985712928 : if (icepack_warnings_aborted(subname)) return
294 :
295 : !-----------------------------------------------------------------
296 : ! Compute the area opening and closing.
297 : !-----------------------------------------------------------------
298 :
299 985712928 : if (closing_flag) then
300 :
301 0 : opning = opening
302 0 : closing_net = closing
303 0 : divu_adv = opening - closing
304 :
305 : else
306 :
307 : call ridge_prep (dt, &
308 : ncat, hin_max, & ! LCOV_EXCL_LINE
309 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
310 : asum, closing_net, & ! LCOV_EXCL_LINE
311 985712928 : divu_adv, opning)
312 :
313 : endif
314 :
315 985712928 : if (icepack_warnings_aborted(subname)) return
316 :
317 : !-----------------------------------------------------------------
318 : ! Compute initial values of conserved quantities.
319 : !-----------------------------------------------------------------
320 :
321 985712928 : if (conserv_check) then
322 :
323 162745968 : do n = 1, ncat
324 139496544 : eicen(n) = c0
325 139496544 : esnon(n) = c0
326 139496544 : sicen(n) = c0
327 139496544 : vbrin(n) = c0
328 :
329 1115972352 : do k = 1, nilyr
330 258241536 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
331 1105596576 : * vicen(n)/real(nilyr,kind=dbl_kind)
332 258241536 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
333 1245093120 : * vicen(n)/real(nilyr,kind=dbl_kind)
334 : enddo
335 278993088 : do k = 1, nslyr
336 36891648 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
337 297438912 : * vsnon(n)/real(nslyr,kind=dbl_kind)
338 : enddo
339 139496544 : vbrin(n) = vicen(n)
340 162745968 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vicen(n)
341 : enddo ! n
342 :
343 : call column_sum (ncat, &
344 23249424 : vicen, vice_init)
345 23249424 : if (icepack_warnings_aborted(subname)) return
346 : call column_sum (ncat, &
347 23249424 : vsnon, vsno_init)
348 23249424 : if (icepack_warnings_aborted(subname)) return
349 : call column_sum (ncat, &
350 23249424 : eicen, eice_init)
351 23249424 : if (icepack_warnings_aborted(subname)) return
352 : call column_sum (ncat, &
353 23249424 : esnon, esno_init)
354 23249424 : if (icepack_warnings_aborted(subname)) return
355 : call column_sum (ncat, &
356 23249424 : sicen, sice_init)
357 23249424 : if (icepack_warnings_aborted(subname)) return
358 : call column_sum (ncat, &
359 23249424 : vbrin, vbri_init)
360 23249424 : if (icepack_warnings_aborted(subname)) return
361 :
362 : endif
363 :
364 985712945 : rdg_iteration: do niter = 1, nitermax
365 :
366 : !-----------------------------------------------------------------
367 : ! Compute the thickness distribution of ridging ice
368 : ! and various quantities associated with the new ridged ice.
369 : !-----------------------------------------------------------------
370 :
371 : call ridge_itd (ncat, aice0, &
372 : aicen, vicen, & ! LCOV_EXCL_LINE
373 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
374 : mu_rdg, & ! LCOV_EXCL_LINE
375 : aksum, apartic, & ! LCOV_EXCL_LINE
376 : hrmin, hrmax, & ! LCOV_EXCL_LINE
377 : hrexp, krdg, & ! LCOV_EXCL_LINE
378 : aparticn, krdgn, & ! LCOV_EXCL_LINE
379 985712945 : mraftn)
380 985712945 : if (icepack_warnings_aborted(subname)) return
381 :
382 : !-----------------------------------------------------------------
383 : ! Redistribute area, volume, and energy.
384 : !-----------------------------------------------------------------
385 :
386 : call ridge_shift (ntrcr, dt, &
387 : ncat, hin_max, & ! LCOV_EXCL_LINE
388 : aicen, trcrn, & ! LCOV_EXCL_LINE
389 : vicen, vsnon, & ! LCOV_EXCL_LINE
390 : aice0, trcr_depend, & ! LCOV_EXCL_LINE
391 : trcr_base, n_trcr_strata,& ! LCOV_EXCL_LINE
392 : nt_strata, krdg_redist, & ! LCOV_EXCL_LINE
393 : aksum, apartic, & ! LCOV_EXCL_LINE
394 : hrmin, hrmax, & ! LCOV_EXCL_LINE
395 : hrexp, krdg, & ! LCOV_EXCL_LINE
396 : closing_net, opning, & ! LCOV_EXCL_LINE
397 : ardg1, ardg2, & ! LCOV_EXCL_LINE
398 : virdg, aopen, & ! LCOV_EXCL_LINE
399 : ardg1n, ardg2n, & ! LCOV_EXCL_LINE
400 : virdgn, & ! LCOV_EXCL_LINE
401 : nslyr, n_aero, & ! LCOV_EXCL_LINE
402 : msnow_mlt, esnow_mlt, & ! LCOV_EXCL_LINE
403 : maero, miso, & ! LCOV_EXCL_LINE
404 : mpond, & ! LCOV_EXCL_LINE
405 985712945 : aredistn, vredistn)
406 985712945 : if (icepack_warnings_aborted(subname)) return
407 :
408 : !-----------------------------------------------------------------
409 : ! Make sure the new area = 1. If not (because the closing
410 : ! and opening rates were reduced above), prepare to ridge again
411 : ! with new rates.
412 : !-----------------------------------------------------------------
413 :
414 985712945 : call asum_ridging (ncat, aicen, aice0, asum)
415 985712945 : if (icepack_warnings_aborted(subname)) return
416 :
417 985712945 : if (abs(asum - c1) < puny) then
418 985712928 : iterate_ridging = .false.
419 985712928 : closing_net = c0
420 985712928 : opning = c0
421 : else
422 17 : iterate_ridging = .true.
423 17 : divu_adv = (c1 - asum) / dt
424 17 : closing_net = max(c0, -divu_adv)
425 17 : opning = max(c0, divu_adv)
426 : endif
427 :
428 : !-----------------------------------------------------------------
429 : ! If done, exit. If not, prepare to ridge again.
430 : !-----------------------------------------------------------------
431 :
432 985712945 : if (iterate_ridging) then
433 17 : write(warnstr,*) subname, 'Repeat ridging, niter =', niter
434 17 : call icepack_warnings_add(warnstr)
435 : else
436 985712928 : exit rdg_iteration
437 : endif
438 :
439 17 : if (niter == nitermax) then
440 0 : write(warnstr,*) ' '
441 0 : call icepack_warnings_add(warnstr)
442 0 : write(warnstr,*) subname, 'Exceeded max number of ridging iterations'
443 0 : call icepack_warnings_add(warnstr)
444 0 : write(warnstr,*) subname, 'max =',nitermax
445 0 : call icepack_warnings_add(warnstr)
446 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
447 0 : call icepack_warnings_add(subname//" ridge_ice: Exceeded max number of ridging iterations" )
448 0 : return
449 : endif
450 :
451 : enddo rdg_iteration ! niter
452 :
453 : !-----------------------------------------------------------------
454 : ! Compute final values of conserved quantities.
455 : ! Check for conservation (allowing for snow thrown into ocean).
456 : !-----------------------------------------------------------------
457 :
458 985712928 : if (conserv_check) then
459 :
460 162745968 : do n = 1, ncat
461 139496544 : eicen(n) = c0
462 139496544 : esnon(n) = c0
463 139496544 : sicen(n) = c0
464 139496544 : vbrin(n) = c0
465 :
466 1115972352 : do k = 1, nilyr
467 258241536 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
468 1105596576 : * vicen(n)/real(nilyr,kind=dbl_kind)
469 258241536 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
470 1245093120 : * vicen(n)/real(nilyr,kind=dbl_kind)
471 : enddo
472 278993088 : do k = 1, nslyr
473 36891648 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
474 297438912 : * vsnon(n)/real(nslyr,kind=dbl_kind)
475 : enddo
476 139496544 : vbrin(n) = vicen(n)
477 162745968 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vbrin(n)
478 : enddo ! n
479 :
480 : call column_sum (ncat, &
481 23249424 : vicen, vice_final)
482 23249424 : if (icepack_warnings_aborted(subname)) return
483 : call column_sum (ncat, &
484 23249424 : vsnon, vsno_final)
485 23249424 : if (icepack_warnings_aborted(subname)) return
486 : call column_sum (ncat, &
487 23249424 : eicen, eice_final)
488 23249424 : if (icepack_warnings_aborted(subname)) return
489 : call column_sum (ncat, &
490 23249424 : esnon, esno_final)
491 23249424 : if (icepack_warnings_aborted(subname)) return
492 : call column_sum (ncat, &
493 23249424 : sicen, sice_final)
494 23249424 : if (icepack_warnings_aborted(subname)) return
495 : call column_sum (ncat, &
496 23249424 : vbrin, vbri_final)
497 23249424 : if (icepack_warnings_aborted(subname)) return
498 :
499 23249424 : vsno_final = vsno_final + msnow_mlt/rhos
500 23249424 : esno_final = esno_final + esnow_mlt
501 :
502 23249424 : fieldid = subname//':vice'
503 : call column_conservation_check (fieldid, &
504 : vice_init, vice_final, & ! LCOV_EXCL_LINE
505 23249424 : puny)
506 23249424 : if (icepack_warnings_aborted(subname)) return
507 23249424 : fieldid = subname//':vsno'
508 : call column_conservation_check (fieldid, &
509 : vsno_init, vsno_final, & ! LCOV_EXCL_LINE
510 23249424 : puny)
511 23249424 : if (icepack_warnings_aborted(subname)) return
512 23249424 : fieldid = subname//':eice'
513 : call column_conservation_check (fieldid, &
514 : eice_init, eice_final, & ! LCOV_EXCL_LINE
515 23249424 : puny*Lfresh*rhoi)
516 23249424 : if (icepack_warnings_aborted(subname)) return
517 23249424 : fieldid = subname//':esno'
518 : call column_conservation_check (fieldid, &
519 : esno_init, esno_final, & ! LCOV_EXCL_LINE
520 23249424 : puny*Lfresh*rhos)
521 23249424 : if (icepack_warnings_aborted(subname)) return
522 23249424 : fieldid = subname//':sice'
523 : call column_conservation_check (fieldid, &
524 : sice_init, sice_final, & ! LCOV_EXCL_LINE
525 23249424 : puny)
526 23249424 : if (icepack_warnings_aborted(subname)) return
527 23249424 : fieldid = subname//':vbrin'
528 : call column_conservation_check (fieldid, &
529 : vbri_init, vbri_final, & ! LCOV_EXCL_LINE
530 23249424 : puny*c10)
531 23249424 : if (icepack_warnings_aborted(subname)) return
532 :
533 : endif ! conserv_check
534 :
535 : !-----------------------------------------------------------------
536 : ! Compute ridging diagnostics.
537 : !-----------------------------------------------------------------
538 :
539 985712928 : dti = c1/dt
540 :
541 985712928 : if (present(dardg1dt)) then
542 985712928 : dardg1dt = ardg1*dti
543 : endif
544 985712928 : if (present(dardg2dt)) then
545 985712928 : dardg2dt = ardg2*dti
546 : endif
547 985712928 : if (present(dvirdgdt)) then
548 985712928 : dvirdgdt = virdg*dti
549 : endif
550 985712928 : if (present(opening)) then
551 985712928 : opening = aopen*dti
552 : endif
553 985712928 : if (present(dardg1ndt)) then
554 5886032400 : do n = 1, ncat
555 5886032400 : dardg1ndt(n) = ardg1n(n)*dti
556 : enddo
557 : endif
558 985712928 : if (present(dardg2ndt)) then
559 5886032400 : do n = 1, ncat
560 5886032400 : dardg2ndt(n) = ardg2n(n)*dti
561 : enddo
562 : endif
563 985712928 : if (present(dvirdgndt)) then
564 5886032400 : do n = 1, ncat
565 5886032400 : dvirdgndt(n) = virdgn(n)*dti
566 : enddo
567 : endif
568 985712928 : if (present(araftn)) then
569 5886032400 : do n = 1, ncat
570 5886032400 : araftn(n) = mraftn(n)*ardg2n(n)
571 : ! araftn(n) = mraftn(n)*ardg1n(n)*p5
572 : enddo
573 : endif
574 985712928 : if (present(vraftn)) then
575 5886032400 : do n = 1, ncat
576 5886032400 : vraftn(n) = mraftn(n)*virdgn(n)
577 : enddo
578 : endif
579 :
580 : !-----------------------------------------------------------------
581 : ! Update fresh water and heat fluxes due to snow melt.
582 : !-----------------------------------------------------------------
583 :
584 : ! use thermodynamic time step (ndtd*dt here) to average properly
585 985712928 : dti = c1/(ndtd*dt)
586 :
587 985712928 : if (present(fresh)) then
588 985712928 : fresh = fresh + msnow_mlt*dti
589 : endif
590 985712928 : if (present(fhocn)) then
591 985712928 : fhocn = fhocn + esnow_mlt*dti
592 : endif
593 985712928 : if (present(faero_ocn)) then
594 1029329616 : do it = 1, n_aero
595 1029329616 : faero_ocn(it) = faero_ocn(it) + maero(it)*dti
596 : enddo
597 : endif
598 985712928 : if (tr_iso) then
599 : ! check size fiso_ocn vs n_iso ???
600 81469056 : do it = 1, n_iso
601 81469056 : fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti
602 : enddo
603 : endif
604 985712928 : if (present(fpond)) then
605 985712928 : fpond = fpond - mpond ! units change later
606 : endif
607 :
608 : !-----------------------------------------------------------------
609 : ! Check for fractional ice area > 1.
610 : !-----------------------------------------------------------------
611 :
612 985712928 : if (abs(asum - c1) > puny) then
613 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
614 0 : call icepack_warnings_add(subname//" total area > 1" )
615 :
616 0 : write(warnstr,*) ' '
617 0 : call icepack_warnings_add(warnstr)
618 0 : write(warnstr,*) subname, 'Ridging error: total area > 1'
619 0 : call icepack_warnings_add(warnstr)
620 0 : write(warnstr,*) subname, 'area:', asum
621 0 : call icepack_warnings_add(warnstr)
622 0 : write(warnstr,*) subname, 'n, aicen:'
623 0 : call icepack_warnings_add(warnstr)
624 0 : write(warnstr,*) subname, 0, aice0
625 0 : call icepack_warnings_add(warnstr)
626 0 : do n = 1, ncat
627 0 : write(warnstr,*) subname, n, aicen(n)
628 0 : call icepack_warnings_add(warnstr)
629 : enddo
630 0 : return
631 : endif
632 :
633 : end subroutine ridge_ice
634 :
635 : !=======================================================================
636 :
637 : ! Find the total area of ice plus open water in each grid cell.
638 : !
639 : ! This is similar to the aggregate_area subroutine except that the
640 : ! total area can be greater than 1, so the open water area is
641 : ! included in the sum instead of being computed as a residual.
642 : !
643 : ! author: William H. Lipscomb, LANL
644 :
645 2209056193 : subroutine asum_ridging (ncat, aicen, aice0, asum)
646 :
647 : integer (kind=int_kind), intent(in) :: &
648 : ncat ! number of thickness categories
649 :
650 : real (kind=dbl_kind), dimension (:), intent(in) :: &
651 : aicen ! concentration of ice in each category
652 :
653 : real (kind=dbl_kind), intent(in) :: &
654 : aice0 ! concentration of open water
655 :
656 : real (kind=dbl_kind), intent(out):: &
657 : asum ! sum of ice and open water area
658 :
659 : ! local variables
660 :
661 : integer (kind=int_kind) :: n
662 :
663 : character(len=*),parameter :: subname='(asum_ridging)'
664 :
665 2209056193 : asum = aice0
666 13216591264 : do n = 1, ncat
667 13216591264 : asum = asum + aicen(n)
668 : enddo
669 :
670 2209056193 : end subroutine asum_ridging
671 :
672 : !=======================================================================
673 :
674 : ! Initialize arrays, compute area of closing and opening
675 : !
676 : ! author: William H. Lipscomb, LANL
677 :
678 985712928 : subroutine ridge_prep (dt, &
679 : ncat, hin_max, & ! LCOV_EXCL_LINE
680 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
681 : asum, closing_net, & ! LCOV_EXCL_LINE
682 : divu_adv, opning)
683 :
684 : integer (kind=int_kind), intent(in) :: &
685 : ncat ! number of thickness categories
686 :
687 : real (kind=dbl_kind), intent(in) :: &
688 : dt ! time step (s)
689 :
690 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
691 : hin_max ! category limits (m)
692 :
693 : real (kind=dbl_kind), intent(in) :: &
694 : rdg_conv , & ! normalized energy dissipation due to convergence (1/s) ! LCOV_EXCL_LINE
695 : rdg_shear ! normalized energy dissipation due to shear (1/s)
696 :
697 : real (kind=dbl_kind), intent(inout):: &
698 : asum ! sum of ice and open water area
699 :
700 : real (kind=dbl_kind), &
701 : intent(out):: & ! LCOV_EXCL_LINE
702 : closing_net, & ! net rate at which area is removed (1/s) ! LCOV_EXCL_LINE
703 : divu_adv , & ! divu as implied by transport scheme (1/s) ! LCOV_EXCL_LINE
704 : opning ! rate of opening due to divergence/shear
705 :
706 : ! local variables
707 :
708 : real (kind=dbl_kind), parameter :: &
709 : big = 1.0e+8_dbl_kind
710 :
711 : character(len=*),parameter :: subname='(ridge_prep)'
712 :
713 : ! Set hin_max(ncat) to a big value to ensure that all ridged ice
714 : ! is thinner than hin_max(ncat).
715 985712928 : hin_max(ncat) = big
716 :
717 : !-----------------------------------------------------------------
718 : ! Compute the net rate of closing due to convergence
719 : ! and shear, based on Flato and Hibler (1995).
720 : !
721 : ! For the elliptical yield curve:
722 : ! rdg_conv = -min (divu, 0)
723 : ! rdg_shear = (1/2) * (Delta - abs(divu))
724 : ! Note that the shear term also accounts for divergence.
725 : !
726 : ! The energy dissipation rate is equal to the net closing rate
727 : ! times the ice strength.
728 : !
729 : ! NOTE: The NET closing rate is equal to the rate that open water
730 : ! area is removed, plus the rate at which ice area is removed by
731 : ! ridging, minus the rate at which area is added in new ridges.
732 : ! The GROSS closing rate is equal to the first two terms (open
733 : ! water closing and thin ice ridging) without the third term
734 : ! (thick, newly ridged ice).
735 : !
736 : ! rdg_conv is calculated differently in EAP (update_ice_rdg) and
737 : ! represents closing_net directly. In that case, rdg_shear=0.
738 : !-----------------------------------------------------------------
739 :
740 985712928 : closing_net = Cs*rdg_shear + rdg_conv
741 :
742 : !-----------------------------------------------------------------
743 : ! Compute divu_adv, the divergence rate given by the transport/
744 : ! advection scheme, which may not be equal to divu as computed
745 : ! from the velocity field.
746 : !
747 : ! If divu_adv < 0, make sure the closing rate is large enough
748 : ! to give asum = 1.0 after ridging.
749 : !-----------------------------------------------------------------
750 :
751 985712928 : divu_adv = (c1-asum) / dt
752 :
753 985712928 : if (divu_adv < c0) closing_net = max(closing_net, -divu_adv)
754 :
755 : !-----------------------------------------------------------------
756 : ! Compute the (non-negative) opening rate that will give
757 : ! asum = 1.0 after ridging.
758 : !-----------------------------------------------------------------
759 :
760 985712928 : opning = closing_net + divu_adv
761 :
762 985712928 : end subroutine ridge_prep
763 :
764 : !=======================================================================
765 :
766 : ! Compute the thickness distribution of the ice and open water
767 : ! participating in ridging and of the resulting ridges.
768 : !
769 : ! This version includes new options for ridging participation and
770 : ! redistribution.
771 : ! The new participation scheme (krdg_partic = 1) improves stability
772 : ! by increasing the time scale for large changes in ice strength.
773 : ! The new exponential redistribution function (krdg_redist = 1) improves
774 : ! agreement between ITDs of modeled and observed ridges.
775 : !
776 : ! author: William H. Lipscomb, LANL
777 : !
778 : ! 2006: Changed subroutine name to ridge_itd
779 : ! Added new options for ridging participation and redistribution.
780 :
781 1223343265 : subroutine ridge_itd (ncat, aice0, &
782 : aicen, vicen, & ! LCOV_EXCL_LINE
783 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
784 : mu_rdg, & ! LCOV_EXCL_LINE
785 : aksum, apartic, & ! LCOV_EXCL_LINE
786 : hrmin, hrmax, & ! LCOV_EXCL_LINE
787 : hrexp, krdg, & ! LCOV_EXCL_LINE
788 : aparticn, krdgn, & ! LCOV_EXCL_LINE
789 1223343265 : mraft)
790 :
791 : integer (kind=int_kind), intent(in) :: &
792 : ncat ! number of thickness categories
793 :
794 : real (kind=dbl_kind), intent(in) :: &
795 : mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) ! LCOV_EXCL_LINE
796 : aice0 ! concentration of open water
797 :
798 : real (kind=dbl_kind), dimension (:), intent(in) :: &
799 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
800 : vicen ! volume per unit area of ice (m)
801 :
802 : integer (kind=int_kind), intent(in) :: &
803 : krdg_partic , & ! selects participation function ! LCOV_EXCL_LINE
804 : krdg_redist ! selects redistribution function
805 :
806 : real (kind=dbl_kind), intent(out):: &
807 : aksum ! ratio of area removed to area ridged
808 :
809 : real (kind=dbl_kind), dimension (0:ncat), intent(out) :: &
810 : apartic ! participation function; fraction of ridging
811 : ! and closing associated w/ category n
812 :
813 : real (kind=dbl_kind), dimension (:), intent(out) :: &
814 : hrmin , & ! minimum ridge thickness ! LCOV_EXCL_LINE
815 : hrmax , & ! maximum ridge thickness (krdg_redist = 0) ! LCOV_EXCL_LINE
816 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1) ! LCOV_EXCL_LINE
817 : krdg ! mean ridge thickness/thickness of ridging ice
818 :
819 : ! diagnostic, category values
820 : real (kind=dbl_kind), dimension(:), intent(out), optional :: &
821 : aparticn, & ! participation function ! LCOV_EXCL_LINE
822 : krdgn ! mean ridge thickness/thickness of ridging ice
823 :
824 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
825 : mraft ! rafting ice mask
826 :
827 : ! local variables
828 :
829 : integer (kind=int_kind) :: &
830 : n ! thickness category index
831 :
832 : real (kind=dbl_kind), parameter :: &
833 : Gstari = c1/Gstar, & ! LCOV_EXCL_LINE
834 : astari = c1/astar
835 :
836 : real (kind=dbl_kind), dimension(-1:ncat) :: &
837 3169108528 : Gsum ! Gsum(n) = sum of areas in categories 0 to n
838 :
839 : real (kind=dbl_kind) :: &
840 121592200 : work ! temporary work variable
841 :
842 : real (kind=dbl_kind) :: &
843 : hi , & ! ice thickness for each cat (m) ! LCOV_EXCL_LINE
844 : hrmean , & ! mean ridge thickness (m) ! LCOV_EXCL_LINE
845 121592200 : xtmp ! temporary variable
846 :
847 : character(len=*),parameter :: subname='(ridge_itd)'
848 :
849 : !-----------------------------------------------------------------
850 : ! Initialize
851 : !-----------------------------------------------------------------
852 :
853 9777245394 : Gsum (-1:ncat) = c0 ! initialize
854 :
855 1223343265 : Gsum (-1) = c0 ! by definition
856 1223343265 : if (aice0 > puny) then
857 1021292875 : Gsum(0) = aice0
858 : else
859 202050390 : Gsum(0) = Gsum(-1)
860 : endif
861 1223343265 : apartic(0) = c0
862 :
863 7330558864 : do n = 1, ncat
864 6107215599 : apartic(n) = c0
865 6107215599 : hrmin (n) = c0
866 6107215599 : hrmax (n) = c0
867 6107215599 : hrexp (n) = c0
868 6107215599 : krdg (n) = c1
869 :
870 : !-----------------------------------------------------------------
871 : ! Compute the thickness distribution of ice participating in ridging.
872 : !-----------------------------------------------------------------
873 :
874 : !-----------------------------------------------------------------
875 : ! First compute the cumulative thickness distribution function Gsum,
876 : ! where Gsum(n) is the fractional area in categories 0 to n.
877 : ! Ignore categories with very small areas.
878 : !-----------------------------------------------------------------
879 :
880 7330558864 : if (aicen(n) > puny) then
881 2359956010 : Gsum(n) = Gsum(n-1) + aicen(n)
882 : else
883 3747259589 : Gsum(n) = Gsum(n-1)
884 : endif
885 : enddo
886 :
887 : ! normalize
888 :
889 1223343265 : if (Gsum(ncat) > c0) then
890 1223343265 : work = c1 / Gsum(ncat)
891 8553902129 : do n = 0, ncat
892 8553902129 : Gsum(n) = Gsum(n) * work
893 : enddo
894 : endif
895 :
896 : !-----------------------------------------------------------------
897 : ! Compute the participation function apartic; this is analogous to
898 : ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
899 : !
900 : ! area lost from category n due to ridging/closing
901 : ! apartic(n) = --------------------------------------------------
902 : ! total area lost due to ridging/closing
903 : !
904 : !-----------------------------------------------------------------
905 :
906 1223343265 : if (krdg_partic == 0) then ! Thornike et al. 1975 formulation
907 :
908 : !-----------------------------------------------------------------
909 : ! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
910 : ! The expressions for apartic are found by integrating b(h)g(h) between
911 : ! the category boundaries.
912 : !-----------------------------------------------------------------
913 :
914 339734496 : do n = 0, ncat
915 339734496 : if (Gsum(n) < Gstar) then
916 2772442 : apartic(n) = Gstari*(Gsum(n ) - Gsum(n-1)) * &
917 34328287 : (c2 - Gstari*(Gsum(n-1) + Gsum(n )))
918 245083342 : elseif (Gsum(n-1) < Gstar) then
919 8434176 : apartic(n) = Gstari*(Gstar - Gsum(n-1)) * &
920 65926176 : (c2 - Gstari*(Gstar + Gsum(n-1)))
921 : endif
922 : enddo ! n
923 :
924 1161634177 : elseif (krdg_partic==1) then ! exponential dependence on G(h)
925 :
926 : !-----------------------------------------------------------------
927 : ! b(h) = exp(-G(h)/astar)
928 : ! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)].
929 : ! The expression for apartic is found by integrating b(h)g(h)
930 : ! between the category boundaries.
931 : !-----------------------------------------------------------------
932 :
933 : ! precompute exponential terms using Gsum as work array
934 1161634177 : xtmp = c1 / (c1 - exp(-astari))
935 1161634177 : Gsum(-1) = exp(-Gsum(-1)*astari) * xtmp
936 8214167633 : do n = 0, ncat
937 7052533456 : Gsum(n) = exp(-Gsum(n)*astari) * xtmp
938 8214167633 : apartic(n) = Gsum(n-1) - Gsum(n)
939 : enddo ! n
940 :
941 : endif ! krdg_partic
942 :
943 : !-----------------------------------------------------------------
944 : ! Compute variables related to ITD of ridged ice:
945 : !
946 : ! krdg = mean ridge thickness / thickness of ridging ice
947 : ! hrmin = min ridge thickness
948 : ! hrmax = max ridge thickness (krdg_redist = 0)
949 : ! hrexp = ridge e-folding scale (krdg_redist = 1)
950 : !----------------------------------------------------------------
951 :
952 1223343265 : if (krdg_redist == 0) then ! Hibler 1980 formulation
953 :
954 : !-----------------------------------------------------------------
955 : ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
956 : !
957 : ! This parameterization is a modified version of Hibler (1980).
958 : ! In the original paper the min ridging thickness is hrmin = 2*hi,
959 : ! and the max thickness is hrmax = 2*sqrt(hi*Hstar).
960 : !
961 : ! Here the min thickness is hrmin = min(2*hi, hi+maxraft),
962 : ! so thick ridging ice is not required to raft.
963 : !
964 : !-----------------------------------------------------------------
965 :
966 278025408 : do n = 1, ncat
967 278025408 : if (aicen(n) > puny) then
968 107568643 : hi = vicen(n) / aicen(n)
969 107568643 : hrmin(n) = min(c2*hi, hi + maxraft)
970 107568643 : hrmax(n) = c2*sqrt(Hstar*hi)
971 107568643 : hrmax(n) = max(hrmax(n), hrmin(n)+puny)
972 107568643 : hrmean = p5 * (hrmin(n) + hrmax(n))
973 107568643 : krdg(n) = hrmean / hi
974 :
975 : ! diagnostic rafting mask not implemented
976 : endif
977 : enddo ! n
978 :
979 : else ! krdg_redist = 1; exponential redistribution
980 :
981 : !-----------------------------------------------------------------
982 : ! The ridge ITD is a negative exponential:
983 : !
984 : ! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin
985 : !
986 : ! where hrmin is the minimum thickness of ridging ice and
987 : ! hrexp is the e-folding thickness.
988 : !
989 : ! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
990 : ! That is, the minimum ridge thickness results from rafting,
991 : ! unless the ice is thicker than maxraft.
992 : !
993 : ! Also, assume that hrexp = mu_rdg*sqrt(hi).
994 : ! The parameter mu_rdg is tuned to give e-folding scales mostly
995 : ! in the range 2-4 m as observed by upward-looking sonar.
996 : !
997 : ! Values of mu_rdg in the right column give ice strengths
998 : ! roughly equal to values of Hstar in the left column
999 : ! (within ~10 kN/m for typical ITDs):
1000 : !
1001 : ! Hstar mu_rdg
1002 : !
1003 : ! 25 3.0
1004 : ! 50 4.0
1005 : ! 75 5.0
1006 : ! 100 6.0
1007 : !-----------------------------------------------------------------
1008 :
1009 7052533456 : do n = 1, ncat
1010 7052533456 : if (aicen(n) > puny) then
1011 2252387367 : hi = vicen(n) / aicen(n)
1012 2252387367 : hi = max(hi,puny)
1013 2252387367 : hrmin(n) = min(c2*hi, hi + maxraft)
1014 2252387367 : hrexp(n) = mu_rdg * sqrt(hi)
1015 2252387367 : krdg(n) = (hrmin(n) + hrexp(n)) / hi
1016 :
1017 : !echmod: check computational efficiency
1018 : ! diagnostic rafting mask
1019 2252387367 : if (present(mraft)) then
1020 1072347513 : mraft(n) = max(c0, sign(c1, hi+maxraft-hrmin(n)))
1021 1072347513 : xtmp = mraft(n)*((c2*hi+hrexp(n))/hi - krdg(n))
1022 1072347513 : mraft(n) = max(c0, sign(c1, puny-abs(xtmp)))
1023 : endif
1024 : endif
1025 : enddo
1026 :
1027 : endif ! krdg_redist
1028 :
1029 : !----------------------------------------------------------------
1030 : ! Compute aksum = net ice area removed / total area participating.
1031 : ! For instance, if a unit area of ice with h = 1 participates in
1032 : ! ridging to form a ridge with a = 1/3 and h = 3, then
1033 : ! aksum = 1 - 1/3 = 2/3.
1034 : !----------------------------------------------------------------
1035 :
1036 1223343265 : aksum = apartic(0) ! area participating = area removed
1037 :
1038 7330558864 : do n = 1, ncat
1039 : ! area participating > area removed
1040 7330558864 : aksum = aksum + apartic(n) * (c1 - c1/krdg(n))
1041 : enddo
1042 :
1043 : ! diagnostics
1044 1223343265 : if (present(aparticn)) then
1045 5886032502 : do n = 1, ncat
1046 5886032502 : aparticn(n) = apartic(n)
1047 : enddo
1048 : endif
1049 1223343265 : if (present(krdgn)) then
1050 5886032502 : do n = 1, ncat
1051 5886032502 : krdgn(n) = krdg(n)
1052 : enddo
1053 : endif
1054 :
1055 1223343265 : end subroutine ridge_itd
1056 :
1057 : !=======================================================================
1058 :
1059 : ! Remove area, volume, and energy from each ridging category
1060 : ! and add to thicker ice categories.
1061 : !
1062 : ! Tracers: Ridging conserves ice volume and therefore conserves volume
1063 : ! tracers. It does not conserve ice area, and therefore a portion of area
1064 : ! tracers are lost (corresponding to the net closing). Area tracers on
1065 : ! ice that participates in ridging are carried onto the resulting ridged
1066 : ! ice (except the portion that are lost due to closing). Therefore,
1067 : ! tracers must be decremented if they are lost to the ocean during ridging
1068 : ! (e.g. snow, ponds) or if they are being carried only on the level ice
1069 : ! area.
1070 : !
1071 : ! author: William H. Lipscomb, LANL
1072 :
1073 985712945 : subroutine ridge_shift (ntrcr, dt, &
1074 : ncat, hin_max, & ! LCOV_EXCL_LINE
1075 : aicen, trcrn, & ! LCOV_EXCL_LINE
1076 : vicen, vsnon, & ! LCOV_EXCL_LINE
1077 : aice0, trcr_depend, & ! LCOV_EXCL_LINE
1078 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
1079 : nt_strata, krdg_redist, & ! LCOV_EXCL_LINE
1080 : aksum, apartic, & ! LCOV_EXCL_LINE
1081 : hrmin, hrmax, & ! LCOV_EXCL_LINE
1082 : hrexp, krdg, & ! LCOV_EXCL_LINE
1083 : closing_net, opning, & ! LCOV_EXCL_LINE
1084 : ardg1, ardg2, & ! LCOV_EXCL_LINE
1085 : virdg, aopen, & ! LCOV_EXCL_LINE
1086 : ardg1nn, ardg2nn, & ! LCOV_EXCL_LINE
1087 : virdgnn, & ! LCOV_EXCL_LINE
1088 : nslyr, n_aero, & ! LCOV_EXCL_LINE
1089 : msnow_mlt, esnow_mlt, & ! LCOV_EXCL_LINE
1090 : maero, miso, & ! LCOV_EXCL_LINE
1091 : mpond, & ! LCOV_EXCL_LINE
1092 985712945 : aredistn, vredistn)
1093 :
1094 : integer (kind=int_kind), intent(in) :: &
1095 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
1096 : nslyr , & ! number of snow layers ! LCOV_EXCL_LINE
1097 : ntrcr , & ! number of tracers in use ! LCOV_EXCL_LINE
1098 : n_aero, & ! number of aerosol tracers ! LCOV_EXCL_LINE
1099 : krdg_redist ! selects redistribution function
1100 :
1101 : real (kind=dbl_kind), intent(in) :: &
1102 : dt ! time step (s)
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 2356095878 : vsnon_init ! snow volume before ridging
1183 :
1184 : real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
1185 12711988415 : atrcrn ! aicen*trcrn
1186 :
1187 : real (kind=dbl_kind), dimension(3) :: &
1188 392740036 : 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 98185009 : 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 98185009 : fvol ! fraction of new ridge volume going to nr
1208 :
1209 : real (kind=dbl_kind) :: &
1210 98185009 : 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 98185009 : wk1 ! work variable
1219 :
1220 : character(len=*),parameter :: subname='(ridge_shift)'
1221 :
1222 5886032502 : do n = 1, ncat
1223 :
1224 : !-----------------------------------------------------------------
1225 : ! Save initial state variables
1226 : !-----------------------------------------------------------------
1227 :
1228 4900319557 : aicen_init(n) = aicen(n)
1229 4900319557 : vicen_init(n) = vicen(n)
1230 4900319557 : vsnon_init(n) = vsnon(n)
1231 :
1232 : !-----------------------------------------------------------------
1233 : ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
1234 : !-----------------------------------------------------------------
1235 :
1236 >17715*10^7 : do it = 1, ntrcr
1237 11145235464 : atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
1238 : + trcr_base(it,2) * vicen(n) & ! LCOV_EXCL_LINE
1239 >17126*10^7 : + trcr_base(it,3) * vsnon(n))
1240 >17616*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1241 35852787560 : do itl = 1, n_trcr_strata(it)
1242 22312784873 : ntr = nt_strata(it,itl)
1243 35852787560 : 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 985712945 : 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 985712945 : if (apartic(0) > c0) then
1264 893122201 : wk1 = apartic(0) * closing_gross * dt
1265 893122201 : if (wk1 > aice0) then
1266 17 : tmpfac = aice0 / wk1
1267 17 : closing_gross = closing_gross * tmpfac
1268 17 : 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 5886032502 : do n = 1, ncat
1277 5886032502 : if (aicen(n) > puny .and. apartic(n) > c0) then
1278 1104833130 : wk1 = apartic(n) * closing_gross * dt
1279 1104833130 : 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 985712945 : aice0 = aice0 - apartic(0)*closing_gross*dt + opning*dt
1292 :
1293 985712945 : 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 985712945 : elseif (aice0 < c0) then ! roundoff error
1301 0 : aice0 = c0
1302 : endif
1303 :
1304 985712945 : 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 5886032502 : do n = 1, ncat
1312 :
1313 : !-----------------------------------------------------------------
1314 : ! Identify grid cells with nonzero ridging
1315 : !-----------------------------------------------------------------
1316 :
1317 482854997 : if (aicen_init(n) > puny .and. apartic(n) > c0 &
1318 5886032502 : .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 1012554047 : ardg1n = apartic(n)*closing_gross*dt
1327 :
1328 1012554047 : 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 1012554047 : ardg1n = min(aicen_init(n), ardg1n)
1337 : endif
1338 :
1339 1012554047 : ardg2n = ardg1n / krdg(n)
1340 1012554047 : 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 1012554047 : virdgn = vicen_init(n) * afrac
1348 1012554047 : vsrdgn = vsnon_init(n) * afrac
1349 :
1350 1012554047 : aicen(n) = aicen(n) - ardg1n
1351 1012554047 : vicen(n) = vicen(n) - virdgn
1352 1012554047 : vsnon(n) = vsnon(n) - vsrdgn
1353 :
1354 : !-----------------------------------------------------------------
1355 : ! Increment ridging diagnostics
1356 : !-----------------------------------------------------------------
1357 :
1358 1012554047 : ardg1 = ardg1 + ardg1n
1359 1012554047 : ardg2 = ardg2 + ardg2n
1360 1012554047 : virdg = virdg + virdgn
1361 :
1362 1012554047 : ardg1nn(n) = ardg1n
1363 1012554047 : ardg2nn(n) = ardg2n
1364 1012554047 : virdgnn(n) = virdgn
1365 :
1366 : !-----------------------------------------------------------------
1367 : ! Place part of the snow and tracer lost by ridging into the ocean.
1368 : !-----------------------------------------------------------------
1369 :
1370 1012554047 : msnow_mlt = msnow_mlt + rhos*vsrdgn*(c1-fsnowrdg)
1371 :
1372 1012554047 : if (tr_aero) then
1373 43502110 : do it = 1, n_aero
1374 2041170 : maero(it) = maero(it) &
1375 : + vsrdgn*(c1-fsnowrdg) & ! LCOV_EXCL_LINE
1376 : *(trcrn(nt_aero +4*(it-1),n) & ! LCOV_EXCL_LINE
1377 45543280 : + trcrn(nt_aero+1+4*(it-1),n))
1378 : enddo
1379 : endif
1380 :
1381 1012554047 : if (tr_iso) then
1382 84420644 : do it = 1, n_iso
1383 588267 : miso(it) = miso(it) + vsrdgn*(c1-fsnowrdg) &
1384 : * (trcrn(nt_isosno+it-1,n) & ! LCOV_EXCL_LINE
1385 85008911 : + trcrn(nt_isoice+it-1,n))
1386 : enddo
1387 : endif
1388 :
1389 1012554047 : if (tr_pond_topo) then
1390 1981166 : mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
1391 15057621 : * 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 1012554047 : dhr = hrmax(n) - hrmin(n)
1400 1012554047 : 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 2090131960 : do k = 1, nslyr
1408 192436354 : esrdgn = vsrdgn * trcrn(nt_qsno+k-1,n) &
1409 1077577913 : / real(nslyr,kind=dbl_kind)
1410 2090131960 : esnow_mlt = esnow_mlt + esrdgn*(c1-fsnowrdg)
1411 : enddo
1412 :
1413 : !-----------------------------------------------------------------
1414 : ! Subtract area- and volume-weighted tracers from category n.
1415 : !-----------------------------------------------------------------
1416 :
1417 30322708303 : do it = 1, ntrcr
1418 :
1419 29310154256 : trfactor(1) = trcr_base(it,1)*ardg1n
1420 29310154256 : trfactor(2) = trcr_base(it,2)*virdgn
1421 29310154256 : trfactor(3) = trcr_base(it,3)*vsrdgn
1422 :
1423 29310154256 : work = c0
1424 >11724*10^7 : do k = 1, 3
1425 >11724*10^7 : work = work + trfactor(k)*trcrn(it,n)
1426 : enddo
1427 29310154256 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1428 6918848156 : do itl = 1, n_trcr_strata(it)
1429 4310944229 : ntr = nt_strata(it,itl)
1430 6918848156 : work = work * trcrn(ntr,n)
1431 : enddo
1432 : endif
1433 30322708303 : 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 6138560969 : do nr = 1, ncat
1442 :
1443 5126006922 : 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 140738115 : if (hrmin(n) >= hin_max(nr) .or. &
1451 4169250 : hrmax(n) <= hin_max(nr-1)) then
1452 28199424 : hL = c0
1453 28199424 : hR = c0
1454 : else
1455 108369441 : hL = max (hrmin(n), hin_max(nr-1))
1456 108369441 : hR = min (hrmax(n), hin_max(nr))
1457 : endif
1458 :
1459 136568865 : farea = (hR-hL) / dhr
1460 136568865 : 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 4989438057 : if (nr < ncat) then
1470 :
1471 4007566743 : hi1 = hrmin(n)
1472 4007566743 : hexp = hrexp(n)
1473 :
1474 4007566743 : if (hi1 >= hin_max(nr)) then
1475 2510632541 : farea = c0
1476 2510632541 : fvol = c0
1477 : else
1478 1496934202 : hL = max (hi1, hin_max(nr-1))
1479 1496934202 : hR = hin_max(nr)
1480 1496934202 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1481 1496934202 : expR_arg = min(((hR-hi1)/hexp),exp_argmax)
1482 1496934202 : expL = exp(-(expL_arg))
1483 1496934202 : expR = exp(-(expR_arg))
1484 1496934202 : farea = expL - expR
1485 : fvol = ((hL + hexp)*expL &
1486 1496934202 : - (hR + hexp)*expR) / (hi1 + hexp)
1487 : endif
1488 :
1489 : else ! nr = ncat
1490 :
1491 981871314 : hi1 = hrmin(n)
1492 981871314 : hexp = hrexp(n)
1493 :
1494 981871314 : hL = max (hi1, hin_max(nr-1))
1495 981871314 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1496 981871314 : expL = exp(-(expL_arg))
1497 981871314 : farea = expL
1498 981871314 : fvol = (hL + hexp)*expL / (hi1 + hexp)
1499 :
1500 : endif ! nr < ncat
1501 :
1502 : ! diagnostics
1503 4989438057 : if (n ==1) then ! only for thinnest ridging ice
1504 991870479 : if (present(aredistn)) then
1505 991870479 : aredistn(nr) = farea*ardg2n
1506 : endif
1507 991870479 : if (present(vredistn)) then
1508 991870479 : 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 5126006922 : aicen(nr) = aicen(nr) + farea*ardg2n
1519 5126006922 : vicen(nr) = vicen(nr) + fvol *virdgn
1520 5126006922 : 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 >15438*10^7 : do it = 1, ntrcr
1534 :
1535 >14824*10^7 : if (it /= nt_alvl .and. it /= nt_vlvl) then
1536 >13820*10^7 : trfactor(1) = trcr_base(it,1)*ardg2n*farea
1537 >13820*10^7 : trfactor(2) = trcr_base(it,2)*virdgn*fvol
1538 >13820*10^7 : trfactor(3) = trcr_base(it,3)*vsrdgn*fvol*fsnowrdg
1539 : else
1540 10045526958 : trfactor(1) = c0
1541 10045526958 : trfactor(2) = c0
1542 10045526958 : trfactor(3) = c0
1543 : endif
1544 :
1545 >14824*10^7 : work = c0
1546 >59299*10^7 : do k = 1, 3
1547 >59299*10^7 : work = work + trfactor(k)*trcrn(it,n)
1548 : enddo
1549 >14824*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1550 35174662192 : do itl = 1, n_trcr_strata(it)
1551 21909955717 : ntr = nt_strata(it,itl)
1552 35174662192 : if (ntr == nt_fbri) then ! brine fraction only
1553 0 : work = work * trcrn(ntr,n)
1554 : else
1555 21909955717 : work = c0
1556 : endif
1557 : enddo
1558 : endif
1559 >15337*10^7 : 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 5886032502 : 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 5383174554 : nt_strata, trcrn(:,n))
1579 5886032502 : 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 272773787 : 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 25534642 : aksum ! ratio of area removed to area ridged
1628 :
1629 : real (kind=dbl_kind), dimension (0:ncat) :: &
1630 721906666 : 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 423598237 : 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 25534642 : 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 272773787 : if (kstrength == 1) then ! Rothrock 1975 formulation
1651 :
1652 : !-----------------------------------------------------------------
1653 : ! Compute thickness distribution of ridging and ridged ice.
1654 : !-----------------------------------------------------------------
1655 :
1656 237630320 : call asum_ridging (ncat, aicen, aice0, asum)
1657 237630320 : 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 237630320 : hrexp, krdg)
1666 237630320 : 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 237630320 : 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 237630320 : elseif (krdg_redist==1) then ! exponential formulation
1686 :
1687 1444526362 : do n = 1, ncat
1688 1444526362 : if (aicen(n) > puny .and. apartic(n) > c0) then
1689 1180039854 : hi = vicen(n) / aicen(n)
1690 115713838 : h2rdg = hrmin(n)*hrmin(n) &
1691 : + c2*hrmin(n)*hrexp(n) & ! LCOV_EXCL_LINE
1692 1180039854 : + c2*hrexp(n)*hrexp(n)
1693 1180039854 : dh2rdg = -hi*hi + h2rdg/krdg(n)
1694 1180039854 : strength = strength + apartic(n) * dh2rdg
1695 : endif
1696 : enddo ! n
1697 :
1698 : endif ! krdg_redist
1699 :
1700 237630320 : 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 35143467 : 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 985712928 : 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 985712928 : flux_bio, closing )
1748 :
1749 : real (kind=dbl_kind), intent(in) :: &
1750 : dt ! time step
1751 :
1752 : integer (kind=int_kind), intent(in) :: &
1753 : ncat , & ! number of thickness categories ! LCOV_EXCL_LINE
1754 : ndtd , & ! number of dynamics supercycles ! LCOV_EXCL_LINE
1755 : nblyr , & ! number of bio layers ! LCOV_EXCL_LINE
1756 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
1757 : nslyr , & ! number of snow layers ! LCOV_EXCL_LINE
1758 : n_aero ! number of aerosol tracers
1759 :
1760 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
1761 : hin_max ! category limits (m)
1762 :
1763 : integer (kind=int_kind), dimension (:), intent(in) :: &
1764 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
1765 : n_trcr_strata ! number of underlying tracer layers
1766 :
1767 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1768 : trcr_base ! = 0 or 1 depending on tracer dependency
1769 : ! argument 2: (1) aice, (2) vice, (3) vsno
1770 :
1771 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
1772 : nt_strata ! indices of underlying tracer layers
1773 :
1774 : real (kind=dbl_kind), intent(inout) :: &
1775 : aice , & ! sea ice concentration ! LCOV_EXCL_LINE
1776 : aice0 , & ! concentration of open water ! LCOV_EXCL_LINE
1777 : rdg_conv , & ! convergence term for ridging (1/s) ! LCOV_EXCL_LINE
1778 : rdg_shear, & ! shear term for ridging (1/s) ! LCOV_EXCL_LINE
1779 : dardg1dt , & ! rate of area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
1780 : dardg2dt , & ! rate of area gain by new ridges (1/s) ! LCOV_EXCL_LINE
1781 : dvirdgdt , & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
1782 : opening , & ! rate of opening due to divergence/shear (1/s) ! LCOV_EXCL_LINE
1783 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
1784 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1785 : fsalt , & ! salt flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1786 : fhocn , & ! net heat flux to ocean (W/m^2) ! LCOV_EXCL_LINE
1787 : fzsal ! zsalinity flux to ocean(kg/m^2/s)
1788 :
1789 : real (kind=dbl_kind), intent(inout), optional :: &
1790 : closing ! rate of closing due to divergence/shear (1/s)
1791 :
1792 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1793 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1794 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1795 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
1796 : dardg1ndt, & ! rate of area loss by ridging ice (1/s) ! LCOV_EXCL_LINE
1797 : dardg2ndt, & ! rate of area gain by new ridges (1/s) ! LCOV_EXCL_LINE
1798 : dvirdgndt, & ! rate of ice volume ridged (m/s) ! LCOV_EXCL_LINE
1799 : aparticn , & ! participation function ! LCOV_EXCL_LINE
1800 : krdgn , & ! mean ridge thickness/thickness of ridging ice ! LCOV_EXCL_LINE
1801 : araftn , & ! rafting ice area ! LCOV_EXCL_LINE
1802 : vraftn , & ! rafting ice volume ! LCOV_EXCL_LINE
1803 : aredistn , & ! redistribution function: fraction of new ridge area ! LCOV_EXCL_LINE
1804 : vredistn , & ! redistribution function: fraction of new ridge volume ! LCOV_EXCL_LINE
1805 : faero_ocn, & ! aerosol flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
1806 : flux_bio ! all bio fluxes to ocean
1807 :
1808 : real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
1809 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
1810 :
1811 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
1812 : trcrn ! tracers
1813 :
1814 : !logical (kind=log_kind), intent(in) :: &
1815 : !tr_pond_topo,& ! if .true., use explicit topography-based ponds ! LCOV_EXCL_LINE
1816 : !tr_aero ,& ! if .true., use aerosol tracers ! LCOV_EXCL_LINE
1817 : !tr_brine !,& ! if .true., brine height differs from ice thickness ! LCOV_EXCL_LINE
1818 : !heat_capacity ! if true, ice has nonzero heat capacity
1819 :
1820 : logical (kind=log_kind), dimension(:), intent(inout) :: &
1821 : first_ice ! true until ice forms
1822 :
1823 : !autodocument_end
1824 :
1825 : ! local variables
1826 :
1827 : real (kind=dbl_kind) :: &
1828 98185008 : dtt ! thermo time step
1829 :
1830 : real (kind=dbl_kind), dimension(:), allocatable :: &
1831 985712928 : l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
1832 :
1833 : real (kind=dbl_kind) :: &
1834 98185008 : l_closing ! local rate of closing due to divergence/shear (1/s)
1835 :
1836 : logical (kind=log_kind) :: &
1837 : l_closing_flag ! flag if closing is passed
1838 :
1839 : character(len=*),parameter :: subname='(icepack_step_ridge)'
1840 :
1841 : !-----------------------------------------------------------------
1842 : ! Identify ice-ocean cells.
1843 : ! Note: We can not limit the loop here using aice>puny because
1844 : ! aice has not yet been updated since the transport (and
1845 : ! it may be out of whack, which the ridging helps fix).-ECH
1846 : !-----------------------------------------------------------------
1847 :
1848 985712928 : if (present(fiso_ocn)) then
1849 985712928 : allocate(l_fiso_ocn(size(fiso_ocn)))
1850 3942851712 : l_fiso_ocn = fiso_ocn
1851 : else
1852 : ! check tr_iso = true ???
1853 0 : allocate(l_fiso_ocn(1))
1854 0 : l_fiso_ocn = c0
1855 : endif
1856 :
1857 985712928 : if (present(closing)) then
1858 0 : l_closing_flag = .true.
1859 0 : l_closing = closing
1860 : else
1861 985712928 : l_closing_flag = .false.
1862 985712928 : l_closing = c0
1863 : endif
1864 :
1865 : call ridge_ice (dt, ndtd, &
1866 : ncat, n_aero, & ! LCOV_EXCL_LINE
1867 : nilyr, nslyr, & ! LCOV_EXCL_LINE
1868 : ntrcr, hin_max, & ! LCOV_EXCL_LINE
1869 : rdg_conv, rdg_shear, & ! LCOV_EXCL_LINE
1870 : aicen, & ! LCOV_EXCL_LINE
1871 : trcrn, & ! LCOV_EXCL_LINE
1872 : vicen, vsnon, & ! LCOV_EXCL_LINE
1873 : aice0, & ! LCOV_EXCL_LINE
1874 : trcr_depend, & ! LCOV_EXCL_LINE
1875 : trcr_base, & ! LCOV_EXCL_LINE
1876 : n_trcr_strata, & ! LCOV_EXCL_LINE
1877 : nt_strata, & ! LCOV_EXCL_LINE
1878 : krdg_partic, krdg_redist, & ! LCOV_EXCL_LINE
1879 : mu_rdg, tr_brine, & ! LCOV_EXCL_LINE
1880 : dardg1dt, dardg2dt, & ! LCOV_EXCL_LINE
1881 : dvirdgdt, opening, & ! LCOV_EXCL_LINE
1882 : fpond, & ! LCOV_EXCL_LINE
1883 : fresh, fhocn, & ! LCOV_EXCL_LINE
1884 : faero_ocn, l_fiso_ocn, & ! LCOV_EXCL_LINE
1885 : aparticn, krdgn, & ! LCOV_EXCL_LINE
1886 : aredistn, vredistn, & ! LCOV_EXCL_LINE
1887 : dardg1ndt, dardg2ndt, & ! LCOV_EXCL_LINE
1888 : dvirdgndt, & ! LCOV_EXCL_LINE
1889 : araftn, vraftn, & ! LCOV_EXCL_LINE
1890 : l_closing_flag, & ! LCOV_EXCL_LINE
1891 985712928 : l_closing )
1892 985712928 : if (icepack_warnings_aborted(subname)) return
1893 :
1894 : !-----------------------------------------------------------------
1895 : ! ITD cleanup: Rebin thickness categories if necessary, and remove
1896 : ! categories with very small areas.
1897 : !-----------------------------------------------------------------
1898 :
1899 985712928 : dtt = dt * ndtd ! for proper averaging over thermo timestep
1900 : call cleanup_itd (dtt, ntrcr, &
1901 : nilyr, nslyr, & ! LCOV_EXCL_LINE
1902 : ncat, hin_max, & ! LCOV_EXCL_LINE
1903 : aicen, trcrn, & ! LCOV_EXCL_LINE
1904 : vicen, vsnon, & ! LCOV_EXCL_LINE
1905 : aice0, aice, & ! LCOV_EXCL_LINE
1906 : n_aero, & ! LCOV_EXCL_LINE
1907 : nbtrcr, nblyr, & ! LCOV_EXCL_LINE
1908 : tr_aero, & ! LCOV_EXCL_LINE
1909 : tr_pond_topo, heat_capacity, & ! LCOV_EXCL_LINE
1910 : first_ice, & ! LCOV_EXCL_LINE
1911 : trcr_depend, trcr_base, & ! LCOV_EXCL_LINE
1912 : n_trcr_strata, nt_strata, & ! LCOV_EXCL_LINE
1913 : fpond, fresh, & ! LCOV_EXCL_LINE
1914 : fsalt, fhocn, & ! LCOV_EXCL_LINE
1915 : faero_ocn, l_fiso_ocn, & ! LCOV_EXCL_LINE
1916 : fzsal, & ! LCOV_EXCL_LINE
1917 985712928 : flux_bio)
1918 985712928 : if (icepack_warnings_aborted(subname)) return
1919 :
1920 3942851712 : if (present(fiso_ocn)) fiso_ocn = l_fiso_ocn
1921 985712928 : deallocate(l_fiso_ocn)
1922 :
1923 985712928 : end subroutine icepack_step_ridge
1924 :
1925 : !=======================================================================
1926 :
1927 : end module icepack_mechred
1928 :
1929 : !=======================================================================
|
