Line data Source code
1 : !=======================================================================
2 :
3 : ! Routines to initialize the ice thickness distribution and
4 : ! utilities to redistribute ice among categories. These routines
5 : ! are not specific to a particular numerical implementation.
6 : !
7 : ! See Bitz, C.M., and W.H. Lipscomb, 1999:
8 : ! An energy-conserving thermodynamic model of sea ice,
9 : ! J. Geophys. Res., 104, 15,669--15,677.
10 : !
11 : ! See Bitz, C.M., M.M. Holland, A.J. Weaver, M. Eby, 2001:
12 : ! Simulating the ice-thickness distribution in a climate model,
13 : ! J. Geophys. Res., 106, 2441--2464.
14 : !
15 : ! authors: C. M. Bitz, UW
16 : ! William H. Lipscomb and Elizabeth C. Hunke, LANL
17 : !
18 : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
19 : !
20 : ! 2004 WHL: Added multiple snow layers, block structure, cleanup_itd
21 : ! 2006 ECH: Added WMO standard ice thickness categories as kcatbound=2
22 : ! Streamlined for efficiency
23 : ! Converted to free source form (F90)
24 : ! 2014 ECH: Converted to column package
25 :
26 : module icepack_itd
27 :
28 : use icepack_kinds
29 : use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny
30 : use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, rhoi
31 : use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall, rhosnew
32 : use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero
33 : use icepack_tracers, only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice
34 : use icepack_tracers, only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, bio_index
35 : use icepack_tracers, only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
36 : use icepack_tracers, only: icepack_compute_tracers
37 : use icepack_parameters, only: skl_bgc, z_tracers, hi_min
38 : use icepack_parameters, only: kcatbound, kitd, saltflux_option, snwgrain, snwredist
39 : use icepack_therm_shared, only: Tmin
40 : use icepack_warnings, only: warnstr, icepack_warnings_add
41 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
42 : use icepack_zbgc_shared,only: zap_small_bgc
43 :
44 : implicit none
45 :
46 : private
47 : public :: aggregate_area, &
48 : shift_ice, & ! LCOV_EXCL_LINE
49 : column_sum, & ! LCOV_EXCL_LINE
50 : column_conservation_check, & ! LCOV_EXCL_LINE
51 : cleanup_itd, & ! LCOV_EXCL_LINE
52 : reduce_area, & ! LCOV_EXCL_LINE
53 : icepack_init_itd, & ! LCOV_EXCL_LINE
54 : icepack_init_itd_hist, & ! LCOV_EXCL_LINE
55 : icepack_aggregate
56 :
57 : !=======================================================================
58 :
59 : contains
60 :
61 : !=======================================================================
62 :
63 : ! Aggregate ice area (but not other state variables) over thickness
64 : ! categories.
65 : !
66 : ! authors: William H. Lipscomb, LANL
67 :
68 4790054376 : subroutine aggregate_area (aicen, aice, aice0)
69 :
70 : real (kind=dbl_kind), dimension(:), intent(in) :: &
71 : aicen ! concentration of ice
72 :
73 : real (kind=dbl_kind), intent(inout) :: &
74 : aice, & ! concentration of ice ! LCOV_EXCL_LINE
75 : aice0 ! concentration of open water
76 :
77 : ! local variables
78 :
79 : integer (kind=int_kind) :: n
80 :
81 : character(len=*),parameter :: subname='(aggregate_area)'
82 :
83 : !-----------------------------------------------------------------
84 : ! Aggregate
85 : !-----------------------------------------------------------------
86 :
87 4790054376 : aice = c0
88 28684096168 : do n = 1, ncat
89 28684096168 : aice = aice + aicen(n)
90 : enddo ! n
91 :
92 : ! open water fraction
93 4790054376 : aice0 = max (c1 - aice, c0)
94 :
95 4790054376 : end subroutine aggregate_area
96 :
97 : !=======================================================================
98 :
99 : ! Rebins thicknesses into defined categories
100 : !
101 : ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
102 :
103 0 : subroutine rebin (trcr_depend, &
104 690510631 : trcr_base, & ! LCOV_EXCL_LINE
105 690510631 : n_trcr_strata, & ! LCOV_EXCL_LINE
106 690510631 : nt_strata, & ! LCOV_EXCL_LINE
107 690510631 : aicen, trcrn, & ! LCOV_EXCL_LINE
108 1381021262 : vicen, vsnon, & ! LCOV_EXCL_LINE
109 690510631 : hin_max, Tf )
110 :
111 : integer (kind=int_kind), dimension (:), intent(in) :: &
112 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
113 : n_trcr_strata ! number of underlying tracer layers
114 :
115 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
116 : trcr_base ! = 0 or 1 depending on tracer dependency
117 : ! argument 2: (1) aice, (2) vice, (3) vsno
118 :
119 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
120 : nt_strata ! indices of underlying tracer layers
121 :
122 : real (kind=dbl_kind), dimension (ncat), intent(inout) :: &
123 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
124 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
125 : vsnon ! volume per unit area of snow (m)
126 :
127 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
128 : trcrn ! ice tracers
129 :
130 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
131 : hin_max ! category limits (m)
132 :
133 : real (kind=dbl_kind), intent(in) :: &
134 : Tf ! freezing temperature
135 :
136 : ! local variables
137 :
138 : integer (kind=int_kind) :: &
139 : n ! category index
140 :
141 : logical (kind=log_kind) :: &
142 : shiftflag ! = .true. if ice must be shifted
143 :
144 : integer (kind=int_kind), dimension (ncat) :: &
145 1381021262 : donor ! donor category index
146 :
147 : real (kind=dbl_kind), dimension (ncat) :: &
148 690510631 : daice , & ! ice area transferred ! LCOV_EXCL_LINE
149 1381021262 : dvice , & ! ice volume transferred ! LCOV_EXCL_LINE
150 690510631 : hicen ! ice thickness for each cat (m)
151 :
152 : character(len=*),parameter :: subname='(rebin)'
153 :
154 : !-----------------------------------------------------------------
155 : ! Initialize
156 : !-----------------------------------------------------------------
157 :
158 4133578164 : do n = 1, ncat
159 3443067533 : donor(n) = 0
160 3443067533 : daice(n) = c0
161 3443067533 : dvice(n) = c0
162 :
163 : !-----------------------------------------------------------------
164 : ! Compute ice thickness.
165 : !-----------------------------------------------------------------
166 4133578164 : if (aicen(n) > puny) then
167 3225665966 : hicen(n) = vicen(n) / aicen(n)
168 : else
169 217401567 : hicen(n) = c0
170 : endif
171 : enddo ! n
172 :
173 : !-----------------------------------------------------------------
174 : ! make sure thickness of cat 1 is at least hin_max(0)
175 : !-----------------------------------------------------------------
176 :
177 690510631 : if (aicen(1) > puny) then
178 688239613 : if (hicen(1) <= hin_max(0) .and. hin_max(0) > c0 ) then
179 2795623 : aicen(1) = vicen(1) / hin_max(0)
180 2795623 : hicen(1) = hin_max(0)
181 : endif
182 : endif
183 :
184 : !-----------------------------------------------------------------
185 : ! If a category thickness is not in bounds, shift the
186 : ! entire area, volume, and energy to the neighboring category
187 : !-----------------------------------------------------------------
188 :
189 : !-----------------------------------------------------------------
190 : ! Move thin categories up
191 : !-----------------------------------------------------------------
192 :
193 3443067533 : do n = 1, ncat-1 ! loop over category boundaries
194 :
195 : !-----------------------------------------------------------------
196 : ! identify thicknesses that are too big
197 : !-----------------------------------------------------------------
198 2752556902 : shiftflag = .false.
199 2752556902 : if (aicen(n) > puny .and. &
200 : hicen(n) > hin_max(n)) then
201 33173 : shiftflag = .true.
202 33173 : donor(n) = n
203 33173 : daice(n) = aicen(n)
204 33173 : dvice(n) = vicen(n)
205 : endif
206 :
207 3443067533 : if (shiftflag) then
208 :
209 : !-----------------------------------------------------------------
210 : ! shift ice between categories
211 : !-----------------------------------------------------------------
212 :
213 : call shift_ice (trcr_depend, &
214 : trcr_base, & ! LCOV_EXCL_LINE
215 : n_trcr_strata, & ! LCOV_EXCL_LINE
216 : nt_strata, & ! LCOV_EXCL_LINE
217 : aicen, trcrn, & ! LCOV_EXCL_LINE
218 : vicen, vsnon, & ! LCOV_EXCL_LINE
219 : hicen, donor, & ! LCOV_EXCL_LINE
220 33173 : daice, dvice, Tf )
221 33173 : if (icepack_warnings_aborted(subname)) return
222 :
223 : !-----------------------------------------------------------------
224 : ! reset shift parameters
225 : !-----------------------------------------------------------------
226 :
227 33173 : donor(n) = 0
228 33173 : daice(n) = c0
229 33173 : dvice(n) = c0
230 :
231 : endif ! shiftflag
232 :
233 : enddo ! n
234 :
235 : !-----------------------------------------------------------------
236 : ! Move thick categories down
237 : !-----------------------------------------------------------------
238 :
239 3443067533 : do n = ncat-1, 1, -1 ! loop over category boundaries
240 :
241 : !-----------------------------------------------------------------
242 : ! identify thicknesses that are too small
243 : !-----------------------------------------------------------------
244 :
245 2752556902 : shiftflag = .false.
246 2752556902 : if (aicen(n+1) > puny .and. &
247 : hicen(n+1) <= hin_max(n)) then
248 31024 : shiftflag = .true.
249 31024 : donor(n) = n+1
250 31024 : daice(n) = aicen(n+1)
251 31024 : dvice(n) = vicen(n+1)
252 : endif
253 :
254 3443067533 : if (shiftflag) then
255 :
256 : !-----------------------------------------------------------------
257 : ! shift ice between categories
258 : !-----------------------------------------------------------------
259 :
260 : call shift_ice (trcr_depend, &
261 : trcr_base, & ! LCOV_EXCL_LINE
262 : n_trcr_strata, & ! LCOV_EXCL_LINE
263 : nt_strata, & ! LCOV_EXCL_LINE
264 : aicen, trcrn, & ! LCOV_EXCL_LINE
265 : vicen, vsnon, & ! LCOV_EXCL_LINE
266 : hicen, donor, & ! LCOV_EXCL_LINE
267 31024 : daice, dvice, Tf )
268 31024 : if (icepack_warnings_aborted(subname)) return
269 :
270 : !-----------------------------------------------------------------
271 : ! reset shift parameters
272 : !-----------------------------------------------------------------
273 :
274 31024 : donor(n) = 0
275 31024 : daice(n) = c0
276 31024 : dvice(n) = c0
277 :
278 : endif ! shiftflag
279 :
280 : enddo ! n
281 :
282 : end subroutine rebin
283 :
284 : !=======================================================================
285 :
286 : ! Reduce area when ice melts for special case of ncat=1
287 : !
288 : ! Use CSM 1.0-like method of reducing ice area
289 : ! when melting occurs: assume only half the ice volume
290 : ! change goes to thickness decrease, the other half
291 : ! to reduction in ice fraction
292 : !
293 : ! authors: C. M. Bitz, UW
294 : ! modified by: Elizabeth C. Hunke, LANL
295 :
296 17292960 : subroutine reduce_area (hin_max, &
297 : aicen, vicen, & ! LCOV_EXCL_LINE
298 : aicen_init,vicen_init)
299 :
300 : real (kind=dbl_kind), intent(in) :: &
301 : hin_max ! lowest category boundary
302 :
303 : real (kind=dbl_kind), intent(inout) :: &
304 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
305 : vicen ! volume per unit area of ice (m)
306 :
307 : real (kind=dbl_kind), intent(in) :: &
308 : aicen_init, & ! old ice area for category 1 (m) ! LCOV_EXCL_LINE
309 : vicen_init ! old ice volume for category 1 (m)
310 :
311 : ! local variables
312 :
313 : real (kind=dbl_kind) :: &
314 : hi0 , & ! initial hi ! LCOV_EXCL_LINE
315 : hi1 , & ! current hi ! LCOV_EXCL_LINE
316 : dhi ! hi1 - hi0
317 :
318 : character(len=*),parameter :: subname='(reduce_area)'
319 :
320 17292960 : hi0 = c0
321 17292960 : if (aicen_init > c0) &
322 3995720 : hi0 = vicen_init / aicen_init
323 :
324 17292960 : hi1 = c0
325 17292960 : if (aicen > c0) &
326 3984632 : hi1 = vicen / aicen
327 :
328 : ! make sure thickness of cat 1 is at least hin_max(0)
329 17292960 : if (hi1 <= hin_max .and. hin_max > c0 ) then
330 0 : aicen = vicen / hin_max
331 0 : hi1 = hin_max
332 : endif
333 :
334 17292960 : if (aicen > c0) then
335 3984632 : dhi = hi1 - hi0
336 3984632 : if (dhi < c0) then
337 1375695 : hi1 = vicen / aicen
338 1375695 : aicen = c2 * vicen / (hi1 + hi0)
339 : endif
340 : endif
341 :
342 17292960 : end subroutine reduce_area
343 :
344 : !=======================================================================
345 :
346 : ! Shift ice across category boundaries, conserving area, volume, and
347 : ! energy.
348 : !
349 : ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
350 :
351 0 : subroutine shift_ice (trcr_depend, &
352 0 : trcr_base, & ! LCOV_EXCL_LINE
353 307353453 : n_trcr_strata, & ! LCOV_EXCL_LINE
354 307353453 : nt_strata, & ! LCOV_EXCL_LINE
355 614706906 : aicen, trcrn, & ! LCOV_EXCL_LINE
356 307353453 : vicen, vsnon, & ! LCOV_EXCL_LINE
357 614706906 : hicen, donor, & ! LCOV_EXCL_LINE
358 307353453 : daice, dvice, Tf )
359 :
360 : integer (kind=int_kind), dimension (:), intent(in) :: &
361 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
362 : n_trcr_strata ! number of underlying tracer layers
363 :
364 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
365 : trcr_base ! = 0 or 1 depending on tracer dependency
366 : ! argument 2: (1) aice, (2) vice, (3) vsno
367 :
368 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
369 : nt_strata ! indices of underlying tracer layers
370 :
371 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
372 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
373 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
374 : vsnon ! volume per unit area of snow (m)
375 :
376 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
377 : trcrn ! ice tracers
378 :
379 : ! NOTE: Third index of donor, daice, dvice should be ncat-1,
380 : ! except that compilers would have trouble when ncat = 1
381 : integer (kind=int_kind), dimension(:), intent(in) :: &
382 : donor ! donor category index
383 :
384 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
385 : daice , & ! ice area transferred across boundary ! LCOV_EXCL_LINE
386 : dvice , & ! ice volume transferred across boundary ! LCOV_EXCL_LINE
387 : hicen ! ice thickness for each cat (m)
388 :
389 : real (kind=dbl_kind), intent(in) :: &
390 : Tf ! freezing temperature
391 :
392 : ! local variables
393 :
394 : integer (kind=int_kind) :: &
395 : n , & ! thickness category index ! LCOV_EXCL_LINE
396 : nr , & ! receiver category ! LCOV_EXCL_LINE
397 : nd , & ! donor category ! LCOV_EXCL_LINE
398 : it , & ! tracer index ! LCOV_EXCL_LINE
399 : ntr , & ! tracer index ! LCOV_EXCL_LINE
400 : itl ! loop index
401 :
402 : real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
403 614706906 : atrcrn ! aicen*trcrn
404 :
405 : real (kind=dbl_kind) :: &
406 : dvsnow , & ! snow volume transferred ! LCOV_EXCL_LINE
407 : datrcr ! aicen*train transferred
408 :
409 : logical (kind=log_kind) :: &
410 : daice_negative , & ! true if daice < -puny ! LCOV_EXCL_LINE
411 : dvice_negative , & ! true if dvice < -puny ! LCOV_EXCL_LINE
412 : daice_greater_aicen, & ! true if daice > aicen ! LCOV_EXCL_LINE
413 : dvice_greater_vicen ! true if dvice > vicen
414 :
415 : real (kind=dbl_kind) :: &
416 : worka, workb
417 :
418 307353453 : real (kind=dbl_kind), dimension(ncat) :: aicen_init
419 614706906 : real (kind=dbl_kind), dimension(ncat) :: vicen_init
420 307353453 : real (kind=dbl_kind), dimension(ncat) :: vsnon_init
421 :
422 : character(len=*),parameter :: subname='(shift_ice)'
423 :
424 : !-----------------------------------------------------------------
425 : ! store initial snow and ice volume
426 : !-----------------------------------------------------------------
427 :
428 1855333400 : aicen_init(:) = aicen(:)
429 1855333400 : vicen_init(:) = vicen(:)
430 1855333400 : vsnon_init(:) = vsnon(:)
431 :
432 : !-----------------------------------------------------------------
433 : ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
434 : !-----------------------------------------------------------------
435 :
436 1855333400 : do n = 1, ncat
437 46181088050 : do it = 1, ntrcr
438 : atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
439 : + trcr_base(it,2) * vicen(n) & ! LCOV_EXCL_LINE
440 44325754650 : + trcr_base(it,3) * vsnon(n))
441 45873734597 : if (n_trcr_strata(it) > 0) then
442 11653722872 : do itl = 1, n_trcr_strata(it)
443 7270525037 : ntr = nt_strata(it,itl)
444 11653722872 : atrcrn(it,n) = atrcrn(it,n) * trcrn(ntr,n)
445 : enddo
446 : endif
447 : enddo ! it
448 : enddo ! n
449 :
450 : !-----------------------------------------------------------------
451 : ! Check for daice or dvice out of range, allowing for roundoff error
452 : !-----------------------------------------------------------------
453 :
454 1547979947 : do n = 1, ncat-1
455 :
456 1240626494 : daice_negative = .false.
457 1240626494 : dvice_negative = .false.
458 1240626494 : daice_greater_aicen = .false.
459 1240626494 : dvice_greater_vicen = .false.
460 :
461 1240626494 : if (donor(n) > 0) then
462 1120354021 : nd = donor(n)
463 :
464 1120354021 : if (daice(n) < c0) then
465 0 : if (daice(n) > -puny*aicen(nd)) then
466 0 : daice(n) = c0 ! shift no ice
467 0 : dvice(n) = c0
468 : else
469 0 : daice_negative = .true.
470 : endif
471 : endif
472 :
473 1120354021 : if (dvice(n) < c0) then
474 0 : if (dvice(n) > -puny*vicen(nd)) then
475 0 : daice(n) = c0 ! shift no ice
476 0 : dvice(n) = c0
477 : else
478 0 : dvice_negative = .true.
479 : endif
480 : endif
481 :
482 1120354021 : if (daice(n) > aicen(nd)*(c1-puny)) then
483 210537 : if (daice(n) < aicen(nd)*(c1+puny)) then
484 210537 : daice(n) = aicen(nd)
485 210537 : dvice(n) = vicen(nd)
486 : else
487 0 : daice_greater_aicen = .true.
488 : endif
489 : endif
490 :
491 1120354021 : if (dvice(n) > vicen(nd)*(c1-puny)) then
492 210537 : if (dvice(n) < vicen(nd)*(c1+puny)) then
493 210537 : daice(n) = aicen(nd)
494 210537 : dvice(n) = vicen(nd)
495 : else
496 0 : dvice_greater_vicen = .true.
497 : endif
498 : endif
499 :
500 : endif ! donor > 0
501 :
502 : !-----------------------------------------------------------------
503 : ! error messages
504 : !-----------------------------------------------------------------
505 :
506 1240626494 : if (daice_negative) then
507 0 : if (donor(n) > 0 .and. &
508 : daice(n) <= -puny*aicen(nd)) then
509 0 : write(warnstr,*) ' '
510 0 : call icepack_warnings_add(warnstr)
511 0 : write(warnstr,*) subname, 'shift_ice: negative daice'
512 0 : call icepack_warnings_add(warnstr)
513 0 : write(warnstr,*) subname, 'boundary, donor cat:', n, nd
514 0 : call icepack_warnings_add(warnstr)
515 0 : write(warnstr,*) subname, 'daice =', daice(n)
516 0 : call icepack_warnings_add(warnstr)
517 0 : write(warnstr,*) subname, 'dvice =', dvice(n)
518 0 : call icepack_warnings_add(warnstr)
519 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
520 0 : call icepack_warnings_add(subname//' shift_ice: negative daice')
521 : endif
522 : endif
523 1240626494 : if (icepack_warnings_aborted(subname)) return
524 :
525 1240626494 : if (dvice_negative) then
526 0 : if (donor(n) > 0 .and. &
527 : dvice(n) <= -puny*vicen(nd)) then
528 0 : write(warnstr,*) ' '
529 0 : call icepack_warnings_add(warnstr)
530 0 : write(warnstr,*) subname, 'shift_ice: negative dvice'
531 0 : call icepack_warnings_add(warnstr)
532 0 : write(warnstr,*) subname, 'boundary, donor cat:', n, nd
533 0 : call icepack_warnings_add(warnstr)
534 0 : write(warnstr,*) subname, 'daice =', daice(n)
535 0 : call icepack_warnings_add(warnstr)
536 0 : write(warnstr,*) subname, 'dvice =', dvice(n)
537 0 : call icepack_warnings_add(warnstr)
538 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
539 0 : call icepack_warnings_add(subname//' shift_ice: negative dvice')
540 : endif
541 : endif
542 1240626494 : if (icepack_warnings_aborted(subname)) return
543 :
544 1240626494 : if (daice_greater_aicen) then
545 0 : if (donor(n) > 0) then
546 0 : nd = donor(n)
547 0 : if (daice(n) >= aicen(nd)*(c1+puny)) then
548 0 : write(warnstr,*) ' '
549 0 : call icepack_warnings_add(warnstr)
550 0 : write(warnstr,*) subname, 'shift_ice: daice > aicen'
551 0 : call icepack_warnings_add(warnstr)
552 0 : write(warnstr,*) subname, 'boundary, donor cat:', n, nd
553 0 : call icepack_warnings_add(warnstr)
554 0 : write(warnstr,*) subname, 'daice =', daice(n)
555 0 : call icepack_warnings_add(warnstr)
556 0 : write(warnstr,*) subname, 'aicen =', aicen(nd)
557 0 : call icepack_warnings_add(warnstr)
558 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
559 0 : call icepack_warnings_add(subname//' shift_ice: daice > aicen')
560 : endif
561 : endif
562 : endif
563 1240626494 : if (icepack_warnings_aborted(subname)) return
564 :
565 1240626494 : if (dvice_greater_vicen) then
566 0 : if (donor(n) > 0) then
567 0 : nd = donor(n)
568 0 : if (dvice(n) >= vicen(nd)*(c1+puny)) then
569 0 : write(warnstr,*) ' '
570 0 : call icepack_warnings_add(warnstr)
571 0 : write(warnstr,*) subname, 'shift_ice: dvice > vicen'
572 0 : call icepack_warnings_add(warnstr)
573 0 : write(warnstr,*) subname, 'boundary, donor cat:', n, nd
574 0 : call icepack_warnings_add(warnstr)
575 0 : write(warnstr,*) subname, 'dvice =', dvice(n)
576 0 : call icepack_warnings_add(warnstr)
577 0 : write(warnstr,*) subname, 'vicen =', vicen(nd)
578 0 : call icepack_warnings_add(warnstr)
579 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
580 0 : call icepack_warnings_add(subname//' shift_ice: dvice > vicen')
581 : endif
582 : endif
583 : endif
584 1547979947 : if (icepack_warnings_aborted(subname)) return
585 : enddo ! boundaries, 1 to ncat-1
586 :
587 : !-----------------------------------------------------------------
588 : ! transfer volume and energy between categories
589 : !-----------------------------------------------------------------
590 :
591 1547979947 : do n = 1, ncat-1
592 :
593 1547979947 : if (daice(n) > c0) then ! daice(n) can be < puny
594 :
595 1059852766 : nd = donor(n)
596 1059852766 : if (nd == n) then
597 378977080 : nr = nd+1
598 : else ! nd = n+1
599 680875686 : nr = n
600 : endif
601 :
602 1059852766 : aicen(nd) = aicen(nd) - daice(n)
603 1059852766 : aicen(nr) = aicen(nr) + daice(n)
604 :
605 1059852766 : vicen(nd) = vicen(nd) - dvice(n)
606 1059852766 : vicen(nr) = vicen(nr) + dvice(n)
607 :
608 1059852766 : worka = daice(n) / aicen_init(nd)
609 1059852766 : dvsnow = vsnon_init(nd) * worka
610 1059852766 : vsnon(nd) = vsnon(nd) - dvsnow
611 1059852766 : vsnon(nr) = vsnon(nr) + dvsnow
612 1059852766 : workb = dvsnow
613 :
614 31107174231 : do it = 1, ntrcr
615 30047321465 : nd = donor(n)
616 30047321465 : if (nd == n) then
617 9914904049 : nr = nd+1
618 : else ! nd = n+1
619 20132417416 : nr = n
620 : endif
621 :
622 : datrcr = trcrn(it,nd)*(trcr_base(it,1) * daice(n) &
623 : + trcr_base(it,2) * dvice(n) & ! LCOV_EXCL_LINE
624 30047321465 : + trcr_base(it,3) * workb)
625 30047321465 : if (n_trcr_strata(it) > 0) then
626 8078772556 : do itl = 1, n_trcr_strata(it)
627 5039702428 : ntr = nt_strata(it,itl)
628 8078772556 : datrcr = datrcr * trcrn(ntr,nd)
629 : enddo
630 : endif
631 :
632 30047321465 : atrcrn(it,nd) = atrcrn(it,nd) - datrcr
633 31107174231 : atrcrn(it,nr) = atrcrn(it,nr) + datrcr
634 :
635 : enddo ! ntrcr
636 : endif ! daice
637 : enddo ! boundaries, 1 to ncat-1
638 :
639 : !-----------------------------------------------------------------
640 : ! Update ice thickness and tracers
641 : !-----------------------------------------------------------------
642 :
643 1855333400 : do n = 1, ncat
644 :
645 1547979947 : if (aicen(n) > puny) then
646 1482437292 : hicen(n) = vicen (n) / aicen(n)
647 : else
648 65542655 : hicen(n) = c0
649 : endif
650 :
651 : !-----------------------------------------------------------------
652 : ! Compute new tracers
653 : !-----------------------------------------------------------------
654 :
655 : call icepack_compute_tracers(trcr_depend, &
656 : atrcrn(:,n), aicen(n), & ! LCOV_EXCL_LINE
657 : vicen(n), vsnon(n), & ! LCOV_EXCL_LINE
658 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
659 1547979947 : nt_strata, trcrn(:,n), Tf)
660 1855333400 : if (icepack_warnings_aborted(subname)) return
661 :
662 : enddo ! ncat
663 :
664 : end subroutine shift_ice
665 :
666 : !=======================================================================
667 :
668 : ! For each grid cell, sum field over all ice categories.
669 : !
670 : ! author: William H. Lipscomb, LANL
671 :
672 415855848 : subroutine column_sum (nsum, xin, xout)
673 :
674 : integer (kind=int_kind), intent(in) :: &
675 : nsum ! number of categories/layers
676 :
677 : real (kind=dbl_kind), dimension (nsum), intent(in) :: &
678 : xin ! input field
679 :
680 : real (kind=dbl_kind), intent(out) :: &
681 : xout ! output field
682 :
683 : ! local variables
684 :
685 : integer (kind=int_kind) :: &
686 : n ! category/layer index
687 :
688 : character(len=*),parameter :: subname='(column_sum)'
689 :
690 415855848 : xout = c0
691 2827103784 : do n = 1, nsum
692 2827103784 : xout = xout + xin(n)
693 : enddo ! n
694 :
695 415855848 : end subroutine column_sum
696 :
697 : !=======================================================================
698 :
699 : ! For each physical grid cell, check that initial and final values
700 : ! of a conserved field are equal to within a small value.
701 : !
702 : ! author: William H. Lipscomb, LANL
703 :
704 165984348 : subroutine column_conservation_check (fieldid, &
705 : x1, x2, & ! LCOV_EXCL_LINE
706 : max_err )
707 :
708 : real (kind=dbl_kind), intent(in) :: &
709 : x1 , & ! initial field ! LCOV_EXCL_LINE
710 : x2 ! final field
711 :
712 : real (kind=dbl_kind), intent(in) :: &
713 : max_err ! max allowed error
714 :
715 : character (len=char_len), intent(in) :: &
716 : fieldid ! field identifier
717 :
718 : character(len=*),parameter :: subname='(column_conservation_check)'
719 :
720 : ! local variables
721 :
722 165984348 : if (abs (x2-x1) > max_err) then
723 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
724 0 : write(warnstr,*) ' '
725 0 : call icepack_warnings_add(warnstr)
726 0 : write(warnstr,*) subname, 'Conservation error: ', trim(fieldid)
727 0 : call icepack_warnings_add(warnstr)
728 0 : write(warnstr,*) subname, ' Initial value = ', x1
729 0 : call icepack_warnings_add(warnstr)
730 0 : write(warnstr,*) subname, ' Final value = ', x2
731 0 : call icepack_warnings_add(warnstr)
732 0 : write(warnstr,*) subname, ' Difference = ', x2 - x1
733 0 : call icepack_warnings_add(warnstr)
734 : endif
735 :
736 165984348 : end subroutine column_conservation_check
737 :
738 : !=======================================================================
739 :
740 : ! Cleanup subroutine that rebins thickness categories if necessary,
741 : ! eliminates very small ice areas while conserving mass and energy
742 : ! and aggregates state variables.
743 : ! It is a good idea to call this subroutine after the thermodynamics
744 : ! (thermo_vertical/thermo_itd) and again after the dynamics
745 : ! (evp/transport/ridging).
746 : !
747 : ! author: William H. Lipscomb, LANL
748 :
749 5981086560 : subroutine cleanup_itd (dt, hin_max, &
750 5981086560 : aicen, trcrn, & ! LCOV_EXCL_LINE
751 2990543280 : vicen, vsnon, & ! LCOV_EXCL_LINE
752 : aice0, aice, & ! LCOV_EXCL_LINE
753 : tr_aero, & ! LCOV_EXCL_LINE
754 : tr_pond_topo, & ! LCOV_EXCL_LINE
755 2990543280 : first_ice, & ! LCOV_EXCL_LINE
756 2990543280 : trcr_depend, trcr_base, & ! LCOV_EXCL_LINE
757 2990543280 : n_trcr_strata,nt_strata, & ! LCOV_EXCL_LINE
758 : fpond, fresh, & ! LCOV_EXCL_LINE
759 : fsalt, fhocn, & ! LCOV_EXCL_LINE
760 2990543280 : faero_ocn, fiso_ocn, & ! LCOV_EXCL_LINE
761 2990543280 : flux_bio, Tf, & ! LCOV_EXCL_LINE
762 : limit_aice, dorebin)
763 :
764 : real (kind=dbl_kind), intent(in) :: &
765 : dt ! time step
766 :
767 : real (kind=dbl_kind), intent(in) :: &
768 : Tf ! Freezing temperature
769 :
770 : real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
771 : hin_max ! category boundaries (m)
772 :
773 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
774 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
775 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
776 : vsnon ! volume per unit area of snow (m)
777 :
778 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
779 : trcrn ! ice tracers
780 :
781 : real (kind=dbl_kind), intent(inout) :: &
782 : aice , & ! total ice concentration ! LCOV_EXCL_LINE
783 : aice0 ! concentration of open water
784 :
785 : integer (kind=int_kind), dimension (:), intent(in) :: &
786 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
787 : n_trcr_strata ! number of underlying tracer layers
788 :
789 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
790 : trcr_base ! = 0 or 1 depending on tracer dependency
791 : ! argument 2: (1) aice, (2) vice, (3) vsno
792 :
793 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
794 : nt_strata ! indices of underlying tracer layers
795 :
796 : logical (kind=log_kind), intent(in) :: &
797 : tr_aero, & ! aerosol flag ! LCOV_EXCL_LINE
798 : tr_pond_topo ! topo pond flag
799 :
800 : logical (kind=log_kind), dimension(ncat), intent(inout) :: &
801 : first_ice ! For bgc and S tracers. set to true if zapping ice.
802 :
803 : ! ice-ocean fluxes (required for strict conservation)
804 :
805 : real (kind=dbl_kind), intent(inout), optional :: &
806 : fpond , & ! fresh water flux to ponds (kg/m^2/s) ! LCOV_EXCL_LINE
807 : fresh , & ! fresh water flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
808 : fsalt , & ! salt flux to ocean (kg/m^2/s) ! LCOV_EXCL_LINE
809 : fhocn ! net heat flux to ocean (W/m^2)
810 :
811 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
812 : flux_bio ! net tracer flux to ocean from biology (mmol/m^2/s)
813 :
814 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
815 : faero_ocn ! aerosol flux to ocean (kg/m^2/s)
816 :
817 : real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
818 : fiso_ocn ! isotope flux to ocean (kg/m^2/s)
819 :
820 : logical (kind=log_kind), intent(in), optional :: &
821 : dorebin, & ! if false, do not call rebin (default true) ! LCOV_EXCL_LINE
822 : limit_aice ! if false, allow aice to be out of bounds
823 : ! may want to allow this for unit tests (default true)
824 :
825 : ! local variables
826 :
827 : integer (kind=int_kind) :: &
828 : n , & ! category index ! LCOV_EXCL_LINE
829 : it ! tracer index
830 :
831 : real (kind=dbl_kind) &
832 : dfpond , & ! zapped pond water flux (kg/m^2/s) ! LCOV_EXCL_LINE
833 : dfresh , & ! zapped fresh water flux (kg/m^2/s) ! LCOV_EXCL_LINE
834 : dfsalt , & ! zapped salt flux (kg/m^2/s) ! LCOV_EXCL_LINE
835 : dfhocn ! zapped energy flux ( W/m^2)
836 :
837 : real (kind=dbl_kind), dimension (n_aero) :: &
838 2990543280 : dfaero_ocn ! zapped aerosol flux (kg/m^2/s)
839 :
840 : real (kind=dbl_kind), dimension (n_iso) :: &
841 5981086560 : dfiso_ocn ! zapped isotope flux (kg/m^2/s)
842 :
843 : real (kind=dbl_kind), dimension (ntrcr) :: &
844 2990543280 : dflux_bio ! zapped biology flux (mmol/m^2/s)
845 :
846 : logical (kind=log_kind) :: &
847 : ldorebin , & ! if true, call rebin ! LCOV_EXCL_LINE
848 : llimit_aice ! if true, check for aice out of bounds
849 :
850 : character(len=*),parameter :: subname='(cleanup_itd)'
851 :
852 : !-----------------------------------------------------------------
853 : ! Initialize
854 : !-----------------------------------------------------------------
855 :
856 2990543280 : if (present(limit_aice)) then
857 0 : llimit_aice = limit_aice
858 : else
859 2990543280 : llimit_aice = .true.
860 : endif
861 :
862 2990543280 : if (present(dorebin)) then
863 1498353240 : ldorebin = dorebin
864 : else
865 1492190040 : ldorebin = .true.
866 : endif
867 :
868 2990543280 : dfpond = c0
869 2990543280 : dfresh = c0
870 2990543280 : dfsalt = c0
871 2990543280 : dfhocn = c0
872 5875487136 : dfaero_ocn(:) = c0
873 3078160944 : dfiso_ocn (:) = c0
874 97478850864 : dflux_bio (:) = c0
875 :
876 : !-----------------------------------------------------------------
877 : ! Compute total ice area.
878 : !-----------------------------------------------------------------
879 :
880 2990543280 : call aggregate_area (aicen, aice, aice0)
881 2990543280 : if (icepack_warnings_aborted(subname)) return
882 :
883 2990543280 : if (llimit_aice) then ! check for aice out of bounds
884 2990543280 : if (aice > c1+puny .or. aice < -puny) then
885 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
886 0 : call icepack_warnings_add(subname//' aggregate ice area out of bounds')
887 0 : write(warnstr,*) subname, 'aice:', aice
888 0 : call icepack_warnings_add(warnstr)
889 0 : do n = 1, ncat
890 0 : write(warnstr,*) subname, 'n, aicen:', n, aicen(n)
891 0 : call icepack_warnings_add(warnstr)
892 : enddo
893 0 : return
894 : endif
895 : endif ! llimit_aice
896 :
897 : !-----------------------------------------------------------------
898 : ! Identify grid cells with ice.
899 : !-----------------------------------------------------------------
900 :
901 2990543280 : if (ldorebin .and. aice > puny) then
902 :
903 : !-----------------------------------------------------------------
904 : ! Make sure ice in each category is within its thickness bounds.
905 : ! NOTE: The rebin subroutine is needed only in the rare cases
906 : ! when the linear_itd subroutine cannot transfer ice
907 : ! correctly (e.g., very fast ice growth).
908 : !-----------------------------------------------------------------
909 :
910 : call rebin (trcr_depend, &
911 : trcr_base, & ! LCOV_EXCL_LINE
912 : n_trcr_strata, & ! LCOV_EXCL_LINE
913 : nt_strata, & ! LCOV_EXCL_LINE
914 : aicen, trcrn, & ! LCOV_EXCL_LINE
915 : vicen, vsnon, & ! LCOV_EXCL_LINE
916 690510631 : hin_max, Tf )
917 690510631 : if (icepack_warnings_aborted(subname)) return
918 :
919 : endif ! aice > puny
920 :
921 : !-----------------------------------------------------------------
922 : ! Zero out ice categories with very small areas.
923 : !-----------------------------------------------------------------
924 :
925 2990543280 : if (llimit_aice) then
926 : call zap_small_areas (dt, &
927 : aice, aice0, & ! LCOV_EXCL_LINE
928 : aicen, trcrn, & ! LCOV_EXCL_LINE
929 : vicen, vsnon, & ! LCOV_EXCL_LINE
930 : dfpond, & ! LCOV_EXCL_LINE
931 : dfresh, dfsalt, & ! LCOV_EXCL_LINE
932 : dfhocn, & ! LCOV_EXCL_LINE
933 : dfaero_ocn, dfiso_ocn, & ! LCOV_EXCL_LINE
934 : tr_aero, & ! LCOV_EXCL_LINE
935 : tr_pond_topo, & ! LCOV_EXCL_LINE
936 : first_ice, & ! LCOV_EXCL_LINE
937 2990543280 : dflux_bio, Tf )
938 :
939 2990543280 : if (icepack_warnings_aborted(subname)) then
940 0 : write(warnstr,*) subname, 'aice:', aice
941 0 : call icepack_warnings_add(warnstr)
942 0 : do n = 1, ncat
943 0 : write(warnstr,*) subname, 'n, aicen:', n, aicen(n)
944 0 : call icepack_warnings_add(warnstr)
945 : enddo
946 0 : return
947 : endif
948 :
949 : endif ! llimit_aice
950 :
951 : !-------------------------------------------------------------------
952 : ! Zap snow that has out of bounds temperatures
953 : !-------------------------------------------------------------------
954 :
955 : call zap_snow_temperature(dt, aicen, &
956 : trcrn, vsnon, & ! LCOV_EXCL_LINE
957 : dfresh, dfhocn, & ! LCOV_EXCL_LINE
958 : dfaero_ocn, tr_aero, & ! LCOV_EXCL_LINE
959 : dfiso_ocn, & ! LCOV_EXCL_LINE
960 2990543280 : dflux_bio)
961 2990543280 : if (icepack_warnings_aborted(subname)) return
962 :
963 : !-------------------------------------------------------------------
964 : ! Update ice-ocean fluxes for strict conservation
965 : !-------------------------------------------------------------------
966 :
967 2990543280 : if (present(fpond)) &
968 2990543280 : fpond = fpond + dfpond
969 2990543280 : if (present(fresh)) &
970 2990543280 : fresh = fresh + dfresh
971 2990543280 : if (present(fsalt)) &
972 2990543280 : fsalt = fsalt + dfsalt
973 2990543280 : if (present(fhocn)) &
974 2990543280 : fhocn = fhocn + dfhocn
975 2990543280 : if (present(faero_ocn)) then
976 5875487136 : do it = 1, n_aero
977 5769887712 : faero_ocn(it) = faero_ocn(it) + dfaero_ocn(it)
978 : enddo
979 : endif
980 2990543280 : if (present(fiso_ocn)) then
981 2986700400 : if (tr_iso) then
982 116823552 : do it = 1, n_iso
983 116823552 : fiso_ocn(it) = fiso_ocn(it) + dfiso_ocn(it)
984 : enddo
985 : endif
986 : endif
987 2990543280 : if (present(flux_bio)) then
988 4900011984 : do it = 1, nbtrcr
989 2015068128 : flux_bio (it) = flux_bio(it) + dflux_bio(it)
990 : enddo
991 : endif
992 :
993 : end subroutine cleanup_itd
994 :
995 : !=======================================================================
996 :
997 : ! For each ice category in each grid cell, remove ice if the fractional
998 : ! area is less than puny.
999 : !
1000 : ! author: William H. Lipscomb, LANL
1001 :
1002 2990543280 : subroutine zap_small_areas (dt, &
1003 : aice, aice0, & ! LCOV_EXCL_LINE
1004 2990543280 : aicen, trcrn, & ! LCOV_EXCL_LINE
1005 5981086560 : vicen, vsnon, & ! LCOV_EXCL_LINE
1006 : dfpond, & ! LCOV_EXCL_LINE
1007 : dfresh, dfsalt, & ! LCOV_EXCL_LINE
1008 : dfhocn, & ! LCOV_EXCL_LINE
1009 2990543280 : dfaero_ocn, dfiso_ocn, & ! LCOV_EXCL_LINE
1010 : tr_aero, & ! LCOV_EXCL_LINE
1011 : tr_pond_topo, & ! LCOV_EXCL_LINE
1012 2990543280 : first_ice, & ! LCOV_EXCL_LINE
1013 2990543280 : dflux_bio, Tf )
1014 :
1015 : real (kind=dbl_kind), intent(in) :: &
1016 : dt ! time step
1017 :
1018 : real (kind=dbl_kind), intent(inout) :: &
1019 : aice , & ! total ice concentration ! LCOV_EXCL_LINE
1020 : aice0 ! concentration of open water
1021 :
1022 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1023 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1024 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1025 : vsnon ! volume per unit area of snow (m)
1026 :
1027 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1028 : trcrn ! ice tracers
1029 :
1030 : real (kind=dbl_kind), intent(inout) :: &
1031 : dfpond , & ! zapped pond water flux (kg/m^2/s) ! LCOV_EXCL_LINE
1032 : dfresh , & ! zapped fresh water flux (kg/m^2/s) ! LCOV_EXCL_LINE
1033 : dfsalt , & ! zapped salt flux (kg/m^2/s) ! LCOV_EXCL_LINE
1034 : dfhocn ! zapped energy flux ( W/m^2)
1035 :
1036 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1037 : dfaero_ocn ! zapped aerosol flux (kg/m^2/s)
1038 :
1039 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1040 : dfiso_ocn ! zapped isotope flux (kg/m^2/s)
1041 :
1042 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1043 : dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s)
1044 :
1045 : real (kind=dbl_kind), intent(in) :: &
1046 : Tf ! Freezing temperature
1047 :
1048 : logical (kind=log_kind), intent(in) :: &
1049 : tr_aero, & ! aerosol flag ! LCOV_EXCL_LINE
1050 : tr_pond_topo ! pond flag
1051 :
1052 : logical (kind=log_kind), dimension (:), intent(inout) :: &
1053 : first_ice ! For bgc tracers. Set to true if zapping ice
1054 :
1055 : ! local variables
1056 :
1057 : integer (kind=int_kind) :: &
1058 : n, k, it, & !counting indices ! LCOV_EXCL_LINE
1059 : blevels
1060 :
1061 : real (kind=dbl_kind) :: xtmp, sicen ! temporary variables
1062 : real (kind=dbl_kind) :: dvssl, dvint ! temporary variables
1063 : real (kind=dbl_kind) , dimension (1):: trcr_skl
1064 2990543280 : real (kind=dbl_kind) , dimension (nblyr+1):: bvol
1065 :
1066 : character(len=*),parameter :: subname='(zap_small_areas)'
1067 :
1068 : !-----------------------------------------------------------------
1069 : ! I. Zap categories with very small areas.
1070 : !-----------------------------------------------------------------
1071 :
1072 17898297984 : do n = 1, ncat
1073 :
1074 : !-----------------------------------------------------------------
1075 : ! Count categories to be zapped.
1076 : !-----------------------------------------------------------------
1077 :
1078 17898297984 : if (aicen(n) < -puny) then
1079 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1080 0 : call icepack_warnings_add(subname//' Zap ice: negative ice area')
1081 0 : return
1082 14907754704 : elseif (abs(aicen(n)) /= c0 .and. &
1083 : abs(aicen(n)) <= puny) then
1084 :
1085 : !-----------------------------------------------------------------
1086 : ! Account for tracers important for conservation
1087 : !-----------------------------------------------------------------
1088 :
1089 17616383 : if (tr_pond_topo) then
1090 : xtmp = aicen(n) &
1091 1044086 : * trcrn(nt_apnd,n) * trcrn(nt_hpnd,n)
1092 1044086 : dfpond = dfpond - xtmp
1093 : endif
1094 :
1095 17616383 : if (tr_aero) then
1096 2529046 : do it = 1, n_aero
1097 : xtmp = (vicen(n)*(trcrn(nt_aero+2+4*(it-1),n) &
1098 1264523 : + trcrn(nt_aero+3+4*(it-1),n)))/dt
1099 2529046 : dfaero_ocn(it) = dfaero_ocn(it) + xtmp
1100 : enddo
1101 : endif
1102 :
1103 17616383 : if (tr_iso) then
1104 373856 : do it = 1, n_iso
1105 280392 : xtmp = vicen(n)*trcrn(nt_isoice+it-1,n)/dt
1106 373856 : dfiso_ocn(it) = dfiso_ocn(it) + xtmp
1107 : enddo
1108 : endif
1109 :
1110 17616383 : if (skl_bgc .and. nbtrcr > 0) then
1111 0 : blevels = 1
1112 0 : bvol(1) = aicen(n)*sk_l
1113 0 : do it = 1, nbtrcr
1114 0 : trcr_skl(1) = trcrn(bio_index(it),n)
1115 : call zap_small_bgc(blevels, dflux_bio(it), &
1116 0 : dt, bvol(1:blevels), trcr_skl(blevels))
1117 0 : if (icepack_warnings_aborted(subname)) return
1118 : enddo
1119 17616383 : elseif (z_tracers .and. nbtrcr > 0) then
1120 1719995 : blevels = nblyr + 1
1121 15479955 : bvol(:) = vicen(n)/real(nblyr,kind=dbl_kind)*trcrn(nt_fbri,n)
1122 1719995 : bvol(1) = p5*bvol(1)
1123 1719995 : bvol(blevels) = p5*bvol(blevels)
1124 33729015 : do it = 1, nbtrcr
1125 : call zap_small_bgc(blevels, dflux_bio(it), &
1126 32009020 : dt, bvol(1:blevels),trcrn(bio_index(it):bio_index(it)+blevels-1,n))
1127 33729015 : if (icepack_warnings_aborted(subname)) return
1128 : enddo
1129 : endif
1130 :
1131 : !-----------------------------------------------------------------
1132 : ! Zap ice energy and use ocean heat to melt ice
1133 : !-----------------------------------------------------------------
1134 :
1135 98029894 : do k = 1, nilyr
1136 : xtmp = trcrn(nt_qice+k-1,n) / dt &
1137 80413511 : * vicen(n)/real(nilyr,kind=dbl_kind) ! < 0
1138 80413511 : dfhocn = dfhocn + xtmp
1139 98029894 : trcrn(nt_qice+k-1,n) = c0
1140 : enddo ! k
1141 :
1142 : !-----------------------------------------------------------------
1143 : ! Zap ice and snow volume, add water and salt to ocean
1144 : !-----------------------------------------------------------------
1145 :
1146 17616383 : xtmp = (rhoi*vicen(n)) / dt
1147 17616383 : dfresh = dfresh + xtmp
1148 :
1149 17616383 : if (saltflux_option == 'prognostic') then
1150 0 : sicen = c0
1151 0 : do k=1,nilyr
1152 0 : sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
1153 : enddo
1154 0 : xtmp = rhoi*vicen(n)*sicen*p001 / dt
1155 : else
1156 17616383 : xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt
1157 : endif
1158 17616383 : dfsalt = dfsalt + xtmp
1159 :
1160 17616383 : aice0 = aice0 + aicen(n)
1161 17616383 : aicen(n) = c0
1162 17616383 : vicen(n) = c0
1163 17616383 : trcrn(nt_Tsfc,n) = Tf
1164 :
1165 : !-----------------------------------------------------------------
1166 : ! Zap snow
1167 : !-----------------------------------------------------------------
1168 : call zap_snow(dt, &
1169 : trcrn(:,n), vsnon(n), & ! LCOV_EXCL_LINE
1170 : dfresh, dfhocn, & ! LCOV_EXCL_LINE
1171 : dfaero_ocn, tr_aero, & ! LCOV_EXCL_LINE
1172 : dfiso_ocn, & ! LCOV_EXCL_LINE
1173 : dflux_bio, & ! LCOV_EXCL_LINE
1174 17616383 : aicen(n))
1175 17616383 : if (icepack_warnings_aborted(subname)) return
1176 :
1177 : !-----------------------------------------------------------------
1178 : ! Zap tracers
1179 : !-----------------------------------------------------------------
1180 :
1181 17616383 : if (ntrcr >= 2) then
1182 710710345 : do it = 2, ntrcr
1183 710710345 : trcrn(it,n) = c0
1184 : enddo
1185 : endif
1186 17616383 : if (tr_brine) trcrn(nt_fbri,n) = c1
1187 17616383 : if (snwredist(1:3) == 'ITD') then
1188 987756 : do k = 1, nslyr
1189 987756 : trcrn(nt_rhos +k-1,n) = rhosnew
1190 : enddo
1191 : endif
1192 17616383 : if (snwgrain) then
1193 10463886 : do k = 1, nslyr
1194 8719905 : trcrn(nt_smice+k-1,n) = rhos
1195 10463886 : trcrn(nt_rsnw +k-1,n) = rsnw_fall
1196 : enddo
1197 : endif
1198 17616383 : first_ice(n) = .true.
1199 :
1200 : endif ! aicen
1201 : enddo ! n
1202 :
1203 : !-----------------------------------------------------------------
1204 : ! II. Count cells with excess ice (aice > c1) due to roundoff errors.
1205 : ! Zap a little ice in each category so that aice = c1.
1206 : !-----------------------------------------------------------------
1207 :
1208 2990543280 : if (aice > (c1+puny)) then
1209 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1210 0 : call icepack_warnings_add(subname//' Zap ice: excess ice area')
1211 0 : return
1212 2990543280 : elseif (aice > c1 .and. aice < (c1+puny)) then
1213 :
1214 34046400 : do n = 1, ncat
1215 :
1216 : !-----------------------------------------------------------------
1217 : ! Account for tracers important for conservation
1218 : !-----------------------------------------------------------------
1219 :
1220 28391883 : if (tr_pond_topo) then
1221 : xtmp = aicen(n) &
1222 : * trcrn(nt_apnd,n) * trcrn(nt_hpnd,n) & ! LCOV_EXCL_LINE
1223 51252 : * (aice-c1)/aice
1224 51252 : dfpond = dfpond - xtmp
1225 : endif
1226 :
1227 28391883 : if (tr_aero) then
1228 235754 : do it = 1, n_aero
1229 : xtmp = (vsnon(n)*(trcrn(nt_aero +4*(it-1),n) &
1230 : + trcrn(nt_aero+1+4*(it-1),n)) & ! LCOV_EXCL_LINE
1231 : + vicen(n)*(trcrn(nt_aero+2+4*(it-1),n) & ! LCOV_EXCL_LINE
1232 : + trcrn(nt_aero+3+4*(it-1),n))) & ! LCOV_EXCL_LINE
1233 117877 : * (aice-c1)/aice / dt
1234 235754 : dfaero_ocn(it) = dfaero_ocn(it) + xtmp
1235 : enddo ! it
1236 : endif
1237 :
1238 28391883 : if (skl_bgc .and. nbtrcr > 0) then
1239 0 : blevels = 1
1240 0 : bvol(1) = (aice-c1)/aice * sk_l
1241 0 : do it = 1, nbtrcr
1242 0 : trcr_skl(1) = trcrn(bio_index(it),n)
1243 : call zap_small_bgc(blevels, dflux_bio(it), &
1244 0 : dt, bvol(1:blevels), trcr_skl(blevels))
1245 : enddo
1246 28391883 : elseif (z_tracers .and. nbtrcr > 0) then
1247 381505 : blevels = nblyr + 1
1248 3433545 : bvol(:) = (aice-c1)/aice*vicen(n)/real(nblyr,kind=dbl_kind)*trcrn(nt_fbri,n)
1249 381505 : bvol(1) = p5*bvol(1)
1250 381505 : bvol(blevels) = p5*bvol(blevels)
1251 3551145 : do it = 1, nbtrcr
1252 : call zap_small_bgc(blevels, dflux_bio(it), &
1253 3169640 : dt, bvol(1:blevels),trcrn(bio_index(it):bio_index(it)+blevels-1,n))
1254 3551145 : if (icepack_warnings_aborted(subname)) return
1255 : enddo
1256 : ! zap snow zaerosols
1257 381505 : dvssl = p5*vsnon(n)/real(nslyr,kind=dbl_kind) ! snow surface layer
1258 381505 : dvint = vsnon(n) - dvssl ! snow interior
1259 :
1260 3551145 : do it = 1, nbtrcr
1261 : xtmp = (trcrn(bio_index(it)+nblyr+1,n)*dvssl + &
1262 3169640 : trcrn(bio_index(it)+nblyr+2,n)*dvint)*(aice-c1)/aice/dt
1263 3551145 : dflux_bio(it) = dflux_bio(it) + xtmp
1264 : enddo ! it
1265 :
1266 : endif
1267 :
1268 28391883 : if (tr_iso) then
1269 1356460 : do it = 1, n_iso
1270 : xtmp = (vsnon(n)*trcrn(nt_isosno+it-1,n) &
1271 : + vicen(n)*trcrn(nt_isoice+it-1,n)) & ! LCOV_EXCL_LINE
1272 1017345 : * (aice-c1)/aice / dt
1273 1356460 : dfiso_ocn(it) = dfiso_ocn(it) + xtmp
1274 : enddo ! it
1275 : endif
1276 :
1277 : !-----------------------------------------------------------------
1278 : ! Zap ice energy and use ocean heat to melt ice
1279 : !-----------------------------------------------------------------
1280 :
1281 225148824 : do k = 1, nilyr
1282 : xtmp = trcrn(nt_qice+k-1,n) &
1283 : * vicen(n)/real(nilyr,kind=dbl_kind) & ! LCOV_EXCL_LINE
1284 196756941 : * (aice-c1)/aice / dt ! < 0
1285 225148824 : dfhocn = dfhocn + xtmp
1286 : enddo ! k
1287 :
1288 : !-----------------------------------------------------------------
1289 : ! Zap snow energy and use ocean heat to melt snow
1290 : !-----------------------------------------------------------------
1291 :
1292 59038618 : do k = 1, nslyr
1293 : xtmp = trcrn(nt_qsno+k-1,n) &
1294 : * vsnon(n)/real(nslyr,kind=dbl_kind) & ! LCOV_EXCL_LINE
1295 30646735 : * (aice-c1)/aice / dt ! < 0
1296 59038618 : dfhocn = dfhocn + xtmp
1297 : enddo ! k
1298 :
1299 : !-----------------------------------------------------------------
1300 : ! Zap ice and snow volume, add water and salt to ocean
1301 : !-----------------------------------------------------------------
1302 :
1303 : xtmp = (rhoi*vicen(n) + rhos*vsnon(n)) &
1304 28391883 : * (aice-c1)/aice / dt
1305 28391883 : dfresh = dfresh + xtmp
1306 :
1307 28391883 : if (saltflux_option == 'prognostic') then
1308 0 : sicen = c0
1309 0 : do k=1,nilyr
1310 0 : sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
1311 : enddo
1312 : xtmp = rhoi*vicen(n)*sicen*p001 &
1313 0 : * (aice-c1)/aice / dt
1314 : else
1315 : xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 &
1316 28391883 : * (aice-c1)/aice / dt
1317 : endif
1318 28391883 : dfsalt = dfsalt + xtmp
1319 :
1320 28391883 : aicen(n) = aicen(n) * (c1/aice)
1321 28391883 : vicen(n) = vicen(n) * (c1/aice)
1322 34046400 : vsnon(n) = vsnon(n) * (c1/aice)
1323 :
1324 : ! Note: Tracers are unchanged.
1325 :
1326 : enddo ! n
1327 :
1328 : !-----------------------------------------------------------------
1329 : ! Correct aice
1330 : !-----------------------------------------------------------------
1331 :
1332 5654517 : aice = c1
1333 5654517 : aice0 = c0
1334 :
1335 : endif ! aice
1336 :
1337 : end subroutine zap_small_areas
1338 :
1339 : !=======================================================================
1340 :
1341 17616383 : subroutine zap_snow(dt, &
1342 17616383 : trcrn, vsnon, & ! LCOV_EXCL_LINE
1343 : dfresh, dfhocn, & ! LCOV_EXCL_LINE
1344 17616383 : dfaero_ocn, tr_aero, & ! LCOV_EXCL_LINE
1345 17616383 : dfiso_ocn, & ! LCOV_EXCL_LINE
1346 17616383 : dflux_bio, aicen)
1347 :
1348 : real (kind=dbl_kind), intent(in) :: &
1349 : dt ! time step
1350 :
1351 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1352 : trcrn ! ice tracers
1353 :
1354 : real (kind=dbl_kind), intent(in) :: &
1355 : aicen ! ice area fraction
1356 :
1357 : real (kind=dbl_kind), intent(inout) :: &
1358 : vsnon , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
1359 : dfresh , & ! zapped fresh water flux (kg/m^2/s) ! LCOV_EXCL_LINE
1360 : dfhocn ! zapped energy flux ( W/m^2)
1361 :
1362 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1363 : dfaero_ocn ! zapped aerosol flux (kg/m^2/s)
1364 :
1365 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1366 : dfiso_ocn ! zapped isotope flux (kg/m^2/s)
1367 :
1368 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1369 : dflux_bio ! zapped bio tracer flux from biology (mmol/m^2/s)
1370 :
1371 : logical (kind=log_kind), intent(in) :: &
1372 : tr_aero ! aerosol flag
1373 :
1374 : ! local variables
1375 :
1376 : integer (kind=int_kind) :: &
1377 : k, it ! counting indices
1378 :
1379 : real (kind=dbl_kind) :: xtmp, dvssl, dvint
1380 :
1381 : character(len=*),parameter :: subname='(zap_snow)'
1382 :
1383 : ! aerosols
1384 17616383 : if (tr_aero) then
1385 2529046 : do it = 1, n_aero
1386 : xtmp = (vsnon*(trcrn(nt_aero +4*(it-1)) &
1387 1264523 : + trcrn(nt_aero+1+4*(it-1))))/dt
1388 2529046 : dfaero_ocn(it) = dfaero_ocn(it) + xtmp
1389 : enddo ! it
1390 : endif ! tr_aero
1391 :
1392 : ! isotopes
1393 17616383 : if (tr_iso) then
1394 373856 : do it = 1, n_iso
1395 280392 : xtmp = vsnon*trcrn(nt_isosno+it-1)/dt
1396 373856 : dfiso_ocn(it) = dfiso_ocn(it) + xtmp
1397 : enddo ! it
1398 : endif ! tr_iso
1399 :
1400 17616383 : if (z_tracers) then
1401 1719995 : dvssl = p5*vsnon/real(nslyr,kind=dbl_kind) ! snow surface layer
1402 1719995 : dvint = vsnon - dvssl ! snow interior
1403 :
1404 33729015 : do it = 1, nbtrcr
1405 : xtmp = (trcrn(bio_index(it)+nblyr+1)*dvssl + &
1406 32009020 : trcrn(bio_index(it)+nblyr+2)*dvint)/dt
1407 33729015 : dflux_bio(it) = dflux_bio(it) + xtmp
1408 : enddo ! it
1409 :
1410 : endif ! z_tracers
1411 :
1412 : ! snow enthalpy tracer
1413 43373850 : do k = 1, nslyr
1414 : xtmp = trcrn(nt_qsno+k-1) / dt &
1415 25757467 : * vsnon/real(nslyr,kind=dbl_kind) ! < 0
1416 25757467 : dfhocn = dfhocn + xtmp
1417 43373850 : trcrn(nt_qsno+k-1) = c0
1418 : enddo ! k
1419 :
1420 : ! snow volume
1421 17616383 : xtmp = (rhos*vsnon) / dt
1422 17616383 : dfresh = dfresh + xtmp
1423 17616383 : vsnon = c0
1424 :
1425 17616383 : end subroutine zap_snow
1426 :
1427 : !=======================================================================
1428 :
1429 2990543280 : subroutine zap_snow_temperature(dt, aicen, &
1430 2990543280 : trcrn, vsnon, & ! LCOV_EXCL_LINE
1431 : dfresh, dfhocn, & ! LCOV_EXCL_LINE
1432 2990543280 : dfaero_ocn, tr_aero, & ! LCOV_EXCL_LINE
1433 2990543280 : dfiso_ocn, & ! LCOV_EXCL_LINE
1434 2990543280 : dflux_bio)
1435 :
1436 : real (kind=dbl_kind), intent(in) :: &
1437 : dt ! time step
1438 :
1439 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1440 : aicen ! concentration of ice
1441 :
1442 : real (kind=dbl_kind), dimension(:), intent(inout) :: &
1443 : vsnon ! volume per unit area of snow (m)
1444 :
1445 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1446 : trcrn ! ice tracers
1447 :
1448 : real (kind=dbl_kind), intent(inout) :: &
1449 : dfresh , & ! zapped fresh water flux (kg/m^2/s) ! LCOV_EXCL_LINE
1450 : dfhocn ! zapped energy flux ( W/m^2)
1451 :
1452 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1453 : dfaero_ocn ! zapped aerosol flux (kg/m^2/s)
1454 :
1455 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1456 : dfiso_ocn ! zapped isotope flux (kg/m^2/s)
1457 :
1458 : real (kind=dbl_kind), dimension (:),intent(inout) :: &
1459 : dflux_bio ! zapped biology flux (mmol/m^2/s)
1460 :
1461 : logical (kind=log_kind), intent(in) :: &
1462 : tr_aero ! aerosol flag
1463 :
1464 : ! local variables
1465 :
1466 : integer (kind=int_kind) :: &
1467 : n, k ! counting indices
1468 :
1469 : real (kind=dbl_kind) :: &
1470 : rnslyr , & ! real(nslyr) ! LCOV_EXCL_LINE
1471 : hsn , & ! snow thickness (m) ! LCOV_EXCL_LINE
1472 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
1473 : zTsn , & ! snow layer temperature (C) ! LCOV_EXCL_LINE
1474 : Tmax ! maximum allowed snow temperature
1475 :
1476 : logical :: &
1477 : l_zap ! logical whether zap snow
1478 :
1479 : character(len=*),parameter :: subname='(zap_snow_temperature)'
1480 :
1481 2990543280 : rnslyr = real(nslyr,kind=dbl_kind)
1482 :
1483 17898297984 : do n = 1, ncat
1484 :
1485 : !-----------------------------------------------------------------
1486 : ! Determine cells to zap
1487 : !-----------------------------------------------------------------
1488 :
1489 14907754704 : l_zap = .false.
1490 :
1491 14907754704 : if (aicen(n) > puny) then
1492 :
1493 : ! snow thickness
1494 3225633860 : hsn = vsnon(n) / aicen(n)
1495 :
1496 : ! check each snow layer - zap all if one is bad
1497 6818286912 : do k = 1, nslyr
1498 :
1499 : ! snow enthalpy and max temperature
1500 3592653052 : if (hsn > hs_min) then
1501 : ! zqsn < 0
1502 2923185352 : zqsn = trcrn(nt_qsno+k-1,n)
1503 2923185352 : Tmax = -zqsn*puny*rnslyr / (rhos*cp_ice*vsnon(n))
1504 : else
1505 669467700 : zqsn = -rhos * Lfresh
1506 669467700 : Tmax = puny
1507 : endif
1508 :
1509 : ! snow temperature
1510 3592653052 : zTsn = (Lfresh + zqsn/rhos)/cp_ice
1511 :
1512 : ! check for zapping
1513 6818286912 : if (zTsn < Tmin .or. zTsn > Tmax) then
1514 0 : l_zap = .true.
1515 0 : write(warnstr,*) subname, "zap_snow_temperature: temperature out of bounds!"
1516 0 : call icepack_warnings_add(warnstr)
1517 0 : write(warnstr,*) subname, "k:" , k
1518 0 : call icepack_warnings_add(warnstr)
1519 0 : write(warnstr,*) subname, "zTsn:", zTsn
1520 0 : call icepack_warnings_add(warnstr)
1521 0 : write(warnstr,*) subname, "Tmin:", Tmin
1522 0 : call icepack_warnings_add(warnstr)
1523 0 : write(warnstr,*) subname, "Tmax:", Tmax
1524 0 : call icepack_warnings_add(warnstr)
1525 0 : write(warnstr,*) subname, "zqsn:", zqsn
1526 0 : call icepack_warnings_add(warnstr)
1527 : endif
1528 :
1529 : enddo ! k
1530 :
1531 : endif ! aicen > puny
1532 :
1533 : !-----------------------------------------------------------------
1534 : ! Zap the cells
1535 : !-----------------------------------------------------------------
1536 14907754704 : if (l_zap) &
1537 : call zap_snow(dt, & ! LCOV_EXCL_LINE
1538 : trcrn(:,n), vsnon(n), & ! LCOV_EXCL_LINE
1539 : dfresh, dfhocn, & ! LCOV_EXCL_LINE
1540 : dfaero_ocn, tr_aero, & ! LCOV_EXCL_LINE
1541 : dfiso_ocn, & ! LCOV_EXCL_LINE
1542 : dflux_bio, & ! LCOV_EXCL_LINE
1543 0 : aicen(n))
1544 17898297984 : if (icepack_warnings_aborted(subname)) return
1545 :
1546 : enddo ! n
1547 :
1548 : end subroutine zap_snow_temperature
1549 :
1550 : !=======================================================================
1551 : !autodocument_start icepack_init_itd
1552 : ! Initialize area fraction and thickness boundaries for the itd model
1553 : !
1554 : ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
1555 : ! C. M. Bitz, UW
1556 :
1557 4222 : subroutine icepack_init_itd(hin_max)
1558 :
1559 : real (kind=dbl_kind), intent(out) :: &
1560 : hin_max(0:ncat) ! category limits (m)
1561 :
1562 : !autodocument_end
1563 :
1564 : ! local variables
1565 :
1566 : integer (kind=int_kind) :: &
1567 : n ! thickness category index
1568 :
1569 : real (kind=dbl_kind) :: &
1570 : cc1, cc2, cc3, & ! parameters for kcatbound = 0 ! LCOV_EXCL_LINE
1571 : x1 , & ! LCOV_EXCL_LINE
1572 : rn , & ! real(n) ! LCOV_EXCL_LINE
1573 : rncat , & ! real(ncat) ! LCOV_EXCL_LINE
1574 : d1 , & ! parameters for kcatbound = 1 (m) ! LCOV_EXCL_LINE
1575 : d2 , & ! ! LCOV_EXCL_LINE
1576 : b1 , & ! parameters for kcatbound = 3 ! LCOV_EXCL_LINE
1577 : b2 , & ! ! LCOV_EXCL_LINE
1578 : b3
1579 :
1580 : real (kind=dbl_kind), dimension(5) :: wmo5 ! data for wmo itd
1581 : real (kind=dbl_kind), dimension(6) :: wmo6 ! data for wmo itd
1582 : real (kind=dbl_kind), dimension(7) :: wmo7 ! data for wmo itd
1583 :
1584 : character(len=*),parameter :: subname='(icepack_init_itd)'
1585 :
1586 : ! thinnest 3 categories combined
1587 : data wmo5 / 0.30_dbl_kind, 0.70_dbl_kind, &
1588 : 1.20_dbl_kind, 2.00_dbl_kind, & ! LCOV_EXCL_LINE
1589 : 999._dbl_kind /
1590 : ! thinnest 2 categories combined
1591 : data wmo6 / 0.15_dbl_kind, &
1592 : 0.30_dbl_kind, 0.70_dbl_kind, & ! LCOV_EXCL_LINE
1593 : 1.20_dbl_kind, 2.00_dbl_kind, & ! LCOV_EXCL_LINE
1594 : 999._dbl_kind /
1595 : !echmod wmo6a
1596 : ! data wmo6 /0.30_dbl_kind, 0.70_dbl_kind, &
1597 : ! 1.20_dbl_kind, 2.00_dbl_kind, & ! LCOV_EXCL_LINE
1598 : ! 4.56729_dbl_kind, & ! LCOV_EXCL_LINE
1599 : ! 999._dbl_kind /
1600 : ! all thickness categories
1601 : data wmo7 / 0.10_dbl_kind, 0.15_dbl_kind, &
1602 : 0.30_dbl_kind, 0.70_dbl_kind, & ! LCOV_EXCL_LINE
1603 : 1.20_dbl_kind, 2.00_dbl_kind, & ! LCOV_EXCL_LINE
1604 : 999._dbl_kind /
1605 :
1606 4222 : rncat = real(ncat, kind=dbl_kind)
1607 4222 : d1 = 3.0_dbl_kind / rncat
1608 4222 : d2 = 0.5_dbl_kind / rncat
1609 4222 : b1 = p1 ! asymptotic category width (m)
1610 4222 : b2 = c3 ! thickness for which participation function is small (m)
1611 4222 : b3 = max(rncat*(rncat-1), c2*b2/b1)
1612 :
1613 : !-----------------------------------------------------------------
1614 : ! Choose category boundaries based on one of four options.
1615 : !
1616 : ! The first formula (kcatbound = 0) was used in Lipscomb (2001)
1617 : ! and in CICE versions 3.0 and 3.1.
1618 : !
1619 : ! The second formula is more user-friendly in the sense that it
1620 : ! is easy to obtain round numbers for category boundaries:
1621 : !
1622 : ! H(n) = n * [d1 + d2*(n-1)]
1623 : !
1624 : ! Default values are d1 = 300/ncat, d2 = 50/ncat.
1625 : ! For ncat = 5, boundaries in cm are 60, 140, 240, 360, which are
1626 : ! close to the standard values given by the first formula.
1627 : ! For ncat = 10, boundaries in cm are 30, 70, 120, 180, 250, 330,
1628 : ! 420, 520, 630.
1629 : !
1630 : ! The third option provides support for World Meteorological
1631 : ! Organization classification based on thickness. The full
1632 : ! WMO thickness distribution is used if ncat = 7; if ncat=5
1633 : ! or ncat = 6, some of the thinner categories are combined.
1634 : ! For ncat = 5, boundaries are 30, 70, 120, 200, >200 cm.
1635 : ! For ncat = 6, boundaries are 15, 30, 70, 120, 200, >200 cm.
1636 : ! For ncat = 7, boundaries are 10, 15, 30, 70, 120, 200, >200 cm.
1637 : !
1638 : ! The fourth formula asymptotes to a particular category width as
1639 : ! the number of categories increases, given by the parameter b1.
1640 : ! The parameter b3 is computed so that the category boundaries
1641 : ! are even numbers.
1642 : !
1643 : ! H(n) = b1 * [n + b3*n*(n+1)/(2*N*(N-1))] for N=ncat
1644 : !
1645 : ! kcatbound=-1 is available only for 1-category runs, with
1646 : ! boundaries 0 and 100 m.
1647 : !-----------------------------------------------------------------
1648 :
1649 4222 : if (kcatbound == -1) then ! single category
1650 222 : hin_max(0) = c0
1651 222 : hin_max(1) = c100
1652 :
1653 4000 : elseif (kcatbound == 0) then ! original scheme
1654 :
1655 3493 : if (kitd == 1) then
1656 : ! linear remapping itd category limits
1657 3395 : cc1 = c3/rncat
1658 3395 : cc2 = c15*cc1
1659 3395 : cc3 = c3
1660 :
1661 3395 : hin_max(0) = c0 ! minimum ice thickness, m
1662 : else
1663 : ! delta function itd category limits
1664 98 : cc1 = max(1.1_dbl_kind/rncat,hi_min)
1665 98 : cc2 = c25*cc1
1666 98 : cc3 = 2.25_dbl_kind
1667 :
1668 : ! hin_max(0) should not be zero
1669 : ! use some caution in making it less than 0.10
1670 98 : hin_max(0) = hi_min ! minimum ice thickness, m
1671 : endif ! kitd
1672 :
1673 20958 : do n = 1, ncat
1674 17465 : x1 = real(n-1,kind=dbl_kind) / rncat
1675 : hin_max(n) = hin_max(n-1) &
1676 20958 : + cc1 + cc2*(c1 + tanh(cc3*(x1-c1)))
1677 : enddo
1678 :
1679 507 : elseif (kcatbound == 1) then ! new scheme
1680 :
1681 122 : hin_max(0) = c0
1682 732 : do n = 1, ncat
1683 610 : rn = real(n, kind=dbl_kind)
1684 732 : hin_max(n) = rn * (d1 + (rn-c1)*d2)
1685 : enddo
1686 :
1687 385 : elseif (kcatbound == 2) then ! WMO standard
1688 :
1689 205 : if (ncat == 5) then
1690 9 : hin_max(0) = c0
1691 54 : do n = 1, ncat
1692 54 : hin_max(n) = wmo5(n)
1693 : enddo
1694 196 : elseif (ncat == 6) then
1695 196 : hin_max(0) = c0
1696 1372 : do n = 1, ncat
1697 1372 : hin_max(n) = wmo6(n)
1698 : enddo
1699 0 : elseif (ncat == 7) then
1700 0 : hin_max(0) = c0
1701 0 : do n = 1, ncat
1702 0 : hin_max(n) = wmo7(n)
1703 : enddo
1704 : else
1705 0 : call icepack_warnings_add(subname//' kcatbound=2 (WMO) must have ncat=5, 6 or 7')
1706 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
1707 0 : return
1708 : endif
1709 :
1710 180 : elseif (kcatbound == 3) then ! asymptotic scheme
1711 :
1712 180 : hin_max(0) = c0
1713 1440 : do n = 1, ncat
1714 1260 : rn = real(n, kind=dbl_kind)
1715 1440 : hin_max(n) = b1 * (rn + b3*rn*(rn+c1)/(c2*rncat*(rncat-c1)))
1716 : enddo
1717 :
1718 : endif ! kcatbound
1719 :
1720 4222 : if (kitd == 1) then
1721 3792 : hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
1722 : endif
1723 :
1724 : end subroutine icepack_init_itd
1725 :
1726 : !=======================================================================
1727 : !autodocument_start icepack_init_itd_hist
1728 : ! Initialize area fraction and thickness boundaries for the itd model
1729 : !
1730 : ! authors: William H. Lipscomb and Elizabeth C. Hunke, LANL
1731 : ! C. M. Bitz, UW
1732 :
1733 375 : subroutine icepack_init_itd_hist (hin_max, c_hi_range)
1734 :
1735 : real (kind=dbl_kind), intent(in) :: &
1736 : hin_max(0:ncat) ! category limits (m)
1737 :
1738 : character (len=35), intent(out) :: &
1739 : c_hi_range(ncat) ! string for history output
1740 :
1741 : !autodocument_end
1742 :
1743 : ! local variables
1744 :
1745 : integer (kind=int_kind) :: &
1746 : n ! thickness category index
1747 :
1748 : character(len=8) :: c_hinmax1,c_hinmax2
1749 : character(len=2) :: c_nc
1750 :
1751 : character(len=*),parameter :: subname='(icepack_init_itd_hist)'
1752 :
1753 375 : write(warnstr,*) ' '
1754 375 : call icepack_warnings_add(warnstr)
1755 375 : write(warnstr,*) subname
1756 375 : call icepack_warnings_add(warnstr)
1757 375 : write(warnstr,*) 'hin_max(n-1) < Cat n < hin_max(n)'
1758 375 : call icepack_warnings_add(warnstr)
1759 2237 : do n = 1, ncat
1760 1862 : write(warnstr,*) hin_max(n-1),' < Cat ',n, ' < ',hin_max(n)
1761 1862 : call icepack_warnings_add(warnstr)
1762 : ! Write integer n to character string
1763 1862 : write (c_nc, '(i2)') n
1764 :
1765 : ! Write hin_max to character string
1766 1862 : write (c_hinmax1, '(f7.3)') hin_max(n-1)
1767 1862 : write (c_hinmax2, '(f7.3)') hin_max(n)
1768 :
1769 : ! Save character string to write to history file
1770 2237 : c_hi_range(n)=c_hinmax1//'m < hi Cat '//c_nc//' < '//c_hinmax2//'m'
1771 : enddo
1772 :
1773 375 : write(warnstr,*) ' '
1774 375 : call icepack_warnings_add(warnstr)
1775 :
1776 375 : end subroutine icepack_init_itd_hist
1777 :
1778 : !=======================================================================
1779 : !autodocument_start icepack_aggregate
1780 : ! Aggregate ice state variables over thickness categories.
1781 : !
1782 : ! authors: C. M. Bitz, UW
1783 : ! W. H. Lipscomb, LANL
1784 :
1785 11112233492 : subroutine icepack_aggregate(aicen, trcrn, &
1786 5556116746 : vicen, vsnon, & ! LCOV_EXCL_LINE
1787 0 : aice, trcr, & ! LCOV_EXCL_LINE
1788 : vice, vsno, & ! LCOV_EXCL_LINE
1789 : aice0, & ! LCOV_EXCL_LINE
1790 5556116746 : trcr_depend, & ! LCOV_EXCL_LINE
1791 5556116746 : trcr_base, & ! LCOV_EXCL_LINE
1792 5556116746 : n_trcr_strata, & ! LCOV_EXCL_LINE
1793 5556116746 : nt_strata, Tf)
1794 :
1795 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1796 : aicen , & ! concentration of ice ! LCOV_EXCL_LINE
1797 : vicen , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1798 : vsnon ! volume per unit area of snow (m)
1799 :
1800 : real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
1801 : trcrn ! ice tracers
1802 :
1803 : integer (kind=int_kind), dimension (:), intent(in) :: &
1804 : trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon ! LCOV_EXCL_LINE
1805 : n_trcr_strata ! number of underlying tracer layers
1806 :
1807 : real (kind=dbl_kind), dimension (:,:), intent(in) :: &
1808 : trcr_base ! = 0 or 1 depending on tracer dependency
1809 : ! argument 2: (1) aice, (2) vice, (3) vsno
1810 :
1811 : integer (kind=int_kind), dimension (:,:), intent(in) :: &
1812 : nt_strata ! indices of underlying tracer layers
1813 :
1814 : real (kind=dbl_kind), intent(out) :: &
1815 : aice , & ! concentration of ice ! LCOV_EXCL_LINE
1816 : vice , & ! volume per unit area of ice (m) ! LCOV_EXCL_LINE
1817 : vsno , & ! volume per unit area of snow (m) ! LCOV_EXCL_LINE
1818 : aice0 ! concentration of open water
1819 :
1820 : real (kind=dbl_kind), dimension (:), intent(out) :: &
1821 : trcr ! ice tracers
1822 :
1823 : real (kind=dbl_kind), intent(in) :: &
1824 : Tf ! freezing temperature
1825 :
1826 : !autodocument_end
1827 :
1828 : ! local variables
1829 :
1830 : integer (kind=int_kind) :: &
1831 : n, it, itl, & ! loop indices ! LCOV_EXCL_LINE
1832 : ntr ! tracer index
1833 :
1834 : real (kind=dbl_kind), dimension (:), allocatable :: &
1835 5556116746 : atrcr ! sum of aicen*trcrn or vicen*trcrn or vsnon*trcrn
1836 :
1837 : real (kind=dbl_kind) :: &
1838 : atrcrn ! category value
1839 :
1840 : character(len=*),parameter :: subname='(icepack_aggregate)'
1841 :
1842 : !-----------------------------------------------------------------
1843 : ! Initialize
1844 : !-----------------------------------------------------------------
1845 :
1846 5556116746 : aice0 = c1
1847 5556116746 : aice = c0
1848 5556116746 : vice = c0
1849 5556116746 : vsno = c0
1850 :
1851 5556116746 : allocate (atrcr(ntrcr))
1852 :
1853 : !-----------------------------------------------------------------
1854 : ! Aggregate
1855 : !-----------------------------------------------------------------
1856 :
1857 >19980*10^7 : atrcr(:) = c0
1858 :
1859 33234597382 : do n = 1, ncat
1860 :
1861 27678480636 : aice = aice + aicen(n)
1862 27678480636 : vice = vice + vicen(n)
1863 27678480636 : vsno = vsno + vsnon(n)
1864 :
1865 >10031*10^8 : do it = 1, ntrcr
1866 : atrcrn = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
1867 : + trcr_base(it,2) * vicen(n) & ! LCOV_EXCL_LINE
1868 >96995*10^7 : + trcr_base(it,3) * vsnon(n))
1869 >96995*10^7 : if (n_trcr_strata(it) > 0) then ! additional tracer layers
1870 >20928*10^7 : do itl = 1, n_trcr_strata(it)
1871 >13059*10^7 : ntr = nt_strata(it,itl)
1872 >20928*10^7 : atrcrn = atrcrn * trcrn(ntr,n)
1873 : enddo
1874 : endif
1875 >99763*10^7 : atrcr(it) = atrcr(it) + atrcrn
1876 : enddo ! ntrcr
1877 : enddo ! ncat
1878 :
1879 : ! Open water fraction
1880 5556116746 : aice0 = max (c1 - aice, c0)
1881 :
1882 : ! Tracers
1883 : call icepack_compute_tracers(trcr_depend, &
1884 : atrcr, aice, & ! LCOV_EXCL_LINE
1885 : vice , vsno, & ! LCOV_EXCL_LINE
1886 : trcr_base, n_trcr_strata, & ! LCOV_EXCL_LINE
1887 5556116746 : nt_strata, trcr, Tf)
1888 5556116746 : if (icepack_warnings_aborted(subname)) return
1889 :
1890 5556116746 : deallocate (atrcr)
1891 :
1892 5556116746 : end subroutine icepack_aggregate
1893 :
1894 : !=======================================================================
1895 :
1896 : end module icepack_itd
1897 :
1898 : !=======================================================================
1899 :
1900 :
1901 :
1902 :
1903 :
1904 :
1905 :
1906 :
1907 :
|
