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

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! Vertical salinity (trcrn(nt_bgc_S)) is solved on the bio grid (bgrid and igrid)
       4             : ! with domain defined by the dynamic brine height (trcrn(nt_fbri) * vicen/aicen).
       5             : ! The CICE Bitz and Lipscomb thermodynamics is solved on the cgrid with height
       6             : ! vicen/aicen.
       7             : ! Gravity drainage is parameterized as nonlinear advection
       8             : ! Flushing is incorporated in the boundary changes and a darcy flow. 
       9             : ! (see Jeffery et al., JGR, 2011).  
      10             : !
      11             : ! authors: Nicole Jeffery, LANL
      12             : !          Elizabeth C. Hunke, LANL
      13             : !
      14             :       module icepack_zsalinity
      15             : 
      16             :       use icepack_kinds
      17             :       use icepack_parameters, only: c0, c1, c2, p001, p5, puny, rhow, depressT, gravit
      18             :       use icepack_parameters, only: rhosi, min_salin, salt_loss
      19             :       use icepack_parameters, only: l_skS, grid_oS, l_sk
      20             :       use icepack_parameters, only: dts_b
      21             :       use icepack_tracers, only: nt_sice
      22             :       use icepack_zbgc_shared, only: remap_zbgc
      23             :       use icepack_zbgc_shared, only: Ra_c, k_o, viscos_dynamic, thinS, Dm, exp_h
      24             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      25             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      26             : 
      27             :       use icepack_brine, only: calculate_drho
      28             :       use icepack_therm_shared, only: calculate_Tin_from_qin
      29             : 
      30             :       implicit none
      31             : 
      32             :       private
      33             :       public :: zsalinity
      34             : 
      35             :       real (kind=dbl_kind), parameter :: & 
      36             :          max_salin = 200.0_dbl_kind, & !(ppt) maximum bulk salinity
      37             :          lapidus_g = 0.3_dbl_kind      ! constant for artificial 
      38             :                                        ! viscosity/diffusion during growth
      39             : !        lapidus_m = 0.007_dbl_kind    ! constant for artificial diffusion during melt
      40             : 
      41             : !=======================================================================
      42             : 
      43             :       contains
      44             : 
      45             : !=======================================================================
      46             : 
      47           0 :       subroutine zsalinity (n_cat,              dt,           &
      48           0 :                             nilyr,              bgrid,        &
      49           0 :                             cgrid,              igrid,        &
      50           0 :                             trcrn_S,            trcrn_q,      &
      51           0 :                             trcrn_Si,           ntrcr,        &
      52             :                             fbri,                             &
      53           0 :                             bSin,               bTin,         &
      54           0 :                             bphin,              iphin,        &
      55           0 :                             ikin,               hbr_old,      &
      56             :                             hbrin,              hin,          &
      57           0 :                             hin_old,            iDin,         &
      58           0 :                             darcy_V,            brine_sal,    &
      59           0 :                             brine_rho,          ibrine_sal,   &
      60           0 :                             ibrine_rho,         dh_direct,    &
      61             :                             Rayleigh_criteria,                &
      62             :                             first_ice,          sss,          &
      63             :                             sst,                dh_top,       &
      64             :                             dh_bot,                           &
      65             :                             fzsal,                            &
      66             :                             fzsal_g,            bphi_min,     &
      67             :                             nblyr,              vicen,        &
      68             :                             aicen,              zsal_tot) 
      69             :              
      70             :       integer (kind=int_kind), intent(in) :: &
      71             :          nilyr         , & ! number of ice layers
      72             :          nblyr         , & ! number of bio layers
      73             :          ntrcr         , & ! number of tracers
      74             :          n_cat             ! category number 
      75             :                     
      76             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
      77             :          bgrid              ! biology nondimensional vertical grid points
      78             : 
      79             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
      80             :          igrid              ! biology vertical interface points
      81             :  
      82             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
      83             :          cgrid              ! CICE vertical coordinate   
      84             : 
      85             :       real (kind=dbl_kind), intent(in) :: &
      86             :          sss           , & ! ocean salinity (ppt)
      87             :          sst           , & ! ocean temperature (oC)
      88             :          hin_old       , & ! old ice thickness (m)
      89             :          dh_top        , & ! brine change in top and bottom for diagnostics (m)
      90             :          dh_bot        , & ! minimum porosity
      91             :          darcy_V       , & ! darcy velocity (m/s)
      92             :          dt            , & ! time step
      93             :          fbri          , & ! ratio of brine height to ice thickness
      94             :          hbr_old       , & ! old brine height  (m)
      95             :          hin           , & ! new ice thickness (m)
      96             :          hbrin         , & ! new brine height  (m)
      97             :          vicen         , & ! ice volume (m)
      98             :          aicen         , & ! ice area (m)
      99             :          bphi_min      , & !
     100             :          dh_direct         ! flooded or runoff amount (m)
     101             : 
     102             :       real (kind=dbl_kind), intent(inout) :: &
     103             :          zsal_tot      , & ! tot salinity (psu*rhosi*total vol ice)
     104             :          fzsal         , & ! total flux of salt out of ice over timestep(kg/m^2/s)
     105             :          fzsal_g           ! gravity drainage flux of salt over timestep(kg/m^2/s)
     106             : 
     107             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     108             :          bTin          , & ! Ice Temperature ^oC (on bio grid)
     109             :          bphin             ! Ice porosity (on bio grid)
     110             : 
     111             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     112             :          bSin          , & ! Ice salinity ppt (on bio  grid)
     113             :          brine_sal     , & ! brine salinity (ppt)
     114             :          brine_rho         ! brine density  (kg/m^3)
     115             : 
     116             :       real (kind=dbl_kind), dimension (nblyr), &
     117             :          intent(inout) :: &
     118             :          trcrn_S           ! salinity tracer ppt (on bio grid)
     119             : 
     120             :       real (kind=dbl_kind), dimension (nilyr), &
     121             :          intent(inout) :: &
     122             :          trcrn_q       , & ! enthalpy tracer 
     123             :          trcrn_Si          ! salinity on CICE grid
     124             : 
     125             :       logical (kind=log_kind), intent(inout) :: &
     126             :          Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c)  was reached 
     127             :       
     128             :       logical (kind=log_kind), intent(in) :: &
     129             :          first_ice         ! for first category ice only .true. 
     130             :                            !initialized values should be used 
     131             : 
     132             :       real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
     133             :          iDin          , & ! Diffusivity on the igrid   (1/s)
     134             :          ikin              ! permeability on the igrid 
     135             : 
     136             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
     137             :          iphin         , & ! porosity on the igrid 
     138             :          ibrine_rho    , & ! brine rho on interface  
     139             :          ibrine_sal        ! brine sal on interface   
     140             :          
     141             :       ! local variables
     142             : 
     143             :       real (kind=dbl_kind) :: &
     144           0 :          fzsaln        , & ! category flux of salt out of ice over timestep(kg/m^2/s)
     145           0 :          fzsaln_g      , & ! category gravity drainage flux of salt over timestep(kg/m^2/s)
     146           0 :          zsal_totn         ! total salt content
     147             : 
     148             :       character(len=*),parameter :: subname='(zsalinity)'
     149             : 
     150             :          call solve_zsalinity     (nilyr, nblyr, n_cat, dt,          &
     151             :                                    bgrid,      cgrid,    igrid,      &
     152             :                                    trcrn_S,            trcrn_q,      &
     153             :                                    trcrn_Si,           ntrcr,        &
     154             :                                    bSin,               bTin,         &
     155             :                                    bphin,              iphin,        &
     156             :                                    ikin,               hbr_old,      &
     157             :                                    hbrin,              hin,          &
     158             :                                    hin_old,            iDin,         &
     159             :                                    darcy_V,            brine_sal,    &
     160             :                                    brine_rho,          ibrine_sal,   &
     161             :                                    ibrine_rho,         dh_direct,    &
     162             :                                    Rayleigh_criteria,                &
     163             :                                    first_ice,          sss,          &
     164             :                                    sst,                dh_top,       &
     165             :                                    dh_bot,                           &
     166             :                                    fzsaln,                           &
     167           0 :                                    fzsaln_g,           bphi_min)
     168           0 :          if (icepack_warnings_aborted(subname)) return
     169             : 
     170           0 :          zsal_totn = c0
     171             : 
     172             :          call column_sum_zsal    (zsal_totn,          nblyr,      &
     173             :                                   vicen,              trcrn_S,    &
     174           0 :                                   fbri)
     175           0 :          if (icepack_warnings_aborted(subname)) return
     176             : 
     177             :          call merge_zsal_fluxes (aicen, &
     178             :                                  zsal_totn,          zsal_tot,       &
     179             :                                  fzsal,              fzsaln,         &
     180           0 :                                  fzsal_g,            fzsaln_g)            
     181           0 :          if (icepack_warnings_aborted(subname)) return
     182             : 
     183             :       end subroutine zsalinity
     184             : 
     185             : !=======================================================================
     186             : !
     187             : ! update vertical salinity 
     188             : ! 
     189           0 :       subroutine solve_zsalinity  (nilyr,              nblyr,        &
     190             :                                    n_cat,              dt,           &
     191           0 :                                    bgrid,    cgrid,    igrid,        &
     192           0 :                                    trcrn_S,            trcrn_q,      &
     193           0 :                                    trcrn_Si,           ntrcr,        &
     194           0 :                                    bSin,               bTin,         &
     195           0 :                                    bphin,              iphin,        &
     196           0 :                                    ikin,               hbr_old,      &
     197             :                                    hbrin,              hin,          &
     198           0 :                                    hin_old,            iDin,         &
     199           0 :                                    darcy_V,            brine_sal,    &
     200           0 :                                    brine_rho,          ibrine_sal,   &
     201           0 :                                    ibrine_rho,         dh_direct,    &
     202             :                                    Rayleigh_criteria,                &
     203             :                                    first_ice,          sss,          &
     204             :                                    sst,                dh_top,       &
     205             :                                    dh_bot,                           &  
     206             :                                    fzsaln,                           &
     207             :                                    fzsaln_g,           bphi_min)
     208             : 
     209             :       integer (kind=int_kind), intent(in) :: &
     210             :          nilyr,          & ! number of ice layers
     211             :          nblyr,          & ! number of bio layers
     212             :          ntrcr,          & ! number of tracers
     213             :          n_cat             ! category number 
     214             :                     
     215             :       real (kind=dbl_kind), intent(in) :: &
     216             :          dt                ! time step
     217             : 
     218             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
     219             :          bgrid              ! biology nondimensional vertical grid points
     220             : 
     221             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
     222             :          igrid              ! biology vertical interface points
     223             :  
     224             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
     225             :          cgrid              ! CICE vertical coordinate   
     226             : 
     227             :       real (kind=dbl_kind), intent(in) :: &
     228             :          sss           , & ! ocean salinity (ppt)
     229             :          sst           , & ! ocean temperature (oC)
     230             :          hin_old       , & ! old ice thickness (m)
     231             :          dh_top        , & ! brine change in top and bottom for diagnostics (m)
     232             :          dh_bot        , &
     233             :          darcy_V    
     234             : 
     235             :       real (kind=dbl_kind), intent(in) :: &
     236             :          hbr_old       , & ! old brine height  (m)
     237             :          hin           , & ! new ice thickness (m)
     238             :          hbrin         , & ! new brine height  (m)
     239             :          bphi_min      , & !
     240             :          dh_direct         ! flooded or runoff amount (m)
     241             :  
     242             :       real (kind=dbl_kind), intent(out) :: &
     243             :          fzsaln        , & ! total flux of salt out of ice over timestep(kg/m^2/s)
     244             :          fzsaln_g          ! gravity drainage flux of salt  over timestep(kg/m^2/s)
     245             : 
     246             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     247             :          bTin          , & ! Ice Temperature ^oC (on bio grid)
     248             :          bphin             ! Ice porosity (on bio grid)
     249             : 
     250             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     251             :          bSin          , & ! Ice salinity ppt (on bio  grid)
     252             :          brine_sal     , & ! brine salinity (ppt)
     253             :          brine_rho         ! brine density  (kg/m^3)
     254             : 
     255             :       real (kind=dbl_kind), dimension (nblyr), &
     256             :          intent(inout) :: &
     257             :          trcrn_S           ! salinity tracer ppt (on bio grid)
     258             : 
     259             :       real (kind=dbl_kind), dimension (nilyr), &
     260             :          intent(inout) :: &
     261             :          trcrn_q       , & ! enthalpy tracer 
     262             :          trcrn_Si          ! salinity on CICE grid
     263             : 
     264             :       logical (kind=log_kind), intent(inout) :: &
     265             :          Rayleigh_criteria ! .true. if minimun ice thickness (Ra_c)  was reached 
     266             :       
     267             :       logical (kind=log_kind), intent(in) :: &
     268             :          first_ice         ! for first category ice only .true. 
     269             :                            !initialized values should be used 
     270             : 
     271             :       real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
     272             :          iDin          , & ! Diffusivity on the igrid   (1/s)
     273             :          ikin              ! permeability on the igrid 
     274             : 
     275             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
     276             :          iphin         , & ! porosity on the igrid 
     277             :          ibrine_rho    , & ! brine rho on interface  
     278             :          ibrine_sal        ! brine sal on interface   
     279             :          
     280             :       ! local variables
     281             : 
     282             :       integer (kind=int_kind) :: &
     283             :          k, nint        ! vertical biology layer index 
     284             : 
     285             :       real (kind=dbl_kind) :: &
     286           0 :          surface_S         ! salinity of ice above hin > hbrin
     287             :       
     288             :       real (kind=dbl_kind), dimension(2) :: &
     289           0 :          S_bot           
     290             : 
     291             :       real (kind=dbl_kind) :: &
     292           0 :          Tmlts         , & ! melting temperature
     293           0 :          dts               ! local timestep (s)
     294             : 
     295             :       logical (kind=log_kind) :: &
     296             :          Rayleigh
     297             :      
     298             :       real (kind=dbl_kind):: &
     299           0 :          Ttemp             ! initial temp profile on the CICE grid
     300             : 
     301             :       real (kind=dbl_kind), dimension (ntrcr+2) :: &
     302           0 :          trtmp0        , & ! temporary, remapped tracers     !need extra 
     303           0 :          trtmp             ! temporary, remapped tracers     !
     304             : 
     305             :       logical (kind=log_kind) :: &
     306             :          cflag
     307             : 
     308             :       character(len=*),parameter :: subname='(solve_zsalinity)'
     309             : 
     310             :       !-----------------------------------------------------------------
     311             :       ! Initialize
     312             :       !-----------------------------------------------------------------
     313             : 
     314           0 :       dts = dts_b
     315           0 :       nint = max(1,INT(dt/dts))  
     316           0 :       dts = dt/nint
     317             : 
     318             :       !----------------------------------------------------------------
     319             :       ! Update boundary conditions
     320             :       !----------------------------------------------------------------
     321             :          
     322           0 :       surface_S = min_salin
     323             : 
     324           0 :       Rayleigh = .true.
     325           0 :       if (n_cat == 1 .AND. hbr_old < Ra_c) then
     326           0 :          Rayleigh = Rayleigh_criteria ! only category 1 ice can be false 
     327             :       endif
     328             : 
     329           0 :       if (dh_bot + darcy_V*dt > c0) then  
     330             :             
     331           0 :          bSin     (nblyr+2) = sss
     332           0 :          bTin     (nblyr+2) = sst
     333           0 :          brine_sal(nblyr+2) = sss 
     334           0 :          brine_rho(nblyr+2) = rhow
     335           0 :          bphin    (nblyr+2) = c1 
     336           0 :          S_bot    (1)       = c0
     337           0 :          S_bot    (2)       = c1  
     338             :    
     339             :       ! bottom melt
     340             :       else  
     341           0 :          bSin (nblyr+2) = bSin(nblyr+1)  
     342           0 :          Tmlts          =  -bSin(nblyr+2)* depressT 
     343           0 :          bTin (nblyr+2) = bTin(nblyr+1)
     344           0 :          bphin(nblyr+2) = iphin(nblyr+1)
     345           0 :          S_bot(1)       = c1
     346           0 :          S_bot(2)       = c0
     347             :       endif
     348             : 
     349           0 :       if (abs(dh_top) > puny .AND. abs(darcy_V) > puny) then 
     350           0 :          bSin(1) = max(min_salin,-(brine_rho(2)*brine_sal(2)/rhosi &
     351             :                  * darcy_V*dt - (dh_top + darcy_V*dt/bphi_min - dh_direct)*min_salin &
     352           0 :                  + max(c0,-dh_direct) * sss  )/dh_top)
     353           0 :          brine_sal(1) = brine_sal(2)
     354           0 :          brine_rho(1) = brine_rho(2)
     355           0 :          bphin(1) = bphi_min
     356             :       else
     357           0 :          bSin(1) = min_salin
     358             :       endif
     359             :             
     360             :       !-----------------------------------------------------------------
     361             :       ! Solve for S using CICE T.  If solve_zsal = .true., then couple back
     362             :       ! to the thermodynamics
     363             :       !-----------------------------------------------------------------
     364             : 
     365             :       call solve_S_dt (cflag,         nblyr,         &
     366             :                        nint         , dts          , &
     367             :                        bSin         , bTin         , &
     368             :                        bphin        , iphin        , &
     369             :                        igrid        , bgrid        , &
     370             :                        ikin         ,                &
     371             :                        hbr_old      , hbrin        , &
     372             :                        hin_old      , &
     373             :                        iDin         , darcy_V      , &
     374             :                        brine_sal    , Rayleigh     , &
     375             :                        first_ice    , sss          , &
     376             :                        dt           , dh_top       , &
     377             :                        dh_bot       , brine_rho    , &
     378             :                        ibrine_sal   , ibrine_rho   , &
     379             :                        fzsaln       , fzsaln_g     , &
     380           0 :                        S_bot        )
     381           0 :       if (icepack_warnings_aborted(subname)) return
     382             :   
     383           0 :       if (n_cat == 1)   Rayleigh_criteria = Rayleigh
     384             : 
     385           0 :       trtmp0(:) = c0
     386           0 :       trtmp (:) = c0
     387             :        
     388           0 :       do k = 1,nblyr                  ! back to bulk quantity 
     389           0 :             trcrn_S(k) =   bSin(k+1) 
     390           0 :             trtmp0(nt_sice+k-1) = trcrn_S(k)
     391             :       enddo           ! k
     392             : 
     393             :       call remap_zbgc   (nilyr, &
     394             :                          nt_sice,          &
     395           0 :                          trtmp0(1:ntrcr),  &
     396           0 :                          trtmp,            &
     397             :                          1,         nblyr, &
     398             :                          hin,       hbrin, &
     399           0 :                          cgrid(2:nilyr+1), &
     400           0 :                          bgrid(2:nblyr+1), &
     401           0 :                          surface_S )
     402           0 :       if (icepack_warnings_aborted(subname)) return
     403             :                
     404           0 :       do k = 1, nilyr
     405           0 :             Tmlts = -trcrn_Si(k)*depressT
     406             :             Ttemp = min(-(min_salin+puny)*depressT, &
     407           0 :                        calculate_Tin_from_qin(trcrn_q(k),Tmlts)) 
     408           0 :             trcrn_Si(k) = min(-Ttemp/depressT, max(min_salin, &
     409           0 :                                  trtmp(nt_sice+k-1)))
     410           0 :             Tmlts = - trcrn_Si(k)*depressT 
     411             :            ! if (cflag)  trcrn_q(k) = calculate_qin_from_Sin(Ttemp,Tmlts)  
     412             :       enddo ! k
     413             : 
     414             :       end subroutine solve_zsalinity
     415             : 
     416             : !=======================================================================
     417             : !
     418             : !  solves salt continuity explicitly using 
     419             : !  Lax-Wendroff-type scheme (MacCormack)
     420             : !  (Mendez-Nunez and Carroll,  Monthly Weather Review, 1993)
     421             : !
     422             : ! authors     Nicole Jeffery, LANL
     423             : !
     424           0 :       subroutine solve_S_dt          (cflag,  nblyr, nint,          &
     425           0 :                                       dts,    bSin,  bTin,          &
     426           0 :                                       bphin,  iphin, igrid,         &
     427           0 :                                       bgrid,  ikin,  hbri_old,      &
     428             :                                       hbrin,  hice_old,      &
     429           0 :                                       iDin,          darcy_V,       &
     430           0 :                                       brine_sal,     Rayleigh,      &
     431             :                                       first_ice,     sss,           &
     432             :                                       dt,            dht,           &
     433           0 :                                       dhb,           brine_rho,     &
     434           0 :                                       ibrine_sal,    ibrine_rho,    &
     435             :                                       fzsaln,        fzsaln_g,      &
     436             :                                       S_bot          )
     437             : 
     438             :       integer (kind=int_kind), intent(in) :: &
     439             :          nblyr            , & ! number of bio layers
     440             :          nint                 ! number of interations
     441             : 
     442             :       logical (kind=log_kind), intent(out) :: &
     443             :          cflag                ! thin or not
     444             : 
     445             :       real (kind=dbl_kind), intent(in) :: &
     446             :          dt               , & ! timestep (s)
     447             :          dts              , & ! local timestep (s)
     448             :          sss              , & ! sea surface salinity
     449             :          dht              , & ! change in the ice top  (positive for melting)
     450             :          dhb              , & ! change in the ice bottom (positive for freezing)
     451             :          hice_old         , & ! old ice thickness (m)
     452             :          hbri_old         , & ! brine thickness (m) 
     453             :          hbrin            , & ! new brine thickness (m)
     454             :          darcy_V              ! Darcy velocity due to a pressure head (m/s) or melt      
     455             : 
     456             :       real (kind=dbl_kind), intent(out) :: &
     457             :          fzsaln           , & ! salt flux +ive to  ocean (kg/m^2/s)
     458             :          fzsaln_g             ! gravity drainage salt flux +ive to ocean (kg/m^2/s)
     459             : 
     460             :       logical (kind=log_kind), intent(inout) :: &
     461             :          Rayleigh             ! if .true. convection is allowed; if .false. not yet
     462             : 
     463             :       logical (kind=log_kind), intent(in) :: &
     464             :          first_ice
     465             : 
     466             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
     467             :          brine_sal        , & ! Internal brine salinity (ppt)
     468             :          brine_rho        , & ! Internal brine density (kg/m^3)
     469             :          bgrid            , & ! biology nondimensional grid layer points 
     470             :          bTin                 ! Temperature of ice layers on bio grid for history (C) 
     471             : 
     472             :       real (kind=dbl_kind), dimension (nblyr+2), intent(inout) :: &
     473             :          bphin            , & ! Porosity of layers
     474             :          bSin                 ! Bulk Salinity (ppt) contains previous timestep
     475             :                               ! and ocean ss
     476             : 
     477             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
     478             :          ibrine_rho       , & ! brine rho on interface 
     479             :          ibrine_sal       , & ! brine sal on interface 
     480             :          igrid                ! biology grid interface points 
     481             : 
     482             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
     483             :          iphin                ! Porosity of layers on interface
     484             : 
     485             :       real (kind=dbl_kind), dimension (nblyr+1), intent(out) :: &
     486             :          iDin             , & ! Diffusivity on the igrid (1/s) with minimum bphi condition
     487             :          ikin                 ! permeability on interface
     488             : 
     489             :       real (kind=dbl_kind), dimension (2), intent(in) :: &
     490             :          S_bot
     491             :        
     492             :       ! local variables
     493             : 
     494             :       integer (kind=int_kind) :: &
     495             :          k, m       ! vertical biology layer index 
     496             : 
     497             :       real (kind=dbl_kind), dimension (nblyr+1) :: &
     498           0 :          iDin_p           , & ! Diffusivity on the igrid (1/s)/bphi^3 
     499           0 :          dSbdx            , & ! gradient of brine rho on grid
     500           0 :          drho             , & ! brine difference rho_a-rho_b  (kg/m^3)
     501             : !        Ci_s             , & !
     502           0 :          Ui_s             , & ! interface function
     503             : !        Vi_s             , & ! for conservation check
     504           0 :          ivel
     505             : 
     506             :       real (kind=dbl_kind), dimension (nblyr+2) :: &
     507           0 :          Din_p            , & ! Diffusivity on the igrid (1/s)/bphi^3  
     508           0 :          Sintemp          , & ! initial salinity
     509           0 :          pre_sin          , & ! estimate of  salinity of layers
     510           0 :          pre_sinb         , & ! estimate of  salinity of layers
     511           0 :          bgrid_temp       , & ! biology nondimensional grid layer points 
     512             :                               ! with boundary values 
     513           0 :          Q_s, C_s         , & ! Functions in continuity equation
     514           0 :          V_s, U_s, F_s    
     515             : 
     516             :       real (kind=dbl_kind) :: &
     517           0 :          dh               , & ! (m) change in hbrine over dts
     518           0 :          dbgrid           , & ! ratio of grid space to spacing across boundary 
     519             :                               ! i.e. 1/nilyr/(dbgrid(2)-dbgrid(1))
     520           0 :          lapidus          , & ! artificial viscosity:  use lapidus_g for growth
     521           0 :          Ssum_old,Ssum_new, & ! depth integrated salt before and after timestep
     522           0 :          fluxcorr,          & ! flux correction to prevent S < min_salin
     523           0 :          Ssum_corr,         & ! numerical boundary flux correction
     524           0 :          fluxb            , & ! bottom, top and total boundary flux (g/kg/m^2)
     525           0 :          fluxg            , & ! bottom, top and total gravity drainage flux (g/kg/m^2)
     526           0 :          fluxm            , & ! bottom, top and total molecular diffusion flux (g/kg/m^2)
     527           0 :          sum_old,sum_new  , & ! integrated salinity at t and t+dt
     528           0 :          dh_dt, dS_dt     , &
     529           0 :          Ssum_tmp
     530             : 
     531             :       real (kind=dbl_kind), dimension (nblyr) :: &
     532           0 :          vel              , & ! advective velocity times dt (m)
     533           0 :          lapidus_diff     , & ! lapidus term and 
     534           0 :          flux_corr        , &
     535           0 :          lapA             , &
     536           0 :          lapB    
     537             : 
     538             :       logical (kind=log_kind) :: &   
     539             :          test_conservation       ! test that salt change is balanced by fluxes 
     540             : 
     541             :       character(len=*),parameter :: subname='(solve_S_dt)'
     542             : 
     543             :       !-----------------------------------------------------------------
     544             :       !  Initialize
     545             :       !-----------------------------------------------------------------
     546             : 
     547           0 :       cflag = .false.
     548           0 :       test_conservation = .false. 
     549           0 :       iDin_p(:) = c0   
     550           0 :       Din_p(:) = c0 
     551           0 :       lapA(:) = c1
     552           0 :       lapB(:) = c1
     553           0 :       lapA(nblyr) = c0
     554           0 :       lapB(1) = c0
     555           0 :       V_s(:) = c0
     556           0 :       U_s(:) = c0
     557           0 :       Q_s(:) = c0
     558           0 :       C_s(:) = c0
     559             : !     Ci_s(:) = c0
     560           0 :       F_s(:) = c0
     561           0 :       Ui_s(:) = c0
     562             : !     Vi_s(:) = c0
     563           0 :       ivel(:) = c0
     564           0 :       vel(:) = c0
     565           0 :       dh = c0
     566           0 :       dbgrid = c2
     567           0 :       fzsaln = c0
     568           0 :       fzsaln_g = c0
     569             : 
     570             :       !-----------------------------------------------------------------
     571             :       ! Find brine density gradient for gravity drainage parameterization
     572             :       !-----------------------------------------------------------------
     573             : 
     574             :       call calculate_drho(nblyr, igrid, bgrid,&
     575           0 :                           brine_rho, ibrine_rho, drho)   
     576           0 :       if (icepack_warnings_aborted(subname)) return
     577             : 
     578             :       !-----------------------------------------------------------------
     579             :       ! Calculate bphi diffusivity on the grid points
     580             :       ! rhosi = 919-974 kg/m^2  set in bio_in
     581             :       ! rhow = 1026.0 density of sea water: uses kinematic viscosity (m^2/s) in Q18
     582             :       ! dynamic viscosity  divided by density = kinematic. 
     583             :       !-----------------------------------------------------------------
     584             : 
     585           0 :       do k = 2, nblyr+1
     586           0 :          iDin_p(k) = k_o*gravit*l_skS/viscos_dynamic*drho(k)/(hbri_old**2)
     587           0 :           Din_p(k) = (iDin_p(k)*(igrid(k)-bgrid(k)) &
     588           0 :                    + iDin_p(k-1)*(bgrid(k)-igrid(k-1)))/(igrid(k)-igrid(k-1))
     589             :       enddo                !k
     590             : 
     591             :       !-----------------------------------------------------------------
     592             :       ! Critical Ra_c value is only for the onset of convection in thinS 
     593             :       ! ice and not throughout, therefore I need a flag to indicate the 
     594             :       ! Ra_c was reached for a particular ice block.
     595             :       ! Using a thickness minimum (Ra_c) for simplicity.
     596             :       !-----------------------------------------------------------------
     597             : 
     598           0 :       bgrid_temp(:) = bgrid(:)
     599           0 :       Din_p(nblyr+2) =  iDin_p(nblyr+1)
     600           0 :       if (.NOT. Rayleigh .AND. hbrin < Ra_c) then
     601           0 :          Din_p(:) =  c0  
     602           0 :          iDin_p(:) = c0  
     603             :       else
     604           0 :          Rayleigh = .true.
     605             :       endif
     606             : 
     607             :       if (hbri_old > thinS .AND. hbrin > thinS .and. &
     608           0 :           hice_old > thinS .AND. .NOT. first_ice) then
     609             : 
     610           0 :          cflag = .true.
     611             : 
     612           0 :          bgrid_temp(1) = c0
     613           0 :          bgrid_temp(nblyr+2) = c1
     614           0 :          dbgrid = igrid(2)/(bgrid_temp(2)-bgrid_temp(1))
     615             : 
     616             :          !-----------------------------------
     617             :          ! surface boundary terms
     618             :          !-----------------------------------
     619             :                              
     620           0 :          lapidus = lapidus_g/real(nblyr,kind=dbl_kind)**2
     621           0 :          ivel(1) = dht/hbri_old
     622           0 :          U_s (1) = ivel(1)/dt*dts 
     623           0 :          Ui_s(1) = U_s(1) 
     624             : !        Ci_s(1) = c0
     625           0 :          F_s (1) = brine_rho(2)*brine_sal(2)/rhosi*darcy_V*dts/hbri_old/bSin(1)
     626             : 
     627             :          !-----------------------------------
     628             :          ! bottom boundary terms
     629             :          !-----------------------------------
     630             : 
     631           0 :          ivel(nblyr+1) =  dhb/hbri_old           
     632           0 :          Ui_s(nblyr+1) = ivel(nblyr+1)/dt*dts  
     633           0 :          dSbdx(nblyr) = (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) &
     634           0 :                       -  ibrine_sal(nblyr)*ibrine_rho(nblyr)) &
     635           0 :                       / (igrid(nblyr+1)-igrid(nblyr))
     636           0 :          C_s(nblyr+1) = Dm/brine_sal(nblyr+1)/brine_rho(nblyr+1)*dts/hbri_old**2 &
     637           0 :                       * (ibrine_sal(nblyr+1)*ibrine_rho(nblyr+1) &
     638           0 :                       -  ibrine_sal(nblyr)*ibrine_rho(nblyr)) &
     639           0 :                       / (igrid(nblyr+1)-igrid(nblyr))
     640           0 :          F_s(nblyr+1) = darcy_V*dts/hbri_old/bphin(nblyr+1)
     641           0 :          F_s(nblyr+2) = darcy_V*dts/hbri_old/bphin(nblyr+2)  
     642           0 :          vel(nblyr) =(bgrid(nblyr+1)*(dhb) -(bgrid(nblyr+1) - c1)*dht)/hbri_old
     643           0 :          U_s(nblyr+1) = vel(nblyr)/dt*dts  
     644           0 :          V_s(nblyr+1) = Din_p(nblyr+1)/rhosi &
     645           0 :                       * (rhosi/brine_sal(nblyr+1)/brine_rho(nblyr+1))**exp_h &
     646           0 :                       * dts*dSbdx(nblyr) 
     647           0 :          dSbdx(nblyr+1) = (brine_sal(nblyr+2)*brine_rho(nblyr+2) &
     648           0 :                         -  brine_sal(nblyr+1)*brine_rho(nblyr+1)) &
     649           0 :                         / (bgrid(nblyr+2)-bgrid(nblyr+1)+ grid_oS/hbri_old )  
     650           0 :          C_s( nblyr+2) = Dm/brine_sal(nblyr+2)/brine_rho(nblyr+2)*dts/hbri_old**2 &
     651           0 :                        * (brine_sal(nblyr+2)*brine_rho(nblyr+2) & 
     652           0 :                        -  brine_sal(nblyr+1)*brine_rho(nblyr+1)) &
     653           0 :                        / (bgrid(nblyr+2)-bgrid(nblyr+1) + grid_oS/hbri_old ) 
     654           0 :          U_s(nblyr+2) = ivel(nblyr+1)/dt*dts 
     655           0 :          V_s(nblyr+2) = Din_p(nblyr+2)/rhosi &
     656           0 :                       * (bphin(nblyr+1)/bSin(nblyr+2))**exp_h &
     657           0 :                       * dts*dSbdx(nblyr+1)
     658             : !        Ci_s(nblyr+1) = C_s(nblyr+2)
     659             : !        Vi_s(nblyr+1) = V_s(nblyr+2) 
     660           0 :          dh = (dhb-dht)/dt*dts
     661             : 
     662           0 :          do k = 2, nblyr  
     663           0 :             ivel(k) =  (igrid(k)*dhb - (igrid(k)-c1)*dht)/hbri_old
     664           0 :             Ui_s(k) = ivel(k)/dt*dts   
     665             : !           Vi_s(k) = iDin_p(k)/rhosi &
     666             : !                   *(rhosi/ibrine_rho(k)/ibrine_sal(k))**exp_h*dts &
     667             : !                   * (brine_sal(k+1)*brine_rho(k+1) &
     668             : !                   -  brine_sal(k)*brine_rho(k)) &
     669             : !                   / (bgrid(k+1)-bgrid(k)) 
     670           0 :             dSbdx(k-1) = (ibrine_sal(k)*ibrine_rho(k) &
     671           0 :                        -  ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1))
     672           0 :             F_s(k) = darcy_V*dts/hbri_old/bphin(k)
     673           0 :             C_s(k) = Dm/brine_sal(k)/brine_rho(k)*dts/hbri_old**2 &
     674           0 :                    * (ibrine_sal(k)*ibrine_rho(k) &
     675           0 :                    -  ibrine_sal(k-1)*ibrine_rho(k-1))/(igrid(k)-igrid(k-1))
     676             : !           Ci_s(k) = Dm/ibrine_sal(k)/ibrine_rho(k)*dts/hbri_old**2 &
     677             : !                   * (brine_sal(k+1)*brine_rho(k+1) &
     678             : !                   -  brine_sal(k)*brine_rho(k))/(bgrid(k+1)-bgrid(k))
     679           0 :             vel(k-1) = (bgrid(k)*(dhb) - (bgrid(k) - c1)* dht)/hbri_old
     680           0 :             U_s(k) = vel(k-1)/dt*dts 
     681           0 :             V_s(k) = Din_p(k)/rhosi &
     682           0 :                    * (rhosi/brine_sal(k)/brine_rho(k))**exp_h*dts*dSbdx(k-1) 
     683           0 :             C_s(2) = c0
     684           0 :             V_s(2) = c0
     685             :          enddo !k
     686             : 
     687             :       !-----------------------------------------------------------------
     688             :       ! Solve
     689             :       !-----------------------------------------------------------------
     690             : 
     691           0 :          do m = 1, nint     
     692             : 
     693           0 :             Sintemp(:) = bSin(:)
     694           0 :             pre_sin(:) = bSin(:)  
     695           0 :             pre_sinb(:) = bSin(:)
     696           0 :             Ssum_old = bSin(nblyr+1)*(igrid(nblyr+1)-igrid(nblyr))
     697             : 
     698             :             ! forward-difference 
     699             : 
     700           0 :             do k = 2, nblyr
     701           0 :                Ssum_old = Ssum_old + bSin(k)*(igrid(k)-igrid(k-1))
     702             : 
     703           0 :                pre_sin(k) =bSin(k) + (Ui_s(k)*(bSin(k+1) - bSin(k)) + &
     704           0 :                  V_s(k+1)*bSin(k+1)**3 - V_s(k)*bSin(k)**3 + &
     705           0 :                  (C_s(k+1)+F_s(k+1))*bSin(k+1)-&
     706           0 :                  (C_s(k)+F_s(k))*bSin(k))/(bgrid_temp(k+1)-bgrid_temp(k)) 
     707             :             enddo    !k
     708             : 
     709           0 :             pre_sin(nblyr+1) = bSin(nblyr+1) + (Ui_s(nblyr+1)*(bSin(nblyr+2) - &
     710           0 :                  bSin(nblyr+1)) +  V_s(nblyr+2)*bSin(nblyr+2)**3 - &
     711           0 :                  V_s(nblyr+1)*bSin(nblyr+1)**3+ (C_s(nblyr+2)+F_s(nblyr+2))*&
     712           0 :                  bSin(nblyr+2)-(C_s(nblyr+1)+F_s(nblyr+1))*bSin(nblyr+1) )/&
     713           0 :                  (bgrid_temp(nblyr+2)- bgrid_temp(nblyr+1))
     714             :               
     715             :             ! backward-difference 
     716             : 
     717           0 :             pre_sinb(2) = p5*(bSin(2) + pre_sin(2) +  (Ui_s(1)&
     718           0 :                   *(pre_sin(2) -pre_sin(1)) + &
     719           0 :                   V_s(2)*pre_sin(2)**3 - &
     720           0 :                   V_s(1)*pre_sin(1)**3 + (C_s(2)+F_s(2))*pre_sin(2)-&
     721           0 :                   (C_s(1)+F_s(1))*pre_sin(1) )/(bgrid_temp(2)-bgrid_temp(1)) )
     722             :               
     723           0 :             do k = nblyr+1, 3, -1  !nblyr+1
     724           0 :                pre_sinb(k) =p5*(bSin(k) + pre_sin(k) +  (Ui_s(k-1)&
     725           0 :                   *(pre_sin(k) - pre_sin(k-1)) + &
     726           0 :                   V_s(k)*pre_sin(k)**3 - &
     727           0 :                   V_s(k-1)*pre_sin(k-1)**3 + (C_s(k)+F_s(k))*pre_sin(k)-&
     728           0 :                   (C_s(k-1)+F_s(k-1))*pre_sin(k-1))/(bgrid_temp(k)-bgrid_temp(k-1)) )
     729             : 
     730           0 :                Q_s(k) = V_s(k)*pre_sin(k)**2 + U_s(k) + C_s(k) + F_s(k) 
     731             :             enddo   !k
     732             : 
     733           0 :             Q_s(2) = V_s(2)*pre_sin(2)**2 + U_s(2) + C_s(2) + F_s(2)
     734           0 :             Q_s(1) = V_s(1)*pre_sin(2)**2 + Ui_s(1) + C_s(1)+ F_s(1)
     735           0 :             Q_s(nblyr+2) = V_s(nblyr+2)*pre_sin(nblyr+1)**2 + & 
     736           0 :             Ui_s(nblyr+1) + C_s(nblyr+2) +  F_s(nblyr+2)
     737             : 
     738             :       !-----------------------------------------------------------------
     739             :       ! Add artificial viscosity   [Lapidus,1967] [Lohner et al, 1985]
     740             :       ! * more for melting ice
     741             :       !-----------------------------------------------------------------
     742             : 
     743           0 :             lapidus_diff(:) = c0
     744           0 :             flux_corr(:) = c0
     745           0 :             Ssum_new = c0
     746           0 :             Ssum_corr = c0
     747           0 :             fluxcorr = c0
     748           0 :             fluxg = c0
     749           0 :             fluxb = c0
     750           0 :             fluxm = c0
     751             : 
     752           0 :             do k = 2, nblyr+1  
     753             : 
     754           0 :                lapidus_diff(k-1) =    lapidus/& ! lapidus/real(nblyr,kind=dbl_kind)**2/&
     755           0 :                   (igrid(k)-igrid(k-1))* &
     756           0 :                   ( lapA(k-1)*ABS(Q_s(k+1)-Q_s(k))*(abs(pre_sinb(k+1))-abs(pre_sinb(k)))/&
     757           0 :                   (bgrid_temp(k+1)-bgrid_temp(k) )**2 - &
     758           0 :                   lapB(k-1)*ABS(Q_s(k)-Q_s(k-1))*(abs(pre_sinb(k))-abs(pre_sinb(k-1)))/&
     759           0 :                   (bgrid_temp(k)-bgrid_temp(k-1))**2)
     760             :                           
     761           0 :                bSin(k) = pre_sinb(k) + lapidus_diff(k-1)
     762             : 
     763           0 :                if (bSin(k) < min_salin) then
     764           0 :                   flux_corr(k-1) = min_salin - bSin(k) !  flux into the ice
     765           0 :                   bSin(k) = min_salin 
     766           0 :                elseif (bSin(k) > -bTin(k)/depressT) then
     767           0 :                   flux_corr(k-1) = bSin(k)+bTin(k)/depressT !  flux into the ice
     768           0 :                   bSin(k) = -bTin(k)/depressT
     769           0 :                elseif (bSin(k) > max_salin) then
     770           0 :                   call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     771           0 :                   call icepack_warnings_add(subname//' bSin(k) > max_salin')
     772             :                endif
     773             :             
     774           0 :                if (k == nblyr+1) bSin(nblyr+2) = S_bot(1)*bSin(nblyr+1) &
     775           0 :                                                   + S_bot(2)*bSin(nblyr+2) 
     776             : 
     777           0 :                Ssum_new = Ssum_new + bSin(k)*(igrid(k)-igrid(k-1))
     778           0 :                fluxcorr = fluxcorr + (flux_corr(k-1) &
     779           0 :                         + lapidus_diff(k-1))*(igrid(k)-igrid(k-1))
     780             : 
     781             :             enddo   !k
     782             : 
     783           0 :             Ssum_tmp = Ssum_old
     784             : 
     785             :             call calc_salt_fluxes (nint, nblyr, &
     786             :                  Ui_s,   dh,dbgrid,hbri_old,Sintemp,    &
     787             :                  pre_sin,   fluxb,fluxg,fluxm,V_s,    &
     788             :                  C_s,   F_s,   Ssum_corr,fzsaln_g,fzsaln, &
     789           0 :                  Ssum_tmp,dts, Ssum_new)
     790           0 :             if (icepack_warnings_aborted(subname)) return
     791             : 
     792           0 :             if (test_conservation) then
     793             :                call check_conserve_salt(nint, m, dt, &
     794             :                                 Ssum_tmp, Ssum_new, Ssum_corr,&
     795             :                                 fluxcorr, fluxb, fluxg, fluxm, &
     796           0 :                                 hbrin, hbri_old)
     797           0 :                call icepack_warnings_add(subname//' check_conserve_salt fails')
     798           0 :                if (icepack_warnings_aborted(subname)) return
     799             :             endif  ! test_conservation
     800             : 
     801             :          enddo !m
     802             : 
     803             :       else    !  add/melt ice only 
     804             : 
     805           0 :          sum_old = c0
     806           0 :          sum_new = c0
     807           0 :          dh_dt = hbrin-hbri_old
     808           0 :          dS_dt = c0
     809           0 :          if (dh_dt > c0) then 
     810           0 :             dS_dt = sss*dh_dt*salt_loss
     811           0 :             do k = 2, nblyr+1 
     812           0 :                bSin(k) = max(min_salin,(bSin(k)*hbri_old + dS_dt)/hbrin)
     813             :             enddo  !k
     814             :          else
     815           0 :             do k = 2, nblyr+1 
     816           0 :                sum_old = sum_old + bSin(k)*hbri_old*(igrid(k)-igrid(k-1))
     817           0 :                bSin(k) = max(min_salin,bSin(k) * (c1-abs(dh_dt)/hbri_old))
     818           0 :                sum_new = sum_new + bSin(k)*hbrin*(igrid(k)-igrid(k-1))
     819             :             enddo  !k
     820             :          endif
     821           0 :          fzsaln = rhosi*(sum_old-sum_new + dS_dt)*p001/dt   !kg/m^2/s
     822           0 :          fzsaln_g = c0
     823             : 
     824             :       endif  ! (hbri_old > thinS .AND. hbrin > thinS &
     825             :              ! .and.  hice_old > thinS .AND. .NOT. first_ice) 
     826             : 
     827             :       !-----------------------------------------------------------------
     828             :       ! Move this to bgc calculation if using tr_salinity
     829             :       ! Calculate bphin, iphin, ikin, iDin and iDin_N 
     830             :       !-----------------------------------------------------------------
     831             : 
     832           0 :       iDin(:) = c0
     833           0 :       iphin(:)  = c1
     834           0 :       ikin(:) = c0    
     835             : 
     836           0 :       do k = 1, nblyr+1
     837           0 :          if (k < nblyr+1) bphin(k+1) = min(c1,max(puny, &
     838           0 :                           bSin(k+1)*rhosi/(brine_sal(k+1)*brine_rho(k+1)))) 
     839           0 :          if (k == 1) then     
     840           0 :             bphin(k) = min(c1,max(puny, bSin(k)*rhosi/(brine_sal(k)*brine_rho(k))))  
     841           0 :             iphin(k) = bphin(2)
     842           0 :          elseif (k == nblyr+1) then
     843           0 :             iphin(nblyr+1) = bphin(nblyr+1)
     844             :          else
     845           0 :             iphin(k) = min(c1, max(c0,(bphin(k+1) - bphin(k))/(bgrid(k+1) &
     846           0 :                            - bgrid(k))*(igrid(k)-bgrid(k)) + bphin(k)))
     847             :          endif 
     848           0 :          ikin(k) = k_o*iphin(k)**exp_h 
     849             :       enddo    !k
     850             : 
     851           0 :       if (cflag) then
     852             :                
     853           0 :          do k = 2, nblyr+1
     854           0 :             iDin(k) =  iphin(k)*Dm/hbri_old**2  
     855           0 :             if (Rayleigh .AND. hbrin .GE. Ra_c) &
     856           0 :             iDin(k) = iDin(k) + l_sk*ikin(k)*gravit/viscos_dynamic &
     857           0 :                          * drho(k)/hbri_old**2 
     858             :          enddo       !k
     859             :       else   !  .not. cflag
     860           0 :          do k = 2, nblyr+1
     861           0 :             iDin(k) = iphin(k)*Dm/hbri_old**2
     862             :          enddo       !k
     863             :       endif
     864             : 
     865             :       end subroutine solve_S_dt
     866             : 
     867             : !=======================================================================
     868             : !
     869             : ! Calculate salt fluxes
     870             : ! 
     871           0 :       subroutine calc_salt_fluxes (mint, nblyr, &
     872           0 :                                    Ui_s,dh,dbgrid,hbri_old,Sintemp,pre_sin,&
     873           0 :                                    fluxb,fluxg,fluxm,V_s,&
     874           0 :                                    C_s,F_s,Ssum_corr,fzsaln_g,fzsaln,Ssum_old, &
     875             : !                                  fluxcorr,dts, Ssum_new)
     876             :                                    dts, Ssum_new)
     877             : 
     878             :       integer(kind=int_kind), intent(in) :: &
     879             :          nblyr,          & ! number of bio layers
     880             :          mint              ! current iteration
     881             : 
     882             :       real (kind=dbl_kind), intent(in) :: &
     883             :          dts      , & !  halodynamic timesteps (s)
     884             :         ! hbrin    , & ! new brine height after all iterations (m)
     885             :          dh       , &  ! (m) change in hbrine over dts
     886             :          dbgrid   , & ! ratio of grid space to spacing across boundary 
     887             :                       ! ie. 1/nilyr/(dbgrid(2)-dbgrid(1))
     888             : !        fluxcorr , & ! flux correction to ensure S >= min_salin
     889             :          hbri_old     ! initial brine height (m)
     890             : 
     891             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: & 
     892             :          Ui_s         ! interface function
     893             : 
     894             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
     895             :          Sintemp  , & ! initial salinity
     896             :          pre_sin  , & ! estimate of salinity of layers
     897             :          C_s      , & ! Functions in continuity equation
     898             :          F_s      , &
     899             :          V_s
     900             : 
     901             :       real (kind=dbl_kind), intent(in) :: &
     902             :          Ssum_old , & ! initial integrated salt content (ppt)/h  
     903             :          Ssum_new     ! next integrated salt content(ppt)/h
     904             : 
     905             :       real (kind=dbl_kind), intent(inout) :: &
     906             :          fzsaln   , & ! total salt flux  out of ice over timestep(kg/m^2/s)
     907             :          fzsaln_g , & ! gravity drainage flux of salt  over timestep(kg/m^2/s)
     908             :          Ssum_corr, & ! boundary flux correction due to numerics
     909             :          fluxb    , & ! total boundary salt flux into the ice (+ into ice)
     910             :          fluxg    , & ! total gravity drainage salt flux into the ice (+ into ice)
     911             :          fluxm        ! total molecular diffusive salt flux into the ice (+ into ice)
     912             : 
     913             :       ! local variables
     914             : 
     915             :       real (kind=dbl_kind) :: &
     916           0 :          Ssum_corr_flux  , & ! numerical boundary flux correction
     917           0 :          fluxb_b, fluxb_t, & ! bottom, top and total boundary flux (g/kg/m^2)
     918           0 :          fluxg_b, fluxg_t, & ! bottom, top and total gravity drainage flux (g/kg/m^2)
     919           0 :          fluxm_b, fluxm_t    ! bottom, top and total molecular diffusion flux (g/kg/m^2)
     920             : 
     921           0 :       real (kind=dbl_kind) :: hin_old, hin_next, dhtmp !, dh
     922             : 
     923             :       character(len=*),parameter :: subname='(calc_salt_fluxes)'
     924             : 
     925           0 :       dhtmp = c1-dh/hbri_old
     926           0 :       hin_next = hbri_old +  real(mint,kind=dbl_kind)*dh
     927           0 :       hin_old  = hbri_old + (real(mint,kind=dbl_kind)-c1)*dh
     928             : 
     929             :       !-----------------------------------------------------------------
     930             :       ! boundary fluxes (positive into the ice)
     931             :       !---------------------------------------------
     932             :       !  without higher order numerics corrections 
     933             :       ! fluxb = (Ui_s(nblyr+1) + F_s(nblyr+2))*Sintemp(nblyr+2)  - (Ui_s(1) + F_s(1))*Sintemp(1)
     934             :       !-----------------------------------------------------------------
     935             : 
     936           0 :       fluxb_b = p5* Ui_s(nblyr+1) * (dhtmp*Sintemp(nblyr+2)*dbgrid &
     937           0 :                                   +        pre_sin(nblyr+1) &
     938           0 :                                   +  dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)) &
     939           0 :               + p5*(F_s(nblyr+2)  *  dhtmp*Sintemp(nblyr+2)*dbgrid &
     940           0 :                   + F_s(nblyr+1)  *       (pre_sin(nblyr+1) &
     941           0 :                                   +  dhtmp*Sintemp(nblyr+1)*(c1-dbgrid)))
     942             : 
     943           0 :       fluxb_t = -p5*Ui_s(1)*(pre_sin(1)*dbgrid +  &
     944           0 :                               dhtmp*Sintemp(2) -  &
     945           0 :                               (dbgrid-c1)*pre_sin(2)) - &
     946           0 :                               p5*(dbgrid*F_s(1)*pre_sin(1) + &
     947           0 :                               F_s(2)*(dhtmp*Sintemp(2) &
     948           0 :                               +(c1-dbgrid)*pre_sin(2)))
     949             : 
     950           0 :       fluxb = fluxb_b + fluxb_t 
     951             : 
     952             :       !-----------------------------------------------------------------
     953             :       ! gravity drainage fluxes (positive into the ice)
     954             :       ! without higher order numerics corrections 
     955             :       ! fluxg =  V_s(nblyr+2)*Sintemp(nblyr+1)**3
     956             :       !-----------------------------------------------------------------
     957             : 
     958             :       fluxg_b =  p5*(dhtmp* dbgrid* &
     959           0 :                                V_s(nblyr+2)*Sintemp(nblyr+1)**3  +  &
     960           0 :                                V_s(nblyr+1)*(pre_sin(nblyr+1)**3 - &
     961             :                                dhtmp*(dbgrid - c1)* &
     962           0 :                                Sintemp(nblyr+1)**3))
     963             :                 
     964           0 :       fluxg_t =  -p5*(dbgrid*V_s(1)*pre_sin(1)**3 + &
     965           0 :                                V_s(2)*(dhtmp*Sintemp(2)**3- &
     966           0 :                                (dbgrid-c1)*pre_sin(2)**3))
     967             :                 
     968           0 :       fluxg =  fluxg_b + fluxg_t
     969             :       
     970             :       !-----------------------------------------------------------------
     971             :       ! diffusion fluxes (positive into the ice)
     972             :       ! without higher order numerics corrections 
     973             :       ! fluxm = C_s(nblyr+2)*Sintemp(nblyr+2)
     974             :       !-----------------------------------------------------------------
     975             : 
     976           0 :       fluxm_b = p5*(dhtmp*C_s(nblyr+2)* Sintemp(nblyr+2)*dbgrid &
     977           0 :                        +  C_s(nblyr+1)*(pre_sin(nblyr+1) &
     978           0 :                   + dhtmp * Sintemp(nblyr+1)*(c1-dbgrid)))
     979             : 
     980           0 :       fluxm_t = -p5 * (C_s(1) *  pre_sin(1)*dbgrid &
     981           0 :               +        C_s(2) * (pre_sin(2)*(c1-dbgrid) + dhtmp*Sintemp(2)))
     982             :            
     983           0 :       fluxm =  fluxm_b + fluxm_t
     984             :               
     985           0 :       Ssum_corr      = (-dh/hbri_old + p5*(dh/hbri_old)**2)*Ssum_old 
     986           0 :       Ssum_corr_flux = dh*Ssum_old/hin_next + Ssum_corr
     987           0 :       Ssum_corr      = Ssum_corr_flux
     988             :           
     989             :       fzsaln_g = fzsaln_g - hin_next *  fluxg_b &
     990           0 :                           * rhosi*p001/dts
     991             :        
     992             :       !approximate fluxes
     993             :       !fzsaln   = fzsaln   - hin_next * (fluxg &
     994             :       !                    + fluxb + fluxm + fluxcorr + Ssum_corr_flux) & 
     995             :       !                    * rhosi*p001/dts  
     996             : 
     997             :       fzsaln   = fzsaln + (Ssum_old*hin_old - Ssum_new*hin_next) &
     998           0 :                           * rhosi*p001/dts  ! positive into the ocean
     999             : 
    1000           0 :       end subroutine calc_salt_fluxes
    1001             : 
    1002             : !=======================================================================
    1003             : ! 
    1004             : ! Test salt conservation:   flux conservative form d(hSin)/dt = -dF(x,Sin)/dx 
    1005             : !  
    1006           0 :       subroutine check_conserve_salt (mmax, mint,     dt, &
    1007             :                                       Ssum_old, Ssum_new, Ssum_corr,        & 
    1008             :                                       fluxcorr, fluxb,    fluxg,     fluxm, &
    1009             :                                       hbrin,    hbri_old)
    1010             : 
    1011             :       integer(kind=int_kind), intent(in) :: &
    1012             :          mint      , & ! current iteration
    1013             :          mmax          ! maximum number of iterations
    1014             : 
    1015             :       real (kind=dbl_kind), intent(in) :: &
    1016             :          dt             , &  ! thermodynamic and halodynamic timesteps (s)
    1017             :          hbrin          , &  ! (m) final brine height 
    1018             :          hbri_old       , &  ! (m) initial brine height
    1019             :          Ssum_old       , &  ! initial integrated salt content
    1020             :          Ssum_new       , &  ! final integrated salt content
    1021             :          fluxcorr       , &  ! flux correction to ensure S >= min_salin
    1022             :          Ssum_corr      , &  ! boundary flux correction due to numerics
    1023             :          fluxb          , &  ! total boundary salt flux into the ice (+ into ice)
    1024             :          fluxg          , &  ! total gravity drainage salt flux into the ice (+ into ice)
    1025             :          fluxm               ! 
    1026             : 
    1027             :      ! local variables
    1028             : 
    1029             :      real (kind=dbl_kind):: &
    1030           0 :          diff2     , & !
    1031           0 :          dsum_flux , & ! salt change in kg/m^2/s
    1032           0 :          flux_tot  , & ! fluxg + fluxb
    1033           0 :          order     , & !
    1034           0 :          dh
    1035             : 
    1036             :      real (kind=dbl_kind), parameter :: &
    1037             :          accuracy = 1.0e-7_dbl_kind ! g/kg/m^2/s difference between boundary fluxes 
    1038             : 
    1039             :      character(len=*),parameter :: subname='(check_conserve_salt)'
    1040             : 
    1041           0 :          dh = (hbrin-hbri_old)/real(mmax,kind=dbl_kind)
    1042             : 
    1043             :          flux_tot = (fluxb + fluxg + fluxm + fluxcorr + Ssum_corr)*&
    1044           0 :                     (hbri_old + (real(mint,kind=dbl_kind))*dh)/dt
    1045             :          dsum_flux =(Ssum_new*(hbri_old + (real(mint,kind=dbl_kind))*dh) - &
    1046             :                      Ssum_old*(hbri_old + (real(mint,kind=dbl_kind)-c1)* &
    1047           0 :                      dh) )/dt
    1048           0 :          order = abs(dh/min(hbri_old,hbrin))
    1049           0 :          if (abs(dsum_flux) > accuracy) then
    1050           0 :            diff2 = abs(dsum_flux - flux_tot)
    1051           0 :            if (diff2 >  puny .AND. diff2 > order ) then 
    1052           0 :               call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1053           0 :               write(warnstr,*) subname, 'Poor salt conservation: check_conserve_salt'
    1054           0 :               call icepack_warnings_add(warnstr)
    1055           0 :               write(warnstr,*) subname, 'mint:', mint
    1056           0 :               call icepack_warnings_add(warnstr)
    1057           0 :               write(warnstr,*) subname, 'Ssum_corr',Ssum_corr
    1058           0 :               call icepack_warnings_add(warnstr)
    1059           0 :               write(warnstr,*) subname, 'fluxb,fluxg,fluxm,flux_tot,fluxcorr:'
    1060           0 :               call icepack_warnings_add(warnstr)
    1061           0 :               write(warnstr,*) subname,  fluxb,fluxg,fluxm,flux_tot,fluxcorr
    1062           0 :               call icepack_warnings_add(warnstr)
    1063           0 :               write(warnstr,*) subname, 'fluxg,',fluxg
    1064           0 :               call icepack_warnings_add(warnstr)
    1065           0 :               write(warnstr,*) subname, 'dsum_flux,',dsum_flux
    1066           0 :               call icepack_warnings_add(warnstr)
    1067           0 :               write(warnstr,*) subname, 'Ssum_new,Ssum_old,hbri_old,dh:'
    1068           0 :               call icepack_warnings_add(warnstr)
    1069           0 :               write(warnstr,*) subname,  Ssum_new,Ssum_old,hbri_old,dh
    1070           0 :               call icepack_warnings_add(warnstr)
    1071           0 :               write(warnstr,*) subname, 'diff2,order,puny',diff2,order,puny
    1072           0 :               call icepack_warnings_add(warnstr)
    1073             :            endif
    1074             :          endif
    1075             : 
    1076           0 :      end subroutine check_conserve_salt
    1077             : 
    1078             : !=======================================================================
    1079             : !
    1080             : ! Aggregate flux information from all ice thickness categories
    1081             : !
    1082           0 :       subroutine merge_zsal_fluxes(aicenS,               &
    1083             :                                    zsal_totn,  zsal_tot, & 
    1084             :                                    fzsal,      fzsaln,   &
    1085             :                                    fzsal_g,    fzsaln_g)
    1086             : 
    1087             :       ! single category fluxes
    1088             :       real (kind=dbl_kind), intent(in):: &          
    1089             :           aicenS   , & ! concentration of ice
    1090             :           fzsaln  , &  ! salt flux                       (kg/m**2/s)
    1091             :           fzsaln_g     ! gravity drainage salt flux      (kg/m**2/s)
    1092             : 
    1093             :       real (kind=dbl_kind), intent(in):: &
    1094             :           zsal_totn   ! tot salinity in category (psu*volume*rhosi)
    1095             : 
    1096             :       real (kind=dbl_kind), intent(inout):: &          
    1097             :           zsal_tot, & ! tot salinity (psu*rhosi*total vol ice)
    1098             :           fzsal   , & ! salt flux                       (kg/m**2/s)
    1099             :           fzsal_g     ! gravity drainage salt flux      (kg/m**2/s)
    1100             : 
    1101             :       character(len=*),parameter :: subname='(merge_zsal_fluxes)'
    1102             : 
    1103             :       !-----------------------------------------------------------------
    1104             :       ! Merge fluxes
    1105             :       !-----------------------------------------------------------------
    1106             : 
    1107           0 :       zsal_tot  = zsal_tot + zsal_totn  ! already *aicenS
    1108             : 
    1109             :       ! ocean tot and gravity drainage salt fluxes
    1110           0 :       fzsal    = fzsal   + fzsaln   * aicenS
    1111           0 :       fzsal_g  = fzsal_g + fzsaln_g * aicenS
    1112             : 
    1113           0 :       end subroutine merge_zsal_fluxes
    1114             : 
    1115             : !==========================================================================
    1116             : !
    1117             : ! For each grid cell, sum field over all ice layers.  "Net" refers to  the column
    1118             : ! integration while "avg"  is normalized by the ice thickness
    1119             : 
    1120           0 :       subroutine column_sum_zsal (zsal_totn, nblyr,   &
    1121           0 :                                   vicenS, trcrn_S, fbri)
    1122             : 
    1123             :       integer (kind=int_kind), intent(in) :: &
    1124             :          nblyr         ! number of layers
    1125             : 
    1126             :       real (kind=dbl_kind),  intent(in):: &          
    1127             :          vicenS    , & ! volume of ice (m)
    1128             :          fbri          ! brine height to ice thickness ratio
    1129             : 
    1130             :       real (kind=dbl_kind), dimension (nblyr), intent(in) :: &
    1131             :          trcrn_S       ! input field
    1132             : 
    1133             :       real (kind=dbl_kind), intent(inout) :: &
    1134             :          zsal_totn    ! avg salinity (psu*rhosi*vol of ice)
    1135             : 
    1136             :       ! local variables
    1137             : 
    1138             :       integer (kind=int_kind) :: &
    1139             :            k           ! layer index
    1140             : 
    1141             :       character(len=*),parameter :: subname='(column_sum_zsal)'
    1142             : 
    1143           0 :       do k = 1, nblyr
    1144             :          zsal_totn = zsal_totn &
    1145           0 :                    + rhosi * trcrn_S(k) &
    1146             :                            * fbri  &
    1147           0 :                            * vicenS/real(nblyr,kind=dbl_kind)
    1148             :       enddo ! k
    1149             : 
    1150           0 :       end subroutine column_sum_zsal
    1151             : 
    1152             : !=======================================================================
    1153             : 
    1154             :       end module icepack_zsalinity
    1155             : 
    1156             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd