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