LCOV - code coverage report
Current view: top level - columnphysics - icepack_zbgc.F90 (source / functions) Hit Total Coverage
Test: 201120-001525:782a1b7d78:3:base,travis,quick Lines: 370 490 75.51 %
Date: 2020-11-19 17:37:37 Functions: 8 8 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : !
       3             : ! Biogeochemistry driver
       4             : !
       5             : ! authors: Nicole Jeffery, LANL
       6             : !          Scott Elliot,   LANL
       7             : !          Elizabeth C. Hunke, LANL
       8             : !
       9             :       module icepack_zbgc
      10             : 
      11             :       use icepack_kinds
      12             :       use icepack_parameters, only: c0, c1, c2, p001, p1, p5, puny
      13             :       use icepack_parameters, only: depressT, rhosi, min_salin, salt_loss
      14             :       use icepack_parameters, only: fr_resp, algal_vel, R_dFe2dust, dustFe_sol, T_max
      15             :       use icepack_parameters, only: op_dep_min, fr_graze_s, fr_graze_e, fr_mort2min, fr_dFe
      16             :       use icepack_parameters, only: k_nitrif, t_iron_conv, max_loss, max_dfe_doc1
      17             :       use icepack_parameters, only: fr_resp_s, y_sk_DMS, t_sk_conv, t_sk_ox
      18             :       use icepack_parameters, only: scale_bgc, ktherm, skl_bgc, solve_zsal
      19             :       use icepack_parameters, only: z_tracers, fsal, conserv_check
      20             : 
      21             :       use icepack_tracers, only: nt_sice, nt_bgc_S, bio_index 
      22             :       use icepack_tracers, only: tr_brine, nt_fbri, nt_qice, nt_Tsfc
      23             :       use icepack_tracers, only: nt_zbgc_frac
      24             :       use icepack_tracers, only: bio_index_o,  bio_index  
      25             : 
      26             :       use icepack_zbgc_shared, only: zbgc_init_frac
      27             :       use icepack_zbgc_shared, only: zbgc_frac_init
      28             :       use icepack_zbgc_shared, only: bgc_tracer_type, remap_zbgc
      29             :       use icepack_zbgc_shared, only: regrid_stationary
      30             :       use icepack_zbgc_shared, only: R_S2N, R_Si2N, R_Fe2C, R_Fe2N, R_Fe2DON, R_Fe2DOC
      31             :       use icepack_zbgc_shared, only: chlabs, alpha2max_low, beta2max
      32             :       use icepack_zbgc_shared, only: mu_max, grow_Tdep, fr_graze
      33             :       use icepack_zbgc_shared, only: mort_pre, mort_Tdep, k_exude
      34             :       use icepack_zbgc_shared, only: K_Nit, K_Am, K_Sil, K_Fe
      35             :       use icepack_zbgc_shared, only: f_don, kn_bac, f_don_Am, f_doc
      36             :       use icepack_zbgc_shared, only: f_exude, k_bac
      37             :       use icepack_zbgc_shared, only: tau_ret, tau_rel
      38             :       use icepack_zbgc_shared, only: R_C2N, R_CHL2N, f_abs_chl, R_C2N_DON
      39             : 
      40             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      41             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      42             : 
      43             :       use icepack_brine, only: preflushing_changes, compute_microS_mushy
      44             :       use icepack_brine, only: update_hbrine, compute_microS 
      45             :       use icepack_algae, only: zbio, sklbio
      46             :       use icepack_therm_shared, only: calculate_Tin_from_qin
      47             :       use icepack_itd, only: column_sum, column_conservation_check
      48             :       use icepack_zsalinity, only: zsalinity
      49             : 
      50             :       implicit none 
      51             : 
      52             :       private
      53             :       public :: add_new_ice_bgc, &
      54             :                 lateral_melt_bgc, &
      55             :                 icepack_init_bgc, &
      56             :                 icepack_init_zbgc, &
      57             :                 icepack_biogeochemistry, &
      58             :                 icepack_load_ocean_bio_array, &
      59             :                 icepack_init_ocean_bio
      60             : 
      61             : !=======================================================================
      62             : 
      63             :       contains
      64             : 
      65             : !=======================================================================
      66             : 
      67             : ! Adjust biogeochemical tracers when new frazil ice forms
      68             : 
      69       66636 :       subroutine add_new_ice_bgc (dt,         nblyr,                &
      70             :                                   ncat,       nilyr,      nltrcr,   &
      71       66636 :                                   bgrid,      cgrid,      igrid,    &
      72      133272 :                                   aicen_init, vicen_init, vi0_init, &
      73       66636 :                                   aicen,      vicen,      vsnon1,   &
      74             :                                   vi0new,                           &
      75       66636 :                                   ntrcr,      trcrn,      nbtrcr,   &
      76       66636 :                                   sss,        ocean_bio,  flux_bio, &
      77             :                                   hsurp)
      78             : 
      79             :       integer (kind=int_kind), intent(in) :: &
      80             :          nblyr   , & ! number of bio layers
      81             :          ncat     , & ! number of thickness categories
      82             :          nilyr    , & ! number of ice layers
      83             :          nltrcr, & ! number of zbgc tracers
      84             :          nbtrcr  , & ! number of biology tracers
      85             :          ntrcr       ! number of tracers in use
      86             : 
      87             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
      88             :          bgrid              ! biology nondimensional vertical grid points
      89             : 
      90             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
      91             :          igrid              ! biology vertical interface points
      92             :  
      93             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
      94             :          cgrid              ! CICE vertical coordinate   
      95             : 
      96             :       real (kind=dbl_kind), intent(in) :: &
      97             :          dt              ! time step (s)
      98             : 
      99             :       real (kind=dbl_kind), dimension (:), &
     100             :          intent(in) :: &
     101             :          aicen_init  , & ! initial concentration of ice
     102             :          vicen_init  , & ! intiial volume per unit area of ice  (m)
     103             :          aicen       , & ! concentration of ice
     104             :          vicen           ! volume per unit area of ice          (m)
     105             : 
     106             :       real (kind=dbl_kind), intent(in) :: &
     107             :          vsnon1          ! category 1 snow volume per unit area (m)
     108             : 
     109             :       real (kind=dbl_kind), dimension (:,:), &
     110             :          intent(inout) :: &
     111             :          trcrn           ! ice tracers
     112             : 
     113             :       real (kind=dbl_kind), intent(in) :: &
     114             :          sss              !sea surface salinity (ppt)
     115             : 
     116             :       real (kind=dbl_kind), intent(in) :: &
     117             :          vi0_init    , & ! volume of new ice added to cat 1 (intial)
     118             :          vi0new          ! volume of new ice added to cat 1
     119             : 
     120             :       real (kind=dbl_kind), intent(in) :: &
     121             :          hsurp           ! thickness of new ice added to each cat
     122             : 
     123             :       real (kind=dbl_kind), dimension (:), &
     124             :          intent(inout) :: &
     125             :          flux_bio   ! tracer flux to ocean from biology (mmol/m^2/s) 
     126             :         
     127             :       real (kind=dbl_kind), dimension (:), &
     128             :          intent(in) :: &
     129             :          ocean_bio       ! ocean concentration of biological tracer
     130             : 
     131             : ! local
     132             : 
     133             :       integer (kind=int_kind) :: &
     134             :          location    , & ! 1 (add frazil to bottom), 0 (add frazil throughout)
     135             :          n           , & ! ice category index
     136             :          k               ! ice layer index
     137             : 
     138             :       real (kind=dbl_kind) :: &
     139       20520 :          vbri1       , & ! starting volume of existing brine
     140       20520 :          vbri_init   , & ! brine volume summed over categories
     141       20520 :          vbri_final      ! brine volume summed over categories
     142             : 
     143             :       real (kind=dbl_kind) :: &
     144       20520 :          vsurp       , & ! volume of new ice added to each cat
     145       20520 :          vtmp            ! total volume of new and old ice
     146             :         
     147             :       real (kind=dbl_kind), dimension (ncat) :: &
     148      189756 :          vbrin           ! trcrn(nt_fbri,n)*vicen(n) 
     149             :        
     150             :       real (kind=dbl_kind) :: &
     151       20520 :          vice_new        ! vicen_init + vsurp
     152             : 
     153             :       real (kind=dbl_kind) :: &
     154       20520 :          Tmlts       ! melting temperature (oC)
     155             : 
     156             :       character (len=char_len) :: &
     157             :          fieldid         ! field identifier
     158             : 
     159             :       character(len=*),parameter :: subname='(add_new_ice_bgc)'
     160             : 
     161             :       !-----------------------------------------------------------------     
     162             :       ! brine
     163             :       !-----------------------------------------------------------------
     164      399816 :       vbrin(:) = c0
     165      399816 :       do n = 1, ncat
     166      333180 :          vbrin(n) = vicen_init(n)
     167      399816 :          if (tr_brine) vbrin(n) =  trcrn(nt_fbri,n)*vicen_init(n)
     168             :       enddo
     169             :      
     170       66636 :       call column_sum (ncat,  vbrin,  vbri_init)
     171       66636 :       if (icepack_warnings_aborted(subname)) return
     172             : 
     173       66636 :       vbri_init = vbri_init + vi0_init
     174     1374588 :       do k = 1, nbtrcr  
     175      385632 :          flux_bio(k) = flux_bio(k) &
     176     1374588 :                             - vi0_init/dt*ocean_bio(k)*zbgc_init_frac(k)
     177             :       enddo
     178             :       !-----------------------------------------------------------------
     179             :       ! Distribute bgc in new ice volume among all ice categories by 
     180             :       ! increasing ice thickness, leaving ice area unchanged.
     181             :       !-----------------------------------------------------------------
     182             : 
     183             :          ! Diffuse_bio handles concentration changes from ice growth/melt
     184             :          ! ice area does not change
     185             :          ! add salt to the bottom , location = 1 
     186             : 
     187       66636 :        vsurp = c0
     188       66636 :        vtmp = c0
     189             : 
     190      399816 :       do n = 1,ncat
     191             :  
     192      399816 :       if (hsurp > c0) then
     193             : 
     194         110 :          vtmp = vbrin(n)
     195         110 :          vsurp = hsurp * aicen_init(n) 
     196         110 :          vbrin(n) = vbrin(n) + vsurp
     197         110 :          vice_new = vicen_init(n) + vsurp
     198         110 :          if (tr_brine .and. vicen(n) > c0) then
     199          22 :             trcrn(nt_fbri,n) = vbrin(n)/vicen(n)
     200          88 :          elseif (tr_brine .and. vicen(n) <= c0) then
     201          88 :             trcrn(nt_fbri,n) = c1
     202             :          endif
     203             : 
     204         110 :          if (nltrcr > 0) then 
     205         110 :             location = 1  
     206             :             call adjust_tracer_profile(nbtrcr,   dt, ntrcr, &
     207          55 :                                        aicen_init(n),       &
     208          55 :                                        vbrin(n),            &
     209             :                                        vice_new,            &
     210           0 :                                        trcrn(:,n),          &
     211             :                                        vtmp,                &
     212             :                                        vsurp,        sss,   &
     213             :                                        nilyr,        nblyr, &
     214             :                                        solve_zsal,   bgrid, & 
     215             :                                        cgrid,               &
     216           0 :                                        ocean_bio,    igrid, &
     217         110 :                                        location)
     218         110 :             if (icepack_warnings_aborted(subname)) return
     219             :          endif       ! nltrcr       
     220             :       endif          ! hsurp > 0
     221             :       enddo          ! n
     222             : 
     223             :       !-----------------------------------------------------------------
     224             :       ! Combine bgc in new ice grown in open water with category 1 ice.
     225             :       !-----------------------------------------------------------------
     226             :        
     227       66636 :       if (vi0new > c0) then
     228             : 
     229         150 :          vbri1    = vbrin(1) 
     230         150 :          vbrin(1) = vbrin(1) + vi0new
     231         150 :          if (tr_brine .and. vicen(1) > c0) then
     232         150 :             trcrn(nt_fbri,1) = vbrin(1)/vicen(1)
     233           0 :          elseif (tr_brine .and. vicen(1) <= c0) then
     234           0 :             trcrn(nt_fbri,1) = c1
     235             :          endif
     236             :        
     237             :       ! Diffuse_bio handles concentration changes from ice growth/melt
     238             :       ! ice area changes
     239             :       ! add salt throughout, location = 0
     240             : 
     241         150 :          if (nltrcr > 0) then 
     242         133 :             location = 0  
     243             :             call adjust_tracer_profile(nbtrcr,  dt,     ntrcr,  &
     244          23 :                                        aicen(1),                &
     245          23 :                                        vbrin(1),                &
     246          23 :                                        vicen(1),                &
     247           0 :                                        trcrn(:,1),              &
     248             :                                        vbri1,                   &
     249             :                                        vi0new,          sss,    &
     250             :                                        nilyr,           nblyr,  &
     251             :                                        solve_zsal,      bgrid,  &
     252             :                                        cgrid,                   &
     253           0 :                                        ocean_bio,       igrid,  &
     254         133 :                                        location)
     255         133 :             if (icepack_warnings_aborted(subname)) return
     256             : 
     257         133 :             if (solve_zsal .and. vsnon1 .le. c0) then
     258           0 :                Tmlts = -trcrn(nt_sice,1)*depressT
     259           0 :                trcrn(nt_Tsfc,1) =  calculate_Tin_from_qin(trcrn(nt_qice,1),Tmlts)
     260             :             endif        ! solve_zsal 
     261             :          endif           ! nltrcr > 0
     262             :       endif              ! vi0new > 0
     263             : 
     264       66636 :       if (tr_brine .and. conserv_check) then
     265           0 :          call column_sum (ncat,   vbrin,  vbri_final)
     266           0 :          if (icepack_warnings_aborted(subname)) return
     267             : 
     268           0 :          fieldid = subname//':vbrin'
     269             :          call column_conservation_check (fieldid,                  &
     270             :                                          vbri_init, vbri_final,    &
     271           0 :                                          puny)
     272           0 :          if (icepack_warnings_aborted(subname)) return
     273             : 
     274             :       endif   ! conserv_check
     275             : 
     276             :       end subroutine add_new_ice_bgc
     277             : 
     278             : !=======================================================================
     279             : 
     280             : ! When sea ice melts laterally, flux bgc to ocean
     281             : 
     282       63007 :       subroutine lateral_melt_bgc (dt,                 &
     283             :                                    ncat,     nblyr,    &
     284       63007 :                                    rside,    vicen,    &
     285       63007 :                                    trcrn,    fzsal,    &
     286       63007 :                                    flux_bio, nbltrcr)
     287             : 
     288             :       integer (kind=int_kind), intent(in) :: &
     289             :          ncat  , & ! number of thickness categories
     290             :          nblyr , & ! number of bio layers
     291             :          nbltrcr   ! number of biology tracers
     292             : 
     293             :       real (kind=dbl_kind), intent(in) :: &
     294             :          dt        ! time step (s)
     295             : 
     296             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
     297             :          vicen     ! volume per unit area of ice          (m)
     298             : 
     299             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
     300             :          trcrn     ! tracer array
     301             : 
     302             :       real (kind=dbl_kind), intent(in) :: &
     303             :          rside     ! fraction of ice that melts laterally
     304             : 
     305             :       real (kind=dbl_kind), intent(inout) :: &
     306             :          fzsal     ! salt flux from layer Salinity (kg/m^2/s)
     307             :   
     308             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
     309             :          flux_bio  ! biology tracer flux from layer bgc (mmol/m^2/s)
     310             : 
     311             :       ! local variables
     312             : 
     313             :       integer (kind=int_kind) :: &
     314             :          k     , & ! layer index
     315             :          m     , & !
     316             :          n         ! category index
     317             : 
     318             :       real (kind=dbl_kind) :: &
     319       17020 :          zspace    ! bio grid spacing
     320             : 
     321             :       character(len=*),parameter :: subname='(lateral_melt_bgc)'
     322             : 
     323       63007 :       zspace = c1/(real(nblyr,kind=dbl_kind))
     324             : 
     325       63007 :       if (solve_zsal) then
     326           0 :          do n = 1, ncat
     327           0 :          do k = 1,nblyr
     328           0 :             fzsal = fzsal + rhosi*trcrn(nt_fbri,n) &
     329           0 :                   * vicen(n)*p001*zspace*trcrn(nt_bgc_S+k-1,n) &
     330           0 :                   * rside/dt
     331             :          enddo
     332             :          enddo
     333             :       endif
     334             : 
     335     1312279 :       do m = 1, nbltrcr
     336     7558639 :          do n = 1, ncat
     337    57466512 :          do k = 1, nblyr+1
     338    13181280 :             flux_bio(m) = flux_bio(m) + trcrn(nt_fbri,n) &
     339    26362560 :                         * vicen(n)*zspace*trcrn(bio_index(m)+k-1,n) &
     340    69398520 :                         * rside/dt
     341             :          enddo
     342             :          enddo
     343             :       enddo
     344             : 
     345       63007 :       end subroutine lateral_melt_bgc 
     346             : 
     347             : !=======================================================================
     348             : !
     349             : ! Add new ice tracers to the ice bottom and adjust the vertical profile 
     350             : !
     351             : ! author: Nicole Jeffery, LANL
     352             : 
     353         243 :       subroutine adjust_tracer_profile (nbtrcr, dt, ntrcr, &
     354             :                                         aicen,      vbrin, &
     355         243 :                                         vicen,      trcrn, &
     356             :                                         vtmp,              &
     357             :                                         vsurp,      sss,   &
     358             :                                         nilyr,      nblyr, &
     359         243 :                                         solve_zsal, bgrid, &
     360         486 :                                         cgrid,      ocean_bio, &
     361         243 :                                         igrid,      location)
     362             : 
     363             :       integer (kind=int_kind), intent(in) :: &
     364             :          location          , & ! 1 (add frazil to bottom), 0 (add frazil throughout)
     365             :          ntrcr             , & ! number of tracers in use
     366             :          nilyr             , & ! number of ice layers
     367             :          nbtrcr            , & ! number of biology tracers
     368             :          nblyr                 ! number of biology layers
     369             : 
     370             :       real (kind=dbl_kind), intent(in) :: &
     371             :          dt              ! timestep (s)
     372             : 
     373             :       real (kind=dbl_kind), intent(in) :: &
     374             :          aicen   , & ! concentration of ice
     375             :          vicen   , & ! volume of ice
     376             :          sss     , & ! ocean salinity (ppt)
     377             :        ! hsurp   , & ! flags new ice added to each cat
     378             :          vsurp   , & ! volume of new ice added to each cat
     379             :          vtmp        ! total volume of new and old ice
     380             : 
     381             :       real (kind=dbl_kind), dimension (nbtrcr), intent(in) :: &
     382             :          ocean_bio
     383             : 
     384             :       real (kind=dbl_kind), intent(in) :: &
     385             :          vbrin       ! fbri*volume per unit area of ice  (m)
     386             :        
     387             :       logical (kind=log_kind), intent(in) :: &
     388             :          solve_zsal 
     389             : 
     390             :       real (kind=dbl_kind), dimension (nblyr+1), intent(in) :: &
     391             :          igrid       ! zbio grid
     392             : 
     393             :       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
     394             :          bgrid       ! zsal grid
     395             : 
     396             :       real (kind=dbl_kind), dimension (nilyr+1), intent(in) :: &
     397             :          cgrid       ! CICE grid
     398             : 
     399             :       real (kind=dbl_kind), dimension (ntrcr), &
     400             :          intent(inout) :: &
     401             :          trcrn       ! ice tracers
     402             :       
     403             :       ! local variables
     404             : 
     405             :       real (kind=dbl_kind), dimension (ntrcr+2) :: &
     406       18692 :          trtmp0, &      ! temporary, remapped tracers
     407       18692 :          trtmp          ! temporary, remapped tracers
     408             :         
     409             :       real (kind=dbl_kind) :: &
     410          78 :          hin     , & ! ice height
     411          78 :          hinS_new, & ! brine height
     412          78 :          temp_S         
     413             : 
     414             :       integer (kind=int_kind) :: &
     415             :          k, m 
     416             : 
     417             :       real (kind=dbl_kind), dimension (nblyr+1) ::  &     
     418        1188 :          C_stationary      ! stationary bulk concentration*h (mmol/m^2)
     419             : 
     420             :       real (kind=dbl_kind), dimension (nblyr) ::  &     
     421         945 :          S_stationary      ! stationary bulk concentration*h (ppt*m)
     422             : 
     423             :       real(kind=dbl_kind) :: &
     424          78 :          top_conc     , & ! salinity or bgc ocean concentration of frazil
     425          78 :          fluxb        , & ! needed for regrid (set to zero here)
     426          78 :          hbri_old     , & ! previous timestep brine height
     427          78 :          hbri             ! brine height 
     428             : 
     429             :       character(len=*),parameter :: subname='(adjust_tracer_profile)'
     430             : 
     431       58787 :       trtmp0(:) = c0
     432       58787 :       trtmp(:) = c0
     433         243 :       fluxb = c0
     434             : 
     435         243 :       if (location == 1 .and. vbrin > c0) then  ! add frazil to bottom
     436             : 
     437          22 :          hbri     = vbrin
     438          22 :          hbri_old = vtmp
     439          22 :          if (solve_zsal) then
     440           0 :             top_conc = sss * salt_loss
     441           0 :             do k = 1, nblyr 
     442           0 :                S_stationary(k) = trcrn(nt_bgc_S+k-1)* hbri_old
     443             :             enddo
     444             :             call regrid_stationary (S_stationary, hbri_old, &
     445             :                                     hbri,         dt,       &
     446             :                                     ntrcr,                  &
     447             :                                     nblyr-1,      top_conc, &
     448           0 :                                     bgrid(2:nblyr+1), fluxb )
     449           0 :             if (icepack_warnings_aborted(subname)) return
     450           0 :             do k = 1, nblyr 
     451           0 :                trcrn(nt_bgc_S+k-1) =  S_stationary(k)/hbri
     452           0 :                trtmp0(nt_sice+k-1) = trcrn(nt_bgc_S+k-1)
     453             :             enddo
     454             :          endif  ! solve_zsal
     455             : 
     456         462 :          do m = 1, nbtrcr
     457         440 :             top_conc = ocean_bio(m)*zbgc_init_frac(m)
     458        3960 :             do k = 1, nblyr+1 
     459        3960 :                C_stationary(k) = trcrn(bio_index(m) + k-1)* hbri_old
     460             :             enddo !k
     461             :             call regrid_stationary (C_stationary, hbri_old, &
     462             :                                     hbri,         dt,       &
     463             :                                     ntrcr,                  &
     464             :                                     nblyr,        top_conc, &
     465         440 :                                     igrid,        fluxb )
     466         440 :             if (icepack_warnings_aborted(subname)) return
     467        3982 :             do k = 1, nblyr+1 
     468        3960 :                trcrn(bio_index(m) + k-1) =  C_stationary(k)/hbri
     469             :             enddo !k                  
     470             :          enddo !m
     471             : 
     472          22 :          if (solve_zsal) then
     473           0 :             if (aicen > c0) then
     474           0 :                hinS_new  = vbrin/aicen
     475           0 :                hin       = vicen/aicen
     476             :             else
     477           0 :                hinS_new  = c0
     478           0 :                hin       = c0
     479             :             endif                   ! aicen
     480           0 :             temp_S    = min_salin   ! bio to cice
     481             :             call remap_zbgc(nilyr,    &
     482             :                             nt_sice,                   &
     483           0 :                             trtmp0(1:ntrcr), trtmp,    &
     484             :                             1,               nblyr,    &
     485             :                             hin,             hinS_new, &
     486           0 :                             cgrid(2:nilyr+1),          &
     487           0 :                             bgrid(2:nblyr+1), temp_S   )
     488           0 :             if (icepack_warnings_aborted(subname)) return
     489          22 :             do k = 1, nilyr
     490           0 :                trcrn(nt_sice+k-1) = trtmp(nt_sice+k-1)   
     491             :             enddo        ! k
     492             :          endif           ! solve_zsal
     493             : 
     494         221 :       elseif (vbrin > c0) then   ! add frazil throughout  location == 0 .and.
     495             : 
     496        1197 :          do k = 1, nblyr+1
     497        1064 :             if (solve_zsal .and. k < nblyr + 1) then
     498           0 :                trcrn(nt_bgc_S+k-1) = (trcrn(nt_bgc_S+k-1) * vtmp &
     499           0 :                                           + sss*salt_loss * vsurp) / vbrin
     500           0 :                trtmp0(nt_sice+k-1) = trcrn(nt_bgc_S+k-1)
     501             :             endif                    ! solve_zsal
     502       21933 :             do m = 1, nbtrcr
     503        9408 :                trcrn(bio_index(m) + k-1) = (trcrn(bio_index(m) + k-1) * vtmp &
     504       28072 :                          + ocean_bio(m)*zbgc_init_frac(m) * vsurp) / vbrin
     505             :             enddo
     506             :          enddo
     507             : 
     508         133 :          if (solve_zsal) then
     509           0 :             if (aicen > c0) then
     510           0 :                hinS_new  = vbrin/aicen
     511           0 :                hin       = vicen/aicen
     512             :             else
     513           0 :                hinS_new  = c0
     514           0 :                hin       = c0
     515             :             endif              !aicen
     516           0 :             temp_S    = min_salin   ! bio to cice
     517             :             call remap_zbgc(nilyr,    &
     518             :                          nt_sice,                   &
     519           0 :                          trtmp0(1:ntrcr), trtmp,    &
     520             :                          1,               nblyr,    &
     521             :                          hin,             hinS_new, &
     522           0 :                          cgrid(2:nilyr+1),          &        
     523           0 :                          bgrid(2:nblyr+1),temp_S    )
     524           0 :             if (icepack_warnings_aborted(subname)) return
     525           0 :             do k = 1, nilyr
     526           0 :                trcrn(nt_sice+k-1) = trtmp(nt_sice+k-1)   
     527             :             enddo        !k
     528             :          endif   ! solve_zsal
     529             : 
     530             :       endif     ! location
     531             : 
     532             :       end subroutine adjust_tracer_profile
     533             : 
     534             : !=======================================================================
     535             : !autodocument_start icepack_init_bgc
     536             : !
     537             : 
     538          20 :       subroutine icepack_init_bgc(ncat, nblyr, nilyr, ntrcr_o, &
     539          20 :          cgrid, igrid, ntrcr, nbtrcr, &
     540          20 :          sicen, trcrn, sss, ocean_bio_all)
     541             : 
     542             :       integer (kind=int_kind), intent(in) :: &
     543             :          ncat  , & ! number of thickness categories
     544             :          nilyr , & ! number of ice layers
     545             :          nblyr , & ! number of bio layers
     546             :          ntrcr_o,& ! number of tracers not including bgc
     547             :          ntrcr , & ! number of tracers in use
     548             :          nbtrcr    ! number of bio tracers in use
     549             :  
     550             :       real (kind=dbl_kind), dimension (nblyr+1), intent(inout) :: &
     551             :          igrid     ! biology vertical interface points
     552             :  
     553             :       real (kind=dbl_kind), dimension (nilyr+1), intent(inout) :: &
     554             :          cgrid     ! CICE vertical coordinate   
     555             : 
     556             :       real (kind=dbl_kind), dimension(nilyr, ncat), intent(in) :: &
     557             :          sicen     ! salinity on the cice grid
     558             : 
     559             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
     560             :          trcrn     ! subset of tracer array (only bgc) 
     561             : 
     562             :       real (kind=dbl_kind), intent(in) :: &
     563             :          sss       ! sea surface salinity (ppt)
     564             : 
     565             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     566             :          ocean_bio_all   ! fixed order, all values even for tracers false
     567             : 
     568             : !autodocument_end
     569             : 
     570             :       ! local variables
     571             : 
     572             :       integer (kind=int_kind) :: &
     573             :          k     , & ! vertical index 
     574             :          n     , & ! category index 
     575             :          mm        ! bio tracer index
     576             : 
     577             :       real (kind=dbl_kind), dimension (ntrcr+2) :: & 
     578        1968 :          trtmp     ! temporary, remapped tracers   
     579             :       
     580             :       character(len=*),parameter :: subname='(icepack_init_bgc)'
     581             : 
     582             :       !-----------------------------------------------------------------------------   
     583             :       !     Skeletal Layer Model
     584             :       !  All bgc tracers are Bulk quantities in units of mmol or mg per m^3
     585             :       !  The skeletal layer model assumes a constant 
     586             :       !  layer depth (sk_l) and porosity (phi_sk)
     587             :       !-----------------------------------------------------------------------------   
     588          20 :          if (skl_bgc) then
     589             :        
     590          24 :             do  n = 1,ncat
     591         344 :             do mm = 1,nbtrcr
     592             :                ! bulk concentration (mmol or mg per m^3, or 10^-3 mmol/m^3)
     593         340 :                trcrn(bio_index(mm)-ntrcr_o, n) = ocean_bio_all(bio_index_o(mm))
     594             :             enddo       ! nbtrcr
     595             :             enddo       ! n 
     596             : 
     597             :       !-----------------------------------------------------------------------------   
     598             :       !    zbgc Model
     599             :       !  All bgc tracers are Bulk quantities in units of mmol or mg per m^3
     600             :       !  The vertical layer model uses prognosed porosity and layer depth
     601             :       !-----------------------------------------------------------------------------   
     602             : 
     603             :          else   ! not skl_bgc
     604             : 
     605          16 :             if (scale_bgc .and. solve_zsal) then ! bulk concentration (mmol or mg per m^3)
     606           0 :                do n = 1,ncat
     607           0 :                do mm = 1,nbtrcr
     608           0 :                   do k = 2, nblyr
     609           0 :                      trcrn(bio_index(mm)+k-1-ntrcr_o,n) = &
     610           0 :                           (p5*(trcrn(nt_bgc_S+k-1-ntrcr_o,n)+ trcrn(nt_bgc_S+k-2-ntrcr_o,n)) &
     611           0 :                          / sss*ocean_bio_all(bio_index_o(mm))) 
     612             :                   enddo  !k
     613           0 :                   trcrn(nt_zbgc_frac-1+mm-ntrcr_o,n) = zbgc_frac_init(mm)
     614           0 :                   trcrn(bio_index(mm)-ntrcr_o,n) = (trcrn(nt_bgc_S-ntrcr_o,n) &
     615           0 :                                          / sss*ocean_bio_all(bio_index_o(mm))) 
     616           0 :                   trcrn(bio_index(mm)+nblyr-ntrcr_o,n) = (trcrn(nt_bgc_S+nblyr-1-ntrcr_o,n) &
     617           0 :                                                / sss*ocean_bio_all(bio_index_o(mm)))
     618           0 :                   trcrn(bio_index(mm)+nblyr+1-ntrcr_o:bio_index(mm)+nblyr+2-ntrcr_o,n) = c0 ! snow
     619             :                enddo ! mm
     620             :                enddo ! n 
     621             :     
     622          16 :             elseif (scale_bgc .and. ktherm == 2) then
     623           0 :                trtmp(:) = c0
     624           0 :                do n = 1,ncat     
     625             :                   call remap_zbgc(nilyr,    &
     626             :                                   1,                          &
     627           0 :                                   sicen(:,n),       trtmp,    &
     628             :                                   0,                nblyr+1,  &
     629             :                                   c1,               c1,       &
     630           0 :                                   cgrid(2:nilyr+1),           &
     631           0 :                                   igrid(1:nblyr+1),           &
     632           0 :                                   sicen(1,n)                  )
     633           0 :                   if (icepack_warnings_aborted(subname)) return
     634             : 
     635           0 :                   do mm = 1,nbtrcr
     636           0 :                   do k = 1, nblyr + 1            
     637           0 :                      trcrn(bio_index(mm)+k-1-ntrcr_o,n) =   &
     638           0 :                           (trtmp(k)/sss*ocean_bio_all(bio_index_o(mm)))
     639           0 :                      trcrn(bio_index(mm)+nblyr+1-ntrcr_o:bio_index(mm)+nblyr+2-ntrcr_o,n) = c0 ! snow
     640             :                   enddo  ! k
     641             :                   enddo  ! mm
     642             :                enddo     ! n 
     643             : 
     644          16 :             elseif (nbtrcr > 0 .and. nt_fbri > 0) then ! not scale_bgc         
     645             :      
     646          96 :                do n = 1,ncat
     647        1616 :                do mm = 1,nbtrcr
     648       13680 :                do k = 1, nblyr+1
     649       17280 :                   trcrn(bio_index(mm)+k-1-ntrcr_o,n) = ocean_bio_all(bio_index_o(mm)) &
     650       23680 :                                              * zbgc_init_frac(mm) 
     651       38000 :                   trcrn(bio_index(mm)+nblyr+1-ntrcr_o:bio_index(mm)+nblyr+2-ntrcr_o,n) = c0 ! snow
     652             :                enddo    ! k
     653        1600 :                trcrn(nt_zbgc_frac-1+mm-ntrcr_o,n) = zbgc_frac_init(mm)
     654             :                enddo    ! mm
     655             :                enddo    ! n 
     656             :               
     657             :             endif  ! scale_bgc
     658             :          endif     ! skl_bgc
     659             : 
     660             :       end subroutine icepack_init_bgc
     661             : 
     662             : !=======================================================================
     663             : !autodocument_start icepack_init_zbgc
     664             : !
     665             : 
     666          90 :       subroutine icepack_init_zbgc ( &
     667         180 :                  R_Si2N_in, R_S2N_in, R_Fe2C_in, R_Fe2N_in, R_C2N_in, R_C2N_DON_in, &
     668         180 :                  R_chl2N_in, F_abs_chl_in, R_Fe2DON_in, R_Fe2DOC_in, chlabs_in, &
     669         180 :                  alpha2max_low_in, beta2max_in, mu_max_in, fr_graze_in, mort_pre_in, &
     670          90 :                  mort_Tdep_in, k_exude_in, K_Nit_in, K_Am_in, K_sil_in, K_Fe_in, &
     671          90 :                  f_don_in, kn_bac_in, f_don_Am_in, f_doc_in, f_exude_in, k_bac_in, &
     672         180 :                  grow_Tdep_in, zbgc_frac_init_in, &
     673          90 :                  zbgc_init_frac_in, tau_ret_in, tau_rel_in, bgc_tracer_type_in, &
     674             :                  fr_resp_in, algal_vel_in, R_dFe2dust_in, dustFe_sol_in, T_max_in, &
     675             :                  op_dep_min_in, fr_graze_s_in, fr_graze_e_in, fr_mort2min_in, fr_dFe_in, &
     676             :                  k_nitrif_in, t_iron_conv_in, max_loss_in, max_dfe_doc1_in, &
     677             :                  fr_resp_s_in, y_sk_DMS_in, t_sk_conv_in, t_sk_ox_in, fsal_in)
     678             : 
     679             :       real (kind=dbl_kind), optional :: R_C2N_in(:)        ! algal C to N (mole/mole)
     680             :       real (kind=dbl_kind), optional :: R_chl2N_in(:)      ! 3 algal chlorophyll to N (mg/mmol)
     681             :       real (kind=dbl_kind), optional :: F_abs_chl_in(:)    ! to scale absorption in Dedd
     682             :       real (kind=dbl_kind), optional :: R_C2N_DON_in(:)    ! increase compare to algal R_Fe2C
     683             :       real (kind=dbl_kind), optional :: R_Si2N_in(:)       ! algal Sil to N (mole/mole) 
     684             :       real (kind=dbl_kind), optional :: R_S2N_in(:)        ! algal S to N (mole/mole)
     685             :       real (kind=dbl_kind), optional :: R_Fe2C_in(:)       ! algal Fe to carbon (umol/mmol)
     686             :       real (kind=dbl_kind), optional :: R_Fe2N_in(:)       ! algal Fe to N (umol/mmol)
     687             :       real (kind=dbl_kind), optional :: R_Fe2DON_in(:)     ! Fe to N of DON (nmol/umol)
     688             :       real (kind=dbl_kind), optional :: R_Fe2DOC_in(:)     ! Fe to C of DOC (nmol/umol)
     689             : 
     690             :       real (kind=dbl_kind), optional :: fr_resp_in         ! frac of algal growth lost due to respiration
     691             :       real (kind=dbl_kind), optional :: algal_vel_in       ! 0.5 cm/d(m/s) Lavoie 2005  1.5 cm/day
     692             :       real (kind=dbl_kind), optional :: R_dFe2dust_in      ! g/g (3.5% content) Tagliabue 2009
     693             :       real (kind=dbl_kind), optional :: dustFe_sol_in      ! solubility fraction
     694             :       real (kind=dbl_kind), optional :: T_max_in           ! maximum temperature (C)
     695             :       real (kind=dbl_kind), optional :: op_dep_min_in      ! Light attenuates for optical depths exceeding min
     696             :       real (kind=dbl_kind), optional :: fr_graze_s_in      ! fraction of grazing spilled or slopped
     697             :       real (kind=dbl_kind), optional :: fr_graze_e_in      ! fraction of assimilation excreted
     698             :       real (kind=dbl_kind), optional :: fr_mort2min_in     ! fractionation of mortality to Am
     699             :       real (kind=dbl_kind), optional :: fr_dFe_in          ! fraction of remineralized nitrogen
     700             :                                                            ! (in units of algal iron)
     701             :       real (kind=dbl_kind), optional :: k_nitrif_in        ! nitrification rate (1/day)
     702             :       real (kind=dbl_kind), optional :: t_iron_conv_in     ! desorption loss pFe to dFe (day)
     703             :       real (kind=dbl_kind), optional :: max_loss_in        ! restrict uptake to % of remaining value
     704             :       real (kind=dbl_kind), optional :: max_dfe_doc1_in    ! max ratio of dFe to saccharides in the ice (nM Fe/muM C)
     705             :       real (kind=dbl_kind), optional :: fr_resp_s_in       ! DMSPd fraction of respiration loss as DMSPd
     706             :       real (kind=dbl_kind), optional :: y_sk_DMS_in        ! fraction conversion given high yield
     707             :       real (kind=dbl_kind), optional :: t_sk_conv_in       ! Stefels conversion time (d)
     708             :       real (kind=dbl_kind), optional :: t_sk_ox_in         ! DMS oxidation time (d)
     709             :       real (kind=dbl_kind), optional :: fsal_in            ! salinity limitation factor (1)
     710             : 
     711             :       real (kind=dbl_kind), optional :: chlabs_in(:)       ! chla absorption 1/m/(mg/m^3)
     712             :       real (kind=dbl_kind), optional :: alpha2max_low_in(:)  ! light limitation (1/(W/m^2))
     713             :       real (kind=dbl_kind), optional :: beta2max_in(:)     ! light inhibition (1/(W/m^2))
     714             :       real (kind=dbl_kind), optional :: mu_max_in(:)       ! maximum growth rate (1/d)
     715             :       real (kind=dbl_kind), optional :: grow_Tdep_in(:)    ! T dependence of growth (1/C)
     716             :       real (kind=dbl_kind), optional :: fr_graze_in(:)     ! fraction of algae grazed
     717             :       real (kind=dbl_kind), optional :: mort_pre_in(:)     ! mortality (1/day)
     718             :       real (kind=dbl_kind), optional :: mort_Tdep_in(:)    ! T dependence of mortality (1/C)
     719             :       real (kind=dbl_kind), optional :: k_exude_in(:)      ! algal carbon  exudation rate (1/d)
     720             :       real (kind=dbl_kind), optional :: K_Nit_in(:)        ! nitrate half saturation (mmol/m^3) 
     721             :       real (kind=dbl_kind), optional :: K_Am_in(:)         ! ammonium half saturation (mmol/m^3) 
     722             :       real (kind=dbl_kind), optional :: K_Sil_in(:)        ! silicon half saturation (mmol/m^3)
     723             :       real (kind=dbl_kind), optional :: K_Fe_in(:)         ! iron half saturation  or micromol/m^3
     724             :       real (kind=dbl_kind), optional :: f_don_in(:)        ! fraction of spilled grazing to DON
     725             :       real (kind=dbl_kind), optional :: kn_bac_in(:)       ! Bacterial degredation of DON (1/d)
     726             :       real (kind=dbl_kind), optional :: f_don_Am_in(:)     ! fraction of remineralized DON to Am
     727             :       real (kind=dbl_kind), optional :: f_doc_in(:)        ! fraction of mort_N that goes to each doc pool
     728             :       real (kind=dbl_kind), optional :: f_exude_in(:)      ! fraction of exuded carbon to each DOC pool
     729             :       real (kind=dbl_kind), optional :: k_bac_in(:)        ! Bacterial degredation of DOC (1/d)    
     730             : 
     731             :       real (kind=dbl_kind), optional :: zbgc_frac_init_in(:)  ! initializes mobile fraction
     732             :       real (kind=dbl_kind), optional :: bgc_tracer_type_in(:) ! described tracer in mobile or stationary phases
     733             :       real (kind=dbl_kind), optional :: zbgc_init_frac_in(:)  ! fraction of ocean tracer  concentration in new ice
     734             :       real (kind=dbl_kind), optional :: tau_ret_in(:)         ! retention timescale  (s), mobile to stationary phase
     735             :       real (kind=dbl_kind), optional :: tau_rel_in(:)         ! release timescale    (s), stationary to mobile phase
     736             : 
     737             : !autodocument_end
     738             : 
     739             :       character(len=*),parameter :: subname='(icepack_init_zbgc)'
     740             : 
     741             :       !--------
     742             : 
     743         225 :       if (present(R_C2N_in))     R_C2N(:)     = R_C2N_in(:)
     744         225 :       if (present(R_chl2N_in))   R_chl2N(:)   = R_chl2N_in(:)
     745         225 :       if (present(F_abs_chl_in)) F_abs_chl(:) = F_abs_chl_in(:)
     746         135 :       if (present(R_C2N_DON_in)) R_C2N_DON(:) = R_C2N_DON_in(:)
     747         225 :       if (present(R_Si2N_in))    R_Si2N(:)    = R_Si2N_in(:)
     748         225 :       if (present(R_S2N_in))     R_S2N(:)     = R_S2N_in(:)
     749         225 :       if (present(R_Fe2C_in))    R_Fe2C(:)    = R_Fe2C_in(:)
     750         225 :       if (present(R_Fe2N_in))    R_Fe2N(:)    = R_Fe2N_in(:)
     751         135 :       if (present(R_Fe2DON_in))  R_Fe2DON(:)  = R_Fe2DON_in(:)
     752         225 :       if (present(R_Fe2DOC_in))  R_Fe2DOC(:)  = R_Fe2DOC_in(:)
     753             : 
     754          90 :       if (present(fr_resp_in))      fr_resp      = fr_resp_in
     755          90 :       if (present(algal_vel_in))    algal_vel    = algal_vel_in
     756          90 :       if (present(R_dFe2dust_in))   R_dFe2dust   = R_dFe2dust_in
     757          90 :       if (present(dustFe_sol_in))   dustFe_sol   = dustFe_sol_in
     758          90 :       if (present(T_max_in))        T_max        = T_max_in
     759          90 :       if (present(op_dep_min_in))   op_dep_min   = op_dep_min_in
     760          90 :       if (present(fr_graze_s_in))   fr_graze_s   = fr_graze_s_in
     761          90 :       if (present(fr_graze_e_in))   fr_graze_e   = fr_graze_e_in
     762          90 :       if (present(fr_mort2min_in))  fr_mort2min  = fr_mort2min_in
     763          90 :       if (present(fr_dFe_in))       fr_dFe       = fr_dFe_in
     764          90 :       if (present(k_nitrif_in))     k_nitrif     = k_nitrif_in
     765          90 :       if (present(t_iron_conv_in))  t_iron_conv  = t_iron_conv_in
     766          90 :       if (present(max_loss_in))     max_loss     = max_loss_in
     767          90 :       if (present(max_dfe_doc1_in)) max_dfe_doc1 = max_dfe_doc1_in
     768          90 :       if (present(fr_resp_s_in))    fr_resp_s    = fr_resp_s_in
     769          90 :       if (present(y_sk_DMS_in))     y_sk_DMS     = y_sk_DMS_in
     770          90 :       if (present(t_sk_conv_in))    t_sk_conv    = t_sk_conv_in
     771          90 :       if (present(t_sk_ox_in))      t_sk_ox      = t_sk_ox_in
     772          90 :       if (present(fsal_in))         fsal         = fsal_in
     773             : 
     774         225 :       if (present(chlabs_in))    chlabs(:)    = chlabs_in(:)
     775         225 :       if (present(alpha2max_low_in)) alpha2max_low(:) = alpha2max_low_in(:)
     776         225 :       if (present(beta2max_in))  beta2max(:)  = beta2max_in(:)
     777         225 :       if (present(mu_max_in))    mu_max(:)    = mu_max_in(:)
     778         225 :       if (present(grow_Tdep_in)) grow_Tdep(:) = grow_Tdep_in(:)
     779         225 :       if (present(fr_graze_in))  fr_graze(:)  = fr_graze_in(:)
     780         225 :       if (present(mort_pre_in))  mort_pre(:)  = mort_pre_in(:)
     781         225 :       if (present(mort_Tdep_in)) mort_Tdep(:) = mort_Tdep_in(:)
     782         225 :       if (present(k_exude_in))   k_exude(:)   = k_exude_in(:)
     783         225 :       if (present(K_Nit_in))     K_Nit(:)     = K_Nit_in(:)
     784         225 :       if (present(K_Am_in))      K_Am(:)      = K_Am_in(:)
     785         225 :       if (present(K_Sil_in))     K_Sil(:)     = K_Sil_in(:)
     786         225 :       if (present(K_Fe_in))      K_Fe(:)      = K_Fe_in(:)
     787         135 :       if (present(f_don_in))     f_don(:)     = f_don_in(:)
     788         135 :       if (present(kn_bac_in))    kn_bac(:)    = kn_bac_in(:)
     789         135 :       if (present(f_don_Am_in))  f_don_Am(:)  = f_don_Am_in(:)
     790          90 :       if (present(f_doc_in))     f_doc(:)     = f_doc_in(:)
     791         225 :       if (present(f_exude_in))   f_exude(:)   = f_exude_in(:)
     792         225 :       if (present(k_bac_in))     k_bac(:)     = k_bac_in(:)
     793             : 
     794        1395 :       if (present(zbgc_frac_init_in))  zbgc_frac_init(:)  = zbgc_frac_init_in(:)
     795        1395 :       if (present(bgc_tracer_type_in)) bgc_tracer_type(:) = bgc_tracer_type_in(:)
     796        1395 :       if (present(zbgc_init_frac_in))  zbgc_init_frac(:)  =  zbgc_init_frac_in(:)
     797        1395 :       if (present(tau_ret_in)) tau_ret(:) = tau_ret_in(:)
     798        1395 :       if (present(tau_rel_in)) tau_rel(:) = tau_rel_in(:)
     799             : 
     800             : 
     801          90 :       end subroutine icepack_init_zbgc
     802             : 
     803             : !=======================================================================
     804             : !autodocument_start icepack_biogeochemistry
     805             : !
     806             : 
     807       88848 :       subroutine icepack_biogeochemistry(dt, &
     808             :                            ntrcr, nbtrcr,  &
     809      266544 :                            upNO, upNH, iDi, iki, zfswin, &
     810       88848 :                            zsal_tot, darcy_V, grow_net,  &
     811      177696 :                            PP_net, hbri,dhbr_bot, dhbr_top, Zoo,&
     812      177696 :                            fbio_snoice, fbio_atmice, ocean_bio, &
     813      355392 :                            first_ice, fswpenln, bphi, bTiz, ice_bio_net,  &
     814       88848 :                            snow_bio_net, fswthrun, Rayleigh_criteria, &
     815       88848 :                            sice_rho, fzsal, fzsal_g, &
     816      177696 :                            bgrid, igrid, icgrid, cgrid,  &
     817             :                            nblyr, nilyr, nslyr, n_algae, n_zaero, ncat, &
     818             :                            n_doc, n_dic,  n_don, n_fed, n_fep,  &
     819       88848 :                            meltbn, melttn, congeln, snoicen, &
     820       88848 :                            sst, sss, fsnow, meltsn, &
     821      266544 :                            hin_old, flux_bio, flux_bio_atm, &
     822      266544 :                            aicen_init, vicen_init, aicen, vicen, vsnon, &
     823      177696 :                            aice0, trcrn, vsnon_init, skl_bgc)
     824             : 
     825             :       real (kind=dbl_kind), intent(in) :: &
     826             :          dt      ! time step
     827             : 
     828             :       integer (kind=int_kind), intent(in) :: &
     829             :          ncat, &
     830             :          nilyr, &
     831             :          nslyr, &
     832             :          nblyr, &
     833             :          ntrcr, &
     834             :          nbtrcr, &
     835             :          n_algae, n_zaero, &
     836             :          n_doc, n_dic,  n_don, n_fed, n_fep
     837             : 
     838             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     839             :          bgrid         , &  ! biology nondimensional vertical grid points
     840             :          igrid         , &  ! biology vertical interface points
     841             :          cgrid         , &  ! CICE vertical coordinate   
     842             :          icgrid        , &  ! interface grid for CICE (shortwave variable)
     843             :          ocean_bio     , &  ! contains all the ocean bgc tracer concentrations
     844             :          fbio_snoice   , &  ! fluxes from snow to ice
     845             :          fbio_atmice   , &  ! fluxes from atm to ice
     846             :          dhbr_top      , &  ! brine top change
     847             :          dhbr_bot      , &  ! brine bottom change
     848             :          darcy_V       , &  ! darcy velocity positive up (m/s)
     849             :          hin_old       , &  ! old ice thickness
     850             :          sice_rho      , &  ! avg sea ice density  (kg/m^3) 
     851             :          ice_bio_net   , &  ! depth integrated tracer (mmol/m^2) 
     852             :          snow_bio_net  , &  ! depth integrated snow tracer (mmol/m^2) 
     853             :          flux_bio     ! all bio fluxes to ocean
     854             : 
     855             :       logical (kind=log_kind), dimension (:), intent(inout) :: &
     856             :          first_ice      ! distinguishes ice that disappears (e.g. melts)
     857             :                         ! and reappears (e.g. transport) in a grid cell
     858             :                         ! during a single time step from ice that was
     859             :                         ! there the entire time step (true until ice forms)
     860             : 
     861             :       real (kind=dbl_kind), dimension (:,:), intent(inout) :: &
     862             :          Zoo            , & ! N losses accumulated in timestep (ie. zooplankton/bacteria)
     863             :                             ! mmol/m^3
     864             :          bphi           , & ! porosity of layers    
     865             :          bTiz           , & ! layer temperatures interpolated on bio grid (C)
     866             :          zfswin         , & ! Shortwave flux into layers interpolated on bio grid  (W/m^2)
     867             :          iDi            , & ! igrid Diffusivity (m^2/s)    
     868             :          iki            , & ! Ice permeability (m^2)  
     869             :          trcrn     ! tracers   
     870             : 
     871             :       real (kind=dbl_kind), intent(inout) :: &
     872             :          grow_net       , & ! Specific growth rate (/s) per grid cell
     873             :          PP_net         , & ! Total production (mg C/m^2/s) per grid cell
     874             :          hbri           , & ! brine height, area-averaged for comparison with hi (m)
     875             :          zsal_tot       , & ! Total ice salinity in per grid cell (g/m^2) 
     876             :          fzsal          , & ! Total flux  of salt to ocean at time step for conservation
     877             :          fzsal_g        , & ! Total gravity drainage flux
     878             :          upNO           , & ! nitrate uptake rate (mmol/m^2/d) times aice
     879             :          upNH         ! ammonium uptake rate (mmol/m^2/d) times aice
     880             : 
     881             :       logical (kind=log_kind), intent(inout) :: &
     882             :          Rayleigh_criteria    ! .true. means Ra_c was reached  
     883             : 
     884             :       real (kind=dbl_kind), dimension (:,:), intent(in) :: &
     885             :          fswpenln        ! visible SW entering ice layers (W m-2)
     886             : 
     887             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     888             :          fswthrun    , & ! SW through ice to ocean            (W/m^2)
     889             :          meltsn      , & ! snow melt in category n (m)
     890             :          melttn      , & ! top melt in category n (m)
     891             :          meltbn      , & ! bottom melt in category n (m)
     892             :          congeln     , & ! congelation ice formation in category n (m)
     893             :          snoicen     , & ! snow-ice formation in category n (m)
     894             :          flux_bio_atm, & ! all bio fluxes to ice from atmosphere  
     895             :          aicen_init  , & ! initial ice concentration, for linear ITD
     896             :          vicen_init  , & ! initial ice volume (m), for linear ITD
     897             :          vsnon_init  , & ! initial snow volume (m), for aerosol 
     898             :          aicen , & ! concentration of ice
     899             :          vicen , & ! volume per unit area of ice          (m)
     900             :          vsnon     ! volume per unit area of snow         (m)
     901             : 
     902             :       real (kind=dbl_kind), intent(in) :: &
     903             :          aice0   , & ! open water area fraction
     904             :          sss     , & ! sea surface salinity (ppt)
     905             :          sst     , & ! sea surface temperature (C)
     906             :          fsnow       ! snowfall rate (kg/m^2 s)
     907             : 
     908             :       logical (kind=log_kind), intent(in) :: &
     909             :          skl_bgc       ! if true, solve skeletal biochemistry
     910             : 
     911             : !autodocument_end
     912             : 
     913             :       ! local variables
     914             : 
     915             :       integer (kind=int_kind) :: &
     916             :          n, mm              ! thickness category index
     917             : 
     918             :       real (kind=dbl_kind) :: &
     919       27360 :          hin         , & ! new ice thickness
     920       27360 :          hsn         , & ! snow thickness  (m)
     921       27360 :          hbr_old     , & ! old brine thickness before growh/melt
     922       27360 :          dhice       , & ! change due to sublimation/condensation (m)
     923       27360 :          kavg        , & ! average ice permeability (m^2)
     924       27360 :          bphi_o      , & ! surface ice porosity 
     925       27360 :          hbrin       , & ! brine height
     926       27360 :          dh_direct       ! surface flooding or runoff
     927             : 
     928             :       real (kind=dbl_kind), dimension (nblyr+2) :: &
     929             :       ! Defined on Bio Grid points
     930      371808 :          bSin        , & ! salinity on the bio grid  (ppt)
     931      371808 :          brine_sal   , & ! brine salinity (ppt)
     932      426528 :          brine_rho       ! brine_density (kg/m^3)
     933             : 
     934             :       real (kind=dbl_kind), dimension (nblyr+1) :: &
     935             :       ! Defined on Bio Grid interfaces
     936      344448 :          iphin       , & ! porosity 
     937      344448 :          ibrine_sal  , & ! brine salinity  (ppt)
     938      344448 :          ibrine_rho  , & ! brine_density (kg/m^3)
     939      310320 :          iTin            ! Temperature on the interface grid (oC)
     940             : 
     941             :       real (kind=dbl_kind) :: & 
     942       27360 :          sloss            ! brine flux contribution from surface runoff (g/m^2)
     943             : 
     944             :       ! for bgc sk
     945             :       real (kind=dbl_kind) :: & 
     946       27360 :          dh_bot_chl  , & ! Chlorophyll may or may not flush
     947       27360 :          dh_top_chl  , & ! Chlorophyll may or may not flush
     948       27360 :          darcy_V_chl     
     949             : 
     950             :       character(len=*),parameter :: subname='(icepack_biogeochemistry)'
     951             : 
     952             : 
     953      533088 :       do n = 1, ncat
     954             : 
     955             :       !-----------------------------------------------------------------
     956             :       ! initialize
     957             :       !-----------------------------------------------------------------
     958      444240 :          hin_old(n) = c0
     959      444240 :          if (aicen_init(n) > puny) then 
     960       56339 :             hin_old(n) = vicen_init(n) &
     961      207881 :                                 / aicen_init(n)
     962             :          else
     963      236359 :             first_ice(n) = .true.
     964      236359 :             if (tr_brine) trcrn(nt_fbri,n) = c1
     965      236359 :             if (z_tracers) then
     966     4654064 :                do mm = 1,nbtrcr
     967     4654064 :                   trcrn(nt_zbgc_frac-1+mm,n) = zbgc_frac_init(mm)
     968             :                enddo
     969             :             endif
     970      236359 :             if (n == 1) Rayleigh_criteria = .false.
     971      236359 :             if (solve_zsal) trcrn(nt_bgc_S:nt_bgc_S+nblyr-1,n) = c0
     972             :          endif
     973             : 
     974      533088 :          if (aicen(n) > puny) then
     975             :           
     976      207881 :             dh_top_chl = c0
     977      207881 :             dh_bot_chl = c0
     978      207881 :             darcy_V_chl= c0
     979     2029244 :             bSin(:)    = c0
     980      207881 :             hsn        = c0
     981      207881 :             hin        = c0
     982      207881 :             hbrin      = c0
     983      207881 :             kavg       = c0
     984      207881 :             bphi_o     = c0
     985      207881 :             sloss      = c0
     986             :   
     987             :       !-----------------------------------------------------------------
     988             :       ! brine dynamics
     989             :       !-----------------------------------------------------------------
     990             : 
     991      207881 :             dhbr_top(n) = c0
     992      207881 :             dhbr_bot(n) = c0
     993             : 
     994      207881 :             if (tr_brine) then 
     995      207881 :                if (trcrn(nt_fbri,n) .le. c0) trcrn(nt_fbri,n) = c1
     996             : 
     997      207881 :                dhice = c0
     998       56339 :                call preflushing_changes  (aicen  (n),   &
     999       56339 :                                  vicen   (n), vsnon  (n),   &
    1000       56339 :                                  meltbn  (n), melttn (n),   &
    1001       56339 :                                  congeln (n), snoicen(n),   &
    1002       56339 :                                  hin_old (n), dhice,        & 
    1003       56339 :                                  trcrn(nt_fbri,n),          &
    1004       56339 :                                  dhbr_top(n), dhbr_bot(n),  &
    1005             :                                  hbr_old,     hin,          &
    1006      207881 :                                  hsn,         first_ice(n)  )
    1007      207881 :                if (icepack_warnings_aborted(subname)) return
    1008             : 
    1009      207881 :                if (solve_zsal)  then  
    1010             : 
    1011             :                   call compute_microS (n,         nilyr,       nblyr,             &
    1012           0 :                                 bgrid,            cgrid,       igrid,             &
    1013           0 :                                 trcrn(1:ntrcr,n), hin_old(n),  hbr_old,           &
    1014           0 :                                 sss,              sst,         bTiz(:,n),         &
    1015           0 :                                 iTin,             bphi(:,n),   kavg,              &
    1016             :                                 bphi_o,           Rayleigh_criteria, &
    1017           0 :                                 first_ice(n),     bSin,        brine_sal,         &
    1018             :                                 brine_rho,        iphin,       ibrine_rho,        &
    1019           0 :                                 ibrine_sal,       sice_rho(n), sloss)
    1020           0 :                   if (icepack_warnings_aborted(subname)) return
    1021             :                else     
    1022             : 
    1023             :                  ! Requires the average ice permeability = kavg(:)
    1024             :                  ! and the surface ice porosity = zphi_o(:)
    1025             :                  ! computed in "compute_microS" or from "thermosaline_vertical"
    1026             : 
    1027     1821363 :                   iDi(:,n) = c0
    1028             : 
    1029             :                   call compute_microS_mushy (nilyr,         nblyr,       &
    1030           0 :                                    bgrid,         cgrid,         igrid,       &
    1031       56339 :                                    trcrn(:,n),    hin_old(n),    hbr_old,     &
    1032           0 :                                    sss,           sst,           bTiz(:,n),   & 
    1033           0 :                                    iTin(:),       bphi(:,n),     kavg,        &
    1034             :                                    bphi_o,        bSin(:),     &
    1035             :                                    brine_sal(:),  brine_rho(:),  iphin(:),    &
    1036       56339 :                                    ibrine_rho(:), ibrine_sal(:), sice_rho(n), &
    1037      320559 :                                    iDi(:,n)       )
    1038      207881 :                   if (icepack_warnings_aborted(subname)) return
    1039             : 
    1040             :                endif ! solve_zsal  
    1041             : 
    1042       56339 :                call update_hbrine (melttn(n),   &
    1043       56339 :                                    meltsn  (n), dt,          &
    1044             :                                    hin,         hsn,         &
    1045       56339 :                                    hin_old (n), hbrin,       &
    1046             : 
    1047             :                                    hbr_old,     &
    1048       56339 :                                    trcrn(nt_fbri,n),         &
    1049       56339 :                                    dhbr_top(n), dhbr_bot(n), &
    1050             :                                    dh_top_chl,  dh_bot_chl,  & 
    1051             :                                    kavg,        bphi_o,      &
    1052       56339 :                                    darcy_V (n), darcy_V_chl, &  
    1053       56339 :                                    bphi(2,n),   aice0,       &
    1054      207881 :                                    dh_direct)
    1055      207881 :                if (icepack_warnings_aborted(subname)) return
    1056             :                
    1057      207881 :                hbri = hbri + hbrin * aicen(n)  
    1058             : 
    1059      207881 :                if (solve_zsal) then 
    1060             : 
    1061             :                   call zsalinity (n,             dt,                  &
    1062           0 :                                   nilyr,         bgrid,               & 
    1063           0 :                                   cgrid,         igrid,               &
    1064           0 :                                   trcrn(nt_bgc_S:nt_bgc_S+nblyr-1,n), &
    1065           0 :                                   trcrn(nt_qice:nt_qice+nilyr-1,n),   &
    1066           0 :                                   trcrn(nt_sice:nt_sice+nilyr-1,n),   &
    1067           0 :                                   ntrcr,         trcrn(nt_fbri,n),    &
    1068           0 :                                   bSin,          bTiz(:,n),           &
    1069           0 :                                   bphi(:,n),     iphin,               &
    1070           0 :                                   iki(:,n),      hbr_old,             &
    1071             :                                   hbrin,         hin,                 &
    1072           0 :                                   hin_old(n),    iDi(:,n),            &
    1073           0 :                                   darcy_V(n),    brine_sal,           & 
    1074             :                                   brine_rho,     ibrine_sal,          & 
    1075             :                                   ibrine_rho,    dh_direct,           &
    1076             :                                   Rayleigh_criteria,                  &
    1077           0 :                                   first_ice(n),  sss,                 &
    1078           0 :                                   sst,           dhbr_top(n),         &
    1079           0 :                                   dhbr_bot(n),                        &
    1080             :                                   fzsal,         fzsal_g,             &
    1081             :                                   bphi_o,        nblyr,               & 
    1082           0 :                                   vicen(n),      aicen_init(n),       &
    1083           0 :                                   zsal_tot) 
    1084           0 :                   if (icepack_warnings_aborted(subname)) return
    1085             : 
    1086             :                endif  ! solve_zsal
    1087             : 
    1088             :             endif ! tr_brine
    1089             : 
    1090             :       !-----------------------------------------------------------------
    1091             :       ! biogeochemistry
    1092             :       !-----------------------------------------------------------------
    1093             : 
    1094      207881 :             if (z_tracers) then 
    1095             :  
    1096             :                call zbio (dt,                    nblyr,                  &
    1097             :                           nslyr,                 nilyr,                  &
    1098       48078 :                           melttn(n),                                     &
    1099       48078 :                           meltsn(n),             meltbn  (n),            &
    1100       48078 :                           congeln(n),            snoicen(n),             & 
    1101             :                           nbtrcr,                fsnow,                  &
    1102           0 :                           ntrcr,                 trcrn(1:ntrcr,n),       &
    1103       48078 :                           bio_index(1:nbtrcr),   aicen_init(n),          &
    1104       48078 :                           vicen_init(n),         vsnon_init(n),          &
    1105       48078 :                           vicen(n),              vsnon(n),               &
    1106       48078 :                           aicen(n),              flux_bio_atm(1:nbtrcr), &
    1107             :                           n,                     n_algae,                &
    1108             :                           n_doc,                 n_dic,                  &
    1109             :                           n_don,                                         &
    1110             :                           n_fed,                 n_fep,                  &
    1111       48078 :                           n_zaero,               first_ice(n),           &
    1112       48078 :                           hin_old(n),            ocean_bio(1:nbtrcr),    &
    1113           0 :                           bphi(:,n),             iphin,                  &     
    1114           0 :                           iDi(:,n),  &
    1115           0 :                           fswpenln(:,n),                                 &
    1116       48078 :                           dhbr_top(n),           dhbr_bot(n),            &
    1117           0 :                           zfswin(:,n),                                   &
    1118             :                           hbrin,                 hbr_old,                &
    1119             : !                         darcy_V(n),            darcy_V_chl,            &
    1120       48078 :                           darcy_V(n),  &
    1121           0 :                           bgrid,   &
    1122           0 :                           igrid,                 icgrid,                 &
    1123             :                           bphi_o,                                        &
    1124             :                           iTin,                   &
    1125           0 :                           Zoo(:,n),                                      &
    1126           0 :                           flux_bio(1:nbtrcr),    dh_direct,              &
    1127             :                           upNO,                  upNH,                   &
    1128           0 :                           fbio_snoice,           fbio_atmice,            &
    1129           0 :                           PP_net,                ice_bio_net (1:nbtrcr), &
    1130      391932 :                           snow_bio_net(1:nbtrcr),grow_net                )
    1131      199620 :                if (icepack_warnings_aborted(subname)) return
    1132             :      
    1133        8261 :             elseif (skl_bgc) then
    1134             : 
    1135             :                call sklbio (dt,                      ntrcr,               &
    1136             :                             nbtrcr,                  n_algae,             &
    1137             :                             n_doc,               &
    1138             :                             n_dic,                   n_don,               &
    1139             :                             n_fed,                   n_fep,               &
    1140           0 :                             flux_bio (1:nbtrcr),     ocean_bio(1:nbtrcr), &
    1141        8261 :                             aicen    (n),        &
    1142        8261 :                             meltbn   (n),            congeln  (n),        &
    1143        8261 :                             fswthrun (n),            first_ice(n),        &
    1144           0 :                             trcrn    (1:ntrcr,n), &
    1145             :                             PP_net,                  upNO,                &
    1146       16522 :                             upNH,                    grow_net             )
    1147        8261 :                if (icepack_warnings_aborted(subname)) return
    1148             : 
    1149             :             endif  ! skl_bgc
    1150             : 
    1151      207881 :             first_ice(n) = .false.
    1152             : 
    1153             :          endif             ! aicen > puny
    1154             :       enddo                ! ncat
    1155             : 
    1156             :       end subroutine icepack_biogeochemistry
    1157             : 
    1158             : !=======================================================================
    1159             : !autodocument_start icepack_load_ocean_bio_array
    1160             : ! basic initialization for ocean_bio_all
    1161             : 
    1162       88868 :       subroutine icepack_load_ocean_bio_array(max_nbtrcr, &
    1163             :           max_algae, max_don, max_doc, max_dic, max_aero, max_fe, &
    1164       88868 :           nit, amm, sil, dmsp, dms, algalN, &
    1165      177736 :           doc, don, dic, fed, fep, zaeros, ocean_bio_all, hum)
    1166             : 
    1167             :       integer (kind=int_kind), intent(in) :: &
    1168             :          max_algae   , & ! maximum number of algal types 
    1169             :          max_dic     , & ! maximum number of dissolved inorganic carbon types 
    1170             :          max_doc     , & ! maximum number of dissolved organic carbon types
    1171             :          max_don     , & ! maximum number of dissolved organic nitrogen types
    1172             :          max_fe      , & ! maximum number of iron types
    1173             :          max_aero    , & ! maximum number of aerosols 
    1174             :          max_nbtrcr      ! maximum number of bio tracers
    1175             : 
    1176             :       real (kind=dbl_kind), intent(in) :: &
    1177             :          nit         , & ! ocean nitrate (mmol/m^3)          
    1178             :          amm         , & ! ammonia/um (mmol/m^3)
    1179             :          sil         , & ! silicate (mmol/m^3)
    1180             :          dmsp        , & ! dmsp (mmol/m^3)
    1181             :          dms         , & ! dms (mmol/m^3)
    1182             :          hum             ! humic material (mmol/m^3)
    1183             : 
    1184             :       real (kind=dbl_kind), dimension (max_algae), intent(in) :: &
    1185             :          algalN          ! ocean algal nitrogen (mmol/m^3) (diatoms, phaeo, pico)
    1186             : 
    1187             :       real (kind=dbl_kind), dimension (max_doc), intent(in) :: &
    1188             :          doc             ! ocean doc (mmol/m^3)  (proteins, EPS, lipid)
    1189             : 
    1190             :       real (kind=dbl_kind), dimension (max_don), intent(in) :: &
    1191             :          don             ! ocean don (mmol/m^3) 
    1192             : 
    1193             :       real (kind=dbl_kind), dimension (max_dic), intent(in) :: &
    1194             :          dic             ! ocean dic (mmol/m^3) 
    1195             : 
    1196             :       real (kind=dbl_kind), dimension (max_fe), intent(in) :: &
    1197             :          fed, fep        ! ocean disolved and particulate fe (nM) 
    1198             : 
    1199             :       real (kind=dbl_kind), dimension (max_aero), intent(in) :: &
    1200             :          zaeros          ! ocean aerosols (mmol/m^3) 
    1201             : 
    1202             :       real (kind=dbl_kind), dimension (max_nbtrcr), intent(inout) :: &
    1203             :          ocean_bio_all   ! fixed order, all values even for tracers false
    1204             : 
    1205             : !autodocument_end
    1206             : 
    1207             :       ! local variables
    1208             : 
    1209             :       integer (kind=int_kind) :: &
    1210             :          k, ks           ! tracer indices
    1211             : 
    1212             :       character(len=*),parameter :: subname='(icepack_load_ocean_bio_array)'
    1213             : 
    1214     2666040 :       ocean_bio_all(:) = c0
    1215             : 
    1216      355472 :       do k = 1, max_algae           
    1217      266604 :          ocean_bio_all(k)      = algalN(k)           ! N
    1218      266604 :          ks = max_algae + max_doc + max_dic + 1
    1219      355472 :          ocean_bio_all(ks + k) = R_chl2N(k)*algalN(k)!chl
    1220             :       enddo   
    1221             : 
    1222       88868 :       ks = max_algae + 1
    1223      355472 :       do k = 1, max_doc
    1224      355472 :          ocean_bio_all(ks + k) = doc(k)              ! doc
    1225             :       enddo  
    1226       88868 :       ks = ks + max_doc
    1227      177736 :       do k = 1, max_dic
    1228      177736 :          ocean_bio_all(ks + k) = dic(k)              ! dic
    1229             :       enddo 
    1230             : 
    1231       88868 :       ks = 2*max_algae + max_doc + max_dic + 7
    1232      177736 :       do k = 1, max_don
    1233      177736 :          ocean_bio_all(ks + k) = don(k)              ! don
    1234             :       enddo  
    1235             : 
    1236       88868 :       ks = max_algae + 1
    1237       88868 :       ocean_bio_all(ks) = nit                        ! nit
    1238             : 
    1239       88868 :       ks = 2*max_algae + max_doc + 2 + max_dic
    1240       88868 :       ocean_bio_all(ks) = amm                        ! Am
    1241       88868 :       ks = ks + 1
    1242       88868 :       ocean_bio_all(ks) = sil                        ! Sil
    1243       88868 :       ks = ks + 1
    1244       27372 :       ocean_bio_all(ks) =  R_S2N(1)*algalN(1) &      ! DMSPp
    1245       27372 :                         +  R_S2N(2)*algalN(2) &
    1246       88868 :                         +  R_S2N(3)*algalN(3) 
    1247       88868 :       ks = ks + 1
    1248       88868 :       ocean_bio_all(ks) = dmsp                       ! DMSPd
    1249       88868 :       ks = ks + 1
    1250       88868 :       ocean_bio_all(ks) = dms                        ! DMS
    1251       88868 :       ks = ks + 1
    1252       88868 :       ocean_bio_all(ks) = nit                        ! PON
    1253       88868 :       ks = 2*max_algae + max_doc + 7 + max_dic + max_don
    1254      266604 :       do k = 1, max_fe
    1255      266604 :          ocean_bio_all(ks + k) = fed(k)              ! fed
    1256             :       enddo  
    1257       88868 :       ks = ks + max_fe
    1258      266604 :       do k = 1, max_fe
    1259      266604 :          ocean_bio_all(ks + k) = fep(k)              ! fep
    1260             :       enddo  
    1261       88868 :       ks = ks + max_fe
    1262      622076 :       do k = 1, max_aero
    1263      622076 :          ocean_bio_all(ks+k) = zaeros(k)             ! zaero
    1264             :       enddo
    1265       88868 :       ks = ks + max_aero + 1 
    1266       88868 :       ocean_bio_all(ks)  = hum                       ! humics
    1267             : 
    1268       88868 :       end subroutine icepack_load_ocean_bio_array
    1269             : 
    1270             : !=======================================================================
    1271             : !autodocument_start icepack_init_ocean_bio
    1272             : !  Initialize ocean concentration
    1273             : 
    1274          20 :       subroutine icepack_init_ocean_bio (amm, dmsp, dms, algalN, doc, dic, don, &
    1275          20 :              fed, fep, hum, nit, sil, zaeros, max_dic, max_don, max_fe, max_aero,&
    1276          20 :              CToN, CToN_DON)
    1277             : 
    1278             :       integer (kind=int_kind), intent(in) :: &
    1279             :         max_dic, &
    1280             :         max_don, &
    1281             :         max_fe, &
    1282             :         max_aero
    1283             : 
    1284             :       real (kind=dbl_kind), intent(out):: &
    1285             :        amm      , & ! ammonium
    1286             :        dmsp     , & ! DMSPp
    1287             :        dms      , & ! DMS
    1288             :        hum      , & ! humic material
    1289             :        nit      , & ! nitrate
    1290             :        sil          ! silicate
    1291             : 
    1292             :       real (kind=dbl_kind), dimension(:), intent(out):: &
    1293             :        algalN   , & ! algae
    1294             :        doc      , & ! DOC
    1295             :        dic      , & ! DIC
    1296             :        don      , & ! DON
    1297             :        fed      , & ! Dissolved Iron
    1298             :        fep      , & ! Particulate Iron
    1299             :        zaeros       ! BC and dust
    1300             : 
    1301             :       real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
    1302             :        CToN     , & ! carbon to nitrogen ratio for algae
    1303             :        CToN_DON     ! nitrogen to carbon ratio for proteins
    1304             : 
    1305             : !autodocument_end
    1306             : 
    1307             :       ! local variables
    1308             : 
    1309             :       integer (kind=int_kind) :: &
    1310             :         k 
    1311             : 
    1312             :       character(len=*),parameter :: subname='(icepack_init_ocean_bio)'
    1313             : 
    1314          20 :        if (present(CToN)) then
    1315           0 :          CToN(1) = R_C2N(1)
    1316           0 :          CToN(2) = R_C2N(2)    
    1317           0 :          CToN(3) = R_C2N(3)     
    1318             :        endif
    1319             : 
    1320          20 :        if (present(CToN_DON)) then
    1321           0 :          CToN_DON(1) = R_C2N_DON(1)
    1322             :        endif
    1323             : 
    1324          20 :        amm  = c1 ! ISPOL < 1 mmol/m^3 
    1325          20 :        dmsp = p1  
    1326          20 :        dms  = p1    
    1327          20 :        algalN(1) = c1  !0.0026_dbl_kind ! ISPOL, Lannuzel 2013(pennate) 
    1328          20 :        algalN(2) = 0.0057_dbl_kind ! ISPOL, Lannuzel 2013(small plankton)
    1329          20 :        algalN(3) = 0.0027_dbl_kind ! ISPOL, Lannuzel 2013(Phaeocystis)
    1330             :                                      ! 0.024_dbl_kind ! 5% of 1 mgchl/m^3 
    1331          20 :        doc(1) = 16.2_dbl_kind ! 18% saccharides
    1332          20 :        doc(2) = 9.0_dbl_kind  ! lipids
    1333          20 :        doc(3) = c1 ! 
    1334          40 :        do k = 1, max_dic
    1335          40 :             dic(k) = c1
    1336             :        enddo  
    1337          40 :        do k = 1, max_don
    1338          40 :             don(k) = 12.9_dbl_kind              
    1339             :             ! 64.3_dbl_kind ! 72% Total DOC~90 mmolC/m^3  ISPOL with N:C of 0.2
    1340             :        enddo  
    1341             :        !ki = 1
    1342             :        !if (trim(fe_data_type) == 'clim') ki = 2
    1343          60 :        do k = 1, max_fe ! ki, max_fe
    1344          40 :             fed(k) = 0.4_dbl_kind ! c1 (nM) Lannuzel2007 DFe, 
    1345             :                                   ! range 0.14-2.6 (nM) van der Merwe 2011
    1346             :                                   ! Tagliabue 2012 (0.4 nM)
    1347          60 :             fep(k) = c2 ! (nM) van der Merwe 2011
    1348             :                         ! (0.6 to 2.9 nM ocean)
    1349             :        enddo 
    1350          20 :        hum  = c1        ! mmol C/m^3
    1351          20 :        nit  = 12.0_dbl_kind
    1352          20 :        sil  = 25.0_dbl_kind
    1353         140 :        do k = 1, max_aero
    1354         140 :          zaeros(k) = c0
    1355             :        enddo
    1356             :  
    1357             : 
    1358          20 :       end subroutine icepack_init_ocean_bio
    1359             : 
    1360             : !=======================================================================
    1361             : 
    1362             :       end module icepack_zbgc
    1363             : 
    1364             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd