Line data Source code
1 : !=======================================================================
2 : !
3 : ! Thermo calculations after call to coupler, related to ITD:
4 : ! ice thickness redistribution, lateral growth and melting.
5 : !
6 : ! NOTE: The thermodynamic calculation is split in two for load balancing.
7 : ! First icepack_therm_vertical computes vertical growth rates and coupler
8 : ! fluxes. Then icepack_therm_itd does thermodynamic calculations not
9 : ! needed for coupling.
10 : !
11 : ! authors William H. Lipscomb, LANL
12 : ! C. M. Bitz, UW
13 : ! Elizabeth C. Hunke, LANL
14 : !
15 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
16 : ! 2004: Block structure added by William Lipscomb
17 : ! 2006: Streamlined for efficiency by Elizabeth Hunke
18 : ! 2014: Column package created by Elizabeth Hunke
19 : !
20 : module icepack_therm_itd
21 :
22 : use icepack_kinds
23 : use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10
24 : use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
25 : use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
26 : use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
27 : use icepack_parameters, only: rhosi, conserv_check
28 : use icepack_parameters, only: kitd, ktherm, heat_capacity
29 : use icepack_parameters, only: z_tracers, solve_zsal
30 :
31 : use icepack_tracers, only: ntrcr, nbtrcr
32 : use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
33 : use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
34 : use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd
35 : use icepack_tracers, only: nt_alvl, nt_vlvl
36 : use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
37 : use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
38 : use icepack_tracers, only: n_aero, n_iso
39 : use icepack_tracers, only: bio_index
40 :
41 : use icepack_warnings, only: warnstr, icepack_warnings_add
42 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
43 :
44 : use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd, get_subdt_fsd
45 : use icepack_itd, only: reduce_area, cleanup_itd
46 : use icepack_itd, only: aggregate_area, shift_ice
47 : use icepack_itd, only: column_sum, column_conservation_check
48 : use icepack_isotope, only: isoice_alpha, isotope_frac_method
49 : use icepack_mushy_physics, only: liquidus_temperature_mush, enthalpy_mush
50 : use icepack_therm_shared, only: hfrazilmin
51 : use icepack_therm_shared, only: hi_min
52 : use icepack_zbgc, only: add_new_ice_bgc
53 : use icepack_zbgc, only: lateral_melt_bgc
54 :
55 : implicit none
56 :
57 : private
58 : public :: linear_itd, &
59 : add_new_ice, &
60 : lateral_melt, &
61 : icepack_step_therm2
62 :
63 : !=======================================================================
64 :
65 : contains
66 :
67 : !=======================================================================
68 : !
69 : ! ITD scheme that shifts ice among categories
70 : !
71 : ! See Lipscomb, W. H. Remapping the thickness distribution in sea
72 : ! ice models. 2001, J. Geophys. Res., Vol 106, 13989--14000.
73 : !
74 : ! Using the thermodynamic "velocities", interpolate to find the
75 : ! velocities in thickness space at the category boundaries, and
76 : ! compute the new locations of the boundaries. Then for each
77 : ! category, compute the thickness distribution function, g(h),
78 : ! between hL and hR, the left and right boundaries of the category.
79 : ! Assume g(h) is a linear polynomial that satisfies two conditions:
80 : !
81 : ! (1) The ice area implied by g(h) equals aicen(n).
82 : ! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
83 : !
84 : ! Given g(h), at each boundary compute the ice area and volume lying
85 : ! between the original and new boundary locations. Transfer area
86 : ! and volume across each boundary in the appropriate direction, thus
87 : ! restoring the original boundaries.
88 : !
89 : ! authors: William H. Lipscomb, LANL
90 : ! Elizabeth C. Hunke, LANL
91 :
92 144085983 : subroutine linear_itd (ncat, hin_max, &
93 : nilyr, nslyr, &
94 144085983 : ntrcr, trcr_depend, &
95 144085983 : trcr_base, n_trcr_strata,&
96 144085983 : nt_strata, &
97 144085983 : aicen_init, vicen_init, &
98 288171966 : aicen, trcrn, &
99 144085983 : vicen, vsnon, &
100 : aice, aice0, &
101 : fpond )
102 :
103 : integer (kind=int_kind), intent(in) :: &
104 : ncat , & ! number of thickness categories
105 : nilyr , & ! number of ice layers
106 : nslyr , & ! number of snow layers
107 : ntrcr ! number of tracers in use
108 :
109 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
110 : hin_max ! category boundaries (m)
111 :
112 : integer (kind=int_kind), dimension (:), intent(in) :: &
113 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
114 : n_trcr_strata ! number of underlying tracer layers
115 :
116 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
117 : trcr_base ! = 0 or 1 depending on tracer dependency
118 : ! argument 2: (1) aice, (2) vice, (3) vsno
119 :
120 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
121 : nt_strata ! indices of underlying tracer layers
122 :
123 : real (kind=dbl_kind), dimension(:), intent(in) :: &
124 : aicen_init, & ! initial ice concentration (before vertical thermo)
125 : vicen_init ! initial ice volume (m)
126 :
127 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
128 : aicen , & ! ice concentration
129 : vicen , & ! volume per unit area of ice (m)
130 : vsnon ! volume per unit area of snow (m)
131 :
132 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
133 : trcrn ! ice tracers
134 :
135 : real (kind=dbl_kind), intent(inout) :: &
136 : aice , & ! concentration of ice
137 : aice0 , & ! concentration of open water
138 : fpond ! fresh water flux to ponds (kg/m^2/s)
139 :
140 : ! local variables
141 :
142 : integer (kind=int_kind) :: &
143 : n, nd , & ! category indices
144 : k ! ice layer index
145 :
146 : real (kind=dbl_kind) :: &
147 7910849 : slope , & ! rate of change of dhice with hice
148 7910849 : dh0 , & ! change in ice thickness at h = 0
149 7910849 : da0 , & ! area melting from category 1
150 7910849 : damax , & ! max allowed reduction in category 1 area
151 7910849 : etamin, etamax,& ! left and right limits of integration
152 7910849 : x1 , & ! etamax - etamin
153 7910849 : x2 , & ! (etamax^2 - etamin^2) / 2
154 7910849 : x3 , & ! (etamax^3 - etamin^3) / 3
155 7910849 : wk1, wk2 ! temporary variables
156 :
157 : real (kind=dbl_kind), dimension(0:ncat) :: &
158 328395313 : hbnew ! new boundary locations
159 :
160 : real (kind=dbl_kind), dimension(ncat) :: &
161 320484464 : g0 , & ! constant coefficient in g(h)
162 320484464 : g1 , & ! linear coefficient in g(h)
163 320484464 : hL , & ! left end of range over which g(h) > 0
164 320484464 : hR ! right end of range over which g(h) > 0
165 :
166 : real (kind=dbl_kind), dimension(ncat) :: &
167 320484464 : hicen , & ! ice thickness for each cat (m)
168 320484464 : hicen_init , & ! initial ice thickness for each cat (pre-thermo)
169 320484464 : dhicen , & ! thickness change for remapping (m)
170 336306162 : daice , & ! ice area transferred across boundary
171 320484464 : dvice ! ice volume transferred across boundary
172 :
173 : real (kind=dbl_kind), dimension (ncat) :: &
174 320484464 : eicen, & ! energy of melting for each ice layer (J/m^2)
175 320484464 : esnon, & ! energy of melting for each snow layer (J/m^2)
176 320484464 : vbrin, & ! ice volume defined by brine height (m)
177 320484464 : sicen ! Bulk salt in h ice (ppt*m)
178 :
179 : real (kind=dbl_kind) :: &
180 7910849 : vice_init, vice_final, & ! ice volume summed over categories
181 7910849 : vsno_init, vsno_final, & ! snow volume summed over categories
182 7910849 : eice_init, eice_final, & ! ice energy summed over categories
183 7910849 : esno_init, esno_final, & ! snow energy summed over categories
184 7910849 : sice_init, sice_final, & ! ice bulk salinity summed over categories
185 7910849 : vbri_init, vbri_final ! briny ice volume summed over categories
186 :
187 : ! NOTE: Third index of donor, daice, dvice should be ncat-1,
188 : ! except that compilers would have trouble when ncat = 1
189 : integer (kind=int_kind), dimension(ncat) :: &
190 151996832 : donor ! donor category index
191 :
192 : logical (kind=log_kind) :: &
193 : remap_flag ! remap ITD if remap_flag is true
194 :
195 : character (len=char_len) :: &
196 : fieldid ! field identifier
197 :
198 : logical (kind=log_kind), parameter :: &
199 : print_diags = .false. ! if true, prints when remap_flag=F
200 :
201 : character(len=*),parameter :: subname='(linear_itd)'
202 :
203 : !-----------------------------------------------------------------
204 : ! Initialize
205 : !-----------------------------------------------------------------
206 :
207 144085983 : hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
208 :
209 869831871 : do n = 1, ncat
210 725745888 : donor(n) = 0
211 725745888 : daice(n) = c0
212 869831871 : dvice(n) = c0
213 : enddo
214 :
215 : !-----------------------------------------------------------------
216 : ! Compute volume and energy sums that linear remapping should
217 : ! conserve.
218 : !-----------------------------------------------------------------
219 :
220 144085983 : if (conserv_check) then
221 :
222 37211811 : do n = 1, ncat
223 :
224 31895838 : eicen(n) = c0
225 31895838 : esnon(n) = c0
226 31895838 : vbrin(n) = c0
227 31895838 : sicen(n) = c0
228 :
229 255166704 : do k = 1, nilyr
230 56204568 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
231 283268988 : * vicen(n)/real(nilyr,kind=dbl_kind)
232 : enddo
233 63791676 : do k = 1, nslyr
234 8029224 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
235 67806288 : * vsnon(n)/real(nslyr,kind=dbl_kind)
236 : enddo
237 :
238 31895838 : if (tr_brine) then
239 0 : vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
240 0 : * vicen(n)/real(nilyr,kind=dbl_kind)
241 : endif
242 :
243 260482677 : do k = 1, nilyr
244 56204568 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
245 283268988 : * vicen(n)/real(nilyr,kind=dbl_kind)
246 : enddo
247 :
248 : enddo ! n
249 :
250 5315973 : call column_sum (ncat, vicen, vice_init)
251 5315973 : if (icepack_warnings_aborted(subname)) return
252 5315973 : call column_sum (ncat, vsnon, vsno_init)
253 5315973 : if (icepack_warnings_aborted(subname)) return
254 5315973 : call column_sum (ncat, eicen, eice_init)
255 5315973 : if (icepack_warnings_aborted(subname)) return
256 5315973 : call column_sum (ncat, esnon, esno_init)
257 5315973 : if (icepack_warnings_aborted(subname)) return
258 5315973 : call column_sum (ncat, sicen, sice_init)
259 5315973 : if (icepack_warnings_aborted(subname)) return
260 5315973 : call column_sum (ncat, vbrin, vbri_init)
261 5315973 : if (icepack_warnings_aborted(subname)) return
262 :
263 : endif ! conserv_check
264 :
265 : !-----------------------------------------------------------------
266 : ! Initialize remapping flag.
267 : ! Remapping is done wherever remap_flag = .true.
268 : ! In rare cases the category boundaries may shift too far for the
269 : ! remapping algorithm to work, and remap_flag is set to .false.
270 : ! In these cases the simpler 'rebin' subroutine will shift ice
271 : ! between categories if needed.
272 : !-----------------------------------------------------------------
273 :
274 144085983 : remap_flag = .true.
275 :
276 : !-----------------------------------------------------------------
277 : ! Compute thickness change in each category.
278 : !-----------------------------------------------------------------
279 :
280 869831871 : do n = 1, ncat
281 :
282 725745888 : if (aicen_init(n) > puny) then
283 687178834 : hicen_init(n) = vicen_init(n) / aicen_init(n)
284 : else
285 38567054 : hicen_init(n) = c0
286 : endif ! aicen_init > puny
287 :
288 869831871 : if (aicen (n) > puny) then
289 687126294 : hicen (n) = vicen(n) / aicen(n)
290 687126294 : dhicen(n) = hicen(n) - hicen_init(n)
291 : else
292 38619594 : hicen (n) = c0
293 38619594 : dhicen(n) = c0
294 : endif ! aicen > puny
295 :
296 : enddo ! n
297 :
298 : !-----------------------------------------------------------------
299 : ! Compute new category boundaries, hbnew, based on changes in
300 : ! ice thickness from vertical thermodynamics.
301 : !-----------------------------------------------------------------
302 :
303 144085983 : hbnew(0) = hin_max(0)
304 :
305 725745888 : do n = 1, ncat-1
306 :
307 581659905 : if (hicen_init(n) > puny .and. &
308 : hicen_init(n+1) > puny) then
309 : ! interpolate between adjacent category growth rates
310 29173336 : slope = (dhicen(n+1) - dhicen(n)) / &
311 571019098 : (hicen_init(n+1) - hicen_init(n))
312 29173336 : hbnew(n) = hin_max(n) + dhicen(n) &
313 541845762 : + slope * (hin_max(n) - hicen_init(n))
314 39814143 : elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0
315 19244842 : hbnew(n) = hin_max(n) + dhicen(n)
316 20569301 : elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0
317 1578952 : hbnew(n) = hin_max(n) + dhicen(n+1)
318 : else
319 18990349 : hbnew(n) = hin_max(n)
320 : endif
321 :
322 : !-----------------------------------------------------------------
323 : ! Check that each boundary lies between adjacent values of hicen.
324 : ! If not, set remap_flag = .false.
325 : ! Write diagnosis outputs if remap_flag was changed to false
326 : !-----------------------------------------------------------------
327 :
328 581659905 : if (aicen(n) > puny .and. hicen(n) >= hbnew(n)) then
329 14 : remap_flag = .false.
330 :
331 14 : if (print_diags) then
332 : write(warnstr,*) subname, 'ITD: hicen(n) > hbnew(n)'
333 : call icepack_warnings_add(warnstr)
334 : write(warnstr,*) subname, 'cat ',n
335 : call icepack_warnings_add(warnstr)
336 : write(warnstr,*) subname, 'hicen(n) =', hicen(n)
337 : call icepack_warnings_add(warnstr)
338 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
339 : call icepack_warnings_add(warnstr)
340 : endif
341 :
342 581659891 : elseif (aicen(n+1) > puny .and. hicen(n+1) <= hbnew(n)) then
343 0 : remap_flag = .false.
344 :
345 : if (print_diags) then
346 : write(warnstr,*) subname, 'ITD: hicen(n+1) < hbnew(n)'
347 : call icepack_warnings_add(warnstr)
348 : write(warnstr,*) subname, 'cat ',n
349 : call icepack_warnings_add(warnstr)
350 : write(warnstr,*) subname, 'hicen(n+1) =', hicen(n+1)
351 : call icepack_warnings_add(warnstr)
352 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
353 : call icepack_warnings_add(warnstr)
354 : endif
355 : endif
356 :
357 : !-----------------------------------------------------------------
358 : ! Check that hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
359 : ! If not, set remap_flag = .false.
360 : ! (In principle we could allow this, but it would make the code
361 : ! more complicated.)
362 : ! Write diagnosis outputs if remap_flag was changed to false
363 : !-----------------------------------------------------------------
364 :
365 581659905 : if (hbnew(n) > hin_max(n+1)) then
366 0 : remap_flag = .false.
367 :
368 : if (print_diags) then
369 : write(warnstr,*) subname, 'ITD hbnew(n) > hin_max(n+1)'
370 : call icepack_warnings_add(warnstr)
371 : write(warnstr,*) subname, 'cat ',n
372 : call icepack_warnings_add(warnstr)
373 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
374 : call icepack_warnings_add(warnstr)
375 : write(warnstr,*) subname, 'hin_max(n+1) =', hin_max(n+1)
376 : call icepack_warnings_add(warnstr)
377 : endif
378 : endif
379 :
380 725745888 : if (hbnew(n) < hin_max(n-1)) then
381 0 : remap_flag = .false.
382 :
383 : if (print_diags) then
384 : write(warnstr,*) subname, 'ITD: hbnew(n) < hin_max(n-1)'
385 : call icepack_warnings_add(warnstr)
386 : write(warnstr,*) subname, 'cat ',n
387 : call icepack_warnings_add(warnstr)
388 : write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
389 : call icepack_warnings_add(warnstr)
390 : write(warnstr,*) subname, 'hin_max(n-1) =', hin_max(n-1)
391 : call icepack_warnings_add(warnstr)
392 : endif
393 : endif
394 :
395 : enddo ! boundaries, 1 to ncat-1
396 :
397 : !-----------------------------------------------------------------
398 : ! Compute hbnew(ncat)
399 : !-----------------------------------------------------------------
400 :
401 144085983 : if (aicen(ncat) > puny) then
402 126088230 : hbnew(ncat) = c3*hicen(ncat) - c2*hbnew(ncat-1)
403 : else
404 17997753 : hbnew(ncat) = hin_max(ncat)
405 : endif
406 144085983 : hbnew(ncat) = max(hbnew(ncat),hin_max(ncat-1))
407 :
408 : !-----------------------------------------------------------------
409 : ! Identify cells where the ITD is to be remapped
410 : !-----------------------------------------------------------------
411 :
412 144085983 : if (remap_flag) then
413 :
414 : !-----------------------------------------------------------------
415 : ! Compute g(h) for category 1 at start of time step
416 : ! (hicen = hicen_init)
417 : !-----------------------------------------------------------------
418 :
419 7910849 : call fit_line(aicen(1), hicen_init(1), &
420 7910849 : hbnew(0), hin_max (1), &
421 7910849 : g0 (1), g1 (1), &
422 144085969 : hL (1), hR (1))
423 144085969 : if (icepack_warnings_aborted(subname)) return
424 :
425 : !-----------------------------------------------------------------
426 : ! Find area lost due to melting of thin (category 1) ice
427 : !-----------------------------------------------------------------
428 :
429 144085969 : if (aicen(1) > puny) then
430 :
431 143701566 : dh0 = dhicen(1)
432 143701566 : if (dh0 < c0) then ! remove area from category 1
433 57723941 : dh0 = min(-dh0,hin_max(1)) ! dh0 --> |dh0|
434 :
435 : !-----------------------------------------------------------------
436 : ! Integrate g(1) from 0 to dh0 to estimate area melted
437 : !-----------------------------------------------------------------
438 :
439 : ! right integration limit (left limit = 0)
440 57723941 : etamax = min(dh0,hR(1)) - hL(1)
441 :
442 57723941 : if (etamax > c0) then
443 53548946 : x1 = etamax
444 53548946 : x2 = p5 * etamax*etamax
445 53548946 : da0 = g1(1)*x2 + g0(1)*x1 ! ice area removed
446 :
447 : ! constrain new thickness <= hicen_init
448 53548946 : damax = aicen(1) * (c1-hicen(1)/hicen_init(1)) ! damax > 0
449 53548946 : da0 = min (da0, damax)
450 :
451 : ! remove area, conserving volume
452 53548946 : hicen(1) = hicen(1) * aicen(1) / (aicen(1)-da0)
453 53548946 : aicen(1) = aicen(1) - da0
454 :
455 53548946 : if (tr_pond_topo) &
456 97208 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
457 715076 : * trcrn(nt_hpnd,1))
458 :
459 : endif ! etamax > 0
460 :
461 : else ! dh0 >= 0
462 85977625 : hbnew(0) = min(dh0,hin_max(1)) ! shift hbnew(0) to right
463 : endif
464 :
465 : endif ! aicen(1) > puny
466 :
467 : !-----------------------------------------------------------------
468 : ! Compute g(h) for each ice thickness category.
469 : !-----------------------------------------------------------------
470 :
471 869831787 : do n = 1, ncat
472 :
473 40223347 : call fit_line(aicen(n), hicen(n), &
474 40223347 : hbnew(n-1), hbnew(n), &
475 40223347 : g0 (n), g1 (n), &
476 765969165 : hL (n), hR (n))
477 725745818 : if (icepack_warnings_aborted(subname)) return
478 :
479 : !-----------------------------------------------------------------
480 : ! Compute area and volume to be shifted across each boundary.
481 : !-----------------------------------------------------------------
482 :
483 725745818 : donor(n) = 0
484 725745818 : daice(n) = c0
485 910055134 : dvice(n) = c0
486 : enddo
487 :
488 725745818 : do n = 1, ncat-1
489 :
490 581659849 : if (hbnew(n) > hin_max(n)) then ! transfer from n to n+1
491 :
492 : ! left and right integration limits in eta space
493 278792114 : etamin = max(hin_max(n), hL(n)) - hL(n)
494 278792114 : etamax = min(hbnew(n), hR(n)) - hL(n)
495 278792114 : donor(n) = n
496 :
497 : else ! hbnew(n) <= hin_max(n); transfer from n+1 to n
498 :
499 : ! left and right integration limits in eta space
500 302867735 : etamin = c0
501 302867735 : etamax = min(hin_max(n), hR(n+1)) - hL(n+1)
502 302867735 : donor(n) = n+1
503 :
504 : endif ! hbnew(n) > hin_max(n)
505 :
506 581659849 : if (etamax > etamin) then
507 470821448 : x1 = etamax - etamin
508 470821448 : wk1 = etamin*etamin
509 470821448 : wk2 = etamax*etamax
510 470821448 : x2 = p5 * (wk2 - wk1)
511 470821448 : wk1 = wk1*etamin
512 470821448 : wk2 = wk2*etamax
513 470821448 : x3 = p333 * (wk2 - wk1)
514 470821448 : nd = donor(n)
515 470821448 : daice(n) = g1(nd)*x2 + g0(nd)*x1
516 470821448 : dvice(n) = g1(nd)*x3 + g0(nd)*x2 + daice(n)*hL(nd)
517 : endif ! etamax > etamin
518 :
519 : ! If daice or dvice is very small, shift no ice.
520 :
521 581659849 : nd = donor(n)
522 :
523 581659849 : if (daice(n) < aicen(nd)*puny) then
524 77877689 : daice(n) = c0
525 77877689 : dvice(n) = c0
526 77877689 : donor(n) = 0
527 : endif
528 :
529 581659849 : if (dvice(n) < vicen(nd)*puny) then
530 77877834 : daice(n) = c0
531 77877834 : dvice(n) = c0
532 77877834 : donor(n) = 0
533 : endif
534 :
535 : ! If daice is close to aicen or dvice is close to vicen,
536 : ! shift entire category
537 :
538 581659849 : if (daice(n) > aicen(nd)*(c1-puny)) then
539 46786 : daice(n) = aicen(nd)
540 46786 : dvice(n) = vicen(nd)
541 : endif
542 :
543 725745818 : if (dvice(n) > vicen(nd)*(c1-puny)) then
544 46786 : daice(n) = aicen(nd)
545 46786 : dvice(n) = vicen(nd)
546 : endif
547 :
548 : enddo ! boundaries, 1 to ncat-1
549 :
550 : !-----------------------------------------------------------------
551 : ! Shift ice between categories as necessary
552 : !-----------------------------------------------------------------
553 :
554 : ! maintain qsno negative definiteness
555 869831787 : do n = 1, ncat
556 1595577605 : do k = nt_qsno, nt_qsno+nslyr-1
557 1451491636 : trcrn(k,n) = trcrn(k,n) + rhos*Lfresh
558 : enddo
559 : enddo
560 :
561 : call shift_ice (ntrcr, ncat, &
562 0 : trcr_depend, &
563 0 : trcr_base, &
564 0 : n_trcr_strata, &
565 0 : nt_strata, &
566 0 : aicen, trcrn, &
567 0 : vicen, vsnon, &
568 0 : hicen, donor, &
569 144085969 : daice, dvice )
570 144085969 : if (icepack_warnings_aborted(subname)) return
571 :
572 : ! maintain qsno negative definiteness
573 869831787 : do n = 1, ncat
574 1595577605 : do k = nt_qsno, nt_qsno+nslyr-1
575 1451491636 : trcrn(k,n) = trcrn(k,n) - rhos*Lfresh
576 : enddo
577 : enddo
578 :
579 : !-----------------------------------------------------------------
580 : ! Make sure hice(1) >= minimum ice thickness hi_min.
581 : !-----------------------------------------------------------------
582 :
583 144085969 : if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then
584 :
585 708777 : da0 = aicen(1) * (c1 - hicen(1)/hi_min)
586 708777 : aicen(1) = aicen(1) - da0
587 708777 : hicen(1) = hi_min
588 :
589 708777 : if (tr_pond_topo) &
590 1708 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
591 13732 : * trcrn(nt_hpnd,1))
592 : endif
593 :
594 : endif ! remap_flag
595 :
596 : !-----------------------------------------------------------------
597 : ! Update fractional ice area in each grid cell.
598 : !-----------------------------------------------------------------
599 :
600 144085983 : call aggregate_area (ncat, aicen, aice, aice0)
601 144085983 : if (icepack_warnings_aborted(subname)) return
602 :
603 : !-----------------------------------------------------------------
604 : ! Check volume and energy conservation.
605 : !-----------------------------------------------------------------
606 :
607 144085983 : if (conserv_check) then
608 :
609 37211811 : do n = 1, ncat
610 :
611 31895838 : eicen(n) = c0
612 31895838 : esnon(n) = c0
613 31895838 : vbrin(n) = c0
614 31895838 : sicen(n) = c0
615 :
616 255166704 : do k = 1, nilyr
617 56204568 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
618 283268988 : * vicen(n)/real(nilyr,kind=dbl_kind)
619 : enddo
620 63791676 : do k = 1, nslyr
621 8029224 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
622 67806288 : * vsnon(n)/real(nslyr,kind=dbl_kind)
623 : enddo
624 :
625 31895838 : if (tr_brine) then
626 0 : vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
627 0 : * vicen(n)/real(nilyr,kind=dbl_kind)
628 : endif
629 :
630 260482677 : do k = 1, nilyr
631 56204568 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
632 283268988 : * vicen(n)/real(nilyr,kind=dbl_kind)
633 : enddo
634 :
635 : enddo ! n
636 :
637 5315973 : call column_sum (ncat, vicen, vice_final)
638 5315973 : if (icepack_warnings_aborted(subname)) return
639 5315973 : call column_sum (ncat, vsnon, vsno_final)
640 5315973 : if (icepack_warnings_aborted(subname)) return
641 5315973 : call column_sum (ncat, eicen, eice_final)
642 5315973 : if (icepack_warnings_aborted(subname)) return
643 5315973 : call column_sum (ncat, esnon, esno_final)
644 5315973 : if (icepack_warnings_aborted(subname)) return
645 5315973 : call column_sum (ncat, sicen, sice_final)
646 5315973 : if (icepack_warnings_aborted(subname)) return
647 5315973 : call column_sum (ncat, vbrin, vbri_final)
648 5315973 : if (icepack_warnings_aborted(subname)) return
649 :
650 5315973 : fieldid = subname//':vice'
651 : call column_conservation_check (fieldid, &
652 : vice_init, vice_final, &
653 5315973 : puny)
654 5315973 : if (icepack_warnings_aborted(subname)) return
655 5315973 : fieldid = subname//':vsno'
656 : call column_conservation_check (fieldid, &
657 : vsno_init, vsno_final, &
658 5315973 : puny)
659 5315973 : if (icepack_warnings_aborted(subname)) return
660 5315973 : fieldid = subname//':eice'
661 : call column_conservation_check (fieldid, &
662 : eice_init, eice_final, &
663 5315973 : puny*Lfresh*rhoi)
664 5315973 : if (icepack_warnings_aborted(subname)) return
665 5315973 : fieldid = subname//':esno'
666 : call column_conservation_check (fieldid, &
667 : esno_init, esno_final, &
668 5315973 : puny*Lfresh*rhos)
669 5315973 : if (icepack_warnings_aborted(subname)) return
670 5315973 : fieldid = subname//':sicen'
671 : call column_conservation_check (fieldid, &
672 : sice_init, sice_final, &
673 5315973 : puny)
674 5315973 : if (icepack_warnings_aborted(subname)) return
675 5315973 : fieldid = subname//':vbrin'
676 : call column_conservation_check (fieldid, &
677 : vbri_init, vbri_final, &
678 5315973 : puny*c10)
679 5315973 : if (icepack_warnings_aborted(subname)) return
680 :
681 : endif ! conservation check
682 :
683 : end subroutine linear_itd
684 :
685 : !=======================================================================
686 : !
687 : ! Fit g(h) with a line, satisfying area and volume constraints.
688 : ! To reduce roundoff errors caused by large values of g0 and g1,
689 : ! we actually compute g(eta), where eta = h - hL, and hL is the
690 : ! left boundary.
691 : !
692 : ! authors: William H. Lipscomb, LANL
693 : ! Elizabeth C. Hunke, LANL
694 :
695 869831787 : subroutine fit_line (aicen, hice, &
696 : hbL, hbR, &
697 : g0, g1, &
698 : hL, hR)
699 :
700 : real (kind=dbl_kind), intent(in) :: &
701 : aicen ! concentration of ice
702 :
703 : real (kind=dbl_kind), intent(in) :: &
704 : hbL, hbR , & ! left and right category boundaries
705 : hice ! ice thickness
706 :
707 : real (kind=dbl_kind), intent(out):: &
708 : g0, g1 , & ! coefficients in linear equation for g(eta)
709 : hL , & ! min value of range over which g(h) > 0
710 : hR ! max value of range over which g(h) > 0
711 :
712 : ! local variables
713 :
714 : real (kind=dbl_kind) :: &
715 48134196 : h13 , & ! hbL + 1/3 * (hbR - hbL)
716 48134196 : h23 , & ! hbL + 2/3 * (hbR - hbL)
717 48134196 : dhr , & ! 1 / (hR - hL)
718 48134196 : wk1, wk2 ! temporary variables
719 :
720 : character(len=*),parameter :: subname='(fit_line)'
721 :
722 : !-----------------------------------------------------------------
723 : ! Compute g0, g1, hL, and hR for each category to be remapped.
724 : !-----------------------------------------------------------------
725 :
726 869831787 : if (aicen > puny .and. hbR - hbL > puny) then
727 :
728 : ! Initialize hL and hR
729 :
730 830821462 : hL = hbL
731 830821462 : hR = hbR
732 :
733 : ! Change hL or hR if hicen(n) falls outside central third of range
734 :
735 830821462 : h13 = p333 * (c2*hL + hR)
736 830821462 : h23 = p333 * (hL + c2*hR)
737 830821462 : if (hice < h13) then
738 158723234 : hR = c3*hice - c2*hL
739 672098228 : elseif (hice > h23) then
740 77224271 : hL = c3*hice - c2*hR
741 : endif
742 :
743 : ! Compute coefficients of g(eta) = g0 + g1*eta
744 :
745 830821462 : dhr = c1 / (hR - hL)
746 830821462 : wk1 = c6 * aicen * dhr
747 830821462 : wk2 = (hice - hL) * dhr
748 830821462 : g0 = wk1 * (p666 - wk2)
749 830821462 : g1 = c2*dhr * wk1 * (wk2 - p5)
750 :
751 : else
752 :
753 39010325 : g0 = c0
754 39010325 : g1 = c0
755 39010325 : hL = c0
756 39010325 : hR = c0
757 :
758 : endif ! aicen > puny
759 :
760 869831787 : end subroutine fit_line
761 :
762 : !=======================================================================
763 : !
764 : ! Given some added new ice to the base of the existing ice, recalculate
765 : ! vertical tracer so that new grid cells are all the same size.
766 : !
767 : ! author: A. K. Turner, LANL
768 : !
769 467164780 : subroutine update_vertical_tracers(nilyr, trc, h1, h2, trc0)
770 :
771 : integer (kind=int_kind), intent(in) :: &
772 : nilyr ! number of ice layers
773 :
774 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
775 : trc ! vertical tracer
776 :
777 : real (kind=dbl_kind), intent(in) :: &
778 : h1, & ! old thickness
779 : h2, & ! new thickness
780 : trc0 ! tracer value of added ice on ice bottom
781 :
782 : ! local variables
783 :
784 1050160028 : real(kind=dbl_kind), dimension(nilyr) :: trc2 ! updated tracer temporary
785 :
786 : ! vertical indices for old and new grid
787 : integer :: k1, k2
788 :
789 : real (kind=dbl_kind) :: &
790 19305078 : z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom
791 19305078 : z2a, z2b, & ! upper, lower boundary of new cell
792 19305078 : overlap , & ! overlap between old and new cell
793 19305078 : rnilyr
794 :
795 : character(len=*),parameter :: subname='(update_vertical_tracers)'
796 :
797 467164780 : rnilyr = real(nilyr,dbl_kind)
798 :
799 : ! loop over new grid cells
800 3737318240 : do k2 = 1, nilyr
801 :
802 : ! initialize new tracer
803 3270153460 : trc2(k2) = c0
804 :
805 : ! calculate upper and lower boundary of new cell
806 3270153460 : z2a = ((k2 - 1) * h2) / rnilyr
807 3270153460 : z2b = (k2 * h2) / rnilyr
808 :
809 : ! loop over old grid cells
810 26161227680 : do k1 = 1, nilyr
811 :
812 : ! calculate upper and lower boundary of old cell
813 22891074220 : z1a = ((k1 - 1) * h1) / rnilyr
814 22891074220 : z1b = (k1 * h1) / rnilyr
815 :
816 : ! calculate overlap between old and new cell
817 22891074220 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
818 :
819 : ! aggregate old grid cell contribution to new cell
820 26161227680 : trc2(k2) = trc2(k2) + overlap * trc(k1)
821 :
822 : enddo ! k1
823 :
824 : ! calculate upper and lower boundary of added new ice at bottom
825 3270153460 : z1a = h1
826 3270153460 : z1b = h2
827 :
828 : ! calculate overlap between added ice and new cell
829 3270153460 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
830 : ! aggregate added ice contribution to new cell
831 3270153460 : trc2(k2) = trc2(k2) + overlap * trc0
832 : ! renormalize new grid cell
833 3737318240 : trc2(k2) = (rnilyr * trc2(k2)) / h2
834 :
835 : enddo ! k2
836 :
837 : ! update vertical tracer array with the adjusted tracer
838 3737318240 : trc = trc2
839 :
840 467164780 : end subroutine update_vertical_tracers
841 :
842 : !=======================================================================
843 : !
844 : ! Given the fraction of ice melting laterally in each grid cell
845 : ! (computed in subroutine frzmlt_bottom_lateral), melt ice.
846 : !
847 : ! author: C. M. Bitz, UW
848 : ! 2003: Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
849 : ! 2016 Lettie Roach, NIWA/VUW, added floe size dependence
850 : !
851 722474472 : subroutine lateral_melt (dt, ncat, &
852 : nilyr, nslyr, &
853 : n_aero, &
854 : fpond, &
855 : fresh, fsalt, &
856 722474472 : fhocn, faero_ocn, &
857 722474472 : fiso_ocn, &
858 : rside, meltl, &
859 : fside, sss, &
860 1444948944 : aicen, vicen, &
861 1444948944 : vsnon, trcrn, &
862 722474472 : fzsal, flux_bio, &
863 : nbtrcr, nblyr, &
864 722474472 : nfsd, d_afsd_latm,&
865 722474472 : floe_rad_c, floe_binwidth)
866 :
867 : real (kind=dbl_kind), intent(in) :: &
868 : dt ! time step (s)
869 :
870 : integer (kind=int_kind), intent(in) :: &
871 : ncat , & ! number of thickness categories
872 : nfsd , & ! number of floe size categories
873 : nilyr , & ! number of ice layers
874 : nblyr , & ! number of bio layers
875 : nslyr , & ! number of snow layers
876 : n_aero , & ! number of aerosol tracers
877 : nbtrcr ! number of bio tracers
878 :
879 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
880 : aicen , & ! concentration of ice
881 : vicen , & ! volume per unit area of ice (m)
882 : vsnon ! volume per unit area of snow (m)
883 :
884 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
885 : trcrn ! tracer array
886 :
887 : real (kind=dbl_kind), intent(in) :: &
888 : rside , & ! fraction of ice that melts laterally
889 : fside ! lateral heat flux (W/m^2)
890 :
891 : real (kind=dbl_kind), intent(inout) :: &
892 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
893 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
894 : fsalt , & ! salt flux to ocean (kg/m^2/s)
895 : fhocn , & ! net heat flux to ocean (W/m^2)
896 : meltl , & ! lateral ice melt (m/step-->cm/day)
897 : fzsal ! salt flux from zsalinity (kg/m2/s)
898 :
899 : real (kind=dbl_kind), dimension (:), intent(in) :: &
900 : floe_rad_c , & ! fsd size bin centre in m (radius)
901 : floe_binwidth ! fsd size bin width in m (radius)
902 :
903 : real (kind=dbl_kind), dimension (:), intent(out) :: &
904 : d_afsd_latm ! change in fsd due to lateral melt (m)
905 :
906 : real (kind=dbl_kind), dimension(nbtrcr), &
907 : intent(inout) :: &
908 : flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s)
909 :
910 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
911 : faero_ocn ! aerosol flux to ocean (kg/m^2/s)
912 :
913 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
914 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
915 :
916 : ! local variables
917 :
918 : integer (kind=int_kind) :: &
919 : n , & ! thickness category index
920 : k , & ! layer index
921 : nsubt ! sub timesteps for FSD tendency
922 :
923 : real (kind=dbl_kind) :: &
924 45106800 : dfhocn , & ! change in fhocn
925 45106800 : dfpond , & ! change in fpond
926 45106800 : dfresh , & ! change in fresh
927 45106800 : dfsalt , & ! change in fsalt
928 45106800 : dvssl , & ! snow surface layer volume
929 45106800 : dvint , & ! snow interior layer
930 90213600 : cat1_arealoss, tmp !
931 :
932 : logical (kind=log_kind) :: &
933 : flag ! .true. if there could be lateral melting
934 :
935 : real (kind=dbl_kind), dimension (ncat) :: &
936 1616729664 : aicen_init, & ! initial area fraction
937 1616729664 : vicen_init, & ! volume per unit area of ice (m)
938 1616729664 : G_radialn , & ! rate of lateral melt (m/s)
939 1616729664 : delta_an , & ! change in the ITD
940 1616729664 : qin , & ! enthalpy
941 1616729664 : rsiden ! delta_an/aicen
942 :
943 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
944 1992135984 : afsdn , & ! floe size distribution tracer
945 1992135984 : afsdn_init ! initial value
946 :
947 : real (kind=dbl_kind), dimension (nfsd) :: &
948 1476652704 : df_flx, & ! finite difference for FSD
949 1688783664 : afsd_tmp, d_afsd_tmp
950 :
951 : real (kind=dbl_kind), dimension(nfsd+1) :: &
952 934605432 : f_flx !
953 :
954 : !echmod - for average qin
955 : real (kind=dbl_kind), intent(in) :: &
956 : sss
957 : real (kind=dbl_kind) :: &
958 45106800 : Ti, Si0, qi0, &
959 45106800 : elapsed_t, & ! FSD subcycling
960 45106800 : subdt ! FSD timestep (s)
961 :
962 : character(len=*), parameter :: subname='(lateral_melt)'
963 :
964 722474472 : flag = .false.
965 722474472 : dfhocn = c0
966 722474472 : dfpond = c0
967 722474472 : dfresh = c0
968 722474472 : dfsalt = c0
969 722474472 : dvssl = c0
970 722474472 : dvint = c0
971 722474472 : cat1_arealoss = c0
972 722474472 : tmp = c0
973 4265674992 : vicen_init = c0
974 4265674992 : G_radialn = c0
975 4265674992 : delta_an = c0
976 4265674992 : qin = c0
977 4265674992 : rsiden = c0
978 8992482552 : afsdn = c1
979 8992482552 : afsdn_init = c0
980 1681670352 : df_flx = c0
981 2404144824 : f_flx = c0
982 :
983 722474472 : if (tr_fsd) then
984 22480848 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
985 22480848 : if (icepack_warnings_aborted(subname)) return
986 :
987 1430896368 : afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
988 134885088 : aicen_init = aicen
989 1430896368 : afsdn_init = afsdn ! for diagnostics
990 281683104 : d_afsd_latm(:) = c0
991 : end if
992 :
993 722474472 : if (tr_fsd .and. fside < c0) then
994 2746914 : flag = .true.
995 :
996 :
997 : ! echmod - using category values would be preferable to the average value
998 : ! Compute average enthalpy of ice (taken from add_new_ice)
999 2746914 : if (sss > c2 * dSin0_frazil) then
1000 2746914 : Si0 = sss - dSin0_frazil
1001 : else
1002 0 : Si0 = sss**2 / (c4*dSin0_frazil)
1003 : endif
1004 2746914 : Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
1005 2746914 : qi0 = enthalpy_mush(Ti, Si0)
1006 :
1007 19228398 : do n = 1, ncat
1008 13734570 : if (ktherm == 2) then ! mushy
1009 109876560 : do k = 1, nilyr
1010 : !qin(n) = qin(n) &
1011 : ! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind)
1012 109876560 : qin(n) = qi0
1013 : enddo
1014 : else
1015 0 : qin(n) = -rhoi*Lfresh
1016 : endif
1017 :
1018 13734570 : if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative
1019 :
1020 16481484 : if (G_radialn(n) < -puny) then
1021 :
1022 :
1023 171416610 : if (any(afsdn(:,n) < c0)) print*,&
1024 0 : 'lateral_melt B afsd < 0',n
1025 :
1026 5064160 : cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
1027 16258420 : * G_radialn(n) / floe_binwidth(1)
1028 :
1029 13726340 : delta_an(n) = c0
1030 171416610 : do k = 1, nfsd
1031 23359150 : delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
1032 171416610 : * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
1033 : end do
1034 :
1035 : ! add negative area loss from fsd
1036 13726340 : delta_an(n) = delta_an(n) - cat1_arealoss
1037 :
1038 13726340 : if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
1039 :
1040 : ! following original code, not necessary for fsd
1041 13726340 : if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
1042 :
1043 13726340 : if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)
1044 :
1045 : end if ! G_radialn
1046 : enddo ! ncat
1047 :
1048 719727558 : else if (rside > c0) then ! original, non-fsd implementation
1049 :
1050 63858080 : flag = .true.
1051 381178080 : rsiden(:) = rside ! initialize
1052 :
1053 : endif
1054 :
1055 722474472 : if (flag) then ! grid cells with lateral melting.
1056 :
1057 397659564 : do n = 1, ncat
1058 :
1059 : !-----------------------------------------------------------------
1060 : ! Melt the ice and increment fluxes.
1061 : !-----------------------------------------------------------------
1062 :
1063 : ! fluxes to coupler
1064 : ! dfresh > 0, dfsalt > 0, dfpond > 0
1065 :
1066 331054570 : dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt
1067 331054570 : dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
1068 331054570 : fresh = fresh + dfresh
1069 331054570 : fsalt = fsalt + dfsalt
1070 :
1071 331054570 : if (tr_pond_topo) then
1072 17929344 : dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
1073 17929344 : fpond = fpond - dfpond
1074 : endif
1075 :
1076 : ! history diagnostics
1077 331054570 : meltl = meltl + vicen(n)*rsiden(n)
1078 :
1079 : ! state variables
1080 331054570 : vicen_init(n) = vicen(n)
1081 331054570 : aicen(n) = aicen(n) * (c1 - rsiden(n))
1082 331054570 : vicen(n) = vicen(n) * (c1 - rsiden(n))
1083 331054570 : vsnon(n) = vsnon(n) * (c1 - rsiden(n))
1084 :
1085 : ! floe size distribution
1086 331054570 : if (tr_fsd) then
1087 13734570 : if (rsiden(n) > puny) then
1088 13220695 : if (aicen(n) > puny) then
1089 :
1090 : ! adaptive subtimestep
1091 13220323 : elapsed_t = c0
1092 165073305 : afsd_tmp(:) = afsdn_init(:,n)
1093 165073305 : d_afsd_tmp(:) = c0
1094 13220323 : nsubt = 0
1095 :
1096 26440646 : DO WHILE (elapsed_t.lt.dt)
1097 :
1098 13220323 : nsubt = nsubt + 1
1099 13220323 : if (nsubt.gt.100) &
1100 0 : print *, 'latm not converging'
1101 :
1102 : ! finite differences
1103 165073305 : df_flx(:) = c0
1104 178293628 : f_flx (:) = c0
1105 151852982 : do k = 2, nfsd
1106 151852982 : f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
1107 : end do
1108 :
1109 165073305 : do k = 1, nfsd
1110 165073305 : df_flx(k) = f_flx(k+1) - f_flx(k)
1111 : end do
1112 :
1113 165073305 : if (abs(sum(df_flx(:))) > puny) &
1114 0 : print*,'sum(df_flx)/=0'
1115 :
1116 : ! this term ensures area conservation
1117 165073305 : tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
1118 :
1119 : ! fsd tendency
1120 165073305 : do k = 1, nfsd
1121 22566626 : d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
1122 165073305 : * (c1/floe_rad_c(k) - tmp)
1123 : end do
1124 :
1125 : ! timestep required for this
1126 13220323 : subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
1127 13220323 : subdt = MIN(subdt, dt)
1128 :
1129 : ! update fsd and elapsed time
1130 165073305 : afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
1131 26440646 : elapsed_t = elapsed_t + subdt
1132 :
1133 :
1134 : END DO
1135 :
1136 165073305 : afsdn(:,n) = afsd_tmp(:)
1137 :
1138 :
1139 : end if ! aicen
1140 : end if ! rside > 0, otherwise do nothing
1141 :
1142 : end if ! tr_fsd
1143 :
1144 : ! fluxes
1145 2648436560 : do k = 1, nilyr
1146 : ! enthalpy tracers do not change (e/v constant)
1147 : ! heat flux to coupler for ice melt (dfhocn < 0)
1148 290280410 : dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
1149 2462522195 : * vicen(n)/real(nilyr,kind=dbl_kind)
1150 2648436560 : fhocn = fhocn + dfhocn
1151 : enddo ! nilyr
1152 :
1153 662109140 : do k = 1, nslyr
1154 : ! heat flux to coupler for snow melt (dfhocn < 0)
1155 41468630 : dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
1156 351788885 : * vsnon(n)/real(nslyr,kind=dbl_kind)
1157 662109140 : fhocn = fhocn + dfhocn
1158 : enddo ! nslyr
1159 :
1160 331054570 : if (tr_aero) then
1161 64834848 : do k = 1, n_aero
1162 2384603 : faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
1163 4769206 : *(trcrn(nt_aero +4*(k-1),n) &
1164 4769206 : + trcrn(nt_aero+1+4*(k-1),n)) &
1165 2384603 : + vicen(n) &
1166 4769206 : *(trcrn(nt_aero+2+4*(k-1),n) &
1167 4769206 : + trcrn(nt_aero+3+4*(k-1),n))) &
1168 74373260 : * rsiden(n) / dt
1169 : enddo ! k
1170 : endif ! tr_aero
1171 :
1172 331054570 : if (tr_iso) then
1173 49253300 : do k = 1, n_iso
1174 382845 : fiso_ocn(k) = fiso_ocn(k) &
1175 765690 : + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
1176 765690 : + vicen(n)*trcrn(nt_isoice+k-1,n)) &
1177 50018990 : * rside / dt
1178 : enddo ! k
1179 : endif ! tr_iso
1180 :
1181 : !-----------------------------------------------------------------
1182 : ! Biogeochemistry
1183 : !-----------------------------------------------------------------
1184 :
1185 397659564 : if (z_tracers) then ! snow tracers
1186 7403785 : dvssl = min(p5*vsnon(n), hs_ssl*aicen(n)) !snow surface layer
1187 7403785 : dvint = vsnon(n)- dvssl !snow interior
1188 148075700 : do k = 1, nbtrcr
1189 253175 : flux_bio(k) = flux_bio(k) &
1190 506350 : + (trcrn(bio_index(k)+nblyr+1,n)*dvssl &
1191 506350 : + trcrn(bio_index(k)+nblyr+2,n)*dvint) &
1192 148582050 : * rsiden(n) / dt
1193 : enddo
1194 : endif
1195 :
1196 : enddo ! n
1197 :
1198 66604994 : if (solve_zsal .or. z_tracers) &
1199 : call lateral_melt_bgc(dt, &
1200 : ncat, nblyr, &
1201 0 : rside, vicen_init, & !echmod: use rsiden
1202 0 : trcrn, fzsal, &
1203 1480757 : flux_bio, nbtrcr)
1204 66604994 : if (icepack_warnings_aborted(subname)) return
1205 :
1206 : endif ! flag
1207 :
1208 722474472 : if (tr_fsd) then
1209 :
1210 22480848 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
1211 22480848 : if (icepack_warnings_aborted(subname)) return
1212 :
1213 : ! diagnostics
1214 281683104 : do k = 1, nfsd
1215 259202256 : d_afsd_latm(k) = c0
1216 1577694384 : do n = 1, ncat
1217 177733200 : d_afsd_latm(k) = d_afsd_latm(k) &
1218 1555213536 : + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
1219 : end do
1220 : end do
1221 : end if
1222 :
1223 : end subroutine lateral_melt
1224 :
1225 : !=======================================================================
1226 : !
1227 : ! Given the volume of new ice grown in open water, compute its area
1228 : ! and thickness and add it to the appropriate category or categories.
1229 : !
1230 : ! NOTE: Usually all the new ice is added to category 1. An exception is
1231 : ! made if there is no open water or if the new ice is too thick
1232 : ! for category 1, in which case ice is distributed evenly over the
1233 : ! entire cell. Subroutine rebin should be called in case the ice
1234 : ! thickness lies outside category bounds after new ice formation.
1235 : !
1236 : ! When ice must be added to categories above category 1, the mushy
1237 : ! formulation (ktherm=2) adds it only to the bottom of the ice. When
1238 : ! added to only category 1, all formulations combine the new ice and
1239 : ! existing ice tracers as bulk quantities.
1240 : !
1241 : ! authors: William H. Lipscomb, LANL
1242 : ! Elizabeth C. Hunke, LANL
1243 : ! Adrian Turner, LANL
1244 : ! Lettie Roach, NIWA/VUW
1245 : !
1246 722474472 : subroutine add_new_ice (ncat, nilyr, &
1247 : nfsd, nblyr, &
1248 : n_aero, dt, &
1249 : ntrcr, nltrcr, &
1250 722474472 : hin_max, ktherm, &
1251 722474472 : aicen, trcrn, &
1252 722474472 : vicen, vsnon1, &
1253 : aice0, aice, &
1254 : frzmlt, frazil, &
1255 : frz_onset, yday, &
1256 : update_ocn_f, &
1257 : fresh, fsalt, &
1258 : Tf, sss, &
1259 722474472 : salinz, phi_init, &
1260 : dSin0_frazil, &
1261 722474472 : bgrid, cgrid, igrid, &
1262 722474472 : nbtrcr, flux_bio, &
1263 722474472 : ocean_bio, fzsal, &
1264 : frazil_diag, &
1265 722474472 : fiso_ocn, &
1266 : HDO_ocn, H2_16O_ocn, &
1267 : H2_18O_ocn, &
1268 : wave_sig_ht, &
1269 722474472 : wave_spectrum, &
1270 722474472 : wavefreq, &
1271 722474472 : dwavefreq, &
1272 722474472 : d_afsd_latg, &
1273 722474472 : d_afsd_newi, &
1274 1444948944 : floe_rad_c, floe_binwidth)
1275 :
1276 : use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
1277 :
1278 : integer (kind=int_kind), intent(in) :: &
1279 : ncat , & ! number of thickness categories
1280 : nfsd , & ! number of floe size categories
1281 : nilyr , & ! number of ice layers
1282 : nblyr , & ! number of bio layers
1283 : ntrcr , & ! number of tracers
1284 : nltrcr, & ! number of zbgc tracers
1285 : n_aero, & ! number of aerosol tracers
1286 : ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
1287 :
1288 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
1289 : hin_max ! category boundaries (m)
1290 :
1291 : real (kind=dbl_kind), intent(in) :: &
1292 : dt , & ! time step (s)
1293 : aice , & ! total concentration of ice
1294 : frzmlt, & ! freezing/melting potential (W/m^2)
1295 : Tf , & ! freezing temperature (C)
1296 : sss , & ! sea surface salinity (ppt)
1297 : vsnon1 ! category 1 snow volume per ice area (m)
1298 :
1299 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1300 : aicen , & ! concentration of ice
1301 : vicen ! volume per unit area of ice (m)
1302 :
1303 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1304 : trcrn ! ice tracers
1305 : ! 1: surface temperature
1306 :
1307 : real (kind=dbl_kind), intent(inout) :: &
1308 : aice0 , & ! concentration of open water
1309 : frazil , & ! frazil ice growth (m/step-->cm/day)
1310 : frazil_diag,& ! frazil ice growth diagnostic (m/step-->cm/day)
1311 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
1312 : fsalt ! salt flux to ocean (kg/m^2/s)
1313 :
1314 : real (kind=dbl_kind), intent(inout), optional :: &
1315 : frz_onset ! day of year that freezing begins (congel or frazil)
1316 :
1317 : real (kind=dbl_kind), intent(in), optional :: &
1318 : yday ! day of year
1319 :
1320 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1321 : salinz ! initial salinity profile
1322 :
1323 : real (kind=dbl_kind), intent(in) :: &
1324 : phi_init , & ! initial frazil liquid fraction
1325 : dSin0_frazil ! initial frazil bulk salinity reduction from sss
1326 :
1327 : logical (kind=log_kind), intent(in) :: &
1328 : update_ocn_f ! if true, update fresh water and salt fluxes
1329 :
1330 : ! BGC
1331 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
1332 : bgrid ! biology nondimensional vertical grid points
1333 :
1334 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
1335 : igrid ! biology vertical interface points
1336 :
1337 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
1338 : cgrid ! CICE vertical coordinate
1339 :
1340 : integer (kind=int_kind), intent(in) :: &
1341 : nbtrcr ! number of biology tracers
1342 :
1343 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1344 : flux_bio ! tracer flux to ocean from biology (mmol/m^2/s)
1345 :
1346 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1347 : ocean_bio ! ocean concentration of biological tracer
1348 :
1349 : ! zsalinity
1350 : real (kind=dbl_kind), intent(inout) :: &
1351 : fzsal ! salt flux to ocean from zsalinity (kg/m^2/s)
1352 :
1353 : ! water isotopes
1354 :
1355 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1356 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
1357 :
1358 : real (kind=dbl_kind), intent(in) :: &
1359 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
1360 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
1361 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
1362 :
1363 : ! floe size distribution
1364 : real (kind=dbl_kind), intent(in) :: &
1365 : wave_sig_ht ! significant height of waves globally (m)
1366 :
1367 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1368 : wave_spectrum ! ocean surface wave spectrum, E(f) (m^2 s)
1369 :
1370 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1371 : wavefreq, & ! wave frequencies (s^-1)
1372 : dwavefreq ! wave frequency bin widths (s^-1)
1373 :
1374 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1375 : floe_rad_c , & ! fsd size bin centre in m (radius)
1376 : floe_binwidth ! fsd size bin width in m (radius)
1377 :
1378 : real (kind=dbl_kind), dimension(ncat) :: & ! for now
1379 : ! change in thickness distribution (area)
1380 1616729664 : d_an_latg , & ! due to fsd lateral growth
1381 1616729664 : d_an_newi ! new ice formation
1382 :
1383 : real (kind=dbl_kind), dimension(:), intent(out) :: &
1384 : ! change in thickness distribution (area)
1385 : d_afsd_latg , & ! due to fsd lateral growth
1386 : d_afsd_newi ! new ice formation
1387 :
1388 : ! local variables
1389 :
1390 : integer (kind=int_kind) :: &
1391 : ncats , & ! max categories to which new ice is added, initially
1392 : n , & ! ice category index
1393 : k , & ! ice layer index
1394 : it ! aerosol tracer index
1395 :
1396 : real (kind=dbl_kind) :: &
1397 45106800 : ai0new , & ! area of new ice added to cat 1
1398 45106800 : vi0new , & ! volume of new ice added to cat 1
1399 : vi0new_lat , & ! volume of new ice added laterally to fsd
1400 45106800 : hsurp , & ! thickness of new ice added to each cat
1401 45106800 : fnew , & ! heat flx to open water for new ice (W/m^2)
1402 45106800 : hi0new , & ! thickness of new ice
1403 45106800 : hi0max , & ! max allowed thickness of new ice
1404 45106800 : vsurp , & ! volume of new ice added to each cat
1405 45106800 : vtmp , & ! total volume of new and old ice
1406 45106800 : area1 , & ! starting fractional area of existing ice
1407 45106800 : alvl , & ! starting level ice area
1408 45106800 : dfresh , & ! change in fresh
1409 45106800 : dfsalt , & ! change in fsalt
1410 45106800 : vi0tmp , & ! frzmlt part of frazil
1411 45106800 : Ti , & ! frazil temperature
1412 45106800 : qi0new , & ! frazil ice enthalpy
1413 45106800 : Si0new , & ! frazil ice bulk salinity
1414 45106800 : vi0_init , & ! volume of new ice
1415 45106800 : vice1 , & ! starting volume of existing ice
1416 45106800 : vice_init, vice_final, & ! ice volume summed over categories
1417 45106800 : eice_init, eice_final ! ice energy summed over categories
1418 :
1419 45106800 : real (kind=dbl_kind) :: frazil_conc
1420 :
1421 : real (kind=dbl_kind), dimension (nilyr) :: &
1422 1659381648 : Sprofile ! salinity profile used for new ice additions
1423 :
1424 : character (len=char_len) :: &
1425 : fieldid ! field identifier
1426 :
1427 : real (kind=dbl_kind), dimension (ncat) :: &
1428 1616729664 : eicen, & ! energy of melting for each ice layer (J/m^2)
1429 1616729664 : aicen_init, & ! fractional area of ice
1430 1616729664 : vicen_init ! volume per unit area of ice (m)
1431 :
1432 : ! floe size distribution
1433 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
1434 2082349584 : afsdn ! floe size distribution tracer (originally areal_mfstd_init)
1435 :
1436 : ! real (kind=dbl_kind), dimension (nfsd) :: &
1437 : ! afsd , & ! fsd tracer for each thickness category
1438 :
1439 : real (kind=dbl_kind), dimension (ncat) :: &
1440 1616729664 : d_an_tot, & ! change in the ITD due to lateral growth and new ice
1441 1616729664 : area2 ! area after lateral growth and before new ice formation
1442 :
1443 : real (kind=dbl_kind), dimension (ncat) :: &
1444 984468792 : vin0new ! volume of new ice added to any thickness cat
1445 :
1446 : real (kind=dbl_kind), dimension (nfsd) :: &
1447 : afsd_ni ! areal mFSTD after new ice added
1448 :
1449 : real (kind=dbl_kind) :: &
1450 : tmp, &
1451 45106800 : latsurf_area, & ! fractional area of ice on sides of floes
1452 45106800 : lead_area , & ! fractional area of ice in lead region
1453 45106800 : G_radial , & ! lateral melt rate (m/s)
1454 45106800 : tot_latg , & ! total fsd lateral growth in open water
1455 45106800 : ai0mod ! ai0new - tot_latg
1456 :
1457 : character(len=*),parameter :: subname='(add_new_ice)'
1458 :
1459 : !-----------------------------------------------------------------
1460 : ! initialize
1461 : !-----------------------------------------------------------------
1462 :
1463 722474472 : hsurp = c0
1464 722474472 : hi0new = c0
1465 722474472 : ai0new = c0
1466 8992482552 : afsdn(:,:) = c0
1467 4265674992 : d_an_latg(:) = c0
1468 4265674992 : d_an_tot(:) = c0
1469 4265674992 : d_an_newi(:) = c0
1470 4265674992 : vin0new(:) = c0
1471 :
1472 722474472 : if (tr_fsd) then
1473 281683104 : d_afsd_latg(:) = c0 ! diagnostics
1474 281683104 : d_afsd_newi(:) = c0
1475 : end if
1476 :
1477 4265674992 : area2(:) = aicen(:)
1478 722474472 : lead_area = c0
1479 722474472 : latsurf_area = c0
1480 722474472 : G_radial = c0
1481 722474472 : tot_latg = c0
1482 722474472 : if (ncat > 1) then
1483 699417192 : hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary
1484 : else
1485 23057280 : hi0max = bignum ! big number
1486 : endif
1487 :
1488 722474472 : if (tr_fsd) then
1489 22480848 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
1490 22480848 : if (icepack_warnings_aborted(subname)) return
1491 : endif
1492 :
1493 4265674992 : do n = 1, ncat
1494 3543200520 : aicen_init(n) = aicen(n)
1495 3543200520 : vicen_init(n) = vicen(n)
1496 3543200520 : eicen(n) = c0
1497 4265674992 : if (tr_fsd) then
1498 1408415520 : do k = 1, nfsd
1499 1408415520 : afsdn(k,n) = trcrn(nt_fsd+k-1,n)
1500 : enddo
1501 : endif
1502 : enddo
1503 :
1504 722474472 : if (conserv_check) then
1505 :
1506 161400960 : do n = 1, ncat
1507 1129806720 : do k = 1, nilyr
1508 242101440 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1509 1227800160 : * vicen(n)/real(nilyr,kind=dbl_kind)
1510 : enddo
1511 : enddo
1512 23057280 : call column_sum (ncat, vicen, vice_init)
1513 23057280 : if (icepack_warnings_aborted(subname)) return
1514 23057280 : call column_sum (ncat, eicen, eice_init)
1515 23057280 : if (icepack_warnings_aborted(subname)) return
1516 :
1517 : endif ! conserv_check
1518 :
1519 : !-----------------------------------------------------------------
1520 : ! Compute average enthalpy of new ice.
1521 : ! Sprofile is the salinity profile used when adding new ice to
1522 : ! all categories, for ktherm/=2, and to category 1 for all ktherm.
1523 : !
1524 : ! NOTE: POP assumes new ice is fresh!
1525 : !-----------------------------------------------------------------
1526 :
1527 722474472 : if (ktherm == 2) then ! mushy
1528 632354952 : if (sss > c2 * dSin0_frazil) then
1529 632354952 : Si0new = sss - dSin0_frazil
1530 : else
1531 0 : Si0new = sss**2 / (c4*dSin0_frazil)
1532 : endif
1533 5058839616 : do k = 1, nilyr
1534 5058839616 : Sprofile(k) = Si0new
1535 : enddo
1536 632354952 : Ti = min(liquidus_temperature_mush(Si0new/phi_init), -p1)
1537 632354952 : qi0new = enthalpy_mush(Ti, Si0new)
1538 : else
1539 440786304 : do k = 1, nilyr
1540 440786304 : Sprofile(k) = salinz(k)
1541 : enddo
1542 90119520 : qi0new = -rhoi*Lfresh
1543 : endif ! ktherm
1544 :
1545 : !-----------------------------------------------------------------
1546 : ! Compute the volume, area, and thickness of new ice.
1547 : !-----------------------------------------------------------------
1548 :
1549 722474472 : fnew = max (frzmlt, c0) ! fnew > 0 iff frzmlt > 0
1550 722474472 : vi0new = -fnew*dt / qi0new ! note sign convention, qi < 0
1551 722474472 : vi0_init = vi0new ! for bgc
1552 :
1553 : ! increment ice volume and energy
1554 722474472 : if (conserv_check) then
1555 23057280 : vice_init = vice_init + vi0new
1556 23057280 : eice_init = eice_init + vi0new*qi0new
1557 : endif
1558 :
1559 : ! history diagnostics
1560 722474472 : frazil = vi0new
1561 722474472 : if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss
1562 :
1563 722474472 : if (present(frz_onset) .and. present(yday)) then
1564 722474472 : if (frazil > puny .and. frz_onset < puny) frz_onset = yday
1565 : endif
1566 :
1567 : !-----------------------------------------------------------------
1568 : ! Update fresh water and salt fluxes.
1569 : !
1570 : ! NOTE: POP assumes fresh water and salt flux due to frzmlt > 0
1571 : ! is NOT included in fluxes fresh and fsalt.
1572 : !-----------------------------------------------------------------
1573 :
1574 722474472 : if (update_ocn_f) then
1575 0 : if (ktherm <= 1) then
1576 0 : dfresh = -rhoi*vi0new/dt
1577 0 : dfsalt = ice_ref_salinity*p001*dfresh
1578 0 : fresh = fresh + dfresh
1579 0 : fsalt = fsalt + dfsalt
1580 : ! elseif (ktherm == 2) the fluxes are added elsewhere
1581 : endif
1582 : else ! update_ocn_f = false
1583 722474472 : if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
1584 632354952 : vi0tmp = fnew*dt / (rhoi*Lfresh)
1585 632354952 : dfresh = -rhoi*(vi0new - vi0tmp)/dt
1586 632354952 : dfsalt = ice_ref_salinity*p001*dfresh
1587 632354952 : fresh = fresh + dfresh
1588 632354952 : fsalt = fsalt + dfsalt
1589 632354952 : frazil_diag = frazil - vi0tmp
1590 : ! elseif ktherm==1 do nothing
1591 : endif
1592 : endif
1593 :
1594 : !-----------------------------------------------------------------
1595 : ! Decide how to distribute the new ice.
1596 : !-----------------------------------------------------------------
1597 :
1598 722474472 : if (vi0new > c0) then
1599 :
1600 121907156 : if (tr_fsd) & ! lateral growth of existing ice
1601 : ! calculate change in conc due to lateral growth
1602 : ! update vi0new, without change to afsdn or aicen
1603 : call fsd_lateral_growth (ncat, nfsd, &
1604 : dt, aice, &
1605 0 : aicen, vicen, &
1606 : vi0new, frazil, &
1607 0 : floe_rad_c, afsdn, &
1608 : lead_area, latsurf_area, &
1609 : G_radial, d_an_latg, &
1610 2218855 : tot_latg)
1611 :
1612 121907156 : if (icepack_warnings_aborted(subname)) return
1613 :
1614 121907156 : ai0mod = aice0
1615 : ! separate frazil ice growth from lateral ice growth
1616 121907156 : if (tr_fsd) ai0mod = aice0-tot_latg
1617 :
1618 : ! new ice area and thickness
1619 : ! hin_max(0) < new ice thickness < hin_max(1)
1620 121907156 : if (ai0mod > puny) then
1621 109985153 : hi0new = max(vi0new/ai0mod, hfrazilmin)
1622 109985153 : if (hi0new > hi0max .and. ai0mod+puny < c1) then
1623 : ! distribute excess volume over all categories (below)
1624 44221539 : hi0new = hi0max
1625 44221539 : ai0new = ai0mod
1626 44221539 : vsurp = vi0new - ai0new*hi0new
1627 44221539 : hsurp = vsurp / aice
1628 44221539 : vi0new = ai0new*hi0new
1629 : else
1630 : ! put ice in a single category, with hsurp = 0
1631 65763614 : ai0new = vi0new/hi0new
1632 : endif
1633 : else ! aice0 < puny
1634 11922003 : hsurp = vi0new/aice ! new thickness in each cat
1635 11922003 : vi0new = c0
1636 : endif ! aice0 > puny
1637 :
1638 : ! volume added to each category from lateral growth
1639 717917821 : do n = 1, ncat
1640 717917821 : if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
1641 : end do
1642 :
1643 : ! combine areal change from new ice growth and lateral growth
1644 121907156 : d_an_newi(1) = ai0new
1645 596010665 : d_an_tot(2:ncat) = d_an_latg(2:ncat)
1646 121907156 : d_an_tot(1) = d_an_latg(1) + d_an_newi(1)
1647 121907156 : if (tr_fsd) then
1648 2218855 : vin0new(1) = vin0new(1) + ai0new*hi0new ! not BFB
1649 : else
1650 119688301 : vin0new(1) = vi0new
1651 : endif
1652 :
1653 : endif ! vi0new > puny
1654 :
1655 : !-----------------------------------------------------------------
1656 : ! Distribute excess ice volume among ice categories by increasing
1657 : ! ice thickness, leaving ice area unchanged.
1658 : !
1659 : ! NOTE: If new ice contains globally conserved tracers
1660 : ! (e.g., isotopes from seawater), code must be added here.
1661 : !
1662 : ! The mushy formulation (ktherm=2) puts the new ice only at the
1663 : ! bottom of existing ice and adjusts the layers accordingly.
1664 : ! The other formulations distribute the new ice throughout the
1665 : ! existing ice column.
1666 : !-----------------------------------------------------------------
1667 :
1668 722474472 : if (hsurp > c0) then ! add ice to all categories
1669 :
1670 337207324 : do n = 1, ncat
1671 :
1672 281063782 : vsurp = hsurp * aicen(n)
1673 :
1674 : ! update ice age due to freezing (new ice age = dt)
1675 281063782 : vtmp = vicen(n) + vsurp
1676 281063782 : if (tr_iage .and. vtmp > puny) &
1677 9512150 : trcrn(nt_iage,n) = &
1678 224590719 : (trcrn(nt_iage,n)*vicen(n) + dt*vsurp) / vtmp
1679 :
1680 281063782 : if (tr_lvl .and. vicen(n) > puny) then
1681 12143910 : trcrn(nt_vlvl,n) = &
1682 12143910 : (trcrn(nt_vlvl,n)*vicen(n) + &
1683 245644799 : trcrn(nt_alvl,n)*vsurp) / vtmp
1684 : endif
1685 :
1686 281063782 : if (tr_aero .and. vtmp > puny) then
1687 8051318 : do it = 1, n_aero
1688 962602 : trcrn(nt_aero+2+4*(it-1),n) = &
1689 4506960 : trcrn(nt_aero+2+4*(it-1),n)*vicen(n) / vtmp
1690 962602 : trcrn(nt_aero+3+4*(it-1),n) = &
1691 8532619 : trcrn(nt_aero+3+4*(it-1),n)*vicen(n) / vtmp
1692 : enddo
1693 : endif
1694 :
1695 281063782 : frazil_conc = c0
1696 281063782 : if (tr_iso .and. vtmp > puny) then
1697 27364508 : do it=1,n_iso
1698 20523381 : if (it==1) &
1699 6841127 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1700 20523381 : if (it==2) &
1701 6841127 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1702 20523381 : if (it==3) &
1703 6841127 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1704 :
1705 : ! dilution and uptake in the ice
1706 333168 : trcrn(nt_isoice+it-1,n) &
1707 333168 : = (trcrn(nt_isoice+it-1,n)*vicen(n) &
1708 : + frazil_conc*rhoi*vsurp) &
1709 20689965 : / vtmp
1710 :
1711 166584 : fiso_ocn(it) = fiso_ocn(it) &
1712 27364508 : - frazil_conc*rhoi*vsurp/dt
1713 : enddo
1714 : endif
1715 :
1716 : ! update category volumes
1717 281063782 : vicen(n) = vtmp
1718 :
1719 337207324 : if (ktherm == 2) then
1720 233766428 : vsurp = hsurp * aicen(n) ! note - save this above?
1721 233766428 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1722 233766428 : if (vicen(n) > c0) then
1723 : call update_vertical_tracers(nilyr, &
1724 0 : trcrn(nt_qice:nt_qice+nilyr-1,n), &
1725 233582390 : vtmp, vicen(n), qi0new)
1726 233582390 : if (icepack_warnings_aborted(subname)) return
1727 : call update_vertical_tracers(nilyr, &
1728 0 : trcrn(nt_sice:nt_sice+nilyr-1,n), &
1729 233582390 : vtmp, vicen(n), Si0new)
1730 233582390 : if (icepack_warnings_aborted(subname)) return
1731 : endif
1732 : else
1733 118748662 : do k = 1, nilyr
1734 : ! factor of nilyr cancels out
1735 71451308 : vsurp = hsurp * aicen(n) ! note - save this above?
1736 71451308 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1737 118748662 : if (vicen(n) > c0) then
1738 : ! enthalpy
1739 13926862 : trcrn(nt_qice+k-1,n) = &
1740 61842775 : (trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n)
1741 : ! salinity
1742 54879344 : if (.not. solve_zsal) &
1743 13926862 : trcrn(nt_sice+k-1,n) = &
1744 61842775 : (trcrn(nt_sice+k-1,n)*vtmp + Sprofile(k)*vsurp) / vicen(n)
1745 : endif
1746 : enddo ! k
1747 : endif ! ktherm
1748 :
1749 : enddo ! n
1750 :
1751 : endif ! hsurp > 0
1752 :
1753 : !-----------------------------------------------------------------
1754 : ! Combine new ice grown in open water with ice categories.
1755 : ! Using the floe size distribution, ice is added laterally to all
1756 : ! categories; otherwise it is added to category 1.
1757 : ! Assume that vsnon and esnon are unchanged.
1758 : ! The mushy formulation assumes salt from frazil is added uniformly
1759 : ! to category 1, while the others use a salinity profile.
1760 : !-----------------------------------------------------------------
1761 :
1762 722474472 : ncats = 1 ! add new ice to category 1 by default
1763 722474472 : if (tr_fsd) ncats = ncat ! add new ice laterally to all categories
1764 :
1765 :
1766 1534872336 : do n = 1, ncats
1767 :
1768 1534872336 : if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n
1769 :
1770 118141680 : area1 = aicen(n) ! save
1771 118141680 : vice1 = vicen(n) ! save
1772 118141680 : area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi
1773 118141680 : aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth
1774 :
1775 118141680 : aice0 = aice0 - d_an_tot(n)
1776 118141680 : vicen(n) = vicen(n) + vin0new(n)
1777 :
1778 118141680 : trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n)
1779 118141680 : trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0)
1780 :
1781 118141680 : if (tr_FY) then
1782 80171377 : trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n)
1783 80171377 : trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1)
1784 : endif
1785 :
1786 118141680 : if (tr_fsd) & ! evolve the floe size distribution
1787 : ! both new frazil ice and lateral growth
1788 : call fsd_add_new_ice (ncat, n, nfsd, &
1789 : dt, ai0new, &
1790 0 : d_an_latg, d_an_newi, &
1791 0 : floe_rad_c, floe_binwidth, &
1792 0 : G_radial, area2, &
1793 : wave_sig_ht, &
1794 0 : wave_spectrum, &
1795 0 : wavefreq, &
1796 0 : dwavefreq, &
1797 0 : d_afsd_latg, &
1798 0 : d_afsd_newi, &
1799 0 : afsdn, aicen_init, &
1800 10275952 : aicen, trcrn)
1801 :
1802 118141680 : if (icepack_warnings_aborted(subname)) return
1803 :
1804 118141680 : if (vicen(n) > puny) then
1805 118141671 : if (tr_iage) &
1806 86133640 : trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n)
1807 :
1808 118141671 : if (tr_aero) then
1809 7821712 : do it = 1, n_aero
1810 547364 : trcrn(nt_aero+2+4*(it-1),n) = &
1811 4184538 : trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n)
1812 547364 : trcrn(nt_aero+3+4*(it-1),n) = &
1813 8095394 : trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n)
1814 : enddo
1815 : endif
1816 :
1817 118141671 : frazil_conc = c0
1818 118141671 : if (tr_iso) then
1819 7754576 : do it=1,n_iso
1820 5815932 : if (it==1) &
1821 1938644 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1822 5815932 : if (it==2) &
1823 1938644 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1824 5815932 : if (it==3) &
1825 1938644 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1826 :
1827 94980 : trcrn(nt_isoice+it-1,1) &
1828 94980 : = (trcrn(nt_isoice+it-1,1)*vice1) &
1829 5910912 : + frazil_conc*rhoi*vi0new/vicen(1)
1830 :
1831 47490 : fiso_ocn(it) = fiso_ocn(it) &
1832 7754576 : - frazil_conc*rhoi*vi0new/dt
1833 : enddo
1834 : endif
1835 :
1836 118141671 : if (tr_lvl) then
1837 102317912 : alvl = trcrn(nt_alvl,n)
1838 7165388 : trcrn(nt_alvl,n) = &
1839 102317912 : (trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n)
1840 7165388 : trcrn(nt_vlvl,n) = &
1841 102317912 : (trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n)
1842 : endif
1843 :
1844 118141671 : if (tr_pond_cesm .or. tr_pond_topo) then
1845 288583 : trcrn(nt_apnd,n) = &
1846 5370123 : trcrn(nt_apnd,n)*area1/aicen(n)
1847 112771548 : elseif (tr_pond_lvl) then
1848 70268978 : if (trcrn(nt_alvl,n) > puny) then
1849 3926775 : trcrn(nt_apnd,n) = &
1850 70206083 : trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n))
1851 : endif
1852 : endif
1853 : endif
1854 :
1855 713416884 : do k = 1, nilyr
1856 713416884 : if (vicen(n) > c0) then
1857 : ! factor of nilyr cancels out
1858 : ! enthalpy
1859 81771478 : trcrn(nt_qice+k-1,n) = &
1860 636160943 : (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n)
1861 : ! salinity
1862 595275204 : if (.NOT. solve_zsal)&
1863 81771478 : trcrn(nt_sice+k-1,n) = &
1864 636160943 : (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n)
1865 : endif
1866 : enddo
1867 :
1868 : endif ! vi0new > 0
1869 :
1870 : enddo ! ncats
1871 :
1872 722474472 : if (conserv_check) then
1873 :
1874 161400960 : do n = 1, ncat
1875 138343680 : eicen(n) = c0
1876 1129806720 : do k = 1, nilyr
1877 242101440 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1878 1227800160 : * vicen(n)/real(nilyr,kind=dbl_kind)
1879 : enddo
1880 : enddo
1881 23057280 : call column_sum (ncat, vicen, vice_final)
1882 23057280 : if (icepack_warnings_aborted(subname)) return
1883 23057280 : call column_sum (ncat, eicen, eice_final)
1884 23057280 : if (icepack_warnings_aborted(subname)) return
1885 :
1886 23057280 : fieldid = subname//':vice'
1887 : call column_conservation_check (fieldid, &
1888 : vice_init, vice_final, &
1889 23057280 : puny)
1890 23057280 : if (icepack_warnings_aborted(subname)) return
1891 23057280 : fieldid = subname//':eice'
1892 : call column_conservation_check (fieldid, &
1893 : eice_init, eice_final, &
1894 23057280 : puny*Lfresh*rhoi)
1895 23057280 : if (icepack_warnings_aborted(subname)) return
1896 :
1897 : endif ! conserv_check
1898 :
1899 : !-----------------------------------------------------------------
1900 : ! Biogeochemistry
1901 : !-----------------------------------------------------------------
1902 722474472 : if (tr_brine .or. nbtrcr > 0) &
1903 : call add_new_ice_bgc(dt, nblyr, &
1904 : ncat, nilyr, nltrcr, &
1905 : bgrid, cgrid, igrid, &
1906 0 : aicen_init, vicen_init, vi0_init, &
1907 0 : aicen, vicen, vsnon1, &
1908 0 : vi0new, ntrcr, trcrn, &
1909 0 : nbtrcr, sss, ocean_bio,&
1910 97529376 : flux_bio, hsurp)
1911 722474472 : if (icepack_warnings_aborted(subname)) return
1912 :
1913 : end subroutine add_new_ice
1914 :
1915 : !=======================================================================
1916 : !autodocument_start icepack_step_therm2
1917 : ! Driver for thermodynamic changes not needed for coupling:
1918 : ! transport in thickness space, lateral growth and melting.
1919 : !
1920 : ! authors: William H. Lipscomb, LANL
1921 : ! Elizabeth C. Hunke, LANL
1922 :
1923 722474472 : subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
1924 : nilyr, nslyr, &
1925 722474472 : hin_max, nblyr, &
1926 722474472 : aicen, &
1927 722474472 : vicen, vsnon, &
1928 722474472 : aicen_init, vicen_init, &
1929 722474472 : trcrn, &
1930 : aice0, aice, &
1931 722474472 : trcr_depend, &
1932 722474472 : trcr_base, n_trcr_strata, &
1933 722474472 : nt_strata, &
1934 : Tf, sss, &
1935 722474472 : salinz, &
1936 : rside, meltl, &
1937 : fside, &
1938 : frzmlt, frazil, &
1939 : frain, fpond, &
1940 : fresh, fsalt, &
1941 : fhocn, update_ocn_f, &
1942 722474472 : bgrid, cgrid, &
1943 722474472 : igrid, faero_ocn, &
1944 722474472 : first_ice, fzsal, &
1945 722474472 : flux_bio, ocean_bio, &
1946 : frazil_diag, &
1947 : frz_onset, yday, &
1948 722474472 : fiso_ocn, HDO_ocn, &
1949 : H2_16O_ocn, H2_18O_ocn, &
1950 : nfsd, wave_sig_ht, &
1951 722474472 : wave_spectrum, &
1952 722474472 : wavefreq, &
1953 722474472 : dwavefreq, &
1954 722474472 : d_afsd_latg, d_afsd_newi, &
1955 722474472 : d_afsd_latm, d_afsd_weld, &
1956 722474472 : floe_rad_c, floe_binwidth)
1957 :
1958 : integer (kind=int_kind), intent(in) :: &
1959 : ncat , & ! number of thickness categories
1960 : nfsd , & ! number of floe size categories
1961 : nltrcr , & ! number of zbgc tracers
1962 : nblyr , & ! number of bio layers
1963 : nilyr , & ! number of ice layers
1964 : nslyr ! number of snow layers
1965 :
1966 : logical (kind=log_kind), intent(in) :: &
1967 : update_ocn_f ! if true, update fresh water and salt fluxes
1968 :
1969 : real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
1970 : hin_max ! category boundaries (m)
1971 :
1972 : real (kind=dbl_kind), intent(in) :: &
1973 : dt , & ! time step
1974 : Tf , & ! freezing temperature (C)
1975 : sss , & ! sea surface salinity (ppt)
1976 : rside , & ! fraction of ice that melts laterally
1977 : fside , & ! lateral heat flux (W/m^2)
1978 : frzmlt , & ! freezing/melting potential (W/m^2)
1979 : wave_sig_ht ! significant height of waves in ice (m)
1980 :
1981 : real (kind=dbl_kind), dimension(:), intent(in) :: &
1982 : wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s)
1983 :
1984 : real(kind=dbl_kind), dimension(:), intent(in) :: &
1985 : wavefreq, & ! wave frequencies (s^-1)
1986 : dwavefreq ! wave frequency bin widths (s^-1)
1987 :
1988 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1989 : floe_rad_c , & ! fsd size bin centre in m (radius)
1990 : floe_binwidth ! fsd size bin width in m (radius)
1991 :
1992 : integer (kind=int_kind), dimension (:), intent(in) :: &
1993 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
1994 : n_trcr_strata ! number of underlying tracer layers
1995 :
1996 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1997 : trcr_base ! = 0 or 1 depending on tracer dependency
1998 : ! argument 2: (1) aice, (2) vice, (3) vsno
1999 :
2000 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
2001 : nt_strata ! indices of underlying tracer layers
2002 :
2003 : real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
2004 : bgrid ! biology nondimensional vertical grid points
2005 :
2006 : real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
2007 : igrid ! biology vertical interface points
2008 :
2009 : real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
2010 : cgrid ! CICE vertical coordinate
2011 :
2012 : real (kind=dbl_kind), dimension(:), intent(in) :: &
2013 : salinz , & ! initial salinity profile
2014 : ocean_bio ! ocean concentration of biological tracer
2015 :
2016 : real (kind=dbl_kind), intent(inout) :: &
2017 : aice , & ! sea ice concentration
2018 : aice0 , & ! concentration of open water
2019 : frain , & ! rainfall rate (kg/m^2 s)
2020 : fpond , & ! fresh water flux to ponds (kg/m^2/s)
2021 : fresh , & ! fresh water flux to ocean (kg/m^2/s)
2022 : fsalt , & ! salt flux to ocean (kg/m^2/s)
2023 : fhocn , & ! net heat flux to ocean (W/m^2)
2024 : fzsal , & ! salt flux to ocean from zsalinity (kg/m^2/s)
2025 : meltl , & ! lateral ice melt (m/step-->cm/day)
2026 : frazil , & ! frazil ice growth (m/step-->cm/day)
2027 : frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day)
2028 :
2029 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
2030 : aicen_init,& ! initial concentration of ice
2031 : vicen_init,& ! initial volume per unit area of ice (m)
2032 : aicen , & ! concentration of ice
2033 : vicen , & ! volume per unit area of ice (m)
2034 : vsnon , & ! volume per unit area of snow (m)
2035 : faero_ocn, & ! aerosol flux to ocean (kg/m^2/s)
2036 : flux_bio ! all bio fluxes to ocean
2037 :
2038 : real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
2039 : trcrn ! tracers
2040 :
2041 : logical (kind=log_kind), dimension(:), intent(inout) :: &
2042 : first_ice ! true until ice forms
2043 :
2044 : real (kind=dbl_kind), dimension(:), intent(out) :: &
2045 : ! change in floe size distribution (area)
2046 : d_afsd_latg , & ! due to fsd lateral growth
2047 : d_afsd_newi , & ! new ice formation
2048 : d_afsd_latm , & ! lateral melt
2049 : d_afsd_weld ! welding
2050 :
2051 : real (kind=dbl_kind), intent(inout), optional :: &
2052 : frz_onset ! day of year that freezing begins (congel or frazil)
2053 :
2054 : real (kind=dbl_kind), intent(in), optional :: &
2055 : yday ! day of year
2056 :
2057 : ! water isotopes
2058 : real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
2059 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
2060 :
2061 : real (kind=dbl_kind), optional, intent(in) :: &
2062 : HDO_ocn , & ! ocean concentration of HDO (kg/kg)
2063 : H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
2064 : H2_18O_ocn ! ocean concentration of H2_18O (kg/kg)
2065 : !autodocument_end
2066 :
2067 : ! local variables
2068 :
2069 : ! water isotopes
2070 : real (kind=dbl_kind), dimension(:), allocatable :: &
2071 722474472 : l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
2072 :
2073 : real (kind=dbl_kind) :: &
2074 45106800 : l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
2075 45106800 : l_H2_16O_ocn , & ! local ocean concentration of H2_16O (kg/kg)
2076 45106800 : l_H2_18O_ocn ! local ocean concentration of H2_18O (kg/kg)
2077 :
2078 : character(len=*),parameter :: subname='(icepack_step_therm2)'
2079 :
2080 : !-----------------------------------------------------------------
2081 : ! Check optional arguments and set local values
2082 : !-----------------------------------------------------------------
2083 :
2084 722474472 : if (present(fiso_ocn)) then
2085 722474472 : allocate(l_fiso_ocn(size(fiso_ocn)))
2086 2889897888 : l_fiso_ocn(:) = fiso_ocn(:)
2087 : else
2088 0 : allocate(l_fiso_ocn(1))
2089 0 : l_fiso_ocn(:) = c0
2090 : endif
2091 722474472 : l_HDO_ocn = c0
2092 722474472 : l_H2_16O_ocn = c0
2093 722474472 : l_H2_18O_ocn = c0
2094 722474472 : if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
2095 722474472 : if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
2096 722474472 : if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
2097 :
2098 : !-----------------------------------------------------------------
2099 : ! Let rain drain through to the ocean.
2100 : !-----------------------------------------------------------------
2101 :
2102 722474472 : fresh = fresh + frain * aice
2103 :
2104 : !-----------------------------------------------------------------
2105 : ! Given thermodynamic growth rates, transport ice between
2106 : ! thickness categories.
2107 : !-----------------------------------------------------------------
2108 :
2109 : ! call ice_timer_start(timer_catconv) ! category conversions
2110 :
2111 : !-----------------------------------------------------------------
2112 : ! Compute fractional ice area in each grid cell.
2113 : !-----------------------------------------------------------------
2114 :
2115 722474472 : call aggregate_area (ncat, aicen, aice, aice0)
2116 722474472 : if (icepack_warnings_aborted(subname)) return
2117 :
2118 722474472 : if (kitd == 1) then
2119 :
2120 : !-----------------------------------------------------------------
2121 : ! Identify grid cells with ice.
2122 : !-----------------------------------------------------------------
2123 :
2124 644544648 : if (aice > puny) then
2125 :
2126 : call linear_itd (ncat, hin_max, &
2127 : nilyr, nslyr, &
2128 0 : ntrcr, trcr_depend, &
2129 0 : trcr_base, &
2130 0 : n_trcr_strata, &
2131 0 : nt_strata, &
2132 0 : aicen_init, &
2133 0 : vicen_init, &
2134 0 : aicen, &
2135 0 : trcrn, &
2136 0 : vicen, &
2137 0 : vsnon, &
2138 : aice , &
2139 : aice0 , &
2140 144085983 : fpond )
2141 144085983 : if (icepack_warnings_aborted(subname)) return
2142 :
2143 : endif ! aice > puny
2144 :
2145 : endif ! kitd = 1
2146 :
2147 : ! call ice_timer_stop(timer_catconv) ! category conversions
2148 :
2149 : !-----------------------------------------------------------------
2150 : ! Add frazil ice growing in leads.
2151 : !-----------------------------------------------------------------
2152 :
2153 : ! identify ice-ocean cells
2154 :
2155 : call add_new_ice (ncat, nilyr, &
2156 : nfsd, nblyr, &
2157 : n_aero, dt, &
2158 : ntrcr, nltrcr, &
2159 : hin_max, ktherm, &
2160 0 : aicen, trcrn, &
2161 0 : vicen, vsnon(1), &
2162 : aice0, aice, &
2163 : frzmlt, frazil, &
2164 : frz_onset, yday, &
2165 : update_ocn_f, &
2166 : fresh, fsalt, &
2167 : Tf, sss, &
2168 0 : salinz, phi_init, &
2169 : dSin0_frazil, bgrid, &
2170 : cgrid, igrid, &
2171 0 : nbtrcr, flux_bio, &
2172 0 : ocean_bio, fzsal, &
2173 : frazil_diag, l_fiso_ocn, &
2174 : l_HDO_ocn, l_H2_16O_ocn, &
2175 : l_H2_18O_ocn, &
2176 : wave_sig_ht, &
2177 0 : wave_spectrum, &
2178 0 : wavefreq, dwavefreq, &
2179 0 : d_afsd_latg, d_afsd_newi, &
2180 722474472 : floe_rad_c, floe_binwidth)
2181 :
2182 722474472 : if (icepack_warnings_aborted(subname)) return
2183 :
2184 : !-----------------------------------------------------------------
2185 : ! Melt ice laterally.
2186 : !-----------------------------------------------------------------
2187 :
2188 : call lateral_melt (dt, ncat, &
2189 : nilyr, nslyr, &
2190 : n_aero, fpond, &
2191 : fresh, fsalt, &
2192 0 : fhocn, faero_ocn, &
2193 : l_fiso_ocn, &
2194 : rside, meltl, &
2195 : fside, sss, &
2196 0 : aicen, vicen, &
2197 0 : vsnon, trcrn, &
2198 0 : fzsal, flux_bio, &
2199 : nbtrcr, nblyr, &
2200 0 : nfsd, d_afsd_latm, &
2201 722474472 : floe_rad_c,floe_binwidth)
2202 722474472 : if (icepack_warnings_aborted(subname)) return
2203 :
2204 : ! Floe welding during freezing conditions
2205 722474472 : if (tr_fsd) &
2206 : call fsd_weld_thermo (ncat, nfsd, &
2207 : dt, frzmlt, &
2208 0 : aicen, trcrn, &
2209 22480848 : d_afsd_weld)
2210 722474472 : if (icepack_warnings_aborted(subname)) return
2211 :
2212 : !-----------------------------------------------------------------
2213 : ! For the special case of a single category, adjust the area and
2214 : ! volume (assuming that half the volume change decreases the
2215 : ! thickness, and the other half decreases the area).
2216 : !-----------------------------------------------------------------
2217 :
2218 : !echmod: test this
2219 722474472 : if (ncat==1) &
2220 2882160 : call reduce_area (hin_max (0), &
2221 2882160 : aicen (1), vicen (1), &
2222 23057280 : aicen_init(1), vicen_init(1))
2223 722474472 : if (icepack_warnings_aborted(subname)) return
2224 :
2225 : !-----------------------------------------------------------------
2226 : ! ITD cleanup: Rebin thickness categories if necessary, and remove
2227 : ! categories with very small areas.
2228 : !-----------------------------------------------------------------
2229 :
2230 : call cleanup_itd (dt, ntrcr, &
2231 : nilyr, nslyr, &
2232 : ncat, hin_max, &
2233 0 : aicen, trcrn(1:ntrcr,:), &
2234 0 : vicen, vsnon, &
2235 : aice0, aice, &
2236 : n_aero, &
2237 : nbtrcr, nblyr, &
2238 : tr_aero, &
2239 : tr_pond_topo, heat_capacity, &
2240 0 : first_ice, &
2241 0 : trcr_depend, trcr_base, &
2242 0 : n_trcr_strata, nt_strata, &
2243 : fpond, fresh, &
2244 : fsalt, fhocn, &
2245 0 : faero_ocn, l_fiso_ocn, &
2246 722474472 : fzsal, flux_bio)
2247 722474472 : if (icepack_warnings_aborted(subname)) return
2248 :
2249 722474472 : if (present(fiso_ocn)) then
2250 2889897888 : fiso_ocn(:) = l_fiso_ocn(:)
2251 : endif
2252 722474472 : deallocate(l_fiso_ocn)
2253 :
2254 722474472 : end subroutine icepack_step_therm2
2255 :
2256 : !=======================================================================
2257 :
2258 : end module icepack_therm_itd
2259 :
2260 : !=======================================================================
|
