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, &
63 : ridge_itd, &
64 : icepack_ice_strength, &
65 : icepack_step_ridge
66 :
67 : real (kind=dbl_kind), parameter :: &
68 : exp_argmax = 100.0_dbl_kind, & ! maximum argument of exponential for underflow
69 : Cs = p25 , & ! fraction of shear energy contrbtng to ridging
70 : fsnowrdg = p5 , & ! snow fraction that survives in ridging
71 : Gstar = p15 , & ! max value of G(h) that participates
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
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 722805000 : subroutine ridge_ice (dt, ndtd, &
94 : ncat, n_aero, &
95 : nilyr, nslyr, &
96 722805000 : ntrcr, hin_max, &
97 : rdg_conv, rdg_shear, &
98 722805000 : aicen, trcrn, &
99 722805000 : vicen, vsnon, &
100 : aice0, &
101 722805000 : trcr_depend, trcr_base, &
102 722805000 : n_trcr_strata, &
103 722805000 : nt_strata, &
104 : krdg_partic, krdg_redist,&
105 : mu_rdg, tr_brine, &
106 : dardg1dt, dardg2dt, &
107 : dvirdgdt, opening, &
108 : fpond, &
109 : fresh, fhocn, &
110 722805000 : faero_ocn, fiso_ocn, &
111 722805000 : aparticn, krdgn, &
112 722805000 : aredistn, vredistn, &
113 722805000 : dardg1ndt, dardg2ndt, &
114 722805000 : dvirdgndt, &
115 722805000 : araftn, vraftn, &
116 : closing_flag,closing )
117 :
118 : integer (kind=int_kind), intent(in) :: &
119 : ndtd , & ! number of dynamics subcycles
120 : ncat , & ! number of thickness categories
121 : nilyr , & ! number of ice layers
122 : nslyr , & ! number of snow layers
123 : n_aero, & ! number of aerosol tracers
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)
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)
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
139 : vicen , & ! volume per unit area of ice (m)
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
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
161 : krdg_redist ! selects redistribution function
162 :
163 : logical (kind=log_kind), intent(in) :: &
164 : closing_flag, &! flag if closing is valid
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)
170 : dardg2dt , & ! rate of fractional area gain by new ridges (1/s)
171 : dvirdgdt , & ! rate of ice volume ridged (m/s)
172 : opening , & ! rate of opening due to divergence/shear (1/s)
173 : closing , & ! rate of closing due to divergence/shear (1/s)
174 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
175 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
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)
180 : dardg2ndt , & ! rate of fractional area gain by new ridges (1/s)
181 : dvirdgndt , & ! rate of ice volume ridged (m/s)
182 : aparticn , & ! participation function
183 : krdgn , & ! mean ridge thickness/thickness of ridging ice
184 : araftn , & ! rafting ice area
185 : vraftn , & ! rafting ice volume
186 : aredistn , & ! redistribution function: fraction of new ridge area
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 1618914432 : eicen ! energy of melting for each ice layer (J/m^2)
199 :
200 : real (kind=dbl_kind), dimension (ncat) :: &
201 1618914432 : esnon, & ! energy of melting for each snow layer (J/m^2)
202 1618914432 : vbrin, & ! ice volume with defined by brine height (m)
203 1618914432 : sicen ! Bulk salt in h ice (ppt*m)
204 :
205 : real (kind=dbl_kind) :: &
206 45487728 : asum , & ! sum of ice and open water area
207 45487728 : aksum , & ! ratio of area removed to area ridged
208 45487728 : msnow_mlt , & ! mass of snow added to ocean (kg m-2)
209 45487728 : esnow_mlt , & ! energy needed to melt snow in ocean (J m-2)
210 45487728 : mpond , & ! mass of pond added to ocean (kg m-2)
211 45487728 : closing_net, & ! net rate at which area is removed (1/s)
212 : ! (ridging ice area - area of new ridges) / dt
213 45487728 : divu_adv , & ! divu as implied by transport scheme (1/s)
214 45487728 : opning , & ! rate of opening due to divergence/shear
215 : ! opning is a local variable;
216 : ! opening is the history diagnostic variable
217 45487728 : ardg1 , & ! fractional area loss by ridging ice
218 45487728 : ardg2 , & ! fractional area gain by new ridges
219 45487728 : virdg , & ! ice volume ridged
220 45487728 : aopen ! area opening due to divergence/shear
221 :
222 : real (kind=dbl_kind), dimension (n_aero) :: &
223 1403196576 : maero ! aerosol mass added to ocean (kg m-2)
224 :
225 : real (kind=dbl_kind), dimension (n_iso) :: &
226 1400698704 : miso ! isotope mass added to ocean (kg m-2)
227 :
228 : real (kind=dbl_kind), dimension (0:ncat) :: &
229 1755377616 : apartic ! participation function; fraction of ridging
230 : ! and closing associated w/ category n
231 :
232 : real (kind=dbl_kind), dimension (ncat) :: &
233 1618914432 : hrmin , & ! minimum ridge thickness
234 1618914432 : hrmax , & ! maximum ridge thickness (krdg_redist = 0)
235 1618914432 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
236 1618914432 : krdg , & ! mean ridge thickness/thickness of ridging ice
237 1618914432 : ardg1n , & ! area of ice ridged
238 1618914432 : ardg2n , & ! area of new ridges
239 1618914432 : virdgn , & ! ridging ice volume
240 987084888 : mraftn ! rafting ice mask
241 :
242 : real (kind=dbl_kind) :: &
243 45487728 : vice_init, vice_final, & ! ice volume summed over categories
244 45487728 : vsno_init, vsno_final, & ! snow volume summed over categories
245 45487728 : eice_init, eice_final, & ! ice energy summed over layers
246 45487728 : vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories
247 45487728 : sice_init ,sice_final, & ! ice bulk salinity summed over categories
248 45487728 : 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
255 : niter , & ! iteration counter
256 : k , & ! vertical index
257 : it ! tracer index
258 :
259 : real (kind=dbl_kind) :: &
260 45487728 : 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 722805000 : msnow_mlt = c0
275 722805000 : esnow_mlt = c0
276 763347384 : maero (:) = c0
277 783906792 : miso (:) = c0
278 722805000 : mpond = c0
279 722805000 : ardg1 = c0
280 722805000 : ardg2 = c0
281 722805000 : virdg = c0
282 4264776000 : ardg1n(:) = c0
283 4264776000 : ardg2n(:) = c0
284 4264776000 : virdgn(:) = c0
285 4264776000 : mraftn(:) = c0
286 722805000 : aopen = c0
287 :
288 : !-----------------------------------------------------------------
289 : ! Compute area of ice plus open water before ridging.
290 : !-----------------------------------------------------------------
291 :
292 722805000 : call asum_ridging (ncat, aicen, aice0, asum)
293 722805000 : if (icepack_warnings_aborted(subname)) return
294 :
295 : !-----------------------------------------------------------------
296 : ! Compute the area opening and closing.
297 : !-----------------------------------------------------------------
298 :
299 722805000 : 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, &
309 : rdg_conv, rdg_shear, &
310 : asum, closing_net, &
311 722805000 : divu_adv, opning)
312 :
313 : endif
314 :
315 722805000 : if (icepack_warnings_aborted(subname)) return
316 :
317 : !-----------------------------------------------------------------
318 : ! Compute initial values of conserved quantities.
319 : !-----------------------------------------------------------------
320 :
321 722805000 : if (conserv_check) then
322 :
323 141225840 : do n = 1, ncat
324 121050720 : eicen(n) = c0
325 121050720 : esnon(n) = c0
326 121050720 : sicen(n) = c0
327 121050720 : vbrin(n) = c0
328 :
329 968405760 : do k = 1, nilyr
330 242101440 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
331 968405760 : * vicen(n)/real(nilyr,kind=dbl_kind)
332 242101440 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
333 1089456480 : * vicen(n)/real(nilyr,kind=dbl_kind)
334 : enddo
335 242101440 : do k = 1, nslyr
336 34585920 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
337 259394400 : * vsnon(n)/real(nslyr,kind=dbl_kind)
338 : enddo
339 121050720 : vbrin(n) = vicen(n)
340 141225840 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vicen(n)
341 : enddo ! n
342 :
343 : call column_sum (ncat, &
344 20175120 : vicen, vice_init)
345 20175120 : if (icepack_warnings_aborted(subname)) return
346 : call column_sum (ncat, &
347 20175120 : vsnon, vsno_init)
348 20175120 : if (icepack_warnings_aborted(subname)) return
349 : call column_sum (ncat, &
350 20175120 : eicen, eice_init)
351 20175120 : if (icepack_warnings_aborted(subname)) return
352 : call column_sum (ncat, &
353 20175120 : esnon, esno_init)
354 20175120 : if (icepack_warnings_aborted(subname)) return
355 : call column_sum (ncat, &
356 20175120 : sicen, sice_init)
357 20175120 : if (icepack_warnings_aborted(subname)) return
358 : call column_sum (ncat, &
359 20175120 : vbrin, vbri_init)
360 20175120 : if (icepack_warnings_aborted(subname)) return
361 :
362 : endif
363 :
364 722805006 : 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 0 : aicen, vicen, &
373 : krdg_partic, krdg_redist, &
374 : mu_rdg, &
375 : aksum, apartic, &
376 0 : hrmin, hrmax, &
377 0 : hrexp, krdg, &
378 0 : aparticn, krdgn, &
379 722805006 : mraftn)
380 722805006 : 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, &
388 0 : aicen, trcrn, &
389 0 : vicen, vsnon, &
390 0 : aice0, trcr_depend, &
391 0 : trcr_base, n_trcr_strata,&
392 0 : nt_strata, krdg_redist, &
393 : aksum, apartic, &
394 0 : hrmin, hrmax, &
395 0 : hrexp, krdg, &
396 : closing_net, opning, &
397 : ardg1, ardg2, &
398 : virdg, aopen, &
399 0 : ardg1n, ardg2n, &
400 0 : virdgn, &
401 : nslyr, n_aero, &
402 : msnow_mlt, esnow_mlt, &
403 0 : maero, miso, &
404 : mpond, &
405 722805006 : aredistn, vredistn)
406 722805006 : 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 722805006 : call asum_ridging (ncat, aicen, aice0, asum)
415 722805006 : if (icepack_warnings_aborted(subname)) return
416 :
417 722805006 : if (abs(asum - c1) < puny) then
418 722805000 : iterate_ridging = .false.
419 722805000 : closing_net = c0
420 722805000 : opning = c0
421 : else
422 6 : iterate_ridging = .true.
423 6 : divu_adv = (c1 - asum) / dt
424 6 : closing_net = max(c0, -divu_adv)
425 6 : opning = max(c0, divu_adv)
426 : endif
427 :
428 : !-----------------------------------------------------------------
429 : ! If done, exit. If not, prepare to ridge again.
430 : !-----------------------------------------------------------------
431 :
432 722805006 : if (iterate_ridging) then
433 6 : write(warnstr,*) subname, 'Repeat ridging, niter =', niter
434 6 : call icepack_warnings_add(warnstr)
435 : else
436 722805000 : exit rdg_iteration
437 : endif
438 :
439 6 : 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 722805000 : if (conserv_check) then
459 :
460 141225840 : do n = 1, ncat
461 121050720 : eicen(n) = c0
462 121050720 : esnon(n) = c0
463 121050720 : sicen(n) = c0
464 121050720 : vbrin(n) = c0
465 :
466 968405760 : do k = 1, nilyr
467 242101440 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
468 968405760 : * vicen(n)/real(nilyr,kind=dbl_kind)
469 242101440 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
470 1089456480 : * vicen(n)/real(nilyr,kind=dbl_kind)
471 : enddo
472 242101440 : do k = 1, nslyr
473 34585920 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
474 259394400 : * vsnon(n)/real(nslyr,kind=dbl_kind)
475 : enddo
476 121050720 : vbrin(n) = vicen(n)
477 141225840 : if (tr_brine) vbrin(n) = trcrn(nt_fbri,n) * vbrin(n)
478 : enddo ! n
479 :
480 : call column_sum (ncat, &
481 20175120 : vicen, vice_final)
482 20175120 : if (icepack_warnings_aborted(subname)) return
483 : call column_sum (ncat, &
484 20175120 : vsnon, vsno_final)
485 20175120 : if (icepack_warnings_aborted(subname)) return
486 : call column_sum (ncat, &
487 20175120 : eicen, eice_final)
488 20175120 : if (icepack_warnings_aborted(subname)) return
489 : call column_sum (ncat, &
490 20175120 : esnon, esno_final)
491 20175120 : if (icepack_warnings_aborted(subname)) return
492 : call column_sum (ncat, &
493 20175120 : sicen, sice_final)
494 20175120 : if (icepack_warnings_aborted(subname)) return
495 : call column_sum (ncat, &
496 20175120 : vbrin, vbri_final)
497 20175120 : if (icepack_warnings_aborted(subname)) return
498 :
499 20175120 : vsno_final = vsno_final + msnow_mlt/rhos
500 20175120 : esno_final = esno_final + esnow_mlt
501 :
502 20175120 : fieldid = subname//':vice'
503 : call column_conservation_check (fieldid, &
504 : vice_init, vice_final, &
505 20175120 : puny)
506 20175120 : if (icepack_warnings_aborted(subname)) return
507 20175120 : fieldid = subname//':vsno'
508 : call column_conservation_check (fieldid, &
509 : vsno_init, vsno_final, &
510 20175120 : puny)
511 20175120 : if (icepack_warnings_aborted(subname)) return
512 20175120 : fieldid = subname//':eice'
513 : call column_conservation_check (fieldid, &
514 : eice_init, eice_final, &
515 20175120 : puny*Lfresh*rhoi)
516 20175120 : if (icepack_warnings_aborted(subname)) return
517 20175120 : fieldid = subname//':esno'
518 : call column_conservation_check (fieldid, &
519 : esno_init, esno_final, &
520 20175120 : puny*Lfresh*rhos)
521 20175120 : if (icepack_warnings_aborted(subname)) return
522 20175120 : fieldid = subname//':sice'
523 : call column_conservation_check (fieldid, &
524 : sice_init, sice_final, &
525 20175120 : puny)
526 20175120 : if (icepack_warnings_aborted(subname)) return
527 20175120 : fieldid = subname//':vbrin'
528 : call column_conservation_check (fieldid, &
529 : vbri_init, vbri_final, &
530 20175120 : puny*c10)
531 20175120 : if (icepack_warnings_aborted(subname)) return
532 :
533 : endif ! conserv_check
534 :
535 : !-----------------------------------------------------------------
536 : ! Compute ridging diagnostics.
537 : !-----------------------------------------------------------------
538 :
539 722805000 : dti = c1/dt
540 :
541 722805000 : if (present(dardg1dt)) then
542 722805000 : dardg1dt = ardg1*dti
543 : endif
544 722805000 : if (present(dardg2dt)) then
545 722805000 : dardg2dt = ardg2*dti
546 : endif
547 722805000 : if (present(dvirdgdt)) then
548 722805000 : dvirdgdt = virdg*dti
549 : endif
550 722805000 : if (present(opening)) then
551 722805000 : opening = aopen*dti
552 : endif
553 722805000 : if (present(dardg1ndt)) then
554 4264776000 : do n = 1, ncat
555 4264776000 : dardg1ndt(n) = ardg1n(n)*dti
556 : enddo
557 : endif
558 722805000 : if (present(dardg2ndt)) then
559 4264776000 : do n = 1, ncat
560 4264776000 : dardg2ndt(n) = ardg2n(n)*dti
561 : enddo
562 : endif
563 722805000 : if (present(dvirdgndt)) then
564 4264776000 : do n = 1, ncat
565 4264776000 : dvirdgndt(n) = virdgn(n)*dti
566 : enddo
567 : endif
568 722805000 : if (present(araftn)) then
569 4264776000 : do n = 1, ncat
570 4264776000 : araftn(n) = mraftn(n)*ardg2n(n)
571 : ! araftn(n) = mraftn(n)*ardg1n(n)*p5
572 : enddo
573 : endif
574 722805000 : if (present(vraftn)) then
575 4264776000 : do n = 1, ncat
576 4264776000 : 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 722805000 : dti = c1/(ndtd*dt)
586 :
587 722805000 : if (present(fresh)) then
588 722805000 : fresh = fresh + msnow_mlt*dti
589 : endif
590 722805000 : if (present(fhocn)) then
591 722805000 : fhocn = fhocn + esnow_mlt*dti
592 : endif
593 722805000 : if (present(faero_ocn)) then
594 763347384 : do it = 1, n_aero
595 763347384 : faero_ocn(it) = faero_ocn(it) + maero(it)*dti
596 : enddo
597 : endif
598 722805000 : 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 722805000 : if (present(fpond)) then
605 722805000 : fpond = fpond - mpond ! units change later
606 : endif
607 :
608 : !-----------------------------------------------------------------
609 : ! Check for fractional ice area > 1.
610 : !-----------------------------------------------------------------
611 :
612 722805000 : 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 1580721216 : 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 1580721216 : asum = aice0
666 9345741534 : do n = 1, ncat
667 9345741534 : asum = asum + aicen(n)
668 : enddo
669 :
670 1580721216 : 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 722805000 : subroutine ridge_prep (dt, &
679 722805000 : ncat, hin_max, &
680 : rdg_conv, rdg_shear, &
681 : asum, closing_net, &
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)
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):: &
702 : closing_net, & ! net rate at which area is removed (1/s)
703 : divu_adv , & ! divu as implied by transport scheme (1/s)
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 722805000 : 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 722805000 : 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 722805000 : divu_adv = (c1-asum) / dt
752 :
753 722805000 : 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 722805000 : opning = closing_net + divu_adv
761 :
762 722805000 : 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 857916216 : subroutine ridge_itd (ncat, aice0, &
782 857916216 : aicen, vicen, &
783 : krdg_partic, krdg_redist, &
784 : mu_rdg, &
785 857916216 : aksum, apartic, &
786 857916216 : hrmin, hrmax, &
787 857916216 : hrexp, krdg, &
788 857916216 : aparticn, krdgn, &
789 857916216 : 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)
796 : aice0 ! concentration of open water
797 :
798 : real (kind=dbl_kind), dimension (:), intent(in) :: &
799 : aicen , & ! concentration of ice
800 : vicen ! volume per unit area of ice (m)
801 :
802 : integer (kind=int_kind), intent(in) :: &
803 : krdg_partic , & ! selects participation function
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
815 : hrmax , & ! maximum ridge thickness (krdg_redist = 0)
816 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
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
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, &
834 : astari = c1/astar
835 :
836 : real (kind=dbl_kind), dimension(-1:ncat) :: &
837 2028224133 : Gsum ! Gsum(n) = sum of areas in categories 0 to n
838 :
839 : real (kind=dbl_kind) :: &
840 53388122 : work ! temporary work variable
841 :
842 : real (kind=dbl_kind) :: &
843 53388122 : hi , & ! ice thickness for each cat (m)
844 53388122 : hrmean , & ! mean ridge thickness (m)
845 53388122 : xtmp ! temporary variable
846 :
847 : character(len=*),parameter :: subname='(ridge_itd)'
848 :
849 : !-----------------------------------------------------------------
850 : ! Initialize
851 : !-----------------------------------------------------------------
852 :
853 6796797966 : Gsum (-1:ncat) = c0 ! initialize
854 :
855 857916216 : Gsum (-1) = c0 ! by definition
856 857916216 : if (aice0 > puny) then
857 715982073 : Gsum(0) = aice0
858 : else
859 141934143 : Gsum(0) = Gsum(-1)
860 : endif
861 857916216 : apartic(0) = c0
862 :
863 5080965534 : do n = 1, ncat
864 4223049318 : apartic(n) = c0
865 4223049318 : hrmin (n) = c0
866 4223049318 : hrmax (n) = c0
867 4223049318 : hrexp (n) = c0
868 4223049318 : 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 5080965534 : if (aicen(n) > puny) then
881 1457091495 : Gsum(n) = Gsum(n-1) + aicen(n)
882 : else
883 2765957823 : Gsum(n) = Gsum(n-1)
884 : endif
885 : enddo
886 :
887 : ! normalize
888 :
889 857916216 : if (Gsum(ncat) > c0) then
890 857916216 : work = c1 / Gsum(ncat)
891 5938881750 : do n = 0, ncat
892 5938881750 : 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 857916216 : 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 2734954 : apartic(n) = Gstari*(Gsum(n ) - Gsum(n-1)) * &
917 34071530 : (c2 - Gstari*(Gsum(n-1) + Gsum(n )))
918 245321355 : 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 796207128 : 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 796207128 : xtmp = c1 / (c1 - exp(-astari))
935 796207128 : Gsum(-1) = exp(-Gsum(-1)*astari) * xtmp
936 5599147254 : do n = 0, ncat
937 4802940126 : Gsum(n) = exp(-Gsum(n)*astari) * xtmp
938 5599147254 : 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 857916216 : 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 105526110 : hi = vicen(n) / aicen(n)
969 105526110 : hrmin(n) = min(c2*hi, hi + maxraft)
970 105526110 : hrmax(n) = c2*sqrt(Hstar*hi)
971 105526110 : hrmax(n) = max(hrmax(n), hrmin(n)+puny)
972 105526110 : hrmean = p5 * (hrmin(n) + hrmax(n))
973 105526110 : 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 4802940126 : do n = 1, ncat
1010 4802940126 : if (aicen(n) > puny) then
1011 1351565385 : hi = vicen(n) / aicen(n)
1012 1351565385 : hi = max(hi,puny)
1013 1351565385 : hrmin(n) = min(c2*hi, hi + maxraft)
1014 1351565385 : hrexp(n) = mu_rdg * sqrt(hi)
1015 1351565385 : krdg(n) = (hrmin(n) + hrexp(n)) / hi
1016 :
1017 : !echmod: check computational efficiency
1018 : ! diagnostic rafting mask
1019 1351565385 : if (present(mraft)) then
1020 686205973 : mraft(n) = max(c0, sign(c1, hi+maxraft-hrmin(n)))
1021 686205973 : xtmp = mraft(n)*((c2*hi+hrexp(n))/hi - krdg(n))
1022 686205973 : 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 857916216 : aksum = apartic(0) ! area participating = area removed
1037 :
1038 5080965534 : do n = 1, ncat
1039 : ! area participating > area removed
1040 5080965534 : aksum = aksum + apartic(n) * (c1 - c1/krdg(n))
1041 : enddo
1042 :
1043 : ! diagnostics
1044 857916216 : if (present(aparticn)) then
1045 4264776036 : do n = 1, ncat
1046 4264776036 : aparticn(n) = apartic(n)
1047 : enddo
1048 : endif
1049 857916216 : if (present(krdgn)) then
1050 4264776036 : do n = 1, ncat
1051 4264776036 : krdgn(n) = krdg(n)
1052 : enddo
1053 : endif
1054 :
1055 857916216 : 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 722805006 : subroutine ridge_shift (ntrcr, dt, &
1074 722805006 : ncat, hin_max, &
1075 1445610012 : aicen, trcrn, &
1076 722805006 : vicen, vsnon, &
1077 722805006 : aice0, trcr_depend, &
1078 722805006 : trcr_base, n_trcr_strata, &
1079 722805006 : nt_strata, krdg_redist, &
1080 722805006 : aksum, apartic, &
1081 722805006 : hrmin, hrmax, &
1082 1445610012 : hrexp, krdg, &
1083 : closing_net, opning, &
1084 : ardg1, ardg2, &
1085 : virdg, aopen, &
1086 722805006 : ardg1nn, ardg2nn, &
1087 722805006 : virdgnn, &
1088 : nslyr, n_aero, &
1089 : msnow_mlt, esnow_mlt, &
1090 722805006 : maero, miso, &
1091 : mpond, &
1092 722805006 : aredistn, vredistn)
1093 :
1094 : integer (kind=int_kind), intent(in) :: &
1095 : ncat , & ! number of thickness categories
1096 : nslyr , & ! number of snow layers
1097 : ntrcr , & ! number of tracers in use
1098 : n_aero, & ! number of aerosol tracers
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
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
1123 : vicen , & ! volume per unit area of ice (m)
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
1138 : hrmax , & ! maximum ridge thickness (krdg_redist = 0)
1139 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
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)
1144 : opning , & ! rate of opening due to divergence/shear (1/s)
1145 : ardg1 , & ! fractional area loss by ridging ice
1146 : ardg2 , & ! fractional area gain by new ridges
1147 : virdg , & ! ice volume ridged (m)
1148 : aopen ! area opened due to divergence/shear
1149 :
1150 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1151 : ardg1nn , & ! area of ice ridged
1152 : ardg2nn , & ! area of new ridges
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)
1157 : esnow_mlt , & ! energy needed to melt snow in ocean (J m-2)
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
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
1174 : k , & ! ice layer index
1175 : it , & ! tracer index
1176 : ntr , & ! tracer index
1177 : itl ! loop index
1178 :
1179 : real (kind=dbl_kind), dimension (ncat) :: &
1180 1709889900 : aicen_init , & ! ice area before ridging
1181 1618914444 : vicen_init , & ! ice volume before ridging
1182 1618914444 : vsnon_init ! snow volume before ridging
1183 :
1184 : real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
1185 5753693454 : atrcrn ! aicen*trcrn
1186 :
1187 : real (kind=dbl_kind), dimension(3) :: &
1188 181950912 : trfactor ! base quantity on which tracers are carried
1189 :
1190 : real (kind=dbl_kind) :: &
1191 45487728 : work , & ! temporary variable
1192 45487728 : expL_arg , & ! temporary exp arg values
1193 45487728 : expR_arg , & ! temporary exp arg values
1194 45487728 : 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 45487728 : afrac , & ! fraction of category area ridged
1200 45487728 : ardg1n , & ! area of ice ridged
1201 45487728 : ardg2n , & ! area of new ridges
1202 45487728 : virdgn , & ! ridging ice volume
1203 45487728 : vsrdgn , & ! ridging snow volume
1204 45487728 : dhr , & ! hrmax - hrmin
1205 45487728 : dhr2 , & ! hrmax^2 - hrmin^2
1206 45487728 : farea , & ! fraction of new ridge area going to nr
1207 45487728 : fvol ! fraction of new ridge volume going to nr
1208 :
1209 : real (kind=dbl_kind) :: &
1210 45487728 : esrdgn ! ridging snow energy
1211 :
1212 : real (kind=dbl_kind) :: &
1213 45487728 : hi1 , & ! thickness of ridging ice
1214 45487728 : hexp , & ! ridge e-folding thickness
1215 45487728 : hL, hR , & ! left and right limits of integration
1216 90975456 : expL, expR , & ! exponentials involving hL, hR
1217 45487728 : tmpfac , & ! factor by which opening/closing rates are cut
1218 45487728 : wk1 ! work variable
1219 :
1220 : character(len=*),parameter :: subname='(ridge_shift)'
1221 :
1222 4264776036 : do n = 1, ncat
1223 :
1224 : !-----------------------------------------------------------------
1225 : ! Save initial state variables
1226 : !-----------------------------------------------------------------
1227 :
1228 3541971030 : aicen_init(n) = aicen(n)
1229 3541971030 : vicen_init(n) = vicen(n)
1230 3541971030 : vsnon_init(n) = vsnon(n)
1231 :
1232 : !-----------------------------------------------------------------
1233 : ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
1234 : !-----------------------------------------------------------------
1235 :
1236 >14233*10^7 : do it = 1, ntrcr
1237 4766608560 : atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
1238 4766608560 : + trcr_base(it,2) * vicen(n) &
1239 >13807*10^7 : + trcr_base(it,3) * vsnon(n))
1240 >14161*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1241 25059782640 : do itl = 1, n_trcr_strata(it)
1242 15576379710 : ntr = nt_strata(it,itl)
1243 25059782640 : 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 722805006 : 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 722805006 : if (apartic(0) > c0) then
1264 651665638 : wk1 = apartic(0) * closing_gross * dt
1265 651665638 : if (wk1 > aice0) then
1266 6 : tmpfac = aice0 / wk1
1267 6 : closing_gross = closing_gross * tmpfac
1268 6 : 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 4264776036 : do n = 1, ncat
1277 4264776036 : if (aicen(n) > puny .and. apartic(n) > c0) then
1278 718451701 : wk1 = apartic(n) * closing_gross * dt
1279 718451701 : 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 722805006 : aice0 = aice0 - apartic(0)*closing_gross*dt + opning*dt
1292 :
1293 722805006 : 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 722805006 : elseif (aice0 < c0) then ! roundoff error
1301 0 : aice0 = c0
1302 : endif
1303 :
1304 722805006 : 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 4264776036 : do n = 1, ncat
1312 :
1313 : !-----------------------------------------------------------------
1314 : ! Identify grid cells with nonzero ridging
1315 : !-----------------------------------------------------------------
1316 :
1317 218792160 : if (aicen_init(n) > puny .and. apartic(n) > c0 &
1318 4264776036 : .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 601147321 : ardg1n = apartic(n)*closing_gross*dt
1327 :
1328 601147321 : 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 601147321 : ardg1n = min(aicen_init(n), ardg1n)
1337 : endif
1338 :
1339 601147321 : ardg2n = ardg1n / krdg(n)
1340 601147321 : 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 601147321 : virdgn = vicen_init(n) * afrac
1348 601147321 : vsrdgn = vsnon_init(n) * afrac
1349 :
1350 601147321 : aicen(n) = aicen(n) - ardg1n
1351 601147321 : vicen(n) = vicen(n) - virdgn
1352 601147321 : vsnon(n) = vsnon(n) - vsrdgn
1353 :
1354 : !-----------------------------------------------------------------
1355 : ! Increment ridging diagnostics
1356 : !-----------------------------------------------------------------
1357 :
1358 601147321 : ardg1 = ardg1 + ardg1n
1359 601147321 : ardg2 = ardg2 + ardg2n
1360 601147321 : virdg = virdg + virdgn
1361 :
1362 601147321 : ardg1nn(n) = ardg1n
1363 601147321 : ardg2nn(n) = ardg2n
1364 601147321 : virdgnn(n) = virdgn
1365 :
1366 : !-----------------------------------------------------------------
1367 : ! Place part of the snow and tracer lost by ridging into the ocean.
1368 : !-----------------------------------------------------------------
1369 :
1370 601147321 : msnow_mlt = msnow_mlt + rhos*vsrdgn*(c1-fsnowrdg)
1371 :
1372 601147321 : if (tr_aero) then
1373 37862990 : do it = 1, n_aero
1374 1819848 : maero(it) = maero(it) &
1375 : + vsrdgn*(c1-fsnowrdg) &
1376 3639696 : *(trcrn(nt_aero +4*(it-1),n) &
1377 39682838 : + trcrn(nt_aero+1+4*(it-1),n))
1378 : enddo
1379 : endif
1380 :
1381 601147321 : if (tr_iso) then
1382 84119616 : do it = 1, n_iso
1383 588897 : miso(it) = miso(it) + vsrdgn*(c1-fsnowrdg) &
1384 1177794 : * (trcrn(nt_isosno+it-1,n) &
1385 84708513 : + trcrn(nt_isoice+it-1,n))
1386 : enddo
1387 : endif
1388 :
1389 601147321 : if (tr_pond_topo) then
1390 1760017 : mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
1391 12320119 : * 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 601147321 : dhr = hrmax(n) - hrmin(n)
1400 601147321 : 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 1202294642 : do k = 1, nslyr
1408 66790842 : esrdgn = vsrdgn * trcrn(nt_qsno+k-1,n) &
1409 601147321 : / real(nslyr,kind=dbl_kind)
1410 1202294642 : esnow_mlt = esnow_mlt + esrdgn*(c1-fsnowrdg)
1411 : enddo
1412 :
1413 : !-----------------------------------------------------------------
1414 : ! Subtract area- and volume-weighted tracers from category n.
1415 : !-----------------------------------------------------------------
1416 :
1417 21605264938 : do it = 1, ntrcr
1418 :
1419 21004117617 : trfactor(1) = trcr_base(it,1)*ardg1n
1420 21004117617 : trfactor(2) = trcr_base(it,2)*virdgn
1421 21004117617 : trfactor(3) = trcr_base(it,3)*vsrdgn
1422 :
1423 21004117617 : work = c0
1424 84016470468 : do k = 1, 3
1425 84016470468 : work = work + trfactor(k)*trcrn(it,n)
1426 : enddo
1427 21004117617 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1428 4433221134 : do itl = 1, n_trcr_strata(it)
1429 2758936791 : ntr = nt_strata(it,itl)
1430 4433221134 : work = work * trcrn(ntr,n)
1431 : enddo
1432 : endif
1433 21605264938 : 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 3602910029 : do nr = 1, ncat
1442 :
1443 3001762708 : 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 140089252 : if (hrmin(n) >= hin_max(nr) .or. &
1451 4148663 : hrmax(n) <= hin_max(nr-1)) then
1452 28300587 : hL = c0
1453 28300587 : hR = c0
1454 : else
1455 107640002 : hL = max (hrmin(n), hin_max(nr-1))
1456 107640002 : hR = min (hrmax(n), hin_max(nr))
1457 : endif
1458 :
1459 135940589 : farea = (hR-hL) / dhr
1460 135940589 : 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 2865822119 : if (nr < ncat) then
1470 :
1471 2295121719 : hi1 = hrmin(n)
1472 2295121719 : hexp = hrexp(n)
1473 :
1474 2295121719 : if (hi1 >= hin_max(nr)) then
1475 1427106041 : farea = c0
1476 1427106041 : fvol = c0
1477 : else
1478 868015678 : hL = max (hi1, hin_max(nr-1))
1479 868015678 : hR = hin_max(nr)
1480 868015678 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1481 868015678 : expR_arg = min(((hR-hi1)/hexp),exp_argmax)
1482 868015678 : expL = exp(-(expL_arg))
1483 868015678 : expR = exp(-(expR_arg))
1484 868015678 : farea = expL - expR
1485 : fvol = ((hL + hexp)*expL &
1486 868015678 : - (hR + hexp)*expR) / (hi1 + hexp)
1487 : endif
1488 :
1489 : else ! nr = ncat
1490 :
1491 570700400 : hi1 = hrmin(n)
1492 570700400 : hexp = hrexp(n)
1493 :
1494 570700400 : hL = max (hi1, hin_max(nr-1))
1495 570700400 : expL_arg = min(((hL-hi1)/hexp),exp_argmax)
1496 570700400 : expL = exp(-(expL_arg))
1497 570700400 : farea = expL
1498 570700400 : fvol = (hL + hexp)*expL / (hi1 + hexp)
1499 :
1500 : endif ! nr < ncat
1501 :
1502 : ! diagnostics
1503 2865822119 : if (n ==1) then ! only for thinnest ridging ice
1504 577367448 : if (present(aredistn)) then
1505 577367448 : aredistn(nr) = farea*ardg2n
1506 : endif
1507 577367448 : if (present(vredistn)) then
1508 577367448 : 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 3001762708 : aicen(nr) = aicen(nr) + farea*ardg2n
1519 3001762708 : vicen(nr) = vicen(nr) + fvol *virdgn
1520 3001762708 : 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 >10857*10^7 : do it = 1, ntrcr
1534 :
1535 >10497*10^7 : if (it /= nt_alvl .and. it /= nt_vlvl) then
1536 99157819068 : trfactor(1) = trcr_base(it,1)*ardg2n*farea
1537 99157819068 : trfactor(2) = trcr_base(it,2)*virdgn*fvol
1538 99157819068 : trfactor(3) = trcr_base(it,3)*vsrdgn*fvol*fsnowrdg
1539 : else
1540 5816277537 : trfactor(1) = c0
1541 5816277537 : trfactor(2) = c0
1542 5816277537 : trfactor(3) = c0
1543 : endif
1544 :
1545 >10497*10^7 : work = c0
1546 >41989*10^7 : do k = 1, 3
1547 >41989*10^7 : work = work + trfactor(k)*trcrn(it,n)
1548 : enddo
1549 >10497*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1550 22215386146 : do itl = 1, n_trcr_strata(it)
1551 13819324193 : ntr = nt_strata(it,itl)
1552 22215386146 : if (ntr == nt_fbri) then ! brine fraction only
1553 0 : work = work * trcrn(ntr,n)
1554 : else
1555 13819324193 : work = c0
1556 : endif
1557 : enddo
1558 : endif
1559 >10797*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 4264776036 : do n = 1, ncat
1574 0 : call icepack_compute_tracers (ntrcr, trcr_depend, &
1575 218792160 : atrcrn(:,n), aicen(n), &
1576 218792160 : vicen(n), vsnon(n), &
1577 0 : trcr_base, n_trcr_strata, &
1578 3760763190 : nt_strata, trcrn(:,n))
1579 4264776036 : 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 169432898 : subroutine icepack_ice_strength (ncat, &
1601 : aice, vice, &
1602 169432898 : aice0, aicen, &
1603 169432898 : vicen, &
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
1611 : vice , & ! volume per unit area of ice (m)
1612 : aice0 ! concentration of open water
1613 :
1614 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1615 : aicen , & ! concentration of ice
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 9993109 : asum , & ! sum of ice and open water area
1627 9993109 : aksum ! ratio of area removed to area ridged
1628 :
1629 : real (kind=dbl_kind), dimension (0:ncat) :: &
1630 406323856 : apartic ! participation function; fraction of ridging
1631 : ! and closing associated w/ category n
1632 :
1633 : real (kind=dbl_kind), dimension (ncat) :: &
1634 376344529 : hrmin , & ! minimum ridge thickness
1635 376344529 : hrmax , & ! maximum ridge thickness (krdg_redist = 0)
1636 376344529 : hrexp , & ! ridge e-folding thickness (krdg_redist = 1)
1637 226897849 : 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 9993109 : hi , & ! ice thickness (m)
1644 9993109 : h2rdg , & ! mean value of h^2 for new ridge
1645 9993109 : 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 169432898 : if (kstrength == 1) then ! Rothrock 1975 formulation
1651 :
1652 : !-----------------------------------------------------------------
1653 : ! Compute thickness distribution of ridging and ridged ice.
1654 : !-----------------------------------------------------------------
1655 :
1656 135111210 : call asum_ridging (ncat, aicen, aice0, asum)
1657 135111210 : if (icepack_warnings_aborted(subname)) return
1658 :
1659 : call ridge_itd (ncat, aice0, &
1660 0 : aicen, vicen, &
1661 : krdg_partic, krdg_redist, &
1662 : mu_rdg, &
1663 : aksum, apartic, &
1664 0 : hrmin, hrmax, &
1665 135111210 : hrexp, krdg)
1666 135111210 : 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 135111210 : 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 135111210 : elseif (krdg_redist==1) then ! exponential formulation
1686 :
1687 816189498 : do n = 1, ncat
1688 816189498 : if (aicen(n) > puny .and. apartic(n) > c0) then
1689 665359412 : hi = vicen(n) / aicen(n)
1690 39310594 : h2rdg = hrmin(n)*hrmin(n) &
1691 39310594 : + c2*hrmin(n)*hrexp(n) &
1692 665359412 : + c2*hrexp(n)*hrexp(n)
1693 665359412 : dh2rdg = -hi*hi + h2rdg/krdg(n)
1694 665359412 : strength = strength + apartic(n) * dh2rdg
1695 : endif
1696 : enddo ! n
1697 :
1698 : endif ! krdg_redist
1699 :
1700 135111210 : 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 34321688 : 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 722805000 : subroutine icepack_step_ridge (dt, ndtd, &
1724 : nilyr, nslyr, &
1725 : nblyr, &
1726 722805000 : ncat, hin_max, &
1727 : rdg_conv, rdg_shear, &
1728 722805000 : aicen, &
1729 722805000 : trcrn, &
1730 722805000 : vicen, vsnon, &
1731 722805000 : aice0, trcr_depend, &
1732 722805000 : trcr_base, n_trcr_strata, &
1733 722805000 : nt_strata, &
1734 : dardg1dt, dardg2dt, &
1735 : dvirdgdt, opening, &
1736 : fpond, &
1737 : fresh, fhocn, &
1738 : n_aero, &
1739 1445610000 : faero_ocn, fiso_ocn, &
1740 722805000 : aparticn, krdgn, &
1741 722805000 : aredistn, vredistn, &
1742 722805000 : dardg1ndt, dardg2ndt, &
1743 722805000 : dvirdgndt, &
1744 722805000 : araftn, vraftn, &
1745 : aice, fsalt, &
1746 722805000 : first_ice, fzsal, &
1747 722805000 : 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
1754 : ndtd , & ! number of dynamics supercycles
1755 : nblyr , & ! number of bio layers
1756 : nilyr , & ! number of ice layers
1757 : nslyr , & ! number of snow layers
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
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
1776 : aice0 , & ! concentration of open water
1777 : rdg_conv , & ! convergence term for ridging (1/s)
1778 : rdg_shear, & ! shear term for ridging (1/s)
1779 : dardg1dt , & ! rate of area loss by ridging ice (1/s)
1780 : dardg2dt , & ! rate of area gain by new ridges (1/s)
1781 : dvirdgdt , & ! rate of ice volume ridged (m/s)
1782 : opening , & ! rate of opening due to divergence/shear (1/s)
1783 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
1784 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
1785 : fsalt , & ! salt flux to ocean (kg/m^2/s)
1786 : fhocn , & ! net heat flux to ocean (W/m^2)
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
1794 : vicen , & ! volume per unit area of ice (m)
1795 : vsnon , & ! volume per unit area of snow (m)
1796 : dardg1ndt, & ! rate of area loss by ridging ice (1/s)
1797 : dardg2ndt, & ! rate of area gain by new ridges (1/s)
1798 : dvirdgndt, & ! rate of ice volume ridged (m/s)
1799 : aparticn , & ! participation function
1800 : krdgn , & ! mean ridge thickness/thickness of ridging ice
1801 : araftn , & ! rafting ice area
1802 : vraftn , & ! rafting ice volume
1803 : aredistn , & ! redistribution function: fraction of new ridge area
1804 : vredistn , & ! redistribution function: fraction of new ridge volume
1805 : faero_ocn, & ! aerosol flux to ocean (kg/m^2/s)
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
1816 : !tr_aero ,& ! if .true., use aerosol tracers
1817 : !tr_brine !,& ! if .true., brine height differs from ice thickness
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 45487728 : dtt ! thermo time step
1829 :
1830 : real (kind=dbl_kind), dimension(:), allocatable :: &
1831 722805000 : l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
1832 :
1833 : real (kind=dbl_kind) :: &
1834 45487728 : 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 722805000 : if (present(fiso_ocn)) then
1849 722805000 : allocate(l_fiso_ocn(size(fiso_ocn)))
1850 2891220000 : 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 722805000 : if (present(closing)) then
1858 0 : l_closing_flag = .true.
1859 0 : l_closing = closing
1860 : else
1861 722805000 : l_closing_flag = .false.
1862 722805000 : l_closing = c0
1863 : endif
1864 :
1865 : call ridge_ice (dt, ndtd, &
1866 : ncat, n_aero, &
1867 : nilyr, nslyr, &
1868 : ntrcr, hin_max, &
1869 : rdg_conv, rdg_shear, &
1870 0 : aicen, &
1871 0 : trcrn, &
1872 0 : vicen, vsnon, &
1873 : aice0, &
1874 0 : trcr_depend, &
1875 0 : trcr_base, &
1876 0 : n_trcr_strata, &
1877 0 : nt_strata, &
1878 : krdg_partic, krdg_redist, &
1879 : mu_rdg, tr_brine, &
1880 : dardg1dt, dardg2dt, &
1881 : dvirdgdt, opening, &
1882 : fpond, &
1883 : fresh, fhocn, &
1884 0 : faero_ocn, l_fiso_ocn, &
1885 0 : aparticn, krdgn, &
1886 0 : aredistn, vredistn, &
1887 0 : dardg1ndt, dardg2ndt, &
1888 0 : dvirdgndt, &
1889 0 : araftn, vraftn, &
1890 : l_closing_flag, &
1891 722805000 : l_closing )
1892 722805000 : 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 722805000 : dtt = dt * ndtd ! for proper averaging over thermo timestep
1900 : call cleanup_itd (dtt, ntrcr, &
1901 : nilyr, nslyr, &
1902 : ncat, hin_max, &
1903 0 : aicen, trcrn, &
1904 0 : vicen, vsnon, &
1905 : aice0, aice, &
1906 : n_aero, &
1907 : nbtrcr, nblyr, &
1908 : tr_aero, &
1909 : tr_pond_topo, heat_capacity, &
1910 0 : first_ice, &
1911 0 : trcr_depend, trcr_base, &
1912 0 : n_trcr_strata, nt_strata, &
1913 : fpond, fresh, &
1914 : fsalt, fhocn, &
1915 0 : faero_ocn, l_fiso_ocn, &
1916 : fzsal, &
1917 722805000 : flux_bio)
1918 722805000 : if (icepack_warnings_aborted(subname)) return
1919 :
1920 2891220000 : if (present(fiso_ocn)) fiso_ocn = l_fiso_ocn
1921 722805000 : deallocate(l_fiso_ocn)
1922 :
1923 722805000 : end subroutine icepack_step_ridge
1924 :
1925 : !=======================================================================
1926 :
1927 : end module icepack_mechred
1928 :
1929 : !=======================================================================
|
