LCOV - code coverage report
Current view: top level - columnphysics - icepack_therm_itd.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 744 846 87.94 %
Date: 2020-07-03 20:11:40 Functions: 6 6 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! Thermo calculations after call to coupler, related to ITD:
       4             : ! ice thickness redistribution, lateral growth and melting.
       5             : !
       6             : ! NOTE: The thermodynamic calculation is split in two for load balancing.
       7             : !       First icepack_therm_vertical computes vertical growth rates and coupler
       8             : !       fluxes.  Then icepack_therm_itd does thermodynamic calculations not
       9             : !       needed for coupling.
      10             : !       
      11             : ! authors William H. Lipscomb, LANL
      12             : !         C. M. Bitz, UW
      13             : !         Elizabeth C. Hunke, LANL
      14             : !
      15             : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
      16             : ! 2004: Block structure added by William Lipscomb  
      17             : ! 2006: Streamlined for efficiency by Elizabeth Hunke
      18             : ! 2014: Column package created by Elizabeth Hunke
      19             : !
      20             :       module icepack_therm_itd
      21             : 
      22             :       use icepack_kinds
      23             :       use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10
      24             :       use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
      25             :       use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
      26             :       use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
      27             :       use icepack_parameters, only: rhosi, conserv_check
      28             :       use icepack_parameters, only: kitd, ktherm, heat_capacity
      29             :       use icepack_parameters, only: z_tracers, solve_zsal
      30             : 
      31             :       use icepack_tracers, only: ntrcr, nbtrcr
      32             :       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
      33             :       use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
      34             :       use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd
      35             :       use icepack_tracers, only: nt_alvl, nt_vlvl
      36             :       use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
      37             :       use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
      38             :       use icepack_tracers, only: n_aero, n_iso
      39             :       use icepack_tracers, only: bio_index
      40             : 
      41             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      42             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      43             : 
      44             :       use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd,  get_subdt_fsd    
      45             :       use icepack_itd, only: reduce_area, cleanup_itd
      46             :       use icepack_itd, only: aggregate_area, shift_ice
      47             :       use icepack_itd, only: column_sum, column_conservation_check
      48             :       use icepack_isotope, only: isoice_alpha, isotope_frac_method
      49             :       use icepack_mushy_physics, only: liquidus_temperature_mush, enthalpy_mush
      50             :       use icepack_therm_shared, only: hfrazilmin
      51             :       use icepack_therm_shared, only: hi_min
      52             :       use icepack_zbgc, only: add_new_ice_bgc
      53             :       use icepack_zbgc, only: lateral_melt_bgc               
      54             :  
      55             :       implicit none
      56             :       
      57             :       private
      58             :       public :: linear_itd, &
      59             :                 add_new_ice, &
      60             :                 lateral_melt, &
      61             :                 icepack_step_therm2
      62             : 
      63             : !=======================================================================
      64             : 
      65             :       contains
      66             : 
      67             : !=======================================================================
      68             : !
      69             : ! ITD scheme that shifts ice among categories
      70             : !
      71             : ! See Lipscomb, W. H.  Remapping the thickness distribution in sea
      72             : !     ice models. 2001, J. Geophys. Res., Vol 106, 13989--14000.
      73             : !
      74             : ! Using the thermodynamic "velocities", interpolate to find the
      75             : ! velocities in thickness space at the category boundaries, and
      76             : ! compute the new locations of the boundaries.  Then for each
      77             : ! category, compute the thickness distribution function,  g(h),
      78             : ! between hL and hR, the left and right boundaries of the category.
      79             : ! Assume g(h) is a linear polynomial that satisfies two conditions:
      80             : !
      81             : ! (1) The ice area implied by g(h) equals aicen(n).
      82             : ! (2) The ice volume implied by g(h) equals aicen(n)*hicen(n).
      83             : !
      84             : ! Given g(h), at each boundary compute the ice area and volume lying
      85             : ! between the original and new boundary locations.  Transfer area
      86             : ! and volume across each boundary in the appropriate direction, thus
      87             : ! restoring the original boundaries.
      88             : !
      89             : ! authors: William H. Lipscomb, LANL
      90             : !          Elizabeth C. Hunke, LANL
      91             : 
      92      852950 :       subroutine linear_itd (ncat,        hin_max,     &
      93             :                              nilyr,       nslyr,       &
      94      852950 :                              ntrcr,       trcr_depend, & 
      95      852950 :                              trcr_base,   n_trcr_strata,&
      96      852950 :                              nt_strata,                &
      97      852950 :                              aicen_init,  vicen_init,  & 
      98     1705900 :                              aicen,       trcrn,       & 
      99      852950 :                              vicen,       vsnon,       & 
     100             :                              aice,        aice0,       & 
     101             :                              fpond                     )
     102             : 
     103             :       integer (kind=int_kind), intent(in) :: &
     104             :          ncat    , & ! number of thickness categories
     105             :          nilyr   , & ! number of ice layers
     106             :          nslyr   , & ! number of snow layers
     107             :          ntrcr       ! number of tracers in use
     108             : 
     109             :       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
     110             :          hin_max      ! category boundaries (m)
     111             : 
     112             :       integer (kind=int_kind), dimension (:), intent(in) :: &
     113             :          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
     114             :          n_trcr_strata  ! number of underlying tracer layers
     115             : 
     116             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
     117             :          trcr_base      ! = 0 or 1 depending on tracer dependency
     118             :                         ! argument 2:  (1) aice, (2) vice, (3) vsno
     119             : 
     120             :       integer (kind=int_kind), dimension (:,:), intent(in) :: &
     121             :          nt_strata      ! indices of underlying tracer layers
     122             : 
     123             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
     124             :          aicen_init, & ! initial ice concentration (before vertical thermo)
     125             :          vicen_init    ! initial ice volume               (m)
     126             : 
     127             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     128             :          aicen  , & ! ice concentration
     129             :          vicen  , & ! volume per unit area of ice      (m)
     130             :          vsnon      ! volume per unit area of snow     (m)
     131             : 
     132             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
     133             :          trcrn     ! ice tracers
     134             : 
     135             :       real (kind=dbl_kind), intent(inout) :: &
     136             :          aice  , & ! concentration of ice
     137             :          aice0 , & ! concentration of open water
     138             :          fpond     ! fresh water flux to ponds (kg/m^2/s)
     139             : 
     140             :       ! local variables
     141             : 
     142             :       integer (kind=int_kind) :: &
     143             :          n, nd        , & ! category indices
     144             :          k                ! ice layer index
     145             : 
     146             :       real (kind=dbl_kind) :: &
     147      291772 :          slope        , & ! rate of change of dhice with hice
     148      291772 :          dh0          , & ! change in ice thickness at h = 0
     149      291772 :          da0          , & ! area melting from category 1
     150      291772 :          damax        , & ! max allowed reduction in category 1 area
     151      291772 :          etamin, etamax,& ! left and right limits of integration
     152      291772 :          x1           , & ! etamax - etamin
     153      291772 :          x2           , & ! (etamax^2 - etamin^2) / 2
     154      291772 :          x3           , & ! (etamax^3 - etamin^3) / 3
     155      291772 :          wk1, wk2         ! temporary variables
     156             : 
     157             :       real (kind=dbl_kind), dimension(0:ncat) :: &
     158     3164760 :          hbnew            ! new boundary locations
     159             : 
     160             :       real (kind=dbl_kind), dimension(ncat) :: &
     161     2872988 :          g0           , & ! constant coefficient in g(h)
     162     2872988 :          g1           , & ! linear coefficient in g(h)
     163     2872988 :          hL           , & ! left end of range over which g(h) > 0
     164     2872988 :          hR               ! right end of range over which g(h) > 0
     165             : 
     166             :       real (kind=dbl_kind), dimension(ncat) :: &
     167     2872988 :          hicen        , & ! ice thickness for each cat     (m)
     168     2872988 :          hicen_init   , & ! initial ice thickness for each cat (pre-thermo)
     169     2872988 :          dhicen       , & ! thickness change for remapping (m)
     170     3456532 :          daice        , & ! ice area transferred across boundary
     171     2872988 :          dvice            ! ice volume transferred across boundary
     172             : 
     173             :       real (kind=dbl_kind), dimension (ncat) :: &
     174     2872988 :          eicen, &     ! energy of melting for each ice layer (J/m^2)
     175     2872988 :          esnon, &     ! energy of melting for each snow layer (J/m^2)
     176     2872988 :          vbrin, &     ! ice volume defined by brine height (m)
     177     2872988 :          sicen        ! Bulk salt in h ice (ppt*m)
     178             : 
     179             :       real (kind=dbl_kind) :: &
     180      291772 :          vice_init, vice_final, & ! ice volume summed over categories
     181      291772 :          vsno_init, vsno_final, & ! snow volume summed over categories
     182      291772 :          eice_init, eice_final, & ! ice energy summed over categories
     183      291772 :          esno_init, esno_final, & ! snow energy summed over categories
     184      291772 :          sice_init, sice_final, & ! ice bulk salinity summed over categories
     185      291772 :          vbri_init, vbri_final    ! briny ice volume summed over categories
     186             : 
     187             :       ! NOTE: Third index of donor, daice, dvice should be ncat-1,
     188             :       !       except that compilers would have trouble when ncat = 1 
     189             :       integer (kind=int_kind), dimension(ncat) :: &
     190     1144722 :          donor            ! donor category index
     191             : 
     192             :       logical (kind=log_kind) :: &
     193             :          remap_flag       ! remap ITD if remap_flag is true
     194             : 
     195             :       character (len=char_len) :: &
     196             :          fieldid           ! field identifier
     197             : 
     198             :       logical (kind=log_kind), parameter :: &
     199             :          print_diags = .false.    ! if true, prints when remap_flag=F
     200             : 
     201             :       character(len=*),parameter :: subname='(linear_itd)'
     202             : 
     203             :       !-----------------------------------------------------------------
     204             :       ! Initialize
     205             :       !-----------------------------------------------------------------
     206             : 
     207      852950 :       hin_max(ncat) = 999.9_dbl_kind ! arbitrary big number
     208             : 
     209     5117700 :       do n = 1, ncat
     210     4264750 :          donor(n) = 0
     211     4264750 :          daice(n) = c0
     212     5117700 :          dvice(n) = c0
     213             :       enddo
     214             : 
     215             :       !-----------------------------------------------------------------
     216             :       ! Compute volume and energy sums that linear remapping should
     217             :       !  conserve.
     218             :       !-----------------------------------------------------------------
     219             : 
     220      852950 :       if (conserv_check) then
     221             : 
     222      289584 :       do n = 1, ncat
     223             : 
     224      241320 :          eicen(n) = c0
     225      241320 :          esnon(n) = c0
     226      241320 :          vbrin(n) = c0
     227      241320 :          sicen(n) = c0
     228             : 
     229     1930560 :          do k = 1, nilyr
     230     1226400 :             eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
     231     2543760 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     232             :          enddo
     233     1447920 :          do k = 1, nslyr
     234      876000 :             esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
     235     1885920 :                      * vsnon(n)/real(nslyr,kind=dbl_kind)
     236             :          enddo
     237             : 
     238      241320 :          if (tr_brine) then
     239           0 :             vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
     240           0 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     241             :          endif
     242             : 
     243     1978824 :          do k = 1, nilyr
     244     1226400 :             sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
     245     2543760 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     246             :          enddo
     247             : 
     248             :       enddo ! n
     249             : 
     250       48264 :       call column_sum (ncat, vicen, vice_init)
     251       48264 :       if (icepack_warnings_aborted(subname)) return
     252       48264 :       call column_sum (ncat, vsnon, vsno_init)
     253       48264 :       if (icepack_warnings_aborted(subname)) return
     254       48264 :       call column_sum (ncat, eicen, eice_init)
     255       48264 :       if (icepack_warnings_aborted(subname)) return
     256       48264 :       call column_sum (ncat, esnon, esno_init)
     257       48264 :       if (icepack_warnings_aborted(subname)) return
     258       48264 :       call column_sum (ncat, sicen, sice_init)
     259       48264 :       if (icepack_warnings_aborted(subname)) return
     260       48264 :       call column_sum (ncat, vbrin, vbri_init)
     261       48264 :       if (icepack_warnings_aborted(subname)) return
     262             : 
     263             :       endif ! conserv_check
     264             : 
     265             :       !-----------------------------------------------------------------
     266             :       ! Initialize remapping flag.
     267             :       ! Remapping is done wherever remap_flag = .true.
     268             :       ! In rare cases the category boundaries may shift too far for the
     269             :       !  remapping algorithm to work, and remap_flag is set to .false.
     270             :       ! In these cases the simpler 'rebin' subroutine will shift ice
     271             :       !  between categories if needed.
     272             :       !-----------------------------------------------------------------
     273             : 
     274      852950 :       remap_flag = .true.
     275             : 
     276             :       !-----------------------------------------------------------------
     277             :       ! Compute thickness change in each category.
     278             :       !-----------------------------------------------------------------
     279             : 
     280     5117700 :       do n = 1, ncat
     281             : 
     282     4264750 :          if (aicen_init(n) > puny) then
     283     3936902 :              hicen_init(n) = vicen_init(n) / aicen_init(n)
     284             :          else
     285      327848 :              hicen_init(n) = c0
     286             :          endif               ! aicen_init > puny
     287             : 
     288     5117700 :          if (aicen (n) > puny) then
     289     3936895 :              hicen (n) = vicen(n) / aicen(n) 
     290     3936895 :              dhicen(n) = hicen(n) - hicen_init(n)
     291             :          else
     292      327855 :              hicen (n) = c0
     293      327855 :              dhicen(n) = c0
     294             :          endif               ! aicen > puny
     295             : 
     296             :       enddo                     ! n
     297             : 
     298             :       !-----------------------------------------------------------------
     299             :       ! Compute new category boundaries, hbnew, based on changes in
     300             :       ! ice thickness from vertical thermodynamics.
     301             :       !-----------------------------------------------------------------
     302             : 
     303      852950 :       hbnew(0) = hin_max(0)
     304             : 
     305     4264750 :       do n = 1, ncat-1
     306             : 
     307     3411800 :          if (hicen_init(n)   > puny .and. &
     308             :              hicen_init(n+1) > puny) then
     309             :              ! interpolate between adjacent category growth rates
     310     1042014 :              slope = (dhicen(n+1) - dhicen(n)) / &
     311     4124798 :                  (hicen_init(n+1) - hicen_init(n))
     312     1042014 :              hbnew(n) = hin_max(n) + dhicen(n) &
     313     3082784 :                       + slope * (hin_max(n) - hicen_init(n))
     314      329016 :          elseif (hicen_init(n) > puny) then ! hicen_init(n+1)=0
     315      103088 :              hbnew(n) = hin_max(n) + dhicen(n)
     316      225928 :          elseif (hicen_init(n+1) > puny) then ! hicen_init(n)=0
     317      102493 :              hbnew(n) = hin_max(n) + dhicen(n+1)
     318             :          else
     319      123435 :              hbnew(n) = hin_max(n)
     320             :          endif
     321             : 
     322             :       !-----------------------------------------------------------------
     323             :       ! Check that each boundary lies between adjacent values of hicen.
     324             :       ! If not, set remap_flag = .false.
     325             :       ! Write diagnosis outputs if remap_flag was changed to false
     326             :       !-----------------------------------------------------------------
     327             : 
     328     3411800 :          if (aicen(n) > puny .and. hicen(n) >= hbnew(n)) then
     329           2 :             remap_flag = .false.
     330             : 
     331           2 :             if (print_diags) then
     332             :                write(warnstr,*) subname, 'ITD: hicen(n) > hbnew(n)'
     333             :                call icepack_warnings_add(warnstr)
     334             :                write(warnstr,*) subname, 'cat ',n
     335             :                call icepack_warnings_add(warnstr)
     336             :                write(warnstr,*) subname, 'hicen(n) =', hicen(n)
     337             :                call icepack_warnings_add(warnstr)
     338             :                write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
     339             :                call icepack_warnings_add(warnstr)
     340             :             endif
     341             : 
     342     3411798 :          elseif (aicen(n+1) > puny .and. hicen(n+1) <= hbnew(n)) then
     343           0 :             remap_flag = .false.
     344             : 
     345             :             if (print_diags) then
     346             :                write(warnstr,*) subname, 'ITD: hicen(n+1) < hbnew(n)'
     347             :                call icepack_warnings_add(warnstr)
     348             :                write(warnstr,*) subname, 'cat ',n
     349             :                call icepack_warnings_add(warnstr)
     350             :                write(warnstr,*) subname, 'hicen(n+1) =', hicen(n+1)
     351             :                call icepack_warnings_add(warnstr)
     352             :                write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
     353             :                call icepack_warnings_add(warnstr)
     354             :             endif
     355             :          endif
     356             : 
     357             :       !-----------------------------------------------------------------
     358             :       ! Check that hbnew(n) lies between hin_max(n-1) and hin_max(n+1).
     359             :       ! If not, set remap_flag = .false.
     360             :       ! (In principle we could allow this, but it would make the code
     361             :       ! more complicated.)
     362             :       ! Write diagnosis outputs if remap_flag was changed to false
     363             :       !-----------------------------------------------------------------
     364             : 
     365     3411800 :          if (hbnew(n) > hin_max(n+1)) then
     366           0 :             remap_flag = .false.
     367             : 
     368             :             if (print_diags) then
     369             :                write(warnstr,*) subname, 'ITD hbnew(n) > hin_max(n+1)'
     370             :                call icepack_warnings_add(warnstr)
     371             :                write(warnstr,*) subname, 'cat ',n
     372             :                call icepack_warnings_add(warnstr)
     373             :                write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
     374             :                call icepack_warnings_add(warnstr)
     375             :                write(warnstr,*) subname, 'hin_max(n+1) =', hin_max(n+1)
     376             :                call icepack_warnings_add(warnstr)
     377             :             endif
     378             :          endif
     379             : 
     380     4264750 :          if (hbnew(n) < hin_max(n-1)) then
     381           0 :             remap_flag = .false.
     382             : 
     383             :             if (print_diags) then
     384             :                write(warnstr,*) subname, 'ITD: hbnew(n) < hin_max(n-1)'
     385             :                call icepack_warnings_add(warnstr)
     386             :                write(warnstr,*) subname, 'cat ',n
     387             :                call icepack_warnings_add(warnstr)
     388             :                write(warnstr,*) subname, 'hbnew(n) =', hbnew(n)
     389             :                call icepack_warnings_add(warnstr)
     390             :                write(warnstr,*) subname, 'hin_max(n-1) =', hin_max(n-1)
     391             :                call icepack_warnings_add(warnstr)
     392             :             endif
     393             :          endif
     394             : 
     395             :       enddo                     ! boundaries, 1 to ncat-1
     396             : 
     397             :       !-----------------------------------------------------------------
     398             :       ! Compute hbnew(ncat)
     399             :       !-----------------------------------------------------------------
     400             : 
     401      852950 :       if (aicen(ncat) > puny) then
     402      751030 :          hbnew(ncat) = c3*hicen(ncat) - c2*hbnew(ncat-1)
     403             :       else
     404      101920 :          hbnew(ncat) = hin_max(ncat)
     405             :       endif
     406      852950 :       hbnew(ncat) = max(hbnew(ncat),hin_max(ncat-1))
     407             : 
     408             :       !-----------------------------------------------------------------
     409             :       ! Identify cells where the ITD is to be remapped
     410             :       !-----------------------------------------------------------------
     411             : 
     412      852950 :       if (remap_flag) then
     413             : 
     414             :       !-----------------------------------------------------------------
     415             :       ! Compute g(h) for category 1 at start of time step
     416             :       ! (hicen = hicen_init)
     417             :       !-----------------------------------------------------------------
     418             : 
     419      291771 :          call fit_line(aicen(1),   hicen_init(1), &
     420      291771 :                        hbnew(0),   hin_max   (1), &
     421      291771 :                        g0   (1),   g1        (1), &
     422      852948 :                        hL   (1),   hR        (1))
     423      852948 :          if (icepack_warnings_aborted(subname)) return
     424             : 
     425             :       !-----------------------------------------------------------------
     426             :       ! Find area lost due to melting of thin (category 1) ice
     427             :       !-----------------------------------------------------------------
     428             : 
     429      852948 :          if (aicen(1) > puny) then
     430             : 
     431      751618 :             dh0 = dhicen(1)
     432      751618 :             if (dh0 < c0) then   ! remove area from category 1
     433      230966 :                dh0 = min(-dh0,hin_max(1))   ! dh0 --> |dh0|
     434             : 
     435             :       !-----------------------------------------------------------------
     436             :       ! Integrate g(1) from 0 to dh0 to estimate area melted
     437             :       !-----------------------------------------------------------------
     438             : 
     439             :                ! right integration limit (left limit = 0)
     440      230966 :                etamax = min(dh0,hR(1)) - hL(1)
     441             : 
     442      230966 :                if (etamax > c0) then
     443      181238 :                   x1 = etamax
     444      181238 :                   x2 = p5 * etamax*etamax
     445      181238 :                   da0 = g1(1)*x2 + g0(1)*x1 ! ice area removed
     446             : 
     447             :                ! constrain new thickness <= hicen_init
     448      181238 :                   damax = aicen(1) * (c1-hicen(1)/hicen_init(1)) ! damax > 0
     449      181238 :                   da0 = min (da0, damax)
     450             : 
     451             :                ! remove area, conserving volume
     452      181238 :                   hicen(1) = hicen(1) * aicen(1) / (aicen(1)-da0)
     453      181238 :                   aicen(1) = aicen(1) - da0
     454             : 
     455      181238 :                   if (tr_pond_topo) &
     456       14468 :                      fpond = fpond - (da0 * trcrn(nt_apnd,1) & 
     457       50232 :                                           * trcrn(nt_hpnd,1))
     458             : 
     459             :                endif            ! etamax > 0
     460             : 
     461             :             else                ! dh0 >= 0
     462      520652 :                hbnew(0) = min(dh0,hin_max(1))  ! shift hbnew(0) to right
     463             :             endif
     464             : 
     465             :          endif                  ! aicen(1) > puny
     466             : 
     467             :       !-----------------------------------------------------------------
     468             :       ! Compute g(h) for each ice thickness category.
     469             :       !-----------------------------------------------------------------
     470             : 
     471     5117688 :          do n = 1, ncat
     472             : 
     473     1458855 :             call fit_line(aicen(n),   hicen(n), &
     474     1458855 :                           hbnew(n-1), hbnew(n), &
     475     1458855 :                           g0   (n),   g1   (n), &
     476     5723595 :                           hL   (n),   hR   (n))
     477     4264740 :             if (icepack_warnings_aborted(subname)) return
     478             : 
     479             :       !-----------------------------------------------------------------
     480             :       ! Compute area and volume to be shifted across each boundary.
     481             :       !-----------------------------------------------------------------
     482             : 
     483     4264740 :             donor(n) = 0
     484     4264740 :             daice(n) = c0
     485     6576543 :             dvice(n) = c0
     486             :          enddo
     487             : 
     488     4264740 :          do n = 1, ncat-1
     489             : 
     490     3411792 :             if (hbnew(n) > hin_max(n)) then ! transfer from n to n+1
     491             : 
     492             :                ! left and right integration limits in eta space
     493     1935580 :                etamin = max(hin_max(n), hL(n)) - hL(n)
     494     1935580 :                etamax = min(hbnew(n),   hR(n)) - hL(n)
     495     1935580 :                donor(n) = n
     496             : 
     497             :             else             ! hbnew(n) <= hin_max(n); transfer from n+1 to n
     498             : 
     499             :                ! left and right integration limits in eta space
     500     1476212 :                etamin = c0
     501     1476212 :                etamax = min(hin_max(n), hR(n+1)) - hL(n+1)
     502     1476212 :                donor(n) = n+1
     503             : 
     504             :             endif            ! hbnew(n) > hin_max(n)
     505             : 
     506     3411792 :             if (etamax > etamin) then
     507     2755027 :                x1  = etamax - etamin
     508     2755027 :                wk1 = etamin*etamin
     509     2755027 :                wk2 = etamax*etamax
     510     2755027 :                x2  = p5 * (wk2 - wk1)
     511     2755027 :                wk1 = wk1*etamin
     512     2755027 :                wk2 = wk2*etamax
     513     2755027 :                x3  = p333 * (wk2 - wk1)
     514     2755027 :                nd  = donor(n)
     515     2755027 :                daice(n) = g1(nd)*x2 + g0(nd)*x1
     516     2755027 :                dvice(n) = g1(nd)*x3 + g0(nd)*x2 + daice(n)*hL(nd)
     517             :             endif               ! etamax > etamin
     518             : 
     519             :             ! If daice or dvice is very small, shift no ice.
     520             : 
     521     3411792 :             nd = donor(n)
     522             : 
     523     3411792 :             if (daice(n) < aicen(nd)*puny) then
     524      424936 :                daice(n) = c0
     525      424936 :                dvice(n) = c0
     526      424936 :                donor(n) = 0
     527             :             endif 
     528             : 
     529     3411792 :             if (dvice(n) < vicen(nd)*puny) then
     530      424936 :                daice(n) = c0
     531      424936 :                dvice(n) = c0
     532      424936 :                donor(n) = 0
     533             :             endif
     534             : 
     535             :             ! If daice is close to aicen or dvice is close to vicen,
     536             :             ! shift entire category
     537             : 
     538     3411792 :             if (daice(n) > aicen(nd)*(c1-puny)) then
     539        3768 :                daice(n) = aicen(nd)
     540        3768 :                dvice(n) = vicen(nd)
     541             :             endif
     542             : 
     543     4264740 :             if (dvice(n) > vicen(nd)*(c1-puny)) then
     544        3768 :                daice(n) = aicen(nd)
     545        3768 :                dvice(n) = vicen(nd)
     546             :             endif
     547             : 
     548             :          enddo                     ! boundaries, 1 to ncat-1
     549             : 
     550             :       !-----------------------------------------------------------------
     551             :       ! Shift ice between categories as necessary
     552             :       !-----------------------------------------------------------------
     553             : 
     554             :          ! maintain qsno negative definiteness
     555     5117688 :          do n = 1, ncat
     556    10347668 :             do k = nt_qsno, nt_qsno+nslyr-1
     557     9494720 :                trcrn(k,n) = trcrn(k,n) + rhos*Lfresh
     558             :             enddo
     559             :          enddo
     560             :  
     561             :          call shift_ice (ntrcr,    ncat,        &
     562           0 :                          trcr_depend,           &
     563           0 :                          trcr_base,             &
     564           0 :                          n_trcr_strata,         &
     565           0 :                          nt_strata,             &
     566           0 :                          aicen,    trcrn,       &
     567           0 :                          vicen,    vsnon,       &
     568           0 :                          hicen,    donor,       &
     569      852948 :                          daice,    dvice        )
     570      852948 :          if (icepack_warnings_aborted(subname)) return
     571             : 
     572             :          ! maintain qsno negative definiteness
     573     5117688 :          do n = 1, ncat
     574    10347668 :             do k = nt_qsno, nt_qsno+nslyr-1
     575     9494720 :                trcrn(k,n) = trcrn(k,n) - rhos*Lfresh
     576             :             enddo
     577             :          enddo
     578             : 
     579             :       !-----------------------------------------------------------------
     580             :       ! Make sure hice(1) >= minimum ice thickness hi_min.
     581             :       !-----------------------------------------------------------------
     582             : 
     583      852948 :          if (hi_min > c0 .and. aicen(1) > puny .and. hicen(1) < hi_min) then
     584             : 
     585         738 :             da0 = aicen(1) * (c1 - hicen(1)/hi_min)
     586         738 :             aicen(1) = aicen(1) - da0
     587         738 :             hicen(1) = hi_min
     588             : 
     589         738 :             if (tr_pond_topo) &
     590           0 :                fpond = fpond - (da0 * trcrn(nt_apnd,1) & 
     591         139 :                                     * trcrn(nt_hpnd,1))
     592             :          endif
     593             : 
     594             :       endif ! remap_flag
     595             : 
     596             :       !-----------------------------------------------------------------
     597             :       ! Update fractional ice area in each grid cell.
     598             :       !-----------------------------------------------------------------
     599             : 
     600      852950 :       call aggregate_area (ncat, aicen, aice, aice0)
     601      852950 :       if (icepack_warnings_aborted(subname)) return
     602             : 
     603             :       !-----------------------------------------------------------------
     604             :       ! Check volume and energy conservation.
     605             :       !-----------------------------------------------------------------
     606             : 
     607      852950 :       if (conserv_check) then
     608             : 
     609      289584 :       do n = 1, ncat
     610             : 
     611      241320 :          eicen(n) = c0
     612      241320 :          esnon(n) = c0
     613      241320 :          vbrin(n) = c0
     614      241320 :          sicen(n) = c0
     615             : 
     616     1930560 :          do k = 1, nilyr
     617     1226400 :             eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
     618     2543760 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     619             :          enddo
     620     1447920 :          do k = 1, nslyr
     621      876000 :             esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
     622     1885920 :                      * vsnon(n)/real(nslyr,kind=dbl_kind)
     623             :          enddo
     624             : 
     625      241320 :          if (tr_brine) then
     626           0 :             vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
     627           0 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     628             :          endif
     629             : 
     630     1978824 :          do k = 1, nilyr
     631     1226400 :             sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
     632     2543760 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     633             :          enddo
     634             : 
     635             :       enddo ! n
     636             : 
     637       48264 :       call column_sum (ncat, vicen, vice_final)
     638       48264 :       if (icepack_warnings_aborted(subname)) return
     639       48264 :       call column_sum (ncat, vsnon, vsno_final)
     640       48264 :       if (icepack_warnings_aborted(subname)) return
     641       48264 :       call column_sum (ncat, eicen, eice_final)
     642       48264 :       if (icepack_warnings_aborted(subname)) return
     643       48264 :       call column_sum (ncat, esnon, esno_final)
     644       48264 :       if (icepack_warnings_aborted(subname)) return
     645       48264 :       call column_sum (ncat, sicen, sice_final)
     646       48264 :       if (icepack_warnings_aborted(subname)) return
     647       48264 :       call column_sum (ncat, vbrin, vbri_final)
     648       48264 :       if (icepack_warnings_aborted(subname)) return
     649             : 
     650       48264 :       fieldid = subname//':vice'
     651             :       call column_conservation_check (fieldid,               &
     652             :                                       vice_init, vice_final, &
     653       48264 :                                       puny)
     654       48264 :       if (icepack_warnings_aborted(subname)) return
     655       48264 :       fieldid = subname//':vsno'
     656             :       call column_conservation_check (fieldid,               &
     657             :                                       vsno_init, vsno_final, &
     658       48264 :                                       puny)
     659       48264 :       if (icepack_warnings_aborted(subname)) return
     660       48264 :       fieldid = subname//':eice'
     661             :       call column_conservation_check (fieldid,               &
     662             :                                       eice_init, eice_final, &
     663       48264 :                                       puny*Lfresh*rhoi)
     664       48264 :       if (icepack_warnings_aborted(subname)) return
     665       48264 :       fieldid = subname//':esno'
     666             :       call column_conservation_check (fieldid,               &
     667             :                                       esno_init, esno_final, &
     668       48264 :                                       puny*Lfresh*rhos)
     669       48264 :       if (icepack_warnings_aborted(subname)) return
     670       48264 :       fieldid = subname//':sicen'
     671             :       call column_conservation_check (fieldid,               &
     672             :                                       sice_init, sice_final, &
     673       48264 :                                       puny)
     674       48264 :       if (icepack_warnings_aborted(subname)) return
     675       48264 :       fieldid = subname//':vbrin'
     676             :       call column_conservation_check (fieldid,               &
     677             :                                       vbri_init, vbri_final, &
     678       48264 :                                       puny*c10)
     679       48264 :       if (icepack_warnings_aborted(subname)) return
     680             : 
     681             :       endif                     ! conservation check
     682             : 
     683             :       end subroutine linear_itd
     684             : 
     685             : !=======================================================================
     686             : !
     687             : ! Fit g(h) with a line, satisfying area and volume constraints.
     688             : ! To reduce roundoff errors caused by large values of g0 and g1,
     689             : ! we actually compute g(eta), where eta = h - hL, and hL is the
     690             : ! left boundary.
     691             : !
     692             : ! authors: William H. Lipscomb, LANL
     693             : !          Elizabeth C. Hunke, LANL
     694             : 
     695     5117688 :       subroutine fit_line (aicen,    hice,            &
     696             :                            hbL,      hbR,             &
     697             :                            g0,       g1,              &
     698             :                            hL,       hR)
     699             : 
     700             :       real (kind=dbl_kind), intent(in) :: &
     701             :          aicen           ! concentration of ice
     702             : 
     703             :       real (kind=dbl_kind), intent(in) :: &
     704             :          hbL, hbR    , & ! left and right category boundaries
     705             :          hice            ! ice thickness
     706             : 
     707             :       real (kind=dbl_kind), intent(out):: &
     708             :          g0, g1      , & ! coefficients in linear equation for g(eta)
     709             :          hL          , & ! min value of range over which g(h) > 0
     710             :          hR              ! max value of range over which g(h) > 0
     711             : 
     712             :       ! local variables
     713             : 
     714             :       real  (kind=dbl_kind) :: &
     715     1750626 :          h13         , & ! hbL + 1/3 * (hbR - hbL)
     716     1750626 :          h23         , & ! hbL + 2/3 * (hbR - hbL)
     717     1750626 :          dhr         , & ! 1 / (hR - hL)
     718     1750626 :          wk1, wk2        ! temporary variables
     719             : 
     720             :       character(len=*),parameter :: subname='(fit_line)'
     721             : 
     722             :       !-----------------------------------------------------------------
     723             :       ! Compute g0, g1, hL, and hR for each category to be remapped.
     724             :       !-----------------------------------------------------------------
     725             : 
     726     5117688 :          if (aicen > puny .and. hbR - hbL > puny) then
     727             : 
     728             :          ! Initialize hL and hR
     729             : 
     730     4688511 :             hL = hbL
     731     4688511 :             hR = hbR
     732             : 
     733             :          ! Change hL or hR if hicen(n) falls outside central third of range
     734             : 
     735     4688511 :             h13 = p333 * (c2*hL + hR)
     736     4688511 :             h23 = p333 * (hL + c2*hR)
     737     4688511 :             if (hice < h13) then
     738      662512 :                hR = c3*hice - c2*hL
     739     4025999 :             elseif (hice > h23) then
     740      579965 :                hL = c3*hice - c2*hR
     741             :             endif
     742             : 
     743             :          ! Compute coefficients of g(eta) = g0 + g1*eta
     744             : 
     745     4688511 :             dhr = c1 / (hR - hL)
     746     4688511 :             wk1 = c6 * aicen * dhr
     747     4688511 :             wk2 = (hice - hL) * dhr
     748     4688511 :             g0 = wk1 * (p666 - wk2)
     749     4688511 :             g1 = c2*dhr * wk1 * (wk2 - p5)
     750             : 
     751             :          else
     752             : 
     753      429177 :             g0 = c0
     754      429177 :             g1 = c0
     755      429177 :             hL = c0
     756      429177 :             hR = c0
     757             : 
     758             :          endif                  ! aicen > puny
     759             : 
     760     5117688 :       end subroutine fit_line
     761             : 
     762             : !=======================================================================
     763             : !
     764             : ! Given some added new ice to the base of the existing ice, recalculate 
     765             : ! vertical tracer so that new grid cells are all the same size. 
     766             : !
     767             : ! author: A. K. Turner, LANL
     768             : !
     769       93904 :       subroutine update_vertical_tracers(nilyr, trc, h1, h2, trc0)
     770             : 
     771             :       integer (kind=int_kind), intent(in) :: &
     772             :          nilyr ! number of ice layers
     773             : 
     774             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     775             :            trc ! vertical tracer
     776             : 
     777             :       real (kind=dbl_kind), intent(in) :: &
     778             :          h1, & ! old thickness
     779             :          h2, & ! new thickness
     780             :          trc0  ! tracer value of added ice on ice bottom
     781             :            
     782             :       ! local variables
     783             : 
     784      452480 :       real(kind=dbl_kind), dimension(nilyr) :: trc2 ! updated tracer temporary
     785             : 
     786             :       ! vertical indices for old and new grid
     787             :       integer :: k1, k2
     788             : 
     789             :       real (kind=dbl_kind) :: &
     790       44112 :          z1a, z1b, & ! upper, lower boundary of old cell/added new ice at bottom
     791       44112 :          z2a, z2b, & ! upper, lower boundary of new cell
     792       44112 :          overlap , & ! overlap between old and new cell
     793       44112 :          rnilyr
     794             : 
     795             :       character(len=*),parameter :: subname='(update_vertical_tracers)'
     796             : 
     797       93904 :         rnilyr = real(nilyr,dbl_kind)
     798             : 
     799             :         ! loop over new grid cells
     800      751232 :         do k2 = 1, nilyr
     801             : 
     802             :            ! initialize new tracer
     803      657328 :            trc2(k2) = c0
     804             : 
     805             :            ! calculate upper and lower boundary of new cell
     806      657328 :            z2a = ((k2 - 1) * h2) / rnilyr
     807      657328 :            z2b = (k2       * h2) / rnilyr
     808             : 
     809             :            ! loop over old grid cells
     810     5258624 :            do k1 = 1, nilyr
     811             : 
     812             :               ! calculate upper and lower boundary of old cell
     813     4601296 :               z1a = ((k1 - 1) * h1) / rnilyr
     814     4601296 :               z1b = (k1       * h1) / rnilyr
     815             :               
     816             :               ! calculate overlap between old and new cell
     817     4601296 :               overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
     818             : 
     819             :               ! aggregate old grid cell contribution to new cell
     820     5258624 :               trc2(k2) = trc2(k2) + overlap * trc(k1)
     821             : 
     822             :            enddo ! k1
     823             : 
     824             :            ! calculate upper and lower boundary of added new ice at bottom
     825      657328 :            z1a = h1
     826      657328 :            z1b = h2
     827             :            
     828             :            ! calculate overlap between added ice and new cell
     829      657328 :            overlap = max(min(z1b, z2b) - max(z1a, z2a), c0)
     830             :            ! aggregate added ice contribution to new cell
     831      657328 :            trc2(k2) = trc2(k2) + overlap * trc0
     832             :            ! renormalize new grid cell
     833      751232 :            trc2(k2) = (rnilyr * trc2(k2)) / h2
     834             : 
     835             :         enddo ! k2
     836             : 
     837             :         ! update vertical tracer array with the adjusted tracer
     838      751232 :         trc = trc2
     839             : 
     840       93904 :       end subroutine update_vertical_tracers
     841             : 
     842             : !=======================================================================
     843             : !
     844             : ! Given the fraction of ice melting laterally in each grid cell
     845             : !  (computed in subroutine frzmlt_bottom_lateral), melt ice.
     846             : !
     847             : ! author: C. M. Bitz, UW
     848             : ! 2003:   Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL
     849             : ! 2016    Lettie Roach, NIWA/VUW, added floe size dependence
     850             : !
     851     1027620 :       subroutine lateral_melt (dt,         ncat,       &
     852             :                                nilyr,      nslyr,      &
     853             :                                n_aero,     &
     854             :                                fpond,      &
     855             :                                fresh,      fsalt,      &
     856     1027620 :                                fhocn,      faero_ocn,  &
     857     1027620 :                                fiso_ocn,               &
     858             :                                rside,      meltl,      &
     859             :                                fside,      sss,        &
     860     2055240 :                                aicen,      vicen,      &
     861     2055240 :                                vsnon,      trcrn,      &
     862     1027620 :                                fzsal,      flux_bio,   &
     863             :                                nbtrcr,     nblyr,      &
     864     1027620 :                                nfsd,       d_afsd_latm,&
     865     1027620 :                                floe_rad_c, floe_binwidth)
     866             : 
     867             :       real (kind=dbl_kind), intent(in) :: &
     868             :          dt        ! time step (s)
     869             : 
     870             :       integer (kind=int_kind), intent(in) :: &
     871             :          ncat    , & ! number of thickness categories
     872             :          nfsd    , & ! number of floe size categories
     873             :          nilyr   , & ! number of ice layers
     874             :          nblyr   , & ! number of bio layers
     875             :          nslyr   , & ! number of snow layers
     876             :          n_aero  , & ! number of aerosol tracers
     877             :          nbtrcr      ! number of bio tracers
     878             : 
     879             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     880             :          aicen   , & ! concentration of ice
     881             :          vicen   , & ! volume per unit area of ice          (m)
     882             :          vsnon       ! volume per unit area of snow         (m)
     883             : 
     884             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
     885             :          trcrn       ! tracer array
     886             : 
     887             :       real (kind=dbl_kind), intent(in) :: &
     888             :          rside   , & ! fraction of ice that melts laterally
     889             :          fside       ! lateral heat flux (W/m^2)
     890             : 
     891             :       real (kind=dbl_kind), intent(inout) :: &
     892             :          fpond     , & ! fresh water flux to ponds (kg/m^2/s)
     893             :          fresh     , & ! fresh water flux to ocean (kg/m^2/s)
     894             :          fsalt     , & ! salt flux to ocean (kg/m^2/s)
     895             :          fhocn     , & ! net heat flux to ocean (W/m^2)
     896             :          meltl     , & ! lateral ice melt         (m/step-->cm/day)
     897             :          fzsal         ! salt flux from zsalinity (kg/m2/s)
     898             :   
     899             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     900             :          floe_rad_c     , & ! fsd size bin centre in m (radius)
     901             :          floe_binwidth      ! fsd size bin width in m (radius)
     902             : 
     903             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
     904             :          d_afsd_latm        ! change in fsd due to lateral melt (m)
     905             : 
     906             :       real (kind=dbl_kind), dimension(nbtrcr), &
     907             :          intent(inout) :: &
     908             :          flux_bio  ! biology tracer flux from layer bgc (mmol/m^2/s)  
     909             : 
     910             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     911             :          faero_ocn     ! aerosol flux to ocean (kg/m^2/s)
     912             : 
     913             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     914             :          fiso_ocn     ! isotope flux to ocean (kg/m^2/s)
     915             : 
     916             :       ! local variables
     917             : 
     918             :       integer (kind=int_kind) :: &
     919             :          n           , & ! thickness category index
     920             :          k           , & ! layer index
     921             :          nsubt           ! sub timesteps for FSD tendency
     922             : 
     923             :       real (kind=dbl_kind) :: &
     924      355716 :          dfhocn  , & ! change in fhocn
     925      355716 :          dfpond  , & ! change in fpond
     926      355716 :          dfresh  , & ! change in fresh
     927      355716 :          dfsalt  , & ! change in fsalt
     928      355716 :          dvssl   , & ! snow surface layer volume
     929      355716 :          dvint   , & ! snow interior layer
     930      711432 :          cat1_arealoss, tmp !
     931             : 
     932             :       logical (kind=log_kind) :: &
     933             :          flag        ! .true. if there could be lateral melting
     934             : 
     935             :       real (kind=dbl_kind), dimension (ncat) :: &
     936     3372984 :          aicen_init, & ! initial area fraction
     937     3372984 :          vicen_init, & ! volume per unit area of ice (m)
     938     3372984 :          G_radialn , & ! rate of lateral melt (m/s)
     939     3372984 :          delta_an  , & ! change in the ITD
     940     3372984 :          qin       , & ! enthalpy
     941     3372984 :          rsiden        ! delta_an/aicen
     942             : 
     943             :       real (kind=dbl_kind), dimension (nfsd,ncat) :: &
     944     6491844 :          afsdn     , & ! floe size distribution tracer
     945     6491844 :          afsdn_init    ! initial value
     946             : 
     947             :       real (kind=dbl_kind), dimension (nfsd) :: &
     948     2344320 :          df_flx, &        ! finite difference for FSD
     949     4056264 :          afsd_tmp, d_afsd_tmp
     950             : 
     951             :       real (kind=dbl_kind), dimension(nfsd+1) :: &
     952     2739564 :          f_flx         !
     953             : 
     954             : !echmod - for average qin
     955             :       real (kind=dbl_kind), intent(in) :: &
     956             :          sss
     957             :       real (kind=dbl_kind) :: &
     958      355716 :          Ti, Si0, qi0, &
     959      355716 :          elapsed_t,    & ! FSD subcycling
     960      355716 :          subdt           ! FSD timestep (s)
     961             : 
     962             :       character(len=*), parameter :: subname='(lateral_melt)'
     963             : 
     964     1027620 :       flag = .false.
     965     1027620 :       dfhocn   = c0
     966     1027620 :       dfpond   = c0
     967     1027620 :       dfresh   = c0
     968     1027620 :       dfsalt   = c0
     969     1027620 :       dvssl    = c0
     970     1027620 :       dvint    = c0
     971     1027620 :       cat1_arealoss  = c0
     972     1027620 :       tmp  = c0
     973     5876136 :       vicen_init = c0
     974     5876136 :       G_radialn  = c0
     975     5876136 :       delta_an   = c0
     976     5876136 :       qin        = c0
     977     5876136 :       rsiden     = c0
     978    14706432 :       afsdn      = c1
     979    14706432 :       afsdn_init = c0
     980     2851596 :       df_flx     = c0
     981     3879216 :       f_flx      = c0
     982             : 
     983     1027620 :       if (tr_fsd) then
     984       98676 :          call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:)) 
     985       98676 :          if (icepack_warnings_aborted(subname)) return
     986             :          
     987     5067216 :          afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:)
     988      592056 :          aicen_init = aicen
     989     5067216 :          afsdn_init = afsdn ! for diagnostics
     990      993708 :          d_afsd_latm(:) = c0
     991             :       end if
     992             : 
     993     1027620 :       if (tr_fsd .and. fside < c0) then
     994       33907 :          flag = .true.
     995             : 
     996             : 
     997             : ! echmod - using category values would be preferable to the average value
     998             :          ! Compute average enthalpy of ice (taken from add_new_ice)
     999       33907 :          if (sss > c2 * dSin0_frazil) then
    1000       33907 :             Si0 = sss - dSin0_frazil
    1001             :          else
    1002           0 :             Si0 = sss**2 / (c4*dSin0_frazil)
    1003             :          endif
    1004       33907 :          Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
    1005       33907 :          qi0 = enthalpy_mush(Ti, Si0)
    1006             : 
    1007      237349 :          do n = 1, ncat
    1008      169535 :             if (ktherm == 2) then  ! mushy
    1009     1356280 :                do k = 1, nilyr
    1010             :                   !qin(n) = qin(n) &
    1011             :                   !       + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind)
    1012     1356280 :                   qin(n) = qi0
    1013             :                enddo
    1014             :             else
    1015           0 :                qin(n) = -rhoi*Lfresh
    1016             :             endif
    1017             : 
    1018      169535 :             if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative
    1019             : 
    1020      203442 :             if (G_radialn(n) < -puny) then
    1021             : 
    1022             :                
    1023     1509815 :                if (any(afsdn(:,n) < c0)) print*,&
    1024           0 :                  'lateral_melt B afsd < 0',n
    1025             : 
    1026      210850 :                cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt &
    1027      269440 :                              * G_radialn(n) / floe_binwidth(1)
    1028             : 
    1029      164015 :                delta_an(n) = c0
    1030     1509815 :                do k = 1, nfsd
    1031      642720 :                   delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) &
    1032     1509815 :                                             * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0
    1033             :                end do
    1034             : 
    1035             :                ! add negative area loss from fsd
    1036      164015 :                delta_an(n) = delta_an(n) - cat1_arealoss
    1037             : 
    1038      164015 :                if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)
    1039             :  
    1040             :                ! following original code, not necessary for fsd
    1041      164015 :                if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)
    1042             : 
    1043      164015 :                if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)
    1044             : 
    1045             :             end if ! G_radialn
    1046             :          enddo ! ncat
    1047             : 
    1048      993713 :       else if (rside > c0) then ! original, non-fsd implementation
    1049             : 
    1050      465887 :          flag = .true.
    1051     2793194 :          rsiden(:) = rside ! initialize
    1052             : 
    1053             :       endif
    1054             : 
    1055     1027620 :       if (flag) then ! grid cells with lateral melting.
    1056             : 
    1057     2996636 :          do n = 1, ncat
    1058             : 
    1059             :       !-----------------------------------------------------------------
    1060             :       ! Melt the ice and increment fluxes.
    1061             :       !-----------------------------------------------------------------
    1062             : 
    1063             :             ! fluxes to coupler
    1064             :             ! dfresh > 0, dfsalt > 0, dfpond > 0
    1065             : 
    1066     2496842 :             dfresh = (rhoi*vicen(n) + rhos*vsnon(n))      * rsiden(n) / dt
    1067     2496842 :             dfsalt =  rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
    1068     2496842 :             fresh  = fresh + dfresh
    1069     2496842 :             fsalt  = fsalt + dfsalt
    1070             : 
    1071     2496842 :             if (tr_pond_topo) then
    1072      329895 :                dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n)
    1073      329895 :                fpond  = fpond - dfpond
    1074             :             endif
    1075             : 
    1076             :             ! history diagnostics
    1077     2496842 :             meltl = meltl + vicen(n)*rsiden(n)
    1078             : 
    1079             :             ! state variables
    1080     2496842 :             vicen_init(n) = vicen(n)
    1081     2496842 :             aicen(n) = aicen(n) * (c1 - rsiden(n))
    1082     2496842 :             vicen(n) = vicen(n) * (c1 - rsiden(n))
    1083     2496842 :             vsnon(n) = vsnon(n) * (c1 - rsiden(n))
    1084             : 
    1085             :             ! floe size distribution
    1086     2496842 :             if (tr_fsd) then
    1087      169535 :                if (rsiden(n) > puny) then
    1088      163187 :                   if (aicen(n) > puny) then
    1089             : 
    1090             :                      ! adaptive subtimestep
    1091      163172 :                      elapsed_t = c0
    1092     1501166 :                      afsd_tmp(:) = afsdn_init(:,n)
    1093     1501166 :                      d_afsd_tmp(:) = c0
    1094      163172 :                      nsubt = 0
    1095             : 
    1096      326344 :                      DO WHILE (elapsed_t.lt.dt)
    1097             : 
    1098      163172 :                          nsubt = nsubt + 1
    1099      163172 :                          if (nsubt.gt.100) &
    1100           0 :                            print *, 'latm not converging'
    1101             :                      
    1102             :                          ! finite differences
    1103     1501166 :                          df_flx(:) = c0
    1104     1664338 :                          f_flx (:) = c0
    1105     1337994 :                          do k = 2, nfsd
    1106     1337994 :                            f_flx(k) =  G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
    1107             :                          end do
    1108             : 
    1109     1501166 :                          do k = 1, nfsd
    1110     1501166 :                           df_flx(k)   = f_flx(k+1) - f_flx(k) 
    1111             :                          end do
    1112             : 
    1113     1501166 :                          if (abs(sum(df_flx(:))) > puny) &
    1114           0 :                            print*,'sum(df_flx)/=0'
    1115             : 
    1116             :                          ! this term ensures area conservation
    1117     1501166 :                          tmp = SUM(afsd_tmp(:)/floe_rad_c(:))
    1118             :                         
    1119             :                          ! fsd tendency
    1120     1501166 :                          do k = 1, nfsd
    1121      639858 :                            d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
    1122     1501166 :                                        * (c1/floe_rad_c(k) - tmp)
    1123             :                          end do
    1124             : 
    1125             :                          ! timestep required for this
    1126      163172 :                          subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
    1127      163172 :                          subdt = MIN(subdt, dt)
    1128             : 
    1129             :                         ! update fsd and elapsed time
    1130     1501166 :                         afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
    1131      326344 :                         elapsed_t = elapsed_t + subdt
    1132             : 
    1133             : 
    1134             :                       END DO
    1135             :  
    1136     1501166 :                      afsdn(:,n) = afsd_tmp(:)
    1137             : 
    1138             :       
    1139             :                   end if ! aicen
    1140             :                end if ! rside > 0, otherwise do nothing
    1141             : 
    1142             :             end if ! tr_fsd
    1143             : 
    1144             :             ! fluxes
    1145    18688354 :             do k = 1, nilyr
    1146             :                ! enthalpy tracers do not change (e/v constant)
    1147             :                ! heat flux to coupler for ice melt (dfhocn < 0)
    1148    11891242 :                dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt &
    1149    22137133 :                       * vicen(n)/real(nilyr,kind=dbl_kind)
    1150    18688354 :                fhocn  = fhocn + dfhocn
    1151             :             enddo                  ! nilyr
    1152             : 
    1153     5958884 :             do k = 1, nslyr
    1154             :                ! heat flux to coupler for snow melt (dfhocn < 0)
    1155     2542802 :                dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt &
    1156     4733443 :                       * vsnon(n)/real(nslyr,kind=dbl_kind)
    1157     5958884 :                fhocn  = fhocn + dfhocn
    1158             :             enddo                  ! nslyr
    1159             : 
    1160     2496842 :             if (tr_aero) then
    1161      177350 :                do k = 1, n_aero
    1162           0 :                   faero_ocn(k) = faero_ocn(k) + (vsnon(n) &
    1163           0 :                                *(trcrn(nt_aero  +4*(k-1),n)   &
    1164           0 :                                + trcrn(nt_aero+1+4*(k-1),n))  &
    1165           0 :                                               +  vicen(n) &
    1166           0 :                                *(trcrn(nt_aero+2+4*(k-1),n)   &
    1167           0 :                                + trcrn(nt_aero+3+4*(k-1),n))) &
    1168      177350 :                                * rsiden(n) / dt
    1169             :                enddo ! k
    1170             :             endif    ! tr_aero
    1171             : 
    1172     2496842 :             if (tr_iso) then
    1173      350660 :                do k = 1, n_iso
    1174           0 :                   fiso_ocn(k) = fiso_ocn(k) &
    1175           0 :                               + (vsnon(n)*trcrn(nt_isosno+k-1,n) &
    1176           0 :                               +  vicen(n)*trcrn(nt_isoice+k-1,n)) &
    1177      350660 :                               * rside / dt
    1178             :                enddo ! k
    1179             :             endif    ! tr_iso
    1180             : 
    1181             :       !-----------------------------------------------------------------
    1182             :       ! Biogeochemistry
    1183             :       !-----------------------------------------------------------------     
    1184             : 
    1185     2996636 :             if (z_tracers) then   ! snow tracers
    1186      315035 :                dvssl  = min(p5*vsnon(n), hs_ssl*aicen(n))       !snow surface layer
    1187      315035 :                dvint  = vsnon(n)- dvssl                         !snow interior
    1188     6561395 :                do k = 1, nbtrcr
    1189     1647660 :                   flux_bio(k) = flux_bio(k) &
    1190     3295320 :                               + (trcrn(bio_index(k)+nblyr+1,n)*dvssl  &
    1191     3295320 :                               +  trcrn(bio_index(k)+nblyr+2,n)*dvint) &
    1192     9856715 :                               * rsiden(n) / dt
    1193             :                enddo
    1194             :             endif
    1195             : 
    1196             :          enddo       ! n
    1197             : 
    1198      499794 :          if (solve_zsal .or. z_tracers) &
    1199             :             call lateral_melt_bgc(dt,                         &
    1200             :                                   ncat,        nblyr,         &
    1201           0 :                                   rside,       vicen_init,    &  !echmod: use rsiden
    1202           0 :                                   trcrn,       fzsal,         &
    1203       63007 :                                   flux_bio,    nbtrcr)
    1204      499794 :             if (icepack_warnings_aborted(subname)) return
    1205             : 
    1206             :       endif          ! flag
    1207             : 
    1208     1027620 :       if (tr_fsd) then
    1209             : 
    1210       98676 :          call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
    1211       98676 :          if (icepack_warnings_aborted(subname)) return
    1212             : 
    1213             :          ! diagnostics
    1214      993708 :          do k = 1, nfsd
    1215      895032 :             d_afsd_latm(k) = c0
    1216     5468868 :             do n = 1, ncat
    1217     1708200 :                d_afsd_latm(k) = d_afsd_latm(k) &
    1218     5370192 :                   + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n)
    1219             :             end do
    1220             :          end do
    1221             :       end if
    1222             : 
    1223             :       end subroutine lateral_melt
    1224             : 
    1225             : !=======================================================================
    1226             : !
    1227             : ! Given the volume of new ice grown in open water, compute its area
    1228             : ! and thickness and add it to the appropriate category or categories.
    1229             : !
    1230             : ! NOTE: Usually all the new ice is added to category 1.  An exception is
    1231             : !       made if there is no open water or if the new ice is too thick
    1232             : !       for category 1, in which case ice is distributed evenly over the
    1233             : !       entire cell.  Subroutine rebin should be called in case the ice
    1234             : !       thickness lies outside category bounds after new ice formation.
    1235             : !
    1236             : ! When ice must be added to categories above category 1, the mushy
    1237             : ! formulation (ktherm=2) adds it only to the bottom of the ice.  When
    1238             : ! added to only category 1, all formulations combine the new ice and
    1239             : ! existing ice tracers as bulk quantities.
    1240             : !
    1241             : ! authors: William H. Lipscomb, LANL
    1242             : !          Elizabeth C. Hunke, LANL
    1243             : !          Adrian Turner, LANL
    1244             : !          Lettie Roach, NIWA/VUW
    1245             : !
    1246     1027620 :       subroutine add_new_ice (ncat,      nilyr,      &
    1247             :                               nfsd,      nblyr,      &
    1248             :                               n_aero,    dt,         &
    1249             :                               ntrcr,     nltrcr,     &
    1250     1027620 :                               hin_max,   ktherm,     &
    1251     1027620 :                               aicen,     trcrn,      &
    1252     1027620 :                               vicen,     vsnon1,     &
    1253             :                               aice0,     aice,       &
    1254             :                               frzmlt,    frazil,     &
    1255             :                               frz_onset, yday,       &
    1256             :                               update_ocn_f,          &
    1257             :                               fresh,     fsalt,      &
    1258             :                               Tf,        sss,        &
    1259     1027620 :                               salinz,    phi_init,   &
    1260             :                               dSin0_frazil,          &
    1261     1027620 :                               bgrid,      cgrid,      igrid,    &
    1262     1027620 :                               nbtrcr,    flux_bio,   &
    1263     1027620 :                               ocean_bio, fzsal,      &
    1264             :                               frazil_diag,           &
    1265     1027620 :                               fiso_ocn,              &
    1266             :                               HDO_ocn, H2_16O_ocn,   &
    1267             :                               H2_18O_ocn,            &
    1268             :                               wave_sig_ht,           &
    1269     1027620 :                               wave_spectrum,         &
    1270     1027620 :                               wavefreq,              &
    1271     1027620 :                               dwavefreq,             &
    1272     1027620 :                               d_afsd_latg,           &
    1273     1027620 :                               d_afsd_newi,           &
    1274     2055240 :                               floe_rad_c, floe_binwidth)
    1275             : 
    1276             :       use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice
    1277             : 
    1278             :       integer (kind=int_kind), intent(in) :: &
    1279             :          ncat  , & ! number of thickness categories
    1280             :          nfsd  , & ! number of floe size categories
    1281             :          nilyr , & ! number of ice layers
    1282             :          nblyr , & ! number of bio layers
    1283             :          ntrcr , & ! number of tracers
    1284             :          nltrcr, & ! number of zbgc tracers
    1285             :          n_aero, & ! number of aerosol tracers
    1286             :          ktherm    ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
    1287             : 
    1288             :       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
    1289             :          hin_max      ! category boundaries (m)
    1290             : 
    1291             :       real (kind=dbl_kind), intent(in) :: &
    1292             :          dt    , & ! time step (s)
    1293             :          aice  , & ! total concentration of ice
    1294             :          frzmlt, & ! freezing/melting potential (W/m^2)
    1295             :          Tf    , & ! freezing temperature (C)
    1296             :          sss   , & ! sea surface salinity (ppt)
    1297             :          vsnon1    ! category 1 snow volume per ice area (m)
    1298             : 
    1299             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1300             :          aicen , & ! concentration of ice
    1301             :          vicen     ! volume per unit area of ice          (m)
    1302             : 
    1303             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
    1304             :          trcrn     ! ice tracers
    1305             :                    ! 1: surface temperature
    1306             : 
    1307             :       real (kind=dbl_kind), intent(inout) :: &
    1308             :          aice0     , & ! concentration of open water
    1309             :          frazil    , & ! frazil ice growth        (m/step-->cm/day)
    1310             :          frazil_diag,& ! frazil ice growth diagnostic (m/step-->cm/day)
    1311             :          fresh     , & ! fresh water flux to ocean (kg/m^2/s)
    1312             :          fsalt         ! salt flux to ocean (kg/m^2/s)
    1313             : 
    1314             :       real (kind=dbl_kind), intent(inout), optional :: &
    1315             :          frz_onset ! day of year that freezing begins (congel or frazil)
    1316             : 
    1317             :       real (kind=dbl_kind), intent(in), optional :: &
    1318             :          yday      ! day of year
    1319             : 
    1320             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    1321             :          salinz     ! initial salinity profile
    1322             : 
    1323             :       real (kind=dbl_kind), intent(in) :: &
    1324             :          phi_init     , & ! initial frazil liquid fraction
    1325             :          dSin0_frazil     ! initial frazil bulk salinity reduction from sss
    1326             : 
    1327             :       logical (kind=log_kind), intent(in) :: &
    1328             :          update_ocn_f ! if true, update fresh water and salt fluxes
    1329             : 
    1330             :       ! BGC
    1331             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
    1332             :          bgrid              ! biology nondimensional vertical grid points
    1333             : 
    1334             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
    1335             :          igrid              ! biology vertical interface points
    1336             :  
    1337             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
    1338             :          cgrid              ! CICE vertical coordinate   
    1339             : 
    1340             :       integer (kind=int_kind), intent(in) :: &
    1341             :          nbtrcr          ! number of biology tracers
    1342             : 
    1343             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1344             :          flux_bio   ! tracer flux to ocean from biology (mmol/m^2/s) 
    1345             :         
    1346             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1347             :          ocean_bio   ! ocean concentration of biological tracer
    1348             : 
    1349             :       ! zsalinity
    1350             :       real (kind=dbl_kind),  intent(inout) :: &
    1351             :          fzsal      ! salt flux to ocean from zsalinity (kg/m^2/s)
    1352             : 
    1353             :       ! water isotopes
    1354             : 
    1355             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1356             :          fiso_ocn       ! isotope flux to ocean  (kg/m^2/s)
    1357             : 
    1358             :       real (kind=dbl_kind), intent(in) :: &
    1359             :          HDO_ocn    , & ! ocean concentration of HDO (kg/kg)
    1360             :          H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
    1361             :          H2_18O_ocn     ! ocean concentration of H2_18O (kg/kg)
    1362             : 
    1363             :       ! floe size distribution
    1364             :       real (kind=dbl_kind), intent(in) :: &
    1365             :          wave_sig_ht    ! significant height of waves globally (m)
    1366             : 
    1367             :       real (kind=dbl_kind), dimension(:), intent(in)  :: &
    1368             :          wave_spectrum  ! ocean surface wave spectrum, E(f) (m^2 s)
    1369             : 
    1370             :       real(kind=dbl_kind), dimension(:), intent(in) :: &
    1371             :          wavefreq,              & ! wave frequencies (s^-1)
    1372             :          dwavefreq                ! wave frequency bin widths (s^-1)
    1373             : 
    1374             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1375             :          floe_rad_c     , & ! fsd size bin centre in m (radius)
    1376             :          floe_binwidth      ! fsd size bin width in m (radius)
    1377             : 
    1378             :       real (kind=dbl_kind), dimension(ncat) :: &  ! for now
    1379             :                             ! change in thickness distribution (area)
    1380     3372984 :          d_an_latg      , & ! due to fsd lateral growth
    1381     3372984 :          d_an_newi          ! new ice formation
    1382             : 
    1383             :       real (kind=dbl_kind), dimension(:), intent(out) :: &
    1384             :                             ! change in thickness distribution (area)
    1385             :          d_afsd_latg    , & ! due to fsd lateral growth
    1386             :          d_afsd_newi        ! new ice formation
    1387             : 
    1388             :       ! local variables
    1389             : 
    1390             :       integer (kind=int_kind) :: &
    1391             :          ncats        , & ! max categories to which new ice is added, initially
    1392             :          n            , & ! ice category index
    1393             :          k            , & ! ice layer index
    1394             :          it               ! aerosol tracer index
    1395             : 
    1396             :       real (kind=dbl_kind) :: &
    1397      355716 :          ai0new       , & ! area of new ice added to cat 1
    1398      355716 :          vi0new       , & ! volume of new ice added to cat 1
    1399             :          vi0new_lat   , & ! volume of new ice added laterally to fsd
    1400      355716 :          hsurp        , & ! thickness of new ice added to each cat
    1401      355716 :          fnew         , & ! heat flx to open water for new ice (W/m^2)
    1402      355716 :          hi0new       , & ! thickness of new ice
    1403      355716 :          hi0max       , & ! max allowed thickness of new ice
    1404      355716 :          vsurp        , & ! volume of new ice added to each cat
    1405      355716 :          vtmp         , & ! total volume of new and old ice
    1406      355716 :          area1        , & ! starting fractional area of existing ice
    1407      355716 :          alvl         , & ! starting level ice area
    1408      355716 :          dfresh       , & ! change in fresh
    1409      355716 :          dfsalt       , & ! change in fsalt
    1410      355716 :          vi0tmp       , & ! frzmlt part of frazil
    1411      355716 :          Ti           , & ! frazil temperature
    1412      355716 :          qi0new       , & ! frazil ice enthalpy
    1413      355716 :          Si0new       , & ! frazil ice bulk salinity
    1414      355716 :          vi0_init     , & ! volume of new ice
    1415      355716 :          vice1        , & ! starting volume of existing ice
    1416      355716 :          vice_init, vice_final, & ! ice volume summed over categories
    1417      355716 :          eice_init, eice_final    ! ice energy summed over categories
    1418             : 
    1419      355716 :       real (kind=dbl_kind) :: frazil_conc
    1420             : 
    1421             :       real (kind=dbl_kind), dimension (nilyr) :: &
    1422     3874176 :          Sprofile         ! salinity profile used for new ice additions
    1423             : 
    1424             :       character (len=char_len) :: &
    1425             :          fieldid           ! field identifier
    1426             : 
    1427             :       real (kind=dbl_kind), dimension (ncat) :: &
    1428     3372984 :          eicen, &     ! energy of melting for each ice layer (J/m^2)
    1429     3372984 :          aicen_init, &    ! fractional area of ice
    1430     3372984 :          vicen_init       ! volume per unit area of ice (m)
    1431             : 
    1432             :       ! floe size distribution
    1433             :       real (kind=dbl_kind), dimension (nfsd,ncat) :: &
    1434     7203276 :          afsdn          ! floe size distribution tracer (originally areal_mfstd_init)
    1435             : 
    1436             : !      real (kind=dbl_kind), dimension (nfsd) :: &
    1437             : !         afsd      , & ! fsd tracer for each thickness category
    1438             : 
    1439             :       real (kind=dbl_kind), dimension (ncat) :: &
    1440     3372984 :          d_an_tot, & ! change in the ITD due to lateral growth and new ice
    1441     3372984 :          area2       ! area after lateral growth and before new ice formation
    1442             : 
    1443             :       real (kind=dbl_kind), dimension (ncat) :: &
    1444     3056796 :          vin0new          ! volume of new ice added to any thickness cat
    1445             : 
    1446             :       real (kind=dbl_kind), dimension (nfsd) :: &
    1447             :          afsd_ni      ! areal mFSTD after new ice added
    1448             : 
    1449             :       real (kind=dbl_kind) :: &
    1450             :          tmp, &
    1451      355716 :          latsurf_area, & ! fractional area of ice on sides of floes
    1452      355716 :          lead_area   , & ! fractional area of ice in lead region
    1453      355716 :          G_radial    , & ! lateral melt rate (m/s)
    1454      355716 :          tot_latg    , & ! total fsd lateral growth in open water
    1455      355716 :          ai0mod          ! ai0new - tot_latg
    1456             : 
    1457             :       character(len=*),parameter :: subname='(add_new_ice)'
    1458             : 
    1459             :       !-----------------------------------------------------------------
    1460             :       ! initialize
    1461             :       !-----------------------------------------------------------------
    1462             : 
    1463     1027620 :       hsurp  = c0
    1464     1027620 :       hi0new = c0
    1465     1027620 :       ai0new = c0
    1466    14706432 :       afsdn(:,:) = c0
    1467     5876136 :       d_an_latg(:) = c0
    1468     5876136 :       d_an_tot(:) = c0
    1469     5876136 :       d_an_newi(:) = c0
    1470     5876136 :       vin0new(:) = c0
    1471             : 
    1472     1027620 :       if (tr_fsd) then
    1473      993708 :           d_afsd_latg(:) = c0    ! diagnostics
    1474      993708 :           d_afsd_newi(:) = c0
    1475             :       end if
    1476             : 
    1477     5876136 :       area2(:) = aicen(:)
    1478     1027620 :       lead_area    = c0
    1479     1027620 :       latsurf_area = c0
    1480     1027620 :       G_radial     = c0
    1481     1027620 :       tot_latg     = c0
    1482     1027620 :       if (ncat > 1) then
    1483      955224 :          hi0max = hin_max(1)*0.9_dbl_kind  ! not too close to boundary
    1484             :       else
    1485       72396 :          hi0max = bignum                   ! big number
    1486             :       endif
    1487             : 
    1488     1027620 :       if (tr_fsd) then
    1489       98676 :          call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:))
    1490       98676 :          if (icepack_warnings_aborted(subname)) return
    1491             :       endif
    1492             : 
    1493     5876136 :       do n = 1, ncat
    1494     4848516 :          aicen_init(n) = aicen(n)
    1495     4848516 :          vicen_init(n) = vicen(n)
    1496     4848516 :          eicen(n) = c0
    1497     5876136 :          if (tr_fsd) then
    1498     4968540 :             do k = 1, nfsd
    1499     4968540 :                afsdn(k,n) = trcrn(nt_fsd+k-1,n)
    1500             :             enddo
    1501             :          endif
    1502             :       enddo
    1503             : 
    1504     1027620 :       if (conserv_check) then
    1505             : 
    1506      434376 :          do n = 1, ncat
    1507     2968236 :          do k = 1, nilyr
    1508     1839600 :             eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
    1509     3815640 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
    1510             :          enddo
    1511             :          enddo
    1512       72396 :          call column_sum (ncat, vicen, vice_init)
    1513       72396 :          if (icepack_warnings_aborted(subname)) return
    1514       72396 :          call column_sum (ncat, eicen, eice_init)
    1515       72396 :          if (icepack_warnings_aborted(subname)) return
    1516             : 
    1517             :       endif ! conserv_check
    1518             : 
    1519             :       !-----------------------------------------------------------------
    1520             :       ! Compute average enthalpy of new ice.
    1521             :       ! Sprofile is the salinity profile used when adding new ice to
    1522             :       ! all categories, for ktherm/=2, and to category 1 for all ktherm.
    1523             :       !
    1524             :       ! NOTE:  POP assumes new ice is fresh!
    1525             :       !-----------------------------------------------------------------
    1526             : 
    1527     1027620 :       if (ktherm == 2) then  ! mushy
    1528      665640 :          if (sss > c2 * dSin0_frazil) then
    1529      665640 :             Si0new = sss - dSin0_frazil
    1530             :          else
    1531           0 :             Si0new = sss**2 / (c4*dSin0_frazil)
    1532             :          endif
    1533     5325120 :          do k = 1, nilyr
    1534     5325120 :             Sprofile(k) = Si0new
    1535             :          enddo
    1536      665640 :          Ti = min(liquidus_temperature_mush(Si0new/phi_init), -p1)
    1537      665640 :          qi0new = enthalpy_mush(Ti, Si0new)
    1538             :       else
    1539     2027088 :          do k = 1, nilyr
    1540     2027088 :             Sprofile(k) = salinz(k)
    1541             :          enddo
    1542      361980 :          qi0new = -rhoi*Lfresh
    1543             :       endif    ! ktherm
    1544             : 
    1545             :       !-----------------------------------------------------------------
    1546             :       ! Compute the volume, area, and thickness of new ice.
    1547             :       !-----------------------------------------------------------------
    1548             : 
    1549     1027620 :       fnew = max (frzmlt, c0)    ! fnew > 0 iff frzmlt > 0
    1550     1027620 :       vi0new = -fnew*dt / qi0new ! note sign convention, qi < 0
    1551     1027620 :       vi0_init = vi0new          ! for bgc
    1552             : 
    1553             :       ! increment ice volume and energy
    1554     1027620 :       if (conserv_check) then
    1555       72396 :          vice_init = vice_init + vi0new
    1556       72396 :          eice_init = eice_init + vi0new*qi0new
    1557             :       endif
    1558             : 
    1559             :       ! history diagnostics
    1560     1027620 :       frazil = vi0new
    1561     1027620 :       if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss
    1562             : 
    1563     1027620 :       if (present(frz_onset) .and. present(yday)) then
    1564     1027620 :          if (frazil > puny .and. frz_onset < puny) frz_onset = yday
    1565             :       endif
    1566             : 
    1567             :       !-----------------------------------------------------------------
    1568             :       ! Update fresh water and salt fluxes.
    1569             :       !
    1570             :       ! NOTE: POP assumes fresh water and salt flux due to frzmlt > 0
    1571             :       !       is NOT included in fluxes fresh and fsalt.
    1572             :       !-----------------------------------------------------------------
    1573             : 
    1574     1027620 :       if (update_ocn_f) then
    1575       72396 :          if (ktherm <= 1) then
    1576       72396 :             dfresh = -rhoi*vi0new/dt 
    1577       72396 :             dfsalt = ice_ref_salinity*p001*dfresh
    1578       72396 :             fresh  = fresh + dfresh
    1579       72396 :             fsalt  = fsalt + dfsalt
    1580             :          ! elseif (ktherm == 2) the fluxes are added elsewhere
    1581             :          endif
    1582             :       else ! update_ocn_f = false
    1583      955224 :          if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
    1584      665640 :             vi0tmp = fnew*dt / (rhoi*Lfresh)
    1585      665640 :             dfresh = -rhoi*(vi0new - vi0tmp)/dt
    1586      665640 :             dfsalt = ice_ref_salinity*p001*dfresh
    1587      665640 :             fresh  = fresh + dfresh
    1588      665640 :             fsalt  = fsalt + dfsalt
    1589      665640 :             frazil_diag = frazil - vi0tmp
    1590             :          ! elseif ktherm==1 do nothing
    1591             :          endif
    1592             :       endif
    1593             : 
    1594             :       !-----------------------------------------------------------------
    1595             :       ! Decide how to distribute the new ice.
    1596             :       !-----------------------------------------------------------------
    1597             : 
    1598     1027620 :       if (vi0new > c0) then
    1599             : 
    1600      497907 :         if (tr_fsd) & ! lateral growth of existing ice
    1601             :             ! calculate change in conc due to lateral growth
    1602             :             ! update vi0new, without change to afsdn or aicen
    1603             :             call fsd_lateral_growth (ncat,       nfsd,         &
    1604             :                                   dt,         aice,         &
    1605           0 :                                   aicen,      vicen,        &
    1606             :                                   vi0new,     frazil,       &
    1607           0 :                                   floe_rad_c, afsdn,        &
    1608             :                                   lead_area,  latsurf_area, &
    1609             :                                   G_radial,   d_an_latg,    &
    1610       61115 :                                   tot_latg)
    1611             : 
    1612      497907 :          if (icepack_warnings_aborted(subname)) return
    1613             : 
    1614      497907 :          ai0mod = aice0
    1615             :          ! separate frazil ice growth from lateral ice growth
    1616      497907 :          if (tr_fsd) ai0mod = aice0-tot_latg
    1617             : 
    1618             :          ! new ice area and thickness
    1619             :          ! hin_max(0) < new ice thickness < hin_max(1)
    1620      497907 :          if (ai0mod > puny) then
    1621      497885 :             hi0new = max(vi0new/ai0mod, hfrazilmin)
    1622      497885 :             if (hi0new > hi0max .and. ai0mod+puny < c1) then
    1623             :                ! distribute excess volume over all categories (below)
    1624        9386 :                hi0new = hi0max
    1625        9386 :                ai0new = ai0mod
    1626        9386 :                vsurp  = vi0new - ai0new*hi0new
    1627        9386 :                hsurp  = vsurp / aice
    1628        9386 :                vi0new = ai0new*hi0new
    1629             :             else
    1630             :                ! put ice in a single category, with hsurp = 0
    1631      488499 :                ai0new = vi0new/hi0new
    1632             :             endif
    1633             :          else                ! aice0 < puny
    1634          22 :             hsurp = vi0new/aice ! new thickness in each cat
    1635          22 :             vi0new = c0
    1636             :          endif               ! aice0 > puny
    1637             : 
    1638             :          ! volume added to each category from lateral growth
    1639     2700018 :          do n = 1, ncat
    1640     2700018 :             if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n)
    1641             :          end do
    1642             : 
    1643             :          ! combine areal change from new ice growth and lateral growth
    1644      497907 :          d_an_newi(1)     = ai0new
    1645     2202111 :          d_an_tot(2:ncat) = d_an_latg(2:ncat)
    1646      497907 :          d_an_tot(1)      = d_an_latg(1) + d_an_newi(1)
    1647      497907 :          if (tr_fsd) then
    1648       61115 :             vin0new(1)    = vin0new(1) + ai0new*hi0new ! not BFB
    1649             :          else
    1650      436792 :             vin0new(1)    = vi0new
    1651             :          endif
    1652             : 
    1653             :       endif                  ! vi0new > puny
    1654             : 
    1655             :       !-----------------------------------------------------------------
    1656             :       ! Distribute excess ice volume among ice categories by increasing
    1657             :       ! ice thickness, leaving ice area unchanged.
    1658             :       !
    1659             :       ! NOTE: If new ice contains globally conserved tracers
    1660             :       !       (e.g., isotopes from seawater), code must be added here.
    1661             :       !
    1662             :       ! The mushy formulation (ktherm=2) puts the new ice only at the
    1663             :       ! bottom of existing ice and adjusts the layers accordingly.
    1664             :       ! The other formulations distribute the new ice throughout the 
    1665             :       ! existing ice column.
    1666             :       !-----------------------------------------------------------------
    1667             : 
    1668     1027620 :       if (hsurp > c0) then   ! add ice to all categories
    1669             : 
    1670       56448 :          do n = 1, ncat
    1671             : 
    1672       47040 :             vsurp = hsurp * aicen(n)
    1673             : 
    1674             :             ! update ice age due to freezing (new ice age = dt)
    1675       47040 :             vtmp = vicen(n) + vsurp
    1676       47040 :             if (tr_iage .and. vtmp > puny) &
    1677           0 :                 trcrn(nt_iage,n) = &
    1678        2655 :                (trcrn(nt_iage,n)*vicen(n) + dt*vsurp) / vtmp
    1679             : 
    1680       47040 :             if (tr_lvl .and. vicen(n) > puny) then
    1681       22056 :                 trcrn(nt_vlvl,n) = &
    1682       22056 :                (trcrn(nt_vlvl,n)*vicen(n) + &
    1683       41642 :                 trcrn(nt_alvl,n)*vsurp) / vtmp
    1684             :             endif
    1685             : 
    1686       47040 :             if (tr_aero .and. vtmp > puny) then
    1687        5310 :                do it = 1, n_aero
    1688           0 :                   trcrn(nt_aero+2+4*(it-1),n) = &
    1689        2655 :                   trcrn(nt_aero+2+4*(it-1),n)*vicen(n) / vtmp
    1690           0 :                   trcrn(nt_aero+3+4*(it-1),n) = &
    1691        5310 :                   trcrn(nt_aero+3+4*(it-1),n)*vicen(n) / vtmp
    1692             :                enddo
    1693             :             endif
    1694             : 
    1695       47040 :            frazil_conc = c0
    1696       47040 :            if (tr_iso .and. vtmp > puny) then
    1697       10620 :              do it=1,n_iso
    1698        7965 :                if (it==1)   &
    1699        2655 :                   frazil_conc = isoice_alpha(c0,'HDO'   ,isotope_frac_method)*HDO_ocn
    1700        7965 :                if (it==2)   &
    1701        2655 :                   frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
    1702        7965 :                if (it==3)   &
    1703        2655 :                   frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
    1704             : 
    1705             :                ! dilution and uptake in the ice
    1706           0 :                trcrn(nt_isoice+it-1,n)  &
    1707           0 :                    = (trcrn(nt_isoice+it-1,n)*vicen(n) &
    1708             :                    + frazil_conc*rhoi*vsurp) &
    1709        7965 :                    / vtmp
    1710             : 
    1711           0 :                fiso_ocn(it) = fiso_ocn(it) &
    1712       10620 :                             - frazil_conc*rhoi*vsurp/dt
    1713             :              enddo
    1714             :             endif
    1715             : 
    1716             :             ! update category volumes
    1717       47040 :             vicen(n) = vtmp
    1718             : 
    1719       56448 :             if (ktherm == 2) then
    1720       47040 :                vsurp = hsurp * aicen(n)  ! note - save this above?
    1721       47040 :                vtmp = vicen(n) - vsurp   ! vicen is the new volume
    1722       47040 :                if (vicen(n) > c0) then
    1723             :                   call update_vertical_tracers(nilyr, &
    1724           0 :                               trcrn(nt_qice:nt_qice+nilyr-1,n), &
    1725       46952 :                               vtmp, vicen(n), qi0new)
    1726       46952 :                   if (icepack_warnings_aborted(subname)) return
    1727             :                   call update_vertical_tracers(nilyr, &
    1728           0 :                               trcrn(nt_sice:nt_sice+nilyr-1,n), &
    1729       46952 :                               vtmp, vicen(n), Si0new)
    1730       46952 :                   if (icepack_warnings_aborted(subname)) return
    1731             :                endif
    1732             :             else
    1733           0 :                do k = 1, nilyr
    1734             :                   ! factor of nilyr cancels out
    1735           0 :                   vsurp = hsurp * aicen(n)  ! note - save this above?
    1736           0 :                   vtmp = vicen(n) - vsurp      ! vicen is the new volume
    1737           0 :                   if (vicen(n) > c0) then
    1738             :                      ! enthalpy
    1739           0 :                      trcrn(nt_qice+k-1,n) = &
    1740           0 :                     (trcrn(nt_qice+k-1,n)*vtmp + qi0new*vsurp) / vicen(n)
    1741             :                      ! salinity
    1742           0 :                      if (.not. solve_zsal) & 
    1743           0 :                      trcrn(nt_sice+k-1,n) = &
    1744           0 :                     (trcrn(nt_sice+k-1,n)*vtmp + Sprofile(k)*vsurp) / vicen(n) 
    1745             :                   endif
    1746             :                enddo               ! k
    1747             :             endif                  ! ktherm
    1748             : 
    1749             :          enddo                     ! n
    1750             : 
    1751             :       endif ! hsurp > 0
    1752             : 
    1753             :       !-----------------------------------------------------------------
    1754             :       ! Combine new ice grown in open water with ice categories.
    1755             :       ! Using the floe size distribution, ice is added laterally to all
    1756             :       ! categories; otherwise it is added to category 1.
    1757             :       ! Assume that vsnon and esnon are unchanged.
    1758             :       ! The mushy formulation assumes salt from frazil is added uniformly
    1759             :       ! to category 1, while the others use a salinity profile.
    1760             :       !-----------------------------------------------------------------
    1761             : 
    1762     1027620 :       ncats = 1                  ! add new ice to category 1 by default
    1763     1027620 :       if (tr_fsd) ncats = ncat   ! add new ice laterally to all categories
    1764             :   
    1765             : 
    1766     2449944 :       do n = 1, ncats
    1767             : 
    1768     2449944 :       if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then  ! add ice to category n
    1769             : 
    1770      741578 :          area1    = aicen(n)   ! save
    1771      741578 :          vice1    = vicen(n)   ! save
    1772      741578 :          area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi
    1773      741578 :          aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth
    1774             :  
    1775      741578 :          aice0    = aice0    - d_an_tot(n)
    1776      741578 :          vicen(n) = vicen(n) + vin0new(n)
    1777             : 
    1778      741578 :          trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n)
    1779      741578 :          trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0)
    1780             : 
    1781      741578 :          if (tr_FY) then
    1782       28576 :             trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n)
    1783       28576 :             trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1)
    1784             :          endif
    1785             : 
    1786      741578 :          if (tr_fsd) & ! evolve the floe size distribution
    1787             :             ! both new frazil ice and lateral growth
    1788             :             call fsd_add_new_ice (ncat, n,    nfsd,          &
    1789             :                                   dt,         ai0new,        &
    1790           0 :                                   d_an_latg,  d_an_newi,     &
    1791           0 :                                   floe_rad_c, floe_binwidth, &
    1792           0 :                                   G_radial,   area2,         &
    1793             :                                   wave_sig_ht,               &
    1794           0 :                                   wave_spectrum,             &
    1795           0 :                                   wavefreq,                  &
    1796           0 :                                   dwavefreq,                 &
    1797           0 :                                   d_afsd_latg,               &
    1798           0 :                                   d_afsd_newi,               &
    1799           0 :                                   afsdn,      aicen_init,    &
    1800      304808 :                                   aicen,      trcrn)
    1801             : 
    1802      741578 :          if (icepack_warnings_aborted(subname)) return     
    1803             :  
    1804      741578 :          if (vicen(n) > puny) then
    1805      741578 :             if (tr_iage) &
    1806       28576 :                trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n)
    1807             : 
    1808      741578 :             if (tr_aero) then
    1809       56748 :                do it = 1, n_aero
    1810           0 :                   trcrn(nt_aero+2+4*(it-1),n) = &
    1811       28374 :                   trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n)
    1812           0 :                   trcrn(nt_aero+3+4*(it-1),n) = &
    1813       56748 :                   trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n)
    1814             :                enddo
    1815             :             endif
    1816             : 
    1817      741578 :            frazil_conc = c0
    1818      741578 :            if (tr_iso) then
    1819      114304 :               do it=1,n_iso
    1820       85728 :                if (it==1)   &
    1821       28576 :                   frazil_conc = isoice_alpha(c0,'HDO'   ,isotope_frac_method)*HDO_ocn
    1822       85728 :                if (it==2)   &
    1823       28576 :                   frazil_conc = isoice_alpha(c0,'H2_16O',isotope_frac_method)*H2_16O_ocn
    1824       85728 :                if (it==3)   &
    1825       28576 :                   frazil_conc = isoice_alpha(c0,'H2_18O',isotope_frac_method)*H2_18O_ocn
    1826             : 
    1827           0 :                 trcrn(nt_isoice+it-1,1) &
    1828           0 :                   = (trcrn(nt_isoice+it-1,1)*vice1) &
    1829       85728 :                   + frazil_conc*rhoi*vi0new/vicen(1)
    1830             : 
    1831           0 :                 fiso_ocn(it) = fiso_ocn(it) &
    1832      114304 :                              - frazil_conc*rhoi*vi0new/dt
    1833             :               enddo
    1834             :            endif
    1835             : 
    1836      741578 :             if (tr_lvl) then
    1837      684814 :                 alvl = trcrn(nt_alvl,n)
    1838      277853 :                 trcrn(nt_alvl,n) = &
    1839      684814 :                (trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n)
    1840      277853 :                 trcrn(nt_vlvl,n) = &
    1841      684814 :                (trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n)
    1842             :             endif
    1843             : 
    1844      741578 :             if (tr_pond_cesm .or. tr_pond_topo) then
    1845           0 :                trcrn(nt_apnd,n) = &
    1846       56764 :                trcrn(nt_apnd,n)*area1/aicen(n)
    1847      684814 :             elseif (tr_pond_lvl) then
    1848      541227 :                if (trcrn(nt_alvl,n) > puny) then
    1849      228816 :                   trcrn(nt_apnd,n) = &
    1850      541227 :                   trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n))
    1851             :                endif
    1852             :             endif
    1853             :          endif
    1854             : 
    1855     5323786 :          do k = 1, nilyr
    1856     5323786 :             if (vicen(n) > c0) then
    1857             :                ! factor of nilyr cancels out
    1858             :                ! enthalpy
    1859     3462526 :                trcrn(nt_qice+k-1,n) = &
    1860     6313471 :               (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n)
    1861             :                ! salinity
    1862     4582208 :                if (.NOT. solve_zsal)&
    1863     3462526 :                trcrn(nt_sice+k-1,n) = &
    1864     6313471 :               (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n)
    1865             :             endif
    1866             :          enddo
    1867             : 
    1868             :       endif ! vi0new > 0
    1869             : 
    1870             :       enddo ! ncats
    1871             : 
    1872     1027620 :       if (conserv_check) then
    1873             : 
    1874      434376 :          do n = 1, ncat
    1875      361980 :             eicen(n) = c0
    1876     2968236 :             do k = 1, nilyr
    1877     1839600 :                eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
    1878     3815640 :                         * vicen(n)/real(nilyr,kind=dbl_kind)
    1879             :             enddo
    1880             :          enddo
    1881       72396 :          call column_sum (ncat, vicen, vice_final)
    1882       72396 :          if (icepack_warnings_aborted(subname)) return
    1883       72396 :          call column_sum (ncat, eicen, eice_final)
    1884       72396 :          if (icepack_warnings_aborted(subname)) return
    1885             : 
    1886       72396 :          fieldid = subname//':vice'
    1887             :          call column_conservation_check (fieldid,               &
    1888             :                                          vice_init, vice_final, &
    1889       72396 :                                          puny)
    1890       72396 :          if (icepack_warnings_aborted(subname)) return
    1891       72396 :          fieldid = subname//':eice'
    1892             :          call column_conservation_check (fieldid,               &
    1893             :                                          eice_init, eice_final, &
    1894       72396 :                                          puny*Lfresh*rhoi)
    1895       72396 :          if (icepack_warnings_aborted(subname)) return
    1896             : 
    1897             :       endif ! conserv_check
    1898             : 
    1899             :       !-----------------------------------------------------------------
    1900             :       ! Biogeochemistry
    1901             :       !-----------------------------------------------------------------     
    1902     1027620 :       if (tr_brine .or. nbtrcr > 0) &
    1903             :          call add_new_ice_bgc(dt,         nblyr,                &
    1904             :                               ncat, nilyr, nltrcr, &
    1905             :                               bgrid,      cgrid,      igrid,    &
    1906           0 :                               aicen_init, vicen_init, vi0_init, &
    1907           0 :                               aicen,      vicen,      vsnon1,   &
    1908           0 :                               vi0new,     ntrcr,      trcrn,    &
    1909           0 :                               nbtrcr,     sss,        ocean_bio,&
    1910       66636 :                               flux_bio,   hsurp)
    1911     1027620 :          if (icepack_warnings_aborted(subname)) return
    1912             : 
    1913             :       end subroutine add_new_ice
    1914             : 
    1915             : !=======================================================================
    1916             : !autodocument_start icepack_step_therm2
    1917             : ! Driver for thermodynamic changes not needed for coupling:
    1918             : ! transport in thickness space, lateral growth and melting.
    1919             : !
    1920             : ! authors: William H. Lipscomb, LANL
    1921             : !          Elizabeth C. Hunke, LANL
    1922             : 
    1923     1027620 :       subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
    1924             :                                      nilyr,        nslyr,         &
    1925     1027620 :                                      hin_max,      nblyr,         &
    1926     1027620 :                                      aicen,                       &
    1927     1027620 :                                      vicen,        vsnon,         &
    1928     1027620 :                                      aicen_init,   vicen_init,    &
    1929     1027620 :                                      trcrn,                       &
    1930             :                                      aice0,        aice,          &
    1931     1027620 :                                      trcr_depend,                 &
    1932     1027620 :                                      trcr_base,    n_trcr_strata, &
    1933     1027620 :                                      nt_strata,                   &
    1934             :                                      Tf,           sss,           &
    1935     1027620 :                                      salinz,                      &
    1936             :                                      rside,        meltl,         &
    1937             :                                      fside,                       &
    1938             :                                      frzmlt,       frazil,        &
    1939             :                                      frain,        fpond,         &
    1940             :                                      fresh,        fsalt,         &
    1941             :                                      fhocn,        update_ocn_f,  &
    1942     1027620 :                                      bgrid,        cgrid,         &
    1943     1027620 :                                      igrid,        faero_ocn,     &
    1944     1027620 :                                      first_ice,    fzsal,         &
    1945     1027620 :                                      flux_bio,     ocean_bio,     &
    1946             :                                      frazil_diag,                 &
    1947             :                                      frz_onset,    yday,          &
    1948     1027620 :                                      fiso_ocn,     HDO_ocn,       &
    1949             :                                      H2_16O_ocn,   H2_18O_ocn,    &
    1950             :                                      nfsd,         wave_sig_ht,   &
    1951     1027620 :                                      wave_spectrum,               &
    1952     1027620 :                                      wavefreq,                    &
    1953     1027620 :                                      dwavefreq,                   &
    1954     1027620 :                                      d_afsd_latg,  d_afsd_newi,   &
    1955     1027620 :                                      d_afsd_latm,  d_afsd_weld,   &
    1956     1027620 :                                      floe_rad_c,   floe_binwidth)
    1957             : 
    1958             :       integer (kind=int_kind), intent(in) :: &
    1959             :          ncat     , & ! number of thickness categories
    1960             :          nfsd     , & ! number of floe size categories
    1961             :          nltrcr   , & ! number of zbgc tracers
    1962             :          nblyr    , & ! number of bio layers
    1963             :          nilyr    , & ! number of ice layers
    1964             :          nslyr        ! number of snow layers
    1965             : 
    1966             :       logical (kind=log_kind), intent(in) :: &
    1967             :          update_ocn_f     ! if true, update fresh water and salt fluxes
    1968             : 
    1969             :       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
    1970             :          hin_max      ! category boundaries (m)
    1971             : 
    1972             :       real (kind=dbl_kind), intent(in) :: &
    1973             :          dt       , & ! time step
    1974             :          Tf       , & ! freezing temperature (C)
    1975             :          sss      , & ! sea surface salinity (ppt)
    1976             :          rside    , & ! fraction of ice that melts laterally
    1977             :          fside    , & ! lateral heat flux (W/m^2)
    1978             :          frzmlt   , & ! freezing/melting potential (W/m^2)
    1979             :          wave_sig_ht ! significant height of waves in ice (m)
    1980             : 
    1981             :       real (kind=dbl_kind), dimension(:), intent(in)  :: &
    1982             :          wave_spectrum  ! ocean surface wave spectrum E(f) (m^2 s)
    1983             : 
    1984             :       real(kind=dbl_kind), dimension(:), intent(in) :: &
    1985             :          wavefreq,              & ! wave frequencies (s^-1)
    1986             :          dwavefreq                ! wave frequency bin widths (s^-1)
    1987             : 
    1988             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1989             :          floe_rad_c     , & ! fsd size bin centre in m (radius)
    1990             :          floe_binwidth      ! fsd size bin width in m (radius)
    1991             : 
    1992             :       integer (kind=int_kind), dimension (:), intent(in) :: &
    1993             :          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
    1994             :          n_trcr_strata  ! number of underlying tracer layers
    1995             : 
    1996             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
    1997             :          trcr_base      ! = 0 or 1 depending on tracer dependency
    1998             :                         ! argument 2:  (1) aice, (2) vice, (3) vsno
    1999             : 
    2000             :       integer (kind=int_kind), dimension (:,:), intent(in) :: &
    2001             :          nt_strata      ! indices of underlying tracer layers
    2002             : 
    2003             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
    2004             :          bgrid              ! biology nondimensional vertical grid points
    2005             : 
    2006             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
    2007             :          igrid              ! biology vertical interface points
    2008             :  
    2009             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
    2010             :          cgrid              ! CICE vertical coordinate   
    2011             : 
    2012             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    2013             :          salinz   , & ! initial salinity profile
    2014             :          ocean_bio    ! ocean concentration of biological tracer
    2015             : 
    2016             :       real (kind=dbl_kind), intent(inout) :: &
    2017             :          aice     , & ! sea ice concentration
    2018             :          aice0    , & ! concentration of open water
    2019             :          frain    , & ! rainfall rate (kg/m^2 s)
    2020             :          fpond    , & ! fresh water flux to ponds (kg/m^2/s)
    2021             :          fresh    , & ! fresh water flux to ocean (kg/m^2/s)
    2022             :          fsalt    , & ! salt flux to ocean (kg/m^2/s)
    2023             :          fhocn    , & ! net heat flux to ocean (W/m^2)
    2024             :          fzsal    , & ! salt flux to ocean from zsalinity (kg/m^2/s)
    2025             :          meltl    , & ! lateral ice melt         (m/step-->cm/day)
    2026             :          frazil   , & ! frazil ice growth        (m/step-->cm/day)
    2027             :          frazil_diag  ! frazil ice growth diagnostic (m/step-->cm/day)
    2028             : 
    2029             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    2030             :          aicen_init,& ! initial concentration of ice
    2031             :          vicen_init,& ! initial volume per unit area of ice          (m)
    2032             :          aicen    , & ! concentration of ice
    2033             :          vicen    , & ! volume per unit area of ice          (m)
    2034             :          vsnon    , & ! volume per unit area of snow         (m)
    2035             :          faero_ocn, & ! aerosol flux to ocean  (kg/m^2/s)
    2036             :          flux_bio     ! all bio fluxes to ocean
    2037             : 
    2038             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    2039             :          trcrn        ! tracers
    2040             :  
    2041             :       logical (kind=log_kind), dimension(:), intent(inout) :: &
    2042             :          first_ice      ! true until ice forms
    2043             : 
    2044             :       real (kind=dbl_kind), dimension(:), intent(out) :: &
    2045             :                             ! change in floe size distribution (area)
    2046             :          d_afsd_latg    , & ! due to fsd lateral growth
    2047             :          d_afsd_newi    , & ! new ice formation
    2048             :          d_afsd_latm    , & ! lateral melt
    2049             :          d_afsd_weld        ! welding
    2050             : 
    2051             :       real (kind=dbl_kind), intent(inout), optional :: &
    2052             :          frz_onset    ! day of year that freezing begins (congel or frazil)
    2053             : 
    2054             :       real (kind=dbl_kind), intent(in), optional :: &
    2055             :          yday         ! day of year
    2056             : 
    2057             :       ! water isotopes
    2058             :       real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
    2059             :          fiso_ocn       ! isotope flux to ocean  (kg/m^2/s)
    2060             : 
    2061             :       real (kind=dbl_kind), optional, intent(in) :: &
    2062             :          HDO_ocn    , & ! ocean concentration of HDO (kg/kg)
    2063             :          H2_16O_ocn , & ! ocean concentration of H2_16O (kg/kg)
    2064             :          H2_18O_ocn     ! ocean concentration of H2_18O (kg/kg)
    2065             : !autodocument_end
    2066             : 
    2067             :       ! local variables
    2068             : 
    2069             :       ! water isotopes
    2070             :       real (kind=dbl_kind), dimension(:), allocatable :: &
    2071     1027620 :          l_fiso_ocn       ! local isotope flux to ocean  (kg/m^2/s)
    2072             : 
    2073             :       real (kind=dbl_kind) :: &
    2074      355716 :          l_HDO_ocn    , & ! local ocean concentration of HDO (kg/kg)
    2075      355716 :          l_H2_16O_ocn , & ! local ocean concentration of H2_16O (kg/kg)
    2076      355716 :          l_H2_18O_ocn     ! local ocean concentration of H2_18O (kg/kg)
    2077             : 
    2078             :       character(len=*),parameter :: subname='(icepack_step_therm2)'
    2079             : 
    2080             :       !-----------------------------------------------------------------
    2081             :       ! Check optional arguments and set local values
    2082             :       !-----------------------------------------------------------------
    2083             : 
    2084     1027620 :       if (present(fiso_ocn)) then
    2085     1027620 :          allocate(l_fiso_ocn(size(fiso_ocn)))
    2086     4110480 :          l_fiso_ocn(:) = fiso_ocn(:)
    2087             :       else
    2088           0 :          allocate(l_fiso_ocn(1))
    2089           0 :          l_fiso_ocn(:) = c0
    2090             :       endif
    2091     1027620 :       l_HDO_ocn = c0
    2092     1027620 :       l_H2_16O_ocn = c0
    2093     1027620 :       l_H2_18O_ocn = c0
    2094     1027620 :       if (present(HDO_ocn)   ) l_HDO_ocn    = HDO_ocn
    2095     1027620 :       if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
    2096     1027620 :       if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
    2097             : 
    2098             :       !-----------------------------------------------------------------
    2099             :       ! Let rain drain through to the ocean.
    2100             :       !-----------------------------------------------------------------
    2101             : 
    2102     1027620 :       fresh  = fresh + frain * aice
    2103             : 
    2104             :       !-----------------------------------------------------------------
    2105             :       ! Given thermodynamic growth rates, transport ice between
    2106             :       ! thickness categories.
    2107             :       !-----------------------------------------------------------------
    2108             : 
    2109             : !      call ice_timer_start(timer_catconv)    ! category conversions
    2110             : 
    2111             :       !-----------------------------------------------------------------
    2112             :       ! Compute fractional ice area in each grid cell.
    2113             :       !-----------------------------------------------------------------
    2114             : 
    2115     1027620 :       call aggregate_area (ncat, aicen, aice, aice0)
    2116     1027620 :       if (icepack_warnings_aborted(subname)) return
    2117             : 
    2118     1027620 :       if (kitd == 1) then
    2119             : 
    2120             :       !-----------------------------------------------------------------
    2121             :       ! Identify grid cells with ice.
    2122             :       !-----------------------------------------------------------------
    2123             : 
    2124      882828 :          if (aice > puny) then
    2125             : 
    2126             :             call linear_itd (ncat,     hin_max,        &
    2127             :                              nilyr,    nslyr,          &
    2128           0 :                              ntrcr,    trcr_depend,    &
    2129           0 :                              trcr_base,        & 
    2130           0 :                              n_trcr_strata,   &
    2131           0 :                              nt_strata,                &
    2132           0 :                              aicen_init,            &
    2133           0 :                              vicen_init,            &
    2134           0 :                              aicen,                 &
    2135           0 :                              trcrn,           & 
    2136           0 :                              vicen,                 &
    2137           0 :                              vsnon,                 &
    2138             :                              aice      ,         &
    2139             :                              aice0     ,         &
    2140      852950 :                              fpond       )
    2141      852950 :             if (icepack_warnings_aborted(subname)) return
    2142             : 
    2143             :          endif ! aice > puny
    2144             : 
    2145             :       endif  ! kitd = 1
    2146             : 
    2147             : !      call ice_timer_stop(timer_catconv)    ! category conversions
    2148             : 
    2149             :       !-----------------------------------------------------------------
    2150             :       ! Add frazil ice growing in leads.
    2151             :       !-----------------------------------------------------------------
    2152             : 
    2153             :       ! identify ice-ocean cells
    2154             : 
    2155             :          call add_new_ice (ncat,          nilyr,        &
    2156             :                            nfsd,          nblyr,        &
    2157             :                            n_aero,        dt,           &
    2158             :                            ntrcr,         nltrcr,       &
    2159             :                            hin_max,       ktherm,       &
    2160           0 :                            aicen,         trcrn,        &
    2161           0 :                            vicen,         vsnon(1),     &
    2162             :                            aice0,         aice,         &
    2163             :                            frzmlt,        frazil,       &
    2164             :                            frz_onset,     yday,         &
    2165             :                            update_ocn_f,                &
    2166             :                            fresh,         fsalt,        &
    2167             :                            Tf,            sss,          &
    2168           0 :                            salinz,        phi_init,     &
    2169             :                            dSin0_frazil,  bgrid,        &
    2170             :                            cgrid,         igrid,        &
    2171           0 :                            nbtrcr,        flux_bio,     &
    2172           0 :                            ocean_bio,     fzsal,        &
    2173             :                            frazil_diag,   l_fiso_ocn,   &
    2174             :                            l_HDO_ocn,     l_H2_16O_ocn, &
    2175             :                            l_H2_18O_ocn,                &
    2176             :                            wave_sig_ht,                 &
    2177           0 :                            wave_spectrum,               &
    2178           0 :                            wavefreq,      dwavefreq,    &
    2179           0 :                            d_afsd_latg,   d_afsd_newi,  &
    2180     1027620 :                            floe_rad_c, floe_binwidth)
    2181             : 
    2182     1027620 :          if (icepack_warnings_aborted(subname)) return
    2183             : 
    2184             :       !-----------------------------------------------------------------
    2185             :       ! Melt ice laterally.
    2186             :       !-----------------------------------------------------------------
    2187             : 
    2188             :       call lateral_melt (dt,        ncat,          &
    2189             :                          nilyr,     nslyr,         &
    2190             :                          n_aero,    fpond,         &
    2191             :                          fresh,     fsalt,         &
    2192           0 :                          fhocn,     faero_ocn,     &
    2193             :                          l_fiso_ocn,               &
    2194             :                          rside,     meltl,         &
    2195             :                          fside,     sss,           &
    2196           0 :                          aicen,     vicen,         &
    2197           0 :                          vsnon,     trcrn,         &
    2198           0 :                          fzsal,     flux_bio,      &
    2199             :                          nbtrcr,    nblyr,         &
    2200           0 :                          nfsd,      d_afsd_latm,   &
    2201     1027620 :                          floe_rad_c,floe_binwidth)
    2202     1027620 :       if (icepack_warnings_aborted(subname)) return
    2203             : 
    2204             :       ! Floe welding during freezing conditions
    2205     1027620 :       if (tr_fsd) &
    2206             :       call fsd_weld_thermo (ncat,  nfsd,   &
    2207             :                             dt,    frzmlt, &
    2208           0 :                             aicen, trcrn,  &
    2209       98676 :                             d_afsd_weld)
    2210     1027620 :       if (icepack_warnings_aborted(subname)) return
    2211             : 
    2212             :       !-----------------------------------------------------------------
    2213             :       ! For the special case of a single category, adjust the area and
    2214             :       ! volume (assuming that half the volume change decreases the
    2215             :       ! thickness, and the other half decreases the area).  
    2216             :       !-----------------------------------------------------------------
    2217             : 
    2218             : !echmod: test this
    2219     1027620 :       if (ncat==1) &
    2220       26280 :          call reduce_area (hin_max   (0),                &
    2221       26280 :                            aicen     (1), vicen     (1), &
    2222       72396 :                            aicen_init(1), vicen_init(1))
    2223     1027620 :          if (icepack_warnings_aborted(subname)) return
    2224             : 
    2225             :       !-----------------------------------------------------------------
    2226             :       ! ITD cleanup: Rebin thickness categories if necessary, and remove
    2227             :       !  categories with very small areas.
    2228             :       !-----------------------------------------------------------------
    2229             : 
    2230             :       call cleanup_itd (dt,                   ntrcr,            &
    2231             :                         nilyr,                nslyr,            &
    2232             :                         ncat,                 hin_max,          &
    2233           0 :                         aicen,                trcrn(1:ntrcr,:), &
    2234           0 :                         vicen,                vsnon,            &
    2235             :                         aice0,                aice,             & 
    2236             :                         n_aero,                                 &
    2237             :                         nbtrcr,               nblyr,            &
    2238             :                         tr_aero,                                &
    2239             :                         tr_pond_topo,         heat_capacity,    &
    2240           0 :                         first_ice,                              &
    2241           0 :                         trcr_depend,          trcr_base,        &
    2242           0 :                         n_trcr_strata,        nt_strata,        &
    2243             :                         fpond,                fresh,            &
    2244             :                         fsalt,                fhocn,            &
    2245           0 :                         faero_ocn,            l_fiso_ocn,       &
    2246     1027620 :                         fzsal,                flux_bio)   
    2247     1027620 :       if (icepack_warnings_aborted(subname)) return
    2248             : 
    2249     1027620 :       if (present(fiso_ocn)) then
    2250     4110480 :          fiso_ocn(:) = l_fiso_ocn(:)
    2251             :       endif
    2252     1027620 :       deallocate(l_fiso_ocn)
    2253             : 
    2254     1027620 :       end subroutine icepack_step_therm2
    2255             : 
    2256             : !=======================================================================
    2257             : 
    2258             :       end module icepack_therm_itd
    2259             : 
    2260             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd