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 865064 : subroutine linear_itd (ncat, hin_max, &
93 : nilyr, nslyr, &
94 865064 : ntrcr, trcr_depend, &
95 865064 : trcr_base, n_trcr_strata,&
96 865064 : nt_strata, &
97 865064 : aicen_init, vicen_init, &
98 1730128 : aicen, trcrn, &
99 865064 : 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 318046 : slope , & ! rate of change of dhice with hice
148 318046 : dh0 , & ! change in ice thickness at h = 0
149 318046 : da0 , & ! area melting from category 1
150 318046 : damax , & ! max allowed reduction in category 1 area
151 318046 : etamin, etamax,& ! left and right limits of integration
152 318046 : x1 , & ! etamax - etamin
153 318046 : x2 , & ! (etamax^2 - etamin^2) / 2
154 318046 : x3 , & ! (etamax^3 - etamin^3) / 3
155 318046 : wk1, wk2 ! temporary variables
156 :
157 : real (kind=dbl_kind), dimension(0:ncat) :: &
158 3320358 : hbnew ! new boundary locations
159 :
160 : real (kind=dbl_kind), dimension(ncat) :: &
161 3002312 : g0 , & ! constant coefficient in g(h)
162 3002312 : g1 , & ! linear coefficient in g(h)
163 3002312 : hL , & ! left end of range over which g(h) > 0
164 3002312 : hR ! right end of range over which g(h) > 0
165 :
166 : real (kind=dbl_kind), dimension(ncat) :: &
167 3002312 : hicen , & ! ice thickness for each cat (m)
168 3002312 : hicen_init , & ! initial ice thickness for each cat (pre-thermo)
169 3002312 : dhicen , & ! thickness change for remapping (m)
170 3638404 : daice , & ! ice area transferred across boundary
171 3002312 : dvice ! ice volume transferred across boundary
172 :
173 : real (kind=dbl_kind), dimension (ncat) :: &
174 3002312 : eicen, & ! energy of melting for each ice layer (J/m^2)
175 3002312 : esnon, & ! energy of melting for each snow layer (J/m^2)
176 3002312 : vbrin, & ! ice volume defined by brine height (m)
177 3002312 : sicen ! Bulk salt in h ice (ppt*m)
178 :
179 : real (kind=dbl_kind) :: &
180 318046 : vice_init, vice_final, & ! ice volume summed over categories
181 318046 : vsno_init, vsno_final, & ! snow volume summed over categories
182 318046 : eice_init, eice_final, & ! ice energy summed over categories
183 318046 : esno_init, esno_final, & ! snow energy summed over categories
184 318046 : sice_init, sice_final, & ! ice bulk salinity summed over categories
185 318046 : 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 1183110 : 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 865064 : hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
208 :
209 5190384 : do n = 1, ncat
210 4325320 : donor(n) = 0
211 4325320 : daice(n) = c0
212 5190384 : dvice(n) = c0
213 : enddo
214 :
215 : !-----------------------------------------------------------------
216 : ! Compute volume and energy sums that linear remapping should
217 : ! conserve.
218 : !-----------------------------------------------------------------
219 :
220 865064 : if (conserv_check) then
221 :
222 289584 : do n = 1, ncat
223 :
224 241320 : eicen(n) = c0
225 241320 : esnon(n) = c0
226 241320 : vbrin(n) = c0
227 241320 : sicen(n) = c0
228 :
229 1930560 : do k = 1, nilyr
230 1226400 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
231 2543760 : * vicen(n)/real(nilyr,kind=dbl_kind)
232 : enddo
233 1447920 : do k = 1, nslyr
234 876000 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
235 1885920 : * vsnon(n)/real(nslyr,kind=dbl_kind)
236 : enddo
237 :
238 241320 : 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 1978824 : do k = 1, nilyr
244 1226400 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
245 2543760 : * vicen(n)/real(nilyr,kind=dbl_kind)
246 : enddo
247 :
248 : enddo ! n
249 :
250 48264 : call column_sum (ncat, vicen, vice_init)
251 48264 : if (icepack_warnings_aborted(subname)) return
252 48264 : call column_sum (ncat, vsnon, vsno_init)
253 48264 : if (icepack_warnings_aborted(subname)) return
254 48264 : call column_sum (ncat, eicen, eice_init)
255 48264 : if (icepack_warnings_aborted(subname)) return
256 48264 : call column_sum (ncat, esnon, esno_init)
257 48264 : if (icepack_warnings_aborted(subname)) return
258 48264 : call column_sum (ncat, sicen, sice_init)
259 48264 : if (icepack_warnings_aborted(subname)) return
260 48264 : call column_sum (ncat, vbrin, vbri_init)
261 48264 : 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 865064 : remap_flag = .true.
275 :
276 : !-----------------------------------------------------------------
277 : ! Compute thickness change in each category.
278 : !-----------------------------------------------------------------
279 :
280 5190384 : do n = 1, ncat
281 :
282 4325320 : if (aicen_init(n) > puny) then
283 3996940 : hicen_init(n) = vicen_init(n) / aicen_init(n)
284 : else
285 328380 : hicen_init(n) = c0
286 : endif ! aicen_init > puny
287 :
288 5190384 : if (aicen (n) > puny) then
289 3996933 : hicen (n) = vicen(n) / aicen(n)
290 3996933 : dhicen(n) = hicen(n) - hicen_init(n)
291 : else
292 328387 : hicen (n) = c0
293 328387 : 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 865064 : hbnew(0) = hin_max(0)
304 :
305 4325320 : do n = 1, ncat-1
306 :
307 3460256 : if (hicen_init(n) > puny .and. &
308 : hicen_init(n+1) > puny) then
309 : ! interpolate between adjacent category growth rates
310 1146380 : slope = (dhicen(n+1) - dhicen(n)) / &
311 4277127 : (hicen_init(n+1) - hicen_init(n))
312 1146380 : hbnew(n) = hin_max(n) + dhicen(n) &
313 3130747 : + slope * (hin_max(n) - hicen_init(n))
314 329509 : elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0
315 102992 : hbnew(n) = hin_max(n) + dhicen(n)
316 226517 : elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0
317 102918 : hbnew(n) = hin_max(n) + dhicen(n+1)
318 : else
319 123599 : 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 3460256 : if (aicen(n) > puny .and. hicen(n) >= hbnew(n)) then
329 2 : remap_flag = .false.
330 :
331 2 : 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 3460254 : 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 3460256 : 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 4325320 : 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 865064 : if (aicen(ncat) > puny) then
402 763201 : hbnew(ncat) = c3*hicen(ncat) - c2*hbnew(ncat-1)
403 : else
404 101863 : hbnew(ncat) = hin_max(ncat)
405 : endif
406 865064 : hbnew(ncat) = max(hbnew(ncat),hin_max(ncat-1))
407 :
408 : !-----------------------------------------------------------------
409 : ! Identify cells where the ITD is to be remapped
410 : !-----------------------------------------------------------------
411 :
412 865064 : 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 318045 : call fit_line(aicen(1), hicen_init(1), &
420 318045 : hbnew(0), hin_max (1), &
421 318045 : g0 (1), g1 (1), &
422 865062 : hL (1), hR (1))
423 865062 : if (icepack_warnings_aborted(subname)) return
424 :
425 : !-----------------------------------------------------------------
426 : ! Find area lost due to melting of thin (category 1) ice
427 : !-----------------------------------------------------------------
428 :
429 865062 : if (aicen(1) > puny) then
430 :
431 763268 : dh0 = dhicen(1)
432 763268 : if (dh0 < c0) then ! remove area from category 1
433 238039 : 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 238039 : etamax = min(dh0,hR(1)) - hL(1)
441 :
442 238039 : if (etamax > c0) then
443 186328 : x1 = etamax
444 186328 : x2 = p5 * etamax*etamax
445 186328 : da0 = g1(1)*x2 + g0(1)*x1 ! ice area removed
446 :
447 : ! constrain new thickness <= hicen_init
448 186328 : damax = aicen(1) * (c1-hicen(1)/hicen_init(1)) ! damax > 0
449 186328 : da0 = min (da0, damax)
450 :
451 : ! remove area, conserving volume
452 186328 : hicen(1) = hicen(1) * aicen(1) / (aicen(1)-da0)
453 186328 : aicen(1) = aicen(1) - da0
454 :
455 186328 : if (tr_pond_topo) &
456 14468 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
457 50232 : * trcrn(nt_hpnd,1))
458 :
459 : endif ! etamax > 0
460 :
461 : else ! dh0 >= 0
462 525229 : 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 5190372 : do n = 1, ncat
472 :
473 1590225 : call fit_line(aicen(n), hicen(n), &
474 1590225 : hbnew(n-1), hbnew(n), &
475 1590225 : g0 (n), g1 (n), &
476 5915535 : hL (n), hR (n))
477 4325310 : if (icepack_warnings_aborted(subname)) return
478 :
479 : !-----------------------------------------------------------------
480 : ! Compute area and volume to be shifted across each boundary.
481 : !-----------------------------------------------------------------
482 :
483 4325310 : donor(n) = 0
484 4325310 : daice(n) = c0
485 6780597 : dvice(n) = c0
486 : enddo
487 :
488 4325310 : do n = 1, ncat-1
489 :
490 3460248 : if (hbnew(n) > hin_max(n)) then ! transfer from n to n+1
491 :
492 : ! left and right integration limits in eta space
493 1957356 : etamin = max(hin_max(n), hL(n)) - hL(n)
494 1957356 : etamax = min(hbnew(n), hR(n)) - hL(n)
495 1957356 : 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 1502892 : etamin = c0
501 1502892 : etamax = min(hin_max(n), hR(n+1)) - hL(n+1)
502 1502892 : donor(n) = n+1
503 :
504 : endif ! hbnew(n) > hin_max(n)
505 :
506 3460248 : if (etamax > etamin) then
507 2800209 : x1 = etamax - etamin
508 2800209 : wk1 = etamin*etamin
509 2800209 : wk2 = etamax*etamax
510 2800209 : x2 = p5 * (wk2 - wk1)
511 2800209 : wk1 = wk1*etamin
512 2800209 : wk2 = wk2*etamax
513 2800209 : x3 = p333 * (wk2 - wk1)
514 2800209 : nd = donor(n)
515 2800209 : daice(n) = g1(nd)*x2 + g0(nd)*x1
516 2800209 : 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 3460248 : nd = donor(n)
522 :
523 3460248 : if (daice(n) < aicen(nd)*puny) then
524 427934 : daice(n) = c0
525 427934 : dvice(n) = c0
526 427934 : donor(n) = 0
527 : endif
528 :
529 3460248 : if (dvice(n) < vicen(nd)*puny) then
530 427934 : daice(n) = c0
531 427934 : dvice(n) = c0
532 427934 : 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 3460248 : if (daice(n) > aicen(nd)*(c1-puny)) then
539 3768 : daice(n) = aicen(nd)
540 3768 : dvice(n) = vicen(nd)
541 : endif
542 :
543 4325310 : if (dvice(n) > vicen(nd)*(c1-puny)) then
544 3768 : daice(n) = aicen(nd)
545 3768 : 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 5190372 : do n = 1, ncat
556 10480922 : do k = nt_qsno, nt_qsno+nslyr-1
557 9615860 : 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 865062 : daice, dvice )
570 865062 : if (icepack_warnings_aborted(subname)) return
571 :
572 : ! maintain qsno negative definiteness
573 5190372 : do n = 1, ncat
574 10480922 : do k = nt_qsno, nt_qsno+nslyr-1
575 9615860 : 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 865062 : if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then
584 :
585 764 : da0 = aicen(1) * (c1 - hicen(1)/hi_min)
586 764 : aicen(1) = aicen(1) - da0
587 764 : hicen(1) = hi_min
588 :
589 764 : if (tr_pond_topo) &
590 0 : fpond = fpond - (da0 * trcrn(nt_apnd,1) &
591 139 : * 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 865064 : call aggregate_area (ncat, aicen, aice, aice0)
601 865064 : if (icepack_warnings_aborted(subname)) return
602 :
603 : !-----------------------------------------------------------------
604 : ! Check volume and energy conservation.
605 : !-----------------------------------------------------------------
606 :
607 865064 : if (conserv_check) then
608 :
609 289584 : do n = 1, ncat
610 :
611 241320 : eicen(n) = c0
612 241320 : esnon(n) = c0
613 241320 : vbrin(n) = c0
614 241320 : sicen(n) = c0
615 :
616 1930560 : do k = 1, nilyr
617 1226400 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
618 2543760 : * vicen(n)/real(nilyr,kind=dbl_kind)
619 : enddo
620 1447920 : do k = 1, nslyr
621 876000 : esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
622 1885920 : * vsnon(n)/real(nslyr,kind=dbl_kind)
623 : enddo
624 :
625 241320 : 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 1978824 : do k = 1, nilyr
631 1226400 : sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
632 2543760 : * vicen(n)/real(nilyr,kind=dbl_kind)
633 : enddo
634 :
635 : enddo ! n
636 :
637 48264 : call column_sum (ncat, vicen, vice_final)
638 48264 : if (icepack_warnings_aborted(subname)) return
639 48264 : call column_sum (ncat, vsnon, vsno_final)
640 48264 : if (icepack_warnings_aborted(subname)) return
641 48264 : call column_sum (ncat, eicen, eice_final)
642 48264 : if (icepack_warnings_aborted(subname)) return
643 48264 : call column_sum (ncat, esnon, esno_final)
644 48264 : if (icepack_warnings_aborted(subname)) return
645 48264 : call column_sum (ncat, sicen, sice_final)
646 48264 : if (icepack_warnings_aborted(subname)) return
647 48264 : call column_sum (ncat, vbrin, vbri_final)
648 48264 : if (icepack_warnings_aborted(subname)) return
649 :
650 48264 : fieldid = subname//':vice'
651 : call column_conservation_check (fieldid, &
652 : vice_init, vice_final, &
653 48264 : puny)
654 48264 : if (icepack_warnings_aborted(subname)) return
655 48264 : fieldid = subname//':vsno'
656 : call column_conservation_check (fieldid, &
657 : vsno_init, vsno_final, &
658 48264 : puny)
659 48264 : if (icepack_warnings_aborted(subname)) return
660 48264 : fieldid = subname//':eice'
661 : call column_conservation_check (fieldid, &
662 : eice_init, eice_final, &
663 48264 : puny*Lfresh*rhoi)
664 48264 : if (icepack_warnings_aborted(subname)) return
665 48264 : fieldid = subname//':esno'
666 : call column_conservation_check (fieldid, &
667 : esno_init, esno_final, &
668 48264 : puny*Lfresh*rhos)
669 48264 : if (icepack_warnings_aborted(subname)) return
670 48264 : fieldid = subname//':sicen'
671 : call column_conservation_check (fieldid, &
672 : sice_init, sice_final, &
673 48264 : puny)
674 48264 : if (icepack_warnings_aborted(subname)) return
675 48264 : fieldid = subname//':vbrin'
676 : call column_conservation_check (fieldid, &
677 : vbri_init, vbri_final, &
678 48264 : puny*c10)
679 48264 : 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 5190372 : 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 1908270 : h13 , & ! hbL + 1/3 * (hbR - hbL)
716 1908270 : h23 , & ! hbL + 2/3 * (hbR - hbL)
717 1908270 : dhr , & ! 1 / (hR - hL)
718 1908270 : 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 5190372 : if (aicen > puny .and. hbR - hbL > puny) then
727 :
728 : ! Initialize hL and hR
729 :
730 4760198 : hL = hbL
731 4760198 : hR = hbR
732 :
733 : ! Change hL or hR if hicen(n) falls outside central third of range
734 :
735 4760198 : h13 = p333 * (c2*hL + hR)
736 4760198 : h23 = p333 * (hL + c2*hR)
737 4760198 : if (hice < h13) then
738 668248 : hR = c3*hice - c2*hL
739 4091950 : elseif (hice > h23) then
740 593672 : hL = c3*hice - c2*hR
741 : endif
742 :
743 : ! Compute coefficients of g(eta) = g0 + g1*eta
744 :
745 4760198 : dhr = c1 / (hR - hL)
746 4760198 : wk1 = c6 * aicen * dhr
747 4760198 : wk2 = (hice - hL) * dhr
748 4760198 : g0 = wk1 * (p666 - wk2)
749 4760198 : g1 = c2*dhr * wk1 * (wk2 - p5)
750 :
751 : else
752 :
753 430174 : g0 = c0
754 430174 : g1 = c0
755 430174 : hL = c0
756 430174 : hR = c0
757 :
758 : endif ! aicen > puny
759 :
760 5190372 : 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 98634 : 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 494100 : 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 49472 : z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom
791 49472 : z2a, z2b, & ! upper, lower boundary of new cell
792 49472 : overlap , & ! overlap between old and new cell
793 49472 : rnilyr
794 :
795 : character(len=*),parameter :: subname='(update_vertical_tracers)'
796 :
797 98634 : rnilyr = real(nilyr,dbl_kind)
798 :
799 : ! loop over new grid cells
800 789072 : do k2 = 1, nilyr
801 :
802 : ! initialize new tracer
803 690438 : trc2(k2) = c0
804 :
805 : ! calculate upper and lower boundary of new cell
806 690438 : z2a = ((k2 - 1) * h2) / rnilyr
807 690438 : z2b = (k2 * h2) / rnilyr
808 :
809 : ! loop over old grid cells
810 5523504 : do k1 = 1, nilyr
811 :
812 : ! calculate upper and lower boundary of old cell
813 4833066 : z1a = ((k1 - 1) * h1) / rnilyr
814 4833066 : z1b = (k1 * h1) / rnilyr
815 :
816 : ! calculate overlap between old and new cell
817 4833066 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
818 :
819 : ! aggregate old grid cell contribution to new cell
820 5523504 : 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 690438 : z1a = h1
826 690438 : z1b = h2
827 :
828 : ! calculate overlap between added ice and new cell
829 690438 : overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
830 : ! aggregate added ice contribution to new cell
831 690438 : trc2(k2) = trc2(k2) + overlap * trc0
832 : ! renormalize new grid cell
833 789072 : trc2(k2) = (rnilyr * trc2(k2)) / h2
834 :
835 : enddo ! k2
836 :
837 : ! update vertical tracer array with the adjusted tracer
838 789072 : trc = trc2
839 :
840 98634 : 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 1038528 : subroutine lateral_melt (dt, ncat, &
852 : nilyr, nslyr, &
853 : n_aero, &
854 : fpond, &
855 : fresh, fsalt, &
856 1038528 : fhocn, faero_ocn, &
857 1038528 : fiso_ocn, &
858 : rside, meltl, &
859 : fside, sss, &
860 2077056 : aicen, vicen, &
861 2077056 : vsnon, trcrn, &
862 1038528 : fzsal, flux_bio, &
863 : nbtrcr, nblyr, &
864 1038528 : nfsd, d_afsd_latm,&
865 1038528 : 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 381996 : dfhocn , & ! change in fhocn
925 381996 : dfpond , & ! change in fpond
926 381996 : dfresh , & ! change in fresh
927 381996 : dfsalt , & ! change in fsalt
928 381996 : dvssl , & ! snow surface layer volume
929 381996 : dvint , & ! snow interior layer
930 763992 : 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 3499920 : aicen_init, & ! initial area fraction
937 3499920 : vicen_init, & ! volume per unit area of ice (m)
938 3499920 : G_radialn , & ! rate of lateral melt (m/s)
939 3499920 : delta_an , & ! change in the ITD
940 3499920 : qin , & ! enthalpy
941 3499920 : rsiden ! delta_an/aicen
942 :
943 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
944 6750180 : afsdn , & ! floe size distribution tracer
945 6750180 : afsdn_init ! initial value
946 :
947 : real (kind=dbl_kind), dimension (nfsd) :: &
948 2366136 : df_flx, & ! finite difference for FSD
949 4183200 : afsd_tmp, d_afsd_tmp
950 :
951 : real (kind=dbl_kind), dimension(nfsd+1) :: &
952 2855592 : f_flx !
953 :
954 : !echmod - for average qin
955 : real (kind=dbl_kind), intent(in) :: &
956 : sss
957 : real (kind=dbl_kind) :: &
958 381996 : Ti, Si0, qi0, &
959 381996 : elapsed_t, & ! FSD subcycling
960 381996 : subdt ! FSD timestep (s)
961 :
962 : character(len=*), parameter :: subname='(lateral_melt)'
963 :
964 1038528 : flag = .false.
965 1038528 : dfhocn = c0
966 1038528 : dfpond = c0
967 1038528 : dfresh = c0
968 1038528 : dfsalt = c0
969 1038528 : dvssl = c0
970 1038528 : dvint = c0
971 1038528 : cat1_arealoss = c0
972 1038528 : tmp = c0
973 5941584 : vicen_init = c0
974 5941584 : G_radialn = c0
975 5941584 : delta_an = c0
976 5941584 : qin = c0
977 5941584 : rsiden = c0
978 13980960 : afsdn = c1
979 13980960 : afsdn_init = c0
980 2704320 : df_flx = c0
981 3742848 : f_flx = c0
982 :
983 1038528 : if (tr_fsd) then
984 83304 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
985 83304 : if (icepack_warnings_aborted(subname)) return
986 :
987 4052664 : afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
988 499824 : aicen_init = aicen
989 4052664 : afsdn_init = afsdn ! for diagnostics
990 793872 : d_afsd_latm(:) = c0
991 : end if
992 :
993 1038528 : if (tr_fsd .and. fside < c0) then
994 31847 : 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 31847 : if (sss > c2 * dSin0_frazil) then
1000 31847 : Si0 = sss - dSin0_frazil
1001 : else
1002 0 : Si0 = sss**2 / (c4*dSin0_frazil)
1003 : endif
1004 31847 : Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
1005 31847 : qi0 = enthalpy_mush(Ti, Si0)
1006 :
1007 222929 : do n = 1, ncat
1008 159235 : if (ktherm == 2) then ! mushy
1009 1273880 : do k = 1, nilyr
1010 : !qin(n) = qin(n) &
1011 : ! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind)
1012 1273880 : qin(n) = qi0
1013 : enddo
1014 : else
1015 0 : qin(n) = -rhoi*Lfresh
1016 : endif
1017 :
1018 159235 : if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative
1019 :
1020 191082 : if (G_radialn(n) < -puny) then
1021 :
1022 :
1023 1385275 : if (any(afsdn(:,n) < c0)) print*,&
1024 0 : 'lateral_melt B afsd < 0',n
1025 :
1026 210850 : cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
1027 259860 : * G_radialn(n) / floe_binwidth(1)
1028 :
1029 154435 : delta_an(n) = c0
1030 1385275 : do k = 1, nfsd
1031 642720 : delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
1032 1385275 : * 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 154435 : delta_an(n) = delta_an(n) - cat1_arealoss
1037 :
1038 154435 : if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
1039 :
1040 : ! following original code, not necessary for fsd
1041 154435 : if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
1042 :
1043 154435 : if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)
1044 :
1045 : end if ! G_radialn
1046 : enddo ! ncat
1047 :
1048 1006681 : else if (rside > c0) then ! original, non-fsd implementation
1049 :
1050 478829 : flag = .true.
1051 2870846 : rsiden(:) = rside ! initialize
1052 :
1053 : endif
1054 :
1055 1038528 : if (flag) then ! grid cells with lateral melting.
1056 :
1057 3061928 : 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 2551252 : dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt
1067 2551252 : dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
1068 2551252 : fresh = fresh + dfresh
1069 2551252 : fsalt = fsalt + dfsalt
1070 :
1071 2551252 : if (tr_pond_topo) then
1072 329895 : dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
1073 329895 : fpond = fpond - dfpond
1074 : endif
1075 :
1076 : ! history diagnostics
1077 2551252 : meltl = meltl + vicen(n)*rsiden(n)
1078 :
1079 : ! state variables
1080 2551252 : vicen_init(n) = vicen(n)
1081 2551252 : aicen(n) = aicen(n) * (c1 - rsiden(n))
1082 2551252 : vicen(n) = vicen(n) * (c1 - rsiden(n))
1083 2551252 : vsnon(n) = vsnon(n) * (c1 - rsiden(n))
1084 :
1085 : ! floe size distribution
1086 2551252 : if (tr_fsd) then
1087 159235 : if (rsiden(n) > puny) then
1088 153793 : if (aicen(n) > puny) then
1089 :
1090 : ! adaptive subtimestep
1091 153783 : elapsed_t = c0
1092 1379109 : afsd_tmp(:) = afsdn_init(:,n)
1093 1379109 : d_afsd_tmp(:) = c0
1094 153783 : nsubt = 0
1095 :
1096 307566 : DO WHILE (elapsed_t.lt.dt)
1097 :
1098 153783 : nsubt = nsubt + 1
1099 153783 : if (nsubt.gt.100) &
1100 0 : print *, 'latm not converging'
1101 :
1102 : ! finite differences
1103 1379109 : df_flx(:) = c0
1104 1532892 : f_flx (:) = c0
1105 1225326 : do k = 2, nfsd
1106 1225326 : f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
1107 : end do
1108 :
1109 1379109 : do k = 1, nfsd
1110 1379109 : df_flx(k) = f_flx(k+1) - f_flx(k)
1111 : end do
1112 :
1113 1379109 : if (abs(sum(df_flx(:))) > puny) &
1114 0 : print*,'sum(df_flx)/=0'
1115 :
1116 : ! this term ensures area conservation
1117 1379109 : tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
1118 :
1119 : ! fsd tendency
1120 1379109 : do k = 1, nfsd
1121 639858 : d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
1122 1379109 : * (c1/floe_rad_c(k) - tmp)
1123 : end do
1124 :
1125 : ! timestep required for this
1126 153783 : subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
1127 153783 : subdt = MIN(subdt, dt)
1128 :
1129 : ! update fsd and elapsed time
1130 1379109 : afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
1131 307566 : elapsed_t = elapsed_t + subdt
1132 :
1133 :
1134 : END DO
1135 :
1136 1379109 : 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 19123634 : do k = 1, nilyr
1146 : ! enthalpy tracers do not change (e/v constant)
1147 : ! heat flux to coupler for ice melt (dfhocn < 0)
1148 12797182 : dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
1149 22970973 : * vicen(n)/real(nilyr,kind=dbl_kind)
1150 19123634 : fhocn = fhocn + dfhocn
1151 : enddo ! nilyr
1152 :
1153 6067704 : do k = 1, nslyr
1154 : ! heat flux to coupler for snow melt (dfhocn < 0)
1155 2672222 : dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
1156 4852563 : * vsnon(n)/real(nslyr,kind=dbl_kind)
1157 6067704 : fhocn = fhocn + dfhocn
1158 : enddo ! nslyr
1159 :
1160 2551252 : if (tr_aero) then
1161 177350 : do k = 1, n_aero
1162 0 : faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
1163 0 : *(trcrn(nt_aero +4*(k-1),n) &
1164 0 : + trcrn(nt_aero+1+4*(k-1),n)) &
1165 0 : + vicen(n) &
1166 0 : *(trcrn(nt_aero+2+4*(k-1),n) &
1167 0 : + trcrn(nt_aero+3+4*(k-1),n))) &
1168 177350 : * rsiden(n) / dt
1169 : enddo ! k
1170 : endif ! tr_aero
1171 :
1172 2551252 : if (tr_iso) then
1173 350660 : do k = 1, n_iso
1174 0 : fiso_ocn(k) = fiso_ocn(k) &
1175 0 : + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
1176 0 : + vicen(n)*trcrn(nt_isoice+k-1,n)) &
1177 350660 : * rside / dt
1178 : enddo ! k
1179 : endif ! tr_iso
1180 :
1181 : !-----------------------------------------------------------------
1182 : ! Biogeochemistry
1183 : !-----------------------------------------------------------------
1184 :
1185 3061928 : if (z_tracers) then ! snow tracers
1186 315035 : dvssl = min(p5*vsnon(n), hs_ssl*aicen(n)) !snow surface layer
1187 315035 : dvint = vsnon(n)- dvssl !snow interior
1188 6561395 : do k = 1, nbtrcr
1189 1647660 : flux_bio(k) = flux_bio(k) &
1190 3295320 : + (trcrn(bio_index(k)+nblyr+1,n)*dvssl &
1191 3295320 : + trcrn(bio_index(k)+nblyr+2,n)*dvint) &
1192 9856715 : * rsiden(n) / dt
1193 : enddo
1194 : endif
1195 :
1196 : enddo ! n
1197 :
1198 510676 : 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 63007 : flux_bio, nbtrcr)
1204 510676 : if (icepack_warnings_aborted(subname)) return
1205 :
1206 : endif ! flag
1207 :
1208 1038528 : if (tr_fsd) then
1209 :
1210 83304 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
1211 83304 : if (icepack_warnings_aborted(subname)) return
1212 :
1213 : ! diagnostics
1214 793872 : do k = 1, nfsd
1215 710568 : d_afsd_latm(k) = c0
1216 4346712 : do n = 1, ncat
1217 1708200 : d_afsd_latm(k) = d_afsd_latm(k) &
1218 4263408 : + 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 1038528 : subroutine add_new_ice (ncat, nilyr, &
1247 : nfsd, nblyr, &
1248 : n_aero, dt, &
1249 : ntrcr, nltrcr, &
1250 1038528 : hin_max, ktherm, &
1251 1038528 : aicen, trcrn, &
1252 1038528 : vicen, vsnon1, &
1253 : aice0, aice, &
1254 : frzmlt, frazil, &
1255 : frz_onset, yday, &
1256 : update_ocn_f, &
1257 : fresh, fsalt, &
1258 : Tf, sss, &
1259 1038528 : salinz, phi_init, &
1260 : dSin0_frazil, &
1261 1038528 : bgrid, cgrid, igrid, &
1262 1038528 : nbtrcr, flux_bio, &
1263 1038528 : ocean_bio, fzsal, &
1264 : frazil_diag, &
1265 1038528 : fiso_ocn, &
1266 : HDO_ocn, H2_16O_ocn, &
1267 : H2_18O_ocn, &
1268 : wave_sig_ht, &
1269 1038528 : wave_spectrum, &
1270 1038528 : wavefreq, &
1271 1038528 : dwavefreq, &
1272 1038528 : d_afsd_latg, &
1273 1038528 : d_afsd_newi, &
1274 2077056 : 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 3499920 : d_an_latg , & ! due to fsd lateral growth
1381 3499920 : 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 381996 : ai0new , & ! area of new ice added to cat 1
1398 381996 : vi0new , & ! volume of new ice added to cat 1
1399 : vi0new_lat , & ! volume of new ice added laterally to fsd
1400 381996 : hsurp , & ! thickness of new ice added to each cat
1401 381996 : fnew , & ! heat flx to open water for new ice (W/m^2)
1402 381996 : hi0new , & ! thickness of new ice
1403 381996 : hi0max , & ! max allowed thickness of new ice
1404 381996 : vsurp , & ! volume of new ice added to each cat
1405 381996 : vtmp , & ! total volume of new and old ice
1406 381996 : area1 , & ! starting fractional area of existing ice
1407 381996 : alvl , & ! starting level ice area
1408 381996 : dfresh , & ! change in fresh
1409 381996 : dfsalt , & ! change in fsalt
1410 381996 : vi0tmp , & ! frzmlt part of frazil
1411 381996 : Ti , & ! frazil temperature
1412 381996 : qi0new , & ! frazil ice enthalpy
1413 381996 : Si0new , & ! frazil ice bulk salinity
1414 381996 : vi0_init , & ! volume of new ice
1415 381996 : vice1 , & ! starting volume of existing ice
1416 381996 : vice_init, vice_final, & ! ice volume summed over categories
1417 381996 : eice_init, eice_final ! ice energy summed over categories
1418 :
1419 381996 : real (kind=dbl_kind) :: frazil_conc
1420 :
1421 : real (kind=dbl_kind), dimension (nilyr) :: &
1422 4053672 : 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 3499920 : eicen, & ! energy of melting for each ice layer (J/m^2)
1429 3499920 : aicen_init, & ! fractional area of ice
1430 3499920 : vicen_init ! volume per unit area of ice (m)
1431 :
1432 : ! floe size distribution
1433 : real (kind=dbl_kind), dimension (nfsd,ncat) :: &
1434 7514172 : 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 3499920 : d_an_tot, & ! change in the ITD due to lateral growth and new ice
1441 3499920 : area2 ! area after lateral growth and before new ice formation
1442 :
1443 : real (kind=dbl_kind), dimension (ncat) :: &
1444 3225384 : 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 381996 : latsurf_area, & ! fractional area of ice on sides of floes
1452 381996 : lead_area , & ! fractional area of ice in lead region
1453 381996 : G_radial , & ! lateral melt rate (m/s)
1454 381996 : tot_latg , & ! total fsd lateral growth in open water
1455 381996 : ai0mod ! ai0new - tot_latg
1456 :
1457 : character(len=*),parameter :: subname='(add_new_ice)'
1458 :
1459 : !-----------------------------------------------------------------
1460 : ! initialize
1461 : !-----------------------------------------------------------------
1462 :
1463 1038528 : hsurp = c0
1464 1038528 : hi0new = c0
1465 1038528 : ai0new = c0
1466 13980960 : afsdn(:,:) = c0
1467 5941584 : d_an_latg(:) = c0
1468 5941584 : d_an_tot(:) = c0
1469 5941584 : d_an_newi(:) = c0
1470 5941584 : vin0new(:) = c0
1471 :
1472 1038528 : if (tr_fsd) then
1473 793872 : d_afsd_latg(:) = c0 ! diagnostics
1474 793872 : d_afsd_newi(:) = c0
1475 : end if
1476 :
1477 5941584 : area2(:) = aicen(:)
1478 1038528 : lead_area = c0
1479 1038528 : latsurf_area = c0
1480 1038528 : G_radial = c0
1481 1038528 : tot_latg = c0
1482 1038528 : if (ncat > 1) then
1483 966132 : hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary
1484 : else
1485 72396 : hi0max = bignum ! big number
1486 : endif
1487 :
1488 1038528 : if (tr_fsd) then
1489 83304 : call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
1490 83304 : if (icepack_warnings_aborted(subname)) return
1491 : endif
1492 :
1493 5941584 : do n = 1, ncat
1494 4903056 : aicen_init(n) = aicen(n)
1495 4903056 : vicen_init(n) = vicen(n)
1496 4903056 : eicen(n) = c0
1497 5941584 : if (tr_fsd) then
1498 3969360 : do k = 1, nfsd
1499 3969360 : afsdn(k,n) = trcrn(nt_fsd+k-1,n)
1500 : enddo
1501 : endif
1502 : enddo
1503 :
1504 1038528 : if (conserv_check) then
1505 :
1506 434376 : do n = 1, ncat
1507 2968236 : do k = 1, nilyr
1508 1839600 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1509 3815640 : * vicen(n)/real(nilyr,kind=dbl_kind)
1510 : enddo
1511 : enddo
1512 72396 : call column_sum (ncat, vicen, vice_init)
1513 72396 : if (icepack_warnings_aborted(subname)) return
1514 72396 : call column_sum (ncat, eicen, eice_init)
1515 72396 : 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 1038528 : if (ktherm == 2) then ! mushy
1528 676548 : if (sss > c2 * dSin0_frazil) then
1529 676548 : Si0new = sss - dSin0_frazil
1530 : else
1531 0 : Si0new = sss**2 / (c4*dSin0_frazil)
1532 : endif
1533 5412384 : do k = 1, nilyr
1534 5412384 : Sprofile(k) = Si0new
1535 : enddo
1536 676548 : Ti = min(liquidus_temperature_mush(Si0new/phi_init), -p1)
1537 676548 : qi0new = enthalpy_mush(Ti, Si0new)
1538 : else
1539 2027088 : do k = 1, nilyr
1540 2027088 : Sprofile(k) = salinz(k)
1541 : enddo
1542 361980 : qi0new = -rhoi*Lfresh
1543 : endif ! ktherm
1544 :
1545 : !-----------------------------------------------------------------
1546 : ! Compute the volume, area, and thickness of new ice.
1547 : !-----------------------------------------------------------------
1548 :
1549 1038528 : fnew = max (frzmlt, c0) ! fnew > 0 iff frzmlt > 0
1550 1038528 : vi0new = -fnew*dt / qi0new ! note sign convention, qi < 0
1551 1038528 : vi0_init = vi0new ! for bgc
1552 :
1553 : ! increment ice volume and energy
1554 1038528 : if (conserv_check) then
1555 72396 : vice_init = vice_init + vi0new
1556 72396 : eice_init = eice_init + vi0new*qi0new
1557 : endif
1558 :
1559 : ! history diagnostics
1560 1038528 : frazil = vi0new
1561 1038528 : if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss
1562 :
1563 1038528 : if (present(frz_onset) .and. present(yday)) then
1564 1038528 : 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 1038528 : if (update_ocn_f) then
1575 72396 : if (ktherm <= 1) then
1576 72396 : dfresh = -rhoi*vi0new/dt
1577 72396 : dfsalt = ice_ref_salinity*p001*dfresh
1578 72396 : fresh = fresh + dfresh
1579 72396 : fsalt = fsalt + dfsalt
1580 : ! elseif (ktherm == 2) the fluxes are added elsewhere
1581 : endif
1582 : else ! update_ocn_f = false
1583 966132 : if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
1584 676548 : vi0tmp = fnew*dt / (rhoi*Lfresh)
1585 676548 : dfresh = -rhoi*(vi0new - vi0tmp)/dt
1586 676548 : dfsalt = ice_ref_salinity*p001*dfresh
1587 676548 : fresh = fresh + dfresh
1588 676548 : fsalt = fsalt + dfsalt
1589 676548 : 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 1038528 : if (vi0new > c0) then
1599 :
1600 499137 : 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 49014 : tot_latg)
1611 :
1612 499137 : if (icepack_warnings_aborted(subname)) return
1613 :
1614 499137 : ai0mod = aice0
1615 : ! separate frazil ice growth from lateral ice growth
1616 499137 : 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 499137 : if (ai0mod > puny) then
1621 499115 : hi0new = max(vi0new/ai0mod, hfrazilmin)
1622 499115 : if (hi0new > hi0max .and. ai0mod+puny < c1) then
1623 : ! distribute excess volume over all categories (below)
1624 9859 : hi0new = hi0max
1625 9859 : ai0new = ai0mod
1626 9859 : vsurp = vi0new - ai0new*hi0new
1627 9859 : hsurp = vsurp / aice
1628 9859 : vi0new = ai0new*hi0new
1629 : else
1630 : ! put ice in a single category, with hsurp = 0
1631 489256 : ai0new = vi0new/hi0new
1632 : endif
1633 : else ! aice0 < puny
1634 22 : hsurp = vi0new/aice ! new thickness in each cat
1635 22 : vi0new = c0
1636 : endif ! aice0 > puny
1637 :
1638 : ! volume added to each category from lateral growth
1639 2707398 : do n = 1, ncat
1640 2707398 : 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 499137 : d_an_newi(1) = ai0new
1645 2208261 : d_an_tot(2:ncat) = d_an_latg(2:ncat)
1646 499137 : d_an_tot(1) = d_an_latg(1) + d_an_newi(1)
1647 499137 : if (tr_fsd) then
1648 49014 : vin0new(1) = vin0new(1) + ai0new*hi0new ! not BFB
1649 : else
1650 450123 : 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 1038528 : if (hsurp > c0) then ! add ice to all categories
1669 :
1670 59286 : do n = 1, ncat
1671 :
1672 49405 : vsurp = hsurp * aicen(n)
1673 :
1674 : ! update ice age due to freezing (new ice age = dt)
1675 49405 : vtmp = vicen(n) + vsurp
1676 49405 : if (tr_iage .and. vtmp > puny) &
1677 0 : trcrn(nt_iage,n) = &
1678 2655 : (trcrn(nt_iage,n)*vicen(n) + dt*vsurp) / vtmp
1679 :
1680 49405 : if (tr_lvl .and. vicen(n) > puny) then
1681 24736 : trcrn(nt_vlvl,n) = &
1682 24736 : (trcrn(nt_vlvl,n)*vicen(n) + &
1683 44007 : trcrn(nt_alvl,n)*vsurp) / vtmp
1684 : endif
1685 :
1686 49405 : if (tr_aero .and. vtmp > puny) then
1687 5310 : do it = 1, n_aero
1688 0 : trcrn(nt_aero+2+4*(it-1),n) = &
1689 2655 : trcrn(nt_aero+2+4*(it-1),n)*vicen(n) / vtmp
1690 0 : trcrn(nt_aero+3+4*(it-1),n) = &
1691 5310 : trcrn(nt_aero+3+4*(it-1),n)*vicen(n) / vtmp
1692 : enddo
1693 : endif
1694 :
1695 49405 : frazil_conc = c0
1696 49405 : if (tr_iso .and. vtmp > puny) then
1697 10620 : do it=1,n_iso
1698 7965 : if (it==1) &
1699 2655 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1700 7965 : if (it==2) &
1701 2655 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1702 7965 : if (it==3) &
1703 2655 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1704 :
1705 : ! dilution and uptake in the ice
1706 0 : trcrn(nt_isoice+it-1,n) &
1707 0 : = (trcrn(nt_isoice+it-1,n)*vicen(n) &
1708 : + frazil_conc*rhoi*vsurp) &
1709 7965 : / vtmp
1710 :
1711 0 : fiso_ocn(it) = fiso_ocn(it) &
1712 10620 : - frazil_conc*rhoi*vsurp/dt
1713 : enddo
1714 : endif
1715 :
1716 : ! update category volumes
1717 49405 : vicen(n) = vtmp
1718 :
1719 59286 : if (ktherm == 2) then
1720 49405 : vsurp = hsurp * aicen(n) ! note - save this above?
1721 49405 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1722 49405 : if (vicen(n) > c0) then
1723 : call update_vertical_tracers(nilyr, &
1724 0 : trcrn(nt_qice:nt_qice+nilyr-1,n), &
1725 49317 : vtmp, vicen(n), qi0new)
1726 49317 : if (icepack_warnings_aborted(subname)) return
1727 : call update_vertical_tracers(nilyr, &
1728 0 : trcrn(nt_sice:nt_sice+nilyr-1,n), &
1729 49317 : vtmp, vicen(n), Si0new)
1730 49317 : if (icepack_warnings_aborted(subname)) return
1731 : endif
1732 : else
1733 0 : do k = 1, nilyr
1734 : ! factor of nilyr cancels out
1735 0 : vsurp = hsurp * aicen(n) ! note - save this above?
1736 0 : vtmp = vicen(n) - vsurp ! vicen is the new volume
1737 0 : if (vicen(n) > c0) then
1738 : ! enthalpy
1739 0 : trcrn(nt_qice+k-1,n) = &
1740 0 : (trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n)
1741 : ! salinity
1742 0 : if (.not. solve_zsal) &
1743 0 : trcrn(nt_sice+k-1,n) = &
1744 0 : (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 1038528 : ncats = 1 ! add new ice to category 1 by default
1763 1038528 : if (tr_fsd) ncats = ncat ! add new ice laterally to all categories
1764 :
1765 :
1766 2410272 : do n = 1, ncats
1767 :
1768 2410272 : if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n
1769 :
1770 694424 : area1 = aicen(n) ! save
1771 694424 : vice1 = vicen(n) ! save
1772 694424 : area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi
1773 694424 : aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth
1774 :
1775 694424 : aice0 = aice0 - d_an_tot(n)
1776 694424 : vicen(n) = vicen(n) + vin0new(n)
1777 :
1778 694424 : trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n)
1779 694424 : trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0)
1780 :
1781 694424 : if (tr_FY) then
1782 28576 : trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n)
1783 28576 : trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1)
1784 : endif
1785 :
1786 694424 : 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 244323 : aicen, trcrn)
1801 :
1802 694424 : if (icepack_warnings_aborted(subname)) return
1803 :
1804 694424 : if (vicen(n) > puny) then
1805 694424 : if (tr_iage) &
1806 28576 : trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n)
1807 :
1808 694424 : if (tr_aero) then
1809 56748 : do it = 1, n_aero
1810 0 : trcrn(nt_aero+2+4*(it-1),n) = &
1811 28374 : trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n)
1812 0 : trcrn(nt_aero+3+4*(it-1),n) = &
1813 56748 : trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n)
1814 : enddo
1815 : endif
1816 :
1817 694424 : frazil_conc = c0
1818 694424 : if (tr_iso) then
1819 114304 : do it=1,n_iso
1820 85728 : if (it==1) &
1821 28576 : frazil_conc = isoice_alpha(c0,'HDO' ,isotope_frac_method)*HDO_ocn
1822 85728 : if (it==2) &
1823 28576 : frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
1824 85728 : if (it==3) &
1825 28576 : frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
1826 :
1827 0 : trcrn(nt_isoice+it-1,1) &
1828 0 : = (trcrn(nt_isoice+it-1,1)*vice1) &
1829 85728 : + frazil_conc*rhoi*vi0new/vicen(1)
1830 :
1831 0 : fiso_ocn(it) = fiso_ocn(it) &
1832 114304 : - frazil_conc*rhoi*vi0new/dt
1833 : enddo
1834 : endif
1835 :
1836 694424 : if (tr_lvl) then
1837 637660 : alvl = trcrn(nt_alvl,n)
1838 291184 : trcrn(nt_alvl,n) = &
1839 637660 : (trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n)
1840 291184 : trcrn(nt_vlvl,n) = &
1841 637660 : (trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n)
1842 : endif
1843 :
1844 694424 : if (tr_pond_cesm .or. tr_pond_topo) then
1845 0 : trcrn(nt_apnd,n) = &
1846 56764 : trcrn(nt_apnd,n)*area1/aicen(n)
1847 637660 : elseif (tr_pond_lvl) then
1848 494073 : if (trcrn(nt_alvl,n) > puny) then
1849 242147 : trcrn(nt_apnd,n) = &
1850 494073 : trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n))
1851 : endif
1852 : endif
1853 : endif
1854 :
1855 4946554 : do k = 1, nilyr
1856 4946554 : if (vicen(n) > c0) then
1857 : ! factor of nilyr cancels out
1858 : ! enthalpy
1859 3649160 : trcrn(nt_qice+k-1,n) = &
1860 6076710 : (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n)
1861 : ! salinity
1862 4252130 : if (.NOT. solve_zsal)&
1863 3649160 : trcrn(nt_sice+k-1,n) = &
1864 6076710 : (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 1038528 : if (conserv_check) then
1873 :
1874 434376 : do n = 1, ncat
1875 361980 : eicen(n) = c0
1876 2968236 : do k = 1, nilyr
1877 1839600 : eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
1878 3815640 : * vicen(n)/real(nilyr,kind=dbl_kind)
1879 : enddo
1880 : enddo
1881 72396 : call column_sum (ncat, vicen, vice_final)
1882 72396 : if (icepack_warnings_aborted(subname)) return
1883 72396 : call column_sum (ncat, eicen, eice_final)
1884 72396 : if (icepack_warnings_aborted(subname)) return
1885 :
1886 72396 : fieldid = subname//':vice'
1887 : call column_conservation_check (fieldid, &
1888 : vice_init, vice_final, &
1889 72396 : puny)
1890 72396 : if (icepack_warnings_aborted(subname)) return
1891 72396 : fieldid = subname//':eice'
1892 : call column_conservation_check (fieldid, &
1893 : eice_init, eice_final, &
1894 72396 : puny*Lfresh*rhoi)
1895 72396 : if (icepack_warnings_aborted(subname)) return
1896 :
1897 : endif ! conserv_check
1898 :
1899 : !-----------------------------------------------------------------
1900 : ! Biogeochemistry
1901 : !-----------------------------------------------------------------
1902 1038528 : 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 66636 : flux_bio, hsurp)
1911 1038528 : 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 1038528 : subroutine icepack_step_therm2 (dt, ncat, nltrcr, &
1924 : nilyr, nslyr, &
1925 1038528 : hin_max, nblyr, &
1926 1038528 : aicen, &
1927 1038528 : vicen, vsnon, &
1928 1038528 : aicen_init, vicen_init, &
1929 1038528 : trcrn, &
1930 : aice0, aice, &
1931 1038528 : trcr_depend, &
1932 1038528 : trcr_base, n_trcr_strata, &
1933 1038528 : nt_strata, &
1934 : Tf, sss, &
1935 1038528 : salinz, &
1936 : rside, meltl, &
1937 : fside, &
1938 : frzmlt, frazil, &
1939 : frain, fpond, &
1940 : fresh, fsalt, &
1941 : fhocn, update_ocn_f, &
1942 1038528 : bgrid, cgrid, &
1943 1038528 : igrid, faero_ocn, &
1944 1038528 : first_ice, fzsal, &
1945 1038528 : flux_bio, ocean_bio, &
1946 : frazil_diag, &
1947 : frz_onset, yday, &
1948 1038528 : fiso_ocn, HDO_ocn, &
1949 : H2_16O_ocn, H2_18O_ocn, &
1950 : nfsd, wave_sig_ht, &
1951 1038528 : wave_spectrum, &
1952 1038528 : wavefreq, &
1953 1038528 : dwavefreq, &
1954 1038528 : d_afsd_latg, d_afsd_newi, &
1955 1038528 : d_afsd_latm, d_afsd_weld, &
1956 1038528 : 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 1038528 : l_fiso_ocn ! local isotope flux to ocean (kg/m^2/s)
2072 :
2073 : real (kind=dbl_kind) :: &
2074 381996 : l_HDO_ocn , & ! local ocean concentration of HDO (kg/kg)
2075 381996 : l_H2_16O_ocn , & ! local ocean concentration of H2_16O (kg/kg)
2076 381996 : 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 1038528 : if (present(fiso_ocn)) then
2085 1038528 : allocate(l_fiso_ocn(size(fiso_ocn)))
2086 4154112 : l_fiso_ocn(:) = fiso_ocn(:)
2087 : else
2088 0 : allocate(l_fiso_ocn(1))
2089 0 : l_fiso_ocn(:) = c0
2090 : endif
2091 1038528 : l_HDO_ocn = c0
2092 1038528 : l_H2_16O_ocn = c0
2093 1038528 : l_H2_18O_ocn = c0
2094 1038528 : if (present(HDO_ocn) ) l_HDO_ocn = HDO_ocn
2095 1038528 : if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
2096 1038528 : 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 1038528 : 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 1038528 : call aggregate_area (ncat, aicen, aice, aice0)
2116 1038528 : if (icepack_warnings_aborted(subname)) return
2117 :
2118 1038528 : if (kitd == 1) then
2119 :
2120 : !-----------------------------------------------------------------
2121 : ! Identify grid cells with ice.
2122 : !-----------------------------------------------------------------
2123 :
2124 893736 : 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 865064 : fpond )
2141 865064 : 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 1038528 : floe_rad_c, floe_binwidth)
2181 :
2182 1038528 : 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 1038528 : floe_rad_c,floe_binwidth)
2202 1038528 : if (icepack_warnings_aborted(subname)) return
2203 :
2204 : ! Floe welding during freezing conditions
2205 1038528 : if (tr_fsd) &
2206 : call fsd_weld_thermo (ncat, nfsd, &
2207 : dt, frzmlt, &
2208 0 : aicen, trcrn, &
2209 83304 : d_afsd_weld)
2210 1038528 : 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 1038528 : if (ncat==1) &
2220 26280 : call reduce_area (hin_max (0), &
2221 26280 : aicen (1), vicen (1), &
2222 72396 : aicen_init(1), vicen_init(1))
2223 1038528 : 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 1038528 : fzsal, flux_bio)
2247 1038528 : if (icepack_warnings_aborted(subname)) return
2248 :
2249 1038528 : if (present(fiso_ocn)) then
2250 4154112 : fiso_ocn(:) = l_fiso_ocn(:)
2251 : endif
2252 1038528 : deallocate(l_fiso_ocn)
2253 :
2254 1038528 : end subroutine icepack_step_therm2
2255 :
2256 : !=======================================================================
2257 :
2258 : end module icepack_therm_itd
2259 :
2260 : !=======================================================================
|
