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