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