LCOV - code coverage report
Current view: top level - columnphysics - icepack_itd.F90 (source / functions) Hit Total Coverage
Test: 200804-015008:4c42a82e3d:3:base,travis,quick Lines: 516 741 69.64 %
Date: 2020-08-03 20:10:57 Functions: 14 14 100.00 %

          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             : 

Generated by: LCOV version 1.14-6-g40580cd