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