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