LCOV - code coverage report
Current view: top level - columnphysics - icepack_mechred.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 473 620 76.29 %
Date: 2020-07-03 20:11:40 Functions: 6 7 85.71 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! Driver for ice mechanical redistribution (ridging)
       4             : !
       5             : ! See these references:
       6             : !
       7             : ! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength
       8             : !  in modeling the thickness distribution of Arctic sea ice,
       9             : !  J. Geophys. Res., 100, 18,611-18,626.
      10             : !
      11             : ! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice
      12             : !  cover, Mon. Wea. Rev., 108, 1943-1973, 1980.
      13             : !
      14             : ! Lipscomb, W. H., E. C. Hunke, W. Maslowski, and J. Jakacki, 2007: 
      15             : !  Improving ridging schemes for high-resolution sea ice models.
      16             : !  J. Geophys. Res. 112, C03S91, doi:10.1029/2005JC003355.
      17             : ! 
      18             : ! Rothrock, D. A., 1975: The energetics of the plastic deformation of
      19             : !  pack ice by ridging, J. Geophys. Res., 80, 4514-4519.
      20             : !
      21             : ! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, 
      22             : !  1975: The thickness distribution of sea ice, J. Geophys. Res., 
      23             : !  80, 4501-4513. 
      24             : !
      25             : ! authors: William H. Lipscomb, LANL
      26             : !          Elizabeth C. Hunke, LANL
      27             : !
      28             : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
      29             : ! 2004: Block structure added by William Lipscomb
      30             : ! 2006: New options for participation and redistribution (WHL)
      31             : ! 2006: Streamlined for efficiency by Elizabeth Hunke
      32             : !       Converted to free source form (F90)
      33             : 
      34             :       module icepack_mechred
      35             : 
      36             :       use icepack_kinds
      37             :       use icepack_parameters,  only: c0, c1, c2, c10, c25, Cf, Cp, Pstar, Cstar
      38             :       use icepack_parameters,  only: p05, p15, p25, p333, p5
      39             :       use icepack_parameters,  only: puny, Lfresh, rhoi, rhos
      40             : 
      41             :       use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
      42             :       use icepack_parameters, only: heat_capacity, conserv_check
      43             : 
      44             :       use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
      45             :       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
      46             :       use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice
      47             :       use icepack_tracers, only: nt_apnd, nt_hpnd
      48             :       use icepack_tracers, only: n_iso
      49             :       use icepack_tracers, only: icepack_compute_tracers
      50             :                            
      51             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      52             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      53             : 
      54             :       use icepack_itd, only: column_sum
      55             :       use icepack_itd, only: column_conservation_check
      56             :       use icepack_itd, only: cleanup_itd
      57             : 
      58             :       implicit none
      59             : 
      60             :       private
      61             :       public :: ridge_ice, &
      62             :                 asum_ridging, &
      63             :                 ridge_itd, &
      64             :                 icepack_ice_strength, &
      65             :                 icepack_step_ridge
      66             : 
      67             :       real (kind=dbl_kind), parameter :: & 
      68             :          exp_argmax = 100.0_dbl_kind, &    ! maximum argument of exponential for underflow
      69             :          Cs = p25         , & ! fraction of shear energy contrbtng to ridging 
      70             :          fsnowrdg = p5    , & ! snow fraction that survives in ridging 
      71             :          Gstar  = p15     , & ! max value of G(h) that participates 
      72             :                               ! (krdg_partic = 0) 
      73             :          astar  = p05     , & ! e-folding scale for G(h) participation 
      74             : !echmod         astar  = p1        , & ! e-folding scale for G(h) participation 
      75             :                               ! (krdg_partic = 1) 
      76             :          maxraft= c1      , & ! max value of hrmin - hi = max thickness 
      77             :                               ! of ice that rafts (m) 
      78             :          Hstar  = c25         ! determines mean thickness of ridged ice (m) 
      79             :                               ! (krdg_redist = 0) 
      80             :                               ! Flato & Hibler (1995) have Hstar = 100 
      81             : 
      82             : !=======================================================================
      83             : 
      84             :       contains
      85             : 
      86             : !=======================================================================
      87             : 
      88             : ! Compute changes in the ice thickness distribution due to divergence
      89             : ! and shear.
      90             : !
      91             : ! author: William H. Lipscomb, LANL
      92             : 
      93     1100016 :       subroutine ridge_ice (dt,          ndtd,       &
      94             :                             ncat,        n_aero,     &
      95             :                             nilyr,       nslyr,      &
      96     1100016 :                             ntrcr,       hin_max,    &
      97             :                             rdg_conv,    rdg_shear,  &
      98     1100016 :                             aicen,       trcrn,      &
      99     1100016 :                             vicen,       vsnon,      &
     100             :                             aice0,                   &
     101     1100016 :                             trcr_depend, trcr_base,  &
     102     1100016 :                             n_trcr_strata,           &
     103     1100016 :                             nt_strata,               &
     104             :                             krdg_partic, krdg_redist,&
     105             :                             mu_rdg,      tr_brine,   &
     106             :                             dardg1dt,    dardg2dt,   &
     107             :                             dvirdgdt,    opening,    &
     108             :                             fpond,                   &
     109             :                             fresh,       fhocn,      &
     110     1100016 :                             faero_ocn,   fiso_ocn,   &
     111     1100016 :                             aparticn,    krdgn,      &
     112     1100016 :                             aredistn,    vredistn,   &
     113     1100016 :                             dardg1ndt,   dardg2ndt,  &
     114     1100016 :                             dvirdgndt,               &
     115     1100016 :                             araftn,      vraftn,     &
     116             :                             closing_flag,closing )
     117             : 
     118             :       integer (kind=int_kind), intent(in) :: &
     119             :          ndtd       , & ! number of dynamics subcycles
     120             :          ncat  , & ! number of thickness categories
     121             :          nilyr , & ! number of ice layers
     122             :          nslyr , & ! number of snow layers
     123             :          n_aero, & ! number of aerosol tracers
     124             :          ntrcr     ! number of tracers in use
     125             : 
     126             :       real (kind=dbl_kind), intent(in) :: &
     127             :          mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) 
     128             :          dt             ! time step
     129             : 
     130             :       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
     131             :          hin_max   ! category limits (m)
     132             : 
     133             :       real (kind=dbl_kind), intent(in) :: &
     134             :          rdg_conv   , & ! normalized energy dissipation due to convergence (1/s)
     135             :          rdg_shear      ! normalized energy dissipation due to shear (1/s)
     136             :  
     137             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     138             :          aicen      , & ! concentration of ice
     139             :          vicen      , & ! volume per unit area of ice          (m)
     140             :          vsnon          ! volume per unit area of snow         (m)
     141             :  
     142             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: & 
     143             :          trcrn          ! ice tracers 
     144             : 
     145             :       real (kind=dbl_kind), intent(inout) :: & 
     146             :          aice0          ! concentration of open water
     147             : 
     148             :       integer (kind=int_kind), dimension (:), intent(in) :: &
     149             :          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
     150             :          n_trcr_strata  ! number of underlying tracer layers
     151             : 
     152             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
     153             :          trcr_base      ! = 0 or 1 depending on tracer dependency
     154             :                         ! argument 2:  (1) aice, (2) vice, (3) vsno
     155             : 
     156             :       integer (kind=int_kind), dimension (:,:), intent(in) :: &
     157             :          nt_strata      ! indices of underlying tracer layers
     158             : 
     159             :       integer (kind=int_kind), intent(in) :: &
     160             :          krdg_partic, & ! selects participation function
     161             :          krdg_redist    ! selects redistribution function
     162             : 
     163             :       logical (kind=log_kind), intent(in) :: &
     164             :          closing_flag, &! flag if closing is valid
     165             :          tr_brine       ! if .true., brine height differs from ice thickness
     166             : 
     167             :       ! optional history fields
     168             :       real (kind=dbl_kind), intent(inout), optional :: &
     169             :          dardg1dt   , & ! rate of fractional area loss by ridging ice (1/s)
     170             :          dardg2dt   , & ! rate of fractional area gain by new ridges (1/s)
     171             :          dvirdgdt   , & ! rate of ice volume ridged (m/s)
     172             :          opening    , & ! rate of opening due to divergence/shear (1/s)
     173             :          closing    , & ! rate of closing due to divergence/shear (1/s)
     174             :          fpond      , & ! fresh water flux to ponds (kg/m^2/s)
     175             :          fresh      , & ! fresh water flux to ocean (kg/m^2/s)
     176             :          fhocn          ! net heat flux to ocean (W/m^2)
     177             : 
     178             :       real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
     179             :          dardg1ndt  , & ! rate of fractional area loss by ridging ice (1/s)
     180             :          dardg2ndt  , & ! rate of fractional area gain by new ridges (1/s)
     181             :          dvirdgndt  , & ! rate of ice volume ridged (m/s)
     182             :          aparticn   , & ! participation function
     183             :          krdgn      , & ! mean ridge thickness/thickness of ridging ice
     184             :          araftn     , & ! rafting ice area
     185             :          vraftn     , & ! rafting ice volume 
     186             :          aredistn   , & ! redistribution function: fraction of new ridge area
     187             :          vredistn       ! redistribution function: fraction of new ridge volume
     188             : 
     189             :       real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
     190             :          faero_ocn      ! aerosol flux to ocean (kg/m^2/s)
     191             : 
     192             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     193             :          fiso_ocn       ! isotope flux to ocean (kg/m^2/s)
     194             : 
     195             :       ! local variables
     196             : 
     197             :       real (kind=dbl_kind), dimension (ncat) :: &
     198     3622896 :          eicen          ! energy of melting for each ice layer (J/m^2)
     199             :  
     200             :       real (kind=dbl_kind), dimension (ncat) :: &
     201     3622896 :          esnon, &       ! energy of melting for each snow layer (J/m^2)
     202     3622896 :          vbrin, &       ! ice volume with defined by brine height (m)
     203     3622896 :          sicen          ! Bulk salt in h ice (ppt*m)
     204             : 
     205             :       real (kind=dbl_kind) :: &
     206      381996 :          asum       , & ! sum of ice and open water area
     207      381996 :          aksum      , & ! ratio of area removed to area ridged
     208      381996 :          msnow_mlt  , & ! mass of snow added to ocean (kg m-2)
     209      381996 :          esnow_mlt  , & ! energy needed to melt snow in ocean (J m-2)
     210      381996 :          mpond      , & ! mass of pond added to ocean (kg m-2)
     211      381996 :          closing_net, & ! net rate at which area is removed    (1/s)
     212             :                         ! (ridging ice area - area of new ridges) / dt
     213      381996 :          divu_adv   , & ! divu as implied by transport scheme  (1/s)
     214      381996 :          opning     , & ! rate of opening due to divergence/shear
     215             :                         ! opning is a local variable;
     216             :                         ! opening is the history diagnostic variable
     217      381996 :          ardg1      , & ! fractional area loss by ridging ice
     218      381996 :          ardg2      , & ! fractional area gain by new ridges
     219      381996 :          virdg      , & ! ice volume ridged
     220      381996 :          aopen          ! area opening due to divergence/shear
     221             : 
     222             :       real (kind=dbl_kind), dimension (n_aero) :: &
     223     2179512 :          maero          ! aerosol mass added to ocean (kg m-2)
     224             : 
     225             :       real (kind=dbl_kind), dimension (n_iso) :: &
     226     1818036 :          miso          ! isotope mass added to ocean (kg m-2)
     227             : 
     228             :       real (kind=dbl_kind), dimension (0:ncat) :: &
     229     4768884 :          apartic        ! participation function; fraction of ridging
     230             :                         ! and closing associated w/ category n
     231             : 
     232             :       real (kind=dbl_kind), dimension (ncat) :: &
     233     3622896 :          hrmin      , & ! minimum ridge thickness
     234     3622896 :          hrmax      , & ! maximum ridge thickness (krdg_redist = 0)
     235     3622896 :          hrexp      , & ! ridge e-folding thickness (krdg_redist = 1) 
     236     3622896 :          krdg       , & ! mean ridge thickness/thickness of ridging ice
     237     3622896 :          ardg1n     , & ! area of ice ridged
     238     3622896 :          ardg2n     , & ! area of new ridges
     239     3622896 :          virdgn     , & ! ridging ice volume
     240     3286872 :          mraftn         ! rafting ice mask 
     241             : 
     242             :       real (kind=dbl_kind) :: &
     243      381996 :          vice_init, vice_final, & ! ice volume summed over categories
     244      381996 :          vsno_init, vsno_final, & ! snow volume summed over categories
     245      381996 :          eice_init, eice_final, & ! ice energy summed over layers
     246      381996 :          vbri_init, vbri_final, & ! ice volume in fbri*vicen summed over categories
     247      381996 :          sice_init ,sice_final, & ! ice bulk salinity summed over categories
     248      381996 :          esno_init, esno_final    ! snow energy summed over layers
     249             : 
     250             :       integer (kind=int_kind), parameter :: &
     251             :          nitermax = 20  ! max number of ridging iterations
     252             : 
     253             :       integer (kind=int_kind) :: &
     254             :          n          , & ! thickness category index
     255             :          niter      , & ! iteration counter
     256             :          k          , & ! vertical index
     257             :          it             ! tracer index
     258             : 
     259             :       real (kind=dbl_kind) :: &
     260      381996 :          dti            ! 1 / dt
     261             : 
     262             :       logical (kind=log_kind) :: &
     263             :          iterate_ridging ! if true, repeat the ridging
     264             : 
     265             :       character (len=char_len) :: &
     266             :          fieldid        ! field identifier
     267             : 
     268             :       character(len=*),parameter :: subname='(ridge_ice)'
     269             : 
     270             :       !-----------------------------------------------------------------
     271             :       ! Initialize
     272             :       !-----------------------------------------------------------------
     273             : 
     274     1100016 :       msnow_mlt = c0
     275     1100016 :       esnow_mlt = c0
     276     2133396 :       maero (:) = c0
     277     1238364 :       miso  (:) = c0
     278     1100016 :       mpond     = c0
     279     1100016 :       ardg1     = c0
     280     1100016 :       ardg2     = c0
     281     1100016 :       virdg     = c0
     282     6310512 :       ardg1n(:) = c0
     283     6310512 :       ardg2n(:) = c0
     284     6310512 :       virdgn(:) = c0
     285     6310512 :       mraftn(:) = c0
     286     1100016 :       aopen     = c0
     287             : 
     288             :       !-----------------------------------------------------------------
     289             :       ! Compute area of ice plus open water before ridging.
     290             :       !-----------------------------------------------------------------
     291             : 
     292     1100016 :       call asum_ridging (ncat, aicen, aice0, asum)
     293     1100016 :       if (icepack_warnings_aborted(subname)) return
     294             : 
     295             :       !-----------------------------------------------------------------
     296             :       ! Compute the area opening and closing.
     297             :       !-----------------------------------------------------------------
     298             : 
     299     1100016 :       if (closing_flag) then
     300             : 
     301      960984 :          opning = opening
     302      960984 :          closing_net = closing
     303      960984 :          divu_adv = opening - closing
     304             : 
     305             :       else
     306             : 
     307             :          call ridge_prep (dt,                      &
     308             :                           ncat,      hin_max,      &
     309             :                           rdg_conv,  rdg_shear,    &
     310             :                           asum,      closing_net,  &
     311      139032 :                           divu_adv,  opning)
     312             : 
     313             :       endif
     314             : 
     315     1100016 :       if (icepack_warnings_aborted(subname)) return
     316             : 
     317             :       !-----------------------------------------------------------------
     318             :       ! Compute initial values of conserved quantities. 
     319             :       !-----------------------------------------------------------------
     320             : 
     321     1100016 :       if (conserv_check) then
     322             : 
     323      434376 :          do n = 1, ncat
     324      361980 :          eicen(n) = c0
     325      361980 :          esnon(n) = c0
     326      361980 :          sicen(n) = c0
     327      361980 :          vbrin(n) = c0
     328             : 
     329     2895840 :          do k = 1, nilyr
     330     1839600 :             eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
     331     3453660 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     332     1839600 :             sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
     333     3815640 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     334             :          enddo
     335     2171880 :          do k = 1, nslyr
     336     1314000 :             esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
     337     2828880 :                      * vsnon(n)/real(nslyr,kind=dbl_kind)
     338             :          enddo
     339      361980 :          vbrin(n) = vicen(n)
     340      434376 :          if (tr_brine) vbrin(n) =  trcrn(nt_fbri,n) * vicen(n)
     341             :          enddo ! n
     342             : 
     343             :          call column_sum (ncat,                     &
     344       72396 :                           vicen, vice_init)
     345       72396 :          if (icepack_warnings_aborted(subname)) return
     346             :          call column_sum (ncat,                     &
     347       72396 :                           vsnon, vsno_init)
     348       72396 :          if (icepack_warnings_aborted(subname)) return
     349             :          call column_sum (ncat,                     &
     350       72396 :                           eicen, eice_init)
     351       72396 :          if (icepack_warnings_aborted(subname)) return
     352             :          call column_sum (ncat,                     &
     353       72396 :                           esnon, esno_init)
     354       72396 :          if (icepack_warnings_aborted(subname)) return
     355             :          call column_sum (ncat,                     &
     356       72396 :                           sicen, sice_init)
     357       72396 :          if (icepack_warnings_aborted(subname)) return
     358             :          call column_sum (ncat,                     &
     359       72396 :                           vbrin, vbri_init)
     360       72396 :          if (icepack_warnings_aborted(subname)) return
     361             : 
     362             :       endif            
     363             : 
     364     2061000 :       rdg_iteration: do niter = 1, nitermax
     365             : 
     366             :       !-----------------------------------------------------------------
     367             :       ! Compute the thickness distribution of ridging ice
     368             :       ! and various quantities associated with the new ridged ice.
     369             :       !-----------------------------------------------------------------
     370             : 
     371             :          call ridge_itd (ncat,        aice0,      &
     372           0 :                          aicen,       vicen,      &
     373             :                          krdg_partic, krdg_redist, &
     374             :                          mu_rdg,                   &
     375             :                          aksum,       apartic,    &
     376           0 :                          hrmin,       hrmax,      &
     377           0 :                          hrexp,       krdg,       &
     378           0 :                          aparticn,    krdgn,      &
     379     2061000 :                          mraftn)    
     380     2061000 :          if (icepack_warnings_aborted(subname)) return
     381             : 
     382             :       !-----------------------------------------------------------------
     383             :       ! Redistribute area, volume, and energy.
     384             :       !-----------------------------------------------------------------
     385             : 
     386             :          call ridge_shift (ntrcr,       dt,          &
     387             :                            ncat,        hin_max,     &
     388           0 :                            aicen,       trcrn,       &
     389           0 :                            vicen,       vsnon,       &
     390           0 :                            aice0,       trcr_depend, &
     391           0 :                            trcr_base,   n_trcr_strata,&
     392           0 :                            nt_strata,   krdg_redist, &
     393             :                            aksum,       apartic,     &
     394           0 :                            hrmin,       hrmax,       &
     395           0 :                            hrexp,       krdg,        &
     396             :                            closing_net, opning,      &
     397             :                            ardg1,       ardg2,       &
     398             :                            virdg,       aopen,       &
     399           0 :                            ardg1n,      ardg2n,      &
     400           0 :                            virdgn,                   &
     401             :                            nslyr,       n_aero,      &
     402             :                            msnow_mlt,   esnow_mlt,   &
     403           0 :                            maero,       miso,        &
     404             :                            mpond,       &
     405     2061000 :                            aredistn,    vredistn)    
     406     2061000 :          if (icepack_warnings_aborted(subname)) return
     407             : 
     408             :       !-----------------------------------------------------------------
     409             :       ! Make sure the new area = 1.  If not (because the closing
     410             :       ! and opening rates were reduced above), prepare to ridge again
     411             :       ! with new rates.
     412             :       !-----------------------------------------------------------------
     413             : 
     414     2061000 :          call asum_ridging (ncat, aicen, aice0, asum)
     415     2061000 :          if (icepack_warnings_aborted(subname)) return
     416             : 
     417     2061000 :          if (abs(asum - c1) < puny) then
     418     1100016 :             iterate_ridging = .false.
     419     1100016 :             closing_net = c0
     420     1100016 :             opning      = c0
     421             :          else
     422      960984 :             iterate_ridging = .true.
     423      960984 :             divu_adv = (c1 - asum) / dt
     424      960984 :             closing_net = max(c0, -divu_adv)
     425      960984 :             opning = max(c0, divu_adv)
     426             :          endif
     427             : 
     428             :       !-----------------------------------------------------------------
     429             :       ! If done, exit.  If not, prepare to ridge again.
     430             :       !-----------------------------------------------------------------
     431             : 
     432     2061000 :          if (iterate_ridging) then
     433      960984 :             write(warnstr,*) subname, 'Repeat ridging, niter =', niter
     434      960984 :             call icepack_warnings_add(warnstr)
     435             :          else
     436     1100016 :             exit rdg_iteration
     437             :          endif
     438             : 
     439      960984 :          if (niter == nitermax) then
     440           0 :             write(warnstr,*) ' '
     441           0 :             call icepack_warnings_add(warnstr)
     442           0 :             write(warnstr,*) subname, 'Exceeded max number of ridging iterations'
     443           0 :             call icepack_warnings_add(warnstr)
     444           0 :             write(warnstr,*) subname, 'max =',nitermax
     445           0 :             call icepack_warnings_add(warnstr)
     446           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     447           0 :             call icepack_warnings_add(subname//" ridge_ice: Exceeded max number of ridging iterations" ) 
     448           0 :             return
     449             :          endif
     450             : 
     451             :       enddo rdg_iteration                    ! niter
     452             : 
     453             :       !-----------------------------------------------------------------
     454             :       ! Compute final values of conserved quantities. 
     455             :       ! Check for conservation (allowing for snow thrown into ocean).
     456             :       !-----------------------------------------------------------------
     457             : 
     458     1100016 :       if (conserv_check) then
     459             : 
     460      434376 :          do n = 1, ncat
     461      361980 :          eicen(n) = c0
     462      361980 :          esnon(n) = c0
     463      361980 :          sicen(n) = c0
     464      361980 :          vbrin(n) = c0
     465             : 
     466     2895840 :          do k = 1, nilyr
     467     1839600 :             eicen(n) = eicen(n) + trcrn(nt_qice+k-1,n) &
     468     3453660 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     469     1839600 :             sicen(n) = sicen(n) + trcrn(nt_sice+k-1,n) &
     470     3815640 :                      * vicen(n)/real(nilyr,kind=dbl_kind)
     471             :          enddo
     472     2171880 :          do k = 1, nslyr
     473     1314000 :             esnon(n) = esnon(n) + trcrn(nt_qsno+k-1,n) &
     474     2828880 :                      * vsnon(n)/real(nslyr,kind=dbl_kind)
     475             :          enddo
     476      361980 :          vbrin(n) =  vicen(n)
     477      434376 :          if (tr_brine)  vbrin(n) =  trcrn(nt_fbri,n) * vbrin(n)
     478             :          enddo ! n
     479             : 
     480             :          call column_sum (ncat,                     &
     481       72396 :                           vicen, vice_final)
     482       72396 :          if (icepack_warnings_aborted(subname)) return
     483             :          call column_sum (ncat,                     &
     484       72396 :                           vsnon, vsno_final)
     485       72396 :          if (icepack_warnings_aborted(subname)) return
     486             :          call column_sum (ncat,                     &
     487       72396 :                           eicen, eice_final)
     488       72396 :          if (icepack_warnings_aborted(subname)) return
     489             :          call column_sum (ncat,                     &
     490       72396 :                           esnon, esno_final)
     491       72396 :          if (icepack_warnings_aborted(subname)) return
     492             :          call column_sum (ncat,                     &
     493       72396 :                           sicen, sice_final)
     494       72396 :          if (icepack_warnings_aborted(subname)) return
     495             :          call column_sum (ncat,                     &
     496       72396 :                           vbrin, vbri_final)
     497       72396 :          if (icepack_warnings_aborted(subname)) return
     498             : 
     499       72396 :          vsno_final = vsno_final + msnow_mlt/rhos
     500       72396 :          esno_final = esno_final + esnow_mlt
     501             : 
     502       72396 :          fieldid = subname//':vice'
     503             :          call column_conservation_check (fieldid,               &
     504             :                                          vice_init, vice_final, &
     505       72396 :                                          puny)
     506       72396 :          if (icepack_warnings_aborted(subname)) return
     507       72396 :          fieldid = subname//':vsno'
     508             :          call column_conservation_check (fieldid,               &
     509             :                                          vsno_init, vsno_final, &
     510       72396 :                                          puny)
     511       72396 :          if (icepack_warnings_aborted(subname)) return
     512       72396 :          fieldid = subname//':eice'
     513             :          call column_conservation_check (fieldid,               &
     514             :                                          eice_init, eice_final, &
     515       72396 :                                          puny*Lfresh*rhoi)
     516       72396 :          if (icepack_warnings_aborted(subname)) return
     517       72396 :          fieldid = subname//':esno'
     518             :          call column_conservation_check (fieldid,               &
     519             :                                          esno_init, esno_final, &
     520       72396 :                                          puny*Lfresh*rhos)
     521       72396 :          if (icepack_warnings_aborted(subname)) return
     522       72396 :          fieldid = subname//':sice'
     523             :          call column_conservation_check (fieldid,               &
     524             :                                          sice_init, sice_final, &
     525       72396 :                                          puny)
     526       72396 :          if (icepack_warnings_aborted(subname)) return
     527       72396 :          fieldid = subname//':vbrin'
     528             :          call column_conservation_check (fieldid,               &
     529             :                                          vbri_init, vbri_final, &
     530       72396 :                                          puny*c10)
     531       72396 :          if (icepack_warnings_aborted(subname)) return
     532             : 
     533             :       endif                     ! conserv_check            
     534             : 
     535             :       !-----------------------------------------------------------------
     536             :       ! Compute ridging diagnostics.
     537             :       !-----------------------------------------------------------------
     538             : 
     539     1100016 :       dti = c1/dt
     540             : 
     541     1100016 :       if (present(dardg1dt)) then
     542     1100016 :          dardg1dt = ardg1*dti
     543             :       endif
     544     1100016 :       if (present(dardg2dt)) then
     545     1100016 :          dardg2dt = ardg2*dti
     546             :       endif
     547     1100016 :       if (present(dvirdgdt)) then
     548     1100016 :          dvirdgdt = virdg*dti
     549             :       endif
     550     1100016 :       if (present(opening)) then
     551     1100016 :          opening = aopen*dti
     552             :       endif
     553     1100016 :       if (present(dardg1ndt)) then
     554     6310512 :          do n = 1, ncat
     555     6310512 :             dardg1ndt(n) = ardg1n(n)*dti
     556             :          enddo
     557             :       endif
     558     1100016 :       if (present(dardg2ndt)) then
     559     6310512 :          do n = 1, ncat
     560     6310512 :             dardg2ndt(n) = ardg2n(n)*dti
     561             :          enddo
     562             :       endif
     563     1100016 :       if (present(dvirdgndt)) then
     564     6310512 :          do n = 1, ncat
     565     6310512 :             dvirdgndt(n) = virdgn(n)*dti
     566             :          enddo
     567             :       endif
     568     1100016 :       if (present(araftn)) then
     569     6310512 :          do n = 1, ncat
     570     6310512 :             araftn(n) = mraftn(n)*ardg2n(n)
     571             : !            araftn(n) = mraftn(n)*ardg1n(n)*p5
     572             :          enddo
     573             :       endif
     574     1100016 :       if (present(vraftn)) then
     575     6310512 :          do n = 1, ncat
     576     6310512 :             vraftn(n) = mraftn(n)*virdgn(n)
     577             :          enddo
     578             :       endif
     579             : 
     580             :       !-----------------------------------------------------------------
     581             :       ! Update fresh water and heat fluxes due to snow melt.
     582             :       !-----------------------------------------------------------------
     583             : 
     584             :       ! use thermodynamic time step (ndtd*dt here) to average properly
     585     1100016 :       dti = c1/(ndtd*dt)
     586             : 
     587     1100016 :       if (present(fresh)) then
     588     1100016 :          fresh = fresh + msnow_mlt*dti
     589             :       endif
     590     1100016 :       if (present(fhocn)) then
     591     1100016 :          fhocn = fhocn + esnow_mlt*dti
     592             :       endif
     593     1100016 :       if (present(faero_ocn)) then
     594     2133396 :          do it = 1, n_aero
     595     2133396 :             faero_ocn(it) = faero_ocn(it) + maero(it)*dti
     596             :          enddo
     597             :       endif
     598     1100016 :       if (tr_iso) then
     599             :          ! check size fiso_ocn vs n_iso ???
     600      184464 :          do it = 1, n_iso
     601      184464 :             fiso_ocn(it) = fiso_ocn(it) + miso(it)*dti
     602             :          enddo
     603             :       endif
     604     1100016 :       if (present(fpond)) then
     605     1100016 :          fpond = fpond - mpond ! units change later
     606             :       endif
     607             : 
     608             :       !-----------------------------------------------------------------
     609             :       ! Check for fractional ice area > 1.
     610             :       !-----------------------------------------------------------------
     611             : 
     612     1100016 :       if (abs(asum - c1) > puny) then
     613           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     614           0 :          call icepack_warnings_add(subname//" total area > 1" ) 
     615             : 
     616           0 :          write(warnstr,*) ' '
     617           0 :          call icepack_warnings_add(warnstr)
     618           0 :          write(warnstr,*) subname, 'Ridging error: total area > 1'
     619           0 :          call icepack_warnings_add(warnstr)
     620           0 :          write(warnstr,*) subname, 'area:', asum
     621           0 :          call icepack_warnings_add(warnstr)
     622           0 :          write(warnstr,*) subname, 'n, aicen:'
     623           0 :          call icepack_warnings_add(warnstr)
     624           0 :          write(warnstr,*) subname,  0, aice0
     625           0 :          call icepack_warnings_add(warnstr)
     626           0 :          do n = 1, ncat
     627           0 :             write(warnstr,*) subname, n, aicen(n)
     628           0 :             call icepack_warnings_add(warnstr)
     629             :          enddo
     630           0 :          return
     631             :       endif
     632             : 
     633             :       end subroutine ridge_ice
     634             : 
     635             : !=======================================================================
     636             : 
     637             : ! Find the total area of ice plus open water in each grid cell.
     638             : !
     639             : ! This is similar to the aggregate_area subroutine except that the
     640             : ! total area can be greater than 1, so the open water area is
     641             : ! included in the sum instead of being computed as a residual.
     642             : !
     643             : ! author: William H. Lipscomb, LANL
     644             : 
     645     3161016 :       subroutine asum_ridging (ncat, aicen, aice0, asum)
     646             : 
     647             :       integer (kind=int_kind), intent(in) :: & 
     648             :          ncat        ! number of thickness categories
     649             : 
     650             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     651             :          aicen          ! concentration of ice in each category
     652             : 
     653             :       real (kind=dbl_kind), intent(in) :: &
     654             :          aice0          ! concentration of open water
     655             : 
     656             :       real (kind=dbl_kind), intent(out):: &
     657             :          asum           ! sum of ice and open water area
     658             : 
     659             :       ! local variables
     660             : 
     661             :       integer (kind=int_kind) :: n
     662             : 
     663             :       character(len=*),parameter :: subname='(asum_ridging)'
     664             : 
     665     3161016 :       asum = aice0
     666    18097344 :       do n = 1, ncat
     667    18097344 :             asum = asum + aicen(n)
     668             :       enddo
     669             : 
     670     3161016 :       end subroutine asum_ridging
     671             : 
     672             : !=======================================================================
     673             : 
     674             : ! Initialize arrays, compute area of closing and opening
     675             : !
     676             : ! author: William H. Lipscomb, LANL
     677             : 
     678      139032 :       subroutine ridge_prep (dt,                          &
     679      139032 :                              ncat,       hin_max,         &
     680             :                              rdg_conv,   rdg_shear,       &
     681             :                              asum,       closing_net,     &
     682             :                              divu_adv,   opning)
     683             : 
     684             :       integer (kind=int_kind), intent(in) :: & 
     685             :          ncat        ! number of thickness categories
     686             : 
     687             :       real (kind=dbl_kind), intent(in) :: &
     688             :          dt             ! time step (s)
     689             : 
     690             :       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
     691             :          hin_max   ! category limits (m)
     692             : 
     693             :       real (kind=dbl_kind), intent(in) :: &
     694             :          rdg_conv   , & ! normalized energy dissipation due to convergence (1/s)
     695             :          rdg_shear      ! normalized energy dissipation due to shear (1/s)
     696             : 
     697             :       real (kind=dbl_kind), intent(inout):: &
     698             :          asum           ! sum of ice and open water area
     699             : 
     700             :       real (kind=dbl_kind), &
     701             :          intent(out):: &
     702             :          closing_net, & ! net rate at which area is removed    (1/s)
     703             :          divu_adv   , & ! divu as implied by transport scheme  (1/s)
     704             :          opning         ! rate of opening due to divergence/shear
     705             : 
     706             :       ! local variables
     707             : 
     708             :       real (kind=dbl_kind), parameter :: &
     709             :          big = 1.0e+8_dbl_kind
     710             : 
     711             :       character(len=*),parameter :: subname='(ridge_prep)'
     712             : 
     713             :       ! Set hin_max(ncat) to a big value to ensure that all ridged ice
     714             :       ! is thinner than hin_max(ncat).
     715      139032 :       hin_max(ncat) = big
     716             : 
     717             :       !-----------------------------------------------------------------
     718             :       ! Compute the net rate of closing due to convergence
     719             :       ! and shear, based on Flato and Hibler (1995).
     720             :       !
     721             :       ! For the elliptical yield curve:
     722             :       !    rdg_conv  = -min (divu, 0)
     723             :       !    rdg_shear = (1/2) * (Delta - abs(divu))
     724             :       ! Note that the shear term also accounts for divergence.
     725             :       !
     726             :       ! The energy dissipation rate is equal to the net closing rate
     727             :       ! times the ice strength.
     728             :       !
     729             :       ! NOTE: The NET closing rate is equal to the rate that open water
     730             :       !  area is removed, plus the rate at which ice area is removed by
     731             :       !  ridging, minus the rate at which area is added in new ridges.
     732             :       !  The GROSS closing rate is equal to the first two terms (open
     733             :       !  water closing and thin ice ridging) without the third term
     734             :       !  (thick, newly ridged ice).
     735             :       !
     736             :       ! rdg_conv is calculated differently in EAP (update_ice_rdg) and 
     737             :       ! represents closing_net directly.  In that case, rdg_shear=0.
     738             :       !-----------------------------------------------------------------
     739             : 
     740      139032 :       closing_net = Cs*rdg_shear + rdg_conv
     741             : 
     742             :       !-----------------------------------------------------------------
     743             :       ! Compute divu_adv, the divergence rate given by the transport/
     744             :       ! advection scheme, which may not be equal to divu as computed
     745             :       ! from the velocity field.
     746             :       !
     747             :       ! If divu_adv < 0, make sure the closing rate is large enough
     748             :       ! to give asum = 1.0 after ridging.
     749             :       !-----------------------------------------------------------------
     750             : 
     751      139032 :       divu_adv = (c1-asum) / dt
     752             : 
     753      139032 :       if (divu_adv < c0) closing_net = max(closing_net, -divu_adv)
     754             : 
     755             :       !-----------------------------------------------------------------
     756             :       ! Compute the (non-negative) opening rate that will give
     757             :       ! asum = 1.0 after ridging.
     758             :       !-----------------------------------------------------------------
     759             : 
     760      139032 :       opning = closing_net + divu_adv
     761             : 
     762      139032 :       end subroutine ridge_prep
     763             : 
     764             : !=======================================================================
     765             : 
     766             : ! Compute the thickness distribution of the ice and open water
     767             : ! participating in ridging and of the resulting ridges.
     768             : !
     769             : ! This version includes new options for ridging participation and
     770             : !  redistribution.
     771             : ! The new participation scheme (krdg_partic = 1) improves stability
     772             : !  by increasing the time scale for large changes in ice strength.
     773             : ! The new exponential redistribution function (krdg_redist = 1) improves 
     774             : !  agreement between ITDs of modeled and observed ridges.   
     775             : !
     776             : ! author: William H. Lipscomb, LANL
     777             : !
     778             : ! 2006: Changed subroutine name to ridge_itd
     779             : !       Added new options for ridging participation and redistribution.  
     780             : 
     781     2061000 :       subroutine ridge_itd (ncat,        aice0,           &
     782     2061000 :                             aicen,       vicen,           &
     783             :                             krdg_partic, krdg_redist,     &
     784             :                             mu_rdg,                       &
     785     2061000 :                             aksum,       apartic,         &
     786     2061000 :                             hrmin,       hrmax,           &
     787     2061000 :                             hrexp,       krdg,            &
     788     2061000 :                             aparticn,    krdgn,           &
     789     2061000 :                             mraft)
     790             : 
     791             :       integer (kind=int_kind), intent(in) :: & 
     792             :          ncat        ! number of thickness categories
     793             : 
     794             :       real (kind=dbl_kind), intent(in) :: &
     795             :          mu_rdg , & ! gives e-folding scale of ridged ice (m^.5) 
     796             :          aice0       ! concentration of open water
     797             : 
     798             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     799             :          aicen   , & ! concentration of ice
     800             :          vicen       ! volume per unit area of ice (m)
     801             : 
     802             :       integer (kind=int_kind), intent(in) :: &
     803             :          krdg_partic  , & ! selects participation function
     804             :          krdg_redist      ! selects redistribution function
     805             : 
     806             :       real (kind=dbl_kind), intent(out):: &
     807             :          aksum       ! ratio of area removed to area ridged
     808             : 
     809             :       real (kind=dbl_kind), dimension (0:ncat), intent(out) :: &
     810             :          apartic     ! participation function; fraction of ridging
     811             :                      ! and closing associated w/ category n
     812             : 
     813             :       real (kind=dbl_kind), dimension (:), intent(out) :: &
     814             :          hrmin   , & ! minimum ridge thickness
     815             :          hrmax   , & ! maximum ridge thickness (krdg_redist = 0)
     816             :          hrexp   , & ! ridge e-folding thickness (krdg_redist = 1) 
     817             :          krdg        ! mean ridge thickness/thickness of ridging ice
     818             : 
     819             :       ! diagnostic, category values
     820             :       real (kind=dbl_kind), dimension(:), intent(out), optional :: &
     821             :          aparticn, & ! participation function
     822             :          krdgn       ! mean ridge thickness/thickness of ridging ice
     823             : 
     824             :       real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
     825             :          mraft       ! rafting ice mask 
     826             : 
     827             :       ! local variables
     828             : 
     829             :       integer (kind=int_kind) :: &
     830             :          n           ! thickness category index
     831             : 
     832             :       real (kind=dbl_kind), parameter :: &
     833             :          Gstari   = c1/Gstar, &
     834             :          astari   = c1/astar
     835             : 
     836             :       real (kind=dbl_kind), dimension(-1:ncat) :: &
     837     8214912 :          Gsum        ! Gsum(n) = sum of areas in categories 0 to n
     838             : 
     839             :       real (kind=dbl_kind) :: &
     840      717192 :          work        ! temporary work variable
     841             : 
     842             :       real (kind=dbl_kind) :: &
     843      717192 :          hi      , & ! ice thickness for each cat (m)
     844      717192 :          hrmean  , & ! mean ridge thickness (m)
     845      717192 :          xtmp        ! temporary variable
     846             : 
     847             :       character(len=*),parameter :: subname='(ridge_itd)'
     848             : 
     849             :       !-----------------------------------------------------------------
     850             :       ! Initialize
     851             :       !-----------------------------------------------------------------
     852             : 
     853    15908832 :       Gsum   (-1:ncat) = c0  ! initialize
     854             : 
     855     2061000 :       Gsum   (-1) = c0     ! by definition
     856     2061000 :       if (aice0 > puny) then
     857     2039124 :          Gsum(0) = aice0
     858             :       else
     859       21876 :          Gsum(0) = Gsum(-1)
     860             :       endif
     861     2061000 :       apartic(0)  = c0
     862             : 
     863    11786832 :       do n = 1, ncat
     864     9725832 :          apartic(n) = c0
     865     9725832 :          hrmin  (n) = c0
     866     9725832 :          hrmax  (n) = c0
     867     9725832 :          hrexp  (n) = c0
     868     9725832 :          krdg   (n) = c1
     869             : 
     870             :       !-----------------------------------------------------------------
     871             :       ! Compute the thickness distribution of ice participating in ridging.
     872             :       !-----------------------------------------------------------------
     873             : 
     874             :       !-----------------------------------------------------------------
     875             :       ! First compute the cumulative thickness distribution function Gsum,
     876             :       !  where Gsum(n) is the fractional area in categories 0 to n.
     877             :       ! Ignore categories with very small areas.
     878             :       !-----------------------------------------------------------------
     879             : 
     880    11786832 :          if (aicen(n) > puny) then
     881     9003348 :             Gsum(n) = Gsum(n-1) + aicen(n)
     882             :          else
     883      722484 :             Gsum(n) = Gsum(n-1)
     884             :          endif
     885             :       enddo
     886             : 
     887             :       ! normalize
     888             : 
     889     2061000 :       if (Gsum(ncat) > c0) then
     890     2061000 :          work = c1 / Gsum(ncat)
     891    13847832 :          do n = 0, ncat
     892    13847832 :             Gsum(n) = Gsum(n) * work
     893             :          enddo
     894             :       endif
     895             : 
     896             :       !-----------------------------------------------------------------
     897             :       ! Compute the participation function apartic; this is analogous to
     898             :       ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
     899             :       !
     900             :       !                area lost from category n due to ridging/closing
     901             :       !  apartic(n) = --------------------------------------------------
     902             :       !                  total area lost due to ridging/closing
     903             :       !
     904             :       !-----------------------------------------------------------------
     905             : 
     906     2061000 :       if (krdg_partic == 0) then  ! Thornike et al. 1975 formulation
     907             : 
     908             :       !-----------------------------------------------------------------
     909             :       ! Assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
     910             :       ! The expressions for apartic are found by integrating b(h)g(h) between
     911             :       ! the category boundaries.
     912             :       !-----------------------------------------------------------------
     913             : 
     914           0 :          do n = 0, ncat
     915           0 :             if (Gsum(n) < Gstar) then
     916           0 :                apartic(n) = Gstari*(Gsum(n  ) - Gsum(n-1)) * &
     917           0 :                       (c2 - Gstari*(Gsum(n-1) + Gsum(n  )))
     918           0 :             elseif (Gsum(n-1) < Gstar) then
     919           0 :                apartic(n) = Gstari*(Gstar - Gsum(n-1)) * &
     920           0 :                       (c2 - Gstari*(Gstar + Gsum(n-1)))
     921             :             endif
     922             :          enddo                  ! n
     923             : 
     924     2061000 :       elseif (krdg_partic==1) then   ! exponential dependence on G(h)
     925             : 
     926             :       !-----------------------------------------------------------------
     927             :       ! b(h) = exp(-G(h)/astar)
     928             :       ! apartic(n) = [exp(-G(n-1)/astar - exp(-G(n)/astar] / [1-exp(-1/astar)]. 
     929             :       ! The expression for apartic is found by integrating b(h)g(h)
     930             :       ! between the category boundaries.
     931             :       !-----------------------------------------------------------------
     932             : 
     933             :          ! precompute exponential terms using Gsum as work array
     934     2061000 :          xtmp = c1 / (c1 - exp(-astari))
     935     2061000 :          Gsum(-1) = exp(-Gsum(-1)*astari) * xtmp
     936    13847832 :          do n = 0, ncat
     937    11786832 :             Gsum(n) = exp(-Gsum(n)*astari) * xtmp
     938    13847832 :             apartic(n) = Gsum(n-1) - Gsum(n)
     939             :          enddo                  ! n
     940             : 
     941             :       endif                     ! krdg_partic
     942             : 
     943             :       !-----------------------------------------------------------------
     944             :       ! Compute variables related to ITD of ridged ice:
     945             :       ! 
     946             :       ! krdg   = mean ridge thickness / thickness of ridging ice
     947             :       ! hrmin  = min ridge thickness
     948             :       ! hrmax  = max ridge thickness (krdg_redist = 0)
     949             :       ! hrexp  = ridge e-folding scale (krdg_redist = 1)
     950             :       !----------------------------------------------------------------
     951             : 
     952     2061000 :       if (krdg_redist == 0) then  ! Hibler 1980 formulation
     953             : 
     954             :       !-----------------------------------------------------------------
     955             :       ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
     956             :       ! 
     957             :       ! This parameterization is a modified version of Hibler (1980). 
     958             :       ! In the original paper the min ridging thickness is hrmin = 2*hi, 
     959             :       !  and the max thickness is hrmax = 2*sqrt(hi*Hstar). 
     960             :       ! 
     961             :       ! Here the min thickness is hrmin = min(2*hi, hi+maxraft), 
     962             :       !  so thick ridging ice is not required to raft. 
     963             :       !
     964             :       !-----------------------------------------------------------------
     965             : 
     966           0 :          do n = 1, ncat
     967           0 :             if (aicen(n) > puny) then 
     968           0 :                hi = vicen(n) / aicen(n) 
     969           0 :                hrmin(n) = min(c2*hi, hi + maxraft) 
     970           0 :                hrmax(n) = c2*sqrt(Hstar*hi) 
     971           0 :                hrmax(n) = max(hrmax(n), hrmin(n)+puny) 
     972           0 :                hrmean = p5 * (hrmin(n) + hrmax(n)) 
     973           0 :                krdg(n) = hrmean / hi 
     974             : 
     975             :                ! diagnostic rafting mask not implemented
     976             :             endif 
     977             :          enddo                  ! n
     978             : 
     979             :       else               ! krdg_redist = 1; exponential redistribution
     980             :  
     981             :       !----------------------------------------------------------------- 
     982             :       ! The ridge ITD is a negative exponential: 
     983             :       ! 
     984             :       !  g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin 
     985             :       ! 
     986             :       ! where hrmin is the minimum thickness of ridging ice and 
     987             :       ! hrexp is the e-folding thickness.
     988             :       ! 
     989             :       ! Here, assume as above that hrmin = min(2*hi, hi+maxraft).
     990             :       ! That is, the minimum ridge thickness results from rafting,
     991             :       !  unless the ice is thicker than maxraft.
     992             :       !
     993             :       ! Also, assume that hrexp = mu_rdg*sqrt(hi).
     994             :       ! The parameter mu_rdg is tuned to give e-folding scales mostly
     995             :       !  in the range 2-4 m as observed by upward-looking sonar.
     996             :       !
     997             :       ! Values of mu_rdg in the right column give ice strengths
     998             :       !  roughly equal to values of Hstar in the left column
     999             :       !  (within ~10 kN/m for typical ITDs):
    1000             :       !
    1001             :       !   Hstar     mu_rdg
    1002             :       !
    1003             :       !     25        3.0
    1004             :       !     50        4.0
    1005             :       !     75        5.0
    1006             :       !    100        6.0
    1007             :       !----------------------------------------------------------------- 
    1008             : 
    1009    11786832 :          do n = 1, ncat
    1010    11786832 :             if (aicen(n) > puny) then
    1011     9003348 :                hi = vicen(n) / aicen(n)
    1012     9003348 :                hi = max(hi,puny)
    1013     9003348 :                hrmin(n) = min(c2*hi, hi + maxraft)
    1014     9003348 :                hrexp(n) = mu_rdg * sqrt(hi)
    1015     9003348 :                krdg(n) = (hrmin(n) + hrexp(n)) / hi
    1016             : 
    1017             :    !echmod:  check computational efficiency
    1018             :                ! diagnostic rafting mask
    1019     9003348 :                if (present(mraft)) then
    1020     9003348 :                   mraft(n) = max(c0, sign(c1, hi+maxraft-hrmin(n)))
    1021     9003348 :                   xtmp = mraft(n)*((c2*hi+hrexp(n))/hi - krdg(n))
    1022     9003348 :                   mraft(n) = max(c0, sign(c1, puny-abs(xtmp)))
    1023             :                endif
    1024             :             endif
    1025             :          enddo
    1026             : 
    1027             :       endif                     ! krdg_redist
    1028             : 
    1029             :       !----------------------------------------------------------------
    1030             :       ! Compute aksum = net ice area removed / total area participating.
    1031             :       ! For instance, if a unit area of ice with h = 1 participates in
    1032             :       !  ridging to form a ridge with a = 1/3 and h = 3, then
    1033             :       !  aksum = 1 - 1/3 = 2/3.
    1034             :       !---------------------------------------------------------------- 
    1035             : 
    1036     2061000 :       aksum = apartic(0) ! area participating = area removed
    1037             : 
    1038    11786832 :       do n = 1, ncat
    1039             :          ! area participating > area removed
    1040    11786832 :          aksum = aksum + apartic(n) * (c1 - c1/krdg(n)) 
    1041             :       enddo
    1042             : 
    1043             :       ! diagnostics
    1044     2061000 :       if (present(aparticn)) then
    1045    11786832 :          do n = 1, ncat
    1046    11786832 :             aparticn(n) = apartic(n)
    1047             :          enddo
    1048             :       endif
    1049     2061000 :       if (present(krdgn)) then
    1050    11786832 :          do n = 1, ncat
    1051    11786832 :             krdgn(n) = krdg(n)
    1052             :          enddo
    1053             :       endif
    1054             : 
    1055     2061000 :       end subroutine ridge_itd
    1056             : 
    1057             : !=======================================================================
    1058             : 
    1059             : ! Remove area, volume, and energy from each ridging category
    1060             : ! and add to thicker ice categories.
    1061             : !
    1062             : ! Tracers:  Ridging conserves ice volume and therefore conserves volume
    1063             : ! tracers. It does not conserve ice area, and therefore a portion of area 
    1064             : ! tracers are lost (corresponding to the net closing).  Area tracers on 
    1065             : ! ice that participates in ridging are carried onto the resulting ridged
    1066             : ! ice (except the portion that are lost due to closing).  Therefore, 
    1067             : ! tracers must be decremented if they are lost to the ocean during ridging
    1068             : ! (e.g. snow, ponds) or if they are being carried only on the level ice 
    1069             : ! area. 
    1070             : !
    1071             : ! author: William H. Lipscomb, LANL
    1072             : 
    1073     2061000 :       subroutine ridge_shift (ntrcr,       dt,              &
    1074     2061000 :                               ncat,        hin_max,         &
    1075     4122000 :                               aicen,       trcrn,           &
    1076     2061000 :                               vicen,       vsnon,           &
    1077     2061000 :                               aice0,       trcr_depend,     &   
    1078     2061000 :                               trcr_base,   n_trcr_strata,   &
    1079     2061000 :                               nt_strata,   krdg_redist,     &
    1080     2061000 :                               aksum,       apartic,         &
    1081     2061000 :                               hrmin,       hrmax,           &
    1082     4122000 :                               hrexp,       krdg,            &
    1083             :                               closing_net, opning,          &
    1084             :                               ardg1,       ardg2,           &
    1085             :                               virdg,       aopen,           &
    1086     2061000 :                               ardg1nn,     ardg2nn,         &
    1087     2061000 :                               virdgnn,                      &
    1088             :                               nslyr,       n_aero,          &
    1089             :                               msnow_mlt,   esnow_mlt,       &
    1090     2061000 :                               maero,       miso,            &
    1091             :                               mpond,           &
    1092     2061000 :                               aredistn,    vredistn)
    1093             : 
    1094             :       integer (kind=int_kind), intent(in) :: & 
    1095             :          ncat  , & ! number of thickness categories
    1096             :          nslyr , & ! number of snow layers
    1097             :          ntrcr , & ! number of tracers in use
    1098             :          n_aero, & ! number of aerosol tracers
    1099             :          krdg_redist      ! selects redistribution function
    1100             : 
    1101             :       real (kind=dbl_kind), intent(in) :: &
    1102             :          dt             ! time step (s)
    1103             : 
    1104             :       integer (kind=int_kind), dimension (:), intent(in) :: &
    1105             :          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
    1106             :          n_trcr_strata  ! number of underlying tracer layers
    1107             : 
    1108             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
    1109             :          trcr_base      ! = 0 or 1 depending on tracer dependency
    1110             :                         ! argument 2:  (1) aice, (2) vice, (3) vsno
    1111             : 
    1112             :       integer (kind=int_kind), dimension (:,:), intent(in) :: &
    1113             :          nt_strata      ! indices of underlying tracer layers
    1114             : 
    1115             :       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
    1116             :          hin_max   ! category limits (m)
    1117             : 
    1118             :       real (kind=dbl_kind), intent(inout) :: &
    1119             :          aice0          ! concentration of open water
    1120             : 
    1121             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1122             :          aicen      , & ! concentration of ice
    1123             :          vicen      , & ! volume per unit area of ice          (m)
    1124             :          vsnon          ! volume per unit area of snow         (m)
    1125             : 
    1126             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
    1127             :          trcrn          ! ice tracers
    1128             : 
    1129             :       real (kind=dbl_kind), intent(in) :: &
    1130             :          aksum          ! ratio of area removed to area ridged
    1131             : 
    1132             :       real (kind=dbl_kind), dimension (0:ncat), intent(in) :: &
    1133             :          apartic        ! participation function; fraction of ridging
    1134             :                         ! and closing associated w/ category n
    1135             : 
    1136             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1137             :          hrmin      , & ! minimum ridge thickness
    1138             :          hrmax      , & ! maximum ridge thickness (krdg_redist = 0)
    1139             :          hrexp      , & ! ridge e-folding thickness (krdg_redist = 1) 
    1140             :          krdg           ! mean ridge thickness/thickness of ridging ice
    1141             : 
    1142             :       real (kind=dbl_kind), intent(inout) :: &
    1143             :          closing_net, & ! net rate at which area is removed    (1/s)
    1144             :          opning     , & ! rate of opening due to divergence/shear (1/s)
    1145             :          ardg1      , & ! fractional area loss by ridging ice
    1146             :          ardg2      , & ! fractional area gain by new ridges
    1147             :          virdg      , & ! ice volume ridged (m)
    1148             :          aopen          ! area opened due to divergence/shear
    1149             : 
    1150             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1151             :          ardg1nn    , & ! area of ice ridged
    1152             :          ardg2nn    , & ! area of new ridges
    1153             :          virdgnn        ! ridging ice volume
    1154             : 
    1155             :       real (kind=dbl_kind), intent(inout) :: &
    1156             :          msnow_mlt  , & ! mass of snow added to ocean (kg m-2)
    1157             :          esnow_mlt  , & ! energy needed to melt snow in ocean (J m-2)
    1158             :          mpond          ! mass of pond added to ocean (kg m-2)
    1159             : 
    1160             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1161             :          maero          ! aerosol mass added to ocean (kg m-2)
    1162             : 
    1163             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1164             :          miso           ! isotope mass added to ocean (kg m-2)
    1165             : 
    1166             :       real (kind=dbl_kind), dimension (:), intent(inout), optional :: &
    1167             :          aredistn   , & ! redistribution function: fraction of new ridge area
    1168             :          vredistn       ! redistribution function: fraction of new ridge volume
    1169             : 
    1170             :       ! local variables
    1171             : 
    1172             :       integer (kind=int_kind) :: &
    1173             :          n, nr      , & ! thickness category indices
    1174             :          k          , & ! ice layer index
    1175             :          it         , & ! tracer index
    1176             :          ntr        , & ! tracer index
    1177             :          itl            ! loop index
    1178             : 
    1179             :       real (kind=dbl_kind), dimension (ncat) :: &
    1180     8214912 :          aicen_init , & ! ice area before ridging
    1181     6780528 :          vicen_init , & ! ice volume before ridging
    1182     6780528 :          vsnon_init     ! snow volume before ridging
    1183             : 
    1184             :       real (kind=dbl_kind), dimension(ntrcr,ncat) :: &
    1185    96799392 :          atrcrn         ! aicen*trcrn
    1186             : 
    1187             :       real (kind=dbl_kind), dimension(3) :: &
    1188     2868768 :          trfactor       ! base quantity on which tracers are carried
    1189             : 
    1190             :       real (kind=dbl_kind) :: &
    1191      717192 :          work       , & ! temporary variable
    1192      717192 :          expL_arg   , & ! temporary exp arg values
    1193      717192 :          expR_arg   , & ! temporary exp arg values
    1194      717192 :          closing_gross  ! rate at which area removed, not counting
    1195             :                         ! area of new ridges
    1196             : 
    1197             : ! ECH note:  the following arrays only need be defined on iridge cells
    1198             :       real (kind=dbl_kind) :: &
    1199      717192 :          afrac      , & ! fraction of category area ridged
    1200      717192 :          ardg1n     , & ! area of ice ridged
    1201      717192 :          ardg2n     , & ! area of new ridges
    1202      717192 :          virdgn     , & ! ridging ice volume
    1203      717192 :          vsrdgn     , & ! ridging snow volume 
    1204      717192 :          dhr        , & ! hrmax - hrmin
    1205      717192 :          dhr2       , & ! hrmax^2 - hrmin^2
    1206      717192 :          farea      , & ! fraction of new ridge area going to nr
    1207      717192 :          fvol           ! fraction of new ridge volume going to nr
    1208             : 
    1209             :       real (kind=dbl_kind) :: &
    1210      717192 :          esrdgn         ! ridging snow energy
    1211             : 
    1212             :       real (kind=dbl_kind) :: &
    1213      717192 :          hi1        , & ! thickness of ridging ice
    1214      717192 :          hexp       , & ! ridge e-folding thickness
    1215      717192 :          hL, hR     , & ! left and right limits of integration
    1216     1434384 :          expL, expR , & ! exponentials involving hL, hR
    1217      717192 :          tmpfac     , & ! factor by which opening/closing rates are cut
    1218      717192 :          wk1            ! work variable
    1219             : 
    1220             :       character(len=*),parameter :: subname='(ridge_shift)'
    1221             : 
    1222    11786832 :       do n = 1, ncat
    1223             : 
    1224             :       !-----------------------------------------------------------------
    1225             :       ! Save initial state variables
    1226             :       !-----------------------------------------------------------------
    1227             : 
    1228     9725832 :          aicen_init(n) = aicen(n)
    1229     9725832 :          vicen_init(n) = vicen(n)
    1230     9725832 :          vsnon_init(n) = vsnon(n)
    1231             : 
    1232             :       !-----------------------------------------------------------------
    1233             :       ! Define variables equal to aicen*trcrn, vicen*trcrn, vsnon*trcrn
    1234             :       !-----------------------------------------------------------------
    1235             : 
    1236   291123864 :          do it = 1, ntrcr
    1237    90645480 :             atrcrn(it,n) = trcrn(it,n)*(trcr_base(it,1) * aicen(n) &
    1238    90645480 :                                       + trcr_base(it,2) * vicen(n) &
    1239   279337032 :                                       + trcr_base(it,3) * vsnon(n))
    1240   289062864 :             if (n_trcr_strata(it) > 0) then    ! additional tracer layers
    1241    60453360 :                do itl = 1, n_trcr_strata(it)
    1242    37075500 :                   ntr = nt_strata(it,itl)
    1243    60453360 :                   atrcrn(it,n) = atrcrn(it,n) * trcrn(ntr,n)
    1244             :                enddo
    1245             :             endif
    1246             :          enddo
    1247             : 
    1248             :       enddo ! ncat
    1249             : 
    1250             :       !-----------------------------------------------------------------
    1251             :       ! Based on the ITD of ridging and ridged ice, convert the net
    1252             :       !  closing rate to a gross closing rate.
    1253             :       ! NOTE: 0 < aksum <= 1
    1254             :       !-----------------------------------------------------------------
    1255             : 
    1256     2061000 :       closing_gross = closing_net / aksum
    1257             : 
    1258             :       !-----------------------------------------------------------------
    1259             :       ! Reduce the closing rate if more than 100% of the open water
    1260             :       ! would be removed.  Reduce the opening rate proportionately.
    1261             :       !-----------------------------------------------------------------
    1262             : 
    1263     2061000 :       if (apartic(0) > c0) then
    1264     2039124 :          wk1 = apartic(0) * closing_gross * dt
    1265     2039124 :          if (wk1 > aice0) then
    1266           0 :             tmpfac = aice0 / wk1
    1267           0 :             closing_gross = closing_gross * tmpfac
    1268           0 :             opning = opning * tmpfac
    1269             :          endif
    1270             :       endif
    1271             : 
    1272             :       !-----------------------------------------------------------------
    1273             :       ! Reduce the closing rate if more than 100% of any ice category
    1274             :       ! would be removed.  Reduce the opening rate proportionately.
    1275             :       !-----------------------------------------------------------------
    1276    11786832 :       do n = 1, ncat
    1277    11786832 :          if (aicen(n) > puny .and. apartic(n) > c0) then
    1278     9003348 :             wk1 = apartic(n) * closing_gross * dt
    1279     9003348 :             if (wk1 > aicen(n)) then
    1280           0 :                tmpfac = aicen(n) / wk1
    1281           0 :                closing_gross = closing_gross * tmpfac
    1282           0 :                opning = opning * tmpfac
    1283             :             endif
    1284             :          endif
    1285             :       enddo                     ! n
    1286             : 
    1287             :       !-----------------------------------------------------------------
    1288             :       ! Compute change in open water area due to closing and opening.
    1289             :       !-----------------------------------------------------------------
    1290             : 
    1291     2061000 :       aice0 = aice0 - apartic(0)*closing_gross*dt + opning*dt
    1292             : 
    1293     2061000 :       if (aice0 < -puny) then
    1294           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1295           0 :          call icepack_warnings_add(subname//' Ridging error: aice0 < 0')
    1296           0 :          write(warnstr,*) subname, 'aice0:', aice0
    1297           0 :          call icepack_warnings_add(warnstr)
    1298           0 :          return
    1299             : 
    1300     2061000 :       elseif (aice0 < c0) then    ! roundoff error
    1301           0 :          aice0 = c0
    1302             :       endif
    1303             : 
    1304     2061000 :       aopen = opning*dt  ! optional diagnostic
    1305             : 
    1306             :       !-----------------------------------------------------------------
    1307             :       ! Compute the area, volume, and energy of ice ridging in each
    1308             :       !  category, along with the area of the resulting ridge.
    1309             :       !-----------------------------------------------------------------
    1310             : 
    1311    11786832 :       do n = 1, ncat
    1312             : 
    1313             :       !-----------------------------------------------------------------
    1314             :       ! Identify grid cells with nonzero ridging
    1315             :       !-----------------------------------------------------------------
    1316             : 
    1317     3375720 :          if (aicen_init(n) > puny .and. apartic(n) > c0 &
    1318    11786832 :                                   .and. closing_gross > c0) then
    1319             : 
    1320             :       !-----------------------------------------------------------------
    1321             :       ! Compute area of ridging ice (ardg1n) and of new ridge (ardg2n).
    1322             :       ! Make sure ridging fraction <=1.  (Roundoff errors can give
    1323             :       !  ardg1 slightly greater than aicen.) 
    1324             :       !-----------------------------------------------------------------
    1325             : 
    1326     6123266 :             ardg1n = apartic(n)*closing_gross*dt
    1327             : 
    1328     6123266 :             if (ardg1n > aicen_init(n) + puny) then
    1329           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1330           0 :                call icepack_warnings_add(subname//' Ridging error: ardg > aicen')
    1331           0 :                write(warnstr,*) subname, 'n, ardg, aicen:', &
    1332           0 :                     n, ardg1n, aicen_init(n)
    1333           0 :                call icepack_warnings_add(warnstr)
    1334           0 :                return
    1335             :             else
    1336     6123266 :                ardg1n = min(aicen_init(n), ardg1n)
    1337             :             endif
    1338             : 
    1339     6123266 :             ardg2n = ardg1n / krdg(n)
    1340     6123266 :             afrac = ardg1n / aicen_init(n)
    1341             : 
    1342             :       !-----------------------------------------------------------------
    1343             :       ! Subtract area, volume, and energy from ridging category n.
    1344             :       ! Note: Tracer values are unchanged.
    1345             :       !-----------------------------------------------------------------
    1346             : 
    1347     6123266 :             virdgn = vicen_init(n) * afrac
    1348     6123266 :             vsrdgn = vsnon_init(n) * afrac
    1349             : 
    1350     6123266 :             aicen(n) = aicen(n) - ardg1n
    1351     6123266 :             vicen(n) = vicen(n) - virdgn           
    1352     6123266 :             vsnon(n) = vsnon(n) - vsrdgn
    1353             : 
    1354             :       !-----------------------------------------------------------------
    1355             :       ! Increment ridging diagnostics
    1356             :       !-----------------------------------------------------------------
    1357             : 
    1358     6123266 :             ardg1 = ardg1 + ardg1n
    1359     6123266 :             ardg2 = ardg2 + ardg2n
    1360     6123266 :             virdg = virdg + virdgn
    1361             : 
    1362     6123266 :             ardg1nn(n) = ardg1n
    1363     6123266 :             ardg2nn(n) = ardg2n
    1364     6123266 :             virdgnn(n) = virdgn
    1365             : 
    1366             :       !-----------------------------------------------------------------
    1367             :       !  Place part of the snow and tracer lost by ridging into the ocean.
    1368             :       !-----------------------------------------------------------------
    1369             : 
    1370     6123266 :             msnow_mlt = msnow_mlt + rhos*vsrdgn*(c1-fsnowrdg)
    1371             : 
    1372     6123266 :             if (tr_aero) then
    1373      679986 :                do it = 1, n_aero
    1374           0 :                   maero(it) = maero(it) &
    1375             :                             + vsrdgn*(c1-fsnowrdg) &
    1376           0 :                             *(trcrn(nt_aero  +4*(it-1),n)   &
    1377      679986 :                             + trcrn(nt_aero+1+4*(it-1),n))
    1378             :                enddo
    1379             :             endif
    1380             : 
    1381     6123266 :             if (tr_iso) then
    1382     1359692 :                do it = 1, n_iso
    1383           0 :                   miso(it) = miso(it) + vsrdgn*(c1-fsnowrdg) &
    1384           0 :                            * (trcrn(nt_isosno+it-1,n) &
    1385     1359692 :                             + trcrn(nt_isoice+it-1,n))
    1386             :                enddo
    1387             :             endif
    1388             : 
    1389     6123266 :             if (tr_pond_topo) then
    1390      117193 :                mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
    1391      674956 :                                       * trcrn(nt_hpnd,n)
    1392             :             endif
    1393             : 
    1394             :       !-----------------------------------------------------------------
    1395             :       ! Compute quantities used to apportion ice among categories
    1396             :       ! in the nr loop below
    1397             :       !-----------------------------------------------------------------
    1398             : 
    1399     6123266 :             dhr  = hrmax(n) - hrmin(n)
    1400     6123266 :             dhr2 = hrmax(n) * hrmax(n) - hrmin(n) * hrmin(n)
    1401             : 
    1402             :       !-----------------------------------------------------------------
    1403             :       ! Increment energy needed to melt snow in ocean.
    1404             :       ! Note that esnow_mlt < 0; the ocean must cool to melt snow.
    1405             :       !-----------------------------------------------------------------
    1406             : 
    1407    13586436 :             do k = 1, nslyr
    1408     5110220 :                esrdgn = vsrdgn * trcrn(nt_qsno+k-1,n) &
    1409     7463170 :                                / real(nslyr,kind=dbl_kind)
    1410    13586436 :                esnow_mlt = esnow_mlt + esrdgn*(c1-fsnowrdg)
    1411             :             enddo
    1412             : 
    1413             :       !-----------------------------------------------------------------
    1414             :       ! Subtract area- and volume-weighted tracers from category n.
    1415             :       !-----------------------------------------------------------------
    1416             : 
    1417   142499980 :             do it = 1, ntrcr
    1418             : 
    1419   136376714 :                trfactor(1) = trcr_base(it,1)*ardg1n
    1420   136376714 :                trfactor(2) = trcr_base(it,2)*virdgn
    1421   136376714 :                trfactor(3) = trcr_base(it,3)*vsrdgn
    1422             : 
    1423   136376714 :                work = c0
    1424   545506856 :                do k = 1, 3
    1425   545506856 :                   work = work + trfactor(k)*trcrn(it,n)
    1426             :                enddo
    1427   136376714 :                if (n_trcr_strata(it) > 0) then    ! additional tracer layers
    1428    39156026 :                   do itl = 1, n_trcr_strata(it)
    1429    24050040 :                      ntr = nt_strata(it,itl)
    1430    39156026 :                      work = work * trcrn(ntr,n)
    1431             :                   enddo
    1432             :                endif
    1433   142499980 :                atrcrn(it,n) = atrcrn(it,n) - work
    1434             : 
    1435             :             enddo                  ! ntrcr
    1436             : 
    1437             :       !-----------------------------------------------------------------
    1438             :       ! Add area, volume, and energy of new ridge to each category nr.
    1439             :       !-----------------------------------------------------------------
    1440             : 
    1441    36314408 :             do nr = 1, ncat
    1442             : 
    1443    30191142 :                if (krdg_redist == 0) then ! Hibler 1980 formulation
    1444             : 
    1445             :       !-----------------------------------------------------------------
    1446             :       ! Compute the fraction of ridged ice area and volume going to
    1447             :       !  thickness category nr.
    1448             :       !-----------------------------------------------------------------
    1449             : 
    1450           0 :                   if (hrmin(n) >= hin_max(nr) .or. &
    1451           0 :                       hrmax(n) <= hin_max(nr-1)) then
    1452           0 :                      hL = c0
    1453           0 :                      hR = c0
    1454             :                   else
    1455           0 :                      hL = max (hrmin(n), hin_max(nr-1))
    1456           0 :                      hR = min (hrmax(n), hin_max(nr))
    1457             :                   endif
    1458             : 
    1459           0 :                   farea = (hR-hL) / dhr
    1460           0 :                   fvol  = (hR*hR - hL*hL) / dhr2
    1461             : 
    1462             :                else         ! krdg_redist = 1; 2005 exponential formulation
    1463             : 
    1464             :       !-----------------------------------------------------------------
    1465             :       ! Compute the fraction of ridged ice area and volume going to
    1466             :       !  thickness category nr.
    1467             :       !-----------------------------------------------------------------
    1468             : 
    1469    30191142 :                   if (nr < ncat) then
    1470             : 
    1471    24067876 :                      hi1  = hrmin(n)
    1472    24067876 :                      hexp = hrexp(n)
    1473             : 
    1474    24067876 :                      if (hi1 >= hin_max(nr)) then
    1475    16379698 :                         farea = c0
    1476    16379698 :                         fvol  = c0
    1477             :                      else
    1478     7688178 :                         hL = max (hi1, hin_max(nr-1))
    1479     7688178 :                         hR = hin_max(nr)
    1480     7688178 :                         expL_arg = min(((hL-hi1)/hexp),exp_argmax)
    1481     7688178 :                         expR_arg = min(((hR-hi1)/hexp),exp_argmax)
    1482     7688178 :                         expL = exp(-(expL_arg))
    1483     7688178 :                         expR = exp(-(expR_arg))
    1484     7688178 :                         farea = expL - expR
    1485             :                         fvol  = ((hL + hexp)*expL  &
    1486     7688178 :                                     - (hR + hexp)*expR) / (hi1 + hexp)
    1487             :                      endif
    1488             : 
    1489             :                   else             ! nr = ncat
    1490             : 
    1491     6123266 :                      hi1  = hrmin(n)
    1492     6123266 :                      hexp = hrexp(n)
    1493             : 
    1494     6123266 :                      hL = max (hi1, hin_max(nr-1))
    1495     6123266 :                      expL_arg = min(((hL-hi1)/hexp),exp_argmax)
    1496     6123266 :                      expL = exp(-(expL_arg))
    1497     6123266 :                      farea = expL
    1498     6123266 :                      fvol  = (hL + hexp)*expL / (hi1 + hexp)
    1499             : 
    1500             :                   endif            ! nr < ncat
    1501             : 
    1502             :                   ! diagnostics
    1503    30191142 :                   if (n ==1) then  ! only for thinnest ridging ice
    1504     5920587 :                      if (present(aredistn)) then
    1505     5920587 :                         aredistn(nr) = farea*ardg2n
    1506             :                      endif
    1507     5920587 :                      if (present(vredistn)) then
    1508     5920587 :                         vredistn(nr) = fvol*virdgn
    1509             :                      endif
    1510             :                   endif
    1511             : 
    1512             :                endif               ! krdg_redist
    1513             : 
    1514             :       !-----------------------------------------------------------------
    1515             :       ! Transfer ice area, ice volume, and snow volume to category nr.
    1516             :       !-----------------------------------------------------------------
    1517             : 
    1518    30191142 :                aicen(nr) = aicen(nr) + farea*ardg2n
    1519    30191142 :                vicen(nr) = vicen(nr) + fvol *virdgn
    1520    30191142 :                vsnon(nr) = vsnon(nr) + fvol *vsrdgn*fsnowrdg
    1521             : 
    1522             :       !-----------------------------------------------------------------
    1523             :       ! Transfer area-weighted and volume-weighted tracers to category nr.
    1524             :       ! Note: The global sum aicen*trcrn of ice area tracers 
    1525             :       !       (trcr_depend = 0) is not conserved by ridging.
    1526             :       !       However, ridging conserves the global sum of volume
    1527             :       !       tracers (trcr_depend = 1 or 2).
    1528             :       ! Tracers associated with level ice, or that are otherwise lost
    1529             :       ! from ridging ice, are not transferred.
    1530             :       ! We assume that all pond water is lost from ridging ice.
    1531             :       !-----------------------------------------------------------------
    1532             : 
    1533   715646850 :                do it = 1, ntrcr
    1534             : 
    1535   679332442 :                   if (it /= nt_alvl .and. it /= nt_vlvl) then
    1536   625749888 :                      trfactor(1) = trcr_base(it,1)*ardg2n*farea
    1537   625749888 :                      trfactor(2) = trcr_base(it,2)*virdgn*fvol
    1538   625749888 :                      trfactor(3) = trcr_base(it,3)*vsrdgn*fvol*fsnowrdg
    1539             :                   else
    1540    53582554 :                      trfactor(1) = c0
    1541    53582554 :                      trfactor(2) = c0
    1542    53582554 :                      trfactor(3) = c0
    1543             :                   endif
    1544             : 
    1545   679332442 :                   work = c0
    1546  2717329768 :                   do k = 1, 3
    1547  2717329768 :                      work = work + trfactor(k)*trcrn(it,n)
    1548             :                   enddo
    1549   679332442 :                   if (n_trcr_strata(it) > 0) then    ! additional tracer layers
    1550   195780130 :                      do itl = 1, n_trcr_strata(it)
    1551   120250200 :                         ntr = nt_strata(it,itl)
    1552   195780130 :                         if (ntr == nt_fbri) then  ! brine fraction only
    1553           0 :                            work = work * trcrn(ntr,n)
    1554             :                         else
    1555   120250200 :                            work = c0
    1556             :                         endif
    1557             :                      enddo
    1558             :                   endif
    1559   709523584 :                   atrcrn(it,nr) = atrcrn(it,nr) + work
    1560             : 
    1561             :                enddo               ! ntrcr
    1562             : 
    1563             :             enddo                  ! nr (new ridges)
    1564             : 
    1565             :          endif                     ! nonzero ridging
    1566             : 
    1567             :       enddo                        ! n (ridging categories)
    1568             : 
    1569             :       !-----------------------------------------------------------------
    1570             :       ! Compute new tracers
    1571             :       !-----------------------------------------------------------------
    1572             : 
    1573    11786832 :       do n = 1, ncat
    1574           0 :          call icepack_compute_tracers (ntrcr,       trcr_depend,   &
    1575     3375720 :                                        atrcrn(:,n), aicen(n),      &
    1576     3375720 :                                        vicen(n),    vsnon(n),      &
    1577           0 :                                        trcr_base,   n_trcr_strata, &
    1578    13101552 :                                        nt_strata,   trcrn(:,n))
    1579    11786832 :          if (icepack_warnings_aborted(subname)) return
    1580             :       enddo
    1581             : 
    1582             :       end subroutine ridge_shift
    1583             : 
    1584             : !=======================================================================
    1585             : !autodocument_start icepack_ice_strength
    1586             : ! Compute the strength of the ice pack, defined as the energy (J m-2)
    1587             : ! dissipated per unit area removed from the ice pack under compression,
    1588             : ! and assumed proportional to the change in potential energy caused
    1589             : ! by ridging.
    1590             : !
    1591             : ! See Rothrock (1975) and Hibler (1980).
    1592             : !
    1593             : ! For simpler strength parameterization, see this reference:
    1594             : ! Hibler, W. D. III, 1979: A dynamic-thermodynamic sea ice model,
    1595             : !  J. Phys. Oceanog., 9, 817-846.
    1596             : !
    1597             : ! authors: William H. Lipscomb, LANL
    1598             : !          Elizabeth C. Hunke, LANL
    1599             : 
    1600           0 :       subroutine icepack_ice_strength (ncat,               &
    1601             :                                       aice,     vice,     &
    1602           0 :                                       aice0,    aicen,    &
    1603           0 :                                       vicen,    &
    1604             :                                       strength)
    1605             : 
    1606             :       integer (kind=int_kind), intent(in) :: & 
    1607             :          ncat       ! number of thickness categories
    1608             : 
    1609             :       real (kind=dbl_kind), intent(in) :: &
    1610             :          aice   , & ! concentration of ice
    1611             :          vice   , & ! volume per unit area of ice  (m)
    1612             :          aice0      ! concentration of open water
    1613             : 
    1614             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
    1615             :          aicen  , & ! concentration of ice
    1616             :          vicen      ! volume per unit area of ice  (m)
    1617             : 
    1618             :       real (kind=dbl_kind), intent(inout) :: &
    1619             :          strength   ! ice strength (N/m)
    1620             : 
    1621             : !autodocument_end
    1622             : 
    1623             :       ! local variables
    1624             : 
    1625             :       real (kind=dbl_kind) :: &
    1626           0 :          asum   , & ! sum of ice and open water area
    1627           0 :          aksum      ! ratio of area removed to area ridged
    1628             : 
    1629             :       real (kind=dbl_kind), dimension (0:ncat) :: &
    1630           0 :          apartic    ! participation function; fraction of ridging
    1631             :                     ! and closing associated w/ category n
    1632             : 
    1633             :       real (kind=dbl_kind), dimension (ncat) :: &
    1634           0 :          hrmin  , & ! minimum ridge thickness
    1635           0 :          hrmax  , & ! maximum ridge thickness (krdg_redist = 0)
    1636           0 :          hrexp  , & ! ridge e-folding thickness (krdg_redist = 1) 
    1637           0 :          krdg       ! mean ridge thickness/thickness of ridging ice
    1638             : 
    1639             :       integer (kind=int_kind) :: &
    1640             :          n          ! thickness category index
    1641             : 
    1642             :       real (kind=dbl_kind) :: &
    1643           0 :          hi     , & ! ice thickness (m)
    1644           0 :          h2rdg  , & ! mean value of h^2 for new ridge
    1645           0 :          dh2rdg     ! change in mean value of h^2 per unit area
    1646             :                     ! consumed by ridging 
    1647             : 
    1648             :       character(len=*),parameter :: subname='(icepack_ice_strength)'
    1649             : 
    1650           0 :       if (kstrength == 1) then  ! Rothrock 1975 formulation
    1651             : 
    1652             :       !-----------------------------------------------------------------
    1653             :       ! Compute thickness distribution of ridging and ridged ice.
    1654             :       !-----------------------------------------------------------------
    1655             : 
    1656           0 :          call asum_ridging (ncat, aicen, aice0, asum)
    1657           0 :          if (icepack_warnings_aborted(subname)) return
    1658             : 
    1659             :          call ridge_itd (ncat,     aice0,      &
    1660           0 :                          aicen,    vicen,      &
    1661             :                          krdg_partic, krdg_redist, &
    1662             :                          mu_rdg,                   &
    1663             :                          aksum,    apartic,    &
    1664           0 :                          hrmin,    hrmax,      &
    1665           0 :                          hrexp,    krdg)   
    1666           0 :          if (icepack_warnings_aborted(subname)) return
    1667             : 
    1668             :       !-----------------------------------------------------------------
    1669             :       ! Compute ice strength based on change in potential energy,
    1670             :       ! as in Rothrock (1975)
    1671             :       !-----------------------------------------------------------------
    1672             : 
    1673           0 :          if (krdg_redist==0) then ! Hibler 1980 formulation
    1674             : 
    1675           0 :             do n = 1, ncat
    1676           0 :                if (aicen(n) > puny .and. apartic(n) > c0)then
    1677           0 :                   hi = vicen(n) / aicen(n)
    1678           0 :                   h2rdg = p333 * (hrmax(n)**3 - hrmin(n)**3)  &
    1679           0 :                                / (hrmax(n) - hrmin(n)) 
    1680           0 :                   dh2rdg = -hi*hi + h2rdg/krdg(n)
    1681           0 :                   strength = strength + apartic(n) * dh2rdg
    1682             :                endif         ! aicen > puny
    1683             :             enddo               ! n
    1684             : 
    1685           0 :          elseif (krdg_redist==1) then ! exponential formulation
    1686             : 
    1687           0 :             do n = 1, ncat
    1688           0 :                if (aicen(n) > puny .and. apartic(n) > c0) then
    1689           0 :                   hi = vicen(n) / aicen(n)
    1690           0 :                   h2rdg =    hrmin(n)*hrmin(n) &
    1691           0 :                         + c2*hrmin(n)*hrexp(n) &
    1692           0 :                         + c2*hrexp(n)*hrexp(n)
    1693           0 :                   dh2rdg = -hi*hi + h2rdg/krdg(n)
    1694           0 :                   strength = strength + apartic(n) * dh2rdg
    1695             :                endif
    1696             :             enddo               ! n
    1697             : 
    1698             :          endif                  ! krdg_redist
    1699             : 
    1700           0 :          strength = Cf * Cp * strength / aksum
    1701             :                        ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow)
    1702             :                        ! Cf accounts for frictional dissipation
    1703             : 
    1704             :       else                      ! kstrength /= 1:  Hibler (1979) form
    1705             : 
    1706             :       !-----------------------------------------------------------------
    1707             :       ! Compute ice strength as in Hibler (1979)
    1708             :       !-----------------------------------------------------------------
    1709             : 
    1710           0 :          strength = Pstar*vice*exp(-Cstar*(c1-aice))
    1711             : 
    1712             :       endif                     ! kstrength
    1713             : 
    1714             :       end subroutine icepack_ice_strength
    1715             : 
    1716             : !=======================================================================
    1717             : !autodocument_start icepack_step_ridge
    1718             : ! Computes sea ice mechanical deformation
    1719             : !
    1720             : ! authors: William H. Lipscomb, LANL
    1721             : !          Elizabeth C. Hunke, LANL
    1722             : 
    1723     1100016 :       subroutine icepack_step_ridge (dt,           ndtd,         &
    1724             :                                     nilyr,        nslyr,         &
    1725             :                                     nblyr,                       &
    1726     1100016 :                                     ncat,         hin_max,       &
    1727             :                                     rdg_conv,     rdg_shear,     &
    1728     1100016 :                                     aicen,                       &
    1729     1100016 :                                     trcrn,                       &
    1730     1100016 :                                     vicen,        vsnon,         &
    1731     1100016 :                                     aice0,        trcr_depend,   &
    1732     1100016 :                                     trcr_base,    n_trcr_strata, &
    1733     1100016 :                                     nt_strata,                   &
    1734             :                                     dardg1dt,     dardg2dt,      &
    1735             :                                     dvirdgdt,     opening,       &
    1736             :                                     fpond,                       &
    1737             :                                     fresh,        fhocn,         &
    1738             :                                     n_aero,                      &
    1739     2200032 :                                     faero_ocn,    fiso_ocn,      &
    1740     1100016 :                                     aparticn,     krdgn,         &
    1741     1100016 :                                     aredistn,     vredistn,      &
    1742     1100016 :                                     dardg1ndt,    dardg2ndt,     &
    1743     1100016 :                                     dvirdgndt,                   &
    1744     1100016 :                                     araftn,       vraftn,        &
    1745             :                                     aice,         fsalt,         &
    1746     1100016 :                                     first_ice,    fzsal,         &
    1747     1100016 :                                     flux_bio,     closing )
    1748             : 
    1749             :       real (kind=dbl_kind), intent(in) :: &
    1750             :          dt           ! time step
    1751             : 
    1752             :       integer (kind=int_kind), intent(in) :: &
    1753             :          ncat  , & ! number of thickness categories
    1754             :          ndtd  , & ! number of dynamics supercycles
    1755             :          nblyr , & ! number of bio layers
    1756             :          nilyr , & ! number of ice layers
    1757             :          nslyr , & ! number of snow layers
    1758             :          n_aero    ! number of aerosol tracers
    1759             : 
    1760             :       real (kind=dbl_kind), dimension(0:ncat), intent(inout) :: &
    1761             :          hin_max   ! category limits (m)
    1762             : 
    1763             :       integer (kind=int_kind), dimension (:), intent(in) :: &
    1764             :          trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon
    1765             :          n_trcr_strata  ! number of underlying tracer layers
    1766             : 
    1767             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
    1768             :          trcr_base      ! = 0 or 1 depending on tracer dependency
    1769             :                         ! argument 2:  (1) aice, (2) vice, (3) vsno
    1770             : 
    1771             :       integer (kind=int_kind), dimension (:,:), intent(in) :: &
    1772             :          nt_strata      ! indices of underlying tracer layers
    1773             : 
    1774             :       real (kind=dbl_kind), intent(inout) :: &
    1775             :          aice     , & ! sea ice concentration
    1776             :          aice0    , & ! concentration of open water
    1777             :          rdg_conv , & ! convergence term for ridging (1/s)
    1778             :          rdg_shear, & ! shear term for ridging (1/s)
    1779             :          dardg1dt , & ! rate of area loss by ridging ice (1/s)
    1780             :          dardg2dt , & ! rate of area gain by new ridges (1/s)
    1781             :          dvirdgdt , & ! rate of ice volume ridged (m/s)
    1782             :          opening  , & ! rate of opening due to divergence/shear (1/s)
    1783             :          fpond    , & ! fresh water flux to ponds (kg/m^2/s)
    1784             :          fresh    , & ! fresh water flux to ocean (kg/m^2/s)
    1785             :          fsalt    , & ! salt flux to ocean (kg/m^2/s)
    1786             :          fhocn    , & ! net heat flux to ocean (W/m^2)
    1787             :          fzsal        ! zsalinity flux to ocean(kg/m^2/s)
    1788             : 
    1789             :       real (kind=dbl_kind), intent(inout), optional :: &
    1790             :          closing      ! rate of closing due to divergence/shear (1/s)
    1791             : 
    1792             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    1793             :          aicen    , & ! concentration of ice
    1794             :          vicen    , & ! volume per unit area of ice          (m)
    1795             :          vsnon    , & ! volume per unit area of snow         (m)
    1796             :          dardg1ndt, & ! rate of area loss by ridging ice (1/s)
    1797             :          dardg2ndt, & ! rate of area gain by new ridges (1/s)
    1798             :          dvirdgndt, & ! rate of ice volume ridged (m/s)
    1799             :          aparticn , & ! participation function
    1800             :          krdgn    , & ! mean ridge thickness/thickness of ridging ice
    1801             :          araftn   , & ! rafting ice area
    1802             :          vraftn   , & ! rafting ice volume 
    1803             :          aredistn , & ! redistribution function: fraction of new ridge area
    1804             :          vredistn , & ! redistribution function: fraction of new ridge volume
    1805             :          faero_ocn, & ! aerosol flux to ocean  (kg/m^2/s)
    1806             :          flux_bio     ! all bio fluxes to ocean
    1807             : 
    1808             :       real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
    1809             :          fiso_ocn     ! isotope flux to ocean  (kg/m^2/s)
    1810             : 
    1811             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    1812             :          trcrn        ! tracers
    1813             : 
    1814             :       !logical (kind=log_kind), intent(in) :: &
    1815             :          !tr_pond_topo,& ! if .true., use explicit topography-based ponds
    1816             :          !tr_aero     ,& ! if .true., use aerosol tracers
    1817             :          !tr_brine    !,& ! if .true., brine height differs from ice thickness
    1818             :          !heat_capacity  ! if true, ice has nonzero heat capacity
    1819             : 
    1820             :       logical (kind=log_kind), dimension(:), intent(inout) :: &
    1821             :          first_ice    ! true until ice forms
    1822             : 
    1823             : !autodocument_end
    1824             : 
    1825             :       ! local variables
    1826             : 
    1827             :       real (kind=dbl_kind) :: &
    1828      381996 :          dtt          ! thermo time step
    1829             : 
    1830             :       real (kind=dbl_kind), dimension(:), allocatable :: &
    1831     1100016 :          l_fiso_ocn     ! local isotope flux to ocean  (kg/m^2/s)
    1832             : 
    1833             :       real (kind=dbl_kind) :: &
    1834      381996 :          l_closing      ! local rate of closing due to divergence/shear (1/s)
    1835             : 
    1836             :       logical (kind=log_kind) :: &
    1837             :          l_closing_flag ! flag if closing is passed
    1838             : 
    1839             :       character(len=*),parameter :: subname='(icepack_step_ridge)'
    1840             : 
    1841             :       !-----------------------------------------------------------------
    1842             :       ! Identify ice-ocean cells.
    1843             :       ! Note:  We can not limit the loop here using aice>puny because
    1844             :       !        aice has not yet been updated since the transport (and
    1845             :       !        it may be out of whack, which the ridging helps fix).-ECH
    1846             :       !-----------------------------------------------------------------
    1847             : 
    1848     1100016 :       if (present(fiso_ocn)) then
    1849     1100016 :          allocate(l_fiso_ocn(size(fiso_ocn)))
    1850     4400064 :          l_fiso_ocn = fiso_ocn
    1851             :       else
    1852             :          ! check tr_iso = true ???
    1853           0 :          allocate(l_fiso_ocn(1))
    1854           0 :          l_fiso_ocn = c0
    1855             :       endif
    1856             : 
    1857     1100016 :       if (present(closing)) then
    1858      960984 :          l_closing_flag = .true.
    1859      960984 :          l_closing = closing
    1860             :       else
    1861      139032 :          l_closing_flag = .false.
    1862      139032 :          l_closing = c0
    1863             :       endif
    1864             : 
    1865             :       call ridge_ice (dt,           ndtd,           &
    1866             :                       ncat,         n_aero,         &
    1867             :                       nilyr,        nslyr,          &
    1868             :                       ntrcr,        hin_max,        &
    1869             :                       rdg_conv,     rdg_shear,      &
    1870           0 :                       aicen,                        &
    1871           0 :                       trcrn,                        &
    1872           0 :                       vicen,        vsnon,          &
    1873             :                       aice0,                        &
    1874           0 :                       trcr_depend,                  &
    1875           0 :                       trcr_base,                    &
    1876           0 :                       n_trcr_strata,                &
    1877           0 :                       nt_strata,                    &
    1878             :                       krdg_partic,  krdg_redist,    &
    1879             :                       mu_rdg,       tr_brine,       &
    1880             :                       dardg1dt,     dardg2dt,       &
    1881             :                       dvirdgdt,     opening,        &
    1882             :                       fpond,                        &
    1883             :                       fresh,        fhocn,          &
    1884           0 :                       faero_ocn,    l_fiso_ocn,     &
    1885           0 :                       aparticn,     krdgn,          &
    1886           0 :                       aredistn,     vredistn,       &
    1887           0 :                       dardg1ndt,    dardg2ndt,      &
    1888           0 :                       dvirdgndt,                    &
    1889           0 :                       araftn,       vraftn,         &
    1890             :                       l_closing_flag,               &
    1891     1100016 :                       l_closing )
    1892     1100016 :       if (icepack_warnings_aborted(subname)) return
    1893             : 
    1894             :       !-----------------------------------------------------------------
    1895             :       ! ITD cleanup: Rebin thickness categories if necessary, and remove
    1896             :       !  categories with very small areas.
    1897             :       !-----------------------------------------------------------------
    1898             : 
    1899     1100016 :       dtt = dt * ndtd  ! for proper averaging over thermo timestep
    1900             :       call cleanup_itd (dtt,                  ntrcr,            &
    1901             :                         nilyr,                nslyr,            &
    1902             :                         ncat,                 hin_max,          &
    1903           0 :                         aicen,                trcrn,            &
    1904           0 :                         vicen,                vsnon,            &
    1905             :                         aice0,                aice,             &          
    1906             :                         n_aero,                                 &
    1907             :                         nbtrcr,               nblyr,            &
    1908             :                         tr_aero,                                &
    1909             :                         tr_pond_topo,         heat_capacity,    &  
    1910           0 :                         first_ice,                              &                
    1911           0 :                         trcr_depend,          trcr_base,        &
    1912           0 :                         n_trcr_strata,        nt_strata,        &
    1913             :                         fpond,                fresh,            &
    1914             :                         fsalt,                fhocn,            &
    1915           0 :                         faero_ocn,            l_fiso_ocn,       &
    1916             :                         fzsal,            &
    1917     1100016 :                         flux_bio)
    1918     1100016 :       if (icepack_warnings_aborted(subname)) return
    1919             : 
    1920     4400064 :       if (present(fiso_ocn)) fiso_ocn = l_fiso_ocn
    1921     1100016 :       deallocate(l_fiso_ocn)
    1922             : 
    1923     1100016 :       end subroutine icepack_step_ridge
    1924             : 
    1925             : !=======================================================================
    1926             : 
    1927             :       end module icepack_mechred
    1928             : 
    1929             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd